Temperature Collocation Algorithm for Fast and Robust Distillation

Ind. Eng. Chem. Res. 2004, 43, 3163-3182

3163

Temperature Collocation Algorithm for Fast and Robust Distillation Design Libin Zhang and Andreas A. Linninger* Laboratory for Product and Process Design, Department of Chemical Engineering, University of Illinois at Chicago, Chicago, Illinois 60607

In this paper, we describe a new minimum bubble-point distance (MIDI) algorithm for assessing the feasibility of a desired distillation specification. The algorithm computes the rectifying and stripping temperature profiles by temperature collocation on finite elements with orthogonal polynomials. We discovered the beneficial use of a dimensionless equilibrium tray temperature as an independent variable. This novel choice is bounded between 0 and 1, improves the numerical quality of the design problem formulation, and is well-behaved even in the vicinity of pinch regions. It also employs the fixed points of column sections as collocation points. Adaptiveelement boundary placement at saddle temperatures can effectively overcome problems with numerical instability near pinch regions. Extensions of our MIDI algorithm to the calculation of the minimum and maximum refluxes are also introduced. We show the application of this algorithm in the separation of quaternary mixtures and provide an outlook of the methodology for optimal column sequencing. Cases studies demonstrate the algorithm’s robustness and reliability. 1. Introduction Distillative separation is among the least expensive methods for separating mixtures that exhibit suitable volatility differences.1 Hence, continuous distillation columns, along with their optimal operation and heat integration, constitute a major engineering activity in refineries and bulk commodity manufacturing. The common practice for distillation design often involves numerous trial-and-error experiments by means of state-of-the-art process flowsheet simulators. However, this time-consuming practice does not guarantee the production of successful designs; it might not provide any information about the feasibility of a particular specification in cases when the efforts do not converge. When the design-by-simulation approach is used, infeasibilities are often discovered only after extensive simulation studies.2 The classical Underwood method for column profile computation is restricted to mixtures of constant relative volatility.3,4 Doherty and co-workers5,6 introduced the boundary value method (BVM), which examines the intersection of rectifying and stripping profiles graphically. However, the BVM is inconvenient for mixtures with more than four components because of the lack of graphical representations for composition trajectories in higher dimensions. Julka and Doherty7,8 extended this methodology to multicomponent systems, employing a geometric theory based on topological concepts. They demonstrated that the feed point and C-1 pinch points lie in one hyperplane satisfying a zero-volume formula at minimum reflux. The rectification body (RB) method9 employed pinch points of the top and bottom sections to approximate the space of reachable column profiles. If the rectification bodies for the rectifying and stripping section just touch each other, the authors conjectured that the minimum reflux ratio is obtained. However, these methods are only accurate for high* To whom correspondence should be addressed. Tel.: (312) 996-2581. Fax: (312) 996-0808. E-mail: [email protected].

purity splits in which the composition profiles approach saddle pinches. They break down for sloppy separations in which the intersection of hyperspaces spanned by the pinch point is a necessary but not a sufficient feasibility condition. In general (sloppy) specifications, only column profiles unambiguously verify the feasibility of a design. Existing approaches avoid calculating the composition profile rigorously because of the huge computational effort. This holds particularly true when several columns are being designed and optimized simultaneously (i.e., in optimimum column sequencing). Orthogonal collocation on finite elements (OCFE) has been shown to reduce the problem size substantially, without significant loss in accuracy compared to the full-order trayby-tray model. Several researchers have explored collocation techniques for distillation column simulation.11-20 An OCFE method was also deployed successfully for minimum-stage-number designs.21-23 Huss and Westerberg24,25 proposed advantageous variable transformations to ensure high accuracy close to the columns’ end sections when using collocation. Despite the attention devoted to continuous distillation simulations, less work has been aimed at developing algorithms to determine whether a given specification, sloppy or sharp, is feasible or not. In this research, we propose a minimum bubble-point distance (MIDI) algorithm for ascertaining the feasibility of any arbitrary design specification for simple continuous distillation columns. Our research objective targets the development of an efficient and reliable computational algorithm for establishing the feasibility or infeasibility of specifications for sharp or sloppy separations of any number of species. The availability of a fast and globally convergent feasibility test would introduce an essential element for automatic computer-aided distillative separation synthesis.26 This work constitutes an important intermediate milestone toward that long-term objective. Outline. A theoretical basis for the novel robust feasibility test algorithm is developed in section 2. Fundamental concepts pertinent to the new approach

10.1021/ie034223k CCC: $27.50 © 2004 American Chemical Society Published on Web 05/15/2004

3164 Ind. Eng. Chem. Res., Vol. 43, No. 12, 2004

such as the attainable temperature window (ATW) and the bubble-point distance (BPD) are introduced. A temperature finite-element collocation method capable of overcoming stationary points in the composition profile is presented. The methodology section also includes numerical results and implementation issues. Section 3 demonstrates the robustness and performance characteristics of our approach for constant relative volatility mixtures, ideal mixtures, and nonideal mixtures with three or four components. Section 4 extends our methodology to ascertain the boundaries of feasible operation, i.e., it describes the robust computation of minimum and maximum reflux ratios. Section 5 addresses problematic design tasks such as the separation of mixtures with tangent pinches and provides suggestions toward novel computational solution options for structural design problems such as optimal column sequencing. The paper closes with conclusions and an outlook for future research. 2. Methodology A performance problem predicting expected distillate and bottoms compositions given the feed and column parameters for an existing column always has a solution that is readily attainable with commercial flowsheet packages (e.g., AspenPlus, HYSYS, Pro/II, etc.). On the other hand, the inverse design problem seeking the column operation and parameters to achieve desired separation targets might not have a solution even for a consistent set of design specifications. In this paper, we shall, without loss of generality, assume given feed compositions and column pressure. A simple column configuration has four design degrees of freedom, typically three product-purity specifications and a desired reflux value, for example. Unfortunately, no generally applicable, globally convergent algorithm exists for determining the feasibility of a given design specification. The classical Underwood method3,4 converges only for feasible designs; it is of limited use in the synthesis of distillation trains because of its extremely nonlinear behavior. Other design methodologies are restricted in either the number of species27 or the properties of the mixture.28 To the best of our knowledge, no algorithm is robust enough for structural flowsheet optimizations, which typically lead to largescale mixed integer nonlinear mathematical programming (MINLP) problems. Critical numerical difficulties in the design problem stem from singularities in the composition profiles near stationary pinch and saddle points that naturally occur even in ideal mixtures. To overcome the existing shortcomings, we propose a combination of two simple but very effective concepts: (1) transformation of column profiles into the space of dimensionless bubble-point temperatures; (2) minimization of a bubble-point distance function between stationary profile nodes in the bubble-point temperature space. This novel approach reduces the dimensionality of the design problem; eliminates singularities encountered in the tray-by-tray approach; extends to any number of species with customary vapor-liquid equilibrium solution models, including constant-relative-volatility, ideal, and nonideal mixtures; and applies to both sharp and sloppy splits. In effect, the column design problem becomes more tractable from a computational point of view. Before the introduction of our new methodology, a discussion of reachable product compositions is in order.

Figure 1. Fixed points in a ternary mixture and residual errors of eqs 1 and 2.

2.1. Reachable Compositions and Their Equilibrium Temperatures. Pinch points delineate the extreme points of all possible column profiles for a given design specification. For a known distillate composition xD for a c-component separation, c stationary points are found without performing tray-by-tray computations by enforcing the pinch equations given by eqs 1 and 2, i.e., simultaneous compound balances and thermodynamics equilibrium.29

Pinch equation for the rectifying section rxi - (r + 1)xiKi + xD,i ) 0

∑xiKi - 1 ) 0

i ) 1, ..., C

(1) (2)

Figure 1a illustrates that the distillate, d, and stable node, γ, span all possible rectifying profiles for a given specification. Rectifying composition profiles start at the distillate and terminate in the stationary pinch point. For high-purity separations, the profiles approach a saddle point dividing the composition profile into two branches that can be reached only after an infinite number of equilibrium trays have been traversed.5 Each pinch point is also associated with a specific bubble-

