Chemical Equilibrium Computation
(Last Revised 11.2006)
? 2006 A. N. Nikitin, A. Yu. Zuev.
Ural State University, Chemical Department
and.nikitin@mail.ru, andrew.zuev@list.ru
Ekaterinburg, Russia.
This worksheet demonstrates calculation of equilibrium for the systems in which independent chemical reactions proceed using total Gibbs Energy minimization method.
1. Introduction
Methane is the main component of natural gas and it is a very important energy source and chemical which can be converted into a variety of other useful chemicals. The efficiency of the energy conversation methods based on combustion of methane and heat released converting is relatively low. Some new methods of an energy conversation are, therefore, still needed especially such of them which allow a direct conversation of energy, e.g. fuel cells. However, the stateoftheart systems of such type do not allow direct use of methane and the preliminary reforming of that to socalled syngase  mixture of and is, therefore, required. Processes of the methane converting to some useful chemicals (e.g. methanol, formaldehyde) include the stage of the methane reforming to syngase as well. The classic technological process of the reforming is based on homogenous interaction of methane with water:
CH4 + H2O = CO + 3H2 (1.1)
CO + H2O = CO2 + H2 (1.2)
But this way of conversation has a relatively low efficiency mostly due to two reasons: (i) severe conditions of reforming (high temperature) and (ii) absorbtion of heat during reforming caused by an endothermical character of these reactions. Alternative way of reforming involving heterogenous partial oxidation of methane can be realized under less severe conditions as compared with those for the former, thus excessive energy required for maintaining of the reforming process is not needed anymore.
Actually the main aim of a reforming process consists in obtaining either hydrogen or syngas, i.e. the mixture of H2 and CO with some ratio between them. The priority goal of this worksheet was to build the general procedure for computation of an equilibrium composition of the gas in the reactor for both homogenous and heterogenous reforming under given conditions. Since the procedure constructed allows to compute equilibrium over whole range of the considered parameters the relationship between the equilibrium composition (in any required units) and the parameters of reforming, e.g. initial composition, temperature or total pressure, can be plotted as well. The method of total Gibbs energy minimization is used in this worksheet since that is the most fundamental thermodynamic method to solve the equilibrium problem.
It is well known that the state of thermodynamic equilibrium is reached under given temperature and pressure when total Gibbs energy of the system amounts to its minimum. The main assumption accepted for the calculations is that the system is considered at sufficient high temperatures in order intermolecular interactions may be ignored, i.e. all the gaseous species in the system are treated as ideal gases. Total Gibbs energy of the system is thus calculated as a sum of Gibbs energies of all compounds involved in equilibrium:
The solving of the equilibrium problem can be, therefore, splitted into two parts:
1. Building the Gibbs energy as a function of required thermodynamic parameters.
2. Minimum of obtained function search under some constrains which makes the system physically valid. The main of them is that a partial pressure of each compound is to be positive.
Note:
This worksheet may take some time to be completely executed.
2. Initialization
> 
with(ScientificConstants):
with(linalg):
with(plots):
with(codegen,optimize,makeproc): 
> 
Constant(R):
R:=GetValue(%):
GetConstant(molar_gas_constant); 

(2.1) 
3. Compounds data
For the calculation three individual properties of the compounds are needed  heat capacity at fixed pressure as a function of temperature (), standard molar entalphy of formation (denoted as in this worksheet) and standard molar entropy ( in this worksheet). Standard molar enthalpy and entropy are defined as lists, but heat capacity is defined as a matrix () in which each string corresponds to given compound and each element of string is coefficient to temperature in some power from the list. Multiplying the matrix by the vector we get the list of heat capacities for compounds considered as a function of temperature (denoted as in this worksheet).
T0 is the standard temperature; all following data retreived from [1] and correspond to standard conditions.
> 
T0:=298.15:
Cp_data:=matrix([
[61.179098,8.01769073E03,0,825569.97,0,1024.0140,0,6604775.7], # T=298..1000K C(s)
[19.825631,3.07772746E03,0,295179.90,1.42984494E06,194.86124,0,0], # T=298..1200K H2(g)
[26.924057,1.69786815E02,0,229329.26,6.76616525E06,79.161659,0,0], # T=298..1000K O2(g)
[304.93108,1.87781656E02,0,4193372.1,0,10143.978,110598.40,0], # T=298..1300K CH4(g)
[90.753584,6.23514260E03,0,1541255.3,0,2589.2549,32063.316,0], # T=298..1700K CO(g)
[103.34460,4.74260628E03,0,40249.044,0,1748.2872,11004.741,0], # T=298..1900K CO2(g)
[25.781640,1.4949716E02,0,27999.319,5.52355894E07,0,1107.2718,0] # T=298..1100K H2O(g)
]):
hf:=[0,0,0,74873.00,110527.00,393522.00,241834.00]: # enthalpy
s:= [5.740000,130.570556,205.037556,186.141556,197.543556,213.795000,188.724556]: # entropy 
> 
r:= [0,1,1/2,2,2,1/2,1,3]: # coefficients in Cp polynome
Cp_sample:= [seq(T^r[i],i=1..nops(r))];
Cp:=multiply(Cp_data,Cp_sample): # heat capacity as function of temperature 

