Optimal Design of Staged Three-Phase Reactive Distillation Columns

Mar 4, 2010 - Therefore, an accurate but compact in size and thus easier to solve .... large set of nonlinear differential algebraic equations (DAEs) ...
1 downloads 0 Views 2MB Size
Ind. Eng. Chem. Res. 2010, 49, 3275–3285

3275

Optimal Design of Staged Three-Phase Reactive Distillation Columns Using Nonequilibrium and Orthogonal Collocation Models Theodoros Damartzis† and Panos Seferlis*,†,‡ Chemical Process Engineering Research Institute, Centre for Research and TechnologysHellas, P.O. Box 60361, 57001 Thermi-Thessaloniki, and Department of Mechanical Engineering, Aristotle UniVersity of Thessaloniki, P.O. Box 484, 54124 Thessaloniki, Greece

Reactive distillation with potential liquid-phase split and subsequent formation of a third phase is a highly complex process system nevertheless quite common in the production of useful solvents and biofuels through esterification. The optimal design of such process systems requires the development and solution of reliable and accurate process models that lead to a computationally demanding mathematical problem. In this work, a nonequilibrium (NEQ) model coupled with the orthogonal collocation on finite elements (OCFE) technique is developed for the simulation and optimal design of three-phase reactive distillation systems. The resulting NEQ/OCFE model combines the predictive accuracy of the NEQ model as well as the model reduction and approximation capabilities of the OCFE formulation. Therefore, an accurate but compact in size and thus easier to solve model that accounts for all the physical phenomena, the interactions among the multiple phases, and the occurring chemical reactions becomes available. The NEQ/OCFE model is enriched with an accurate prediction procedure for the identification and tracking of the phase transition boundaries (point of transition for a single liquid phase to two liquid-phase regimes and vice versa) inside the column for varying operating conditions. The model is validated using experimental results and utilized in the optimal design and the dynamic simulation of a staged reactive distillation column for the production of butyl acetate via the esterification reaction of butanol with acetic acid. The optimal column configuration defined by the number of stages in each column section, the feed strategy (single feed or multiple feed points), the location of the feed stages, and the operating conditions are calculated through a rigorous design optimization procedure for tight production purity specifications. At optimal conditions the specific column appears to have both twoand three-phase regions that are separated by the feed stage. 1. Introduction Reactive distillation (RD) is a relatively new and attractive enhancement of the classical distillation processes. The incorporation of chemical reactions and product separation into one single process may result in significant economical benefits due to a drastic reduction of investment costs and higher product yields when compared to conventional reaction distillation processes especially for chemical equilibrium controlled reactions with substantial volatility differences among reactants and products. Reactive distillation has been mainly used for the industrial production of ether and acetate solvents such as TAME (tert-amyl methyl ether), MTBE (methyl tert-butyl ether), and esters such as methyl, ethyl, and butyl acetate. The achieved process intensification results in higher overall efficiency, but in the meantime process complexity increases significantly. Indeed, reactive distillation is a complex unit operation, where simultaneous mass and energy transfer, coupled with phase equilibrium and reactive phenomena, occurs. The reactive separation system becomes even more perplexed when the formation of a second liquid phase is favored by thermodynamics in some sections of the distillation column. The splitting of a single liquid phase into two distinct liquid phases is a phenomenon that cannot be neglected in design, simulation, and control of the reactive separation system as it may greatly affect the product stream quality and subsequently the overall column performance. On several occasions, the addition of a * To whom correspondence should be addressed. Tel.: +30 2310 994229. E-mail: [email protected]. † Centre for Research and TechnologysHellas. ‡ Aristotle University of Thessaloniki.

component in the reactive mixture aims to create such thermodynamic conditions that favor the formation of a second liquid phase for the efficient separation of azeotropic mixtures. Phase transition phenomena greatly complicate the associated process modeling equations due to the higher degree of nonlinearity that is introduced. Furthermore, phase transition causes discontinuities in the model due to the spontaneous appearance and disappearance of the second liquid phase, thus increasing model complexity and consequently solution effort. Another challenge relates to the fact that the boundaries between the two-phase and three-phase regions in the column are not known a priori as these are solely determined by the prevailing conditions within the column that become available only after a solution has been obtained. This implies that an alternation of modeling relations that correspond to two- and three-phase regimes for successive column sections may be required. Several researchers have studied extensively the optimal design of three-phase reactive distillation systems. More specifically, Venimadhavan et al.1 and Steinigeweg and Gmehling2 using simplified process models proposed optimal column configurations and operating conditions for a high-purity butyl acetate column. Steyer et al.3 proposed a novel reactive distillation process for the production of cyclohexanol from cyclohexene and water taking into consideration liquid-phase splitting. Khaledi and Bishnoi4 and Khaledi and Bishnoi5 used an equilibrium model to describe three-phase reactive distillation and examined various reactive systems. Shyamsundar and Rangaiah6 used an equilibrium model coupled with a fictitious parameter (τ method) that enabled the tractability of phase transition boundaries through suitable model structure adjustments. Bruggemann et al.7 and Radulescu et al.8 investigated

10.1021/ie901260b  2010 American Chemical Society Published on Web 03/04/2010

3276

Ind. Eng. Chem. Res., Vol. 49, No. 7, 2010

the dynamic features of three-phase reactive distillation systems using equilibrium models. In addition to modeling, design optimization has also been examined for reactive distillation columns using equilibrium models. Gangadwala et al.9 and Gangadwala and Kienle10 examined various column configurations and case studies for the optimal design of butyl acetate production via catalytic two-phase reactive distillation using mixed-integer nonlinear programming. However, in the aforementioned modeling and design attempts equilibrium models have been employed that require all phases present in a column tray to be in thermodynamic equilibrium. While the equilibrium model provides greater simplicity, it does not provide a realistic representation of the phenomena occurring during reactive distillation, leading to discrepancies between simulated predictions and experimental data. Such a process can be described in a more accurate way by introducing a nonequilibrium (NEQ) or rate-based model.11,12 Lao and Taylor11 developed an NEQ model for a single distillation stage which was later expanded for an entire column by Higler et al.12 Higler et al.12 provide the most comprehensive NEQ model for conventional three-phase distillation, which was tested for its accuracy of predictions in a number of relevant separation systems. While these two works consider the general case of rate-based mass and energy transfer among all three phases, Eckert and Vanek13 used an NEQ model assuming that the two liquid phases are in equilibrium with each other. This last assumption is an interesting physical and mathematical simplification that applies in cases where strong agitation of the liquid phases is present. Evidently, the employment of NEQ models further increases the complexity and size of the resulting process model. Unfortunately, such highly complex models consist of a very large set of nonlinear differential algebraic equations (DAEs) that require significant computational effort usually forbidding real time applications. In addition, it is necessary to develop a technique that would enable the identification of the phase transition boundaries and also trace the phase discontinuities inside the column as the operating conditions change in design optimization problems or dynamic simulations. Such requirements necessitate that the developed process model possesses the flexibility to adapt in terms of model structure as phase transition points emerge or disappear within the column domain. Approximating techniques such as orthogonal collocation on finite elements (OCFE) provide a solid solution to the significant reduction in terms of the total number of modeling equations without compromising the resolution of modeling details. OCFE is a technique that transforms the column in a continuous analogue, thus eliminating the need for integer design variables associated with the existence of column trays, thus tackling the discrete character of staged distillation columns. The OCFE formulation approximates the temperature and composition profiles in staged or packed columns by continuous polynomial functions of position inside the column. Consequently, the OCFE formulation uses a smaller set of equations, proportional to the degree of the approximating polynomials, to describe the reactive distillation process than other conventional modeling methods and in this way reduces the computational effort for the solution of the model. Moreover, the approximation of the concentration and temperature profiles as continuous functions of position in the column for the OCFE formulation eliminates the need for the use of integer variables for the representation of stages as column size becomes a continuously varying variable. Therefore, the optimal design of staged separation units

