© 2022 Horia Metiu

WorkBookT.13:  Thermochemistry

WorkBookT.13.1  Calculate the heat of the reaction (Subscript[N, 2]) +3(Subscript[H, 2]) ->2 (Subscript[NH, 3])

Calculate the heat of the reaction
	(Subscript[N, 2]) +3(Subscript[H, 2]) ->2 (Subscript[NH, 3])
at
	T = 623.15 K and p = 394 atm
when you know that the heat of reaction at
	(Subscript[T, i]) = 298.15 K and  (Subscript[p, i]) = 1 atm 
is equal to 
	ΔH(298 K, 1 atm) = - 92.22 kJ

Cell 1:  Explain the procedure

The reaction takes place at constant pressure therefore the heat of the reaction is equal to the enthalpy of reaction.

I want to calculate the change of the enthalpy of reaction in the transformation
	 { T = 298.1, p = 1,(Subscript[v, i])} -> {T = 623.15 K, p = 394 atm, (Subscript[v, f]) }
Enthalpy is a function of state and its change in a transformation is independent of path.  Therefore I can chose a convenient path consisting of two segments:

Segment 1:  { T = 298.15, p = 1 atm,(Subscript[v, i])} -> {T = 623.15 K, p = 1 atm, (Subscript[v, 1]) }
Segment 2: { T = 623.15 K , p = 1,(Subscript[v, 1])} -> {T = 623.15 K, p = 394 atm, (Subscript[v, f]) }
 • I denote  (Subscript[ΔH, 1]) the enthalpy change along segment 1
 • I denote  (Subscript[ΔH, 2]) the enthalpy change along segment 2
 • The enthalpy change in the transformation
 	{ T=298.15 K , p=1} -> {T=623.15 K, p=394 atm } 
 is
 	ΔH = (Subscript[ΔH, 1]) + (Subscript[ΔH, 2])

The change of reaction enthalpy on the first segment
	{ T = 298.1, p = 1,(Subscript[v, i])} -> {T = 623.15 K, p = 1 atm, (Subscript[v, 1]) }
 is 
	(Subscript[ΔH, 1]) ≡ ΔH(298.15 ,p=1 ) +(∫_( 298.15)^( 623.15)) (Subscript[Δc, p])(T p=1  )dT 		(1)
with
	(Subscript[Δc, p])(T,p) = (∑_( α=1)^( c))(Subscript[ν, α])(Subscript[c, p])((Subscript[A, α]); T, p)						(2)
• (Subscript[ν, α])  is the stoichiometric coefficient of the compound (Subscript[A, α]).  
• The stoichiometric coefficients of the products are positive, those of the reactants are negative.
• I want ΔH in kilocalorie and therefore I must use (Subscript[c, p]) in kilocalories/ mol K

The change of reaction enthalpy along the second segment
	{ T = 623.15 K , p = 1,(Subscript[v, 1])} -> {T = 623.15 K, p = 394 atm, (Subscript[v, f]) }
is
	(Subscript[ΔH, 2])≡(∑_(α=1)^c)(Subscript[ν, α])I((Subscript[A, α]); (Subscript[T, f]))							(3)
with
 	I((Subscript[A, α]); (Subscript[T, f])) = (24.2173/1000)*(∫_(Subscript[v, 1])^( (Subscript[v, f])))(Subscript[{v (Subscript[( ((∂p)/(∂v))), T])  + (Subscript[T(((∂p)/(∂T))), v])}, (Subscript[A, α]), (Subscript[T, f])])dv	(4)
• To evaluate the integral one uses the equation of state of compound (Subscript[A, α])
•  (Subscript[v, f]) is the molar volume in the state  {T = 623.15 K, p = 394 atm} and  
• (Subscript[v, 1]) is the molar volume of the state {T = 623.15 K, p = 1 atm}
• These equations were derived in Chapter 13, §12. 
• The factor multiplying the integral converts liter atm to kilocalories
• I want all thermal units in kcal

