Chemical Equilibrium Calculations Applied to the Haber Process ​
​Authors: Housam Binous and Ahmed Bellagi

Reference : Housam Binous and Ahmed Bellagi, Calculations of complex chemical reaction equilibria using stoichiometric and non-stoichiometric approaches in combination with arc-length continuation, Engineering Reports, 4(10) , 2022.
Thermodynamic data obtained from the NIST WebBook : https://webbook.nist.gov/
In[]:=
(*-----------------HYDROGENSECTION------------------------------------------*)​​WBH2={33.066178,-11.363417,11.432816,-2.772874,-0.158558,-9.980797,172.707974,0.0};(*ShomateEquationParametersforH2*)​​t=Tn/1000;​​THH2={t,t^2/2,t^3/3,t^4/4,-1/t,1,0,-1};​​TSH2={Log[t],t,t^2/2,t^3/3,-1/(2*t^2),0,1,0};​​​​HH2=Sum[(WBH2*THH2)[[i]],{i,1,8}];​​SH2=Sum[(WBH2*TSH2/1000)[[i]],{i,1,8}];​​H298H2=0.000;​​GjH2[Tn_]:=H298H2+HH2-Tn*SH2;​​​​(*-----------------NITROGENSECTION------------------------------------------*)​​WBN2={19.50583,19.88705,-8.598535,1.369784,0.527601,-4.935202,212.3900,0.0};(*ShomateEquationParametersforN2*)​​t=Tn/1000;​​THN2={t,t^2/2,t^3/3,t^4/4,-1/t,1,0,-1};​​TSN2={Log[t],t,t^2/2,t^3/3,-1/(2*t^2),0,1,0};​​​​HN2=Sum[(WBN2*THN2)[[i]],{i,1,8}];​​SN2=Sum[(WBN2*TSN2/1000)[[i]],{i,1,8}];​​H298N2=0.000;​​GjN2[Tn_]:=H298N2+HN2-Tn*SN2;​​​​(*-----------------AMMONIASECTION------------------------------------------*)​​WBNH3={19.99563,49.77119,-15.37599,1.921168,0.189174,-53.30667,203.8591,-45.89806};(*ShomateEquationParametersforNH3*)​​t=Tn/1000;​​THNH3={t,t^2/2,t^3/3,t^4/4,-1/t,1,0,-1};​​TSNH3={Log[t],t,t^2/2,t^3/3,-1/(2*t^2),0,1,0};​​HNH3=Sum[(WBNH3*THNH3)[[i]],{i,1,8}];​​SNH3=Sum[(WBNH3*TSNH3/1000)[[i]],{i,1,8}];​​H298NH3=-45.940;​​GjNH3[Tn_]:=H298NH3+HNH3-Tn*SNH3;​​​​​​G[x_,P_,1]:=GjH2[x]+RxLog[P/P1];​​G[x_,P_,2]:=GjN2[x]+RxLog[P/P1];G[x_,P_,3]:=GjNH3[x]+RxLog[P/P1];
In[]:=
P1=1;R=8.31410^-3;
Mole fractions
In[]:=
y[i_]:=n[i]/Sum[n[j],{j,1,3}];
Total Gibbs free energy assuming contribution of Residual Gibbs free energy is negligible
In[]:=
Gtot[x_,P_]:=Sum[n[j]*(G[x,P,j]+RxLog[n[j]/Sum[n[k],{k,1,3}]]),{j,1,3}];
Building the Lagrangian, for a stoichiometric mixture of hydrogen and nitrogen, in order to minimize Gtot under the equality constraints (i.e., the
material balance equations).
In[]:=
Lagragian[x_,P_,l1_,l2_,n1_,n2_,n3_]:=Gtot[x,P]+l1(2n[1]+3n[3]-2×0.75)+l2(2n[2]+n[3]-2×0.25)/.n[1]n1/.n[2]n2/.n[3]n3
In[]:=
​​f1[l1_,l2_,n1_,n2_,n3_,T_]:=Evaluate[Re[D[Lagragian[T,P,l1,l2,n1,n2,n3],l1]]];​​​​f2[l1_,l2_,n1_,n2_,n3_,T_]:=Evaluate[Re[D[Lagragian[T,P,l1,l2,n1,n2,n3],l2]]];​​​​f3[l1_,l2_,n1_,n2_,n3_,T_]:=Evaluate[Re[D[Lagragian[T,P,l1,l2,n1,n2,n3],n1]]];​​​​f4[l1_,l2_,n1_,n2_,n3_,T_]:=Evaluate[Re[D[Lagragian[T,P,l1,l2,n1,n2,n3],n2]]];​​​​f5[l1_,l2_,n1_,n2_,n3_,T_]:=Evaluate[Re[D[Lagragian[T,P,l1,l2,n1,n2,n3],n3]]];
Equilibrium composition at 600 bar and 800 K
In[]:=
P=600;
In[]:=
FindRoot[{f1[l1,l2,n1,n2,n3,Tn]0,​​f2[l1,l2,n1,n2,n3,Tn]0,​​f3[l1,l2,n1,n2,n3,Tn]0,​​f4[l1,l2,n1,n2,n3,Tn]0,​​f5[l1,l2,n1,n2,n3,Tn]0}/.Tn800,{l1,1},{l2,1},{n1,0.05},{n2,0.01},{n3,0.5}]
Out[]=
{l137.3045,l265.3731,n10.410118,n20.136706,n30.226588}