Mathematical and numerical analysis of classes of property models

Fluid Phase Equilibria 250 (2006) 1–32

Mathematical and numerical analysis of classes of property models Rafiqul Gani a,∗ , Nuria Muro-Su˜ne´ a , Mauricio Sales-Cruz a , Claude Leibovici b , John P. O’Connell c a

CAPEC, Department of Chemical Engineering, Technical University of Denmark, DK-2800 Lyngby, Denmark b CFL Consultant, Pau, France c Department of Chemical Engineering, University of Virginia, Charlottesville, VA 22911, USA Received 4 April 2006; received in revised form 24 August 2006; accepted 14 September 2006 Available online 19 September 2006

Abstract A general analysis is presented to focus on the mathematical and numerical elements of property models for stand-alone results as well as within other calculations. This includes listing all of the equations that constitute a property model and articulating the sets of known and unknown variables for given problems. Ordering of equations is done to determine if a direct solution is feasible where the incidence matrix of unknown variables and equations can be written in tridiagonal form or if simultaneous or iterative solution of multiple equations is required. Example pure compound and mixture property models are: cubic equations of state; more advanced non-cubic equations of state; activity coefficient models; and models for polymer and electrolyte systems. Uses of this analysis are discussed for property calculations as a part of other calculations, such as for process models and parameter regression. © 2006 Elsevier B.V. All rights reserved. Keywords: Property models; Mathematical analysis; Numerical analysis; Thermodynamic consistency

1. Introduction Property models are commonly used for synthesis, design, control and analysis of chemical processes and products. In many cases they are used as stand-alone (when property values needed for a specific problem are not available through a database, they need to be calculated through a model) or as part of a bigger model (when the property model is a subset of the total model equations as in a process simulation model). In either case, they are used in the so-called service role [1] where, given the values of a set of intensive variables (such as temperature, pressure and/or composition) and the chemical compound identities, the required properties are calculated. In this way, the property models are usually treated as black-box or procedure equations. That is, values for the known variables are transferred through an interface to a black-box or procedure, which returns calculated property values. Then, the model calculations form an inner calculation loop that is repeated for every request for property values by an outer calculation loop, which occurs when any of the intensive variables change. Although attempts have been made in the past to analyze the property model equations [2–4], the focus was on how to reduce the computational time per call of the property model rather than thorough analysis of the property model equations themselves. Other attempts have been made, for example, by Michelsen [5,6] and Hendriks [7] to reduce the number of equations to be solved per call of inner-loop property models, and Leibovici [8] to order the property model equations, e.g., cubic equations of state, into variant and invariant sets depending on the outer loop dependence on intensive variables. From the points of view of more effective application and education, as well as a better understanding of property models, we believe a thorough mathematical and numerical analysis could be of value. For reliable and efficient usage, in addition to their principles, approximations, limitations and values, property models should be appreciated in terms of the equations that need to be ∗

Corresponding author. Tel.: +45 45252882; fax: +45 45932906. E-mail address: [email protected] (R. Gani).

0378-3812/$ – see front matter © 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.fluid.2006.09.010

2

R. Gani et al. / Fluid Phase Equilibria 250 (2006) 1–32

solved, the model parameters that need to be supplied, and the model (internal) variables that need to be calculated. This information can also help users to implement property models on their own for service, advice and solve roles [1]. Finally, such analysis should be valuable for utilizing third party software as well as implementation of property model exchange, e.g., CAPE-OPEN [9], when potential users know more about model content. The objective of this paper is to present a systematic mathematical and numerical analysis that can be applied to a wide range of property models (also known as constitutive models). The analysis provides for each model a list of the set of independent (model) equations and their corresponding variables together with a classification of the variables as known and unknown variables. The known variables are further classified in terms of those fixed by the system, those fixed by the problem, those set by the property model and finally, those that are regressed or retrieved from model parameter tables. The set of unknown variables and the corresponding set of independent model equations are presented through an ordered incidence matrix of equations and (unknown) variables, which highlights the structure that elucidates issues in solving the equations. In particular, if the structure is lower tridiagonal, the solution approach is simpler than when the structure is not lower tridiagonal, in which case, simultaneous solution (or an extra iteration-loop) is necessary. The property models analyzed in this paper are cubic equations of state [10], advanced non-cubic equations of state: CPA [11–13]; PC-SAFT [14,15]; and SAFT [16], models for polymer systems: GC-Flory [17]; and models for electrolyte systems: UNIQUAC-Elec: [18,19]. Finally, these aspects of property models are discussed in connection with parameter regression and in the common process engineering calculations of two-phase separations and two-phase equilibrium. 2. Property model analysis The systematic analysis performed for each property model consists of the following steps: Step 1. List the independent set of equations representing the model for the desired property and all the variables found in the equations; classify the variables as scalars, vectors and matrices. Step 2. Determine the degrees of freedom (DOF) by subtracting the number of equations from the number of variables. Based on the DOF, select the variables that need to be specified and classify them as those fixed by the problem, fixed by the system, fixed by the property model and adjustable (regressed) model parameters. The remaining variables would typically be the remaining state conditions, properties and intermediate values. Step 3. Establish an incidence matrix of equations and all variables (except those fixed by the system and the property model). If the case has relatively few variables, this may be done visually by putting equations in columns and variables in rows and putting a * on the column–row index where a variable appears in an equation. Step 4. Eliminate the columns for variables that have been selected as specified variables (i.e., variables set by the problem). The incidence matrix must then become square. Step 5. Order the equations so that a lower tridiagonal form will appear, if possible. (Standard equation ordering techniques exist in computer-aided modeling tools [24] to obtain the ordered equation set.) a. If the incidence matrix shows a lower tridiagonal form, then all the equations can be solved sequentially (one unknown for each equation) corresponding to the order giving the lower tridiagonal form. Expect some of the equations to be non-linear, requiring iterative solution. b. If there are elements in the upper tridiagonal portion of the matrix, equations will need to be solved simultaneously and/or iteratively. Note that steps 3–5 may also be combined by directly generating the ordered incidence for the square system. In the text below, some common property models are analyzed according to the systematic analysis given above. We do not give the model derivations since these can be found in the original references. Only the set of equations representing the property model are given and analyzed. In the model equations and their analysis below, a vector will be represented by an underlined variable name (for example, x), a matrix will be represented by an underlined bold variable name (for example, u). An element of a vector or a matrix will be indicated by subscripts (for example, xi or uij ). In this case, unless otherwise indicated, variables will be single-valued. The notation we have used is typically that of the papers cited. This means that the same quantity may have different symbols for different models and vice versa. We have denoted in the text those variables that might not be obvious, but the reader should consult the original references for symbols of uncertain meaning. 2.1. Cubic equations of state: SRK EOS We consider the SRK cubic equation state in its most commonly used form [10]. First we will consider the pure compound property, followed by the compound properties in a mixture.

R. Gani et al. / Fluid Phase Equilibria 250 (2006) 1–32

3

Table 1a List of property model equations for pure compound property calculation with the SRK EOS Equations

Number of equations

RT ai − V − bi V (V + bi ) mi = 0.48 + 1.574ωi − 0.176ωi2 ψB RTci bi = Pci P=



ai = ψA

R2 (Tci ) Pci

 2

1 + mi

1 − T

(1)

1

(2)

1

(3)

1

(4)

1

1/2 2

Tci

Total number of equations = 4 Total number of variables = 12 [R, T, P, V, ψA , ψB , ai , bi , mi , ωi , Tci , Pci ] Table 1b Incidence matrix for SRK EOS for pure components Equations

1 2 3 4

Variables ai

bi

*

*

mi

P

V

T

*

*

*

* * *

*

*

Table 1c Ordered incidence matrix for calculating the molar volume of a pure compound at specified T, P with the SRK EOS Equations

Variables bi

3 2 4 1

mi

ai

V

* *

*

* * * *

2.1.1. Pure compound property The SRK EOS relates the pressure (P), molar volume (V), and temperature (T) of a single pure compound and is written as: P=

RT a − V − b V (V + b)

(1)

In the above equation, R is the universal gas constant while a and b are model parameters. Therefore, given values of R, a, b, T, and P, Eq. (1) can be used to calculate V. However, the value of the parameter a depends on T as well as the model and properties of the compound, while parameter b depends only on the model and the compound properties. Step 1: The independent equations and variables are listed in Table 1a. Step 2: From Table 1, the total number of equations is 4 and the total number of variables is 12, meaning that the DOF is 8. This means that eight variables need to be specified. These are the following: three set by the acentric factor, critical temperature and critical pressure of the system compound (ωi , Tci , Pci ), three set by the property model (R, ψA , ψB ) and two set by the problem (any two from P, V, T). Step 3: An incidence matrix for all the equations and variables minus those set by the system and property model is given in Table 1b. Table 1d Ordered incidence matrix for calculating pure compound volume at specified P, V with SRK EOS Equations

Unknown variables bi

3 2 4 1

mi

ai

T

* *

* *

* * * *

4

R. Gani et al. / Fluid Phase Equilibria 250 (2006) 1–32

Table 2a List of property model equations for a mixture property calculation with the SRK EOS Equations

Number of equations

RT a − V −b V (V + b) mi = 0.48 + 1.574ωi − 0.176ωi2 ,

P=

bi =

ψB RTci , Pci



αi = ψA

i = 1, NC

i = 1, NC

R2 (Tci ) Pci

 2

1 + mi

1 − T Tci

1/2 2 ,

i = 1, NC

b = Σi xi bi aij = (αi αj )1/2 (1 − kij ),

i = 1, NC; j = i, NC

a = Σi Σj xi xj aij

(1)

1

(2)

NC

(3)

NC

(4)

NC

(5)

1

(6)

NC(NC + 1)/2

(7)

1

Total number of equations = 3 + 3NC + NC*(NC + 1)/2 Total number of variables = 8 [R, T, P, V, ψA , ψB , a, b] + 7NC [x- , a- , b- , m - , ω- , T- c , P- c ] + 2NC ∗ (NC + 1)/2 [a- , k- ] = 8 + 7NC + 2NC ∗ (NC + 1)/2

Step 4. To make the matrix square, two columns must be removed by specifying two variables from those set by the problem (any two from P, V, T). There are three possible problem cases: (a) specify T, P and find V; (b) specify T, V and find P, and (c) specify P, V and find T. In all cases, two of the last three columns of the matrix are eliminated. Step 5. The incidence matrix for case (a) is shown in Table 1c after ordering is done to achieve a lower tridiagonal form. Problem (b) gives the same result when P is substituted for V in the matrix. However, problem (c) does not give a lower tridiagonal form, as shown in Table 1d. This means that, for problem (c), Eqs. (1) and (4) need to be solved simultaneously for T. 2.1.2. Bulk and compound properties in a mixture Bulk properties are for the whole of a single-phase mixture, such as the molar volume, V, or residual enthalpy at specified T, P, and x- (compound compositions). Compound properties of a single-phase mixture are associated with variables such as the partial molar volume or fugacity coefficient of each compound present in the mixture at the specified T, P, and x- . 2.1.2.1. Analysis for bulk property of the mixture. This case is similar to that in Section 2.1.1. Here we will focus on obtaining the molar volume, V, at a given T and P of a mixture of NC compounds with composition vector x- . Step 1: The list of property model equations for the cubic SRK EOS is given in Table 2a. Step 2: From Table 2a, the DOF is 5 + 4NC + NC*(NC + 1)/2. The 5 + 4NC + NC*(NC + 1)/2 variables that need to be specified are classified as, NC + 2 variables fixed by the problem [T, P, x- ], 3NC fixed by the system compounds [ω- , T- c , P- c ], three fixed by the property model [R, ψA , ψB ] and NC*(NC + 1)/2 variables [k- ] are regressed (or retrieved from model parameter tables) model parameters (with all kii = 0). For the binary system with NC = 2, this makes 12 model equations with 16 variables that must be specified. In this case, they are all constants. The remaining 12 variables [V, a, b, α1 , α2 , a11 , a12 , a22 , b1 , b2 , m1 , m2 ] are the unknown variables for this problem. Steps 3–5: The incidence matrix for V can always be arranged in tridiagonal form. For illustrative purposes only the binary case is shown in Table 2b. As for the pure component case (Table 1c), the incidence matrix for calculating a mixture T, given V, P, and x- , will again not show a tridiagonal form. 2.1.2.2. Analysis for compound properties in a mixture. Consider now a NC compounds mixture at given T, P and x- , whose phase (liquid or vapor) is also known. We are interested in calculating the component fugacity coefficients for each compound i, ϕi , of the mixture at the specified condition. Step 1: The model equations needed for using the cubic SRK EOS are those listed in Table 2c. Step 2: From Table 2c, the DOF is: 4NC + NC*NC + 5. The 4NC + NC*NC + 5 variables that need to be specified are classified as, NC + 2 variables fixed by the problem [T, P, x- ], 3NC fixed by the system compounds [ω- , T- c , P- c ], three fixed by the property model [R, ψA , ψB ] and NC*NC variables [k- ] are regressed (or retrieved from model parameter tables with all kii = 0 and kij = kji ). For the binary system with NC = 2, this makes 29 model equations and 46 variables, out of which 17 must be specified. The remaining 29 variables [Z, A, B, a, b, T- R , m - , α- , a- , b- , a- , b- , a- r , b- r , ϕ- ] are the unknown variables for this problem. Note that because of symmetry, only NC*(NC + 1)/2 elements of the matrices a- and b- are unique. Step 3–5: The ordering of the 8NC + 2NC*NC + 5 equations is highlighted in Table 2d through the incidence matrix. Use of this and other models involving the computation of multiphase equilibrium calculations is highlighted in Section 3. Example 1 in Appendix B highlights the implementation of the above calculation order for a binary mixture of n-decane/ethane.

