A Mixed Integer Nonlinear Programming Formulation for Optimal

Ind. Eng. Chem. Res. 2006, 45, 1373-1388

1373

A Mixed Integer Nonlinear Programming Formulation for Optimal Design of a Catalytic Distillation Column Based on a Generic Nonequilibrium Model Jorge M. Go´ mez,† Jean-Michel Reneaume,*,† Michel Roques,† Michel Meyer,‡ and Xuan Meyer‡ Laboratoire de Thermique, Energe´ tique et Proce´ de´ s (LaTEP), ENSGTI, Rue Jules Ferry, BP 75 11, 64075 Pau Cedex, France, and Laboratoire de Ge´ nie Chimique (LGC), ENSIACET, 118 Route de Narbonne, 31077 Toulouse Cedex 4, France

The objective of this contribution is to propose a mixed integer nonlinear programming (MINLP) formulation for optimal design of a catalytic distillation column based on a generic nonequilibrium (NEQ) model. The use of this NEQ model presents two main advantages: (i) the computation of tray efficiencies is entirely avoided and (ii) the geometrical parameters of the column’s hardware can be optimized. The minimization of the total annualized cost is submitted to three sets of constraints: the model equations, the product specification, and the tray hydraulic equations. The solution strategy for the optimization uses a combination of simulated annealing and sequential quadratic programming. Catalytic distillation of ethyl tert-butyl ether (ETBE) is considered as an illustrative example. The results of the optimization are discussed. Pre- and postoptimal sensitivity analysis is also performed. 1. Inroduction Reactive distillation (RD) is a multifunctional reactor concept which combines a distillative separation with a chemical reaction. More specifically, catalytic distillation (CD) is a process where chemical reaction occurs over a solid catalyst and the ensuing mixture of reactants and products is separated by fractional distillation. The combination of distillation and chemical reaction may result in some advantages over the conventional process (reactor + distillation column), namely, investment cost reduction, high conversion, improved selectivity, or heat integration. The research on the CD domain has focused principally on three aspects: design, experimental study, and simulation. The available approaches for the design of RD or CD can be classified into three main groups:1 those based on graphical/ topological considerations, those based on optimization techniques, and those from heuristic/evolutionary considerations. The second group of design methods, based on optimization techniques, can be classified into three different methodologies (problem formulations): mixed integer non linear programming (MINLP), orthogonal collocation on finite elements (OCFE), and mixed integer dynamic optimization (MIDO). The MINLP optimization of a RD column has already been studied by several authors. Tables 1 and 2 give a concise overview of the different works. Ciric and Gu2 proposed a MINLP optimization for the synthesis of reactive distillation columns with a tray-by-tray model. In this case, the MINLP problem is solved using the generalized bender decomposition (GBD) algorithm. The technique is illustrated with the RD of the ethylene glycol. In 2001, Poth et al.3 presented a MINLP optimization of a kinetically controlled reactive distillation process, based on the concepts of equilibrium stage and tray efficiencies. The outer-approximation/equality-relaxation/augmented-penalty (OA/ER/AP) algorithm is used. Cardoso et al.4 used a simulated annealing-based algorithm to solve the MINLP * To whom correspondence should be addressed. Tel.: +33 559407817. Fax: +33 559407801. E-mail: jean-michel.reneaume@ univ-pau.fr. † LaTEP, ENSGTI. ‡ LGC, ENSIACET.

problem for a reactive distillation column. The considered reaction is the production of ethylene glycol. In 2003, Poth et al.5 proposed an extension of their previous works: the MINLP problem is solved using the general algebraic modeling systems (GAMS) in combination with external functions. Sand et al.6 presented the solution of the integer master problem (IP) by means of three approaches: branch and bound (BB), complete enumeration of the subspace of binary variables, and an interval reduction algorithm. The nonlinear programming (NLP) subproblem is solved by the generalized reduced gradient method, and the column is designed to produce methyl tert-butyl ether (MTBE). In all the previous works, the authors used the equilibrium (EQ) model: component material balances, phase equilibrium equations, summation equations and heat balances (MESH equations). The main drawback with these equations is that the column hardware parameters (the diameter and the tray configuration) are not taken into account directly. Thus, the authors generally used the concept of tray efficiency which is very difficult to estimate for multicomponent reactive systems. The LaTEP laboratory (Pau, France) and the LGC laboratory (Toulouse, France) have been working since 1999 on the reactive distillation issue. Our works have covered the three previously mentioned aspects: design, experimental study, and simulation.9-11 The experimental studies are supported by a tray-catalytic distillation column provided by the French Institute of Petroleum (Institut Franç ais du Pe´trole) at Pau and by a section of a packed column at Toulouse, allowing experimental studies at the pilot plant level. In this contribution, we propose a MINLP formulation for optimal design of a catalytic distillation column based on a generic nonequilibrium (rate based) model. The first objective is to avoid the use of tray efficiency. The second objective is to optimize the geometrical parameters of the column’s hardware. This approach is supported by fundamental mass and heat transfer principles. Of course, the main drawback is that the model is much larger and much more complex. In the following section, the nonequilibrium model is described. In section 3, we present the optimization problem formulation: optimization variables and constraints. The fourth section is dedicated to the two stage solution strategy: the nonlinear problem (NLP) subproblem and the MINLP master

10.1021/ie0504506 CCC: $33.50 © 2006 American Chemical Society Published on Web 01/21/2006

1374

Ind. Eng. Chem. Res., Vol. 45, No. 4, 2006

Table 1. Different Optimal Designs of a Catalytic Distillation Column author

Ciric and Gu, 19942

Frey and Stichmair, 20007

Cardoso et al., 20004

Poth et al., 20013

type of model phase reaction kinetic thermodynamics algorithm example

equilibrium liquid homogeneous ideal GBD ethylene glycol

equilibrium liquid heterogeneous wilson OA/ER/AP methyl acetate

equilibrium liquid homogeneous ideal/wilson SA/N methyl acetate

equilibrium liquid heterogeneous wilson OA/ER/AP MTBE

Table 2. Different Optimal Designs of a Catalytic Distillation Column author

Jackson and Grossmann, 20018

Poth et al., 20035

Sand et al., 20046

this article

type of model phase reaction kinetic thermodynamics algorithm example

equilibrium liquid homogeneous ideal OA ethylene glycol

equilibrium liquid heterogeneous wilson OA/ER/AP methyl acetate

equilibrium liquid heterogeneous wilson BB/GRG MTBE

nonequilibrium liquid heterogeneous UNIFAC SQP/SA ETBE

problem. The example and the results are described in sections 5 and 6. Section 7 is dedicated to a brief comparison of the best optimal solutions using the NEQ and the EQ models. Finally, conclusions and future works are discussed. 2. The Model According to our pilot plant design, the model relies on separation stages (trays) and reactive stages. The separation stages are cross-flow sieve trays, of a single pass (Figure 1). Inside the reaction stages, one finds the catalyst pellet in direct contact with the liquid flow, but not with the vapor flow (Figure 2). The stages (reaction or separation) are numbered from the top to the bottom, including the condenser (subscript 0) and the reboiler (subscript nt + 1). 2.1. Separation Stages. For modeling purposes, there are two different available approaches: the equilibrium stage model and

the nonequilibrium stage model. The EQ model includes the assumption that the streams leaving the stages are at the physical equilibrium. Component material balances, phase equilibrium equations, summation equations, and heat balances for each stage (MESH equations) are solved to give composition, temperature, and flow profiles. The classical way of dealing with equilibrium in nonreactive systems is to incorporate the tray efficiency concept. The prediction and scaling-up of this tray efficiency, however, is very difficult and unreliable in the case of multicomponent separation with reactions.12 On the other hand, the NEQ models are commonly based on the two film theory (Figure 3). The equations are the component material balances and energy balances for each phase together with mass and energy transfer rate equations and equilibrium equations at

Figure 2. Schematic representation of a reactive stage.

Figure 1. Sieve-plate diagram.

Figure 3. Diagram of a nonequilibrium stage.

Ind. Eng. Chem. Res., Vol. 45, No. 4, 2006 1375

the interface. In this model, it is the rate of mass and heat transfer, and not the physical equilibrium, that often limit the separation. The NEQ model can be classified, depending on the way heat and mass transfer rates are calculated, into differential and integral NEQ models. According to the differential model, Maxwell Stefan equations, for example, are integrated on the thickness of the films. On the other hand, the integral model uses mass and heat transfer coefficients to determine the flux at the interface. It is not necessary to determinate the thickness of the film in this case. In this work, the nonequilibrium model for the separation stages is based on the works by Krishnamurthy and Taylor,13-15 who presented a nonequilibrium model for a nonreactive distillation column. This model uses the two film theory and the heat and mass transfer coefficients to determine the flux at the interface (integral NEQ model). For optimization purposes, the integral NEQ model is a good trade off between the complexity of the differential NEQ model and the equilibrium model. 2.1.1. Equations of the Integral Nonequilibrium Model. The assumptions considered in our model are the following: steady state operation, vapor and liquid bulks on each side of the interface are perfectly mixed, there is no accumulation of mass and heat at the interface, the liquid-vapor interface is uniform, the interface is at thermodynamic equilibrium, and the feeds are saturated liquid. For the estimation of the thermodynamic properties, the model uses the BibPhy FORTRAN library, from ProSim S. A. The nonequilibrium model equations for the separation section are as follows, where subscript i refers to the component (i ) 1,..., nc) and the subscript j refers to the nonequilibrium stage (j ) 1,..., nt): Liquid bulk, energy balance l LjHlj - Lj-1Hj-1 - lj - hljaj(Tij - Tlj) - FljHLF j )0

