Optimization of reactive distillation processes with simulated annealing

20 Pages·1904·0.37 MB·English

ChemicalEngineeringScience55(2000)5059}5078
Optimization of reactive distillation processes
with simulated annealing
M. F. Cardoso, R. L. Salcedo*, S. Feyo de Azevedo, D. Barbosa
DepartamentodeEngenhariaQun&mica,FaculdadedeEngenhariadaUniversidadedoPorto,RuadosBragas,4099Porto,Portugal
Received6August1999;accepted4April2000
Abstract
A simulated annealing-based algorithm (MSIMPSA) suitable for the optimization of mixed integer non-linear programming
(MINLP)problemswasappliedtothesynthesisofanon-equilibriumreactivedistillationcolumn.Asimulationmodelbasedonan
extensionofconventionaldistillationisproposedforthesimulationstepoftheoptimizationproblem.Inthecaseofidealvapor}liquid
equilibrium,thesimulationresultsaresimilartothoseobtainedbyCiricandGu(1994,AIChEJournal,40(9),1479)usingtheGAMS
environment and to those obtained with the AspenPlus modular simulator. The optimization results are also similar to those
previouslyreportedandsimilartothoseusinganadaptiverandomsearchalgorithm(MSGA).Theoptimizationswerealsoperformed
withnon-idealvapor}liquidequilibrium,consideringeitherdistributedfeedandreactiontraysorsinglefeedandreactiontray.The
resultsshowthattheoptimizedobjectivefunctionvaluesareverysimilar,andmostlyindependentofthenumberoftraysandofthe
reactiondistribution.Itisshownthattheproposedsimulation/optimizationequation-orientedenvironmentsarecapableofproviding
optimized solutions which are close to the global optimum, and reveal its adequacy for the optimization of reactive distillation
problemsencounteredinchemicalengineeringpractice. ( 2000ElsevierScienceLtd.Allrightsreserved.
Keywords: Simulatedannealing;Reactivedistillation;Simulation;Optimization
1. Introduction studies in reactive distillation. Afterwards, rigorous
mathematicalmodels for computer simulation began to
Reactive distillation columns can be attractive when- bedeveloped.Theseareessentiallyextensionsofmethods
ever conversions are limited by unfavorable reaction available for conventional distillation (Suzuki, Yagi,
equilibrium and when selectivity can be increased by Kamatsu & Hirata, 1971; Hunek, Foldes & Sawinsky,
allowing simultaneous reaction and separation in the 1979; Murthy, 1984; Chang & Seader, 1988; Doherty
sameprocessingunit.Thesearethetwomainadvantages & Buzad, 1992).
of reactive distillation over more conventional alterna- Mostof the studies in reactivedistillation are simula-
tives (Barbosa & Doherty, 1988a,b; Okasinski tion studies, viz. the resulting composition pro"les are
&Doherty,1998),wherebythetotalannualizedproduc- determined after speci"cation of feed composition and
tion cost can be decreased by an order of magnitude quality, column pressure, re#ux or reboiler ratio, total
(Doherty & Buzad, 1992). number of stages, feed plate location and liquid-phase
Reactivedistillationprocesses havereceivedattention volumes on each stage. There are several simulation
sincethebeginningof thecentury,whenthe"rstpropo- algorithms available (Buzad & Doherty, 1995), some
sal for the use of a distillation column as a reactor was implemented in commercial simulators, such as Aspen-
considered.In the early1920s a numberof patentswere Plus (Venkataraman, Chan & Boston, 1990). Abufares
registeredforcarryingoutesteri"cationreactionsindis- and Douglas (1995) developed a dynamic mathematical
tillation columns (Backhaus, 1921a}c 1922a}e 1923a,b). modelforamethyltertiarybutylether(MTBE)catalytic
Mostof the work preceeding1970 refer to experimental distillation process, which was implemented in the dy-
namic simulator SpeedUp. Recently, Huss, Chen, Mal-
one and Doherty (1999) proposed a general framework
basedonconceptualdesignforsolvingbothequilibrium
*Correspondingauthor.Fax:[phone].
E-mailaddress:[email protected](R.L.Salcedo). and rate-based reactive distillation processes. Also,
0009-2509/00/$-seefrontmatter ( 2000ElsevierScienceLtd.Allrightsreserved.
PII: S0009-2509(00)00119-6
5060 M.F.Cardosoetal./ChemicalEngineeringScience55(2000)5059}5078
Schenk, Gani, Bogle and Pistikopoulos (1999) have de- MINLPproblemusingageneralizedBendersdecompo-
"ned a general modeling framework that allows both sition algorithm (Geo!rion, 1972).
rate-basedandequilibrium-basedmodels,steadystateor Buzad and Doherty (1995) developed a procedure
dynamic,to be formulated and solved.Jimoh, Arellano- which determines the optimum value for the total
Garcia, Bock and Wozny (1999) have experimentally amountofliquidholdupfora"xedproductcomposition.
validatedtheirsteady-stateanddynamicmodelforreact- Theprocedureisa"xedpointmethodwhichdividesthe
ive distillation using a pilot-scale column and investi- compositionspace into regions of feasibility, and is pre-
gated the transesteri"cation of methyl myristate with sentedasanalternativeto methodsdevelopedforreact-
isopropanol to methanol and isopropyl myristate. ive distillation columns operating close to reaction
Inthesimulationofadistillationcolumntheresulting equilibriumconditions.Theprocedurewasdevelopedfor
separation is evaluated for a speci"ed column, while in distillationprocesseswiththreecomponentswherereac-
a design problem the desired separationis speci"edand tions such as 2C HA#B occur.
the column parameters are then evaluated (Buzad Pekkanen(1995)presentsadesignalgorithmfora&lo-
& Doherty, 1995). While it has been relatively easy to cal optimization’approach for reactivedistillation. This
extend simulation methods developed for non-reactive approach is characterized as a &stage by stage speci"ca-
distillation, it has not been easy to extend non-reactive tion’, where the design procedure starts from both col-
distillationdesign techniquesto reactivedistillationcol- umn ends and design speci"cations are made at each
umns. stage as the calculation progresses. The column is thus
Thereareseveralreasonsforthedi$cultyindevelop- notoptimizedasawholesincenotalltheparametersare
ing methods for the design of reactive distillation pro- simultaneouslyvaried.
cesses, as referred by Ciric and Gu (1994). In OkasinskiandDoherty(1998)haveextendedthe"xed
conventional distillation design, the major design vari- point design method for systems with isomolar or non-
ablesarethenumberoftraysinthecolumn,thefeedtray isomolar liquid-phase reaction, non-ideal vapor}liquid
location and the re#ux ratio. It is usually reasonable to equilibrium and allowing for a distribution of liquid
assumeconstantmolar#owrateineachcolumnsection, holdups on the reactive stages.
andthateachfeedstreamisintroducedinasingletray.In Hussetal.(1999)giveemphasistoconceptualcolumn
opposition, in reactive distillation processes the molar designthroughgeometricmethods,andhaveshownthat
#owratescanonlybesetconstantifthefollowingcondi- equilibriumreactivedesignscanbethestartingpointfor
tionsarealsoveri"ed:heatofreactionnegligibleandthe rate-based designs.
total number of moles constant. The liquid holdup on Inthiswork,weproposeanintegratedapproachtothe
eachtraycannotbeignored,evenforsteady-stateopera- simulation/optimizationof reactivedistillationcolumns,
tion,sinceitisherethatreactionstakeplaceanditisthis byadoptingtheMINLPcolumnmodelofCiricandGu
volume that determines the extent of the reactions. (1994).Thesimulationmodelisbasedonanextensionof
Barbosa and Doherty (1988a,b), Ciric and Gu (1994), conventionaldistillationcolumnsandtheoptimizationis
Buzad and Doherty (1995), Pekkanen (1995), Okasinski performedbyasimulatedannealing-basedMINLPalgo-
andDoherty(1998)andHussetal.(1999)aresomeofthe rithm (the MSIMPSA algorithm of Cardoso, 1998).
researchers who developed methods for the design of Resultsarealsocomparedwiththose obtainedusingan
reactive distillation columns. adaptive random search algorithm (the MSGA algo-
Barbosa and Doherty (1988a,b) developed techniques rithm of Salcedo, 1992). The proposed method can be
based on methods applicable to the design of conven- appliedtosingleormultiplereactions,idealornon-ideal
tional distillation columns. However, these techniques vapor}liquid equilibrium and distributed or single-
were developed for systems near chemical equilibrium, staged reaction zones.
with constant molar #ow and which involve a single
reaction.Whenanyoneoftheseassumptionsisviolated,
the method is no longer valid. 2. Simulated annealing and the MSIMPSA algorithm
CiricandGu(1994)proposedaninterestingtechnique
forthedesignofreactivedistillationcolumns,byformu- Simulatedannealingisawell-establishedtechniquefor
latingasynthesisproblemwhichcanbeappliedtomore the optimization of combinatorial problems, i.e. for the
genericsituations,whenthereismorethanonechemical optimization of large-scale functions that may assume
reaction,orwhenreactiveequilibriumorconstantmolar several distinct discrete con"gurations (Metropolis,
#ow cannot be assumed. With their column model they Rosenbluth,Rosenbluth,Teller&Teller,1953;Kirkpat-
have built a mixed integer non-linear programming rick, Gellatt & Vechi, 1983). It may also be used for
(MINLP) problem. The solution of this problem yields continuous(NLP)ormixed-integer(MINLP)problems.
the optimum number of trays, the optimal feed tray For a review of the available algorithms, the reader
location(s)andcomposition(s),andthecompositionpro- shouldrefertotheavailableliterature(Cardoso,Salcedo
"les within the column. These authors solved the & Azevedo, 1994, 1996).
M.F.Cardosoetal./ChemicalEngineeringScience55(2000)5059}5078 5061
Brie#y, the Metropolis algorithm simulates the phys- NLP and combinatorial problems. When applied to
icalprocessofannealing,i.e.meltingasolidbyincreasing NLP problems, it performs as the SIMPSA algorithm
itstemperature,followedbyslowcoolingandcrystalliza- (Cardosoetal.,1996)andwhenappliedtocombinatorial
tion into a minimum free energy state. From an optim- problems it performs as the original Metropolis algo-
ization point of view, simulated annealing explores the rithm (Cardoso et al., 1994).
keyfeatureof thephysicalannealingprocessofgenerat-
ing transitions to higher energy states, applying to the
new states an acceptance/rejection probability criterion
3. Synthesis of a reactive distillation column for the
whichshouldnaturallybecomemoreandmorestringent
production of ethylene glycol
with the progress of the optimization. Thus, simulated
annealingmaybeviewedasarandomizationdevicethat
3.1. Problem dexnition
allowssomewrong-waymovementsduringthecourseof
the optimization, through an adaptive acceptance/rejec-
Thesynthesisproblemfordesigningareactivedistilla-
tion criterion (Schoen, 1991). It is recognized that these
tioncolumnis statedby Ciric andGu (1994) asfollows.
arepowerfultechniqueswhichmighthelpinthelocation
Given:
of near-optimum solutions, despite the drawback that
theymay not be rigorous(exceptin an asymptoticway) f a set of chemical species i"1,2,I,
and may be computationally expensive (Grossmann f a set of desired products and production rate,
& Daichendt, 1994). f a set of chemical reactions j"1,2,R, the
The SIMPSA solver (Cardoso et al., 1996) is stoichiometric coe$cients for all species, l , and rate
ij
a simulated annealing-based algorithm that was de- expressions r "f(x,„),
j i
velopedforsolvingconstrainedNLPproblems.Itisbuilt f theenthalpyofvaporizationandvapor}liquidequilib-
onaschemeproposedbyPressandTeukolsky(1991)for rium data,
unconstrained continuous optimization, and combines f the feed streams composition,
the simplex method of Nelder and Mead (1965) with f the cost parameters,
simulatedannealing.Theroleofthenon-linearsimplexis
to generate continuous system con"gurations, while the the goal is to determine the optimal number of trays,
role of simulated annealing is to allow for wrong-way holdup per tray, re#ux ratio, condenser and reboiler
movements, simultaneously providing (asymptotic) con- duties and feed tray locations. Fig. 1 shows a schematic
vergence to the global optimum. The algorithm showed diagram of the process.
good robustness, viz. good insensitivity to the starting Ethyleneglycol(C2H6O2) is producedfromthereac-
point,andreliabilityinattainingtheglobaloptimum,for tion of ethylene oxide (C2H4O) and water:
a number of di$cult NLP problems described in the
C H O#H OPC H O . (R1)
literature(Cardosoetal.,1996).Itmayalsobeusedwith 2 4 2 2 6 2
a faster non-equilibrium variant of simulated annealing
However, the ethylene glycol produced can further
(Cardoso, Salcedo & Azevedo, 1994).
react with ethylene oxide to produce the unwanted by-
The MSIMPSA algorithm (Cardoso et al., 1997) was
product diethylene glycol (C H O ):
4 10 3
developedasanextensionoftheSIMPSAsolvertodeal
with MINLP problems. The continuous non-linear sol- C H O#C H O PC H O . (R2)
2 4 2 6 2 4 10 3
ver is used in an inner loop to update the continuous
parameters,whiletheMetropolisalgorithmisusedinan Both reactions are highly exothermic and occur at
outer loop to update the complete set of decision vari- moderatetemperatures,allowingproductionviaareact-
ables.Thestochasticallygeneratingschemeproposedby ive distillation column.
Dolan, Cummings and Levan (1989) is employed each According to Ciric and Gu (1994) there are two
timeadi!erentsetofdiscretevariablesisneeded,where- reasons for producing ethylene glycol via reactive
bya"rstrandomnumberchoosesthediscretevariableto distillation.Firstly,the largedi!erenceinvolatilitiesbe-
bechangedandasecondonechoosesifthechangeisto tween ethylene oxide and ethylene glycol will lead to
beupwards ordownwardsby one unit.The MSIMPSA a rapid separation of these two components in the col-
algorithm is simple to use, since it neither requires the umn,improvingtheoverallselectivity.Secondly,partof
preliminary identi"cation of integer candidates for the theheatrequiredfortheseparationisobtainedfromthe
globaloptimumnortheidenti"cationoffeasibleorgood heat of reaction, which allows the reduction of energy
starting points. It was tested with several di$cult costs.
MINLP functions published in the literature, showing Data for this problem can be found in Appendix
agoodrobustnessin attainingthe globaloptimum.The A*Table6(reactiondata),Table7(vapor}liquidequi-
MSIMPSAalgorithmhasawideapplicability,sinceitis librium constants, ideal VLE), Table 8 (cost data) and
applicable not only to MINLP problems, but also to Table 9 (thermodynamic data, non-ideal VLE). The
5062 M.F.Cardosoetal./ChemicalEngineeringScience55(2000)5059}5078
maximum value of the liquid holdup, and 20 for the
maximumnumberoftrays.Althoughtheselimitsarenot
strictly necessary for the optimization step, they were
initially enforced for comparing our results with those
from Ciric and Gu (1994). These constraints are
< )F , k"1,2,N,
k .!9
‚ )F , k"1,2,N,
k .!9
C
+ F )F , k"1,2,N,
ik .!9
i/1
= )= , k"1,2,N,
k .!9
aswellasthemaximumnumberoftrays,N .Thefeed
.!9
tothedistillationcolumnineachtrayiscomposedofone
feed stream of ethylene oxide and one feed stream of
water. Obviously, some of these feeds may have a null
value. This problem can be considered as a fairly large
highly non-linear non-convex problem (function of
N ),sincethematerialbalancescontainbiandtrilinear
.!9
terms and the reaction terms, the vapor}liquid equilib-
rium evaluation and the objective function are highly
non-linear.
For the optimization problem under study, there are
(3N #1) decision variables. One is integer corre-
.!9
sponding to the maximum number of trays, N ,
.!9
and3N decisionvariablesarecontinuous:theboil-up
.!9
fraction, b, N feed #ow rates corresponding to one
.!9
of the components, ethylene oxide or water, N !1
.!9
feed #ow rates corresponding to the other component
Fig. 1. Schematic diagram of a reactive distillation column for the (all but the feed #ow rate into the "rst tray next to the
productionofethyleneglycol.
reboiler), and N liquid holdups. For a maximum
.!9
numberof20traysthereare341variablesand61degrees
offreedom,whereoneisinteger.Theintegervariablecan
berepresentedbysixbinaryvariablesaccordingtoabi-
enthalpyofvaporizationofthemixtureisassumedtobe
nary expansion (Cardoso et al., 1997), since with the
constant all over the column with the value of
binaryexpansionlargervariationsoftheoriginaldiscrete
40]103 J mol~1(Ciric,1995).Themathematicalformu-
variable are allowed at each step of the optimization
lation of this problem as formulated by Ciric and Gu
process.
(1994) is completely described in Appendix B. The con-
ForthesimulationofacolumnwithN trays,allthe
.!9
stant C used for evaluating the column diameter (Eq.
D 3N continuous variables are required; however, for
.!9
(B.24)) has the value 0.01331 (dimensionless). a column with N trays, where N(N , only 3N con-
.!9
The objective function is the minimization of the
tinuous variables are involved while the other
total annualized cost which is composed of two basic 3(N !N) are simply ignored.
.!9
terms: annual operating costs and the annualized
investment. The raw material, steam and cooling water
set the annual operating costs. The costs of the 3.2. Simulation algorithm
columnshell,trays,reboilerandcondensercorrespondto
the annualized investment. To evaluate the objective The simulation of reactive distillation processes
function, the operating conditions determined by the involves the simultaneous solution of material and en-
column simulation for each decision vector have to be ergy balances and stoichiometric relationships, and
known. this corresponds to the solution of a considerable large
set of non-linear equations. The calculation proced-
3.1.1. Constraints ures reported for solving sets of non-linear equations
Ciric (1995) has set as 1000 kmol h~1 the maximum can be divided broadly into four categories (Chang
valueforthevaporandliquid#owrates,14.16 m3forthe & Seader, 1988; Venkataraman et al., 1990; Biegler,
M.F.Cardosoetal./ChemicalEngineeringScience55(2000)5059}5078 5063
Grossmann & Westerberg, 1997; Lee & Dudukovic, Newton or quasi-Newton methods converge quickly
1998): from suitable starting guesses (Murthy, 1984; Chang
&Seader,1988).Whenthestartingpointis farfromthe
(i) methods using equation decomposition (or tear- solution,thesemethodscanconvergetoimpossiblephys-
ing/partitioning); icalconditionsormaynotconvergeatall.Bastos(1987)
(ii) relaxationtechniques; hasdevelopedanalgorithmforthesimulationofconven-
(iii) methodsincorporatingNewtonorquasi-Newtonal- tionaldistillationcolumns thatmakes use of a modi"ed
gorithms; Newton}Raphsonmethod. As a "rst step the algorithm
(iv) homotopy-continuationmethods. applies a relaxation technique to improve the starting
guess for the composition pro"le and #ow rates in the
Decomposition methods allow the identi"cation of column, which are needed by the Newton}Raphson
blockswhichneedthesimultaneoussolutionandblocks method.
which can be solved sequentially. For example, Ledet Homotopy-continuationmethodshavetheadvantage
andHimmelbleau(1970)proposeamethodthatidenti"- of forcing the desired solution by tracking a homotopy
esaminimumsetofrecyclevariables,whichmaythenbe curveregardlessofthechoiceoftheinitialestimates.Lee
usedwithdirectsubstitution.Thesemethodsarefastand andDudukovic1998have foundtheseto besuperiorto
e$cient in what concerns the use of computer storage Newton-basedmethodsforsolvingthenon-linearsystem
space.However,whendi!erencesinboilingpointvalues ofequationsarisinginreactivedistillation,despitealon-
between components are large, when kinetics are com- ger computational burden.
plex,orwhenliquidsolutionsarehighlynon-ideal,these In this work, we propose a procedure for the simula-
direct methods su!er from poor convergence character- tionofreactivedistillationcolumnsbasedonthemethod
istics. developed by Bastos (1987) for conventional distillation
Withrelaxationtechniques,theliquid-phasecomposi- columns (see Appendix C). Table 1 shows the variables’
tions are computed based on non-steady-state material identi"cation,bothcalculatedandspeci"edbythesimu-
balances (Bastos, 1987), which in subsequent iterations lation algorithm.
proceed towards the steady-state solution. Relaxation
techniques are reliable but can be slow specially as the 3.2.1. General structure of the simulation problem
solution is approached (Bastos, 1987; Chang & Seader, For the reactive distillation column speci"ed as
1988; Venkataraman et al., 1990). above, we propose a simulation model based on the
Table1
Variablescalculatedandspeci"ed*simulationalgorithm
Variable Number Speci"ed Calculated
Numberoftrays 1 1
Compositions
Liquid C(N#1) C(N#1)
Distillate C C
Flowrates
Liquid N#1 N#1
Vapor N#1 N#1
Distillate 1 1
Bottomcomponents C 1 C!1
Feeds NNF NNF!1 1
Columntemperatures N P N
Boil-upfraction 1 1
Liquidholdups N N
Duties
Condenser 1 1
Reboiler 1 1
Extentofreactions NR NR
Vapor}liquidequilibrium NC NC
constants
Columndimensions
Diameter 1 1
Height 1 1
Total N(NF#2C#R#4)#3C#9 N(NF#1)#3 N(2C#R#3)#3C#6
5064 M.F.Cardosoetal./ChemicalEngineeringScience55(2000)5059}5078
following steps: The simulation is aborted when convergence on the
torn variable is obtained but not on the global material
(i) An initial estimate of the molar feed #ow rate not
balance,orwhen theNewton}Raphsonmethodexceeds
speci"ed by the decision vector, viz. a torn (recycle)
a pre-speci"ed maximum number of iterations. When
variable,isspeci"ed.WehavechosenF asthetorn
21 oneofthesesituationsisveri"ed,thecolumnisinfeasible
variable, although other recycle variables may
andtheobjectivefunctionvalueissimplypenalizedwith
equally work. This torn variable was found to be
a large constant value. This simple penalizingscheme is
appropriatefromaninformation#owpointofview,
possibleduetothestochasticnatureoftheoptimization
andneed not be a true recycle variable (from a ma-
algorithms employed in the present work.
terial #ow point of view).
The simulation algorithm was written in Fortran 77,
(ii) An initial estimate of the composition pro"le in the
and all runs were performed with double precision on
column is next obtained.
aHP730Workstation,runningcompileroptimizedFor-
(iii) This estimate is improved by a relaxation method.
tran77 code.The averagetimespentin eachsimulation
(iv) Thematerialbalanceequationsarethensolvedwith
variedfrom0.1sforatwo-traycolumntoabout1.5sfor
the Newton}Raphson method, and a new composi-
a column with 20 trays.
tion pro"le is computed.
(v) The torn variable is re-evaluated and the equations
describingthe reactive distillationprocess are again 3.2.2. Details on initialization
solved (Fig. 2). The starting guess of the torn variable F should be
21
(vi) The algorithm ends up when the global material largeenoughtoaccomodatethedesiredproductionrate,
balance to the column is veri"ed (after convergence considering that conversion is not complete. When the
on the torn variable) or when a non-plausible col- simulationstarts,thisvariableisestimatedsuchthatthe
umn is obtained. sumofallthefeed#owratesofthiscomponentisatleast
50% larger than the speci"ed production #ow rate of
ethyleneglycol.Thisvaluewassetbytakingintoaccount
several items: the stoichiometry of the ethylene glycol
productionreaction(reaction(R1)),theexistenceofasec-
ondary reaction that produces diethylene glycol from
ethyleneglycol(reaction(R2));andthatnotallthewater
fed to the column reacts, some exiting the column with
the products.
Forthesubsequentstepofthesimulationofareactive
or non-reactive distillation column, a starting guess for
thecompositionpro"leisalsonecessary.Fornon-react-
ive distillation the guessed initial composition pro"le is
generally based on feed compositions, considering the
same liquid composition all over the column (Bastos,
1987):
+N F
x " j/1 ij , i"1, 2,C; j"1,2,N. (1)
ij +C +N F
i/1 j/1 ij
This procedure is sometimes used also with reactive
distillationprocesses(Murthy,1984).However,thestart-
ingpro"lecanbetoofarfromthecorrectone,drivingthe
simulation to impossible physical conditions or to non-
convergenceof the iterative process. This can occur, for
instance, when one of the feed species corresponds to
a reactant with a high consumption rate.
Venkataramanetal.(1990)initializethecolumncom-
position pro"le by #ashing the feed at the average col-
umnpressure,takingthedi!erentreactionsintoaccount.
Schenk et al. (1999) solve "rstly a simple equilibrium
model, and next more complex simulations using as
initial guesses the results from the previous simulations.
Fig.2. Simulationalgorithm. The same approach is adopted by Huss et al. (1999) in
M.F.Cardosoetal./ChemicalEngineeringScience55(2000)5059}5078 5065
theirconceptualdesignofnon-equilibriumreactivemod- This normalization is a common feature of algo-
els.LeeandDudukovic(1998)assumetheliquidphaseto rithms employed to solve distillation models (Lee
beanidealsolutionandthestartingtemperaturepro"le & Dudukovic, 1998).
to be linear. (iv) For each one of the other trays, with the estimated
Mostoftheauthors,however,donotreferanymethod compositionsof thetopandbottomtrays,thecom-
to generate the initial conditions necessary to start the positions of components 2}4 are computed using
iterative process. Venkataraman et al. (1990), who in- linear interpolation, and the composition of com-
corporatedtheiralgorithmforthesimulationofreactive ponent1is setnullforalltrays.Ifontraykthereis
distillation columns in the AspenPlus simulator, refer liquid holdup for reaction,the compositionof com-
that a correct initialization by the user can improve the ponent1isrecalculatedusingtheprocedurefortray
success of the starting point. N.
Analgorithmfor the initializationof the composition
With the starting guess for the composition pro"le, the
pro"le for the case study was developed taking into
temperaturepro"lecanbecomputed,asdescribedbelow,
accountthetworeactionsandthespeci"edproduct#ow
as well as all the other variables of the process.
rate.Fromthestoichiometryofthereactions,thefollow-
It should be noted that it is not bene"cial to use the
ing relations are obtained:
previous composition and temperature pro"les as the
N initial point for subsequent iterations, for two reasons:
P "+ F !P !2P , (2)
1 1k 3 4 the optimization algorithm may produce at each iter-
k/1
ationacompletelydi!erentnumberofstages,viz.Ndoes
N notchangesequentially;andthesimplexmoves,even at
P "+ F !P !P . (3)
2 2k 3 4 each "xed value for N, may produce a distributed feed
k/1
which is widely di!erent in composition and location.
The P #ow rate is estimated with Eq. (2) if the total Thus, with the MSIMPSA or MSGA algorithms, it is
4
ethylene oxide feed #ow rate is less than the total water also not bene"cial to start the optimization from
feed#owrate,otherwiseitisestimatedwithEq.(3).This a single-tray column.
estimateisperformedbysupposingthatthereactantwith
thelowertotalfeed#owratehasaconsumptionratenear 3.2.3. Temperature proxle
100%.Theotherreactantproduct#owrateisestimated The temperature on each tray is computed with the
withtheremainderequation.Thecompositioninthetray corresponding stoichiometric vapor phase (Eq. (A.12)).
next to the reboiler and the global extent of reactions Thesolution of this equationis obtainedwith the New-
(R1)and(R2)cannowbecomputedwiththeestimatesof ton}Raphsonmethod.
all product #ow rates, P }P .
1 4
The initial composition pro"le is estimated based on 3.2.4. Liquid and vapor molar yow rates
the following assumptions: Theliquidmolar#owrateo!the"rsttrayiscomputed
with the reboiler material balance for the species corre-
(i) The extent of reactions (R1) and (R2) is the same in
spondingtothespeci"edproduct(Eq.(B.14)).Thevapor
alltrayswithliquidholdup,andnullinalltheother
#owrateintothe"rsttray,< ,isthe‚ fractionvaporiz-
0 1
trays.
ed in the reboiler (Eq. (B.19)). Knowing < , the vapor
0
(ii) The starting guess for the compositionof the liquid
#ow rates o! all the trays of the column are recursively
phaseonthelasttray,N,is95%ofcomponent2and
computed from stage 1 to stage N, using the energy
5% of component 3.
balance over each tray (Eq. (B.13)).
(iii) The composition of component 1 on tray N is re-
Theliquidmolar#owrateso!eachtrayofthecolumn
evaluated if liquid holdup for reaction exists. This
are recursivelycomputedfrom stage N to stage 2, using
composition value must be correctly estimated be-
theglobalmaterialbalanceovereachtray(Eqs.(B.9)and
cause if it is null, then the reaction term is also null
(B.10)). The liquid molar #ow rate into the top tray,
giving rise to infeasible columns. The correction of
‚ is computed with Eq. (B.16).
N‘1
this composition is performed taking into account
the extension value for reaction (R1):
3.2.5. Relaxation method
m "= kr x x , (4) The use of the Newton}Raphson method requires
1N N 1N 1N 2N
asuitablestartingguessforthecompositionpro"le.The
where = is the liquid holdup, kr is the reaction initial starting guess for the composition pro"le is sub-
N 1N
constant, and x and x are the composition of sequently improved by the application of a relaxation
1N 2N
components 1 and 2, in tray N. The composition methodduringapre-speci"ednumberof iterations.The
valuesfor all components are then normalizedsuch relaxationmethodusedis describedinAppendixC,and
thatthestoichiometricequationfortrayNisobeyed. is based on a method described by Bastos (1987) for
5066 M.F.Cardosoetal./ChemicalEngineeringScience55(2000)5059}5078
non-reactive systems. The extension to reactive systems The convergence of the recycle variable, R , is veri"ed
F
can be seen, for example, in Murthy (1984). when
KRn!Rn~1K
F F (e , (6)
3.2.6. The Newton}Raphson method RnF R
The Newton}Raphson algorithm is used to compute
where e and e are appropriate values. The simulation
the compositionpro"le by solving the material balance. C R
process ends up with a feasible reactive distillation col-
The application of this method was based on the algo-
umn, which obeys the global material balance to within
rithm proposed by Bastos (1987), but considering that
10~3(e inEq.(5)).Toavoidprematureconvergenceon
the chemical reaction may be distributed over several C
thetornvariablewithoutconvergenceontheglobalmass
column trays, such as described in Appendix C. With
balance, the convergence criterion for the torn variable
non-ideal VLE, it was also considered that the reaction
was set to a much lower value, 10~10 (e in Eq. (6)).
might be circumscribed to a single tray (Okasinski R
& Doherty, 1998).
Asthecompositionofcomponent1intheliquidphase 3.2.8. Simulation results
isnear zero,theapplicationof therelaxationmethodto The simulation method developed was validated by
this component wouldproduce large errors,not only in performing the simulation of the distillation column re-
thiscompositionbutalsointhecompositionoftheother portedbyCiricandGu(1994),correspondingtothebest
components.Sinceareasonableestimateofthecomposi- con"guration obtained by these authors. For a better
tionofcomponent1isneededin thetrayshavingliquid evaluation of the simulation algorithm proposed in the
holdup,theprocedureadoptedis similartothatused in present work, the same column was simulated in the
the initialization of the composition pro"le, which re- AspenPlusmodularenvironment.Twootherdistillation
quiresthereactionextensionineachtray.Theextension columns obtained during the optimization using the
valuesforthetworeactionswere setequalto thevalues MSIMPSA algorithm were simulated with the Aspen-
previously obtained in the initialization of the column Plus simulator. One was obtained including the con-
composition pro"le. straintsonthemaximumvaluesofvaporandliquid#ow
The composition values obtained in the liquid phase ratesasproposedbyCiricandGu(1994),whiletheother
are then normalized. With the composition pro"le, the was obtained by relaxing these constraints.
temperature pro"le along the column and all the other The simulationof the distillation column reported by
operating variables are calculated (except the reaction Ciric and Gu (1994) corresponds to the water and ethy-
rateswhicharesetconstantduringtheapplicationofthe leneglycolcompositionandtemperaturepro"lesshown
relaxationmethod).Themaximumnumberof iterations in Figs. 3 and 4. These "gures show that the pro"les
was set at 20, after a preliminary study. obtained with the algorithm developed in this work are
The Newton}Raphsonmethod can also be solved for very close to the pro"les obtained by Ciric and Gu, as
componentmolar#ows,insteadofmolarfractions.How- expected since the column model is the same. In what
ever, similar numericalproblems may occur, since com- concerns the liquid and vapor molar #ows, column
ponent #ows may be very small and become negative. height, condenser and reboiler duties, the results ob-
For example, with the data of Ciric and Gu (1994) and tainedbythesetwomethodsareagainthesame,andvery
the MSIMPSA algorithm, while molar fractions vary
from 1.3]10~11 to 7.7]10~5, the corresponding com-
ponent molar #ows vary from 3.5]10~9 to 2.1]10~2.
Also, molar fractions are easier to normalize.
3.2.7. Convergence criteria
The simulation ends up with a feasible distillation
column when the global material balance of all
components is obeyed. If convergence on the recycle
variable is obtained without convergence on the global
materialbalance,thedistillationcolumnisnotafeasible
one.Theglobalmaterialbalanceisassumedtobeveri"-
ed when
S C C N A R B D
+ + Fik#+ tijmjk !Pi 2(ec. Fig. 3. Water composition and ethylene glycol pro"les for the best
(5)
i/1 k/1 j/1 solutionobtainedbyCiricandGu(1994).
M.F.Cardosoetal./ChemicalEngineeringScience55(2000)5059}5078 5067
Fig.4. Temperaturepro"leforthebestsolutionobtainedbyCiricand Fig.5. Waterandethyleneglycolcompositionpro"lesforasolution
obtainedwithactiveconstraintsonthemaximummolar#owrates.
Gu(1994).
close to those obtained with the AspenPlus simulator.
However,atthebottom,theliquidandvapormolar#ow
ratesobtainedwithAspenPlusdi!erfromtheothertwo
simulatorsby 2%, and at the top, by 30%. The column
diametersobtainedbyCiricandGu(1994),bythesimu-
lationmethoddevelopedinthisworkandbytheAspen-
Plussimulatorareveryclosetoeachother,respectively,
1.30, 1.34 and 1.38 m.
The di!erencesbetween the composition andtemper-
ature pro"les obtained by the AspenPlus simulator and
the other two simulators can be due not only to the
di!erent methods used to evaluate the physical proper-
ties, but also to the di!erent assumptions made when
developingthemathematicalmodelforthecolumn.The
energy balance in the model proposed by Ciric and Gu
(1994)andadoptedin this work,Eq. (B.13),will onlybe Fig. 6. Temperature pro"le for a solution obtained with active con-
straintsonthemaximummolar#owrates.
validwithnegligibleenthalpiesoftheliquidstreamsand
constantenthalpyofvaporization.TheAspenPlussimu-
latormostprobablydoesnotfollowanyoftheseassump-
tions,sincetheenthalpiesoftheliquidstreamscomputed
by this simulatorare of the same order of magnitudeas rates. It can be seen that the pro"les obtained with the
the enthalpies of the vapor streams, certainly not negli- proposed simulation method and with the AspenPlus
gible.IntheCiricandGumodel,theenergybalanceused simulator are similar. Also, the column diameter is the
makes the vapor stream remain constant in the bottom same, 1.68 m.
trayswherethereactiontermsarenull.Thisisthereason Thecolumnspeci"cationobtainedwhenthemaximum
why the application of this model gives constant vapor liquid and vapor molar #ow rates are relaxed is com-
andliquid#owratesforthe"rstfourtrays,wherethereis paredwiththeresultsobtainedwiththeAspenPlussimu-
no reaction. On the contrary the AspenPlus simulator latorinFigs.7and8.Itcanbeseenthatthepro"lesare
givesrise to di!erent #ow rates, which changesthe "nal now nearly the same all over the column.
simulationresult.Inspiteofthedi!erent#owratesinthe The pro"les obtained with the proposed simulation
toptrays,thepro"lesobtainedwiththedi!erentmodels algorithmand withthe AspenPlussimulatorare similar
aresimilartoeachother,probablybecauseinthosetrays forthelasttwocon"gurations.Thiscanbeduetothefact
thereis reactionwhichis themainfactor in establishing that in both con"gurations the reaction occurs in all
the tray compositions. trays.
Figs. 5 and 6 show the pro"les corresponding to Inconclusion,thevalidationofthesimulationmethod
a simulation obtained with the MSIMPSA algorithm developed was performed by comparing the results ob-
with active constraints on the maximum molar #ow tainedwiththosefromtheAspenPlussimulator,forthree
5068 M.F.Cardosoetal./ChemicalEngineeringScience55(2000)5059}5078
3.3. Optimization of the production of ethylene glycol via
reactive distillation
3.3.1. Ideal vapor}liquid equilibrium
Thesynthesisofanon-equilibriumreactivedistillation
column for the production of ethylene glycol, as de-
scribedabove,wasproposedandsolvedbyCiricandGu
(1994). These authors used a modi"ed GBD algorithm
(Geo!rion, 1972).
Ciric and Gu (1994) state that the GBD algorithm
simplyrequiresselectinganinitialsolutionvectorofthe
discretevariables,viz.thenumberof theoreticaltrays in
thedistillationcolumn.IntheGBDalgorithm,amaster
problemselectsthenumberoftrayswhileaprimalprob-
lemoptimizesthecontinuousvariables.Forhighlynon-
Fig.7. Waterandethyleneglycolcompositionpro"lesforasolution linear MINLPs, it is standard practice to solve a third
obtainedwithoutconstraintsonthemaximummolar#owrates.
optimization problem, between the master and the pri-
mal, to "nd a feasible set of continuous variables by
minimizing the total constraint violation.
TheMSIMPSAalgorithm,developedforthesolution
of MINLP problems was applied to this example. With
theMSIMPSAalgorithm,thestartingsolutionvectorsof
both continuous and discrete variables are randomly
generated and correspond most probably to infeasible
points. However, since this algorithm is an infeasible
pathmethod,itdoesnotrequireneitheraninitialfeasible
pointnorsolvingintermediateoptimizationproblemsfor
feasibility. For highly constrained and non-linear NLP
or MINLP problems, a dynamic penalizing scheme is
available for "nding feasible solution sets, while main-
taining convergence to the global optimum (Cardoso et
al., 1997). This dynamic penalizing scheme roughly
doubles the CPU time, but it was not needed in the
present case.
Fig.8. Temperaturepro"leforasolutionobtainedwithoutconstraints The MSIMPSA algorithm adopts the stochastic
onthemaximummolar#owrates.
scheme proposed by Dolan et al. (1989) for the genera-
tion of discrete con"gurations. Since this scheme only
allows a one unit variation up or down to the next
di!erent con"gurations of a reactive distillation column discretevalue,largejumpsinthediscretesearchspaceare
for the production of ethylene glycol. All column speci- notallowedfromoneiterationtothenext,andthismay
"cationscorrespondtofeasiblecolumndesigns.Eventhe weaken the optimization process. To allow larger vari-
con"guration obtained with the MSIMPSA algorithm ations,thenumberoftraysinthecolumnwasrepresent-
withouttakingintoaccountconstraintsontheallowable ed by a binary expansion (Cardoso et al., 1997).
liquidandvapormolar#owrates,whichare50%larger Firstly,the MSIMPSAalgorithmwas appliedinclud-
thanthelimitsimposedbyCiricandGu,correspondsto ing the additional constraints on maximum molar #ow
a feasible column. The AspenPlus simulator, with the rates as stated by Ciric (1995), and again an optimum
same con"guration, could also end up with a feasible con"gurationwith10theoreticaltrayswasobtained.The
column. This agrees with the heuristic rules given by best solution found by the MSIMPSA algorithm has
Douglas(1988),whichstatethatifthereisneither#ood- a total cost of 15.26]106 US$ yr~1, similar to the cost
ing nor liquid entrainment there are no limits to the found by Ciric and Gu of 15.69]106 US$ yr~1. The
molar#owrateswithinthecolumn.Theonlyconstraints solution vectors are however di!erent, and the solution
statedbythoserulesareamaximumlimitof53minthe obtainedwith the MSIMPSA algorithm corresponds to
column height, a minimum limit of 0.45 m in the dia- activeconstraintsonthemaximummolar#owrates.The
meter, and that the column diameter should be lower twosolutionsarecomparedinTable2andinFigs.9and
than5% of the columnheight.Noneoftheselimits was 10.Itcanbeseenthatwiththesolutionobtainedinthis
violated in any of the above cases. work there is reaction in all trays and the feed is more