Gibbs Energy Minimization Using Simulated Annealing for Two-phase Equilibrium Calculations in ...

14 Pages·2008·0.38 MB·English

A.BONILLA-PETRICIOLETetal.,Gibbs EnergyMinimizationUsing Simulated…,Chem. Biochem.Eng. Q.22 (3)285–298 (2008) 285
GibbsEnergyMinimizationUsingSimulatedAnnealing
forTwo-phaseEquilibriumCalculationsinReactiveSystems
A. Bonilla-Petriciolet,a,* U. I. Bravo-Sánchez,a F. Castillo-Borja,a S. Frausto-Hernández,a
and J. G. Segovia-Hernándezb
aInstituto Tecnológico de Aguascalientes, Departamento de Ingeniería Química
Av. López Mateos 1801, CP 20256, Aguascalientes, Ags. México Original scientific paper
bUniversidad de Guanajuato, Facultad de Ciencias Químicas, Received: April 11, 2007
Noria Alta s/n, C.P. 36050, Guanajuato, Gto. México Accepted: February 28, 2008
Phaseequilibriumcalculationsinsystemssubjecttochemicalreactionsareinvolved
in the design, synthesis and optimization of reactive separation processes. Until now,
several methods have been developed to perform simultaneously physical and chemical
equilibrium calculations. However, published methods may face numerical difficulties
such as variable initialization dependence, divergence and convergence to trivial solu-
tions or unstable equilibrium states. Besides, these methods generally use conventional
composition variables and reactions extents as unknowns which directly affect the nu-
mericalimplementation,reliabilityand efficiency of solving strategies. The objective of
this work is to introduce and test an alternative approach to perform Gibbs energy
minimization in phase equilibrium problems for reactive systems. Specifically, we have
employedthe transformed compositionvariables of Ung and Doherty and the stochastic
optimization method Simulated Annealing for two-phase equilibrium calculations in re-
acting systems. Performance of this strategy has been tested using several benchmark
problems and results show that proposed approach is generally suitable for the global
minimization of transformed Gibbs energy in reactive systems with two-phase equilib-
rium.
Key words:
Global optimization,Gibbs energyminimization,simulatedannealing, chemicalequilib-
rium, phase equilibrium
Introduction consequence, there are frequently computational
difficulties in these calculations and published
Phase equilibrium calculations in reacting sys-
methods could not be reliable generally. Specifi-
tems are often required for designing processes of
cally, there are initialization troubles, the presence
the chemical, petrochemical and metallurgical in-
of trivial solutions or local minimums of Gibbs en-
dustries. During the last years, there has been a
ergy and numerical methods may also present slow
growing interest for developing new numerical
convergence or divergence. It is important to note
tools to model the thermodynamic behavior of mix-
that few global solving methods have been pro-
tures under physical and chemical equilibrium.1
posed for this area.2,7,9,10 Specifically, McDonald
Until now, several methods have been proposed to
and Floudas2 proposed a deterministic global opti-
perform reactive phase equilibrium calculations.1–12
mization method which guarantees the global
Principally, these methods have been developed to
minimization of Gibbs energy in reacting mixtures
model reactive distillation process and they are
using solution models. Also, Jalali and Seader7
based on equilibrium constant (K-value) method or
have successfully applied a nonlinear optimization
Gibbs energy minimization techniques.11 Also, they
based on homotopy continuation and Lagrange
canbeclassedaseitherstoichiometricornonstoichio-
function while Burgos-Solórzano et al.10 have used
metric, depending on the way in which the elemen-
a reliable stability analysis, based on interval math-
tal abundance constraints are used in the algo-
ematics, to validate the results of phase and chemi-
rithm.13
cal equilibrium calculations.
Strong interactions among components, phases
and reactions may cause that this thermodynamic The stochastic optimization methods such as
problem exhibits highly nonlinear behavior. By Simulated Annealing (SA), Tabu Search (TS) or
Genetic Algorithm (GA) have not been extensively
studied in thermodynamic calculations for reacting
*Corresponding author:
mixtures. These methods have a great potential in
tel:[phone] ext.127, e-mail:[email protected]
286 A.BONILLA-PETRICIOLETetal.,Gibbs EnergyMinimizationUsing Simulated…,Chem. Biochem.Eng. Q.22 (3)285–298 (2008)
this area because they are reliable numerical tools namic properties of non reacting mixtures can be
thatcanbeappliedformultivariableandnonconvex extended to reacting systems and, by consequence,
problemsand,by consequence, they aresuitablefor the non-reactive phase equilibrium algorithms can
performing phase equilibrium calculations in reac- be easily modified to account for the equilibrium
tive systems. For example, Reynolds et al.14 have reactions.9 Also, reactive phase diagrams in trans-
outlined a general procedure for the global mini- formed coordinates look very similar to the non-re-
mization of Gibbs energy in non-reactive and reac- active ones in the standard mole fraction coordi-
tive mixtures using SA method and a nonstoichio- nates. The main difference is in the shape of the
metric approach. Unfortunately, these authors have composition space, which depends on the type of
reported results for only non-reactive systems. reaction in the system and could change if different
On the other hand, we consider convenient to reference components are selected.9 Transformed
remarkthat mostof published methods use reaction amount fractions X are defined as
i
extentsandmolefractionsasindependentvariables.
x (cid:2)v N(cid:2)1x
For example,in this category the methodsproposed i i ref
X (cid:1) i(cid:1)1,(cid:3),c(cid:2)R (1)
by McDonald and Floudas2 and Jalali and Seader7 i 1(cid:2)v N(cid:2)1x
TOT ref
areclassed.Asuitablechoicetoreducetheproblem
dimensionality, and to favor the numerical perfor- wherecisthenumberofcomponents,Risthenum-
mance of solution algorithms, consists in using ber of independent reactions, x is the column
ref
techniques of variable transformation.4,15 Only a vector of R reference component mole fractions, v
i
few methods have used these transformed variables is the row vector of stoichiometric number of com-
with the aim of improving the numerical behavior ponent i for each reaction, v is a row vector
TOT
(efficiency andreliability) ofsolving strategies.3,4,8,9 where each element corresponds to the sum of the
Thesekindsofalgorithmsarevery attractiveforthe stoichiometric coefficients for all components that
simulationofseparationprocessandfavorthestudy participate in reaction r, and N is a square matrix
of complex multireactive multicomponent systems. formed from the stoichiometric coefficients of the
Moreover, these approaches, in combination with referencecomponentsintheRreactions.3,15Itisim-
stochastic optimization techniques, can be used to portant to note that eq. (1) provides constant values
develop reliable strategies for phase equilibrium of tranformed mole fractions in single-phase reac-
calculations in reactive systems. tions, or constant overall mole fractions in hetero-
The objective of this work is to introduce an geneous reactions:
alternative approach for performing two-phase
z (cid:2)v N(cid:2)1z
equilibrium calculations in reactive systems. We use i i ref
Z (cid:1) i(cid:1)1,(cid:3),c(cid:2)R.
the transformed variables of Ung and Doherty3,15,16 i 1(cid:2)v N(cid:2)1z
TOT ref
and the stochastic optimization method Simulated
Annealing for the global minimization of Gibbs en-
In multiple-phase reactions, X -values in coexistent
ergy in reacting mixtures with two-phase equilib- i
phases are variable and subject to phase equilibria
rium. Numerical performance of this strategy is
requirements.
tested using several multicomponent reactive mix-
Thereferencemolefractionsx arecalculated
tures and our results show that it is generally a suit- ref
using eq. (1) and the equilibrium constants for each
ablemethodtoperformthiskind ofthermodynamic
reaction Kr by solving a system of R nonlinear
calculations. eq
equations given by
c
Thermodynamic problem statement Kr (cid:1)(cid:4)avir r(cid:1)1,(cid:3),R (2)
eq i
i(cid:1)1
Ung and Doherty3,15,16 proposed the use of
transformed composition variables with the objec- wherea istheactivityofcomponentiandvr isthe
i i
tive of developing a simpler thermodynamic frame- stoichiometric number of component i in the reac-
work for modeling reactive systems. These trans- tion r, respectively. When we know the R reference
formed variables depend only on the initial compo- mole fractions, for a set of c(cid:2)R specified trans-
sition of each independent chemicalspecies and are formed variables X , the corresponding mole frac-
i
constant as the reactions proceed. They also restrict tions of c(cid:2)R non-reference components are calcu-
the solution space to the compositions that satisfy lated using eq. (1). In this work, we used the bisec-
stoichiometry requirements and reduce the dimen- tion method to find the mole fraction of reference
sion of the composition space by the number of in- component in all examples with only one reaction.
dependent reactions. These characteristics allow Our experience indicates that bisection method
that all of the procedures used to obtain thermody- works very well for this purpose.17 On the other
A.BONILLA-PETRICIOLETetal.,Gibbs EnergyMinimizationUsing Simulated…,Chem. Biochem.Eng. Q.22 (3)285–298 (2008) 287
hand, for multireactive mixtures, the Newton (cid:16)
(cid:12)
method can be used to find the mole fractions of F j (cid:1)1 (6)
reference components.3,15 j(cid:1)(cid:17)
Classical thermodynamics indicates that mini-
c(cid:2)R
mizationofGibbsenergyisanaturalcourseforcal- (cid:12)X j (cid:1)1 j(cid:1)(cid:17),(cid:3),(cid:16) (7)
i
culating the equilibrium state of a system. Gibbs
i(cid:1)1
energy minimization algorithm was introduced by
White et al.18 and, for reactive mixtures, minimiz- We note thatF j is defined as
ing Gibbs energy is equivalent to minimizing the
(1(cid:2)v N(cid:2)1x j )
transformed Gibbs energy.19 In a multicomponent F j (cid:1)(cid:18) j TOT ref j(cid:1)(cid:17),(cid:3),(cid:16) (8)
multireaction mixture, with c components and R in- (1(cid:2)v N(cid:2)1z )
TOT ref
dependent chemical reactions, the dimensionless
transformed molar Gibbs free energy of mixing can being(cid:18) j themolefractionofphasejwhosefeasible
be defined as9
domain is (0, 1) and, in accordance with,16 is subject
(cid:16)
(cid:5)g(cid:1) g(cid:1)(cid:2)g(cid:1) c(cid:12)(cid:2)R (cid:6)(cid:8)(cid:5)(cid:13)(cid:1) {X}(cid:9)(cid:11) to (cid:12)(cid:18) j (cid:1)1(cid:20)v (cid:19) where (cid:19) is the column vector
(cid:1) 0 (cid:1) X (cid:8) i (cid:11)(cid:1) TOT
R T R T i(cid:7) R T (cid:10) j(cid:1)(cid:17)
g g i(cid:1)1 g
(3) ofRextentsofreactionand z isthecolumnvector
ref
c(cid:2)R c(cid:2)R (cid:6)x (cid:15)(cid:1) {X}(cid:9) of R reference component mole fractions associated
(cid:1)(cid:12)X ln(x (cid:14) {X})(cid:1)(cid:12)X ln(cid:8)(cid:8) i i (cid:11)(cid:11) to the transformed feed Z, respectively.
i i i i (cid:7) (cid:15) (cid:10)
i(cid:1)1 i(cid:1)1 i (cid:5)g(cid:1)
At equilibrium, must be at the global
R T
whereR isthegasuniversalconstant,Tisthetem- g
g minimum which is a necessary and sufficient condi-
(cid:5)(cid:13)(cid:1) (cid:13)(cid:1) {X}(cid:2)(cid:13)0
perature, i (cid:1) i i is the transformed tion for a thermodynamically stable state.3,9,15,19
R T R T However, the global minimization of transformed
g g
chemical potential of component i, g(cid:1) is the pure Gibbs energy is a challenging optimization problem
component free energy, g(cid:1) is the transf0ormed molar due to the objective function is nonconvex and may
Gibbs free energy, (cid:15) is the fugacity coefficient of havemultiplelocaloptimumsevenfortwo-phasere-
pure component, (cid:15)(cid:1) iis the fugacity coefficient of acting systems modeled with simple thermodynamic
componentiinthemiixtureand(cid:14) istheactivityco- equations. Considering this fact, local optimization
efficient of component i, respectiively. Note that (cid:15)(cid:1) methods are not suitable to solve this problem and a
and(cid:14) arefunctionsofthetransformedcompositioni robustoptimizationstrategymustbeused.Inthefol-
variabiles X. Also, calculations of pure component lowing section we describe the optimization proce-
free energies are avoided, which do not influence dure used to perform the global minimization of
the equilibrium and stability calculations, if instead transformed Gibbs energy for two-phase equilibrium
of g(cid:1) we use the transformed molar Gibbs energy of calculations in reactive systems.
(cid:5)g(cid:1)
mixing .9
R T
g Optimization approach
For a feed with a global transformed composi-
(cid:5)g(cid:1) In this paper we have extended the optimiza-
tion Z that splits in (cid:16) phases, is given by
tion strategy proposed by Rangaiah,20 which was
R T
g originally developed to non-reactive mixtures, to
perform an unconstrained Gibbs energy minimiza-
(cid:5)g(cid:1) (cid:1)(cid:12)(cid:16) F j(cid:12)c(cid:2)R X j(cid:6)(cid:8)(cid:5)(cid:13)(cid:1) ij(cid:9)(cid:11) (4) tion in two-phase reactive systems. Using this ap-
(cid:8) (cid:11)
R T i(cid:7)R T(cid:10) proach, the transformed Gibbs energy minimization
g j(cid:1)(cid:17) i(cid:1)1 g
problem can be simplified by using a set of new
The equilibrium transformed mole fractions variables instead of the transformed mole numbers
X j forallphasesmustsatisfythematerialbalance or fractions as decision variables in the optimiza-
i tion strategy. The introduction of these variables
eliminates the restrictions imposed by material bal-
(cid:16)
Z (cid:2)(cid:12)F jX j (cid:1)0 i(cid:1)1,(cid:3),c(cid:2)R (5) ances,reducesproblemdimensionalityandtheopti-
i i
mization problem is transformed to an uncon-
j(cid:1)(cid:17)
strained one. This approach is more suitable than
where F j is the transformed amount fraction for the Lagrange multiplier formulation due to the sig-
phase j. All transformed variables are subject to the nificant reduction of problem dimensionality. So,
following restrictions for two-phase reacting systems where all trans-
288 A.BONILLA-PETRICIOLETetal.,Gibbs EnergyMinimizationUsing Simulated…,Chem. Biochem.Eng. Q.22 (3)285–298 (2008)
formed mole fractions have values in the range point is accepted, the algorithm moves on from that
X (cid:21)(0,1) or its phase equilibrium region satisfies point. If the trial is rejected, another point is chosen
i
this restriction, real variables (cid:22) in the range [0, 1] for a trial evaluation. After having adjusted each el-
i
are defined and employed as optimization variables ement of VM periodically, half of all function eval-
using the next equations uationsinthatdirectionareacceptable.Thetemper-
ature reduction factor RT is used to decrease T
n(cid:1)(cid:17)(cid:1)(cid:22) n(cid:1)T i(cid:1)1,(cid:3),c(cid:2)R (9) employing T j(cid:20)1(cid:1)RT(cid:25)T j where j is the iteratioSAn
i i i SA SA
counter. Thus, asT declines, downhill moves are
SA
n(cid:1)(cid:23)(cid:1)n(cid:1)T (cid:2)n(cid:1)(cid:17) i(cid:1)1,(cid:3),c(cid:2)R (10) less likely to be accepted and the percentage of re-
i i i
jections rises. As T declines, VM falls and the
SA
wheren(cid:1)T (cid:1)n(cid:1)(cid:17)(cid:20)n(cid:1)(cid:23) isthetransformedmolenumber method focuses upon the most promising area for
i i i
of component i in the feed, n(cid:1) j (cid:1)n (cid:2)vTN(cid:2)1n is optimization. If the final function values from the
i i i ref
the transformed mole number of component i at last temperatures (NEPS = 4) differ from the corre-
phase j and (cid:22) is the optimization variable of the sponding value at the current temperature by less
i
unconstrained optimization problem, respectively. than a tolerance value (EPS) and the final function
Thebenefitofthismodificationisthatalltrialcom- value at the current temperature differs from the
positions will satisfy the material balance which al- current optimal function value by less than this tol-
lowsthe easy application of optimization strategies. erance value, algorithm execution terminates.
It is important to note that this formulation can be Corana et al.30 provide a full description of this al-
used if n(cid:1)T (cid:24)0 for all i(cid:1)1,(cid:3)c(cid:2)R. gorithm and we used the subroutine implemented
i
by Goffe et al.32 Fig. 1 shows the flow diagram of
Transformed Gibbs energy function is mini-
this algorithm.
mized using the simulated annealing (SA) stochas-
tic optimization method. SA is a robust numerical
tool that presents a reasonable computational effort
in the optimization of multivariable functions; it is
applicable to ill-structure or unknown structure
problems, it requires only calculations of the objec-
tive function and can be used with all thermody-
namic models.21,22 In fact, SAhas the attributes of a
good numerical optimization method if is properly
implemented: generality, acceptable computational
time,reliability and easeof use.23–25 Itisconsidered
as global optimization strategy and is one of the
most used stochastic methods in engineering appli-
cations.26 Specifically, in the field of thermodynam-
ics,thismethodhasbeensuccessfullyusedinphase
stability and equilibrium calculations14,17,20–22,27 and
nonlinear parameter estimation.28,29 This work in-
troduces the application of this optimization strat-
egy in the minimization of transformed Gibbs en-
ergy in reactive mixtures.
SA is a generalization of a Monte Carlo
method and its concept is based on the thermody-
namic process of cooling of molten metals to
achievethelowestenergystate.Generally,itcanlo-
cate the global optimum independently of initial
guesses if the values for its algorithm parameters
are properly selected. We have used the algorithm
proposed by Corana et al.30 because of its proved
reliability in thermodynamic calculations.17,20–22,29
In this algorithm, a trial point, randomly chosen
within the step length VM (a vector of length n op-
timization variables) is the starting point. The func-
tion evaluated at this trial point is compared to its
value atthe initial point. The Metropolis criterion,31
with aparametertermedannealing temperatureT ,
SA
is used to accept or reject the trial point. If the trial Fig. 1 – Flow diagram of SA optimization algorithm29
A.BONILLA-PETRICIOLETetal.,Gibbs EnergyMinimizationUsing Simulated…,Chem. Biochem.Eng. Q.22 (3)285–298 (2008) 289
Numerical performance (reliability and effi- used in severalapplications of other stochasticopti-
ciency) of SA method is significantly affected by mization methods.33–35 We must remark that the
the cooling schedule which is linked to the parame- global optimumOBJ for all exampleswasdeter-
min
ters T0 ,RT and NT. These parameters require minedby solving theequality oftransformedchem-
SA
pre-calibration for new problems. We have tuned ical potentials using several initial estimations. The
these parameters by performing several flash cal- stable two-phase solutions obtained by this proce-
culations for the reactive mixture: acetic acid dure were considered as the corresponding global
+ n-butanol (cid:1) water + n-butyl acetate at 25 °C. optimum OBJ values.
min
Thermodynamic properties of this mixture were On the other hand, the numerical performance
modeled with UNIQUAC equation using the model of our optimization strategy is characterized using
parameters reported by Wasylkiewicz and Ung.9 widespread criterions reported in the literature of
Based on these calculations (not reported in this stochastic optimization methods;25,33,36,37 a) success
paper) and our numerical experience in phase sta- rate (SR) of finding the global minimum given as
bility analysis for several reactive mixtures,17 we percent of calculations performed that satisfies eq.
have defined that SA works acceptably well using (11), b) mean total number of function evaluations
thefollowing conditions:T0 =10,RT=0.85, NT= (NFEV)during the optimization procedure, c)mean
SA
5·(c–R) where NT is the iteration number before computational time and d) mean absolute percent-
annealing temperature reduction. age deviation of the calculated compositions from
the known compositions at global minimum
Results and discussion 100 (cid:12)2 c(cid:12)(cid:2)R Ximjin (cid:2)XCijalc
AAD(cid:1) (12)
2(c(cid:2)R) X min
We have tested the numerical performance of i(cid:1)1 j(cid:1)1 ij
proposed method using several examples with dif-
where X min is the global optimum value for the
ferent dimensionality and thermodynamic models.
ij
Most of our examples are standard benchmark in transformed composition of component i at phase j
the literature of reactive distillation process, predic- and XCalc is the calculated value for the trans-
ij
tion of reactive homogeneous and heterogeneous
formedcompositionofcomponentiatphasejusing
azeotropes, multiple steady states in reactive sepa-
the stochastic method, respectively. Criteria a) and
ration units and phase stability analysis. All exam-
d) are used to characterize the reliability of optimi-
ples were solved several times using different feed
zation strategy while remaining ones are indicators
conditions (temperature,pressure,chemicalequilib-
of algorithm efficiency.25 Mean values of NFEV,
rium constants or feed compositions). We have as-
AAD and computational time are calculated consid-
sumed that all chemicalreactions are reversible and
ering only the successful calculations. All calcula-
they occur in both phases. For variable transforma-
tions were performed on a Processor Intel Pentium
tion,inallexampleswehaveselectedarbitrarilythe
M 1.73 GHz with 504 MB of RAM using Fortran
reference components. Also, tested conditions were
4.0 software.
arbitrary but we consider that they are sufficient to
demonstrate the numerical performance of pro- Phase stability analysis of calculated equilib-
posed strategy. rium compositions was performed using the Reac-
To determine the computational behavior of tive Tangent Plane Distance Function (RTPDF)
our algorithm, 100 runs are performed for each ex- which is given by9,19
ampleusing random initial values for decision vari-
c(cid:2)R
ables (cid:22) andrandomnumberseedsforSA.Wehave (cid:12)
i RTPDF(cid:1) X ((cid:13)(cid:1) (cid:2)(cid:13)(cid:1) ) (13)
defined a tolerance value of ESP = 1·10–6 for the i i X i Z
i(cid:1)1
convergence of SA method. Since stochastic opti-
mization methods do not provide an accurate solu-
where the global minimum of RTPDF is (cid:27) 0 for
tion of the global optimum,25,33 we have considered
stable mixtures. This function was globally mi-
that the global minimization of transformed Gibbs
nimized also using SA optimization method.
energy is successful upon satisfying the condition
Bonilla-Petriciolet et al.17 have reported that SA is
robust to perform phase stability analysis in react-
|OBJ (cid:2)OBJ |(cid:26)|OBJ |(cid:25)10(cid:2)4(cid:20)10(cid:2)5 (11)
calc min min ing and non-reacting mixtures. All reported solu-
tions in this work are thermodynamically stable.
where OBJ is the global minimum of the trans-
min
formed Gibbs energy and OBJ is the calculated Finally, the overall method is outlined as fol-
calc
value for this thermodynamic function with the op- lows for aflash calculation in amulticomponent re-
timization method, respectively. Eq. (11) has been active system:
290 A.BONILLA-PETRICIOLETetal.,Gibbs EnergyMinimizationUsing Simulated…,Chem. Biochem.Eng. Q.22 (3)285–298 (2008)
1. Input n0, equilibrium constants and model where X (cid:21)(0,1).Margules solution modelwith the
i i
parameters. data reported by Ung and Doherty15is used to cal-
culate thermodynamic properties
n0
2. Calculate z by z (cid:1) i and Z using eq.
i i (cid:12)c n0 i Gex (cid:1)1(cid:12)(cid:12)A x x (16)
i R 2 ij i j
i(cid:1)1 g i j
(1). Do a phase stability analysis using Z and eq.
with A (cid:1)478.6,A (cid:1)1074.484and A (cid:1)626.9,re-
(13). Continue with step 3 if the global optimum of 12 13 23
spectively.Wehavestudied thissystemat323.15 K
RTPDF < 0, otherwise propose other n0 and repeat
i and phase equilibrium calculations are performed
step 2.
for different values of the reaction equilibrium con-
(cid:12)c stantK (cid:21)(2.25,30) and feedcompositions. Based
3. Calculate n(cid:1)T (cid:1) n0(cid:2)v N(cid:2)1n0 . eq
i TOT ref on our formulation, we have only two optimization
i(cid:1)1 variables for this example. Table 1 shows the re-
4. Guess (cid:22) (cid:21)(0,1) for i(cid:1)1,(cid:3),c(cid:2)R. sults of equilibrium calculations for this reacting
i
5.Calculaten(cid:1)(cid:17) andn(cid:1)(cid:23) usingeq.(9)and(10). mixture. For all calculations performed, the pro-
i i posed optimization strategy is very reliable to find
n(cid:1) j n(cid:1) j the equilibrium compositions corresponding to
6. Calculate X j (cid:1) i and F j (cid:1) where
i n(cid:1) j n(cid:1)T the global minimum of transformed Gibbs energy.
With respect to efficiency of SA method, mean
c(cid:2)R
n(cid:1)T (cid:1)n(cid:1)(cid:17)(cid:20)n(cid:1)(cid:23) and n(cid:1) j (cid:1)(cid:12)n(cid:1) j. NFEV ranged from 40 068 to 41 489 where this
i
numerical effort is equivalent to 5 s of compu-
i(cid:1)1
7. Define x j as function of X j and x j . tational time. For this mixture, AAD is around
i i ref 6.1·10–03 %.
8. Introduce x j in the equations of equilibrium
i
constants for each reaction Kr .
eq
9. Determine x j by solving the system of R Table 1 – Numerical performance of simulated annealing
ref in the unconstrained global minimization of Transformed
nonlinear equations formed with Kr . This step is
performed for both phases to determeqine their com- Gibbsenergyforthereactingmixture A1(cid:20)A2(cid:1)A3at323.15K
(Margules solution model)
position in conventional mole fractions x.
Feed Numerical
10. Calculate transformed Gibbs energy using Equilibrium conditions
conditions performance1
eq. (4).
transformed
11. Check for minimum; proceed with SA
Gibbs
method with new (cid:22) (cid:21)(0,1) and repeat steps 5–11
energy
i Z K X(cid:17) X(cid:23) NFEV SR, %
until satisfy the convergence criterion of SA 1 eq 1 1 (cid:5)g(cid:1)
method.
RT
g
This algorithm has been applied for all calcula-
tions performed in this paper. In the following text, 0.408 2.25 0.3921 0.4239 –0.4918 40068±1584 100
wedescribetheresultsobtained for severalreacting
0.5 2.5 0.4863 0.7400 –0.5237 40785±1458 100
systems using this strategy.
0.7 3 0.4936 0.8069 –0.4267 40582±1437 100
Example 1. Our first example is a hypothetical
0.55 4 0.5000 0.8606 –0.6761 40977±1381 100
reacting ternary mixture A (cid:20)A (cid:1) A with liquid –
1 2 3
liquid equilibrium. This system was introduced by 0.776 5 0.5030 0.8851 –0.4343 40748±1269 100
UngandDoherty15intheirseriesofpublicationsre-
0.63297.5 0.5065 0.9124 –0.7844 41100±1375 100
lated to their transformed variables for reacting
mixtures. Transformed mole fractions, considering 0.75 10 0.5080 0.9243 –0.6246 41073±1225 100
third component as reference component, are given
0.887 15 0.5094 0.9354 –0.3576 41105±1169 100
by
0.535 20 0.5101 0.9406 –1.4197 41489±1238 100
x (cid:20)x
X (cid:1) 1 3 (14) 0.706 30 0.5107 0.9457 –1.0340 40598±1471 100
1 1(cid:20)x
3 1NFEVisthemeantotalnumberof function evaluations involved in
theglobalminimizationoftransformedGibbsenergyandSRisthesuc-
x (cid:20)x
X (cid:1) 2 3 (cid:1)1(cid:2)X (15) cessratetofindthetransformedphaseequilibriumcompositions.100
2 1(cid:20)x 1 trials performedwith randominitial values for optimization variables
3 and randomnumberseeds.
A.BONILLA-PETRICIOLETetal.,Gibbs EnergyMinimizationUsing Simulated…,Chem. Biochem.Eng. Q.22 (3)285–298 (2008) 291
Example 2. Butyl acetate is widely used as sol-
vent in coating, adhesives and paint industries and it
can be produced via the liquid-phase reaction of
butanol and acetic acid in the presence of a suitable
acidiccatalyst.11Thisreactionisgivenbyaceticacid
(1) + n-butanol (2) (cid:1) water (3) + n-butyl acetate (4).
Khaledi and Bishnoi11 have modeled the production
of butyl acetate via reactive distillation process. In
this work, we have calculated several tie-lines for
this mixture at 25°C. UNIQUAC model is used
to predict thermodynamic properties with the pa-
rameters reported by Wasylkiewicz and Ung.9 The
chemical equilibrium constant is calculated using
450 Fig. 2 – Calculated tie-lines in transformed mole fractions
lnK (cid:1) (cid:20)0.8 where T is given in K. Trans- for the reaction (1)acetic acid + (2)n-butanol (cid:1) (3)water +
eq T (4) n-butyl acetate at 25 °C. UNIQUAC solution model.
formed compositions are defined using n-butyl ace-
tate as reference component and are given by
transformed Gibbs energy. Specifically, an incre-
X (cid:1)x (cid:20)x (17) ment onT0 and NT provides a better performance,
1 1 4 SA
of course, at the cost of a greater computational ef-
X (cid:1)x (cid:20)x (18)
2 2 4 fort. Also, we can improve the performance of SA
X (cid:1)x (cid:2)x (cid:1)1(cid:2)X (cid:2)X (19) method in these difficult problems using a proper
3 3 4 1 2 initialization strategy. For example, we can use the
where X ,X (cid:21)(0,1) and X (cid:21)((cid:2)1,1). Trans- results of phase stability as initial values for Gibbs
1 2 3
formed Gibbs energy is minimized considering energy minimization.4 Finally, with respect to effi-
threeoptimizationvariables.Calculatedtie-linesfor
this system are shown in Fig. 2 and the numerical
performanceof optimization strategy appears in Ta- Table 2 – Numericalperformanceofsimulatedannealingin
the unconstrained global minimization of transformed Gibbs
ble 2. In this table, we also report the slopes of
energyforthereactingmixture(1)aceticacid+(2)n-butanol(cid:1)
tie-lines X(cid:28) which are calculated using
12 (3)water+(4)n-butylacetateat25°C(UNIQUACsolutionmodel)
X (cid:17)(cid:2)X (cid:23) Equilibrium Numerical
j j
X(cid:28) (cid:1) j(cid:1)2,(cid:3),c(cid:2)R (20) conditions performance1
1j X (cid:17)(cid:2)X (cid:23)
1 1 transformed
Gibbs energy
c(cid:2)R
where(cid:12)X(cid:28) (cid:1)(cid:2)1.Ourresultsshowthatproposed Z X1(cid:28)2 (cid:5)g(cid:1) NFEV SR, %
1j
RT
j(cid:1)2 g
algorithm is reliable to find the equilibrium compo-
(0.01, 0.4, 0.59) 46.0948 –0.2217 93757±2543 100
sitions and it generally exhibits a 100 % success
rate. Only for two feeds, it shows a poor perfor- (0.1, 0.2, 0.7) 2.6796 –0.4533 94477±2839 100
mance. In these cases, SA method converged to
(0.15, 0.5, 0.35) 3.7574 –0.8034 94357±2744 100
trivial solutions (localoptimums)whereitsvalueof
transformed Gibbs energy is very near to the global (0.2, 0.3, 0.5) 1.8920 –0.7867 94537±2491 100
one. For example, at Z(0.394, 0.274, 0.332), the
(0.3, 0.3, 0.4) 1.3410 –0.9934 94273±2655 100
trivial solution (Z (cid:1)X (cid:17)(cid:1)X (cid:23)) has a transformed
i i i (0.3, 0.4, 0.3) 1.6227 –1.1063 93985±2454 100
Gibbs energy value of –1.115 167 while the global
optimum is equal to –1.115 171. We note that these (0.397, 0.294,
1.0649 –1.1496 91297±3413 0
0.309)
local optimums satisfy the condition given by eq.
(11). However, the phase stability analysis helped (0.394, 0.274,
1.0323 –1.1152 91789±3620 3
us to identify these failures of the proposed optimi- 0.332)
zation method. As indicated by Burgos-Solozarno
(0.3, 0.15, 0.55) 0.9692 –0.7847 94237±2435 100
et al.,10 phase stability analysis is a fundamental
procedure which must be used to validate the re- (0.27, 0.1, 0.63) 0.9176 –0.6617 93961±2244 100
sults of any phase equilibrium calculation. On the 1NFEVisthemeantotalnumberof function evaluations involved in
other hand, for these difficult cases,it is convenient theglobalminimizationoftransformedGibbsenergyandSRisthesuc-
to modify the cooling schedule of SA algorithm to cessratetofindthetransformedphaseequilibriumcompositions.100
trials performedwith randominitial values for optimization variables
favor the convergence to the global optimum of
and randomnumberseeds.
292 A.BONILLA-PETRICIOLETetal.,Gibbs EnergyMinimizationUsing Simulated…,Chem. Biochem.Eng. Q.22 (3)285–298 (2008)
ciency, the computational effort in terms of mean
function evaluations ranged from 91 297 to 94 537
while mean computation time is around 40 seconds
and AAD is lower than 4.5·10–02 % for all tested
conditions.
Example 3. Third example is the reaction of
isobutene (1), methanol (2) and methyl tert-butyl
ether (3) with n-butane (4) as an inert. MTBE is an
important industrial chemical and has been used an
anti-knock agent to replace tetra-ethyl lead in gaso-
line. Several simulation and experimental researches
have been performed to study the MTBE reactive
Fig. 3 – Calculated tie-lines in transformed mole fractions
distillation process.38 For example, Okansinski and
for the reaction (1) isobutene + (2) methanol (cid:1) (3) methyl
Doherty39 have studied the effect of the reaction tert-butyl ether and(4)n-butane asinert at 100 °Cand10.13
equilibrium constant on the existence and location bar. Wilson model and Antoine equation.
of reactive homogeneous azeotropes in this mix-
ture. In this work, we have considered the vapor –
lished methods cannot handle the presence of inert
liquid equilibrium for this reaction at p = 10.1325
components.3 The inert components do not partici-
barand (cid:29)=100 °C.Transformedmolefractions are
pate in any of the reactions, but have an influence
defined using MTBE as reference component and
on the phase equilibrium. Based on these results, it
are given by
appears that proposed strategy is robust even when
x (cid:20)x inert components are considered.
X (cid:1) 1 3 (21)
1 1(cid:20)x
3
Table 3 – Numerical performance of simulated annealing
x (cid:20)x
X (cid:1) 2 3 (22) intheunconstrainedglobalminimizationoftransformedGibbs
2 1(cid:20)x energyforthereactingmixture(1)isobutene+(2)methanol(cid:1)
3
(3)methyl tert-butyl etherand(4)n-butaneasinertat100°C
x and 10.13 bar. (Wilson solution model and ideal gas).
X (cid:1) 4 (cid:1)1(cid:2)X (cid:2)X (23)
4 1(cid:20)x 1 2 Equilibrium Numerical
3 conditions performance1
where all transformed mole fractions ranged from 0
transformed
to 1. Wilson solution model and Antoine equation Gibbs energy
are used to calculated thermodynamic properties. Z X1(cid:28)2 (cid:5)g(cid:1) NFEV SR, %
Model parameters are taken from Maier et al.40 and RT
g
thereactionequilibriumconstantiscalculatedusing
(cid:5)G0 /R (cid:1)(cid:2)4205.05(cid:20)10.0982T(cid:2)0.2667T lnT (0.278, 0.365, 1.2106 –1.3846 94285±2756 100
rxs g 0.357)
whereTisgiveninK.Foridealgasbehavior,trans-
formed Gibbs energy is given by41 (0.3, 0.3, 0.4) 1.2028 –1.4343 94381±2415 100
(cid:5)g(cid:1) c(cid:2)R (cid:6)x p(cid:9) (0.35, 0.25, 0.4) 1.9706 –1.4894 94177±2658 100
R T (cid:30)(cid:12)Xi ln(cid:7)(cid:8)(cid:8)pisat(cid:10)(cid:11)(cid:11) (24) (0.4, 0.25, 0.35) 3.4749 –1.5332 93949±2966 100
g i(cid:1)1 i
(0.5, 0.3, 0.2) –5.7654 –1.5961 94729±2721 100
where psat is the vapor pressure of pure component
i (0.7, 0.25, 0.05) –1.1764 –1.3969 94177±2620 100
i.Forthisreactingmixture,wealsohavethreeopti-
mizationvariables.Calculatedtie-linesforthismix- (0.6713, 0.316,
–1.0483 –1.3990 94201±2499 100
0.0127)
ture appear in Fig. 3 and details of equilibrium cal-
culations are reported in Table 3. For all calcula- (0.15, 0.6, 0.25) 10.7574 –0.9409 94429±2519 100
tions performed, SA method shows a 100 % reli-
(0.00251,0.42965,
ability where phase equilibrium compositions are –11621.85 –0.4909 94813±2360 100
0.56784)
located without numerical problems. In the other
hand, mean value of NFEV ranged from 93 949 to (0.6, 0.3, 0.1) –1.5955 –1.5417 94141±2888 100
94 813 for all calculations performed while mean 1NFEVisthemeantotalnumberof function evaluations involved in
computation time is around of 85 s. For this react- the minimizationof transformed Gibbs energy and SR is the success
ing mixture, ADD ranged from 6.2·10–03 to ratetofindthetransformedphaseequilibriumcompositions.100trials
performed with random initial values for optimization variables and
7.48·10–02%.Itisimportanttonotethatsomepub- randomnumberseeds.
A.BONILLA-PETRICIOLETetal.,Gibbs EnergyMinimizationUsing Simulated…,Chem. Biochem.Eng. Q.22 (3)285–298 (2008) 293
Example4.Thisexampleisalsoahypothetical
reacting ternary mixture A (cid:20)A (cid:1) A with liquid –
1 2 3
liquid equilibrium. We have considered a reaction
equilibrium constant independent of temperature
and thermodynamic properties are calculated using
the Margules solution model where
gE
(cid:1)3.6x x (cid:20)2.4x x (cid:20)2.3x x (25)
R T 1 2 1 3 2 3
g
Consideringthirdcomponentasreferencecom-
ponent, transformed compositions are given by eqs.
(14) and (15). Iglesias-Silva et al.41 used this sys-
tem to introduce the concept of equal area rule for
multiphase equilibrium in reacting mixtures based
on transformed variables. Based on their paper, this
mixture shows three-phase equilibrium at K (cid:1)
eq
0.980019 and Z (cid:21)(0.218336, 0.815122). We have
1
performedphaseequilibriumcalculationsfordiffer-
ent values of K and feed compositions that show
eq
a transformed Gibbs function with multiple local
(cid:5)g(cid:1)
minima.Withillustrativepurposes, -surfacefor
R T
g
some tested conditions appears in Fig. 4 where the
tangent planes for stable and unstable equilibriums
are indicated. Fig. 4 shows that this mixture can
present unstable two-phase equilibrium states
where a local optimization method can be easily
trapped. For example, at Z(0.6, 0.4) and K =
eq
0.9825 this system has two unstable phase equilib-
rium states: X (cid:17)(0.3677,0.6323), X (cid:23)(0.8104,0.1896)
Table 4 – Numerical performance of simulated annealing
intheunconstrainedglobalminimizationoftransformedGibbs
energy for the reacting mixture A (cid:20)A (cid:1) A (Margules solu-
1 2 3
gE
tion model (cid:1)3.6x x (cid:20)2.4x x (cid:20)2.3x x )
RT 1 2 1 3 2 3
g
Feed condi- Numerical
Equilibrium conditions
tions performance1
transformed
Gibbs energy Fig. 4– Transformed Gibbs energy surface of a hypotheti-
Z1 Keq X1(cid:17) X1(cid:23) (cid:5)g(cid:1) NFEV SR, % cal reacting mixture A (cid:20)A (cid:1) A (Margules solution model
1 2 3
RT gE
g (cid:1)3.6x x(cid:20)2.4x x(cid:20)2.3x x )a)K =0.95,b)K =0.9825
RT 1 2 1 3 2 3 eq eq
0.6 0.98250.4845 0.8150 –0.1445 37772±724 100 g
and c) K = 1.0
eq
0.35 0.985 0.2184 0.4824 –0.1495 39857±376 100
(cid:5)g(cid:1)
0.75 0.995 0.4891 0.8147 –0.1439 39996±659 100 with = –0.1425 and X (cid:17)(0.2197, 0.7803),
R T
0.5 0.975 0.2156 0.8167 –0.1448 41105±1202 100 g (cid:5)g(cid:1)
X (cid:23)(0.8143,0.1857) with = –0.1442. So, this
R T
0.63 0.95 0.2037 0.8240 –0.1389 41553±723 100
g
example is a good choice to test the reliability of
0.559 1.0 0.4907 0.8146 –0.1507 39905±1169 100 the proposed optimization strategy. Two optimiza-
1NFEVisthemeantotalnumberof function evaluations involved in tion variables are considered in the minimization of
the minimizationof transformed Gibbs energy and SR is the success transformed Gibbs energy and results of phase
ratetofindthetransformedphaseequilibriumcompositions.100trials equilibrium calculations are reported in Table 4.
performed with random initial values for optimization variables and
randomnumberseeds. Again, for all tested conditions, SAmethod is capa-
294 A.BONILLA-PETRICIOLETetal.,Gibbs EnergyMinimizationUsing Simulated…,Chem. Biochem.Eng. Q.22 (3)285–298 (2008)
ble of finding the global minimum of transformed
Gibbs energy without numerical problems. On the
other hand, NFEV ranged from 37 772 to 41 553
while computation time is around 2.5 s. ADD is
lower than 0.02 % for all cases.
Example 5. Tert-amyl methyl ether (TAME) is
an important chemical for gasoline and is produced
by liquid-phase synthesis from methanol and
iso-amylenes catalyzed by a sulfonic acid ion ex-
change resin.38 In this reaction, five components
take part: methanol, 2-methyl-1-butene, 2-methyl-
-2-butene, TAME and n-pentane as inert compo-
nent. In this work, we have considered the lumped
single reaction which can be written as: 2-methyl- Fig. 5 – Calculated tie-lines in transformed mole fractions
-1-butene (1) +2-methyl-2-butene (2) +2 methanol forthereaction(1)2-methyl-1-butene+(2)2-methyl-2-butene+
+(3)2methanol (cid:1)(4)2tert-amyl methyl ether at 335Kand
(3) (cid:1) 2 TAME (4) with n-pentane (5) as an inert
1.52 bar. Wilson model and ideal gas.
solvent.12,38 In first instance, VLE of this reaction
without n-pentane as inert is studied. Transformed
mole fractions, considering TAME as reference Transformed Gibbs energy is minimized con-
component, are given by sidering four optimization variables. Phase equilib-
rium calculations are performed at 335 K and 1.52
x (cid:20)0.5x
X (cid:1) 1 4 (26) bar where the details of calculations and numerical
1 1(cid:20)x
4
x (cid:20)0.5x
X (cid:1) 2 4 (27) Table 5 – Numerical performance of simulated annealing
2 1(cid:20)x
4 intheunconstrainedglobalminimizationoftransformedGibbs
energy for the reacting mixture (1) 2-methyl-1-butene + (2)
x (cid:20)x
X (cid:1) 3 4 (cid:1)1(cid:2)X (cid:2)X (28) 2-methyl-2-butene + (3) 2 methanol (cid:1) (4) 2 tert-amyl methyl
3 1(cid:20)x 1 2 ether at 335 K and 1.52 bar (Wilson model and ideal gas)
4
Equilibrium Numerical
All transformed fractions ranged from 0 to 1.
conditions performance1
Wilsonandidealgasmodelshavebeenusedtocalcu-
late thermodynamic properties of this mixture where transformed
thermodynamic parameters are taken from Chen et Gibbsenergy
al.38 Reactionequilibriumconstantiscalculatedusing Z X1(cid:28)2 (cid:5)g(cid:1) NFEV SR, %
K (cid:1)1.057(cid:25)10(cid:2)04 e4273.5/T where T is given in K. RT
eq g
Phase equilibrium calculations are performed for sev-
(0.3, 0.15, 0.55) –0.2072 –1.0868 94453±2354 100
eral feeds at 335 K and 1.52 bar. Details of calcula-
tions are reported in Table 5 and calculated tie-lines (0.32, 0.2, 0.48) –0.2800 –1.2179 94393±2554 100
reported in transformed mole fractions appear in Fig.
(0.354, 0.183,
5. Only for onefeed, theSA method shows asuccess –0.2856 –1.2264 93089±2881 63
0.463)
rate lower than 100 %. For remaining calculations,
this method is very reliable to find the global mini- (0.2, 0.07, 0.73) –0.0076 –0.7003 94345±2935 100
mumoftransformedGibbsenergy.Withrespecttoef-
(0.15, 0.02, 0.83) 0.0064 –0.3950 94429±2519 100
ficiency, NFEV ranged from 93 781 to 94 873 where
mean computational time is around 83 s and AAD is (0.27, 0.3, 0.43) 0.8089 –1.2913 93781±2205 100
lower than 1.5·10–02 % for all tested conditions.
(0.2, 0.35, 0.45) –3.6767 –1.2218 94573±2549 100
In our second scenario for this example, we
have considered the presence of n-pentane as inert. (0.1, 0.35, 0.5) –8.5301 –0.9623 94477±2335 100
Transformedmolefractionsaredefinedbyeqs.(26)
(0.05, 0.3, 0.65) –162.6184 –0.7089 94873±2418 100
– (27) and
(0.025, 0.3,
x (cid:20)x 327.5080 –0.5787 94117±2614 100
X (cid:1) 3 4 (29) 0.675)
3 1(cid:20)x
4 1NFEVisthemeantotalnumberof function evaluations involved in
theglobalminimizationoftransformedGibbsenergyandSRisthesuc-
x
X (cid:1) 5 (cid:1)1(cid:2)X (cid:2)X (cid:2)X (30) cessratetofindthetransformedphaseequilibriumcompositions.100
5 1(cid:20)x 1 2 3 trials performedwith randominitial values for optimization variables
4 and randomnumberseeds.