(3.1) 
In order to compute total Gibbs energy we need to introduce a chemical potential of each compound as
where the standard part is given by:
It is necessary to note that the chemical potentials of carbon and oxygen are fixed at given temperature since carbon is solid and oxygen partial pressure is regarded by us as external controlled parameter.
> 
var_p:= [0,1,0,1,1,1,1];
names:= ["C(s)","H2","O2","CH4", "CO", "CO2","H2O"]; # names of compounds for pretty plots 

(3.2) 
> 
for i from 1 to nops(var_p) do
mu[i]:=hf[i]+int(Cp[i],T=T0..T,'continuous')T*(s[i]+int(Cp[i]/T,T=T0..T,'continuous'))
end do:
mu:=convert(mu,list): 
4. General procedures used to solve the problem
This section contains all general procedures which will be used further in the present worksheet. They are the "makers" or the "constructors"of the structures needed for calculations and Gibbs energy building using given matrix of the reactions in the system, initial quantities of species, temperature, and pressure as well. Minimization tool and a simple procedure of accuracy estimation are also given here.
> 
SetReaction:= proc(U,IC)
global n,Ngas,Nxi,NullNxi,x,usedp,usedpn;
local i;
Nxi:= rowdim(U):
NullNxi:= [0 $ Nxi]:
n:= convert(multiply(['xi[i]' $ 'i'=1..Nxi],U),list)+IC; # quantities
Ngas:= sum('var_p[i]*n[i]','i'=1..coldim(U)): # total quantity of gaseous species
x:= zip ((x,y) > x*y/Ngas, n, var_p): # molar ratios
print(x);
usedp:=[]:
for i from 1 to nops(var_p) do
if n[i]<>0 then
usedp:=[op(usedp),i]
end if
end do:
usedpn:= nops(usedp)
end proc: 
Total Gibbs energy is given as
Since we are interested to compute equilibrium properties under different conditions we will consider Gibbs energy as a function of initial composition, temperature, pressure, and other controled parameters. In this worksheet we will build a procedure which builds in turn Gibbs energy as a function of given parameters.
> 
SetTP:= proc(TEL,PEL)
global mu,p,G,F,N,rho,hf,Cp,s;
local GG,PEl,i,P;
description "Sets ii pressure to value pp and refreshes Gibbs Energy Func and gradient of it";
p:= map(y>y*P,x):
for i from 3 to nargs do
p[args[i][1]]:=args[i][2]
end do;
PEl:= PELmultiply(map(y>(y+1) mod 2,var_p),p):
simplify(subs([T=TEL,P=PEl],sum('n[usedp[i]]*(mu[usedp[i]]+R*T*ln(p[usedp[i]]))','i'=1..usedpn))):
G:= optimize(unapply(%,xi)): # Gibbs Energy function
F:= optimize(unapply(simplify([seq(diff(G(xi),xi[k]),k=1..Nxi)]),xi)): # gradient of Gibbs Energy
rho:=optimize(unapply(simplify(subs(P=PEl,[seq(p[usedp[i]],i=1..usedpn)])),xi)) # pressure as function of XIs
end proc: 
As noted above we have inequiality constrains in our system which provide that all pressures of species in the system accept only positive values. So in minimization procedure in which the minimum is iteratively searched and we need some procedure which gives us only phisycally valid values of variables. We need some procedure to determine a termination point of our calculations. Here a simple procedure is build which determines whether there is a minimum in the vicinity of given point within some accuracy or not. It can be shown that our model with gaseous species has only one minimum in the domain of physically valid values of parameters except oxygen which has a fixed potential.
> 
GetValidPnt:= proc(pnt,delta,dir)
local pp,jk,t,d;
description "Returns valid point in dir direction and delta increment.";
t:=true;
d:=delta;
while t do
t:=false;
pp:=pnt;
pp[dir]:= pp[dir]+d;
for jk from 1 to usedpn do
if rho(pp)[jk]<=0 then t:=true end if
end do;
d:=d/2
end do;
return(pp)
end proc:
CheckPrecision:= proc(pnt,eps)
local j,i,epsi,lpnt,hpnt,R;
global F;
epsi:= eps/2;
for i from 1 to Nxi do
lpnt:= GetValidPnt(pnt,eps[i],i);
hpnt:= GetValidPnt(pnt,eps[i],i);
R:= zip((x,y) > sign(x)*sign(y), F(lpnt),F(hpnt));
if sum('R[j]','j'=1..Nxi)=Nxi then return(true) else return(false) end if
end do
end proc: 
Here the main procedure which actually solves the problem  the minimization procedure is introduced. The method used is based on steepest descent method but only one coordinate per step is changed. At each step the procedure moves along the coordinate on which absolute value of descending is maximal. The minimization method used is simple but not efficient. This method is convenient in our case since we need to check positiveness of species partial pressures at each step. So we can control descent along each coordinate independently.
> 
PL:=300: # number of iterations between two stopchecks
CalcPnt:= proc(pnt,d,epsilon,MC)
local pnt0,pnt1,F0,F1,delta,cv,evn,Fs,mx,jk,zeta,i;
global G,F;
delta:=d; pnt0:= pnt;
evn:= NullNxi:
F0:=F(pnt0): Fs:=map(abs,F0);
mx:=1:
for jk from 2 to Nxi do
if Fs[jk]>Fs[mx] then mx:=jk end if
end do:
zeta:= sign(F0[mx]):
for cv from 1 to MC do
if (cv mod PL)=0 then
if CheckPrecision(pnt0,epsilon) then break end if;
if delta[mx]<10^(1Digits) then ERROR("delta is too small, please increase Digits",delta); break end if
end if;
pnt1:=pnt0;
pnt1[mx]:=pnt1[mx]zeta*delta[mx]:
F1:=F(pnt1):
i:=0; # checking quantities  they must be positive
for jk from 1 to usedpn do
if rho(pnt1)[jk]<=0 then i:=i+1 end if
end do;
if i<>0 then
delta[mx]:= delta[mx]/2;
evn[mx]:= 0;
next
end if;
if sign(F0[mx])*sign(F1[mx])<0 then
delta[mx]:=delta[mx]/2;
evn[mx]:=0;
next
end if;
if evn[mx]= 3 then
delta[mx]:=delta[mx]*2;
evn[mx]:=0
end if;
evn[mx]:= evn[mx]+1;
pnt0[mx]:=pnt1[mx];
F0:=F1;
Fs:=map(abs,F0):
mx:=1:
for jk from 2 to Nxi do
if Fs[jk]>Fs[mx] then mx:=jk end if
end do:
zeta:= sign(F0[mx]):
end do:
if cv=MC then WARNING("Number of iterations exceeded given limit "  MC  ", INTERRUPTED!") end if;
return([pnt0,rho(pnt0)])
end proc: 
5. Classic method of the methane reforming (2 reactions) and "BlueBird"
The matrix hereafter corresponds to a matrix of chemical reactions. Each row of this matrix corresponds to given independent reaction and each column of the matrix corresponds to stoichometric coefficients of given component in the reactions regarded. IC hereafter is the list of initial quantities of species.
> 
U:= matrix([[0,3,0,1,1,0,1], # reaction coefficients
[0,1,0,0,1,1,1]]);
IC:=[0,0,0,1,0,0,1]; # initial quantities of compounds
SetReaction(U,IC): 

