Membrane Parallel Hybrid: A

Jan 20, 2010 - of olefins from paraffins, namely the retrofit of a C3-splitter and the retrofit of a C2-splitter through parallel hybridization with f...
4 downloads 0 Views 252KB Size
Ind. Eng. Chem. Res. 2010, 49, 2295–2305

2295

Energy Requirement of a Distillation/Membrane Parallel Hybrid: A Thermodynamic Approach Etienne Ayotte-Sauve´, Mikhail Sorin,* and Fernand Rheault CanmetENERGY, Natural Resources Canada, 1615 Lionel-Boulet, P.O. Box 4800, Varennes, Quebec, J3X 1S6, Canada

This paper presents a new thermodynamic approach, based on the notion of power of separation, for the retrofit problem of finding the minimal energy requirement of an existing binary distillation column when coupling it in parallel with a membrane unit. A new geometric interpretation of this concept, which is supported by a rigorous mathematical proof, is introduced. From it results an efficient and accurate shortcut method to tackle the aforementioned problem. Numerical examples are considered for the energy intensive separation of olefins from paraffins, namely the retrofit of a C3-splitter and the retrofit of a C2-splitter through parallel hybridization with facilitated transport membranes. The results of the proposed shortcut method are compared to those obtained via the nonlinear programming optimization solvers GAMS-CONOPT and GAMS-CoinIpopt for a superstructure based problem formulation, which act as a reference. In both case studies, the shortcut method results in a significant reduction in problem size and in the number of solver iterations, while yielding only a small error on the minimal energy requirement of the column within the hybrid and on the hybrid architecture (i.e., the position of side-streams along the column). It could therefore be used to carry out a rapid screening of alternatives (e.g., different membrane technologies) in order to evaluate potential energy improvements. The proposed approach could also provide mathematical programming algorithms with a good initial guess of the solution. 1. Introduction In recent years, hybrid processes composed of a distillation column coupled with another separation unit (membrane, pressure-swing adsorption, etc.) have received increased attention by both the industry and scientific community.1–7 For the separation of binary zeotropic close-boiling mixtures (e.g., propane-propylene or ethane-ethylene), distillation alone often requires a large number of trays (sometimes over 240 trays6). This number can be reduced, often by a significant percentage and with a noteworthy reduction in total annualized cost,6 by coupling the column with a membrane unit. However, many of today’s industrial sites already carry out the separation of such mixtures with distillation alone. Such columns often have high energy requirements (e.g., ethane-ethylene cryogenic distillation). As a consequence, the engineer is faced with the retrofit task of augmenting an existing distillation column with another separation unit in order to diminish the column’s energy consumption. In principle, the approach proposed here applies to any type of additional separation unit. Nevertheless, the scope is restricted to the case where this additional unit is a membrane. Hybrid processes involving a distillation column and a membrane can be realized in many different ways depending on the properties of the mixture to be separated, the process targets, etc. For example, the separation of binary azeotropic mixtures can be carried out by a serial arrangement where the membrane is located on the top (or bottom) of the distillation column in order to bypass the azeotropic point.5 For the separation of binary zeotropic close-boiling mixtures (e.g., propane-propylene or ethane-ethylene), one particular arrangement, termed “parallel”, is recommended8 and was investigated by many authors (e.g., Moganti et al.,3 Stephan et al.,4 and Kookos6). A parallel hybrid process involving a distillation column and a membrane unit * Corresponding author: Tel.: 450-652-3513. Fax: 450-652-5918. E-mail: [email protected]. 10.1021/ie9009155

consists in taking the membrane unit’s feed as a side-stream drawn from the column and feeding its products back into the column. This paper considers the problem of determining the minimal energy requirement attainable by an existing binary distillation column while coupling it in parallel with a membrane. Note that the energy consumption associated with the membrane’s compressors is generally small when compared to that of the distillation column. For example, Kookos6 presents a case study where the electricity consumption for the membrane’s compressor is roughly 1.1% of the energy associated with the steam used by the distillation column. Therefore, the column’s energy consumption can be viewed as an approximation of the energy requirements for the whole hybrid process. Of all binary mixtures to be treated by a distillation/membrane parallel hybrid, close-boiling zeotropic mixtures have arguably received the most attention from both the scientific and industrial communities thus far. Given that a column separating such a mixture involves many trays, there exists a large number of possible parallel hybrid architectures. For example, a column with 140 trays generates more than 5 × 107 possible parallel design arrangements with a single membrane unit. As was remarked by Kookos6 and other authors, the combinatorial nature of the problem and the nonconvexities present in the model equations may give rise to varying results when it is formulated as a mathematical programming problem and solved using a gradient based algorithm (e.g., the general reduced gradient method in GAMS-CONOPT9) supplied with various initial guesses of the hybrid’s architecture and operating conditions. This paper presents a new thermodynamically based approach to tackle the aforementioned problem. It is based on the concept of power of separation introduced in the works of Sorin et al.10 and Sorin and Rheault11 and later modified in the work of Ayotte-Sauve´ et al.12 and Sorin et al.13 The new approach has its roots in a new geometric interpretation of the power of

Published 2010 by the American Chemical Society Published on Web 01/20/2010

2296

Ind. Eng. Chem. Res., Vol. 49, No. 5, 2010

separation of a column section which will be presented in section 2. Not only does this geometric result yield new insights into the thermodynamic limitations imposed on both processes, but it also gives rise to a shortcut method which accurately predicts both the column’s minimal energy consumption within the hybrid and the optimal position of the input and output streams of the membrane along the column as well as the column’s feed position (which is allowed to vary). In the next section, the definition and basic properties of the power of separation are introduced. A mathematical expression linking the power of separation of each column section and the areas of the regions between the column operating lines and the diagonal in the McCabe-Thiele diagram is found. A shortcut method for the retrofit design of parallel hybrid processes separating a binary mixture is defined (section 3.2). The results obtained via this shortcut method and via a mathematical programming problem formulation (section 3.1), which involves a superstructure of all possible hybrid architectures and is termed the superstructure approach for short, are presented in section 3.4 for two case studies: the retrofit of a C3-splitter (propanepropylene) and the retrofit of a C2-splitter (ethane-ethylene) through parallel hybridization with facilitated transport membranes. This comparative study aims at showing that the shortcut method is both efficient and accurate with respect to the superstructure approach and can therefore be used to carry out a rapid screening of different scenarios (e.g., different membrane technologies) and to supply mathematical programming approaches with a good initial guess. 2. Power of Separation: Definition and Properties 2.1. Definition. In this section, the notion of power of separation is defined and its relation with the concept of compositional exergy is introduced. Consider a continuous steady-state separation process (involving no chemical reactions) with input streams Inj (j ) 1, ..., p) and output streams Outk (k ) 1, ..., q) separating a mixture with c components (i ) 1, ...,c). The power of separation of the process (Psep) is defined as Psep ) RgT0[

∑ Out z

k i,k

ln(zi,k)-

i,k

∑ In z

j i,j

ln(zi,j)]

