Modeling and Simulation of Reactive Distillation Operations

---

3188

Ind. Eng. Chem. Res. 1997, 36, 3188-3197

Modeling and Simulation of Reactive Distillation Operations P. A. Pilavachi Faculty of Applied Sciences, Universite´ Libre de Bruxelles, 50, Av. Fr. Roosevelt, 1050 Brussels, Belgium

M. Schenk, E. Perez-Cisneros, and R. Gani* Department of Chemical Engineering, Technical University of Denmark, 2800 Lyngby, Denmark

Important aspects related to modeling and simulation of reactive distillation processes are presented. Reactive distillation processes are system specific and are subject to the sensitivity of the model parameters. The sensitive model parameters have been identified as those belonging to the models describing the physical and/or chemical equilibrium of the reactive system. The influence of the sensitive model parameters on simulation/design is highlighted through a systematic analysis of the models typically employed for steady-state and dynamic simulations of reactive distillation operations. For reliable and consistent simulation and design of reactive distillation operations, a necessary first step is a systematic analysis of the model parameters and the design/operational variables. Validated numerical results from test problems involving two reactive systems are presented. Introduction Separation by reactive distillation has become a feasible alternative when separating liquid mixtures involving azeotropes or isomers. Separation of these mixtures by conventional distillation requires highenergy costs, while separation with solvent-based distillation requires a suitable solvent, which in some cases may not be available. This solvent may also contaminate the product. Another important advantage of reactive distillation is that it is able to increase the product yield, by allowing separation of the reaction products from the reactants in the same unit operation. Examples of application of reactive distillation can be found in esterification processes (Agreda et al., 1990; Sawistowski and Pilavachi, 1988), separation of isomers (Saito et al., 1971), and in the manufacture of “antiknock” enhancers such as methyl tert-butyl ether (MTBE) and tert-amyl methyl ether (TAME) (e.g., De Garmo et al., 1992; Ignatius et al., 1995). Although the use of reactive distillation as a separation technique has been known for a long time, it has only been shown recently that we can solve separation problems which were previously considered difficult or unsolvable. Despite the recent advances by Doherty and co-workers (Barbosa and Doherty, 1988; Doherty and Buzad, 1992; Ung and Doherty, 1995), it can be claimed that there are still no generally accepted design methods for reactive distillation. One important reason is the difficulty associated with the modeling of these processes and the parametric uncertainty of the model parameters (Gani and Constantinou, 1996). As a result, reliable and consistent models (for property prediction, reaction kinetics, and process simulation) are not available, and the dependence of the model parameters is unknown. For nonreactive systems, the question of model-parameter uncertainty has been studied in detail by Laroche et al. (1992) and Whiting et al. (1993). It is obvious that without appropriate “tools”, the design and operation of any new process would be subject to unacceptable uncertainties. Modeling of reactive separation processes requires simultaneous prediction of reaction and separation. Several alternatives are possible. For example, treat * Author to whom all correspondence should be addressed. S0888-5885(96)00640-9 CCC: $14.00

the process as a nonequilibrium system (kinetically controlled reaction plus separation through mass/heat transfer), or treat the process as only in physical equilibrium with the kinetically controlled reaction, or treat the process as an equilibrium process (chemical and physical equilibrium). Even though no reactive distillation process will ever achieve total equilibrium, nevertheless, the equilibrium-based model provides the theoretical limits of the achievable separation. Models based only on physical equilibrium and not on chemical equilibrium are a more practical approach to the simulation of industrial reactive distillation operations. The nonequilibrium-based (or rate-based) models appear to be the most correct and appropriate models for simulation of any reactive distillation operation. The key advantage of the rate-based models is that the use of empirical factors such as efficiencies (for tray columns) and HEPTs (for packed columns) is completely avoided. The values of the efficiencies are uncertain due to the effect of the separation coupled with the chemical reaction. In rate-based models, this effect is considered explicitly in the mass and heat balance equations for the liquid and vapor phases. Although rate-based models for distillation are now quite common (see, for example, Seader, 1989; Taylor and Lucia, 1995), such models for dynamic simulation of reactive distillation operations are still in the development stage. A commercial version of a generalized steady-state rate-based model, RATEFRAC, capable of handling reactive systems, has been described by Sivasubramanian and Boston (1990). Presently, the main drawback of these models is that they need a large number of parameters (some of which may be difficult to obtain) for modeling the mass/heat-transfer rates, for example, model parameters for many of the needed transport properties. In addition, the rate-based models require the specification of design parameters such as packing type and shape, column diameter, etc.; design of reactive distillation columns through rate-based models is not simple. Therefore, in this paper, the models based on chemical and/or physical equilibrium are only considered. The objective of this paper is to investigate, in a systematic manner, the important modeling and simulation aspects with respect to correct and consistent simulation of reactive distillation operations. First, an © 1997 American Chemical Society

Ind. Eng. Chem. Res., Vol. 36, No. 8, 1997 3189