Ind. Eng. Chem. Res., Vol. 43, No. 12, 2004 3165

stripping and rectifying sections, we can already exclude many design specifications without actually performing tray-by-tray computations as expressed in remark 1. This property can also be very valuable in assessing the separability of mixtures with azerotropes. Remark 1. Attainable Temperature Windows. If the attainable temperature window (ATW) of the rectifying and stripping section is empty, the given design specification is infeasible. 2.2. Bubble-Point Distance (BPD). In cases of overlapping ATWs, more detailed quantitative analysis is required. Numbers of trays or column heights are unsuitable independent variables in numerical computations because of their indefinite value range, namely, ]0, ∞[. We propose an affine mapping of the column height into a dimensionless bubble-point temperature, Θ. The mapping establishes a bijective relationship between height and bubble-point temperature. We also convert the tray-by-tray difference equations into ordinary differential equations as proposed by Doherty.5 The composition trajectories expressed as functions of the new independent variable, Θ, do not approach infinity hampering the traditional tray-by-tray methods. Using the dimensionless temperature instead of the column height, we arrive at a system of differential-algebraic equations. Equations 3 determine the composition profiles, xi(Θ), as a function of the equilibrium constant, Ki; the distillate specification, xD,i; and the reflux, r. The new temperature transformation reduces the nonlinearity of the profile equations significantly, because the use of temperature as an independent variable makes temperature-dependent expressions explicit. Equation 4 represents the stripping profile in term of the bottoms composition, xB,i, and the reboiler ratio, s. The detailed derivations of continuous rectifying and stripping profiles as well as the total derivative of the equilibrium constant, dKj/dΘ, in terms of temperature and composition for a general nonideal mixture model are provided in the appendices.

Continuous rectifying profile equations

Figure 2. Attainable temperature windows in (top) composition and (bottom) temperature space.

point temperature as illustrated in Figure 1b. Consequently, realizable rectifying sections exhibit tray equilibrium temperatures within the interval delimited by the temperatures of the distillate, TD, and the stationary pinch point, TP1. Correspondingly, the bubble points of the bottoms, TB, and the stationary stripping pinch, TP2, bracket the temperature range for all possible stripping sections. For a simple column configuration (i.e., with no side trays or multiple feeds), the stripping and rectifying profiles must intersect at the liquid tray composition of the feed plate. Hence, the temperature interval between min(TP1, TB) and max(TP2, TD) constitutes the attainable temperature window (ATW) for the given column design specification. The ATW marks all possible equilibrium tray temperatures that can be reached by the column sections, as illustrated in Figure 2 (i.e., boiling points of mixtures on an equilibrium tray). If the ATW associated with a design specification is empty (i.e., TP1′ < TP2′,), then one can conclude unequivocally that the specification is infeasible (e.g., Figure 2b). By comparing the temperature ranges in the

dxRi

( ∑[(

) ∑( ) )]

r+1 R 1 yi + xD,i r r

xRi -

)-

dΘ

C

r+1 R 1 yj + xD,j Kj r r R R F [x (Θ),Θ,p]

C

xRj -

j)1

j)1

dKj xRj ) dΘ

i ) 1, ..., C (3)

Continuous stripping profile equations dxSi dΘ

)-

( ∑[(

s

ySi - xSi +

s+1

C

j)1

s

s+1

ySj - xSj +

1

xB,i

s+1 1

S

C

j)1

dKj xSj ) dΘ

xB,j Kj

s+1 F [x (Θ),Θ,p] S

) ∑( ) )]

i ) 1, ..., C (4)

The bubble-point distance (BPD) is defined as the Euclidean distance between two equilibrium compositions belonging to the stripping and rectifying profiles with the same bubble-point temperature. Figure 3 depicts two instances of bubble-point distances, dA and dB, for two discrete bubble-point temperatures (TA and TB, respectively), as well as the traces of the bubblepoint surfaces in a nonideal ternary mixture. For a


Figure 3. Bubble-point distance and boiling-point trajectories in a nonideal mixture.

Figure 5. Illustration of shortcut feasibility for a ternary mixture.

Figure 4. Bubble-point distance function (BPD) for a feasible design (+, lower curve) and an infeasible design (*, upper curve) in a ternary separation.

design specification to be feasible, there must exist a pair of compositions with BPD ) 0. If the globally minimum BPD for a specified separation target is greater than zero, then this target is guaranteed infeasible. Figure 4 contrasts a feasible solution and an infeasible design of a ternary separation problem. The infeasible specification has a minimum BPD that is greater than zero. Our feasibility test also exploits remark 2 for ternary mixtures: Remark 2. Shortcut Infeasibility. For a feasible ternary column to exist, it is necessary and sufficient that the vectors P1-Xr(TP2) and P2-Xs(TP1) intersect, where Xr(TP2) is the point on the rectifying profile corresponding to the temperature of the stable stripping pinch P2 and

Xs(TP1) is the point on the stripping profile corresponding to the temperature of the stable rectifying pinch P1. Intuitively, we can then restate remark 2 as follows: for column profiles to be feasible, the linear approximation of the rectifying and stripping profiles taken between the two pinch-point temperatures must intersect as depicted in Figure 5. The validity of remark 2 becomes clearer when one considers strict temperature monotonicity of the composition profiles; this property is observed in industrially relevant columns. Exceptions to this rule are possible but require special heat effects such as can occur in reactive or extractive distillations. These situations are not the objective of this work. Hence, profiles are assumed not to intersect the same bubble-point surface twice even in highly nonideal solutions. For mixtures of four or more species, the generalization of remark 2 to hypersurfaces is necessary, but not sufficient. BPD for More than Four Species. The design specifications for the separation of a quaternary mixture do not fully determine the products. In this case, the feasibility test and determination of the remaining


stripping profiles as formulated in eqs 5-9 below. This problem can also be considered a “dynamic” optimization problem in which the dimensionless bubble-point temperature, Θ, assumes the role of the independent variable “time”. The equality constraints of eq 6 correspond to the rectifying and stripping profiles of eqs 3 and 4. Equality 7 represents the pinch equations yielding the dimensionless temperature boundary for this dynamic optimization problem.30 Inequality constraint 8 ensures that the dimensionless temperature, Θ, remains within the bounds of the attainable temperature window (ATW). This mathematical program always has a solution, even for infeasible separation targets. The final output of the feasibility test returns the actual column profile or impossibility of realizing the desired separation with simple distillation within the precision limits of the algorithm. Figure 6. Bubble-point distance surface for a feasible separation of a quaternary mixture. Table 1. Pseudocode for the Feasibility Test (MIDI Algorithm) 1. 2. 3. 4. 5.

Compute the stable pinch points in the rectifying and stripping sections, P1 and P2 (eqs 1 and 2) Compute the points Xs(P1) and Xr(P2) corresponding to the pinch point on the opposite profile. IF the shortcut is infeasible, THEN the column specifications are infeasible. Stop. Perform global minimization of the BPD using orthogonal collocation on finite elements. IF min ||BPD|| < , THEN column specs are feasible. Stop. ELSE design specifications are infeasible. Stop.

compositional degrees of freedom must be performed simultaneously. Figure 6 depicts a feasible solution for a quaternary problem in a two-dimensional problem space consisting of the bubble-point temperature and remaining compositional degree of freedom, i.e., distillate composition, xD,3. Again, the minimum BPD is zero for a feasible specification. It is worth noting the smooth BPD functions in both cases; this property of the bubblepoint distance function is essential for numerical algorithms requiring gradient information. Furthermore, the BPD surface, as well as all projections of the BPD for a given xD,3 value exhibits only one minimum. 2.3. Minimum Bubble-Point Distance (MIDI) Algorithm. Feasible separations targets lead to intersecting stripping and rectifying profiles, whereas these profiles miss each other in infeasible designs. Enforcing the equality of the feed tray composition in the top and bottom sections is not a good basis for a numerical algorithm, as it cannot be converged for infeasible specifications. To ensure numerical convergence for any situation, we propose to search for the globally minimum BPD instead. The information flow for the new feasibility test based on the minimum bubble-point distance (MIDI algorithm) is summarized in Table 1. The first step searches for the stable pinch points for the stripping and rectifying profile. Steps 2 and 3 examines the shortcut infeasibility (remark 2), which involves the liquid tray compositions corresponding to the bubble points of the stationary nodes X(Θ1) and X(Θ2) of the opposite column section. For ternary mixtures, the shortcut infeasibility of remark 2 can detect many infeasible specifications without the need for detailed profile computations. Step 4 entails a nonlinear optimization problem seeking the global minimum distance between the rectifying and