(1)

i,j

where Rg is the ideal gas constant; T0 is the temperature of the surroundings; Inj (respectively Outk) is the total molar flow rate of input j (respectively, output k); zi,j (respectively zi,k) is the molar fraction of component i in the stream Inj (respectively, Outk). The notion of power of separation for component i (Psep,i) is defined mathematically as the term in the preceding equation corresponding to component i: Psep,i ) RT0[

∑ Out z

k i,k

ln(zi,k) -

k

∑ In z

j i,j

ln(zi,j)]

(2)

j

The power of separation of a process can be viewed as part of its compositional exergy. To see this, recall that by comparing the process to a reversible process carrying out the same separation, the first and second laws of thermodynamics give an expression for the compositional exergy (Ex) of the process14,15 Ex ) RT0[

∑ i,k

Outkzi,k ln(γi,kzi,k)-



Injzi,j ln(γi,jzi,j)]

i,j

(3) where γi,j (respectively, γi,k) is the activity coefficient corresponding to component i in the stream Inj (respectively, Outk),

calculated at temperature T0. Some authors refer to Ex as the minimal (rate of) work of separation of the process.16 Using basic properties of the logarithm, the right-hand side of the preceding equation can be split into two parts: Ex ) Psep + RT0[

∑ Out z

k i,k

ln(γi,k)-

i,k

∑ In z

j i,j

ln(γi,j)]

i,j

(4) Therefore, the compositional exergy of the process can be written as the sum of two terms: the power of separation of the process and the excess (rate of) change in Gibbs free energy between the input and output streams. Note that for an ideal mixture, i.e. a mixture for which all corresponding activity coefficients are equal to one, Psep ) Ex. An application of the concept of power of separation for multicomponent zeotropic and azeotropic distillation is introduced in the work of AyotteSauve´ et al.12 and further examples are given in the work of Sorin et al.13 The next section describes the properties of the power of separation in the context of binary distillation and introduces new geometric results. 2.2. Properties. Additivity. Consider a (steady state) process (A) composed of two (steady state) subprocesses (B and C) such that the output of each subprocess is either sent as an input for the other subprocess or sent to the environment. It follows from the definition of Psep that the power of separation of A (Psep,A) is the sum of the power of separation of B (Psep,B) and C (Psep,C). In other terms, the power of separation is an additive function. Note that we assume that a stream going from B to C exits B with the same composition and flowrate as when it enters C (and vice versa), but not necessarily with the same temperature or pressure. Geometric Properties. In the present paper, the geometrical results of Ayotte-Sauve´ et al.12 are further developed in the context of the separation of a binary mixture by a distillation column. By transforming the x-y coordinates of the traditional McCabe-Thiele diagram for such a column, a link naturally arises between the power of separation of each column section and the areas of the regions between the corresponding operating lines and the diagonal. Remark. The following result is independent of the nature of the column’s trays (i.e., it is valid for equilibrium trays or for trays with a Murphree efficiency different from unity), it is independent of the column’s complexity (i.e., it is valid for a column with or without side-streams) and is independent of the mixture assumptions (i.e., it is valid for ideal, nonideal, zeotropic and azeotropic binary mixtures). Consider a section of a distillation column, i.e. a maximal set of consecutive trays for which the inlet and outlet total liquid section denote the (respectively vapor) flow rates are equal. Let Psep,i power of separation of the section (for component i of the mixture) and let A denote the area between the image of the operating line and the image of the diagonal by the transformation (xi,yi) |f (xi,ln(yi)). One can then show that the dimensionless power of separation (for component i) is equal to the area A in the modified McCabe-Thiele diagram (with coordinates xi-ln(yi)), that is section Psep,i )A RgT0L

(5)

where L is the liquid flow rate in the section. A mathematical proof of this result is given in Appendix A. Figure 1 gives a visual representation of these results for the case of a singlefeed binary distillation column separating a hexane-p-xylene

Ind. Eng. Chem. Res., Vol. 49, No. 5, 2010

2297

Figure 1. Power of separation (of the light component) for a single-feed binary column.

mixture (example 3.1 from the work of Doherty and Malone17). A direct corollary of this result and the additivity property of the power of separation is that the power of separation of the whole column (for component i) is equal to the sum of the areas (As) in the xi-ln(yi) diagram corresponding to each section (s) weighed by the liquid flow rate in each section (Ls) and the reboiler and condenser’s power of separation (for component i). These results gives a link between the area (A) associated to a column section and its power of separation. Note that an increase in the area (A) would result in an increase in the number of distillation trays in the section.3 As a consequence, the previous results give an intuitive link between the number of trays in a column section and its power of separation. 3. Optimization of a Distillation Column/Membrane Parallel Hybrid Separating a Binary Mixture The problem considered in the present section is the minimization of the energy consumption of an existing (adiabatic) binary distillation column by coupling it in parallel with a membrane. Recall that a parallel hybrid process involving a distillation column and another separation unit consists in taking the additional separation unit’s feed as a side-stream drawn from the column and feeding its products back into the column. To be more precise, the following problem is treated: For a parallel hybrid configuration separating a binary mixture with a fixed number of distillation trays (N0) and a giVen separation task, that is, giVen F, xF, q, xB, and xD, find the membrane surface area, the position of the membrane’s feed and product streams along the column, the position of the column feed, and the hybrid profile (i.e., the flow rates and composition of all streams) such that the column’s reflux ratio is minimal. The preceding problem formulation assumes that the pressures at both sides of the membrane are given. It also assumes that upper bounds on the membrane surface area and on the membrane’s input flow rate are given. An upper bound on the membrane surface area can either be defined via investment cost constraints or by other considerations, such as the variation of the membrane’s selectivity with respect to its surface area.3 On the other hand, an upper bound on the input flow rate of the membrane can be obtained either by considering the variation of its separation efficiency with respect to its input flow rate3 or by considering other constraints, such as hydrodynamic limits for the trays of the existing distillation column. For example, if too much liquid is drawn from the column to feed the membrane,