analysis of the modeling options for steady-state and dynamic simulations of chemical and/or physical equilibrium is presented. Next, a procedure for the analysis of the model parameters is provided. This analysis includes model/parameter validation, sensitivity analysis of the model parameters, and their effects on the design of reactive distillation operations. Finally, through numerical examples, the model-parameter-simulation/ design relationships are highlighted through two wellknown reactive systems (the MTBE process and an esterification process).

γp,i ) f(xp,i,T,πactivity)

(5)

where eqs 3 and 4 are considered when the φ-φ phase equilibrium approach is applied, while eqs 3 and 5 are for the γ-φ approach. The functionality of the equilibrium constant and the extent of the reaction and reaction rate are

Ka ) f(T,πpure)

(6)

ξ ) f(φ,γ,xp,i,yp,i,T,P)

(7)

ra ) f(φ,γ,Ka,xp,i,yp,i)

(8)

Modeling of Reactive Distillation In this section, only the chemical and/or physical equilibrium-based models will be described. Typically, the mathematical models are represented by a set of differential and algebraic equations for the dynamic model and a set of nonlinear algebraic equations for the steady-state model. The model equations can be decomposed into subsets of balance (mass and energy) equations, equilibrium equations (chemical and/or physical equilibrium computations), hydraulic equations (relationships between column geometry and liquid/ vapor flows), and defined equations (reaction kinetics, controller equations, constraint equations, etc.). In this work, steady-state and dynamic models for distillation columns (Gani et al., 1986) have been updated to also handle reactive distillation. Thus, most of the model equations and the solution procedure have been described elsewhere and need not be given here. In this section, the equilibrium equations and the defined equations are briefly described together with the calculation procedure. Note that each term of the balance equations is computed through the other subset of equations. The complete set of equations is given in the Appendix. Equilibrium Equations. Physical Equilibrium (PE). In the case of only physical equilibrium, as proposed by Gani et al. (1986), the phase equilibrium calculations involve the prediction of the bubble-point temperature and the corresponding vapor phase in equilibrium, given the pressure and the liquid composition (for each plate). The reaction term in the balance equation is determined, in this case, through the kinetic equations (see “defined equations” section in the Appendix). In this case, the parameters (π) that may affect the simulation results are the phase equilibrium model parameters and the coefficients for the kinetic expression. It should be noted that in this case, the phase equilibrium model parameters may be heavily correlated to the reaction kinetic parameters. That is, if the phase equilibrium model parameters were estimated for a given reaction kinetics expression, the simulation results may not match the experimental data unless the corresponding sets of model parameters are used. The equilibrium relations are

yp,i ) Kp,ixp,i Kp,i )

L φp,i V φp,i

psi γp,i )

V Pφp,i

(1) (2)

The nonideality of the phases is represented as V φp,i ) f(yp,i,T,P,πEOS) L φp,i

) f(xp,i,T,P,πEOS)

Chemical and Physical Equilibrium (CPE). In the case of simultaneous chemical and physical equilibrium, in order to compute the reactive bubble-point temperatures, simultaneous chemical and physical equilibrium must be considered. It can be noted in the above equations (eqs 3-5) that computation of chemical and/ or physical equilibrium requires the use of a thermodynamic model to supply the fugacity coefficients (φ) or activity coefficients (γ) for the mixture components in their respective phases. From the same equations, it can also be seen that chemical equilibrium is a function of the model parameters (π). In addition, eq 6 shows that the equilibrium constant may also be affected by another set of model parameters (πpure) which are related to the estimation of pure-component properties. One way of combining CPE in reactive distillation modeling is to employ the chemical model approach recently proposed by Pe´rez-Cisneros et al. (1997). This approach has been derived from the chemical model theory, where the equations of chemical equilibrium with any appropriate physical model yielding the chemical potentials are incorporated into an element-based model (called the chemical model). The solution of the chemical model equations together with the conditions of equilibrium provides the element composition for the reactive system. Indeed, the use of the chemical model permits the statement of the combined chemical and physical equilibrium problem to be formally identical to the physical equilibrium problem for a mixture of elements. That is, the relation between the elemental composition vector b and the elemental chemical potential vector λ (which are the Lagrange multipliers used in the minimization of the Gibbs free energy) is completely equivalent to the relation between the component molar composition vector n and the component chemical potential vector µ. Also, the transformation of the chemical physical equilibrium problem for a multicomponent system into a physical equilibrium problem for an element system permits the definition of element azeotrope, which is identical (in principle) to the azeotropes found in nonreactive systems. Thus, with the element balance approach, the conditions for the existence of element azeotropes are the point where the compositions of element j (j ) 1, 2, ..., M) in all coexisting phases are identical, or equivalently, the element relative volatility of element j (j ) 1, 2, ..., M) with respect to element R is unity. The element relative volatility of element j with respect to reference element R is defined as (for a two-phase system)

(3) (4)

Rj,R )

Kj KR

(9)

3190 Ind. Eng. Chem. Res., Vol. 36, No. 8, 1997