can be formulated as a conventional nonlinear program, thus enormously facilitating the solution effort. In this work, the NEQ/OCFE model developed by Dalaouti and Seferlis14 for the optimal design of two-phase columns is expanded to efficiently model three-phase reactive distillation units. The model combines the accuracy and validity of a ratebased model with the approximating properties of the OCFE formulation. Phase boundaries are predicted and traced by adaptively placing a finite element breakpoint at the phase transition points as proposed by Swartz and Stewart.15 Phase splitting is monitored through liquid-liquid (L-L) flash calculations for the liquid phases. The model is then employed for the optimal design and process optimization of a tray reactive distillation column for the production of butyl acetate. 2. Model Formulation The NEQ model for a distillation column as presented by Higler et al.12 involves the description of simultaneous mass and heat transfer among the existing phases, phase equilibrium at the interface, and chemical reactions. For the description of mass transfer near the phase interfaces, the thin film model is used. According to the latter, the resistance to the mass transfer is solely located in two thin films adjacent to each interface. Thermodynamic equilibrium occurs only at the interface, while the bulk of the phases are in contact only with the corresponding films. Mass transfer through the thin films is modeled with the use of the Maxwell-Stefan equations for multicomponent mixtures. Hydrodynamic conditions determine the interfacial areas, film thickness, and mass transfer coefficients. In its most general case, three-phase distillation involves one gaseous phase and two liquid phases in contact with each other, leading to three interfaces and six thin films (two for each interface). According to Higler et al.,12 the second liquid phase can be considered either as immiscible with or dispersed in the other liquid phase. In the present work it is assumed that the second liquid phase is formed as a dispersed phase inside the bulk of the continuous liquid phase. Therefore, the second liquid phase is not in direct contact with the gaseous phase, and thus, only two interfaces are distinct; they are gas-liquid I and liquid I-liquid II. Possible contact of the gaseous and dispersed liquid phases would be incidental, and the influence it may have on the overall column behavior can be considered negligible due to the short contact time.13 The NEQ model has been proved to be more rigorous and accurate in the prediction of the column behavior than the conventional equilibrium model.12,16 The coupling of the NEQ model with the OCFE formulation creates an NEQ/OCFE model which combines the accuracy of the rate-based NEQ model with the model reduction properties of the OCFE formulation. In this work the generic OCFE formulation developed by Seferlis and Hrymak17 and Seferlis and Grievink,18 extended for NEQ modeling relations by Dalaouti and Seferlis,14 and finally adapted to accommodate the effect of a third phase (dispersed liquid phase) and equipped with a robust algorithm for the tracking of the phase transition boundary has been employed. Following the work of Dalaouti and Seferlis,14 the distillation column is divided into sections defined as the regions between two streams entering or leaving the column. Each section is divided into finite elements (FEs) of variable size, with each element including a specified number of collocation points (CPs) determined by the degree of polynomial approximation. The mass and energy balances are satisfied exactly only at the position of the designated CP. The collocation points are selected as the roots of discrete Hahn family orthogonal polynomials.19

Ind. Eng. Chem. Res., Vol. 49, No. 7, 2010

3277

the following form for the liquid and gaseous phases: s - sk s - sk k)0,k*j j n

WLj (s) )



s - sk s - sk k)1,k*j j

j ) 0, ..., n

(1)

j ) 1, ..., n + 1

(2)

n+1

WGj (s) )



Considering the thin film model which assumes that mass and heat transfer is concentrated in thin films adjacent to the phase interface, the model representation of each collocation point or discrete stage for a three-phase system is shown in Figure 2. Consequently, the mass balances for each phase are as follows: liquid I phase Figure 1. Schematic of the NEQ/OCFE formulation in a reactive distillation column.

Hahn orthogonal polynomials would place a collocation point at the exact location of a column tray in the limiting case that the number of collocation points equals the number of actual trays in the column, so that the full-order tray-by-tray model is recovered. Streams entering or leaving the column are modeled as discrete stages to isolate discontinuities in the concentration and temperature profiles. The column condenser, the reboiler, and the decanter are also modeled as separate equilibrium stages. In the general case, where the number of collocation points is much smaller than the number of column trays, the latter can be extracted by the means of calculating the lengths of the finite elements in the column. The element length is a design variable that represents the number of stages that the element accounts for. In this way, the total number of stages in the column equals the sum of the employed element lengths including the discrete stages such as feed stages, the reboiler, and the condenser. A schematic representation of the distillation column with one side feed stream and a decanter at the top is depicted in Figure 1. The following assumptions are made for the development of the process model: (1) one-dimensional mass and heat transport across phase interfaces, (2) thermodynamic equilibrium at phase interfaces, (3) no axial dispersion along the distillation column, (4) no entrainment of the liquid phases in the vapor phase, (5) complete mixing of bulk phases, (6) no contact between the gas and dispersed liquid phase, (7) negligible heat transfer through the films, (8) negligible pressure drop inside the column. Component molar flow rates and stream molar enthalpies are approximated by Lagrange interpolating polynomials that take

dmLI i (sj) ) L˜Ii (sj - 1) - L˜Ii (sj) + (φLI(sj) RLIb i (sj) + dt NGLb (sj)aGL - NLLIb (sj)aLL)Acol∆h (3) i i i ) 1, ..., NC

j ) 1, ..., n

liquid II (dispersed liquid) phase dmLII i (sj) ) L˜IIi (sj - 1) - L˜IIi (sj) + (φLII(sj) RLIIb i (sj) + dt LL col NLIIb i (sj)a )A ∆h (4) i ) 1, ..., NC

j ) 1, ..., n

gas phase dmGi (sj) ˜ i(sj - 1) - G ˜ i(sj) + (φG(sj) RGb )G i (sj) dt Gb Ni (sj)aGL)Acol∆h (5) i ) 1, ..., NC j ) 1, ..., n In the dynamic mass balances above, mi is the component molar accumulation which is computed from the relations