Figure 2. Superstructure for parallel hybrid configurations.

then the reflux ratio of the column may be reduced lower than a threshold value, under which blowing occurs.18 To tackle the optimization problem, an approach based on the geometric interpretation of the power of separation presented in the previous section is introduced. In order to show that the proposed shortcut method is efficient and predicts accurately the optimal reflux ratio and an optimal hybrid configuration, its results are compared to those of a mathematical programming problem formulation obtained by writing the model equations for the column, here chosen as the MESH equations (see, for example, the work of Henley and Seader16), and the membrane model equations in the context of a superstructure representing all hybrid configurations considered. The latter problem formulation defines a nonlinear programming problem (NLP). Note that the superstructure presented here is a modification to the retrofit case of the superstructure proposed by Kookos.6 The following two sections present the superstructure problem formulation and the shortcut method, respectively. 3.1. Superstructure Problem Formulation. In order to formulate the mathematical programming problem, the superstructure presented in Figure 2 is used. The distillation trays are numbered from bottom to top j ) 1, ..., N, including the condenser and reboiler. Figure 3 gives a summary of the notation. Since the column feed position is allowed to vary, the feed stream is split into N - 2 streams and Fj denotes the substream entering tray j. It is assumed that part of the liquid stream flowing out of tray j (Mj) is drawn out of the column to be eventually fed to the membrane. All streams Mj are mixed, forming a stream M, which is totally evaporated and fed to the membrane. The permeate and retentate products of the membrane are totally condensed, split, and then fed on each column tray (streams Pj and Rj, respectively). Even if the problem at hand involves a binary mixture, the distillation model will be presented for a c-component mixture

2298

Ind. Eng. Chem. Res., Vol. 49, No. 5, 2010 N-1

M)

∑M

j

(16)

j)2

N-1

MyM )

∑Mx

j j

(17)

j)2

N-1

P)

∑P

j

(18)

∑R

(19)

j)2

N-1

Figure 3. Distillation tray within the superstructure.

R)

(i ) 1, ..., c). The column is modeled using the MESH equations. These equations involve equilibrium relations coupled with material and energy balances on each tray. The internal pressure of the column is assumed to be constant, and all trays are taken to be at equilibrium. The equations modeling the column inside the hybrid are given below for a partial condenser and a partial reboiler. These relations can be easily modified for the case of a column equipped with a total condenser, a total reboiler or vapor side-streams. Equilibrium Relations. yi,j ) Ki,jxi,j,

for i ) 1, ..., c and j ) 1, ..., N

(6)

Mole Fraction Summations.

∑x )∑y i,j

i

i,j)1,

for j ) 1, ..., N

(7)

i

Material Balances. Condenser. VN-1yN-1 ) LNxN + VNyN

(8)

Column Tray. FjxF + Vj-1yj-1 + PjyP + RjyR + Lj+1xj+1 ) (Lj + Mj)xj + Vjyj, for j ) 2, ..., N - 1 (9) Reboiler. L2x2 ) V1y1 + L1x1

(10)

Energy Balances. Condenser. V VN-1HN-1 ) LNHLN + VNHVN + Qc

(11)

Column Tray. V L FjHF + Vj-1Hj-1 + PjHLP + RjHLR + Lj+1Hj+1 ) (Lj + Mj)HLj + VjHVj , for j ) 2, ..., N - 1

(12)

Reboiler. L2HL2 + Qr ) V1HV1 + L1HL1

(13)

Definition of Reflux Ratio. rVN ) LN

(14)

Mixer and Splitter Equations. N-1

F)

∑F

j

j)2

(15)

j

j)2

For both case studies treated in this paper, facilitated transport membranes with silver ion as the carrier were considered. A brief model description and model equations for each case study are presented in Appendix B. The hybrid model is complete when the vapor-liquid equilibrium (VLE) relations defining the K-values and molar enthalpies are added. The optimization problem can then be stated: Minimize r subject to eqs 6-19, the membrane model equations (see Appendix B), the K-value and molar enthalpy definitions and bound constraints on the variables. 3.2. Shortcut Method: Basic Principles. Before stating the algorithmic steps involved in the shortcut method for the retrofit optimization problem formulated in section 3.1, the basis of the procedure is explained intuitively in the present section. It involves the interplay between the retrofit problem and its “dual” problem, i.e. the design problem, and their interpretation via the geometric results linking the power of separation of column sections and areas in the McCabe-Thiele diagram (section 2.2). In particular, the retrofit problem can be reduced to finding the reflux ratio for which the minimal number of trays for the design problem corresponds to the actual number of trays in the existing column. The following paragraphs give precise meaning to the previous statements. Note that as in section 2.2, it is assumed that the streams F, P, and R enter the column with no exergy losses due to mixing, i.e. F (P and R, respectively) enters the column where the composition of the internal liquid stream is equal to the composition of F (P and R, respectively). 3.2.1. Dual Problem. The dual problem to the retrofit problem (section 3.1) consists in minimizing the number of distillation trays within a distillation/membrane parallel hybrid for a given value of the reflux ratio and for a given separation task. In light of the geometric results of section 2.2, consider the areas between the various operating lines and the diagonal in the McCabe-Thiele diagram of such a process (Figure 4). Note that the operating lines for each column section within the hybrid are in bold, while the operating lines of a column alone (i.e., without a membrane)scarrying out the same separation as the hybrid with the same reflux ratiosare the bold lines in sections 1 and 5 extended by the dotted lines in the remaining sections (2-4). The area associated with the column within the hybrid (in gray) corresponds to the column’s power of separation within the hybrid. On the other hand, the area of the triangle whose boundary is defined by the column’s operating lines (no membrane) and the diagonal corresponds to the power of separation of a column alone carrying out the same separation task as the hybrid (and with the same reflux ratio). Since the column alone and the hybrid have the same input and output streams (F, B, and D), their power of separation is the same. It follows that the area of the triangle corresponds to the power of separation of the hybrid. The difference between the area

Ind. Eng. Chem. Res., Vol. 49, No. 5, 2010

2299

Figure 4. Example of a McCabe-Thiele diagram for a parallel hybrid separating an ideal binary mixture.