l Ni,jH h i,j ∑ i)1

(3)

Ljxi,j - Lj-1xi,j-1 - fli,j - Ni,j ) 0

(4)

(5)

Liquid bulk, component transfer rate equations nc-1

i ) 1, ..., nc - 1 (6)

nc

-1)0

energy balance component material balance material balance component transfer rate

Vapor Bulk vapor bulk temp vapor bulk comp vapor bulk flow rate vapor comp at the interface

nt n t nc nt ncnt

Interface liquid comp at the interface interface temp

ncnt nt

Total

5nt + 5ncnt

equilibrium equality of energy flux at the interface

Vapor bulk, energy balance v v v v i VjHvj - Vj+1Hj+1 - Fvj HVF j + j + hj aj(Tj - Tj) ) 0 (8)

where vj is the energy transfer rate due to the mass transfer. It is calculated by the following equation: nc

vj

)

v Ni,jH h i,j ∑ i)1

(9)

Vapor bulk, component material balances

Vjyi,j - Vj+1yi,j+1 + Ni,j ) 0

(10)

Vapor bulk, material balance

Vj - Vj+1 + Nt,j ) 0

(11)

Vapor bulk, component transfer rate equations nc - 1

v i aj(yk,j - yk,j ))0 ∑ khik,j

i ) 1, ..., nc - 1 (12)

i yi,j -1)0 ∑ i)1

(7)

(13)

Interface, equilibrium equations

(14)

Interface, equality of energy flux at the interface

vj - lj + hvj aj(Tvj - Tij) - hljaj(Tij - Tlj) ) 0

Lj - Lj-1 - Flj - Nt,j

∑ i)1

nt n t nc nt ncnt

Ki,jxii,j - yii,j ) 0 i ) 1, ..., nc

Liquid bulk, material balance

i xi,j

Liquid Bulk liquid bulk temp liquid bulk comp liquid bulk flow rate component transfer rate

nc

Liquid bulk, component material balances

l i aj(xk,j - xk,j) ) 0 ∑ khik,j

energy balance component material balance material balance component transfer rate

(2)

∫Ni,j daj

k)1

number

k)1

The net loss or gain of component i on stage j due to interface transport is defined as

Ni,j - Nt,jxi,j -

variable

Ni,j - Nt,jyi,j -

nc

Ni,j )

equation

(1)

where lj is the energy transfer rate due to the mass transfer. It is calculated by the following equation:

lj )

Table 3. Nonequilibrium Model: Equations and Variables

(15)

Table 3 presents a concise classification of the variables and the equations of the model. 2.1.2. Mass Transfer Coefficient for Muticomponent Systems. In eqs 6 and 12, the mass transfer coefficients for multicomponent systems are calculated directly from the matrix function of inverted binary mass transfer coefficients as follows:14

[khvj ] ) [Rvj ]-1

(16)

[khlj] ) [Rlj]-1[Γj]

(17)

1376


The elements of the matrix [R] are estimated as follows:

yi,j

Rii,j )

nc

ki.nc,jaj

Ril,j ) -yl,j +

yk,j

∑ k* ik

+

k)1

2.1.3. Heat Transfer Coefficients. In the absence of specific methods to evaluate individual phase heat transfer coefficients (eqs 8 and 15), the Chilton-Colburn analogy is used:14

(18)

ik,jaj

(

v hvj ) (kvav)j(Cp,av )j(Le)j2/3

)

1 1 kil,jaj ki.nc,jaj

(i* l)

(19)

The previous equations are written for the vapor phase (composition y) but can be used for the liquid phase (composition x). Concerning the vapor phase, the binary mass transfer coefficients are evaluated according to the AIChE method:14

kvik,jaj )

(

0.776 + 4.567hw - 0.238Fs +

)

87.319Qlj Vj w

(Scvj )-0.5 v Scvj ) µvj /Fvj Dik,j

(20) (21)

For the estimation of diffusion coefficients in vapor mixtures, the correlation of Fuller16 is used:

Dvik,j ) 1.013 × 10-2(Tvj )1.75

x{Mi + Mk}/MiMk 3

3

P{xϑi + xϑk}2

(22)

In this expression, ϑ refers to the molecular diffusion volume (m3/mol). Concerning the liquid phase, the binary mass transfer coefficients are also evaluated using the AIChE correlation:14

klik,jaj

) 19

700(Dlik,j)0.5(0.4Fs tlj

)

+

0.17)tljLj

hlzw/Qlj

(23) (24)

The clear liquid height (hl) is obtained from the Benett equation.16 The diffusion coefficients in concentrate liquid mixtures are estimated with the correlation of Vignes:16 o xi,j o xk,j Dlik,j ) (Dik,j ) (Dki,j )

(25)

Estimation of diffusion coefficients in dilute liquid mixtures is given by the correlation of Wilke and Chang:16 o ) (7.4 × Dik,j

(φkMk)1/2Tlj -8 10 ) l 0.6 µk,jVi

(26)

( )

Γik,j ) δik,j + xi,j

∂ ln γi,j ∂xk,j

The heat transfer coefficient of the liquid phase is arbitrarily set to 1000hvj in order to keep the liquid saturated, a condition observed in all of the experiments.14 2.2. Reaction Stages. As was said before, according to our pilot plant configuration, there is no vapor-liquid interface in the reactive stages (Figure 2). When a reactive stage is considered, a pseudohomogeneous reaction term is introduced in the mass balances of the liquid bulk. These stages are considered to be stirred cell reactors. For programming conveniences, the reactive stage model is formulated in such a way that the number of variables and equations remains the same as that of the separation stages. The equations for each reactive section j are as follows: Liquid bulk, energy balance l + mcat,j Rj∆Hr,j - FljHLF LjHlj - Lj-1Hj-1 j )0

(27)

As noted above (see eqs 20 and 24), a nonequilibrium model cannot proceed without knowledge of the type of columns and the internal layout (hw and w). Tray type and mechanical layout data, for example, are needed in order to calculated the mass transfer coefficients for each tray.

(29)

(reference for enthalpies, Href ) 0 at 398.15 K, ideal gas, 101325 Pa) Liquid bulk, component material balances

Ljxi,j - Lj-1xi,j-1 - mcat,jRi,j - fli,j ) 0

(30)

Liquid bulk, material balance

Lj - Lj-1 - Flj - mcat,jRj ) 0

(31)

Liquid bulk, component transfer rate equations

xi,j - xii,j ) 0

(32)

Vapor bulk, energy balance v )0 Tvj - Tj+1

(33)

Vapor bulk, component material balances

yi,j - yi,j+1 ) 0

(34)

Vapor bulk, material balance

Vj - Vj+1 ) 0

(35)

Vapor bulk, component transfer rate equations

Ni,j ) 0

In the correlation of Wilke and Chang, the following specific notations are used: V volume molar (cm3 mol); µ viscosity (mPa s). The elements of the matrix of thermodynamic factors are defined by

(28)

(36)

Interface, equilibrium equations

yi,j - yii,j ) 0

(37)

Interface, equality of energy flux at the interface

Tvj - Tij ) 0

(38)

2.3. Condenser and the Reboiler. Nonequilibrium models of multicomponent condensers and reboilers have been devised. However, these models need to solve differential equations rather than algebraic equations. We have not attempted to include this yet. The model considers a total condenser and a partial reboiler. Both the condenser and the reboiler are

Ind. Eng. Chem. Res., Vol. 45, No. 4, 2006 1377 Table 4. Condenser: Equations and Variables equation

variable

number

total material balance component material balance energy balance equilibrium mole fraction summation

distillate (D0) liquid comp condenser duty vapor comp temp Total

1 nc 1 nc 1 3 + 2nc

Table 5. Reboiler: Equations and Variables equation

variable

number

total material balance component material balance energy balance equilibrium mole fraction summation

vapor flow rate vapor comp bottoms flow rate liquid comp temp Total

1 nc 1 nc 1 3 + 2nc

nonreactive stages. The MESH equations, with some modifications for the condenser, are used to describe theses stages. At the level of the column model, both the reflux ratio (TR) and the reboiler duty (QB) are fixed parameters. Tables 4 and 5 present the variable and equation classification. The equations for these sections are as follows (subscript 0 refers to the total condenser and subscript nt + 1 refers to the reboiler): Total material balance