Cell 2: Evaluate the change of the reaction enthalpy with temperature

The change of reaction enthalpy along segment 1
	{ T = 298.1, p = 1,(Subscript[v, i])} -> {T = 623.15 K, p = 1 atm, (Subscript[v, 1]) }
is 
	ΔH1(623.15, p= 1)= ΔH(298.15,p=1) + (∫_298.15^623.15)(Subscript[Δc, p])(T, p = 1)dT

Cell 2.1:  Input the data

I have the following data
	ΔH(T= 298.15 K,p=1 atm) = - 92.22 kJ   
	 
	(Subscript[C, p])((Subscript[N, 2])(g))  =  27.565 + 5.230 (10^(-3))T - 0.04 (10^(-7))(T^2)         J/K mol
	
	(Subscript[C, p])((Subscript[H, 2])(g))  =  +29.894 - 0.836 (10^(-3)) T + 20.17 (10^(-7 ))(T^2)   J/K mol
	
	(Subscript[C, p])((Subscript[NH, 3])(g))= 25.895 + 32.999 (10^(-3))T - 30.46 (10^(-7))(T^2)     J/K mol
	
I will use Benedict-Webb-Rubin equation of state with the parameters given in Appendix 6.

Input the data.  ΔH(298.15 K, p = 1atm) is given in kJ and the (Subscript[c, p])-s in J/mol K.  Convert units of the data to kcal and kcal/mol K.  The factor in purple color performs the conversion.

ΔH298 = - 92.22*0.239 (*kcal*);

cpN2[T_] := (0.239/1000)*(27.565 + 5.23*T*10^-3 - 0.04*T^2*10^-7) (*kcal/mol K*)

cpH2[T_] :=  (0.239/1000)*(29.894 - 0.836*T*10^-3 + 20.17*T^2*10^-7)  (*kcal/mol K*)

cpNH3[T_] :=  (0.239/1000)*(25.895 + 32.999*T*10^-3 - 30.46*T^2*10^-7 )(*kcal/mol K*)

Cell 2.2 : Calculate (Subscript[Δc, p])

Define some lists needed later. 
Make a list with the names of the compounds in the order they appear in the reaction
	(Subscript[N, 2]) + 3(Subscript[H, 2]) = 2 (Subscript[NH, 3])

compounds={N2,H2,NH3};

Make a list with the stoichiometric coefficients with the appropriate signs (minus for reactants plus for products) in the order they appear in the reaction (Subscript[N, 2]) + 3(Subscript[H, 2]) = 2(Subscript[NH, 3])

stoList={-1.,-3.,2.};

Make a list with the heat capacities at constant pressure for each compound,  in the order in which the compounds appear in the reaction (and in the list compounds)

cpList={cpN2[T],cpH2[T],cpNH3[T]}  (*kcal/mol K*)

Calculate (Subscript[Δc, p])(T, p=1atm) using Eq. (2). 
For the specific reaction discussed here (Subscript[Δc, p]) is:
	(Subscript[Δc, p])(T)  =  2 (Subscript[c, p])((Subscript[NH, 3]);T) - (Subscript[c, p])((Subscript[N, 2]);T)  - 3 (Subscript[c, p])((Subscript[H, 2]);T)    at p = 1 atm

Clear[Δcp]
	Δcp[T_]= Plus@@(stoList*cpList) //Expand  (*in kcal/mol K*)

Mathematica explanation:

stoList*cpList 
multiplies each (Subscript[c, p ])in cpList with the corresponding stoichiometric coefficient in stoList to give

intermediateList=stoList*cpList

Plus@@AnyList 
adds all the terms in the list AnyList.

Plus@@intermediateList

Calculate (Subscript[Δc, p])(400) to check that the order of magnitude is right.

Δcp[400] (*kcal/mol K*)

Cell 2.3:  Define the function ΔH1[T] that gives  ΔH(T, p = 1 atm).