associated with the hybrid and the area associated with the column within the hybrid, i.e. the white area between the dotted lines and the column operating lines in sections 2-4, corresponds to the membrane’s power of separation. The additivity of these areas translates into the additivity of the power of hybrid column membrane column ) Psep + Psep , where Psep is the separation, i.e. Psep power of separation of the column within the hybrid. The preceding intuitive geometric interpretation generates a heuristic which characterizes the minimal number of trays for a distillation/membrane parallel hybrid with fixed separation task and fixed reflux ratio. The number of distillation trays in a parallel hybrid (separating a binary mixture) is linked with the area between the operating line and the diagonal in the McCabe-Thiele diagram of the process (i.e., the gray area in Figure 4). In fact, minimizing this area gives a good approximation of the minimal number of distillation trays necessary to carry out the separation task, a fact which was investigated and validated by Moganti et al.3 and Stephan et al.4 Since this area corresponds to Pcolumn , one could be led to think that minimizing sep may minimize the number of distillation trays (for a given Pcolumn sep separation task and a given reflux ratio). The inlet and outlet streams of the hybrid having fixed flow rates and compositions, the power of separation of the hybrid is fixed. It then follows from the additivity of the power of separation that minimizing column membrane is equivalent to maximizing Psep . Psep One also notes that at minimum reflux, a binary single-feed (adiabatic) distillation column alone either presents a feed pinch or a tangent pinch (at composition x ) xpinch), for which a neighborhood contains an infinite number of trays. In order to minimize the number of distillation trays by adding a membrane, it could then seem natural to ask that the membrane be located at that pinch point. The preceding arguments thus give rise to the following heuristic: Heuristic. Consider a distillation column/membrane parallel hybrid separating a binary mixture, with F, xF, xB, xD, and r fixed. A “good” approximation of the minimal number of trays membrane for in the distillation column is given by maximizing Psep a value of xM corresponding to the pinch-point of the distillation column. 3.2.2. Retrofit Problem. Adding a membrane to an existing column (with a total of N0 trays) changes its McCabe-Thiele diagram (Figure 4). Since the area between the operating lines and the diagonal is reduced (see gray area in Figure 4), the number of distillation trays calculated via the staircase

McCabe-Thiele construction with the new operating lines (and the same reflux ratio as the column alone) is smaller than N0. In order to meet the requirement of a fixed number of trays, this area must be changed by specifying new values for the degrees of freedom of the hybrid system, which in the current context can be chosen as: the column’s reflux ratio and the membrane characteristics (e.g., surface area, input flow rate, and input composition). Remember that since it is assumed that the position of column side-streams is determined by their composition (i.e., no exergy losses due to the mixing of side-streams with internal column streams), no additional degrees of freedom are required. Aiming at minimizing the column’s reflux ratio, choosing the hybrid’s degrees of freedom is a nontrivial task. The following paragraphs present a method to determine these values. Consider the retrofit problem of minimizing an existing column’s reflux ratio (with a total of N0 fixed trays and a fixed separation task) by augmenting the column with a membrane unit located in parallel (section 3.1). As a first step, the membrane’s characteristics can be determined according to the heuristic presented in section 3.2.1, that is by maximizing its power of separation (with input composition corresponding to the column’s pinch point). Afterward, there remains only one degree of freedom to determine: the column’s reflux ratio. As was previously hinted, this can be achieved by exploiting the link between the retrofit problem and the design problem. To be more precise, the minimal reflux ratio (rmin) for the retrofit problem (with a total of N0 fixed trays, section 3.1) can be obtained by finding the reflux ratio (r) for which the minimal number of trays (Nmin ) Nmin(r)) obtained for the design problem (section 3.2.1) verifies the equality Nmin ) N0. Note that a sufficient condition for the previous statement to be true is that the function Nmin ) Nmin(r) is strictly decreasing (see Figure 6 of section 3.4). In that case, one has rmin ) r. In the following section, an algorithm based on the heuristic presented in the section 3.2.1 to solve the design problem is presented. Finding the value of r for which Nmin(r) ) N0 can be achieved with an appropriate root finding algorithm (e.g., Newton’s method). 3.3. Shortcut Method: Algorithm Steps. The present section introduces a shortcut method based on the notion of power of separation to tackle the problem presented in section 3.1. 1. Determination of the Side-Stream Characteristics (M, P, R, yM, yP, yR). In order to determine the values of the

2300

Ind. Eng. Chem. Res., Vol. 49, No. 5, 2010

unknowns M, P, R, yM, yP, yR, and S, solve the following optimization problem: membrane maximize Psep ) RgT0[P(yP ln(yP) + (1 - yP) ln(1 yP)) + R(yR ln(yR) + (1 - yR) ln(1 - yR)) - M(yM ln(yM) + (1 - yM) ln(1 - yM))] (20)

subject to: (a) the membrane model equations (Appendix B); (b) yM ) xpinch; (c) 0 e M e MU and 0 e S e SU, where MU and SU are upper bounds on the membrane feed flow rate and surface area respectively (see section 3.4) and xpinch is the column’s pinch composition at minimal reflux. In the remaining steps of the algorithm, the side-stream characteristics M, P, R, yM, yP, and yR assume the values obtained in step 1. 2. Initialization of the Reflux Ratio (r). Initialize the column’s reflux ratio (r). For example, the value of the base case reflux ratio for the column alone could be taken. 3. Calculation of the Number of Trays (N(r)). Calculate the number of distillation trays N ) N(r) for the given reflux ratio r and the membrane characteristics obtained in step 1. 4. Updating the Reflux Ratio (r). If |N - N0| < tolerance, where N is the number of trays calculated in step 3 and N0 is the number of fixed distillation trays, then the procedure stops and the current reflux ratio value is chosen as an approximation of rmin (see Figure 6). If |N - N0| g tolerance, then change the value of r according to a root finding algorithm (e.g., Newton’s algorithm) and go to step 3. In the above procedure, the power of separation of the membrane is maximized only once (step 1). For the membrane specifications obtained in step 1 and a given reflux ratio, the number of distillation trays is calculated (step 3), for example by an iterative procedure assuming the column verifies the constant molar overflow (CMO) hypothesis and that the streams F, P, and R enter the column with no exergy losses due to mixing, i.e. F (P and R, respectively) enters the column where the composition of the internal liquid stream is equal to the composition of F (P and R, respectively). If the CMO hypothesis is not verified, the Ponchon-Savarit procedure18 can be used to calculate the total number of trays instead of the McCabe-Thiele approach. This step can be naturally adapted to the case of vapor side-streams. The resulting number of trays (N) must be equal to the fixed number of trays in the column (N0). If equality is not met up to a specified tolerance, then the reflux ratio value is changed (step 4) and the number of trays is recalculated. In the next section, results of the shortcut method are compared to that of the superstructure formulation, which acts as a reference, based on two case studies for the separation of olefins from paraffins: the retrofit of a C3-splitter and the retrofit of a C2-splitter through parallel hybridization with facilitated transport membranes. 3.4. Numerical Results. Case Study 1: Retrofit of a C3-Splitter. The following case study was taken from the work of Moganti et al.3 and Stephan et al.4 It consists in the separation of a mixture of propane and propylene by a distillation column/ facilitated transport membrane parallel hybrid to produce high (>99%) purity propylene. For such a separation task, Gottschlich and Roberts19 showed that a hybrid process can reach a higher thermodynamic efficiency and a lower processing cost than either a distillation column or a membrane unit taken alone. The liquid and vapor phase of the mixture are assumed ideal and the VLE is described by the Soave-Redlich-Kwong (SRK) model.20 The column’s internal pressure is assumed constant (17.21 bar); the condenser and reboiler are partial and equilibrium trays are considered. Table 1 gives the column base case

