Efficient Optimization-Based Design of Distillation ... - ACS Publications

Mar 6, 2006 - 1.2. MINLP Optimization. The RBM is not applicable for detailed distillation column design, since the minimum energy demand is calculate...
2 downloads 0 Views 373KB Size
Download PDF
8492

Ind. Eng. Chem. Res. 2006, 45, 8492-8502

Efficient Optimization-Based Design of Distillation Columns for Homogenous Azeotropic Mixtures Sven Kossack, Korbinian Kraemer, and Wolfgang Marquardt* Lehrstuhl fu¨r Prozesstechnik, RWTH Aachen UniVersity, Turmstrasse 46, D-52056 Aachen, Germany

Design of distillation columns using mathematical programming has been a research topic for over a decade. However, a lack of robustness and high sensitivity toward the initial solution have prevented the widespread use. To overcome the numerical problems of the optimization algorithms, a combination of shortcut calculations with the rectification body method (RBM) and rigorous optimization is proposed. All potential splits can be generated and ranked based on the minimum energy demand calculated in the shortcut step. The least attractive split alternatives can be discarded to reduce the complexity of the rigorous model, while the results of the shortcut calculations for the promising alternatives can be used further for the initialization and bounding of the rigorous optimization. The rigorous optimization problems are then optimized as successive relaxed mixedinteger nonlinear programming (SR-MINLP) problems: instead of a combinatorial expensive MINLP problem, a succession of NLP problems with decreasing bounds is solved. An integer solution is achieved in all case studies considered. The robustness of the algorithm is improved, and the calculation times are reduced compared to existing methods. The design methodology is reliable and yields good results. It is illustrated with several case studies involving azeotropic mixtures. 1. Introduction The purification of reaction products in a chemical plant is typically accomplished through distillation. Especially for largescale separation of nonideal mixtures, distillation is still the method of choice, as stated by Widagdo and Seider.1 Evaluating and synthesizing the best possible distillation column setup has, therefore, generated a lot of scientific interest. Most of this previous research work, however, is related to nonazeotropic or even ideal mixtures. For the optimal design of distillation columns, a number of procedures have been suggested, which can be grouped into three main areas2: heuristics, shortcut calculations, and rigorous optimization. While heuristics, such as those given by Douglas,3 certainly are a valuable design technique, they cannot quantify the effort required for a desired separation task. Therefore, they do not replace quantitative studies. The separation cost is, in this case, usually determined by simulation studies, where the different feasible alternatives are optimized by successive simulation experiments. Shortcut procedures can be used to determine the minimum energy demand of a previously specified separation for zeotropic4 and azeotropic5-8 mixtures. They also allow an approximate evaluation of the feasibility of a desired split. The calculated minimum energy demand from these methods allows for a fast screening and ranking of the different design alternatives. These methods, however, usually calculate no detailed information on the optimal column design. This optimal design of a distillation column can be computed by rigorously optimizing a column model. Here, the design variables are usually optimized with mixed-integer nonlinear programming (MINLP) techniques,9-12 for which, in the general case, large superstructures10 need to be defined. These superstructures need to account for all possible splits and column configurations. The optimization runs are computationally * Corresponding author. Tel.: +49/241/80-94668. Fax: +49/241/ 80-92326. E-mail: [email protected].

expensive, and the quality of the final solution depends strongly11 on the initial values because of the nonlinearity and nonconvexity of the problem. These drawbacks of the rigorous MINLP optimizations have so far prohibited their use in industrial practice. A combination of shortcut calculations and rigorous MINLP optimization is proposed in this contribution in order to develop a robust method for distillation column design. This leads to a sequential procedure. In the first step, possible splits are generated and evaluated with the help of shortcut calculations. The results are then exported to a rigorous optimization model, where they are used for an efficient initialization and tight bounding of the optimization variables.13 This handover can be automated, so that no manual intervention is necessary once the shortcut calculations are completed. In a series of papers, Barttfeld and co-workers12,14,16 propose a similar method. These authors use the theory of reversible separation15 to generate initial values for the rigorous optimization. Their method is mainly limited by the drawbacks of this so-called “preferred separation”, because, for azeotropic mixtures, usually nonsharp splits are generated. An extension to the much more common sharp splits is not trivial. Further, their proposed initialization procedure cannot consistently handle separations that involve a tangent pinch situation. The rectification body method (RBM)7 is an extension of the preferred separation for sharp splits. Since the RBM is also applicable to splits in which a tangent pinch occurs, this versatile shortcut method is used to initialize and bound a rigorous optimization in this contribution. It is observed that the robustness of the rigorous MINLP optimization is greatly enhanced by the good initial values and tight bounds generated from the RBM. The required solution time is further reduced by using an efficient solution algorithm for the MINLP model. Only single columns with two products are considered here. The method can, however, be extended to more general cases and, in particular, to column sequences. 1.1. Rectification Body Method. The rectification body method (RBM) is used in the shortcut step to compute the minimum energy demand of the selected separation and to

10.1021/ie060117h CCC: $33.50 © 2006 American Chemical Society Published on Web 06/03/2006

Ind. Eng. Chem. Res., Vol. 45, No. 25, 2006 8493