Clear[ΔH1]
ΔH1[T_]:= (ΔH298 + Integrate[Δcp[x],{x,298.15,T}])//Expand  (*kcal*)
ΔH1[T]

Calculate ΔH(T = 623.15 K , p = 1 atm):

ΔH1[623.15] (*kcal*)

Recall the heat of reaction at T = 298.15 K

ΔH298

• At 298 K and p = 1 atm the the heat of reaction is -22.04 kcal.  
• At T = 623.15 K and p = 1 atm the heat of reaction is -25.07  kcal. 
• The change in the heat of reaction with temperature is fairly substantial

Cell 2.4:  Look how ΔH(T, p=1atm) changes with temperature.

TableForm[Table[{i,Evaluate[ΔH1[i]]},{i,300,1000,100}], 
TableHeadings->{None,{"T K", "ΔH1(T,1atm) kJ"}}]

Plot ΔH1( T, p=1atm).

Plot[ΔH1[T],{T,298,1000},
	AxesLabel->{"T  K","ΔH1/kJ"},
	PlotLabel->"(Subscript[N, 2]) + 3 (Subscript[H, 2]) ⟶ (Subscript[NH, 3]) "]

If I assumed that (Subscript[c, p]) of the compounds participating in the reaction are independent of temperature, I would 
have obtained that 
	ΔH(T,1atm)= ΔH(298.15 K, 1 atm) + (Subscript[Δc, p])* (T-298.15). 
This equation predicts a linear dependence of the heat for reaction on temperature.  If this were true, the plot would have been a straight line.  

Assuming that (Subscript[c, p]) is temperature independent is not a good approximation. It is reasonable though to make this assumption if the final  temperature is not very different from the initial temperature.

Cell 3: The change of enthalpy with pressure

Cell 3.1 The procedure

I want to calculate the change of reaction enthalpy in the transformation
Segment 2:
	{ T = 623.15 K , p = 1,(Subscript[v, 1])} -> {T = 623.15 K, p = 394 atm, (Subscript[v, f]) }
This is given by
	ΔH2≡(∑_(α=1)^c)(Subscript[ν, α])I((Subscript[A, α]); (Subscript[T, f]))