Figure 2. Schematic stage/collocation point representation of the NEQ/OCFE model.

LI LI mLI i (sj) ) φ (sj) d (sj)

i ) 1, ..., NC

LIi (sj) LIt (sj)

Acol∆h

j ) 1, ..., n

(6)

3278

Ind. Eng. Chem. Res., Vol. 49, No. 7, 2010 LII LII mLII i (sj) ) φ (sj) d (sj)

i ) 1, ..., NC

LIIi (sj) LIIt (sj)

Acol∆h

(7)

j ) 1, ..., n

mGi (sj) ) φG(sj) dG(sj) i ) 1, ..., NC

∂cGLf ∂NGLf i (sj) i (sj) + - RGLf i (sj) ) 0 GLf ∂t ∂η i ) 1, ..., NC

Gi(sj) col A ∆h Gt(sj)

(8)

0 < ηGLf e δGLf(sj)

boundary conditions GLf NLIb i (sj) ) Ni (sj)

j ) 1, ..., n

j ) 1, ..., n

(16)

where φ stands for the phase volumetric holdup fraction and d stands for the phase molar density at each collocation point. The total molar flow rates and molar fractions for each phase are computed from the component flow rates as follows:

|

GLf xLIb i (sj) ) xi (sj)

ηGLf)δGLf(sj)

|

ηGLf)δGLf(sj)

(17) i ) 1, ..., NC

j ) 1, ..., n

L-L interface, liquid I side film NC

∑ L (s ) ) L (s ) I i

I t

j

j

xI,Lb i (sj) )

i)1

i ) 1, ..., NC

LIi (sj) LIIt (sj)

j ) 1, ..., n

NC

∑ L (s ) ) L (s ) II i

II t

j

j

xII,Lb (sj) ) i

i)1

i ) 1, ..., NC

Gi(sj) yGb i (sj) ) Gt(sj)

∑ G (s ) ) G (s ) j

t

LIIi (sj) LIIt (sj)

j

i)1

i ) 1, ..., NC

(10)

(18)

0 < ηLLIf e δLLIf(sj)

boundary conditions

|

LLIf xLIb (sj) i (sj) ) xi

ηLLIIf)0

|

i ) 1, ..., NC

j ) 1, ..., n

ηLLIf)0

(19)

(11) j ) 1, ..., n

L-L interface, liquid II side film

dU(sj) ˜ LI(sj - 1) + L˜IIt (sj - 1) H ˜ LII(sj - 1) + ) L˜It (sj - 1) H dt ˜ t(sj + 1) H ˜ G(sj + 1) - L˜It (sj) H ˜ LI(sj) - L˜IIt (sj) H ˜ LII(sj) G G ˜ t(sj) H ˜ (sj) + Q(sj) (12) G j ) 1, ..., n

(sj) ∂cLLIIf ∂NLLIIf (sj) i i + - RLLIIf (sj) ) 0 i LLIIf ∂t ∂η i ) 1, ..., NC

j ) 1, ..., n

NC

LI LII LII G G {mLI i (sj) ui (sj) + mi (sj) ui (sj) + mi (sj) ui (sj)}

(20)

0 < ηLLIIf e δLLIIf(sj)

boundary conditions LLIIf NLIIb (sj) i (sj) ) Ni

In eq 12 Q(sj) denotes the net heat rate exchanged between the ˜ the stream molar enthalpy, and column and its surroundings, H U(sj) the total molar energy accumulation, which is calculated from the relation



j ) 1, ..., n

LLIIf (sj) NLIb i (sj) ) Ni

Under the assumption of negligible heat transfer through the films, the dynamic energy balance attains the form

U(sj) )

(sj) ∂cLLIf ∂NLLIf (sj) i i + - RLLIf (sj) ) 0 i LLIf ∂t ∂η i ) 1, ..., NC

j ) 1, ..., n

NC

i

(9)

|

ηLLIIf)δLLII(sj)

LLIIf xLIIb (sj) i (sj) ) xi

|

ηLLIIf)δLLIIf(sj)

(21) i ) 1, ..., NC

j ) 1, ..., n

Thermodynamic equilibrium is valid at the phase interfaces (G-LI and LI-LII) and is described by the following equations:

i)1

(13) j ) 1, ..., n The dynamic mass balances in the gas and liquid films at the two possible interfaces (G-LI and LI-LII) along with the boundary conditions are as follows:

GL GL yGL i (sj) ) Ki (sj) xi (sj)

II,LL xI,LL (sj) ) KLL (sj) i i (sj) xi

i ) 1, ..., NC

i ) 1, ..., NC

G-L interface, gas side film ∂cGf ∂NGf i (sj) i (sj) + - RGf(sj) ) 0 ∂t ∂ηGf i ) 1, ..., NC

j ) 1, ..., n

(14)

0 < ηGf e δGf(sj)

boundary conditions Gf NGb i (sj) ) Ni (sj)| ηGf)0

i ) 1, ..., NC G-L interface, liquid side film

Gf yGb i (sj) ) yi (sj)| ηGf)0

j ) 1, ..., n

(15)

j ) 1, ..., n (22)

j ) 1, ..., n (23)

and xII,LL represent the molar fractions In eqs 22 and 23, xI,LL i i of the first and second liquid phases on the L-L interface, respectively. The boundary conditions that arise from the equilibrium equations at the interfaces are GL GLf NGf i (sj)| ηGf)δGf(sj) ) Ni (sj) ) Ni (sj)| ηGLf)0

i ) 1, ..., NC

j ) 1, ..., n

(24)

Ind. Eng. Chem. Res., Vol. 49, No. 7, 2010

)

yGf i (sj)| ηGf)δGf(sj)

yGL i (sj)

xGL i (sj)

)

xGLf i (sj)| ηGLf)0

(25) i ) 1, ..., NC

i ) 1, ..., NC

j ) 1, ..., n

The reaction terms Ri used in the above equations represent the total rate of reaction per component i and are computed from the following equation: NR

i ) 1, ..., NC

j ) 1, ..., n

r

(28) rr are the rates of the reactions taking place, and νi,r are the corresponding stoichiometric coefficients. The component molar flux terms Ni at the collocation points are calculated through the multicomponent Maxwell-Stefan equations. Considering one-dimensional transport normal to each interface, the Maxwell-Stefan equations are applied at each collocation point as follows: gas phase ∂yGf i (sj) ∂ηGf

NC

Gf Gf Gf yGf k (sj) Ni (sj) - yi (sj) Nk (sj)

k)1,k*i

(P(sj)/RTGf(sj))ÐGik



)-

(29) i ) 1, ..., NC

j ) 1, ..., n

0 < ηGf e δGf(sj)

liquid phases NC-1

∑Γ

i,k(sj)

k)1

∂xGLf i (sj)

)

∂ηGLf NC

GLf GLf GLf xGLf k (sj) Ni (sj) - xi (sj) Nk (sj)

k)1,k*i