initialize the MINLP optimization. The main ideas of the RBM will be briefly reviewed here; a more detailed presentation has been given by Bausa et al.7 According to Levy et al.,5 the tray-to-tray profiles of the rectifying and stripping section need to intersect for a feasible column design. The point of intersection belongs to both column sections and is located on the feed tray. Their boundary value method (BVM) calculates the respective profiles for an estimated energy demand by starting at the respective products. When the two profiles intersect, a feasible design has been found. This design method, however, is very sensitive toward trace impurities in the product compositions; therefore, a number of trayto-tray calculations are typically necessary to find the minimum energy demand of a separation. The BVM, therefore, cannot be automated and can only be used for ternary or quaternary mixtures, since the intersection of the profiles has to be checked visually. On the basis of the BVM, the zero-volume criterion (ZVC),6 an algebraic criterion for minimum reflux based on the alignment of pinch points, has been proposed. Since the fulfillment of the criterion can be checked algebraically, it can be applied to any number of components. The method is, however, only applicable to direct and indirect cuts, since the topology of the split has to be known a priori.7 The RBM uses the pinch points to approximate the column profile. The location of the pinch points is independent of the trace components and is mainly determined by the thermodynamic behavior of the mixture under consideration. After the pinch points have been calculated, the rectification bodies are determined by an automatic procedure. The resulting rectification bodies can be seen as a linear approximation of the envelope in the composition space of all possible column profiles. When the rectification bodies of the stripping and rectifying section intersect in a single point, the minimum energy of the separation has been found. The algorithm is summarized as follows. For a given feed and given products, the pinch lines are calculated. The locations of the pinch points on these pinch lines depend on the separation energy, which is varied automatically until the rectification bodies intersect in a single point. Since the intersection of the rectification bodies is checked with an algebraic criterion, no visual inspection is necessary and the method can be used for n-component mixtures. The RBM can be used for the generation of feasible splits and for a quick economic evaluation of the separation. Furthermore, the occurrence of a tangent pinch can be detected and handled automatically. More recently the RBM has been extended to complex columns17,18 and to extractive columns.19 1.2. MINLP Optimization. The RBM is not applicable for detailed distillation column design, since the minimum energy demand is calculated for a column with an infinite number of trays. After a first screening of design alternatives, the optimal tradeoff between investment costsswhich mainly depend on the size and number of trayssand the operating costsswhich are defined mainly by the cost of heating and cooling of the process streamsshas to be found. The optimal column design is typically computed through an optimization with an economic objective function.10,11,14 Since the number of trays and the feed tray location are discrete variables, an MINLP problem is solved in the design optimization. MINLP optimization of distillation columns using rigorous models has been introduced by Viswanathan and Grossmann.9 Since then, a number of extensions to the models and algorithms

Figure 1. Balance envelope for tray-to-tray and pinch calculations.

have been presented; however, mostly ideal mixtures have been considered. The most commonly used solution methods for these problems are decomposition algorithms such as outer approximation,20 where NLP and MILP problems are solved alternately, and branch-and-bound techniques21 with NLP subproblems. More recent modeling and solution techniques include general disjunctive programming22 or a continuous reformulation,23 where the discrete variable set in the MINLP is replaced by continuous variables which are forced to take discrete values by special constraints. Reliability of the MINLP optimization and convergence to good local optima or even to the global optimum still remain open issues. These are addressed in this paper through a novel solution approach, where a succession of relaxed MINLP problems with decreasing bounds is solved (see Section 3.2). 1.3. Overview on the Paper. In the following section, the models used in the shortcut calculations and in the rigorous MINLP optimization are discussed. The optimization procedure with successive relaxed MINLP (SR-MINLP) and its initialization are presented in Section 3. The method is illustrated through case studies with azeotropic mixtures in Section 4. In the concluding Section 5 of the paper, further possible extensions are outlined. 2. Column Models In this section, a short derivation of the models used in the shortcut step as well as in the rigorous optimization is given. The underlying assumptions and the relations between both model types are briefly outlined. 2.1. Shortcut Model. After the condenser or reboiler duties and the products of the column have been specified, the profiles of each column section can be solved if one starts at the respective column end. The balance envelope for the calculation of a tray-to-tray profile of the distillate section is shown in Figure 1. The equations for this section are24 v 0 ) -Lnhln + Vn+1hn+1 - DhD + QD, n ∈ N

(1)

0 ) -Lnxn,i + Vn+1yn+1,i - DxD,i, n ∈ N, i ∈ I

(2)

0 ) yn,i - Kn,i(x,y,T,p)xn,i, n ∈ N, i ∈ I

(3)

C

xn,i - 1 ) 0, ∑ i)1

n ∈ N, i ∈ I

(4)

n ∈ N, i ∈ I

(5)

C

yn,i - 1 ) 0, ∑ i)1

8494

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

where n ∈ N denotes the tray number and i ∈ I denotes the component index. Pinch points are defined as the fixed points of the tray-totray maps. They can be calculated from the pinch equations

0 ) -Lhl + Vhv - DhD + QD

(6)

0 ) -Lxi + Vyi - DxD,i, i ∈ I

(7)

0 ) yi - Ki(x,y,T,p)xi, i ∈ I

(8)

C

xi - 1 ) 0, ∑ i)1

i∈I

(9)

i∈I

(10)

C

yi - 1 ) 0, ∑ i)1

which can be derived from eqs 1-5 for a column with a large number of trays. The distribution coefficient Ki in eqs 3 and 8 is, at moderate pressures, defined by

γip0i , i∈I Ki ) pφ′′i

(11)

The liquid activity coefficient γi is calculated in this work by means of either the Wilson or the UNIQUAC model. The extended Antoine equation25 is used for the calculation of the pure-component vapor pressure p0i . Nonidealities in the vapor phase can be included by a model for the vapor-phase fugacity φ′′n,i; this is, however, not necessary for the mixtures considered here. The enthalpies needed in eq 6 are modeled by rigorous models such as the DIPPR equations. The rectification bodies are calculated from the pinch equations 6-10 in a very efficient way: the energy balance eq 6 is decoupled from eqs 7-10. The solution of this reduced set of equations are the pinch lines, on which the pinch points are located. The locations of the pinch points on these lines depend on QD. This parameter is varied in the intersection algorithm (see also Section 1.1) until the resulting bodies intersect in a single point. This resembles the intersection of the tray-to-tray profiles of the rectifying and stripping section, since the rectification bodies approximate the column profile. The energy determined in the intersection algorithm corresponds to the minimum energy demand of the specified separation. The pinch points are the fixed points of the tray-to-tray calculation. No separation progress is possible, once a pinch point is touched. These points are, thus, never reached in a column with a finite number of trays. The pinch points are, therefore, extreme points of a hypothetical column profile and bound a column profile for a column with a finite number of trays. 2.2. Rigorous Optimization Model. The RBM can be used to gain insight into feasible splits and to determine the minimum energy demand for a distillation column with an infinite number of trays. A second step is necessary to determine the optimal number of trays and the optimal feed tray location. Here, numerical optimization of a rigorous column model is used to compute the most cost-effective design of a distillation column. For a given product purity, the total annualized cost TAC, which combines investment cost Cinv and operational cost Copta on an annual basis, is minimized. The objective of the optimization is stated as