where
	I((Subscript[A, α])) = (24.2173/1000) (∫_( (Subscript[v, 1]))^(Subscript[v, f]))(Subscript[{v(Subscript[(((∂p)/(∂v))), T])+((Subscript[T(((∂p)/(∂T))), v])}, (Subscript[A, α])])dv   kcal
• The subscript (Subscript[A, α]) reminds me that I must use the equation of state for the compound (Subscript[A, α]).
• The factor 24.2173 converts liter atm to calories
• Division by 1000 gives I((Subscript[A, α])) in kilocalories
• ΔH2 is the change of the heat of reaction when the pressure changes from 1atm to 394 atm, and the temperature is held constant T = 623.15 K

• To evaluate I((Subscript[A, α]))  I calculate the derivative of pressure with v and T by using the equation of state  for the compound (Subscript[A, α]). 
• The limits of integration are 
	• the molar volume (Subscript[v, 1]) of the state { T = 623.15 K , p = 1,(Subscript[v, 1])} 
	• the molar volume (Subscript[v, f]) in the final state { p = 394 atm, T = 298.15 K,(Subscript[v, f])). 
	• These two volumes are obtained from the equation of state.

Cell 3.2  Define the Benedict-Webb-Rubin equation of state

The Benedict-Webb-Rubin equation is
p = R T ρ + (B R T-A-(Subscript[C, 0])T-2)(ρ^2) + (ρ^3)(b R T-a)+ a α (ρ^6)  + 
	(ρ^3)c T-2 (1 + γ (ρ^2)) exp[-γ (ρ^2)]	
with
	ρ = 1/v

Define the function pBWR which gives the pressure according to the equation given above.

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

Input the coefficients for the BWR equation (from Appendix 6)

H2Coefs={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} ;
N2Coefs={A->1.1925,  B->0.0458,
			                                            C0->5.8891*10^3, a->0.0149,b->1.98154*10^-3,  c->5.48064*10^2,α->2.91545*10^-4,γ->7.5*10^-3 } ;
NH3Coefs={A->3.78928,  B->5.16461*10^-2,
			                                            C0->1.78567*10^5, a->0.10354,b->7.19561*10^-4,  c->1.57536*10^2,α->4.651890*10^-6,γ->1.980*10^-2} ;

Cell 3.3 Calculate the needed molar  volumes for each compound

Use the function volumeG from WB3:

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

The transformation is
	{ T = 623.15 K , p = 1,(Subscript[v, 1])} -> {T = 623.15 K, p = 394 atm, (Subscript[v, f]) }

The volumes (Subscript[v, 1]) are

v1N2=volumeG[623.15,1,pBWR,N2Coefs]

v1H2=volumeG[623.15,1,pBWR,H2Coefs]

v1NH3=volumeG[623.15,1,pBWR,NH3Coefs]

At this high temperature the gases are nearly ideal and their molar volumes are nearly equal.

The final molar volumes are:

vfN2=volumeG[623.15,394,pBWR,N2Coefs]

vfH2=volumeG[623.15,394,pBWR,H2Coefs]

vfNH3=volumeG[623.15,394,pBWR,NH3Coefs]

Cell 3.4  Calculate the integrals  I((Subscript[A, α]))

The path is 
	{ T = 623.15 K , p = 1,(Subscript[v, 1])} -> {T = 623.15 K, p = 394 atm, (Subscript[v, f]) }
I need to evaluate 
	I((Subscript[A, α]))  = (∫_(Subscript[v, 1])^(Subscript[v, f])) {v (Subscript[(((∂p)/(∂v))), T])+T (Subscript[(((∂p)/(∂T))), v])}dv    for T = 623.15 K .  
There is one integral for each compound.

The integral for (Subscript[N, 2]) . The factor 24.2173/1000 is needed to  convert liter atm to kcal (see Appendix  7).

IN2=(24.2173/1000) * Integrate[v*D[pBWR[623.15,v,N2Coefs],v] + T D[pBWR[T, v, N2Coefs],T]/.T->623.15,{v,v1N2,vfN2},GenerateConditions->False] (*kcal*)

The integral for (Subscript[H, 2])

IH2 =(24.2173/1000) * Integrate[v*D[pBWR[623.15,v,H2Coefs],v] + T D[pBWR[T, v, H2Coefs],T]/.T->623.15,{v,v1H2,vfH2}, GenerateConditions->False]  (*kJ*)

The integral for (Subscript[NH, 3])

INH3 =(24.2173/1000)* Integrate[v*D[pBWR[623.15,v,NH3Coefs],v] + 
T D[pBWR[T, v, NH3Coefs],T]/.T->623.15,
{v,v1NH3,vfNH3}, GenerateConditions->False]  (*kJ*)

Cell 3.5  Calculate (Subscript[ΔH, 2]) (use Eq. 4)

ΔH2= 2 INH3 -IN2- 3 IH2  (*kJ*)

The effect of pressure change on the heat of reaction is fairly small but it cannot be neglected. Keeping the reactor at constant temperature is essential and we need to know the heat of reaction accurately.  Moreover, the equilibrium constant (which is the quantity that ultimately a chemist wants) depends exponentially on ΔH.

Cell 3.5  The heat of reaction at T = 623.15, p=394.8 atm

The heat of reaction is 
	ΔH = ΔH1[623.15] +ΔH2

{ΔH298,ΔH1[623.15],ΔH2,ΔH=ΔH1[623.15]+ΔH2}

WorkBookT.13.2  Optional Get the heat of one reaction from the heats of other reactions.

I know the standard heat of reaction for the following reactions 
	(Subscript[C, 3])(Subscript[H, 6]) + (Subscript[H, 2 ]) =    (Subscript[C, 3])(Subscript[H, 8])                       	(1)
	(Subscript[C, 3])(Subscript[H, 8]) + 5 (Subscript[O, 2]) = 3 C(Subscript[O, 2]) + 4 (Subscript[H, 2])O		(2)
	(Subscript[H, 2])O = (Subscript[H, 2])+(1/2)(Subscript[O, 2])					 (3)
 Use this information to calculate the heat of the reaction
 
	 (Subscript[C, 3])(Subscript[H, 6])+ (9/2)(Subscript[O, 2])=3 (Subscript[CO, 2])+ 3 (Subscript[H, 2])O		          (4)

Remove["Global`*"];

Procedure: I follow the instructions given in  §35 in Supplement 13.1
1) Make several lists:
a) The list comp contains the names of all the compounds that enter in the reactions.
b) Make the lists nu[[r]] which contain the stoichiometric coefficients of the  compounds in the reaction r. Include only the reactions whose heat is known (Eqs. 1-3). 
c) Make the list nur={nu[1], nu[2], ... nu[(Subscript[n, c])]} where (Subscript[n, c])is the number of compounds in the list comp.  nur[[r, c]] is the stoichiometric coefficient of the compound c in the reaction r.
d) Make a list nuR  with the stoichiometric coefficients of the compounds in the list comp for the reaction whose heat you want to calculate (Eq. 4)
2)  Set up the equations for the coefficients (Subscript[α, c]),  c=1, 2, ... (Subscript[n, c]).  (Subscript[n, c]) is the number of compounds in the list comp.  The equations are (see Eq. 13.94, in Complement 13. 1)
	       		(∑_(r=1)^(Subscript[n, r]))(Subscript[ν, c r]) (Subscript[α, r])=(Subscript[ν, c])             			(5)
 or in the notation introduced at 1) 
	       		       (∑_(r=1)^(Subscript[n, r])) nur[[c,r]]*α[r] 			(6)