Table 1. C3-Splitter Base Case Data data

value

F (saturated liquid) xF yD xB r N0 VLE model

2.78 (mol/s) 0.44 0.99 0.04 18 137 SRK

data used in the calculations. It is assumed that the membrane’s permeate is recompressed to the internal pressure of the column and the membrane’s feed is drawn from the column as a saturated liquid to be totally evaporated before being fed in the membrane. The saturated vapor products of the membrane are then totally condensed before entering the column. An upper bound of 90 m2 was taken for the membrane surface area based on a study of the membrane’s selectivity as a function of its surface area.3 An upper bound of 2 mol/s was chosen for the membrane’s input flow rate. This bound corresponds to roughly 10% of the liquid flow rate inside the rectifying section of the column operating at its base case regime. More precise hydrodynamic constraints could be taken into account in order to determine this upper bound;18 such a task is outside the scope of this text. The superstructure formulation given in section 3.1 was treated in the GAMS environment with the use of the CONOPT solver, which implements a generalized reduced gradient algorithm suited for large-scale problems.9 On the other hand, for the shortcut method, step 1 (section 3.2) was tackled with GAMS-CONOPT while step 4 was implemented with Newton’s algorithm. Figure 5 shows the results of the shortcut method and of the superstructure approach. Note that the GAMSCONOPT algorithm was run with various initial points for the superstructure approach (see below for details) and the best result obtained is shown in Figure 5a. The shortcut method’s results were independent of the initial guess given to GAMSCONOPT (section 3.3, step 1 of the algorithm). While examining the results, one notes that the minimal reflux ratio estimation by the shortcut method is only 2.4% higher than the one obtained via the superstructure formulation, the latter corresponds to a reduction of 29.56% of the column’s reflux ratio with respect to the base case. The authors believe the difference in the results of both methods can be attributed, in part, to a combination of the following factors: the shortcut method assumes that side-streams and the corresponding column internal streams have the same composition, as a consequence nonintegral tray numbers in each column section are obtained for the shortcut method (these were rounded to the nearest integer in Figure 5b); the heuristic given in section 3.2.1 offers an approximation of optimality, not an exact characterizations partly because it may not be possible to have a tray with liquid output composition matching exactly the pinch composition. Not only do both methods give similar estimates of the minimal reflux ratio, there is also a noticeable similarity between the hybrid architectures (i.e., the side-stream positions along the column). In particular, while no conditions were imposed on the molar fraction of propylene in the membrane’s feed in the superstructure approach, the results show that the stream M should be located near the feed tray in order to minimize the column’s reflux ratio. In both cases, the membrane’s feed flow rate reaches its upper bound. In terms of the power of separation of the membrane, membrane ) 909 the superstructure approach yields a value of Psep membrane ) 913 W while the shortcut approach gives a value of Psep

Ind. Eng. Chem. Res., Vol. 49, No. 5, 2010

2301

Figure 5. Optimization results for the retrofit of a C3-splitter: (a) reference case (best result for superstructure approach) and (b) shortcut method.

W. These result, coupled with the fact that the superstructure approach gives roughly a value of yM equal to xF, confirms the intuition (section 3.2) that the membrane’s power of separation should be near its maximum when the column’s reflux ratio is minimal (for a given separation task). Consider both resolution methods. The superstructure formulation necessitates an initial guess for the hybrid configuration and operating conditions. To be specific, in the current case study, initial guesses consisted in specifying values for the following: the reflux ratio (and hence liquid and vapor flow rates inside the column); the position of the column’s feed and that of the membrane’s feed, permeate, and retentate streams; the flow rates of the membrane feed, permeate, and retentate streams; rough estimates of the liquid and vapor composition profiles in the column (generated randomly via a uniform probability distribution in the interval [0.85, 0.95]). If we assume that the initial hybrid configuration is such that the column feed, the membrane’s feed, permeate, and retentate streams are not allowed to be split into substreams (or mix; see Figure 2), then there are over 5 × 107 such possible architectures encoded in the superstructure (Figure 2)sfor the present case study. In general, if side-stream splitting (or mixing) is not allowed and N denotes the number of distillation trays (including condenser and reboiler), then there are (N - 2)(N - 3)2(N 4)/6 such configurations, assuming the permeate (respectively retentate) is located strictly above (respectively below) the membrane’s feed along the column. Therefore, the number of possible choices for initial hybrid architecture in the present case study is large enough to warrant multiple trials. The results obtained via GAMS-CONOPT for the superstructure formulation (section 3.1) varied with the choice of the initial hybrid configuration. In other terms, for given values of the initial estimates for the reflux ratio, the membrane input and output flow rates, and the liquid and vapor composition profiles, different results were obtained when the position of the column feed or the membrane feed, permeate, or retentate streams were changed in the initial guess. Only the best resultsin terms of the column’s reflux ratiosis shown in Figure 5a. On the other

Table 2. Efficiency of the Shortcut Method for the Retrofit of a C3-Splitter

number of equations number of variables number of iterations

short-cut method

superstructure approach

5 8 27 CONOPT iterations + 19 Newton iterations

558 1097 575 CONOPT iterations