min TAC ) fcfLCinv + Copta,

(12)

with the capital charge factor fc, which annualizes one-time investment cost and other cost factors such as repair and maintenance, insurance, and overhead cost as well as interests and depreciation, and the Lang factor fL, which accounts for the direct cost associated with the installation of the equipment, including support structures, insulation, instrumentation, piping, etc. A value of 4 is assigned to the Lang factor fL and a value of 0.25 is assigned to the capital charge factor fc. The annual operation time ta is assumed to be 8000 h/a. The investment cost is calculated from the number of trays, the column diameter, and the heat exchanger areas, while the operating cost is determined by the energy cost. The details of the economic evaluation can be found in the Appendix. The optimization is constrained by a column model with a variable number of trays. The two different column superstructures used in this paper are shown in Figure 2. One column end, either the reboiler or the condenser, is fixed. For a column with a fixed reboiler (Figure 2a), all trays above the variable condenser are nonexistent, whereas in the case of a column with a fixed condenser (Figure 2b), all trays below the variable reboiler are nonexistent. Modeling the column and optimizing the model can be done with less nonlinear equations, if the purity constraint for a product is applied at a fixed tray. Hence, the superstructure is chosen such that the purity constraint for the key component of the separation is placed at the fixed column end. In the following, the mathematical representation of a distillation column is outlined for a superstructure with a fixed reboiler and variable condenser (see Figure 2a). The reformulation into a model with a fixed condenser tray and variable reboiler location (Figure 2b) is straightforward. The trays of the column are modeled as equilibrium trays as shown in Figure 3. In this model, the feed can be introduced on each tray and each tray can be cooled or heated by a heat exchanger. The products can leave the column at any tray. Accordingly, the material balance for a tray is given by

0 ) Fnzn,i + Ln-1xn-1,i + Vn+1yn+1,i - Lnxn,i - Vnyn,i Plnxn,i - Pvnyn,i, n ∈ N, i ∈ I (13) The tray energy balance is defined by l v + Vn+1hn+1 - Lnhln - Vnhvn - Plnhln 0 ) Fnhf + Ln-1hn-1

Pvnhvn + Qn, n ∈ N (14) The vapor-liquid equilibrium (VLE) (eq 3) and the closure relations for the liquid (eq 4) and vapor phase (eq 5) complete the equilibrium tray model. The overall column mass balance and the overall column energy balance are not needed to develop a valid column model, since the tray-by-tray material and energy balances are formulated consistently. It has been observed, however, that numerical stability is increased when the overall column material balance Nmax

Nmax

Nmax

∑ Fn ) n)1 ∑ Pln + n)1 ∑ Pvn n)1

(15)

is included in the optimization model as a redundant equation. The feed flowrate F ) ∑nFn and its compositions zi, i ∈ I, are known from the design problem specification. To meet the desired product specifications, a product purity constraint is

Ind. Eng. Chem. Res., Vol. 45, No. 25, 2006 8495

with F being the specified feed flow rate. The number of trays of the optimized column is then given by Nmax

NT ) Nmax -

Figure 2. Column superstructures: (a) column with a fixed reboiler (product purity constraint for the bottom product) and (b) column with a fixed condenser (product purity constraint for the distillate).

Figure 3. Single tray superstructure.

introduced. The liquid distillate product flow Nmax - 1

∑ n)1

Pln ) D

(16)

and the bottom product purity

xn,i g xb,i, n ) Nmax, i ∈ I

(17)

are set in the case where the model has a fixed reboiler, as considered here. In our work, only products leaving as a saturated liquid are considered. Therefore, we have

∀n ∈ N

Pvn ) 0,

(18)

Feasible product specifications can be determined, within the accuracy of the shortcut calculation, in the shortcut step. The discrete decisions regarding the existence of a certain tray or the feed tray location are modeled with binary variables and special constraints. These constraints involving binary variables for the location of the feed bfn and for the location of the condenser bcn are

Fn e

bfnF,

n∈N

(19)

Nmax

∑

bfn

) 1,

bfn

∈ {0, 1}

(20)

n)1

Qn e bcnQmax, n ∈ N

(21)

Nmax

∑ bcn ) 1,

n)1

bcn ∈ {0, 1}

(22)

bcnn - 1 ∑ n)1

(23)

Qmax in eq 21 needs to be a valid upper bound to the condenser duty. The exact value of Qmax is not known a priori. Therefore, it needs to be estimated, resulting in a weak relaxation of the optimization problem. As a consequence of this weak relaxation, the binary variables bcn are detached from the condenser duty Qn and from the liquid distillate product flows Pln when the relaxed optimization problem is solved. A value of zero for the binary variable bcn is equivalent to Qn ) 0, but a relaxed value of 0.3 for bcn is not equivalent to 0.3Qc, since Qc is overestimated by Qmax in eq 21. In a first step of the solution algorithm, the binary variables bfn and bcn are relaxed. Typically, the solution of the NLP resulting from the relaxed MINLP yields noninteger values for these variables. The value of the objective function is, therefore, a lower bound for the optimal integer solution. The locally optimal integer values for the binary variables are determined either through a branch-and-bound tree evaluation or through outer approximation iterations. According to our experience and as mentioned by Viswanathan and Grossmann9 and later by Du¨nnebier and Pantelides,11 the distribution of a single column feed is optimal when the entire feed is introduced on a single tray. The optimal solution of the relaxed problem will, therefore, determine an integer solution for the relaxed binary variable bfn. This optimal feed tray location minimizes the total number of trays.26 This can be shown in the binary case with a McCabeThiele diagram, where a feed with a fixed composition and thermal state is assigned a single optimal feed location. The optimal feed stage is located at the point where the operating line of the rectifying section crosses the so-called “q-line”, which denotes the thermal state of the feed. This single optimal feed point minimizes the total number of stages of the column, since here the switching of the active operating line maximizes the orthogonal distance to the VLE curve. Around the feed pinch zone in a multicomponent mixture, the situation is comparable; the graphical interpretation is, however, not extensible for mixtures with more than two components. Still the feed tray has to belong to both column sections for a feasible column. To maximize the number of effective trays in each section, the feed tray has to be placed farthest away from each column end and should, thus, not be distributed between several trays. Equations 19 and 20 can, therefore, be rewritten as