(

V1 - L0 1 +

)

1 )0 TR

Lnt - Vnt+1 - Lnt+1 ) 0

(39) (40)

Component material balances

xi,0 - yi,1 ) 0

(41)

Lntxi,nt - Lnt+1xi,nt+1 - Vnt+1yi,nt+1 ) 0

(42)

Energy balance

(

V1Hv1 - L0 1 +

)

1 l H - QC ) 0 TR 0

LntHnl t - Lnt+1Hnl t+1 - Vnt+1Hnvt+1 + QB ) 0

(43) (44)

Equilibrium equations

yi,0 - Ki,0xi,0 ) 0

(45)

yi,nt+1 - Ki,nt+1xi,nt+1 ) 0

(46)

Model fraction summation nc

(yi,0 - xi,0) ) 0 ∑ i)1

(47)

(yi,n +1 - xi,n +1) ) 0 ∑ i)1

(48)

nc

t

t

2.4. Conclusions. The elaborated model is an integral nonequilibrium stage model (integral rate based model), which uses mass and heat transfer coefficients to determine the flux at the interface. The column hardware design (column diameter D h ; tray configuration hw, w, z) has a direct influence on the mass and heat transfer coefficients (eqs 20 and 23). Thus, the

computation of stage efficiencies has been entirely avoided. Furthermore, the dimension of the trays can be optimized. For a system of nc components and a column of nt nonequilibrium sections (reaction + separation), there is a total of (6 + 4nc + 5nt + 5ncnt) equations and the same number of variables. Considering separation sections, the presented model can be easily extended to a packaged column. The packing is divided into a number of sections, each of which is modeled as a nonequilibrium stage. Mass transfer coefficients and interfacial areas can be calculated using suitable models: Bravo and Fair,17 Onda,18 etc. One can also consider sections with simultaneous separation and reaction: the catalytic pellets act as packing elements. In this case, the reaction section of the presented model disappears and the separation section of the model is modified in order to introduce a pseudohomogeneous reaction term. However, this structure is far from our pilot plant configuration. 3. Optimization Problem Formulation There are two major formulations for mathematical representation of problems involving discrete and continuous variables:19 mixed integer nonlinear programming (MINLP) and general disjunctive programming (GDP) when the logic is represented through disjunctions and propositions.8 In this study, the optimal design of a catalytic distillation column is stated as a MINLP problem: the annualized cost is minimized according to continuous and discrete variables. This minimization is submitted to a set of nonlinear constraints. In the present work, we consider the following hypothesis. The catalytic distillation column has two inlet streams (feeds). This hybrid column is composed of three sections: a rectification section (subscript 1), a reactive section between the feeds (subscript 2), and a stripping section (subscript 3). The total feed flow rates F, partial feed flow rates (feed composition)fi, feed temperatures TLF, operational pressure P, hole diameter dh, and pitch of perforations ph are fixed parameters. 3.1. Optimization Variables. The variables to be optimized, can be classified into three groups: (i) two operational variables, reflux ratio TR and reboiler duty QB; (ii) (6 + 4nc + 5nt + 5ncnt) model variables (see section 2), liquid bulk temperature Tlj, liquid bulk compositions xi,j, liquid bulk flow rate Lj, component transfer rates Ni,j, vapor bulk temperature Tvj , vapor bulk compositions yi,j, vapor bulk flow rate Vj, vapor compositions at the interface yii,j, liquid compositions at the interface xii,j, interface temperature Tij, distillate flow rate Do, condenser duty QC, and bottoms flow rate Lnt+1; (iii) (10 + 2nr2) design variables. This third group of variables can be classified into two subgroups. We first have a set of (4 + nr2) integer variables: number of stages in each section ns,j (j ) 1, 2, 3), number of reactive stages in the second section nr2, and location of the catalytic stages in the second section cl2. The term cl2 is an integer vector with nr2 elements: cl2,i (i ) 1, nr2) refers to the location of catalyst stage i in the second section. Then, the cl2,i value lies between (ns1) and (ns1 + ns2): the catalyst stage is located between the two inlet streams. We then have a set of (6 + nr2) continuous variables: weir height hw, column diameter D h , hole area Ah, active area Aa, liquid flow path length z, weir length w, and catalyst load per stage mcat. The variable mcat is a vector of nr2 elements since this is the value of the number of reactive stages in section 2. One should note that the number of integer variables is a function of the optimization variable nr2 which is the number of reactive stages in section 2. This can be easily handled using

1378


the simulated annealing algorithm (see next section). All these variables are bounded with physical criteria. Bounds are given for each specific example. 3.2. Optimization Constraints. The optimization is submitted to three sets of constraints: the model equations (see section 2), the product specification, and the tray hydrodynamic equations. 3.2.1. Product Specification. The most straightforward way is to specify the purity of one of the products in the bottom or in the top. In our case, we will consider the composition of ethyl tert-butyl ether (ETBE) in the residue. The conversion rate is not constrained since the role of the optimization is to increase it implicitly.

xE,nt+1 - E g 0

[] dh ph

Ah ) 0.907Aa

2

(50)

Aa ) At - 2Ada

uvN,max ) Csbf

x

Fl - Fv Fv

(60)

The capacity parameter is evaluated by means of the Kister and Hass correlation:

Csbf ) 0.37

( ) ()() dh2σ

0.125

Fl

Fv Fl

0.1

s hcl

)

hcl )

( )

996 0.157d0.833p-0.791 l -0.59 Q Fl 1 + 1.04 × 10-4 p-1.791 w

()

p)

0.5(1-0.91dh/p)

Ah Aa

(52)

(53)

The liquid flow path length and the weir length are proportional to the column diameter:

z ) 0.6D h

(54)

w ) 0.8D h

(55)

3.2.2.2. Undesirable Effects. To have a representative evaluation of the hydraulic behavior of the column in both rectification and stripping sections, the following equations and inequalities are calculated on the first nonequilibrium stage (subscript 1) and on the last one (subscript nt). Entrainment flooding occurs when the upward vapor velocity is high enough to suspend a liquid droplet. The minimum column diameter recommended to avoid this problem is calculated using the correlations of Kister and Hass:20

D h gD h min D h min )

(

4At,min π

(56)

)

(62)

(63)

Downflow flooding occurs when columns flood because of their inability to handle large quantities of liquid. To avoid this problem, one must ensure that excessive backup does not occur. Then, a pressure-balance equation is written:21

∆P + ∆Pq g(Fl - Fv)

(64)

In the previous inequality, the right-hand side is equal to the height of liquid in the downcomer (hdc in Figure 1). The clear liquid height hl is the sum of the weir height hw and the weir height crest how (without hydraulic gradient across plate height):

hl ) hw + how

()

(61)

The clear liquid height at the froth-to-spray transition is given by

where β is determined as

w β ) 2 arcsin D h

0.5

(51)

D h 2 π β - sin(β) Ada ) 0.5 2 180

( )(

(59)

v uN,max

(hw + s) g hl +

The active area is determined as

(58)

Qv

Aa,min )

(49)

3.2.2. Tray Hydrodynamic Equations. Geometrical parameters of the trays are submitted to a set of constraints in order to (i) ensure feasibility of the design from a geometrical point of view and (ii) avoid three undesirable effects on the tower: entrainment flooding, down-flow flooding, and weepingdumping. 3.2.2.1. Geometrical Relations. The holes are located in the corners of equilateral triangles. The ratio between the pitch of perforation (distance between the centers, ph) and the hole diameter (dh) varies from 2.5 to 5. As was said before, dh and ph are, in our case, fixed parameters. For such an arrangement, the hole area is given by

Aa,min 0.7

At,min )

1/2

(57)

(65)

where the weir height crest how may be calculated from the Francis weir equation for a segmental weir:21

how ) 0.6

() Ql w

2/3

(66)

Pressure drop across downcomer is given by

[( ) ( ) ]

∆Pq ) 1.62Fl

Ql 2 Qv + Sq Sq′

2

(67)

The pressure drop for the gas, ∆P, is the sum of the effects of the pressure drop across the tray ∆Ps, the pressure drop due to the superficial tension ∆Pσ, and that caused by the presence of liquid ∆PL.

∆P ) ∆Ps + ∆Pσ + ∆PL

(68)

Pressure drop across the tray

∆Ps )

()

1 Qv 2 FG(1 - p2) 2 A 2K0 h

(69)


plifications,4 the objective function is as follows (assuming that D h 1.066 = D h ):

where

(

K0 ) 880.6 - 67.7

()

()

( ))

dh dh 2 dh + 7.32 - 0.338 ep ep ep

3

× 10-3 (70)

4σl dh

l (fi,j C$,i) - Ln +1V$,n +1 + CRQB + ∑ ∑ j)1 i)1

( (

t

nt

h 1.55 CconQC + AFCCAT + CTD

Pressure drop due to the superficial tension

∆Pσ )

{

n t nc

Cost ) min Co +

CSHD h s0 +

(71)

nt

t

(

∑ s + 1.27 J)1

D h2

)