Feasibility Test: Mathematical Problem Formulation min ||BPD|| s.t.

(5)

dx(Θ) ) F[x(Θ),Θ,p] dΘ

(6)

G[x(Θ),Θ,p] ) 0

(7)

Θl e Θ e Θu

(8)

0 e x(Θ) e 1

(9)

Θ,p

2.4. Finite-Element Approximation of Column Profiles. As stated above, the mathematical program of eqs 5-9 constitutes a dynamic optimization problem with differential algebraic (DAE) constraints.31 For its solution, we approximated the differential equations with finite temperature elements and collocation on orthogonal polynomials.23,32,33 A detailed discussion of orthogonal collocation on finite elements (OCFE) for discretizing differential and partial differential equations can be found in the literature.22,34 This method is also referred to as spline collocation.35 The continuous differential equations are converted into a set of nonlinear algebraic equations in terms of the unknown weights, i.e., unknown compositions. The OCFE discretization renders algebraic expressions of the variable derivatives for each collocation node as described in eq 10. With the help of gradient expression in eq 10, the differential equations become algebraic as shown in eq 11.

Discretization of derivatives of node j in element [i]

|

dx ) AΘ[i]x j dΘ Θ)Θ[i]j

(10)

Transformation of profiles in eqs 6 via point collocation discretization of node j in element [i] AΘ[i]x ) F(x,Θ,p) j

(11)

[i] , represent The elements of the variable vector xj,c the unknown compositions for each species, i.e., the equilibrium tray composition of each species c at node j in element i corresponding to temperature, Θ[i] j . The finite elements decompose the temperature range into segments; the nodes of the orthogonal collocation poly-

3168 Ind. Eng. Chem. Res., Vol. 43, No. 12, 2004 Chart 1

nomials further divide each element according to a predetermined number of roots. Within each element, the variable derivatives are expressed in terms of a derivative matrix, A, of the collocation polynomials and the unknown liquid tray compositions, x, each one corresponding to the specific node temperature. Zero-order continuity is ensured from element to element. The detailed expression resulting from the OCFE discretization on dimensionless temperature using Lagrangian and Legendre base polynomials for the rectifying and stripping profiles are given in eqs 12-17. The scalar entries, λk(Θ[i] j ), of the derivative matrix, A, in eq 12 represent the first derivative of the kth collocation polynomial, lk, evaluated at node j in element [i], i.e., [∂lk(Θ[i] j )]/∂Θ. The elements of the right-hand-side vec[i] tor, gj,c , involve the expressions and bc, where ξ[i] j represents the right-hand sides of the profile equations, i.e., eq 14 for the rectifying section and eq 16 for the stripping profile. The term bc accounts for the product purity, i.e., eqs 15 and 17. Substitution of the approximating functions into eqs 3 and 4 yields a corresponding set of residual functions that can be solved to obtain the desired liquid concentration trajectories in each column section. It is important to stress that the solvability of the resulting design equations is substantially enhanced because of the explicit expressions in dimensionless temperature as well as the boundedness of all dependent and independent variables between 0 and 1. [i] [i] [i] gj,c ) ξ[i] j xj,c + bcξj

(13)

Right-hand-side expression for the rectifying profile

(

1-

ξ[i] j ) -

C

∑

m)1

r + 1 [i] Kj,c r

)

× r + 1 1 [i] [i] [i] [i] xj,m xk,m + xD,m Kj,m Kj,m r r C dK[i] j,m [i] (14) xj,m m)1 dΘ

[(

∑

1 bc ) xD,c r

(

) ]

)

(15)

Right-hand-side expression for the stripping profile

ξ[i] j ) -

(

s

C

∑

m)1

[(

s

)

[i] -1 Kj,c

s+1

[i] [i] [i] xj,m - xj,m + Kj,m

s+1

1

[i] xB,m Kj,m s+1 C dK[i] j,m [i] (16) xj,m m)1 dΘ

∑

bc )

1 x s + 1 B,c

) ]

×

(

)

(17)

The implementation of the MIDI algorithm is discussed in the next subsection. 2.5. Implementation of the MIDI Algorithm. Pinch Equations. The solution of the pinch equations for the computation of all stationary points in step 1 can be executed with well-known homotopy continuation methods.8,36 Determining all pinch points in nonideal mixtures is a multivariable algebraic problem because of the dependence of equilibrium constants on composition as well as temperature. Different versions of the MIDI algorithm identify zeros by continuation or interval search. The computation of stationary points for ideal mixtures simplifies to a one-dimensional zero search, because the equilibrium constant depends only on the temperature. Bisection and bounded one-dimensional Newton methods are successful for ideal mixtures. For constant-R mixtures, the problem becomes a simple polynomial root finding. This can be addressed by the classical Newton-Horner scheme.37 The MIDI algorithm also benefits from the fact that unstable nodes lying outside the temperature window spanned by the product compositions are not needed. Choices for Node Placements. The solution accuracy benefits from node placement in regions with steep gradients. It is worthwhile mentioning the advantage of element boundary placement at saddle-point temperatures. We observed that the saddle-point temperature nearly coincides with the maximum curvature of the composition profile in most cases. This aspect holds particularly true in sharp-split separations, as demonstrated in Figure 7. In consequence, it is advisable to place an element boundary at the saddle-point

Ind. Eng. Chem. Res., Vol. 43, No. 12, 2004 3169 Table 2. Initialization Methods Tested for the Case Study A initialization method

speed to solve (s)

performance

finite-difference Underwood ODE

0.10 0.82 1.26

occasional convergence failure always convergent always convergent

Table 3. Comparison of Performancea of Different Numerical Method

example ternary mixtureb quaternary mixturec min/max reflux ratiod column sequencee

golden Newtonhomotopy section Raphson continuation search method method rSQP C C C NA

C F C C

C F C C

C F C C

a C ) converges, F ) fails occasionally. b Ternary mixture example from Figure 9. c Quaternary mixture example from Figure 12. d Min/max reflux ratio example from Figures 16 and 17. e Column sequence example from Figure 22.

Figure 7. Saddle temperature coincident with region of maximum profile curvature.

temperature. Numerical problems in the column profile computation exhibiting long, flat profile sections can be circumvented by reducing the number of elements. Initialization. The influence of the initial values on the convergence and computational effort required to find an accurate solution was investigated. Three options were explored: (a) Underwood initialization, (b) ODE initialization, and (c) finite-difference initialization. For the Underwood initialization, the equilibrium relation is approximated by a geometric mean of the relative volatilities at the bottoms and distillate temperatures. The initial weights at specified node temperatures were extracted by means of the Underwood method using a pseudo-temperature (described in more detail in section 3). In the ODE initialization, we directly integrate the profile with a fourth-order explicit RungeKutta method37 under the ideal-mixture assumptions expressed in eqs 3 and 4. This initialization technique already yields the correct result for ideal mixtures. The finite-difference initialization obtains an approximate composition at the next node by employing a one-stage explicit Euler step from the previous node. For all zeotropic mixtures, Underwood initialization is recommended, because of its balance between speed and reliability. The ODE initialization proved to be the most