Fn e cfnF, cfn ∈R, n ∈ N

(24)

Nmax

∑ cfn ) 1,

cfn ∈ R

(25)

n)1

The continuous positive variables cfn denote the location of the feed. An integer value is obtained at the optimum. In the case studies of Section 4, we have observed that the local optima of the column optimization model with relaxed condenser location (0 e bcn e 1) are located at integer points. This is enforced by augmenting constraint 21 with the constraint

Pln e bcnD, n ∈ N

(26)

8496

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

relating the binary variables bcn to the distillate product flow rate D. This constraint is tighter than eq 21, because the upper bound, D, is either specified or defined by the overall column material balance. A relaxed value of 0.2 of the binary variable is equivalent to one-fifth of the total product flow. The detachment of the relaxed binary variables bcn denoting the condenser tray and the distillate product flows Pln is, thus, not possible. Hence, the total number of stages is now connected to the stage where the product of that respective column end is drawn off. A distributed product draw-off cannot be optimal, since the product purity is maximized at the extreme column ends. A distributed product draw-off would lead to a mixing of streams of different purity, which would then require a higher reboiler energy or an increased number of stages. Since integer solutions are reached for the relaxed binary variables, eqs 21-23 and 26 are reformulated with continuous variables ccn to result in

Qn e ccnQmax, ccn ∈R, n ∈ N

(27)

Nmax

ccn ) 1, ∑ n)1

ccn ∈R

Pln e ccnPlnmax ) ccnD , ccn ∈R, n ∈ N

(28) (29)

Nmax

NT ) Nmax -

ccnn - 1, ∑ n)1

ccn ∈R, n ∈ N

(30)

Integer values are obtained for the continuous variables for the location of the condenser ccn at the optimum. The rigorous optimization problem now only contains continuous variables and is, thus, no longer a MINLP problem. As a consequence, no computationally expensive tree search or decomposition algorithm is needed. Furthermore, NLP solvers can be used for the optimization of this problem. Note that no nonlinear and nonconvex expressions forcing the continuous decision variables to take on binary values have been introduced, which would have been the case for a continuously reformulated MINLP.23 Here, a tight model formulation identifies the local optima of the optimization problem at integer points due to the nature of the problem. The continuous column model consisting of eqs 3-5, 11-18, 24, 25, and 27-30 exhibits considerably shorter solution times than the equivalent MINLP model as demonstrated in the case studies in Section 4. The difference between the solution times is expected to become even larger, when combinatorially more demanding problems such as complex columns or column sequences are optimized. 3. Solution Procedure The solution procedure for the optimization of distillation columns in this work can be broadly divided into RBM calculation, initialization, and rigorous optimization. The RBM has already been discussed in the previous Sections 1.1 and 2.1. In the following Section 3.1, the initialization of the rigorous optimization with information from the RBM calculation is shown. The novel solution procedure for the rigorous optimization is presented in Section 3.2. For an illustration of the initialization algorithm, a ternary mixture containing acetone/ methanol/ethanol as presented in the first case study (see Section 4.1 for further details) is considered. Note that the RBM and

Figure 4. Rectification bodies and linearized column composition profile (chain dotted line) used for the initialization.

the optimization model are applicable to mixtures with an arbitrary number of components. 3.1. Initialization of Rigorous Optimization. The pinch points calculated with the RBM can be used to approximate the initial profile for the rigorous model of the distillation column. The compositions and temperatures of the pinch points as well as the minimum energy demand are exported to GAMS, where the rigorous optimizations are performed. A number of different initialization strategies can be thought of if the pinch points are known. A very simple procedure, where the tray compositions and temperatures are initialized with a piecewise linear function through the pinch points (shown in Figure 4), works sufficiently well. As in the RBM, the pinch points are connected in the order of increasing stability.7 Still, reasonable estimates are required for the initial feed tray and the location of the pinch points. These trays need to be assigned by the designer. After the initialization with the linear column profile, an initial optimization of the column is performed, as shown in Figure 5. In this initial optimization, the feed tray location is optimized by minimizing the reboiler duty

MinQB cfn

(31)

The number of trays is fixed to the maximum number of trays by assigning the condenser to the first tray:

ccn ) 1, n ) 1 ccn ) 0, 2 e n e Nmax

(32)

All other column model eqs 3-5, 11, 13-18, 24, 25, and 27-30 are included in this optimization. Note that the product specification eq 17 is also included in this model. The initial feed tray optimization derives a feasible column which meets the purity specifications. This increases the robustness of the optimization of a column with a variable number of trays. In particular, the impact of the initial choice of the feed tray and the trays associated with the pinch points on the quality of the solution of the final optimization problem is minimized. The bounding of the variables, especially of the tray temperatures and enthalpies, is as critical as the initialization. These key variables greatly influence calculation of the physical

Ind. Eng. Chem. Res., Vol. 45, No. 25, 2006 8497

Figure 6. Big-M values and objective function value in the successive SR-MINLP optimization procedure.

The relaxation of the continuous optimization model is realized by reformulating the constraint in eq 29 with a big-M relaxation to

1 ccn D e Pln e ccnMD, n ∈ N M Figure 5. Solution strategy for the MINLP optimization.

properties. The bounds for the temperature and enthalpy are set to the maxima and minima of the considered system. Note that tighter temperature and enthalpy bounds yield suboptimal solutions, since the temperatures and enthalpies on nonexisting trays can exceed the expected distillate and bottom values during the solution process. The minimum energy demand of the separation as determined by the RBM can be used as an excellent initial value for the reboiler duty in the rigorous optimization. 3.2. Successive Optimization Procedure for the Rigorous Continuous Model. Since the model used for the feed tray optimization in Section 3.1 only differs in two aspects, namely, in the evaluation of the objective function eq 31 and in the fixed number of trays eq 32, from the full optimization model, all final values of the variables are used to initialize the rigorous column optimization. It should be noted that this includes the feed tray location calculated in the previous step, which is not fixed in the full optimization model. The rigorous continuous optimization model as stated in Section 2.2, consisting of eqs 3-5, 11-18, 24, 25, and 27-30, can be solved with common NLP solvers such as SNOPT27 or CONOPT.28 Integer values for the decision variables cfn (feed tray) and ccn (number of trays) are computed in the locally optimal solutions for the case studies presented in Section 4. However, the obtained solutions proved to be local optima of low quality: they are considerably worse than the global optimal solution, which can be identified by enumeration. The lower solution quality can be attributed to the tight model formulation, which aggravates the impact of the initial values on the optimization solution. This influence of the initial guess on the final result is reduced by a successive solution procedure. A series of successive relaxed MINLP (SR-MINLP) is solved with decreasing bounds, representing a tightening model formulation. The idea is to further relax the continuous problem in a first solution step in order to create initial values for a second solution step with tighter bounds. The bounds are then further tightened in successive steps to yield the tight model formulation defined in Section 2.2, which exhibits integer solutions for the decision variables.

