1. science
2. physics
3. thermodynamics

Haber Process and Steam-Coal Gasification:
Two Standard Thermodynamic Problems
Elucidated Using Two Distinct Approaches
Housam Binous*1, Ahmed Aheed1,2 and Mohammad M. Hossain1,2
1
Department of Chemical Engineering, 2KACST-TIC for CCS
King Fahd University of Petroleum & Minerals, Dhahran 31261, KSA
*Corresponding
author
Key words: Gibbs free energy minimization, reaction coordinates, equilibrium
constants, Haber process, steam-coal gasification, MATHEMATICA©
Abstract
The present study shows how the direct Gibbs free energy minimization technique is
sometimes superior to the reaction coordinates – equilibrium constants method when
the thermodynamic analysis of complex systems is performed. In this respect, the
above two methods are applied for two different processes to determine their
equilibrium compositions: (1) the Haber process: a simple problem consisting of a
single gas-phase reaction where both methods are expected to give the same result,
and (2) steam-coal gasification process which is more complex since it involves
numerous reactions and solid-phase chemical species (i.e., coal or carbon). Thus, in
the second case the reaction coordinates – equilibrium constants method fails to
provide correct predictions of the equilibrium composition. In addition, the authors
give the optimum conditions for the two processes: high pressure for Haber process
1
and high temperature for steam-coal gasification in agreement with LeChatelier
principle. In order to deal with deviation from the ideal-case assumption, the PengRobinson Equation of State (PR EOS) is implemented for both techniques and both
case studies. All computations and figures are generated using a single computer
algebra, MATHEMATICA©.
Introduction
Thermodynamic analysis calculates the equilibrium distribution of different chemical
species involved in specific system under certain conditions [1] (i.e., temperature,
pressure and/or initial reactants’ mole fractions). Such system is subject to chemical
reactions occurring among the involved chemical species in presence of a catalyst.
Thus, it is very important to know the equilibrium conversions (or mole fractions),
that can be achieved, in order to evaluate the catalyst efficiency and overall
performance.
To perform equilibrium conversions calculations in reactive systems, the general idea
is that the system attains the equilibrium when its total Gibbs free energy is at its
global minimum value [2]. To achieve this, two common methods are usually
employed: (1) reaction coordinates – equilibrium constants method and (2) direct total
Gibbs free energy minimization technique [1,2].
Generally, the direct Gibbs free energy minimization technique is the preferred
method if multiple reactions occur simultaneously. Indeed, this method does not
require neither a chemical reaction(s) set nor a knowledge of the equilibrium
constants. The only thing that is needed is a list of all chemical species involved in the
system as well as a method of computation of their Gibbs free energy. Moreover,
2
reaction coordinates – equilibrium constants method fails to handle solid reactants and
products [1]. Here, we apply both methods for two different case studies: (1) Haber
process and (2) steam-coal gasification. However, instructors can choose other
examples of reaction sets such as biomass gasification, ethane steam cracking,
catalytic cracking of crude oil, and methanol synthesis from syngas. The present
Gibbs free energy minimization approach has been introduced in two graduate level
chemical engineering courses at King Fahd University of Petroleum & Minerals
(KFUPM): (1) CHE 530 or Advanced Reaction Engineering and (2) OGSF 508 or
Industrial Chemical Reactions - Catalysts - R&D Trends. This teaching helps
students: (1) understand the thermodynamic characteristics of industrial chemical
reactions, (2) find the feasibility of a reaction, and (3) determine the displacement of
chemical equilibria as a function of operating conditions (i.e., temperature and
pressure) and reactants initial composition. Students appreciate especially the use of
Mathematica© because it make difficult mathematical problems easily tractable and
allow them to focus on the chemical engineering side of the problem at hand.
Several pedagogical papers have dealt with thermodynamics [3–6] and chemical
reaction engineering [7, 8]. Most of the pedagogical articles considered reaction
coordinates – equilibrium constants methods and/or ideal gas-phase assumption (e.g.,
Ref. 6) to estimate the equilibrium product compositions.
To the best of our
knowledge, only one article has applied the minimization of the Gibbs free energy
technique to determine the equilibrium mole fractions of the coal methanation process
[9]. Although, this article have not compared the prediction between the two after
mentioned approaches. Each method has several advantages over the other. In some
cases, the Gibbs free energy minimization approach shows clear superiority in
accurate prediction. Especially, for the complex reaction system where defining the
3
exact reaction sets is difficult (or in some cases almost impossible). The
thermodynamic equilibrium modeling using Gibbs free energy minimization
technique has become a very popular approach for studying new reactions [10-15].
Keeping the above in mind, the present paper reports a comparison of the two
approaches by considering two case studies.
After a short description of the PR EOS, we provide insights about both resolution
approaches. Then, the two case studies are presented. Finally, we conclude the paper
by giving details about the pedagogical value of the presented work.
The Peng-Robinson Equation of State
Even though the well-known Peng-Robinson EOS is especially suitable for
hydrocarbon systems, it can also handle other common systems (except in case of
sulfur present) containing water, carbon dioxide and monoxide, nitrogen and oxygen
among other chemical species [14]. This EOS can describe both liquid-phase and
vapor-phase behavior. It relates molar volume to pressure and temperature as follows
[15]:
𝑃=
RT
a α(T)
−
V − b V(V + b) + b(V − b)
(1)
Where:
𝑎 = 0.45724
𝑅 2 𝑇𝑐2
𝑅 𝑇𝑐
, 𝑏 = 0.0778
𝑃𝑐
𝑃𝑐
4
(2)
𝛼(𝑇) = [1 + (0.37464 + 1.54226 𝜔 − 0.26992 𝜔
2)
(1 −
1 2
𝑇𝑟2 )]
(3)
and
𝑇𝑟 =
𝑇
𝑇𝑐
(4)
Note that the above equation is valid only for pure component. In the case of
mixtures, the following mixing and combining rules should be used [1]:
𝑃=
𝑅𝑇
𝑎𝑚𝑖𝑥 𝛼(𝑇)
−
𝑉𝑚𝑖𝑥 − 𝑏𝑚𝑖𝑥 𝑉𝑚𝑖𝑥 (𝑉𝑚𝑖𝑥 + 𝑏𝑚𝑖𝑥 ) + 𝑏𝑚𝑖𝑥 (𝑉𝑚𝑖𝑥 − 𝑏𝑚𝑖𝑥 )
(5)
Where 𝑎𝑚𝑖𝑥 and 𝑏𝑚𝑖𝑥 are related to individual components ones by:
𝑛
𝑏𝑚𝑖𝑥 = ∑ 𝑦𝑖 𝑏𝑖 ,
𝑖=1
𝑛
𝑛
𝑎𝑚𝑖𝑥 = ∑ ∑ 𝑦𝑖 𝑦𝑗 𝑎𝑖𝑗
, 𝑎𝑖𝑗 = √𝑎𝑖 𝑎𝑗(1 − 𝑘𝑖𝑗 )
(6)
𝑖=1 𝑗=1
0.45724𝛼𝑖 (𝑇)𝑅2 𝑇𝑐 2𝑖
0.0778 𝑅𝑇𝑐 𝑖
𝑎𝑖 =
, 𝑏𝑖 =
𝑃𝑐 𝑖
𝑃𝑐 𝑖
(7)
2
𝛼𝑖 (𝑇) = [1 + 𝑘𝑖 (1 − √𝑇𝑟 𝑖 )] , 𝑘𝑖 = 0.37464 + 1.54226 𝜔𝑖 − 0.26992 𝜔2𝑖
(8)
And
𝑇
𝑇𝑟 𝑖 = 𝑇
𝑐𝑖
5
(9)
Here, 𝑘𝑖𝑗 , the binary interaction parameters, can be set equal to zero if the molecules
are almost same in size and charge, or if the system approaches ideality, otherwise
they have to be determined either experimentally or by using approximation models
[2].
Reaction Coordinates
Reaction coordinate – equilibrium constant method considers that at equilibrium
(minimum Gibbs free energy point) the differential (or change) of total Gibbs free
energy with respect to the number of moles (of each chemical species involved in the
system) is zero. However, this is true only for constant temperature and pressure
systems. This method can also be applied for flow systems since each point (specific
temperature and pressure) can be dealt with as closed system and because the Gibbs
free energy is a state function which means it does not depend on the path of the
process [1].
The steps of this method are as follow:
a) Selection of complete and independent chemical reactions set. Completeness
means that the chemical reactions’ set should include all chemical compounds
involved in the system. On the other hand, independence requires that each
chemical reaction has its own direction (combinations of reactants and products)
and is written in its overall form. To reduce complexity of this step, very common
and good simplification can be made by ignoring a reaction with either a relatively
high or a relatively small equilibrium constant.
6
Any mistake made in the reactions set selection will result in very different results
from what must be. Completeness problems result from the absence of some
chemical species from the calculations. While independency problems leads to
either no solution or an infinite number of solutions for the studied system.
b) Relating mole fraction of different species to the reaction coordinates of the
selected chemical reactions.
𝑦𝑖 =
𝑛𝑖 0 + ∑𝑗 𝜈𝑖,𝑗 𝜀𝑗
𝑛0 + ∑𝑗 𝜈𝑗 𝜀𝑗
(10)
c) Calculating equilibrium constants for those reactions
𝛥𝐺𝑗 𝑜
𝐾𝑗 = exp (−
)
𝑅𝑇
(11)
d) Relating equilibrium constants to mole fractions, then to reaction coordinates.
−𝜈𝑗
̂ )𝜈𝑖,𝑗 = ( 𝑃 )
∏(𝑦𝑖 𝜙
𝑜
𝑖
𝑃
𝑖
𝐾𝑗
(12)
e) Solving the resulting system of non-linear algebraic equations to evaluate the
values of reaction coordinates. Generally, this step done by using iterative method
due to presence of fugacity coefficients (calculation of fugacity coefficients needs
knowing of the mole fractions and vice versa as stated by the above equation, so
trial and error procedure is needed). A good initial guess for this trial and error
method considers all fugacity coefficients to be equal to unity (i.e., the ideal gas
mixture assumption). In the present work, we make use of the built-in command
of MATHEMATICA© called FindRoot to solve simultaneously the corresponding
system of nonlinear algebraic equations. In addition, since we seek solution for
7
various values of parameters such as temperature, pressure and/or initial reactant
mole fractions, the problem is readily tractable using arc-length continuation. A
short description of this method is given in the following section of the present
paper.
f) The final step is the calculation of the involved chemical species mole fractions
from the found reaction coordinates (back calculations). Note that in several
instances there are more than one solution of such system of nonlinear algebraic
equations. The constraints below allow the determination of the appropriate
solution:
∑ 𝑦𝑖 = 1, 0 ≤ 𝑦𝑖 ≤ 1 ∀𝑖
𝑖
(13)
The advantages of the reaction coordinates – equilibrium constants method resides in
constructing a solvable algebraic system by shifting the system's degree of freedom to
a zero value even though it may not seem to be the case when the set of mole fractions
relations used as is.
Arc-length Continuation
Suppose one has to solve a system of nonlinear algebraic equations which solution
depends on a parameter called 𝛼 in the text thereafter. This parameter must be set first
in order to get the system of equations’ solution. Generally getting the solution of
suchlike systems as a continuous function of this parameter is not easy task because of
nonlinearity and solutions’ multiplicity. Thus, it is advisable to use a parametric
method to see how the solution changes with 𝛼 (especially if one solution is correct
and the others are rejected). An excellent approach consists in the use of the concept
8
of arc-length continuation by letting all of the system's unknowns and the parameter
functions in the arc-length [16], labeled 𝑠 in the text thereafter. In such case, one has
to introduce an additional differential equation called the arc-length auxiliary
equation. Note that the resulting differential-algebraic system needs an initial
condition to be solved. The true solution obtained at a specific value of the parameter,
𝛼, should be used for 𝑠 = 0. All of the above computations are readily achieved using
the built-in MATHEMATICA© commands called FindRoot and NDSolve.
In vector form the system of nonlinear algebraic equation is written as follows [16]:
𝐟(𝐮; α) = 𝟎
(14)
where the components of the vectors u and f are:
𝐟 = {f1 , f𝟐 , f𝟑 , … , fN } ,
𝐮 = {u1 , u𝟐 , … , uN }
(15)
Suppose the system has a solution 𝐮0 when α=α0 , or 𝐟(𝐮𝟎 ; α0 ) = 𝟎, now by letting
𝐮 = 𝐮(s) and α = 𝛼(s) and taking total derivative of f:
d𝐮
dα
d𝐟 = 𝐟𝐮 (𝐮; α) ds + 𝐟α (𝐮; α) ds = 𝟎
du
the unknowns number becomes N+1 as follow: { ds1 ,
du2
ds
,…,
(16)
duN dα
ds
, ds }, but the total
differential has just N independent equation, so the N+1st equation is :
2
d𝐮 d𝐮
dα
∙
+( ) =1
ds ds
ds
In the present study, in order to obtain the equilibrium mole fraction distribution using
the reaction coordinates – equilibrium constants method, one has to solve a system of
nonlinear algebraic equations where either the temperature, the pressure and/or the
initial reactants number of moles is a regarded as the parameter, 𝛼. This makes the
arc-length continuation technique a method of choice when one tackles such
problems.
9
(17)
Gibbs-Free Energy Minimization
Direct total Gibbs free energy minimization can be performed by using any normal
optimization method [2]. Here, the objective function that should be minimized is the
total Gibbs free energy, and since it depends on each chemical species number of
moles, the number of moles must be taken as decision variables. It should be noted
that the number of moles are subject to constraints steaming from the various atomic
mass balance equations. It is in these above constraints where the initial reactants
mole fractions come into play.
So the general form of our optimization problem is [1]:
𝑚𝑖𝑛 𝐺 𝑡𝑜𝑡𝑎𝑙 ({𝑛𝑖 })
Subjected to:
∑ 𝑛𝑖 𝑎𝑖𝑘 = 𝐴𝑘
𝑖
(18)
So the optimum (minimum) point can be obtained analytically (by using Lagrangian
multipliers since all of constraints are of equality ones), or numerically by using any
suitable
algorithm
[17]
or
the
MATHEMATICA© built-in
command
called
FindMinimum.
In case of using the Lagrangian multipliers, the equations to be solved are [2]:
∑ 𝑛𝑖 𝑎𝑖𝑘 = 𝐴𝑘
(𝑘 = 1, 2, 3, … , 𝑤)
𝑖
̂𝑃
𝑦𝜙
𝐺𝑜𝑖 + 𝑅𝑇𝑙𝑛 ( 𝑖 𝑜𝑖 ) + ∑ 𝜆𝑘 𝑎𝑖𝑘 = 0
𝑃
𝑘
10
(𝑖 = 1, 2, 3, … , 𝑁)
(19)
(20)
Where 𝜆𝑘 are the Lagrangian multipliers. Note also that there is a need of an iterative
method (due to fugacity coefficients/number of moles dependency) when using this
method. Moreover, this method gives all of the critical points (maxima, minima and
saddle points) which adds another step to the calculations (i.e., the determination the
global minimum). Because of this, using one of the numerical minimization
algorithms designed to converge to the minimum point directly is the preferred
approach. Regarding numerical algorithms, any interior or exterior point algorithm
can be applied. However, in case of using interior point algorithm, one must insure
that the initial guess is feasible. Thus, to overcome this difficulty, it is suitable to use
exterior point algorithm.
Total Gibbs free energy for all system is mainly the summation of each single-phase
total Gibbs free energy existing in that system [2]. In this respect, it can be divided
into two main parts:
i.
Gibbs free energy for gas and liquid phases (because two of them calculated by
applying the PR EOS and taking different roots for the molar volume) and
ii.
Gibbs free energy of solid phase
a. Gas and Liquid phase mixture Gibbs free energy:
Firstly, the liquid-phase mixture can be dealt with as if it is gas-phase with change the
value of molar volume only. This the great advantage of the cubic equations of state
such as the PR EOS. It should be pointed out that separate calculations of each
species’ fugacity can be avoided by using the concept of mixing and combining rules
(𝑎𝑚𝑖𝑥 and 𝑏𝑚𝑖𝑥 ) and using the residual property relation [14]:
11
𝑖𝑔
𝐺𝑔𝑎𝑠,𝑚𝑖𝑥 = 𝐺𝑚𝑖𝑥 + 𝐺𝑅𝑚𝑖𝑥
(21)
The ideal gas mixture Gibbs free energy is given by:
𝑛
𝑖𝑔
𝐺𝑚𝑖𝑥
=
𝑛
𝑖𝑔
∑ 𝑦𝑖 𝐺𝑖
𝑖=1
+ 𝑅 𝑇 ∑ 𝑦𝑖 𝑙𝑛𝑦𝑖
(22)
𝑖=1
The ideal gas Gibbs free energy of pure substance is given by:
𝑖𝑔
𝑇
𝐺𝑖 = ∫
𝑇
𝐶𝑝,𝑖
𝑑𝑇 +
298.15 𝑇
𝐶𝑝,𝑖 𝑑𝑇 − 𝑇 ∫
298.15
𝑅 𝑇𝑙𝑛
𝑃
− 𝑆𝑖 298.15 (𝑇 − 298.15) + 𝐺𝑖 298.15
𝑃1
(23)
where 298.15 K is the reference temperature while 𝑃1 is the reference pressure.
One can notice that, all of the terms on the right hand side of the equation (23) are
functions in temperature except one, which has also pressure-dependency.
Regarding residual Gibbs free energy, it is the very same formula for pure substances
that is used but with some modifications made in the EOS constants by using mixing
and combining rules shown below.
For pure components and by using Peng-Robinson EOS [14]:
𝐺𝑅
𝑃𝑏
𝑎
𝑉 + (1 + √2)𝑏
= 𝑍 − 1 − ln (𝑍 − ) −
ln (
)
𝑅𝑇
𝑅𝑇
2√2𝑏𝑅𝑇
𝑉 + (1 − √2)𝑏
Thus, for mixtures, we find what follows:
𝐺𝑅𝑚𝑖𝑥
𝑃𝑏𝑚𝑖𝑥
𝑎𝑚𝑖𝑥
𝑉𝑚𝑖𝑥 + (1 + √2)𝑏𝑚𝑖𝑥
= 𝑍𝑚𝑖𝑥 − 1 − ln (𝑍𝑚𝑖𝑥 −
ln (
)−
)
𝑅𝑇
𝑅𝑇
2√2𝑏𝑚𝑖𝑥 𝑅𝑇
𝑉𝑚𝑖𝑥 + (1 − √2)𝑏𝑚𝑖𝑥
12
(24)
(25)
𝑍𝑚𝑖𝑥 =
𝑃𝑉𝑚𝑖𝑥
𝑅𝑇
(26)
where 𝑉𝑚𝑖𝑥 is determined from the roots of the cubic PR EOS. It is worth mentioning
at this point that the relation between Gibbs free energy and number of moles clearly
comes from the mole fractions, which appear in molar Gibbs free energy:
𝑦𝑖 =
𝑛𝑖
𝑛𝑇
(27)
On the other hand, in the reaction coordinates – equilibrium constants method, one
needs to calculate fugacity coefficients from the PR EOS as follows [14, 18]:
̂ = exp ((𝑧 − 1)
𝜙
𝑖
𝐵𝑖
− 𝑙𝑛(𝑧 − 𝐵)
𝐵
(28)
−
2
𝐴
(
2√2𝐵
∑𝑛𝑐
𝑗=1 𝑦𝑖
𝐴𝑖𝑗
𝐴
𝐴𝑖 =
−
𝐵𝑖
𝑧 + (1 + √2)𝐵
) 𝑙𝑛 (
))
𝐵
𝑧 + (1 − √2)𝐵
𝑎𝑖 𝑃
𝑅𝑇
2
, 𝐵𝑖 =
𝑏𝑖 𝑃
𝑅𝑇
𝐴𝑖𝑗 = √𝐴𝑖 𝐴𝑗 (1 − 𝑘𝑖𝑗 )
It should be clear that there are three values for z or v, which result from solving the
PR EOS at specific temperatures and pressures. Either two of them have complex
conjugate values (because the polynomial coefficients are real ones), in such case,
there is only one phase existing in the system which is either the vapor or liquid
phase. The other case where all values are real, in such case the largest value
corresponds to the vapor-phase molar volume while the smallest one is its liquidphase counterpart and this can be proved by stability tests [2].
13
(29)
(30)
b. Solid Phase Gibbs free energy:
Fortunately, total Gibbs free energy of solid phase can be set equal to the sum of each
total Gibbs free energies for all components existing in such phase due to negligible
solid-solid interactions.
𝑛𝑡𝑜𝑡𝑎𝑙,𝑠𝑜𝑙𝑖𝑑 𝐺𝑠𝑜𝑙𝑖𝑑,𝑚𝑖𝑥 = ∑ 𝑛𝑖,𝑠𝑜𝑙𝑖𝑑 𝐺𝑖,𝑠𝑜𝑙𝑖𝑑
(31)
𝑖
In this respect, Gibbs free energies of solid species (reactants and products such as
carbon or coal in the second example) are calculated and modeled as a function of
temperature only (since solid phase is not affected by pressure remarkably).
Thus the total Gibbs free energy of the system is:
𝐺𝑡𝑜𝑡𝑎𝑙 = 𝑛𝑡𝑜𝑡𝑎𝑙,𝑔𝑎𝑠 𝐺𝑔𝑎𝑠,𝑚𝑖𝑥 + 𝑛𝑡𝑜𝑡𝑎𝑙,𝑠𝑜𝑙𝑖𝑑 𝐺𝑠𝑜𝑙𝑖𝑑,𝑚𝑖𝑥
(32)
The two Case Studies
a. Haber process:
The Haber process produces ammonia gas by using hydrogen and nitrogen gases as
main reactants and iron as catalyst according to the following reaction:
3𝐻2 + 𝑁2 ⇆ 2𝑁𝐻3
14
(r.1)
This chemical process is easy enough (i.e., a single reaction in the gas-phase) to be
completely studied and understood thermodynamically. A single gas-phase reaction
implies neither dependency nor completeness problems are expected when the
reaction coordinates – equilibrium constants method is utilized. In addition, no solid
material is involved in such system, so results obtained by performing the reaction
coordinates – equilibrium constants method should coincide with those obtained by
using direct Gibbs free energy minimization technique.
a.1 The reaction coordinates – equilibrium constants method applied to
the Haber process
Firstly, by simple considerations, all of mole fractions are related to the extent of the
reaction as follow:
𝑦𝐻2 =
𝑛𝐻2 0 − 3𝜀
(33)
𝑛𝑁2 0 − 𝜀
(34)
𝑛𝑁𝐻3 0 + 2𝜀
(35)
𝑛0 − 2𝜀
𝑦𝑁2 =
𝑦𝑁𝐻3 =
𝑛0 − 2𝜀
𝑛0 − 2𝜀
Since hydrogen and nitrogen are elements, their Gibbs of formation energies are
zeros [1], so Gibbs energy change for this reaction is simply that one for ammonia
formation (under same conditions). This allows us to write an expression for the
single equilibrium constant involved in the process:
15
𝐾 = exp (−
𝛥𝐺𝑁𝐻3 𝑓
𝑅𝑇
)
(36)
Now, by substituting mole fractions in terms of reaction extent or coordinate:
2
𝑛𝑁𝐻3 0 + 2𝜀
̂
(
𝑛0 − 2𝜀 𝜙𝑁𝐻3 )
𝛥𝐺𝑁𝐻3𝑓
𝑃
= ( ) exp (−
)
1
𝑃𝑜
𝑅𝑇
2
3
𝑛𝐻2 0 − 3𝜀
𝑛𝑁2 0 − 𝜀
̂
̂
∗
𝜙
(
(
)
𝐻2
𝑛0 − 2𝜀
𝑛0 − 2𝜀 ∗ 𝜙𝑁2 )
(37)
This single equation is easily tractable using a trial and error method due to fugacity
coefficients dependency on mole fractions. In the present work, we choose to solve
simultaneously the system nonlinear algebraic of equations consisting of the above
equation, the three definitions of fugacity coefficients and the cubic PR EOS.
Moreover, we apply arc-length continuation in order to determine the solutions at
various pressures and temperatures (see Figure 1 with a block diagram explaining the
solution steps using MATHEMATICA©).
a.2. Direct Gibbs free energy minimization applied to the Haber process
The problem statement when this method is used is the following:
𝑚𝑖𝑛 𝐺 𝑡𝑜𝑡𝑎𝑙 (𝑛𝐻2 , 𝑛𝑁2 , 𝑛𝑁𝐻3 )
Subjected to:
𝑛𝐻2
𝐴
2
0
3
[
] [ 𝑛𝑁 ] = [ 𝐻 ]
𝐴𝑁
0 2 1 𝑛𝑁𝐻2
3
𝑛𝐻2 , 𝑛𝑁2 and 𝑛𝑁𝐻3 ≥ 0
16
(38)
Where:
2𝑛𝐻2 0 + 3 𝑛𝑁𝐻3 0
𝐴𝐻
]
]=[
𝐴𝑁
2 𝑛𝑁2 0 + 𝑛𝑁𝐻3 0
[
(39)
Since the reaction is taking place in the gas-phase, we have:
𝐺𝑡𝑜𝑡𝑎𝑙 = 𝑛𝑡𝑜𝑡𝑎𝑙,𝑔𝑎𝑠 𝐺𝑔𝑎𝑠,𝑚𝑖𝑥
Figure 2 gives the block diagram for the solution steps involved in this calculation
using MATHEMATICA©. Figure 3 shows the mole fractions of hydrogen, nitrogen and
ammonia. The solid lines represents the mole fractions calculated by equilibrium
constants – reaction coordinates method using arc-length continuation while the
rhombus shape points refers to the results obtained by applying Gibbs free energy
minimization. As expected and explained previously, it is clear that both approaches
give the same results. Using the built-in MATHEMATICA© command, Manipulate,
it is possible for user of the program to vary the temperature and observe its effect
instantaneously on the graph. As expected, it is found that ammonia synthesis is
favored by high pressures (i.e., up to a few hundred bars) as well as low temperatures.
b. Steam-Coal Gasification
Steam-coal gasification process is much more complicated than the Haber process.
All limitations of the reaction coordinates – equilibrium constants method apply.
Indeed, there is a plethora of reactions and a solid-phase component is involved in the
reaction scheme.
Table 1 shows the chemical reactions taking place in the steam-coal gasification
process [19].
17
(40)
b.1. The reaction coordinates – equilibrium constants method applied to
the Steam-Coal Gasification
If all reactions are taken into account, a dependency problem will arise as explained
below. One can naively consider the system of equations given below:
̂ 𝑦 𝜙
̂
𝑦𝐶𝑂 𝜙
𝐶𝑂 𝐻2 𝐻2
̂
𝑦𝐻2𝑂 𝜙
𝐻2 𝑂
= 𝐾2
(41)
2
̂
𝑦2𝐶𝑂 𝜙
𝐶𝑂
̂
𝑦𝐶𝑂2 𝜙
𝐶𝑂2
̂
𝑦𝐶𝐻4 𝜙
𝐶𝐻4
2
̂
𝑦2𝐻2 𝜙
𝐻2
2
= 𝐾3
(42)
= 𝐾4
(43)
2
̂ 𝑦2 𝜙
̂
𝑦2𝐻2 𝜙
𝐻2 𝐶𝑂 𝐶𝑂
̂ 𝑦 𝜙
̂
𝑦𝐶𝑂2 𝜙
𝐶𝑂2 𝐶𝐻4 𝐶𝐻4
= 𝐾5
(45)
3
̂ 𝑦 𝜙
̂
𝑦3𝐻2 𝜙
𝐻2 𝐶𝑂 𝐶𝑂
̂
̂
𝑦𝐻2𝑂 𝜙
𝐻2 𝑂 𝑦𝐶𝐻4 𝜙𝐶𝐻4
̂ 𝑦 𝜙
̂
𝑦𝐶𝑂2 𝜙
𝐶𝑂2 𝐻2 𝐻2
̂ 𝑦
̂
𝑦𝐶𝑂 𝜙
𝐶𝑂 𝐻 𝑂 𝜙𝐻2 𝑂
2
18
= 𝐾6
(46)
= 𝐾7
(47)
However, this set of equations has no solution as a result of the conflict created by
substituting equilibrium relations for chemical reactions (r.2) and (r.4) in the relation
for reaction (r.6), which does not lead to
𝐾2
𝐾4
= 𝐾6 as expected. In addition, one can
note that this system consists of 5 unknowns (i.e., the five equilibrium mole fractions)
while there are six equations. Thus, the problem is over-determined (i.e., we have a
negative degree of freedom) and not all the equations are independent from each
other. Thus, a carefully selected chemical reaction set must be used, one can propose
the following set:
𝐶 + 2𝐻2 𝑂 ⇆ 𝐶𝑂2 + 2𝐻2
(r.8)
𝐶 + 𝐻2 𝑂 ⇆ 𝐶𝑂 + 𝐻2
(r.9)
𝐶 + 2𝐻2 ⇆ 𝐶𝐻4
(r.10)
According to this proposed system, we have what follows:
2
̂ 𝑦2 𝜙
̂
𝑦𝐶𝑂2 𝜙
𝐶𝑂2 𝐻2 𝐻2
2
̂
(𝑦𝐻 𝑂 𝜙
𝐻 𝑂)
2
= 𝐾8
(48)
2
̂ 𝑦 𝜙
̂
𝑦𝐶𝑂 𝜙
𝐶𝑂 𝐻2 𝐻2
̂
𝑦𝐻2𝑂 𝜙
𝐻2 𝑂
̂
𝑦𝐶𝐻4 𝜙
𝐶𝐻4
2
̂
𝑦2𝐻2 𝜙
𝐻2
19
= 𝐾9
(49)
= 𝐾10
(50)
The dependency versus temperature of 𝐾8 , 𝐾9 and 𝐾10 are shown in Figure 4. It is
clear from the trend of these curves that reactions (r.8) and (r.9) are endothermic but
reaction (r.10) is exothermic (i.e., 𝐾10 for reaction (r.10) decreases as 𝑇 increases).
Also, the equilibrium constant, 𝐾9 , for reaction (r.9) is larger than 𝐾8 for reaction (r.8)
at higher temperatures. Thus, one expects methane, water, and carbon dioxide to be
depleted at higher temperatures in favor of carbon monoxide and hydrogen. Also, a
higher molar ratio 𝑟 of water to carbon shifts reaction (r.8) to the product side, with
more carbon dioxide being formed.
The relations among mole fractions and reactions' extents are:
𝑛𝐻2𝑂 0 − 2𝜀8 − 𝜀9
(51)
𝑛𝐻2 0 + 2𝜀8 + 𝜀9 − 2𝜀10
(52)
𝑦𝐻 2 𝑂 =
𝑦𝐻2 =
𝑛0 + 𝜀8 + 𝜀9 − 𝜀10
𝑛0 + 𝜀8 + 𝜀9 − 𝜀10
𝑛𝐶𝑂2 0 + 𝜀8
𝑦𝐶𝑂2 =
𝑛0 + 𝜀8 + 𝜀9 − 𝜀10
𝑦𝐶𝑂 =
𝑛𝐶𝑂 0 + 𝜀9
𝑛0 + 𝜀8 + 𝜀9 − 𝜀10
𝑦𝐶𝐻4 =
𝑛𝐶𝐻4 0 + 𝜀10
𝑛0 + 𝜀8 + 𝜀9 − 𝜀10
Now, by substituting equations 51 through 55 into equations 48 through 50, the
system will be have zero degree of freedom (three unknowns and three independent
equations). The unknowns are now the three extents of reactions. This may lead to
more than one solution and one has to select the solution that provides mole fractions
laying in the interval [0,1]. Results obtained by applying such approach and arc20
(53)
(54)
(55)
length continuation are displayed in Figure 5. Using the built-in MATHEMATICA©
command, Manipulate, it is possible for user of the program to vary the pressure
and observe it effect instantaneously on the graph. As expected, it is found that
syngas (i.e., hydrogen and carbon monoxide) production is favored by low pressures
(i.e., atmospheric pressure) as well as high temperatures. In Figure 6 the extents of
reaction for the three reactions (i.e. 𝜀8 , 𝜀9 and 𝜀10 ) are plotted.
b2. Direct Gibbs free energy minimization applied to the Steam-Coal
Gasification
Based on the above description, one has to solve the following optimization problem:
𝑚𝑖𝑛 𝐺 𝑡𝑜𝑡𝑎𝑙 (𝑛𝐻2 𝑂 , 𝑛𝐻2 , 𝑛𝐶𝑂2 , 𝑛𝐶𝑂 , 𝑛𝐶𝐻4 , 𝑛𝐶 )
Subjected to:
2
[1
0
𝑛𝐻2𝑂
𝑛𝐻
𝐴𝐻
2 0 0 4 0 𝑛 2
𝐶𝑂2
= [𝐴𝑂 ]
0 2 1 0 0] 𝑛
𝐶𝑂
𝐴𝐶
0 1 1 1 1 𝑛
𝐶𝐻4
[ 𝑛𝐶 ]
(51)
𝑛𝐻2𝑂 , 𝑛𝐻2 , 𝑛𝐶𝑂2 , 𝑛𝐶𝑂 , 𝑛𝐶𝐻4 and 𝑛𝐶 ≥ 0
where:
2 𝑛𝐻2𝑂 0 + 2𝑛𝐻2 0 + 4𝑛𝐶𝐻4 0
𝐴𝐻
[𝐴𝑂 ] = [ 𝑛𝐻2 𝑂 0 + 2 𝑛𝐶𝑂2 0 + 𝑛𝐶𝑂 0 ]
𝐴𝐶
𝑛𝐶𝑂2 + 𝑛𝐶𝑂 0 + 𝑛𝐶𝐻4 + 𝑛𝐶 0
0
0
21
(52)
It is obvious that there is only one component in this system appears in solid state
which is carbon namely. Thus, we have:
𝐺𝑡𝑜𝑡𝑎𝑙 = 𝑛𝑡𝑜𝑡𝑎𝑙,𝑔𝑎𝑠 𝐺𝑔𝑎𝑠,𝑚𝑖𝑥 + 𝑛𝐶 𝐺𝐶
Figure 7 displays the results of such minimization procedure (i.e., the water-free
equilibrium mole fractions, for methane, carbon monoxide, hydrogen and carbon
dioxide, versus temperature). Again Using the built-in MATHEMATICA© command,
Manipulate, it is possible for user of the program to vary either the pressure or the
steam/coal initial ratio and observe their effect instantaneously on the graph. This is a
very useful pedagogical feature that instructors can use in the classroom. Indeed,
recalculation for different parameters and display of the results takes only a few
seconds. As expected, it is found that syngas (i.e., hydrogen and carbon monoxide)
production is favored by low pressures (i.e., atmospheric pressure) as well as high
temperatures. The result using this methodology is far more accurate than the
previous one (i.e., using reaction coordinates – equilibrium constants method).
Indeed, the calculation using the Gibbs free energy takes into account more accurately
the solid-phase reactant (i.e., coal).
Conclusion
From the present study the following list of conclusion is drawn:
1. The reaction coordinates – equilibrium constants method requires clearly
defined reactions incorporating all the species involved in the reactions.
22
(53)
2. The Gibbs free energy minimization technique does not require the
stoichiometric reactions. This approach relies on thermodynamic databases
that contain the values of the standard Gibbs energy of the components.
3. The reaction coordinates – equilibrium constants method always considers
the solid phase reactant activity as unity.
4. The Gibbs free minimization approach can take into account solid reactants
such as coal. Therefore, this method represents the system more accurately.
It is worth, at this point, stressing the importance of this present pedagogical
contribution to the teaching of chemical engineering thermodynamics. Indeed, we find
equilibrium mole fractions for two important processes that are often taught in senior
undergraduate and/or first-year graduate level courses. We use a common technique
that is ubiquitous in most chemical engineering thermodynamic courses: the reaction
coordinates – equilibrium constants method. However, we go beyond what is usually
considered by studying high-pressure effects through the incorporation of the fugacity
coefficients using the PR EOS. In addition, we use the Gibbs free minimization
technique, which is less commonly taught in undergraduate courses despite its
versatility. Indeed, we see in our second case study (i.e., the steam-coal gasification
process) that it is easily implemented in MATHEMATICA© using the built-in command
called FindMinimum. In addition, this method can handle solid-phase reactants such
as coal and requires no prior knowledge of the numerous chemical reactions taking
place in the process. Finally, all the MATHEMATICA© codes, used in the presented
theoretical calculations, are available upon request from the corresponding author or
at
the
Wolfram
Demonstration
(http://demonstrations.wolfram.com/index.html).
23
Project
Acknowledgement
The research team acknowledges the financial support provided by King AbdulAziz
City for Science and Technology (KACST) to this research under KACST-TIC for
CCS project no 03. The team also acknowledges the facilities and support provided by
KFUPM.
References
[1]
J. M. Smith, H. C. Van Ness, and M. M. Abbot, Introduction to chemical
engineering thermodynamics (7th ed.), McGraw-Hill Higher Education, New
York, 2005.
[2]
J. W. Tester and M. Modell, Thermodynamics and Its Applications (3rd ed.),
Prentice Hall, Upper Saddle River, NJ, 1996.
[3]
R. Baur, J. Bailey, B. Brol, A. Tatusko, and R. Taylor, Maple and the art of
thermodynamics, Comput Appl Eng Educ, 6 (1998), pp. 223–234.
[4]
Y. Liu, Development of instructional courseware in thermodynamics
education, Comput Appl Eng Educ, 19 (2011), pp. 115–124.
[5]
F. Cruz‐Peragón, J. M. Palomar, E. Torres‐Jimenez, and R. Dorado,
Spreadsheet for teaching reciprocating engine cycles, Comput Appl Eng Educ,
20 (2012), pp. 681–691.
[6]
K. T. Klemola, “Chemical reaction equilibrium calculation task for chemical
engineering undergraduates – simulating Fritz Haber’s ammonia synthesis with
thermodynamic software,” Chem Eng Educ, 48 (2014), pp. 115–120.
24
[7]
L. M. Porto and R. H. Ogeda, Java applets for chemical reaction engineering,
Comput Appl Eng Educ, 6 (1998), pp. 67–77.
[8]
R. Rodrigues, R. P. Soares, and A. R. Secchi, Teaching chemical reaction
engineering using EMSO simulator, Comput Appl Eng Educ, 18 (2010), pp.
607–618.
[9]
A. L. Myers, “Computation of Multiple Reaction Equilibria,” Chem Eng Educ,
25 (1991), pp. 112–116.
[10] J. Mazumder and H. I. de Lasa, “Fluidizable La2O3 promoted Ni/γ-Al2O3
catalyst for steam gasification of biomass: Effect of catalyst preparation
conditions,” Appl. Catal. B Environ, 168 (2015), pp. 250–265.
[11] Y. Cao, W-P. Pan, “Investigation of Chemical Looping Combustion by Solid
Fuels. 1. Process Analysis”, Energy & Fuels, 20 (2006), pp. 1836-1844.
[12] M. Asadullah, “Barriers of commercial power generation using biomass
gasification gas: A review,” Renew. Sustain. Energy Rev., 29 (2014), pp. 201215.
[13] H. I. de Lasa, E. Salaices, J. Mazumder, and R. A. Lucky, “Catalytic Steam
Gasification of Biomass: Catalysts, Thermodynamics and Kinetics, Chem. Rev.
111(2011) 5404-5433.
[14] S. I. Sandler, Chemical and engineering thermodynamics (3 rd ed.), Wiley, NY,
1999.
25
[15] D. Y. Peng and D. B. Robinson, “A New Two-Constant Equation of State,”
Indust. and Eng. Chemistry: Fundamentals, 15 (1976), pp. 59–64.
[16] H. Binous and A. A. Shaikh, “Introduction of the arc-length continuation
technique in the chemical engineering graduate program at KFUPM,” Comput
Appl Eng Educ, 23 (2015), pp. 344–351.
[17] M. S. Bazaraa, H. D. Sherali, and C. M. Shetty, Nonlinear Programming:
Theory and Algorithms (3rd ed.), John Wiley & Sons, New York, 2006.
[18] H. Binous, “Applications of the Peng-Robinson Equation of State Using
Mathematica,” Chem Eng Educ, 42 (2008), pp. 47–51,
[19] E. Salaices, B. Serrano, and H. De Lasa, “Biomass Catalytic Steam
Gasification Thermodynamics Analysis and Reaction Experiments in a CREC
Riser Simulator,” Ind. Eng. Chem. Res., 49 (2010), pp. 6834–6844.
26
Nomenclature
𝑎𝑖𝑘
Number of atoms of kth element involved in ith component
𝐴𝑘
Total number of atomic masses of kth element in the system
𝐶𝑝,𝑖
ith component constant pressure heat capacity (J/mol*K)
𝐺𝑖 298.15
Gibbs free energy of ith component at 298.15 K (J/mol)
𝐺𝑖𝑜
Standard Gibbs free energy for ith component (J/mol)
𝑖𝑔
𝐺𝑚𝑖𝑥
Ideal gas Gibbs free energy of the mixture (J/mol)
𝑅
𝐺𝑚𝑖𝑥
Residual Gibbs free energy of the mixture (J/mol)
𝛥𝐺𝑗 𝑜
jth chemical reaction Gibbs free energy (J/mol)
𝑘𝑖𝑗
Binary interaction parameter between ith and jth components
𝐾𝑗
jth chemical reaction equilibrium constant
𝑛𝑇
Total number of moles in the system
𝑛𝑖
Number of moles of ith component
𝑛0
Initial total number of moles
𝑛𝑖 0
Initial number of moles of ith component
𝑃
Pressure (bar)
𝑃𝑐
Critical pressure (bar)
𝑃𝑜
Standard pressure (bar)
𝑅
Universal gas constant (J/mol*K)
𝑆𝑖 298.15
ith component entropy at 298.15 K (J/mol*K)
𝑇
Temperature (K)
𝑇𝑐
Critical temperature (K)
𝑇𝑟
Reduced temperature
𝑉
Molar volume (cm3/mol)
𝑦𝑖
Mole fraction of component i
27
𝑍
Compressibility factor
Greek Letters
𝜀𝑗
Extent of jth chemical reaction
𝜆𝑘
kth Lagrange multiplier
𝜈𝑖,𝑗
Stoichiometric coefficient of ith component in jth chemical reaction
𝜙̂𝑖
ith component fugacity coefficient
𝜔
Acentric factor
Abbreviations
EOS
Equation of State
PR
Peng-Robinson
28
Figure 1 – Algorithm for the solution of the Haber process using
the arc-length continuation technique and the reaction coordinate method
29
Figure 2 – Algorithm for the solution of the Haber process
using the Gibbs free energy minimization technique
30
Figure 3 – Comparison of hydrogen, nitrogen and ammonia mole fractions calculated
by Gibbs free energy minimization technique and
the reaction coordinates – equilibrium constants method.
31
Figure 4 – Equilibrium constants for reactions r.8, r.9 and r.10
32
Figure 5 – Equilibrium mole fractions of Methane, carbon monoxide, hydrogen, and
carbon dioxide on a water-free basis calculated by reaction coordinates – equilibrium
constants method and arc-length continuation.
33
Figure 6 – Extent of reactions for reactions r.8, r.9 and r.10 calculated by reaction
coordinates – equilibrium constants method and arc-length continuation.
34
Figure 7 – Equilibrium mole fractions of Methane, carbon monoxide, hydrogen, and
carbon dioxide on a water-free basis calculated
by direct Gibbs free energy minimization.
35
∆𝐆𝐟𝐨 (𝟐𝟗𝟖𝐊)
∆𝐇𝐨𝐟 (𝟐𝟗𝟖𝐊)
[𝐤𝐉/𝐦𝐨𝐥𝐞]
[𝐤𝐉/𝐦𝐨𝐥𝐞]
Chemical Reaction
K (800℃)
𝐂 + 𝐇𝟐 𝐎 ⇆ 𝐇𝟐 + 𝐂𝐎
(r.2)
89.824
130.414
7.0401
𝐂 + 𝐂𝐎𝟐 ⇆ 𝟐𝐂𝐎
(r.3)
118.362
172.615
6.499
𝐂 + 𝟐𝐇𝟐 ⇆ 𝐂𝐇𝟒
(r.4)
-50.273
-74.9
0.049
𝐂𝐇𝟒 + 𝐂𝐎𝟐 ⇆ 𝟐𝐂𝐎 + 𝟐𝐇𝟐
(r.5)
168.635
123.76
132.013
𝐂𝐇𝟒 + 𝐇𝟐 𝐎 ⇆ 𝐂𝐎 + 𝟑𝐇𝟐
(r.6)
140.098
205.31
169.182
𝐂𝐎 + 𝐇𝟐 𝐎 ⇆ 𝐇𝟐 + 𝐂𝐎𝟐
(r.7)
-28.538
-42.2
1.0051
Table 1 – List of the reactions involved in the steam-coal gasification process
36