where the element separation factor, Kj, is

Kj )

WVj WLj

(10)

It can be noted that the above conditions are identical to the set of conditions that define the existence of azeotropes in nonreactive systems. A more detailed description of the element base chemical model for CPE can be found in Pe´rez-Cisneros et al. (1997). By using the chemical model approach, eqs 6 and 7 are included in the solution of the physical equilibrium problem (eqs 1-5) for a chemical system defined in terms of elements. At the physical and chemical equilibrium condition, the extent of reaction is calculated (eq 7). The reaction term in the mass and energy balance equations, in this case, is represented as a function of the extent of reaction. Thus, by simply adding the reactive bubble-point temperature algorithm of Pe´rez-Cisneros et al. (1997) and eq 7 into the set of model equations for a nonreactive distillation model, a CPE-based reactive distillation model is generated. Solution of the Model Equations. For steady-state simulation, the model equations represent a set of nonlinear equations consisting of the mass and energy balance equations, the equilibrium equations, and the defined equations. This set of equations is solved typically with the Naphatali-Sandholm (1971) method. For dynamic simulation, the model equations represent a set of differential and algebraic equations (DAEs) consisting of the balance equations, equilibrium equations, hydraulic equations, and defined equations. Perregaard et al. (1992) have discussed various simulation strategies for solving sets of DAEs. Note that, in principle, the dynamic model can also be used for steady-state simulation (Gani and Cameron, 1989). Model/Parameter Analysis In this paper, the effect of the model parameters (π) in the γ-φ expression for the equilibrium constant is investigated within the context of simulation (steadystate and dynamic) and design of reactive distillation operations. Through this analysis, a set of parameters which are sensitive for a given reactive system is identified. At the same time, the choice of the model parameters with respect to validation of the model is investigated. The analysis is performed in three levels. In the first (model validation) level, the reliability and consistency of the “tools” to be employed are established through a systematic analysis of the suitability of the models (and their parameters) for the estimation of the relevant pure component and mixture properties. The prediction of CPE as well as steady-state and dynamic simulations needs to be validated before a design/analysis study can be initiated. The objective of level 1 is to identify the set of sensitive properties for any reactive distillation problem, to select the appropriate models, and, finally, to validate the models and their parameters. In the second (relationships) level, the relationship between model parameters and design objectives and the relationship between the reaction kinetics, physical properties, steady-state operation, and transient behavior are established. For example, how are the chemical equilibrium constants and the activity coefficient model parameters related to the conditions of operation of reactive distillation columns? The objective of level 2

is to determine the sensitivities of the selected model parameters to the prediction of CPE and the variables related to the condition of operation of the process. Also, the influence of the phase diagrams and the location of the element azeotropes on the design and operation of reactive distillation columns are established. In the third and final level, the determined relationships from the lower levels are employed to analyze their effect in the design and operation of reactive distillation processes. Computational Tools Required. The proposed analysis approach is largely dependent on computational tools such as estimation methods for purecomponent properties (heats of formation, vapor pressures, enthalpies, etc.), mixture properties (fugacity coefficients, activity coefficients, densities, etc.), and physical and/or chemical equilibria. These estimation methods need to be employed in the calculation procedures that generate phase diagrams and locate the existence of reactive/nonreactive azeotropes and in the process simulation models (steady-state and dynamic). In this work, for the analysis of CPE and its relation to design/analysis, a new computational procedure employing the element balance approach (Pe´rez-Cisneros et al., 1997) has been applied. The advantage of this approach over other approaches (e.g., Ung and Doherty, 1995) is that the method of solution is similar to the well-known nonreactive physical equilibrium problems. Also, the model variables (element fractions) can be directly employed in phase diagrams and for design/ analysis of reactive distillation. Finally, this approach reduces the dimension of the problem, thereby allowing the representation of many multicomponent reactive systems in two- or three-dimensional phase diagrams. Results and Discussion The proposed analysis approach has been tested through a number of reactive systems (for example, MTBE process, esterification of alcohols, separation of xylenes, production of ethanol from ethylene and water, production of ethylene glycol, and more). In this paper, the application results are highlighted through two of the reactive systems studied. They are the well-known MTBE process (the production of methyl tert-butyl ether from the reaction of isobutene and methanol with the presence of one inert component) and the esterification process (the production of ethyl acetate and water from ethanol and acetic acid). Application results for the other reactive systems can be obtained from us. Level 1: Validation of the Computational Tools. Steady-state simulation results for the two selected reactive systems have been reported in the literature. For the MTBE process, it was tried to reproduce the results of Jacobs and Krishna (1993) and Hauan et al. (1995). Both sets of authors have reported the existence of multiple solutions. For the esterification process, experimental data reported by Suzuki et al. (1971) have been successfully reproduced by steady-state and dynamic simulations. Set of Sensitive Properties. For both reactive systems, the following set of properties have been found to be highly sensitive to the simulation results: purecomponent Gibbs free energies, temperature dependence of the chemical equilibrium constant, liquid-phase activity coefficients (because the reactions are assumed to occur only in the liquid phase which represents a nonideal mixture), and, in the case of the esterification example, the effect of dimerization in the vapor-phase