(33)

with M being a sufficiently large constant. The relation between the continuous decision variables ccn and the distillate product flows Pln is relaxed for M > 1. The relaxed continuous solution is a lower bound to the tightest formulation (M ) 1), which determines an integer solution. To relax the model further, the constraint defined in eq 27 connecting the variables ccn with the condenser duty Qn is omitted in the relaxed steps of the optimization procedure. This additional constraint is added again in the final step of the optimization. The optimization problem defined by eqs 3-5, 11-18, 24, 25, 28, 30, and 33 is then solved in successive optimization steps, where the big-M values M in eq 33 are subsequently tightened as illustrated in Figure 6. The big-M value is set to a rather large value in the first step of the sequence to facilitate convergence of the optimization problem and to generate a favorable noninteger starting point for the subsequent steps. The bounds are then tightened by decreasing the big-M values in the successive steps. The value of the objective function rises with tighter constraints, as can be seen in Figure 6. The successive relaxation is not sensitive to the choice of the big-M constants as long as the difference between successive values is not too large. In the penultimate step of the sequence, the big-M value is set to 1, which results in the tightest possible model formulation. The constraint defined in eq 33 is now equivalent to the constraint defined in eq 29. The local optima of the decision variables are located on integer points. Constraint 27 is added to the model in a last step, forcing the condenser to be located on the same tray as the distillate product. Note that the optimal feed tray and the optimal number of column trays does not change after constraint 27 has been included in the last step. 4. Illustrative Case Studies The rigorous optimization of single columns is illustrated with three case studies. In all case studies, the feed flow rate is set to 10 mol/s and a saturated liquid feed is considered. An overview of the feed compositions for the case studies presented is given in Table 1. In addition, the parameter sources and the models used to calculate vapor-liquid equilibrium in eq 3 are given. The thermodynamic relations for the calculation of the vapor-liquid equilibrium introduce highly nonlinear equations to the optimization problem. To reduce the size of the optimiza-

8498

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

Figure 7. Acetone/methanol/ethanol mixture: (a) residue curves and (b) mass balance line, rectification bodies, and column profile. Table 1. Case Studies Considered in This Paper components

feed composition

thermodynamic model

parameter source

acetone/methanol/ethanol methanol/ethanol/2-propanol/water acetone/ethanol/water

0.35/0.15/0.5 0.4/0.2/0.3/0.1 0.0564/0.0637/0.8799

UNIQUAC UNIQUAC Wilson

ASPEN (2004) ASPEN (2004) Knapp and Doherty (1992)

Table 2. Optimization Specifications and Results for the Case Studies mixture

purity specification

QB,min (kW) RBM

number of trays (feed tray)

QB (kW) opt.

TAC

total CPU time (s)

AME (Sec. 4.1) MEPW (Sec. 4.2) AEW (Sec. 4.3)

xB,ethanol g 0.999 xD,methanol g 0.999 xB,water g 0.9999

479.25 540.3 174.52

48 (18) 60 (28) 41 (39)

521.96 638.82 211.42

EUR 211927 EUR 253454 EUR 94249

26.26 86.81 27.7

tion problem in the optimization software and to enhance the flexibility to choose more complicated activity coefficient models (such as UNIQUAC), external functions are used in GAMS to calculate these properties and the required derivatives. While the robustness increases, the communication overhead between the solver and the external function raises the required solution time. The pressure of the columns is fixed at atmospheric pressure (1.013 bar) although pressure loss correlations can be introduced into the model, since the tray pressure is defined as a variable. The case studies presented in this section are solved in seven successive steps with the big-M values M of eq 33 being gradually reduced: The big-M values are fixed at 200, 30, 10, 5, 2, 1.5, and 1 in successive steps, as shown in Figure 6. The successive solution steps are computationally easy to solve (the solution times for each step can also be found in Figure 6), since the optimization models in each step are identical, except for the tightened big-M values. The successive steps are, thus, initialized close to the optimal solution of each optimization step by the solution of the previous step. With this successive solution procedure, the global optimum is found in the first case study, as can be shown by enumeration of all possible solutions. All other reported optima can be expected to be of a very good quality; however, it cannot be guaranteed that the global optimum is found with this solution procedure, since it relies on local optimization of nonconvex problems. The solution algorithm as well as the presented case studies have been implemented in GAMS 21.7.29 The optimizations are

solved with the NLP solver SNOPT 6.2-1 on a workstation with a 3.06 GHz CPU and 3.75 GB RAM. The solution times for the case studies proved to be very short, as documented in Table 2, where an overview of the calculated results is given. 4.1. Ternary Azeotropic Mixture. A ternary mixture of acetone/methanol/ethanol is considered first. As shown in Figure 7a, this mixture exhibits a minimum-boiling binary azeotrope. This azeotrope is the light boiler or unstable node of the residue curve map (RCM). Ethanol is the heavy boiler, corresponding to a stable node of the RCM. All residue curves originate, therefore, at the acetone/methanol azeotrope and end at the pure ethanol vertex. The residue curve originating from an infinitesimal distance of the azeotrope is shown as a solid line in Figure 7a. All residue curves above this line pass the acetone vertex, which represents a saddle of the residue curve map. For all residue curves below this line, the saddle is the methanol vertex. Ethanol is to be separated from acetone and methanol from the feed with a specified composition of xacetone ) 0.35, xmethanol ) 0.15, and xethanol ) 0.5. The preferred distillation products12 for this feed would result in a nonsharp split between a binary acetone/methanol mixture and a binary methanol/ethanol mixture. The distillate has, in the case of the preferred separation, a higher concentration of acetone than the azeotropic composition. In the specification of the sharp separation, however, a concentration lower than the azeotropic composition is required. Since all ethanol has to be separated from the feed, a bottom product purity of xethanol g 0.999 and a bottom product flow rate of 5 mol/s are specified. The minimum reboiler demand as calculated by the RBM is 479.25 kW for this split. As can be

