A disjunctive programming approach for the optimal design of reactive distillation columns

Computers and Chemical Engineering 25 (2001) 1661– 1673 www.elsevier.com/locate/compchemeng

A disjunctive programming approach for the optimal design of reactive distillation columns Jennifer R. Jackson, Ignacio E. Grossmann * Department of Chemical Engineering, Carnegie Mellon Uni6ersity, Pittsburgh, PA15213, USA Received 31 January 2001; received in revised form 3 August 2001; accepted 6 August 2001

Abstract A generalized disjunctive programming formulation is presented for the optimal design of reactive distillation columns using tray-by-tray, phase equilibrium and kinetic based models. The proposed formulation uses disjunctions for conditional trays to apply the MESH and reaction kinetics equations for only the selected trays in order to reduce the size of the nonlinear programming subproblems. Solution of the model yields the optimal feed tray locations, number of trays, reaction zones, and operating and design parameters. The disjunctive program is solved using a logic-based outer-approximation algorithm where the MILP master problem is based on the big-M formulation of disjunctions, and where a special initialization scheme is used to reduce the number of initial NLP subproblems that need to be solved. Two examples are presented that include reactive distillation for the metathesis reaction of 2-pentene and for the production of ethylene glycol. The results show that the proposed method can effectively handle these difficult nonlinear optimization problems. © 2001 Elsevier Science Ltd. All rights reserved. Keywords: Reactive distillation; Disjunctive programming; Optimal distillation column design

Nomenclature The continuous 6ariables for the model constraints are defined as follows: F ik DIS BOT NT STGk LIQk VAPk L ik V ik Hk x ik y ik i R QB

molar feed flow of component i on tray k total molar flow of distillate total molar flow of bottoms total number of trays counter for the existence of a tray total molar flow of liquid out of tray k total molar flow of vapor out of tray k molar liquid flow of component i out of tray k molar vapor flow of component i out of tray k height of tray k liquid mole fraction of component i out of tray k vapor mole fraction of component i out of tray k boil-up fraction reflux fraction reboiler heat load

* Corresponding author. Tel.: + 1-412-268-2230; fax: + 1-412-268-7139. E-mail address: [email protected] (I.E. Grossmann). 0098-1354/01/$ - see front matter © 2001 Elsevier Science Ltd. All rights reserved. PII: S 0 0 9 8 - 1 3 5 4 ( 0 1 ) 0 0 7 3 0 - X

J.R. Jackson, I.E. Grossmann / Computers and Chemical Engineering 25 (2001) 1661–1673

1662

QC ZT R lm ZT C lm HLik HVik H ivap Tk PRjk DHjk K ik xjk Wk

condenser heat load log-mean temperature difference in reboiler log-mean temperature difference in condenser liquid molar enthalpy out of tray k vapor molar enthalpy out of tray k heat of vaporization for component i temperature out of tray k production rate of reaction j on tray k heat of reaction of reaction j on tray k vapor –liquid equilibrium constant for component i on tray k reaction rate of reaction j on tray k liquid holdup on tray k

1. Introduction Process intensification is an area of growing interest that is defined as any chemical engineering development that leads to a substantially smaller, cleaner, and more energy efficient technology (Stankiewicz & Moulijn, 2000). An excellent example of process intensification is reactive distillation, which provides an alternative to conventional processing schemes that include for instance a reactor followed by a distillation column in chemical production by combining reaction and separation (Fig. 1). The reduction in the number of processing units and the direct heat integration between reaction and separation can reduce capital investment as well as utility costs. Increased overall conversion, as well as improved selectivity in competing reactions, can be achieved in reactive distillation by the continuous removal of products from the reaction zone of equilibrium limited reactions (Ciric & Gu, 1994). Other process improvements can be realized such as reducing byproduct formation, improving raw material usage, and overcoming chemical

Fig. 1. Reactor/distillation column integrated by reactive distillation.

equilibrium limitations (Okasinski & Doherty, 1998). Examples of industrial reactive distillation processes include the methyl acetate process developed by Agreda and Partin (1984), Agreda, Partin, and Heise (1990) and commercialized by Eastman Chemical Company, and the methyl tert-butyl ether (MTBE) process invented by Smith (1981) and Smith (1990). Design issues for reactive distillation systems are significantly more complex than those involved in ordinary distillation (Doherty & Buzad, 1992). Catalyst selection, liquid holdup on each tray, and position of feeds become important design considerations. Reaction often occurs in the liquid holdup so that the reaction volume is a major design parameter, and constant molar overflow cannot be assumed. Also, a single feed may not be appropriate and a distributed feed must be considered. Research work in reactive distillation design has focused on either equilibrium chemical reactions or kinetically controlled reactions. The most popular focus has been the design of reactive distillation systems that are near chemical equilibrium. Much of this work utilizes geometric and graphical methods. Barbosa and Doherty (1988) developed an approach that utilizes residue curve maps to analyze ideal and nonideal mixtures. Blessling, Schembecker, and Simmrock (1997) presented a design method using reactive distillation line diagrams based on component concentration profiles. Lee and Westerberg (2000) proposed a graphical method for the visualization of tray-by-tray calculations to gain design insights into binary reactive distillation. The only optimization based design method for equilibrium limited reactions to date has been work by Frey and Stichlmair (2000), which combines simulation with a rigorous mixed-integer nonlinear programming (MINLP) model.