Ind. Eng. Chem. Res., Vol. 36, No. 8, 1997 3191 Table 1. CPE Modeling Alternatives for Esterification of Acetic Acid with Ethanol (γ-φ Approach) liquid-phase model

source for liquid-phase model parameters

vapor-phase model

UNIFAC UNIFAC UNIFAC UNIQUAC UNIQUAC Margules-modified expression WILSON

Fredenslund et al. (1977) Fredenslund et al. (1977) Fredenslund et al. (1977) Kang et al. (1992) Kang et al. (1992) Suzuki et al. (1970) Suzuki et al. (1970)

ideal SRK virial equation ideal SRK ideal SRK

source for vapor-phase model parameters binary interaction coefficients set to zero Hayden and O’Connell (1975) binary interaction coeffs set to zero binary interaction coeffs set to zero

Table 2. CPE Modeling Alternatives for the MTBE System (γ-φ Approach) liquid-phase model

source for liquid-phase model parameters

vapor-phase model

source for vapor-phase model parameters

UNIFAC UNIFAC UNIQUAC

Fredenslund et al. (1977) Fredenslund et al. (1977) Rehfinger and Hoffmann (1990)

ideal SRK SRK

binary interaction coeffs set to zero binary interaction coeffs set to zero

Table 3. Kinetic Expressions for the Two Case Studies case studies MTBE

kinetic expression

ra ) k

(

aIB 1 aMTBE 2 aMeOH Ka a MeOH

ra ) KDCAACE - KICEACW

)

equilib const

source for model parameters

Ka ) f(T)

Rehfinger and Hoffmann (1990) Izarraraz et al. (1980)

esterification

ra ) k

( ) xAA,p

2

Ka )

KD KI Suzuki et al. (1971)

NC

∑F

i,pxi,p

i)1

fugacity coefficients. The criterion for selection of these sensitive properties has been to note the changes in simulation (steady-state) results due to the changes in the predicted property values. Note that a study of the model parametric uncertainty has not been made here. The effect of the model parameters corresponding to the prediction of these properties needs, therefore, to be investigated. Since different choices of the model parameters are available, a careful selection needs to be made. The various chemical and/or physical equilibrium modeling alternatives considered in this paper for the two examples are given in Tables 1 and 2 in terms of the liquid-phase and vapor-phase models together with the sources for their model parameters. In the case of only physical equilibrium, the corresponding kinetic expressions and the sources for their model parameters are given in Table 3. Selection/Validation of Models. Steady-state simulations for the MTBE process have been performed with ASPEN-PLUS (1995) for the model combinations given in Table 2. With the combination of the UNIQUAC model for the liquid-phase activity coefficients, the Redlich-Kwong EOS for the vapor phase, the Rehfinger-Hoffmann (1990) expression for the temperaturedependent equilibrium constant term, and the nonideal expression for the stoichiometry-based equilibrium constant term, it was possible to match the results (including multiple solutions) reported by Jacobs and Krishna (1993) and those of Hauan et al. (1995). Figure 1a shows the variation of fractional conversion of isobutene as a function of the methanol feed location. Changing the UNIQUAC model with the UNIFAC model (for the liquid-phase activity coefficient), although the highconversion solutions were found to be approximately similar, the low conversion solutions (between feed stages 8 and 10) were found to be somewhat different (also shown in Figure 1a). The temperature profile for the multiple solutions (high conversion and low conver-

sion) is shown in Figure 1b. In Figure 2, the variation of the fractional conversion of isobutene with methanol feed location is shown as a function of temperature dependence of the equilibrium constant. Three cases are considered: the Rehfinger-Hoffmann expression (1990); the Colombo et al. (1983) expression; and the expression based on the Gibbs free-energy data from the DIPPR Data Bank (Daubert and Danner, 1986). The temperature dependence of the equilibrium constant for the three expressions is shown in Figure 3. It can be noted from Figure 2 that while the Rehfinger-Hoffmann expression yielded multiple solutions, with the Colombo et al. expression or the DIPPR Data Bank based expression, multiple solutions could not be found. From Figure 3, it can be noted that the differences between the three expressions increases with increasing temperature, while from Figure 2, it can be noted that the differences between the different solutions (considering only the high conversion solution) lies in stages 10-14. From Figure 1b, however, it can be noted that these stages correspond to the “higher” temperatures where the differences between the expressions for the temperature dependence of the equilibrium constant increases as the temperature increases. In Figure 4, composition-phase diagrams computed with the Rehfinger-Hoffmann expression confirm the existence of an element azeotrope (as defined by Pe´rez-Cisneros et al. (1997)). Using the Colombo expression, however, the element azeotrope could not be found. Therefore, for the MTBE process, it appears that the choice of the temperature-dependent equilibrium constant expression is very important, as it affects the simulation results significantly. The esterification steady-state and dynamic simulations have been performed with ASPEN-PLUS and DYNSIM (Gani et al., 1992), respectively, using various model combinations. The steady-state simulation results obtained with ASPEN-PLUS are shown in Figure