hand, such considerations do not come into account when considering the shortcut method, since the hybrid architecture is imposed as being that for which the membrane’s power of separation is maximal (and for which the membrane’s feed corresponds to the column’s feed pinch). In terms of efficiency, both approaches are compared in Table 2. The differences in problem size and number of iterations of the optimization algorithm show that the shortcut method is indeed efficient. This section concludes by showing that the power of separation approach compares favorably with respect to the “minimum area method”.3 Figure 6 shows results of both methods for the design problem (section 3.2.1) for various reflux ratios. The approximation of the minimal number of trays by the power of separation method is better than that offered by the minimum area method. To be more precise, the mean difference in Nmin estimation between both methods is about six trays. The minimum area method therefore overestimates the minimum reflux ratio for the retrofit problem by 7.6% (see Figure 6 and Table 3), while the power of separation approach gives an error of 2.4%. The hybrid architecture given by the power of separation shortcut method predicts with an error of at most four trays the position of all side-streams corresponding to the input and outputs of the membrane along the column, while the minimal area method gives a mean error of 12 trays (with errors ranging between 8 and 20 trays). Remark. For propane-propylene at 17.21 bar and according to the SRK model, the molar heats of vaporization of the mixture components are roughly the same (their relative difference is smaller than 0.8%). Therefore, the CMO assumption was

2302

Ind. Eng. Chem. Res., Vol. 49, No. 5, 2010

Figure 6. Design and retrofit problems for the power of separation method and the minimum area method (C3-splitter case study). Table 3. Power of Separation Method vs Minimum Area Method for the Retrofit of a C3-Splitter

method reference (superstructure approach) Psep method min area method

rmin approximation

maximal error on side-stream position (in number of trays)

12.6790

0

12.9822 13.6387

4 20

assumed to be valid and the following results for the shortcut method were obtained by calculating tray numbers (step 3, section 3.3) with the McCabe-Thiele iterative procedure. Nevertheless, the shortcut method was also implemented by taking into account energy balances thus discarding the CMO hypothesis and the McCabe-Thiele method in favor of the Ponchon-Savarit method. The minimal reflux ratios of the C3splitter inside the hybrid predicted by the shortcut method by assuming the CMO hypothesis and by discarding this assumption were the same up to two decimals. Since step 1 of the shortcut method does not depend on the way tray numbers are calculated, both versions of the shortcut method yielded the same membrane characteristics; tests showed for the present case study that these two methods yielded hybrid architectures differing by at most one tray. Step 4 of the shortcut method (section 3.3) also took under 30 Newton iterations to converge for both versions (using the base case reflux ratio as a starting value). In conclusion, for this case study there is no considerable gain in choosing the Ponchon-Savarit procedure over the McCabe-Thiele method to count tray numbers in each column section in step 3 of the shortcut method (section 3.3). Case Study 2: Retrofit of a C2-Splitter. In many petrochemical industries, the separation of a mixture of ethane and ethylene is carried out by cryogenic distillation. This process having high energy requirements, coupling it with a membrane offers a way to diminish its energy consumption. The stream data for the distillation column considered here was taken from King.18 The base case reflux ratio and the (fixed) number of equilibrium trays were obtained by multiplying the column’s minimum reflux ratio by a factor of roughly 1.12. The resulting C2-splitter base case data is given in Table 4. The liquid and

Table 4. C2-Splitter Base Case Data data

value

F (saturated liquid) xF yD xB r N0 VLE model

76.239 (mol/s) 0.40455 0.9982 0.016 5.9553 60 SRK

vapor phase of the mixture are assumed ideal and the VLE is described by the SRK model.20 The column’s internal pressure is assumed constant (20.03 bar); the condenser and reboiler are partial and equilibrium trays are considered. As for the previous case study, it is assumed that the membrane’s permeate is recompressed to the internal pressure of the column and the membrane’s feed is drawn from the column as a saturated liquid to be totally evaporated before being fed in the membrane. The saturated vapor products of the membrane are then totally condensed before entering the column. Since no a priori study of the membrane’s selectivity as a function of its surface area was available, a “large” upper bound (SU ) 105 m2) was chosen for preliminary optimization tests. As the results given below confirm, this bound was more than sufficient, i.e. the optimization results yielded a membrane surface area strictly smaller than SU. An upper bound of 30 mol/s was chosen for the membrane’s input flow rate. This bound corresponds to roughly 10% of the liquid flow rate inside the rectifying section of the column operating at its base case regime. Preliminary numerical tests for the superstructure formulation did not give satisfactory results with the version of GAMSCONOPT used in the previous case study. As a consequence, both the superstructure formulation (section 3.1) and step 1 of the shortcut method (section 3.2) were treated with the GAMSCoinIpopt solver, which produced better results. Remark. For ethane-ethylene at 20.03 bar and according to the SRK model, the molar heats of vaporization of the mixture components have a relative difference of roughly 7.8%. Assuming CMO, the minimal reflux ratio of the C2-splitter inside

Ind. Eng. Chem. Res., Vol. 49, No. 5, 2010

2303

Figure 7. Optimization results for the retrofit of a C2-splitter: (a) reference case and (b) shortcut method.

the hybrid obtained with the shortcut method is 5.64% below that obtained via the superstructure approach. The results shown below were obtained not by assuming CMO, but by taking into account the energy balances on each tray via the Ponchon-Savarit method. As will be seen, a considerable improvement is gained. Figure 7 shows the results of the shortcut method and of the superstructure approach. The GAMS-CoinIpopt algorithm was run with various initial points for the superstructure approach. All initial points tested yielded the hybrid shown in Figure 7a. The shortcut method’s results were also independent of the initial guess given to GAMS-CoinIpopt (section 3.3, step 1 of the algorithm). Note that the minimal reflux ratio estimation by the shortcut method is only 1.52% higher than the one obtained via the superstructure formulation, the latter which corresponds to a reduction of 9.35% of the column’s base case reflux ratio (i.e., without membrane). As in the previous case study, the shortcut method gives a good approximation of the best result obtained via the superstructure approach in terms of the objective function value (i.e., in terms of the value of the reflux ratio) and in terms of the position of the streams F, M, P, and R along the column. In both cases, the membrane’s feed flow rate reaches its upper bound. The power of separation of the membrane given by the shortcut method (respectively superstructure approach) is membrane membrane ) 22.9 kW (respectively Psep ) 24.7 kW). These Psep result, coupled with the fact that the superstructure approach gives a difference between yM and xF that is less than 0.04, confirms the intuition (section 3.2) that the membrane’s power of separation is near its maximum when the column’s reflux ratio is minimal (for a given separation task). As in the previous case study, the differences in problem size and number of iterations of the optimization algorithm between the shortcut method and the superstructure approach show that the shortcut method is indeed efficient (Table 5). The power of separation approach compares favorably with respect to the “minimum area method”3 for the current case study (Table 6). The minimum area method underestimates the rmin value given by the superstructure approach by 4.39% while

