- Email: [email protected]
Simulation of continuous reactive distillation by a homotopy-continuation method
Compur.them. Engng,Vol. 12, No. 12, PP. 1243-1255, 1988 Printedin Great Britain.All rightsreserved
009% I354/88 S3.00 + 0.00 Copyright6 1988 PergamonPressplc
SIMULATION OF CONTINUOUS REACTIVE BY A HOMOTOPY-CONTINUATION
DISTILLATION METHOD
Y. A. CHANG and J. D. SEADER~ Department of Chemical Engineering, University of Utah, Salt Lake City, UT 84112, U.S.A. (Received 25 March 1988; final revision received 5 May 1988; received for publication 19 May 1988) Abstract-Reactive distillation was studied theoretically for the esterification of acetic acid with ethanol to produce ethyl acetate and water, using second-order reversible kinetics. Because of the high degree of interaction among the iteration variables, convergence was difficult-to-impossible to achieve with Newton’s method. Consequently, a robust homotopy-continuation method was applied to solve the simultaneous nonlinear equations used to model the reactive-distillation system.Case studiesfor several different specifications were conducted to show how design variables affect the reactive-distillation system.
INTRODUCllON
Because of the potential importance of reactive distillation, a number of publications, both patents and articles, have considered the problems of its application and performance. Prior to the advent of digital computers, the literature dealt mainly with proposed applications and simplified calculational procedures. Backhaus (1921), in perhaps the first reference to reactive distillation, published a patent describing an esterification process carried out in a distillation column. Leyes and Othmer (1945), as well as Bermans et al. (1948), treated an esterification reaction in a distillation column and utilized experimental vapor-liquid equilibrium data and kinetic and chemical equilibrium relations to get numerical solutions. Belck (1955) made assumptions to simplify calculations for hypothetical reversible reactions carried out in a single distillation column, Since 1970, calculation procedures suitable for implementation on digital computers, for rigorously solving the mass, equilibrium, and energy balances in multicomponent, multistage reactive-distillation columns have been reported. The procedures, which involve the solution of large sets of combined nonlinear and linear equations, as well as applications, are summarized in Table 1. From a numerical calculation procedure, these previously reported studies can be divided into three categories: (1) methods using tear variables (bubble-point method); (2) relaxation techniques (dynamic approach); and (3) methods incorporating the Newton-Raphson method. For the tearing method, it is not necessary to compute derivatives, and the tridiagonal matrix algorithm can be implemented to decrease computer storage space and increase computation efficiency. However, when differences in boiling point between components is too large, when the kinetics are comtAuthor to whom all correspondence should be addressed.
plex, or when the liquid solutions are highly nonideal, the tearing method is very difficult to converge. The relaxation method has better stability, but may require many more iterations to obtain a solution. Although the Newton-Raphson method converges quickly from a suitable starting guess and can be applied more widely to different specifications, it needs large storage and involves complicated derivative operations. Convergence is not guaranteed for any of the procedures used in previous studies. Table 1 shows that the following reaction processes have been reported in the literature: esterification, saponification and the manufacture of propylene oxide from chlorohydrin. The esterification of acetic acid with ethanol is the most widely studied case. Not included in Table 1 are the required specifications, which are essentially the same for all references listed. The user must specify the number of stages, feedplate location, feed condition, reflux ratio, and distillate or bottoms flowrate. In this study, it was desired to develop and demonstrate a robust numerical procedure for reactive distillation, involving the use of homotopycontinuation in the manner of Wayburn and Seader (1984). The procedure was designed to permit great flexibility with respect to column configuration (including interlinked systems), specifications and choice of thermodynamic properties.
MATHEMATICAL
MODEL
Consider an N-stage column separating C components, as shown in Fig. 1, where stage 1 is a total condenser that produces a saturated liquid and stage N is a partial reboiler. Furthermore, let chemical reactions in the liquid phase occur at each stage. Assume that: (1) the process has reached a continuous steady state; (2) each stage is a perfectly mixed stirred-tank reactor (thus, the liquid composition at each stage is homogeneous and equal to
1243
Carra er u/. (1979)
Combine tearing method (choose the profile of liquid mol fraction of [OH-] as tearing variables) and Newton-Raphson method
rithm
Reaction type not stated
(i) propylene chlorohydrin t 0.5 calcium hydroxide = propylene oxide t 0.5 calcium chloride (ii) propylene oxide _t water = propandiol Vapor pressure, Antoine equation, activity coefficient by van Laar equation, salting Evaluate new profik of liq uid mol fraction of [OH-]
Not stated
effect
(i) AtB=CtD (ii) normal alcohol t carboxylic acid + eater t water (iii) saponification, acetate t steam = alcohol t acetic acid Reaction type not stated Not stated
Combine Newton-Raphson method and banded structure of Jacobian matrix algo-
AcOH t EtOH = AcOEt t H,O Second-order Reversible
K is a function of temperature and liquid composition
vergence
Kaibel er al. (1979)
method
Multi 8-n method of con-
Newton-Raphson
Komatsu and Holland (1977)
position
Second-order Reversible
Reaction rate constant Arrhenius equation
(i) AcOH t EtOH = AcOEt t H,O (ii) AcryOH t EtOH = AcryOEt t Hz0 Second-order Reversible
AcOH + EtOH = AcOEt t H,O
Relaxation factor and corm section for liquid com-
Relaxation technique
Komatsu (1977)
Vapor pressure, Antoine equation, activity coefficient by NRTL equation
r =~,~~l~~~lI~l-~,M~~l~
Reaction-rate quation
A+B=C+D
AcOH t EtOH = AcOEt t H,O Second-order Irreversible
Reaction and its type
K is a function of temperature and liquid
Relaxation factor
Relaxation technique
Jelinek and Hlavacck (1976)
Not stated clearly
equation
Case 2: vapor pressure, Antoine equation, activity coetiicient by modified Margules
Case I: K is a function of temperature only
Thermodynamic correlation for vapor-liquid equilibrium
systems
composition
Damping factor
method
NewtonRaphson
Nelson (1971)
Convergence method Normalize the liquid cornposition and use modified Mullet’s method for convergence of temperature
procedure
Choose temperature as tear variable and use the tridiagonal matrix algorithm for the solution of the linearized material balance equations
Calculation
Suzuki er al. (1971)
Reference
Table I. Rigorous methods for computing reactive distillation
Homotopy-continuation
for reactive distillation
1245
the composition of the liquid leaving the stageCSTR assumption); and (3) the vapor and liquid leaving any stage are in physical equilibrium. The extent of liquid-phase reaction at each stage is governed by reversible and/or irreversible kinetic rate expressions. Figure 2 shows a schematic representation of a general reactive-distillation stage. The three types of functions which describe the physical and chemical processes on stage i are: 1. Component
material balances:
M,=(l
-ts,)L),,-t(l -f;j-gij-
+si)fij-ui,,J-~l,_,J 2iRYN.s ui c v,mr’in “-=I
j=l,...,C
(1)
where the last term accounts for chemical reactions. 2. Energy balances: Ei = (1 + S,)H,
+ (1 + si)hi -ha-,
Hi, L -
b-j -
&,
-
Qi,
(2)
where unlike the component material balances, it is not necessary to take the heat of reaction into account because enthalpies are referred to the elements rather than to the components. 3. Equilibrium relationships: Q, = K,, iii 2, -
vy ,
j-1
9--m, C
(3)
where Q,J is derived from the definition of Key: K,j = ti x,j
or
Kijx,
-
y, = 0,
(4)
and x and y represent component mol fractions. Specifications for equations (l-3) include the total number of stages, stage locations of all feeds, sidestreams, and heat exchangers; all stage pressures and liquid-phase holdup volumes; and complete specification of each fresh feed. This leaves N specifications that must be made before solving the N(2C + 1) equations for the corresponding N(2C + 1) iteration variables. If the N specifications are N stage heattransfer rates, then the equations are solved for 2CN component molar flows and N stage temperatures. However, because condenser and reboiler duties are strongly dependent and can generally be specified independently over only narrow unknown ranges, other specifications are more desirable. The total condenser stage is a special case which is not represented by equation (3). Let the component molar flows of the liquid distillate be called If the reflux ratio is specified, the “I.1 3U1\2, . . . 3 h.c-
Y. A. CHANG
1246
and J. D. SFZADER
Fig. 1. Configuration for reactive-distillation column.
top-stage
energy
balance
equation
is changed to a
For the case of all r, = 0, Naphtali and Sandholm (1971) showed that ordering the above equations as follows for each stage: Ei, M,,, . . . , M ,,=, and differentiating the equations in the Q,,, ...qQt.,; order: Ti_ ,,I, -,,, ,.‘.,~i-,,c,~i.,,..-,~i,c, Ti;.,4.L,.-., T.+, leads to a block-tridiagonal 6.,90,+1.19....~i+1.e, Jacobian matrix with minimum bandwidth. In the procedure applied here, the block-reduction scheme of Waybum (1983) was used to deal with the matrix, making the problem of bandwidth irrelevant. To avoid unnecessary pivoting, Waybum placed the derivative of the energy balance equation with respect to temperature on the diagonal. The equations were ordered M,,, , . . . , M,,,, E,, Q;,, , . . . , Q,,, and then differentiated in the order of Naphtali and Sandholm (1971). To achieve flexibility with respect to specifications, the method of Ferraris (1981), as implemented by Waybum and Seader (1984) was used. This involved adding specification equations as rows at the bottom of the standard set and adding the additional unknowns as columns at the right side of the Jacobian matrix. The total condenser, modeled by the reflux-ratio equation (5), the bubble-point equation (6) and mol-fraction identities (7) were considered as standard equations. Chemical-reaction rates, rin. in the liquid phase were modeled by reversible power-law kinetic rate equations. Vapor molar holdup and vapor-phase chemical reactions were neglected. Because the material balance equations (1) are in the form outminus-in, the term represents the moles of component j produced in the liquid phase. Therefore, v,,,, the stoichiometric coefficient for component j in the n th chemical reaction, is positive if component j is created by the nth reaction. The reaction-rate equations were written in terms of concentrations of the relevant components:
reflux-ratio equation: E,=L,-RV,,
(5)
where R is the reflux ratio. The distillate composition, D1.1r”,,2, . , qE, is the same composition as in the liquid leaving the condenser and, therefore, the equilibrium relations are not independent. the equilibrium equation for component 1 is replaced by the bubblepoint equation:
s iY
ij
t
Pi
Q,., = 2 Wi,p4.p)~ - vl.p* I p-1
(6)
where component 1 must be present in the distillate in nontrivial amounts. Like the energy balance equation, the other equilibrium equations, Q,*-Q1=, reflux are substituted for by the component equations; viz.
I i-l./
' ij
Ti
7
'ij
Y.
r+rj
l ij
Sij
Fig. 2. General reactive-distillation stage.
1 ij
Homotopy-continuation for reactive distillation
1247
of temperature only. By making initial guesses for the R, values, the modified Wang-Henke linear equations could be readily solved for initial composition profiles. SIMULATION
=xA,exp P
(8)
where C,= i,/(Livk) is the concentration of component q. The variable A, is the Arrhenius preexponential factor, E, is the energy of activation, R is the universal gas constant and 7’, is the temperature at stage i. This reaction model can deal with either irreversible or reversible reactions. For example, for the reversible reaction: A+B=C+D
(9) where p = 1 indicates the forward reaction, p = 2 is the reverse reaction, /cl is positive and kz is negative. To achieve robustness in convergence, the modeling equations, call them f(x) = 0, were embedded into a Newton homotopy: H(x, 1) =f(x)
- (1 - ~)f(x’)
= 0,
(10)
where t is a scalar homotopy parameter, which is varied from 0 to 1 as continuation is used to track the path from the starting guess x0 (at t = 0) to the solution at t = 1. Tracking was by means of differential arclength continuation using an Euler predictor and a Newton corrector as described by Waybum and Seader (1984). Initially, a maximum step size equivalent to moving from t = 0 to 2 = 1 was attempted. This is Newton’s method, which was iterated with and without a line search. If it failed to converge, for reasons discussed later, continuation was employed using the step-size algorithm of den Heijer and Rheinboldt (198 1), which was used by Lin, Seader, and Waybum (1987). If turning points occurred along the path, they were tracked by the method of Chan (1984), as described by Lin et al. (1987). For providing the starting guesses for the iteration variables (unknowns), a modification of the first step in the Wang and Henke (1966) method was used. For each stage, equation (1) is written in mol-fraction form as: Mj=(l
+
si)KYt,+ -
C1
Li_,,ixi_
+si)Llxij-
,.i-
C+l.jYt+l.j
Fjzii -
R,=
0.
(11)
By combining equations (4) and (11) and selecting T, and Vi as tear variables, the resulting NC equations
are linear in the unknown liquid mol-fractions provided that K values are composition independent. Since this was not the case for the system studied here, the method of Suzuki er al. (1971) was used to estimate starting guesses for K-values as a function
.STUDIES
All the simulation case studies reported here were made on the four-component reactive system involving acetic acid (l), ethanol (2), water (3) and ethyl acetate (4). with the sole reaction being: acetic acid + ethanol - ethyl acetate + water. Since at least 1921, this system has been considered to be suitable for reactive distillation. Because of the non-ideal nature of the liquid-phase mixture caused by the polar molecules, it was necessary to give careful attention to the thermodynamic correlations used to compute vapor-liquid equilibrium. In particular, acetic acid tends to form a dimer and a trimer in the vapor phase and it is well-known that binary mixtures of ethyl acetate-water, acetic acid-water, and ethyl alcohol-water, as well as the ternary mixture of ethyl acetate-acetic acid-water, all form azeotropes. K-values were computed from Kj = yj Pg/ P where yj was computed from a modified Margules equation developed especially by Suzuki et al. (1970) for the system studied here. Suzuki et al. (1970) also present equations for correcting the vapor-phase composition for the dimer and trimer of acetic acid as well as the constants for each component in the Antoine vapor pressure expression. For vapor and liquid enthalpy, the correlations presented by Izarraraz et al. (1980), which are referenced to the elements, were used. Details of the thermodynamic and kinetic equations, together with references to listings of the constants in the equations, are given in the Appendix. The first simulation, Case 1, was made using the same specifications as Suzuki et al. (1971) and Izarraraz et al. (1980), as shown in Fig. 3. There are 13 theoretical stages, all at 1 atm pressure, counting the partial reboiler and the total condenser. The feed rate is 107.6lbmol h-‘, with a liquid distillate of 20.8 lbmol h-‘, giving a bottoms of 86.8 lbmol h-l. The feed is preheated to its bubble point at feed tray pressure and the reflux ratio is 10. The feed is fed to the sixth stage from the top and has the following mol-fraction composition: acetic acid z, = 0.4962, ethanol z2 = 0.4808, water z, = 0.0229 and ethyl acetate z, = 0.0. Holdup volumes are 266.95 and 80.085 ft3, respectively, for the reboiler and each of the stages l-l 2. Although these specifications are the same as those of Suzuki et al. (1971) and Izarraraz et ai. (1980), the results of this study as discussed below, were different, but much closer in agreement to Suzuki et al. (1971). Comparisons with Izarraraz er al. (1980) are not included here because they used considerably different vapor-liquid equilibrium ratios as compared to those used here and by Suzuki et al. (1971).
Y. A. CHANG
1248
P=lotm
Retlfm Ratio-10.0 Hddups : Reboller=2&5.9Bcuf Stage 1-12 [email protected] /each 5
Saturated
Uquld
,brn0l AcOH EtOH W&m AcoEi
9
I
I 6
53Ao 51.74 2.46 0.00 107.60
12
E+ 13
1
Bottoms
Fig.
3.
B&B LbmoC h-’
Specifications for Case
1
Table 2 compares the computed temperature, liquid flow and vapor flow profiles with those of Suzuki et al. (1971). For stages above the feed, the profiles computed in this study agree well with those of Suzuki et al. However, as the reboiler stage is approached, the profiles begin to deviate. The deviations are caused not by differences in K-values but by differences in the reaction-rate expressions. As shown in Table 3, K-values for each component at each stage are almost identical based on liquid compositions from Suzuki ef al. (1971) who did not allow for the backward reaction, while this study did. As seen in Table 2, the resulting reaction-rate profiles are quite different. Small negative production rates of ethyl acetate were computed here for the upper four stages, while Suzuki et al. show zero rates there. For most of the lower nine stages, the reaction rate is greater than that of Suzuki et al. due to the use of Table
2. Comparison
of Case
and J. D. SEADER
different kinetic constants for the forward reaction as obtained from Izarraraz et al. (1980). For Case 1, the extent of reaction and separation achieved are quite unfavorable. The extent of reaction, based on the limiting reactant, is only 28.3% and the distillate contains appreciable ethanol. Accordingly, attention was directed toward changes in column configuration and operating conditions which might increase conversion and improve the separation. The most common industrial schemes are to use an excess of one reactant (usually the least expensive one) and to recycle the unconverted material. If the recycle configuration proves to be uneconomical, it is possible to use a countercurrently operated sequence of continuous stirred-tank reactors or to remove continuously one or more of the products from the reaction zone. In such devices, the reaction is thus driven toward completion. In some situations, a staged distillation column is a very efficient reactor. By use of the technology of reactive distillation, it is possible to get high reaction conversion and separation simultaneously. For example, consider a reversible reaction of the type, A + B = C + D, where B and D are more volatile than A and C. Liquid rich in A is fed to a stage between the top and the middle portion of the column. Vapor rich in B is fed to a stage between the middle portion and the bottom of the column. Reactant B is absorbed into the liquid phase where the reaction takes place and product D is stripped from the liquid phase and carried out the top of the column. In the case of esterification, the order of volatility is ethyl acetate, ethanol, water and acetic acid. Thus, it is expected that the middle portion of the distillation tower is the chief reaction zone. The rectifying section fractionates the ethyl acetate out of the acetic acid, and the stripping section removes alcohol from water. Ideally, the ethyl acetate is the distillate and water is the bottoms product. However, over a wide range of composition, as shown in Table 3, ethanol and water do not differ greatly in volatility, making it difficult to move ethanol up the column. Case 2 was designed to gain a better understanding of reactive distillation for the acetic acid-ethanol system. The only difference from the specifications of I results with Suzuki
et al. (1971) Suzuki
This study stage No. 1 2 3 4 5 6 7 8 9 10 I1 12 13
(Ibm&‘) 162.8 163.1 163.5 164.2 166.3 176.1 176.5 177.0 177.6 178.3 180.0 185.4 199.1
208.0 206.6 204.5 201.7 198.5 303.7 302.4 300.9 299.3 297.8 297.2 303.3 86.8
(lbmCL) 20.8 228.8 227.4 225.3 222.5 219.3 216.9 215.6 214.1 212.5 211.0 210.4 216.5
Reaction rate (lbmol h- ‘) -0.019 -0.017 -0.016 -0.013 0.05 0.56 0.60 0.65 0.72 0.80 0.96 I .43 8.03
(a,) 163.0 163.3 163.7 164.4 166.5 176.2 176.6 177.1 177.7 178.4 180.0 184.9 197.9
(lbm:
h-‘)
208.0 206.6 204.6 201.8 198.7 303.8 303.0 302.1 301.3 300.5 300.6 306.5 86.8
er al. (1971) (lbm:h-‘) 20.8 228.8 227.4 225.4 222.6 219.5 217.0 216.2 215.3 214.5 213.7 213.8 219.7
Reaction rate (Ibmol hW’) 0.00 0.00 0.00 0.00 0.00 0 27 0.29 0.30 0.32 0.35 0.45 0.94 11.66
Homotopy-continuation for reactive distillation Table Acetic
10 II 12 13 Case Case
Table
No. 1
EtOH
-
-
f 8 1: 11 12 13
11,12-i=
0.9120 0.8965 0.8865 0.9061 1.0797 1.0813 .085 1 1.0918 1.1042 1.1357 1.2383 1.5011
I
1
2
acetate 2
053
0.93 0.9 1 0.88 0.89
0.92 0.9 1 0.88 0.89
1.12 1.06 0.91 0.99
t.11 .05 0.91 0.98
1.05 1.09 1.15 .24
0.05 0.13 0.13 0.13 0.14 0.14
0.05 0.13 0.13 0.13 0.14 0.14
0.90 1.07 I .08 .oa
0.90 .07 1.08 1.08 1.09 1.10
0.88 1.07 1.07 1.07 1.07 1.07
0.88 1 .Ol 1.07 1 .O7 1.07 1.07
1.42 1 .I32 1.86 1.92 1.98
1.05 1.09 1.15 1.24 1.42 1 .I32 1.86 1.92 1.98
0.14 0.20 0.35
0.15 0.20 0.35
1.13 1.23 1.48
1.13
1.08 1.14 1.34
1.08 1.14 I.35
2.06 2.16 2.29 2.45
2.06 2.16 2.28 2.44
eauilibrium
AcOH
Ethyl
I
t I.09
1.10
t
1.23 1.48
I
1 = Suzuki er al. (1971). 2 = This study.
Vapo-liauid
0.036 1 0.0363 0.0380 0.0491 0.1313 0.1322 0.1339 0.1366 0.1414 0.1551 0.2077 0.3650
2 3 4 5
CACE.
4.
(K-values)
water
0.04 0.04 0.04
Case 1 was the feed composition, which was changed to an equimolar bubble-point mixture of acetic acid and ethanol. It is noted above that ethanol and ethyl acetate are more volatile than acetic acid and water. Thus, the reaction should occur mainly in the middle portion of the distillation tower. By examining the vapor-liquid equilibrium data shown in Table 4, the volatility of acetic acid appears to be suitable. It is by far the least volatile component in the system, with K values much smaller than for the other three components. Unfortunately, the ratios of KEIOH to KHlo vary only from 0.87 to 1.06, which are too close to 1. Ethanol is not much more volatile than water and, in fact, less volatile in the upper stages of the column. The results of Case 2 are listed in Table 5. The temperature profile of the column is far from linear, with two temperature jumps evident. One occurs at the feed point between stages 5 and 6 and the other appears at the reboiler. The change in temperature for the top five stages is only 3.S”F. The change in temperature for stages 6-l 1 is just 0.8”F. Between the 12th stage and the reboiler, the temperature difference is 13.8”F. Since esterification is an endothermic reaction, equilibrium shifts at higher temperatures toward the products. As seen in Table 5, the column temperature profile causes the chemical reaction to take place in the lower portion instead of the middle part of the distillation tower. The volatility of acetic acid coincides with the requirement that ideally it should be less volatile than
Sync
ratios
2
-
9
equilibrium
Ethanol
2
1
I 2 3 4
8
Vapor-liquid
1
stage
5 6 7
3.
acid
1249
4
data
for
Case
2
Water
AcOEt
1.0439 0.9691 0.9021 0.8823 .0730 1.0754 1.0777 1.0794 1.0814 1.0961 1.1795 1.4135
1.1001 1.1590 I .2496 I.4161 1.8123 1.8598 1.9178 1.9871 2.0697 2.1720 2.3012 2~4652
-
I
P 8
B
5 = t
B 3
;;‘
I
Y. A. CHANG and J. D. SEADER
1250
Table
Case
6. Specification
Differences from Case k
compare with Casek
3
2
Two
liquid
4
3
5
4
One Liquid feed one vapor feed Feed locations
6
4
Feed locations
7
2
Feed compositions
8
2
Feed compositions
feds
9
6
Feed compositions
10
6
Feed compositions
II
9
Holdup
12
9
Distillate
13
2
Number Holdup
rate of stages
ethanol and ethyl acetate. But the K values of acetic acid are too small compared to those of the other components. Most acetic acid tends to flow downwards and appears in the lower portion of the column only; very little exists in the upper and middie parts. In the middle portion, the liquid mol-fractions of acetic acid are not more than 0.2. Therefore, the production rate of ethyl acetate is limited by the deficiency of the reactant acetic acid. The liquid composition of ethanol increases from the condenser to a maximum at the 5th stage. Then it decreases from the 5th stage to the reboiler. The distribution of ethanol is favored for the reaction taking place in the middle portion of the column. Since there is a strong interaction between the reaction and separation, each specification of the system exerts a complex effect that is difficult to predict. Accordingly, as listed in Table 6, a number of variations to the specifications of Case 2 were studied in an attempt to achieve a higher conversion and a better separation. Resulting product distributions for all cases are given in Table 7. In Case 3, two separate saturated liquid feeds are used. Acetic acid is fed to stage 6 and ethanol is fed to stage 12. The separation and the reaction rate do not improve much. The big difference between Case 2 and Case 3 is the compositions of ethanol and water in the distillate. Of 20.8 lbmol h-l distillate, there are 9.86293 lbmol h-’ ethanol and 0.85561 lbmol h-i water in Case 2, while in Case 3 there are 6.62495 lbmol h-’ ethanol and 3.11422 lbmol h-r water. However, the total amount of ethanol and water just changes slightly. In the ideal case, it is preferred to have two feeds, one liquid and the other vapor. In Case 4, saturated liquid acetic acid is fed to stage 6 and saturated vapor enters at stage 11. Compared with Case 3, the two-
for all -
Specification
changes
Liquid AcOH -6th stage Liquid EtOH + 12th stage Liquid AcOH +6th stage Vapor EtOH - 11 th stage Liquid AcOH + 5th stage Vapor EtOH - 9th stage Liquid AcOH -4th stage Vapor EtOH - 10th stage 58.8 lbmol h- ’ AcOH 48.8 lbmol h-’ EtOH 48.8 lbmol h-’ AcOH 58.8 lbmol h-’ EtOH 58.8 lbmol h ’ AcOH 48.8 lbmol h - ’ EtOH 48.8 lbmol h ’ AcOH 58.8 lbmol h-’ EtOH stages l-l 3 = 94.459 fts/cach stage Distillate = 16.00 lbmol h-’ Reflw. ratio = 13.00 Number of stages = 24 Holdups: stage 1 80.085 ft3 stage 2-23 40.0425 ft’/each stage 24 266.95 ft”
stage
liquid-feed case, the reaction rate in the middle portion of the column is a little greater than for Case 3, but the separation is not improved. In Cases 5 and 6, the two feed stages are moved upward from Case 4. Comparing these three cases in Table 7, the production rate of ethyl acetate in Case 6 is about 10% more than that of the other two cases. As mentioned earlier, acetic acid has a strong tendency to flow downwards. If the feed stage for acetic acid is moved upwards, mol-fractions of acetic acid in the upper portion should increase, causing the forward reaction of esterification to increase. In the upper part of the column, either the negative values of production rate of ethyl acetate decrease or the positive values increase. Therefore, the total production rate of ethyl acetate increases but not significantly. One-feed cases, Cases 7 and 8, were made for comparing the effect of feed composition changes with Case 2. Two-feed cases, Cases 9 and 10, were chosen to compare with Case 6. In Table 6, feeds in Cases 7 and 9 consist of an excess of acetic acid. Flowrates of feeds in Cases 8 and 10 consist of an excess of ethanol. Conclusions are the same for the one-feed and the two-feed cases. The greater the acetic acid-to-ethanol ratio, the more ethyl acetate is created and vice versa. Since most of the reaction takes place at the higher temperature in the lower portion of the column, an increase in the amount of acetic acid causes the reaction rate to increase. Because volatilities of ethanol and acetic acid differ, ethanol is distributed more evenly than acetic acid. In the main reaction zone, the increased concentration of ethanol does not compensate for the loss of acetic acid. From these cases, distillate compositions do not change significantly with feed composition change. When the feed flowrate of ethanol was changed from
Homotopycootinuation Table Component
AcOH EtOH Water AcOEt
7. Final
product
Feed
53.4000 51.7400 2.4600 0.0000 107.600
AcOH EtOH Water AcOEt
53.8000 53.8000 O.lYJOCt O.OOC+l 107.600
AcOH EtOH Water AcOEt
53.8000 53.8000 0.0000 0.0000 107.600
AcOH EtOH Water AcOEt
53.8ooO 53.8IXtCI 0.0000 o.oooo 107.6OfJ
AcOH EtOH Water AcOEt
53.8000 53.8000 o.oooo 0.0000 107.600
AcOH EtOH Water AcOEt
53.8000 53.8000 0.0000 O.OODO to7600
AcOH EtOH Water AcOEt
58.8000 48.8000 0.0000 0.0000 107.6000
AcOH EtOH Water AcOEt
48.8000 58.8Oaa 0.0000 0.0000 107.6LM
AcOH EtOH Water AcOEt
58.8000 48.8000 O.OOOfJ O.OOOCl 107.600
AcOH EtOH Water AcOEt
48.8000 58.8000 0.0000 0.0000 107.600
AcOH EtOH Water AcOEt
distribution Distillate (Ibmol
case 1 0.0019 9.1843 1.3121 IO.3018
13 cases
BOttOll-lS
RePcti0tl
h-‘)
39.6524 28.8100 14.8937 3.4440
20.800
86.800
Case 2 0.0011 9.8629 0.8556 10.0804
40.2491 30.3873 12.6942 3.4694
20.800
86.800
Case 3 0.0064 6.6250 3.1142 I 1.0544
39.8220 33.2034 10.8574 2.9172
20.800 Case 4 o.On5 1 7.1773 2.6834 10.9343
39.3644 32.1921 11.7472 3.4%3
20.800
86.800 39.0127 30.9404 12.7823 3.4440
20.800
86.800
Case 6 0.0057 7.2238 2.4363 11.1342
38.8323 31.6141 12.5258 3.8279
20.800
86.800
Case 7 0.0014 9.2747 1.0383 10.4856 20.800 Case 8 0.0009 10.5905 0.7217 9.4869 20.800 Case 9 0.0067 6.7776 2.6414 11.3744 20.800 Case to 0.0050 7.7564 2.2740 10.7647 20.800
20.800
- 13.746 - 13.746 13.746 13.746
- l3SM - 13.550 13.550 13.550
- 13.972 - 13.972 13.912 13.972
86.800
Case 5 0.0033 8.0456 2.0017 10.3018
case 11 58.8OOoO 0.0081 48.80000 6.8466 O.OOODO 2.8459 O.OOOOO 11.0994 107.600
for
for reactive distillation
44.5078 25.2345 t3.2525 3.8052
- 14.431 - 14.431 14.431 14.431
- 14.784 - 14.784 t4.784 14.784
- 14.962 - 14.962 14.962 14.962
- 14.291 - 14.291 14.291 14.291
86.800 36.2464 35.6568 11.8310 3.0658
- 12.553 --12.553 12.553 12.553
86.800 43.285 1 26.5143 12.8668 4.1338
- 15.508 - 15.508 15.508 15.54!8
AcOH EtOH Water AcOEt
58.8000 48.8ooO 0.0000 0.0000 107.600
AcOH EtOH Water AcOEt
53.8000 53.8000 0.0000 0.0000 107.cwoct
- 14.219 - 14.219 14.219 14.219
86.800 45.4601 28.6216 to.4859 2.2324 86.800
- 13.332 - 13.332 13.332 13.332
Case 12 0.0051 5.0929 1.9198 8.9823
44.8661 29.7783 12.0090 4.9466
16.000
91.600
Case I3 0.0005 9.3724 0.4013 1 t .0259
39.9264 30.554s 13.4718 2.8473
20.800
- 13.929 - 13.929 13.929 13.929
- 13.873 - 13.873 13.873 13.873
86.800
to 58.8 lbmol h-l, ethanol appeared in the distillate in almost the same amount, with the excess appearing in the bottoms product. In Cases I-10, the volume holdup ratio for the reboiler to each stage was 10 to 3. Approximately half of the ethyl acetate was produced in the reboifer. In Case 11, the holdup of each stage is specified to be equal to that of the reboiler (94.459 ft’), with the total holdup of the column unchanged. Because the temperature and the concentrations of components (especially acetic acid) in the reboiler is advantageous for the forward reaction, the total rate decreases. Due to the change of holdup, the reaction is distributed more uniformly over the column. The separation is slightly improved, with a 13% lower production rate of 13.332 Ibmol h-’ of ethyl acetate. In Cases l-11, the overall reaction rate is computed to be approx. 14 lbmol h-‘. With a distillate specification of 20.8 lbmol h-‘, the purity of the resulting ethyl acetate distillate is poor. In Case 12, the distillate is reduced to 16.0 lbmol h-’ in an effort to increase purity. At the same time, the reflux ratio is changed to 13.0 to keep the vapor and liquid traffic unchanged. These changes improve the purity of ethyl acetate in the distillate by less than 5% even though the overall reaction rate is maintained at 141bmol h-l. In Case 13, the number of stages between the condenser and the reboiler is doubled, but the holdup of each stage is reduced to 50% of the holdup per stage so as to maintain the total holdup of the column. The distribution and the overall reaction rate remain almost the same as Case 2. The amount of 48.8
Table
8. Continuation Homotopy 0.0000 0.0660
86.800 39.6524 36.8242 11.9454 3.4547
1251
2 3 4 5 6 7 8 9 10 II I2 13 14 Converged
0.1382 0.2 150 0.2881 0.3617 0.4358 0.5105 0.5860 0.6621 0.7391 0.8 170 0.8959 0.9757 I .0533 l.OOoo
path
for case
1
Euclidean 195.97 182.93 168.75 153.80 139.46 125.04 110.51 95.85 81.07 66.14 51.05 35.78 20.33 4.71 10.70 0.01022
Y.
1252
A. CHANGand J. D. SEADER
ethyl acetate in the distillate is increased by 10% from 10.08 to 11.03 lbmol h-l. NUMERICAL
ROBUSTNESS
AND
EFFICIENCY
Newton’s method and Newton’s method with a golden-section line search are economical techniques nonlinear equations. solving systems of for methods can be timeHomotopy-continuation consuming procedures, but are more reliable. When used properly, homotopy-continuation can reach a solution with a probability of one. In this study, Newton’s method and Newton’s method with a golden-section line search routinely failed with the most likely reasons for failure being inadequate starting guesses and the small acetic acid flowrates in the upper portion of the column. The temperature and reaction-rate profiles were found to have a great influence on liquid-phase composition. Thus, it was important to provide reasonable starting guesses. However, if the starting values were. located in improper regions, Newton’s method and Newton’s method with a golden-section line search failed. Due to the strong tendency of acetic acid to flow downwards, very little acetic acid exists in the vapor phase of the upper portion of the column. With Newton’s method, especially in early iterations, negative values of acetic acid vapor flowrates in the upper portion of the column were generated. When this occurred, the mapping Function of Naphtali and Sandholm (1971) was used to transform negative values to positive values. But Newton’s method and Newton’s method with a golden-section line search still diverged or stopped. Only homotopy-continuation was successful for all cases studied. As discussed by Waybum and Seader (1984), differental arclength homotopy continuation is not an efficient as Newton’s method. However, continuation can result in a converged solution when Newton’s method fails. A typical continuation path for this study is that of Case 1, which is given in Table 8. Each of the tist 13 steps consists of one Euler predictor and four Newton correctors. The Euclidian norm of the Functions [equations (l-3)] decreases from 195.97 at t = 0 to 10.70 at t = 1.0533 after the 14th predictor step. Seventeen Newton corrections are then used at t = 1.0 to reduce the norm From 0.86332 to an acceptable value of 0.01022. Thus, a total of 14 Euler predictors and 69 Newton correctors were used. The computations were performed on a UNIVAC 1170 and on a DEC MicroVAX II. On the MicroVAX, the executable module resulting from compilation and linking of the FORTRAN code consisted of 418 blocks (0.204 MBytes). The CPU time for Case 1, which was typical of all cases, was 20.6min on the MicroVAX. CONCLUSIONS
Modifications and extensions of the homotopy continuation program reported by Waybum and
Seader (1984) made that code suitable For studying reactive-distillation systems involving reversible, kinetically controlled reactions. In all cases, the calculations converged only when applying the homotopy-continuation method. For all the sets of specifications studied, a nearly pure ethyl acetate distillate was not obtained and a high conversion was not possible. The main problems were: (1) unfavorable total reaction conversion; (2) vapor-liquid equilibria; and (3) temperature profile. It was observed that K-values of ethanol and water were close, making separation of these two components very difficult. The volatility of acetic acid caused it to flow downwards and esterification occurred mainly in the lower part of the column. In the upper portion of the column, K-values of ethyl acetate were too close to ethanol and water, making the separation of ethyl acetate From ethanol and water in these stages very difficult. The temperature distribution in the column under different sets of specifications did not vary appreciably. Because the reaction conversion and K-values are strongly related to temperature, it was difficult to cause the temperature profile, reaction-rate profile, and K-values to change independently of each other. It is concluded that the application of reactive distillation to the esterification of ethanol with acetic acid is not technically Favorable in a two-product distillation column for the conditions imposed in this study. However, reactive distillation may be suitable for other reactive systems. In particular, the saponification reaction between propylene chlorohydrin and lime, as reported by Carra et al. (1979), to Form propylene oxide appears to be more suitable. The components involved have more Favorable volatility and reaction-rate characteristics than the system studied here. NOMENCLATURE A, B, C, ~1,b, c, a, y = Constants
A- =
y A,, A,._ . . . . A,,.” = AT, B,“, C,” =
C = Cj = C, = d = E, = Ep =
F, =
in the Bcnedict-WebbRubin equation of state The Arrhenius oreexoonential factor for the pth te& a coefficients for the Regression modified Mar&es equation Antoine constants of component j Number of components Molar heat capacity Concentration of the reference component for the qth reaction on stage i Molar density The discrepancy (residual) of the energy equation for stage i Energy of activation for the pth term
c f,,
,=I
f = Function fu = Liquid component molar flow for component j in the fresh feed to stage i
G, = i: g,, ,=I
Homotopy-continuation g, = Vapor component molar flow for component j in the fresh feed to stage i H = Homotopy function ZZ, = Enthalpy of the primary vapor leaving stage i; i.e. not including the vapor sidestream ZZG,= Enthalpy of Gi ffui = Molar enthalpy of the vapor leaving stage i IT;, = Virtual value of the partial molar enthalpy of component j in the vapor phase leaving stage i Hz. = Enthalpy of pure component j in the perfect gas state at temperature I: = Heat of formation of component j at AZ%, 298 K h, = Enthalpy of the primary liquid leaving stage i; i.e. not including the sidestream or entrained liquid h, = Liquid enthalpy of F, h, = Molar enthalpy of the liquid leaving stage i h;, ., = Virtual value of the nartial molar enthalpy of componen‘t j in the liquid phase leaving plate i K = Vapor-liquid equilibrium ratio Kti = Vapor-liquid equilibrium ratio for component j on stage i k, = Reaction coefficient for the pth term L, = i l,i ,=1 i, = Moles of component j in the primary liquid leaving stage i ; i.e. not including sidestream liquid or entrained liquid M,, = Discrepancy (residual) of the component material balance for component j and stage i N = Number of stages including condensers and reboilers NRXNS = Number of chemica1 reactions au = Exponent in chemical kinetics rate expression P = Pressure P, = Pressure on stage i P,” = Vapor pressure of pure component j Qi = Heat duty at stage i Q,_ = Discrepancy (residual) of the equilibrium relation for component j and stage i R = Reflux ratio; universal gas constant; molar reaction rate R,j = Initial guess for production rate due to chemical reaction of component j at stage i r, = dC,jdL, reaction rate for reaction n (with respect to the concentration of the reference component for the qth reaction) on stage i T = temperature; when used as a superscript denotes matrix transpose Tb = Normal boiling point temperature T, = Critical temperature T, = Temperature on stage i t = Time; temperature; homotopy parameter S, = Ratio of vapor drawoff to primary vapor (vapor not withdrawn or entrained on stage i s, = Ratio of liquid drawoff to primary liquid (liquid not withdrawn or entrained) on stage i un = Liquid holdup on stage i
for reactive distillation
1253
q = i:
=
cly
uij = Moles of component j in the primary vapor leaving stage i ; i.e. not including sidestream vapor nf = Liquid mol volume at stage i W, = F, + Gi x = Unknown x0 = Initial guess of unknown x, = Mol-fraction in liquid of component j x,/ = Mol-fraction in liquid of component j at stage i y, = Mol-fraction in vapor of component j vi, = Mel-fraction in vapor of component j at stage i Greek 2, = Molar latent heat of vaporization for component j at the temperature of stage i v,+,= Stoichiometric coefficient for component j in chemical reaction n R = Ri, - w;.. molar enthalpy departure function REFERENCES
Backhaus A. A., U.S. Patent 1400849 (1921). Benedict M., G. B. Webb and L. C. Rubin, An empirical equation for thermodynamic properties of light hydrocarbons and their mixtures. J. Chem. Phys. 10,747 (1942). Belck L. H., Continuous reactions in distillation equipment AZChE Ji. 1, 467 (1955). Bermans S., H. Isbenjian, A. Sedoff and D. F. Othmer, Esterification-continuous production of dibutyl phthalate in a distillation column. Znd. Engng Chem. 40, 2139 (1948). Carra S., M. Morbidelli, E. Santacesarla and G. Buzzi, Synthesis of propylene oxide from propylene chlorohydrins-II. Chem. Engng Sci. 34, 1133 (1979). Chan T. F., Deflation techniques and block-elimination algorithm for solving bordered singular systems. SIAM J. Sci. Srarist. Comput. 5, 121 (1984). Ferraris G. B., Interlinked, multistaged separations with nonstandard specifications solved by the NewtonRaphson method. AIChE JI. 27, 163 (1981). den Heijer C. and W. C. Rheinboldt, On steplength algorithm for a class of continuation methods. SIAM J. Namer. Analysis IS, 925 (1981). Izarraraz A., G. W. Bentzen. R. G. Anthony and C. D. Holland, Solve more distillation problems. l?Zydrocarbon Process 59, No. 6, 195 (1980). Jelinek J. and V. Hlavacek, Steady state countercurrent equilibrium stage separation with chemical reaction by relaxation method. Chem. Engng Common. 2, 79 (1976). Jing W., Q. Wang, W. Ling, X. Wang and Y. Yang, The simulation study of a reactor-three phases azeotropic distillation tower system. ZRSI. Chem. Engng Symp. Ser. 92, 675 (1985). Kaibel G., H. H. Mayer and B. Seid, Reaction in distillation columns. Ger. Chem. Engng 2, 180 (1979). Kinoshita M., I. Hashimoto and T. Takamatsu, A new simulation procedure for multicomponent distillation column processing non-ideal solutions or reactive solutions. J. Chem. Engng Japan 16, 370 (1983). Komatsu H., Application of the relaxation method for solving distillation. J. Chem. Engng Japan 10, 200 (1977). Komatsu H. and C. D. Holland, A new method of convergence for solving reacting distillation problems. J. Chem. Engng Japan 10, 292 (1977). Leyes C. E. and D. F. Othmer, Continuous esterification of
Y. A. CHANG
1254
and J. D. SEADER
butanol and acetic acid, kinetic and distillation consideration. Trans. AfChE 41, 157 (1945). Lin W., J. D. Seader and T. L. Waybum, Computing multiple solutions to systems of interlinked separation columns. AZChE JI. 33, 886 (1987). Naphatali L. M. and D. P. Sandholm, Multicomponent separation calculations by linearization. AZChE JI. 17,
where Z$=u;+n,
(A6)
Enthalpies for the ideal-gas state were computed
from:
AZChE
JI. II,
205 (1965).
Suzuki I., H. Komatsu and H. Hirata, Formulation and prediction of quaternary vapor-liquid equilibria accompanied by esterification. J. Chem. Engng JQPU~S3, 152 (1970). Suzuki I., H. Yagi, H. Komatsu and M. Hirata, Calculation of multicomponent distillation accompanied by a chemical reaction. J. Chem. Engng J~pm 4, 26 (1971). Wang J. C. and G. E. Henke, Tridiagonal matrix for distillation. ZZydrocarbon Process 45, No. 9, 169 (1966). Watson K. M., Thermodynamics of the liquid state, generalized prediction of properties. Znd. Engng Chem. 35, 398 (1943). Wayburn T. L., Modeling interlinked distillation columns by differential homotopy-continuation. Ph.D. Dissertation, University of Utah (1983). Wayburn T. L. and J. D. Seader, Solutions of systems columns by homotopyof interlinked distillation continuation methods. Proc. 2nd Znt. Conf. Foundation Computer-Aided Process Design (A. W. Westerberg and H. H. Chien, Eds). CACHE (1984). Xu X. E. and H. F. Chen, Simulation of distillation processes with reactions in series. Paper No. 5Of, 1986 Spring National Mtg, AZChE, New Orleans (1986). APPENDIX
A
H;
+ A+:
+ A&
+ ASx2xq + A6+xq 10xxx I
+A,x,x:+A
+
A,3~q~,~g
+ A,x,x:
Ci=ap+&T+y,Tz
C;dT,
(A7) (A8)
c
l-Tr2
AZ]=&
n
1
d
l-T,.,
(A9)
’
where the subscripts 1 and 2 refer to temperatures 1 and 2, and T, = T/T,. A common choice for n is 0.38 and for subscript 1 is the normal boiling point. Constants for equation (AS) are given by Izarrarax et Q/. (1980) in their Table 3.5. The Benedict et OZ. (1942) equation of state was selected to derive the enthalpy departure function: B,RT
P = RTd +
-A,,
- 2
d* + (bRT
(
-
a)d’
+ aad
>
+$[(l
(AlO)
+ydz)exp(-yd*)].
By introducing dimensionless groups into the BWR equation, the reduced equation can be written as: P,=T,d+d2x
I , (All) 1
T,B’-A’-$
A’ = A; + Q’d(l
- a’d”),
(Al2)
+ yd2)exp(-yd*),
(-414)
B’=B;,+b’d
(‘413)
C’ = C; -
c’d(l
P,=$
+ A,,x,x:
+ A,,x,x:
T, = f
c
c
(Al)
where x, is the liquid mol-fraction of component j, y, represents the activity coefficient, and A,-A,, are constants determined from quaternary equilibrium data by Suzuki et al. (1970) and given in their Table 3. The expression for the activity coefficient of the remaining components can be obtained by rotating the subscripts: I+ 2 + 3 * 4 4 1. The vapor pressure of pure liquid is calculated by using the Antoine equation: logP,o=A;--
298
using the expressions of Ixarraraz et ai. (1980) in their Tables 3 and 4. The latent heat of vaporixation decreases with increasing temperature and is zero at the critical point. A widely used correlation between A,, and T, is the Watson (1943) relation:
2 1+A II XXX 2 3 4+A 12XXX 3 . ,
+ A,,x&
+
where dimensionless groups are:
+ A,x,x, + A,x,x:
= AH&,,
where
K-Values Suzuki et al. (1970) rearranged the Margules equation as a polynomial series in mol-fractions of the components in the mixture: logy, = A,$
s r/
148 (1971).
Nelson P. A., Countercurrent equilibrium stage separation with reaction. AZChE JI. 17, 1043 (1971). Su G. J. and D. S. Viswanath, Part II. Generalized Benedict-Webb-Rubin equation of state for real gases.
(A5)
4) = R,, - rl,(Ti)
BO
A,=_
c;=_
Q’=-
c’=_
T+C,O
where the constants are given by Suzuki et al. (1970) in their Table 1. Vapor and liquid enthalpies The vapor and liquid enthalpies were calculated, respectively, by use of the virtual values of the partial molar enthalpies: c (A3) (A4)
CCIP, R=T4’c UP: R3T3’e CPZ R3TS’c
b’=’
AoPc R2T,2 bP f R2T2 c QEP,5 __ =R6Tz
”
Qa y’=-
YPf R2T2’c
Using a least-squares technique, a set of parameters for the generalized BWR equation was determined by Su and Viswamath (1965) who recommended the following generalixed constants: 47;= 0.04407059,
A; = 0.24180980, a’a’=0.11369478
x lO-3,
C; = 0.21217005,
b’=0.03715171, c’ = 0.06448001
B; = 0.07643101,
and
y’=O.O6.
The coel%ients for the BWR equation of state for the components of interest are given by Ixarraraz et al. (1980)
Homotopy-continuation
for reactive distillation
in their Table 3.6. In the mixture, each of the eight BWR constants are related to composition by mixing-rule equations of the form:
Values of /3 for each BWR constant are given Izarraraz et al. (1980) in their Table 3.7. From thermodynamics, the enthalpy departure function, 0, can then be obtained from the equation of state: R==
ew(-_yd*) 2
1255
+ yd*exp(-yd*)
. I
Kinetics The
reaction-rate expression for the reversible reaction: AcOH
+ EtOHuAcOEt
A
+ H,O
B
C D was obtained from Izarraraz ef ai. (1980): rA=k,C,C,-kkC,C,,
+ RT(& The resuit is: f-l=(B,,RT-2ZA0-~)d+(26R;-3a)dz 6aads t---+-3 5
ca2 T
I -
exp(--yd*) yd2
-
I>.
(A16)
(A17)
(‘4lW
where: kA = 29OOOexp(-7150/T),
6419)
k; = 7380exp(-7150/T),
(A20)
with kA and k_A in 1 gmol-’ min-’ and Tin “K. For comparison,, Suzuki ef al. (1970) used just the forward rate portion of equation (A18) with log,,k,,
=
-2710/T
+ 3.70.
6421)
009% I354/88 S3.00 + 0.00 Copyright6 1988 PergamonPressplc
SIMULATION OF CONTINUOUS REACTIVE BY A HOMOTOPY-CONTINUATION
DISTILLATION METHOD
Y. A. CHANG and J. D. SEADER~ Department of Chemical Engineering, University of Utah, Salt Lake City, UT 84112, U.S.A. (Received 25 March 1988; final revision received 5 May 1988; received for publication 19 May 1988) Abstract-Reactive distillation was studied theoretically for the esterification of acetic acid with ethanol to produce ethyl acetate and water, using second-order reversible kinetics. Because of the high degree of interaction among the iteration variables, convergence was difficult-to-impossible to achieve with Newton’s method. Consequently, a robust homotopy-continuation method was applied to solve the simultaneous nonlinear equations used to model the reactive-distillation system.Case studiesfor several different specifications were conducted to show how design variables affect the reactive-distillation system.
INTRODUCllON
Because of the potential importance of reactive distillation, a number of publications, both patents and articles, have considered the problems of its application and performance. Prior to the advent of digital computers, the literature dealt mainly with proposed applications and simplified calculational procedures. Backhaus (1921), in perhaps the first reference to reactive distillation, published a patent describing an esterification process carried out in a distillation column. Leyes and Othmer (1945), as well as Bermans et al. (1948), treated an esterification reaction in a distillation column and utilized experimental vapor-liquid equilibrium data and kinetic and chemical equilibrium relations to get numerical solutions. Belck (1955) made assumptions to simplify calculations for hypothetical reversible reactions carried out in a single distillation column, Since 1970, calculation procedures suitable for implementation on digital computers, for rigorously solving the mass, equilibrium, and energy balances in multicomponent, multistage reactive-distillation columns have been reported. The procedures, which involve the solution of large sets of combined nonlinear and linear equations, as well as applications, are summarized in Table 1. From a numerical calculation procedure, these previously reported studies can be divided into three categories: (1) methods using tear variables (bubble-point method); (2) relaxation techniques (dynamic approach); and (3) methods incorporating the Newton-Raphson method. For the tearing method, it is not necessary to compute derivatives, and the tridiagonal matrix algorithm can be implemented to decrease computer storage space and increase computation efficiency. However, when differences in boiling point between components is too large, when the kinetics are comtAuthor to whom all correspondence should be addressed.
plex, or when the liquid solutions are highly nonideal, the tearing method is very difficult to converge. The relaxation method has better stability, but may require many more iterations to obtain a solution. Although the Newton-Raphson method converges quickly from a suitable starting guess and can be applied more widely to different specifications, it needs large storage and involves complicated derivative operations. Convergence is not guaranteed for any of the procedures used in previous studies. Table 1 shows that the following reaction processes have been reported in the literature: esterification, saponification and the manufacture of propylene oxide from chlorohydrin. The esterification of acetic acid with ethanol is the most widely studied case. Not included in Table 1 are the required specifications, which are essentially the same for all references listed. The user must specify the number of stages, feedplate location, feed condition, reflux ratio, and distillate or bottoms flowrate. In this study, it was desired to develop and demonstrate a robust numerical procedure for reactive distillation, involving the use of homotopycontinuation in the manner of Wayburn and Seader (1984). The procedure was designed to permit great flexibility with respect to column configuration (including interlinked systems), specifications and choice of thermodynamic properties.
MATHEMATICAL
MODEL
Consider an N-stage column separating C components, as shown in Fig. 1, where stage 1 is a total condenser that produces a saturated liquid and stage N is a partial reboiler. Furthermore, let chemical reactions in the liquid phase occur at each stage. Assume that: (1) the process has reached a continuous steady state; (2) each stage is a perfectly mixed stirred-tank reactor (thus, the liquid composition at each stage is homogeneous and equal to
1243
Carra er u/. (1979)
Combine tearing method (choose the profile of liquid mol fraction of [OH-] as tearing variables) and Newton-Raphson method
rithm
Reaction type not stated
(i) propylene chlorohydrin t 0.5 calcium hydroxide = propylene oxide t 0.5 calcium chloride (ii) propylene oxide _t water = propandiol Vapor pressure, Antoine equation, activity coefficient by van Laar equation, salting Evaluate new profik of liq uid mol fraction of [OH-]
Not stated
effect
(i) AtB=CtD (ii) normal alcohol t carboxylic acid + eater t water (iii) saponification, acetate t steam = alcohol t acetic acid Reaction type not stated Not stated
Combine Newton-Raphson method and banded structure of Jacobian matrix algo-
AcOH t EtOH = AcOEt t H,O Second-order Reversible
K is a function of temperature and liquid composition
vergence
Kaibel er al. (1979)
method
Multi 8-n method of con-
Newton-Raphson
Komatsu and Holland (1977)
position
Second-order Reversible
Reaction rate constant Arrhenius equation
(i) AcOH t EtOH = AcOEt t H,O (ii) AcryOH t EtOH = AcryOEt t Hz0 Second-order Reversible
AcOH + EtOH = AcOEt t H,O
Relaxation factor and corm section for liquid com-
Relaxation technique
Komatsu (1977)
Vapor pressure, Antoine equation, activity coefficient by NRTL equation
r =~,~~l~~~lI~l-~,M~~l~
Reaction-rate quation
A+B=C+D
AcOH t EtOH = AcOEt t H,O Second-order Irreversible
Reaction and its type
K is a function of temperature and liquid
Relaxation factor
Relaxation technique
Jelinek and Hlavacck (1976)
Not stated clearly
equation
Case 2: vapor pressure, Antoine equation, activity coetiicient by modified Margules
Case I: K is a function of temperature only
Thermodynamic correlation for vapor-liquid equilibrium
systems
composition
Damping factor
method
NewtonRaphson
Nelson (1971)
Convergence method Normalize the liquid cornposition and use modified Mullet’s method for convergence of temperature
procedure
Choose temperature as tear variable and use the tridiagonal matrix algorithm for the solution of the linearized material balance equations
Calculation
Suzuki er al. (1971)
Reference
Table I. Rigorous methods for computing reactive distillation
Homotopy-continuation
for reactive distillation
1245
the composition of the liquid leaving the stageCSTR assumption); and (3) the vapor and liquid leaving any stage are in physical equilibrium. The extent of liquid-phase reaction at each stage is governed by reversible and/or irreversible kinetic rate expressions. Figure 2 shows a schematic representation of a general reactive-distillation stage. The three types of functions which describe the physical and chemical processes on stage i are: 1. Component
material balances:
M,=(l
-ts,)L),,-t(l -f;j-gij-
+si)fij-ui,,J-~l,_,J 2iRYN.s ui c v,mr’in “-=I
j=l,...,C
(1)
where the last term accounts for chemical reactions. 2. Energy balances: Ei = (1 + S,)H,
+ (1 + si)hi -ha-,
Hi, L -
b-j -
&,
-
Qi,
(2)
where unlike the component material balances, it is not necessary to take the heat of reaction into account because enthalpies are referred to the elements rather than to the components. 3. Equilibrium relationships: Q, = K,, iii 2, -
vy ,
j-1
9--m, C
(3)
where Q,J is derived from the definition of Key: K,j = ti x,j
or
Kijx,
-
y, = 0,
(4)
and x and y represent component mol fractions. Specifications for equations (l-3) include the total number of stages, stage locations of all feeds, sidestreams, and heat exchangers; all stage pressures and liquid-phase holdup volumes; and complete specification of each fresh feed. This leaves N specifications that must be made before solving the N(2C + 1) equations for the corresponding N(2C + 1) iteration variables. If the N specifications are N stage heattransfer rates, then the equations are solved for 2CN component molar flows and N stage temperatures. However, because condenser and reboiler duties are strongly dependent and can generally be specified independently over only narrow unknown ranges, other specifications are more desirable. The total condenser stage is a special case which is not represented by equation (3). Let the component molar flows of the liquid distillate be called If the reflux ratio is specified, the “I.1 3U1\2, . . . 3 h.c-
Y. A. CHANG
1246
and J. D. SFZADER
Fig. 1. Configuration for reactive-distillation column.
top-stage
energy
balance
equation
is changed to a
For the case of all r, = 0, Naphtali and Sandholm (1971) showed that ordering the above equations as follows for each stage: Ei, M,,, . . . , M ,,=, and differentiating the equations in the Q,,, ...qQt.,; order: Ti_ ,,I, -,,, ,.‘.,~i-,,c,~i.,,..-,~i,c, Ti;.,4.L,.-., T.+, leads to a block-tridiagonal 6.,90,+1.19....~i+1.e, Jacobian matrix with minimum bandwidth. In the procedure applied here, the block-reduction scheme of Waybum (1983) was used to deal with the matrix, making the problem of bandwidth irrelevant. To avoid unnecessary pivoting, Waybum placed the derivative of the energy balance equation with respect to temperature on the diagonal. The equations were ordered M,,, , . . . , M,,,, E,, Q;,, , . . . , Q,,, and then differentiated in the order of Naphtali and Sandholm (1971). To achieve flexibility with respect to specifications, the method of Ferraris (1981), as implemented by Waybum and Seader (1984) was used. This involved adding specification equations as rows at the bottom of the standard set and adding the additional unknowns as columns at the right side of the Jacobian matrix. The total condenser, modeled by the reflux-ratio equation (5), the bubble-point equation (6) and mol-fraction identities (7) were considered as standard equations. Chemical-reaction rates, rin. in the liquid phase were modeled by reversible power-law kinetic rate equations. Vapor molar holdup and vapor-phase chemical reactions were neglected. Because the material balance equations (1) are in the form outminus-in, the term represents the moles of component j produced in the liquid phase. Therefore, v,,,, the stoichiometric coefficient for component j in the n th chemical reaction, is positive if component j is created by the nth reaction. The reaction-rate equations were written in terms of concentrations of the relevant components:
reflux-ratio equation: E,=L,-RV,,
(5)
where R is the reflux ratio. The distillate composition, D1.1r”,,2, . , qE, is the same composition as in the liquid leaving the condenser and, therefore, the equilibrium relations are not independent. the equilibrium equation for component 1 is replaced by the bubblepoint equation:
s iY
ij
t
Pi
Q,., = 2 Wi,p4.p)~ - vl.p* I p-1
(6)
where component 1 must be present in the distillate in nontrivial amounts. Like the energy balance equation, the other equilibrium equations, Q,*-Q1=, reflux are substituted for by the component equations; viz.
I i-l./
' ij
Ti
7
'ij
Y.
r+rj
l ij
Sij
Fig. 2. General reactive-distillation stage.
1 ij
Homotopy-continuation for reactive distillation
1247
of temperature only. By making initial guesses for the R, values, the modified Wang-Henke linear equations could be readily solved for initial composition profiles. SIMULATION
=xA,exp P
(8)
where C,= i,/(Livk) is the concentration of component q. The variable A, is the Arrhenius preexponential factor, E, is the energy of activation, R is the universal gas constant and 7’, is the temperature at stage i. This reaction model can deal with either irreversible or reversible reactions. For example, for the reversible reaction: A+B=C+D
(9) where p = 1 indicates the forward reaction, p = 2 is the reverse reaction, /cl is positive and kz is negative. To achieve robustness in convergence, the modeling equations, call them f(x) = 0, were embedded into a Newton homotopy: H(x, 1) =f(x)
- (1 - ~)f(x’)
= 0,
(10)
where t is a scalar homotopy parameter, which is varied from 0 to 1 as continuation is used to track the path from the starting guess x0 (at t = 0) to the solution at t = 1. Tracking was by means of differential arclength continuation using an Euler predictor and a Newton corrector as described by Waybum and Seader (1984). Initially, a maximum step size equivalent to moving from t = 0 to 2 = 1 was attempted. This is Newton’s method, which was iterated with and without a line search. If it failed to converge, for reasons discussed later, continuation was employed using the step-size algorithm of den Heijer and Rheinboldt (198 1), which was used by Lin, Seader, and Waybum (1987). If turning points occurred along the path, they were tracked by the method of Chan (1984), as described by Lin et al. (1987). For providing the starting guesses for the iteration variables (unknowns), a modification of the first step in the Wang and Henke (1966) method was used. For each stage, equation (1) is written in mol-fraction form as: Mj=(l
+
si)KYt,+ -
C1
Li_,,ixi_
+si)Llxij-
,.i-
C+l.jYt+l.j
Fjzii -
R,=
0.
(11)
By combining equations (4) and (11) and selecting T, and Vi as tear variables, the resulting NC equations
are linear in the unknown liquid mol-fractions provided that K values are composition independent. Since this was not the case for the system studied here, the method of Suzuki er al. (1971) was used to estimate starting guesses for K-values as a function
.STUDIES
All the simulation case studies reported here were made on the four-component reactive system involving acetic acid (l), ethanol (2), water (3) and ethyl acetate (4). with the sole reaction being: acetic acid + ethanol - ethyl acetate + water. Since at least 1921, this system has been considered to be suitable for reactive distillation. Because of the non-ideal nature of the liquid-phase mixture caused by the polar molecules, it was necessary to give careful attention to the thermodynamic correlations used to compute vapor-liquid equilibrium. In particular, acetic acid tends to form a dimer and a trimer in the vapor phase and it is well-known that binary mixtures of ethyl acetate-water, acetic acid-water, and ethyl alcohol-water, as well as the ternary mixture of ethyl acetate-acetic acid-water, all form azeotropes. K-values were computed from Kj = yj Pg/ P where yj was computed from a modified Margules equation developed especially by Suzuki et al. (1970) for the system studied here. Suzuki et al. (1970) also present equations for correcting the vapor-phase composition for the dimer and trimer of acetic acid as well as the constants for each component in the Antoine vapor pressure expression. For vapor and liquid enthalpy, the correlations presented by Izarraraz et al. (1980), which are referenced to the elements, were used. Details of the thermodynamic and kinetic equations, together with references to listings of the constants in the equations, are given in the Appendix. The first simulation, Case 1, was made using the same specifications as Suzuki et al. (1971) and Izarraraz et al. (1980), as shown in Fig. 3. There are 13 theoretical stages, all at 1 atm pressure, counting the partial reboiler and the total condenser. The feed rate is 107.6lbmol h-‘, with a liquid distillate of 20.8 lbmol h-‘, giving a bottoms of 86.8 lbmol h-l. The feed is preheated to its bubble point at feed tray pressure and the reflux ratio is 10. The feed is fed to the sixth stage from the top and has the following mol-fraction composition: acetic acid z, = 0.4962, ethanol z2 = 0.4808, water z, = 0.0229 and ethyl acetate z, = 0.0. Holdup volumes are 266.95 and 80.085 ft3, respectively, for the reboiler and each of the stages l-l 2. Although these specifications are the same as those of Suzuki et al. (1971) and Izarraraz et ai. (1980), the results of this study as discussed below, were different, but much closer in agreement to Suzuki et al. (1971). Comparisons with Izarraraz er al. (1980) are not included here because they used considerably different vapor-liquid equilibrium ratios as compared to those used here and by Suzuki et al. (1971).
Y. A. CHANG
1248
P=lotm
Retlfm Ratio-10.0 Hddups : Reboller=2&5.9Bcuf Stage 1-12 [email protected] /each 5
Saturated
Uquld
,brn0l AcOH EtOH W&m AcoEi
9
I
I 6
53Ao 51.74 2.46 0.00 107.60
12
E+ 13
1
Bottoms
Fig.
3.
B&B LbmoC h-’
Specifications for Case
1
Table 2 compares the computed temperature, liquid flow and vapor flow profiles with those of Suzuki et al. (1971). For stages above the feed, the profiles computed in this study agree well with those of Suzuki et al. However, as the reboiler stage is approached, the profiles begin to deviate. The deviations are caused not by differences in K-values but by differences in the reaction-rate expressions. As shown in Table 3, K-values for each component at each stage are almost identical based on liquid compositions from Suzuki ef al. (1971) who did not allow for the backward reaction, while this study did. As seen in Table 2, the resulting reaction-rate profiles are quite different. Small negative production rates of ethyl acetate were computed here for the upper four stages, while Suzuki et al. show zero rates there. For most of the lower nine stages, the reaction rate is greater than that of Suzuki et al. due to the use of Table
2. Comparison
of Case
and J. D. SEADER
different kinetic constants for the forward reaction as obtained from Izarraraz et al. (1980). For Case 1, the extent of reaction and separation achieved are quite unfavorable. The extent of reaction, based on the limiting reactant, is only 28.3% and the distillate contains appreciable ethanol. Accordingly, attention was directed toward changes in column configuration and operating conditions which might increase conversion and improve the separation. The most common industrial schemes are to use an excess of one reactant (usually the least expensive one) and to recycle the unconverted material. If the recycle configuration proves to be uneconomical, it is possible to use a countercurrently operated sequence of continuous stirred-tank reactors or to remove continuously one or more of the products from the reaction zone. In such devices, the reaction is thus driven toward completion. In some situations, a staged distillation column is a very efficient reactor. By use of the technology of reactive distillation, it is possible to get high reaction conversion and separation simultaneously. For example, consider a reversible reaction of the type, A + B = C + D, where B and D are more volatile than A and C. Liquid rich in A is fed to a stage between the top and the middle portion of the column. Vapor rich in B is fed to a stage between the middle portion and the bottom of the column. Reactant B is absorbed into the liquid phase where the reaction takes place and product D is stripped from the liquid phase and carried out the top of the column. In the case of esterification, the order of volatility is ethyl acetate, ethanol, water and acetic acid. Thus, it is expected that the middle portion of the distillation tower is the chief reaction zone. The rectifying section fractionates the ethyl acetate out of the acetic acid, and the stripping section removes alcohol from water. Ideally, the ethyl acetate is the distillate and water is the bottoms product. However, over a wide range of composition, as shown in Table 3, ethanol and water do not differ greatly in volatility, making it difficult to move ethanol up the column. Case 2 was designed to gain a better understanding of reactive distillation for the acetic acid-ethanol system. The only difference from the specifications of I results with Suzuki
et al. (1971) Suzuki
This study stage No. 1 2 3 4 5 6 7 8 9 10 I1 12 13
(Ibm&‘) 162.8 163.1 163.5 164.2 166.3 176.1 176.5 177.0 177.6 178.3 180.0 185.4 199.1
208.0 206.6 204.5 201.7 198.5 303.7 302.4 300.9 299.3 297.8 297.2 303.3 86.8
(lbmCL) 20.8 228.8 227.4 225.3 222.5 219.3 216.9 215.6 214.1 212.5 211.0 210.4 216.5
Reaction rate (lbmol h- ‘) -0.019 -0.017 -0.016 -0.013 0.05 0.56 0.60 0.65 0.72 0.80 0.96 I .43 8.03
(a,) 163.0 163.3 163.7 164.4 166.5 176.2 176.6 177.1 177.7 178.4 180.0 184.9 197.9
(lbm:
h-‘)
208.0 206.6 204.6 201.8 198.7 303.8 303.0 302.1 301.3 300.5 300.6 306.5 86.8
er al. (1971) (lbm:h-‘) 20.8 228.8 227.4 225.4 222.6 219.5 217.0 216.2 215.3 214.5 213.7 213.8 219.7
Reaction rate (Ibmol hW’) 0.00 0.00 0.00 0.00 0.00 0 27 0.29 0.30 0.32 0.35 0.45 0.94 11.66
Homotopy-continuation for reactive distillation Table Acetic
10 II 12 13 Case Case
Table
No. 1
EtOH
-
-
f 8 1: 11 12 13
11,12-i=
0.9120 0.8965 0.8865 0.9061 1.0797 1.0813 .085 1 1.0918 1.1042 1.1357 1.2383 1.5011
I
1
2
acetate 2
053
0.93 0.9 1 0.88 0.89
0.92 0.9 1 0.88 0.89
1.12 1.06 0.91 0.99
t.11 .05 0.91 0.98
1.05 1.09 1.15 .24
0.05 0.13 0.13 0.13 0.14 0.14
0.05 0.13 0.13 0.13 0.14 0.14
0.90 1.07 I .08 .oa
0.90 .07 1.08 1.08 1.09 1.10
0.88 1.07 1.07 1.07 1.07 1.07
0.88 1 .Ol 1.07 1 .O7 1.07 1.07
1.42 1 .I32 1.86 1.92 1.98
1.05 1.09 1.15 1.24 1.42 1 .I32 1.86 1.92 1.98
0.14 0.20 0.35
0.15 0.20 0.35
1.13 1.23 1.48
1.13
1.08 1.14 1.34
1.08 1.14 I.35
2.06 2.16 2.29 2.45
2.06 2.16 2.28 2.44
eauilibrium
AcOH
Ethyl
I
t I.09
1.10
t
1.23 1.48
I
1 = Suzuki er al. (1971). 2 = This study.
Vapo-liauid
0.036 1 0.0363 0.0380 0.0491 0.1313 0.1322 0.1339 0.1366 0.1414 0.1551 0.2077 0.3650
2 3 4 5
CACE.
4.
(K-values)
water
0.04 0.04 0.04
Case 1 was the feed composition, which was changed to an equimolar bubble-point mixture of acetic acid and ethanol. It is noted above that ethanol and ethyl acetate are more volatile than acetic acid and water. Thus, the reaction should occur mainly in the middle portion of the distillation tower. By examining the vapor-liquid equilibrium data shown in Table 4, the volatility of acetic acid appears to be suitable. It is by far the least volatile component in the system, with K values much smaller than for the other three components. Unfortunately, the ratios of KEIOH to KHlo vary only from 0.87 to 1.06, which are too close to 1. Ethanol is not much more volatile than water and, in fact, less volatile in the upper stages of the column. The results of Case 2 are listed in Table 5. The temperature profile of the column is far from linear, with two temperature jumps evident. One occurs at the feed point between stages 5 and 6 and the other appears at the reboiler. The change in temperature for the top five stages is only 3.S”F. The change in temperature for stages 6-l 1 is just 0.8”F. Between the 12th stage and the reboiler, the temperature difference is 13.8”F. Since esterification is an endothermic reaction, equilibrium shifts at higher temperatures toward the products. As seen in Table 5, the column temperature profile causes the chemical reaction to take place in the lower portion instead of the middle part of the distillation tower. The volatility of acetic acid coincides with the requirement that ideally it should be less volatile than
Sync
ratios
2
-
9
equilibrium
Ethanol
2
1
I 2 3 4
8
Vapor-liquid
1
stage
5 6 7
3.
acid
1249
4
data
for
Case
2
Water
AcOEt
1.0439 0.9691 0.9021 0.8823 .0730 1.0754 1.0777 1.0794 1.0814 1.0961 1.1795 1.4135
1.1001 1.1590 I .2496 I.4161 1.8123 1.8598 1.9178 1.9871 2.0697 2.1720 2.3012 2~4652
-
I
P 8
B
5 = t
B 3
;;‘
I
Y. A. CHANG and J. D. SEADER
1250
Table
Case
6. Specification
Differences from Case k
compare with Casek
3
2
Two
liquid
4
3
5
4
One Liquid feed one vapor feed Feed locations
6
4
Feed locations
7
2
Feed compositions
8
2
Feed compositions
feds
9
6
Feed compositions
10
6
Feed compositions
II
9
Holdup
12
9
Distillate
13
2
Number Holdup
rate of stages
ethanol and ethyl acetate. But the K values of acetic acid are too small compared to those of the other components. Most acetic acid tends to flow downwards and appears in the lower portion of the column only; very little exists in the upper and middie parts. In the middle portion, the liquid mol-fractions of acetic acid are not more than 0.2. Therefore, the production rate of ethyl acetate is limited by the deficiency of the reactant acetic acid. The liquid composition of ethanol increases from the condenser to a maximum at the 5th stage. Then it decreases from the 5th stage to the reboiler. The distribution of ethanol is favored for the reaction taking place in the middle portion of the column. Since there is a strong interaction between the reaction and separation, each specification of the system exerts a complex effect that is difficult to predict. Accordingly, as listed in Table 6, a number of variations to the specifications of Case 2 were studied in an attempt to achieve a higher conversion and a better separation. Resulting product distributions for all cases are given in Table 7. In Case 3, two separate saturated liquid feeds are used. Acetic acid is fed to stage 6 and ethanol is fed to stage 12. The separation and the reaction rate do not improve much. The big difference between Case 2 and Case 3 is the compositions of ethanol and water in the distillate. Of 20.8 lbmol h-l distillate, there are 9.86293 lbmol h-’ ethanol and 0.85561 lbmol h-i water in Case 2, while in Case 3 there are 6.62495 lbmol h-’ ethanol and 3.11422 lbmol h-r water. However, the total amount of ethanol and water just changes slightly. In the ideal case, it is preferred to have two feeds, one liquid and the other vapor. In Case 4, saturated liquid acetic acid is fed to stage 6 and saturated vapor enters at stage 11. Compared with Case 3, the two-
for all -
Specification
changes
Liquid AcOH -6th stage Liquid EtOH + 12th stage Liquid AcOH +6th stage Vapor EtOH - 11 th stage Liquid AcOH + 5th stage Vapor EtOH - 9th stage Liquid AcOH -4th stage Vapor EtOH - 10th stage 58.8 lbmol h- ’ AcOH 48.8 lbmol h-’ EtOH 48.8 lbmol h-’ AcOH 58.8 lbmol h-’ EtOH 58.8 lbmol h ’ AcOH 48.8 lbmol h - ’ EtOH 48.8 lbmol h ’ AcOH 58.8 lbmol h-’ EtOH stages l-l 3 = 94.459 fts/cach stage Distillate = 16.00 lbmol h-’ Reflw. ratio = 13.00 Number of stages = 24 Holdups: stage 1 80.085 ft3 stage 2-23 40.0425 ft’/each stage 24 266.95 ft”
stage
liquid-feed case, the reaction rate in the middle portion of the column is a little greater than for Case 3, but the separation is not improved. In Cases 5 and 6, the two feed stages are moved upward from Case 4. Comparing these three cases in Table 7, the production rate of ethyl acetate in Case 6 is about 10% more than that of the other two cases. As mentioned earlier, acetic acid has a strong tendency to flow downwards. If the feed stage for acetic acid is moved upwards, mol-fractions of acetic acid in the upper portion should increase, causing the forward reaction of esterification to increase. In the upper part of the column, either the negative values of production rate of ethyl acetate decrease or the positive values increase. Therefore, the total production rate of ethyl acetate increases but not significantly. One-feed cases, Cases 7 and 8, were made for comparing the effect of feed composition changes with Case 2. Two-feed cases, Cases 9 and 10, were chosen to compare with Case 6. In Table 6, feeds in Cases 7 and 9 consist of an excess of acetic acid. Flowrates of feeds in Cases 8 and 10 consist of an excess of ethanol. Conclusions are the same for the one-feed and the two-feed cases. The greater the acetic acid-to-ethanol ratio, the more ethyl acetate is created and vice versa. Since most of the reaction takes place at the higher temperature in the lower portion of the column, an increase in the amount of acetic acid causes the reaction rate to increase. Because volatilities of ethanol and acetic acid differ, ethanol is distributed more evenly than acetic acid. In the main reaction zone, the increased concentration of ethanol does not compensate for the loss of acetic acid. From these cases, distillate compositions do not change significantly with feed composition change. When the feed flowrate of ethanol was changed from
Homotopycootinuation Table Component
AcOH EtOH Water AcOEt
7. Final
product
Feed
53.4000 51.7400 2.4600 0.0000 107.600
AcOH EtOH Water AcOEt
53.8000 53.8000 O.lYJOCt O.OOC+l 107.600
AcOH EtOH Water AcOEt
53.8000 53.8000 0.0000 0.0000 107.600
AcOH EtOH Water AcOEt
53.8ooO 53.8IXtCI 0.0000 o.oooo 107.6OfJ
AcOH EtOH Water AcOEt
53.8000 53.8000 o.oooo 0.0000 107.600
AcOH EtOH Water AcOEt
53.8000 53.8000 0.0000 O.OODO to7600
AcOH EtOH Water AcOEt
58.8000 48.8000 0.0000 0.0000 107.6000
AcOH EtOH Water AcOEt
48.8000 58.8Oaa 0.0000 0.0000 107.6LM
AcOH EtOH Water AcOEt
58.8000 48.8000 O.OOOfJ O.OOOCl 107.600
AcOH EtOH Water AcOEt
48.8000 58.8000 0.0000 0.0000 107.600
AcOH EtOH Water AcOEt
distribution Distillate (Ibmol
case 1 0.0019 9.1843 1.3121 IO.3018
13 cases
BOttOll-lS
RePcti0tl
h-‘)
39.6524 28.8100 14.8937 3.4440
20.800
86.800
Case 2 0.0011 9.8629 0.8556 10.0804
40.2491 30.3873 12.6942 3.4694
20.800
86.800
Case 3 0.0064 6.6250 3.1142 I 1.0544
39.8220 33.2034 10.8574 2.9172
20.800 Case 4 o.On5 1 7.1773 2.6834 10.9343
39.3644 32.1921 11.7472 3.4%3
20.800
86.800 39.0127 30.9404 12.7823 3.4440
20.800
86.800
Case 6 0.0057 7.2238 2.4363 11.1342
38.8323 31.6141 12.5258 3.8279
20.800
86.800
Case 7 0.0014 9.2747 1.0383 10.4856 20.800 Case 8 0.0009 10.5905 0.7217 9.4869 20.800 Case 9 0.0067 6.7776 2.6414 11.3744 20.800 Case to 0.0050 7.7564 2.2740 10.7647 20.800
20.800
- 13.746 - 13.746 13.746 13.746
- l3SM - 13.550 13.550 13.550
- 13.972 - 13.972 13.912 13.972
86.800
Case 5 0.0033 8.0456 2.0017 10.3018
case 11 58.8OOoO 0.0081 48.80000 6.8466 O.OOODO 2.8459 O.OOOOO 11.0994 107.600
for
for reactive distillation
44.5078 25.2345 t3.2525 3.8052
- 14.431 - 14.431 14.431 14.431
- 14.784 - 14.784 t4.784 14.784
- 14.962 - 14.962 14.962 14.962
- 14.291 - 14.291 14.291 14.291
86.800 36.2464 35.6568 11.8310 3.0658
- 12.553 --12.553 12.553 12.553
86.800 43.285 1 26.5143 12.8668 4.1338
- 15.508 - 15.508 15.508 15.54!8
AcOH EtOH Water AcOEt
58.8000 48.8ooO 0.0000 0.0000 107.600
AcOH EtOH Water AcOEt
53.8000 53.8000 0.0000 0.0000 107.cwoct
- 14.219 - 14.219 14.219 14.219
86.800 45.4601 28.6216 to.4859 2.2324 86.800
- 13.332 - 13.332 13.332 13.332
Case 12 0.0051 5.0929 1.9198 8.9823
44.8661 29.7783 12.0090 4.9466
16.000
91.600
Case I3 0.0005 9.3724 0.4013 1 t .0259
39.9264 30.554s 13.4718 2.8473
20.800
- 13.929 - 13.929 13.929 13.929
- 13.873 - 13.873 13.873 13.873
86.800
to 58.8 lbmol h-l, ethanol appeared in the distillate in almost the same amount, with the excess appearing in the bottoms product. In Cases I-10, the volume holdup ratio for the reboiler to each stage was 10 to 3. Approximately half of the ethyl acetate was produced in the reboifer. In Case 11, the holdup of each stage is specified to be equal to that of the reboiler (94.459 ft’), with the total holdup of the column unchanged. Because the temperature and the concentrations of components (especially acetic acid) in the reboiler is advantageous for the forward reaction, the total rate decreases. Due to the change of holdup, the reaction is distributed more uniformly over the column. The separation is slightly improved, with a 13% lower production rate of 13.332 Ibmol h-’ of ethyl acetate. In Cases l-11, the overall reaction rate is computed to be approx. 14 lbmol h-‘. With a distillate specification of 20.8 lbmol h-‘, the purity of the resulting ethyl acetate distillate is poor. In Case 12, the distillate is reduced to 16.0 lbmol h-’ in an effort to increase purity. At the same time, the reflux ratio is changed to 13.0 to keep the vapor and liquid traffic unchanged. These changes improve the purity of ethyl acetate in the distillate by less than 5% even though the overall reaction rate is maintained at 141bmol h-l. In Case 13, the number of stages between the condenser and the reboiler is doubled, but the holdup of each stage is reduced to 50% of the holdup per stage so as to maintain the total holdup of the column. The distribution and the overall reaction rate remain almost the same as Case 2. The amount of 48.8
Table
8. Continuation Homotopy 0.0000 0.0660
86.800 39.6524 36.8242 11.9454 3.4547
1251
2 3 4 5 6 7 8 9 10 II I2 13 14 Converged
0.1382 0.2 150 0.2881 0.3617 0.4358 0.5105 0.5860 0.6621 0.7391 0.8 170 0.8959 0.9757 I .0533 l.OOoo
path
for case
1
Euclidean 195.97 182.93 168.75 153.80 139.46 125.04 110.51 95.85 81.07 66.14 51.05 35.78 20.33 4.71 10.70 0.01022
Y.
1252
A. CHANGand J. D. SEADER
ethyl acetate in the distillate is increased by 10% from 10.08 to 11.03 lbmol h-l. NUMERICAL
ROBUSTNESS
AND
EFFICIENCY
Newton’s method and Newton’s method with a golden-section line search are economical techniques nonlinear equations. solving systems of for methods can be timeHomotopy-continuation consuming procedures, but are more reliable. When used properly, homotopy-continuation can reach a solution with a probability of one. In this study, Newton’s method and Newton’s method with a golden-section line search routinely failed with the most likely reasons for failure being inadequate starting guesses and the small acetic acid flowrates in the upper portion of the column. The temperature and reaction-rate profiles were found to have a great influence on liquid-phase composition. Thus, it was important to provide reasonable starting guesses. However, if the starting values were. located in improper regions, Newton’s method and Newton’s method with a golden-section line search failed. Due to the strong tendency of acetic acid to flow downwards, very little acetic acid exists in the vapor phase of the upper portion of the column. With Newton’s method, especially in early iterations, negative values of acetic acid vapor flowrates in the upper portion of the column were generated. When this occurred, the mapping Function of Naphtali and Sandholm (1971) was used to transform negative values to positive values. But Newton’s method and Newton’s method with a golden-section line search still diverged or stopped. Only homotopy-continuation was successful for all cases studied. As discussed by Waybum and Seader (1984), differental arclength homotopy continuation is not an efficient as Newton’s method. However, continuation can result in a converged solution when Newton’s method fails. A typical continuation path for this study is that of Case 1, which is given in Table 8. Each of the tist 13 steps consists of one Euler predictor and four Newton correctors. The Euclidian norm of the Functions [equations (l-3)] decreases from 195.97 at t = 0 to 10.70 at t = 1.0533 after the 14th predictor step. Seventeen Newton corrections are then used at t = 1.0 to reduce the norm From 0.86332 to an acceptable value of 0.01022. Thus, a total of 14 Euler predictors and 69 Newton correctors were used. The computations were performed on a UNIVAC 1170 and on a DEC MicroVAX II. On the MicroVAX, the executable module resulting from compilation and linking of the FORTRAN code consisted of 418 blocks (0.204 MBytes). The CPU time for Case 1, which was typical of all cases, was 20.6min on the MicroVAX. CONCLUSIONS
Modifications and extensions of the homotopy continuation program reported by Waybum and
Seader (1984) made that code suitable For studying reactive-distillation systems involving reversible, kinetically controlled reactions. In all cases, the calculations converged only when applying the homotopy-continuation method. For all the sets of specifications studied, a nearly pure ethyl acetate distillate was not obtained and a high conversion was not possible. The main problems were: (1) unfavorable total reaction conversion; (2) vapor-liquid equilibria; and (3) temperature profile. It was observed that K-values of ethanol and water were close, making separation of these two components very difficult. The volatility of acetic acid caused it to flow downwards and esterification occurred mainly in the lower part of the column. In the upper portion of the column, K-values of ethyl acetate were too close to ethanol and water, making the separation of ethyl acetate From ethanol and water in these stages very difficult. The temperature distribution in the column under different sets of specifications did not vary appreciably. Because the reaction conversion and K-values are strongly related to temperature, it was difficult to cause the temperature profile, reaction-rate profile, and K-values to change independently of each other. It is concluded that the application of reactive distillation to the esterification of ethanol with acetic acid is not technically Favorable in a two-product distillation column for the conditions imposed in this study. However, reactive distillation may be suitable for other reactive systems. In particular, the saponification reaction between propylene chlorohydrin and lime, as reported by Carra et al. (1979), to Form propylene oxide appears to be more suitable. The components involved have more Favorable volatility and reaction-rate characteristics than the system studied here. NOMENCLATURE A, B, C, ~1,b, c, a, y = Constants
A- =
y A,, A,._ . . . . A,,.” = AT, B,“, C,” =
C = Cj = C, = d = E, = Ep =
F, =
in the Bcnedict-WebbRubin equation of state The Arrhenius oreexoonential factor for the pth te& a coefficients for the Regression modified Mar&es equation Antoine constants of component j Number of components Molar heat capacity Concentration of the reference component for the qth reaction on stage i Molar density The discrepancy (residual) of the energy equation for stage i Energy of activation for the pth term
c f,,
,=I
f = Function fu = Liquid component molar flow for component j in the fresh feed to stage i
G, = i: g,, ,=I
Homotopy-continuation g, = Vapor component molar flow for component j in the fresh feed to stage i H = Homotopy function ZZ, = Enthalpy of the primary vapor leaving stage i; i.e. not including the vapor sidestream ZZG,= Enthalpy of Gi ffui = Molar enthalpy of the vapor leaving stage i IT;, = Virtual value of the partial molar enthalpy of component j in the vapor phase leaving stage i Hz. = Enthalpy of pure component j in the perfect gas state at temperature I: = Heat of formation of component j at AZ%, 298 K h, = Enthalpy of the primary liquid leaving stage i; i.e. not including the sidestream or entrained liquid h, = Liquid enthalpy of F, h, = Molar enthalpy of the liquid leaving stage i h;, ., = Virtual value of the nartial molar enthalpy of componen‘t j in the liquid phase leaving plate i K = Vapor-liquid equilibrium ratio Kti = Vapor-liquid equilibrium ratio for component j on stage i k, = Reaction coefficient for the pth term L, = i l,i ,=1 i, = Moles of component j in the primary liquid leaving stage i ; i.e. not including sidestream liquid or entrained liquid M,, = Discrepancy (residual) of the component material balance for component j and stage i N = Number of stages including condensers and reboilers NRXNS = Number of chemica1 reactions au = Exponent in chemical kinetics rate expression P = Pressure P, = Pressure on stage i P,” = Vapor pressure of pure component j Qi = Heat duty at stage i Q,_ = Discrepancy (residual) of the equilibrium relation for component j and stage i R = Reflux ratio; universal gas constant; molar reaction rate R,j = Initial guess for production rate due to chemical reaction of component j at stage i r, = dC,jdL, reaction rate for reaction n (with respect to the concentration of the reference component for the qth reaction) on stage i T = temperature; when used as a superscript denotes matrix transpose Tb = Normal boiling point temperature T, = Critical temperature T, = Temperature on stage i t = Time; temperature; homotopy parameter S, = Ratio of vapor drawoff to primary vapor (vapor not withdrawn or entrained on stage i s, = Ratio of liquid drawoff to primary liquid (liquid not withdrawn or entrained) on stage i un = Liquid holdup on stage i
for reactive distillation
1253
q = i:
=
cly
uij = Moles of component j in the primary vapor leaving stage i ; i.e. not including sidestream vapor nf = Liquid mol volume at stage i W, = F, + Gi x = Unknown x0 = Initial guess of unknown x, = Mol-fraction in liquid of component j x,/ = Mol-fraction in liquid of component j at stage i y, = Mol-fraction in vapor of component j vi, = Mel-fraction in vapor of component j at stage i Greek 2, = Molar latent heat of vaporization for component j at the temperature of stage i v,+,= Stoichiometric coefficient for component j in chemical reaction n R = Ri, - w;.. molar enthalpy departure function REFERENCES
Backhaus A. A., U.S. Patent 1400849 (1921). Benedict M., G. B. Webb and L. C. Rubin, An empirical equation for thermodynamic properties of light hydrocarbons and their mixtures. J. Chem. Phys. 10,747 (1942). Belck L. H., Continuous reactions in distillation equipment AZChE Ji. 1, 467 (1955). Bermans S., H. Isbenjian, A. Sedoff and D. F. Othmer, Esterification-continuous production of dibutyl phthalate in a distillation column. Znd. Engng Chem. 40, 2139 (1948). Carra S., M. Morbidelli, E. Santacesarla and G. Buzzi, Synthesis of propylene oxide from propylene chlorohydrins-II. Chem. Engng Sci. 34, 1133 (1979). Chan T. F., Deflation techniques and block-elimination algorithm for solving bordered singular systems. SIAM J. Sci. Srarist. Comput. 5, 121 (1984). Ferraris G. B., Interlinked, multistaged separations with nonstandard specifications solved by the NewtonRaphson method. AIChE JI. 27, 163 (1981). den Heijer C. and W. C. Rheinboldt, On steplength algorithm for a class of continuation methods. SIAM J. Namer. Analysis IS, 925 (1981). Izarraraz A., G. W. Bentzen. R. G. Anthony and C. D. Holland, Solve more distillation problems. l?Zydrocarbon Process 59, No. 6, 195 (1980). Jelinek J. and V. Hlavacek, Steady state countercurrent equilibrium stage separation with chemical reaction by relaxation method. Chem. Engng Common. 2, 79 (1976). Jing W., Q. Wang, W. Ling, X. Wang and Y. Yang, The simulation study of a reactor-three phases azeotropic distillation tower system. ZRSI. Chem. Engng Symp. Ser. 92, 675 (1985). Kaibel G., H. H. Mayer and B. Seid, Reaction in distillation columns. Ger. Chem. Engng 2, 180 (1979). Kinoshita M., I. Hashimoto and T. Takamatsu, A new simulation procedure for multicomponent distillation column processing non-ideal solutions or reactive solutions. J. Chem. Engng Japan 16, 370 (1983). Komatsu H., Application of the relaxation method for solving distillation. J. Chem. Engng Japan 10, 200 (1977). Komatsu H. and C. D. Holland, A new method of convergence for solving reacting distillation problems. J. Chem. Engng Japan 10, 292 (1977). Leyes C. E. and D. F. Othmer, Continuous esterification of
Y. A. CHANG
1254
and J. D. SEADER
butanol and acetic acid, kinetic and distillation consideration. Trans. AfChE 41, 157 (1945). Lin W., J. D. Seader and T. L. Waybum, Computing multiple solutions to systems of interlinked separation columns. AZChE JI. 33, 886 (1987). Naphatali L. M. and D. P. Sandholm, Multicomponent separation calculations by linearization. AZChE JI. 17,
where Z$=u;+n,
(A6)
Enthalpies for the ideal-gas state were computed
from:
AZChE
JI. II,
205 (1965).
Suzuki I., H. Komatsu and H. Hirata, Formulation and prediction of quaternary vapor-liquid equilibria accompanied by esterification. J. Chem. Engng JQPU~S3, 152 (1970). Suzuki I., H. Yagi, H. Komatsu and M. Hirata, Calculation of multicomponent distillation accompanied by a chemical reaction. J. Chem. Engng J~pm 4, 26 (1971). Wang J. C. and G. E. Henke, Tridiagonal matrix for distillation. ZZydrocarbon Process 45, No. 9, 169 (1966). Watson K. M., Thermodynamics of the liquid state, generalized prediction of properties. Znd. Engng Chem. 35, 398 (1943). Wayburn T. L., Modeling interlinked distillation columns by differential homotopy-continuation. Ph.D. Dissertation, University of Utah (1983). Wayburn T. L. and J. D. Seader, Solutions of systems columns by homotopyof interlinked distillation continuation methods. Proc. 2nd Znt. Conf. Foundation Computer-Aided Process Design (A. W. Westerberg and H. H. Chien, Eds). CACHE (1984). Xu X. E. and H. F. Chen, Simulation of distillation processes with reactions in series. Paper No. 5Of, 1986 Spring National Mtg, AZChE, New Orleans (1986). APPENDIX
A
H;
+ A+:
+ A&
+ ASx2xq + A6+xq 10xxx I
+A,x,x:+A
+
A,3~q~,~g
+ A,x,x:
Ci=ap+&T+y,Tz
C;dT,
(A7) (A8)
c
l-Tr2
AZ]=&
n
1
d
l-T,.,
(A9)
’
where the subscripts 1 and 2 refer to temperatures 1 and 2, and T, = T/T,. A common choice for n is 0.38 and for subscript 1 is the normal boiling point. Constants for equation (AS) are given by Izarrarax et Q/. (1980) in their Table 3.5. The Benedict et OZ. (1942) equation of state was selected to derive the enthalpy departure function: B,RT
P = RTd +
-A,,
- 2
d* + (bRT
(
-
a)d’
+ aad
>
+$[(l
(AlO)
+ydz)exp(-yd*)].
By introducing dimensionless groups into the BWR equation, the reduced equation can be written as: P,=T,d+d2x
I , (All) 1
T,B’-A’-$
A’ = A; + Q’d(l
- a’d”),
(Al2)
+ yd2)exp(-yd*),
(-414)
B’=B;,+b’d
(‘413)
C’ = C; -
c’d(l
P,=$
+ A,,x,x:
+ A,,x,x:
T, = f
c
c
(Al)
where x, is the liquid mol-fraction of component j, y, represents the activity coefficient, and A,-A,, are constants determined from quaternary equilibrium data by Suzuki et al. (1970) and given in their Table 3. The expression for the activity coefficient of the remaining components can be obtained by rotating the subscripts: I+ 2 + 3 * 4 4 1. The vapor pressure of pure liquid is calculated by using the Antoine equation: logP,o=A;--
298
using the expressions of Ixarraraz et ai. (1980) in their Tables 3 and 4. The latent heat of vaporixation decreases with increasing temperature and is zero at the critical point. A widely used correlation between A,, and T, is the Watson (1943) relation:
2 1+A II XXX 2 3 4+A 12XXX 3 . ,
+ A,,x&
+
where dimensionless groups are:
+ A,x,x, + A,x,x:
= AH&,,
where
K-Values Suzuki et al. (1970) rearranged the Margules equation as a polynomial series in mol-fractions of the components in the mixture: logy, = A,$
s r/
148 (1971).
Nelson P. A., Countercurrent equilibrium stage separation with reaction. AZChE JI. 17, 1043 (1971). Su G. J. and D. S. Viswanath, Part II. Generalized Benedict-Webb-Rubin equation of state for real gases.
(A5)
4) = R,, - rl,(Ti)
BO
A,=_
c;=_
Q’=-
c’=_
T+C,O
where the constants are given by Suzuki et al. (1970) in their Table 1. Vapor and liquid enthalpies The vapor and liquid enthalpies were calculated, respectively, by use of the virtual values of the partial molar enthalpies: c (A3) (A4)
CCIP, R=T4’c UP: R3T3’e CPZ R3TS’c
b’=’
AoPc R2T,2 bP f R2T2 c QEP,5 __ =R6Tz
”
Qa y’=-
YPf R2T2’c
Using a least-squares technique, a set of parameters for the generalized BWR equation was determined by Su and Viswamath (1965) who recommended the following generalixed constants: 47;= 0.04407059,
A; = 0.24180980, a’a’=0.11369478
x lO-3,
C; = 0.21217005,
b’=0.03715171, c’ = 0.06448001
B; = 0.07643101,
and
y’=O.O6.
The coel%ients for the BWR equation of state for the components of interest are given by Ixarraraz et al. (1980)
Homotopy-continuation
for reactive distillation
in their Table 3.6. In the mixture, each of the eight BWR constants are related to composition by mixing-rule equations of the form:
Values of /3 for each BWR constant are given Izarraraz et al. (1980) in their Table 3.7. From thermodynamics, the enthalpy departure function, 0, can then be obtained from the equation of state: R==
ew(-_yd*) 2
1255
+ yd*exp(-yd*)
. I
Kinetics The
reaction-rate expression for the reversible reaction: AcOH
+ EtOHuAcOEt
A
+ H,O
B
C D was obtained from Izarraraz ef ai. (1980): rA=k,C,C,-kkC,C,,
+ RT(& The resuit is: f-l=(B,,RT-2ZA0-~)d+(26R;-3a)dz 6aads t---+-3 5
ca2 T
I -
exp(--yd*) yd2
-
I>.
(A16)
(A17)
(‘4lW
where: kA = 29OOOexp(-7150/T),
6419)
k; = 7380exp(-7150/T),
(A20)
with kA and k_A in 1 gmol-’ min-’ and Tin “K. For comparison,, Suzuki ef al. (1970) used just the forward rate portion of equation (A18) with log,,k,,
=
-2710/T
+ 3.70.
6421)