+

)) }

W s + 1.27 J)1 D h2

∑

W

0.802

(79)

The exchanger investment costs are supposed to be linear functions of the heat duty:

Pressure drop due to the liquid

∆PL ) faFlghl

(72)

CB ) Cr1 + Cr2QB

(80)

CC ) Cc1 + Cc2QC

(81)

with By lumping constants,

fa ) 0.981e(-0.411Fa) if Fa e 1

(73)

Co ) AF(Cr1 + Cc1)

(82)

fa ) 0.718e(-0.079Fa) if Fa > 1

(74)

CR ) CH + AFCr2

(83)

Ccon ) CW + AFCc2

(84)

The Fa factor for flow through holes is calculated as follows:

Qv Fa ) xFv Aa

(75)

The last undesirable effect, weeping-dumping, occurs if the vapor velocity through the holes is too small. Then, the liquid will drain through them. To avoid this problem, a minimum gas velocity through the holes is calculated: v uvh g uh,min

(76)

CT and CSH lump the constants for the installed costs of the trays and of the shell. An interesting article by Ciric and Gu2 gives the values for the constants used in eq 79 (considering the Marshall and Swift index of 1994):

Co: $10,000 year-1 CR: $146.8 kW-1‚year-1 Ccon: $24.5 kW-1‚year-1 CT: $15.7 year-1

The minimum gas velocity (weep point) through the holes below which excessive weeping is likely (Locket equation):21 v ) uh,min

0.68 ( 0.12

xFv/(Flgfahl)

(77)

3.3. Objective Function. The objective function to be minimized is the sum of two terms: (i) the annualized capital cost (ACC) and (ii) the annual operating costs (AOP). The ACC includes the installed cost of: shell, trays, the reboiler, the condenser, and catalyst, with a 5 year payback. The AOP considers the consumption of raw material, steam, and cooling water. The value of the different products is also included in the AOP. n t nc

l (fi,j C$,i) - Ln +1V$,n +1 + CHQB + ∑ ∑ j)1 i)1

Cost ) min {

t

t

CWQC + AF(CCS + CCi + CB + CC + CCAT)} (78)

CSH: $222 year-1 3.4. Conclusion. The MINLP optimization problem for optimal design of catalytic distillation columns has been stated. The geometrical parameters of the trays have direct incidence on the mass and heat transfer that, in turn, have a strong influence on the production and, consequently, on the cost function. For a system of nc components, and a column of nt sections (reaction + separation), there are (18 + 5nt + 4nc + 5ntnc + 2nr2) variables to be optimized, including (4 + nr2) integer variables. There are (10 + 5nt + 4nc + 5ntnc) equality constraints: eqs 50, 51, 54, 55, and the model equations. There also are 4 inequality constraints: relations 49, 56, 64, and 76. Please note that the number of variables and constraints is a 3 function of the optimization variables: nt )∑j)1 ns,j. Furthermore, this MINLP problem is nonconvex. Thus, particular attention has to be paid to the selection of the solution strategy and the optimization algorithms. 4. Solution Strategy

The different expressions to calculate the column investment costs (shell, trays, condenser, and reboiler) are from ref 22. Substituting these expressions, and after some algebraic sim-

4.1. Introduction. For such a complex, nonconvex MINLP model, solution strategy is a critical point. Optimization techniques can be classified into two main groups: deterministic

1380


Figure 4. Solution strategy.

and stochastic algorithms.23 Deterministic algorithms include local or global optimization techniques. In a first attempt to solve the formulated problem, we have chosen a two stage strategy that combines both stochastic and deterministic algorithms. The main reasons are as follows: Local deterministic MINLP optimization techniques include branch and bound,24 generalized bender’s decomposition (GBD),25 outer approximation,26 and extended cutting plane.27 With such algorithms, depending on the initial point, the probability to reach the global optimum is quite weak. Global deterministic MINLP methods28,29 are very promising techniques: Rbranch and bounds (RBB),30 branch and reduce algorithm,23 symbolic manipulation algorithm,31 interval subdivision,32 etc. Deterministic MINLP algorithms for global optimization take advantage of the mathematical structure of the problem but are generally restricted to particular classes of models. Furthermore, only medium scale problems can be considered. Stochastic methods can increase the probability to achieve the global optimum: random search (RS),33 simulated annealing (SA),4 and genetic algorithm (GA).34 For large scale MINLP problems, pure stochastic techniques result in a prohibitive computational load. For all these reasons, we have implemented a two stage solution strategy that combines a stochastic and a local deterministic algorithm (Figure 4). A nonlinear programming subproblem is formulated and solved using a deterministic algorithm: successive quadratic programming (SQP). Most of the continuous variables are considered at this level. Integer and remaining continuous variables are optimized at the level of the MINLP master problem. This medium scale MINLP problem is solved using a stochastic algorithm: simulated annealing (SA). In the literature, there are works with this combination of algorithms (a stochastic and a local deterministic algorithm). An example of this kind of optimization structure is presented by Athier et al.35 They propose an approach for the heat exchanger networks (HENs) problem with simulated annealing at an upper level and a lower level within a nonlinear programming (NLP/SQP) frame to optimize the operation variables. Another example in the literature is proposed by Montastruc et al.36 To solve a solid-liquid equilibrium optimization problem, they use an optimization technique based on

the successive use of a genetic algorithm (GA) and of a classical sequential quadratic programming (SQP) method. 4.2. NLP Subproblem. At this level, the objective is to minimize the cost function. Optimization variables are: model variables (6 + 4nc + 5nt + 5ncnt), operational variables (2), and the catalyst load per stage (nr2) which are part of the design variables. This minimization is submitted to the model constraints (6 + 4nc + 5nt + 5ncnt equality equations) and to the product specification (1 inequality equation). Here, the other variables, especially integer variables, are fixed parameters. This problem results in a nonlinear programming model. Basically, there are three different methods to solve the NLP problems: generalized reduce gradient, sequential quadratic programming, and adaptations of successive linear programming (cutting plane method, etc). Due to the size of the NLP subproblem (6 + 4nc + 5nt + 5ncnt equalities) and to the degree of freedom (2 + nr2), we have implemented a variation of the SQP method: the reduced hessian successive quadratic programming (rSQP). This algorithm has been shown to be successful in solving problems with large numbers of equations but relatively few degrees of freedom.37 More precisely, the solver QPKWIK, based upon the dual algorithm of Goldfard and Idnani and developed by Schmid, Lorenz, and Biegler at Carnegie Mellon University,38 is used. We have experienced that initialization and bounding of the variables is an essential part for a successful optimization. For that reason, at the initialization step, the catalytic distillation model is solved using the Newton-Raphson method. Then, at the initial point of the optimization procedure, the model equations are satisfied. Among all the tested initialization strategies, this one has been proven to be the best one. Furthermore, variables are bounded with physical criteria. 4.3. MINLP Master Problem. At this level, the design variables are optimized (except the catalyst loads): 4 + nr2 integer variables (number of reactive stages in the second section nr2, number of stages per section ns,j)1-3, location of the reactive stages in the second section cl2,j)1-nr2); 6 continuous variables (weir height hw, column diameter D h , hole area Ah, active area Aa, liquid flow path length z, weir length w). Optimization constraints include the geometrical equations and the hydrodynamic relations. The objective function is the penalized cost function.


To solve the MINLP master problem, we have selected the simulated annealing algorithm. The name of the method is derived from the simulation of thermal annealing of critically heated solids. A slow and controlled cooling of a heated solid ensures proper solidification with a highly ordered crystalline state that corresponds to the lowest internal energy. Rapid cooling, on the other hand, causes defects inside the material. It is recognized that SA is a good technique which might help in the location of a near-optimum solution. The role of the master problem is to randomly propose a set of design variables, allowing wrong way movements to increase the probability of convergence toward the global optimum or, at least, the location of a near-optimum solution. All changes that favorably affect the objective function are accepted, and all changes that unfavorably affect the objective function are accepted with some finite probability. This probability depends on the parameter called “temperature”, which decreases periodically during the optimization. The accepting or rejecting is determined by the Metropolis criterion with a pseudorandom number sequence function. This process continues until some predetermined stopping criterion has been met: no more improvement in the objective function or maximum number of move attempts (sufficiently low temperature). 4.4. Conclusion. The proposed strategy is a compromise between stochastic and deterministic algorithms. We do not claim that a global optimum will be reached. Such strategy can just improve the probability to achieve a global optimum in a reasonable computation time. One must keep in mind that, at the NLP subproblem level, the SQP may be trapped in a local minimum. Care has been shown to the initialization step in order to improve the robustness of the NLP subproblem solution. One should note that analytical derivatives are used at the SQP level in order to reduce the computation load. In this work, continuous design variables are considered in the MINLP master problem since analytical derivatives are not easily calculable: simulated annealing is a direct optimization technique (derivatives are not required).