(5.1) 
Temperature and pressure of the system are also collected in the set.
> 
TE:=950; PE:=3;
SetTP(TE,PE): 

(5.2) 
Now we have the explicit expression of Gibbs energy as a function of only .
The obtained function can be plotted. Let us call this plot as "BlueBird" of equilibrium.
> 
plot3d(G([s1,s2]),s1=0..1,s2=0..0.5,orientation=[125,71],shading=ZHUE,style=patchnogrid,axes=boxed,labels=["xi[1]","xi[2]","G"], title=cat("BlueBird; G(xi[1],xi[2]); [T,P]=[",convert(TE,string),",",convert(PE,string),"]"), titlefont=[COURIER,OBLIQUE,14]); 
Minimization procedure which works iteratively needs an initial point which given here as . Here and mean a level of accuracy of the minimization result and an initial step respectively.
> 
TE:=1000; PE:=1;
SetTP(TE,PE):
pnt0:=[.5,.25]; # initial point
epsilon:= [1e7 $ Nxi]: # precision
delta:=[1e5 $ Nxi]: # initial step 

(5.3) 
Let us define now a list of the temperatures at which the equilibruim will be calculated. In the plots presented below we can see that reaction (2) with formation of has a maximum and one (2) does not have. Therefore, if the aim is to obtain as much as possible then higher temperature for the reforming process has to be chosen.
> 
j:=12:
T_int:=[seq(700+25*i,i=0..j)];
DA:=array(0..j,[]):
for i from 0 to j do
SetTP(T_int[i+1],1):
DA[i]:=CalcPnt(pnt0,delta,epsilon,5000);
end do:
DL:=convert(DA,list): 

