Exothermic Isomerization in One-Stage Reactive ... - ACS Publications

Exothermic Isomerization in One-Stage Reactive Distillation: Steady-State Behaviour. Floris Sebastiaan Vos, and Costin Sorin Bildea*. Delft University...
0 downloads 0 Views 226KB Size
Ind. Eng. Chem. Res. 2007, 46, 203-210

203

Exothermic Isomerization in One-Stage Reactive Distillation: Steady-State Behaviour Floris Sebastiaan Vos and Costin Sorin Bildea* Delft UniVersity of Technology, Julianalaan 136, 2628 BL Delft, The Netherlands

The steady-state behavior of one-stage reactive distillation processes is analyzed, for an exothermic isomerization reaction with first-order kinetics and a light-boiling reactant. The dimensionless model distinguishes between state variables, control parameters, and system properties. Combination of state multiplicity and feasibility boundaries corresponding to total reflux or total reboil leads to complex steady-state behavior. Singularity theory is used to divide the Damko¨hler number-heat input space into regions with qualitatively different behaviors. Forty-five bifurcation diagrams are presented, exhibiting a maximum of three steady states, five feasibility boundaries, and three disconnected solution branches. Introduction Reactive distillation processes exhibit strongly nonlinear behavior. Nonlinearity might lead to start-up complications or undesired transitions between steady states,1,7 oscillatory behavior,2 and difficult control.3 State multiplicity has been found in a variety of reactive distillation processes, for example, synthesis of methyl tert-butyl ether (MTBE),1,2,4,7-9 ethylene glycol,5,6,19 ethyl tert-butyl ether (ETBE),7 and tert-amyl methyl ether (TAME).1,8,9 Gu¨ttinger and Morari,10 Malone and Doherty,11 Taylor and Krishna,12 and Almeida-Rivera et al.13 review the relevant literature. Several model-systems were considered for systematic investigation of multiple steady-state occurrences in reactive distillation. Gu¨ttinger and Morari10,14 assumed chemical and phase equilibrium in an infinitely long column operated at infinite reflux and applied the ∞/∞ analysis. The same type of behavior was found in an isothermal continuous stirred tank reactor (CSTR) and reactive distillation units having different numbers of stages.1 Rodriguez et al.15,16 analyzed the simplest reactive separation process, the isobaric adiabatic reactive flash. They showed that multiplicity is caused by the interaction between separation and reaction in systems in which the activation energy and the gradient of boiling temperature with composition are sufficiently large. The same system was analyzed by Lakerveld et al.,17 who presented 25 bifurcation diagrams exhibiting a maximum of three steady states and five feasibility boundaries. The two-phase reactor under boiling conditions18 adds more complexity by including a condenser, although the system operates at total reflux. Three necessary conditions for the existence of multiple steady states have been identified: the reactant has to be the light-boiling component, the difference in boiling points of the reactant and the product has to be sufficiently large, and the reaction mechanism has to include some form of self-inhibition. In this paper, we continue the investigation of state multiplicity in model-systems for reactive distillation. We extend the reactive-flash model by including the condenser and the reboiler. More complex behavior is found, because of the feedback effect introduced by the reflux and boilup. We note that one-stage systems behave similarly to real-life, multistage reactive distillation columns, as shown by the study of methyl formate and ethylene glycol synthesis.19 In our investigation, we use dimensionless equations which distinguish the state variables, * To whom correspondence should be addressed. Tel.: +31 15 278 6076. Fax: +31 15 278 5006. E-mail: [email protected].

Figure 1. One-stage reactive distillation.