Table 6. Thermodynamic, Physical, and Mass Transfer Models equilibrium ratio equation of state liquid molar volume binary mass transfer coefficients enthalpy liquid activity coefficients vapor pressure heat transfer coefficients liquid binary diffusion coefficients (liquid mixtures) liquid binary diffusion coefficients (dilute liquids) vapor binary diffusion coefficients reference for enthalpy calculation

(CH2)CdCH2 + C2H5OH S (CH3)3COC2H5

(85)

The main side reactions are the dimerization of isobutylene to form diisobutylene (DIB) and the hydratation of isobutylene to form isobutyl alcohol, if there is water in the system. The

Fuller Href ) 0 at 398.15 K, ideal gas, 101 325 Pa

(CH3)2CdCH2 + (CH3)CdCH2 S [(CH3)2CdCH2]2 (CH3)2CdCH2 + H2O S (CH3)3COH

(86) (87)

The preferred catalyst A15 is a macroreticular sulfonic resin manufactured by Rohm and Hass. It is a copolymer of styrene and divinilbenzene sulfonated which has a surface area of 45 m2‚g-1 and an ion exchange capacity of 4.8 mequiv of H+‚g-1 dry.43 The catalyst is susceptible to both deactivation and poisoning. Deactivation occurs over a much longer period, and it is accelerated by thermal degradation (hot spots). The kinetics used corresponds to an expression based on the Langmuir-Hinsherwood-Hougen-Watson (LHHW) model, which involves two adsorbed ethanol sites for the reaction in the rate-limited step. The liquid phase is strongly nonideal and has to be taken into consideration by using activities (γx) instead of molar fractions (x):44

RETBE,j )

The example given below is the cost optimization of ETBE production, as a main product. A minimum 95% molar purity of ETBE is required. The significance of this reaction resides in the growing worldwide consumption of the ETBE. The use of the latter is in line with an endeavor to improve the air quality in certain sensitive areas, by using oxygenated compounds as a substitute to other substances which cause noxious emissions. The example is first described. In the second step, a sensitivity analysis is performed. Then, optimization results are discussed. 5.1. Reaction. This reaction has been treated by several researchers: in 1997, Sneesby39 modeled a catalytic distillation column with this reaction, whereas Thiel40 worked on systems of discontinuous distillation by means of residual curves. AlArfaj and Luyben41 carried out a study on distillation control, and Jhon and Lee42 investigated the dynamic simulation of a catalytic column with the reaction of ETBE production. The formation of ETBE from ethanol and isobutene is an acid catalyzed (Amberlyst 15, A15), reversible etherification, moderately exothermic reaction, in the liquid phase.

Wilke-Chang

first side reaction can be minimized by using excess ethanol, and the latter is neglected in this example because it occurs only in the presence of water.

(

κrate,j(γx)2EtOH,j (γx)IB,j 5. ETBE Illustrative Example

Gamma-Phi Soave-Redlich-Kwong Rackett AIChE ideal UNIFAC Antoine Chilton-Colburn Vignes

(γx)ETBE,j κETBE,j(γx)EtOH,j

)

(1 + κA,j (γx)EtOH,j)3

(88)

Ln(κETBE,j) ) 10.387 + 4060.59/Tlj - 2.89055 ln(Tlj) 0.01915144(Tlj) + 5.28586 × 10-5(Tlj)2 - 5.32977 × 10-8(Tlj)3 (89) Ln(κA,j) ) -1.0707 +

1323.1 Tlj 3

(90) l

κrate,j ) 2.0606 × 1012 e(-60.4×10 /RTj)

(91)

5.2. Thermodynamic Properties. The thermodynamic equilibrium at the vapor-liquid interface is described as follows: i yi,j

sat γi,jPsat i φi Ki,j ) i ) φi,jP xi,j

(92)

The modified UNIFAC (Dortmund) model for the prediction of the activity coefficients is used. The Soave-Redlich-Kwong (SRK) equation is used to calculate the vapor fugacity coefficients. Thermodynamic properties are shown in Table 6. 5.3. Economical Parameters. The raw material costs are the following: $15 kmol-1 for ethanol; $8.25 kmol-1 for butenes

1382


Table 7. Initial Configurations variable

IF

INF

column diameter (m) reflux ratio reboiler duty (W) weir height (m) weir length (m) active area (m2) flow path length (m) hole area (m) number of reactive sections reactive section location number of trays feed tray positions ETBE purity annualized cost ($‚year-1)

10-2

12 × 10-2 10 10 000 5.0 × 10-2 9.0 × 10-2 5.0 × 10-3 8.0 × 10-2 7.0 × 10-4 1 15 30 10/2 29%

8.8 × 7.93 103 59 1.9 × 10-2 7.0 × 10-2 4.3 × 10-3 5.3 × 10-2 5.5 × 10-4 6 20/24/25/29/30/36 48 19/36 95% 10,642

Figure 5. Effect of the reflux ratio on the bottom ETBE purity and objective function.

Table 8. Column Feed feed

temp (K)

flow rate (mol‚s-1)

XEtOH (molar)

XIB (molar)

XETBE (molar)

Xn-B (molar)

1 2

400 320

0.02853 0.091135

1.0 0.0

0.0 0.3

0.0 0.0

0.0 0.7

(iso + normal); $7.7 kg-1 for catalyst.41 The ETBE value is $25.3 kmol-1. 5.4. Initial Configurations. To check the global convergence of the proposed solution strategy, we have tested several initial points. In this article, we present the main results considering two initial configurations (Table 7). We first consider a feasible (IF) configuration: all optimization constraints (model equations, purity, and hydrodynamic constraints) are satisfied. This configuration has been obtained by trial and error, using simple simulations, without any optimization algorithm. This is a long and tedious work. A second initial point is considered. This configuration is said infeasible (INF) since optimization constraints are not satisfied. In the feasible configuration, on each of the six reactive sections, 400 g of catalyst is introduced. The total amount of catalyst in the reactive sections is 2400 g. For the other initialization, only 150 g of catalyst is used. The following parameters are fixed: operational pressure (950 000 Pa), hole pitch (0.008 m), hole diameter (0.003 m), tray spacing (0.15 m), number of feeds (2), feed compositions, and feed temperatures. Table 8 gives the summary of the feed column specifications. 5.5. Sensitivity Analysis. Let us consider the feasible initial configuration. We now perform a sensitivity analysis in order to better understand the influence of the different optimization variables (both operational and design variables) on the objective function and on the ETBE purity. Variables are varied separately. The limitations introduced by the constraints (mainly hydraulic constraints) are discussed. Figure 5 shows the ETBE purity in the bottom product and the objective function at different reflux ratios. The ETBE purity is non-monotonic with respect to the reflux, and it shows the presence of an optimum ratio: when the reflux ratio is too high, the recirculation of unreacted heavy reactant does not occur in the column; when the reflux ratio is too small, unreacted light reactant is not recycled. In both cases, the conversion rate decreases. Thus, ETBE purity presents a maximum. The annualized cost is a decreasing function of the reflux ratio: annualized cost is a function of the bottom total flow rate. If the reflux ratio increases, the bottom flow increases. Since incomes are increasing, the annualized cost decreases. A symmetric situation occurs with the bottom product ETBE purity and the objective function at different reboiler duties (Figure 6). Both operational variables (reflux ratio and reboiler

Figure 6. Effect of the reboiler duty on the bottom ETBE purity and objective function.

Figure 7. Effect of the diameter on the bottom ETBE purity and objective function.

duty) have maximum peaks close to 5 and 7000, respectively. They are strongly limited by the bottom ETBE purity (95% minimum). In conclusion, Figures 5 and 6 show how sensitive the ETBE purity and the objective function are to the operational variables. According to the objective function formulation (eq 79), the annualized cost increases with the column diameter (Figure 7). When the diameter increases, the vapor velocity decreases. Thus, the duration of the contact between the vapor and the liquid is longer and, consequently, the mass transfer and ETBE purity is better. When the diameter is too small, the following constraints are violated: first, entrainment flooding and, next, ETBE purity and downflow flooding. For higher values of the diameter, weeping-dumping occurs. ETBE purity and annualized cost also increase with the weir height (Figure 8): mass transfer is more efficient for high values of the weir height (eqs 20 and 23). The minimum weir height is limited by the downflow flooding problems (eq 64), and the maximum is limited by the downflow flooding and also the weeping-dumping problems (eq 77): the minimum gas velocity grows with the weir height. The other analyses show that there is no guarantee that the bottom ETBE purity can be achieved in the column unless there


Figure 10. Exponential temperature schedule for the simulated annealing. Figure 8. Effect of the weir height on the bottom ETBE purity and objective function.

Figure 9. Effect of the number of rectification trays on the bottom ETBE purity and objective function.

