Calculation of High-Pressure 
Chemical Equilibrium
Author's Data: Housam Binous
Department of Chemical Engineering
National Institute of Applied Sciences and Technology
Tunis, TUNISIA
Email: binoushousam@yahoo.com

Acknowledgement :
This problem was presented in the excellent book by S. I. Sandler, Chemical and Engineering Thermodynamics, 3rd Edition, Wiley 1999 (Illustration 9.1-8 page 654).

Off[General::"spell1"]

Equilibrium constant versus temperature for the reaction:

			(1/2)(Subscript[N, 2])+(3/2)(Subscript[H, 2])=(Subscript[NH, 3])

Ka298=Exp[16.5 1000/(8.314 298.15)]

ΔH[T_]:=-46100+(T-298.15) Δa+Δb/2 (T^2-298.15^2)+Δc/3 (T^3-298.15^3)+Δd/4 (T^4-298.15^4)

Δa=24.619-0.5 27.318-1.5 26.879

Δb=(3.75-0.5 (0.623)-1.5 (0.435)) 10^-2

Δc=(-0.138+0.5 (0.095)+1.5 (0.033)) 10^-5

Δd=(0.5 (2.871)+1.5 (0.87)) 10^-9

Ka[T_]:=Ka298 Exp[Δa/8.314 Log[T/298.15]+Δb/(2 8.314) (T-298.15)+Δc/(6 8.314) (T^2-298.15^2)+Δd/(12 8.314) (T^3-298.15^3)+1/8.314 (46100+(298.15) Δa+Δb/2 (298.15^2)+Δc/3 (298.15^3)+Δd/4 (298.15^4)) (1/T-1/298.15)]

Ka[500]

Ka[600]

Ka[700]

Ka[800]

Reduced Pressure and Temperature

Tre[i_]:=T/Tc[i]

Pre[i_]:=P/Pc[i]

Fugacity coefficients using the Peng Robinson Equation of State. Binary interaction parameters for the Peng Robinson EOS were taken from HYSYS 3.2.

ap[i_]:=0.45724 a[i]  R^2 Tc[i]^2/Pc[i]

Ap[i_]:=ap[i] P/(R T)^2;

bp[i_]:=0.07780 R Tc[i]/Pc[i]

Bp[i_]:=P bp[i]/(R T)

(Subscript[ϕ, v])[i_]:=Exp[((Subscript[Z, v])-1) Bp[i]/(Subscript[B, v])-Log[(Subscript[Z, v])-(Subscript[B, v])]-(Subscript[A, v])/(2 (Sqrt[2])(Subscript[B, v])) (2 (∑_(j=1)^3)y[j] Ab[i,j]/(Subscript[A, v])-Bp[i]/(Subscript[B, v])) Log[((Subscript[Z, v])+(1+(Sqrt[2]))(Subscript[B, v]))/((Subscript[Z, v])+(1-(Sqrt[2]))(Subscript[B, v]))]]

a[i_]:=(1+m[i] (1-Tre[i]^0.5))^2

m[i_]:=0.37464+1.54226 ω[i]-0.26992 ω[i]^2

(Subscript[A, v])=(∑_(i=1)^3)(∑_(j=1)^3)y[i] y[j] Ab[i,j];

Ab[i_,j_]=(Ap[i] Ap[j])^0.5(1-k[i,j]);

(Subscript[B, v])=(∑_(i=1)^3)y[i] Bp[i];

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;

Gas constant

R=8.314;

Critical Pressure and Temperature and acentric factor for Hydrogen, Nitrogen and Ammonia

Pc[1]=13.16;Tc[1]=-239.71+273.15;ω[1]=0.0;

Pc[2]=33.94;Tc[2]=-146.96+273.15;ω[2]=0.0206;

Pc[3]=112.77;Tc[3]=132.40+273.15;ω[3]=0.2582;

The correction factor (Subscript[K, v]) and the extent of reaction X versus P for T=500K or 900R

T=500;

P=0.001;For[i=1,i<17,{sol900=FindRoot[{Ka[T](0.5 (1-X))^0.5 (1.5 (1-X))^1.5 (Subscript[ϕ, v])[1]^1.5 (Subscript[ϕ, v])[2]^0.5==  X (2-X)/ P (Subscript[ϕ, v])[3],y[3](2-X)==  X,
y[2](2-X)==  0.5(1-X),
y[1](2-X)== 1.5 (1-X),
(Subscript[Z, v])^3+(-1+(Subscript[B, v]))(Subscript[Z, v])^2+(Subscript[Z, v]) ((Subscript[A, v])-3(Subscript[B, v])^2-2(Subscript[B, v]))-(Subscript[A, v]) (Subscript[B, v])+(Subscript[B, v])^2+(Subscript[B, v])^3==0},
{y[2],0.2},{y[1],0.1},{y[3],0.4},{X,0.95},{(Subscript[Z, v]),0.95},MaxIterations-> 10000],Pp[i]=P,Kv[i]=(Subscript[ϕ, v])[3]/((Subscript[ϕ, v])[1]^1.5 (Subscript[ϕ, v])[2]^0.5)/.sol900,X[i]=X/.sol900,Print[i," ",P],P=P+100,i=i+1}]

tbl900Kv=Table[{Pp[i],Chop[Kv[i]]},{i,1,16}];

pltkv900=ListPlot[tbl900Kv,Joined-> True, PlotStyle-> Hue[0.1],AspectRatio->1];

tbl900X=Table[{Pp[i],Chop[X[i]]},{i,1,16}];

pltX900=ListPlot[tbl900X,Joined-> True, PlotStyle-> Hue[0.1],AspectRatio->1];

The correction factor (Subscript[K, v]) and the extent of reaction X versus P for T=600K or 1080R