Ind. Eng. Chem. Res., Vol. 45, No. 25, 2006 8499

Figure 8. Results for the acetone/methanol/ethanol case study.

seen in Figure 7b, the column composition profile calculated in ASPEN for a fairly large column (60 trays overall and feed introduced on tray 30) with the minimum energy is approximated very well by the rectification bodies. The profile passes through both products as well as the feed pinch. It should be noted that the occurrence of a feed pinch is not necessary for the RBM (as opposed to the reversible distillation method), since tangent pinches can be detected and handled appropriately by the RBM, as will be demonstrated in Section 4.3. The initialization and optimization procedure as outlined in the previous sections is applied. The results of the optimization steps are presented in Figure 6. As can be seen in Figure 8, the identified final solution is also the global optimal solution. It is identified after 26 s of CPU time. The optimal total tray number is 48, and the feed is introduced on tray 18 of the optimized column. The optimal reboiler duty is 521.96 kW. The total annualized cost, representing the optimal value of the objective function, is EUR 211 927. When the model is solved as a MINLP with the branch-and-bound solver SBB30 and the NLP solver SNOPT starting from the same initial conditions, the global optimum is found; however, the solution time increases to 240 s, which is almost 10 times the time required for the successive solution procedure. All additional results in Figure 8 are calculated by fixing the total number of trays and only solving for the optimal feed location. In the vicinity of the global optimum, only small differences between the local optimal solution exists; therefore, any solution between 45 total trays (TAC ) EUR 212 663) and 51 total trays (TAC ) EUR 212 385) seems to be an acceptable result. While the reboiler and condenser duties increase with lower tray numbers as expected, investment cost first decreases with lower tray numbers (until about 41 total trays) and then starts to increase again. The increased vapor load in the column, due to a larger reboiler energy, causes an increase in the tray diameter. The individual trays are, therefore, more expensive. 4.2. Quaternary Example. The RBM and the optimization model presented in the previous sections are applicable to mixtures with an arbitrary number of components. A quaternary azeotropic mixture is chosen in this second case study to highlight this fact. A mixture of methanol, ethanol, 2-propanol, and water is considered. The feed composition is given by xmethanol ) 0.4, xethanol ) 0.2, x2-propanol ) 0.3, and xwater ) 0.1. A study of the residue curve map reveals that only methanol can be recovered with high purity in a single column. The optimal column design for this separation, which is equivalent to a direct split, is determined. The RBM calculates for this separation an energy demand of 540.3 kW. The pinch points are exported to GAMS for the initialization procedure. The distillate product flow is set to 4

Figure 9. Acetone/ethanol/water: rectification bodies (solid lines) and final composition profile (chain dotted line).

mol/s. Since the purity constraint of xmethanol g 0.999 is set for the distillate product, a superstructure with a fixed condenser, as shown in Figure 2b, is used. The solution procedure presented in the previous section is followed, and an integer solution is found in 87 s. This seems to be the global optimum, since no better solution is found for different initial values; however, a complete enumeration of all possible solutions has not been done. Note that, because of the nonconvex model and the use of a local NLP solver, global optimality cannot be guaranteed. The optimal column consists of 60 active trays. The optimal feed is at tray 28. A value of EUR 253 454 for the objective function and a reboiler duty of 638.82 kW is determined in the optimization. 4.3. Mixture Exhibiting a Tangent Pinch. As mentioned in Section 1, the initialization procedure proposed by Barttfeld and Aguirre12 based on reversible separation theory cannot consistently handle systems in which a tangent pinch occurs. The initialization procedure proposed in this paper does not have this limitation. An example for a separation, in which a tangent pinch occurs, is the first column of a pressure-swing process for the separation of ethanol and water. As proposed by Knapp and Doherty,31 acetone is added to the mixture to enhance the pressure-swing area. Here, only the first column, which is operated at atmospheric pressure, is considered. The product of this column is water at very high purity with the specification of xwater g 0.9999 at the bottom. Since the bottom flowrate is set to 8.717 mol/s, the top product still contains some water; the split in this column is, therefore, a nonsharp split. As can be seen in Figure 9, the rectification body in the rectifying section of the column reduces to a straight line. Other calculated pinch points in this section, which cannot be reached by a column profile, are discarded.32 The actual column profile does not resemble a line but rather shows a curvature between the end points of the rectification body, as can be seen in Figure 9. The initialization routine as well as the optimization procedure can still be applied for this column. A sensitivity toward the initial choices for the tray numbers associated with the pinch points is observed. An integer solution is found within 28 s CPU time. In this case, the optimal column has 43 trays and the feed is introduced on tray 39. The value of the objective function at the optimum

8500

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

is EUR 94 249, and a reboiler duty of 211.42 kW is determined by the optimization.

with c1 ) EUR4100, c2 ) EUR1800, and c3 ) EUR3100 taken from Brusis.33 The height of the column Hcol is computed from

Hcol ) NTHtray + H0,