are a sufficient number of reactive sections, with a sufficient amount of catalyst. A minimum number of stripping and rectifying stages is also required. However, Figure 9 shows that an increment in the rectification section above the optimum is not necessarily beneficial for the reactive section conditions. Too much rectification separation can result in isobutylene loss via the distillate, and the overall effect is usually to decrease the isobutylene conversion rate and then the corresponding product purity. From a quantitative point of view, paying attention to the annualized cost scale, one should note that the objective function is not very sensitive to this optimization variable. 6. Optimization Results and Discussions. Both optimizations were simultaneously solved on a workstation with two 1.50 GHz CPUs and 2096 MB RAM. The average CPU time per iteration was approximately 360 s, and the maximum number of NLP evaluations (number of move

attempts) was fixed at 10 000. The optimization duration is between 12.5 and 16 days. In this example, we use an exponential temperature schedule determined by the parameter R, equal to 0.8. The initial temperature is equal to 10 000, evaluated according to the method proposed by Pibouleau et al.45 with an initial acceptation rate between 0.94 and 0.96. The number of moves per temperature level (the process of changing the system state) is equal to 100 (Figure 10). The main results concerning the geometrical characteristics for both optimizations are presented in Table 9. The final configurations are very close. Comparing the two optimal objective functions, one should note that the optimal solution is slightly better in the IF case than in the INF case. The difference is 1.13% between the objective functions. It should be pointed out, however, that the agreement between the two solutions is good for almost all the variables. There are only two variables with significant differences: the number of rectification and stripping stages. However, the total number of trays is different only by three trays. The reason that the results from the optimization are not significantly dependent on the number of trays is directly related to the cost function. On the basis of the sensitivity analyses discussed earlier, it can be seen that the objective function is not very sensitive to the number of trays in both sections. However, this does not mean that any total number of trays will be accepted as solutions because it has to satisfy the purity constraint (95% minimum, ETBE). The final solutions have 400 g of catalyst on each reactive section. This is the maximum quantity possible with the mechanical configuration of the column: the catalyst weight on each reactive section was limited in such a way that the weir height cannot exceed 10 cm to prevent excessive pressure drop. The saturation of theses variables indicates that the position of the reactive section is adapted to the reaction condition: a bad

Table 9. Final Configurations variable

IF

INF

search range

column diameter (m) reflux ratio reboiler duty (W) weir height (m) weir length (m) active area (m2) flow path length (m) hole area (m) total number of trays number of trays in the rectification section: ns1 number of trays in the stripping section: ns3 number of trays in the reactive section: ns2 number of reactive stages in section 2: nr2 reactive stage location: cl2

9.3 × 10-2 5.07 6816 3.2 × 10-2 7.4 × 10-2 4.9 × 10-3 5.6 × 10-2 6.2 × 10-4 43 11 8 24 13 12/13/14/16/19/21/22/23/ 27/28/30/33/35 95.00% 98.16% 9,407

9.5 × 10-2 5.15 6974 2.3 × 10-2 7.6 × 10-2 5.1 × 10-3 5.7 × 10-2 6.5 × 10-4 46 7 10 24 12 7/10/11/14/17/20/22/ 24/25/27/29/30 95.00% 97.85% 9,519

(8.0-12) × 10-2 2-15 2000-15000 (0.8-10) × 10-2 w ) 0.8D h Aa ) At - 2Ada z ) 0.6D h Ah ) 0.907Aa[dh/ ph]2 30-55 5-30 5-30 4-25 1-13 between feed 1 and feed 2

ETBE purity isobutene conversion annualized cost($‚year-1)

min (95%)

1384


increment, the optimization increases the column diameter. This column diameter is limited by the weeping problem and by the increment of the cost function. Analyzing more in-depth the best optimal solution, one could note that the limiting constraints are the ETBE bottom purity and the downflow flooding: these constraints are saturated. The minimum diameter and maximum weir height recommended by the hydraulic equations are respectively 7.1 and 0.33 cm. A postoptimal sensitivity analysis is also performed (Table 10). Small changes in the optimization variables do not improve the solution: the objective function is increased or a constraint is violated. We can conclude that the proposed best solution is a local minimum. Of course, it is impossible to determine if it is a global optimum. 7. Comparison between NEQ and EQ Optimal Solutions

Figure 11. (a) Objective function trajectory and (b) initial objective function trajectory.

location makes the optimizer (at the NLP subproblem level) move the catalyst load toward the minimum allowed value, which is fixed at 100 g for this example. Figure 11a shows, for the two MINLP runs, the trajectories for the best objective function values. For the IF case, the operational cost is decreased by 11.60%. In the case of infeasible initialization (INF), at the beginning, the objective function is strongly penalized (1 × 1018) because the initial configuration does not satisfy the purity constraint and the hydrodynamic restrictions (see Figure 11b). The operational cost is decreased by 13.95% between the first feasible solution ($11,063 year-1) and the final solution ($9,519 year-1). Comparing the initial configuration (feasible initialization, Table 7) and the optimal solution (Table 9), the tendency of the optimization is to simultaneously decrease the values of the reflux ratio and the reboiler duty, to reduce the operating costs but respecting the product specification. The operation is highly sensitive to these two variables. Sneesby et al.46 have shown that only a few combinations of the reflux ratio and the reboiler duty will produce the desired results of a high isobutylene conversion and a high ether purity at the same time. For these new conditions, the weir height is increased to compensate the mass transfer loss by the reduction of the reflux rate. This increment is limited by a hydraulic restriction (downflow flooding) and by a small increment of the cost function. To avoid downflow flooding problems due to the weir height

In this section, we propose a brief comparison of the optimizations based on the equilibrium model and the nonequilibrium model. For sake of comparison, both optimizations are performed using the same solution strategy and the same algorithms. The objective function and the various numerical parameters (such as the maximal number of iterations, initializations, penalizations, etc.) remain the same. The equilibrium model used in this study follows the conventional approach for a tray column which assumes physical and mechanical equilibrium. In a first attempt, Murphree tray efficiencies were assumed to be 1.0. Table 11 and Figure 12 show the differences between both optimizations. As an expected result, the cost is smaller and the isobutene conversion is greater using the equilibrium model. One can note that, in this case, optimal structures (number of separation sections, number of catalytic sections) are quite similar. However, there are significant differences concerning operating and design parameters: reflux, reboiler duty, and column diameter. Again, a detailed comparison of the two approaches is beyond the objectives of this contribution. From a more general point of view, the question to know if rate based models are more reliable than equilibrium models is a quasi-philosophical discussion. From our point of view, this discussion is still open, but let us consider the following points. The mentioned differences can be a consequence of the efficiency value. Then, the problem is how to compute the tray efficiencies: (i) constant values (1.0, 0.7, ...); (ii) empirical correlation (Mac Farland, etc.); (iii) mechanistic models (AIChE, etc.). It is very difficult to determine, a priori, the correct constant values of the tray efficiencies. The use of empirical correlations is said to be unreliable for general multicomponent, reactive systems.47 Various authors have proposed to use the AIChE model in order

Table 10. Postoptimal Sensitivity Analysis ETBE purity

annualized cost ($‚year-1)

5.00 5.07 5.10 6800 6816 6850 9.0 × 10-2

0.967 10 0.950 00 0.944 40 0.946 20 0.950 00 0.957 00 0.949 00

9,861 9,407 9,280 9,316 9,407 9,584 9,302

9.3 × 10-2 9.5 × 10-2 3.0 × 10-2 3.2 × 10-2 3.5 × 10-2

0.950 00 0.950 30 0.949 40 0.950 00 0.950 00

9,407 9,413 9,400 9,407 9,414

variable reflux ratio reboiler duty (W)

column diameter (m) weir height (m)

violated constraint

ETBE purity ETBE purity downflow flooding, ETBE purity ETBE purity downflow flooding


Figure 12. Final configurations: NEQ and EQ optimal solutions. Table 11. Comparison between the Equilibrium and the Nonequilibrium Optimal Solution NEQ model final config

variable column diameter (m) reflux ratio reboiler duty (W) weir height (m) weir length (m) active area (m2) flow path length (m) total number of trays number of trays in the rectification section: ns1 number of trays in the stripping section: ns3 number of trays in the reactive section: ns2 number of reactive stages in section 2: nr2 reactive stage location: cl2 ETBE purity isobutene conversion annualized cost ($‚year-1)

EQ model final config

9.3 × 10-2 5.07 6816 3.2 × 10-2 7.4 × 10-2 4.9 × 10-3 5.6 × 10-2 43 11

8.4 × 10-2 4.5 6473

8

7

24

24

13

11

12/13/14/16/19/21/ 22/2327/28/30/33/35 95.00% 98.16% 9,407

16/18/20/21/22/ 24/26/31/34/36 95.00% 99.77% 8,970