R. Gani et al. / Fluid Phase Equilibria 250 (2006) 1–32

5

Table 2b Incidence matrix of SRK EOS property model equations for a binary mixture (bulk property calculation) m1

m2

b1

b2

α1

α2

a11

a12

a22

b

a

V

* *

*

*

2

* *

3

* 4

*

* *

* * *

6

* *

*

* 5

*

*

*

*

7 1

*

*

* *

Table 2c List of SRK EOS model equations for calculation of fugacity coefficients Equations ln φi = bri (Z − 1) − ln(Z − B) −



ari =

0.5 (α0.5 i Tci )/Pci

j

Tci /Pci



bri =

0.5 (xj α0.5 j Tcj /Pcj )

J

(xj Tcj /Pcj )

,

,



A 0.5 (2ari − bri ) ln 1 + B Z

i = 1, NC

i = 1, NC

Z3 − Z2 + Z(A − B − B2 ) − AB = 0 aP

A= B= a= b=

(RT )2 bP RT

i

xi

j



i

xj aij ,

i = 1, NC; j = 1, NC

i = 1, NC; j = 1, NC

xi xj bij ,

aij = (ai aj )1/2 (1 − kij ), bij =

bi + bj , 2



ai = ψA αi (T ) bi = ψB

ci

Pci

i = 1, NC; j = 1, NC

i = 1, NC; j = 1, NC

,

i = 1, NC

,

i = 1, NC 2

0.5 αi (T ) = [1 + mi (1 − TRi )] ,

i = 1, NC

mi = 0.48 + 1.574ωi − 0.176ωi2 , TRi =

T , Tci

Number of equations ,

i = 1, NC

(1)

NC

(2)

NC

(3)

NC

(4)

1

(5)

1

(6)

1

(7)

1

(8)

1

(9)

NC*NC

(10)

NC*NC

(11)

NC

(12)

NC

(13)

NC

(14)

NC

(15)

NC

 2

R2 Tci Pci

 RT 

 B

i = 1, NC

i = 1, NC

Total number of equations = 5 + 8NC + 2NC*NC Total number of variables = 10 [T, P, R, ψA , ψB , a, b, A, B, Z] + 12NC [x- , T- c , P- c , ω- , T- R , α- , m - , a- , b- , a- r , b- r , ϕ- ] + 3NC ∗ NC[a- , b- , k- ] = 10 + 12NC + 3NC ∗ NC

6

R. Gani et al. / Fluid Phase Equilibria 250 (2006) 1–32

Table 2d Incidence matrix of SRK EOS property model equations Equations

Unknown variables Tr

15 14

*

13 12

*

m

α

b-

b-

a-

a

* *

*

b

B

* *

*

A

Z

* *

*

b- r

a- r

ϕ -

* *

* *

11 10

*

* *

9 7

* *

8 6

*

5 4a

* *

3 2

* *

1 a

a-

* *

*

*

*

*

*

Requires the solution of a cubic equation in Z.

2.2. Advanced non-cubic equations of state The analysis given above is now repeated for the CPA EOS [11–13], the PC-SAFT EOS [14,15] and the SAFT EOS [16]. For the CPA EOS, the incidence matrix for the ordered property model equations is illustrated for the calculation of a pure component property (fugacity coefficient) and for the component fugacity coefficients in a mixture, while for the PC-SAFT and SAFT EOS, only the case of component fugacity coefficients in a mixture is described. 2.2.1. CPA EOS model analysis 2.2.1.1. Analysis for pure component property. The CPA EOS model is analyzed for the calculation of the fugacity coefficient of a pure associating compound at given T and P. The compound has NS unique sites for association with symmetric parameters εAB = εBA giving the variable AB = BA . Step 1: The CPA EOS model equations [11] are listed in Table 3a. Note that since the calculation of ZSRK requires the same set of equations as listed in Table 2c and since their incidence matrix shows a lower tridiagonal form, they can also be represented by a single (procedure) equation (Eq. (3) in Table 3a). Note that Eqs. (6) and (7) of Table 3a can be reduced to the following forms for a pure associating compound with two active sites (NS = 2) using the 2B association model (see Appendix A). XA =

−1 + (1 + 4ρ )1/2 2ρ 2

−(XA ) (XA + XA ρ(∂ /∂ρ)) ∂XA = 2 ∂ρ 1 + ρ (XA ) In the above equations, XA = XB . Note that we have used the simplified version of CPA in Eq. (9) [11]. Step 2: From Table 3a, the DOF is 6 + NS/2. The 6 + NS/2 variables that need to be specified are classified as two variables fixed by the problem [T, P], one variable fixed by the property model [R], 1 + NS/2 variables that are regressed (or retrieved from model parameter tables) parameters [ε- AB , β] and two variables [a, b] that are internal parameters of the SRK EOS (see Section 2.1.2). The unknown variables are therefore [Z, V, ZSRK , Zassoc , ρ, g, ∂g/∂ρ, y, ϕ, Ar /RT, ASRK /RT, Aassoc / RT, AB , ∂ AB /∂ρ, XA- , ∂XA- /∂ρ, XB- , ∂XB- /∂ρ], with the two internal parameters for the SRK EOS calculated or supplied through the SRK EOS model equations. Step 3–5: The ordered incidence matrix for the above CPA model consisting of 6 + 3NS equations is highlighted in Table 3b for NS = 2. It can be noted that a lower tridiagonal structure is not obtained for this model. Note, however, that we are not claiming that Table 3b shows the optimal order of the CPA EOS model equations.

R. Gani et al. / Fluid Phase Equilibria 250 (2006) 1–32

7

Table 3a The CPA EOS model equations Equations

Number of equations

Z = ZSRK + Zassoc ZRT V = P V a SRK Z = − V −b RT (V + b)

(1)

1

(2)

1

(3)

1

(4)

1

(5)

1

(6)

NS

∂XA −(XA ) (XA AB + XA ρ(∂ /∂ρ)) = 2 ∂ρ 1 + ρ (XA )

(7)

NS

∂ AB ∂g = exp ∂ρ ∂ρ

(8)

NS/2

(9)

1

(10)

1

(11)

NS/2

(12)

1

(13)

1

(14)

1

(15)

1

(16)

1

 1 NS

Zassoc = ρ

−

XA

 1 ∂XA

2

∂ρ

A

ρ=

1 V



XA =



−1

NS

1+ρ

XB AB

,

A = 1, NS

B 2





εAB RT

1 1 − 1.9y ∂g 1.9b/4 = ∂ρ (1 − 1.9y)2





− 1 bβ

g=





AB = g exp y=



 − 1 bβ

b 4V

ln ϕ = Ar

εAB RT

Ar RT

+ Z − 1 − ln Z

ASRK

Aassoc

= + RT RT RT     ASRK b a b − = −ln 1 − ln 1 + RT V RTb V



1 1 Aassoc = ln XA − XA + RT 2 2 NS



A

Total number of equations = 12 + 3NS Total number of variables = 18 [Z, V, ZSRK , T, R, Zassoc , a, b, P, ρ, β, g, ∂g/∂ρ, y, ϕ, Ar /RT, ASRK /RT, Aassoc /RT ] + AB AB AB A A B B 3NS/2 [ - /∂ρ, ε- ] + 2NS [X - , ∂X - /∂ρ, X - , ∂X - /∂ρ] - , ∂

2.2.1.2. Analysis for compound properties in a mixture. Consider now a mixture of NC compounds mixture at a given T, P and x- , whose phase (liquid or vapor) is also known. We wish to calculate the component fugacity coefficients for each compound i of the mixture at the specified conditions. Of the NC compounds, NCA compounds are associating and NCNA compounds are non-associating, while NS is the number of active sites for an associating compound (here assumed to be the same for all compounds). Step 1: The model equations needed for using the CPA EOS are those listed in Table 3c. Eq. (19) (from Table 3c) for an associating compound with two active sites A and B (NS = 2), that is, for the 2B association model (see Appendix A) when XBi = XAi becomes (Eq. (19a) with reference to Table 3c) XAi =

−1 + (1 + 4(xi /V ) Ai Bj )1/2 2(xi /V ) Ai Bj

Note that for one associating compound with two active sites, i = j = 1 here and that for Eq. (16) (Table 3c), the simplified CPA form is used [11]. Step 2: From Table 3c, the DOF is 6 + 5NC + 5NCA + NC*NC + NS*NCA, where the variables to be specified are classified as: 3 + 2NC variables fixed by the problem [T, P, n, x- , n- ], 3NC variables fixed by the system [P- c , T- c , ω- ], three variables fixed by the property model [R, Ωa , Ωb ] and 5NCA [b- , a- 0 , c- 1 , εA , βA ] plus NC*NC [k- ] (where all kii = 0 and kij = kji ), plus NS*NCA [n- A ] variables that must be retrieved from the model parameter tables or regressed.

8

R. Gani et al. / Fluid Phase Equilibria 250 (2006) 1–32

Table 3b Incidence matrix of CPA property model equations for a pure property calculation (note that XA = XB ) Equations

Unknown variables y

g

12 9

* *

*

10 8

*

11 5

dg/dρ

d AB /dρ

AB

P

XA

dXA /dρ

ZSRK

Z

* *

ASRK /RT

Ar /RT

ϕ

* *

*

*

*

* *

4 3

*

* *

* *

*

*

*

*

* *

1 2

* *

*

*

* * *

* *

*

* *

*

Aassoc /RT

* *

15 16

V *

*

6 7

14 13

Zassoc

*

Steps 3–5: The CPA model equations used for computation of component fugacity coefficients in a mixture, represented by 11 + 9NC + NCA + 3NCNA + 2NC*NC + NS*NCA + 3NCA*NCA equations, are ordered as shown in Table 3d. Again, it can be seen that a lower-tridiagonal form is not obtained, meaning that a sub-set of equations will need to be solved simultaneously (or an iterative scheme needs to be derived for their solution). Example 2 in Appendix B highlights the implementation of the above calculation order for a binary mixture of methanol/methane. 2.2.2. PC-SAFT EOS model analysis The PC-SAFT EOS model [14,15] is analyzed for the calculation of the fugacity coefficient of compound k in a NC component mixture at given T, P and x- . Step 1: The PC-SAFT EOS model equations, as described by Gross and Sadowski [14,15], are listed in Table 4a. The right hand side for Eq. (33) in the above table (Table 4a) is given by:

   ∂˜ahs ζ0,xk 1 3(ζ1,xk ζ2 + ζ1 ζ2,xk ) 3ζ1 ζ2 ζ3,xk 3ζ22 ζ2,xk ζ 3 ζ3,x (3ζ3 − 1) hs a˜ + =− + + + 2 2 k 2 2 ∂xk ζ0 ζ0 1 − ζ3 (1 − ζ3 ) ζ3 (1 − ζ3 ) ζ2 (1 − ζ3 )3 T,v,xj =k

(ζ0 − ζ23 /ζ32 )ζ3,xk 3ζ22 ζ2,xk ζ3 − 2ζ23 ζ3,xk + − ζ0,xk ln(1 − ζ3 ) + (1 − ζ3 ) ζ33 Step 2: The DOF is found to be: 4NC + 3(Ja + Jb ) + NC*NC + 2, where the variables to be specified are classified as, NC + 2 variables fixed by the problem [T, P, x- ], 3NC variables fixed by the system [m - , σ- , (ε- /k)], 3(Ja + Jb ) variables fixed by the property model [a- 0 , a- 1 , a- 2 , b- 0 , b- 1 , b- 2 ], and NC*NC [k- ] parameters regressed or retrieved from model parameter tables (with kii = 0 and kij = kji ). Steps 3–5: The set of NC*(3NC + 14) + Nz*(NC + 1) + Ja (NC + 1) + Jb (NC + 1) + 20 unknown variables are listed in Table 4b, which also shows the ordered incidence (square) matrix of the model equations for the calculation of the compound fugacity coefficients in a mixture of NC compounds. It can be noted that a lower tridiagonal form is not obtained. Example 3 in Appendix B highlights the implementation of the above calculation order for a binary mixture of n-decane/ethane. 2.2.3. SAFT EOS The SAFT EOS model [16] is analyzed for the calculation of the fugacity coefficient of compound k in a NC component mixture at given T, P and x- . Of the NC compounds, NCA compounds are associating while NS is the number of active sites for an associating compound. Step 1: The SAFT EOS model equations are listed in Table 4c. Note that the equation numbers used in column 1 of Table 4c are the same as those used by Huang and Radosz [16].

R. Gani et al. / Fluid Phase Equilibria 250 (2006) 1–32

9

Table 3c List of property model equations for mixture property calculation with CPA Equations