as liquid and vapor compositions, from system properties, as relative volatility, heat of reaction, or activation energy, and control parameters, as feed rate, reflux, or reboiler duty. The system properties are chosen such that most of the possible behavior patterns are uncovered. As a result, we show 45 bifurcation diagrams that are qualitatively different. This rich behavior arises from the combination of state multiplicity and feasibility boundaries corresponding to total reflux and total reboil. This article is organized as follows. We start by presenting the mathematical model of the one-stage reactive distillation process and introducing the singular points where the number of steady states changes. Then, a dimensionless model is derived and discussed. In the following section, we choose a reaction system characterized by small heat of reaction, moderate activation energy, and high relative volatility of the reactant. We illustrate the qualitative changes undergone by the bifurcation diagrams when the control parameters are modified, and we achieve the classification of the steady-state behavior by computing relevant codimension-two varieties. Finally, we discuss the influence of system properties. The article ends with conclusions. One-Stage Reactive Distillation Model Figure 1 shows the flowsheet of a one-stage reactive distillation process consisting of reactive tray, total condenser, and reboiler. Pure reactant is fed to the reactive tray, with the flow rate F, preheated to temperature T0 and containing the vapor fraction . The first-order irreversible reaction A f B takes place only on the reactive tray, in the liquid phase. The

10.1021/ie060393w CCC: $37.00 © 2007 American Chemical Society Published on Web 12/08/2006

204

Ind. Eng. Chem. Res., Vol. 46, No. 1, 2007

molar holdup of the reactive tray, nt, is constant and assumed to be well-mixed. The condenser duty QC is adjusted by a pressure controller which manipulates the cooling water flow rate. From the condenser drum, liquid is returned as reflux with flow rate L0. In the reboiler, condensing steam provides heat at the rate QR. Two level controllers manipulate the flow rates of distillate and bottom products, D and B, respectively. The reflux rate L0 and reboiler duty QR are set by means of two flow controllers. zF ) 1, x2, and xD are the reactant mole fractions in the feed and two product streams, respectively. Additionally, we assume ideal gas and liquid behavior, vapor pressures following the Antoine equation, and the same heat capacity and heat of vaporization for both components. The steady-state behavior of the process is described by total and component mole balances eqs 1-6, energy balances eqs 7 and 8, and equilibrium relations eqs 9-13.

Total mole balance: V1 - D - L0 ) 0

(1)

F + V2 + L0 - V1 - L1 ) 0

(2)

L1 - V2 - B ) 0

(3)

∆vapH Cp for T0 ) Tbp,A (liquid/vapor feed) (14b)

T0* ) Tbp,A + 

gas

T0* ) Tbp,A +

∆vapH Cp + (T - Tbp,A) Cp Cp 0 for T0 > Tbp,A (only vapor in feed) (14c)

and the vapor pressure is calculated according to Antoine’s law:

(

psat i ) exp Ai -

∆vapH RT

)

(15)

The model unknowns are the compositions xD, x1, x2, y1, and y2; the temperatures TD, T1, and T2; and the flow rates D, B, L1, V1, and V2. The validity of the one-stage reactive distillation model is limited by the boundaries D ) 0 (total reflux) and B ) 0 (total reboil). The set of model parameters can be divided into system-specific (k0, Ea/R, ∆rH, ∆vapH, Ai, Cp), control (L0, T0*, F, QR, P) and design (nt) parameters. The choice of system properties (Table 1) deserves a special

Component mole balance: x D - y1 ) 0

(4)

( )

Ea zFF + y2V2 + xDL0 - y1V1 - x1L1 - x1n k0 exp )0 RT1 (5) t

x1L1 - y2V2 - x2B ) 0

(6)

Table 1. System Properties system property

dimensional variable

pre-exponential factor Arrhenius temperature Antoine coefficients

k0 ) 5.88 ×1010 s-1 TA ) 7 300 K A1 ) 9.784 A2 ) 7.704 ∆vapH ) 20 112 J/mol Cp ) 83.7 J/mol/K ∆rH ) -2 092 J/mol

heat of vaporization heat capacity heat of reaction

dimensionless variable γ ) 30 a1 ) 9.784 a2 ) 7.704 λ)1 β ) 0.1

Energy balance: FCpT/0 + L0CpT1 + V2(CpT2 + ∆vapH) + x1ntk0 Ea r exp ∆ H - V1(CpT1 + ∆vapH) - L1CpT1 ) 0 (7) RT1

( )

L1Cp(T1 - T2) + QR - V2∆vapH ) 0

(8)

Equilibrium relations: sat xDpsat A (T)TD) + (1 - xD)pB (T)TD) - p ) 0

(9)

sat x1psat A (T)T1) + (1 - x1)pB (T)T1) - p ) 0

