© 2022 Horia Metiu

WorkBookT.19 Calculate ((ΔG)^(0 ))and the equilibrium constant K

WorkBookT.19.1. Use ΔG of formation to calculate the equilibrium constant K

I use as an example for the reaction 
	(Subscript[C, graphite])+ 2(Subscript[H, 2])O(g) = CO(g) + 2(Subscript[H, 2])(g)						(1)
at temperatures of 
	298 K,  500 K,  1000 K  and  2000 K. 
Use data from Appendix 10 which gives ΔG of formation (denoted (Subscript[ΔG, f]^0))

Plan:  I will use the equation 
	- RT lnK(T, p =1bar) =((ΔG)^0)(T, p = 1 bar),		(2)
to calculate the equilibrium constant K at the temperatures requested in the problem. 
To calculate ((ΔG)^0)at T  and p=1 bar, I use: 
	((ΔG)^0)(T, p=1 bar) =(∑_(i=1)^c)ν(i) (Subscript[ΔG, f]^0) (i; T, p=1bar)	(3)
 (Subscript[ΔG, f]^0)(i; T,p) is ΔG of formation of the pure compound i

Remove["Global`*"] ;
Needs["Utilities`CleanSlate`"];
CleanSlate[];

Cell 1:  Input and organize the data.

I need information about four compounds. I make a list of the compounds participating in the reaction
	(Subscript[C, graphite])+ 2(Subscript[H, 2])O(g) = CO(g) + 2(Subscript[H, 2])(g)	
The terms in the list appear in the order they appear in the reaction

compounds={"(Subscript[C, gr])","(Subscript[H, 2])O(g)","(Subscript[CO, 2])(g)","(Subscript[H, 2])(g)"}

Make a list with the stoichiometric coefficients ν(i). The coefficients appear in the same order as the compounds in the list compounds. For example, -2 is the stoichiometric coefficient of (Subscript[H, 2])O(g).

nu={-1,-2,1,2}

Make a list of the temperatures you work at:

TList={298,500,1000,2000}

Input the heat of formation at the temperatures of interest. This is a list of lists. First list has the (Subscript[ΔG, f]^(0 )) of the four compounds (in order they appear in the list compounds) at T = 298K; the second list contains  (Subscript[ΔG, f]^(0 ))for the four compounds (in order they appear in the list compounds) at T = 500K; etc.

ΔGf={
		{0.,-228.582,-394.389,0.(*(Subscript[ΔG, f]^0)(298,1 bar) for graphite, (Subscript[H, 2])O(g),(Subscript[CO, 2])(g) and (Subscript[H, 2])(g)*)}  ,
		{0.,-219.051,-394.939,0.(*(Subscript[ΔG, f]^0)(500,1 bar) for (Subscript[graphiteH, 2])O(g),CO(g) and (Subscript[H, 2])(g)*)}  ,{0.,-192.590,-395.886,0.(*(Subscript[ΔG, f]^0)(1000,1 bar) for graphite,(Subscript[H, 2])O(g),CO(g) and (Subscript[H, 2])(g)*)},
		{0.,-135.528,-396.333,0.(*(Subscript[ΔG, f]^0)(2000,1 bar) for graphite,(Subscript[H, 2])O(g),CO(g) and (Subscript[H, 2])(g)*)}
      } (*KJ/mol*)

I make a Table with the data to check easily that I made the right input.

Prepare a list with the data in a form that can be nicely displayed by TableForm.  Add the temperature to the  list containing the values of Δ(Subscript[G, f]^0).  I do this by using the syntax Prepend[x,a] which adds a in  the first position of the list x.

toTable=MapThread[Prepend ,{ΔGf,TList}]

The table below gives the values of (Subscript[ΔG, f]^0) for the compounds participating in the reaction at the temperatures of interest.  It is curious that (Subscript[ΔG, f]^0) for (Subscript[CO, 2]) is almost temperature independent while (Subscript[ΔG, f]^0) for water changes quite a bit. The reason is that to form water one converts 1.5 mole of gas into one mole and this transformation causes a large change in entropy.  Because ΔG =Δ H - T ΔS the change of entropy is large this causes a large dependence of ΔG on T.