(5.4) 
> 
xis:= map(y > y[1],DL):
ps:= map(y > y[2],DL): 
> 
for i from 1 to Nxi do
xi_T[i]:= zip((x,y)>[x,y[i]],T_int,xis)
end do:
xi_T:=convert(xi_T,list):
plot(xi_T,legend=map(y>convert(y,string),['xi[i]' $ 'i'=1..Nxi]), labels=["T","xi_i"], title="xi[i] (T)", titlefont=[COURIER,OBLIQUE,14]); 
> 
for i from 1 to usedpn do
p_T[i]:= zip((x,y) > [x,y[i]],T_int,ps)
end do:
p_T:=convert(p_T,list):
plot(p_T,legend=map(y>names[y],usedp), labels=["T","p[i]"], title="p[i] (T)", titlefont=[COURIER,OBLIQUE,14]); 
6. Partial oxydation of methane under fixed oxygen partial pressure (3 reactions)
The main idea of this section is to consider the participation of molecular oxygen in the reforming system. Since the number of independent reactions (degree of freedom) in the system is given by equation , where is the number of compounds in the system and is the number of chemical elements comprising these compounds, three independent reactions have to be regarded. Let us consider the following set of reactions:
CH4 + H2O = CO + 3H2 (6.1)
CO + H2O = CO2 + H2 (6.2)
H2 + 1/2O2 = H2O (6.3)
> 
U:= matrix([ [0,3,0,1,1,0,1], # reaction coefficients
[0,1,0,0,1,1,1],
[0,1,1/2,0,0,0,1]]);
IC:=[0,0,2,1,0,0,1]; # initial quantities of compounds
SetReaction(U,IC): 

(6.1) 
> 
TE:=1000; PE:=1;
epsilon:= [1e7 $ Nxi]: # precision
delta:=[1e5 $ Nxi]: # initial step
pnt0:= [.5,.25,.2]; # initial point 

(6.2) 
As mentioned above oxygen partial pressure on membrane is regarded as an external controlled parameter. In order to calculate the equilibrium composition depending on given oxygen partial pressure we should define the range of that.
> 
ABint:=[24,15,46]: # [(range in log10 form),number of points]
lgpARR:= evalf([seq(ABint[2](ABint[2]ABint[1])/(ABint[3]1)*(i1),i=1..ABint[3])]):
pARR:= map(y > evalf(10^y),lgpARR): 
> 
i_p:=3: # number of points for interpolation the next one
DA:= array(1..ABint[3],[]):
for i from 1 to ABint[3] do
SetTP(TE,PE,[3,pARR[i]]): # for gases only
DA[i]:= CalcPnt(pnt0,delta,epsilon,5000);
if i>i_p then
for k from 1 to Nxi do
pnt0[k]:=subs(v=lgpARR[i],interp([seq(lgpARR[iz],z=1..i_p)],[seq(DA[iz][1][k],z=1..i_p)],v)) # interpolation
end do;
else
pnt0:=DA[i][1]
end if
end do:
LA:= convert(DA,list):
if nops(LA)<>ABint[3] then ERROR("Some points were NOT calculated. Please check this!!!") end if; 
> 
xis:= map(y > y[1],LA):
ps:= map(y > y[2],LA):
lgpARR:= map(y > log10(y[2]),ps): 
> 
xi_lgP:='xi_lgP':
for i from 1 to Nxi do
xi_lgP[i]:= zip((x,y) > [x,y[i]],lgpARR,xis)
end do:
xi_lgP:=convert(xi_lgP,list):
plot(xi_lgP,legend=map(y>convert(y,string),['xi[i]' $ 'i'=1..Nxi]), labels=["lg(pO2)","xi[i]"], title="xi[i] (lg(p[3]))", titlefont=[COURIER,OBLIQUE,14]); 
> 
p_lgP:='p_lgP':
for i from 1 to usedpn do
p_lgP[i]:= zip((x,y) > [x,y[i]],lgpARR,ps)
end do:
p_lgP:= convert(p_lgP,list):
plot(p_lgP,legend=map(y>names[y],usedp), labels=["lg(p)","p[i]"], title="p[i] (lg(pO2))", titlefont=[COURIER,OBLIQUE,14]); 
7. Partial oxydation of methane with methane cracking CH4 = C + 2H2 (4 reactions)
By analogy with previous section we introduce new reaction in the system in order to take methane cracking, i.e. solid carbon formation, into account. The chemical potential of solid carbon is defined by only its standard part. It is, therefore, useful to redefine the set of reactions in following way:
C + 2H2 = CH4 (7.1)
C + 1/2O2 = CO (7.2)
C + O2 = CO2 (7.3)
H2 + 1/2O2 = H2O (7.4)
> 
U:= matrix([ [1,2,0,1,0,0,0], # reaction coefficients
[1,0,1/2,0,1,0,0],
[1,0,1,0,0,1,0],
[0,1,1/2,0,0,0,1]]);
IC:=[0,0,2,1,0,0,1]; # initial quantities
SetReaction(U,IC): 