Cell 1: Input the data you need

Make a list of all compounds involved in the reaction

comp={"C3H6", "H2", "C3H8", "O2", "CO2", "H2O"}

Make a list of  the stoichiometric coefficients of each reaction.

The stoichiometric coefficients for reaction 1
     C3H6+H2=C3H8

nu[1]={-1,-1,1,0,0,0}

nu[1].comp

This is saying that the reaction is  C3H6 + H2 = C3H8

The stoichiometric coefficients for reaction 2
C3H8 + 5O2 = 3 CO2 +4H2O 
are

nu[2]={0,0,-1,-5,3,4}

nu[2].comp

The stoichiometric coefficients for reaction 3
(Subscript[H, 2])O=(Subscript[H, 2]) +(1/2) (Subscript[O, 2])

nu[3]={0,1,0,1/2,0,-1}

nu[3].comp

The stoichiometric coefficients for the reaction of interest (whose heat we don't know)
(Subscript[C, 3])(Subscript[H, 6]) +9/2  (Subscript[O, 2]) =3 (Subscript[CO, 2]) + 3 (Subscript[H, 2])O

nuR={-1,0,0,-9/2,3,3}

Make the list nur containing the stoichiometric coefficients of all the reactions whose heat of reaction we know

nur=Table[nu[i],{i,1,3}]  (*we have three reactions whose heat we know*)

Check whether this works: the stoichiometric coefficient of the third compound (in the list compounds) for the reaction 1 is equal to 1. That of the fourth compound in reaction two is -5.

{nur[[1,3]],nur[[2,4]]}

Cell 2:   Write a function that given the stoichiometric coefficients of a reaction and the list of compounds returns the reaction. This is used to check that the stoichiometric coefficients are typed correctly.

You don't need to understand how this function works. The arguments are the list comp containing the names of the compounds involved in all the reactions and the list nu that contains the stoichiometric coefficients of the reaction of interest.

Clear[reaction];
reaction[comp_List,nu_List]:=Module[{sright,sleft,i},
 sleft=0;sright=0;
 For[i=1,i<=Length[nu],i++,   
     If[nu[[i]]>0,sright=sright+nu[[i]]*comp[[i]],
     sleft=sleft-nu[[i]]*comp[[i]]
         ]  ];
        Print[sleft," = ",sright]]

Test the input.

Table[reaction[comp,nu[i]],{i,1,4}];

Cell 3:  Write a function equations[nu, nur,nbReactions,nbCompounds] that will produce a list of the equations for the coefficients a[i], when given the stoichiometric coefficients nu of the reactions for which I know the heat of reaction, the stoichiometric coefficients nur for the reaction whose heat I want to calculate, he number of reactions nbReactions and the number of compounds nbCompounds involved in the reactions..

I'll write down the steps that go into the calculation of Eq. 4, 
	  	(∑_(r=1)^(Subscript[n, r]))(Subscript[ν, c r]) (Subscript[α, r])=(Subscript[ν, c]) ,     
 so that you can understand how the function works.

Set nbReactions equal to 3 and create the list unKnowns ={(Subscript[α, 1]), (Subscript[α, 2]), (Subscript[α, 3])}

nbReactions=3;
unKnowns=Table[α[i],{i,1,nbReactions}]

The list containing the stoichiometric coefficients (Subscript[⋁, c r]) is nur

nur

The list with the stoichiometric coefficients of the reaction whose heat I want to get is nuR

nuR

Create a list containing the sums

	(∑_(r=1)^(Subscript[n, r]))(Subscript[ν, c r]) (Subscript[α, r])
	
The sum is over the reactions and produces one term for each component

unKnowns.nur

Mathematica explanation: here is how the dot operation works:

vec=Table[a[i],{i,1,3}]

matr=Table[M[i][j],{i,1,3},{j,1,6}]

vec.matr

Create a list with the equations 
		(∑_(r=1)^(Subscript[n, r]))(Subscript[ν, c r]) (Subscript[α, r])=(Subscript[ν, c])

MapThread[Equal[#1,#2]&,{unKnowns.nur,nuR}]

Mathematica explanation:   you could have achieved the same result by threading the Equal function over the List head (see below the syntax)

Thread[Equal[unKnowns.nur,nuR]]

Remember that the FullForm of == is

a==b//FullForm

Here is how thread works:  it switches the head f with the head List

Thread[f[{a,b}]]

Now put all these steps into one function:

Clear[equations];
equations[nuR_List,nur_List,nbReactions_Integer,nbCompounds_Integer]:=Module[{unKnowns },
Clear[α];
unKnowns=Table[α[i],{i,1,nbReactions}];
MapThread[Equal[#1,#2]&,{unKnowns.nur,nuR}] ]

Check that it works:

eq=equations[nuR,nur,3,6]

Display the equations in a form that is easy to read:

eq//TableForm

These equations are the same as the ones derived in §38  of Complement 13.1

Cell 4:  Pick up the equations you want to use and put them in eqlist

I have three reactions with known heat of reaction. Therefore I need three equations (there are three unknown coefficients: a[1], a[2] and a[3]). I pick up the first three equations

eqList={eq[[1]],eq[[2]],eq[[3]]}

Write a function slv[eqlist, heatList, nr] that will solve the equations and give the heat of the reaction of interest. Here eqlist is the list of the equations, heatList  is the list of the known heats of reaction and nr is the number of equations picked for solving.  The number of equations is equal to the length of the heatList.

Clear[slv];
slv[eqlist_List,heatList_List, nr_Integer]:=
Module[{tb,cond},
Clear[α];
tb=Table[α[j],{j,1,nr}];
(*solve for α[1],α[2],α[3]*)
cond=Solve[eqlist,tb]//Flatten; 
(*calculate the heat of reaction*)
 tb.heatList/.cond ]

The result: the heat of the reaction is the sum of the heats of the three reactions whose heat of reaction I know.

slv[eqList,{ΔH1,ΔH2,ΔH3},3]

To get  a number you need to replace ΔH1, ΔH2 and ΔH3 with their numerical values