High-Pressure Chemical Equilibrium Calculations​​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;​​​​R=8.31410^-3;​​​​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;
Mole fractions
In[]:=
y[i_]:=n[i]/Sum[n[j],{j,1,3}];
Total Gibbs free energy including GR the Residual Gibbs free energy contribution (obtained using the Peng-Robinson Equation of State)
In[]:=
k[1,1]=0;k[2,2]=0;k[3,3]=0;​​k[2,1]=k[1,2];k[3,1]=k[1,3];k[2,3]=k[3,2];​​k[1,2]=-0.036;k[1,3]=0;k[3,2]=0.222;​​​​P1=1;Rm=8.31410^-3;​​​​m[i_]:=0.37464+1.54226ω[i]-0.26992ω[i]^2;​​​​Pc[1]=13.16;Tc[1]=33.44;ω[1]=-0.1201;(*hydrogen*)​​Pc[2]=33.94;Tc[2]=126.194;ω[2]=0.04;(*nitrogen*)​​Pc[3]=112.8;Tc[3]=405.549;ω[3]=0.255;(*ammonia*)​​​​Tre[i_,T_]:=T/Tc[i];​​ap[i_,T_]:=0.45724a[i,T]Rm^2Tc[i]^2/Pc[i];​​Ap[i_,P_,T_]:=ap[i,T]P/(RmT)^2;bp[i_,T_]:=0.07780RmTc[i]/Pc[i];​​Bp[i_,P_,T_]:=Pbp[i,T]/(RmT);​​a[i_,T_]:=(1+m[i](1-Tre[i,T]^0.5))^2;​​Av[P_,T_]:=
3
∑
i=1
3
∑
j=1
y[i]y[j]Ab[i,j,P,T];​​Ab[i_,j_,P_,T_]:=(Ap[i,P,T]Ap[j,P,T])^0.5(1-k[i,j]);​​Bv[P_,T_]:=
3
∑
i=1
y[i]Bp[i,P,T];​​GR[T_,P_]:=RmT(Sum[n[i],{i,1,3}])*(Zv-1-Log[(Zv-(Bv[P,T]))]-Av[P,T]/2^1.5/Bv[P,T]*Log[((Zv+(Bv[P,T])(1+2^.5))/(Zv+(Bv[P,T])(1-2^.5)))]);​​
In[]:=
Gtot[x_,P_]:=Sum[n[j]*(G[x,P,j]+RxLog[n[j]/Sum[n[k],{k,1,3}]]),{j,1,3}]+GR[x,P];
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_,Zv_,l1_,l2_,l3_,n1_,n2_,n3_]:=Gtot[x,P]+l1(Zv^3+(-1+Bv[P,T])Zv^2+Zv(Av[P,T]-3Bv[P,T]^2-2Bv[P,T])-Av[P,T]Bv[P,T]+Bv[P,T]^2+Bv[P,T]^3)+l2(2n[1]+3n[3]-1.5)+l3(2n[2]+n[3]-0.5)/.n[1]n1/.n[2]n2/.n[3]n3
In[]:=
​​f1[Zv_,l1_,l2_,l3_,n1_,n2_,n3_,T_]:=Evaluate[Re[D[Lagragian[T,P,Zv,l1,l2,l3,n1,n2,n3],Zv]]];​​​​f2[Zv_,l1_,l2_,l3_,n1_,n2_,n3_,T_]:=Evaluate[Re[D[Lagragian[T,P,Zv,l1,l2,l3,n1,n2,n3],l1]]];​​​​f3[Zv_,l1_,l2_,l3_,n1_,n2_,n3_,T_]:=Evaluate[Re[D[Lagragian[T,P,Zv,l1,l2,l3,n1,n2,n3],l2]]];​​​​f4[Zv_,l1_,l2_,l3_,n1_,n2_,n3_,T_]:=Evaluate[Re[D[Lagragian[T,P,Zv,l1,l2,l3,n1,n2,n3],l3]]];​​​​f5[Zv_,l1_,l2_,l3_,n1_,n2_,n3_,T_]:=Evaluate[Re[D[Lagragian[T,P,Zv,l1,l2,l3,n1,n2,n3],n1]]];​​​​f6[Zv_,l1_,l2_,l3_,n1_,n2_,n3_,T_]:=Evaluate[Re[D[Lagragian[T,P,Zv,l1,l2,l3,n1,n2,n3],n2]]];​​​​f7[Zv_,l1_,l2_,l3_,n1_,n2_,n3_,T_]:=Evaluate[Re[D[Lagragian[T,P,Zv,l1,l2,l3,n1,n2,n3],n3]]];
Equilibrium composition at 600 bar and 800 K
In[]:=
P=600;
In[]:=
FindRoot[{f1[Zv,l1,l2,l3,n1,n2,n3,Tn]0,​​f2[Zv,l1,l2,l3,n1,n2,n3,Tn]0,​​f3[Zv,l1,l2,l3,n1,n2,n3,Tn]0,​​f4[Zv,l1,l2,l3,n1,n2,n3,Tn]0,​​f5[Zv,l1,l2,l3,n1,n2,n3,Tn]0,​​f6[Zv,l1,l2,l3,n1,n2,n3,Tn]0,​​f7[Zv,l1,l2,l3,n1,n2,n3,Tn]0}/.Tn800,{Zv,0.9},{l1,1},{l2,1},{l3,1},{n1,0.05},{n2,0.01},{n3,0.5}]
Out[]=
{Zv1.14907,l11.50466×
-16
10
,l236.9846,l364.8994,n10.366855,n20.122285,n30.25543}