ct(sj)ÐLik



-

i ) 1, ..., NC NC-1



j ) 1, ..., n

∂xLLIf (sj) i Γi,k(sj) LLIf

(30)

0 < ηGLf e δGLf(sj)

Equations 1-33 form the NEQ/OCFE model for an element where three phases are present. An element with two phases would involve only balances without the terms associated with the third phase. Connecting equations between successive elements or between an element and a discrete stage (e.g., feed stage) are reported elsewhere.14,17,18 One major feature of three-phase distillation units is that the liquid-phase split regions (i.e., regions in which two liquid phases are present) are not known before a solution of the modeling equations is obtained. Unless the correct set of modeling equations corresponding to two- or three-phase regimes is used, the simulated results will not represent accurately the column behavior and occasionally fail to reach a feasible solution. The NEQ/OCFE model formulation is therefore properly modified to automatically detect the presence of a second liquid phase and use a suitable set of modeling equations for the column sections. The point of appearance or disappearance of the second liquid phase causes a discontinuity in the concentration profiles that would otherwise jeopardize the accuracy of the approximating interpolation polynomials. Such a situation requires the proper adaptation of the set of the modeling equations to account for changes in the phase distribution. According to the proposed methodology,15 an element breakpoint is adaptively placed at the phase transition boundary. The identification of the discontinuity point (liquidphase split or complete mixing of two liquid phases) is performed with the aid of L-L flash calculations applied at the prevailing conditions of the liquid phase at the element breakpoint. The L-L flash calculations determine whether the formation of the second liquid phase is feasible and subsequently specify whether the element breakpoint is a phase transition boundary. As the element breakpoint may be set free to move in the OCFE model formulation, the determination of the exact location of the breakpoint inside the column is achieved through the satisfaction of the bubble point calculation for the appearing second liquid phase. Similarly, in the case of disappearance of a liquid phase the element breakpoint is placed at the bubble point of the liquid-liquid equilibrium. The L-L flash calculations involve the following equations: LIIb Φ(s0) xLIb i (s0) + (1 - Φ(s0))xi (s0) ) xi(s0) i ) 1, ..., NC

)

(34)

∂η

k)1

-

NC

LLIf xLLIf (sj) - xLLIf (sj) NLLIf k (sj) Ni i k (sj)

k)1,k*i

ct(sj)ÐLik



i ) 1, ..., NC NC-1



T˜(sj),P(sj),xk(sj),k*i)1,...,NC-1

(26)

xII,LL (sj) ) xLLIIf (sj)| ηLLIIf)0 i i (27)

i ) 1, ..., NC

i,rrr(sj)

|

(33)

j ) 1, ..., n

xLLIf (sj)| ηLLIf)δLLIF(sj) ) xI,LL (sj) i i

∑ν

∂ ln γi(sj) ∂xk(sj)

j ) 1, ..., n

LLIIf NLLIf (sj)| ηLLIf)δLLIf(sj) ) NLL (sj)| ηLLIIf)0 i i (sj) ) Ni

Ri(sj) )

Γi,k(sj) ) δi,k + xi(sj)

3279

j ) 1, ..., n

∂xLLIIf (sj) i Γi,k(sj) LLIIf

NC

(31)

0 < ηLLIf e δLLIf(sj)

)

∂η

k)1

-

NC

xLLIIf (sj) NLLIIf (sj) - xLLIIf (sj) NLLIIf (sj) k i i k

k)1,k*i

ct(sj)ÐLik



i ) 1, ..., NC where

j ) 1, ..., n

(32)

0 < ηLLIIf e δLLIIf(sj)

∑ (x

LIb i (s0)

- xLIIb i (s0)) ) 0

i ) 1, ..., NC

(35)

i)1

LIb LIIb LIIb γLIb i xi (s0) ) γi xi (s0)

i ) 1, ..., NC

(36)

where γi is the activity coefficient of each component and Φ is the phase split fraction (Φ ) 1 implies a single liquid phase, whereas 0 < Φ < 1 implies two liquid phases). The L-L flash calculations are performed at the top of the finite element (s0 corresponds to the interpolation point for the liquid phase at the top of each element). Element breakpoints are therefore placed at suitable positions so that the condition Φ(s0) ) 1.0 holds.

3280

Ind. Eng. Chem. Res., Vol. 49, No. 7, 2010

3. Optimal Design Problem: Solution Algorithm 3.1. Optimal Design Problem Definition. For staged separation processes such as reactive distillation, the problem of optimal design involves the determination of the column configuration which can be expressed through the location of the feed stages, the number of stages in each column section, and the stage holdups or catalyst loads and the operating conditions such as condenser and reboiler heat duties and feed stream policy (ratio of reactants) that optimize a criterion relative to the economic feasibility of the process while satisfying the process operational and safety constraints and product specifications. The most common economic criterion is the minimization of the total annualized cost of the column that includes the process investment and operating costs. The problem statement for the optimal design of a staged reactive distillation column can therefore be formulated as follows: Given (1) a set of components and chemical reactions, (2) a flowsheet configuration, (3) the condition of the fresh feed streams, (4) a set of product specifications as well as safety and operating constraints, and (5) economic data (prices for products and reactants, investment cost data), calculate (1) the total number of stages for each column section, (2) the location of the feed stages, (3) the stage liquid-phase holdups or catalyst load for every column section, and (4) the column operating conditions that minimize the annualized process investment and operating costs subject to the steady-state NEQ/OCFE process model. The process design as derived by the above formulation ensures that all safety regulations as well as all operating, product, environmental, and economical constraints are satisfied by the optimal design under nominal model parameter values. The formulation of the design optimization problem enables the performance of extensive sensitivity analysis of the optimal solution under model or process uncertainty.18 3.2. Model Solution Algorithm. The numerical solution of the mathematical model for a three-phase distillation column is a highly demanding task, not only because of the possible phase transition along the column that results in a varying model structure but also because of the complexity caused by the large set of nonlinear differential algebraic equations. Attaining an accurate solution by solving such a large and complex set requires consistent initial values that satisfy both the differential and algebraic equations. However, as the model is formulated in a time-varying fashion the steady-state solutions can be easily obtained by integration using a consistent initial point. Moreover, as discussed in section 2, the phase splitting distribution along the column length is generally not known a priori. As a result, modeling the entire column as a three-phase column could lead to an infeasible solution as there may be sections within the column where the second liquid phase is not thermodynamically stable. To ensure that a feasible solution will be attained, the column is first modeled as a two-phase column. Liquid-liquid flash calculations are then performed in an offline fashion at the conditions obtained from the solution of the two-phase model at each collocation point. The liquid-liquid flash calculations clearly indicate whether a second liquid phase is present without interfering with the column balance equations of the two-phase model. Consequently, the sections with two or three phases are determined and finite elements with two- or three-phase modeling equations, respectively, are employed accordingly for an overall accurate representation. The predictions of the offline L-L flash calculations are utilized as starting values for the two liquid-phase concentrations in the three-phase model equations. The exact location of the phase transition inside the column is then