(10)

x1psat A (T)T1) - y1p ) 0

(11)

sat x2psat A (T)T2) + (1 - x2)pB (T)T2) - p ) 0

(12)

x2psat A (T)T2) - y2p ) 0

(13)

where T0* is defined as

T0* ) T0

for T0 < Tbp,A (only liquid in feed) (14a)

discussion. Because we aim to reveal the possible behavior, we choose parameter values that enhance the complexity. The number of behavior patterns is higher for systems with large heat of reaction, large activation energy, and large relative volatility.15-17 For isomerization reactions of A f B type, low heat of reaction is expected. In practice, reactive distillation units having few trays could be used for systems with large relative volatility. Consequently, we accept the system properties used by Rodriguez et al.,15,16 including the large relative volatility R ≈ 8. For simple reaction-rate law, there can only be one steady state if the reactant is the heavy component.15,16 Therefore, we consider a light-boiling reactant. From the set of model parameters, we choose F as a distinguished one and use it to plot (Figure 2) the dependence of the reactant concentration in the bottom stream, x2, for a vapor-liquid feed ( ) 0.15, T0* ) 280 K) and different values of the reboiler duty QR. The plots, called bifurcation diagrams, are qualitatively different. For QR ) 50 × 106 W and small F, the entire feed is vaporized (B < 0) and the numerical solution has no physical significance. By increasing F, one solution becomes feasible when the total-reboil boundary (B ) 0) is crossed. The bifurcation diagram computed for QR ) 35 × 106 W has different features: the feasible solution appears also at the total-reboil boundary (B ) 0), but the diagram exhibits two turning points, also known as limits or folds. For increasing values of F, the number of steady states found are zero (small

Ind. Eng. Chem. Res., Vol. 46, No. 1, 2007 205

x1l1 - y2V2 - x2b ) 0 Energy balance:

(21)

( )

γθ1 1 + θ1 V1(λθ1 + 1) - l1λθ1 ) 0 (22)

fλθ/0 + l0λθ1 + V2(λθ2 + 1) + x1Daβ exp

l1λ(θ1 - θ2) + q - V2 ) 0

(23)

Equilibrium relations: sat xDpsat A (θ)θ)D + (1 - xD)pB (θ)θ)D - p ) 0

(24)

sat x1psat A (θ)θ1) + (1 - x1)pB (θ)θ1) - p ) 0

(25)

x1psat A (θ)θ1) - y1p ) 0

(26)

sat x2psat A (θ)θ2) + (1 - x2)pB (θ)θ2) - p ) 0

(27)

x2psat A (θ)θ2) - y2p ) 0

(28)

Figure 2. Composition of the bottom product x2 versus feed flow rate F, for different values of the heat duty: (9) total reflux feasibility boundary (D ) 0); (0) total reboil feasibility boundary (B ) 0); nt ) 5000 mol, T0 ) 280 K, L0 ) 1500 mol/s, and p ) 1 bar.

F), one (after passing the B ) 0 point), three (at the right of the first turning point), and one (at large values of F). Consequently, we will denote this diagram by 0-1-3-1. Reducing the reboiler duty to QR ) 32 × 106 W, two total-reflux feasibility boundaries (D ) 0) appear, and the bifurcation diagram consists of two disconnected branches. In this case, the multiplicity pattern is 0-1-3-2-1. Finally, the bifurcation diagram corresponding to QR ) 30 × 106 W shows only a single steady state, which appears at a total-reflux boundary. For each plot in Figure 2, the turning points and feasibility boundaries correspond to certain values of one parameter, the feed flow rate F. We say that these are codimension-one bifurcation or singular points. The projection of the singular points onto the parameter space is called the bifurcation set. When the singular set is crossed, the number of steady states changes by two (at folds) or one (at feasibility limits). We also note that the type of bifurcation diagram depends upon the number and relative location of turning points and feasibility boundaries. In this article, we will divide the space of model parameters into regions with qualitatively different bifurcation diagrams. However, aiming at general results, we will first write the model into a dimensionless form. This ensures that the number of parameters is kept to a minimum, while the type of behavior predicted by the model is preserved. Equations 1-13 can be rewritten in terms of dimensionless variables and parameters as follows:

where

θ0 * ) θ0 θ0* )

for θ0 < θbp,A ) 0 (only liquid in feed) (29a)

 λ

for θ0 ) θbp,A ) 0 (liquid/vapor feed) (29b)

gas

θ/0

1 Cp ) + θ λ Cp 0 for θ0 > θbp,A ) 0 (only vapor in feed) (29c)

and the vapor pressure is given by

(

psat i (θ)θj) ) exp ai

V1 - d - l0 ) 0

(16)

state variables: xi, θi )

Ti - Tref (i ) D, 1, 2), yk, Tref

Vk L1 D B (k ) 1, 2), d ) ,b) , l1 ) Fref Fref Fref Fref

Ea -∆rH ;γ) ; CpTref RTref CpTref A1 Cp ) vap ; ai ) Ai; b1 ) 1; b2 ) λ) RA1 ∆ H A2

system properties: β )

f + V 2 + l 0 - V1 - l 1 ) 0

(17)

l 1 - V2 - b ) 0

(18)

control parameters: f ) Component mole balance: xD - y 1 ) 0

zFf + y2V2 + xDl0 - y1V1 - x1l1 - x1Da exp

(30)

The model contains the following:

Vk )

Total mole balance:

)

1 - bi + θj 1 + θj

(19)

( )

γθ1 )0 1 + θ1 (20)

L0 λQR F ; l0 ) ;q) ; Fref Fref FrefCpTref

Da )

k(Tref)nt / T0/ - Tref p ; θ0 ) ;p) Fref Tref pref

The boiling point of component A is used as the reference temperature. For an existing design, the nominal value of the feed flow rate and pressure would be the usual choices of the reference flow rate and pressure, respectively. Note that the

206

Ind. Eng. Chem. Res., Vol. 46, No. 1, 2007

dimensionless model is written such that the operational parameters F, L0, QR, ntk(Tref), T0*, and P enter only in one dimensionless parameter (f, l0, q, Da, θ0*, and p, respectively). Classification Of The Steady-State Behavior As previously mentioned, we choose the feed flow rate f as the distinguished parameter. We will divide the Da - q space into regions with different x2 vs f bifurcation diagrams. Singularity theory20 states that the qualitative features of the bifurcation diagram may change only when the parameter set crosses the hysteresis, isola, or double-limit varieties. Only the hysteresis variety is found for the one-stage reactive distillation. When the hysteresis variety is crossed, the number of possible steady states changes by two, as two turning points appear or disappear. The defining conditions of the hysteresis variety are presented, for example, by Kuznetsov.21 When feasibility boundaries exist, the bifurcation diagram may also change at special sets of parameters.22 In this work, we considered the following: • The boundary-limit set (BL) , a turning point occurs at a feasibility boundary. • The boundary-tangent set (BT), the bifurcation diagram is tangent to the feasibility boundary. After crossing this set, the bifurcation diagram changes by appearance/disappearance of two feasibility boundaries. • The cross-and-limit set (CL), the position of one turning point relative to one solution located at the feasibility boundary changes. • The corner set (C), the bifurcation parameter and the state variable are, simultaneously, at their boundaries. • The double-cross set (DC), the relative locations of two feasibility boundaries changes. The above-mentioned varieties are said to be of codimension two because their defining conditions can be solved for the state variables and two parameters (for example, f and q). When labeling codimension-two singular points, we will use the subscripts d and b to denote the feasibility boundary involved. A second, numeral subscript will distinguish between the different varieties of the same type. For tracing the bifurcation diagrams and the loci of singular points, we used a continuation algorithm based on local parametrization23 with the Newton method as corrector. The derivatives needed for calculating the singular points were obtained using the symbolic toolbox from Matlab. Starting values for continuation of codimension-two points were taken from the codimension-one loci. The corner set at f ) d ) b ) 0 was a practical starting value for continuation of the feasibility boundaries. Boundary limit points on b ) 0 or d ) 0 loci provided good initial guesses for the continuation of folds. Singular Set at Large Da. In this section, we discuss the locus of codimension-one singular point at large Damko¨hler number, Da ) 0.3, where the one-stage reactive distillation behaves similarly to the reactive flash.17 This similarity is explained by the dominant role played by the reaction (compared to separation), which applies in both cases. Figure 3 presents the projection of codimension-one singular point in the f - q space (bifurcation set). Because f - x bifurcation diagrams are obtained for a given value of the q parameter, the number and relative location of singular points can be obtained by the intersection of q ) constant horizontal lines with the locus of codimension-one singular points. At large q, there is only one intersection with the boundary b ) 0, and the bifurcation diagram is similar to the case QR ) 50 × 106 W in Figure 2 or diagrams (1) and (3) in Figure 4.