Table 5. Efficiency of the Shortcut Method for the Retrofit of a C2-Splitter

number of equations number of variables number of iterations

short-cut method

superstructure approach

5 8 39 CoinIpopt iterations + 25 Newton iterations

374 604 422 CoinIpopt iterations

Table 6. Power of Separation Method vs Minimum Area Method for the Retrofit of a C2-Splitter

method reference (superstructure approach) Psep method min area method

rmin approximation

maximal error on side-stream position (in number of trays)

5.3985

0

5.4805 5.1615

2 5

the power of separation shortcut method overestimates it by 1.52%. The new approach also is more accurate in its prediction of the hybrid architecture. 4. Conclusion An efficient and accurate shortcut method for determining the minimum reflux ratio and the side-stream positions for an existing binary distillation column when augmenting it with a membrane unit in parallel is presented. It is based on a new geometric result linking the power of separation (and thus the compositional exergy) of a column section with the area between the operating line and diagonal in the McCabe-Thiele diagram of the process. In order to show that the new shortcut method is indeed both efficient and accurate, it is compared to a superstructure approach for the energy intensive separation of olefins from paraffins. Two case studies are presented: the retrofit of a C3-splitter and the retrofit of a C2-splitter through parallel hybridization with facilitated transport membranes. For the C3-splitter (respectively C2-splitter) case study, the shortcut method gives an approximation of the minimal reflux

2304

Ind. Eng. Chem. Res., Vol. 49, No. 5, 2010

ratio of the column within the hybrid which is only 2.4% (respectively 1.52%) higher than the value obtained via the superstructure formulation and predicts with an error of at most four trays (respectively, two trays) the optimal position of all side-streams corresponding to the input and outputs of the membrane along the column as well as the column’s feed position. While the GAMS-CONOPT (respectively GAMSCoinIpopt) solver used to tackle the superstructure problem formulation for the C3-splitter (respectively C2-splitter) case study takes roughly 570 iterations (respectively 420 iterations) to convergesinvolving more than 550 equations (respectively 370 equations)sthe shortcut method takes only 27 CONOPT (respectively 39 CoinIpopt) iterations (involving only the four membrane model equations) and 19 (respectively 25) Newton iterations. The proposed approach can therefore be used to do a rapid and reliable screening of different alternatives (e.g., different membrane technologies) and can provide mathematical programming approaches with a good initial guess of the optimum. Since the new method applies in principle to parallel hybrids between a distillation column and any additional separation unit, future work includes applying it in the context where the additional unit is a pervaporation membrane or a pressure-swing adsorption process. Acknowledgment The authors would like to thank the Office of Energy Research and Development (OERD) of Natural Resources Canada for its funding and support of this study. They would also like to offer their gratitude to Forouzan Sadeghi who supplied the membrane model for the second case study. Appendix A: Geometric Properties of the Power of Separation This appendix contains an elementary mathematical proof of eq 5 (section 2.2). To find the link between the area (A) and the power of separation of the section for component i (Psection sep,i ), note that by definition one has section ) RgT0[Lxout ln(xout) + Vyout ln(yout) - Lxin ln(xin) Psep,i Vyin ln(yin)]

where x and y refer to liquid (respectively vapor) molar fractions of component i and the subscripts in and out refer to the corresponding inlet or outlet streams for the section. Note that the molar fractions in the liquid and vapor crossing at a given height in the section are related through the operating line y)

L x+b V

where L and V are the total liquid and vapor flow rates, and b is determined by the flow rates (of component i) in the sidestreams and product streams going through the mass balance envelope used to obtain the operating line. Transforming the coordinates in the McCabe-Thiele diagram from x-y to x-ln(y), the area between the operating line and the diagonal (y ) x) becomes A)



xin

xout

[ln( VL x + b) - ln(x)] dx

It will be assumed that the composition and flow rate of all streams into and out of the column are known. It follows that b is constant (i.e., b does not depend on x). By definition of a

column section, L and V are also constant. Using these observations and basic properties of the logarithm, one obtains

∫ ) ∫

A)

xin

xout xin xout

)

ln

( VL + bx )dx

ln (Lx + bV)dx -

[ ( (

)



xin

xout

ln (Vx)dx

(

)

Lxin + bV Lxout + bV 1 Lxin ln - Lxout ln + L Vxin Vxout bV + Lxin bV ln bV + Lxout

)

The operating line equation gives the following relations yin )

L x +b V out

yout )

L x +b V in

To conclude, these relations can then be used to rewrite the section . area in the x-ln(y) diagram as A ) [1/(RgT0L)]Psep,i Appendix B: Membrane Models This appendix presents the membrane models used for both case studies. For the separation of propane and propylene, a facilitated transport membrane with silver ion as the carrier is considered. Quoting Moganti et al.:3 “Here, propylene reacts with the silver ion to form a silver-propylene complex which diffuses across the membrane. Hence, the net transport of propylene is due to a combination of solute diffusion and the diffusion of the complex.” The membrane model studied was previously considered by Gottschlich and Roberts19 and later by Moganti et al.3 and Stephan et al.4 Instead of repeating the modeling hypothesis and the list of coefficients used in the equations (e.g., Henry constants, diffusivity, etc.), the reader is referred to the aforementioned articles. The two mass balances and two transport relations which constitute the membrane model are given below. Note that the feed-side pressure of the membrane (PR) is fixed at 17.21 bar and that the permeate-side pressure (PP) is fixed at 3.44 bar (see the aforementioned references and the first case study in section 3.4 of the present paper). The membrane feed, permeate, and retentate are all saturated vapor streams. In the following model equations, the symbol y denotes the molar fraction propylene. M)P+R MyM ) PyP + RyR

(

PyP ) S(PRyR - PPyP) 0.00195 + 0.254 (1 + 0.479PPyP)(1 + 0.479PRyR)

)

