© 2022 Horia Metiu

WorkBookT.20  Calculate the change of K and composition with T and p

WorkBookT.20.1. Find the dependence of K on temperature

Calculate the equilibrium constant for the reaction 
	(Subscript[CH, 4])+(Subscript[H, 2])O  =  CO + 3 (Subscript[H, 2])					
at a temperatures 298 K, 500 K, 1000K and 1500 K and pressure of 1 bar. Assume that the gases form an ideal mixture.

Cell 1. Outline of the method of calculation

• Calculate K(T, 1 bar) from  Eq. 20.7 of the textbook:
	ln K(T, 1bar ) = ln K(298.15, 1bar) + 
	(∫_298.15^T)((ΔH)^0)(T, 1bar)/ (R (T^2))dT				(1)

• To evaluate the right-hand-side of  Eq.1 I need to know 
	((ΔH)^0)(T, 1bar) and K(298.15, 1bar)

• Calculate ((ΔH)^0)(T, 1bar)  from (see Chapter 13,  Eq. 13.21)
((ΔH)^0)(T, 1bar) =((ΔH)^0)(298.15K, 1bar) +(∫_298^T)(Subscript[Δc, p]^0)(T,1bar)dT
where
      (Subscript[Δc, p]^0)(T, 1bar) = (∑_(i=1)^c)ν(i)(Subscript[c, p]^0)(i;T, 1bar)

• To calculate ΔH(298.15 K, 1bar) I use
     ΔH(298K, 1bar) = (∑_(i=1)^c)ν(i) (Subscript[ΔH, f])(i; 298.15 K, 1 bar)		 
where 
	(Subscript[ΔH, f])(i; 298.15K, 1bar) 
is the heat of formation of compound i.

• Calculate  K(298.15K, 1 bar)  from
	- RT ln[K(298.15 K, 1bar)] = ((ΔG)^0)(298.15 K, 1bar)	
with 
((ΔG)^0)(298.15 K,1bar) =(∑_(i=1)^c)ν(i) (Subscript[ΔG, f])(i; 298.15 K, 1bar)    
Here (Subscript[ΔG, f]^( ))(i; 298.15 K, 1bar) is the ΔG of formation of compound i.

• This collection of six equations allows me to calculate 
	K(T, 1bar) at any temperature T and pressure of 1 bar.

Cell 2:  Input data:

I need 
	(Subscript[c, p])(T) for all compounds 
	(Subscript[ΔH, f]^0)(298.15K,1bar) for all compounds
	(Subscript[ΔG, f]^0)(298.15K,1bar) for all compounds.
These data are available in Appendices 8.2 and 9 or in the NIST data base

Off[Remove::rmnsm]
Remove["Global`*"] ;
Needs["Utilities`CleanSlate`"];
CleanSlate[];

Make a list with the compounds participating in the reaction.

compounds={"CH4","H2O","H2","CO"}

Define below a list with the stoichiometric coefficients (the place of the stoichiometric coefficient for a compound i in this list is the same as the place of the compound i in the list compounds; for example, the stoichiometric coefficient of (Subscript[H, 2]) is in the third place in list nu because the symbol (Subscript[H, 2]) is in the third place in the list compounds.  This rule is followed for all other lists).

nu={-1,-1,3,1}

Input Δ(Subscript[G, f]^0), the Gibbs free energy of formation for the compounds at 298.15 K and p = 1 bar; the first value is for (Subscript[CH, 4]), next for (Subscript[H, 2])O etc.  Be careful to use Δ(Subscript[G, f]^0) for gaseous water, not for liquid.

ΔGf= {-50.72,-228.572,0,-137.168}(*kj/mol*)

Input Δ(Subscript[H, f]^0) of formation for the compounds, at 298.15 K and 1 bar; the heat of formation of (Subscript[CH, 4]) is given first, etc. Use the heat of formation of the gaseous water.

ΔHf = {-74.81, -241.818,0,-110.525}(*kJ/mol K*)

Make a list of the temperatures at which you want to perform calculations.

temp={298.15,500,1000,1500}