of a catalytic distillation column, based on a steady state nonequilibrium (rate based) model, has been investigated. It should be recalled that the aim of a nonequilibrium model is to describe a real column and, therefore, it requires information about its design and internal configuration. This information is needed to calculate mass transfer coefficients, avoiding the use of Murphree tray efficiencies which are very difficult to calculate in the case of a multicomponent separation. The optimized mechanical tray configuration satisfies the hydraulic constraints in order to guarantee the correct design of the column. The column model is first described, next the optimization problem is formulated, and finally the solution strategy and the selected algorithms are discussed. We consider the ETBE catalytic distillation as an illustrative example. This is a quite challenging example because of nonideality and reaction exothermicity. Because of the complex interactions between the different optimization variables and their influence on the objective function and the limiting phenomena (optimization constraints), it is quite impossible to predict the great trends of the design. Nevertheless, we propose a sensitivity analysis in order to better understand the global influence of the different variables. Next, optimization results are discussed in connection with the sensitivity analysis. These results are also analyzed performing a postoptimal sensitivity analysis. Future developments will include (i) integration of the proposed optimization techniques and the feasibility analysis,10 as a preliminary step, to improve the initialization and variable search range definition, (ii) optimization of the operational pressure, (iii) a more extensive comparison between equilibrium and nonequilibrium models, including the study of the influence of tray efficiency, and (iv) generalization of the proposed approach to a packed column. Nomenclature

44 13

to compute the efficiency value.16 Krishnamurthy13 has reported that it is more complicated to use the AIChE relation to compute efficiencies than to perform mass balances on each phase, as done in the rate based model. Moreover, the use of tray efficiency does not allow one to describe the heat transfer phenomenon.48 For all these reasons, we consider, at the moment, that the use of nonequilibrium model should be preferred for MINLP optimization since (i) such optimization is possible as shown in this contribution and (ii) the rate based model seems to be “physically more realistic”.15 8. Conclusions By using a combination of simulated annealing (SA) and sequential quadratic programming (SQP), the optimal design

Aa ) active area (m2) Ada ) area under the downcomer (m2) AF ) annualizing factor (year-1) Ah ) hole area (m2) At ) total column cross-sectional area (m2) a ) interfacial area (m2) CB ) installed cost of reboiler ($) CC ) installed cost of condenser ($) CCAT ) cost of the catalyst ($) CCi ) installed cost of trays and reactive sections ($) CCS ) installed cost of shell ($) CH ) cost of heating ($‚W-1‚year-1) Cp ) molar heat capacity (J‚mol-1‚K) Csbf ) capacity parameter (m‚s-1) CW ) cost of cooling ($‚W-1‚year-1) C$ ) raw material cost ($‚mol-1) cl2 ) catalytic stage location in section 2 (-) D ) Fick diffusivity in the vapor phase (m2‚s) D0 ) distillate (mol‚s-1) D h ) column diameter (m) dh ) hole diameter (m2) E ) referring of the product, molar fraction of the principal product (-) ep ) tray thickness (m) F ) total feed flow rate of liquid (mol‚s-1) Fa ) F-factor for flow through holes (m‚s-1‚(kg‚m-3)0.5) Fs ) superficial factor (m) f ) feed flow rate of a component (mol‚s-1)

1386


fa ) aeration factor (-) g ) acceleration of gravity (m‚s2) H ) molar enthalpy (J‚mol-1) H h ) partial molar enthalpy of a component (J‚mol-1) h ) heat transfer coefficient (W‚m-2‚K) hcl ) clear liquid height at the froth-to-spray transition (m) hdc ) height in downcomer (m) hl ) clear liquid height (m) hw ) weir height (m) how ) weir height crest (m) K ) equilibrium ratio (-) [k] ) matrix of mass transfer coefficient in binary mixture (mol‚m-2‚s-1) [kh] ) matrix of multicomponent low flux mass transfer coefficients (mol‚m-2‚s-1) K0 ) Economopoulos coefficient (-) k ) mass transfer coefficient in binary mixture (mol‚m-2‚s-1) κA ) adsorption equilibrium constant (-) kETBE ) equilibrium constant, kinetic (-) κrate ) reaction rate constant (mol‚kg-1‚s-1) kh ) multicomponent low flux mass transfer coefficients (mol‚m-2‚s-1) L ) total liquid flow (mol‚s-1) Le ) Lewis number (-) M ) molar mass (g‚mol-1) mcat ) catalytic mass (kg) N ) mass transfer rate (mol‚s-1) N ) molar flux of a component (mol‚m-2‚s-1) nc ) total number of components nf ) total number of feed trays nr2 ) number of reaction stages in section 2 ns ) number of stages nt ) total number of separation sections and reaction sections P ) operational pressure (Pa) ph ) pitch of perforations (m) Q ) volumetric flow rates (m3‚s-1) QB ) condenser duty (W) QC ) condenser duty (W) R ) rate of reaction (mol‚kg-cat-1‚s-1) [R] ) matrix function of inverted binary mass transfer coefficients (m2‚s‚mol-1) Sc ) Schmidt number (-) Sq ) vertical area between tray and downcomer apron (m2) Sq′ ) horizontal area between downcomer apron and inlet weir (m2) s ) tray spacing (m) s0 ) extra column height (m) T ) temperature (K) TR ) reflux ratio (-) t ) residence time (s) uh ) superficial velocity, based on the hole area (m‚s-1) uN ) superficial velocity, based on the active area (m‚s-1) V ) total vapor flow rate (mol‚s-1) V$ ) product value ($‚mol-1) WK ) holdup volume on tray (m3) w ) weir length (m) x ) molar fraction liquid phase (-) y ) molar fraction vapor phase (-) z ) liquid flow path length (m) Greek Letters β ) angle formed by the semicircle of the downcomer γ ) activity coefficient (-) [Γ] ) matrix of thermodynamic factor (-)

δ ) Kronecker delta (-) ) energy transfer due to the mass transfer (J‚s-1) p ) fraction hole area based on active area (-) µ ) viscosity (Pa‚s-1) F ) density (kg‚m-3) σ ) surface tension (N‚m-1) sat φi,j ) fugacity coefficient of saturated vapor, for pure species i (-) φi,j ) fugacity coefficient of i at the vapor phase (-) φ ) association factor for the solvent. (-) ∆Hr ) heat of reaction (J‚mol-1) ∆P ) pressure drop (Pa) ∆PL ) pressure drop due to the liquid (Pa) ∆Ps ) pressure drop across the tray (Pa) ∆Pq ) pressure drop across downcomer (Pa) ∆Pσ ) pressure drop due to the superficial tension (Pa) Subscripts av ) average B ) referring to reboiler E ) principal product i, j, k ) component indices. j ) referring to stage or section number max ) maximum min ) minimum nc ) total number of components nt ) total number of stages or sections nt + 1 ) referring to reboiler P ) constant pressure T ) constant temperature t ) referring to total mixture ∑ ) constant all the other species 0 ) referring to condenser Superscripts i ) referring to interface LF ) referring to liquid feed l ) referring to liquid phase nc ) total number of components nt ) total number of stages or sections sat ) referring to saturation VF ) referring to vapor feed v ) referring to vapor phase o ) referring to infinitely low concentration Miscellaneous-Abbreviations ACC ) annualized capital cost ($‚year-1) AOP ) annual operating cost ($‚year-1) A15 ) Amberlyst 15 BB ) branch and bound CD ) catalytic distillation CDC ) catalytic distillation column cl ) catalyst locations (-) EQ ) equilibrium model ETBE ) ethyl tert-butyl ether ftl ) feed tray locations (-) GAMS ) general algebraic modeling systems GBD ) generalized bender decomposition GDP ) general disjunctive programming GRG ) reduced gradient generalized IF ) feasible initialization INF ) infeasible initialization. IP ) integer master problem LaTEP ) Laboratoire de Thermique, Energe´tique et Proce´de´s LGC ) Laboratoire de Genie Chimique


