3672
Ind. Eng. Chem. Res. 2004, 43, 3672-3684
Robust Dynamic Simulation of Three-Phase Reactive Batch Distillation Columns Stefan Bru 1 ggemann, Jan Oldenburg, Ping Zhang,† and Wolfgang Marquardt* Lehrstuhl fu¨ r Prozesstechnik, RWTH Aachen,Templergraben 55, D-52056 Aachen, Germany
The dynamic simulation of heterogeneous distillation columns is a demanding task because the number of liquid phases can vary over time on each tray. The success of the simulation depends on the rapid, reliable, and robust assessment of the number of phases and the correct initialization of the compositions of the liquid phases. In this work, a flexible simulation strategy is presented that decouples the phase split and flash algorithms from the process model. Therefore, no discontinuous model switching occurs in the process model, and the column can easily be simulated using standard tools employing a simultaneous equation-oriented simulation strategy. The phase state is assessed with a rapid and reliable homotopy-continuation algorithm that is well-suited for dynamic simulation because it decomposes the phase stability problem into two steps to minimize the necessary computational effort. The simulation strategy is applied to the production of butyl acetate from butanol and acetic acid in a three-phase reactive batch distillation column as an illustrative example. Introduction Batch processes have been used in the chemical industry since the beginning of chemical processing. They have been the focus of renewed interest in recent years, especially for the production of low-volume, highvalue specialty chemicals. The high flexibility of batch units also makes them attractive for seasonal production campaigns. Processes frequently consist of at least one batch reactor and at least one batch distillation column. For equilibrium-limited reactions such as the esterification of methanol to methyl acetate, it has been demonstrated in continuous processing that the application of process intensification results in a significant improvement and greatly reduces process cost.1 The combination of reaction and separation in reactive distillation2 has the potential of achieving essentially complete conversion of the reactants because the products are continuously removed from the reaction zone, thus shifting the reaction equilibrium toward the forward reaction. Venimadhavan et al.3 showed that reactive distillation in batch columns can also be applied successfully to the esterification of butanol to butyl acetate, a powerful solvent that is mainly used in coatings. Because of the existence of several low-boiling azeotropes, butyl acetate cannot be recovered directly as the top product of the reactive column. However, the production of high-purity butyl acetate can be facilitated by using a clever reflux policy that exploits the appearance of a miscibility gap at the condenser and in the upper part of the column. Although extremely beneficial for the process itself, the appearance of a second liquid phase makes the dynamic simulation of the column a much more difficult task. The main obstacles here lie * To whom correspondence should be addressed. Tel.: +49241-8096712. Fax: +49-241-8092326. E-mail: marquardt@ lpt.rwth-aachen.de. † Present address: Institut fu¨r Verfahrenstechnik/TVT, Martin-Luther-Universita¨t Halle-Wittenberg, FB Ingenieurwissenschaften, D-06099 Halle (Saale), Germany.
in the rapid, robust, and reliable determination of the phase state on each tray over time; the calculation of the compositions on the heterogeneous trays, and the determination of the switches in the process model associated with a change in the phase state on some tray. In this work, the dynamic simulation of heterogeneous reactive batch distillation using a standard equationoriented dynamic simulation package is addressed. First, a homotopy-continuation method for the rapid and reliable assessment of phase stability at each point in time is reviewed. As explained in the subsequent section, this augmented flash with a variable number of phases is connected to the column model as an external procedure. The computationally beneficial properties of the decoupling of the DAE column model and the phase stability test are discussed in detail. An outline of the software implementation is given next. The new simulation approach is then applied to the production of butyl acetate in a heterogeneous batch distillation column. A case study shows how the dynamic simulation can be employed to find an improved operational strategy for this process. Finally, the consequences of a failure of the phase test are discussed. Determination of Phase Stability in Simulation In the past three decades, several algorithms for the determination of phase stability have been under research. Most algorithms for steady-state or dynamic simulation of heterogeneous columns are based on local solution methods for the determination of phase stability. In 1990, Cairns and Furzer4 presented a review of algorithmic approaches in use at that time. They also suggested a new approach based on the NapthaliSandholm algorithm. In this approach, the phase state is determined by a local analysis of the tangent plane criterion of Michelsen,5 which, in addition, also provides a first initialization for the compositions of the heterogeneous liquid phases. Usually, multiple minima of the tangent plane function exist, and therefore, local minimization might lead to erroneous conclusions depending
10.1021/ie034045v CCC: $27.50 © 2004 American Chemical Society Published on Web 01/15/2004
Ind. Eng. Chem. Res., Vol. 43, No. 14, 2004 3673
on the initialization. Another approach was presented by Buzzi Ferraris and Morbidelli,6 who used a series of local flash calculations to solve the phase stability problem. A modification of this algorithm was used by Kingsley and Lucia7 for the global steady-state optimization of a heterogeneous distillation process with a tunneling algorithm. A similar approach was employed by Eckert and Kubicek8 for the dynamic simulation of a distillation column. Landwehr et al.9 introduced a phase test based on the evaluation of isoactivity constraints and additional mass balances. Finally, the τ-method defines slack variables for the closure relation that allow a temporary deviation of the simulation from physical reality in favor of more robust convergence properties.10 In so far as their modes of operation are documented, all state-of-the-art commercial steady-state or dynamic process simulators are also based on similar local phase stability tests. However, dynamic simulation studies of heterogeneous batch distillation using commercial software are scarce. A notable exception is the paper of Rodriguez-Donis et al.11 The low computational demand of the local solution methods makes them well-suited for steady-state or dynamic simulation. However, such methods are known to fail in some cases because good initial values for the heterogeneous flash calculation cannot always be obtained.12 To overcome this shortcoming, global methods based on either global optimization (McDonald and Floudas,13 Yushan and Zhizhong14) or interval methods (Tessier et al.15) exploit sufficient stability criteria, such as the tangent plane criterion of Michelsen.5 Although these global approaches guarantee the detection of any phase split, they have a high computational demand and require specialized model formulations and solvers. These deficiencies currently prohibit the use of such methods for routine simulation and optimization purposes. In the past several years, homotopy-continuation methods have been developed that promise a good compromise between the computational efficiency of the local algorithms and the reliability of the global approaches. The methods proposed by Sun and Seider16 and Wasylkiewicz et al.17 aim at finding all solutions to the tangent plane criterion. Then, the compositions of the equilibrium phases are determined by solving the heterogeneous flash problem either by performing a second continuation step or by applying a local solver that is initialized with the global solution of the tangent plane criterion. In both cases, a two-step algorithm must be solved for each test mixture. Recently, Bausa and Marquardt12 presented a novel hybrid homotopy-continuation approach that incorporates a priori knowledge about the phase diagram and decomposes the phase stability problem into two steps. In the preprocessing step, all heterogeneous regions of the system are divided into convex regions. For each region, the overall, liquid-phase, and vapor-phase compositions, as well as the phase fraction, pressure, and temperature, of one reference state, uref ) [xrefT, xIrefT, I T T T xII 0 , yref , φref, pref, Tref] , are stored, where xref, xref, II xref, and yref are the component concentration vectors of dimension C in the bulk liquid, both in the liquid phases and in the vapor phase. Here, C denotes the total number of components contained in the mixture. This step needs to be performed only once for a given pressure p and temperature T. It can be carried out using any local or global method for the detection of
phase stability. Usually, a local search in all binary subsystems, as suggested by Bausa and Marquardt,12 is sufficient. The actual vapor-liquid-liquid flash calculation at bubble-point condition on each stage is performed by solving the set of flash equations F(u) ) 0 with
Fi(u) ) (1 - φ)xIi + φxII i - xi ) 0, i ) 1, ..., C FC+i(u) ) yi - Ki(xI,y,p,T)xIi ) 0, i ) 1, ..., C II
F2C+i(u) ) yi - Ki(x
,y,p,T)xII i
) 0, i ) 1, ..., C
(1) (2) (3)
C
F3C+1(u) )
(yk - xkI) ) 0 ∑ k)1
(4)
C
F3C+2(u) )
(yk - xkII) ) 0 ∑ k)1
(5)
for u ) [xT, (xI)T, (xII)T, yT, φ, T, p]T where the overall liquid composition x and the pressure p define the test mixture on some tray of the column model. However, the solution of eqs 1-5 guarantees a reliable indication of the phase situation only if good initial values for u are always present. This is ensured by defining the homotopy
Hi(u,λ)) (1 - φ)xIi + φxII i - [λxi + (1 - λ)xref,i], i ) 1, ..., C (6) that replaces the mass balance of eq 1 of the heterogeneous flash. At λ ) 0, the results of the homotopy u ) uref are already known from the preprocessing step. A predictor/corrector continuation algorithm is used to obtain the results at λ ) 1. Note that such a homotopy run is necessary for all reference states that were found during the preprocessing step. Assuming at most two liquid phases, a nontrivial solution with positive phase fractions φ and (1 - φ) indicates a phase split of the test mixture x. Otherwise, the mixture is homogeneous, and the equilibrium vaporphase composition y and the bubble-point temperature T can easily be calculated by solving the homogeneous flash equations
0 ) yi - Ki(x,y,p,T)xi, i ) 1, ..., C
(7)
C
0)1-
∑ yi
(8)
k)1
with a suitable initialization strategy. To produce on-specification products, the operational policy for a column often requires the appearance of a phase split on some prespecified tray. For example, some processes require a heterogeneous mixture at the top of the column because only an aqueous phase or an organic phase is to be withdrawn as the distillate product. This phase split can appear or disappear depending on the choice of the reflux ratio. Correspondingly, the set of equations describing this tray can change. To safely detect and handle such explicit model switching, therefore, it is not only required to know whether a phase split occurs or not, but also helpful to know the distance of the current operating point to the phase boundary. An approximate measure Ψ for the distance to the phase boundary can be obtained directly
3674 Ind. Eng. Chem. Res., Vol. 43, No. 14, 2004
Figure 1. Calculation of the phase state and distance to the phase boundary Ψ.
from the result of the phase stability algorithm. The use of Ψ as a switching criterion for the reflux policy is discussed in one of the following sections. The determination of Ψ is illustrated in Figure 1 for three representative test mixtures. Test mixture xA is heterogeneous (vapor-liquid-liquid equilibrium; VLLE) and splits into the liquid phases xIA and xII A . These liquid-phase compositions are determined by continuation in the heterogeneous domain starting from the reference solution xref. The distance ΨA of xA to the phase boundary can then be calculated from the tie line as the distance of the closest equilibrium composition on the binodal curve. Test mixture xB is homogeneous but still located below the tangent to the critical point of the miscibility gap. The continuation is first performed in the heterogeneous region until the binodal curve is reached. Then, the continuation is continued in the homogeneous region. However, according to the concept of the negative flash, it is still possible to find liquid-phase compositions xI and xII that fulfill eqs 1-5. We therefore speak of negative vapor-liquid-liquid equilibrium (NVLLE). Two equilibrium phase compositions xIB and xII B are found for bulk composition xB with a phase ratio φB > 1. The distance ΨB to the boundary is calculated analogously to the heterogeneous case. The last case is shown by test mixture xC, which is homogeneous and above the tangent to the critical point. A negative vapor-liquid-liquid equilibrium solution can no longer be found after the homogeneous continuation path crosses the tangent to the critical point. In this case, the distance to the miscibility gap is approximated by the distance ΨC from the test mixture composition xC to the critical point xcrit. Table 1 summarizes the equations for the calculation of Ψ for all three phase state situations. The decomposition of the homotopy phase stability test into preprocessing and flash steps makes this algorithm well-suited for dynamic simulation. This is especially true for the simulation of distillation columns, given that the temperature is at the boiling point on every stage and the pressure range along the column is small or can even be neglected completely for simulation. Therefore, the heterogeneous regions do not change qualitatively, and one application of the preprocessing step at boiling conditions is sufficient for the whole simulation, as the dynamic model contains only the numerically less intensive flash calculation. It should
Table 1. Calculation of the Distance to the Phase Boundary Ψ type of equilibrium VLLE
sign of Ψ
definition of Ψ
Ψg0
If φ < 0.5, then φmin ) φ; else φmin ) 1 - φ
x
NVLLE
Ψ