Table 1. Reaction Kinetic Parameters2 reaction

k0 (mol s-1)

E (kJ mol-1)

esterification hydrolysis

6.1084 × 104 9.8420 × 104

56.67 67.66

determined through the adaptive element length selection performed during the solution of the modeling equations. The main aim in the solution procedure is to provide sufficient flexibility to the overall model through the selection of an appropriate model structural combination of two- or three-phase finite elements for an accurate representation of the distillation column. In conclusion, the solution algorithm can be summarized in the following steps: Step 1. Model the column entirely as a two-phase column. Step 2. Employ the LL flash equations at the element boundaries to determine whether the formation of the third phase is possible. If the LL flash calculations are in agreement with the modeling assumption regarding the phase split, then STOP. Otherwise, proceed to Step 3. Step 3. Model those elements with potential phase split (0 < Φ(sj) < 1) as three phase elements/sections. Step 4. Determine the phase transition boundaries through adaptive placement of the element breakpoints between two successive elements with alternating two/ three phase modeling equations (condition Φ(s0) ) 1.0) and simultaneously solve the optimal design problem. Go to Step 2. The exact satisfaction of the condition Φ(s0) ) 1.0 can be relaxed to |Φ(s0)| e ε, where ε denotes a desired tolerance, to facilitate convergence. The solution of the dynamic model equations is performed using gPROMS,20 an efficient modeling environment that enables a fast and reliable convergence provided that a feasible initial point that satisfies the modeling equations is available. 4. Design Case Study: Butyl Acetate Production via Reactive Distillation 4.1. Process Description. Butyl acetate is an organic compound used mainly in the lacquer and color industries as a solvent and in the food industry as a synthetic food flavoring and is produced by the esterification reaction of butanol with acetic acid. The design of a tray reactive distillation column for the simultaneous production and separation of n-butyl acetate (BuAc) is investigated. 1-Butanol (BuOH) and acetic acid (AcOOH) are fed into the column at a ratio to be determined by the design optimization problem. The esterification-hydrolysis reactions take place in the liquid phase within the column according to the following scheme: BuOH + AcOOH T BuAc + H2O The reaction requires a strongly acidic catalyst. Homogenous acidic catalysts however have been proven to be less attractive as the corrosion problems in the column arising from their use and the formation of the poisonous side product dibutyl ether (DBE) need to be taken into consideration.9 More recent studies2,9 favor heterogeneous catalysts in the form of ion exchange resins such as Amberlyst-15. The resins are loaded in the column inside the catalytic reaction zone at the beginning of the process. The reaction constants of both reactions are considered to be Arrhenius type, following the basic form k ) k0e-E/RT

(37)

Ind. Eng. Chem. Res., Vol. 49, No. 7, 2010

Figure 3. Experimental validation of the model.

The values of k0 and E have been taken from Steinigeweg and Gmehling2 for both esterification and hydrolysis reactions and are given in Table 1. Reaction and separation proceed simultaneously inside the column, separating the products from the reactants and driving the esterification reaction to the right-hand side, thus leading to higher product yields. An equilibrium reboiler ensures the recycling of the gaseous phase back into the column, while a total condenser coupled with a decanter condenses and separates the two liquid phases of the top product stream. The overhead stream leaving the condenser is composed by water as well as butyl acetate. In the decanter this stream is separated into an aqueous phase and an organic phase due to the large miscibility gap between butyl acetate and water.8 The aqueous phase is driven away, while the organic one, composed primarily of butyl acetate, is recycled back to the top of the column. The bottom stream is almost pure butyl acetate with small traces of water, and unreacted acetic acid and butanol. The binary diffusion coefficients used in the multicomponent Maxwell-Stefan equations have been calculated using the Chapman-Enskog-Wilke-Lee model21 for the gaseous phase, while the liquid diffusion coefficients have been estimated using the method proposed by Siddiqi and Lucas.22 Thermodynamic properties, such as molar enthalpies, vapor pressures, and activity coefficients, for each component have been estimated using the UNIFAC23,24 property method, a quite reasonable choice as reported elsewhere.2 4.2. Experimental Model Validation. The proposed process model has been validated with experimental data reported in the work of Steinigeweg and Gmehling.2 In this work, the production of butyl acetate is performed in a laboratory scale packed reactive distillation column. The reactants are fed into the column at two different stages. The reactive zone utilized a strongly acidic ion exchange resin (Amberlyst-15) as the catalyst. The experimental conditions have been utilized for the NEQ/OCFE simulations. According to the experimental data and the simulated results there is no appearance of a second liquid phase within the column except the decanter at the top of the column. The column is simulated with 28 stages as indicated by Steinigeweg and Gmehling.2 Figure 3 compares the experimental data2 for the liquid composition with the simulation results obtained from the NEQ/OCFE model along the entire column. A total number of 11 collocation points distributed among 5 finite elements have been utilized in the process model. The symbols on the simulation lines correspond to the location of the collocation points and discrete stages in

3281

the NEQ/OCFE model along the column. The agreement of the experimental data and the simulated results at the two ends of the column is very good, with some discrepancies in the concentration profiles present in the reactive zone between the two feeds. However, similar deviations are also reported by Steinigeweg and Gmehling2 when comparing the experimental data to the simulation results from an equilibrium stage-bystage model, indicating that the simulated points are within the bounds of the experimental error. Therefore, a satisfactory match of the NEQ/OCFE process model has been obtained. 4.3. Design Optimization. The column is divided into two sections as shown in Figure 1, namely, the rectifying and stripping sections, which represent the section above the feed stage up to the condenser and the section below the feed stage down to the reboiler, respectively. The feed stage location is indirectly determined from the calculation of the sizes of the rectifying and stripping sections. For the particular column under study, one finite element is used in the rectifying section and four finite elements are used in the stripping section. Four collocation points are placed within each finite element. The liquid holdups along the column are assumed to be varying between sections but maintaining a constant value per section due to construction constraints. The vapor-phase holdups as well as the pressure drop inside the column and gas phase reactions are considered negligible. The NEQ/OCFE model is employed in the rectifying and stripping sections, while the feed stage is modeled as a discrete nonequilibrium stage. The condenser, the reboiler, and the decanter are modeled as single equilibrium stages. The decanter separates the liquid phase into an organicrich phase that is recycled back into the column and an aqueousrich phase that is removed as the top product. The partial differential equations referring to the mass transfer between the boundary films have been discretized using orthogonal collocation on finite elements. Specifically, an approximating scheme with four finite elements and second-order polynomials has been selected for the balance equations in the liquid and gas boundary films (eqs 14-33). Optimization (design) variables in this process are the number of stages per section, as represented by the lengths of the associated elements in each section, and the column operating conditions such as the condenser and reboiler duties, BuOH to AcOOH feed flow rate ratio, and liquid holdups per section. The problem has been formulated into the gPROMS optimizer along with given constraints such as the desirable purity of the product streams and an objective function representing the total annual cost for the process. The objective function consists of three terms representing the capital cost25 for the shell, the trays, the condenser, and the reboiler, the energy cost,25 and the material loss cost represented as the differential cost26 for useful product and reactant lost in the overhead stream as shown in top the following relation: Ltop BuAc and LBuOH represent the liquid molar