5. Summary and Conclusion In this contribution, the optimal design of distillation columns is calculated through a combination of a shortcut model and a rigorous model. The RBM is used to identify feasible products and to determine the minimum energy demand of the separation. Calculation results from the initial shortcut step are then used to initialize and bound the rigorous model. The integration of the rigorous optimization with the shortcut calculations increases the robustness in two ways: it is not necessary to define a complicated superstructure for the rigorous optimization, since the split generation and evaluation can be done in the shortcut step. Only the most promising alternatives need to be optimized in the rigorous optimization step, which can be initialized and bounded very well using the results of the shortcut step. The optimization gives reliably good results. The initialization of the rigorous optimizations with the rectification body method does not suffer from the limitations of a previously presented initialization strategy12 based on reversible distillation theory. Mixtures which exhibit a tangent pinch can be handled consistently. Furthermore, the product purities can be set to values above 99.9%, even for azeotropic mixtures. The solution of the MINLP problems is done in a novel way. Instead of separating the MI from the NLP part, as is usually done in MINLP solvers, a succession of relaxed MINLP problems is solved with successively tightened bounds. The optimal solutions of the column model are located on integer values for all case studies considered. If local optima are not located on integer values, continuous reformulation23 could be employed to enforce an integer solution. This will retain the advantage that NLP solvers can still be used to solve these problems at the expense of additional nonconvexity. The successive relaxed integer optimization (SR-MINLP) is a very fast and robust solution technique, which will be further explored in our future research. To start the SR-MINLP optimization directly from the shortcut calculations, the initialization and bounding procedure will be automated. This integrated method is also suited for the design of column sequences. The extension to a series of columns without recycle streams seems to be straightforward. If column sequences with recycles are considered, an intermediate optimization with shortcut methods as suggested by Bru¨ggemann32 should provide additional robustness. In addition, the successive solution procedure of relaxed MINLP optimization models is expected to reduce the combinatorial complexity when flowsheets with more than one column are considered.

with Htray ) 0.5 m being the height of a column tray. The estimated clearance H0 ) 4 m accounts for liquid distributors, etc. The column diameter Dcoln at tray n can be computed from

Dcoln )

xx

RT

4Vn

Mi ∑ iy n,i

πF

p

Dcol ) Dcoln)Nmax ) Dcoln ∀n

(A1)

Cs ) c1DcolHcol

(A2)

Cint ) c2DcolHtrayNT

(A3)

Chex ) c3(Areb + Acon)0.65,

(A4)

(A6)

(A7)

While the actual diameter of the column might differ slightly from the diameter at the chosen tray, this simplification should not affect the quality of the solution in a noticeable way. The heat exchange areas for the reboiler and condenser can be obtained from

∑

Qn)Nmax n*Nmax , Acon ) Areb ) ∆Trebk Qn∆Tconk

(A8)

Since the temperature of the bottom product and the temperature of the steam are fixed, the temperature difference at the reboiler ∆Treb is also constant. The temperature difference in the condenser, ∆Tcon, can be calculated from the logarithmic mean temperature. Fixed values are assumed here, since the temperature of the reboiler or condenser is set by the product specifications of the problem. The mean heat transfer coefficient k is assumed to be 800 W/(m2 K) for both heat exchangers.33 The operational costs for heating and cooling add up to

Cop ) ccwmcw + cstmst

(A9)

The cost for cooling water at 298 K and reboiler steam at 3 bar is fixed at ccw ) 0.05 EUR/t and cst ) 10 EUR/t, respectively (cf. ref 33). The mass flow of cooling water mcw and mst can be derived from

∑

mcw ) -

n*Nmax

Qncpw∆Tcw Qn)Nmax

Cinv ) Cs + Cint + Chex

, n∈N

where the gas load F is set to a value of 2xPa. In general the column diameter should be chosen to be the diameter at the widest tray; however, such a formulation proved to add to the nonconvexity of the problem and led to suboptimal solutions. Therefore, the column diameter is fixed to the diameter of the reboiler tray

Appendix: Economic evaluation The investment cost for column shell Cs, column internals Cint, and heat exchanger equipment Chex can be estimated by the correlations3

(A5)

mst )

rst

(A10)

(A11)

The heat capacity of the cooling water at 298 K is cpw ) 4.183 kJ/(kg K), and the heat of vaporization of the steam at 3 bar is rst ) 2163.2 kJ/kg. ∆Tcw is the maximum temperature difference of the cooling utility. Acknowledgment The authors would like to thank their former colleague Stefan Bru¨ggemann for fruitful discussions on the subject of this work. The financial support from Deutsche Forschungsgemeinschaft

Ind. Eng. Chem. Res., Vol. 45, No. 25, 2006 8501

under the Gottfried Wilhelm Leibniz program is gratefully acknowledged. Nomenclature A ) area (m2) B ) bottoms flow (mol/s) b ) binary variable C ) cost (EUR) C˙ ) cost per hour (EUR/h) c ) continuous variable/heat capacity [kJ/(kg K)]/cost parameter D ) distillate flow (mol/s)/diameter (m) F ) feed stream (mol/s)/gas load (xPa) f ) function of/factor for objective function H ) height (m) h ) enthalpy (kJ/kmol) K ) equilibrium constant k ) heat transfer coefficient (W/m2‚K) L ) liquid stream (mol/s) M ) Big-M constant/molar mass (kg/kmol) m˘ ) mass stream (kg/s) N,n ) tray number NT ) total number of trays Q ) energy (W) P ) product stream (mol/s) p ) pressure (Pa) R ) gas constant (kJ/kmol‚K) r ) reflux ratio/heat of condensation (kJ/kmol) T ) temperature (K) t ) time (h) TAC ) total annualized cost (EUR) V ) vapor stream (mol/s) x ) liquid composition y ) vapor composition z ) total composition Greek Letters γ ) activity coefficient ∆ ) difference φ ) fugacity coefficient Subscripts 0 ) column without internals a ) annualized b ) bottom col ) column con ) condenser cw ) cooling water d ) distillate hex ) heat exchanger i ) component i int ) internals inv ) investment n ) tray number max ) maximum min ) minimum op ) operating p ) constant pressure r ) reflux reb ) reboiler s ) shell st ) steam tray ) tray w ) water