T=600;

P=0.001;For[i=1,i<17,{sol1080=FindRoot[{Ka[T](0.5 (1-X))^0.5 (1.5 (1-X))^1.5 (Subscript[ϕ, v])[1]^1.5 (Subscript[ϕ, v])[2]^0.5==  X (2-X)/ P (Subscript[ϕ, v])[3],y[3](2-X)==  X,
y[2](2-X)==  0.5(1-X),
y[1](2-X)== 1.5 (1-X),
(Subscript[Z, v])^3+(-1+(Subscript[B, v]))(Subscript[Z, v])^2+(Subscript[Z, v]) ((Subscript[A, v])-3(Subscript[B, v])^2-2(Subscript[B, v]))-(Subscript[A, v]) (Subscript[B, v])+(Subscript[B, v])^2+(Subscript[B, v])^3==0},
{y[2],0.2},{y[1],0.1},{y[3],0.4},{X,0.6},{(Subscript[Z, v]),0.95},MaxIterations-> 10000],Pp[i]=P,Kv[i]=(Subscript[ϕ, v])[3]/((Subscript[ϕ, v])[1]^1.5 (Subscript[ϕ, v])[2]^0.5)/.sol1080,X[i]=X/.sol1080,Print[i," ",P],P=P+100,i=i+1}]

tbl1080Kv=Table[{Pp[i],Chop[Kv[i]]},{i,1,16}];

pltkv1080=ListPlot[tbl1080Kv,Joined-> True, PlotStyle-> Hue[0.3],AspectRatio->1];

tbl1080X=Table[{Pp[i],Chop[X[i]]},{i,1,16}];

pltX1080=ListPlot[tbl1080X,Joined-> True, PlotStyle-> Hue[0.3],AspectRatio->1];

The correction factor (Subscript[K, v]) and the extent of reaction X versus P for T=700K or 1260R

T=700;

P=0.001;For[i=1,i<17,{sol1260=FindRoot[{Ka[T](0.5 (1-X))^0.5 (1.5 (1-X))^1.5 (Subscript[ϕ, v])[1]^1.5 (Subscript[ϕ, v])[2]^0.5==  X (2-X)/ P (Subscript[ϕ, v])[3],y[3](2-X)==  X,
y[2](2-X)==  0.5(1-X),
y[1](2-X)== 1.5 (1-X),
(Subscript[Z, v])^3+(-1+(Subscript[B, v]))(Subscript[Z, v])^2+(Subscript[Z, v]) ((Subscript[A, v])-3(Subscript[B, v])^2-2(Subscript[B, v]))-(Subscript[A, v]) (Subscript[B, v])+(Subscript[B, v])^2+(Subscript[B, v])^3==0},
{y[2],0.2},{y[1],0.1},{y[3],0.4},{X,0.6},{(Subscript[Z, v]),0.95},MaxIterations-> 10000],Pp[i]=P,Kv[i]=(Subscript[ϕ, v])[3]/((Subscript[ϕ, v])[1]^1.5 (Subscript[ϕ, v])[2]^0.5)/.sol1260,X[i]=X/.sol1260,Print[i," ",P],P=P+100,i=i+1}]

tbl1260Kv=Table[{Pp[i],Chop[Kv[i]]},{i,1,16}];

pltkv1260=ListPlot[tbl1260Kv,Joined-> True, PlotStyle-> Hue[0.5],AspectRatio->1];

tbl1260X=Table[{Pp[i],Chop[X[i]]},{i,1,16}];

pltX1260=ListPlot[tbl1260X,Joined-> True, PlotStyle-> Hue[0.5],AspectRatio->1];

The correction factor (Subscript[K, v]) and the extent of reaction X versus P for T=800K or 1440R

T=800;

P=0.001;For[i=1,i<17,{sol1440=FindRoot[{Ka[T](0.5 (1-X))^0.5 (1.5 (1-X))^1.5 (Subscript[ϕ, v])[1]^1.5 (Subscript[ϕ, v])[2]^0.5==  X (2-X)/ P (Subscript[ϕ, v])[3],y[3](2-X)==  X,
y[2](2-X)==  0.5(1-X),
y[1](2-X)== 1.5 (1-X),
(Subscript[Z, v])^3+(-1+(Subscript[B, v]))(Subscript[Z, v])^2+(Subscript[Z, v]) ((Subscript[A, v])-3(Subscript[B, v])^2-2(Subscript[B, v]))-(Subscript[A, v]) (Subscript[B, v])+(Subscript[B, v])^2+(Subscript[B, v])^3==0},
{y[2],0.2},{y[1],0.1},{y[3],0.4},{X,0.6},{(Subscript[Z, v]),0.95},MaxIterations-> 10000],Pp[i]=P,Kv[i]=(Subscript[ϕ, v])[3]/((Subscript[ϕ, v])[1]^1.5 (Subscript[ϕ, v])[2]^0.5)/.sol1440,X[i]=X/.sol1440,Print[i," ",P],P=P+100,i=i+1}]

tbl1440Kv=Table[{Pp[i],Chop[Kv[i]]},{i,1,16}];

pltkv1440=ListPlot[tbl1440Kv,Joined-> True, PlotStyle-> Hue[0.7],AspectRatio->1];

tbl1440X=Table[{Pp[i],Chop[X[i]]},{i,1,16}];

pltX1440=ListPlot[tbl1440X,Joined-> True, PlotStyle-> Hue[0.7],AspectRatio->1];

Plots of the correction factor (Subscript[K, v]) and the extent of reaction X versus P for T from 500K to 800K

Show[pltX900,pltX1080,pltX1260,pltX1440,PlotRange->All]

Show[pltkv900,pltkv1080,pltkv1260,pltkv1440,PlotRange->All]