Input data for the heat capacity (Subscript[c, p]) of the compounds. The first expression in the list cp is (Subscript[c, p]) of (Subscript[CH, 4]), next is (Subscript[c, p]) of (Subscript[H, 2])O(g),...  The order is that of the compounds in the list compounds.  The data is from Appendix 8.2

cpCH4=4.171+14.4500 T 10^-3+0.2670 10^-6-1.722 10^-9 T^3 (*cal/mol K*)

cpH2O= 7.219 +2.3740 10^-3 T+0.2670 10^-6 T^2 (*cal/mol K*)

cpH2=6.947 - 0.2000 10^-3 T +0.4808 10^-6 T^2(*cal/mol K*)

cpCO=6.342+1.836 10^-3 T-0.2801 10^-6 T^2 (*cal/mol K*)

compounds
cp={cpCH4,cpH2O,cpH2,cpCO}(*cal/mol K*)

Display the heat capacities to check that the input is correct

TableForm[Transpose[{compounds,cp }]](*cal/mol K*)

Convert heat capacity to kJ/mol K

cp=cp*4.184/1000 //N(*kJ/mol K*)

Cell 3:  Define a function that gives ΔH(T, 1 bar)

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

Calculate ΔH(298K, 1bar)

ΔH0298 = nu.ΔHf (*kJ*)

The reaction 
	(Subscript[CH, 4])+(Subscript[H, 2])O  =  CO + 3 (Subscript[H, 2])		
is very endothermic, which we should expect because methane and water are very stable molecules.

Calculate (Subscript[Δc, p]) from Eq. 3:
	(Subscript[Δc, p]^0)(T, 1bar) = (∑_( i=1)^( c))ν(i)(Subscript[c, p]^0)(i;T,1bar)

nu.cp//N//Expand (*kJ/mol K*)

Calculate ((ΔH)^0)(T, 1 bar) from Eq. 2:
((ΔH)^0)(T, 1bar) =((ΔH)^0)(298.15K, 1bar) +(∫_298^T)(Subscript[Δc, p]^0)(T,1bar)dT

Clear[z,ΔH0,T]
ΔH0[z_]:=(ΔH0298+Integrate[nu.cp,{T,298,z}]) (*kJ*)

Look at the result:

ΔH0[T]

Calculate ((ΔH)^0) at 1000 K:

ΔH0[1000]

The reaction is very endothermic at high temperature also.

Cell 4:  Calculate ((ΔG)^0)(298.15 K,1 atm)

Use 
	((ΔG)^0)(298.15 K, 1bar) = (∑_( i=1)^( c))ν(i) (Subscript[ΔG, f]^0)(i; 298.15 K, 1bar)
to calculate 
	((ΔG)^0)(298.15 K, 1bar)

ΔG0298 = nu.ΔGf (*kJ*)

This means that no products will be formed at 298 K

Cell 5:  Calculate ln K(298.15 K, 1 bar)

Use
	-RT lnK(298.15 K, 1bar) = ((ΔG)^0)(298.15 K, 1bar)	
to calculate 
	lnK(298.15 K, 1bar) =  -((ΔG)^0)(298.15 K, 1bar)/RT 
Since the unit of ((ΔG)^0) is  kJ, use 
	R = 8.134 (10^(-3))kJ/mol K.

lnK298=-ΔG0298/(8.314*10^-3*298.15)

Cell 6: Calculate ln K(T, 1 bar)

Use the equation  
ln K(T, 1 ) = ln K(298.15, 1) +(∫_( 298.15)^( T))((ΔH)^0)(T, 1)/ (R (T^2))dT 
to calculate the equilibrium constant at a temperature T.

lnK[x_]:=lnK298+Integrate[ΔH0[y]/(R y^2),{y,298.15,x},GenerateConditions->False]

Show the result. Use R = 8.314 (10^(-3))kJ/K mol

lnK[T]/.R->8.314*10^-3//N//Expand

Calculate the equilibrium constant, as a function of T:

Clear[y, K]
 K[x_]:= Module[{R=8.314*10^-3, expr},
expr=lnK298+Integrate[ΔH0[y]/(R y^2),{y,298.15,x},GenerateConditions->False];
Exp[expr] ]

Take a look at it:

K[T]

Calculate a list of equilibrium constants for the temperatures in the list temp

Klist=K/@temp

No product is made at 298 K and 500 K.  A lot of product is made at 1000 K and there is almost no reactant left at T = 1500 K

Cell 7: Display the results of all these calculations

toTable=.
toTable=Transpose[{temp,K/@temp}];
TableForm[toTable,TableHeadings->{None,{"T K","K"}}]

WorkBookT.20.2. Calculate the dependence of K on pressure

Cell 1. Explain the problem and the procedure

I start with the equation
	 (Subscript[(((∂ln[K(T, p)])/(∂p))), T, n]) = -((((ΔV)^0)(T, p))/(R T))
Integrating this leads to 
	ln[K(T,p)]=ln[K(298, 1bar)] = -(∑_( i=1)^( 4))(((ν(i))/(R (Subscript[T, 0]))))(Subscript[I, v])(i; T,p) 
with 
	 (Subscript[I, v])(i; T,p) =(∫_(Subscript[v, 1])^(Subscript[v, 2]))v(i; T,p) (Subscript[(((∂p)/(∂v))), T])dv				 
	 (Subscript[v, 1]) = (v^0)(T,1bar)							 
and 
	(Subscript[v, 2])  = (v^0)(T, p)							 
	
are the molar volumes of the gases participating in the reaction.  They are calculated from the equation of state.   In all calculations T= 1000K

I have already calculated, in WorkBookT.20.1, that 
	K(1000K, 1bar) = 25.93
• All I need is good equations of state for the gases participating in the reaction. 
• For methane, CO and (Subscript[H, 2]) I use the Benedict-Webb-Rubin (BWR) equation. 
• The parameters are given in Appendix 6.  
• For (Subscript[H, 2])O I use the van der Waals equation with the parameters given in Appendix 4.

Cell 2: Data

Off[Remove::rmnsm]
Remove["Global`*"] ;
Needs["Utilities`CleanSlate`"];
CleanSlate[];

Define the equations of state. For (Subscript[H, 2])O(g) I use the van der Waals equation

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

The parameters for water are

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

Define a function 
	pBWR[T, v, parameters] 
which calculates the pressure by using the Benedict-Web-Rubin equation (Eq. 4):
	p = R T ρ + ( B R T-A-(Subscript[C, 0])T-2) ρ2 + (b R T-a) ρ3 + a α ρ6  + c (ρ^3) (T^(-2)) (1 + γ ρ2) exp[-γ (ρ^2)] 
Here ρ = 1/v where v is the molar volume.

Clear[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];
SetAttributes[pBWR,HoldAll];

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};
paramBWRCH4={A->1.79894,B->4.5625*10^-2,C0->3.18382*10^4,a->0.04352,b->2.52033*10^-3,c->3.5878*10^3,α->3.30*10^-4,  γ->1.805*10^-2};
paramBWRCO={A->1.03115,B->0.04,C0->31.124*10^4,a->3.665*10^-2,b->2.6316*10^-3,c->1.04*10^3,α->1.350*10^-4,  γ->0.006};

The list of the pressures for which the equilibrium constant will be calculated

pListBar={10,50,100};(*bar*)

The equations of state I use require the pressure to be in atm. Because of this I create the list pListAtm that give the above pressures in atm

pListAtm=0.987*pListBar

The equilibrium constant at 1000K and p = 1 bar (from WorkBookT.20.1);

K1000=25.93;

Cell 3: Calculate the integrals   (Subscript[I, v])(CO; (Subscript[T, 0]),p)

I need to evaluate
	 (Subscript[I, v])(i; T,p) =(∫_(Subscript[v, 1])^(Subscript[v, 2]))v(i; T,p) (Subscript[(((∂p)/(∂v))), T])dv				 
	 (Subscript[v, 1]) = (v^0)(T,1bar)							 
and 
	(Subscript[v, 2])  = (v^0)(T, p)		