3192 Ind. Eng. Chem. Res., Vol. 36, No. 8, 1997

Figure 2. Fractional conversion of isobutene as a function of the methanol feed location using three different expressions for Ka(T).

Figure 1. (a, top) Comparison between different thermodynamic models for the MTBE case. (b, bottom) Temperature profiles corresponding to high and low conversions for the MTBE system.

5, and, the corresponding steady-states obtained through dynamic simulations with DYNSIM are shown in Figure 6. From these figures, an important point to note is that the experimental data of Suzuki et al. (1971) are quite difficult to match unless the “fitted” activity coefficient model (modified Margules in Figure 6c) proposed by Suzuki et al. (1970) is employed together with the corresponding expressions for the reaction kinetics. With ASPEN-PLUS, the best match of the Suzuki et al. (1971) data could be achieved with the combination of UNIFAC for the liquid phase, the Virial equation for the vapor phase (taking into account association and dimerization), and the Izarraraz et al. (1980) rate of reaction expression (Figure 5a). The CPE program package did not find any reactive azeotrope for this system. Level 2: Relationships. Model/Parameter Sensitivities. For the MTBE process, the activity coefficient model parameters have been found to be sensitive

Figure 3. Comparison between different data for the equilibrium constant.

to the low conversion solution. Figure 7a shows the equilibrium liquid lines using the UNIQUAC model with the parameters reported by Rehfinger and Hoffmann (1990) and the Wilson model with the parameters reported by Ung and Doherty (1995). The main differences in the behavior with respect to the models used are only noted when the composition of element A (isobutene) is small, and they disappear as it increases. However, the differences in the behavior with respect to the choice of the temperature-dependent expression (or the parameters for this expression) for the equilibrium constant appear to be very important for the lowconversion region. Figure 7b shows that this difference is larger when the amount of inert component (1-butene) is small. Note that the low conversion corresponds to higher temperatures and that the deviation between the

Ind. Eng. Chem. Res., Vol. 36, No. 8, 1997 3193

Figure 4. Element azeotrope for the MTBE system with 1-butene as the inert component (properties and model parameters from Rehfinger and Hoffmann (1990)).

Rehfinger and Hoffmann and the Colombo et al. expressions increases as the temperature increases. For this reason, only the low-conversion multiple solutions are not obtained when the Colombo et al. expression is used. For the esterification process, the sensitivity of the model parameters is also very high. If the reaction is described by a correlated reaction kinetic expression, then the corresponding activity coefficient model parameters used together with this correlation must be employed in the simulations. That is, the activity coefficient model parameters estimated purely from vapor-liquid equilibrium data, when employed in simulations involving reacting systems and physical equilibrium, may not produce acceptable results (see Figure 5). Some model combination may, however, provide a good match of the experimental data even though the physical equilibrium model parameters were estimated from physical equilibrium data only. For example, in the esterification system, if dimerization of the acetic acid in the vapor phase is considered together with the UNIFAC model for the liquid phase, a good match of the experimental data is possible (see Figure 5). It can be noted in Figure 5b that there are some stages where the rate of generation of ethyl acetate is approximately zero. Therefore, it is likely that these stages are in chemical equilibrium. Subsequent computations of the rate of change of Gibbs energy (from the known liquidphase compositions) have also confirmed that equilibrium has indeed been attained in these stages. Influence of Phase Diagrams. The phase diagrams of the type shown in Figure 8 (temperaturecomposition diagrams) provide information related to the location of the feed and the influence of the inert composition (e.g., for the MTBE process). The existence of element azeotropes increases the importance of the correct determination of phase diagrams. Figure 9 shows the location of the stage compositions for the high- and low-conversion solutions when the methanol feed stage is 10. It can be noted that the highconversion-solution stage compositions are located in one corner of the triangular diagram (close to the pure inert component 1-butene), while the low conversion solution is located close to pure methanol. Level 3: Analysis for Design and Operation. Models/Parameters and Design. Since the activity

Figure 5. (a, top) Liquid composition profiles of ethyl acetate for different thermodynamic models (using ASPEN-PLUS). (b, middle) Rate of generation of ethyl acetate for different thermodynamic models (using ASPEN-PLUS). (c, bottom) Activity coefficient profiles of acetic acid for different thermodynamic models (using ASPEN-PLUS).

3194 Ind. Eng. Chem. Res., Vol. 36, No. 8, 1997

Figure 7. (a, top) Comparison of the effect of different thermodynamic models for the liquid phase on the computation of CPE. (b, bottom) Comparison of the effect of different chemical equilibrium constant expressions (different thermochemical data) for the liquid phase on the computation of CPE.