J.R. Jackson, I.E. Grossmann / Computers and Chemical Engineering 25 (2001) 1661–1673

When reaction equilibrium cannot be assumed, parameters such as reaction kinetics, heat of reaction, and liquid holdup must be considered as important design factors. Several methods for kinetically controlled reactive distillation design have been reported over the past several years. Buzad and Doherty (1994) and Buzad and Doherty (1995) developed design methods for ternary reactive distillation systems. They first proposed an approach using a fixed-point method that set certain design parameters based on fixed points of condenser and reboiler composition profiles. Later they presented a design procedure in which chemical reactions are handled implicitly through the Damko¨ hler number (Da), which is the ratio of the liquid residence time to the reaction time. Both methods assume the same liquid holdup on each stage, which may yield suboptimal designs. Okasinski and Doherty (1998) proposed a method to generate feasible column designs over a range of design conditions. Their approach uses both the Damko¨ hler number and the fixed points of both material and energy balances. While these methods provide design insights, they do not provide an optimum column design, nor do they provide many of the essential design and operating parameters. Several optimization methods for the design of nonequilibrium reactive distillation columns have been proposed that utilize mathematical programming models. Ciric and Gu (1994) were first to present a rigorous, tray-by-tray model. Their MINLP model, which is solved with Generalized Benders Decomposition, is similar to the work by Viswanathan and Grossmann (1990) and Viswanathan and Grossmann (1993) for optimizing a conventional distillation column, but ignores effects of liquid enthalpies. The number of trays is optimized in this model by multiplying related constraints with a binary variable that represents the existence of a tray. This introduces bilinearities which complicates the solution and exhibits poor numerical behavior. The other two optimization methods employ a network of modules or blocks as the superstructure representation. Papalexandri and Pistikopoulos (1996) characterize the reactive distillation column as a set of mass and heat transfer modules that integrate a given network. An MINLP model is then formulated from the superstructure that contains all possible connections between modules, and where discrete decisions are modeled with big-M constraints. Smith and Pantelides (1995) uses building blocks of reactive flash vessels networked together to form the column. By using an enumerative procedure for the optimization, the resulting nonlinear programming (NLP) models are each formulated for a fixed number of blocks (trays). Both these methods emphasize a framework where unit operations and structures are not prepos-

1663

tulated (e.g. a full column with trays). However, due to the complexity of reactive distillation, the modules/ blocks basically correspond to trays in the column, and therefore the models are similar to conventional rigorous tray-by-tray models. The objective of this paper is to develop a new optimization model for the rigorous design of kinetically controlled reactive distillation columns. The proposed model is based on the Generalized Disjunctive Programming (GDP) framework developed by Yeomans and Grossmann (2000) for the design and synthesis of distillation columns. The nonlinear trayby-tray model is described, as well as the solution algorithm that is applied to two different problems.

2. Problem statement Given is a set of chemical species that are to be processed by reactive distillation. A subset of reactants and a subset of products are specified. Either the feed rates along with purity requirements or the required production rates must be specified. Also given is a set of chemical reactions with stoichiometric coefficients and rate expressions, vapor–liquid equilibrium data, heats of vaporization, and a maximum number of trays for the column. All reactions are assumed to take place in the liquid phase, which is homogeneous. The problem then consists of determining the optimal number of trays and feed tray locations, the holdup per tray, the reflux and boil-up, and the condenser and reboiler duties. The objective is to minimize the total annualized investment costs plus operating costs.

3. Problem formulation

3.1. Reacti6e distillation column superstructure A reactive distillation column is modeled on a trayby-tray basis, where each tray represents a stage where phase equilibrium and reaction take place, only phase equilibrium takes place, or neither. The simplest form of a reactive distillation column contains a reboiler, a condenser, and a feed tray with some liquid holdup. This case is similar to a conventional reactor with two outlets, and liquid and vapor recycle. Our representation consists of two heat exchangers, the reboiler and condenser which are always active, and a set of conditional trays which can be either active or inactive (Fig. 2). Consider a conditional stage that can correspond to a phase equilibrium tray, a tray with equilibrium and reaction, or no tray at all. In our model we use disjunctions with

1664

J.R. Jackson, I.E. Grossmann / Computers and Chemical Engineering 25 (2001) 1661–1673

