Ind. Eng. Chem. Res. 2008, 47, 51-65
51
PROCESS DESIGN AND CONTROL Process Design Approach for Reactive Distillation Based on Economics, Exergy, and Responsiveness Optimization Cristhian P. Almeida-Rivera* Process Science Department, UnileVer Food and Health Research Institute, OliVier Van Noortlaan 120, 3130 AC Vlaardingen, The Netherlands
Johan Grievink Process Systems Engineering Group, Delft UniVersity of Technology, Julianalaan 136, 2628 BL Delft, The Netherlands
The interactions between economic performance, thermodynamic efficiency, and responsiveness in reactive distillation process design are explored in this contribution. This motivation is derived from taking a sustainable life span perspective, where economics and avoidance of potential losses of resources (i.e., mass, energy, and exergy) in process operation over the process life span are taken into consideration. The approach to reactive distillation column design involves defining a generic lumped reactive distillation volume element and the development of rigorous dynamic models for the behavior of the element. As an extension of conventional existing models, this model, which is used as a standard building block, includes entropy and exergy computations. Three objective functions are formulated, which account for the process performance regarding economy, exergy loss, and responsiveness. A fundamental understanding of the strengths and shortcomings of this approach is developed by addressing various case studies: (i) steady-state simulation of a classical MTBE design based on economics only, useful as a reference case; the multiobjective optimization of a MTBE reactive distillation column with respect to (ii) economics and thermodynamic efficiency; and (iii) economics, thermodynamic efficiency, and responsiveness. Structural differences and dynamic responses are identified between the optimized and classical designs to stress the importance of considering exergy and responsiveness criteria. Thus, the classical (economics-driven) design is compared to a green (economicsand exergy-driven bioptimized) design. The green design produces less entropy than the classical one, while being marginally less attractive from an economic standpoint. The bioptimized design has an improved closedloop performance compared to the classical design, keeping the same control structure and settings. It is concluded that incorporating both economic- and exergy-related objectives in the unit design results in a process with better closed-loop performance and reduced exergy loss. Introduction Advances in computing and information technology allow chemical engineers to solve complex design problems arising from a need to improve sustainability of the biosphere and human society. In such a context, the performance of a chemical plant needs to be optimized over its manufacturing life span, while accounting for the use of multiple resources in the design and manufacturing stages. Since such resources are of a different nature (e.g., capital, raw materials, and labor) with different degrees of depletion and rates of replenishment, the performance of a chemical plant is characterized by multiple objectives with tradeoffs.1 An integrated unit operation like reactive distillation (RD) can offer significant benefits for resource utilization, when it is accompanied with the avoidance of chemical equilibrium limitations, enhancement of conversion, selectivity, and yield, and the reduction of both operational and capital costs. The nonlinear coupling of chemical reaction, phase equilibria, and * To whom correspondence should be addressed. Tel: +31 10 460 6504. Fax: +31 10 460 5025. E-mail: Cristhian.Almeida-Rivera@ Unilever.com.
transport phenomena can, however, lead to undesired, systemdependent features, such as multiple steady states2-17 and reactive azeotropes.18-25 The results of a case study on the interactions between economics, exergy efficiency, and process responsiveness in the conceptual design of RD processes are discussed in this contribution. In the wider context of conceptual process design, sustainability issues like feedstock selection, (re)use of catalyst, and solvents are of key importance. In our perspective of life span, we do not cover those aspects as we focus exclusively on the minimization of loss of resources in the operational phase. Currently, conceptual design in chemical engineering is predominantly driven by (steady-state) economic considerations, postponing to latter engineering stages any issue related to exergy analysis, responsiveness, and control performance. The importance of incorporating exergy efficiency in the early stages of the design cycle is stressed by the increasing depletion of nonrenewable energy resources and the relatively low second law efficiencies (5-20%) for the majority of chemical processes. In spite of these facts, exergy analysis is in most of the cases performed as a follow-up task to process synthesis.26-29 In
10.1021/ie061521j CCC: $40.75 © 2008 American Chemical Society Published on Web 11/30/2007
52 Ind. Eng. Chem. Res., Vol. 47, No. 1, 2008
addition to economic and exergy considerations,30 the performance optimization of a chemical plant should include controllability and process responsiveness. These aspects are important to maintain product quality and to avoid excessive use of utilities for control actions and reprocessing of off-spec material. The focus of this contribution is the importance of including interactions between economics, exergy efficiency, and process responsiveness in RD conceptual design. Aiming at these goals, a generic compartment model is defined, which accounts for mass, energy, and momentum balances and includes entropy and exergy computations. This behavioral model is integrated to create a RD column structure. The design is cast into a nonlinear optimization problem, where degrees of freedom in the design model are fixed by optimization of multiple objectives and applying constraints, reflecting heuristics and engineering knowledge. At this column level, three performance criteria are defined for the optimization of the unit design. These embrace the unit performance regarding economics, thermodynamic efficiency, and process responsiveness. The economic performance index is given by the total annualized cost, whereas the thermodynamic efficiency index accounts for the process irreversibilities and is expressed in terms of the entropy production rate by chemical reaction, heat, and mass transfer. The responsiveness index is similar to the response time defined by Luyben et al.31 The three objective functions are solved simultaneously within a single multiobjective optimization problem. At this level of understanding, it is believed that a multiobjective optimization gives more insight into trade-offs between objectives than combining the three objectives into a single monetary objective function using ad hoc weighting factors. Comparing the design based on the 3-fold performance criteria with a reference economic-driven design allows us to identify structural differences. The response of the optimized design in the presence of a deterministic disturbance scenario is validated using dynamic simulation. The thrust of this contribution is the optimization-based design approach. Multiple objective functions are considered and embedded in a large-scale optimization problem. Generic Lumped Reactive Distillation Volume Element The starting point for a fundamental understanding of the strengths and shortcomings of the proposed multiobjective design approach is the definition of a generic lumped RD volume element GLRDVE. This element (or compartment) is chosen as a building block in the design approach and selected to represent the governing phenomena occurring in RD processing. This behavioral model accounts for the following: (i) mass, energy, and momentum balances, together with entropy and exergy computations; (ii) thermodynamic equations of state for the specific heat and specific density and for the chemical potential per species; (iii) thermodynamic phase equilibrium conditions at the interphase; (iv) the constitutive rate equations for chemical reactions; (v) the rate equations for interphase transport of species mass and energy; (vi) the rate equation for heat transfer with an external solid surface; (vii) the hydrodynamics-based rate equations for convective mass and energy flows with the external world and the pressure drop; and (viii) the equations for the interphase contact surfaces and for the void fraction of the vapor phase in the liquid phase. The following attributes are introduced for the system proposed in Figure 1: (i) two internally (well-mixed) homogeneous fluid phases; and (ii) chemical reactions in the liquid phase only. Thus, mass transfer between phases and heat exchange are allowed in the compartment, and chemical reaction
Figure 1. Schematic representation of the GLRDVE. Feed and side drawoff streams and heat exchange are allowed.
is triggered by the presence of an appropriate catalyst. Each GLRDVE has its own complexities due to the many interacting phenomena. However, this does not account for the interactions between the stages in a real multistage column. The latter class of interactions could give rise to nonlinear dynamic behavior and poor controllability if an unfortunate combination of reflux ratio and number of (non)reactive trays is chosen. Issues related to connecting compartments and superstructure description are addressed in Integration of Volume Elements to a Column Structure, Definition of Optimization Criteria, and Optimization Approach. Mass/Energy/Entropy/Momentum Balances. Component, energy, and entropy balances are developed in the GLRDVE. The entropy and, in particular, the entropy production rate are of paramount importance to determine exergy losses and the responsiveness of the process. The balances lead to the following set of differential algebraic equations, (R) ∂(F(R) i ×V )
j)nst
)
∂t
Fj,(R) - Ψi × aΨ + δ(R) × ri × VLRX ∑ i j)1 i ∈ Znc (1)
∂U(R) ∂t
j)nst i)nc
)
Fj,(R) × hj,(R) + Ψe × aΨ +.... + ∑ ∑ i i j)1 i)1 j)nst i)nc
Fj,(R) × gj,(R) + δ(R) × Q(R) ∑ ∑ i i j)1 i)1
(2)
where Q(R) is assumed to only occur in the liquid phase (i.e., Q(R) ) QL). The Kronecker delta function defined as
δ(R) :)
{
0 (if (R) * L) 1 (if (R) ) L)
(3)
The energy flux associated with the interphase mass and heat transfer Ψe is composed of a measurable heat flux and the partial molar enthalpies brought by the component fluxes. Namely, i)nc
Ψe ) -
(R′-R) Ψi × h(R) ∑ i +q i)1
(4)
The entropy balance in the transient state at a macroscopic level is given by the net amount of entropy transferred into the volume element per unit time and the amount of entropy used/ produced in the control volume per unit time. Thus, the entropy balance for a given phase (R) can be stated as,
∂S(R) ∂t
j)nst i)nc
)
× Sj,(R) + Σ(R) ∑ ∑ Fj,(R) i i V j)1 i)1
(5)
where Σ(R) V includes exclusively interfacial diffusional and chemical reaction contributions.
Ind. Eng. Chem. Res., Vol. 47, No. 1, 2008 53
Alternatively, the entropy balance can be expressed in terms of the entropy density s(R) and the fluxes of entropy density,32,33
∂s(R) (R) ) -∇J(R) s + σV ∂t
(6)
For stream j, the flux of entropy J js is related to the total flux j according to the following of energy J je and mass fluxes J m,i 32 expression,
) Jj,(R) s
1 T
i)nc
∑ i)1
× Jj,(R) e
µj,(R) i T
j,(R) × Jm,i
j ∈ Znst
(7)
The fluxes of mass and energy required in expression 7 can be obtained from the equation of continuity per volume unit and energy balance, respectively, for an open unsteady-state multicomponent system.34 The entropy production rate at macroscopic level Σ(R) V (eq 5) can be obtained after integration of the local entropy production rate density σ(R) V over the volume element in phase (R). For the sake of simplification, the following assumptions are adopted: (i) entropy contributions of all involving phases in the GLRDVE are combined in a single (R) variable (Σ(*) V ) f (ΣV )); (ii) entropy is produced due to mass and heat diffusion in the liquid-vapor interphase, chemical reaction in the bulk of the liquid, and heat transfer from external heat sources/sinks and the surrounding environment to the liquid phase only; and (iii) the effects of electric field, external forces, and diffusive momentum fluxes are neglected. The following expression for the integral entropy production term can be then derived,32,35,36
∫ ∫ ∫V σ
(*)
dV )
Σ(*) V
1
)
L
T
ALi × Li × VLRX + JLq × ∑ i)1
Jm,i × aΨ × hi × Xe + ∑ i)1
i)nc
(-
Jm,i × aΨ × Xm,i) + QL × XHX ∑ i)1
(8)
The entropy contribution due to chemical reaction contains the affinity of chemical reaction and its degree of advancement, which are defined respectively as, i)nc
ALi ) -
∑ i)1
νLi,j × µLi
j ∈ Znrx
(9)
j)nrx
rLi )
parameter
driving force
mass transfer interphase energy transfer chemical reaction external heat exchange
in Xm,i ) -(1/T)∆µi,T ) R × ln(yin i × P /yi × P) Xe ) ∆(1/TV) ) 1/2((1/Tout,V) - (1/Tin,V)) Xrx,j ) (ALj /T) XHX ) (1/T) - (2/(Tu,in + Tu,out))
a
i ∈ Znc; j ∈ Znrx .
For component i the driving force Xm,i is related to the chemical potential gradient between the incoming and outgoing vapor streams.37,38 The integrated flux associated with Xm,i is calculated by mass balance over the vapor phase,28,37 j)nst
Ji × aΨ )
j Fj,V ∑ total × yi j)1
i ∈ Znc
(11)
This integrated flux expression is consistent with the component balance given by eq 1 if Ψi equals Jm,i. Moreover, the total streamflow is related to the component flow according to, j j,V Fj,V total × yi ) Fi
i ∈ Z nc
(12)
The driving force Xe is estimated as suggested by Koeijer and Rivero.28 The measurable heat flux through the interphase associated with Xq is derived from the energy balance in the vapor side,37 j)nst
Jq × aΨ )
j,V Fj,V ∑ total × hi j)1
(13)
i)nrx
i)nc
a Ψ × Xq +
Table 1. Expressions To Estimate Driving Forces in the GLRDVE of Figure 1a
νLi,j × Lj ∑ j)1
i ∈ Z nc
(10)
To estimate the driving forces and fluxes occurring in the volume element, additional assumptions are introduced: (i) the GLRDVE is composed of incoming and outgoing streams of liquid and vapor phases; (ii) the vapor phase behaves ideally; (iii) the flux of mass and heat through the interphase is positive from the liquid side to the vapor side; (iv) dynamic changes in gas phase hold-ups of species mass and energy are ignored; and (v) Gibbs free energies at liquid and vapor sides are equal. The underlying reasoning behind the inclusion of these assumptions is a more transparent derivation of the driving forces and flux expressions, but without compromising the validity of the approach. The expressions used to estimate the averaged driving forces are summarized in Table 1.
In this expression, the enthalpy reference state is calculated at the average temperature at which components are transferred through the interphase. The expression suggested by Koeijer et al.39 and Koeijer37 has been used to estimate the XHX average driving force. In the development of the GLRDVE, a reduced form of the momentum balance is considered. We assume a pseudo steady state and negligible viscous contribution. Moreover, the gradient of pressure is assumed to be much larger than the rate of momentum gain by convention of the fluid. These simplifications applied to the momentum balance expression given elsewhere,34 result in the expression,
∇P(R) ) F(R) ‚g(R)
(14)
Availability and Availability Function. Two additional variables of key importance for this thermodynamic-based approach are the availability A and the availability function B. In this context, availability is regarded as a measure of the maximum amount of useful energy that can be extracted when the system is brought to equilibrium with the dead-state surroundings.35 Thus, the availability of the system in phase (R) and state (T1,P1) is given by the expression,
A(R) ) (h(R) - T0 × S(R))T1,P1 - (h(R) - T0 × S(R))T0,P0 (15) where the dead-state conditions are those at which T0 ) 298.15 K and P0 ) 1 atm. The availability function B at state given by (T,P) is computed according to the following equality, (R) ) (h(R) - T0 × S(R))T,P BT,P
(16)
54 Ind. Eng. Chem. Res., Vol. 47, No. 1, 2008
Contrary to the availability, the availability function of the system at a given state (T1,P1) can be referenced to any other state (T2,P2),
∆B
(R)
) (h
(R)
- T0 × S )T1,P1 - (h (R)
(R)
- T0 × S )T2,P2 (17) (R)
According to their definitions, A and B increase with increasing enthalpy and decrease with increasing entropy. Moreover, the numerical values of their variations coincide (i.e., ∆B(R) ) ∆A(R)). Based on the combination of energy and entropy balances and according to the definition of availability an exergy balance can be obtained for a given phase (R) in the absence of external forces (see Appendix A in Luyben et al.31),
∂B(R) × M(R)
j)nst
)
∂t
∑ j)1
Fj,(R) total
×B
j,(R)
+
( ) 1-
T0
TQ
Q(R) - T0 × Σ(R) V (18)
Design Decision Variables. According to a degree of freedom analysis performed in the GLRDVE the following design variables are to be specified for a given operating conditions and feed specifications: (a) reaction volume V(R) RX
V(R) RX
)
{
V(R) mcat × Fm,cat
-1
phase homogeneously (for liquid ) catalyzed reaction(s) for heterogeneously catalyzed ( ) reaction(s)
(b) external and heating/cooling flow rate Q(R), and (c) the interfacial area for heat and mass transfer aΨ. Connecting streams and associated flow rates Fj,(R) are i assumed to be given. Operational variables such as reflux ratio and pressure are of key relevance in RD and should be chosen with thoroughness. In our approach, these variables are specified at a higher decision making layer (cf. Volume Elements to a Column Structure, Definition of Optimization Criteria, and Optimization Approach), aiming at operating windows where steady-state multiplicity and reactive azeotropy phenomena are not likely to occur. The perspective adopted for selecting the design decision variables is task-oriented rather than the conventionally used equipment-oriented design approach. First, the task-based concept would allow us to derive the phase volume V(R) from the reaction time (i.e., residence time in element) to meet a specified degree of advancement of reaction; second, the amount of heat to be added or removed allows one to compute the heat exchange surface and coolant flows; and third, the interfacial areas for species mass and heat transfer may be translated into the mechanical design of a tray and its prevailing hydrodynamic regime such that the required areas can be met within constraints. The main benefit of this functional, task-oriented perspective involves the process design unconstrained by equipment geometry. Thus, task-oriented design prevails against equipmentoriented design. Economics/Exergy Efficiency and Responsiveness Objective Functions. A set of three performance criteria indicators is defined for the GLRDVE: the economic performance index fecon, the sustainability performance index fexergy, and the responsiveness performance index fresponse. The economic performance indicator fecon represents the total annualized cost TAC of the GLRDVE, including capital
investment and operational costs of the process, j)nsti)nc
fecon ≡ TAC )
Fj,(R) ci + Fcuccu + Fhuchu + ∑ ∑ i j)1 i)1 cunit+catalyst
R ∈ Znp (19)
The exergy efficiency fexergy is accounting for the process irreversibilities and is expressed in terms of efficiency of the second law of thermodynamics, ηII,
fexergy ≡ ηII )
Wideal Wideal - T0ΣV
(20)
where Wideal (Wideal < 0) is defined as by Seider et al.35 and fixed for a given process.37 From a sustainability perspective, the desired goal is the enhancement of ηII. According to eq 20 the maximization of fexergy is equivalent to the minimization of ΣV.37,40 In this contribution, we adopt this second approach and compute T0ΣV for a fixed T0 value. The responsiveness index fresponse is approximated to the settling time defined by Luyben et al.,31
fresponse ≡ τ ) max x
∆B(x) T0∆ΣV(x)
(21)
At steady-state conditions B(x) is at its minimum value and increases when the system is perturbed by a disturbance. Since a similar situation holds for ΣV(x), the response time τ measures how fast a small variation ∆B(x) would dissipate due to an increased ∆ΣV(x).31 In the previous definition of fresponse, Luyben et al.31 assumed that both the availability and entropy production rate are exclusively functions of the vector of state variables x. Although this assumption is fundamentally correct, it somehow disguises the way of estimating the effect of a small disturbance or a change in an input on the system. Therefore, we propose a subtle modification of the arguments of B and ΣV, by which a second argument is included, representing the systemic view on a dynamic system. Thus, the availability and entropy production rate are regarded as functions of an internal state vector x and a vector of disturbances or manipulated external excitations u,
exergy: B(x,u)
(22a)
entropy production rate: ΣV(x,u)
(22b)
Combining the definitions of eqs 22 with eq 21 leads to a more generalized expression for fresponse,
fresponse ≡ τ )
∆B(x,u) T0∆ΣV(x,u)
(23)
The previous expression could be further elaborated by defining the values of state variables x and disturbances u at each state. Hence,
with
∆B(x,u) ) B(x,u) - B(x0,u0)
(24)
∆ΣV(x,u) ) ΣV(x,u) - ΣV(x0,u0)
(25)
Ind. Eng. Chem. Res., Vol. 47, No. 1, 2008 55
u ) u0 + ∆u
(26)
x ) x(u)
(27)
ΣV(x,u) g ΣV(x0,u0)
(28)
An additional feature of Luyben’s original definition of fresponse, not shown in expression 23, is related to the maximization for the time constant τ over all possible steady states surrounding the reference steady state. In this way, one looks for the worst case over all possible small-scale disturbances and input change scenarios. This optimization-driven approach justifies Luyben’s approach of including a single argument for ∆B and ∆ΣV. In our approach, such an maximization approach is not considered due to its inherent high computational effort. However, a rather formal definition of fresponse is adopted (eq 23) together with a single input excitation. By not searching all possible state trajectories, the generality of the results might be questionable and the external excitation is a fairly conservative one. However, engineering judgment can be exercised in selecting the disturbance variable u and the magnitude of the disturbance ∆u. Integration of Volume Elements to a Column Structure, Definition of Optimization Criteria, and Optimization Approach Creating a Column Structure from the Generic Volume Element. The integration of compartments to create a column structure introduces additional degrees of freedom, which should be conveniently specified (e.g., reflux ratio). In the more general context of process synthesis, the issue of coupling compartments to a process unit has been frequently addressed41-43 and is not, therefore, explicitly considered in this research. Briefly speaking, the compartment integration is carried out by introducing a superstructure and finding the optimal configuration using hybrid optimization techniques over the superstructure. The description of the superstructure and the rules for connecting compartments with each other and with the outer world are embedded in the structure model. Extensive information on structure models can be found in the research work of, for example, Ismail et al.41 and Papalexandri and Pistikopoulos.42 The rules of connectivity used in the RD superstructure state the following: (i) the outgoing liquid at a given compartment enters the compartment located below; (ii) the outgoing gas at a given compartment enters the compartment located above; (iii) the liquid stream leaving the condenser is partially recycled back to the top compartment and withdrawn from the unit as top product; (iv) the liquid leaving the bottom compartment enters the reboiler; (v) the streams leaving the reboiler are a vapor stream (entering the bottom compartment) and a liquid stream (withdrawn from the unit as a bottom product); and (vi) side feed streams and side withdrawal streams may occur. The connectivity rules, however, impose some configuration limitations: (i) heat integration over volume elements is not considered; (ii) side reactors are not present; and (iii) no parallel structures are considered. During the development of the RD unit, the following set of simplifications are made: (i) the liquid and vapor streams leaving a tray are in thermodynamic equilibrium; (ii) both liquid and vapor hold-ups are fully considered; (iii) the liquid and vapor on each tray are perfectly well mixed (i.e., no radial gradient in temperature and composition and the liquid on the tray has a concentration equal to that of the liquid leaving the tray); (iv) the reactions are taking place in the liquid phase and
enhanced by a solid catalyst; (v) the reactions occur only on the trays containing catalyst; (vi) vapor and liquid entrainment is neglected; (vii) the column is equipped with a condenser, reflux drum, and reboiler; (viii) the location(s) of the feed stream(s) and product stream(s) are known a priori; (ix) the liquid hold-up on a tray is related to the liquid outflow according to Francis weir formula; (x) the flow of vapor leaving a tray is determined by the pressure drop between two consecutive trays; (xi) the tray pressure drop includes contributions of dry pressure drop, liquid height, and residual head; (xii) the physical properties of the involved streams are computed as a function of temperature or phase composition; (xiii) the heat of reaction is omitted in the energy balance, as the heat of formation at standard conditions (1.033 × 105 Pa/298 K) is used as reference in the enthalpy computations; (xiv) the reboiler, condenser, and reflux drum are nonreactive units; (xv) the tray sizing is estimated by rules-of-thumb criteria;44 and (xvi) the column diameter is larger than the minimum flooding value. The compartments representing the behavior of the condenser and reboiler are, in principle, identical to the GLRDVE, provided the following characteristics are present: (i) no chemical reaction takes place due to the absence of catalyst; (ii) heat exchange takes place in both phases; (iii) the condenser is a total condenser; and (iv) the reboiler is a partial reboiler. In addition to the specification of the reflux ratio, the compartment integration to the superstructure involves the following features of the structure model: (i) all compartments can have a different catalyst load; and (ii) the column sections can have different diameters. Development of Three Criteria for Optimization of Column Design. Stepping out from the GLRDVE case to the integrated RD unit is intuitive for fecon and fexergy. The variables ∆B(x,u) and ∆ΣV(x,u), on the other hand, need to be globally evaluated over all the structure for the responsiveness criterion fresponse. The resulting performance criteria are given by the following set of expressions for the kth compartment, k)nt j)nst i)nc
fecon )
k k k [Fj,k ∑ ∑ ∑ i ci + Fcuccu + Fhuchu + cunit+catalyst] k)1 j)1 i)1
(29a)
fexergy )
(Wideal)
total
(Wideal)total - T0(ΣV)total
fresponse )
(∆B(x,u))total T0(∆Σ V(x,u))total
(29b)
(29c)
The validity of an additivity property (i.e., the objective functions per compartment can be added up to a corresponding objective function over the entire structure) is implicitly assumed in these definitions. An indication of the implications of process responsiveness for design is revealed by the definition of fresponse. According to eq 30, a short responsiveness time is desired and can be attained by a large change in the entropy production rate for a given change in stored exergy. Using an ideal heat exchanger as a case study, for instance, Almeida-Rivera45 and Meeuse46 demonstrated that better process responsiveness is attained at smaller driving forces and at small hold-up values. Moreover, Almeida-Rivera45 suggested that as a heuristic for heat-transfer processes the exchange surface per unit mass should be maximized, while driving force is kept small at constant transfer duty.
56 Ind. Eng. Chem. Res., Vol. 47, No. 1, 2008 Table 2. Utopia Point in Optimization Problems with More Than One Objective Function Utopia point U in a RΦ-space is given by the Φ-objective function values that solve independently the Φ scalar optimization problems.62 For instance, in an optimization problem involving two objective functions the mathematical formulation of the utopia point U is min f1(x) ) flow 1 (x*) x∈Γ
min f2(x) ) x∈Γ
flow 2 (x*)
and
f2(x*) ) fup 2 (x*)
(43)
and
f1(x*) )
(44)
fup 1 (x*)
f low i (x*)
and f up i (x*) are the lower and upper limits of function fi, respectively, and the set of constraints
where the objective defines the space of design variables x. The coordinates of the utopia point are then low given by (f low 1 , f2 )
Approach To Design Optimization and References Cases. Two approaches are adopted to solve the optimization problem involving more than one objective function, (a) Each objective function is used independently, disregarding the other two; this approach localizes the coordinates of the utopia point (cf. Table 2) and, subsequently, reduces the computational effort by limiting the search space, (b) The interactions and trade-offs among the objectives are accounted; in this approach, a multiobjective problem is formulated by using a weighted p-norm for the vector of objective functions, as given by Clark and Westerberg,47
min F ≡ [ D∈Γ
∑i (wifi)p]1/p
p ∈ [1, ∞)
(1) Steady-State Entropy Production Profile in a MTBE Reactive Distillation Column. The esterification of isobutene and methanol is used as a case study in this application. MTBE is produced from methanol (MeOH) and isobutene (iC4) in an equilibrium-limited, liquid-phase reaction48,49 catalyzed heterogeneously by bentonites, zeolites, or strong acidic macroporous ion-exchange resins or homogeneously by sulfuric acid according to the following stoichiometry,
C4H8 + CH3OH h C5H12O isobutene methanol MTBE
(31)
FCC is a commonly used feedstock for the MTBE production.49 This feedstock stream contains a complete range of butanes and butenes. For the sake of simplicity and provided that the reaction is highly selective for iC4, all the organic components in the iC4 stream are lumped together as inert nC4.5 The chemical reaction is represented by a pseudohomogeneous kinetics as reported in the literature,17,50,51
(
aMeOH
-
iC4 MeOH MTBE nC4
aMTBE
)
Keq × aMeOH2
(32)
In this study, we use the γ - φ thermodynamic formulation.52 The vapor phase is assumed to behave ideally, whereas the liquid phase is modeled using the Wilson activity coefficient equation.53 The interaction parameters for the MTBE system at a given operation condition are listed in Table 3. The reaction equilibrium and kinetics of MTBE synthesis have been the subject of intensive research during the past
MeOH
0.00000 0.74200 -0.2413 0.00000
iC4 MeOH MTBE nC4
MTBE
aij -0.74200 0.00000 -0.98330 -0.81492 bij (K) -85.5447 0.0000 204.5029 -192.4019
0.0000 -1296.719 -136.6574 0.0000
nC4
0.24130 0.98330 0.0000 0.00000 30.2477 -746.3971 0.0000 0.0000
0.00000 0.81492 0.00000 0.00000 0.0000 -1149.280 0.0000 0.0000
j)nc a ln(Λ ) ) a + b /T; S ij ij ij W,i ) ∑j)1 xj × Λij; ln(γi) ) 1 - ln(SW,i) k)nc ∑k)1 (xk × Λki/SW,k).
Table 4. Equilibrium Constant and Kinetic Constant for the Synthesis of MTBE58
(
ln Keq ) ln Keq,0 + Rk
)
()
1 1 T + βk × ln + γk × (T - T0) + δk × T T0 T0
(T2 - T02) + k × (T3 - T03) + φk × (T4 - T04) parameters
validity
101;
[ (
validity
(45)
Rk ) -1.49277 × βk ) -7.74002 × γk ) 5.07563 × 10-1; δ ) -9.12739 × 10-4; -6 k ) 1.10649 × 10 ; φk ) -6.27996 × 10-10; T0 ) 298.15 K; Keq,0 ) 284 heterogeneous catalyst (Amberlyst 15) 103;
kf(T) ) kf(T*) × exp parameters
Applications
aiC4
iC4
(30)
where a norm of p e 2 is used to localize noninferior points in the nonconvex regions.47 This theory will now be applied to the multiobjective design cases for a methyl tert-butyl ether (MTBE) column.
grx(x,T) )
Table 3. Wilson Interaction Parametersa for the System iC4-MeOH-MTBE-nC4 at 11 × 105 Pa5,74
Ea 1 1 R T T*
)]
(46)
kf (T*) ) 15.5 × 10-3 mol × (s × eq[H+])-1; Ea ) 92400 J × mol-1; T* ) 333.15 K heterogeneous catalyst (Amberlyst 15)
Table 5. Nominal Values in the MTBE Synthesis Nominal Values in the MTBE Synthesisa parameter
value
number of trays feed 1 location flow rate composition temperature feed 2 location flow rate composition temperature operating pressure reflux ratio mass of catalyst/tray extent of reaction column diameter condenser area reboiler area
15 [2 rectifying, 8 reactive, 5 stripping] iC4 + nC4, 10th tray (g) 455 mol s-1 iC4, 0.37; nC4, 0.63 350 K MeOH, 9th tray (l) 168 mol s-1 MeOH, 1.00 320 K 11 × 105 Pa (top) 7 204.1 kg ) 1000 eq[H+] 150 mol s-1 3.573 m 4.460 × 101 m2 2.456 × 102 m2
a Operational conditions are those reported in Hauan et al.,8 Jacobs and Krishna,4 and Seader and Henley.59
decade.54-58 In this study, we use the model reported by Thiel et al.58 and listed in Table 4. The GLRDVEs are integrated to a column structure in such a way that the superstructure is comparable to the unit studied elsewhere4,8,59 and specified in Table 5. In this column structure, the catalyst load is evenly distributed along the reactive trays4,8,59 and all the sections of the column have a common value. The structure to be further considered during the course of this contribution is depicted schematically in Figure 2. To prepare for the optimization task and to reduce the problem complexity, we identify the dominant contributions to entropy production rate over the MTBE column. The expressions to compute all contributions in the GLRDVE are the terms of eq
Ind. Eng. Chem. Res., Vol. 47, No. 1, 2008 57 Table 6. Optimized Design of a RD Column for MTBE Synthesis Based on Economic Performance
Figure 2. Schematic representation of a RD column as the integration of GLRDVEs. Numbers represent tray/compartment location.
8. These expressions are applicable for all reactive and nonreactive compartments. The reboiler and condenser models assume that the entropy contribution by external heat transfer is dominant. The reference design (E-design) is the column structure obtained by steady-state optimization of the spatial design variables with respect to the unit’s economic performance. For this optimization problem, the following formulation is considered,
Function: s.t.
min
Dcol,Areb,Acond∈Γ
fecon ) TAC
(33)
gx(‚) ) 0 hx(‚) e 0 D D e xMTBE,spec xMTBE B B xMTBE g xMTBE,spec
where the path constraints hx(‚) must be satisfied at all times during the system operation.60,61 The relevant parameters and schematic representation of the optimized design are given in Table 6 and Figure 2, respectively. The entropy profiles along the column for this reference design are presented in Figure 3 and calculated according to eq 8. The following results are of interest, (a) The entropy profiles show that energy transfer contributes largely to the overall entropy production and particularly in the stripping section of the column, followed by mass diffusion. In the reactive zone, however, the dominant source of entropy production is chemical reaction. In the rectifying section, both mass and energy-transfer contributions are equally important, (b) The entropy contribution originated from chemical reaction shows negative values at the bottom of the reactive section. This fact may be explained by possible inconsistency between the reaction kinetics model and the affinity expression. Although MTBE decomposition takes place at that column location, due to high MTBE concentration, the reaction affinity is excepted also to change sign exactly at the same point in the (P,T,x) space. These facts result always in a nonnegative entropy production contribution. Moreover, second law requires only that the sum of all contributions to the entropy production is larger than zero, (c) At the feed stages (9th-10th), the conjugated driving forces are considerably larger than elsewhere due to the
parameter
value
number of trays feed 1 location flow rate composition temperature feed 2 location flow rate composition temperature operating pressure reflux ratio mass of catalyst/tray bottom flow rate column diameter condenser area reboiler area operational constraints
15 [2 rectifying, 8 reactive, 5 stripping]
product specifications
iC4 + nC4, 10th tray (g) 455 mol s-1 iC4, 0.37; nC4, 0.63 350 K MeOH, 9th tray (l) 168 mol s-1 MeOH, 1.00 320 K 11 × 105 Pa (top) 7 204.1 kg )1000 eq[H+] 197 mol s-1 7.132 m 1.059 × 103 m2 6.93 × 102 m2 P e Pcond Treb < Tc,iC4 xDMTBE e 0.1% xBMTBE g 99%
difference of chemical potential and temperature between the entering streams and the mixture coexisting on those stages, (d) External heat transfer plays a dominant role in the reboiler and condenser, exclusively. In terms of entropy production, summarizing, the separation function dominates over the reaction function for this MTBE RD column configuration. (2) Biobjective Optimization of a MTBE Reactive Distillation Column. A multiobjective optimization problem is formulated for the RD column with respect to economic performance and exergy efficiency (see Table 7). The formulation includes the balances described by eqs 1-5 and 18, the criteria definitions given in eqs 29, and the optimization formulation described by eq 30. Constraints imposed by operating conditions and product specifications are included. For the sake of building up the optimization complexity, the first approximation of this approach omits the response time constant as objective function. The resulting optimization problem is formulated as follows and is implemented and solved in gPROMS/gOPT,
Figure 3. Entropy production rate profile for a 15-stage RD column for MTBE synthesis. Stages are numbered from top to bottom; reboiler and condenser drum are not depicted; MeOH feed stream is fed at 9th tray and iC4 stream at 10th tray. The operational and design variables are listed in Table 6.
58 Ind. Eng. Chem. Res., Vol. 47, No. 1, 2008 Table 7. Optimized Design of a RD Column for MTBE Synthesis Based on Economic Performance and Exergy Efficiencya parameter
value
number of trays feed 1 location flow rate composition temperature feed 2 location flow rate composition temperature operating pressure reflux ratio mass of catalyst (eq[H+])
15 [2 rectifying, 8 reactive, 5 stripping] iC4 + nC4, 10th tray (g) 455 mol s-1 iC4, 0.37; nC4, 0.63 350 K MeOH, 9th tray (l) 168 mol s-1 MeOH, 1.00 320 K 11 × 105 Pa (top) 7 m1:2 ) 0; m3 ) 771.9 m4 ) 826.122; m5 ) 762.405 m6 ) 688.412; m7 ) 652.292 m8 ) 785.991; m9 ) 1359.42 m10 ) 982.448; m11:15 ) 0 197 mol s-1
bottom flow rate column diameter rectifying reactive stripping condenser area reboiler area a
Figure 4. Pareto optimal curve f ′econ versus f ′exergy. Each point of the curve represents an optimal solution in terms of the design variables [mcat, Q, aΨ] and weighting factors w; axes are normalized so that utopia point coordinates are as follows: U ) (0,0); original coordinates of U are: U0 ) (5.251 × 104 J × (s × K)-1,2.082 × 108 euro × y-1).
6.95 m 6.95 m 6.33 m 1.059 × 103 m2 6.93 × 102 m2
w ) [0.3 0.7]T.
min J(fecon(x),fexergy(x),w)
(34)
D∈Γ
s.t. gx(‚) ) 0 hx(‚) e 0 where the set of design variables D include variables as ([mcat Q aΨ]). Operational variables such as reflux ratio and condenser pressure have been fixed based on the explanation given in Design Decision Variables at appropriate values as listed in Table 6. The inequality constraints hx(‚) are defined by operating conditions and product specifications (cf. bottom section of Table 6). The planar view of the two normalized objective functions at different weighting factors is depicted in Figure 4. The normalization of the objectives requires a priori knowledge of the minimum (fmin) and maximum (fmax) values for each objective. These values are found by solving the optimization problem at different weighting factors. Each objective value is then scaled based on the corresponding utopia coordinate and for a given weighting factor. Thus, the normalized objective functions (f′) are computed according to the expression,
f′i )
hf i - hf min i hf max - hf min i i
( )
hf i ) wi ×
fi
fUi
i ∈ Znp
deviations between optimal designs and the utopia point (up to 60%) for the pair f ′econ - f ′exergy.66 The decision maker should have in mind the importance of compromising both objectives at early stages of the design cycle. Significant structural differences exist between the economicdriven (E) and exergy-driven (X) designs, which are respectively given as
E; w ) [1 0 ]T (35)
where the objective function hfi is scaled with respect to the corresponding coordinates of the utopia point, p
Figure 5. Normalized catalyst distribution in the synthesis of MTBE with respect to economic performance and exergy efficiency.
(36)
The curve of Figure 4 is termed Pareto curve and clearly shows that one objective function can only be improved at the expense of the other objective function.1,62 The economicsexergy trade-off has been extensively studied in energy-related literature27,63,64 and in the engineering branch of exergoeconomics.65 In our RD case study, this trade-off leads to significant
and
X; w ) [0 1 ]T As depicted in Figure 5, for instance, the catalyst load along the reactive trays differs among the proposed designs when compared to the reference case. As expected, the optimal design based on economic performance (E) coincides with the nominal case (R). The catalyst load of the nominal case is used to normalize the plot. When exergy is fully considered as an objective function, the last reactive tray vanishes together with the sixth reactive tray. Furthermore, the catalyst load of the seventh reactive tray is considerably higher than that of the nominal case. This output suggests an alternative design, where a side reactor with high catalyst load is placed at the seventh
Ind. Eng. Chem. Res., Vol. 47, No. 1, 2008 59 Table 8. Entropy Produced in Classical and Green Designsa w [1.0 [0.3 [0.0
0.0]T 0.7]T 1.0]T
fecon (euro y-1)
fexergy (J (s K)-1)
2.082 × 2.085 × 108 2.086 × 108
5.604 × 104 5.305 × 104 5.251 × 104
108
a Classical design corresponds to w ) [1 0]T; green design corresponds to w ) [0.3 0.7]T.
Figure 6. Entropy production rate profile for an optimal design of a MTBE RD column based on exergy efficiency (X-design). Stages are numbered from top to bottom; MeOH feed stream is fed at 9th tray and iC4 stream at 10th tray.
tray. The catalyst load on the remaining reactive trays decreases when exergy is taken into account. The design with equally weighted objectives (Q; w ) [0.5 0.5]T) presents an interesting catalyst load profile, specially at the ninth tray. The rather large catalyst load might be explained by the reduced catalyst load on the last reactive tray and the constraint of targeting the product specifications. The reduced catalyst load on the ninth tray diminishes the MTBE decomposition at high MTBE composition and the product synthesis at low MTBE composition. As a constraint in the product specification is included in the problem formulation, the required extent of reaction is accommodated only by having the ninth tray fully loaded with catalyst. Since MeOH is fed at the ninth tray, larger holdup of reactant is present and only converted into products by an increased amount of catalyst. The effect of catalyst amount on the performance of a RD column has also been addressed,67 aiming at further internal heat integration between the reaction and separation tasks. When the entropy profiles of the X-design (cf. Figure 6) are compared to those of the E-design, it can be realized that the mayor contribution of entropy is still coming from the heattransfer process in the stripping section. However, mass diffusion becomes more relevant than in the E-design. In the reactive zone and in contrast to the results of Figure 3, the exergy by chemical reaction is unevenly produced along the reactive trays. Thus, the ninth tray largely contributes to the overall exergy produced within the reactive section. A fully catalyst load on that tray intuitively explains this behavior. According to these last results, for the exergy-driven design (X-design), it seems more attractive to reach high conversions on one or a few trays and low conversion on surrounding trays. Apparently, this finding seems to disagree with the theorem of equipartition of driving forces.68,69 As demonstrated by Johannessen and Rosjorde40 and mentioned by Jimenez et al.,63 this theorem is a good approximation to the state of minimum entropy production in an optimally controlled system, only if the system has sufficient freedom. For instance, in conventional distillation equipartition of driving forces results in a diabatic separation column with minimum total entropy production at large areas and at large number of trays.40,63 (3) Triobjective Optimization of a MTBE Reactive Distillation Column: A Sensitivity-Based Approach. For the sake of simplification, the following strategy is adopted: (i) an optimized design based in economics and exergy efficiency is considered; (ii) a disturbance scenario is defined; and (iii) the
responsiveness index of the optimized design is estimated at the given disturbance scenario. Optimized Design. According to applications 1 and 2, the optimized designs based on fecon and fexergy present an important tradeoff between both objectives. This tradeoff is dependent on the vector of weighting factors (w). The ultimate decision on w is made by the designer, according to the requirements of the design problem. In our case, and only for tutorial purposes, we choose the vector
w ) [0.3 0.7 ]T over all available designs. The optimized design variables for this case are listed in Table 8. Disturbance Scenario. The selection of a disturbance scenario is carried out by analyzing qualitatively the entropy profiles of Figure 3 (E-design) and Figure 6 (X-design) and by focusing on process disturbances occurring in industrial practice. An upstream variation of the MeOH feed flow rate up to -10% is chosen as disturbance scenario.
FMeOH ) FSS MeOH × (1 - Rτ(t)) (for 0 e Rτ(t) e 0.1, t ∈[0,tf]) (37) The disturbance scenario is chosen so that any steady-state transition is avoided. Control Structure. To mitigate any effect caused by a disturbance scenario, a control structure is proposed for the RD unit shown in Figure 2. The following assumptions are considered: (i) four control valves are associated with the RD unit; (ii) the unit is operated at constant reflux ratio RR () FR,L/FD,L ) 7); and (iii) the throughput of the unit is set. Based on a rigorous degree of freedom analysis45 a 4 × 4 SISO control structure is suggested. The controlled variables are chosen aiming at the satisfaction of product specifications and protection of the unit. The pairing of controlled-manipulated variables, on the other hand, tries to minimize cross-interactions between manipulations of input variables on the controlled output variables and looks for a dominant 1-1 correspondence between input and output variables. Thus, two control valves are used to control the liquid levels in the reflux drum and column bottom. As the reflux ratio is larger than a value of 4 (RR > 4), the reflux drum level is controlled with the liquid reflux or drum liquid outflow (Richardson’s rule31). The level at the column bottom might be by convention controlled with the reboiler liquid outflow. The third control valve is used to manipulate the flow rate of cooling utility to control the column pressure. The remaining control degree of freedom is used to control the purity of the MTBE product by manipulating the reboiler heat duty. The composition of MTBE at the bottom stream is conveniently inferred from temperature measurements at a given tray. Steady-state open-loop simulations are carried out to assess the effect of varying manipulated variable (i.e., reboiler duty Qreb). The simulation results indicate45 that the largest temperature changes are occurring either at the rectifying section, the top section of the reactive section or at the reboiler. By selecting the reboiler temperature to infer the MTBE
60 Ind. Eng. Chem. Res., Vol. 47, No. 1, 2008
Figure 7. Schematic representation of a MTBE RD column with a 4 × 4 SISO control structure. Numbers represent tray location.
composition in the bottom stream, any hydraulic lag is effectively avoided. This single-end point composition control structure is originated from the assumption of invariant reflux ratio. As mentioned by Luyben et al.,31 this control structure is more advantageous from both operational (i.e., easy tuning) and economical (i.e., comparable energy demand) points of view than the dual composition control. Figure 7 depicts the proposed control structure. The controllers are tuned by solving an optimization problem in the time domain, aiming at a minimization of the integral absolute error of the control loops,70
Function: s.t.
min
Ki,Bi,τc,i,SPi
F)
fx(‚) ) 0
∫0t [(xb - xjb)2 + (xd - xjd)2] dt f
t ∈ [t0,tf]
gx(‚) ) 0
t ∈ [t0,tf]
hx(‚) e 0
t ∈ [t0,tf]
dx(‚) ) 0
t ∈ [t0,tf]
ix(‚) e 0
t ∈ [t0,tf]
max Tmin reb e Treb e Treb
Responsiveness Index. For each MeOH feed rate, the variables B(x,u) and ΣV(x,u) are estimated in the corresponding steady state and the responsiveness index is computed. The responsiveness index is treated as a function of the magnitude of the MeOH feed perturbation, fresponse(FMeOH). Specifically, this MeOH flow rate disturbance results in a variation of the mass-transfer driving forces for each component along the unit. From this point on, the analysis is focused on the top tray (1st), the 9th tray, and the bottom tray (15th) of the RD column specified in Table 8. Considering these tray locations, we provide a generalized overview on the unit’s responsiveness in the rectifying, reactive and stripping sections. The 9th tray is chosen because it is on that particular feeding tray where MeOH is fed to the column and, therefore, where a more pronounced effect of feed flow rate disturbance is expected. The corresponding mass-transfer and chemical reaction driving forces are plotted for the 1st, 9th, and 15th trays in Figure 8. It can be seen from these plots that both MeOH masstransfer driving forces increase with larger disturbances, whereas chemical reaction driving force decreases. The response time of the MTBE RD column as a function of Xm,MeOH (1st and 15th trays) and Xrx (9th tray) is depicted in Figure 9. The variations of B(x,u) and ΣV(x,u) are computed between the initial reference state and between consecutive states. The response time becomes larger at small Xm,MeOH-values in the rectifying and stripping sections, whereas small driving forces lead to short response times in the reactive section. From a responsiveness perspective, therefore, small driving forces are desirable for reactive sections and large driving forces are desirable for separation (nonreactive) sections. The first statement suggests that in the case of reactive sections driving forces should be manipulated to control the unit while flux is specified (i.e., a certain conversion rate). To explain this statement. we require the introduction of an approximated relation between fluxes and forces for a kth compartment,31,73
Jk ≈ R × Xk T0 × ∆ΣV ≈ 2 × T0
∫V dV ∑ R × Xk × ∆Xk
(39) (40)
k
t ∈ [t0,tf]
out out,max Tout,min water e Twater e Twater
t ∈ [t0,tf]
out out,max Tout,min steam e Tsteam e Tsteam
t ∈ [t0,tf]
(38)
where the objective function F accounts for deviations from product and top specifications (xj) and the constraint paths set bounds for the temperature in the reboiler and the outlet temperatures of the cold and hot utilities. Alternative control configurations are proposed for the MTBE synthesis by Wang et al.71 and Sneesby et al.72 addressing switchover to different steady states. As the focus of this contribution is the dynamic behavior around a known steady state, no steady-state transition is considered by constraining the disturbances’ variability.
where R is a design factor, which Meeuse called the thermodynamic design factor.46 At small gradient X, the design factor R should be large enough to obtain the required flux. In the reactive section, small reaction driving forces require large reaction volumes on the trays. Small variations in the MeOH feed flow rate produce large changes in the affinity and subsequently large entropy production rate values. As the stored exergy changes only modestly with a small flow disturbance, the overall effect enhances the process responsiveness (i.e., small response times). The second statement can be explained by following the same reasoning. In a distillation (separation) process, the driving force (top-bottom) is specified and the flux remains to influence. Thus, in nonreactive sections, flux should be manipulated to control the unit. The results of these observations fully agree with the analysis given by Luyben et al.31 and Meeuse46 for a reactor and distillation column. Meaningfulness and Validity of Results. In the current approach, a continuous range of disturbances for one input variable is considered for the computation of the responsiveness time constant. In comparison with the generalized definition of
Ind. Eng. Chem. Res., Vol. 47, No. 1, 2008 61
Figure 8. Driving forces as a function of the MeOH feed flow rate. MeOH feed flow rate varies from the nominal value (168 mol s-1) to 150 mol s-1; the optimal design variables and sizing parameters corresponds to a fecon - fexergy design, with w ) [0.3 0.7]T.
Figure 9. Response time as a function of the MeOH feed flow rate. MeOH feed flow rate varies from the nominal value (168 mol s-1) to 150 mol s-1; the optimal design variables and sizing parameters corresponds to a fecon - fexergy design, with w ) [0.3 0.7]T.
fresponse by Luyben et al.,31 this is a very specific case. Their (Luyben’s et al.) definition demands an optimization of the states x to find the worst case situation over a domain of states. Therefore, there is a latent risk that the presented results are
not optimized over the feasible domain of states for normal operation. Our approach might be locally optimized for the given operational window, but might have a low degree of validity to other cases.
62 Ind. Eng. Chem. Res., Vol. 47, No. 1, 2008
Comparison between Classical and Green Designs In the course of this contribution, the economic-driven design (classical design) was implicitly compared to a new bioptimized design (green design) for the case of a MTBE RD column. A set of indicators can be used to compare both designs: (i) exergy loss; and (ii) closed-loop control performance. The exergy loss (or entropy produced) can be extracted from the Pareto curve presented in Figure 4. The relevant data for the classical and green designs and for the fully exergy-driven design, i.e.,
w ) [0 1 ]T are listed in Table 8. As was expected, the green design produces less entropy than the classical one; however, economically it is less attractive than the conventional design, as it presents a higher fecon value. The increase in costs is rather marginal (∼0.2%) and within the uncertainty range of the correlation to predict investments costs for units and catalyst. A tendency toward increasing cost can be suggested though,66 it being not significant for this particular case. A more significant change is observed for fexergy (∼6%). A more revealing concept to assess the impact of making a process greener is the shadow price λp. It is a cost penalty function for changing entropy production and given by
λp )
∂fecon ∂fexergy
(41)
The shadow prices can be computed directly out from the data in Table 8 or from the Pareto curve in Figure 4. It was found that, when compared to the classical case, the green design presents a cost penalty of approximately -2.02, whereas the X-design presents a value of approximately -4.37 [euro × K × MJ-1]. These values confirm the tradeoff between both objectives and suggest that the greener the design the larger the involved cost penalty. The closed-loop control performance of both designs is now analyzed. The disturbance scenario under consideration represents a step variation in the flow rate of methanol during the time horizon of 6 × 104 s and given by,
FMeOH )
{
FSS MeOH 1.1 ×
FSS MeOH
(
)
if t ∈ [0 s, 1 × 104 s] or t ∈ [4 × 104 s, 6 × 104 s] (42) 4 4 (if t ∈ [1 × 10 s, 4 × 10 s])
The dynamic response of the classic and green designs are depicted in Figure 10. In this figure, the MTBE composition at the bottom stream is presented as a function of time in the presence of the MeOH feed disturbance. The bioptimized design (i.e., green) has an improved closed-loop performance compared to the classical design. The deviations from xMTBE set point are less pronounced in the green design, which smoothly returns back to the initial steady state when the disturbance disappears. The set point of the green design is slightly lower (0.05%) than the one of the classic design. This negligible variation in set points does not affect the quality of the analysis results. It is suggested from our results that incorporating both economic and exergy-related objectives in the unit design results in a process with better closed-loop performance, provided that control structure and tuning are kept invariant.A well-tuned control can compensate for any change in the process operability feature.
Figure 10. Time variation of MTBE product stream for the classic and green designs in the presence of a MeOH feed flow rate disturbance. Green design corresponds to w ) [0.3 0.7]T; column specifications for classic and green designs are given in Tables 6 and 8, respectively.
Concluding Remarks The starting point of this contribution was the definition of a generic building block for process synthesis (GLRDVE). Being embedded in a superstructure, each GLRDVE block has its own complexities due to the many interacting phenomena. However, the modeling framework of the GLRDVE does not account for the interactions between the stages in a real multistage column. To overcome this limitation, operational variables such as reflux ratio and condenser pressure were specified at a higher decisionmaking layer aiming at operating windows where steady-state multiplicity and reactive azeotropy phenomena are not likely to occur. This a priori specification might be considered a drawback of the proposed approach, but it is certainly a topic to be further researched. The sources and the distribution of the entropy production rate in a complete MTBE column were then determined and analyzed. It is believed that insight into such a distribution is helpful for assessing which column-related design decision variables (e.g., feed tray location, reflux ratio, heat load distribution using diabatic operation, catalyst load, and mass/heat-transfer area) are key to optimize the column performance. The perspective adopted for selecting the design decision variables was task-oriented rather than the conventionally used equipment-oriented design approach. The chosen design decision variables included the reaction volume, the amount of heat to be added or removed, and the interfacial areas for species mass and heat transfer. Analyzing the steady-state entropy production profiles in a MTBE RD column showed that, in terms of entropy production, the separation function dominates over the reaction function. A biobjective function was then studied regarding economic and exergy performance. It is suggested that designs based exclusively on economic or exergy performances differ significantly in their structure, showing an existing trade-off. The Pareto curve provides sensitivity information on the cost effects of improving the exergy efficiency. In addition to economics and exergy objectives, this optimization-based approach covers responsiveness aspects. To define a responsiveness criterion, the availability and entropy production rate were regarded as functions of an internal state vector x and a vector of disturbances or manipulated external excitations u. In our approach, the maximization approach suggested by Luyben et al.31 was not considered due to its inherent need for high computational effort. However, a rather formal definition of fresponse was adopted together with a single input excitation. By not searching over all possible state trajectories,
Ind. Eng. Chem. Res., Vol. 47, No. 1, 2008 63
the generality of the results might be questionable and the external excitation is a fairly conservative one. By analyzing the responsiveness criteria of the bioptimized design, it is reconfirmed that small driving forces are desirable for reactive sections and large driving forces are desirable for separation (nonreactive) sections.31,46 The first statement suggests that, in the case of reactive sections, driving forces should be manipulated to control the unit. The second statement suggests that in nonreactive sections flux should be manipulated to control the unit. Finally, the economic-driven design (classical design) was compared to the new bioptimized design (green design) for the case of the MTBE RD column. The green design produced less entropy than the classical one, while being marginally less attractive from an economic standpoint. The bioptimized design had an improved closed-loop performance compared to the classical design, keeping the same control structure and settings. It is concluded that incorporating both economic and exergyrelated objectives in the unit design results in a process with better closed-loop performance and reduced exergy loss. In this contribution, we have focused on the performance of RD process designs, paying less explicit attention to the industrial conventional reaction-separation processing route. Therefore, a comparison between RD and non-RD processes with regard to entropy production and the design implications for the principle of equipartition of forces are recommended topics to be further researched. Nomenclature List of Symbols A ) area B ) controller’s bias Dcol ) column diameter j,V ) F j,V total ) total flow of stream j in the vapor phase; Ftotal i)nc Σ i)1 Fij,V F j,(R) ) convective flow of component i in phase R in the i incoming or outgoing stream F j,k i ) convective flow of component i in compartment k in the incoming or outgoing stream j for the exchange with the external world J (R) e ) total energy flux in phase (R) J Lq ) measurable heat flux into the liquid phase J (R) s ) entropy flux in phase (R) (R) J m,i ) mass flux of component i in phase (R) K ) controller’s gain Keq ) temperature-dependent chemical equilibrium constant M(R) ) molar holdup in phase (R) P ) absolute pressure Q(R) ) external heat exchange rate in phase (R) RR ) reflux ratio SW ) variable in Wilson model Se ) entropy produced by interfacial energy transfer S j,(R) ) molar entropy of component i in stream j and phase i (R) Sm ) entropy produced by interfacial mass diffusion Srx ) entropy produced by chemical reaction T(R) ) absolute temperature in phase (R) V ) volume of compartment V LRX ) reaction volume in the liquid phase Wideal ) minimum work required by the process Xq ) conjugated driving forces for heat exchange between phases Xe ) conjugated driving forces for interfacial heat diffusion
Xm,j ) conjugated driving forces for interfacial mass diffusion of component j XHX ) conjugated driving forces for heat exchange with the external world Xrx,j ) conjugated driving forces for the chemical reaction of reaction j hf ) scaled objective function xj ) product specification E ) economic-based optimal design Q ) optimal design based equally on economic performance and exergy efficiency U ) Utopia point X ) exergy-based optimal design A Li ) affinity of chemical reaction i in the liquid phase D ) set of design or decision variables R ) universal gas constant A ) availability B ) availability function B(x,u)total ) availability over all the superstructure aΨ ) mass and heat transfer area ai ) activity of component i aij, b ij ) Wilson model parameters c ) investment cost of unit and catalyst; cost of component or utility per molar unit dx(‚) ) set of disturbances f ) objective function f′ ) normalized objective function fecon ) economic performance indicator fexergy ) sustainability performance index fresponse ) responsiveness performance index fx(‚) ) set of differential expressions ) external force exerted on component i per molar unit in gj,(R) i stream j and phase (R) gx(‚) ) set of algebraic expressions grx(‚) ) chemical reaction driving force based on activities hj,(R) ) molar enthalpy of component i in stream j and phase i (R) hx(‚) ) set of inequalities and path constraints ix(‚) ) set of initial conditions m ) mass nc ) number of components no ) number of objective functions np ) number of phases nt ) number of compartments in the integrated structure nrx ) number of chemical reactions nst ) number of incoming and outgoing streams p ) norm q(R′ - R) ) heat flux from phase (R′) to phase (R) ri(R) ) net change of the amount of component i in phase (R) due to all reactions s(R) ) entropy density in phase (R) t ) time u ) internal energy density of the volume element w ) vector of weighting factors x ji ) liquid molar fraction of component i in stream j y ji ) vapor molar fraction of component i in stream j R ) nominal design U(R) ) internal energy of phase (R) or compartment Acronyms nC4 ) n-butane iC4 ) isobutene GLRDVE ) generic lumped reactive distillation volume element MeOH ) methanol
64 Ind. Eng. Chem. Res., Vol. 47, No. 1, 2008
MTBE ) methyl tert-butyl ether RD ) reactive distillation SP ) set point TAC ) total annualized cost F ) multiobjective function Greek Symbols R ) thermodynamic design factor δ(R) ) Kronecker delta function ηII ) efficiency of the 2nd law of thermodynamics Γ ) set of constraints γi ) liquid activity coefficient for component i λp ) shadow price Λij ) Wilson model binary interaction parameter for the pair i, j ∈ Znc j,(R) µi ) chemical potential of component i in stream j and phase (R) ν (R) i,j ) stoichiometric coefficient of component i in reaction j and phase (R) Ψe ) energy flux associated with the interfacial mass and heat transfer Ψi ) interfacial molar flux of component i ) molar density of component i in phase (R) F(R) i (R) ΣV ) total entropy production rate in phase (R) ΣV(x,u) ) total rate of entropy production over all superstructure σ (R) V ) entropy production rate per unit of volume in phase (R) τ ) settling or response time τc ) controller’s response time Li ) degree of advancement of chemical reaction i in liquid phase Subscripts and Superscripts 0 ) initial value; infinitely large sink with similar surrounding conditions U ) Utopia point S(R) ) entropy of compartment in phase (R) B ) bottom stream cu ) cool utility hu ) hot utility D ) distillate stream e ) energy-based f ) final HX ) external heat exchange m ) mass-based Q ) heat source or heat sink s ) entropy-based cat ) catalyst in ) inlet conditions low ) lower bound max ) maximal value min ) minimal value out ) outlet conditions reb ) reboiler spec ) specification value SS ) steady state total ) total value or combined stream u ) utility up ) upper bound Vectors g(R) ) external forces in phase (R) u ) disturbances or manipulated excitations
x ) state variables Literature Cited (1) Hakanen, J.; Hakal, J.; Manninen, J. An Integrated Multiobjective Design Tool for Process Design. Appl. Therm. Eng. 2006, 26, 1393. (2) Rodriguez, I. E.; Zheng, A.; Malone, M. F. Parametric Dependence of Solution Multiplicity in Reactive Flashes. Chem. Eng. Sci. 2004, 59, 1589. (3) Lakerveld, R.; Bildea, C. S.; Almeida-Rivera, C. P. Exothermic Isomerization Reaction in Reactive Flash: Steady-State Behavior. Ind. Eng. Chem. Res. 2005, 44, 3815. (4) Jacobs, R.; Krishna, R. Multiple Solutions in Reactive Distillation for Methyl tert-Butyl Ether Synthesis. Ind. Eng. Chem. Res. 1993, 32, 1706. (5) Gu¨ttinger, T. E. Multiple Steady States in Azeotropic and Reactive Distillation. Ph.D. Dissertation, Swiss Federal Institute of Technology, Zu¨rich, 1998. (6) Kienle, A.; Groebel, M.; Gilles, E. D. Multiple Steady States in Binary Distillation. Theoretical and Experimental Results. Chem. Eng. Sci. 1995, 50, 2691. (7) Ciric, A.; Miao, P. Steady-State Multiplicites in a Ethylene Glycol Reactive Distillation Column. Ind. Eng. Chem. Res. 1994, 33, 2738. (8) Hauan, S.; Hertzberg, T.; Lien, K. Multiplicity in Reactive Distillation of MTBE. Comput. Chem. Eng. 1995, 19, S327. (9) Hauan, S.; Hertzberg, T.; Lien, K. Why Methyl tert-Butyl Ether Production by Reactive Distillation May Yield Multiple Solutions. Ind. Eng. Chem. Res. 1995, 34, 987. (10) Gu¨ttinger, T. E.; Morari, M. Predicting Multiple Steady States in Distillation: Singularity Analysis and Reactive Systems. Comput. Chem. Eng. 1997, 21, S995. (11) Gu¨ttinger, T. E.; Morari, M. Predicting Multiple Steady States in Equilibrium Reactive Distillation. 2. Analysis of Hybrid Systems. Ind. Eng. Chem. Res. 1999, 38, 1649. (12) Gu¨ttinger, T. E.; Morari, M. Predicting Multiple Steady States in Equilibrium Reactive Distillation. 1. Analysis of Nonhybrid Systems. Ind. Eng. Chem. Res. 1999, 38, 1633. (13) Higler, A.; Taylor, R.; Krishna, R. Nonequilibrium Modelling of Reactive Distillation: Multiple Steady States in MTBE Synthesis. Chem. Eng. Sci. 1999, 54, 1389. (14) Chen, F.; Huss, R. S.; Doherty, M.; Malone, M. F. Multiple Steady States in Reactive Distillation: Kinetic Effects. Comput. Chem. Eng. 2002, 26, 81. (15) Sneesby, M. G.; Tade, M. O.; Smith, T. N. Mechanistic Interpretation of Multiplicity in Hybrid Reactive Distillation: Physically Realizable Cases. Ind. Eng. Chem. Res. 1998, 37, 4424. (16) Sneesby, M. G.; Tade, M. O. R. S. Multiplicity and Pseudomultiplicity in MTBE and ETBE Reactive Distillation. Chem. Eng. Res. Des. 1998, 76, 524. (17) Mohl, K. D.; Kienle, A.; Gilles, E. D.; Rapmund, P.; Sundmacher, K.; Hoffmann, U. Steady-state Multiplicities in Reactive Distillation Columns for the Production of Fuel Ethers MTBE and TAME: Theoretical Analysis and Experimental Validation. Chem. Eng. Sci. 1999, 54, 1029. (18) Malone, M. F.; Doherty, M. F. Reactive Distillation. Ind. Eng. Chem. Res. 2000, 39, 3953. (19) Barbosa, D.; Doherty, M. Design and Minimum-reflux Calculations for Single-feed Multicomponent Reactive Distillation Columns. Chem. Eng. Sci. 1988, 43, 1523. (20) Song, W.; Huss, R. S.; Doherty, M.; Malone, M. F. Discovery of Reactive Azeotropes. Nature 1997, 388, 561. (21) Frey, T.; Stichlmair, J. Reactive Azeotropes in Kinetically Controlled Reactive Distillation. Trans. Inst. Chem. Eng., Part A 1999, 77, 613. (22) Harding, S. T.; Floudas, C. A. Locating All Heterogeneous and Reactive Azeotropes in Multicomponent Mixtures. Ind. Eng. Chem. Res. 2000, 39, 1576. (23) Maier, R.; Brennecke, J.; Stadtherr, M. Reliable Computation of Reactive Azeotropes. Comput. Chem. Eng. 2000, 24, 1851. (24) Frey, T.; Stichlmair, J. Thermodynamic Fundamentals of Reactive Distillation. Chem. Eng. Technol. 1999, 22, 11. (25) Ung, S.; Doherty, M. Necessary and Sufficient Conditions for Reactive Azeotropes in Multireaction Mixtures. AIChE J. 1995, 41, 2383. (26) Apaiah, R.; Kinnemann, A.; vander Kooi, H. J. Exergy Analysis: A Tool to Study the Sustainability of Food Supply Chains. Food Res. Int. 2006, 39, 1. (27) Rivero, R.; Garcia, M.; Urquiza, J. Simulation, Exergy Analysis and Application of Diabatic Distillation to a Tertiary Amyl Methyl Ether Production Unit of a Crude Oil Refinery. Energy 2004, 29, 467. (28) Koeijer, G.; Rivero, R. Entropy Production and Exergy Loss in Experimental Distillation Columns. Chem. Eng. Sci. 2003, 58, 1587.
Ind. Eng. Chem. Res., Vol. 47, No. 1, 2008 65 (29) Takeshita, K.; Ishida, M. Optimum Design of Multi-stage Isotope Separation Process by Exergy Analysis. Energy 2006, 31, 3097. (30) Iwakabe, K.; Nakaiwa, M.; Huang, K.; Nakanishi, T.; Rosjorde, A.; Ohmori, T.; Endo, A.; Yamamoto, T. Energy Saving in Multicomponent Separation Using an Internally Heat-integrated Distillation Column (HIDiC). Appl. Therm. Eng. 2006, 26, 1362. (31) Luyben, W. L.; Tyreus, B. D.; Luyben, M. L. Plantwide Process Control; McGraw-Hill: New York, 1999. (32) Kjelstrup, S.; Bedeaux, D. Elements of IrreVersible Thermodynamics; International Centre for Applied Thermodynamics: Istanbul, 2001. (33) Marquardt, W. Modeling and Analysis of Chemical Process Systems. Presented at Delft University of Technology, Delft, The Netherlands, July 2004. (34) Bird, B.; Stewart, W.; Lightfoot, E. Transport Phenomena; John Wiley and Sons, Inc.: New York, 2003. (35) Seider, W.; Seader, J. D.; Lewin, D. Process Design Principles; John Wiley and Sons, Inc.: New York, 1999. (36) Bedeaux, D.; Kjelstrup, S. Irreversible Thermodynamics. A Tool to Describe Phase Transitions Far from Global Equilibrium. Chem. Eng. Sci. 2004, 59, 109. (37) Koeijer, G. Energy Efficient Operation of Distillation Columns and a Reactor Applying Irreversible Thermodynamics. Ph.D. Dissertation, Norwegian University of Science and Technology, Trondheim, Norway, 2002. (38) Kjelstrup, S. Norwegian University of Science and Technology at Trondheim, Trondheim, Norway. Personal communication, 2005. (39) Koeijer, G.; Kjelstrup, S.; vander Kooi, H. J.; Gros, B.; Knoche, K. F.; Andersen, T. R. Positioning Heat Exchangers in Binary Tray Distillation Using Isoforce Operation. Energy ConVers. Manage. 2002, 43, 1571. (40) Johannessen, E.; Rosjorde, A. Equipartition of Entropy Production as an Approximation to the State of Minimum Entropy Production in Diabatic Distillation. Energy 2007, 32, 467. (41) Ismail, S.; Proios, P.; Pistikopoulos, E. Modular Synthesis Framework for Combined Separation/Reaction Systems. AIChE J. 2001, 47, 629. (42) Papalexandri, K.; Pistikopoulos, E. Generalized Modular Representation Framework for Process Synthesis. AIChE J. 1996, 42, 1010. (43) Bermingham, S. A Design Procedure and Predictive Models for Solution Crystallisation Processes. Ph.D. Dissertation, Delft University of Technology, Delft, The Netherlands, 2003. (44) Coulson, J. M.; Richardson, J. F. Chemical Engineering Design; Chemical Engineering, Vol. 6; Butterworth-Heinemann: Oxford, Boston, 1997. (45) Almeida-Rivera, C. P. Designing Reactive Distillation Processes with Improved Efficiency. Ph.D. Dissertation, Delft University of Technology, Delft, The Netherlands, 2005. (46) Meeuse, M. On the Design of Chemical Processes with Improved Controllability Characteristics. Ph.D. Dissertation, Delft University of Technology, Delft, The Netherlands, 2003. (47) Clark, P.; Westerberg, A. Optimization for Design Problems Having More than One Objective. Comput. Chem. Eng. 1983, 7, 259. (48) Jimenez, L.; Wanhschafft, O.; Julka, V. Analysis of Residue Curve Maps of Reactive and Extractive Distillation Units. Comput. Chem. Eng. 2001, 25, 635. (49) Peters, U.; Nierlich, F.; Schulte-Ko¨rne, E.; Sakuth, M. Methyl tertButyl Ether. In Ullmann’s Encyclopedia of Industrial Chemistry; John Wiley and Sons, Inc.: New York, 2000. (50) Sundmacher, K.; Hoffmann, U. Multiple Reactions in Catalytic Distillation Processes for the Production of the Fuel Oxygenates MTBE and TAME: Analysis by Rigourous Model and Experimental Validation. Chem. Eng. Sci. 1999, 54, 2839. (51) Schrans, S. M.; Wolf, S.; Baur, R. Dynamic Simulation of Reactive Distillation: An MTBE Case Study. Comput. Chem. Eng. 1996, 20, S1619. (52) Malanowski, S.; Anderko, A. Modelling Phase Equilibria. Thermodynamic Background and Practical Tools; John Wiley and Sons, Inc.: New York, 1992.
(53) Kooijman, H.; Taylor, R. The ChemSep Book; Library books on demand: Norderstedt, 2000. (54) Colombo, F.; Cori, L.; Dalioro, L.; Delogu, P. Equilibrium Constant for the Methyl tert-Butyl Ether Liquid Phase Synthesis by Use of UNIFAC. Ind. Eng. Chem. Fundam. 1983, 22, 219. (55) Rehfinger, A.; Hoffmann, U. Kinetics of Methyl Tertiary Butyl Ether Liquid Phase Synthesis Catalyzed by Ion Exchange ResinsI. Intrinsic Rate Expression in Liquid Phase Activities. Chem. Eng. Sci. 1990, 45, 1605. (56) Zhang, T.; Datta, R. Integral Analysis of Methyl tert-Butyl Ether Synthesis Kinetics. Ind. Eng. Chem. Res. 1995, 34, 730. (57) Izquierdo, J.; Cunill, F.; Vila, M.; Tejero, J.; Iborra, M. Equilibrium Constants for Methyl tert-Butyl Ether Liquid-phase Synthesis. J. Chem. Eng. Data 1992, 37, 339. (58) Thiel, C.; Sundmacher, K.; Hoffmann, U. Residue Curve Maps for Heterogeneously Catalysed Reactive Distillation of Fuel Ethers MTBE and TAME. Chem. Eng. Sci. 2002, 52, 993. (59) Seader, J. D.; Henley, E. J. Separation Process Principles; John Wiley and Sons, Inc.: New York, 1998. (60) Bansal, V.; Perkins, J.; Pistikopoulos, E.; Ross, R.; van Schijndel, J. Simultaneous Design and Control Optimisation under Uncertainty. Comput. Chem. Eng. 2000, 24, 261. (61) Process Systems Enterprise. gPROMS AdVanced User Guide; Process Systems Enterprise: London, 1999. (62) Lim, Y. I.; Floquer, P.; Joulia, X.; Kim, S. D. Multiobjective Optimization in Terms of Economics and Potential Environment Impact for Process Design and Analysis in a Chemical Process Simulator. Ind. Eng. Chem. Res. 1999, 38, 4729. (63) Jimenez, E.; Salamon, P.; Rivero, R.; Rendon, C.; Hoffman, K.; Schaller, M.; Andresen, B. Optimization of a Diabatic Column with Sequential Heat Exchangers. Ind. Eng. Chem. Res. 2004, 43, 7566. (64) Larsson, M.; Wang, C.; Dahl, J. Development of a Method for Analysing Energy, Environmental and Economic Efficiency for an Integrated Steel Plant. Appl. Therm. Eng. 2006, 26, 1353. (65) Tsatsaronis, G. Definition and Nomenclature in Exergy Analysis and Exergoeconomics. Energy 2007, 32, 249. (66) Almeida-Rivera, C. P.; Swinkels, P. L. J.; Grievink, J. Economics and Exergy Efficiency in the Conceptual Design of Reactive Distillation Processes. AIChE Symp. Ser. 2004, 2004, 237. (67) Huang, K.; Iwakabe, K.; Nakaiwa, M.; Tsutsumi, A. Towards Further Internal Heat Integration of Reactive Distillation Columns. Part I: The Design Principle. Chem. Eng. Sci. 2005, 60, 4901. (68) Nummedal, L.; Kjelstrup, S. Equipartition of Forces as a Lower Bound on the Entropy Production in Heat Exchange. Int. J. Heat Mass Transfer 2001, 44, 2827. (69) Sauar, E.; Siragusa, G.; Andresen, B. Equal Thermodynamic Distance and Equipartition of Forces Principles Applied to Binary Distillation. J. Phys. Chem. A 2001, 105, 2312. (70) Georgiadis, M.; Schenk, M.; Pistikopoulos, E.; Gani, R. The Interactions of Design, Control and Operability in Reactive Distillation Systems. Comput. Chem. Eng. 2002, 26, 735. (71) Wang, S. J.; Wong, D. S. H.; Lee, E. K. Effect of Interaction Multiplicity on Control System Design for a MTBE Reactive Distillation Column. J. Proc. Control 2003, 13, 503. (72) Sneesby, M. G.; Tade, M. O.; Smith, T. N. Steady-state Transitions in the Reactive Distillation of MTBE. Comput. Chem. Eng. 1998, 22, 879. (73) Hasse, R. Thermodynamics of IrreVersible Processes; AddisonWesley Publishing Co.: Darmstadt, Germany, 1969. (74) Ung, S.; Doherty, M. Vapor-Liquid Phase Equilibrium in Systems with Multiple Chemical Reactors. Chem. Eng. Sci. 1995, 50, 23.
ReceiVed for reView November 28, 2006 ReVised manuscript receiVed September 10, 2007 Accepted September 11, 2007 IE061521J