Figure 6. Comparison of (a, top) temperature profiles and (b, middle) composition profiles calculated with ASPEN-PLUS and DYNSIM (in steady state, t ) 20 h) and (c, bottom) ethyl acetate composition profiles for different thermodynamic models for the esterification process (using DYNSIM).

coefficient model parameters are sensitive to simulation, they also have significance with respect to design. For example, for the esterification process, the reboiler heat duties needed to obtain the same product specification (amount of ethyl acetate in the top product) using two different model combinations (UNIFAC with SRK EOS for the vapor phase as opposed to UNIFAC with the Virial equation for the vapor phase) differ by more than 50%. Depending on which combination is correct and which combination is used, the design results may lead to wasted energy or off-specification products. The same behavior was also found with DYNSIM when the corresponding activity coefficient model parameters were not used together with the known kinetic expression. Since some stages (1-5) are in chemical equilibrium (see level 2), dividing the column into two reactive zones (a chemical equilibrium zone and a kinetically controlled zone) and repeating the simulation should give approximately the same results as those shown in Figure 5a. Therefore, it is possible to identify the

Ind. Eng. Chem. Res., Vol. 36, No. 8, 1997 3195

diagram) for the corresponding steady-state solution. One point to note, however, is that element azeotropes do not have the same significance as nonreactive azeotropes. This is because the element azeotrope point can be pass through a nonreactive stage. The two multiple solutions for the MTBE case show large differences in the total liquid and vapor flow rates (9302.4 as opposed to 13 910.4 kmol/h). Thus, it may be possible to obtain these multiple solutions with steady-state simulation models which do not consider equipment sizing parameters. With dynamic simulation, however, the flexibility of the equipment design parameters will determine if all the multiple solutions can be obtained. This indicates that even if multiple solutions exist, the equipment design may exclude some of these solutions. The opposite, equipment design introducing new solutions, has not been investigated in this work. Conclusions Figure 8. Reactive phase diagram for the MTBE system at different temperatures (properties and model parameters from Colombo et al. (1983)).

Figure 9. Location of the reactive stage compositions (element mass fractions) corresponding to low and high conversions (simulation results for methanol feed stage 10; properties and model parameters from Colombo et. al. (1983)).

equilibrium behavior of a section of the reactive column through the proposed analysis. CPE and Design/Operation. The existence of element azeotropes is closely related to the final steadystate condition of operation. If element azeotropes are present (as in the case of the MTBE process), the initial condition of the column stage compositions and the feed location determine which of the multiple solutions the column will arrive at. If element azeotropes are not present (as in the case of the esterification process), it is the feed location and not the initial condition that determines the final steady-state condition of operation. In this work, we have only investigated the effects of feed location and initial conditions. Other variables may also lead to the same observations described above. In the case of systems with element azeotropes, however, the feed location and initial condition must match (i.e., they must lie on the correct side of the phase

Through two well-known reactive distillation processes, we have highlighted the results of a systematic analysis related to design and operation of these processes. The analysis shows that there exist a set of sensitive properties and that the selected models for their predictions play an important role in the design/ operation of the distillation columns. The analysis also illustrates how design decisions related to feed location, product specification, and initial conditions can be made and points out the dangers related to design/operation of the reactive distillation columns due to incorrect choice of the models/parameters. Therefore, the basic steps of the proposed analysis approach are also valid for other reactive systems. Multiple solutions appear to be related to the existence of element azeotropes, and the design/operation of the reactive distillation columns are closely related to the phase diagrams of the reactive system. Also, element azeotropes do not pose the same separation limits as nonreactive azeotropes. It should be noted that the above conclusions are not based on the results from the two examples highlighted in this paper but a larger number of reactive systems that have been studied. Current and future work is focusing on the development of a general reactive distillation model and the development of an algorithm which combines design and operation through the element-based phase diagrams. More work is needed to establish the relationships between equipment design parameters and multiple solutions. Also, further investigations are needed to determine the differences (if any) between the number of steady states that can be obtained through steady-state simulation and through dynamic simulation. The results of this investigation will be useful for control/operation of reactive distillation columns. Nomenclature a ) activity b ) element composition vector D ) mass or energy balances E ) total energy holdup on the tray (kJ) Ea ) activation energy (J/kmol) EL ) liquid entrainment flow rate (kmol/h) EV ) vapor entrainment flow rate (kmol/h) F ) feed flow rate (kmol/h) HF ) feed enthalpy (kJ/kmol) H ) vapor enthalpy (kJ/kmol) h ) liquid enthalpy (kJ/kmol) K ) physical equilibrium constant