accurate, but it is also most expensive in computational time. The finite-difference initialization is fastest, but it often leads to insufficiently accurate profiles, potentially causing numerical problems in the subsequent collocation procedure. This initialization method can be applied when many elements and nodes are used, so that the temperature gap is small. The computational performance and reliability for these three different initialization methods are reported in the Table 2. We found the Underwood initialization using the geometricmean constant relative volatility to exhibit the best overall result. Computation of Minimum Distance. The numerical solution of the NLP in eqs 5-9 still poses a formidable challenge because the interesting solutions are often sought at extreme concentrations. Several gradient-based and informed-search techniques were investigated for performance and robustness. For the solution of the nonlinear profile equations, we tested the Newton-Raphson method with Armijo line search,37 the autoscaled Newton solver using dynamic column and row scaling,38 the inexact Newton method based on Krylov subspace iterations,39 the Broyden method37 and the trust region solver (MINPACK).40 The homotopy method41 was deployed for the fixed-point equations (eqs 1 and 2). For the optimization problem of the minimum BPD, the golden section search, the Newton-Raphson method,37 and the rSQP method,42,43 in combination with the Newton-type method for the OCFE, were deployed in the case studies. We investigated the Newton-type and reduced successive quadratic programming (rSQP) methods for the simultaneous solution of the dynamic optimization problem as presented in eqs 5-9. For ternary mixtures, the direct dynamic programming succeeded, but this approach often failed in quaternary problems. The performance results are summarized in Table 3. We discovered that the BPD distance function admits multiple extrema whenever saddle-point temperatures fall within the ATW. To locate multiple local minima, MIDI launches multiple searches for each temperature subinterval created by saddle-point temperatures within the ATW. In the case of multiple local convergence, the smallest distance is declared to be the “globally” minimum BPD. In our experience with design problems in ternary and quaternary mixtures, MIDI proved reliable,44 because it warps the design problem into a


relatively small space in terms of the independent variables (i.e., Θ, compositional degrees of freedom). Alternatively global search algorithm based on interval arithmetic could be used. All solvers were implemented in C or Fortran on a Pentium II 450-MHz computer. For industrial users, a user-friendly graphical interface was developed. Users can input the specifications and select from among all solution algorithms described above. We also provide graphical interfaces to track the convergence trajectory of the numerical procedures.45 In a preliminary study, the gradient information in the optimization problems was obtained by numerical differentiation. Gradient information obtained by numerical differences can fail in quaternary separation problems because the gradients are often fairly flat along the free composition coordinates (e.g., xD,3). This fact is quite clear in Figure 6. Inaccurate gradient information often leads to near singular matrices, causing the mathematical program to diverge. These problems can be circumvented by using first-order sensitivity information, which is superior to numerical differences, as will be shown in case study E.46 3. Design Applications and Validation This section demonstrates the application of the new feasibility test for a wide range of zeotropic mixtures for both sloppy and sharp splits. For constant-R mixtures, the MIDI algorithm provides a more robust and faster alternative to the Underwood method. The MIDI algorithm extends seamlessly to ideal and nonideal mixtures. The examples described here characterize and validate the robustness and performance of the novel feasibility test. 3.1. Case Study A: Constant Relative Volatility Mixtures. For mixtures with constant relative volatilities, the Underwood method can be used to assess feasible column specifications. Attacking design problems with the Underwood method requires the simultaneous solution of the Underwood roots, an n-dimensional root-finding problem, together with the search for the numbers of trays in the rectifying and stripping section for profile intersection.3,4 For more than four species with sharp desired product splits, this problem is often numerically unstable and has no solution for infeasible column specifications. As a consequence, infeasibility, divergence due to sharp split specifications, and insufficiently accurate initialization are numerically indistinguishable. Hence, the Underwood method does not constitute a robust feasibility test for computational purposes in separation synthesis. As a solution, we advocate the modified MIDI algorithm. The MIDI algorithm requires the introduction of an artificial bubble-point temperature, θ′, for mixtures with constant relative volatility by eq 21. For converting the constant-R equilibrium condition into the format required by the MIDI algorithm, we define the purecomponent vapor pressure in relation to an artificial bubble-point temperature, θ′. The pseudo-vapor pressure, Ψi, which is linear in the bubble-point temperature θ′, as expressed in eq 20, is compatible with the equilibrium condition for constant relative volatility in eq 18, as well as with the ideal solution model in eq 19. In effect, one obtains linear pressure-temperature expressions with a temperature scale chosen between the boiling points of the lightest and heaviest com-

pounds as a reference. With this modification, the MIDI feasibility test is adapted for constant relative volatility mixtures.

yC )

RCxC C

(18)

Rixi ∑ i)1 ψC x P C

(19)

P [(R - 1)θ′ + 1] R1 1

(20)

yC ) ψC ) 1

C

[(R1 - 1)θ′ + 1]

R1

Rixi - RC ) 0 ∑ i)1

(21)

A computation experiment documented in Figure 8 applied the feasibility test to 10 000 randomly chosen column design problems for a ternary mixture with constant relative volatility. All design problems converged within a few iterations; no convergence failure was reported. The entire execution time amounted to 39.4 CPU s on a Pentium II 450-MHz computer for 10 000 design computations! Figure 8 also visualizes all feasible designs in a composition triangle and lists the number of infeasible designs. More design experiments for mixture with R ranging from 1 to 1000 are documented elsewhere.44 The algorithm was found to work robustly and reliably in all constant-R mixtures. Failures were observed only when all species had extremely close volatilities (e.g., RA ) 1.1, RB ) 1.05, RC ) 1). The results of the ternary case matched the analytical results obtained by Doherty and Malone.29 These results provide evidence of the robustness and speed of the new approach. 3.2. Case Study B: Ideal Mixtures. The MIDI algorithm for mixtures with ideal vapor-liquid behavior was formulated using eqs 11-17. The equilibrium constant, Ki, is a function of the pure-component vapor pressure only, which are conveniently modeled by the Antoine equation.47 It is recommended that analytical derivatives available at specialized high-performance physical property servers48 be used. The simulation results of Figure 9 depict product compositions and corresponding column profiles obtained by the new feasibility test. The rectifying and stripping profiles in Figure 9a corresponding to sloppy separation targets run far from the saddle point as expected. In sharp separations, as depicted in Figure 9b, the profiles approach the saddle pinch requiring a huge number of trays but only a few temperature nodes. These results demonstrate the speed and reliability of the feasibility test based on the novel idea of temperature collocation for ideal mixtures in sharp and sloppy splits. Figure 10 reports typical results of a simulation experiment with 10 000 randomly chosen column design problems. Comparison of the feasible designs with the feasibility boundaries established by Doherty and Malone29 shows good agreement. More computational results show stable and rapid performance for a wide range of vapor pressure differences. Numerical failures occurred when attempts were made to design separations for extremely narrow boilers such as isomers.


Figure 8. Feasible designs (right) filtered from 10 000 random specifications (left) (R1 ) 6.35, R2 ) 2.47, R3 ) 1).

However, for narrow boilers, distillation is, in all likelihood, not an economical process option. 3.3. Case Study C: Nonideal Zeotropic Mixtures. The MIDI algorithm extends easily to nonideal mixtures. The vapor-liquid data models employed the Antoine equation for the vapor pressure and the twoparameter Wilson equation for liquid-phase activities.47 Analytical derivatives of the activity coefficients with respect to temperature were computed by means of a high-speed thermodynamic property server.48 The collocation equations for nonideal mixtures are analogous to those presented for ideal solutions, except for expressions modeling the equilibrium constant, Ki. They introduce the species activities and their temperature dependence in the profile equations as quantified in eq 22.

(

s 1 dpi ) γi + psi dΘ P dΘ

dKi

(∑

∂γi dxj ∂γi ‚ + ∂Θ j)1 ∂xj dΘ C

))

(22)

When eq 22 is used to include the composition dependence of the equilibrium constant, all equations presented for ideal mixtures extend to nonideal mixtures. Figure 11 reports all feasible specifications upon solution of 10 000 feasibility tests for the separation of the nonideal mixture of acetaldehyde, methanol, and water. The execution time of 1000 s to perform 10 000 random nonideal column design problems is still reasonably fast. Extensive design problems for the separation of ideal and nonideal solutions at different feed compositions and a wide range of volatilities are discussed elsewhere. 44 Whereas feasible regions can be determined analytically for ternary mixture by computations based on pinch regions,27,49 extensions to higher dimensions do not exist. The MIDI algorithm does not suffer from this limitation.