Figure 3. f - q bifurcation set; T stands for turning point; d ) 0 and b ) 0 are the feasibility boundaries; H stands for the cusp; BT stands for boundary tangent; CL stands for cross-and-limit; and BL stands for boundary limit. Qualitatively different bifurcation diagrams corresponding to regions (1)-(10) are shown in Figure 4: Da ) 0.3; θ/0 ) 0.13245; l0 ) 0.3; and p ) 1. Physical properties are given in Table 1.

Figure 4. x2 - f bifurcation diagrams at large Da. Diagrams (1)-(10) correspond to regions (1)-(10) in Figures 3 and 5.

When q is decreased, a narrow region with three b ) 0 boundaries is encountered. This region is bounded by the local extremes of the q vs f plot, which are the boundary-tangent points BTb1 and BTb2. Diagram (2) is representative for the behavior in this region. Decreasing q further, two turning points are met below the cusp H. State multiplicity becomes possible, as shown in diagram (4) of Figure 4. Note that the hysteresis represents the locus of cusp points. Many changes occur at lower q, a region that is detailed in a close-up. First, the boundarytangent set BTd1 is crossed to region (5), leading to two d ) 0 boundaries. Immediately, the relative locations of a turning point and a d ) 0 feasibility boundary change at the cross-and-limit point CL, leading to diagram (6). Further decreasing q to region (7), the d ) 0 and b ) 0 boundaries move toward f ) 0 and disappear at corner points (not shown). One turning point leaves the feasible domain at the boundary limit BLd1, and the maximum number of steady states becomes two, as in diagram

Ind. Eng. Chem. Res., Vol. 46, No. 1, 2007 207

Figure 5. Steady-state classification of the one-stage reactive distillation model: θ/0 ) 0.13245; l0 ) 0.3; and p ) 1. Physical properties are given in Table 1. H stands for hysteresis; C stands for corner set; BT stands for boundary tangent set; CL stands for cross-and-limit set; BL stands for boundary limit set; and DC stands for double-cross. Bifurcation diagrams corresponding to regions (1)-(10) and (11)-(45) are shown in Figures 4 and 6, respectively.

(8). Finally, the other turning point becomes unfeasible, and the state multiplicity disappears in region (10). We also note a very narrow range of q values, bounded by the boundary-tangent sets BTd2 and BTd3, where three d ) 0 boundaries exist (diagram 9). Steady-State Classification. The steady-state classification is achieved by computing the codimension-two points for different values of an additional parameter, the Damko¨hler number. In this way, the Da - q space is divided into several regions (Figure 5). The different f - x2 bifurcation diagrams existing in each region are shown in Figures 4 and 6. We found a total of 45 different bifurcation diagrams. We are quite confident that these cover the behavior possible for the given set of parameters, although this statement is difficult to prove. In general, at large q all diagrams start at a b ) 0 boundary (similar to region 1). In contrast, diagrams at small q contain a d ) 0 boundary (as in region 10). For intermediate values of the reboiler duty, the steady-state behavior becomes very complex because of the existence of up to five feasibility

Figure 6. x2 - f bifurcation diagrams at small Da. Diagrams (11)-(45) correspond to regions (11)-(45) in Figure 5.

boundaries in one diagram. As mentioned in the previous section, the behavior at large Da is similar to the one found for the reactive flash.17 In addition to this, we should remark the occurrence of two double-cross sets at lower Da values. At DCdd, the relative location of two d ) 0 limits changes. Pair of diagrams 12 and 13 is representative for the change in the type of bifurcation diagram that occurs here. At DCbd, d ) 0 and b ) 0 feasibility limits exchange their position, as shown by diagrams 15 and 16. We also note diagrams 19-28, which

