620
Ind. Eng. Chem. Res. 1993,32, 620-633
Simulation, Design, and Analysis of Azeotropic Distillation Operationsf Bjarne S. Bossen, Sten Bay Jargensen, and Rafiqul Ganr Department of Chemical Engineering, Technical University of Denmark, 2800 Lyngby, Denmark
The computational tools needed for simulation, design, and analysis of azeotropic distillation operations are described. These tools include simple methods to identify the existence of binary and ternary azeotropes and to classify ternary mixtures as homogeneous or heterogeneous. The tools also include more complex methods to compute the phase diagram (or a heterogeneous liquid boiling surface), predict liquid-vapor phase equilibrium, and/or predict liquid-liquid-vapor phase equilibrium for simulations of batch and continuous distillation column operations. Important new features of these tools are the incorporation of a fast and efficient method for test of phase stability in simulation of distillation operations, the ability to handle a large range of mixtures (including mixtures with supercritical compounds), and the ability for computations covering wide ranges of temperature and pressure. On the basis of these tools, simple and consistent design algorithms are developed. The applicability of the design algorithms is verified through process simulation and analysis of the predicted behavior and data from the open literature. Conditions are given for examples illustrating (when and how possible distillation boundaries can be crossed) how multiple steady states can be obtained. Finally, the effect of changes in operating conditions on the dynamic behavior of the azeotropic distillation columns and the sensitivity of design to the prediction of phase equilibria are presented.
Introduction Azeotropicdistillation operationsrepresent complexbut interesting behavior which is also difficult to compute or simulate. Accurate and reliable design of these processes therefore needs accurate and reliable computational tools. Doherty and co-workers [Doherty and Perkins, 1978;Van Dongen and Doherty, 1985; Pham and Doherty, 1990a,b; Rovaglio and Doherty, 1990 (to name a few)] in their pioneering work in the area of azeotropic distillation have proposed several new concepts, design algorithms, and analysis of the behavior of the process under different operating conditions. These papers have also firmly established the use of ternary diagrams for design of the azeotropic distillation operations and have introduced the concepts of distillation boundaries and residue curves. Other workers, for example, Prokopakis and Seider (1983a,b),Kovach I11 and Seider (1987a,b),and Andersen et al. (19891, have also presented results that have contributed to a better understanding of the azeotropic distillation operations. Several workers (Magnussen et al., 1979; Kovach I11 and Seider, 1987a; Rovaglio and Doherty, 1990) have reported multiplicities in the numerical solution of the azeotropicdistillationcolumn model equations. In spite of the advances in recent years in the computational techniques and methods of solution, questions still remain on the consistency and reliability of the existing methods and algorithms. Questions regarding the reliability of predicted phase equilibria, which is a true distillation boundary, the existence of multiple states, and many more still remain unclear. Also, inconsistencies in the presently available design techniques have been reported (Laroche et al., 1990) recently. In this paper, an attempt is made to address some of the issues that still remain unclear within the context of a more clear understanding of azeotropic distillation operations. The paper shows why and how predicted phase diagrams can be used for design of azeotropic distillation
* Author to whom all correspondence should be addressed. +Presented a t the AIChE Annual Meeting, Los Angeles, November 1991; Session 60, paper no. 60d.
operations. It proposes a definition for a true distillation boundary and shows analytically and numerically that the so-called distillation boundaries (separatrices) can be crossed. Consequently, these boundaries should be more correctly labeled as possible distillation boundaries. Finally, the paper presents dynamic and steady-state simulation results and highlights the dangersof inaccurate/ inconsistent prediction of properties with respect to design and operation of heterogeneous azeotropic distillation columns. This paper therefore briefly describes the important features of some of the computational tools (for example, computation of ternary phase diagrams, fast and simple phase stability tests for the liquid phase, steady-state and dynamic distillation models, etc.) that are important for studies related to desigdanalysis of azeotropic distillation operations. The reliability of predictions of phase equilibria is discussed, and comparisons with experimental data are presented. Brief descriptions of the model equations and the method of solution are given, and new features are highlighted. For several design/analysis problems, the proper use of these tools is demonstrated and simple and reliable design algorithms are proposed. The simulation models are then used to verify the results from the design/analysis problems. At the same time, limitations of the proposed algorithms and the currently available computational tools are analyzed with respect to accurate/consistent prediction of properties.
Computational Tools In this section, a brief description of the computational tools needed for design/analysis of azeotropic distillation operations is given. The tools include computer programs for determination of binary and ternary azeotropes, computation of phase diagrams and residue curves (separatrices), and steady-state and dynamic simulation of azeotropic distillation columns. The model equations and the computationalmethods related to these toolsare briefly described below. The models and methods are valid for ternary mixtures and are not restricted to only the 7-4 approach for computation of phase equilibria. Two types of models (requiring the 7-4 and 4-4 approaches for
08SS-5885/93/2632-0620$04.00/00 1993 American Chemical Society
Ind. Eng. Chem. Res., Vol. 32, No. 4, 1993 621 computation of phase equilibrium) are considered in this paper. The first type of models are those based on UNIFAC (Fredenslund, 1989). These models are used for computation of the liquid-phase activity coefficients at low pressures. The second type of models combine activity coefficient models with equations of state (for example, the MHVB model (Dah1et al., 1991)). With the MHVB model, it is possible to predict phase equilibria for mixtures containing supercritical compounds and for wide ranges of temperatures and pressures. Note also that since with the MHV2 model the 4-4 approach is used for computation of phase equilibria, in principle any other equation of state can also be used. All previous work in this area has employed only 7-4 approach for computation of phase equilibria (for example, Pham and Doherty (1990a)). Determination of Binary Azeotropes. Unless two azeotropes (one low boiling and one high boiling) exist, the existence of an azeotrope can be detected by relative volatility defined as
For an ideal vapor phase and unity Pointing factor, eq 1 can be written as Y.
P?
Note that eq 2 is only valid at low pressures (with the 7-4 approach) and, at higher pressures, the 4-4 approach is preferred and the assumption of vapor phase ideality is not needed (eq 1is used). The value of a becomes unity at the azeotrope, and since eq 1 represents a continuous function, a has to be greater than unity on one side of the azeotrope and less than unity on the other side. Thus, the values of a at x = 0 and x = 1detect the existence of an azeotrope. If both values of a are greater than or smaller than unity, no azeotrope exists. If an azeotrope exists, the value of x corresponding to a = 1 can easily be determined. In this work, a "secant-like" method has been employed and the temperature of the binary azeotrope is determined by solving eq 3 together with eq 2 (7-4 approach) or eq 4 together witheq 1(4-4 approach). With
(4)
eqs 2 and 3, any gE-based model can be used, while with eqs 1 and 4, any equation of state can be used. Phase Diagrams. In this work, the ternary phase diagrams are described by a set of tie lines in a liquid boiling surface. The calculations are started on the binary lines of a ternary mixture. If all the binary lines are found to be stable (with respect to liquid-liquid equilibria), the system is considered to be homogeneous. Otherwise, the system is considered to be heterogeneous. Phase diagrams showing the heterogeneous liquid boiling surfaces are needed only for heterogeneous systems. The equations needed to compute each tie line in the heterogeneous liquid boiling surface are given below (fori = 1,2, ...,C; C indicates the total number of components). In the discussion below, all derivations are shown for the 7-4 approach. The corresponding equations for the 4-4 approach are given
in the Appendix. 1 1 0 = niyi - n i2y i2
(5)
O = ~ n ~ - N ,
(7)
i
0=Xnf-N, i
o = (niVni)- nipec
(9)
Note that in eq 9, for each tie line, the mole number for component 2 can be specified for liquid phase 1 or liquid phase 2. This will result in a set of seven equations with seven unknown variables (ni,nf, and 7'). If the system is heterogeneous, eqs 4-8 are solved successively with increasing values of nipec to determine the tie lines connecting the two liquid phases within the phase diagram. Nt in eqs 6-8 indicate total mole number. It should be noted that this formulation allows the computation of phase diagrams at low pressures (with gE-based models) and at high pressures (see the Appendix) with the MHVB model or other equations of state. The stability of the phase is checked with the method of Michelsen (1982). The Newton-Raphson method with analytical derivatives is employed to solve the set of equations (eqs 5-9) at different values (moles) of component 2 in liquid phase 1or 2. The tie line from the previously estimated point serves as initial estimate for the next point. The calculations start from a binary line and are continued until the critical point (point at which the two liquid phases have the same composition) is found. Note that each overall liquid composition also signifies the vapor composition in equilibrium with the two liquid phases. Thus the equilibrium vapor line corresponding to the two liquid phases within the phase boundary is indicated by the line starting from the binary azeotrope and ending on a point corresponding to the critical point. Note that each tie line or surface point represents also a distinct temperature and that the vapor line may end inside or outside the heterogeneous liquid boiling surface. Phase Stability. For the calculation of the separatrices and residue curves, it is necessary to know whether the system is homogeneous or heterogeneous. In the case of a heterogeneous system, it needs to be known whether any point on the curve is in the heterogeneous or homogeneous region. Typically, a check of phase stability through an analysis of the Gibbs energy (Michelsen, 1982) will provide the needed information. This however is computationally expensive when repetitive calculations are involved, and therefore, a faster and simpler method applicable only for ternary systems has been developed. First the phase diagram for the ternary system is determined. If the system is heterogeneous, then there will be a heterogeneous liquid boiling surface. Any composition within the surface will therefore have two liquid phases, and for all other cases, there will be only one liquid phase. The advantage of this method is that the liquid boiling surface can easily be precomputed, and once it is available, the check for phase stability is made by locating the liquid composition on the ternary diagram. If the liquid-phase composition is located within the phase boundary, it is immediately evident that there are two liquid phases. Also, the corresponding tie line passing through the located liquid composition point gives very good estimates for liquid-liquid phase split calculations.
622 Ind. Eng. Chem. Res., Vol. 32, No. 4, 1993
This procedure is new and is highly suitable not only for computations of residue curves but also for simulation of continuous and batch distillation operations (as shown below). With the inclusion of the 4-4 approach for computation of phase equilibria and equations of state, it is possible to study, as well, the effect of pressure on the heterogeneous liquid boiling surface. This feature is important in supercritical extraction processes. Separatrix and Residue Curves. The separatrices are a subset of residue curves, and in a ternary phase diagram, they indicate what may be labeled as possible distillation boundaries for the ternary system (Doherty and Perkins, 1978; Pham and Doherty, 1990a,b). The residue curves are calculated by simulating a simple batch distillation with no reflux and where the vapor in equilibrium with the remaining liquid is continuously removed. The following assumptions are employed in this work to model such a process: (a) The total vapor flow out of the system is constant. (b) The vapor phase is ideal. (c) There is equilibrium between liquid phase@) and vapor phase. The process can be represented by a single ordinary differential equation (Doherty and Perkins, 19781, dnpldt = -ui (10) In the above equation, np is the total liquid holdup of component i and ui is the vapor flow rate of component i. In addition to eq 10, the following algebraic equations are needed together with eqs 5 and 6 (note that there are C components in the mixture and two liquid phases may or may not exist),
case) or T (for homogeneous case) have been determined. For the heterogeneous case, eqs 5, 6, and 11 are solved with the Newton-Raphson method with analytical derivatives. By using the predicted values at the last time step as the initial estimates, convergence is achieved very rapidly (for most cases, between one and two iterations only are needed to achieve a convergence criteria of lV9). For the homogeneous case, eq 6 is solved for T (superscript 1 is changed to 0)and v is determined from eq 12. The calculations are started for either case from a specified overall liquid composition. Simulation Model. The dynamic distillation model of Gani et al. (1986) has been extended to handle azeotropic distillation operations with a stability check and computation of liquid-liquid-vapor phase equilibria in the column section. A decanter model for computation of liquid-liquid phase split has also been developed and added. For the computation of phase equilibria, the following equations are employed for the distillation and decanter models. Note that only the algebraic equations related to the computation of phase equilibria are given below. The ordinary differential equations (ODES) representing the mass and energy balance equations for a distillation column can be found in Gani et al. (1986). NJ
(13)
o = np - C n ; I
nfy!CV, 0 = v;-
In eq 11, superscript j indicates phase j . If there are two liquid phases, there are 4C f 1 equations (eqs 5,6,10,11, and 12) and 4C 1 variables (no, nl,n2,v, 7‘). Equations 5,6,and 11 describe a bubble point condition combined with a liquid-liquid phase split while eq 12 relates the bubble point condition to the molar flow rates. nofrom the 4C + 1 unknown variables are considered as the differential variables (that is, the integration routine determines their values at any time t ) , and the rest are considered as algebraic variables. In the homogeneous case, j = 1, and hence eqs 5 and 11 are not used and the superscript 1’s in eqs 6 and 12 are changed to 0’s leaving 2C + 1 equations and 2C + 1 variables. The differential-algebraic equation set obtained for the homogeneous or heterogeneous system is solved using the BDFSH package (Srarensen, 1991) and is based on the backward differentiation formula (BDF). Since the liquidliquid-vapor phase computation is very sensitive with respect to convergence, in this work the algebraic equations (eqs 5,6,11,12)are solved separately as a procedure (that is, solve the algebraic equations separately from the ordinary differential equations). Note that this simulation strategy has the advantage that it avoids the redefinition of the number of equations when liquid-liquid phase split is encountered on any plate. Examination of the equations reveal that eqs 5,6, and 11 need to be solved simultaneously and eq 12 can be solved after n1and T (for heterogeneous
+
In eq 13, nb is the total molar holdup of component i on plate p at time t (value of this variable is known a t time t). Equation 13 therefore indicates the distribution of component i in the different phases. Equations 14 and 15 are related to the isofugacitycriteria for phase equilibrium. Equations 13-15 are used for the distillation column and the decanter (if it is operating at the bubble point temperature of the mixture). If the decanter is operating at a subcooled temperature, eq 15 is omitted from the equation set and the temperature is specified or found from an energy balance over the distillate side of the condenser. Note that if the system on plate p is homogeneous, only eq 14 is needed thereby pointing to the need for a test for phase stability. In order to determine when the liquid-liquid phase split in the distillation should be considered in the calculation method and also to determine if the decanter is operating in the heterogeneous region, the method for phase stability test described above is employed for each plate and the decanter at every instant of time. If the decanter is operating at the bubble point, the check for phase stability is made by comparison between the overall liquid composition and a computed (a priori) heterogeneous liquid boiling surface. If the decanter is operating a t a constant temperature, the comparison should not be made with the heterogeneous liquid boiling surface but with a constant-temperature binodal curve. The computation of the constant-temperature binodal curve is a special condition of the method described above for the compu-
Ind. Eng. Chem. Res., Vol. 32, No. 4,1993 623 V
F1
T
+ B1
' K ! q ? h +- B
SQnetoh
OB 2070 MJ/h
Figure 1. Flowsheetfor azeotropicdistillation column with decanter. DIST, distillation column; DEC, decanter; C, condenser; REB, reboiler.
tation of the phase diagram. For low-pressureoperations, the effect of pressure is neglected since thegE-basedmodels are used. For high-pressureoperations however, the effect of pressure is taken into account by precomputing a series of liquid boiling surfaces at the pressures of interest. The interpolation procedure to obtain the estimates for the liquid-liquid phase split calculations becomes more complex in this case. In addition to the phase equilibrium relations (eqs 1315), the dynamic decanter model has the following differential-algebraic system of equations.
(16)
In the above equations, F is the molar flow rate into the decanter, F, is the molar flow rate out of the decanter corresponding to a full tank, z is the tank height, A is the cross-sectional area of the tank, FUand FLare the molar flow rate out of the decanter of the upper and lower streams, respectively, FTis the total molar flow rate out of the decanter corresponding to the actual liquid height (h)in the tank, H is the total molar holdup, 9 is the ratio between the two liquid phases (molar holdup) present in the decanter, and p is the molar density. In this model, eq 16 represents the component mass balance for the decanter, eqs 17-19 are "defined" equations relating the output streams to the input streams, and eq 20 determines the height of the total liquid phase. In order to solve eqs 1620 (together with eqs 13-15), three variables need to be specified (two of them must be F, and z ) in addition to F ,x:, and P. The remaining variable can be any one from 9,Fu, and FL. In Figure 1, a typical distillation column configuration with a decanter is shown. The remaining set of equations for the distillation column model (ODES from mass and energy balance, algebraic equations from plate hydraulics, etc.) can be found in Gani et al. (1986). The pressure on each plate is determined (if the dynamics of the energy holdup is considered) or is specified if the dynamics of the energy
holdup is assumed to be negligible. Alternatively, the dynamics of the energy holdup is included only at the top and the bottom plates to get a variation of pressure with time (a uniform pressure profile is determined from the top and bottom plate pressures). In the azeotropic distillation systems studied in this work, however, the pressure is specified equal on all plates since the columns operate at nearly atmospheric pressures. For supercritical extraction processes, the effect of pressure is important and, therefore, the simulation model considers the variation of pressure with time. The set of differential and algebraic equations (DAEs) are again solved with the BDFSH package (Sprrensen, 1991). Note that the system of DAEs representing this process has index 1, and it satisfies the criteria recently proposed by Gani and Cameron (1992)to avoid higher index problems. The algebraicequations are decomposed into two subsets: one subset consists of "explicit algebraic equations" (EAEs) with only one unknown variable per equation and the other subset consists of "implicit algebraic equations" (IAEs) with more than one unknown variable per equation. The EAEs are solved separately at every function evaluation (determination of the right-hand sides of the ODES and the IAEs). The IAEs may be solved together with the ODES as a DAE system or separately, as procedures (that is, solved separately every time the right-hand sides of the ODESare evaluated). In this work, the IAEs have been solved as procedures. Note that Perregaard et al. (1991) have shown that, in many distillation problems, solving the set of equations as a DAE system is not necessarily more efficient than solving as ODES and procedures.
Design/Analysis The tools that have been described above have been used to generate design procedures and for analysis of azeotropic distillation operations (use of phase diagrams rather than residue curves is encouraged). First an analysis of residue curves is given, then a simple design procedure (which should be used only as a first estimate) is proposed, and, finally, some implications of actual distillation boundaries and multiple feeds on the design and operation of azeotropic distillation systems are discussed. Analysis of Residue Curves. Doherty and Perkins (1978)have pointed out that the distillation boundaries (computed for batch distillation) or separatrices can be crossed by simple continuous distillation only if they are stable, meaning that the residue curves converge into it. A separatrix is unstable if all residue curves diverge from it. Residue curves are therefore computed to determine the stability of a separatrix. Pham and Doherty (1990b) have shown how these properties can be used for design and analysis of continuous heterogeneous azeotropic distillation operations. In this paper, it will be shown when, why, and how the distillation boundaries (computed for batch distillation) can be crossed (in simple continuous distillation or batch distillation with reflux) and that the crossing of a separatrix is not governed by its stability but rather by the location of the initial state of the mixture and the shape of the curve which consequently is labeled as a possible distillation boundary. We start with a definition of the shape of the separatrix. For example, a separatrix can be convex, concave, or a straight line. In order to determine the shape of the separatrices, we note that they describe curves in the ternary composition diagram given by (note that, for this analysis, the solvent
624 Ind. Eng. Chem. Res., Vol. 32, No. 4, 1993 BENZP(E
ETWOL
(== K)
1.0 A Azeotropes Ternary Azeotrope at 337.6 K
0.9
4\
0.5
4
Rddue-aurvo -___--Companding distillate-curve A z h ~ s
0.4 0.3 0.2 0.1 (3731 K)O.O WATER
0.1 0 2 0.30 . 4 0.5 0.60.7 0.80.9 1.0 (353.2 K)
Figure 2. Schematic mass balance on the ternary phase diagram for the system shown in Figure 1.
is component 1 and the solute with which it is partially miscible is component 3)
0.0 (=J 00.b 0.'1 0.b 0,.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
Residue-curve _ _ - -Corresponding distillate-curve ....._ tie-lines - - Liquid-vapour Moss balance line for plate 2
0.45 1 0.40
= 0.30 = 0.35
x; = f ( x 4
(21)
since x i can be determined when xp and x i are known. We can, therefore, define the shape of the separatrices (or residue curves) as concave, if
(22)
0.25 1 0.20
=
0.15 1
-
0.10
~,,,,1,,,,1,,,,1,,,,,,*,,,,,,,,,,~,, 0.41
0.42
0.43
0.44
0.45
0.46
0.47
0.48
Figure 3. (a, top) Residue-curve map for the system chloroformbenzene-acetone. (b, bottom) Mass balances based on residue curve.
convex, if
d2x:,0 -=o
straight, if
dxy
According to this definition, the separatrix between regions I and I1 in Figure 2 is concave (with respect to component 1) and region I is designated as the concave side while region I1 is designated as the convex side. In order to understand when and how these separatrices can be crossed, we have to look at the residue curves together with the corresponding vapor lines (distillate curves). Figure 3a shows a separatrix and a residue curve together with the corresponding vapor composition or distillate curve (the dashed curve). In Figure 3b, a part of the residue curve and a part of the distillate curve are magnified. If we now assume that the residue curves and the corresponding distillate curves can represent the behavior of continuous distillation columns (for example, a column operating at total reflux and with a single feed), then any point on these curves indicate the equilibrium stage conditions. Li and Vi, i = 1, 2, 3, therefore indicate the liquid and vapor compositions on stages 1, 2, and 3, respectively (see Figure 3b). The steady-state mass balance for stage 2 is therefore given by (assuming no feed enters or side streams leave from these stages)
+ v, - L, - v,
(25) This mass balance is represented by point 0 in Figure 3b. Note that point 0 here is equidistant from points LZand V:! and from points L1 and Vs, meaning that their ratios should be close to unity. Such a condition can indeed be satisfied by a column operating under continuous equimolar overflow conditions. In Figure 3, we can also note that 0 = L,
the vapor lines lie on the convex side of the residue curves. This will always be the case since the tangency rule which states that the slope of the tie line between the liquid and vapor composition is equal to the slope of the residue curve at the liquid terminus of the tie line is satisfied. For an overall liquid composition xoin equilibrium with a vapor of composition yo and an infinitesimal amount of vapor removed 6V, the composition of the remaining liquid x is obtained from eq 26. =
x ~ -Hyp 6V
(26) H-6V A consequence of satisfying the tangency rule and the use of eq 26 means that two residue curves cannot intersect since this will mean two different tangents at the same point and, therefore, two overall vapor compositions in equilibrium with the same liquid. Therefore,the following theorem is proposed: Theorem 1: Two residue curves on the convex side of a separatrix cannot intersect, and therefore, it is not possible to cross the separatrix from the convex side. Now let us consider the situations when the separatrices can be crossed. In Figure 4a, the heterogeneous liquid boiling surface and a separatrix are shown for the ethanolwater-benzene system. Point L in Figure 4a indicates a point on a residue curve which is on the concave side of the separatrix and very close to it. The dashed curve is the corresponding distillate curve, and according to the tangency rule, it must lie on the convex side of the separatrix. Now if we have a composition on the top plate in a distillation column close to point L, we will have a corresponding equilibrium vapor composition in region xi
Ind. Eng. Chem. Res., Vol. 32, No. 4, 1993 625 ETHWOL
(311.4
K)
- Seporatrix
-.....Corresponding distillate;cuwe
Heterogeneous liquid boiling surface Vapour line A Critical point 0 Azeotropes Ternary Azeotrope at 337.2 K
,G 9 0.7
0.6
4;
0.3 o.2 0.1
1
0.0 (373.1 K I O . 0 0.1 0.2 0.30.40.5 0.6 0.7 0.80.91.0 (353.2 WATER
K)
BENZENE
J,?, ETHANOL
(351.4
K)
,e 0.8
_.._._ Liquid
composition profile Heterogeneous liquid boiling surface Vopour line A Critical point 0 Azeotropes Ternory Azeotrope at 337.2 K
0.7
0.6
0.3 0.1
j
0.0 (373.1
KIO.0 0.1 0 . 2 0.30.40.50.60.7 0.80.9 1.0 (353.2 K) BENZENE
WATER
J*?
,4.
Figure 4. (a, top) Separatrix and corresponding distillate curve for the system benzene-ethanol-water (batch distillation). Point L indicates starting liquid composition. (b, bottom) Simulated (continuous distillation operation) steady-state liquid compositionprofile superimposed on the ternary diagram for the case with one feed stream at point L (top plate of a 20-plate column).
11. When this vapor is condensed in a total condenser, the reflux liquid stream from the condenser is in a distillation region other than the initial composition on the top plate (in agreement with Pham and Doherty (1990b)). If we assume that the initial column composition profile follows a residue curve, then for a top plate composition at point L the liquid on the plate below would also be lying close to the separatrix. Therefore, the vapor from the plate below has a composition in region 11. The consequence of this is that both the liquid stream (that is reflux) and the vapor stream (from plate below) entering the top plate have a composition in region 11. Hence, the composition on the top plate will have to move into region 11,crossing the possible separation boundary. Since the separatrix is unstable, the liquid leaving the top plate will move further away from the separatrix and into region 11. The liquid leaving the top plate will then move the composition on the plate below closer to the separatrix (or even cross it), thus driving the vapor from the second plate further into region 11. If no feed stream is introduced further down the column, the result will be that the entire profiles moves to the other side of the separatrix. Thus, the condensation of the top liquid starts a force that drives the profile to the other side of the separatrix. If we consider a feed (of composition L)to the column
on plate NP, the location of the correspondingplate liquid composition will depend on the liquid and vapor streams entering the plate. Also, one can say that the liquid composition profiles moves from one residue curve to another. Therefore, the following theorem is proposed Theorem 2: All curved separatrices can be crossed from the concave side of the separatrix if it is possible to obtain a liquid composition profile close to the separatrix and within the concave region. Figure 4 shows the verification by simulation of the crossing of the possible distillation boundary described above. In Figure 4a, the location of the vapor line (distillate curve) corresponding to a liquid composition close to the separatrix is shown for a batch distillation operation. In Figure 4b, the simulated liquid composition profile is shown for a continuous distillation operation. In this case, a feed stream of 1 kmollh with a composition given by point L ( x i = 0.35,0.525,0.125; i = 1,2,3) is fed to the top plate of a 20-plate column. The steady state is obtained from the simulation in terms of the distillate curve (dashed line). This steady state clearly verifies the possibility of crossingthe separatrix by continuous distillation. It should be noted also that because the distillate line (vapor line) always lies on the convex side of the liquid composition line, it is not possible to cross the separatrix from the convex side. This also has been verified by simulation. Finally, based on the analysis given above, theorem 3 is proposed: Theorem 3: A true distillation boundary must always be a straight line. The proof follows from the above since all curved boundaries can be crossed from the concave side. A Simple Design Procedure. The first decision with respect to azeotropic distillation (continuous operation) is related to the type of the system-whether the system is homogeneous or heterogeneous, and where the binary azeotrope(s)and the ternary azeotrope(s)are located. From the information on the type of the system and the location of the azeotropic points, a preliminary design of the separation system can be obtained. Simple design procedures serve to provide first estimates and limitingvalues. Since all curved possible distillation boundaries can be crossed from the concave side to the convex side (note that the equilibrium vapor line is always on the convex side), in order to obtain the limiting values, it is simpler to consider a true distillation boundary (see theorem 3 above). Consider the following step by step procedure (assume that at this stage the system type has been identified and the locationsof the azeotropicpoints have been determined by the tools described above). This procedure does not use residue curves but restricts the possible separation by employingtheorem 3. Therefore,at this stage of the design procedure, all distillation boundaries are assumed to be straight lines. Note that this does not mean the design is looking for systems that will provide straight line distillation boundaries. Assume for simplicity that there is only one binary azeotropic point per binary pair and only one ternary azeotrope (if it exists) for the ternary system. Assume also that solute A is to be separated from solute B and the entrainer (or solvent) is component S. Note also that, for a ternary system, a two-column configuration is needed where solute B can be obtained as the bottom or top product from column 1,solute A can be obtained as the bottom or top product from column 2, and the solvent (with some solutes) is recycled back from column 2 to column 1 (note also that two liquid phases may exist in both columns). The temperatures of the
626 Ind. Eng. Chem. Res., Vol. 32, No. 4, 1993
binary and ternary azeotropes and the pure component boiling points indicate the type (bottom or top) of products. Step 1. Determination of approximate possible distillation boundaries: Draw lines by simply joining the ternary azeotropicpoint (if it exists) with each of the binary azeotropic points. If all the three binary azeotropes exist, this will result in three regions. In the case of a homogeneous system with a single binary azeotrope, a ternary azeotropewill usually not exist. In this case,join the binary azeotropic point with the pure component point of the third component. The pure component boiling points (temperatures) and the temperatures of the binary (also ternary azeotropic points) need to be considered to define the regions when more than one binary azeotrope exists. In this case, the procedure given by Foucher et al. (1990) can be employed. Note that the total number of regions can be less than, equal to, or greater than the number of components in the system and that the type (light or heavy) of product is governed by the temperatures of the pure components and the azeotropic points. Step 2. Overall feed compositions: The overall feed composition for column 1may lie inthe region bracketing the binary azeotrope between A and B, the ternary azeotropic point (if it exists), and the binary azeotropic point between A and S (if it exists) or in the region bracketing the binary azeotrope between A and B, the ternary azeotropic point (if it exists), and the binary azeotrope between B and S (if it exists). Note that, in the case of former, this will mean a crossing of the distillation boundary. The overall feed composition for column 2 should lie in the region bracketing the binary azeotrope between A and B, the ternary azeotropic point (if it exists), and the binary azeotrope between A and S (if it exists). Step 3. Product compositions: From the first column, the limiting product composition can be pure solute B. From the second column, the limiting product composition can be pure solute A. The location of the products (top or bottom product) depend on the pure component boiling temperatures and the azeotropic point temperatures. Step 4. Feed composition to the decanter (only for heterogeneous systems): The liquid composition entering the decanter must lie within the heterogeneous liquid boiling surface (if at boiling point-the case of subcooling can be considered at a later stage of the design procedure). Typically, the lighter product (or the solvent-rich phase) from the decanter is sent back to the column as reflux and the heavier product is sent as feed to the other column. For this step, the phase diagram for the heterogeneous system is needed. Note that, at this stage of the design process, crossing of distillation boundaries is not being considered (except in the case of overall feed composition for column l-see step 2) and, therefore, the separatrices are assumed to be straight lines (see theorem 3). The actual boundaries are determined later using the method outlined earlier (see the section on computational tools) for a more accurate design. A typical schematic mass balance for the azeotropic distillation operation is shown in Figure 2 where the actual curved separatrices are shown. In Figure 2, solute A is water, solute B is ethanol, and solvent S is benzene. Point B is the bottom product, point F1 is the ethanol-water feed mixture, point D1 is the entrainer feed (pure solvent), and point V is the vapor from the top plate which splits into two liquid phases in the decanter. The three so-calleddistillation regions are indicated by regions I, 11, and 111. Point FO, indicating the overall feed composition is obtained from the intersection of the lines
D1 V
1
1I I -
DlST
~
REC
D
F1
Bl ETHWOL (351.4 KI
Ternary Azeotrope at 337.8 K
(573.1
K)O.O 0.1 0.20.30.40.5 0.60.7 0.80.91.0 (353.2
WATER
K)
BENZENE J*?.J
4
Figure 5. (a, top) Column sequence for a two-column azeotropic distillation operation. (b, bottom) Schematic mass balances for azeotropic distillation operation.
D-B and F1-D1. Note that FO is in region I while B is in region 11, indicating that a crossing of the possible distillation boundary must take place. The step by step procedure given above is simple and needs only the information related to the phase behavior of the system. If more accurate analysis is needed, the simple distillation boundaries proposed by Wahnschafft et al. (1992) could also be computed (the proof given above in this paper for crossingof possible distillation boundaries corresponds to the definition of distillation boundaries given by Doherty and Perkins (1978)). Note, however, that computation of any type of distillation boundaries also needs the same phase behavior information used above (see Computational Tools). Design of Distillation Operations with Multiple Feeds. Figure 5a shows the schematic mass balance for the system shown in Figure 5b (this is similar to that reported by Pham and Doherty (1990b)). In Figure 5, B1 is assumed to be pure solute A (water) and point D1 is determined from the intersection of the lines through B1-D and F1-L (notethat F1, FO, D1, and L must lie no a straight line). Since F1 does not contain the entrainer (solvent S or benzene), it can, theoretically, be anywhere on the solute B-solute A (ethanol-water) binary line but, actually, lies below the binary azeotropic point. The location of the points D1 and FO is therefore dependent on the location (or specification) of points B1 and V and the shape of the
Ind. Eng. Chem. Res., Vol. 32, No. 4, 1993 627
Uquld mole fmctlon of benzene _--mole fmctlon of ahanol - - Uquld Uquld mda fmctlon of water
Uquid mole fmctlon of benzene Uquid mole fmctlon of ethanol Uquid mole fmctlon of water
.- - - - - - -
-
\
$
\
\\ \
0.7
II \
\
r
*z 0.6 0.5
d 0.3
= 0.2 /
/ 1 1 1 1 1 1 1 1 1 1 1 1 1 , 1 1 1 1 1 1 1 1 1 ~ 1 ~
l l l l l l l l l l l l l l l , l l l l l l i ~ l ~
1 3 5 7 9 11 1315171921232527 Plate number ETHANOL
(351.4 K)
1.0
A
1 3 5 7 9 11 1315171921232527 Plate number ETWOL (351 4 K)
_____. Liquid cornposition profile
Heterogeneous liquid --- --- Vomur line
boiling surfoce
_.._._ Liquid composition
profile Heterogeneous liquid boiling ourfoce Vopour line A Critical point a Azeotropes Ternory Azeotrope at 337.2 K
A Cdtical point
,,\
a Azeotropes
%
0.8
,5\
Ternory Azeotrope at 337.2 K
4 \\
3
0.7
0.6 0'7
0.6 0.5
,,' -
-\&
*
0.4
0.0 (373.1 '00.00.1 WATER
0.2 0.30.4 0.5 0.6 0.7 0.8 0.9 1.0 (353.2 342.4
K
K)
BENZENE
WATER
342.4 K
BENZENE
Figure 6. Simulated multiple steady state 1 composition profiles (a, top) and liquid composition profile superimposed on the ternary diagram (b, bottom).
Figure 7. Simulated multiple steady state 2 composition profiles (a, top) and liquid composition profile superimposed on the ternary diagram (b, bottom).
separatrix. If the separatrix is concave with respect to the solvent (as in Figure 5b), for all points of D1 in region I, FO is also located in region 1. If, however, the separatrix was convex with respect to the solvent, for all points of FO in region 11, D1 will also be located in region 11. Also, if V moves toward the phase boundary, D1 moves away from the separatrix in the case of a concave separatrix with respect to the solvent and toward the separatrix for a convex separatrix. Furthermore, since B1, D, and D1 lie in a straight line, D1 moves away from the separatrix if the solvent content in B1 increases for a concave separatrix (with respect to the solvent). The opposite is true if the solute B (ethanol) content of B1 increases. Therefore, in order to minimize the loss of solvent, D1 must be located nearer the separatrix in the case of the concave separatrix. Since this will cause loss of solute B (compound to be recovered), B1 should be located as near as possible to pure solute A. Crossing of the possible distillation boundary may occur in both columns. For the recovery column, crossingoccurs when D1 is in region I1 or 111,since D is in region I. Note however that the crossing is possible only when the separatrix is concave with respect to the solvent and not convex. Thus, location of D1 in region I1or I11represents an infeasible design for a convex separatrix. For the azeotropic column, crossing occurs when FO is in region I, since B is in region 11. Since crossing is not possible in the case of a convex separatrix, the location of FO in region I in the case of a convex separatrix represents an infeasible design. Since the location of FO and D1 may lead to infeasible design, instead of determining D1 and FO from specified B1, F1, and V, a safer approach may be to simply choose a point on the extension of B1-D for the location of D1. A consistent location for FO and F1 can then be obtained
by extending the line L-D1. This procedure may also be used for a single-column configuration (that is, without the solvent recovery column). In this case, the slope of the line D1-L line can be adjusted by the location of the solvent makeup. This procedure has been used to obtain the preliminary design of the azeotropic column (column 1 in Figure 5a; details given in Figure 1). Steady-state and dynamic responses for this column have been computed with the computational toolsdescribed earlier. Also, simulations have confirmed that it is indeed possible to obtain a product B in region I1 with FO located in region I (for the ethanol-water-benzene system shown in Figure 5b). Steady-state and Dynamic Simulation. Simulations have been made for the dehydration of ethanol in a 27plate column with a feed tray on plate 23 (numbered from the bottom) of a binary mixture of ethanol and water using benzene as the entrainer. The results have shown that it is possible to obtain ethanol of 99.997% purity from a feed stream of 89% ethanol and 11%water. With the same specifications of input streams, reflux ratio in the decanter, and reboiler heat duty, it was possible to obtain four different steady states. Three of the steady states (feasible steady states) corresponded to the desired separation, and the fourth (infeasible steady state) corresponded to no separation due to loss of phase interface in the decanter. The composition profiles are shown in Figures 6a-9a. From Figures 6a, 7a, and 8a, it can be noted that the difference between the steady states is the location of the front where the ethanol concentration decreases and the benzene concentration increases. The two product streams leaving the column in the three cases have exactly the same compositionsand flow rates. When the profiles are superimposed on the ternary phase diagrams (Figures 6b, 7b, and 8b), we can observe that,
628 Ind. Eng. Chem. Res., Vol. 32,No. 4,1993 Uquld mole f m d o n of benzene Uquid mole fmdlon of ethanol Uquld mole f m d o n of water
-Uquld mole f m d o n of mole fmctlon of - - Uquld Uquld mole fmctlon of
1 .o
'-.
benzene ethanol water
\ \
I I \ \
P 0.6
f-2
e
d =
0.5 0.4 0.3 0.2
/
0.0
i l l l , l l l i , i , l l l l i l l l l l i ~ i ~
ETWOL (351.4 K)
1.0
1 3 5 7 9 1 1 1315171921232527
1 3 5 7 9 111315171921232527
Plate number
Plate number
7
......... Liquid composition profile Heterogeneous liquid boiling --- --- VoDour line
ETHANOL
(351.4
surfoce
...._. Liquid
composition profile Heterogeneous liquid boiling surfoce Vapour line A Criticol point Azeotropes Ternory Azeotrope at 337.2 K
K)
A Crkicol point
,5\
o Azeolropes
3
f
Ternory Azeotrope at 337.2 K
0.8
0.7
0.7
o,2 0.1 0.0 (J73(
K)o.o0.1 0.20.30.40.50.6 0.70.80.9 10 (3532
WATER
342.4 K
0.0 K)
BENZENE
Figure 8. Simulated multiple steady state 1 composition profiles (a, top) and liquid composition profile superimposed on the ternary diagram (b, bottom).
except for the top section of the column, the three profiles are almost identical. This shows that the composition profile up through the column and that the actual location of the front have importance. Another important result is that, in all cases, the possible distillation boundaries have been crossed (from region I11 to region 11-in terms of feed plate location), indicating the effect of the reflux from the decanter on the column. The reason that there is no separation for the fourth steady state is that the decanter operates in the homogeneous region and therefore becomes a simple tankmixer. As a result, too much water is recycled into the column and the desired separation is not achieved. The first steady state was reached through dynamic simulation from an initial condition equivalent to an estimate of the expected steady state. The steady-state option of the simulator was then used to verify this steady state. The two other feasible steady states were obtained by manipulating the reboiler heat duty. First the reboiler heat duty was changed by +5%. This change was maintained until the front was near the top of the column. Then the reboiler heat duty was returned to its original value. This gave the second steady state. The third steady state was obtained in a similar manner but with a -5% change in the reboiler heat duty. The fourth steady state was obtained in a similar manner as the second, but the reboiler heat duty was kept at the increased value for a longer period. This pushed the front so far up the column that the vapor leaving from the top plate could no longer be condensed into a heterogeneous liquid. To check for the stability of the steady states, dynamic simulations for an additional 50 h were carried out. At this point, the eigenvalues for the Jacobian matrices were evaluated. Negative eigenvalues (real part) were found for the three feasible steady states. Note that the Jacobian matrix is
4
\\
1
(373.1 K)O.O
0.1 0.20.30.40.50.6 0.7 0.80.9 1.0 (353.2
WATER
342.4 K
K)
BENZENE
Figure 9. Simulated multiple steady state 2 composition profiles (a, top) and liquid composition profile superimposed on the ternary diagram (b, bottom).
fmctlon of ethanol on top plote --- Mole Mole fmctlon of ethanol on feed plate --
.------.
Mole fmctlon of ethanol on bottom plote
---___
,---_ \
0.0
B 0
2
4
6
8 10 12 14 16 18 20
Time (hours)
Figure 10. Dynamicresponseto +5% positive step change in reboiler heat duty of the steady state in Figure 7a.
evaluated with respect to the ODESonly. The steady state here is defined as the state where the residuals for all the ODES and IAEs are less than 10-6. Note that since the EAEs are solved analytically, they therefore have exact solutions. Dynamic simulations have also been made for step disturbances in the entrainer feed stream and the reboiler heat duty with the four steady states as the reference point. Except for the homogeneous steady state (where the disturbances had no effect), the other three steady states all showed the same behavior (leading to different steady states). The common behavior was a movement of the front upward in response to increased reboiler duty and decreased entrainer feed and a front movement downward for decreased reboiler duty and increased entrainer feed. The transient responses for the liquid mole fraction of
Ind. Eng. Chem. Res., Vol. 32, No. 4, 1993 629 Ethanol
Ethanol
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0.0 Water
Ethanol
E tl
IO1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0.0
0.0
WaterO
0
Watero.O
0.2
0.4
0.6
0.8
tiU
2
Figure 11. COz-ethanol-water predicted (MHVBmodel) and experimental ternary liquid-liquid equilibrium at T = 304.2 K and P = (a, upper left) 69, (b, upper right) 67.8, (c, lower left) 66.9, and (d, lower right) 66.8 atm.
ethanol on the top plate, the feed plate, and the bottom plate are shown in Figure 10 for an increase in the reboiler heat duty. Finally, a word about the dynamic simulator. Because of the proposed phase stability test procedure, it has been possible to run all the simulations at very low computational times. Also, there have not been any numerical problems due to failures in the convergence procedures for phase equilibria computations. In fact, the analytical derivatives and the Newton method employed to solve the phase equilibria problem are very robust and efficient. For most computations (afterthe fist function evaluation), between one and three (mostly one) iterations are needed. Reliability of the Results. The simulated results have been compared with those reported by Magnussen et al. (1979),Prokopakis and Seider (1983b),and Rovaglio and Doherty (1990). Also, the phase diagrams have been compared with published data (Gmehling et al., 1981). The comparison of results with those from Magnussen et al. (1979) and Prokopakis and Seider (198313) showed some differences (Pedersen, 1990). These differences can be attributed to the fact that Magnussen et al. and
Table I. Comparison of Experimental and Computed Binary Azeotropes (Azeotropic Composition and TemDerature at 1 atm) exptl system variable0 value benzene (1)-methanol ~ ( 1 ) 0.3969 T 330.65 heptane (1)-methanol r(1) 0.2316 T 332.25 ethyl acetate (1)-water r(1) 0.6885 T 343.53 benzene (1)-2-propanol r(1) 0.5720 T 344.35 water (1)-2-propanol r(1) 0.3112 T 351.47 benzene (1)-water r(1) 0.7003 T 343.65
computed value mod. UNIFAC MHVB 0.3862 0.3888 330.33 331.18 0.2577 0.2483 331.74 331.68 0.6775 0.6865 344.26 343.98 0.6024 0.6128 344.49 344.98 0.3062 0.3087 353.41 352.97 0.7019 0.7000 342.50 342.30
nr(l), mole fraction of component 1 at azeotropic point; T, temperature (K) at the azeotropic point.
Prokopakis and Seider used not a decanter model in their simulation but, rather, a constant recycle stream (to
630 Ind. Eng. Chem. Res., Vol. 32, No. 4, 1993 ETHANOL (508.2 K)
ethanol
1 .o Heterogeneous liquid boiling surface Vopour line A Criticol point
O" 0.6
1
\
-
,\
\ (54659 K)O.O WATER
MHVI, MHVZ
***** Azeotrope
0.8
0.8
0.5 0.4
--- vapor%
Hetero rneous li uid boiling surface UNIFAC Vapor 4ine U N I F A ~(VLE) M k e o t r o o e UNIFAC (ME) Hater oneour liquid bbiling surface MhV2
0.9
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
0.0 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
water
benzene
(261.8 K)
Heterogonoous liquid boiling surface --- _-A _Vapour line Critical point
CARBON DIOXIDE
Azeotropes
(416.6 K)
WATER
BENZENE
00 0.9 LO2 10 Figure 12. COz-ethanol-water ternary heterogenous liquid boiling surface at P = 60 atm computed with MHVB model (b, bottom) showing a magnified portion of the COz corner).
represent the decanter). The reflux composition obtained from simulation in this work was found to be x(benzene) = 0.74, %(ethanol)= 0.22, and x(water) = 0.04 while the values used in the two cited papers were x(benzene) = 0.54, x(ethano1) = 0.31, and x(water) = 0.15. Comparison with the results of Rovaglio and Doherty (1990) shows very good agreement even though Rovaglio and Doherty used the UNIQUAC model for computation of liquid activities, while in this work, modified UNIFAC and UNIFACWLE) were used (results with UNIFAC(VLE)were closer to those obtained by Rovaglio and Doherty). It can be concluded, therefore, that the two simulation models yield reliable results and can be used in simulation and design of azeotropic distillation operations (at least for the ethanol-water-benzene system). One must however be careful that multiple steady states are not obtained due to numerical artifacts. Also, all steady states must satisfy the global component and energy balances up to a specified accuracy. Even though we have obtained good agreement with reported results by other workers, the problem of reliability of phase equilibrium predictions of liquid-liquid-vapor systems needs to be addressed. Different thermodynamic models may predict different phase behavior. In this paper therefore, the performances of the UNIFAC-based models and the MHVB model are compared with experimental data (when available) as well as with each other (to study the effect of predicted phase behavior on design). Table I compares the computed binary azeotropes with the modified UNIFAC model and
JQo mu"AN (423.4 K) 6
Heterogeneous liquid boiling Vapour line Critical point Azeotropes
0.6
0.5 0.4
31 \
\
0.1
0.0 (452.7
\
,,
U0.0 0.1 0.20.30.4 0.50.6 0.70.80.9 1.0 (453.0 K)
WATER
BENZENE
?'e.!
Figure 13. Benzene-ethanol-water ternary phase diagrams predicted with MHV2 at P = (a, top) 1,(b, middle) 5, and (c, bottom) 10 atm. Table 11. Design Specifications for the UNIFAC(VLE) and MHV2 Models That Correspond to Similar Computed Profiles for Continuous Distillation Operations
model UNIFAC(VLE) MHV2
Fl/Dl 24.8 13.4
D/B 0.45 0.65
the MHVB model with experimental data reported in Pham and Doherty (1989a). It can be seen that both models
Ind. Eng. Chem. Res., Vol. 32,No. 4, 1993 631 WATER (373.1
1 .o
G
0.9 0.8 0.7
0.6
0.5 0.4
0.3 0.2
0.1
0.0 (337-7 G0.0 0.1 YEl"0L
0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 2-PROPANOL
WATER P3.1 G
1.0
0.9 0.8 0.7
0.6
0.5 0.4
0.3 0.2
0.1
0.0
1 -0
(337.7 QO.0 MEl"OL
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
(355.7
next two examples. Figure 13 shows the ternary phase diagrams for the ethanol-water-benzene system at 1,5, and 10 atm. At 1 atm, the predicted phase diagrams (Figure 13a) by the two models (UNIFAC and MHV2) compare apparently fairly well. However, when the distillation operation was simulated using the MHV2 model with the same specificationsas for UNIFAC, it was not possible to obtain the same steady state. The benzene makeup rate and the reflux ratio had to be varied in order to obtain the same composition profiles. The reason for this difference is that the computed tie lines are different. Even though the difference is not very large, it is sufficient to change the decanter and top stage conditions of the column (composition as well as temperature). Thus one must take care to check that even if the phase boundary is predicted accurately, the tie lines may have different slopes. This would not only change the results quantitatively, but also will affect design decisions (for example, there may or may not be liquid-phase split in the column section). Table I1 gives the values of FUD1 and D/B that are needed for the UNIFAC model and the MHVB model to obtain approximately the same composition and temperature profiles. The second example shows the qualitative effect on design as a consequence of the use of different thermodynamic models. Figure 14 shows, in ternary diagrams, the location of binary azeotropes, the possible distillation boundaries, and some residue curves. It can be noted that while Figure 14a (obtained with UNIFAC(VLE)) shows a convex separatrix (from region I), Figure 14b (obtained with modified UNIFAC) shows a concave separatrix from the same region. Thus in the case of Figure 14a, it is not possible to cross from region I to region 11, while in the case of Figure 14b, it is possible to do so.
G
2-PROPANOL
Figure 14. (a, top) Residue curve map for the system 2-propAno1water-methanol simulated with modified UNIFAC. (b, bottom) Residue curve map for the system 2-propanol-water-methanol Simulated with UNIFACNLE).
predict these values to acceptable accuracies. For the MHVB model, Figure 11 compares the computed and experimental (reported by Dah1et al. (1992))liquid-phase splits for the system carbon dioxide-ethanol-water. It can be seen that the three-phase regions are predicted by the MHVZ model with more accuracy than the phase compositions, indicating that the effect on design would be quantitative rather than qualitative. It should be noted however that these are very difficult calculations and very few models can match the predictions of the MHVZ model. For the case of the ethanol-water-benzene system, however, it is shown below that the predicted phase boundary and the phase compositionsaffect the design qualitatively as well as quantitatively. From a design point of view, it is also important whether the vapor line is inside the heterogeneous region or outside. The computed heterogeneous liquid boiling surface for carbon dioxide-ethanolwater at 60 atm is shown in Figure 12. Figure 12bmagnifies the carbon dioxide end of the ternary phase diagram. It can be seen that the vapor line is actually outside the heterogeneous region, indicating that a decanter is not needed. By reducing the pressure, however, it is possible to bring the vapor line inside the heterogeneous region. The operating pressure is therefore an important variable in the design of supercritical extraction processes. The effect on design due to predictions of phase behavior by different thermodynamic models is analyzed in the
Conclusions Computational tools that are needed for solving problems related to simulation and design/analysisof azeotropic distillation operations have been developed, and with examples, their use has been demonstrated. It has been shown how ternary phase diagrams (with the 7-6and 6-6 approaches) can be computed and how they can be used for design of the azeotropic separation process. The main features of the computational tools, a fast stability test procedure to be used in distillation computations and the ability to study systems which operate at high pressures, have been highlighted with examples. The concept of possible distillation boundaries has been analyzed, and it has been shown that true boundaries must be straight lines; otherwise, they can be crossed from the concave side. A new definition of a distillation boundary has recently been proposed by Wahnschafft et al. (1992). Since the equilibrium vapor lines for these boundaries also lie on the convex side, it remains to be seen if they can also be crossed. Simpledesign algorithms (based on distillation boundaries) for continuous distillation column operation therefore should use the definition of a true distillation boundary. A simple design algorithm which requires only the use of phase diagrams and not residue curves has been proposed. The design/analysis results have been illustrated with simulations. With the simulation model, multiple steady states have been found. Some typical conditions which produce these steady states have been identified (the only differences between the different steady states is the location of the "front" and the liquid composition curve of the top few plates in the column). This result leads to the speculation that the number of multiple (feasible)
632 Ind. Eng. Chem. Res., Vol. 32, No. 4, 1993
steady states is related to the number of possible locations for the “front”. This number is also related to the number of plates in the column. Current work is under way to study the phenomenon in greater detail. It has been shown, also, that infeasible steady states can be obtained. The effectsof changes in the design variableshave been studied, and it has been noted that the three feasible steady states show similar behavior. It is hoped that the results and analysis presented in this paper will contribute to a better understanding of the behavior of azeotropic distillation operations. Also, although the results obtained with the computational tools seem to agree well with other workers, there is a danger of making significant errors for other systems. To reduce the possibilities for design errors (qualitative as well as quantitative), the phase behavior should be predicted accurately. The comparisons of simulated distillation operations and the corresponding predicted phase diagrams have indicated potential problem areas. Further work therefore is needed not only in improving the reliability of predicted phase equilibria, but also in improving and extending the areas of application of the simulation models (for example, distillation operations involving supercritical gases and polar compounds). Currently, work is being done to address the above-mentioned issues as well as study, in greater detail, the causes and effects of the multiple steady states.
Acknowledgment The authors would like to thank Professor W. D. Seider, Professor M. L. Michelsen, and Dr. S. Widgado for previewing the manuscript and for their helpful suggestions and comments. The authors would also like to thank Allan Hansen for incorporating the MHV2 model in the computational tools. Nomenclature A = cross-sectional area of the decanter, m2 B = bottom product flow rate, kmollh D = distillate product flow rate, kmol/h E = entrainer feed (makeup) flow rate, kmol/h F = molar flow rate into the decanter, kmol/h FL = molar flow rate out of the decanter (heavier phase), kmoV h Fs = total molar flow rate out of the decanter corresponding to a full tank, kmol/h FT = total molar flow rate out of the decanter corresponding to the actual liquid height, kmol/h FU = molar flow rate out of the decanter (lighter phase), kmol/h F1 = feed flow rate (mixture to be separated), kmol/h Fo= overall feed flow rate for azeotropic distillation column (DIST), kmol/h h = actual liquid height in the decanter, m H = total molar holdup, kmol Kip = equilibrium constant for component i on plate p L , = liquid flow rate from plate p , kmol/h ni = molar holdup of component i, kmol P = pressure, atm Q b = reboiler heat duty, MJ/h Qc = condenser heat duty, MJ/h t = independent variable (time), h T = temperature, K u, = molar volume of component i, rnVkmol V,, = vapor flow rate from plate p , kmol/h VT = total volume, m3 x , = liquid-phase mole fraction of component i yi = vapor-phase mole fraction of component i z = decanter tank height, m
Greek Letters a = relative volatility y = activity coefficient
Z = molar density, m3/kmol = ratio between two liquid phases (see eq 19)
9
Subscripts
i = component i j = component j or phase j k = component k
P = Plate P 1 = component 1 or plate 1 2 = component 2 or plate 2 3 = plate 3 Superscripts j = phase j
L = liquid phase o = total (before phase split) s = vapor pressure V = vapor phase 1 = phase 1 2 = phase 2
Appendix Computationof VLLE with the &#Approach. The equations given below and the method of solution concern the computation of phase diagrams for ternary systems where azeotropes may be present and liquid-liquid phase split may occur. The algorithm described is related to the computation of one point of the phase diagram. For a ternary system, the isofugacity criteria for equilibrium is given by The fugacity,
e,for component i in phase k is given by = %“fP
where 23 indicates the composition of component i in phase k. The set of equations to be solved simultaneously therefore are the following:
Fli = x iL1$iL1 -xiL2$i L2 F2, = x:&
- yi$,”
~ 3 = 1 - ~ y - ~ y - ~ y ~4=1-~?2-~?-~? F5 = 1- y1-y2 - y3 F6 = x ; -~( x~y V~x y ) In the above equations, x and y indicate the component liquid and vapor mole numbers. Note that since the basis of the computations is 1 mol, mole fraction and mole number are the same. The above equation system shows that we have 10 equations and 12 unknown variables (T, P, xL1,xL2,y, xipec). Since there are three phases and there are three components, we have two degrees of freedom ( P or T and xipecare specified). Note that since we are computing the entire phase boundary, it is convenient to choose the pressure (or temperature) and one of the compositions as the specified variables. Note also that the equations for the fugacity coefficients are not shown since at a given T,P, and zk the fugacity coefficient of component i in phase k can be computed through an appropriate thermodynamic model (for example, any equation of state). In this work, the MHV2 model has been used.
Ind. Eng. Chem. Res., Vol. 32,No. 4, 1993 633 In order to employ the Newton-Raphson method, derivatives with respect to the unknown variables are needed. Expressions for these derivatives will be supplied to the reader upon request. The structure of the Jacobian matrix is shown below (*denotes a nonzero number and blank denotes a zero). It can also be noted that derivatives of the fugacity coefficients with respect to temperature and component mole numbers (in a given phase) are needed (which are available for the thermodynamic models used in this work). All the elements of the Jacobian matrix are evaluated analytically. x+' xi' 4' x y @ x3L.' y1 y2 y3 T(or P) * t * * F11 * * * F12 t * * * Fl3 * * * F21 * * F22 F23 F3 F4 F5 F6
* -1
-1
*
t
*
*
*
*
-1
-1
-1
*
-1
-1 -1 (0)
-1
-1
O(-1)
Pattern of the Jacobian matrix for VLLE computations (for function F6,one set of values are only used-either the values outside the parentheses or those inside the parentheses). Only the computation of the first point on the phase diagram is sensitive to the initial estimate. At all other points, the value from the previous point serves as a very good initial estimate. For the first point, since we start from the solvent-solute (A) binary line, the computation is not difficult unless we have a pressure close to the mixture critical point. Thus, for most calculations, this algorithm has been found to be robust and efficient. At the solution, the stability of the phases are checked to ensure that all three phases are stable.
Literature Cited Andersen, H. W.; Laroche, L.; Morari, M. Effect of Design on the Control of HomogenousAzeotropicDistillation. Presented at the AIChE Annual Meeting, San Francisco, 1989; paper no. 27D. Dahl, S.; Fredenslund, Aa.; Rasmussen, P. The MHVB model-A UNIFAC Based Model for Prediction of Gas Solubility and VaporLiquid Equilibria at Low and High Pressures. Znd. Eng. Chem. Res. 1991,30, 1936-1944. Dahl, S.;Dunalewica,A.; Fredenslund, Aa.; Rasmussen, P. The MHVB Model: Prediction of Phase Equilibria at Sub- and Supercritical Conditions. J. Supercrit. Fluids 1992,5, 42-47. Doherty, M. F.; Perkins, J. D. On the Dynamics of Distillation Processes-I. The Simple Distillation of Multicomponent Nonreacting, Homogenous Liquid Mixtures. Chem. Eng. Sci. 1978, 33,281-301.
Foucher, E. R.; Doherty, M. F.; Malone, M. F. Automatic Screening of Entrainers for HomogeneousAzeotropicDistillation. Znd.Eng. Chem. Res. 1991,30,760-772. Fredenslund, Aa. UNIFAC and Related Group-Contribution Prediction of Phase Equilibria. Fluid Phase Equilib. 1989,52,135150. Gani, R.; and Cameron, I. T. Modelling for Dynamic Simulation of ChemicalProcesses-The INDEX Problem. Chem.Eng. Sci. 1992, 45 (5), 1311-1315. Gani, R.; Ruiz, C. A.; Cameron, I. T. A Generalized Dynamic Model for Distillation Columns-I Model Description and Applications. Comput. Chem. Eng. 1986,10,181-198. Gmehling, J.; Onken, U.; Arlt, W. Vapor-Liquid Equilibrium Data Collection. Aqueous-Organic Systems. Chemistry Data Series; DECHEMA Frankfurt, 1981; Vol. 1, Part la, pp 572-577. Kovach 111, J. W.; Seider, W. D. Heterogenous Azeotropic Distillation: Experimental and Simulation Resulta. AZChE J. 1987a,33, 1300-1314. Kovach 111, J. W.; Seider, W. D. Heterogenous Azeotropic Distillation: Homotopy-Continuation Method. Comput. Chem. Eng. 1987b, 11, 593-605. Laroche, L.; Andersen, H. W.; Morari, M. Selecting the Right Entrainer for Homogenous Azeotropic Distillation. Presented at the AIChE Annual Meeting, Chicago, 1990; paper no. 215d. Magnussen, T. M.; Michelsen, M. L.; Fredenslund, Aa. Azeotropic Distillation Using UNIFAC. Znst. Chem. Eng. Symp. Ser. 1979, NO.56, 4.211-42/19. Michelsen, M. L. The Isothermal Flash Problem. Part I. Stability. Fluid Phase Equilib. 1982,9, 1-9. Pedersen, B. S. Static and Dynamic Behaviour of Multicomponent Distillation Columns. M.Sc. Dissertation, The Technical University of Denmark, Lyngby, Denmark, 1990. Perregaard, J.; Pedersen, B. S.; Gani, R. Steady State and Dynamic Simulation of Complex Chemical Processes. Trans Znst. Chem. Eng. 1992, 70,99-109. Pham, H. N.; Doherty, M. F. Design and Synthesis of Heterogeneous AzeotropicDistillations-I. Column Sequences. Chem.Eng. Sci. 1990a, 45 (7), 1823-1836. Pham, H. N.; Doherty, M. F. Design and Synthesis of Heterogeneous Azeotropic Distillations-111. Column Sequences. Chem. Eng. Sci. 1990b, 45 (7), 1845-1854. Prokopakis, G. J.;Seider, W. D. Feasible Specifications in Azeotropic Distillation. AZChE J. 1983a, 29, 49-60. Prokopakis, G. J.; Seider, W. D. Dynamic Simulation of Azeotropic Distillation. AZChE J. 1983b, 29, 1017-1029. Rovaglio, M.; Doherty, M. F. Dynamics of Heterogenous Azeotropic Distillation Columns. AZChE J. 1990, 36 (l),39-52. Ssrensen, E. L. Dynamic Simulation of Chemical Processes with Accurate Thermodynamic Models. Ph.D. Dissertation, The Technical University of Denmark, Lyngby, Denmark, 1991. Van Dongen, D. B.; Doherty, M. F. Design and Synthesis of HomogeneousAzeotropic Distillations. I. Problem Formulation for a Single Column. Ind. Eng. Chem. Fundam. 1985, 24,454463. Wahnschafft, 0. M.; Koehler, J. W.; Blass, E.; Westerberg, A. W. The Product Composition Regions of Single-Feed Azeotropic Distillation Columns. Znd. Eng. Chem. Res. 1992,31,2345-2344.
Received for review November 23, 1992 Accepted December 21, 1992