3.4. Case Study D: Quaternary Mixtures. So far, we have shown the utility of the feasibility test for ternary mixtures. In this section, the method is extended to four components, a problem for which analytic methods for determining feasible regions do not exist. In mixtures with four or more components, the design problem involves the search for one or more compositional degrees of freedom (DOF) of the distillate or bottom product.7 The problem becomes a two-dimensional search in the quaternary case, e.g., temperature and one compositional degree of freedom. The mathematical formulation is identical to the problem in eq 5, when one compositional DOF is considered as an additional unknown parameter p, e.g., the distillate purity of a species, xD,3, without loss of generality. The problem is again a dynamic optimization problem in which temperature plays the role of “time”.50 In nearideal quaternary mixtures, there are four fixed points with two unstable nodes in each column section. The bubble-point temperatures corresponding to the two saddles are chosen as element boundaries. Figure 12 displays typical results for quaternary separation problems. These diagrams are designed to document the applicability of the novel approach to quaternary mixture of constant R, ideal mixtures and nonideal mixtures. Other experiments cover a variety of problems including different feeds, direct and indirect splits, and different production specifications. Exhaustive simulation results for a variety of nonideal mixtures as well as a description of the software providing an easy-touse interface for industrial users are documented elsewhere.44 The previous sections have demonstrated the suitability of the MIDI algorithm as a reliable and fast feasibility test for ternary and quaternary mixtures. The next section extends the proposed method for the computation of the feasible operational range, i.e., minimum and maximum reflux.


Figure 10. Experiments with 10 000 random designs for assessing the feasibility of separating an ideal mixture of pentane, hexane, and heptane.

Figure 9. Feasible composition profiles of an ideal mixture: (a) sloppy split specification, (b) sharp split specification.

4. Feasible Range of Column Operation: Minimum and Maximum Reflux During the early stages of distillation design, it is often useful to ascertain the range of refluxes capable of realizing desired product purities. Existing graphical methods for determining the minimum reflux are limited to binary and ternary mixtures.51-53 Shortcut methods apply only to constant relative volatility mixtures3,4,54,55 and/or sharp split assumptions.10 In general, there exist definite lower and upper reflux bounds for feasible column operation. It is possible to relate the problem of feasible refluxes to the feasibility test discussed previously, when considering three product specifications as given and the reflux ratio as an open variable. Then, the problem can be restated as the search for the interval [rmin, rmax] for which the feasibility test is positive.

Figure 11. Experiments with 10 000 random designs for assessing the feasibility of separating a nonideal mixture of acetaldehyde, methanol, and water.

4.1. Minimum Reflux Ratio. The compact mathematical formulation of our new feasibility test allows for the accurate computation of the minimum feasible reflux ratios for desired product purities. A new method based on the feasibility test for minimum reflux ratios makes use of remark 3: Remark 3. Minimum Reflux Ratio, rmin. A reflux specification leading to a zero bubble-point distance (BPD) between the rectifying and stripping profiles at the stationary pinch point is the minimum.


Figure 12. Feasible specification of quaternary mixtures. (a) Constant-R mixture [A(1), B(2), C(3), D(4)] (RA ) 8, RB ) 4, RC ) 2, RD ) 1. (b) Ideal mixture [pentane (1), hexane (2), heptane (3), octane (4)]. (c) Nonideal mixture [methanol (1), ethanol (2), n-propanol (3), acetic acid (4)].

For a ternary mixture, three possible pinch scenarios are depicted in Figure 13. The desired value corresponds to the smallest of the three possible cases. We propose a two-level optimization algorithm for calculating the minimum reflux ratio, as shown in Figure 14. The inner loop executes the feasibility test using the MIDI algorithm as discussed in the previous sections. The outer loop modifies the reflux ratio until

the limit of feasible operation is identified within desired accuracy. A modified golden section method was employed to endow the algorithm with robustness.37 Figure 15 depicts the information flowsheet of the outer-loop optimizer for updating the reflux ratio. We bracket the reflux ratio between 0 and 1 using the dimensionless reflux F ) r/(r + 1). The initial guesses are typically given at 0.01, 0.44, 0.66, and 0.999 as the starting


Figure 14. Overview of the extended MIDI algorithm for computing minimum and maximum reflux ratios.

Figure 13. Three different pinch scenarios at minimum reflux of a ternary constant-R mixture.

interval. The search produces a series of contracting bounds until rmin is found. The MIDI algorithm enforces

the feasibility test for each specification in the inner loop. Two branches can be met depending on whether a feasible specification was found in the interval. Branch 1 determines whether none of the four values in the current interval is feasible. It compares the signs of the gradients of the minimum distances to determine the next step according to two possible scenarios. (1a) When the signs of the gradients are the same, the procedure terminates and concludes that the specifications are infeasible. (1b) If the signs differ, a new search is launched within the subinterval corresponding to refluxes with sign changes as the new upper and lower boundaries. The range they span is divided consistent with the golden section rule. Branch 2 sets the smallest feasible reflux as the upper bound and subdivides the adjacent interval according to the golden cut. Finally, the minimum reflux ratio is obtained when all boundaries are within a small interval of length . rmin is equal to the last upper bound, which is a feasible reflux. Typical results for two ideal ternary mixtures of pentane, hexane, and heptane and methanol, ethanol, and n-propanol for sloppy and sharp separations are depicted in Figure 16. Table 4 contrasts the outcomes of other minimum-reflux algorithms for a variety of separation problems to those of our temperature collocation approach. The numerical results of all methods are almost identical for the sharp-split specification, as expected. Figure 16a and b shows the correct minimumreflux profiles obtained with the MIDI algorithm in all circumstances. These diagrams also explain the deviations befalling pinch alignment methods, as the pinch points and the feed points for sloppy specifications are


Figure 15. Detailed information flow diagram for computing the minimum reflux ratio based on the MIDI algorithm.

not collinear at minimum reflux. This shortcoming also occurs for pinch alignment in nonideal mixtures, for which MIDI is still accurate. We found previous methods to exhibit errors amounting to as much as 20% compared to the correct minimum reflux obtained with the MIDI algorithm. 4.2. Maximum Reflux Ratio. The maximum reflux is the largest reflux ratio for which a simple column specification can be realized. Unless the distillate and bottoms happen to lie on the same residue curve, which would be an unusual coincidence, a given separation target is realizable only up to a particular reflux limit, rmax. This explanation also makes clear why total reflux does not always render the best product splits. A reflux r is at its feasible maximum for a given separation target if the rectifying profile goes through the bottoms product or if the stripping profile penetrates the opposite liquid equilibrium point of a tie line through the distillate composition. This is expressed in remark 4. Remark 4. Maximum Reflux Ratio, rmax. A reflux leading to a zero bubble-point distance (BPD) between the rectifying and stripping profiles at the bottoms product temperature (bubble point) or distillation product temperature (dew point) is the maximum. Figure 17 shows the solution to the maximum reflux problem for a direct or indirect split of the ideal mixture pentane, hexane, and pentane obtained by the MIDI algorithm. The temperature collocation algorithm of this paper provides the first computational procedure for computing the maximum reflux ratio without resorting to graphical means. We should also note that the distillate composition is in phase equilibrium with the liquid composition on the top stage. The stripping profile passes through a liquid composition point in equilibrium

with the desired distillate composition xD, as depicted in Figure 17b. Figures 18 and 19 document the successful extension of the new approach to the calculation of the minimum and maximum reflux ratios for a nonideal quaternary separation of methanol (1), ethanol (2), n-propanol (3), and acetic acid (4). The activities were modeled with a two-parameter Wilson solution model with Wohl’s expansion.47 The minimum reflux ratio is reached when the minimum BPD equals zero at the stationary node, i.e., the stable node pinch on the opposite column profile (Figures 18a and 19a). For the maximum reflux, the composition profile passes through the product (Figures 18b and 19b). Although graphical representations are not available for mixtures with more than four components, the computational approaches described in this paper extend to multicomponent problems. 5. Limitations and Future Research Directions This section points out difficult situations related to tangent pinches and describes methods for overcoming them. Furthermore, we suggest an extension of the novel methodology for separation synthesis problems. 5.1. Validation against Flowsheet Simulator. To assess the accuracy of the approach against a commercial flowsheet simulator, we validated the results of temperature collocation against HYSYS. Figure 20 compares the column profiles obtained by the MIDI algorithm and with the commercial flowsheet simulator HYSYS. Even though we used identical vapor-liquid data sets, as expected, differences remain for two reasons: (a) The OCFE formulation solves continuous profile equations. It is well-known that this formulation