208

Ind. Eng. Chem. Res., Vol. 46, No. 1, 2007

Figure 8. Locus of the cusp points in the β - R space, for different values of the activation energy γ; Da ) 0.3; q ) 0.35. Physical properties are given in Table 1. The dashed line represents the locus of cusp points occurring at the feasibility limit b ) 0.

Figure 7. f - R bifurcation set: T stands for turning point; d ) 0 and b ) 0 are the feasibility boundaries; H stands for the cusp; BT stands for boundary tangent; CL stands for cross-and-limit; and BL stands boundary limit; θ/0 ) 0.13245; l0 ) 0.3; and p ) 1. Physical properties are given in Table 1.

consist of three disconnected branches. It should be remarked that, in many of the bifurcation diagrams, the solution branch that corresponds to high product concentration in the bottoms stream (small x2) has a limited extent or is interrupted by several unfeasibility domains. This domain of small concentrations is important, because it corresponds to high-purity product. Discussion When drawing Figures 3-6, a high value for the relative volatility R was used, in order to be able to reveal the whole complexity as previously explained. The effect of lower relative volatility is investigated in Figure 7, which shows the bifurcation set in f - R coordinates. In the first diagram, the cusp H is located at R ≈ 3.5. Therefore, state multiplicity is possible even at low R values, although the absence of the d ) 0 feasibility boundary simplifies the behavior. In the second diagram, obtained at lower Da number, the cusp is located at a higher value of R. Below the cusp, the bending of the b ) 0 locus leads to the occurrence of boundary tangent points. The result is bifurcation diagrams consisting of two disconnected branches because of three feasibility boundaries. These can be better seen

in the last diagram, obtained at even lower Da value. Also, we note that the loci of singular points are confined to a narrower range of f values. Figure 8 shows the location of the cusp point in the heat of reaction (β)-relative volatility (R) space, for several values of the activation energy γ. The dashed line corresponds to the cusp point occurring at the b ) 0 feasibility boundary. Each γ ) constant line divides the β - R plane into two regions: at large values of R, state multiplicity occurs, while the low-R region corresponds to state unicity. Figure 8 demonstrates that multiple steady states are possible for small values of the relative volatility, if the heats of reaction and/or activation energy are sufficiently large. It should be remarked that the maximum values of the dimensionless variables used in Figure 8 correspond to dimensional values of practical significance, Ea ) 105 J/mol and ∆rH ) 2.09 × 105 J/mol. Note that, even if the cusp is unfeasible (dotted line, b < 0 region), there is a region in the Da - q space where the fold points are feasible (for example, region 29 of Figure 5). Although the feasibility boundaries b ) 0 and d ) 0 are still present at low R, the extent of the parameter ranges for which the bifurcation diagrams show three or five feasibility boundaries is very limited. Moreover, we repeated some of the calculations for θ0 ) 0 ( ) 0, saturated liquid) and θ0 ) -0.132 ( ) -0.15, subcooled liquid) and observed only a minor influence on the multiplicity behavior. The calculations presented here were also performed for other values of the reflux l0. We observed similar behavioral trends, which seem to be governed by net amount of heat supplied, q - l0. Comparing the one-tray reactive distillation and the reactive flash, we can make the following observations. The corner set is a codimension-two variety that is important for steady-state classification. Above and below the corner set, the bifurcation diagrams start with a b ) 0 and a d ) 0 point, respectively. For the reactive flash, the corner set occurs for heat input q ) 0. For the one-tray column, the corner set is located at q ) l0, which signifies that the heat input in the reboiler is balanced by the heat removed in the condenser. At large values of Da, the type of bifurcation diagrams and the order in which they appear for different values of q is the same for both systems. For the reactive flash, diagrams with five feasibility boundaries were found only at large values of the activation energy and heat of reaction. For the one-tray column, such boundaries occur also for slightly exothermal reactions with small activation energy. We would like to make clear that there are few reactions that would be carried on in a reactive distillation column with few trays. One candidate could be butene dimerization,24 a very exothermal reaction that involves one reactant and one product,