(7.1) 
> 
TE:=1000; PE:=1;
epsilon:= [1e7 $ Nxi]: # precision
delta:=[1e5 $ Nxi]: # initial step
pnt0:= [.5,.2,.05,.5]; # initial point 

(7.2) 
> 
ABint:=[22,26,41]; # [(interval in log10 form),number of points]
lgpARR:= evalf([seq(ABint[2](ABint[2]ABint[1])/(ABint[3]1)*(i1),i=1..ABint[3])]):
pARR:= map(y > evalf(10^y),lgpARR): 

(7.3) 
> 
i_p:=3: # number of points for interpolation the next one
DA:= array(1..ABint[3],[]):
for i from 1 to ABint[3] do
SetTP(TE,PE+1,[1,1],[3,pARR[i]]): # for gases and solids
DA[i]:= CalcPnt(pnt0,delta,epsilon,5000);
if i>i_p then
for k from 1 to Nxi do
pnt0[k]:=subs(v=lgpARR[i],interp([seq(lgpARR[iz],z=1..i_p)],[seq(DA[iz][1][k],z=1..i_p)],v)) # interpolation
end do;
else
pnt0:=DA[i][1]
end if
end do:
LA:= convert(DA,list): 
> 
if nops(LA)<>ABint[3] then ERROR("Some points were NOT calculated. Please inspect this!!!") end if; 
> 
xis:= map(y > y[1],LA):
ps:= map(y > y[2],LA):
lgpARR:= map(y > log10(y[3]),ps): 
It can be shown by computations without any restrictions with respect to oxygen partial pressure that the equilibrium in the reactions (7.2) and (7.3) cannot be reached at partial pressures higher than that of 1e22 atm. The computations were, therefore, carried out at oxygen partial pressures lower than this bound.
> 
xi_lgP:='xi_lgP':
for i from 1 to Nxi do
xi_lgP[i]:= zip((x,y) > [x,y[i]],lgpARR,xis)
end do:
xi_lgP:=convert(xi_lgP,list):
plot(xi_lgP,legend=map(y>convert(y,string),['xi[i]' $ 'i'=1..Nxi]), labels=["lg(p[3])","xi[i]"], title="xi[i] (lg(p[3]))", titlefont=[COURIER,OBLIQUE,14]); 
> 
p_lgP:='p_lgP':
for i from 1 to usedpn do
p_lgP[i]:= zip((x,y) > [x,y[i]],lgpARR,ps)
end do:
p_lgP:= convert(p_lgP,list):
plot(p_lgP,legend=map(y>names[y],usedp), labels=["lg(p)","p[i]"], title="p[i] (lg(p[3]))", titlefont=[COURIER,OBLIQUE,14]); 
8. Conclusions
In this worksheet we showed the application of total Gibbs energy minimization method for calculation of chemical equilibrium in the process of methane reforming. It is obvious that presented method can be applied for any chemical system with arbitrary number of independent both homogenous and heterogeneous reactions if appropriate thermodynamic data are available.
References:
[1] D.R. Stull and H. Prophet. "JANAF Thermochemical Tables", U.S. Department of Commerce, Washington, 1985. Cp Fitted by CRCT, Montreal.
Legal Notice: The copyright for this application is owned by the author(s). Neither Maplesoft nor the author are responsible for any errors contained within and are not liable for any damages resulting from the use of this material. This application is intended for noncommercial, nonprofit use only. Contact the author for permission if you wish to use this application in forprofit activities.