Figure 16. Minimum reflux ratio of an ideal mixture: (a) sloppy split, (b) sharp split.

Figure 17. Maximum reflux ratio of an ideal mixture: (a) pass the bottom production, (b) pass the distillation liquid composition.

derived from a first-order profile approximation deviates slightly from the tray-by-tray computations.5 (b) The OCFE model does not account for the effects of heat (i.e., the assumption of constant molar overflow). In light of the inherent model differences, the level of disagreement between the continuous collocation method and the trayby-tray model is judged to be acceptable for practical design purposes. More comparison and validation studies are recorded elsewhere.44 5.2. Tangent Pinch. Tangent pinches can occur in nonideal mixtures. In a separation problem with tangent pinches, the composition profiles exhibit additional stationary nodes, i.e., tangent pinch nodes. The existence of a tangent pinch can be detected by bifurcation analysis, a detailed discussion of which is beyond the scope of this paper.56,57 In a separation controlled by a tangent pinch, the composition profile terminates at the

tangent pinch point. Such a situation is depicted in Figure 21a, for which the MIDI algorithm accurately renders the shorter profile between the points xD and T. In near-tangent-pinch situations, as depicted in Figure 21b, the MIDI algorithm stably traverses the problematic column region. Standard tray-by-tray methods would fail because of the requirement for infinitely many trays, or they would simply produce inaccurate profiles because of the accumulation of round-off errors. The MIDI algorithm behaves well in actual tangent pinch control. In near-tangent-pinch situations, the MIDI algorithm passes the region with similar robustness as shown in the vicinity of saddle points. It should be noted that the feasibility test only aims to establish whether a design is realizable. Even when feasible, a column operation might be found to be uneconomical. Such questions require economic tradeoffs;

Ind. Eng. Chem. Res., Vol. 43, No. 12, 2004 3177 Table 4. Comparison of Minimum Reflux Ratios (rmin) and Relative Errors (er) Obtained by Different Methods split type

Underwood method rmin er (%)

zero-volume method rmin er (%)

RBM9 rmin

er (%)

rmin

BVM6 er (%)

MIDI rmin

er (%)

Mixturea

sloppy sharp

c

1.63 2.844

2.97 0.2

1.22 2.851

Constant Relative Volatility 27.4 1.22 27.4 0.03 2.851 0.03

sloppy sharp

0.99 1.54

25.6 1.9

sloppy sharp

0.64 0.43

31 2.3

1.68 2.85

1 0

1.663 2.85

0 0

25.6 1.9

1.33 1.54

2.2 1.9

1.36 1.51

0 0

Nonideal Mixturec 0.64 31 0.43 2.3

0.93 0.43

8.6 2.3

1.01 0.42

0 0

Ideal Mixtureb 0.99 1.54

a Constant-relative-volatility mixture’s specifications from Figure 13a and b. b Ideal mixture’s specifications from Figure 16a and b. Nonideal mixture’s specificatiosn from Figure 11.

Figure 18. Limiting refluxes for a direct split of a quaternary mixture.

Figure 19. Limiting refluxes for an indirect split of a nonideal quaternary mixture.

a simple example for the advantageous use of MIDI for a column design optimization problems is presented in the next subsection. 5.3. Column Sequencing. The speed of the algorithm holds high promise for its use in computer-aided

separation synthesis and optimization design problems. With rapid iterations and reliable detections of infeasibility, it is tempting to solve the problem of optimal column sequences for the achievement of desired products in closed form. The reduction in problem size


Figure 20. Composition profile for an OCFE model and flowsheet simulator (HYSYS).

afforded by temperature collocation on finite elements speeds convergence; its compactness and scaled dimensionless variable space enhance the stability of the optimization procedure.18 The performance of the new approach within the column optimization problem also benefits from the ready availability of sensitivity information offered by the BPD formulation. First-order sensitivity information requires only a simple differentiation of the OCFE model equations with respect to desired design variables, e.g., reflux. This is easily accomplished thanks to the compact expressions in terms of the orthogonal polynomials. Hence, the gradients required for the Hessians can be computed with little extra effort and superior precision over expensive and inaccurate numerical differentiations. In our preliminary studies of optimum column sequencing problems, we found a 40%

Figure 21. Composition profiles showing the presence of a tangent pinch.

reduction in CPU time when using sensitivity information rather than numerical differentiation. Even larger gains can be expected in more complex synthesis problems, i.e., more species and longer separation trains. MIDI plus sensitivity information also converged from initial values where regular SQP failed. Case Study E: Column Sequence for Minimum Operating Cost. We studied the optimal sequencing of a pentane-heptane-hexane separation task with the proposed feasibility test. The objective of this optimiza-


Figure 22. Direct and indirect sequences and profiles for the minimum-cost sequence for the desired product recovery of pentane, hexane, and heptane.

tion approach was to identify a minimum operating cost separation train for predefined purification targets (i.e., pentane in P1 ) 99.5%, heptane in P2 ) 96.0%, and hexane in P3 ) 84.5%) such as are often required for high-purity pharmaceutical products or solvent recovery tasks. A very simple performance model according to the marginal cost approach was adopted.58 Accordingly, operating cost were specified as functions of only the vapor rate. The optimum sequencing task consists of finding the configuration with the smallest total vapor rate, Q1 + Q2. Two configurations are possible, namely, a direct sequence and an indirect sequence, as depicted in Figure 22. The solution of the optimization problem outlined in eqs 23-27 is equivalent to finding the minimum reflux for each of the columns, where BPD1 and BPD2 are the bubble-point distances of the first and second columns, respectively, between the rectifying and stripping profiles at the corresponding temperatures Θ1 and Θ2.

min (Q1 + Q2)

Θ1,Θ2,r1,r2

(23)

s.t. ||BPD1|| e

(24)

||BPD2|| e

(25)

Θl1 e Θ1 e Θu1

(26)

Θl2 e Θ2 e Θu2

(27)

The cheapest sequence for the separation of a 100 kmol/h feed mixture employs a direct split in the first column for the removal of pentane with a minimum reflux ratio of r ) 1.53. In the second binary column, hexane and heptane are separated, leading to a total vapor rate of 63.71 kmol/h. The best indirect sequence requires a total vapor rate of 79.3 kmol/h. The optimal column configuration operates at minimum refluxes as can be inferred from Figure 22. Compared to the

marginal cost method, which is limited to constant relative volatilities, our method is able to produce results that extend to ideal and nonideal mixtures. Conclusions and Significance In this paper, we propose a novel fast and reliable algorithm for testing the feasibility of a distillative separation target. The nucleus of the new approach is the dimensionless bubble-point temperature as an independent variable, instead of the commonly used column heights or number of trays. The novel MIDI algorithm was shown herein to be effective for a wide variety of separation problems involving both sharp and sloppy specifications. The temperature collocation technique offers critical advantages over the regular tray-by-tray models, which can lock up in singular points (i.e., saddle and pinch points) requiring an infinite number of trays to be overcome. The temperature collocation method also dramatically reduces the problem size, thereby providing a computationally efficient feasibility test. Another advantage of temperature collocation lies in the improved numerical scale. Even independent variables are dimensionless, bounded within known limits. Better problem scaling was found to lead to superior numerical stability. The advantages of the MIDI algorithm were fully exploited in minimum- and maximum reflux ratio computations. Promising opportunities in column sequencing in which an inexpensive feasibility tests was needed in each iteration were briefly demonstrated. In the future, we will extend the new methodology to azeotropic mixtures and apply the MIDI algorithm to complex column configurations. We will explore the potential of the MIDI algorithm in conjunction with structural optimization techniques for the synthesis of distillative networks with desired performances. Acknowledgment Acknowledgment is made to the donors of the Petroleum Research Fund, administered by the ACS, for partial support of this research (ACS-PRF 35702-G9).