Number of equations

Z = ZSRK + Zassoc

(1)

1

ZRT V = P

(2)

1

(3)

1

(4)

NCA

(5)

NCNA

(6)

NCNA

(7)

NCNA

(8)

NC

(9)

NC*NC

ZSRK =

a V − V −b RT (V + b)

ai = a0,i (1 + c1,i (1 − 2 Ωa αi R2 Tc,i

ai =

Pc,i

αi = [1 + (0.480



2

Tr,i )) ,

i = 1, NCA

i = NCA + 1, NC

,

1/2 2 + 1.574ωi − 0.176ωi2 )(1 − Tr,i )] ,

i = NCA + 1, NC

Ωb RTc,i bi = , i = NCA + 1, NCNA Pc,i T , i = 1, NC Tr,i = T √ c,i aij = ai aj (1 − kij ), i = 1, NC; j = 1, NC a=

NC NC

(10)

1

(11)

1

(12)

NC*NC

(13)

1

(14)

1

(15)

1

(16)

1

(17)

1

(18)

1

i = 1, NCA

(19)

NS*NCA

i = 1, NCA; j = 1, NCA

(20)

NCA*NCA

(21)

NCA*NCA

(22)

NCA*NCA

(23)

NC

(24)

NC

xi xj aij i=1 j=1

b=

NC NC

xi xj bij i=1 j=1

bij =

1 (bi + bj ), 2

Zassoc = −

1 2



i = 1, NC; j = 1, NC

1+ρ

∂ ln g ∂ρ

NC NC  i=1

1 V 1 ∂g ∂(ln g) = ∂ρ g ∂r 1 g= 1 − 1.9y b y= 4V ∂g 1.9b/4 = ∂ρ (1 − 1.9y)2

xi nA (1 − XAi )

Ai

ρ=

⎡

XAi = ⎣1 + ρ

NC NS



βij =

xj XBj Ai Bj ⎦

j=1 Bj

Ai Bj = g exp εij =

⎤−1

 εAi Bj 

1 i (ε + εj ), 2



RT



,

− 1 bij βAi Bj ,

i = 1, NCA; j = 1, NCA

βi βj , i = 1, NCA; j = 1, NCA μi − ln Z, i = 1, NC ln ϕi = RT   μi Ar ∂ , i = 1, NC = RT ∂ni RT      Aassoc  Ar ASRK ∂ ∂ ∂ = + , ∂ni RT ∂ni RT ∂ni RT μSRK ∂ i = ∂ni RT



ASRK RT



= −ln(1 − b/V ) +

2 nj aij ln(1 + b/V ), RTb

i = 1, NC

a n + V −b RTb

(25)

 ln(1 + b/V ) b

−

NC 2 1 V +b

n

NC

 nj bij − b

NC

j=1

NC

−

j=1

i = 1, NC

(26)

10

R. Gani et al. / Fluid Phase Equilibria 250 (2006) 1–32

Table 3c (Continued) Equations μassoc ∂ i = RT ∂ni

Number of equations

 Aassoc 

1 ∂g ∂ ln g = , ∂ni g ∂ni

RT

=

NS

n ∂ ln g xj nBj (1 − XBj ), 2 ∂ni NC

nA ln XAi −

Ai

i = 1, NC

∂g(ρ) 1.9(∂y/∂ni ) = , i = 1, NC ∂ni (1 − 1.9y)2   ∂y 1 bi n − b = + b , i = 1, NC ∂ni 4V n

NS

i = 1, NC

(27)

NC

(28)

NC

(29)

NC

(30)

NC

j=1 Bj

Total number of equations = 11 + 9NC + NCA + 3NCNA + 2NC*NC + NS*NCA + 3NCA*NCA Total number of variables = 17 [Z, V, ZSRK , T, R, P, Zassoc , Ωa , Ωb , a, b, ρ, ∂ ln g/∂ρ, g, ∂g/∂ρ, y, n] + 6NCA [a- , b- , a- 0 , c- 1 , εA- , βA- ] + 3NCNA [a- , b- , α- ] + 14NC [P- c , T- c , ω- , T- r , x- , ϕ, μ/RT, ∂Ares /∂n- , ∂ASRK /∂n- , ∂Aassoc /∂n- , n- , ∂ ln g/∂n- , ∂g/∂n- , ∂y/∂n- ] + 3NC ∗ NC[k- , b- , a- ] + 2NS ∗ - NCA[n- A , XA ] + 3NCA ∗ NCA[ΔAB , εAB , βAB ]

Step 2: The DOF is found to be: 42 + 4NC + 3NCA + NC*NC, where the variables to be specified are classified as: 2 + NC [T, P, X -] variables fixed by the problem, NCA [M ] variables fixed by the system, 41 variables fixed by the property model (D , τ, N , CI, R) Av 00 A A and 4NC [u- 0 , e- /k, m - , v- ] + 2NCA [ε , κ ] + NC ∗ NC[k- ] variables that must be retrieved or regressed (with all kii = 0 and kij = kji ). Steps 3–5: For the calculation of the specified problem in Step 1, the ordered incidence matrix for the model equations from Table 4c are highlighted in Table 4d . It is clear that a lower tridiagonal form is not obtained. Example 4 in Appendix B highlights the implementation of the above calculation order for a binary mixture of methane/hexadecane. 2.3. Activity coefficient models Unlike the EOS property models, the activity coefficient models, as their name suggests, give directly the activity coefficient γ i for compound i in NC compounds mixture. Two-activity coefficient models are highlighted below, the Wilson model [20] and the UNIFAC group contribution model [21]. The structure of the most other well-known activity coefficient models would be similar to these. 2.3.1. Wilson activity coefficient model The Wilson [20] model is analyzed below for the computation of activity coefficients of the compounds 1 and 2 in a binary mixture. That is, for a given T, x- , values of γ i (i = 1, 2) are calculated. Step 1: The Wilson model equations are listed in Table 5a. Step 2: From the model equations listed in Table 5a, the DOF is 9, where the variables to be specified are classified as three variables fixed by the problem [T, x- ] two variables fixed by the system compounds [V- ] and four variables [A - ] that must be retrieved from model parameter tables or regressed (with A11 = A22 = 0). Steps 3–5: The ordered incidence matrix is shown in Table 5b for the eight unknown variables. It can be noted that a tridiagonal structure is obtained. Consider a binary equimolar mixture of methanol–water at 300 K. According to the above analysis, in order to calculate ln γ 1 (methanol) and ln γ 2 (water), we need to specify three variables fixed by the problem [T = 300, x1 = x2 = 0.5]; two variables fixed by the system compounds [V1 = 0.0805, V2 = 0.0181] and four variables [A12 = −22.6004, A21 = 298.3496, with A11 = A22 = 0] that are retrieved from model parameter tables. For these specified variable values, solution of the Wilson model equations according to Tables 5a and 5b will give ln γ 1 = 0.0491 and ln γ 2 = 0.1471. 2.3.2. UNIFAC-VLE (original) activity coefficient model The UNIFAC-VLE model [21] is analyzed for the computation of the liquid phase activity coefficients of compound i in a NC compound mixture at specified T and x- . The system compounds have a total of NMG main groups and NSG subgroups. Step 1: The UNIFAC-VLE model equations, as described by Fredenslund et al. [21], are listed in Table 5c. Step 2: Based on the equations listed in Table 5c, the DOF for the UNIFAC-VLE model is 1 + 3NC + 2NMG*NMG, where the variables to be specified are classified as: 1 + NC [T, x- ] variables fixed by the problem, NSG*NC [v- ] variables fixed by the system and 2NSG [Q, R ] plus NMG*NMG [a- ] parameters that must be retrieved from the model parameter tables or regressed (with all aii = 0). - Steps 3–5: The ordered incidence matrix for the set of unknown variables is shown in Table 5d. It can be noted that a lower tridiagonal structure is obtained.

Table 3d Incidence matrix of CPA model equations (note that XAi = XBi ) Equations

Unknown variables B

b-

7 12

* *

*

11 17

*

*

Y

* *

*

*

* *

15 13

g

dg/dρ

d ln g/dρ

ρ

ε- AB

βAB -

T- r

a-

α-

a- 0

a-

a

AB Δ -

XAi

V

dy/dn-

dg/dn-

* *

*

d ln g/dn-

d(Aassoc /RT )/dn-

* *

*

d(ASRK /RT )/dn-

d(Ares /RT )/dn-

μ/RT -

Φ -

* *

*

* *

*

* *

* * * * *

6 5

*

9 10

* * * *

20 19

*

*

*

* *

* *

*

*

* *

*

3 13

*

*

* *

*

* *

1 2 *

* *

* *

*

* *

* *

*

28 27

24 23

Z

*

8 4

26 25

Zassoc

*

21 22

30 29

ZSRK

R. Gani et al. / Fluid Phase Equilibria 250 (2006) 1–32

16 18

b

*

* *

*

*

*

*

* *

* *

* *

*

11

12

R. Gani et al. / Fluid Phase Equilibria 250 (2006) 1–32

Table 4a The PC-SAFT EOS model equations and variables Equations

Number of equations

μk − ln(Z), kT hc Z = 1 + Z + Zdisp

ln φk =

a˜

res

= a˜ + a˜ hc

a˜ hc = m˜ ¯ ahs − m ¯ =

disp





i

k = 1, NC

i

xi mi ,

(mi − 1) ln(giihs ),

i = 1, NC

i = 1, NC 1

a˜ hs =

ζ0 [3ζ1 ζ2 /(1 − ζ3 ) + ζ23 /ζ3 (1 − ζ3 )2 + (ζ23 /ζ33 − ζ0 ) ln(1 − ζ3 )]

