Calculation of High-Pressure Chemical Equilibrium​
​ Authors: Housam Binous and Ahmed Bellagi

H . Binous, Applications of the Peng - Robinson Equation of State using Mathematica®, Chemical Engineering Education, 42 (1), pp . 47 - 51, 2008.
​
Acknowledgement :
​​
This problem was presented by S. I. Sandler, Chemical and Engineering Thermodynamics, 3rd Edition, Wiley 1999
​
(Illustration 9.1-8 page 654).
​
In[]:=
Off[General::"spell1"]
Equilibrium constant versus temperature for the reaction:​
1
2
N
2
+
3
2
H
2
=
NH
3
​
In[]:=
Ka298=Exp[16.51000/(8.314×298.15)]
Out[]=
777.742
In[]:=
Δ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)
In[]:=
Δa=24.619-0.5×27.318-1.5×26.879
Out[]=
-29.3585
In[]:=
Δb=(3.75-0.5(0.623)-1.5(0.435))10^-2
Out[]=
0.02786
In[]:=
Δc=(-0.138+0.5(0.095)+1.5(0.033))10^-5
Out[]=
-4.1×
-7
10
In[]:=
Δd=(0.5(2.871)+1.5(0.87))10^-9
Out[]=
2.7405×
-9
10
In[]:=
Ka[T_]:=Ka298Exp[Δa/8.314Log[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)]
In[]:=
Ka[500]
Out[]=
0.32826
In[]:=
Ka[600]
Out[]=
0.04347
In[]:=
Ka[700]
Out[]=
0.00990121
In[]:=
Ka[800]
Out[]=
0.00320049
Reduced Pressure and Temperature
In[]:=
Tre[i_]:=T/Tc[i]
In[]:=
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.
In[]:=
ap[i_]:=0.45724a[i]R^2Tc[i]^2/Pc[i]
In[]:=
Ap[i_]:=ap[i]P/(RT)^2;
In[]:=
bp[i_]:=0.07780RTc[i]/Pc[i]
In[]:=
Bp[i_]:=Pbp[i]/(RT)
In[]:=
ϕ
v
[i_]:=Exp(
Z
v
-1)Bp[i]
B
v
-Log[
Z
v
-
B
v
]-
A
v
2
2
B
v
2
3
∑
j=1
y[j]Ab[i,j]
A
v
-Bp[i]
B
v
Log
Z
v
+1+
2

B
v

Z
v
+1-
2

B
v

In[]:=
a[i_]:=(1+m[i](1-Tre[i]^0.5))^2
In[]:=
m[i_]:=0.37464+1.54226ω[i]-0.26992ω[i]^2
In[]:=
A
v
=
3
∑
i=1
3
∑
j=1
y[i]y[j]Ab[i,j];
In[]:=
Ab[i_,j_]=(Ap[i]Ap[j])^0.5(1-k[i,j]);
In[]:=
B
v
=
3
∑
i=1
y[i]Bp[i];
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;
Gas constant
In[]:=
R=8.314;
Critical Pressure and Temperature and acentric factor for Hydrogen, Nitrogen and Ammonia
In[]:=
Pc[1]=13.16;Tc[1]=-239.71+273.15;ω[1]=0.0;
In[]:=
Pc[2]=33.94;Tc[2]=-146.96+273.15;ω[2]=0.0206;
In[]:=
Pc[3]=112.77;Tc[3]=132.40+273.15;ω[3]=0.2582;
The correction factor
K
v
and the extent of reaction X versus P for T=500K or 900R
In[]:=
T=500;
In[]:=
P=0.001;For[i=1,i<17,{sol900=FindRoot[{Ka[T](0.5(1-X))^0.5(1.5(1-X))^1.5
ϕ
v
[1]^1.5
ϕ
v
[2]^0.5X(2-X)/P
ϕ
v
[3],y[3](2-X)X,​​y[2](2-X)0.5(1-X),​​y[1](2-X)1.5(1-X),​​
Z
v
^3+(-1+
B
v
)
Z
v
^2+
Z
v
(
A
v
-3
B
v
^2-2
B
v
)-
A
v
B
v
+
B
v
^2+
B
v
^30},​​{y[2],0.2},{y[1],0.1},{y[3],0.4},{X,0.95},{
Z
v
,0.95},MaxIterations10000],Pp[i]=P,Kv[i]=
ϕ
v
[3]/(
ϕ
v
[1]^1.5
ϕ
v
[2]^0.5)/.sol900,X[i]=X/.sol900,Print[i," ",P],P=P+100,i=i+1}]
1 0.001
2 100.001
3 200.001
4 300.001
5 400.001
6 500.001
7 600.001
8 700.001
9 800.001
10 900.001
11 1000.
12 1100.
13 1200.
14 1300.
15 1400.
16 1500.
In[]:=
tbl900Kv=Table[{Pp[i],Chop[Kv[i]]},{i,1,16}];
In[]:=
pltkv900=ListPlot[tbl900Kv,JoinedTrue,PlotStyleHue[0.1],AspectRatio1];
In[]:=
tbl900X=Table[{Pp[i],Chop[X[i]]},{i,1,16}];
In[]:=
pltX900=ListPlot[tbl900X,JoinedTrue,PlotStyleHue[0.1],AspectRatio1];
The correction factor
K
v
and the extent of reaction X versus P for T=600K or 1080R
In[]:=
T=600;
In[]:=
P=0.001;For[i=1,i<17,{sol1080=FindRoot[{Ka[T](0.5(1-X))^0.5(1.5(1-X))^1.5
ϕ
v
[1]^1.5
ϕ
v
[2]^0.5X(2-X)/P
ϕ
v
[3],y[3](2-X)X,​​y[2](2-X)0.5(1-X),​​y[1](2-X)1.5(1-X),​​
Z
v
^3+(-1+
B
v
)
Z
v
^2+
Z
v
(
A
v
-3
B
v
^2-2
B
v
)-
A
v
B
v
+
B
v
^2+
B
v
^30},​​{y[2],0.2},{y[1],0.1},{y[3],0.4},{X,0.6},{
Z
v
,0.95},MaxIterations10000],Pp[i]=P,Kv[i]=
ϕ
v
[3]/(
ϕ
v
[1]^1.5
ϕ
v
[2]^0.5)/.sol1080,X[i]=X/.sol1080,Print[i," ",P],P=P+100,i=i+1}]
1 0.001
2 100.001
3 200.001
4 300.001
5 400.001
6 500.001
7 600.001
8 700.001
9 800.001
10 900.001
11 1000.
12 1100.
13 1200.
14 1300.
15 1400.
16 1500.