Ind. Eng. Chem. Res., Vol. 46, No. 1, 2007 209

with very high relative volatility of the reactant. The possibility of state multiplicity in reactive distillation is well-known. The results presented in this paper show that regions of unfeasibility are equally important, occurring for quite a large range of model parameters. In practice, the d < 0 model prediction means that no liquid is left in the condenser drum and the column “dries”. Similarly, b < 0 means no liquid left in the reboiler. These unfeasible regions tend to occur at high purity of the bottom product stream, which is of interest in practice. This behavior must be taken into account during design by considering not only the nominal operating point but also changes of the operating parameters or inaccuracies of the kinetic and thermodynamic data. During operation, one could guard against losing the level in the reboiler or condenser drum by implementing adequate control, for example, override schemes. Conclusions The steady-state behavior of a one-stage reactive distillation process carrying an exothermic isomerization reaction with firstorder kinetics and light-boiling reactant is analyzed. Steadystate classification in the Damko¨hler number-reboiler duty space is achieved by calculating relevant singularities. At large Da values, reaction dominates and the behavior is similar to the reactive flash.17 At intermediate and small Da values, the separation becomes important. Because of the feedback effect introduced by reflux and boilup, more complex behavior was observed, such as bifurcation diagrams with up to three disconnected solution branches and five total reboil and total reflux feasibility boundaries. In general, high heat of reaction, large relative volatility, and large activation energy enlarge the multiplicity region. For quite a large range of the control parameters, the solution branch corresponding to high product concentration is narrow or interrupted by unfeasibility domains. These conclusions can be extended to reactive distillation columns containing several trays. Acknowledgment Financial support of the European Commission (Research Training Network MRTN-CT-2004-512233, “Towards Knowledge-Based Processing Systems”) is gratefully acknowledged. Nomenclature Ai ) Antoine coefficient of component i ai ) dimensionless Antoine coefficient of component i B ) bottom product flow rate (mol s-1) b ) dimensionless bottom product flow rate bi ) coefficient in dimensionless Antoine equation of component i Cp ) heat capacity of liquid (J mol-1 K-1) gas Cp ) heat capacity of vapor (J mol-1 K-1) D ) distillate flow rate (mol s-1) Da ) Damko¨hler number Ea ) activation energy, (J mol-1) F ) feed flow rate (mol s-1) f ) dimensionless feed flow rate P ) pressure (atm) p ) dimensionless pressure k0 ) pre-exponential factor (s-1) L ) liquid flow rate (mol s-1) l ) dimensionless liquid flow rate nt ) molar hold up (mol) Q ) reboiler duty (W)