LHHW ) Langmuir-Hinshelwood-Hougen-Watson model MIDO ) mixed integer dynamic optimization MINLP ) mixed integer nonlinear programming MTBE ) methyl tert-butyl ether NEQ ) nonequilibrium model nt ) total number of stages or sections OA/ER/AP ) outer-approximation/equality-relaxation/augmentedpenalty OCFE ) orthogonal collocation on finite elements RD ) reactive distillation rSQP ) reduced Hessian successive quadratic programming SA/N ) simulated annealing/Newton SQP ) sequential quadratic programming SRK ) Soave-Redlich-Kwong Literature Cited (1) Almeida-Rivera, C. P.; Swinkels, P. L. J.; Grievink, J. Designing Reactive Distillation Processes: Present and Future. Comput. Chem. Eng. 2004, 28, 1997. (2) Ciric, A. R.; Gu, D. Synthesis of Non Equilibrium Reactive Distillation Processes by MINLP Optimization. AIChE J. 1994, 40 (9), 1479. (3) Poth, N.; Frey, Th.; Stichlmair, J. MINLP Optimization of Kinetically Controlled Reactive Distillation Processes. 34th Eur. Symp. Work. Party Comput. Aided Process Eng.-11 2001, Supplementary Proceeding Volume, 79. (4) Cardoso, M. F.; Salcedo, R. L.; Feyo de Azevedo, S.; Barbosa, D. Optimization of Reactive Distillation Processes with Simulated Annealing. Chem. Eng. Sci. 2000, 55, 5059. (5) Poth, N.; Brusis, D.; Stichlmair, J. Rigorous Optimization of Reactive Distillation in GAMS with the Use of External Functions. Eur. Symp. Comput. Aided Process Eng.-13 2003, 869. (6) Sand, G.; Barkmann, S.; Engell, S.; Schembecker, G. Structuring of Reactive Distillation Columns for Non-Ideal Mixtures Using MINLP Techniques. Eur. Symp. Comput. Aided Process Eng.-14 2004, 493. (7) Frey, T.; Stichlmair, J. MINLP Optimization of Reactive Distillation Columns. Eur. Symp. Comput. Aided Process Eng.-10 2000, 115. (8) Jackson, J. K.; Grossmann, I. A Disjunctive Programming Approach for the Optimal Design of Reactive Distillation Columns. Comput. Chem. Eng. 2001, 25, 1661. (9) Rouzineau, D.; Prevost, M.; Meyer, M. Evaluation of Coupled Reactive Distillation Performances by Means of a Rigorous Simulation Procedure. Comput.-Aided Chem. Eng. 2001, 9, 267. (10) Thery, R. Analyse de Feasibilite´ Synthe`se et Conception de Proce´de´s de Distillation Reáctive. Ph.D. Thesis, INP., Toulouse, France, 2002. (11) Druart, F. Contribution a` la Distillation Catalytique: Mode´lisation et Etude Expe´rimentale. Ph.D. Thesis, Universite´ de Pau et des Pays de l’Adour, France, 2002. (12) Higler, A. P.; Taylor, R.; Krishna, R. Nonequilibrium Modelling of Reactive Distillation: Multiple Steady States in MTBE Synthesis. Chem. Eng. Sci. 1999, 54, 1389. (13) Krishnamurthy, R.; Taylor, R. A Non Equilibrium Stage Model of Multicomponent Separation Processes. Part I: Model Description and Method of Solution. AIChE J. 1985, 31 (3), 449. (14) Krishnamurthy, R.; Taylor, R. A Non Equilibrium Stage Model of Multicomponent Separation Processes. Part II: Comparison with Experiment. AIChE J. 1985, 31 (3), 456. (15) Krishnamurthy, R.; Taylor, R. A Non Equilibrium Stage Model of Multicomponent Separation Processes. Part III: The Influence of Unequal Component - Efficiencies in Process Design Problems. AIChE J. 1985, 31 (12), 1973 (16) Taylor, R.; Krishna, R. Multicomponent Mass Transfer; John Wiley and Sons: New York, 1993. (17) Bravo, J. L.; Fair, J. R. Generalized Correlation for Mass Transfer in Packed Distillation Columns. Ind. Eng. Chem. Process Des. DeV. 1982, 21, 162. (18) Onda, K.; Takehuchi, H.; Okumoto, Y. Mass Transfer Coefficients Between Gas and Liquid Phases in Packed Columns. J. Chem. Eng. Jpn. 1968, 1, 56. (19) Battfeld, M., Aguirre, P.; Grossmann, I. Alternative Representation and Formulation for the Economic Optimization of Multicomponent Distillation Columns. Comput. Chem. Eng. 2003, 27, 363.

(20) Perry, H.; Green, D. Perry’s Chemical Engineers’ Handbook; McGraw Hill: New York, 1997. (21) Cicile, J. C. Techniques de l’Inge´ nieur, No d’impression; 600409, J7.12, Paris, 1996. (22) Douglas, J. M. Conceptual Design of Chemical Processes; McGraw Hill: New York, 1998. (23) Ryoo, H. S.; Sahinidis, N. V. Global Optimization of Non Convex NLPs and MINLP with Applications in Process Design. Comput. Chem. Eng. 1995, 19, 551. (24) Beale, E. M. L. Integer Programming: the State of the Art in Numerical Analysis; Academic Press: London, 1977. (25) Geoffrion, A. M. Generalized Bender Decomposition. J. Optim. Theory Appl. 1972, 10, 237. (26) Duran, A. M.; Grossmann, I. E. A Mixed-Integer Nonlinear Programming Algorithm for Process Systems Synthesis. AIChE J. 1986, 32, 592. (27) Westerlund, T.; Petterson, F. An Extended Cutting Plane Method for Solving Convex MINLP Problems. Comput. Chem. Eng. 1995, 19, S131. (28) Floudas, C. A.; Akrotiriankis, I. G.; Caratzoulas, S.; Meyer, C. A.; Kallrath, J. Global Optimization in the 21st Century: Advances and Challenges. Eur. Symp. Comput. Aided Process Eng.-14 2004, 23. (29) Grossmann, I. E. Global Optimization in Engineering Design; Kluwer Academic Publishers: Dordrecht, 1996. (30) Adjiman, C. S.; Androulakis, I. P.; Maranas, C. D.; Floudas, C. A. A Global Optimization Methods RBB, for Process Design. Comput. Chem. Eng. 1996, 20, S419. (31) Smith, E. M. B.; Pantelides, C. C. Global Optimisation of General Process Models. In Global Optimization in Engineering Design; Grossmann, I. E., Ed.; Kluwer Academic Publishers: Dordrecht, 1996; p 355. (32) Byrne, R. P.; Bogle, I. D. L. Solving Nonconvex Process Optimisation Problems Using Interval Subdivision Algorithms. In Global Optimization in Engineering Design; Grossmann, I. E., Ed.; Kluwer Academic Publishers: Dordrecht, 1996; p 155. (33) Salcedo, R. L.; Goncalves, M. J.; Feyo de Azevedo, S. An Improved Random Search Algorithm for Non Linear Optimization. Comput. Chem. Eng. 1990, 14, 1111. (34) Kian, H. L.; Sorensen, E. Simultaneous Optimal Design and Operation of Multipurpose Bath Distillation Columns. Chem. Eng. Process. 2004, 43, 273. (35) Athier, G.; Floquet, P.; Pibouleau, L.; Domenech, S. Process Optimizations by Simulated Annealing and NLP Procedures. Application to Heat Exchange Network Synthesis. Comput. Chem. Eng. 1997, 21 (Suppl.), S475. (36) Montastruc, L.; Azzaro-Pantel, C.; Pibouleau, L.; Domenech, S. Use of Genetic Algorithms and Gradient Based Optimization Techniques for Calcium Phosphate Precipitation. Chem. Eng. Process. 2004, 43, 1289. (37) Ternet, D. J.; Biegler, L. T. Interior-Point Methods for Reduced Hessian Successive Quadratic Programming. Comput. Chem. Eng. 1999, 23, 859. (38) Schmid, C.; Biegler, L. T. Quadratic Programming Methods for Reduced Hessian SQP. Comput. Chem. Eng. 1994, 8, 817. (39) Sneesby, M. G.; Tade´, M. O.; Datta, R.; Smith, T. N. ETBE Synthesis Via Reactive Distillation. 1. Steady-State Simulation and Design Aspects. Ind. Eng. Chem. Res. 1997, 36, 1855. (40) Thiel, C.; Sundmancher, K.; Hoffmann, U. Synthesis of ETBE: Residual Curve Maps for the Heterogeneously Catalyzed Reactive Distillation Process. Chem. Eng. J. 1996, 66, 181. (41) Al-Arfaj, M.; Luyben, W. L. Control Study of Ethyl tert-Butyl Ether Reactive Distillation. Ind. Eng. Chem. Res. 2002, 41, 3784. (42) Jhon, Y.; Lee, T. Dynamic Simulation for Reactive Distillation with ETBE Synthesis. Sep. Purif. Technol. 2003, 31, 301. (43) Franç oisse, O.; Thyrion, F. C. Kinetics and Mechanism of Ethyl tert-Butyl Ether Liquid-Phase Synthesis. Chem. Eng. Process. 1991, 30, 141. (44) Jensen, K. L.; Datta, R. Ethers from Ethanol. 1. Equilibrium Thermodynamics Analysis of the Liquid-Phase Ethyl tert-Butyl Ether Reaction. Ind. Eng. Chem. Res. 1995, 34, 392. (45) Pibouleau, L.; Domenech, S.; Davin, C.; Azzaro-Pantel, Davin. Expe´rimentations Nume´riques sur les Variantes et Parame`tres de la Me`thode du Recuit Simule´. Chem. Eng. 2005, 105, 117. (46) Sneesby, M. G.; Tade´, M. O.; Datta, R.; Smith, T. N. Two - Point Control of a Reactive Distillation Column for Composition and Conversion. J. Process Control 1999, 9, 19.

1388


(47) Lee, J. H.; Dudukovic, M. P. A Comparison of the Equilibrium and Nonequilibrium Models for a Multicomponent Reactive Distillation Column. Comput. Chem. Eng. 1998, 23, 159. (48) King, C. J. Separation Processes; McGraw Hill Inc.: New York, 1980.

ReceiVed for reView April 14, 2005 ReVised manuscript receiVed September 5, 2005 Accepted November 30, 2005 IE0504506

A Mixed Integer Nonlinear Programming Formulation for Optimal

Recommend Documents