The authors gratefully acknowledge Professor Lorenz T. Biegler from Carnegie Mellon University for providing the source code of his rSQP algorithm.

C

{Kj[T,x(T)]xj}) ∑ j)1

d(

dx1

Nomenclature C ) number of components C1, C2 ) operating costs of cooling water BPD ) bubble-point distance K ) number of nodes in each element Ki ) equilibrium constant of component i n ) order of polynomial NB ) tray number of the rectifying profile NT ) tray number of the stripping profile NS ) number of elements in each profile psi ) vapor pressure of component i, kPa P ) total pressure, kPa p ) parameter Q1, Q2 ) vapor rates, kmol/h r ) reflux ratio s ) reboiler ratio T ) temperature, K xi ) liquid mole fraction of component i yi ) vapor mole fraction of component i

)0

(A-4)

Equation A-4 can be expanded further to yield C

∑ j)1 C

∑ j)1

{

{

}

dKj[T,x(T)] dxj xj + Kj )0 dx1 dx1

dKj[T,x(T)] dT dT

xj + Kj

dx1

}

dxj

(A-5)

)0

dx1

(A-6)

Rearrangement and collection of items in eq A-6 produces eq A-7 and an expression for dx1/dT in eq A-8. C

∑ j)1

{

dKj[T,x(T)] dT dT

dx1 C

Greek Letters Ri ) relative volatility of component i Θ ) dimensionless temperature, Θ ) (T - Tmin)/(Tmax Tmin) θ′ ) pseudo temperature γi ) liquid activity coefficient of component i ψi ) pseudo-vapor pressure of component i, kPa Subscripts B ) bottoms product c ) component index D ) distillation product i, j ) indices for nodes and elements p ) pinch point or fixed point

dx1

)-

∑ j)1

}

xj ) -

{

dxj

C

Kj ∑ dx j)1

dKj[T,x(T)] xj dT

dT

C

}

(A-8)

dxj

Kj ∑ dx j)1

(A-7)

1

1

The derivatives dxj/dx1 can be eliminated with the column height h to give the desired relationship for dx1/ dT as shown in eq A-9.

dx1

)

dT C

Superscripts l ) lower bound R ) rectifying section S ) stripping section u ) upper bound

-

∑ j)1

{

dKj[T,x(T)] xj dT

C

∑ j)1

Kj

}

C

)-

dxj dh

∑ dx1 j)1

}

[

dKj[T,x(T)] xj dT

dh

C

∑ j)1

dh dx1

Kj

dxj dh (A-9)

Appendix I Differential Equation Model with Temperature as the Independent Variable. The continuous differential equation models for both the rectifying section and the stripping section of the column with column height as the independent variable are given in eqs A-1 and A-2.5

dxi r+1 1 ) xi yi + xi,D dh r r

Inserting eq A-1 into A-9 provides the temperature dependence of the liquid composition profiles, xi, for all species, as desired.

dxi dT

(

) - xi -

∑

(A-1)

C

dxi s 1 ) y - xi + x dh s + 1 i s + 1 i,B

Equation A-3 is the bubble-point equation with equilibrium constants Kj, and compositions xj. C

{Kj[T,x(T)]xj} ) 1 ∑ j)1

∑ j)1

(A-2)

(A-3)

Total differentiation of eq A-3 with respect to x1 gives

)

r+1 1 yi + xi,D × r r C dK [T,x(T)] j xj dT j)1

[(

{

}

)]

r+1 1 xj yj + xj,D Kj r r

(A-10)

The derivation can be repeated in analogy for the stripping profile. In eq A-10, the differential of the equilibrium constant with respect to temperature, dKi/ dT is refined further as explained in Appendix II. Appendix II The temperature differential of the equilibrium constant, Ki, with respect to T includes the differential


pressure and differentials of the activity coefficients, as in eq A-12.

Ki )

psi γi[T,x(T)] P

(A-11)

(

)

dγi[T,x(T)] dKi 1 dpsi ) γ + psi dT P dT i dT

(A-12)

The differential of the liquid activity is a function of both the temperature and the composition for nonideal mixtures. The differential of the activity coefficient with respect to temperature follows from the chain rule.

dγi dT

C

)

∂γi dxj

∑ j)1 ∂x

+

dT

j

∂γi

(A-13)

∂T

Hence, the desired expression for Ki, dependent on both temperature, T, and composition, x, follows in eq A-14.

[

s 1 dpi ) γi + psi dT P dT

dKi

(

C

∂γi dxj

∑ j)1 ∂x

j

dT

+

)]

∂γi ∂T

(A-14)

Literature Cited (1) Widago, S.; Seider, W. Azetropic Distillation. AIChE J. 1996, 42, 96. (2) Wasylkiewicz, S. K.; Kobylka, L. C.; Castillo, F. J. L. Optimal Design of Complex Azetropic Distillation Columns. Chem. Eng. J. 2000, 79, 219. (3) Underwood, A. J. V. Fractional Distillation of Ternary Mixtures. Part I. J. Inst. Pet. 1946, 31, 111. (4) Underwood, A. J. V. Fractional Distillation of Ternary Mixtures. Part II. J. Inst. Pet. 1946, 31, 598. (5) Van Dongen, D. B.; Doherty, M. F. Design and Synthesis of Homogeneous Azeotropic Distillations. 1. Problem Formulation for a Single Column. Ind. Eng. Chem. Fundam. 1985, 24, 454. (6) Levy, S. G.; Van Dongen, D. B.; Doherty, M. F. Design and Synthesis of Homogeneous Azeotropic Distilations. 2. Minimum Reflux Calculations for Nonideal and Azeotripic Columns. Ind. Eng. Chem. Fundam. 1985, 24, 463. (7) Julka, V.; Doherty, M. F. Geometric Behavior and Minimum Flows for Nonideal Multicomponent Distillation. Chem. Eng. Sci. 1990, 45, 1801. (8) Julka, V.; Doherty, M. F. Geometric Nonlinear Analysis of Multicomponent Nonideal Distillation: Simple Computer-Aided Design Procedure. Chem. Eng. Sci. 1993, 48, 1367. (9) Bausa, J.; Watzdorf, R. V.; Marquardt, W. Shortcut Methods for Nonideal Multicomponent Distillation: 1. Simple Columns. AIChE J. 1998, 44, 2181. (10) Koehler, J.; Aguirre, P.; Blass, E. Minimum Reflux Calculations for Nonideal Mixtures Using the Reversible Distillation Model. Chem. Eng. Sci. 1991, 46, 3007. (11) Cho, Y. S.; Joseph, B. Reduced-Order Steady-State and Dynamic Models for Separation Processes I. Development of the Model Reduction Procedure. AIChE J. 1983, 29, 261. (12) Cho, Y. S.; Joseph, B. Reduced-Order Steady-State and Dynamic Models for Separation Processes II. Application to Nonlinear Multicomponent Systems. AIChE J. 1983, 29, 270. (13) Cho, Y. S.; Joseph, B. Reduced-Order Steady-State and Dynamic Models for Separation Columns III: Application to Columns with Multiple Feeds and Sidestreams. Comput. Chem. Eng. 1984, 8, 81. (14) Srivastava, R. K.; Joseph, B. Simulation of Packed-Bed Separation Processes Using Orthogonal Collocation. Comput. Chem. Eng. 1984, 8, 43. (15) Srivastava, R. K.; Joseph, B. Reduced Order Models for Staged Separation Columns. V. Selection of Collocation Points. Comput. Chem. Eng. 1985, 9, 601. (16) Srivastava, R. K.; Joseph, B. Reduced Order Models for Staged Separation Columns. IV. Treatment of Columns with