are a subset of trays kTC. Therefore, TC is given by TC =NCT @ NRT @ TM. Let FEEDi represent the feed of flowrate component i, PRODi the minimum product flowrate of species i, and  i the minimum purity of species i in the bottoms and/or distillate. We assume specifications in the reboiler only to simplify the presentation of the model. The trays of the column (n= 1, 2, …, N) are numbered from bottom to top, so that the reboiler is tray number 1 and the condenser is tray number N. The objective function in Eq. (1) minimizes the total annualized cost ANNCOST of equipment and utilities. ANNCOST consists of the utility costs along with the annualized column cost f(HT, DIA, AR, AC), which is a function of column height HT which includes reaction holdup, column diameter DIA, reboiler area AR, and condenser area AC. min ANNCOST =f(HT, DIA, AR, AC)+ CSQB + CCWQC.

Fig. 2. Conditional tray framework for reactive distillation column.

two terms, one for existing trays and one for nonexisting trays. An existing conditional tray is an equilibrium stage where reaction can possibly take place. Mass-transfer on each of these trays is modeled with the MESH equations, which include material balance, vapor–liquid equilibrium, mole fraction summation, and enthalpy balance constraints. Feed is allowed on all existing trays. To consider reaction, constraints that account for reaction kinetics and thermodynamics are included. Production is a function of reaction volume so that reaction takes place only when there is positive liquid holdup. A nonexisting or inactive tray is one where no mass transfer and no reaction take place, and consists solely of a simple input– output operation. Mass and energy balances that model these trays equate inlet and outlet flows and enthalpies for liquid and vapor streams. No feed is allowed on inactive trays.

3.2. GDP model In this section we present the detailed Generalized Disjunctive Programming model for the reactive distillation column representation. C is the set of components i, and R is the set of reactions j in the system. TC is the set of all trays k in the column. NCT is the tray k  TC corresponding to the column condenser. NRT is the tray k  TC corresponding to the reboiler. TM is the set of conditional trays which

(1)

Areas are calculated using the Chen approximation (Chen, 1987) for the log-mean temperature differences, heat duties, and overall heat transfer coefficients U R and U C for the reboiler and condenser, respectively, and are substituted in the objective function Eq. (1). The global constraints include the following overall column mass balance, the specifications for minimum production and purity requirements, PRODi,  i, the distribution of feeds, the calculation of the number of trays, NT, column diameter based on flooding correlations, column height, reboiler and condenser areas, and the summation of mole fractions: B i ] PRODi x ik ]  i

(2)

iC

iC, kNRT

FEEDi = %k  TMF ik

iC

NT = %k  TCSTGk DIA ] f(VAPk )

kNRT @ NCT @ TM

HT = %k  TMHk AR = QR/U RZT R lm AC = QC/U CZT C lm %i  Cx ik = 1

kTC

%

kTC

i iC k

y =1

It should be noted that the last two equations can also be written in the more common form as Si  Cx ik − Si  Cy ik = 0. However, in our experience we

J.R. Jackson, I.E. Grossmann / Computers and Chemical Engineering 25 (2001) 1661–1673

found better numerical performance with the two equations which is probably due to the particular pivoting strategy used in CONOPT. The reboiler and condenser trays are permanent stages and are fixed in the column configuration. Therefore, their tray constraints are valid for any column configuration. The constraints for the reboiler and condenser are given in Eqs. (3) and (4), respectively. They include mass and energy balances, the calculation of the

1665

vapor–liquid equilibria, equilibrium constant K ik, kinetic rate xjk, production rate of reaction j, and tray height. The production rate, PRjk, is computed from the liquid holdup, Wk, and kinetic rate relations. Tray height, Hk, is the minimum tray spacing H0 plus additional spacing required to accommodate the liquid holdup, Wk, needed for reaction. The stage counter STGk is set to one. For simplicity in notation we omit the condition i C in the equations below.

F ik +L ik + 1 −B i − V ik +%j  Rw ijPRjk =0  à à V ik =VAPky ik à à à B i =(1−i)L ik à à B i =BOTx ik à à à i BOT =%i  CB à à à y ik = K ikx ik à à i K k =f(Tk ) à à à QB =i%i  C LIQkH ivap à à Ìk NRT i H vap = f(Tk ) à à à %i  CF ikSHFi + %i  CL ik + 1HLik + 1 −L ikHLik −V ikHVik + iL ikHVik − %j  R DHjkPRjk = 0 à à à i HLk =f(Tk ) à à à HVik =f(Tk ) à à DHjk =f(Tk ) à à à xjk = f(x ik, Tk ) à à PRjk =Wkxjk à à Hk =H0 + (1.27Wk )/DIA2 à à ŠDT R =f(T ) lm k L ik =LIQkx ik

STGk =1

(3)

1666

J.R. Jackson, I.E. Grossmann / Computers and Chemical Engineering 25 (2001) 1661–1673