3196 Ind. Eng. Chem. Res., Vol. 36, No. 8, 1997 K ) element separation factor Ka ) chemical equilibrium constant (evaluated at plate temperature) KD ) forward kinetic rate constant KI ) reverse kinetic rate constant k ) kinetic rate constant (k ) f(T)) L ) liquid flow rate (kmol/h) M ) total molar hold-up (kmol) m ) molar hold-up (kmol) n ) component composition vector NR ) number of chemical reactions NC ) number of components ps ) saturation pressure (kPa) P ) tray pressure (kPa) PL ) liquid product side stream (kmol/h) PV ) vapor product side stream (kmol/h) Q ) heat removed or added (J/h) QR ) heat of reaction (J/h) Rp,i ) reaction term (kmol/h) R ) ideal gas constant (J/(kmol K)) ra ) reaction rate (kmol/h) t ) time (h) T ) plate temperature (K) T0 ) reference temperature (K) V ) vapor flow rate (kmol/h) Wm ) elemental fractions x ) liquid composition (mole fraction) y ) vapor composition (mole fraction) z ) feed composition (mole fraction)

Mass balance around plate p for component i

Dp,i ) Fpzp,i + Vp+1yp+1,i + Lp-1xp-1,i + EVp-1yp-1,i + ELp+1xp+1,i - (PVp + Vp + EVp)yp,i - (PLp + Lp + NR

ELp)xp,i + Rp,i + Total energy balance

DH p ) FpHFp + Vp+1Hp+1 + Lp-1hp-1 + EVp-1Hp-1 + ELp+1Hp+1 - (PVp + Vp + EVp)Hp - (PLp + Lp + ELp)hp + QRp + Qp Steady state

Subscripts AA ) acetic acid Al ) chemical element Al (isobutene) Bl ) chemical element Bl (methanol) Cl ) chemical element Cl (1-butene) E ) ethanol EA ) ethyl acetate EOS ) equation of state i ) component number IB ) isobutene k ) reaction number MeOH ) methanol MTBE ) methyl tert-butyl ether 0 ) reaction constant (evaluated at T0) p ) plate number (top ) 1, bottom ) np) W ) water Superscripts R, β ) kinetic power law coefficient ν ) stoichiometric coefficients L ) liquid phase H ) enthalpy balance V ) vapor phase

Appendix. Reactive Distillation Model Equations Reaction Distillation Model. (1) Balance Equations.

DH p ) 0

Dp,i ) 0 Dynamic

Dp,i )

dmp,i dt

DH p )

dEp dt

(2) Physical Equilibrium.

yp,i ) Kp,ixp,i Kp,i ) f(xp,iyp,i,T,P) (3) Hydraulic Correlations.

Greek Letters R ) element relative volatility γ ) activity coefficients λ ) element potential vector ξ ) extent of reaction π ) parameters F ) molar density (kmol/m3) φ ) fugacity coefficients µ ) chemical potential vector

∑ (ξk,pvk,p)

k)1

∆Pp ) Pp - Pp-1 ELp ) f(∆Pp,Mp,plate geometry) Vp ) f(∆Pp,Mp,plate geometry) Lp ) f(Mp,plate geometry) EVp ) f(∆Pp,Mp,plate geometry) (4) Defined Equations. Kinetic expressions NC

rak,p ) KDk,p

(γixp,iFp) ∏ i)1

NC

Ri,k

- KIk,p

(γixp,iFp)β ∏ i)1

( ) ( )

KDk,p ) KD0k,p exp -

EaD,k RTp

EaI,k

KIk,p ) KI0k,p exp -

RTp

Equilibrium constant expression

Ka )

∏(xiγi)ν

i

Rate of reaction expression NR

Rp,i )

∑ (rak,pνi,k)

k)1

0 ) Kap,k -

∏iap,iν k

k ) 1, ..., NR

Heat of reaction NR

QR p )

∑ (Qk,prak,p)

k)1

i,k

Ind. Eng. Chem. Res., Vol. 36, No. 8, 1997 3197