P(1 - yP) ) 0.000608S(PR(1 - yR) - PP(1 - yP)) For the separation of ethane and ethylene, a bulk flow liquid membrane (BFLM) with silver nitrate as a carrier was considered.21 In such a membrane, a carrier solution is supplied to the feed side and is able to permeate to other side of the membrane. The transferred carrier solution is then recycled back to the feed side. An idealized permeation model proposed by Teramoto et al.21 was used in order to predict permeances. In this model, a value of 4 × 10-5 m/s was considered for the volumetric flux of the carrier solution. The reader can find all

Ind. Eng. Chem. Res., Vol. 49, No. 5, 2010

physical properties of the membrane in the aforementioned paper. In the following model equations, the symbol y denotes the molar fraction of ethylene. Note that the feed-side pressure of the membrane (PR) is fixed at 20.03 bar and that the permeateside pressure (PP) is fixed at 3.4535 bar. M)P+R MyM ) PyP + RyR

(

PyP ) 2.16771 × 10-3S(PRyR - PPyP) 1 + 473.846 (1 + 0.308PPyP)(1 + 0.308PRyR) -5

)

P(1 - yP) ) 7.12 × 10 S(PR(1 - yR) - PP(1 - yP)) Nomenclature Symbols A ) area between the operating line and the diagonal in the x-ln(y) diagram B ) bottoms flow rate [mol/s] D ) distillate flow rate [mol/s] F ) column feed flow rate [mol/s] K ) vector of equilibrium ratios M ) membrane feed flow rate [mol/s] N ) total number of trays (including condenser and reboiler) P ) membrane permeate flow rate [mol/s] PP ) membrane’s permeate side pressure [bar] PR ) membrane’s feed side pressure [bar] Psep ) power of separation [J/s] q ) thermal quality of column feed R ) membrane retentate flow rate [mol/s] Rg ) ideal gas constant [J/(K mol)] r ) reflux ratio S ) membrane surface area [m2] x ) liquid phase molar fraction y ) vapor phase molar fraction Subscripts and Superscripts i ) referring to component i of the mixture in ) referring to the input of a section j ) referring to distillation tray j out ) referring to the output of a section AbbreViations CMO ) constant molar overflow NLP ) nonlinear programming SRK ) Soave-Redlich-Kwong VLE ) vapor-liquid equilibrium

Literature Cited (1) Davis, J. C.; Valus, R. J.; Eshraghi, R.; Velikoff, A. E. Facilitated transport membrane hybrid systems for olefin purification. Sep. Sci. Technol. 1993, 28, 463–476.

2305

(2) Ghosh, T. K.; Lin, H. D.; Hines, A. L. Hybrid adsorption-distillation process for separating propane and propylene. Ind. Eng. Chem. Res. 1993, 32, 2390–2399. (3) Moganti, S.; Noble, R. D.; Koval, C. A. Analysis of a membrane/ distillation column hybrid process. J. Membr. Sci. 1994, 93, 31–44. (4) Stephan, W.; Noble, R. D.; Koval, C. A. Design methodology for a membrane/distillation column hybrid process. J. Membr. Sci. 1995, 99, 259– 272. (5) Bausa, J.; Marquardt, W. Shortcut design methods for hybrid membrane/distillation processes for the separation of non-ideal multicomponent mixtures. Ind. Eng. Chem. Res. 2000, 39, 1658–1672. (6) Kookos, I. Optimal design of membrane/distillation column hybrid processes. Ind. Eng. Chem. Res. 2003, 42, 1731–1738. (7) Lucia, A.; Amale, A.; Taylor, R. Energy Efficient Hybrid Separation Processes. Ind. Eng. Chem. Res. 2006, 45, 8319–8328. (8) Kreis, P.; Go`rak, A. Process analysis of hybrid separation processes combination of distillation and pervaporation. Chem. Eng. Res. Des. 2006, 84 (A7), 595–600. (9) Brooke, A.; Kendrick, D.; Meeraus, A.; Raman, R. GAMS, A user’s guide; GAMS Development Corporation: WA, 1998. (10) Sorin, M.; Jedrzejak, S.; Bouchard, C. On maximum power of reverse osmosis separation processes. Desalination 2006, 190, 212–220. (11) Sorin, M.; Rheault, F. Thermodynamically guided intensification of separation processes. Appl. Therm. Eng. 2007, 27, 1191–1197. (12) Ayotte-Sauve´, E.; Sorin, M.; Rheault, F.; Sadeghi, F. Power of separation and energy requirements in distillation. 8th World Congress of Chemical Engineering Montreal; Montreal, Canada, Aug 23-27, 2009. (13) Sorin, M.; Ayotte-Sauve´, E.; Sadeghi, F.; Rheault, F. On productivity, selectivity and useful exergy effect of separation processes. Proceedings of 24th European symposium on applied thermodynamics; Santiago de Compostela, Spain, June 27-July 1, 2009. (14) Szargut, J.; Morris, D. R.; Steward, F. R. Exergy Analysis of Thermal, Chemical and Metallurgical Processes; Hemisphere Publishing Corporation: New York, 1988. (15) Rivero, R.; Rendon, C.; Monroy, L. The Exergy of Crude Oil Mixtures and Petroleum Fractions: Calculation and Application. Int. J. Appl. Thermodyn. 1999, 2, 115–123. (16) Henley, E. J.; Seader, J. D. Equilibrium-Stage Separation Operations in Chemical Engineering; John Wiley & Sons, Inc.: New York, 1981. (17) Doherty, M. F.; Malone, M. F. Conceptual design of distillation systems; McGraw-Hill: Boston, 2001. (18) King, C. J. Separation Processes, second edition; McGraw-Hill: New-York, 1980. (19) Gottschlich, D. E.; Roberts, D. L. Energy minimization of separation processes using conVentional/membrane hybrid systems; Final Report under DOE Contract No. DE91-004710, 1990. (20) Prausnitz, J. M.; Lichtenthaler, R. N.; de Azevedo, E. G. Molecular thermodynamics of fluid-phase equilibria; PTR Prentice Hall: New Jersey, 1986. (21) Teramoto, M.; Takeuchi, N.; Maki, T.; Matsuyama, H. Ethylene/ Ethane separation by facilitated transport membrane accompanied by permeation of aqueous silver nitrate solution. Sep. Purif. Technol. 2002, 28, 117–124.

ReceiVed for reView June 3, 2009 ReVised manuscript receiVed September 22, 2009 Accepted December 28, 2009 IE9009155