F ik +V ik − 1 −L ik −D i +%j  Rw ijPRjk = 0  à à i i V k =VAPky k à à à D i =(1 −R)Vk à à i i D =DISy k à à à DIS = %i  CD i à à à i i i y k =K kx k à à K ik =f(Tk ) à à à QC =%i  CV ikH vap,i à k à ÌkNCT H vap,i =f(Tk ) à k à à %i  CF ikSHFi + %i  CRVikHLik +V ik − 1HVik − 1 −L ikHLik − V ikHVik − %j  R DHjk PRjk = 0 à à à HLik =f(Tk ) à à i à HVk =f(Tk ) à à DHjk =f(Tk ) à à à xjk =f(x ik, Tk ) à à PRjk =Wkxjk à à Hk =H0 +(1.27Wk )/DIA2 à à ŠDT C =f(T ) lm k L ik =LIQkx ik

(4)

STGk =1

The following mass and energy balances, along with the heat of reaction calculation, are always valid for

conditional trays (see Fig. 2) regardless of whether or not a tray exists:

J.R. Jackson, I.E. Grossmann / Computers and Chemical Engineering 25 (2001) 1661–1673

 à F +L +V − L −V + % w PRjk =0 à jR à à LIQk = % L ik à iC à à i VAPk = % V k à iC Ìk TM % F ikSHF i + L ik + 1HL ik + 1 +V ik − 1HV ik − 1 −L ikHL ik − V ikHV ik − % DHjkPRjk = 0 à à iC jR à i HL k =f(Tk ) à à i HV k =f(Tk ) à à DHjk =f(Tk ) à Ši k

i k+1

i k−1

i k

i k

1667

i j

The aforementioned constraints involve only continuous variables and hold for any column configuration. The disjunctions in Eq. (6) are required for the discrete decision to enforce mass transfer and possible reaction in an existing tray. We introduce the Boolean variable Zk to represent this discrete choice. When a tray is selected, Zk takes the value of true and the vapor and liquid component mass balances are applied along with constraints that calculate the vapor– liquid equilibria, equilibrium constant, kinetic rate and production rate of reaction j, and tray height. The stage counter is set to one. If Zk takes the value of false, the composition and the liquid and vapor flows entering the tray are set equal to those exiting the tray. Similarly, the temperature is set equal to the temperature of the tray above. Feeds, liquid holdup, height, and the stage counter are set to zero. Æ ¬Z Ç k ÃÃ Æ Ç Ã Zk à à à x ik = x ik + 1 à à L ik =LIQkx ik Ãà à à à i i à à à y k = y k − 1 Ãà i i V k =VAPky k à à à L i = L i Ãà k+1 Ãà à à à k y ik =K ikx ik i à à ÃV k = V ik − 1 à à à Ìk TM à à –Ã K ik =f(Tk ) F ik = 0 à à à à à xjk =f(x ik, Tk ) Ãà à à à Tk = Tk + 1 à à à à à PRjk =Wkxjk à à à Wk =0 à à Ãà ÃHk =H0 + 1.27Wk /DIA2 à à à à à Hk =0 à à STGk =1 ÃÃ È É Ã È STGk =0 É Å

(6) Multiple solutions with the same objective function are possible because different combinations of conditional trays can be deactivated for the same total number of trays. A logic constraint is introduced to eliminate these duplicate solutions by enforcing the selected trays to be

(5)