Literature Cited Agreda, V. H.; Partin, L. R.; Heise, W. H. High-Purity Methyl Acetate via Reactive Distillation. Chem. Eng. Prog. 1990, 86, 40. Aspen Technology Inc. The ASPEN -PLUS User’s Manual, 1995. Barbosa, D.; Doherty, M. F. Design of Multicomponent Reactive Distillation Columns. Ind. Chem. Eng. Symp. Ser. 1988, 194, B369. Colombo, F.; Cori, L.; Dalloro, L.; Delogu, P. Equilibrium Constant for the Methyl Tertiary-Butyl Ether Liquid Phase Synthesis by Using UNIFAC. Ind. Eng. Chem. Fundam. 1983, 22, 219. Daubert, T. E.; Danner, R. P. DIPPR Data Compilation; AIChE: New York, 1986. De Garmo, J. L.; Parulekar, V. N.; Pinjala, V. Consider Reactive Distillation. Chem. Eng. Prog. 1992, 88, 43. Doherty, M. F.; Buzad, G. Reactive Distillation by Design. Trans. Inst. Chem. Eng. 1992, 70A, 448. Fredenslund, Aa.; Gmehling, J.; Rasmussen, P. Vapour-Liquid Equilibria Using UNIFAC; Elsevier: Amsterdam 1977. Gani, R.; Cameron, I. T. Extension of Dynamic Models of Distillation Columns to Steady State Simulation. Comp. Chem. Eng. 1989, 13, 271. Gani, R.; Constantinou, L. Molecular Structure Based Estimation of Properties for Process Design. Fluid Phase Equilibr. 1996, 116, 75. Gani, R.; Ruiz, C. A.; Cameron, I. T. A Generalized Model for Distillation Columns-I. Comp. Chem. Eng. 1986, 10, 181. Gani, R.; Sørensen, E. L.; Perregaard, J. Design and Analysis of Chemical Processes through DYNSIM. Ind. Eng. Chem. Res. 1992, 31, 244. Hauan, S.; Hertzberg, T.; Lien, K. M. Why Methyl Tertiary-Butyl Ether Production by Reactive Distillation May Yield Multiple Solutions. Ind. Eng. Chem. Res. 1995, 34, 987. Hayden, J. C.; O’Connell, J. P. A Generalised Method for Predicting Second Virial Coefficients. Ind. Eng. Chem. Process. Des. Dev. 1975, 14, 209. Ignatius, J.; Ja¨rvelin, H.; Lindquist, P. Use TAME and Heavier Ethers to Improve Gasoline Properties. Hydrocarbon Process. 1995, 2, 51. Izarraraz, A.; Bentzen, W.; Anthony, R. G.; Holland, C. D. Solve More Distillation Problems. Hydrocarbon Process. 1980, 59, 195. Jacobs, R.; Krishna, R. Multiple Solutions in Reactive Distillation for Methyl Tertiary-Butyl Ether Synthesis. Ind. Eng. Chem. Res. 1993, 32, 1706. Kang, Y. W.; Lee, Y. Y.; Lee, W. K. Vapour-Liquid Equilibria with Chemical Reaction. Equilibrium-Systems Containing Acetic Acid, Ethyl Alcohol, Water and Ethyl Acetate. J. Chem. Eng. Japan. 1992, 25, 649.

Laroche, L.; Bekiaris, N.; Andersen, H. W.; Morari, M. The Curious Behavior of Homogeneous Azeotropic Distillation-Implications for Entrainer Selection. AIChE J. 1992, 38, 1309. Naphatali, L. M.; Sandholm, D. P. Multicomponent Separation Calculations by Linearization. AIChE J. 1971, 17, 148. Pe´rez-Cisneros, E. S.; Gani, R.; Michelsen, M. L. Reactive Separation Systems. Part I: Computation of Physical and Chemical Equilibrium. Chem. Eng. Sci. 1997, 52, 527. Perregaard, J.; Pedersen, B. S.; Gani, R. Steady State and Dynamic Simulation of Complex Chemical Processes. Trans. Inst. Chem. Eng. 1992, 70A, 99. Rehfinger, A.; Hoffmann, U. Kinetics of the Methyl Tertiary-Butyl Ether Liquid Phase Synthesis Catalysed by Ion Exchange Resin-I. Intrinsic Rate Expression in the Liquid Phase Activities. Chem. Eng. Sci. 1990, 45, 1605. Saito, S.; Michishta, T.; Maeda, S. Separation of Meta- and ParaXylene Mixture by Distillation Accompanied by Chemical Reactions. J. Chem. Eng. Jpn. 1971, 4, 37. Sawistowski, H.; Pilavachi, P. A. Performance of Esterification in a Reaction-Distillation Column. Chem. Eng. Sci. 1988, 43, 355. Seader, J. D. the Rate-based Approach for Modeling Staged Separations. Chem. Eng. Res. 1989, 41. Sivasubramanian, M. S.; Boston, J. F. The Heat and Mass Transfer Rate-Based Approach for Modeling Multicomponent Separation Processes. Comp. Appl. Chem. Eng. 1990, 331. Suzuki, I.; Komatsu, H.; Hirata, M. Formulation and Prediction of Quaternary Vapour-Liquid Equilibria Accompanied by Esterification. J. Chem. Eng. Jpn. 1970, 3, 152. Suzuki, I.; Yagi, H.; Komatsu, H.; Hirata, M. Calculation of Multicomponent Distillation Accompanied by a Chemical Reaction. J. Chem. Eng. Jpn. 1971, 4, 26. Taylor, R.; Lucia, A. Modelling and Analysis of Multicomponent Separation Processes. AIChE Symp. Ser. 1995, 91, 19. Ung, S.; Doherty M. F. Vapour-Liquid Phase Equilibrium in System with Multiple Chemical Reactions. Chem. Eng. Sci. 1995, 50, 23. Whiting, W. B.; Tong, T.-M.; Reed, M. E. Effect of Uncertainties in Thermodynamic Data and Model Parameters on Calculated Process Performance. Ind. Eng. Chem. Res. 1993, 32, 1367.

Received for review October 11, 1996 Revised manuscript received March 24, 1997 Accepted March 25, 1997X IE9606404 X Abstract published in Advance ACS Abstracts, June 15, 1997.