for each reaction participant

Use the function volumeG to get the molar volumes

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}] ];

Use the Benedict-Webb-Rubin equation of state  to calculate  
	v( (Subscript[CO, 2]); 1000K, p) for p =10 atm, 50 atm and 100 atm

v2CO=volumeG[1000,#,pBWR,paramBWRCO]&/@pListAtm

Calculate v( (Subscript[CO, 2]); 1000K, p=1bar ). Remember that 1bar = 0.987 atm and that the pressure in the argument of the function VolumeG is in atmosphere

v1CO=volumeG[1000,0.987 (*atm*),pBWR,paramBWRCO] (*liter/mol*)

Calculate the integral  (Subscript[I, v])((Subscript[CO, 2]); 1000K, p)  by using its definition  
 (Subscript[I, v])(i; (Subscript[T, 0]),p) =(∫_( (Subscript[v, 1]))^( (Subscript[v, 2]))) v(i; (Subscript[T, 0]),p) (Subscript[(((∂p)/(∂v))), T])dv			 
Calculate it for three pressures.  Note the  factor 101.311*(10^(-3)) that converts liter atm into kJ.

IVCOList=Module[{integrand,int2,int},
 integrand=v D[pBWR[1000,v,paramBWRCO],v];
int=Integrate[integrand,{v,v1CO,#}]&/@v2CO;(*liter atm*)
101.311*10^-3 *int (*converted liter atm to kJ*)]

Cell 4:  Calculate (Subscript[I, v])((Subscript[H, 2]); (Subscript[T, 0])=1000,p)

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

Get the final molar volume for (Subscript[H, 2]) for three pressures and T = 1000

Get the final molar volume for (Subscript[H, 2]) for three pressures.

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

Get the initial volume for (Subscript[H, 2]) for three pressures.

v1H2=volumeG[1000,0.987 (*atm*),pBWR,paramBWRH2]

IVH2List=Module[{int1,int2,int},
 int1=v D[pBWR[1000,v,paramBWRH2],v];
int=Integrate[int1,{v,v1H2,#}]&/@v2H2  (*liter atm*);
101.311*10^-3 *int (*kJ*)]

Cell 5:  Calculate (Subscript[I, v])((Subscript[CH, 4]); (Subscript[T, 0])=1000,p)

Calculate 
	 (Subscript[I, v])(i; (Subscript[T, 0]),p) =(∫_(Subscript[v, 1])^(Subscript[v, 2]))v(i; (Subscript[T, 0]),p) (Subscript[(((∂p)/(∂v))), T])dv	
for methane, at the pressures in pListAtm

v2CH4=volumeG[1000,#,pBWR,paramBWRCH4]&/@pListAtm (*liter/mol*)

v1CH4=volumeG[1000,0.987 (*atm*),pBWR,paramBWRCH4] (*liter/mol*)

IVCH4List=Module[{int1,int},
 int1=v D[pBWR[1000,v,paramBWRCH4],v];
int=Integrate[int1,{v,v1CH4,#}]&/@v2CH4  (*liter atm*);
101.311*10^-3 *int (*kJ*)]

Cell 6:  Calculate (Subscript[I, v])((Subscript[H, 2])O; (Subscript[T, 0])=1000,p) for (Subscript[H, 2])O(g)

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

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

v1H2O=volumeG[1000,0.987 (*atm*),pvdw,paramvdwH2O]

IVH2OList=Module[{int1,int2,int},
 int1=v D[pvdw[1000,v,paramvdwH2O],v];
int=Integrate[int1,{v,v1H2O,#}]&/@v2H2O (*liter atm*);
101.311*10^-3 *int (*kJ*)]//Chop

Cell 7:  Display the molar volumes and  (Subscript[I, v])(i; 1000,p)

Add the molar volumes at p=1atm to the beginning of the list containing the volumes at p=10, 50, 100 atm. Use the function Prepend.  Also add the pressure p=1 atm to the list pList containing 10, 50, 100.

vCO=Prepend[v2CO,v1CO];
vH2=Prepend[v2H2,v1H2];
vH2O=Prepend[v2H2O,v1H2O];
vCH4=Prepend[v2CH4,v1CH4];
p=Prepend[pListBar,1];

Make a table with the molar volumes at different pressures:

toTable=MapThread[{#1,#2,#3,#4,#5}&,{p,vCO,vH2,vH2O,vCH4}];TableForm[toTable,TableHeadings->{None,{"p bar","vCO l/mol","(Subscript[vH, 2]) l/mol","(Subscript[vH, 2])O l/mol","(Subscript[vCH, 4]) l/mol"}}]

Display the values of (Subscript[I, v])(i; 1000K, p) for the three pressures.  Make a list in which the (Subscript[I, v]) for the same pressure and different compounds are stored in one sublist.

toTable=MapThread[{#1,#2,#3,#4,#5}&,{pListBar,IVCOList,IVH2List,IVH2OList, IVCH4List}];TableForm[toTable,TableHeadings->{None,{"p bar","(Subscript[I, v])(CO)kJ","(Subscript[I, v])((Subscript[H, 2])) kJ","(Subscript[I, v])((Subscript[H, 2])O) kJ","(Subscript[I, v])((Subscript[CH, 4])) kJ"}}]

Cell 8:  Calculate  the equilibrium constants

Define the stoichiometric coefficients ν(i)

compounds={"(Subscript[CH, 4])","(Subscript[H, 2])O","CO","(Subscript[H, 2])"};
nu={-1,-2,1,3};

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

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

in which the formulae appear in the list compounds

IVList=MapThread[{#1,#2,#3,#4}&,{IVCH4List, IVH2OList,IVCOList,IVH2List}] ;
toTable=MapThread[Prepend[#1,#2]&,{IVList,pListAtm}]; TableForm[toTable, TableHeadings->{None,{"p atm","(Subscript[I, v])((Subscript[CH, 4])) kJ","(Subscript[I, v])((Subscript[H, 2])O) kJ","(Subscript[I, v])(CO) kJ","(Subscript[I, v])((Subscript[H, 2])) kJ"}}]

Take the scalar product of nu with the elements of the list IVList to get  
		(∑_( i=1)^( 4))ν(i) (Subscript[I, v])(i; 1000K,p) ; 
add ln[K(T, 1bar)  to it.
To calculate ln[K(T,p] use
ln[K(T,p)]=ln[K(T, 1bar)] - (∑_( i=1)^4)(((ν(i))/(R T )))(Subscript[I, v])(i; T,p) 	
The terms (Subscript[I, v]) were calculated in kJ.  Therefore I must use 
	R = 8.314 (10^(-3))  kJ 
to make the sum dimensionless.  Therefore this equation becomes
	ln[K(T,p)]=ln[K(T, 1bar)] - 
	(∑_( i=1)^4)(((ν(i))/(8.314 (10^(-3)) T )))(Subscript[I, v])(i; T,p) 
In this equation T = 1000 K

lnKp=Log[K1000]- (1/(1000*8.314*10^-3))*(nu.#&/@IVList)  (*this is Eq. 1 with R=8.314*10^-3 kJ/Kmol and T=1000K*)

Make a list of equilibrium constants (use K = exp[ln[K]]  )

KpList=Exp[lnKp]

This list contains the values of K(1000K, p)  for p = 10, 50  and 100 bar, respectively.   The equilibrium constant decreases when pressure is increased.  This means that the equilibrium is displaced to the left (more methane is present) when the pressure increases. This is in agreement with Le Chatelier's principle.

Display the results. Add p = 1bar to the pressure list and K(1000K, 1 bar) to the list of equilibrium constants.  The reaction is 
	(Subscript[CH, 4])(g) + (Subscript[H, 2])O(g) = CO(g) + 3 (Subscript[H, 2])(g)

The results shown below illustrate one frequent dilemma: the equilibrium is favorable at low temperature while the rate at low temperature tends to be low.

TableForm[Transpose[{ Prepend[pListBar,1],Prepend[KpList,K1000]}],
TableHeadings->{None,{"p bar","K"}}]