activated above the reboiler1 (Eq. (7)). Zk [ Zk − 1

kTM.

(7)

4. Special cases The model can easily be modified to handle special cases. We can account for the addition or removal of heat on each tray by adding the heat duties of hot and cold heat exchangers to the individual tray energy balances as given by the following equation: %i  CF ikSHF + LIQk + 1HLk + 1 + VAPk − 1HVk − 1

(8)

− LIQk HLk − VAPk HVk − %j  R DHjk PRjk − QBk + QCk = 0, where QBk, QCk, are nonnegative variables representing the heat added/removed in tray k. The model can also optimize the location of a single feed tray by introducing binary variables uk to represent feed F ik to each tray k, and by adding to the feed summation equation in (Eq. (2)) the following constraints: F ik − F iUPuk 5 0 %k  TCuk = 1,

i C, kTC,

(9)

uk = 0, 1, kTC

where F iUP is the upper bound of FEEDi if it is a variable, and equal to FEEDi if fixed. The optimization of the feed tray locations, single or distributed feed case, can be improved by including cost penalties in the objective for added investment costs for feed piping. The model can accommodate these costs with the addition of 0–1 variables (uk ) to represent the 1

Trays are numbered from bottom to top of column.

J.R. Jackson, I.E. Grossmann / Computers and Chemical Engineering 25 (2001) 1661–1673

1668

existence of each feed tray, and by adding the fixed cost term k  TC h ku k to Eq. (1). 5. Solution algorithm The problem given by Eqs. (1)– (7) corresponds to a Generalized Disjunctive Programming (GDP) problem (Raman & Grossmann, 1994). The solution method used is very similar to that of Yeomans and Grossmann (2000). The proposed GDP model is solved with a modification of the logic-based outer approximation (OA) algorithm of Tu¨ rkay and Grossmann (1996). This decomposition algorithm iterates between reduced NLP subproblems and an MILP master problem. The NLP subproblem contains only the equations for the terms in disjunctions that are true and corresponds to a column with the selected trays. The solution of the NLP subproblem provides an upper bound for the objective. The master problem is a linear approximation of the problem and provides a new combination of discrete variables Table 1 Kinetics and physical properties for pentene/butene/hexene system Rate (/h)

0.5kf (x 2C5−(xC4xC6/Keq))

kf (/h) Keq

213216 e(−27633 [kJ/kmol]/RT) exp(−DG/RT)

corresponding to a column with a different number of trays. The MILP master problem is constructed by applying big-M constraints to the disjunctions, instead of using the convex hull formulation proposed in the Tu¨ rkay and Grossmann algorithm. As discussed in Yeomans and Grossmann (2000), the use of big-M transformations is justified because: (1) the convex hull formulation does not yield a tight lower bound due to the looseness of the linearized constraints in each sub-region; and (2) due to the linkage between variables on each tray with other trays there are many disaggregated variables introduced, which greatly increases the size of the MILP. The linearization of the bilinear terms in the disjunctive model is performed by replacing these terms with convex envelopes in the MILP master problem (McCormick, 1976; Quesada & Grossmann, 1995). Initialization becomes a nontrivial issue as solutions are very different from conventional distillation, and are consequently difficult to predict. All conditional trays can be chosen to be active for the first NLP. This would provide a solution for the maximum number of trays, as well as give the necessary values for the linear approximations of the MILP, however, this NLP itself might be difficult to solve. There are several initialization schemes that may be applied to help alleviate this problem. For instance, Fletcher and Morton (2000) have suggested a scheme for initializing NLP models for distillation columns. However, this method does not specifically address the case of reactive distillation. In this paper we have used the initialization scheme described in the appendix. An upper bound for the maximum number of trays (nmax) to be considered for the column design must be selected. The following algorithm can be applied for selecting nmax: 2.1. Fix all trays greater or equal to nmax to be active. 2.2. Increment nmax by Dn. 2.3. Return to step 1. Else: optimal solution found . 6. Numerical examples

Fig. 3. Optimal column configuration for Section 6.1.

The proposed reactive distillation column model Eqs. (1)–(7) was tested with two simultaneous reaction/separation problems. The first example involves separating butene and hexene produced from pentene. The second example deals with the production of ethylene glycol from ethylene oxide and water. Ideal vapor–liquid equilibrium on each tray is assumed, although the proposed model does not exclude the possibility of using nonlinear thermodynamic models. Both examples were solved on a 500 MHz Pentium III PC, with 256 MB RAM. The models and solution algorithm were coded in the GAMS modeling environ-

J.R. Jackson, I.E. Grossmann / Computers and Chemical Engineering 25 (2001) 1661–1673

1669

Table 2 Computational results for Section 6.1 Model description Discrete variables Continuous variables Constraints OA iterations NLP initialization and subproblems (s) MILP master problems (s) Total CPU time (s) Optimal solution Objective (1000-$ per year) No. of trays Feed locations Reaction zone Column diameter (m) Column height (m) Reboiler duty (kW) Condenser duty (kW) Boilup fraction Reflux fraction Total liquid holdup (L)

30 876 862 3 56 133 189 605.47 27 Trays 12, 13, 17 and 19 Trays 5–25 1.22 36.7 596 408 0.735 0.516 23,742

ment (Brooke, Kendrick, Meerhaus, & Raman, 1997). The solver CONOPT 2.0 was used for the NLP subproblems and the solver CPLEX 6.6.1 was used for the MILP master problems.

6.1. Example 1: reacti6e distillation of pentene, butene, and hexene For the first example, we consider the metathesis reaction of 2-pentene to form 2-butene and 3-hexene studied by Okasinski and Doherty (1998). The reaction, shown in Eq. (10), takes place at atmospheric pressure and can be represented with good accuracy using ideal vapor–liquid equilibrium (Reid, Prausmitz, & Poling, 1987; Yaws, 1991).

Fig. 5. Concentration profile of chemical components in Section 6.1.

Fig. 6. Temperature profile of column in Section 6.1.

2C5H10 U C4H8 + C6H12

(10)

The reaction kinetics and physical property data are given in Table 1. The heat of reaction DH is calculated from the heats of formation of the reaction components. A saturated liquid stream of 120 kmol/h pure pentene is fed to the column. There are 30 conditional trays in the column model giving an upper bound of 32 trays in the column. The purity requirement for both butene and hexene was 95%. Cost data for the column, trays, and utilities were obtained from Peters and Timmerhaus (1991). The

Fig. 4. Liquid holdup profile and feeds for Section 6.1.

J.R. Jackson, I.E. Grossmann / Computers and Chemical Engineering 25 (2001) 1661–1673

1670

Table 3 Reaction data for ethylene gycol system Reaction

Rate (mol/cm3s)

DH (kJ/mol)

i

3.15×109 exp(−9547/T (K)) xC2H4O xH2O

−80

ii

6.3×109 exp(−9547/T (K)) xC2H4O xC2H6O2

−13.1

Table 4 Physical property data for ethylene glycol system Component

K

Ethylene oxide 71.9 exp{5.72[(T(K)−469)/(T(K)−35.9)]} Water 221.2 exp{6.31[(T(K)−647)/(T(K)−52.9)]} Ethylene glycol 77 exp{9.94[(T(K)−645)/(T(K)−71.4)]} Diethylene glycol 47 exp{10.42[(T(K)−681)/(T(K)−80.6)]}

Table 5 Cost data for ethylene glycol system Column cost ($ per year) Heating utility cost Cooling utility cost Ethylene oxide feedstock cost Water feedstock cost Downstream separation cost

15.7 DIA1.55 HT +222DIA1.066 HT0.802 $146.8/kW year $24.5/kW year $43.70/kmol $21.90/kmol $0.15/kmol water in effluent

objective was to minimize the total annualized cost. Because the problem is nonconvex, the stopping criterion for the logic-based OA algorithm was no improvement in the NLP solution (see Viswanathan & Grossmann, 1990). The optimal column configuration, which has a total annualized cost of $6.055×105 per year, is shown in Fig. 3. The computational results and optimal design characteristics are given in Table 2. Again, trays are numbered from bottom to top, so that the reboiler is tray number 1. The liquid holdup profile as well as the feed tray locations is shown in Fig. 4. The holdup per tray stays around 1000 l throughout the reaction zone, however, the holdup values on the four feed trays hit the specified upper bound of 2124 l (75 ft3). The component concentration profiles (Fig. 5) show that the reactant pentene stays at its maximum concentration throughout the reaction zone in the column where it is being converted to the products. This is expected as higher concentrations of butene and hexene occur at the top and bottom of the column, respectively. The temperature profile in Fig. 6 illustrates the near isothermal behavior in the range of the feed trays. This is due to temperature increases caused by the saturated liquid feed, as heat effects from this slightly endothermic reaction are negligible. Many design methods for reactive distillation consider only a column with a single feed. The example was re-optimized with a single feed where the number of trays was fixed at 27 trays (original optimal solution). The optimal column configuration has the all pentene fed into tray 17 with an annualized cost of $618,000 per year, $13,000 more per year than with the distributed feed model. Although the cost penalty is not large in this case, the example illustrates the cost benefits that can be realized by allowing distributed feeds.

6.2. Example 2: production of ethylene glycol from ethylene oxide and water 6ia reacti6e distillation In this example, we examine the reaction of ethylene oxide and water to form ethylene glycol at atmospheric pressure (Eq. (11)). Unwanted diethylene glycol byproduct is produced when ethylene oxide further reacts with ethylene glycol as shown in Eq. (12).

Fig. 7. Optimal column configuration for Section 6.2.

(i)

C2H4O + H2O “ C2H6O2

(11)

(ii)

C2H4O + C2H6O2 “ C4H10O3

(12)

Since this system was reviewed and solved by Ciric and Gu (1994), Smith (1996) and Papalexandri and Pistikopoulos (1996), data was used from these papers whenever possible for comparison purposes. Reaction, physical property, and costing data of the ethylene glycol system are found in Tables 3–5, respectively. Our model optimizes the saturated liquid feeds of ethylene oxide and water. There are 38 conditional trays

J.R. Jackson, I.E. Grossmann / Computers and Chemical Engineering 25 (2001) 1661–1673 Table 6 Computational results for Section 6.2 Model description Discrete variables Continuous variables Constraints OA iterations NLP initialization and subproblems (s) MILP master problems (s) Total CPU time Optimal solution Objective (million-$ per year) No. of Trays Feed locations Reaction Zone Column diameter (m) Column height (m) Reboiler duty (MW) Condenser duty (MW) Boilup fraction Reflux fraction Total liquid holdup (l)

38 1705 1593 5 1270 811 55 min 32 s 16.6895 31 Trays 1 and 11–29 Trays 5–31 1.99 19.5 6.523 6.993 0.934 1.000 1832

1671

in the column model giving an upper bound of 40 trays in the column. A 25-kmol/h-production rate of ethylene glycol is required out of the column bottoms. The objective was to minimize the total annualized cost, and no NLP improvement was the stopping criterion for termination. The optimal column configuration, which has a total annualized cost of $16.90× 106 per year, is shown in Fig. 7. The computational results and optimal design characteristics are given in Table 6. Note that in this example the computational requirements are substantially higher, and that more than 50% of the time is spent in solving the NLP subproblems. Fig. 8 shows the profile of the liquid holdup and the feed tray locations. Of the total ethylene oxide fed to the column, 25% is fed to tray 1 while the remaining 75% is distributed about evenly on trays 11–29 (Figs. 7 and 8). The concentration profiles for each component are shown in Fig. 9. The temperature profile shown in Fig. 10 illustrates the near isothermal behavior in the reaction zone of the column due to reaction heat effects. The exothermic reaction contributes to the temperature increase which leads to the isothermal behavior as opposed to the expected monotonic temperature decrease up the column.

7. Conclusions

Fig. 8. Liquid holdup profile and feeds for Section 6.2.

Fig. 9. Concentration profiles of chemical components in Section 6.2.

A new disjunctive programming model has been proposed for the optimal design of reactive distillation columns which allows the systematic selection of number of trays, feed tray locations, locations of reaction trays, and other major decision parameters. The conditional tray framework effectively handles the disappearance of trays without the introduction of bilinearities, and reduces the size of the NLP subproblems to be solved. Two example problems have been solved in this paper with the proposed method: the metathesis reaction of 2-pentene, and the production of ethylene glycol. The numerical results have shown that although the model involves many nonlinear, nonconvex constraints, the proposed method aids significantly in its convergence. Multiple steady states in reactive distillation systems have been confirmed experimentally (Bravo, Pyhalathi & Jaervelin, 1993; Mohl, Kienle, Gilles, Rapmund, Sundmacher, & Hoffmann, 1999; Rapmund, Sundmacher, & Hoffmann, 1998), and the nature of these systems leads to nonlinear, nonconvex models, which yield local solutions. The identification of different steady states and the assurance of global optimality are both important points for future research, but are beyond the scope of the work for this paper.

1672

J.R. Jackson, I.E. Grossmann / Computers and Chemical Engineering 25 (2001) 1661–1673

Fig. 10. Temperature profile of column in Section 6.2.

Acknowledgements The authors acknowledge the financial support provided by the National Science Foundation Graduate Fellowship and the Center for Advanced Process Decision-making at Carnegie Mellon University.

Appendix A. Initialization of NLP subproblems A simulation can be run first with a different objective function that will provide better initial values for the NLP. Another option is to run multiple NLPs with subsets of trays selected to be active so that all conditional trays will be active a least once. The resulting NLPs are reduced in size, which can make them easier to solve. The following describes the initialization scheme used in Sections 6.1 and 6.2. Section 6.1 was initialized with a simulation run in GAMS. Here, the NLP was solved with all equipment costing constraints removed and the objective function only considered utility costs. The initial values for the simulation were set by the following. The liquid and vapor mole fractions were set to 0.5, and the reflux and boilup fractions were set to 0.6. Liquid and vapor flows were set assuming equimolar overflow and set to 60 and 100 kmol/h for the component and total flows, respectively. Fixing the condenser and reboiler values to their dew point and bubble point temperature, respectively, and using linear interpolation for the rest of the trays initialized temperature. Compositions for these calculations were 0.5 for pentene, 0.5 for butene and 0.0 for hexene in the condenser, and 0.0 for butene and 0.5 for hexene in the reboiler. The following variables were specified to ensure nonzero starting points: DHk, QB, QC, Wk, and xk. After the simulation run, the NLP with the full set of constraints and the original objective

function was run. All conditional trays were chosen to be active for these initial runs. In Section 6.2, the initial NLP, with all conditional trays chosen to be active, was run directly. The liquid and vapor mole fractions were set to 0.5, and the reflux and boilup fractions were set to 0.9. Liquid and vapor flows were set assuming equimolar overflow and set to 100 and 300 kmol/h for the component and total flows, respectively. The temperature for all trays was set to 380 K, which is the approximate average of the lower and upper bounds which are based on the boiling and critical points of the lightest component (ethylene oxide). A feed of 3 kmol/h was used for both feed components on all trays. The following variables were specified to ensure nonzero starting points: B EO, B W, D EO, D W, K ik, x 1k, x 2k, Wk, QB, QC, and DIA. Note that enthalpies, heat loads, and costs for both examples were scaled so that their values remained in the range of O(10 − 1)–O(10). Temperature was also scaled for Section 6.2. References Agreda, V. H., Partin, L. R. (1984). Reactive distillation process for the production of methyl acetate. U.S. Patent No. 4,435,595. Agreda, V. H., Partin, L. R., & Heise, W. H. (1990). High-purity methyl acetate via reactive distillation. Chemical Engineering Progress, 86 (2), 40 – 46. Barbosa, D., & Doherty, M. F. (1988). Simple distillation of homogeneous reactive mixtures. Chemical Engineering Science, 43 (3), 541 – 550. Blessling, B., Schembecker, G., & Simmrock, K. H. (1997). Design of processes with reactive distillation line diagrams. Industrial Engineering and Chemical Research, 36, 3032 – 3042. Bravo, J. L., Pyhalathi, A., & Jaervelin, H. (1993). Investigations in a catalytic distillation pilot plant: vapor/liquid equilibrium, kinetics and mass transfer issues. Industrial Engineering and Chemical Research, 32, 2220 – 2225. Brooke, A., Kendrick, D., Meerhaus, A., Raman, R. (1997). GAMS

J.R. Jackson, I.E. Grossmann / Computers and Chemical Engineering 25 (2001) 1661–1673 Release 2.25 Language Guide, Version 92; Gams Development Corporation. Buzad, G., & Doherty, M. F. (1994). Design of three-component kinetically controlled reactive distillation columns using fixedpoint methods. Chemical Engineering Science, 49 (12), 1947 – 1963. Buzad, G., & Doherty, M. F. (1995). New tools for the design of kinetically controlled reactive distillation columns for ternary mixtures. Computers and Chemical Engineering, 19 (4), 395 – 408. Chen, J. J. (1987). Letter to the editor: comments on improvement on a replacement for the logarithmic mean. Chemical Engineering Science, 42, 2488. Ciric, A. R., & Gu, D. (1994). Synthesis of nonequilibrium reactive distillation processes by MINLP optimization. American Institute of Chemical Engineering Journal, 40, 1479–1487. Doherty, M. F., & Buzad, G. (1992). Reactive distillation by design. Chemical Engineering Research Design, 70 (A5), 448 –458. Fletcher, R., & Morton, W. (2000). Initialising distillation column models. Computers and Chemical Engineering, 23, 1811 –1824. Frey, Th., & Stichlmair, J. (2000). In S. Pierucci, MINLP Optimization of Reacti6e Distillation Columns. European Symposium on Computer Aided Process Engineering, vol. 10 (pp. 115 –120). Amsterdam: Elsevier Science B.V. Lee, J.W., Hauen, S., Westerberg, A. W. (2000). Graphical methods for reaction distribution in a reactive distillation column. American Institute of Chemical Engineering Journal 46 (6), 1218 – 1233. McCormick, G. P. (1976). Computability of global solutions to factorable nonconvex programs —Part I — convex underestimating problems. Mathematical Programming, 10, 147–175. Mohl, K. D., Kienle, A., Gilles, E. D., Rapmund, P., Sundmacher, K., & Hoffmann, U. (1999). Steady-state multiplicities in reactive distillation columns for the production of fuel ethers MTBE and TAME: theoretical analysis and experimental verification. Chemical Engineering Science, 54, 1029–1043. Okasinski, M. J., & Doherty, M. F. (1998). Design method for kinetically controlled, staged reactive distillation columns. Industrial Engineering and Chemical Research, 37, 2821–2834. Papalexandri, K. P., & Pistikopoulos, E. N. (1996). Generalized modular representation framework for process synthesis. American Institute of Chemical Engineering Journal, 42 (4), 1010 – 1032. Peters, M. S., & Timmerhaus, K. D. (1991). Plant Design and

1673

Economics for Chemical Engineers (4th). New York: McGrawHill. Quesada, I., & Grossmann, I. E. (1995). Global optimization of bilinear process networks with multicomponent flows. Computers and Chemical Engineering, 12, 1219 – 1242. Raman, R., & Grossmann, I. E. (1994). Modeling and computational techniques for logic based integer programming. Computers and Chemical Engineering, 18 (7), 563 – 578. Rapmund, P., Sundmacher, K., & Hoffmann, U. (1998). Multiple steady states in a reactive distillation column for the production of the fuel ether TAME II experimental validation. Chemical Engineering and Technology, 21, 136 – 139. Reid, R. C., Prausmitz, J. M., & Poling, B. E. (1987). Properties of Gases and Liquids (4th). New York: McGraw Hill. Smith, E. M. B., & Pantelides, C. (1995). Design of reaction/separation networks using detailed models. Computers and Chemical Engineering, S19, S83– S88. Smith, L. A. (1981). Catalytic distillation process. U.S. Patent No. 4,307,254. Smith, L. A. (1990). Method for the preparation of methyl tert-butyl ether. U.S. Patent No. 4,978,807. Stankiewicz, A. I., & Moulijn, J. A. (2000). Process intensification: transforming chemical engineering. Chemical Engineering Progress, 96 (1), 22 – 34. Tu¨ rkay, M., & Grossmann, I. E. (1996). Logic-based MINLP algorithms for the optimal sythesis of process networks. Computers and Chemical Engineering, 20, 959 – 978. Viswanathan, J., & Grossmann, I. E. (1990). Combined-penalty function and outer-approximation method for MINLP optimization. Computers and Chemical Engineering, 14, 769 – 782. Viswanathan, J., & Grossmann, I. E. (1993). Optimal feed locations and number of trays for distillation columns with multiple feeds. Industrial Engineering and Chemical Research, 32 (11), 2942 –2949. Yaws, C. L. (1991). Calculate liquid heat capacity. Hydrocarbon Processing, 70 (12), 73 – 77. Yeomans, H., & Grossmann, I. E. (2000). Disjunctive programming models for the optimal design of distillation columns and separation sequences. Industrial Engineering and Chemical Research, 39 (6), 1637 – 1648.