TableForm[ toTable ,TableHeadings->{None,{"T",compounds}//Flatten}]

Cell 2:  Calculate ((ΔG)^0) ( T,p=1bar) at the temperatures in TList

I  evaluate 
	((ΔG)^0) =  (∑_( i=1)^c)ν(i) (Subscript[ΔG, f])(i)  
at various temperatures.This means that I take the scalar product of the list nu with each list in  ΔGf. I get a list of ((ΔG)^0) for the reaction (Subscript[C, graphite])+ 2(Subscript[H, 2])O(g) = CO(g) + 2(Subscript[H, 2])(g)  at the four temperatures for which I have data.

ΔG0= Dot[nu,#]&/@ΔGf (*KJ/mole*)

Mathematica explanation:  This looks so complicated because I want to compute ((ΔG)^0)at all temperatures given  in the list TList.
The function Dot[nu, #]&  calculates the scalar product of nu with the variable #  
The operator /@  Maps the function Dot[nu,#]& into each sublist of the list ΔGf.  The  first sublist contains the ΔG-s of formation at the first temperature in TList, the second sublist contains the ΔG-s of formation at the second temperature in TList, etc.  The output is a list with  ((ΔG)^0)at the temperatures in TList.

The list ΔG0 contains  Δ(G^0) for the reaction at the temperatures 298, 500, 1000, 2000 in kJ/mol.
Make a table with the results

TableForm[MapThread[{#1,#2}&,{TList,ΔG0}],TableHeadings->{None,{"T K","((ΔG)^0) kJ"}}]

Remember now the equation 
	((ΔG)^0) = - R T ln(K) 
and the fact that 
	K = ((((x((Subscript[H, 2])))^2)x(CO))/((x((Subscript[H, 2])O))^2))
These two equations are telling me that at 298 K and 500 K the reaction 
	(Subscript[C, graphite])+ 2(Subscript[H, 2])O(g) = CO(g) + 2(Subscript[H, 2])(g) 
produces practically no CO and no hydrogen. At 1000 and 2000 K the equilibrium favors the formation of the products 
As a general rule:	
 A very negative value for ((ΔG)^0) for a reaction 
 means 
 a large value of the equilibrium constant, 
 which means
 a lot of product is formed

Cell 3:  Calculate the equilibrium constants K at the temperatures in the list

Evaluate
	 ln K = - Δ(G^0)/R T  
Since ((ΔG)^0) is in kJ/mol I need to use R = 0.008314 kJ/ mol K

Calculate ln(K) from
	ln(K) = - ((ΔG)^0)/RT

lnK=MapThread[-#1/(0.008314*#2)&,{ΔG0,TList}]

Calculate K from ln(K)  by using exp[ln[x]] = x

K=Exp[lnK]

Show the results in a table.

TableForm[MapThread[{#1,#2,#3,#4}&,{TList,ΔG0,lnK,K}],
		TableHeadings->{None,{"T K","((ΔG)^0)(T)","ln(K(T)","K(T)"}}]

• The equilibrium constant grows rather rapidly with the temperature because it depends exponentially on ((ΔG)^0). 
• When T=1000 K there is quite a bit of CO and (Subscript[H, 2]) formed (at equilibrium).

WorkBookT.19.2. Temperature dependence of  ((ΔG)^0) and K

Cell 1:  Outline the equations used in the calculations

In the previous  workbook I happened to have data for ΔG of formation at the temperatures for which I wanted to calculate (Subscript[ΔG, 0]) and K.  For some compounds I have heat of formation under standard conditions only (i.e. T = 298.15 K and p = 1 atm). Nevertheless, I want to calculate (Subscript[ΔG, 0]) at other temperatures and p = 1 atm

I will use as example the reaction
		(Subscript[C, gr])(s) + 2(Subscript[H, 2])O(g) ->(Subscript[CO, 2])(g)+2(Subscript[H, 2])(g)					  
 I want to calculate (Subscript[ΔG, 0]) and K at the temperatures
 	T = 298K, 500K, 1000K and 2000K and p=1 atm

• I will use the equation
	((ΔG)^0)(T,p=1) = ((ΔH)^0)(T,p=1) - T ((ΔS)^0)(T,p=1)				  (1)		
• Calculate (Subscript[ΔH, 0])(T, p = 1) from
	 ΔH(T , 1bar) = ΔH(298.15 , 1bar) + (∑_( i=1)^( c))ν(i) (∫_( 298)^T)(Subscript[c, p])(i; T) dT	  (2)
• Calculate ΔH(298.15 , 1bar) from
	ΔH(298.15 , 1bar) = (∑_( i=1)^( c))ν(i)  (Subscript[ΔH, f])(i;298.15, 1bar)			  (3)
where (Subscript[ΔH, f])(i;298.15, 1bar) is the heat of formation of compound i at standard conditions
• Calculate the entropy of reaction from
	ΔS(T, 1bar) = ΔS(298K, 1bar) + (Underoverscript[∑ , (i=1), ( c)])ν(i) (∫_298^T)(((Subscript[c, p])(i; T))/T) dT			  (4)
• Calculate  ΔS(298.15, 1 bar)  from
	ΔS(298K, 1bar) = (∑_( i=1)^c)ν(i) (s^0)(i; 298 K, 1 bar) 				   (5)
where the molar entropies (s^0)(i; 298 K, 1 bar) are given in the NIST WebBook
• Calculate ((ΔG)^0)(T,p=1 ) with
	((ΔG)^0)(T,p=1) = ((ΔH)^0)(T,p=1) - T ((ΔS)^0)(T,p=1)				   (6)
• Calculate the equilibrium constant K from
	  K(T,1bar)  = Exp[ - ((ΔG)^0)(T,1bar)/R T]  					    (7)

Cell 2:  Input the data:

Make a  list named compounds with the compounds in the order they appear in the reaction:

Remove["Global`*"] ;
Needs["Utilities`CleanSlate`"];
CleanSlate[];

compounds={"(Subscript[C, gr])","(Subscript[H, 2])O(g)","(Subscript[CO, 2])(g)","(Subscript[H, 2])(g)"}

Make a list of stoichiometric coefficients:

nu={-1,-2,1,2}

Make a list of heats of formation under standard conditions. From thermodynamic tables I get
	(Subscript[ΔH, f])((Subscript[C, gr]); 298.15 K,1 atm)  = 0.
	(Subscript[ΔH, f])((Subscript[H, 2])O(g); 298.15 K, 1 atm)   = -57.798 Kcal/mol
	(Subscript[ΔH, f])((Subscript[CO, 2])(g); 298.15K,1 atm)   = -94.051 Kcal/mol
	(Subscript[ΔH, f])((Subscript[H, 2])(g); 298.15 K, 1 atm)  =0

ΔHf={0,-57.798, -94.051,0} (*kcal/mol*);

Convert the heats of formation to kJ:

ΔHf=ΔHf*4.184 (*kJ/mol*)

The formation of (Subscript[H, 2])O(g) and (Subscript[CO, 2])(g) is extremely exothermic

Input  the equations for heat capacity dependence on temperature.

(Subscript[c, p])((Subscript[C, gr]); T) = 2.673 + 2.617 10-3 T - 0.1169*(10^6)(T^(-2)) ,
(Subscript[c, p])((Subscript[H, 2])O(g);T) =  7.219 + 2.374 (10^(-3)) T  +0.267 (10^(-6)) (T^2) 
(Subscript[c, p])((Subscript[CO, 2])(g);T) =  5.152 + 15.224 (10^(-3)) T - 9.681 (10^(-6)) (T^2) +2.313 (10^(-9)) T3, 
(Subscript[c, p])((Subscript[H, 2])(g);T) = 6.947 - 0.2 (10^(-3)) T + 0.4808 (10^(-6)) (T^2)

Clear [T]
cpList={2.673+2.617*10^-3 T-0.1169*10^6 T^-2,
7.219+2.374*10^-3 T+0.267*10^-6 T^2,
5.152+15.224*10^-3 T -9.681*10^-6 T^2+2.313*10^-9 T^3, 
6.947-0.2*10^-3 T+0.4808*10^-6 T^2};  (*cal/mol K*)

Check that (Subscript[c, p]) for water gas has a reasonable numerical value ( in cal/mole K)

cpList[[2]]/.T->500 (*cal/mole K*)

This is of the same order of  magnitude as the values given in Appendix 9.

TableForm[Transpose[{compounds, cpList}],
TableHeadings->{None,{"compounds", "(Subscript[c, p]^0) cal/mol K"}}]

The heat capacities for the four compounds participating in the reaction were given in cal/mole.  Convert them to kiloJoule/mole K

cpList=cpList*4.184/1000//Expand;(*kJ/mole K*)

cpList[[2]]/.T->500

I need the entropies of the compounds at T = 298.15 K and p = 1 bar.  I get them from the table given in Appendix 9:

I need the entropies of all the compounds at T =298.15 K and p = 1 bar
 I get them from the thermodynamic table given in Appendix 9.
	(s^0)((Subscript[C, gr]); 298K, 1bar) =  5.74 J/mol K
	(s^0)((Subscript[H, 2])O(g); 298K, 1bar)  =  188.826  
	(s^0)((Subscript[CO, 2])(g); 298 K, 1 bar)  =  213.74  J/mol K  
          (s^0)((Subscript[H, 2])(g); 298 K, 1bar)  =  130.684 J/mol K
Make a list with these values converted to kJ/mol K:

S0List={5.74,188.826,213.74,130.684}/1000 (*kJ/mol K*)

Cell 3:  Calculate the heat of reaction ΔH(T, 1 bar)

ΔH(T K, 1bar) = ΔH(298 K, 1bar) + (∑_( i=1)^c)ν(i) (∫_( 298)^T)(Subscript[c, p])(i; T) dT	
with
 	ΔH(298 K, 1bar) = (∑_(i=1)^c)ν(i) (Subscript[ΔH, f]^0)(298 K, 1 bar)

Calculate  ΔH(298 K, 1bar)

ΔH298=nu.ΔHf (*kJ/mol*)

TList={298.15,500,1000,2000};

Calculate the heat of reaction ΔHT at the temperatures contained in TList using	
 	ΔH(T K, 1bar) = ΔH(298 K, 1bar) + (∑_( i=1)^c)ν(i) (∫_( 298)^( T)) (Subscript[c, p])(i; T) dT

ΔHT=ΔH298 + Integrate[nu.cpList,{T,298,#}]&/@TList  (*kJ*)

Display the result

TableForm[Transpose[{TList,ΔHT}],
TableHeadings-> {None,{"T K","((ΔH)^0)(T,p=1 atm) kJ"}}]

The heat of the reaction 
	(Subscript[C, gr])(s) + 2(Subscript[H, 2])O(g) ->(Subscript[CO, 2])(g)+2(Subscript[H, 2])(g)
does not change much with the temperature.  This is because the number of moles of reactant is equal to the number of moles of product and
	(Subscript[Δc, p]) = (∑_(α=1)^c)ν(α) (Subscript[c, p])(α)
is small which means that the integral in 
	ΔH(T K, 1bar) = ΔH(298 K, 1bar)+ (∑_( i=1)^c)ν(i) (∫_298^('T))(Subscript[c, p])(i; T) dT
is small.

Cell 4: Calculate the entropy of reaction ΔS(T, 1 bar)

Cell 4: Calculate the entropy of reaction at the temperatures in TList

The entropy of reaction at temperature T is 
	ΔS(T, 1atm) =((ΔS)^0)(298,1) +(∫_298^T)(∑_( i=1)^( 4))((ν(i)(Subscript[c, p])(i;T,1))/T)dT
with
 	ΔS(298K, 1bar) = (∑_( i=1)^c)ν(i) (s^0)(i; 298 K, 1 bar)

ΔS298=(nu.S0List) (*kJ/mol K*)

Evaluate
	ΔS(T,1atm) =((ΔS)^0)(298,1) +(∫_298^T)(∑_( i=1)^( 4))((ν(i)(Subscript[c, p])(i;T,1))/T)dT
for the four temperatures given in TList

ΔST=ΔS298+( Integrate[nu.cpList/T,{T,298,#},GenerateConditions->False])&/@TList (*kJ/mole K*)

The Mathematica syntax is the same as the one used to calculate the heats of reaction

TableForm[Transpose[{TList,ΔST}],
 TableHeadings->{None,{"T K","ΔS kj/mol K"}}]

Cell 5:  Calculate ((ΔG)^0)(T, 1 bar)

Evaluate
	 ((ΔG)^0)=((ΔH)^0) -T ((ΔS)^0)
at the four temperatures in TList

ΔGT= MapThread[(#1 -#2*#3)&,{ΔHT,TList,ΔST}]

Display the results

toTable=Transpose[{TList,ΔHT,ΔST,ΔGT}]

TableForm[toTable,TableHeadings->{None,{"T K","((ΔH)^0)(T,1bar) kJ","((ΔS)^0)(T,1bar) kJ/K","((ΔG)^0)(T,1bar) kJ"}}]

In the WorkBookT19.1  I obtained (by a different method)
		({{"T K", "ΔG kJ"}, {298, 62.77499999999998`}, {500, 43.162999999999954`}, {1000, -10.706000000000017`}, {2000, -125.27700000000004`}})

If the data was perfect the two calculations would have given identical results.

Cell 6. Calculate the equilibrium constant at the temperatures in the list TList.

Make a list with the equilibrium constants at the temperatures in TList by using
	((ΔG)^0) = - R T ln(K) 
which gives
	K = exp[- ((ΔG)^0)/ R T]

EquilConstList=MapThread[Exp[-#1/(0.008314 #2)]&,{ΔGT,TList}]

Display T, ((ΔG)^0) and K

TableForm[Transpose[{TList,ΔGT,EquilConstList}],TableHeadings->{None,{"T K","((ΔG)^0)(T) kJ","K"}}]

In WorkBookT. 19.1 I obtained (by a different method):
	({{"T K", "\!\(ΔG\^0\)(T)", "ln(K(T)", "K(T)"}, {298, 62.77499999999998`, -25.33730603994555`, 9.91169205535952`*^-12}, {500, 43.162999999999954`, -10.383209044984353`, 0.000030947787582715855`}, {1000, -10.706000000000017`, 1.2877074813567497`, 3.6244678638682513`}, {2000, -125.27700000000004`, 7.5340991099350525`, 1870.7582518489032`}}),

The two methods give slightly different values for the equilibrium constants. If both data sets used in these calculations were perfect the results of these calculations would have been identical. However the difference is largest at 2000 K but since both values are very large it makes no difference: both calculations predict that at equilibrium the system consist in nearly 100% products.

WorkBook T19.3. Calculate the equilibrium constant at different pressures

Cell 1. The method of calculation

Cell 1.1.  Specify the problem

I will perform the calculation for the reaction
	(Subscript[C, gr])(s) + 2(Subscript[H, 2])O(g) ->(Subscript[CO, 2])(g)+2(Subscript[H, 2])(g)					 	(1)	
at
	 T = 1000 K  
 and the pressures
	p = 1, 10, 50, 100 atm.

I will use the equation
	((ΔG)^0)(1000K, p)= ((ΔH)^0)(1000K,p)  -1000 ((ΔS)^0)(1000K,p)	        (2)
This requires me to calculate the change of enthalpy and of the entropy of reaction when T is constant and p changes.

To calculate the change in the reaction enthalpy or  reaction entropy in the transformation
	{T = 298 K, p = 1 atm} -> {T = 1000 K, p}
I need to pick a path.  I pick one having two segments:
Path 1:  {T = 298 K, p = 1 atm} -> {T = 1000 K, p = 1}
Path 2:  {T = 1000 K, p = 1}->  {T = 1000 K, p}
The pressure p in the last state is one of the values 
	p = 1, 10, 50,100 bar.

Cell 1.2. The change in the reaction enthalpy

Path 1:  To calculate how the enthalpy of reaction changes along the path
	{T = 298 K, p = 1 atm} -> {T = 1000 K, p = 1}
I use
((ΔH)^0)(1000 K,1 bar) = ((ΔH)^0)(298 K, 1bar) + (∑_( i=1)^c)ν(i) (∫_298^('1000))(Subscript[c, p]^0)(i; T) dT    (3)
I have calculated ((ΔH)^0)(1000 K,1 bar) in the previous WorkBook

Path 2: The change of the enthalpy of reaction along the path
	{T = 1000 K, p = 1}->  {T = 1000 K, p}
I use
	((ΔH)^0)(1000 K, p bar)=((ΔH)^0)(1000 K, 1bar) + (∑_( i=1)^c)ν(i)  (Subscript[I, H])(i; 1000,p) 
with 
	(Subscript[I, H])(i; 1000K,p) = (∫_( (Subscript[v, 1]))^( (Subscript[v, 2])))(Subscript[{T (Subscript[(((∂p)/(∂T))), v])  + v (Subscript[(((∂p)/(∂v))), T])}, (Subscript[A, i])])dv	          (4) 
• The notation  (Subscript[{...}, (Subscript[A, i])]) reminds me that when evaluating (Subscript[I, H])(i;1000K,p) for compound i I use the equation of state for that compound.  
• After evaluating the derivatives present in the integrand, I set T = 1000 and perform the integrals.   
•  (Subscript[v, 1]) is the molar volume at pressure 1 atm and T = 1000 
• (Subscript[v, 2]) is the molar volume at T = 1000 K and pressure p.
• (Subscript[v, 1]) and (Subscript[v, 2])  are calculated as usual from the equation of state.

• The only input needed for the calculation of ((ΔH)^0)(1000 K, p ) are the equations of state for all the gases participating in the reaction.  
• The properties of solid graphite change so little with pressure that its enthalpy can be assumed to be independent of pressure.   
• I will use the Benedict-Webb-Rubin equation of state with the parameters given in Appendix 6

Cell 1.3. The change in the reaction entropy

Path 1:  To calculate how the entropy of reaction changes along the path
	{T = 298 K, p = 1 atm} -> {T = 1000 K, p = 1}
I use
((ΔS)^0)(1000K, p = 1) = ((ΔS)^0)(298K, 1atm) + (∑_( i=1)^( c))ν(i) (∫_( 298)^( 1000))(((Subscript[c, p]^0)(i; T))/T) dT	  (3)
I have calculated ((ΔS)^0)(1000 K,1 bar) in the previous WorkBook

Path 2: The change of the entropy of reaction along the path
	{T = 1000 K, p = 1,(Subscript[v, 1])}->  {T = 1000 K, p,(Subscript[v, 2])}
I use
	((ΔS)^0)(1000K, p) = ((ΔS)^0)(1000 K, 1atm)   + (∑_( i=1)^( c))ν(i) (Subscript[I, S])(i; 1000K, p)   
with 
	(Subscript[I, S])(i; 1000K,p)=(∫_( (Subscript[v, 1]))^( (Subscript[v, 2])))(Subscript[{(Subscript[(((∂p)/(∂T))), v])}, (Subscript[A, i])])dv          (4) 
• The notation  (Subscript[{...}, (Subscript[A, i])]) reminds me that when evaluating (Subscript[I, S])(i;1000K,p) for compound i I use the equation of state for that compound.  
• After evaluating the derivatives present in the integrand, I set T = 1000 and perform the integrals.   
•  (Subscript[v, 1]) is the molar volume at pressure 1 atm and T = 1000 
• (Subscript[v, 2]) is the molar volume at T = 1000 K and pressure p.
• (Subscript[v, 1]) and (Subscript[v, 2])  are calculated as usual from the equation of state.

• The only input needed for the calculation of ((ΔS)^0)(1000 K, p bar) are the equations of state for all the gases participating in the reaction.  
• The properties of solid graphite change so little with pressure that its enthalpy can be assumed to be independent of pressure.   
• I will use the Benedict-Webb-Rubin equation of state with the parameters given in Appendix 6

Cell 1.4. The calculation of (Subscript[ΔG, 0])(1000, p) and of K(1000,p)

I use 
	((ΔG)^0)(1000K, p)= ((ΔH)^0)(1000K,p)  +1000 ((ΔS)^0)(1000K,p)		(2)
and 
	K(1000K, p) = Exp {-((((ΔG)^0)(1000,p))/(1000 R))}						(15)

Cell 2. Input data

In WorkBookT. 19.2 I have evaluated  ((ΔH)^0)(1000 K, 1bar) and I found that

ΔH1=93.819 (*kJ*)(*this is ((ΔH)^0)(1000 K, 1bar)*)

ΔS1=0.111(*kJ/mol K*)(*this is ((ΔS)^0)(1000 K, 1bar) *)

This means that I already have the change of enthalpy and the entropy of reaction along path1, which is
	{T = 298 K, p = 1 atm} -> {T = 1000 K, p = 1}
I need to calculate how the enthalpy and the entropy of reaction change along path2, which is
	{T = 1000 K, p = 1}->  {T = 1000 K, p}

Make a list with the pressures of the final states

pList={10,50, 100}

Use the Benedict-Webb-Rubin equation of state which is (Eq. 2.11 in texbook) 
	p=R Tρ + (B R T-A-(ρ^2)(Subscript[C, 0])T-2  + (ρ^3)(b R T-a)+ a α (ρ^6)  + 
	(ρ^3)c T-2 (1 + γ (ρ^2)) exp[-γ (ρ^2)]					
with 
	ρ = 1/v

Define pBWR which gives pressure according to Benedict-Webb-Rubin

ClearAll[pBWR]
pBWR[T_,v_,parameters_]:=Module[
{R=0.08205,p1,p2,p3,p4,p5,ρ},
Clear[A,B,C0,a,b,c,α,γ];
 ρ=1/v;
p1=R T ρ;
p2=(B R T-A-C0 T^-2)ρ^2;
p3=(b R T-a)ρ^3;
p4=a α ρ^6;
p5= c T^-2 ρ^3 (1+γ ρ^2) Exp[-γ ρ^2];
(p1+p2+p3+p4+p5)/.parameters]

Input the parameters for pBWR.  For (Subscript[CO, 2]) (see Appendix 6)

paramBWRCO2={A->2.51604,B->4.48842*10^-2,C0->1.474405*10^5,a->1.3681*10^-1,b->4.12381*10^-3,c->1.49183*10^4,α->8.4658*10^-5,  γ->5.253*10^-3}

Let's see whether the function works properly

pBWR[298,24,paramBWRCO2]

The pressure should be close to 1 atm and it is.

For (Subscript[H, 2]) I use the Benedict-Webb-Rubin equation with the parameters given below (Appendix 6)

paramBWRH2={A->9.7319*10^-2,B->1.8041*10^-2,C0->3.8914*10^2,a->-9.2211*10^-3,b->1.7976*10^-4,c->-2.4613*10^2,α->-0.34215*10^-5,  γ->1.89*10^-3}

pBWR[298,24,paramBWRH2]

I don’t have the BWR equation for water so I am going to use the van der Waals equation. This would not work well at pressure of 100 atm.  You cannot do this in your professional life.

Define the function pwdw[T,v,parameters] that gives the pressure

ClearAll[pvdw]
pvdw[T_,v_,parameters_]:=Module[{R=0.08205 (*liter atm/mol*) },
   Clear[a,b];
		(R T/(v-b)-a/v^2)/.parameters]

Define a replacement rule with the parameters for the van der Waals equation for water

paramvdwH2O={a->5.536*0.987 (*(liter^2)atm/(mol^2)*), b->0.03049 (*liter/mmol*)}

Cell 3. Calculate the heat of reaction at p = 10, 50, 100

Cell 3.1. The plan of calculation

• Calculate (Subscript[I, H]) for each gas, for the specified pressures, from 
		(Subscript[I, H])(i; T, p) = (∫_( (Subscript[v , 1]))^( (Subscript[v, 2]))){T (Subscript[(((∂p)/(∂T))), v])  + v (Subscript[(((∂p)/(∂v))), T])}dv			
• Use these values in 
	ΔH(1000 K, p) = 101.035 kJ +  (∑_(i=1)^c)ν(i)(Subscript[I, H])(i; 1000K,p)
to calculate the heat of reaction at T=1000K and the specified pressures.

Since this calculation is for the transformation (path 2)
	{T = 1000 K, p = 1}->  {T = 1000 K, p}
the molar volume (Subscript[v, 1]) is the molar volume at T = 1000 and p = 1 
the molar volume (Subscript[v, 2]) is the molar volume at T = 1000 and p = 10, 50 or 100 atm

I need to calculate the molar volumes (Subscript[v, 1])=v( 1000K, 1 atm) and (Subscript[v, 2]) =v(1000K, p), for the specified values of p.  Use the function VolumeG defined below:

Cell 3.2.  Calculate (Subscript[I, H])(1000, p) for (Subscript[CO, 2])

Calculate the change of molar volume of (Subscript[CO, 2]) in the transformation

volumeG[T_,p_,f_,parameters_]:=
	Module[{R=0.082 (*liter atm/mol*),eq,sol},
Clear[a,b,c,A,B,C0,γ];
	eq=p==f[T,v ,parameters];
		v/.FindRoot[eq,{v,R T/p}] ];

Calculate the molar volumes for T=1000K and the pressures in pList

v2CO2=volumeG[1000,#,pBWR,paramBWRCO2]&/@pList

The volumes go down with the pressure, as they should. The units are liter/mole

Calculate the molar volume of (Subscript[CO, 2]) for T = 1000 K and p = 1 atm

v1CO2=volumeG[1000,1,pBWR,paramBWRCO2]

Calculate the integral  (Subscript[I, H])((Subscript[CO, 2]); 1000K, p) for the specified pressures. Note the  factor 101.311*(10^(-3))that converts liter atm into kJ.

IHCO2List=Module[{integrand1,integrand2,integral},
 integrand1=v D[pBWR[1000,v,paramBWRCO2],v];
integrand2=T D[pBWR[T,v,paramBWRCO2],T]/.T->1000 ;
integral=Integrate[integrand1+integrand2,{v,v1CO2,#},GenerateConditions->False]&/@v2CO2  (*liter atm*);
101.311*10^-3 *integral (* converted liter atm to kJ*)]

Cell 3.3. Calculate (Subscript[I, H])(1000,p) for (Subscript[H, 2])

This calculation is very similar to that for (Subscript[CO, 2])

v2H2=volumeG[1000,#,pBWR,paramBWRH2]&/@pList

v1H2=volumeG[1000,1,pBWR,paramBWRH2]

IHH2List=Module[{int1,int2,int},
 int1=v D[pBWR[1000,v,paramBWRH2],v];
int2=T D[pBWR[T,v,paramBWRH2],T]/.T->1000 ;
int=Integrate[int1+int2,{v,v1H2,#},GenerateConditions->False]&/@v2H2  (*liter atm*);
101.311*10^-3 *int (*kJ*)(* converted liter atm to kJ*)]

Cell 3.4. Calculate (Subscript[I, H])(1000,p)  for (Subscript[H, 2])O

This is similar to the calculations for (Subscript[CO, 2]) and (Subscript[H, 2])  but I use the van der Waals equation rather then Benedict-Webb-Rubin

paramvdwH2O

v2H2O=volumeG[1000,#,pvdw,paramvdwH2O]&/@pList

v1H2O=volumeG[1000,1,pvdw,paramvdwH2O]

IHH2OList=Module[{integrand1,integrand2,integral},
 integrand1=v D[pvdw[1000,v,paramvdwH2O],v];
integrand2=T D[pvdw[T,v,paramvdwH2O],T]/.T->1000 ;
integral=Integrate[integrand1+integrand2,{v,v1H2O,#},GenerateConditions->False]&/@v2H2O (*liter atm*);
101.311*10^-3 *integral (*kJ*)]//Chop(* converted liter atm to kJ*)

Cell 3.5. Display the results for molar volumes and for (Subscript[I, H])

Make lists with the quantities you want to display:

vCO2=Prepend[v2CO2,v1CO2]

vH2=Prepend[v2H2,v1H2]

vH2O=Prepend[v2H2O,v1H2O]

p=Prepend[pList,1]

Organize the lists in the format needed for using TableForm:

toTable=MapThread[{#1,#2,#3,#4}&,{p,vCO2,vH2,vH2O}]

Print["
Molar volumes at different pressures and T=1000 K"]

TableForm[toTable,TableHeadings->{None,{"p atm","(Subscript[v, (Subscript[CO, 2])]) L/mol","(Subscript[v, (Subscript[H, 2])]) L/mol","(Subscript[v, (Subscript[H, 2]) O]) L/mol"}}]

Display the values of (Subscript[I, H])(i; 1000K, p) for the three pressures.  Organize the lists in the format required by TableForm

toTable=MapThread[{#1,#2,#3,#4}&,{pList,IHCO2List,IHH2List,IHH2OList}]

As you can see in the table below the effect of pressure is small.

TableForm[toTable,TableHeadings->{None,{"p atm","(Subscript[I, H])((Subscript[CO, 2]))kJ","(Subscript[I, H])((Subscript[H, 2])) kJ","(Subscript[I, H])((Subscript[H, 2])O) kJ"}}]

Cell 3.6.  Calculate ΔH(1000K, p ).

Use
	ΔH(1000 K, p) = 101.035 kJ +  (∑_(i=1)^c)ν(i)(Subscript[I, H])(i; 1000K,p) 
where 
	101.035 kJ is the reaction enthalpy at 1000 K and p = 1

Define the stoichiometric coefficients (Subscript[ν, i]).  I don’t need the stoichiometric coefficient for graphite because its properties change very little with pressure.  I also don’t have an equation of state for it.

compounds={"(Subscript[CO, 2])","(Subscript[H, 2])","(Subscript[H, 2])O"}

nu={1,2,-2}

Make a lists of the (Subscript[I, H])(i; 1000K, p) which contain the values of (Subscript[I, H]) at the same pressure. Put them in a list in the order 

{(Subscript[I, H])((Subscript[CO, 2]); 1000K, p), (Subscript[I, H])((Subscript[H, 2]); 1000K, p), (Subscript[I, H])((Subscript[H, 2])O; 1000K, p)},

which is the same as the order of the chemical formulae  in the list compounds. Each list (in the list IHList) is (Subscript[I, H]) for one of the compounds

IHList=MapThread[{#1,#2,#3}&,{IHCO2List,IHH2List,IHH2OList}]

Take the scalar product of nu with the elements of this list to get  
	(∑_(i=1)^3)ν(i) (Subscript[I, H])(i; 1000K,p) ; 
add ΔH1 to it to calculate ΔH(1000K, p)

ΔH=ΔH1+nu.#&/@IHList

This list contains the values of the heat of reaction ΔH(1000K, p) for p = 10 atm, 50 atm and 100atm, respectively.  For example  ΔH(1000K, 50 atm) = 94. 90143 kJ.  
 • We have  ΔH(1000, p = 1) = 101.035 kJ.  The change with pressure is small but not negligible. 
 • The equilibrium constant K depends exponentially on ΔH.  Therefore, small changes in ΔH could result in fair changes in K
• Keep in mind that van der Waals equation is very inaccurate at p = 100 atm and this affects the high pressure result.

Cell 4: Calculate the entropy of reaction at p = 10, 50, 100

Cell 4.1.  The equations used  in calculations

Use
	ΔS(1000K, p) = ΔS(1000K, p=1) +  (∑_(i=1)^c)ν(i) (Subscript[I, S])(i; 1000K, p)	
with
	(Subscript[I, S])(i; 1000K, p) =(∫_(Subscript[v, 1])^(Subscript[v, 2]))(Subscript[(((∂p)/(∂T))), v])dv. 					
In the previous WorkBook I calculated that
	ΔS(1000K, p = 1) - ΔS(298, 1)  = 0.111437 kJ/mol K

I use the Benedict-Webb-Rubin equation of state for (Subscript[H, 2]) and (Subscript[CO, 2]) and the van der Waals equation for (Subscript[H, 2])O(g).  These equation of state and the appropriate parameters should be in memory.

I have already calculated the molar volumes when I calculated the change of enthalpy with pressure.  At this point they should be in the memory

Cell 4.2. Calculate (Subscript[I, S]) ((Subscript[CO, 2]); T = 1000 K, p) for p = 10, 50, 100 atm

Use 
	(Subscript[I, S])((Subscript[CO, 2]) ; 1000K, p) = (∫_(Subscript[v, 1])^(Subscript[v, 2]))(Subscript[(((∂p)/(∂T))), v])dv
Since the path 2 is 
	Path 2:  {T = 1000 K, p = 1,(Subscript[v, 1])}->  {T = 1000 K, p,(Subscript[v, 2])}
(Subscript[v, 1]) is the molar volume at T = 1000 K and p = 1 atm
(Subscript[v, 2]) is the molar volume at T = 1000 K and p = 10, 50 or 100 atm

Use Eq. 14 to calculate the integral  (Subscript[I, S])((Subscript[CO, 2]); 1000K, p) for the specified pressures. Note the  factor 101.311*(10^(-3)) that converts liter atm into kJ. Remember that the molar volumes for p = 1 atm are stored in v1CO2 and the ones for the other pressures are in v2CO2.

The output is (Subscript[I, S])((Subscript[CO, 2]); T= 1000 K, p) for p = 10, 50, 100 atm

ISCO2List=Module[{integrand, integral},
integrand= D[pBWR[T,v,paramBWRCO2],T]/.T->1000 ;
integral=Integrate[integrand,{v,v1CO2,#}]&/@v2CO2  (*liter atm*);
101.311*10^-3 *integral (*kJ*)]

Cell 4.4: Calculate (Subscript[I, S ])for (Subscript[H, 2])(g) at the specified pressures

This calculation is very similar to that for (Subscript[CO, 2])

ISH2List=Module[{integrand, integral},
integrand= D[pBWR[T,v,paramBWRH2],T]/.T->1000 ;
integral=Integrate[integrand,{v,v1H2,#}]&/@v2H2  (*liter atm*);
101.311*10^-3 *integral (*kJ*)]

Cell 4.5: Calculate (Subscript[I, S]) for (Subscript[H, 2])O at the specified pressures

This is similar to the calculations for (Subscript[CO, 2]) and (Subscript[H, 2]) but I use the van der Waals equation rather then Benedict-Webb-Rubin

ISH2OList=Module[{integrand, integral},
integrand= D[pvdw[T,v,paramvdwH2O],T]/.T->1000 ;
integral=Integrate[integrand,{v,v1H2O,#}]&/@v2H2O (*liter atm*);
101.311*10^-3 *integral//Chop (*kJ*)]

Cell 4.6: Display the results for the (Subscript[I, S])

Make a lists of the (Subscript[I, S])(i; 1000K, p) which contain the values of (Subscript[I, S]) at the same pressure. Put them in the order 
	{(Subscript[I, S])((Subscript[CO, 2]); 1000K, p), (Subscript[I, S])((Subscript[H, 2]); 1000K, p), (Subscript[I, S])((Subscript[H, 2]); 1000K, p)},
 in which the formulae appear in the list compounds

ISList=MapThread[{#1,#2,#3}&,{ISCO2List,ISH2List,ISH2OList}]

Display the results for (Subscript[I, S ])(1000K, p). Prepare the lists you want to use in TableForm.  I add the pressure to the list of (Subscript[I, S])-s

toTable=MapThread[Prepend[#1,#2]&,{ISList,pList}]

TableForm[toTable,TableHeadings->{None,{"p atm","(Subscript[I, S])((Subscript[CO, 2]))kJ","(Subscript[I, S])((Subscript[H, 2])) kJ","(Subscript[I, S])((Subscript[H, 2])O) kJ"}}]

Cell 4.7: Use these values of (Subscript[I, s])(i; 1000K, p) to calculate ΔS(1000K, p )

Use 
	ΔS(1000 K, p) = 0.111  kJ +  (∑_( i=1)^c)ν(i)(Subscript[I, S])(i; 1000K,p)

Define the stoichiometric coefficients ν(i)

compounds={"(Subscript[CO, 2])","(Subscript[H, 2])","(Subscript[H, 2])O"}

nu={1,2,-2}

Take the scalar product of nu (the list of stoichiometric coefficients) with the elements of the list ISList to get 
	 (∑_( i=1)^( 3))ν(i) (Subscript[I, S])(i; 1000K,p) ; 
add ΔS1 to it to calculate  ΔS(1000K, p)

ΔS=ΔS1+nu.#&/@ISList (*kJ*)

This list contains the values of ((ΔS)^0)(1000K, p) for p = 10 atm, 50 atm and 100atm, respectively.  Thus ((ΔS)^0)(1000K, 50 atm) = 102.118 kJ.  The change of ΔS with pressure is small. We have seen this to be the case when we studied thermochemistry.

Cell 5.  Calculate ((ΔG)^0)  and K for p = 10, 50, 100

Use

	((ΔG)^0)(1000K, p) =((ΔH)^0)(1000K, p) -1000*((ΔS)^0)(1000K, p) 
and
	K=exp[-(((ΔG)^0)/(R T))]

ΔG0p=ΔH-1000*ΔS

Calculate the equilibrium constant at 1000 K and various pressures from 
	K(1000,p)=Exp[-((ΔG)^0)(1000K, p)/ 1000 R]

K =Exp[-ΔG0p/(1000*8.314*10^-3)]

I used the following Mathematica syntax

Exp[-a {1,2,3}/b]

Display the results.  The pressures are in pList, the heats of reaction in ΔH, the entropies in ΔS, the ((ΔG)^0) in ΔG0p and the equilibrium constants in K . Prepare a list that can be used in TableForm.

toTable=MapThread[{#1,#2,#3, #4,#5}&,{pList,ΔH,ΔS,ΔG0p,Kp}]

Display the results for 
	(Subscript[C, gr])(s) + 2(Subscript[H, 2])O(g) ->(Subscript[CO, 2])(g)+2(Subscript[H, 2])(g)					 (1)

TableForm[toTable,TableHeadings->{None,{"p atm","ΔH kJ","ΔS kJ","((ΔG)^0) kJ","K"}}]

The equilibrium constant is fairly sensitive to pressure. As the pressure is increased the amount of products formed decreases (K becomes smaller).  This is in agreement with Le Chatelier principle.