flow rates (mol/h) of BuAc and BuOH in the top water-rich product stream, respectively, while D and H are the column diameter and height. The last term in eq 39 represents the cost due to the loss of some amounts of products (BuAc) and raw materials (BuOH) in the distillate stream with the coefficients chosen in such a way to enforce a high conversion of reactants. The column diameter obeys the following flooding constraint:25

3282

Ind. Eng. Chem. Res., Vol. 49, No. 7, 2010

De

[(



)( 0.5

NC

yiMwiRgTr

i)1

Pr Grtot

]

0.5

4Grtot 3.14 × 1.5 × 3600

)

Table 2. Optimization Results

(39)

where Tr, Pr, and denote the reboiler temperature, pressure, and total gas streamflow, respectively. According to the solution algorithm presented in section 3.2, the column is initially modeled using a two-phase model. After convergence to an optimal design point, an L-L flash is performed at the condition of all collocation in the column to identify possible three-phase regions in the column. On the basis of the L-L flash results, the rich in water rectifying section exhibited liquid-phase splitting, and therefore, the three-phase NEQ/OCFE model is utilized in that section. The solution to the optimal design problem is then obtained for the hybrid three/ two-phase model. In this case, the element size in the rectifying section not only defines the number of stages in the respective section but also designates the boundary between the threephase rectifying section and the two-phase feed stage. The optimal design results for solely two-phase and hybrid two/ three-phase models are summarized in Table 2. Optimization results show that there are some significant differences in the optimal column configurations obtained by the two models. With reference to the same BuAc purity at the bottom stream, the solely two-phase model predicts a smaller column in both diameter and number of stages and a lower bottom flow rate than the hybrid two/three-phase model. The minimization of acetic acid in the top stream is directly involved in the objective function, and according to the results, the hybrid two/three-phase model predicts a 3 times lower concentration than the solely two-phase model. It is therefore evident that the inclusion of the third phase formation in the mathematical model not only alters the model structure, but also results in significant changes in the model predictions and optimal column design configurations. Such an observation implies that failure to trace and properly model the three-phase regions in the column as these spontaneously form due to the prevailing operating conditions in the column would result in erroneous column design solutions and unreliable model predictions. The remaining design parameters for the separation system such as the feed stage location, column section liquid holdups, reflux and reboil ratios, condenser and reboiler heat duties, alcohol to acid feed ratio, and total column stages as resulted from the solution of the design optimization problem are shown in Table 3. Simulation results show that at the optimum operating conditions the column has both two- and three-phase regions. In particular, the rectifying section is entirely a three-phase section, while the stripping section has only two phases. In this way, the phase transition boundary is determined by the location of the feed stage, which in place has been determined by the adaptive placement of the rectifying section breakpoint. Thermodynamic conditions indicate that the third phase appears because of the miscibility gap between butyl acetate and water. It has been noted from the simulation results that the threephase region tends to disappear near the feed stage due to the elevated concentrations of the incoming through the feed stream 1-butanol and acetic acid. The second, water-rich liquid phase is hence formed when a certain composition ratio between butyl acetate and water is reached near the top section where the reactant concentrations are low. The column composition profiles at the optimal operating point are shown in Figure 4. Stage 1 represents the column condenser, while the last stage represents the reboiler. The butyl acetate concentration maintains an increasing trend in the

model

three-phase NEQ/OCFE

two-phase NEQ/OCFE

column diameter (m) total no. of stages feed stage condenser duty (kW) reboiler duty (kW) feed ratio (BuOH/AcOOH) bottom flow (kmol/h) bottom composition (mole fraction), BuOH/AcOOH/BuAc/H2O top composition (mole fraction), BuOH/AcOOH/BuAc/H2O total annual cost ($)

1.317 60 10 -1559 1420 1.005 51.46 0.0035/0.0015/0.98/0.015

1.213 43 4 -1368 1142 1.004 42.43 0.007/0/0.98/0.013

0/0.002/0.001/0.997

0/0.006/0.001/0.993

1617.64K

1240.01K

Table 3. Design Conditions for the Butyl Acetate Reactive Distillation Column feed stream temperature (K) flow rate (kmol/h), BuOH/AcOOH length of elements in rectifying/ stripping section liquid holdups (m3) condenser rectifying section feed stage stripping section reboiler reflux ratio reboil ratio

291 34.992/34.809 (BuOH/ AcOOH ) 1.005) 8.541/12 0.676 0.040 0.250 0.197 0.593 5.52 2.52

stripping section due to the combined effect of the esterification reaction and the separation process. The concentration of BuAc drops drastically near the top of the column as the water concentration increases, but then it rises sharply shortly after as the water-rich phase is removed in the decanter and the reflux stream enters at the top of the rectifying section. The BuAc and water compositions reach the lowest level at the feed stage (stage 10). The profile for AcOOH has been omitted as it is almost identical to the BuOH profile. The rectifying section also contains the third phase (water-rich liquid phase), which is composed almost entirely of H2O. The maximum attainable product purity with this NEQ/OCFE model was 98% BuAc as the bottom product and almost pure water (99.7%) as the top product. Figures 5 and 6 depict the temperature and reaction rate (esterification) profiles throughout the column, respectively. The

Figure 4. Component composition profiles along the column.

Ind. Eng. Chem. Res., Vol. 49, No. 7, 2010

3283

Figure 5. Temperature profile along the column. Figure 7. Liquid flow rate profile. Table 4. Comparison of Design Optimization Results between Single- and Double-Feed-Stage Configurations

Figure 6. Reaction rate profile along the column.

first temperature peak around the feed stage can be attributed to the esterification reaction exothermal nature and the increased reaction rate due to the high composition of BuOH and AcOOH near the feed stage region. The temperature profile remains almost uniform in the stripping section before reaching another peak at the reboiler. The reaction rates exhibit a sharp peak at the feed stage due to the introduction of the reactants into the column followed by a large decrease and an almost steady profile for the largest part of the stripping section. The aqueousphase reaction rates on the other hand remain low because of the very low concentrations of reactants in the rectifying section (three-phase region). The liquid flow rate profile for both liquid phases inside the column is shown in Figure 7. As can be observed, the waterrich liquid phase completely disappears after the feed stage. At the left side of the graph (rectifying section) three phases exist in the column, while at the right side (stripping section) only two phases are stable. The phase boundary is decided through the optimization process by means of determining the section boundaries and therefore the sizes (lengths) of the finite elements. The addition of a second decanter aiming at removing the water-rich phase through a side stream in the rectifying section does not improve the economics of the design, and therefore, it has not been further investigated. 4.4. Double-Feed-Stage Configuration. Another column configuration involves the feeding of the reactants at different stages. Therefore, a scenario has been formulated using a column