q ) dimensionless reboiler duty R ) gas constant (J mol-1 K-1) T ) temperature (K) V ) vapor flow rate (mol s-1) V ) dimensionless vapor flow rate x ) liquid molar fraction of component A y ) vapor molar fraction of component A z ) molar fraction of component A ∆rH ) heat of reaction (J mol-1) ∆vapH ) heat of vaporization (J mol-1) Greek Letters R ) relative volatility β ) dimensionless heat of reaction γ ) dimensionless activation energy  ) vapor fraction in feed λ ) dimensionless heat capacity θ ) dimensionless temperature Subscripts 0, 1, 2, D ) feed, reactive tray, reboiler, distillate b ) total reboil d ) total reflux F ) feed bp ) boiling point A, B ) chemical species ref ) reference Superscripts sat ) saturated * ) feed condition Literature Cited (1) Mohl, K.; Kienle, A.; Gilles, E.; Rapmund, P.; Sundmacher, K.; Hoffmann, U. Steady-state multiplicities in reactive distillation columns for the production of fuel ethers MTBE and TAME: theoretical analysis and experimental verification. Chem. Eng. Sci. 1999, 54, 1029. (2) Hauan, S.; Hertzberg, T.; Lien, K. M. Multiplicity in reactive distillation of MTBE. Comput. Chem. Eng. 1997, 21, 1117. (3) Al-Arfaj, M. A.; Luyben, W. L. Comparative control study of ideal and methyl acetate reactive distillation. Chem. Eng. Sci. 2002, 57, 5039. (4) Sneesby, M. G.; Tade´, M. O.; Smith, T. N. Steady-state transitions in the reactive distillation of MTBE. Comput. Chem. Eng. 1998, 22, 879. (5) Ciric, A. R.; Miao, P. Steady state multiplicities in an ethylene glycol reactive distillation column. Ind. Eng. Chem. Res. 1994, 33, 2738. (6) Kumar, A.; Daoutidis, P. Modeling, analysis and control of ethylene glycol reactive distillation column. AIChE J. 1997, 45, 51. (7) Sneesby, M. G.; Tade´, M. O.; Smith, T. N. Multiplicity and pseudomultiplicity in MTBE and ETBE reactive distillation. Chem. Eng. Res. Des. 1998, 76 (A2), 525. (8) Baur, R.; Taylor, R.; Krishna, R. Dynamic behaviour of reactive distillation columns described by a nonequilibrium stage model. Chem. Eng. Sci. 2001, 56, 2085. (9) Chen, F.; Huss, R. S.; Doherty, M. F.; Malone, M. F. Multiple steady states in reactive distillation: Kinetic effects. Comput. Chem. Eng. 2002, 26, 81. (10) Gu¨ttinger, T. E.; Morari, M. Predicting multiple steady states in equilibrium reactive distillation. 1. Analysis of nonhybrid systems. Ind. Eng. Chem. Res. 1999, 38, 1633. (11) Malone, M. F.; Doherty, M. F. Reactive Distillation. Ind. Eng. Chem. Res. 2000, 39, 3953. (12) Taylor, R.; Krishna, R. Modelling reactive distillation. Chem. Eng. Sci. 2000, 55, 5183. (13) Almeida-Rivera, C. P.; Swinkels, P. L. J.; Grievink, J. Designing reactive distillation processes: present and future. Comput. Chem. Eng. 2004, 28, 1997. (14) Gu¨ttinger, T. E.; Morari, M. Predicting multiple steady states in equilibrium reactive distillation. 2. Analysis of hybrid systems. Ind. Eng. Chem. Res. 1999, 38, 1649.

210

Ind. Eng. Chem. Res., Vol. 46, No. 1, 2007

(15) Rodriguez, I. E.; Zheng, A.; Malone, M. F. The stability of a reactive flash. Chem. Eng. Sci. 2001, 56, 4737. (16) Rodriguez, I. E.; Zheng, A.; Malone, M. F. Parametric dependence of solution multiplicity in reactive flash. Chem. Eng. Sci. 2004, 59, 1589. (17) Lakerveld, R.; Bildea, C. S.; Almeida-Rivera, C. P. Exothermic isomerization reaction in a reactive flash: Steady-state behaviour. Ind. Eng. Chem. Res. 2005, 44, 3815. (18) Waschler, R.; Pushpavanam S.; Kienle, A. Multiple steady states in two-phase reactors under boiling conditions. Chem. Eng. Sci. 2003, 58, 2203. (19) Gehrke, V.; Marquardt, W. A. singularity theory approach to the study of reactive distillation. Comput. Chem. Eng. 1997, 21, S1001. (20) Golubitsky, M.; Schaeffer, D. G. Singularities and groups in bifurcation theory; Springer-Verlag: New York, 1985.

(21) Kuznetsov, Y. A. Elements of applied bifurcation theory; SpringerVerlag: New-York, 1998. (22) Balakotaiah, V.; Luss, D. Global analysis of the multiplicity features of multi-reaction lumped parameter systems. Chem. Eng. Sci. 1984, 39, 865. (23) Seydel, R.; Hlavacek, V. Role of continuation in engineering analysis. Chem. Eng. Sci. 1987, 42, 1281. (24) Honkela, M.; Krause, O. Kinetic Modeling of the Dimerization of Isobutene. Ind. Eng. Chem. Res. 2004, 43, 3251.

ReceiVed for reView March 29, 2006 ReVised manuscript receiVed October 2, 2006 Accepted October 17, 2006 IE060393W