Multiple Feeds and Sidestreams via Spline Fitting. Comput. Chem. Eng. 1987, 11, 159. (17) Srivastava, R. K.; Joseph, B. Reduced Order Models for Staged Separation Columns. VI. Columns with Steep and Flat Composition Profiles. Comput. Chem. Eng. 1987, 11, 165. (18) Seferlis, P.; Hrymak, A. N. Optimization of Distillation Units Using Collocation Model. AIChE J. 1994, 40, 813. (19) Seferlis, P.; Hrymak, A. N. Adaptive Collocation of Finite Elements Models for the Optimization of Multi-Stage Distillation Unit. Chem. Eng. Sci. 1994, 49, 1369. (20) Stewart, W. E.; Levien, K. L.; Morari, M. Simulation of Fractionation by Orthogonal Collocation. Chem. Eng. Sci. 1985, 40, 409. (21) Swartz, C. L. E.; Stewart, W. E. A Collocation Approach to Distillation Column Design. AIChE J. 1986, 32, 1832. (22) Swartz, C. L. E.; Stewart, W. E. Finite-Element SteadyState Simulation of Multiphase Distillation. AIChE J. 1987, 33, 1977. (23) Cuthrell, J. E.; Biegler, L. T. On the Optimization of Differential, Algebraic Process Systems. AIChE J. 1987, 33, 1257. (24) Huss, R. S.; Westerberg, A. W. Collocation Methods for Distillation Design I: Model Description and Testing. Ind. Eng. Chem. Res. 1996, 35, 1603. (25) Huss, R. S.; Westerberg, A. W. Collocation Methods for Distillation Design II: Application for Distillation. Ind. Eng. Chem. Res. 1996, 35, 1611. (26) Chakraborty A.; Linninger, A. A. Plant-Wide Waste Management. 1. Synthesis and Multi-Objective Design. Ind. Eng. Chem. Res. 2002, 41, 4591. (27) Fidkowski, Z. T.; Doherty M. F.; Malone, M. F. Feasibility of Separations for Distillation of Nonideal Ternary Mixtures. AIChE J. 1993, 39, 1303. (28) Gani, R.; Bek-Petersen, E. Simple New Algorithm for Distillation Column Design. AIChE J. 2000, 46, 1271. (29) Doherty, M. F.; Malone, M. F. Conceptual Design of Distillation Systems; McGraw-Hill: New York, 2001. (30) Feehery, W. F. Dynamic Optimization with Path Constraints. Ph.D. Dissertation, MIT, Cambridge, MA, 1998. (31) Chowdhry, C.; Krendl, H.; Linninger, A. A. A Symbolic Numeric Index Analysis Algorithm for Differential Algebraic Equations Model Simplification for Dynamic Systems. Ind. Eng. Chem. Res., in press. (32) Biegler, L. T. Solution of Dynamic Optimization Problems by Successive Quadratic Programming and Orthogonal Collocation. Comput. Chem. Eng. 1984, 8, 243. (33) Cervantes, A.; Biegler, L. T. Large Scale DAE Optimization Using Simultaneous Nonlinear Programming Formulations. AIChE J. 1998, 44, 1038. (34) Renfro, J. G. Computational Studies in the Optimization of Systems Described by Differential Algebraic Equations. Ph.D. Dissertation, University of Houston, Houston, TX, 1986. (35) Villadsen, J. V.; Michelsen, M. L. Solution of Differential Equations by Polynomial Approximation; Prentice Hall: Upper Saddle River, NJ, 1978. (36) Fidkowski, Z. T.; Malone, M. F.; Doherty, M. F. Computing Azeotropes in Multicomponent Mixtures. Comput. Chem. Eng. 1993, 17, 1141. (37) Press, W. H.; Teukolsky, S. A.; Flannery, B. P. Numerical Recipes in C; Cambridge University Press: New York, 1995. (38) Linninger, A. A. Lecture Notes of Numerical Methods in Chemical Engineering; UIC-LPPD-090302; University of Illinois at Chicago: Chicago, IL, 2002. (39) Pernice, M.; Walker, H. F. NITSOL: A Newton Iterative Solver for Nonlinear Systems. SIAM J. Sci. Comput. 1998, 19, 302. (40) More, J. J.; Sorenson, D. C.; Garbow, B. S.; Hillstrom, K. E. The MINPACK Project. In Sources and Development of Mathematical Software; Wayne R. C., Ed.; Prentice Hall: Upper Saddle River, NJ, 1984; p 88. (41) Watson, L. T.; Billups, S. C.; Morgan, A. P. ALGORITHM 652: HOMPACK: A Suite of Codes for Globally Convergent Homotopy Algorithm. ACM TOMS 1987, 13, 281. (42) Schmid, C.; Biegler, L. T. Quadratic Programming Algorithms for Reduced Hessian SQP. Comput. Chem. Eng. 1994, 18, 817. (43) Biegler, L. T.; Grossmann I. E.; Westerberg A. W. Systematic Methods of Chemical Process Design; Prentice Hall: Upper Saddle River, NJ, 1997.

3182 Ind. Eng. Chem. Res., Vol. 43, No. 12, 2004 (44) Mark, N. P.; Zhang, L.; Linninger, A. A. Robust and Reliable Column Design Method; REU Report UIC-LPPD-073003; University of Illinois at Chicago: Chicago, IL, 2003. (45) Linninger, A. A.; Zhang, L. The Solver Manual; UIC-LPPD031003; University of Illinois at Chicago: Chicago, IL, 2003. (46) Tanartkit, P.; Biegler, L. T. Stable Decomposition for Dynamic Optimization. Ind. Eng. Chem. Res. 1995, 34, 1253. (47) Gmehling J.; Onken, U. Vapor-Liquid Equilibrium Data Collection; DECHEMA Chemistry Data Series; DECHEMA: Frankfurt, Germany, 1977. (48) Gregor, H.; Linninger, A. A. TD2000 Property Database; UIC-LPPD-081901; University of Illinois at Chicago: Chicago, IL, 2001. (49) Wahnschaft, O.; Koehler, J.; Blass, E.; Westerberg, A. The Product Composition Region of Single-Feed Azeotropic Distillation Columns. Ind. Eng. Chem. Res. 1992, 31, 2345. (50) Vassiliadis, V. Computational Solution of Dynamic Optimization Problems with General Differential-Algebraic Constraints. Ph.D. Dissertation, Imperial College, London, 1993. (51) McCalbe, W. L.; Thiele, E. W. Graphical Design of Fractionating Columns. Ind. Eng. Chem. 1925, 17, 605. (52) Ponchon, M. Graphical Study of Distillation. Technol. Mod. 1921, 13, 20.

(53) Reyes, J. A.; Gomez, A.; Marcilla, A. Graphical Concepts to Orient the Minimum Reflux Ratio Calculation on Ternary Mixtures Distillation. Ind. Eng. Chem. Res. 2000, 39, 3912. (54) Colburn, A. P. The Calculation of Minimum Reflux Ratio in the Distillation of Multicomponent Mixtures. Trans. Am. Inst. Chem. Eng. 1941, 37, 805. (55) Gilliland, E. R. Multicomponent Rectification: Estimate of the Number of Theoretical Plates as a Function of the Reflux Ratio. Ind. Eng. Chem. 1940, 32, 1220. (56) Levy, S. G.; Doherty, M. F. A Simple Exact Method for Calculating Tangent Pinch Points in Multicomponent Nonideal Mixtures by Bifurcation Theory. Chem. Eng. Sci. 1986, 41, 3155. (57) Fidkowski, Z. T.; Doherty M. F.; Malone, M. F. Nonideal Multicomponent Distillation: Use of Bifurcation Theory for Design. AIChE J. 1991, 37, 1761. (58) Modi, A. K.; Westerberg, A. W. Distillation Column Sequencing Using Marginal Price. Ind. Eng. Chem. Eng. 1992, 31, 839.

Received for review October 30, 2003 Revised manuscript received March 18, 2004 Accepted March 22, 2004 IE034223K

Temperature Collocation Algorithm for Fast and Robust Distillation

Recommend Documents