configuration

single feed

double feed

column diameter (m) total no. of stages feed stage location finite element lengths condenser duty (kW) reboiler duty (kW) feed ratio (BuOH/AcOOH) bottom flow (kmol/h) bottom composition (mole fraction), BuOH/AcOOH/BuAc/H2O top composition (mole fraction), BuOH/AcOOH/BuAc/H2O total annual cost ($)

1.317 60 10 8.541/12/12/12/ 12 -1559 1,420 1.005 51.46 0.0035/0.0015/0.98/0.015

1.180 58 4 (AcOOH)/9 (BuOH) 2/3.177/16/16/16 -1344 931.9 1.058 52.15 0.002/0.005/0.98/0.013

0/0.002/0.001/0.997

0/0.003/0.001/0.996

1617.64K

1409.54K

configuration with two feed stages, where butanol is fed at the lower feed location and acetic acid is fed at the upper feed location. Butanol is lighter than acetic acid and prefers the vapor phase, which is a disadvantage since the reaction takes place primarily in the liquid phase. However, feeding at different stages distributes the reactants in a larger section of the column, thus enhancing the reaction rates. The two-feed-stage configuration is optimized in a fashion similar to that of the singlefeed configuration without however considering the additional investment cost of the second feed stage in the objective function value (e.g., piping and fitting costs). Five finite elements are used for the column representation, one for the section above the upper feed stage, one for the section between the two feed stages, and three for the section below the lower feed stage. The results from the design optimization are shown in Table 4 and compared to the single-feed configuration. The optimization results reveal that the rectifying section is a three-phase region whereas the sections below the upper feed stage (acetic acid feed) as well as the feed stages themselves are two-phase regions in the double-feed-stage configuration. The results show that although the total number of stages is almost the same, the column’s diameter and the necessary reboiler duty are smaller in this case due to the larger product bottom stream, which leads to a smaller total annualized cost for this column. In conclusion, the maximization of butyl acetate production is favored by a double-feed configuration and a placement of the feed stages near the top of the column. 4.5. Dynamic Simulation Studies. Modeling eqs 1-33 for the reactive three-phase column include the dynamic components, thus allowing the investigation of the dynamic behavior

3284

Ind. Eng. Chem. Res., Vol. 49, No. 7, 2010

calculations have been performed at the collocation point prevailing conditions to ensure the utilization of the correct modeling structure (three-phase vs two-phase systems) along the column. Obviously, the element breakpoints have been kept fixed during the dynamic transition. However, the adaptive placement of the element breakpoints to follow possible phase boundary drift during dynamic transition still remains a challenging research problem. 5. Conclusions

Figure 8. Dynamic response to a 6% decrease in the feed AcOOH concentration.

Figure 9. Dynamic response to a 5% increase in the feed H2O concentration.

of the column around the optimal operating point. During the dynamic simulation studies, the same assumptions regarding the negligible pressure drop and vapor holdups apply in the column. The design parameters for the staged column such as the condenser and reboiler duties, section liquid holdups, feed stage location, and number of stages per column section are taken from the design optimization problem (Tables 2 and 3). The liquid holdup per stage indicates the volume of the reacting mixture for a fixed catalyst density per volume of stage capacity that also considers the porosity of the catalyst substrate. Dynamic simulations enable the study of the influence of disturbances in the product specifications. During dynamic transition the phase boundary may migrate within the column. Figures 8 and 9 show the dynamic behavior of the process for a 6% decrease in the acetic acid flow rate in the feed stream and a 5% increase in the water content of the feed stream, respectively. In Figure 8, the butyl acetate product purity (BuAc mole fraction) gradually decreases before it reaches a final steady state of 0.964. Such behavior is attributed to the decrease in the esterification reaction rate due to the reduced amount of acetic acid in the feed. In the second case, the water entering the column acts as a hindrance for the process as the final product purity is decreased by 0.8%. An increased amount of water in the column would drive the reaction to the left (hydrolysis), reducing the effectiveness of the process for butyl acetate production. In both cases the column exhibits three phases in the section above the feed stage and two phases in the section below the feed stage. L-L flash

In this work, a model for the steady-state optimal design of a reactive distillation column with possible liquid-phase splitting is developed. The model is a combination of rate-based NEQ balance equations and the OCFE formulation that also includes a liquid-liquid flash at the element breakpoints that aims to adaptively place the element breakpoints at the exact location of the phase transition point. The model takes into consideration all the mass and energy transfer phenomena as well as thermodynamic equilibrium and chemical reactions while maintaining a compact structure that facilitates the solution of the mathematical program. The model has first been validated against experimental results from the literature and then employed for the design optimization of a tray reactive distillation column for the production of butyl acetate. At the optimal design, the column appears to have three phases in the rectifying section whereas the stripping section is an entirely two-phase region. In that case, the phase transition boundary is determined by the adaptive placement of the rectifying section element breakpoint. As design variables for the process optimization, the number of stages per column section, the feed tray location, the reactant feed ratio, the condenser and reboiler heat duties, and the liquid molar holdups per column section have been selected. Moreover, the influence of the column feed strategy (single-feed vs doublefeed configuration) on the optimal solution is investigated. A 98% butyl acetate stream has been achieved as the bottom product, while the top stream consists of 99.7% water. The column design obtained using a process model that ignores the formation of a second liquid phase fails to meet the purity constraints and product specifications of the actual threephase system. Therefore, the modeling of the second liquid phase is imperative in providing an accurate representation of the column steady-state and dynamic behavior. Furthermore, OCFE formulation successfully manages to reduce the model size, remove the need for integer-valued design variables, and provide a flexible model structure that can efficiently handle the discontinuous nature of three-phase systems while maintaining a detailed representation of the column phenomena using NEQ modeling equations. Application of the NEQ/OCFE model formulation in the study of the dynamic behavior of a threephase system with an automatic phase transition monitoring scheme remains however a challenging field of research. Nomenclature aint ) specific interfacial area (m2/m3) Acol ) column stage cross section (m2) c ) molar concentration (kmol/m3) d ) molar density (kg/m3) D ) column diameter (m) Ð ) Maxwell-Stefan diffusion coefficient (m2/s) E ) reaction activation energy (J/kmol) G ) gas molar flow rate (kmol/h) h ) column height (m)

