Ind. Eng. Chem. Res. 2005, 44, 6845-6855
6845
A Method for Steady-State Simulation of Multistage Three-Phase Separation Columns Rahman Khaledi and P. R. Bishnoi* Department of Chemical and Petroleum Engineering, The University of Calgary, 2500 University Drive NW, Calgary, Alberta, Canada T2N 1N4
The development of a new algorithm for steady-state modeling of multistage three-phase separation columns is presented. In the developed algorithm, prior knowledge of phase pattern in the column is not required. The phase pattern determination and equilibrium calculations are carried out simultaneously. The stability of a phase is determined by coupling a stability equation with the model equations of the column. The resulting nonlinear model equations of the system are solved using an “inside-out” method. The ability and flexibility of the algorithm for the simulation of three-phase separation columns are shown by solving a variety of threephase distillation column examples. The algorithm is also applicable to two-phase separation columns. Introduction Three-phase distillation is a separation process in which vapor-liquid-liquid phases coexist in the distillation column. Three-phase distillation columns are common in chemical engineering processes. It is not unusual to have liquid-phase splitting on one or more stages within the distillation column where a highly nonideal feedstock is being processed. The formation of a second liquid phase in the column may enhance the separation process. In some cases, separation of the multicomponent mixture is either difficult or impossible by conventional distillation. Examples of such cases are separation of close-boiling mixtures and azeotropic mixtures. Heterogeneous azeotropic distillation enhances the separation process of such mixtures.1 An example of heterogeneous azeotropic distillation is separation of ethanol-water binary azeotropic mixture using benzene as entrainer.1 On the other hand, the formation of a second liquid phase may cause some operational problems. The presence of the second liquid phase leads to lower tray efficiency due to less efficient contact among the phases. There is also the possibility of cyclic buildup of the heavier liquid phase in the inactive regions on the trays due to settling of the heavier liquid phase.2 An example of such a case is a 1-butanol dehydration process in which a water-rich liquid phase forms toward the bottom of the column where the concentration of 1-butanol is higher.3 In this case, side withdrawal of the second liquid phase can prevent the above-mentioned operational problems. Therefore, simulation of three-phase distillation columns has a great importance for design engineers. It enables them to predict the behavior of three-phase distillation columns for their design and optimization. The simulation of three-phase distillation columns is more difficult than the simulation of two-phase columns. This difficulty is due to the fact that the number of phases present on each stage is not known until the calculations are performed. * To whom correspondence should be addressed. Fax: +(403) 282-3945. Tel: +(403)220-6695. E-mail: bishnoi@ ucalgary.ca.
In dealing with the presence of two liquid phases, there are three general approaches for simulating threephase column: (1) prespecified phase pattern methods; (2) methods involving a stability test for the stages in each calculation iteration; (3) simultaneous phase pattern identification and phase equilibrium calculations methods with specified maximum number of expected phases. In the first approach, the number of phases present at equilibrium on each stage is assumed “a priori”, and the model equations are solved with this predefined phase pattern. In these methods, if the assumed number of phases is less than what is actually present in the system, the method may converge to a pseudosolution.4 The appearance and disappearance of the second liquid phase may also lead to instability problems in the solution method. Most of the recent methods in the literature are based on the second approach. In this approach, first a phase pattern is assumed on each stage of the column. Subsequently, a phase stability test is performed for each stage and at each outer loop iteration to ensure the correctness of the assumed phase pattern. The process of phase stability test and search for other possible phases in these methods make the procedure computationally time-consuming. These methods also suffer from instability when the appearance or disappearance of a phase in the inner loop is encountered. Cairns and Furzer5 proposed a Naphtali-Sandholm6 based algorithm for the simulation of three-phase distillation columns. They formulated their model equations based on a single liquid phase with an average liquid mole fractions. The Michelsen7 phase stability method was performed on each stage to decide if the liquidphase splitting occurs. Georgoulaki and Korchinsky8 utilized a mixed K-value approach to handle the two liquid phases as a single liquid phase. They used the stability method developed by Pham and Doherty9 to determine whether the liquid phase on each stage is stable or not. In the third approach, the phase stability and equilibrium calculations are performed simultaneously for the system. Shyamsundar and Rangaiah10 proposed a method for the simulation and optimization of three-
10.1021/ie050088v CCC: $30.25 © 2005 American Chemical Society Published on Web 07/20/2005
6846
Ind. Eng. Chem. Res., Vol. 44, No. 17, 2005
phase distillation columns. Their method is based on the simultaneous phase stability and equilibrium calculation method developed by Han and Rangaiah,11 known as the τ-method. In this method, a phase characteristic variable, τ, is introduced into the mole fraction summation of each phase that allows the mole fraction summation to be less than unity for the nonexistent phase. Their approach for the solution of the model equations includes simultaneous solution of the governing equations and optimization for the phase pattern identification. In the present work, a new method for the simulation of multistage multiphase separation processes is developed. In this method, a prior knowledge of the phase pattern is not required and the stability of the phases is determined simultaneously during the calculations. This is accomplished by including a phase stability equation, developed by Gupta et al.,12 in governing equations. By simultaneous solution of these equations, the equilibrium composition and phase pattern are determined. The method is illustrated by presenting results for five three-phase distillation columns. Phase Stability Method. Gupta et al.12 developed a method for simultaneous phase equilibrium and stability calculations for multiphase reacting and nonreacting systems. This method is based on the minimization of the total Gibbs free energy of a system subject to material balance and nonnegative moles constraints. Solution of the minimization problem results in a nonlinear equation that describe the stability of a phase, i.e.,
R k θk ) 0
(k ) 1,...,π; k * r)
(1)
subject to
Rk g 0
and
θk g 0
(2)
where Rk is the phase fraction and θk is the phase stability variable. The phase stability variable, θk, is related to the fugacities of component i in phase k and reference phase r by the following equation:
θk ) ln(fˆi,k/fˆi,r)
(i ) 1,...,Nc; k ) 1,...,π; k * r)
(3)
It is noted that, at the solution point, θk is independent of the component i in the phase k. Equation 1 is used to describe the stability condition of a phase. It states that, in a system at equilibrium, for any phase k, either the phase fraction, Rk, or the phase stability variable, θk, is zero, or both are zero. In other words, during the phase equilibrium and stability computations, when the phase fraction, Rk, becomes zero, the stability variable, θk, for that phase should become either a positive value for an unstable phase or zero for an incipient phase to satisfy the phase stability equation (eq 1). Similarly, the phase fraction, Rk, should be either a positive value for a stable phase or zero for an incipient phase when the phase stability variable is zero to satisfy the phase stability equation (eq 1). Thus, the phase stability equation conveniently allows the appearance and disappearance of phases during the calculations. Equation 1 is modified to eq 4 to overcome some numerical difficulties that are discussed in detail by Gupta et al.12
Rkθk )0 Rk + θk
(k ) 1,...,π; k * r)
(4)
Figure 1. Schematic of a stage in multiphase column (k ) 1, vapor; k ) 2, ...π, liquid).
In the next section, the model equations for multiphase multistage separation processes are developed. Equation 4 is included in the model equations to handle the phase pattern determination simultaneously during the computations by specifying the maximum number of expected phases on each stage. Model Equations for Multiphase Distillation Columns. Figure 1 shows the schematic of a typical stage in a multiphase multicomponent separation column. The process is steady state with no chemical reaction involved. The phases are completely separated from each other and are at thermodynamic equilibrium when leaving the stage. A vapor and two liquid side streams may be withdrawn from each stage. The stage number increases from the top to the bottom of the column. The first stage is a combination of condenser and decanter, and the last stage is a reboiler. If the condenser does not exist, as in the case of stripping or absorption columns, the first stage would be the first tray in the column. Similarly, for absorption or stripping columns, where the reboiler does not exist, Ns would be the last tray in the column. If on any stage, side withdrawal, feed, or energy input do not exist, it is made equal to zero. The governing equations for the equilibrium stage are the component mass balance, energy balance, mole fraction summation, phase fraction summation, and thermodynamic phase stability and equilibrium equations. These equations are written as discrepancy equations as follow:
material balance equation π
δijM ) Fjzij -
∑ nijk + (1 - wrj+1,1)ni,j+1,1 +
k)1 π
∑ (1 - wrj-1,k)ni,j-1,k k)2
(i ) 1,...,Nc; j ) 1,...,Ns) (5)
energy balance equation Nc π
δjE ) Hfj + Qj -
∑ ∑ nijkhh ijk + i)1k)1 Nc
ni,j+1,1h h i,j+1,1 + ∑ i)1
(1 - wrj+1,1) Nc π
∑ ∑ (1 - wrj-1,k)ni,j-1,khh i,j-1,k i)1k)2
(j ) 1,...,Ns) (6)
In eqs 5 and 6, the terms involving j - 1 and j + 1, are dropped from the equations for the first and the last stage in the column.
Ind. Eng. Chem. Res., Vol. 44, No. 17, 2005 6847
mole fraction summation equation Nc
δjkS )
(xijk - xijr) (k ) 1,...,π; k * r; j ) 1,...,Ns) ∑ i)1
(7)
In eq 7, xijr refers to mole fraction in the reference phase. Since the vapor phase always exists in the distillation, absorption, or reboiled absorption columns, it is taken as reference phase.
zero with the disappearance of some of the phases and does not cause numerical problem during the calculations. The component molar flow rate, nijk, is related to mij by introducing variable σijk as
nijk ) σijkmij where
σijk )
phase fraction summation equation δj
)
∑ Rjk - 1
(j ) 1,...,Ns)
(8)
k)1
phase stability equation The modified form of phase stability equation (eq 4) suggested by Gupta et al.12 is written for stage j in the form of a discrepancy equation as,
δjkSt )
Rjkθjk (j ) 1,...,Ns; k ) 1,...,π; k * r) Rjk + θjk
Nc
Rjk )
(k ) 1,...,π; k * r)
Nc π
(15)
or Nc
Rjk )
(11)
Equation 11 gives the composition of all stable and unstable phases. For a stable phase, the phase stability variable, θk, is zero and eq 11 reduces to the well-known equilibrium relation. When the phase stability variable, θk, is a positive value, eq 11 gives the composition of a phase that is unstable (i.e., absent) and is not at equilibrium with the other existing phases. The model equations can be expressed in terms of mij, total moles of component i leaving the stage j,
nijk ∑ i)1 Nc
(16)
mlj ∑ l)1
(10)
The relation of the phase stability variable with the fugacities of components in mixture in different phases is given by eq 3. The fugacity coefficients in eq 10 can be expressed in terms of mole fractions and fugacities. By combining eq 3 with eq 10, the following relation for phase equilibrium is obtained:
xijk ) xijrKijkeθjk
nijk ∑ i)1
∑ ∑ nljp l)1p)1
The K-value is defined as the ratio of fugacity coefficients of reference phase and that of phase k by the following relation:
(k ) 1,...,π; k * r)
(14)
The variable σijk represents the fraction of component i leaving stage j through phase k and may be called the component phase fraction. Expressing the governing equations in terms of mij and σijk simplifies the form of equations and their derivatives. The phase fraction Rjk and mole fraction xijk are related to the component molar flow with the following relations:
It is noted that the reference phase always exists and the phase stability variable for this phase is always equal to zero. Hence, the phase stability equation for this phase is excluded from model equations.
φˆ ijr Kijk ) φˆ ijk
π
k)1
(9)
thermodynamic phase equilibrium equation
nijk
∑ nijk
π
PS
(13)
xijk )
nijk Nc
(17)
nljk ∑ l)1 Nc Substituting the ∑i)1 nijk term from eq 16 into eq 17 gives
xijk )
nijk Nc
Rjk
(18)
mij ∑ i)1
or Nc
nijk ) xijkRjk
mlj ∑ l)1
(19)
π
mij )
∑ nijk
(12)
k)1
When a phase disappears, the component molar flow rate of that phase, nijk, becomes zero. The advantage of the variable mij is that even with zero component molar flow rate, nijk, the total moles of component i leaving the stage, mij, is not zero. Therefore, mij does not become
Equation 13 also can be substituted into eq 18 to give the mole fraction relation as
xijk )
σijkmij Nc
Rjk
mij ∑ i)1
(20)
6848
Ind. Eng. Chem. Res., Vol. 44, No. 17, 2005
Substitution of eq 11 into eq 19 and then substituting its result into eq 14 gives the component phase fraction as
RjkeθjkKijk
σijk )
(21)
π
∑ Rjpe
θjp
Kijp
p)1
Equation 21 implicitly incorporates the equilibrium relation given by eq 11. Equation 21 when substituted in eq 20 gives the phase composition as
mijKijkeθjk
xijk )
Nc
(22)
π
mlj ∑ Rjpe ∑ l)1 p)1
θjp
Kijp
In fact, eq 22 is a form of equilibrium eq 11 that is used for updating phase compositions in the outer loop as discussed later. Substitution of eq 13 into eqs 5 and 6 and also eq 20 into eq 7 forms a new set of governing equations given below:
material balance equation
∑ (1 - wrj-1,k)σi,j-1,kmij-1
k)2
(i ) 1,...,Nc; j ) 1,...,Ns) (23) energy balance equation
∑ ∑ σijkmijhh ijk + i)1k)1 σi,j+1,1mi,j+1h h ij+11 + ∑ i)1
(j ) 1,...,Ns) (24) mole fraction summation equation 1 Nc σijk σijr mij δjkS ) Nc Rjr i)1 Rjk mlj
)
(k ) 1,...,π; k * r; j ) 1,...,Ns) (25) phase fraction summation equation π
δj
)
∑ Rjk - 1
(j ) 1,...,Ns)
(26)
k)1
phase stability equation δjkSt )
Rjkθjk Rjk + θjk
∑ ∑σi,Ns,kmi,Ns - B
(30)
k)2i)1
∑ ∑ (1 - wrj-1,k)σi,j-1,kmi,j-1hh ij-1,k i)1k)2
PS
(29)
π Nc
Nc π
∑ l)1
δCond.E ) R1,1 ) 0
δReb.E )
Nc
(
(28)
For specified bottom product rate, the last stage (reboiler) energy balance equation is replaced by
Nc π
∑
∑ wr1,kR1,k - R1,1
where RVL is top product vapor-to-liquid ratio. In the case of total condenser, this ratio is zero and eq 28 reduces to
π
(1 - wrj+1,1)
π
δCond.E ) RVL
k)2
δijM ) -mij + Fij + (1 - wrj+1,1)σi,j+1,1mij+1 +
δjE ) Hfj + Qj -
liquid phases during computations. They are also general and can model separation processes such as distillation, stripping, absorption, and reboiled absorption. The total number of these equations is Ns(Nc + 2π), namely, Nc‚Ns material balance, Ns energy balance, Ns(π - 1) mole fraction summation, Ns phase fraction summation, and Ns(π - 1) phase stability. Considering specified feed condition, side heat duties, side withdrawal ratios, and known stage pressures, it is required to solve the governing equations for Ns stage temperatures, Tj, Nc‚Ns overall component stage flow rates, mij, Nsπ phase fractions, Rjk, and Ns(π - 1) phase stability variables, θjk. The total number of these unknown variables is Ns(Nc + 2π). It is noted that the phase stability variable, θjk, corresponding to the reference phase (vapor) is always zero. Typically, in the distillation columns the bottom product rate and reflux ratios are specified while condenser and reboiler duties are calculated. The energy balance equation then can be replaced with the new specification equation. For specified reflux ratios, the first stage (condenser) energy balance equation is replaced by
(k ) 1,...,π; j ) 1,...,Ns) (27)
In the above equations, σijk is calculated using the phase equilibrium relation given by eq 21. Equations 23-27 are model equations for the steady-state simulation of multicomponent, multiphase, and multistage separation processes. The advantage of the above set of equations is that it can handle two- and three-phase distillation and easily allows for appearance and disappearance of
Once the solution is found, the reboiler and condenser duties can be calculated using energy balance eq 24. Method of Solution. The flow diagram of the solution algorithm used to solve the system of equation is shown in Figure 2. The procedure can be compared with equation tearing methods in which the equations are divided in two sets and solved in an inner and an outer loop. In the inner loop, material balance, energy balance, mole fraction summation, phase fraction summation, and stability equations (eqs 23-27) are solved simultaneously. An approach similar to the “inside-out” method of Boston and Sullivan13 is adopted for solution of these governing equations. The thermodynamic properties, approximated as simple functions of temperature, are used in the inner loop for computations. When the inner loop is converged, the parameters of the approximate models are updated using the rigorous thermodynamic models in the outer loop and the calculations are repeated. The approximate thermodynamic model for the K-value and the partial molar residual enthalpy as simple functions of temperature are given by the following equations:
(
ln Kik(T) ) Aik - Bik
1 1 T Tb
h h ikR ) Cik + Dik(T - Tb)
)
(31) (32)
where Tb is a base temperature and subscript k repre-
Ind. Eng. Chem. Res., Vol. 44, No. 17, 2005 6849
the inverse of temperatures, 1/Tj, the phase stability variables, θjk, and the phase fractions, Rjk. For all the stages, the independent variables may be represented in the vector form as
X ) {x1,x2,...,xj,...,xNs}T
(39)
where
xj )
{
}
1 m1,j,m2,j,...,mNc,j, ,θj,2,...,θj,k,...,θj,π,Rj,1,...,Rj,k,...,Rj,π Tj (40)
Figure 2. Flowchart of inside-out algorithm.
sents the phase number. It is assumed that the parameters of the approximate models, Aik, Bik, Cik, and Dik, are independent of temperature and composition and are kept constant in the inner loop. These parameters are related to rigorous thermodynamic properties by the following equations:
Aik ) lnKik(Tb)
)
(34)
h ikR(Tb) Cik ) h
(35)
Dik ) (∂h h ikR/∂T)T)Tb
(36)
Bik ) -
(
(33)
∂ ln Kik ∂(1/T)
T)Tb
It can be seen from eq 37 that the model equations are grouped by the stage in the Naphtali-Sandholm form.6 This grouping results in a block tridiagonal form for the Jacobian matrix of the model equations. In the inner loop, the set of Ns(Nc + 2π) independent nonlinear equations, F, is solved simultaneously for Ns(Nc + 2π) independent variables, X, using the Powell method14,15 or a combination of the Powell and homotopy continuation method.16,17 These methods provide a much larger domain of convergence than the Newton-Raphson method.18 To solve the model equations, first an attempt is made by using the Powell method. If the Powell method is unable to converge with the provided initial estimates, the homotopy continuation method, which is more robust, is used in the first few outer loop iterations to provide a better initial estimate for the Powell method. This procedure utilizes the robustness of homotopy continuation and the speed of the Powell method in solving the model equations. Once the inner loop is converged, the component compositions are updated in the outer loop using eq 22. With the updated compositions and temperatures, the parameters of approximate thermodynamic models are updated using the rigorous thermodynamic models in the outer loop. The inner loop is converged when the Euclidean norm of residual functions in eqs 23-27 is less that the specified tolerance, in. It is not necessary to converge the inner loop to a very small value of in in the initial outer loop iterations. Beginning from a relatively large value, the inner loop tolerance, in, can be gradually tightened as the outer loop iterations increase. This policy reduces the unnecessary inner loop computations in the initial outer loop iterations. An adjustable convergence tolerance given by eq 41 is used as inner loop convergence criteria.19
The set of inner loop nonlinear equations can be represented in the vector form as
in(i+1) ) in(i)/x10
F ) {f1,f2,...,fj,...,fNs}T
This equation is applied for every outer loop iteration provided that the number of inner loop iterations in the previous outer loop iteration is less than or equal to five. If the number of inner loop iterations is more than five, the in is kept constant until the number of inner loop iterations drops below six. In this work, the in value starts from 0.1 and is allowed to decrease to a minimum value of 10-6. The outer loop is considered to be converged when the convergence tolerance criteria out, which is based on composition change between two consecutive outer loop iterations, is less than or equal to the specified value of 10-8, that is,
(37)
where fj is the residual or discrepancy of equations for the jth stage, that is
fj ) {δ1,jM,δ2,jM,...,δNc,jM,δjE,δj,2S,...,δj,kS,...,δj,πS, δjPS,δj,2St,...,δj,kSt,...,δj,πSt}T (38) The independent variables in the inner loop are grouped by stage. For each stage, the independent variables consist of the overall component molar flow rates, mij,
(41)
6850
Ind. Eng. Chem. Res., Vol. 44, No. 17, 2005
x∑∑∑ Nc Ns π
out )
(∆xijk)
i)1 j)1k)1
Nc‚Ns‚π
2
e 10-8
(42)
Initialization Method. Two-phase distillation columns are initialized by assuming a linear stage temperature profile between the condenser bubble point and reboiler dew point temperatures evaluated at average composition of all feeds introduced to the column. A linear pressure profile is considered in the column. Vapor- and liquid-phase fractions are initialized by assuming constant molal overflow in the column, and initial phase stability variables for all phases are set equal to zero. The initial values for component overall flow rates, mij, are calculated by solving the tridiagonal system of linear equations formed by component mass balance, eq 23. The component phase fractions, σijk, in eq 23 are calculated using eq 21 with ideal equilibrium ratios (K-value) evaluated at stage temperatures and pressures. The initial phase compositions then are calculated from eq 22. Subsequently, the parameters for the approximate thermodynamic models are calculated at the initial stage temperatures and compositions using rigorous thermodynamic models. The three-phase column is initialized first by performing one iteration of two-phase distillation calculations. Subsequently, a liquid-liquid flash calculation is performed on the generated liquid on each stage. The initial K-values for liquid-liquid flash calculation are generated using a noniterative stability analysis method.20 The resulting compositions and phase fractions from the liquid-liquid flash calculation are utilized as initial values for the two liquid phases in the column. If the performed liquid-liquid flash calculations do not find a second liquid on any of the stages, composition and phase fraction of a nearest stage that has two liquid phases is used to initialize composition and phase fraction of this stage. These compositions are used to initialize the parameters of approximate thermodynamic models for the first outer loop iteration in the three-phase column calculations. Simulation Results. Although the developed algorithm was used to simulate a number of two- and threephase separation columns,21 only a selected number of three-phase distillation examples are reported here. Table 1 gives column specifications for five selected
three-phase distillation columns. For all the examples, the vapor-phase thermodynamic properties are calculated using the Peng-Robinson22 equation of state. The ideal gas heat capacities and the standard state fugacities are taken from Prausnitz et al.23 The feed and operating pressures for all the three-phase distillation examples, given in Table 1, are 1 atm. The simulation results for these examples are discussed in detail to illustrate the ability and performance of the algorithm in handling different situations in three-phase separation processes. The computer program for this algorithm is written in the C++ language. The compiler used is Visual C++ version 5.0. The program is run on a personal computer with a Pentium 4 CPU and a 2.4-GHz processor. The required numbers of outer loop iterations and CPU times to converge the algorithm for the examples are given in Table 2. For some of the examples, the calculation is started with the homotopy method. As can be seen from the table, the algorithm was able to converge in a reasonable number of iterations and in a short computational time for all the examples. Example 1: 1-Butanol-Water-1-Propanol. Example 1 is separation of 1-butanol, water, and 1-propanol mixture. This example was first studied by Block and Hegner3 and then tested by many other researchers as a benchmark example.2,4,10,24-26 The mixture contains two partially immiscible componentsswater and 1-butanol. 1-Propanol is soluble in the other two components and distributes between the immiscible phases. The liquid phase splits into organic and aqueous liquid phases in the stripping section of the column where the concentration of 1-butanol is increasing toward the bottom stages. For this problem, the NRTL activity model with parameters taken from Gmehling and Onken27 is used to predict the liquid-phase thermodynamic properties. The simulation results are presented in Figures 3-6 and Table 3. A good agreement is found between the results obtained in this work and those obtained by Block and Hegner.3 The source of thermodynamic data used by Block and Hegner3 is not reported. Hence, slight differences in results are likely due to the use of different thermodynamic model parameters. A single liquid phase is found on stages 1-7, while two liquid phases (aqueous and organic phases) are present on stages 8-12. As it can be seen from Figure 7, the stability variable for the second liquid phase on
Table 1. Details of Column Specifications for Selected Three-Phase Distillation Examplesa feed conditions example
components
1
1-butanol water 1-propanol acetone chloroform water acetonitrile trichloroethylene water 1-butanol water 1-butyl acetate ethanol benzene water
2 3 4 5
a
compositionb FI FII 0.13 0.65 0.22 0.167 0.333 0.5 0.45 0.52 0.029 0.24 0.3 0.46 0.2184 0.7419 0.0397
0.6 0.2 0.2 0.6848 0.0002 0.3150
0.89 0.0 0.11
column specifications rate (kmol/h)
condition
no. of stage
F5: 50
liq at 363.15 K
12
R1 ) 3.0 R2 ) 3.0
21
F4: 60 F5: 100
sat vap sat liq
11
R1 ) 2.5 R2 ) 2.5
120
F2: 104.9 F5: 100
sat liq sat liq
13
R1 ) 0.0(aq) R2 ) 1.5
66.1
F2: 50
sat liq
7
R1 ) total R2 ) 0.0 (aq)
34.84
F1: 45.32 F5: 100
sat liq sat liq
28
R1 ) 0.0 R2 ) 4.0 (org)
70
Feed and column pressure, 1 atm. b FI, first feed from top; FII, second feed from top.
reflux ratio
B (kmol/h)
Ind. Eng. Chem. Res., Vol. 44, No. 17, 2005 6851
Figure 3. Stage temperature profile in example 1 (line, this work; points, ref 3).
Figure 4. Vapor and liquid flow rates in example 1 (line, this work; points, ref 3). Table 2. Number of Iterations and Computational Time for Examples of Table 1
Figure 5. Vapor composition profile in example 1 (line, this work; points, ref 3).
Figure 6. Average liquid composition profile in example 1 (line, this work; points, ref 3).
outer loop iteration example
homotopy
Powell
total
CPU time (s)
1 2 3 4 5
1 0 1 2 0
9 26 51 7 12
10 26 52 9 12
3.5 3.6 7.9 2.0 3.3
stages 1-7 is a positive value, which means that the second liquid phase cannot exist on these stages. The phase fraction of the second liquid phase is zero for these stages. On the other hand, for stages 8-12, the stability variable is zero and the phase fraction on these stages is a positive value. Figure 8 shows variation of stability variable and phase fraction for the second liquid phase at the eighth stage during the calculations. The outer loop iteration is initialized with positive phase fractions and zero stability variables for the second liquid phase according to the method discussed previously. As can be seen from Figure 8, the phase fraction of the second liquid phase is a positive value and the phase exists until the fourth iteration, when the phase disappears. The phase stability factor in the fourth iteration is a small positive value that indicates the absence of the phase. The phase fraction remains zero and the phase stability variable remains a positive value until the eighth iteration, when the phase fraction becomes a positive value and the phase stability variable becomes zero. This indicates appearance of the second liquid phase on the stage. The variation of the phase fraction and the phase stability variable at the eighth stage illustrates the ability of the algorithm in allowing for appearance and disappearance of a phase during the computations.
Figure 7. Phase fraction and stability variable profiles of liquid II (aqueous) in example 1.
In the proposed algorithm, the composition of an unstable phase is calculated along with the composition of the stable phases in each iteration. Although the phase fraction of an unstable phase is zero, its composition is used to calculate its thermodynamic properties. Thus, the computation for a phase can move smoothly from its unstable state to a stable state without the need for composition reinitialization and reconstruction of model equations. A similar transition takes place when a stable phase becomes unstable during the computations. Stage compositions for both liquid phases are compared with the results obtained by Block and Hegner3 and given in Table 3. This example clearly demonstrates the ability of the proposed algorithm in handling appearance and disappearance of phases during computations. Example 2: Acetone-Chloroform-Water. Example 2 is distillation of acetone, chloroform, and water
6852
Ind. Eng. Chem. Res., Vol. 44, No. 17, 2005
Table 3. Stage Compositions for Liquid Phases in Example 1 (xI, Liquid I; xII, Liquid II Composition) stage no.
liquid compsn
1
xI xII xI xII xI xII xI xII xI xII xI xII xI xII xI xII xI xII xI xII xI xII xI xII
2 3 4 5 6 7 8 9 10 11 12
1-butanol this work ref 3
this work
water ref 3
1-propanol this work ref 3
0.029 483
0.030 840
0.611 181
0.619 607
0.359 336
0.349 552
0.058 348
0.057 763
0.643 561
0.647 114
0.298 091
0.295 123
0.087 781
0.084 724
0.674 330
0.672 762
0.237 890
0.242 514
0.112 134
0.107 482
0.695 780
0.690 801
0.192 086
0.201 717
0.135 812
0.130 386
0.695 832
0.689 119
0.168 356
0.180 495
0.152 878
0.145 868
0.711 568
0.702 242
0.135 553
0.151 890
0.169 493
0.161 569
0.723 130
0.712 946
0.107 377
0.125 485
0.201 068 0.025 721 0.230 953 0.024 279 0.257 495 0.023 507 0.280 356 0.022 944 0.299 643 0.022 475
0.187 443 0.025 090 0.216 519 0.023 784 0.243 862 0.022 770 0.268 828 0.021 992 0.291 195 0.021 393
0.709 388 0.951 956 0.696 736 0.959 885 0.686 077 0.965 217 0.677 234 0.969 178 0.669 972 0.972 213
0.705 573 0.947 484 0.693 227 0.955 993 0.682 274 0.962 448 0.672 697 0.967 355 0.664 390 0.971 128
0.089 543 0.022 323 0.072 311 0.015 836 0.056 428 0.011 276 0.042 410 0.007 878 0.030 385 0.005 312
0.106 984 0.027 422 0.090 254 0.020 223 0.073 864 0.014 782 0.058 475 0.010 650 0.044 416 0.007 478
system. This problem was previously studied by Boston and Shah,28 Cairns and Furzer,26 Georgoulaki and Korchinsky,8 and Shyamsundar and Rangaiah.10 The UNIQUAC model with parameters taken from Sørensen and Arlt29 is used to model the liquid-phase thermodynamic properties in the present work. Simulation results show that the phase stability variable for the second liquid phase is zero on all stages and two immiscible liquid phases are found on all stages. Similar results were found by other authors. The results are
shown in Figures 9-12. It can be seen from Figure 9 that two liquid phases are present on all the stages. Figures 10 and 11 present the mole fraction profiles for the first (organic) and second (aqueous) liquid phases, respectively. Average liquid mole fractions are also plotted in Figure 12. A good agreement is found between average liquid mole fraction obtained in this work and those obtained by Shyamsundar and Rangaiah,10 except for a slight difference (less than 0.02) in acetone and chloroform average mole fraction on the first stage.
Figure 8. Variation of phase fraction and stability variable of liquid II (aqueous) with iteration number at eighth stage in example 1.
Figure 10. Organic-phase (liquid I) composition profile in example 2.
Figure 9. Phase fraction profile in example 2.
Figure 11. Aqueous-phase (liquid II) composition profile in example 2.
Ind. Eng. Chem. Res., Vol. 44, No. 17, 2005 6853
Figure 12. Average liquid composition profile in example 2 (line, this work; points, ref 10). Figure 14. Composition profile of liquid I in example 4.
Figure 13. Liquid-phase composition profiles in example 3 (lines, liquid I; points, liquid II).
Example 3: Acetonitrile-TrichloroethyleneWater. This example is heterogeneous azeotropic distillation of acetonitrile and water mixture using trichloroethylene as entrainer. This system originally was studied by Pratt30 and then simulated by Cairns and Furzer.26 The UNIQUAC model with interaction parameters taken from Sørensen and Arlt29 is used to model the liquid-phase thermodynamic properties. Two liquid phases are found on the first stage (condenserdecanter). The aqueous phase consists mainly of water and is withdrawn completely from the column, and the organic phase is refluxed partially to the column. A highly concentrated acetonitrile stream is taken as product from the bottom of the column. Figure 13 presents the liquid mole fraction profiles for this example. In some cases during three-phase distillation calculations, the composition of the second liquid phase approaches the composition of the first liquid phase, and consequently, the K-values for both liquid phases become identical. In this case, mole fraction summation equations of the two liquids become similar in the systems of equations and results in singular Jacobian. To avoid the problem of Jacobian singularity, the two liquid phases are collapsed by adding up the phase fraction of the two liquid phases and the phase number of that stage reduces to two. After the three-phase distillation calculation is converged, a phase stability test is performed to ensure the correctness of number of phases on that stage. If a second liquid is found, the calculation is resumed with the composition found from the stability test as the initial estimate for the second liquid phase. However, for all the examples studied in this work, none of the collapsed phases was found to
Figure 15. Liquid-phase composition profiles in example 5 for solution S1, the high ethanol concentration product (lines, liquid I; points, liquid II).
have more than one liquid phase present after the stability test. Example 3 is one of the cases in which for all the stages except for the first stage (where two liquid phases are present) the algorithm found identical liquid-phase compositions and collapsed the two liquids. Example 4: 1-Butanol-Water-1-Butyl Acetate. This example is dehydration of 1-butanol and 1-butyl acetate mixture that was studied by Block and Hegner,3 Cairns and Furzer,26 and recently by Shyamsundar and Rangaiah.10 The NRTL activity model with interaction coefficients taken from Gmehling and Onken27 is used to predict the liquid thermodynamic properties in this work. Two liquid phases are present only on the first stage (condenser-decanter). The first liquid phase (organic) is totally refluxed to the column while the second liquid phase (aqueous) is withdrawn completely from the decanter. The aqueous phase withdrawn from the decanter contains mainly water (1-butanol, 0.0089 mol %; water, 0.9894 mol %; 1-butyl acetate, 0.0017 mol %) whereas bottom product contains virtually no water. Mole fraction profile for the first liquid is plotted in Figure 14. Example 5: Ethanol-Benzene-Water. This example is the heterogeneous azeotropic distillation of ethanol and water mixture using benzene as entrainer. This system was studied by Magnussen et al.,31 Prokopakis and Seider,32 Kovach and Seider,33 and Georgoulaki and Korchinsky.8 The UNIQUAC model with interaction parameters taken from Gmehling and Onken27 reported by Prokopakis and Seider32 is used to calculate the liquid-phase thermodynamic properties.
6854
Ind. Eng. Chem. Res., Vol. 44, No. 17, 2005
Figure 16. Liquid-phase composition profiles in example 5 for solution S2, the low ethanol concentration product (lines, liquid I; points, liquid II).
Figure 18. Liquid-phase initial composition profiles for example 5 (lines, IS1, initial composition for solution S1; lines with points, IS2, initial composition for solution S2).
profiles for the same column specifications are obtained for the S1 and S2 solutions. The initial temperature profile that results in solution S2, is lower than the initial temperature profile that results in solution S1. Similarly, the initial composition profile for the ethanol that results in solution S1 is higher than the initial composition profile that results in solution S2, as shown in Figure 18. Conclusions Figure 17. Stage temperature profile in example 5 (S1, temperature profile for solution S1; S2, temperature profile for solution S2; IS1, initial temperature profile for solution S1; IS2, initial temperature profile for solution S2).
The simulation results for this example show liquidphase splitting on the first stage (condenser-decanter). The aqueous phase consists mainly of water and is withdrawn completely from the column, and the organic phase contains mainly benzene and is partially refluxed to the column. A highly concentrated ethanol stream is taken as product from the bottom of the column. The concentration profile is shown in Figure 15. For this system, Magnussen et al.31 reported multiple steady-state solutions. Other authors32,33 also observed multiplicity in solutions for this system. We found two steady-state solutions for the specifications given in Table 1. These solutions are similar to the solutions found by Magnussen et al.31 The first solution denoted by S1 gives a high ethanol recovery with bottom product containing a high concentration of ethanol. The second solution, denoted by S2, gives a lower concentration of ethanol in the bottom product (∼93 mol %). The concentration profiles for S1 and S2 solutions are presented in Figures 15 and 16, respectively. The first solution, S1, was found simply by solving the problem for the specifications given in Table 1 for this example. The second solution, S2, was found first by solving the column with bottom product flow rate of B ) 100 kmol/h (instead of the specified value of B ) 70 kmol/h) and then changing the bottom product flow rate to its specified value of B ) 70 kmol/h. In this way, the initial composition and temperature profile for starting the computations would be different from the initial temperature and concentration profiles that result in solution S1. The solution and initial temperature profiles for the cases S1 and S2 is given in Figure 17. As can be seen from this figure, two different final temperature
An algorithm for simulating three-phase distillation columns and two-phase multistage separation processes was developed. In this algorithm, prior knowledge of phase pattern in the column is not required and the phase pattern determination and equilibrium calculation are performed simultaneously. An “inside-out” method with approximate thermodynamic properties model was adapted to solve the resulting nonlinear model equations of the system. The algorithm was successfully applied to simulate a variety of two-phase21 and three-phase distillation columns, and the computations converged very fast for all the examples. The simulation results were compared with literature data where available. The ability of the algorithm in handling the appearance and disappearance of phases during computations was well demonstrated. The extension of the algorithm for the simulation of two- and three-phase reactive distillation columns will be presented in a later publication. Acknowledgment The financial support for this work was provided by the National Science and Engineering Research Council (NSERC) of Canada. Nomenclature A ) parameter in eq 31 B ) parameter in eq 31 or bottom product rate (kmol/h) C ) parameter in eq 32 D ) parameter in eq 32 F ) model equations vector F ) feed molar rate (kmol/h) f ) model equation functions ˆf ) fugacity (kPa) h h ) partial molar enthalpy (kJ/kmol) Hf ) feed enthalpy (kJ) K ) equilibrium ratio
Ind. Eng. Chem. Res., Vol. 44, No. 17, 2005 6855 m ) total molar flow rate (kmol/h) n ) molar flow rate (kmol/h) Nc ) number of components Ns ) number of stages Q ) stage energy input (kJ) RVL ) vapor liquid ratio T ) temperature (K) wr ) side-withdrawal ratio x ) mole fraction or stage-independent variable vector X ) variables vector z ) feed mole fraction Greek Symbols R ) phase fraction δ ) discrepancy or function residual ) convergence tolerance φˆ ) fugacity coefficient π ) number of phases θ ) phase stability variable σ ) component phase fraction Subscripts and Superscripts i ) component index j ) stage number k ) phase index R ) residual properties
Literature Cited (1) Holland, C. D. Fundamentals of Multicomponent Distillation; McGraw-Hill Inc.: New York, 1981. (2) Ross, B. A.; Seider, W. D. Simulation of Three-Phase Distillation Towers. Comput. Chem. Eng. 1980, 5, 7. (3) Block, U.; Hegner, B. Development and Application of a Simulation Model for Three-Phase Distillation. AICHE J. 1976, 22, 582. (4) Buzzi Ferraris, G.; Morbidelli, M. Distillation Models for Two Partially Immiscible Liquids. AICHE J. 1981, 27, 881. (5) Cairns, B. P.; Furzer, I. A. Multicomponent Three-Phase Azeotropic Distillation. 1. Extensive Experimental Data and Simulation Results. Ind. Eng. Chem. Res. 1990, 29, 1349. (6) Naphtali, L. M.; Sandholm, D. P. Multicomponent Separation Calculations by Linearization. AICHE J. 1971, 17, 148. (7) Michelsen, M. The Isothermal Flash Problem. Part I. Stability. Fluid Phase Equilib. 1982, 9, 1. (8) Georgoulaki, A.; Korchinsky, W. J. Simulation of Heterogeneous Azeotropic Distillation: An Improved Algorithm Using Modified UNIFAC Thermodynamic Predictions and the NaphtaliSandholm Matrix Method. Trans. Inst. Chem. Eng. 1997, 75, 101. (9) Pham, H. N.; Doherty, M. F. Design and Synthesis of Heterogeneous Azeotropic Distillations. I. Heterogeneous Phase Diagrams. Chem. Eng. Sci. 1990, 45, 1823. (10) Shyamsundar, V.; Rangaiah, G. P. A Method for Simulation and Optimization of Multiphase Distillation. Comput. Chem. Eng. 2000, 24, 23. (11) Han, G.; Rangaiah, G. P. A Method for Multiphase Equilibrium Calculations. Comput. Chem. Eng. 1998, 22, 897. (12) Gupta, A. K.; Bishnoi, P. R.; Kalogerakis, N. A Method for the Simultaneous Phase Equilibria and Stability Calculations for Multiphase Reacting and Non-Reacting Systems. Fluid Phase Equilib. 1991, 63, 65. (13) Boston, J. F.; Sullivan, S. L. A New Class of Solution Methods for Multicomponent, Multistage Separation Processes. Can. J. Chem. Eng. 1974, 52, 52.
(14) Powell, M. J. D. A Hybrid Method for Non-Linear Equations. In Numerical Methods for Nonlinear Algebraic Equations; Rabinowitz, P., Ed.; Gordon and Breach: New York, 1970. (15) Garbow, B. S.; Hillstrom, K. E.; More, J. J. Hybrid Subroutine. Minpack Library, Netlib; Argonne National Laboratory, 1980. (16) Rheinboldt, W.; Burkardt, J. A Program for a Locally Parameterized Continuation Process. ACM Trans. Math. Soft. 1983, 9, 236. (17) Rheinboldt, W.; Burkardt, J. Pitcon 6.1 Subroutine. Conti Library, Netlib; University of Pittsburgh, 1991. (18) Wayburn, T. L.; Seader, J. D. Homotopy Continuation Methods for Computer-Aided Process Design. Comput. Chem. Eng. 1987, 11, 7. (19) Saeger, R. B.; Bishnoi, P. R. A Modified “Inside-Out” Algorithm for Simulation of Multistage Multicomponent Separation Processes Using the UNIFAC Group-Contribution Method. Can. J. Chem. Eng. 1986, 64, 759. (20) Trebble, M. A. A Preliminary Evaluation of Two and Three Phase Flash Initiation Procedures. Fluid Phase Equilib. 1989, 53, 113. (21) Khaledi, R. Simulation of Multiphase Reactive and NonReactive Separation Processes. Ph.D. Thesis, University of Calgary, Calgary, Alberta, Canada, 2005. (22) Peng, D. Y.; Robinson, D. B. A New Two-Constant Equation of State. Ind. Eng. Chem. Fundam. 1976, 15, 59. (23) Prausnitz, J. M.; Anderson, T. F.; Grens, E. A.; Eckerts, C. A.; Hsieh, R.; O’Connell, J. P. Computer Calculations for Multicomponent Vapor-Liquid and Liquid-Liquid Equilibria; Prentice-Hall: Englewood Cliffs, NJ, 1980. (24) Kinoshita, M.; Hashimoto, I.; Takamatsu, T. A Simulation Procedure for Multicomponent Distillation Column within which Three Phases of Vapor and Two Partially Immiscible Liquids are Present. J. Chem. Eng. Jpn. 1983, 16, 513. (25) Swartz, C. L. E.; Stewart, W. E. Finite-Element SteadyState Simulation of Multiphase Distillation. AICHE J. 1987, 33, 1977. (26) Cairns, B. P.; Furzer, I. A. Multicomponent Three-Phase Azeotropic Distillation. 3. Modern Thermodynamic Models and Multiple Solutions. Ind. Eng. Chem. Res. 1990, 29, 1383. (27) Gmehling, J.; Onken, U. Vapor-Liquid Equilibrium Data Collection; Chemistry Data Series Vol. I; Dechema: Flushing, NY, 1977; Part 1. (28) Boston, J. F.; Shah, V. B. An Algorithm for Rigorous Distillation Calculation with Two Liquid Phases. Paper Presented at AICHE Meeting, Houston, TX, April 1979. (29) Sørensen, J. M.; Arlt, W. Liquid-Liquid Equilibrium Data Collection (Ternary Systems); Chemistry Data Series Vol. V; Dechema: Flushing, NY, 1979; Part 2. (30) Pratt, H. R. C. Continuous Purification and Azeotropic Dehydration of Acetonitrile Produced by the Catalytic Acetic AcidAmmonia Reaction. Trans. Inst. Chem. Eng. 1947, 25, 43. (31) Magnussen, T.; Michelsen, M. L.; Fredenslund, A. Azeotropic Distillaton Using UNIFAC. Inst. Chem. Eng. Symp. Ser. No. 56, 3rd Inter. Symp. on Distillation, ICE, Rugby, Warwickshire, England, 1979. (32) Prokopakis, G. J.; Seider, W. D. Feasible Specifications in Azeotropic Distillation. AICHE J. 1983, 29, 49. (33) Kovach, J. W., III; Seider, W. D. Heterogeneous Azeotropic Distillation. Homotopy-Continuation Methods. Comput. Chem. Eng. 1987, 11, 593.
Received for review January 21, 2005 Revised manuscript received May 13, 2005 Accepted June 8, 2005 IE050088V