1 + (di /2)(3ζ3 /(1 − ζ3 )2 + (di /2)2 (2ζ22 /(1 − ζ3 )3 ), (1 − ζ3 ) π ζn = ρ xi mi din , i = 1, NC; n = 0, . . . , 3 6 i   ε  r di = σi 1 − 0.12 exp , i = 1, NC kT giihs =

a˜ disp = −2πρI1 m2 εσ 3 − πρmC ¯ 1 I2 m2 ε2 σ 3



C1 =

1 + m(8η ¯ − 2η2 ) (1 − η)4

+

(1 − m)(20η ¯ − 27η2 + 12η3 − 2η4 )

i = 1, NC

−1

[(1 − η)(2 − η)]2 ε  ij m203 = m2 εσ 3 = xi xj mi mj σij3 , i = 1, NC; j = 1, NC kT i j  ε 2 ij m223 = m2 ε2 σ 3 = xi xj mi mj σij3 , i = 1, NC; j = 1, NC kT i j 1 σij = (σi + σj ), i = 1, NC; j = 1, NC 2 εij = kT I1 = I2 =

 ε   ε 1/2 j i

j

i = 1, NC; j = 1, NC

1

(5)

1

(6)

1

(7)

NC

(8)

Nz = 4

(9)

NC

(10)

1

(11)

1

(12)

1

(13)

1

(14)

NC*NC

(15)

NC*NC 1

bj η,

j = 0, . . . , 6

(17)

1

j = 0, . . . , 6

(18)

Ja = 7

j = 0, . . . , 6

(19)

Jb = 7

(20)

1

(21)

1

(22)

1

(23)

1

i

(mi − 1)(giihs )

−1

ρ(∂ ln giihs /∂ρ)

3ζ 3 − ζ3 ζ23 ζ3 3ζ1 ζ2 + + 2 2 1 − ζ3 ζ0 (1 − ζ3 ) ζ0 (1 − ζ3 )3

 hs

Zdisp = −

1

(16)

Zhc = mZ ¯ hs −

∂ ln gii ∂ρ

(3) (4)

j = 0, . . . , 6

P = ZkTρ(1010 )

ρ

1

aj η,

¯ − 1) ¯ − 1)(m ¯ − 2) a1j (m a2j (m , + m ¯ m ¯2 ¯ − 1) ¯ − 1)(m ¯ − 2) b1j (m b2j (m bj = b0j + , + 2 m ¯ m ¯  6   −1 η xi mi di3 , i = 1, NC ρ= π i



(2)

kT

aj = a0j +

Zhs =

NC

kT



j

(1 − kij ),

(1)

=

ζ3 di + 1 − ζ3 2





3ζ2 (1 − ζ3 )

2

+

6ζ2 ζ3 (1 − ζ3 )



  2  d i

4ζ22

2

(1 − ζ3 )

3

2πρ∂(ηI1 ) C1 ∂(ηI2 ) − πρm ¯ + C2 ηI2 m2 ε2 σ 3 ∂η ∂ηm2 εσ 3

3

+

6ζ22 ζ3 (1 − ζ3 )4

 , i = 1, NC(24)

NC

(25)

1

∂(ηI1 ) = aj (j + 1)ηj , ∂η j

j = 0, . . . , 6

(26)

1

∂(ηI2 ) = bj (j + 1)ηj , ∂η j

j = 0, . . . , 6

(27)

1

R. Gani et al. / Fluid Phase Equilibria 250 (2006) 1–32

13

Table 4a (Continued ) Equations

Number of equations

 C2 = −C12

2 + 20η + 8) m(−4η ¯

(1 − η)5

μk = a˜ res + (Z + 1) + kT ζn,xk =

π mk (dk3 ), 6

 ∂˜ares    

∂xk

T,v,xj =k

∂˜ahs ∂xk ∂giihs ∂xk



=



∂˜ahc ∂xk

[(1 − η)(2 − η)]3

T,v,xj =k



∂˜ahc ∂xk



 + T,v,xj =k



T,v,xj =k

∂˜ahs ∂xk



= f (˜ahs , ζn , ζn,xk ),



T,v,xj =k

= T,v,xj =k



−

  a˜ res  xj

j

xk



ζ3,xk (1 − ζ3 )2

+

∂˜adisp xk

T,v,xj =k

,

k = 1, NC

2

3

k = 1, NC



223 2 2 3 (m - )x = (m ε σ )xk = 2mk

,

k = 1, NC

T,v,xj =k

−



i

T,v,xj =k

−1 xi (mi − 1)(giihs )



∂giihs ∂xk

NC

(30)

Nz*NC

(31)

NC

(32)

NC

(33)

NC

(34)

NC*NC

 hs − (mk − 1) ln(gkk ),

k = 1, NC

T,v,xj =k

k = 1, NC

(1 − ζ3 )2

+

6ζ2 ζ3,xk (1 − ζ3 )2

  2  4ζ2 ζ2,xk di +

2

(1 − ζ3 )3

+

6ζ22 ζ3,xk



(1 − ζ3 )4

j

xj mj

j

ε  kj

, i = 1, NC; k = 1, NC

kT



j



j

NC (35)

3 σkj ,

3 xj mj (εkj /kT )2 σkj ,

j = 1, NC; k = 1, NC j = 1, NC; k = 1, NC

C1,xk = C2 ζ3,xk − C12 {mk (8η − 2η2 )/(1 − η)4 − mk (20η − 27η2 + 12η3 − 2η4 )/[(1 − η)(2 − η)]2 }

I2,xk =

(29)

= −2πρ[I1,xk m2 εσ 3 + I1 (m2 εσ 3 )xk ] − πρ{[mk C1 I2 + mC ¯ 1,xk I2 + mC ¯ 1 I2,xk ]m2 ε2 σ 3 T,v,xj =k 2 2 3

)x = (m εσ )xk = 2mk

I1,xk =

1



 d   3ζ 2,xk i 2

(28)



n = 0, . . . , 3

+ mC ¯ 1 I2 (m ε σ )xk , (m -

∂xk

= mk a˜ − m ¯

∂˜adisp ∂xk

203

3 + 12η2 − 48η + 40) (1 − m)(2η ¯

 ∂˜ares 

hc



+

−1

, k = 1, NC

(36)

NC

(37)

NC

(38)

NC

[aj jζ3,xk ηj−1 + aj,xk ηj ],

j = 0, . . . , 6; k = 1, NC

(39)

NC

[bj jζ3,xk ηj−1 + bj,xk ηj ],

j = 0, . . . , 6; k = 1, NC

(40)

NC

aj,xk =

a1j mk a2j mk , + 2 m ¯2 m ¯ (3 − 4/m) ¯

j = 0, . . . , 6; k = 1, NC

(41)

Ja *NC

bj,xk =

b1j mk b2j mk , + 2 m ¯2 m ¯ (3 − 4/m) ¯

j = 0, . . . , 6; k = 1, NC

(42)

Jb *NC

Total number of equations = NC*(3NC + 14) + Nz*(NC + 1) + Ja *(NC + 1) + Jb *(NC + 1) + 19 Total number of variables = 21 [T, η, Z, a˜ res , a˜ hc , m, ¯ a˜ hs , a˜ disp , C1 , m203 , m223 , I1 , I2 , ρ, P, Zhc , Zhs , Zdisp , ∂(ηI1 )/∂η, ∂(ηI2 )/∂h, C2 ] + hs hs 203 223 /kT ), (∂˜ares /∂x- )T,v,xj =k , (∂˜ahc /∂x- )T,v,xj =k , (∂˜ahs /∂x- )T,v,xj =k , (∂˜adisp /∂x- )T,v,xj =k , (m 18NC [x- , m - )x , (m - )x , - , σ- , ε- /k, φ- , g- , d- , ρ(∂ ln g /∂ρ), (μ hs C- 1,x , I- 1,x , I- 2,x ] + 3NC ∗ NC[σ- , (ε- /kT ), (∂g /∂x- )T,v,x =k ] + Nz ∗ (NC + 1)[ζ , ζ x ] + Ja (NC + 4)[a- x , a- 0 , a- 1 , a- 2 , a- ] + Jb (NC + 4)[bx , b- 0 , b- 1 , b- 2 , b- ] + j - NC ∗ NC[k- ] = NC(3NC + 18) + Nz ∗ (NC + 1) + Ja (NC + 4) + Jb (NC + 4) + NC ∗ NC + 22 Note that Nz defines the dimensions of variables [ζ , ζ x ], while Ja and Jb define the dimensions of the variables [a- x , a- 0 , a- 1 , a- 2 , a- , b- x , b- 0 , b- 1 , b- 2 , b- ]. The vector g - contains only the diagonal elements of the matrix g.

Example 5 in Appendix B highlights the implementation of the above calculation order for a binary mixture of diethylamine/nheptane. 2.4. GC-Flory EOS model for polymer systems The GC-Flory model [17] is a group contribution based model for the calculation of phase equilibrium involving polymers in solution. The model is analyzed in terms of computation of activity coefficients of compound (polymer) i, in the polymer system. Note that the GC-Flory model has been derived as an EOS. Step 1: The GC-Flory [17] model equations are listed in Table 6a. Step 2: The DOF is found to be: 4NC + NMG*NMG + 6 + 5NSG, where the variables to be specified are classified as: 3NC + 2 variables fixed by the problem [T, x- , n, n- , w - ], 4 + NC variables fixed by the property model [T0 , z, R, P, P- ] and 5NSG , Q ] plus NMG*NMG [ε ] variables that must be retrieved from the model parameter tables or regressed. [C- T0 , C- T , C- 0 , R - -

14

Table 4b Incidence matrix of PC-SAFT EOS model equations Equations Unknown variables 203 ) (m223 ) ζ I disp hs hs hc res d- σ- (ε- /kT ) m ¯ m203 m223 a- b- I2 C1 ρ ζ ghs C2 ∂(ηI1 )/∂η ∂(ηI2 )/∂η Zdisp ρ(∂ ln ghs /∂ρ) Zhs Zhc Z P η I1 a˜ disp a˜ hs a˜ hc a˜ res a- x b- x (m /kT ln φ x x - x - 1,x I- 2,x C - 1,x ∂˜a /∂x- ∂g- /∂x- ∂˜a /∂x- ∂˜a /∂x- ∂˜a /∂x- μ - -

9 14

* *

15 5

*

12 13

* * * *

* * * * *

17 11

*

* * * *

20 8

* *

7 28

*

* * * * * * *

*

*

*

*

*

*

*

*

* *

* *

*

* *

*

* *

* *

* * *

* * *

6 4

*

*

*

* * *

*

*

* *

3 41

* *

* *

* *

*

*

42 36

*

*

*

* *

*

* *

*

*

39 40

* * *

38 35

29 1

* *

* *

*

16 10

32 31

* *

2 21

34 33

*

*

23 22

37 30

*

* * *

26 27 25 24

* *

*

R. Gani et al. / Fluid Phase Equilibria 250 (2006) 1–32

18 19

* * *

*

* *

* *

* *

*

* *

*

*

*

*

* *

*

* *

*

* *

* *

*

* *

* * *

*

* *

*

* *

* *

*

*

* * *

* * *

*

R. Gani et al. / Fluid Phase Equilibria 250 (2006) 1–32

15

Table 4c The SAFT EOS model equations and variables Equations Z =1+ρ

Number of equations

 ∂˜a  ∂ρ

 ∂˜a  ∂ρ

 ∂˜a 

=

∂ρ

T,X

i=1

+

∂βi

i=1

i

∂ρ



∂ ares ∂ρ RT

= T,β



E

T,β

jDij Gi

i

T,ρ,βj=`ı



πD 6

=

 ∂β  ∂ρ

1

(2)

1

(3)

1

(4)

1

T,X

∂˜a ∂A ∂˜a ∂B ∂˜a ∂C ∂˜a ∂D ∂˜a ∂E ∂˜a ∂F ∂˜a ∂G ∂˜a ∂H + + + + + + + ∂A ∂ρ ∂B ∂ρ ∂C ∂ρ ∂D ∂ρ ∂E ∂ρ ∂F ∂ρ ∂G ∂ρ ∂H ∂ρ

= T,X

4 9  πD 

6τ

8   ∂˜a

 ∂β  T,ρ,βj=`ı

 ∂˜a  ∂ρ

+ T,β

8   ∂˜a

∂βi

(1) T,X

P ρ= ZRT

2C3 /D2 A + (3BC/D) − (C3 /D2 ) (3BC/D)ζ − (C3 /D2 ) + + 2 1−ζ (1 − ζ) (1 − ζ)3

 ζ j−1 τ

 1

(5)

i=1 j=1

A=

NC

Xi mi (di )0

(6)

1

Xi mi (di )1

(7)

1

(8)

1

(9)

1

(10)

1

(11)

1

(12)

1

(13)

NC*NC

(14)

NC

(15)

NC*NC

(16)

NC

(17)

NC

(18)

NC

∂˜a = −ln(1 − ζ) ∂A

(19)

1

∂˜a (3C/D)ζ = ∂B 1−ζ

(20)

1

(3B/D)ζ − (3C2 /D2 ) ∂˜a 3C2 3C2 /D2 = + 2 ln(1 − ζ) + ∂C 1−ζ D (1 − ζ)2

(21)

1

i=1

B=

NC i=1 NC

C=



Xi mi (di )2

i=1

D=

NC i=1

E=m=

Xi mi (di )3 NC

Xi mi i=1

G=

u kT

NC NC

u = kT

i=1

j=1

Xi Xj mi mj [uij /kT ](v0 )ij

NC NC i=1

j=1

Xi Xj mi mj (v0 )ij

uij = (1 − kij )(ui uj )1/2 , k



ui = u0i 1 + (v0 )ij = v0i =

1 2



ei 1 , k T 1/3

1/3



di = σi 1 − CI exp

 6τ πNAv

v00 i

3 ,

i = 1, NC; j = 1, NC

i = 1, NC



σi =

i = 1, NC

+ (v0 )j ]

[(v0 )i

πNAv 3 d , 6τ i

i = 1, NC; j = 1, NC

1/3 ,

3u0 − i kT

 ,

i = 1, NC

i = 1, NC

16

R. Gani et al. / Fluid Phase Equilibria 250 (2006) 1–32

Table 4c (Continued ) Equations (A/D)ζ + (C3 /D3 )(2 − ζ) (3BC/D2 )ζ 2 − (C3 /D3 )(2 + ζ) ∂˜a (2C3 /D3 )ζ = + + ∂D 1−ζ (1 − ζ)2 (1 − ζ)3 2C3 ln(1 − ζ) + D3

−

 πρ  4

9

E

6τ

jDij Gi

Number of equations 1

 ζ j−1

(22)

τ

i=1 j=1

 j

∂˜a ζ Dij Gi = ∂E τ 4

9

(23)

1

(24)

1

(25)

1

=1

(26)

1

=0

(27)

1

=0

(28)

1

=0

(29)

1

=0

(30)

1

=0

(31)

1

(32)

1

i=1 j=1

∂˜a =1 ∂F

 j

∂˜a ζ iDij Gi−1 =E ∂G τ 4

9

i=1 j=1

∂˜a ∂H ∂A ∂ρ ∂B ∂ρ ∂C ∂ρ ∂D ∂ρ ∂E ∂ρ

∂F 1 ∂gii = Xi (1 − mi ) ∂ρ gii ∂ρ NC



i=1

di dj 3di dj ξ2 1 + +2 1 − ξ3 di + dj (1 − ξ3 )2 di + dj i = 1, NC; j = 1, NC ∂G =0 ∂ρ

2

gij (dij )seg ≈ gij (dij )hs =

ξ22 (1 − ξ3 )3

NC*NC

, (33) (34)

1

(35)

1

(36)

NCA*NS

(37)

NCA*NCA

(38)

NC*NC

εAi Bj =

(39)

NCA*NCA

κAi Bj

(40)

NCA*NCA

∂H Xi = ∂ρ NCA



NS  1

XAi

i=1

XAi = ⎣1 + ρ

NCA NS i=1

∂ρ

 ∂XAi ∂ρ

⎤−1 Xj XBj Ai Bj ⎦

Bj



Ai Bj = gij (dij )seg exp

∂XAi

1 2

Ai

⎡

σij =

−

σi + σj , 2

 εAi Bj  kT

i = 1, NCA

,



− 1 (σij )3 κAi Bj ,

i = 1, NCA; j = 1, NCA

i = 1, NC; j = 1, NC

εAi + εBj , i = 1, NCA; j = 1, NCA 2 √ A B = κ i κ j , i = 1, NCA; j = 1, NCA

⎡ ⎤ NCA NCA NCA NS NS NS B A B j i j ∂X ∂ ⎦, = −(XAi )2 ⎣ Xj XBj Ai Bj + ρ Xj Xj XBj Ai Bj + ρ ∂ρ

j=1 Bj

i = 1, NCA ∂ Ai Bj ∂ρ

∂gij = ∂ρ

= (σij )3 κ

π 6

 Ai Bj ∂gij ∂ρ

exp

 εAi Bj  kT

 3di dj + di + dj (1 − ζ3 )2 D

i = 1, NC; j = 1, NC



−1 ,

(1 − ζ3 )2

j=1 Bj

(41)



C

NCA*NS

∂ρ

j=1 Bj

i = 1, NCA; j = 1, NCA

+

2Dζ2 (1 − ζ3 )3



 +2

3di dj di + dj

2 

(42)

2Cζ2 (1 − ζ3 )3

+

2D(ζ2 )2 (1 − ζ3 )4

NCA*NCA

 , (43)

NC*NC

R. Gani et al. / Fluid Phase Equilibria 250 (2006) 1–32

17

Table 4c (Continued ) Equations

Number of equations

πNAv ρ Xi mi (di )k , 6 NC

ξk =

 ∂(N a˜ ) 

k = 1, 4

(44)

4

(45)

NC

(46)

NC

i=1

ln ϕi =

∂ni

 ∂(N a˜ )  ∂ni

+ Z − 1 − ln Z,

i = 1, NC

ρ,T,nj=`ı

= a˜ +

 ∂(˜a)  ∂Xi

ρ,T,nj=`ı

− ρ,T,Xk=`ı

  NC ∂(˜a) Xj

∂Xj



j=1 C3 /D2

(3BC/D)ζ − (C3 /D2 ) ares = + + RT 1−ζ (1 − ζ)2

a˜ =

+F +E

4 9

Dij Gi

 ζ j τ

i = 1, NC

,



ρ,T,Xk=j

C3 − A ln(1 − ζ) D2

1

+H

(47)

i=1 j=1

achain Xi (1 − mi ) ln gii (di ) = RT NC

F=

i=1 NC

H=

aassoc Xi = RT



NS 

i=1

 ∂(˜a)  ∂Xi

= ρ,T,Xk=`ı

∂A = mk , ∂Xk

ln XAi −

1 Ai X 2

(48)

1

(49)

1

(50)

NC

(51)

NC

(52)

NC

(53)

NC

(54)

NC

(55)

NC

(56)

NC

(57)

NC

(58)

NC

(59)

NCA*NS*NC

(60)

NCA*NCA*NC



+

1 Mi 2

Ai

8 ∂˜a ∂βl l=1

∂βl ∂Xi

,

βl = parameters A to H; i = 1, NC

k = 1, NC

∂B = mk dk , k = 1, NC ∂Xk ∂C = mk (dk )2 , k = 1, NC ∂Xk ∂D = mk (dk )3 , k = 1, NC ∂Xk ∂E = mk , k = 1, NC ∂Xk

∂F 1 ∂gii = (1 − mk ) ln gkk + Xi (1 − mi ) , ∂Xk gii ∂Xk NC

∂G = ∂Xk ∂H = ∂Xk ∂XAi ∂Xk

2

NC



j=1

NS 

Xj mk mj [ukj /kT ](v0 )kj − (G) 2

NC NC i=1

ln X

Ak

⎡

Ak

= −(XAi )2 ⎣ρ

XAk − 2



NS

6

1 + Mk 2

+



NC

Xi



mk (dk )3 (1 − ζ3 )

2

+

kT



Xj



−1 ,

mk (dk )2 (1 − ζ3 )

NS  1

XAi

NC NS

 εAi Bj 

3di dj di + dj

,

Ak

∂XBj

j=1 Bj



Xj mk mj (v0 )kj

Xi Xj mi mj (v0 )ij

XBk Ai Bk + ρ

Bk

 πρ 

j=1



j=1

i=1

∂gij ∂ Ai Bj = (σij )3 κAi Bj exp ∂Xk ∂Xk ∂gij = ∂Xk

k = 1, NC

  NC

i=1

2

∂Xk

k = 1, NC



1 ∂XAi − + 2 ∂Xk

Ai Bj + ρ

NC NS

k = 1, NC

,

⎤

Xj XBj

∂ Ai Bj ∂Xk

j=1 Bj

⎦ , k = 1, NC

i = 1, NCA; j = 1, NCA; k = 1, NC +

2mk (dk )3 ζ2 (1 − ζ3 )

3



 +2

3di dj di + dj

2 

2mk (dk )2 ζ2 (1 − ζ3 )

3

+

2mk (dk )3 (ζ2 )

 2

NC*NC*NC

,

(1 − ζ3 )4

i = 1, NC; j = 1, NC; k = 1, NC (61) Total number of equations = 35 + 15NC + 5NC*NC + 2NCA*NS + 4NCA*NCA + NCA*NCA*NC + NCA*NS*NC + NC*NC*NC



Total number of variables = 21 Z, r, P, R, T, t, CI, NAv , 16

 ∂˜a

∂˜a ∂˜a ∂˜a ∂˜a ∂˜a ∂A , ∂B , ∂C , ∂D , ∂E , ∂F

, A, B, C, D, E, F, G, H, a˜ ,

∂˜a ∂˜a ∂A ∂B ∂C ∂D ∂E ∂F ∂G , ∂H , ∂ρ , ∂ρ , ∂ρ , ∂ρ , ∂ρ , ∂ρ



,

 ∂H

∂G ∂ρ , ∂ρ

 ∂˜a 

∂ρ T,X ,

8  ∂˜a  

i=1 ∂β T,ρ,βj=`ı

 ∂βi  ∂ρ

T,X

,

 ∂˜a   ∂ρ T,β

+ 4[ξ] + 36[D] + 19NC u- 0 , ε- /k, v- 0 , v- 00 , σ- , d- , m - , ϕ- , - ,X

+

 ∂(N a˜ ) 

, , ∂A , ∂B , ∂C , ∂D , ∂E , ∂F , ∂G , ∂H + 3NCA [Mi , εAi , κAi ] + 2NCA ∗ NS[XAi , ∂XAi /∂r] + 6NC ∗ NC[u- /k, v- 0 , g, σ- , k- , ∂g/∂ρ] + 4NCA ∗ ∂X - ρ,T,X- k=i ∂X- ∂X- ∂X- ∂X- ∂X- ∂X- ∂X- ∂X] NCA[ΔAB , εAB , κAB , ∂ΔAB /∂ρ] + NCA ∗ NCA ∗ NC[∂ΔAB /∂X] + NCA ∗ NS ∗ NC[∂XA /∂X] + NC ∗ NC ∗ NC[∂g/∂X - Note: ζ = ζ 3 , from Eq. (44).

 ∂(˜a) 

,

u kT

∂X -

ρ,T,nj=i

18

Table 4d Incidence matrix of SAFT EOS model equations Equations

Unknown variables σ* *

d* *

v- 0

* *

u-

u- /k

* *

*

* *

v- 0

u/k

* *

*

* * * *

A

B

C

D

E

F

G

ζ3

ζ2

∂a/∂A

∂a/∂B

∂a/∂C

∂a/∂D

∂a/∂E

∂a/∂F

∂a/∂G

∂a/∂H

g

∂g/∂ρ, ∂gij /∂ρ -

εAB

σ-

κAB

AB

XA-

* *

*

*

*

* * * * * *

*

* *

* *

*

* *

* * *

* * *

*

* *

* * * * *

*

*

* * * * * *

*

* *

* *

*

* *

*

* *

* * *

*

* * *

*

*

*

*

*

*

*

*

*

* *

*

* *

*

*

*

*

*

*

*

* * *

*

* *

*

* *

* * * * *

*

* * *

*

*

*

*

*

*

*

*

*

*

*

*

*

*

R. Gani et al. / Fluid Phase Equilibria 250 (2006) 1–32

18 17 16 14 13 15 12 6 7 8 9 10 11 44 44b 19 20 21 22 23 24 25 26 33 43 39 38 40 37 36 41 42 27 28 29 30 31 32 34 35 5 4 3 1 2 48 49 61 60 59 51 52 53 54 55 56 57 58 50 47 46 45

Equations Unknown variables A ∂XA- /∂ρ ∂ AB /∂ρ ∂A/∂ρ ∂B/∂ρ ∂C/∂ρ ∂D/∂ρ ∂E/∂ρ ∂F/∂ρ ∂G/∂ρ ∂H/∂ρ ∂a/∂ρ (∂a/∂β∂β/∂ρ)T,β ∂a/∂ρ Z ρ F H ∂g/∂X ∂ AB /∂X - ∂X - /∂X - ∂A/∂X - ∂B/∂X - ∂C/∂X - ∂D/∂X - ∂E/∂X - ∂F/∂X - ∂G/∂X - ∂H/∂X - ∂a/∂X - a ∂Na/∂n- ϕ- -

18 17 16 14 13 15 12 6 7 8 9 10 11 44

*

44b 19

*

20 21 *

R. Gani et al. / Fluid Phase Equilibria 250 (2006) 1–32

22 23 24 25 26 33 43 39 38 40 37 36 41 42

* *

* *

27 28

* *

29 30

* *

31 32 34 35 5 4 3 1

* * * *

* * *

*

*

*

*

*

*

*

* *

*

* *

2 48

* * * * *

49 61

*

*

60 59

*

* *

* *

*

51 52

* *

53 54

* *

55 56

* *

57 58

* *

50 47

* *

*

*

*

*

*

*

*

*

* *

* *

*

* * *

*

19

46 45

*

20

R. Gani et al. / Fluid Phase Equilibria 250 (2006) 1–32

Table 5a The Wilson model equations for calculation of liquid phase activity coefficients Model equations δ1−1 = δ1−2

V 

 −A

1

1−1

Number of equations



exp V1 T V   −A  2 1−2 = exp V1 T

δ2−1 = δ2−2 =

V  1

V2

V  2

V2

exp exp

 −A



T



2−1

 −A

2−2

1 1 1 1

T

E1 = x1 δ1−1 + x2 δ1−2

1

E2 = x1 δ2−1 + x2 δ2−2

1

ln γ1 = 1 − ln(E1 ) − (x1 δ1−1 + x2 δ2−1 )

1

ln γ2 = 1 − ln(E2 ) − (x1 δ1−2 + x2 δ2−2 )

1

Total number of equations = 8 Total number of variables = 17 [T, δ1−1 , δ1−2 , δ2−1 , δ2−2 , V1 , V2 , A1−1 , A2−1 , A1−2 , A2−2 , E1 , E2 , x1 , x2 , γ 1 , γ 2 ]

Steps 3–5: For the calculation of the liquid phase activity coefficients, the ordered incidence matrix for the set of unknown variables are highlighted in Table 6b. Note that a lower tridiagonal structure could not be found. Generally, it seems that if one adds a term to the EOS, which is dependent on density or volume, the solution is not straightforward. The GC-Flory model uses the EOS to get the molar volumes (pure compound and mixture) but not at the system pressure. Note that the other terms do not have pressure as a variable. The EOS is used to obtain the smallest pressure at which a liquid-like root is obtained. Example 6 in Appendix B highlights the implementation of the above calculation order for a binary mixture of methyl ethyl ketone/polybutadiene. 2.5. Extended UNIQUAC model for electrolyte systems The extended UNIQUAC model [18,19] is analyzed for activity coefficients of the solvent and ionic species in a solution at a specified T and solvent and salt composition, z- . The system has of NSOLV solvents and NC − NSOLV salt compounds. The salts fully dissociate via stoichiometric coefficients, νij , to form NI ions, each having an ionic charge of Zi . The total number of species in the system is NSP = NSOLV + NI. There are parameters for all pure and binary species interactions. The electrostatic contributions to the activity coefficients are identified by the Debye–H¨uckel model. Step 1: The extended UNIQUAC model equations given by Rasmussen and co-workers [18,19] are listed in Table 7a. Step 2: The DOF is found to be: 4 + 3NSP + NC*(NI + 1) + 2NSP*NSP, where the variables to be specified are classified as: 1 + NC [T, z- ] variables fixed by the problem, NI [Z - ] plus NSOLV [M - ] + NC ∗ NI [ν- ] + one [e] variables fixed by the system, two [z, b] variables fixed by the property model and 2NSP [r- , q] plus 2NSP*NSP [u0 , ut ] variables that must be retrieved from the model parameter tables or regressed. In principle, some of the charges, Z - , could be determined from the stoichiometric coefficients, [ν- ], and the requirement of charge neutrality for each of the NC − NSOLV salts, NC i=1 νij Zi = 0. In practice all of them are input for the system compounds with the user determining charge neutrality. Table 5b Wilson model equations ordered into a lower tri-diagonal form Equations

Unknown variables δ1−1

1 2

δ1−2

δ2−1

δ2−2

E1

E2

ln γ 1

ln γ 2

* *

3 4

* *

5 6

*

7 8

*

* *

* *

* *

* *

*

* *

*

R. Gani et al. / Fluid Phase Equilibria 250 (2006) 1–32

21

Table 5c List of the UNIFAC-VLE model equations for calculation of liquid phase activity coefficients Equations

Number of equations

ln γi = ln γiC + ln γiR ,

i = 1, NC

ln γiC = 1 − Ji + ln Ji − 5qi

1 − J

i

 J  i

+ ln , i = 1, NC Li Li  s   ϑk ski ki ln γiR = qi (1 − ln Li ) − − Gki ln , i = 1, NC; k = 1, NSG ηk ηk k qi



, i = 1, NC; j = 1, NC (xj qj ) ri Ji =  , i = 1, NC; j = 1, NC (x r ) j j j Li =

(1)

NC

(2)

NC

(3)

NC

(4)

NC

(5)

NC

(6)

NSG

(7)

NSG*NC

(8)

NMG*NMG

(9)

NSG

j



ηk =

(xi ski ),



i = 1, NC; k = 1, NSG

i

ski =

τmn = exp ϑk =



i = 1, NC; k = 1, NSG; m = 1, NMG

(Gmi τmk ),

m

 −a  mn

,

T

(xi Gki ),

m = 1, NMG; n = 1, NMG

i = 1, NC; k = 1, NSG

i

Gki = qi =



l

i = 1, NC; l = 1, NSG; k = 1, NSG

(vli Ql ),

(10)

NSG*NC

(vli Ql ),

l = 1, NSG; i = 1, NC

(11)

NC

(vli Rl ),

l = 1, NSG; i = 1, NC

(12)

NC

l

ri =

l

Total number of equations = 7NC + 2NSG + 2NSG*NC + NMG*NMG C R Total number of variables = 1 [T ] + 8NC [x- , L , R, η, ϑ] + 3NSG ∗ NC[s- , G - , J- , q- , r- , γ- , γ- , γ- ] + 4NSG [Q - , v- ] + 2NMG ∗ NMG[a- , τ- ] - - - Table 5d Incidence matrix for the UNIFAC-VLE model equations for calculation of activity coefficient for compound i at T, x- and model parameters Equations

Unknown variables r-

12 11 10 9 8 7 6 5 4 3 2 1

q -

G -

ϑ-

* *

*

τ-

s-

η -

* *

*

J-

L -

γR -

* * *

*

γC -

γ -

* *

*

* *

* *

* *

* * * *

*

*

*

* *

*

Step 3–5: For the calculation of the specified problem in Step 1, the ordered incidence matrix for the set of unknown variables are highlighted in Table 7b. It can be noted that a lower tridiagonal structure is obtained. Example 7 in Appendix B highlights the implementation of the above calculation order for a binary mixture of water/sodium chloride. 3. Use of property models within another model Three examples, which represent a wide range of problems where property models are used within other models, are given in this section. The first problem deals with the saturation point calculation (where an equilibrium condition needs to be satisfied in addition to the solution of the property model equations), the second problem deals with the regression of property model parameters and the third problem deals with a two-phase flash separation problem (where the property model equations, together with the equilibrium condition and mass balance equations need to be solved). These methods have been described in a somewhat different fashion by

22

R. Gani et al. / Fluid Phase Equilibria 250 (2006) 1–32

Table 6a List of the GC-Flory EoS model equations for calculation of liquid phase activity coefficients Equations

Number of equations

ln Ωi = ln Ωicomb + ln Ωifv + ln Ωiattr , ϕi ϕi ln Ωicomb = ln + 1 − , i = 1, wi xi xi v∗i Φi = , i = 1, NC NC x v∗ j=1 j j

i = 1, NC NC



SG

v∗i = (1.448)(15.17)

i = 1, NC

vn Rn ,



n=1

ln Ωifv = 3(1 + Ci ) ln

Ci =

SG



vn CT0 ,n + CT,n

n=1

ln Ωiattr

1/3

v˜ i − 1 v˜ 1/3 − 1



1 = zqi 2

 − Ci ln

1

1 − T T0

v˜ i , v˜

i = 1, NC

SG 

+

SG

SG

m=1

n=1

Rm

Cn0

i = 1, NC

,

1 θj exp(− εji /RT ) − [εii (˜v) − εii (˜vi )] + 1 − ln RT NC

NC

θ exp(− εji /RT ) εji

j=1 j NC θ k=1 k



i = 1, NC

vn Qn ,

NC

(2)

NC

(3)

NC

(4)

NC

(5)

NC

(6)

NC

(7)

NC

(8)

NC

(9)

NC*NC

Rn

j=1

qi =

(1)

 ,

exp(− εki /RT )

i = 1, NC

n=1

εji (˜v) = εii (˜v) = εji (˜vi ) =

ε0ji v˜

ε0ii , v˜ 0 εji v˜ i



ε0ji = θn(i) =

i = 1, NC; j = 1, NC − 1

,

i = 1, NC i = 1, NC; j = 1, NC

,

(i) θm



m (i) vn Qn

qi

(j)

θn εnm ,

i = 1, NC; j = 1, NC

NC

(11)

NC*NC

(12)

NC*NC

(13)

NC*NSG

(14)

NMG*NMG

(15)

NC

(16)

NC*NC

(17)

NC

(18)

1

(19)

1

(20)

1

(21)

NC

(22)

1

n

vn Qn , qi

=

i = 1, NC; j = 1, NSG

εnm = −[εmm εnn ]1/2 + εnm , x qi

i

θi =

(10)

m = 1, NMG; n = 1, NMG

i = 1, NC

,

xj qj j

εji = εji − εii = εji (˜v) − εii (˜v), 1/3

Pi =

nRT (˜vi

1/3 v∗i v˜ i (˜vi

+ Ci ) − 1)

1 ni Ci = xi Ci n

C=

v∗ =

NC

NC

i=1

i=1

NC

xj v∗j

j=1

E

(˜vi ) =

NC 1

2



E

(˜v) =

NC 1

2 i=1

 zqi ni

NC

εii (vi ) +

zqi ni

i=1 attr

i = 1, NC

Eattr (˜v) nRT (˜v1/3 + C) + v∗ v˜ v∗ v˜ (˜v1/3 − 1)

P=

attr

Eattr (˜vi ) , v∗i v˜ i

+

i = 1, NC; j = 1, NC

εii (v) +

θ exp(− εji /RT ) εji

j=1 j NC θ k=1 k



NC 

, exp(− εki /RT )

θ exp(− εji /RT ) εji

j=1 j NC θ k=1 k

 i = 1, NC



exp(− εki /RT )

Total number of equations = 12NC + 4NC*NC + NC*NSG + NMG*NMG + 4 comb fv attr ∗ attr ,Φ Total number of variables = 10 [˜v, P, Eattr (˜v), C, v∗ , T0 , T, z, R, n] + 16NC [Ω - ,Ω - , v- , Ω - , C- , Ω - , v˜- , q- , ε- (˜v), θ- , P- , E - (˜vi ), x- , n- , w -] 0 0 + 5NSG [R ] + 4NC ∗ NC[ε- (˜v), ε- (˜v- ), ε- , Δε] + NC ∗ NSG[θ- ] + 2NMG ∗ NMG[Δε, ε- ] - , C- T0 , C- T , C- , Q -

R. Gani et al. / Fluid Phase Equilibria 250 (2006) 1–32

23

Table 6b Incidence matrix of GC-Flory EoS model equations Equations

4 3 8 13 15 14 12 16 11 10 6 21 9 17 20 19 22 18 2 5 7 1

Unknown variables v∗i

ϕi

* *

*

(i)

qi

θn

* * *

*

θi

εji

ε0ji

* *

*

εji

εji (˜vi )

εii (˜v)

Ci

Eattr (˜vi )

εji (˜v)

v˜ i

v∗

C

Eattr (˜v)

v˜

Ωicomb

Ωifv

Ωiattr

Ωi

* *

*

* *

* * *

*

*

*

* *

* *

*

*

*

*

*

*

*

* *

*

*

* * *

* *

*

*

*

* *

*

* *

*

*

* * *

*

*

*

*

*

*

* *

*

Seader and Henley [22] and by Westerberg et al. [23], where the property model related equations are represented only by one equation (as a procedure equation). This is valid as long as the property model equations can be ordered in a lower triangular form. 3.1. Saturation point calculation Considering a closed two-phase system where the phase equilibrium (isofugacity criteria) is represented by the gamma–phi approach (at low enough pressure to ignore the Poynting correction), xi γi Pisat ϕisatV = yi PϕiV ,

i = 1, NC

(2)

or, the phase equilibrium (isofugacity criteria) is represented by the phi–phi approach (at any temperature and pressure with the same fugacity coefficient model for both phases), xi ϕiL = yi ϕiV ,

i = 1, NC

(3)

From the above equations, it becomes clear that in order to predict VLE, we need to predict first γ i using a model for the liquid phase activity coefficients and Pisat using a model for pure compound vapor pressures and a model for ϕiV and ϕisatV using an EOS model for vapor phase fugacity coefficients or an EOS model for both ϕiL and ϕiV . We will now analyze the phase equilibrium model equations by inserting property models and illustrating the forms of the ordered set of equations displayed in an incidence matrix. This analysis is more detailed than those of Seader and Henley [22] and of Westerberg et al. [23]. Consider a liquid that is non-ideal, that is, γ i = 1 in equilibrium with an ideal vapor (ϕiV = 1 = ϕisatV ). Using the UNIFAC-VLE [22] model (Section 2.3.2) for the liquid phase activity coefficients (γ i ) and the Antoine model for component vapor pressures ln Pisat =

Ai − Bi , (T + Ci )

i = 1, NC

(4)

We need to add the following two equations to satisfy the condition of phase equilibrium (at least for a local solution), Residy = Σi yi − 1 = 0,

i = 1, NC

(5)

Note that if Eqs. (2) and (4) are inserted into Eq. (5), we will have one implicit equation with one unknown variable, T. However, we prefer to add one more equation that relates the sum in Eq. (5) to T. T = f (T, Σi yi , Residy ),

i = 1, NC

(6)

The total number of equations for this problem is given by the model equations listed in Table 5b plus Eqs. (2), (4)–(6). All these equations become the ordered set in the incidence matrix of Table 8a (note that we omit the known variables A -,B - , C- , x- and P from

24

R. Gani et al. / Fluid Phase Equilibria 250 (2006) 1–32

Table 7a The extended UNIQUAC model equations and variables Equation

Number of equations Debye–H¨uckel

ln γw = ln γwcomb + ln ␥residual + ln γw w ln γi∗

=

ln γicomb Φi =

− ln γicomb,∞

ln γicomb = ln

Φ  xi

xi ri

xi qi

i = 1, NSP

(2)

NI

(3)

NSP

,

i = 1, NSP

(5)

NSP

(6)

NSP

(7)

NSP*NSP

(8)

NSP*NSP

(9)

NSOLV

NSP



T

k=1

∗,Debye–H¨uckel

θk ψik

NSP

θψ l=1 l lk

,

i = 1, NSP

k = 1, NSP; i = 1, NSP

,

 2A b3

= −Zi2

NSP k=1

uki = u0ki + utki (T − 298.15), = Mw

−

θk ψki

 u −u  ki ii

ψki = exp −

ln γi

i = 1, NI

NSOLV

NSP

= qi 1 − ln

Debye–H¨uckel



Φi +1− , θi

(1)

(4)



ln γw

θi

∗,Debye–H¨uckel + ln γi ,

i = 1, NSP

xq i=1 i i

ln γiresidual

− ln γiresidual,∞

,

xr i=1 i i

NSP

θi =

w = 1, NSOLV

,

 Φ  i

Φi z +1− − qi ln xi 2

i

NSP

+ ln γiresidual

k = 1, NSP; i = 1, NSP

1 − 2 ln(1 + bI 1/2 ) , 1 + bI 1/2

l + bI 1/2 −

AI 1/2 , 1 + bI 1/2

1 mi Zi2 2



w = 1, NSOLV

i = 1, NI

(10)

NI

(11)

1

(12)

NI

(13)

1

(14)

NI

(15)

NI

(16)

NI

(17)

NSOLV

NI

I=



i=1



NSOLV NI xs + x s=1 i=1 i NSOLV xs Ms s=1





mj = xj

j = 1, NI

,

A = 1.131 + 1.335e − 3(T − 273.15) + 1.164e − 5(T − 273.15)2 r   r q   ri z ri qw i i w ln γicomb,∞ = ln +1− − qi ln +1− , rw rw 2 rw qi rw qi ln γiresidual,∞ = qi [1 − ln ψwi − ψiw ],

i = 1, NI

i = 1, NI

NC

xi =

NSOLV w=1

xw =

ν z j=1 ij j NI zw + k=1

NSOLV w=1



zw zw +

NC

NI

x j=1 j

,

ν z j=1 kj j ,

i = 1, NI

w = 1, NSOLV

Total number of equations = 2 + 4NSP + 6NI + 3NSOLV + 2NSP*NSP comb , γ residual ] + 6NI [Z, γ ∗ , γ Debye–H¨uckel , γ comb,∞ , γ residual,∞ , m] Total number of variables = 6 [T, z, A, e, b, I] + 7NSP [x- , Φ - , θ- , r- , q- , γ- - Debye – H¨ u ckel ] + NC[z- ] + NC ∗ NI[ν- ] + 4NSP ∗ NSP[ψ, u- , u- 0 , u- t ] + 3NSOLV [M - , γ- , γ-

the incidence matrix). Eqs. (2), (4)–(6) add 2NC + 2 equations and 2NC [y, P- sat ] plus 2 [Residy , T] new variables. As shown in Table 8a, a lower tridiagonal structure is not obtained. The upper off-diagonal variable is the saturation temperature T. This means that in every step of an iterative solution, the value for T needs to be found. Consequently, all other related variables and equations within the loop, will need to be solved again. Let us now consider the use of the phi–phi approach with Eq. (3) representing the isofugacity criteria using an EOS model that allows for both liquid and vapor phases to be non-ideal (ϕiL (T, P, x- ) = 1 = ϕiV (T, P, y) along with Eq. (3) and Eqs. (5) and (6). Again, we specify x- and P. It is assumed that all the EOS model parameters are available. As shown in Table 8b, we have 3NC + 2 L V equations with 3NC + 2 unknown variables [ϕ , ϕ , y, Residy , T ], which are listed in the last five columns of row 1 in the incidence - - matrix. Now, the incidence matrix shows additional unknown variables (the vapor compositions) on the upper tri-diagonal part of the incidence matrix. This means that in addition to T, the y also will be found at each iteration. A detailed model analysis of an EOS in terms of calculation of the fugacity coefficient has already been given previously for several EOS models and therefore is not repeated here. Note that the property model is within the iterated set of equations. The consequence is that those EOS models (such as the PC-SAFT, CPA, GC-Flory, etc.), which are not lower tridiagonal in structure, will have additional iteration loops. For these models, a simultaneous equation solution mode may be more efficient.

R. Gani et al. / Fluid Phase Equilibria 250 (2006) 1–32

25

Table 7b Incidence matrix of the extended UNIQUAC model equations Equations

Unknown variables x-

16/17 8 7 5 4 13 12 11 14 15 10 3 6 2 9 1

U-

ψ -

* *

*

θ-

A

Φ -

m -

I

* *

*

γ comb,∞ -

γ ∗,DH -

γ residual,∞ -

γ residual -

γ∗ -

*

* *

*

*

*

γ comb -

γ Debye–H¨uckel -

γ -

* *

*

* * * *

* *

* * * *

*

*

*

*

* * *

*

*

*

3.2. Regression of property model parameters Regression of property model parameters involves those variables that cannot be retrieved from property model parameter tables in order to use a property model for a calculation. Consider, for example, Table 8a (where the UNIFAC model is used for calculation of saturation temperatures). If ND experimental VLE data sets of P, Texp , x- and yexp are available, a property model parameter regression problem can be formulated, such as exp 2

Minimize Fobj =

exp 2

w1 [Σk (Tk − Tk ) ] w2 [Σk {(y1 )k − (y1k )} ] + , ND ND

k = 1, ND

(7)

Table 8a Incidence matrix for VLE saturation temperature calculation model (gamma–phi approach with ideal vapor phase) Equations

A -

τ-

a- = 0) τ- = f (a- , T ) r- = f (ν- , R -) q = f (ν- , Q) G - = f (ν, Q) J- = f (x- , r- ) L - = f (x- , q- ) ln γ C = f (J- , L - , z, q- ) ϑ- = f (x- , G ) s- = f (τ- , G -) η = f (x- , s- ) ln γ R = f (L - , q- , G - , ϑ- , s- , η- ) ln γ = ln γ C + ln γ R 2 4 5 6

* *

*

f (a0 , T, T

r-

q -

G -

J-

ln γ C -

L -

s-

ϑ-

η -

ln γ R -

ln γ -

P- sat

y -

Residy

T * *

* * * *

* * * * *

*

*

*

* *

* *

*

* * *

* *

* *

* *

* *

* *

* * *

* *

*

Table 8b Incidence matrix for VLE calculation model (phi–phi approach with an EOS) Equations

ϕL (EOS) ϕV (EOS) 3 5 6

Known variables

Unknown variables

x-

P

ϕL -

*

*

*

* *

ϕV -

y -

*

*

*

* *

*

*

*

Residy

T * *

*

26

R. Gani et al. / Fluid Phase Equilibria 250 (2006) 1–32

Table 9a Incidence matrix for a PT-flash model with gamma–phi approach Equations

Known variables F

Unknown variables P

z-

4 12 9 2 5 β= 11

*

T

P- sat

* *

*

* *

* *

γ -

x-

*

*

* *

* *

y -

* *

V

L

* * *

V F

β

* *

* *

*

subject to the optimization variables, a- 0 (matrix of UNIFAC group intraction parameters), and the constraints of isofugacity, Eqs. (2). The data can be weighted by the factors w1 and w2 . The constraints are typically those of Table 8a, but alternative objective functions to Eq. (7) require the same treatment. Thus the NMG*NMG variables a- 0 , which were considered known in Table 8a, are not available. Adding Eq. (7) and a- 0 into the calculations leads to an extra iteration loop. In the inner-loop, the saturation point calculation is performed for each data-point while in the outer-loop, the estimates for the parameters a- 0 are generated based on the values of the objective function. The result is the addition of ND*NC + 1 equations and NMG*NMG unknown variables. Since this is an optimization problem, the only requirement is that the total number of equations be greater than the total number of unknown variables. An alternative regression is to obtain the first derivatives of Eq. (7) with respect to each of the parameters a- 0 , setting these derivatives to zero (condition for local mimima) and solving this set of NMG*NMG equations for the NMG*NMG variables a- 0 together with the ND*NC equations of saturation point calculation model. Either way, parameter regression will add upper off-diagonal elements to the incidence matrix structure, requiring additional iteration loops (in a sequential calculation) or a simultaneous solution approach (where, depending on the number and location of upper off-diagonal elements, a sparse solver may be useful). 3.3. Two-phase separation model Now let us consider a typical two-phase PT-flash calculation involving vapor and liquid phases by displaying the incidence matrices for the gamma–phi and phi–phi approaches. In both cases, we must introduce mass balance equations. Fzi = Vyi + Lxi ,

i = 1, NC

(8)

Eq. (8) can be rearranged to the following form, where β = V/F. zi = βyi + (1 − β)xi ,

i = 1, NC

(9)

In addition to Eq. (5), the liquid compositions x are constrained to sum to unity. Residx = Σi xi − 1,

i = 1, NC

(10)

Since the Residx and Residy both must be zero at solution of the model equations, this means that inserting Eqs. (5) and (10) into Eq. (8) results in the total mass balance equation, F =V +L

(11)

For the PT-flash problem described above, the model equations are those listed in Fig. 2 plus Eqs. (8)–(10) as shown in Table 9a where we represent the activity coefficients by a generic expression since they have a tridiagonal form. γ i = f (T, x- ) (12) The result is 4NC + 3 equations in 5NC + 6 variables. We specify NC + 3 of these: F, z- , T, and P. For purposes of illustration, let us consider a binary mixture. There are 11 unknown variables: β, L, V, x- , y, γ , and P- sat to be solved - by the 11 equations of the PT-flash model represented by Eqs. (2), (4), (5), (9), (11), (12) plus β = V/F. From the above analysis, it is clear that a sub-set of equations (Eqs. (12), (9), (2), (5)) need to be solved simultaneously for γ , x- , y and β. Also, by assuming values for the two variables lying in the upper tridiagonal (x- and β), an iterative solution scheme can be generated. Note that in this case not all the property models belong to the equations within the repetitive cycle of equations (or the set of simultaneous equations). The same analysis as above can also be shown using an EOS for the phi–phi approach where the fugacity coefficient model equations are used twice. The two-phase flash model is given by Eqs. (3), (5), (8), (11) and β = V/F. The unknown variables are, ϕL , ϕV , x- , y, β, V and L. The incidence matrix in Table 9b shows that initial estimates are necessary for x- , y and β. Also, while only -

R. Gani et al. / Fluid Phase Equilibria 250 (2006) 1–32

27

Table 9b Incidence matrix for a PT-flash model with phi–phi approach Equations

Known variables F

ϕL (EOS) ϕV (EOS) 9 3 5 β=

z-

Unknown variables P

T

ϕL -

*

*

*

*

*

* *

ϕV -

x-

* *

* *

* *

*

11

L

*

* V F

V

β

* *

* *

y -

* *

*

*

*

*

the equations for the fugacity coefficient model plus Eqs. (3), (5), (9) need to be solved simultaneously, all of the property model equations belong to the set of equations within the repetitive cycle of equations (or the set of simultaneous equations). This means much greater computational effort. 4. Conclusions General analyses of the mathematical and numerical elements of several well-known property models have been presented together with numerical verification of the solution of the model equations. The analysis is based on a systematic methodology for organizing the property model equations into ordered, square (incidence) matrices of equations and unknown variables. In this paper, the total set of equations and variables are given together with an ordered set of independent equations and a corresponding sample calculation a variety of property models. The approach can provide a thorough understanding of the property models in terms of the variables set by the problem, the variables set by the system compounds, and the variables set by the property model or found by parameter table retrieval or regression. In addition, the analysis identifies a sequencing of the property model equations, often into lower tridiagonal form, which suggests efficiencies in solving the equations. It is found that simple property model equations such as the cubic EOS or UNIFAC can often be ordered into lower tridiagonal form, complex property model equations such as the SAFT and CPA EOS cannot, demonstrating the need for simultaneous solution (or an additional iteration-loop) solution methods. We believe this exposition will help readers develop their own calculation routines and appreciate more fully the nature of property model equations and what information, including data, is needed to apply them. Our intention in current and future work is to investigate how property model equations might be modified, simplified and/or manipulated for specific applications and more useful roles in product and process synthesis, design and analysis. Appendix A. Association models See Tables A1 and A2. Table A1 Monomer fractions for different bonding types [16] Type

approximations

1

AA = 0

2A

AA = AB = BB = 0

2B

XA approximations

XA =

−1+(1+4ρ )1/2 2ρ

X A = XB

=

−1+(1+8ρ )1/2 4ρ

AA = BB = 0, AB = 0

XA = X B

=

−1+(1+4ρ )1/2 2ρ

3A

AA = AB = BB = AC = BC = CC = 0

X A = X B = XC

=

−1+(1+12ρ )1/2 6ρ

3B

AA = AB = BB = CC = 0, AC = BC = 0

XA = XB , XC = 2XA − 1

=

−(1−ρ )+((1+ρ )2 +4ρ ) 6ρ

4A

AA = AB = BB = AC = BC = CC = AD = BD = CD = DD = 0

X A = X B = XC = XD

=

−1+(1+16ρ )1/2 8ρ

4B

AA = AB = BB = AC = BC = CC = DD = 0, AD = BD = CD = 0

XA = XB = XC , XD = 3XA − 2

=

−(1−2ρ )+((1+2ρ )2 +4ρ ) 6ρ

4C

AA = AB = BB = CC = CD = DD = 0, AC = AD = BC = BD = 0

XA = XB = XC = XD

=

−1+(1+8ρ )1/2 4ρ

1/2

1/2

28

R. Gani et al. / Fluid Phase Equilibria 250 (2006) 1–32

Table A2 Types of bonding in real associating fluids [16] Species

Rigorous type

Assigned type

Acid Alkanol

1 3B

1 2B

Water Amines tertiary

4C 1

3B Non self-associating

Amines secondary Amines primary

2B 3B

2B 3B

Ammonia

4B

3B

Appendix B. Test examples This appendix contains sample solutions of property model equations as shown in the incidence matrix tables. The cases and results given here have been solved through ICAS-MoT [24] and then creating a COM-object. Such objects are available from R. Gani as COM-objects that can be run from EXCEL since once the incidence matrix has been made, the equations can be solved through any appropriate external software. We have given relatively simple problems that have been reported in the open literature so that the numerical solutions can be straightforwardly be verified (Table B1). B.1. Equation of state: Soave–Redlich–Kwong (SRK) For the binary system n-decane(1)/ethane at 373 K and 1 atm, find ϕ1 and ϕ2 when x1 = 0.3 and x2 = 0.7. B.1.1. Input of fixed variables System: n-decane(1)/ethane(2); variables fixed by the problem (NC + 2 = 4): [T, P, x1 , x2 ] = [373.0, 1.0, 0.3, 0.7]; variables fixed by the system (3NC = 6): [ωi , Tci , Pci ] i 1 2

Tc (K)

Pc (atm)

ω

617.700 305.320

20.8240 48.0830

0.4923 0.0995

Variables fixed by the property model (3): [ψA , ψB , R] = [0.42747, 0.08664, 82.0575]; variables regressed or retrieved from model parameter tables (NC*NC = 4) [k11 = 0 = k22 ; k12 = k21 = 0]. B.1.2. Desired output using equations in Table 2c as in the incidence matrix of Table 2d Variables calculated: [ϕ1 , ϕ2 , Z] = [0.9428, 1.0012, 0.9831]. B.2. Equation of state: CPA EOS For the binary system methanol/methane at 373 K and 1 atm, find ϕ1 and ϕ2 when x1 = 0.5 and x2 = 0.5. Table B1 Problem descriptions for the tested property models Binary system

Fugacity coefficient ϕ1

ϕ2

Example 1: SRK Example 2: CPA

n-Decane/ethane Methanol/methane

0.9428 0.9638

1.0012 1.0108

Example 3: PC-SAFT Example 4: SAFT

n-Decane/ethane Methane/hexadecane

0.63198 1.4963

1.0160 9.996 × 10−8

Binary system

Activity coefficient γ1

γ2

Example 5: UNIFAC Example 6: GC-Flory

Diethylamine/n-heptane Methyl ethyl ketone/polybutadiene

1.1330 0.0256

1.047 1.000

Example 7: extended UNIQUAC

Water/sodium chloride

1.002

0.2953 (Na+ ), 1.4482 (Cl− )

R. Gani et al. / Fluid Phase Equilibria 250 (2006) 1–32

29

B.2.1. Input of fixed variables System: methanol(1)/methane(2); variables fixed by the problem (3 + 2NC = ): [T, P, n, x- , n- ] = [373.0, 1.0, 1.0, 0.5, 0.5, 0.5, 0.5]; variables fixed by the system (3NC = 6): [T- c , P- c , ω- ]. - i 1 2

Tci (K)

Pci (atm)

ωi

512.640 190.564

79.911 45.389

0.5640 0.0115

Variables fixed by the property model (3) [R, Ωa , Ωb ] = [82.0575, 0.42747, 0.08664]; variables regressed or retrieved from model parameter tables (5NCA + NC ∗ NC + NS ∗ NCA = 10) [a- 0 , b- , C11 , εAB , β, k- , n- A ]. a01 (atm cm6 /mol)

b1 (cm3 /mol)

C11

3 εAB 1 (atm cm /mol)

β1

4.0001 × 106

30.980

0.431

242690

0.0161

k11 = k12 = k21 = k22 = 0; nA1 = 1 B.2.2. Output using equations in Table 3c as in the incidence matrix of Table 3d Variables calculated: [ϕ, Z] = [0.9638, 1.0108, 0.9870]. B.3. Equation of state: PC-SAFT EOS For the binary system n-decane/ethane at 511 K and 20 bar, find ϕ1 and ϕ2 when x1 = 0.3 and x2 = 0.7. B.3.1. Input of fixed variables System: n-decane(1)/ethane(2); variables fixed by the problem: (NC + 2 = 4) [T, P, x1 , x2 ] = [511.0, 20.0, 0.3, 0.7]; Variables fixed by the system (3NC = 6) [m - , σ- , ε- /k]. Name

mi

˚ σ i (A)

εr /k (K)

n-Decane Ethane

4.6627 1.6069

3.8384 3.5206

243.87 191.42

Variables fixed by the property model (3Ja + 3Jb = 42) [a- , b- ]. i

a0i

a1i

a2i

0 1

0.910563 0.636128

−0.308402 0.186053

−0.090615 0.452784

2 3

2.686134 −26.547362

−2.503005 21.419794

0.596270 −1.724183

4 5

97.759209 −159.591541

−65.255885 83.318680

−4.130211 13.776632

6

91.297774

−33.746923

−33.746923

i

b0i

b1i

b2i

0 1

0.724095 2.238279

−0.575550 0.699510

0.097688 −0.255757

2 3

−4.002585 −21.003577

3.892567 −17.215472

−9.155856 20.642076

4 5

26.855641 206.551338

192.672264 −161.826462

−38.804430 93.626774

6

−355.602356

−165.207693

−29.666906

Variables regressed or retrieved from model parameter tables (NC*NC = 4) [k11 = 0 = k22 ; k12 = k21 = 0]. B.3.2. Output using equations in Table 4a as in the incidence matrix of Table 4b Variables calculated: [ϕ, Z] = [0.63198, 1.0160, 0.871]. -

30

R. Gani et al. / Fluid Phase Equilibria 250 (2006) 1–32

B.4. Equation of state: SAFT EOS For the binary system methane/n-hexadecane at 300 K and 13 MPa, find ϕ1 and ϕ2 when x1 = 0.7 and x2 = 0.3. B.4.1. Input of fixed variables System: methane(1)/hexadecane(2); variables fixed by the problem (2 + NC = 4): [T, P, x- ] = [300.0, 13.0, 0.7, 0.3]; variables fixed by the system (NCA = 0); variables fixed by the property model (=40) [D - , τ, NAv , CI, R]. i

D0i

D1i

D2i

D3i

0 1

−8.8043 4.1646270

2.9396 −6.0865383

−2.8225 4.7600148

0.3400 −3.1875014

2 3

−48.203555 140.43620

40.137956 −76.230797

11.257177 −66.382743

12.231796 −12.110681

4 5

−195.23339 113.51500

−133.70055 860.25349

69.248785 0.0

0.0 0.0

6 7

0.0 0.0

−1535.3224 1221.4261

0.0 0.0

0.0 0.0

8

0.0

−409.10539

0.0

0.0

τ = 0.7405; NAv = 6.022136736 E +023; CI = 0.12; R = 8.31441.

0 Variables regressed or retrieved from model parameter tables (4NC + 2NCA + NC*NC = 12) [v- 00 , m - , u- /k, k- ].

Name

ν00 (mL/mol)

m

u0 /k (K)

Methane Hexadecane

21.576 12.3

1.0 11.209

190.29 210.65

[k11 = k22 = 0; k12 = k21 = 0.118]. B.4.2. Output using equations in Table 4c as in the incidence matrix of Table 4d Variables calculated: [ϕ1 , ϕ2 , Z] = [1.4963, 9.996 × 10−8 , 0.6438]. B.5. Excess Gibbs energy model: UNIFAC-VLE For the binary system diethylamine/n-heptane at 308.15 K, find γ 1 and γ 2 when x1 = 0.4 and x2 = 0.6. B.5.1. Input of fixed variables System: diethylamine(1)/n-heptane(2); variables fixed by the problem (1 + NC = 3): [T, x] = [308.15, 0.4, 0.6]; variables fixed by the system (NSG*NC = 6) [ν- ]. Name

Group formula

Diethylamine n-Heptane

CH3 CH2 NH CH2 CH3 CH3 (CH2 )5 (CH3 )

Subgroup i

ν1i

ν2i

CH3 CH2

2 1

2 5

CH2 NH

1

0

Variables regressed or retrieved from model parameter tables (2NSG + NMG*NMG/2 = 8) [Q, R , a]. - - Subgroup

Main group

Ri

Qi

CH3 CH2

CH2 CH2

0.9011 0.6744

0.8480 0.5400

CH2 NH

CNH

1.2070

0.936

aCH2 –CNH = 65.33 K; aCNH–CH2 = 255.7 K.

R. Gani et al. / Fluid Phase Equilibria 250 (2006) 1–32

31

B.5.2. Output using equations in Table 5c as in the incidence matrix of Table 5d Variables calculated: [γ 1 , γ 2 ] = [1.1330, 1.047]. B.6. Equation of state: GC-Flory EOS For the binary system methylethylketone (MEK)/polybutadiene (PBD) at 373 K, find γ 1 and γ 2 when x1 = 1 × 10−7 and x2 ∼ 1. B.6.1. Input of fixed variables −7 System: MEK (1)/PBD (2); variables fixed by the problem [3NC + 2 = 8]: [T, x- , n, n- , w - ] = [373.0, 1 × 10 , −7 −10 0.9999999, 1, 1 × 10 , 0.9999999, 3 × 10 , 1.0]; variables fixed by the property model [4 + NC = 6]: [T0 , z, R, P, P- ] = (i) (i) [298.15, 10.0, 1.9872, 0.0, 0.0, 0.0]; variables fixed by the system (NC*NSG = 8): v1 , v2 . Compound

Group formula

MEK PBD

CH3 CH2 CO CH3 [CH2 CH CH CH2 ]

(i)

(i)

Subgroup

v1

v2

CH3 CH2

1 1

0 2

CH3 CO CH CH

1 0

0 1

Variables regressed or retrieved from model parameter tables [=29]: [C- T0 , C- T , C- 0 , ε- , R ] -,Q Subgroup

C- T0 ,n

CT,n

Cn0

CH3 CH2

−0.0738 0.108

−3.57 −3.57

0 0

−4.117 −15.01

0 0

CH3 CO CH = CH

0.126 0.1762

εmn (cal/q-unit)

CH2

CH2 C O

C C

CH2 CH2 C O

−544.0 110.0

110.0 −1080.0

−16.2 30.0

−16.2

30.0

−630.0

C C

Subgroup

Main group

Ri

Qi

CH3 CH2

CH2 CH2

0.9011 0.6744

0.848 0.540

CH3 CO CH CH

CH2 C O C C

1.6742 1.1167

1.488 0.867

Output: variables calculated: [γ 1 , γ 2 ] = [0.0256, 1.0]. B.7. Activity coefficient model: extended UNIQUAC model For the binary system water(1)/sodium chloride(2) at 298.15 K, find the activity coefficients of all neutral and ionic species γ i for i = 1–3 using the extended UNIQUAC model for z1 = 0.98 and z2 = 0.02. B.7.1. Input of fixed variables System: water(1)/sodium chloride(2); variables fixed by the problem [1 + NC = 3]: [T, z- ] = [298.15 K, 0.98, 0.02]; variables fixed by the system [=33] [M - , z, b, x- , uT ] = [0.01801528, 10.0, 1.5].

32

R. Gani et al. / Fluid Phase Equilibria 250 (2006) 1–32

i

Species

xi

0 1

H2 O Na+

0.9608 0.0196

2

Cl−

0.0196

uT0i 0 0.4872 14.631

uT1i

uT2i

0 0

8.5455 20.278

0

13.628

uT3i

uT4i

0.4872 0 15.635

14.631 15.635 14.436

Variables regressed or retrieved from model parameter tables [=39] [u- 0 , r- , q, Z ] - i

u00i 0 1

0 733.286

2

1523.39

i

u01i

u02i

100,000 1 × 1010

600.495 1398.14

1 × 1010

1895.52

ri 0 1

0.92 1.4034

2

10.386

qi

u03i

u04i

733.286 0

1523.39 1443.23

1443.23

2214.81

Zi

1.4 1.199

– 1

10.197

−1

Output: variables calculated: [γ 0 , γ 1 , γ 2 ] = [1.002, 0.2953, 1.4482]. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24]

R. Gani, J.P. O’Connell, Comp. Chem. Eng. 25 (2000) 3–14. S. Macchietto, E.H. Chimowitz, T.F. Andersen, L.F. Stutzman, IEC Proc. Des. Dev. 25 (1986) 674–682. J. Perregaard, E.L. Sorensen, Comp. Chem. Eng. 16S (1992) 247–254. S. Storen, T. Hertzberg, Comp. Chem. Eng. 21 (Suppl. S) (1997) S709–S714. M.L. Michelsen, IEC Res. 24 (1986) 184–188. M.L. Michelsen, Fluid Phase Equil. 60 (1990) 213–219. E. Hendriks, IEC Res. 27 (1988) 1732–1736. C.F. Leibovici, Fluid Phase Equil. 84 (1993) 1–8. M. Pons, The CAPE-OPEN interface specification for physical properties databanks, in: Proceedings of the 16th Conference on Thermophysical Properties, London, UK, September, 2002. G. Soave, Chem. Eng. Sci. 27 (1972) 1197–1203. G.M. Kontogeorgis, E.C. Voutsas, I.V. Yakoumis, D. Tassios, IEC Res. 35 (1996) 4310–4318. G.M. Kontogeorgis, I.V. Yakoumis, H. Meijer, E.M. Hendriks, Fluid Phase Equil. 158–160 (1999) 201–211. G.M. Kontogeorgis, Association models—the CPA equation of state, in: G.M. Kontogeorgis, R. Gani (Eds.), Computer Aided Property Estimation for Product and Process Design, CACE-19, Elsevier B.V., 2004, pp. 113–142 (Chapter 6). J. Gross, G. Sadowski, IEC Res. 40 (2001) 1244–1260. J. Gross, G. Sadowski, IEC Res. 41 (2002) 1084–1093. S.H. Huang, M. Radosz, IEC Res. 30 (1991) 1994–2005. G. Bogdanic, Aa. Fredenslund, IEC Res. 33 (1994) 1331–1333. H. Nicolaison, P. Rasmussen, J.M. Sorensen, Chem. Eng. Sci. 48 (1993) 3149–3158. K. Thomsen, P. Rasmussen, R. Gani, Chem. Eng. Sci. 51 (1996) 3675–3683. G.M. Wilson, J. Am. Chem. Soc. 86 (1964) 127–130. Aa. Fredenslund, P. Rasmussen, J. Gmehling, Vapor–Liquid Equilibria Using UNIFAC, Elsevier Scientific Publishing Company, 1977. J.D. Seader, E.J. Henley, Separation Process Principles, Wiley, New York, USA, 1998, pp. 253–263 (Chapter 5). A.W. Westerberg, H.P. Hutchison, R.L. Motard, P. Winter, Process Flowsheeting, Cambridge University Press, Cambridge, UK, 1979, pp. 27–101 (Chapter 3). M. Sales-Cruz, R. Gani, A modelling tool for different stages of the process life, in: S.P. Asprey, S. Macchietto (Eds.), Computer-Aided Chemical Engineering, vol. 16: Dynamic Model Development, Elsevier, Amsterdam, 2003, pp. 209–249.