Ind. Eng. Chem. Res., Vol. 49, No. 7, 2010 H ) stream molar enthalpy (J/kmol) k0 ) pre-exponential kinetic parameter k ) kinetic constant (kmol/s) K ) equilibrium K value L ) liquid molar flow rate (kmol/h) m ) component molar holdup (kmol) Mw ) component molecular weight n ) number of collocation points in a finite element N ) molar flux (kmol/(m2 s)) NC ) number of components NE ) number of finite elements in a column section NR ) number of chemical reactions P ) pressure (Pa) Q ) rate of heat loss/heat duty (J/s) r ) rate of reaction (kmol/(m3 s)) R ) component rate of reaction (kmol/(m3 s)) Rg ) ideal gas constant (8.314 J/(mol K)) s ) column position coordinate t ) time (h) T ) temperature (K) u ) component internal energy (J) U ) stream internal energy (J) Vb ) component molar volume (m3/kmol) W ) Langrange interpolating polynomial x ) liquid-phase mole fraction y ) gas-phase mole fraction Greek Letters γ ) activity coefficient Γ ) thermodynamic factor δ ) film thickness (m) δi,k ) Kronecker δ, 1 if i ) k, 0 if i * k ∆h ) stage height (m) η ) film coordinate µ ) component dynamic viscosity (Pa s) ν ) stoichiometric coefficient σ ) characteristic molecule diameter (Å) φ ) phase volumetric holdup fraction (m3(phase)/m3(stage)) Φ ) phase split fraction Ω ) collision integral Subscripts t ) total Superscripts G ) gas phase Gb ) gas bulk region Gf ) gas film in gas-liquid I interface GL ) gas-liquid interface GLf ) liquid film in gas-liquid I interface LI ) first liquid phase (continuous) LIb ) first liquid phase (continuous) bulk region LLIf ) liquid I film in liquid-liquid interface LII ) second liquid phase (dispersed) LIIb ) second liquid phase (dispersed) bulk region LLIIf ) liquid II film in liquid-liquid interface LL ) liquid-liquid interface

3285

Literature Cited (1) Venimadhavan, G.; Malone, M. F.; Doherty, M. F. A Novel Distillate Policy for Batch Distillation with Application to the Production of Butyl Acetate. Ind. Eng. Chem. Res. 1999, 38, 714–722. (2) Steinigeweg, S.; Gmehling, J. n-Butyl Acetate Synthesis via Reactive Distillation: Thermodynamic Aspects, Reaction Kinetics, Pilot-Plant Experiments and Simulation Studies. Ind. Eng. Chem. Res. 2002, 41, 5483–5490. (3) Steyer, F.; Qi, Z.; Sundmacher, K. Synthesis of Cyclohexanol by Three-Phase Reactive Distillation: Influence of Kinetics on Phase Equilibria. Chem. Eng. Sci. 2002, 57, 1511–1520. (4) Khaledi, R.; Bishnoi, P. R. A Method for Steady-State Simulation of Multistage Three-Phase Separation Columns. Ind. Eng. Chem. Res. 2005, 44, 6956–6855. (5) Khaledi, R.; Bishnoi, P. R. A Method for Modeling Two- and ThreePhase Distillation Columns. Ind. Eng. Chem. Res. 2006, 45, 6007–6020. (6) Shyamsundar, V.; Rangaiah, G. P. A Method for Simulation and Optimization of Multiphase Distillation. Comput. Chem. Eng. 2000, 24, 23–37. (7) Bruggemann, S.; Oldenburg, J.; Zhang, P.; Marquardt, W. Robust Dynamic Simulation of Three-Phase Reactive Batch Distillation Columns. Ind. Eng. Chem. Res. 2004, 43, 3672–3684. (8) Radulescu, G.; Gangadwala, J.; Paraschiv, N.; Kienle, A.; Sundmacher, K. Dynamics of Reactive Distillation Processes with Potential Liquid Phase Splitting Based on Equilibrium Staged Models. Comput. Chem. Eng. 2009, 33, 551–574. (9) Gangadwala, J.; Kienle, A.; Stein, E.; Mahajani, S. Production of Butyl Acetate by Catalytic Distillation: Process Design Studies. Ind. Eng. Chem. Res. 2004, 43, 136–143. (10) Gangadwala, J.; Kienle, A. MINLP Optimization of Butyl Acetate Synthesis. Chem. Eng. Process. 2007, 46, 107–118. (11) Lao, M.; Taylor, R. Modeling Mass Transfer in Three-Phase Distillation. Ind. Eng. Chem. Res. 1994, 33, 2637–2650. (12) Higler, A.; Chande, R.; Taylor, R.; Baur, R.; Krishna, R. Nonequilibrium Modelling of Three-Phase Distillation. Comput. Chem. Eng. 2004, 28, 2021–2036. (13) Eckert, E.; Vanek, T. Some Aspects of Rate-Based Modeling and Simulation of Three-Phase Distillation Columns. Comput. Chem. Eng. 2001, 25, 603–612. (14) Dalaouti, N.; Seferlis, P. A Unified Modeling Framework for the Optimal Design and Dynamic Simulation of Staged Reactive Distillation Processes. Comput. Chem. Eng. 2006, 30, 1264–1277. (15) Swartz, C. L. E.; Stewart, W. E. Finite-Element Steady State Simulation of Multiphase Distillation. AIChE J. 1987, 33 (12), 1977–1985. (16) Klo¨ker, M.; Kenig, E. Y.; Hoffmann, A.; Kreis, P.; Go´rak, A. RateBased Modeling and Simulation of Reactive Separations in Gas/VapourLiquid Systems. Chem. Eng. Process.: Process Intensif. 2005, 44 (6), 617– 629. (17) Seferlis, P.; Hrymak, A. N. Optimization of Distillation Units Using Collocation Models. AIChE J. 1994, 40 (5), 813–825. (18) Seferlis, P.; Grievink, J. Optimal Design and Sensitivity Analysis of Reactive Distillation Units Using Collocation Models. Ind. Eng. Chem. Res. 2001, 40, 1673–1685. (19) Stewart, W. E.; Levien, K. L.; Morari, M. Simulation of Fractionation by Orthogonal Collocation. Chem. Eng. Sci. 1985, 40, 409. (20) Process Systems Enterprise Ltd. gPROMS Introductory and AdVanced User’s Guide, R.2.3; London, 2003. (21) Reid, R. C.; Prausnitz, J. M.; Poling, B. E. The Properties of Gases and Liquids; McGraw-Hill: New York, 1987. (22) Siddiqi, M. A.; Lucas, K. Correlations for Prediction of Diffusion in Liquids. Can. J. Chem. Eng. 1986, 64, 839. (23) Fredenslund A.; Gmehling J.; Rasmussen P. Vapor-Liquid Equilibria Using UNIFAC; Elsevier: Amsterdam, 1977. (24) Nova´k J. P.; Matousˇ J.; Pick J. Liquid-Liquid Equilibria: Studies in Modern Thermodynamics 7; Elsevier: Amsterdam, 1987. (25) Douglas J. M. Conceptual Design of Chemical Processes; McGrawHill: Singapore, 1988. (26) Edgar T. F.; Himmelblau D. M. Optimization of Chemical Processes; McGraw-Hill: New York, 2001.

ReceiVed for reView August 10, 2009 ReVised manuscript receiVed January 24, 2010 Accepted February 4, 2010 IE901260B