Superscripts 0 ) pure component c ) (partial) condenser f ) (partial) feed l ) liquid low ) lower bound up ) upper bound v ) vapor Literature Cited (1) Widagdo, S.; Seider, W. D. Azeotropic distillation. AIChE J. 1996, 42, 96. (2) Wahnschafft, O. M. Synthese Von Trennsystemen fu¨r azeotrope Gemische unter besonderer Beru¨cksichtigung Von Rektifikationsprozessen, Fortschr.-Ber. VDI: Reihe 3, Nr. 311; VDI Verlag: Duesseldorf, Germany, 1993. (3) Douglas, J. M. Conceptual Design of Chemical Processes; McGrawHill: New York, 1988. (4) Underwood, A. Fractional distillation of multicomponent mixtures. Chem. Eng. Prog. 1948, 44, 603. (5) Levy, S.; Van Dongen, D.; Doherty, M. Design and synthesis of homogeneous azeotropic distillations. 2. Minimum reflux calcualtions for nonideal and azeotropic columns. Ind. Eng. Chem. Fundam. 1985, 24, 463. (6) Julka, V.; Doherty, M. Geometric behavior and minimum flows for nonideal multicomponent distillation. Chem. Eng. Sci. 1990, 45, 1801. (7) Bausa, J.; von Watzdorf, R.; Marquardt, W. Shortcut methods for nonideal multicomponent distillation. 1. Simple columns. AIChE J. 1998, 44, 2181. (8) Liu, G.; Jobson, M.; Smith, R.; Wahnschafft, O. Shortcut design method for columns separating azeotropic mixtures. Ind. Eng. Chem. Res. 2004, 43, 3908. (9) Viswanathan, J.; Grossmann, I. Optimal Feed Locations and Number of Trays for Distillation Columns with Multiple Feeds. Ind. Eng. Chem. Res. 1993, 32, 2942. (10) Bauer, M.; Stichlmair, J. Design and economic optimization of azeotropic distillation processes using mixed-integer nonlinear programming. Comput. Chem. Eng. 1998, 22, 1271. (11) Du¨nnebier, G.; Pantelides, C. Optimal design of thermally coupled distillation columns. Ind. Eng. Chem. Res. 1999, 38, 162. (12) Barttfeld, M.; Aguirre, P. Optimal synthesis of multicomponent zeotropic distillation processes. 1. Preprocessing phase and rigorous optimization for a single unit. Ind. Eng. Chem. Res. 2002, 41, 5298. (13) Bru¨ggemann, S.; Marquardt, W. Shortcut methods for nonideal multicomponent distillation: ProViding initials and bounds for rigorous optimization; Technical Report; Lehrstuhl fu¨r Prozesstechnik: RWTH Aachen, Germany, 2001; http://www.lpt.rwth-aachen.de/Publication/abstract. php?Nummer)LPT-2001-16. (14) Barttfeld, M.; Aguirre, P. Optimal synthesis of multicomponent zeotropic distillation processes. 2. Preprocessing phase and rigorous optimization of efficient sequences. Ind. Eng. Chem. Res. 2003, 42, 3441. (15) Koehler, J.; Aguirre, P.; Blass, E. Minimum Reflux Calculations for Nonideal Mixtures Using the Reversible Distillation Model. Chem. Eng. Sci. 1991, 46, 3007. (16) Barttfeld, M.; Aguirre, P.; Grossmann, I. Alternative representations and formulations for the economic optimization of multicomponent distillation columns. Comput. Chem. Eng. 2003, 27, 363. (17) von Watzdorf, R.; Bausa, J.; Marquardt, W. Shortcut methods for nonideal multicomponent distillation. 2. Complex columns. AIChE J. 1999, 45, 1615. (18) Bru¨ggemann, S.; Marquardt, W. Rapid screening of design alternatives for nonideal multiproduct distillation processes. Comput. Chem. Eng. 2005, 29, 165. (19) Bru¨ggemann, S.; Marquardt, W. Shortcut methods for nonideal multicomponent distillation. 3. Extractive distillation columns. AIChE J. 2004, 50, 1129. (20) Duran, M.; Grossmann, I. An outer-approximation algorithm for a class of mixed-integer nonlinear programs. Math. Programming 1986, 36, 307. (21) Gupta, O.; Ravindran, V. Branch and bound experiments in convex nonlinear programming. Manage. Sci. 1985, 31, 1533. (22) Tu¨rkay, M.; Grossmann, I. Logic-based MINLP algorithms for the optimal synthesis of process networks. Comput. Chem. Eng. 1996, 20, 959. (23) Stein, O.; Oldenburg, J.; Marquardt, W. Continuous reformulation of discrete continuous optimization problems. Comput. Chem. Eng. 2004, 28, 1951.

8502

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

(24) Bausa, J. Shortcut-methods for Conceptual Design and Thermodynamic Analysis of Distillation Processes, Fortschr.-Ber. VDI: Reihe 3, Nr. 692; VDI Verlag: Duesseldorf, Germany, 2001. (25) ASPEN Technology. ASPEN Physical Property System 12.1; ASPEN Technology, Inc.: Cambridge, MA, 2003. (26) Vogelpohl, A. Die optimale Lage des Zulaufbodens bei der Vielstoffrektifikaion. Chem.-Ing.-Tech. 1975, 47, 895. (27) Gill, P.; Murray, W.; Saunders, M.; Drud, A.; Kalvelagen, E. SNOPT: An SQP algorithm for large-scale constrained optimization. In GAMSsThe SolVer Manuals; GAMS Development Corporation: Washington, DC, 2005; p 393. (28) Drud, A. CONOPT. In GAMSsThe SolVer Manuals; GAMS Development Corporation: Washington, DC, 2005; pp 39-82. (29) Brooke, A.; Kendrick, D.; Meeraus, A.; Raman, R. GAMSsA Users Guide; GAMS Development Corporation: Washington, DC, 2005. (30) Drud, A. SBB. In GAMSsThe SolVer Manuals; GAMS Development Corporation: Washington, DC, 2005; p 377.

(31) Knapp, J.; Doherty, M. A new pressure-swing-distillation process for separating homogeneous azeotropic mixtures. Ind. Eng. Chem. Res. 1992, 31, 346. (32) Bru¨ggemann, S. Rapid Screening of Conceptual Design AlternatiVes for Distillation Processes, Fortschr.-Ber. VDI: Reihe 3, Nr. 841; VDI Verlag: Duesseldorf, Germany, 2005. (33) Brusis, D. Synthesis and Optimisation of Distillation Processes with MINLP Techniques, Fortschr.-Ber. VDI, Reihe 3 Nr. 797; VDI Verlag: Duesseldorf, Germany, 2003.

ReceiVed for reView January 27, 2006 ReVised manuscript receiVed April 12, 2006 Accepted April 13, 2006 IE060117H