This book is written in LaTeX, a complete typesetting language, and
set in the standard Times-Roman 11 point font. It is also provided
with ChemSep in ASCII text form (file CHEMSEP.TXT) for
online reference which was generated with a LaTeX to ASCII converter.
The conversion is limited, with the result that the ASCII text file
contains some unconverted LaTeX formatting. A PostScript file
(BOOK.PS) also generated by LaTeX can be downloaded from our ftp site.
Harry Kooijman
Ross Taylor
We would like to express our appreciation to Professor Hans Wesselingh (now at the University of Groningen, the Netherlands) who initially promoted the project and made various resources available and encouraged us by letting students use the program for their course work. This has been an indispensable source of feedback that has helped us to improve the program. We also like to thank Peter Verheijen for his enthusiasm and contributions in the early years of the project. Also, various students have worked on projects to check and improve the programs and documentation, which was very helpful. Finally, ChemSep owes its very existence to the Internet which enabled the authors to keep in touch and continue development while living on different continents.
F(x) = 0where F is a vector consisting of all the equations to be solved and x is, again, the vector of variables. A Taylor series expansion of the function vector around the point x_o which the functions are evaluated gives (ignoring second and higher order terms):
F(x) = F(x ) + J(x - x )
o o
where J is the Jacobian matrix of partial derivatives of F
with respect to the independent variables x:
J = (d F / d x ) ij i jIf x is the actual solution to the system of equations, then F(x)=0 and we can rewrite the above equation as:
J(x - x ) = - F(x )
o o
This linear system of equations may be solved for a new
estimate of the vector x. If the new vector, x, obtained in this way
does not actually satisfy the set of equations, F, then the procedure
can be repeated using the calculated x as a new x. The entire
procedure is summarized below.
An important drawback of the Newton's method can be its sensitivity to the initial guess, x_o, since quadratic convergence is only achieved ``close" to the solution. In order to obtain convergence, Newton's method requires that reasonable initial estimates be provided for all independent variables. It is obviously impractical to expect the user of a SC method to guess this number of quantities. Thus, the designer of a computer code implementing a SC method must provide one or more methods of generating initial estimates of all the unknown variables. Several techniques have been developed to improve the convergence away from the solution and to prevent the method from taking a step in a wrong direction. The simplest and most common technique is to ``damp" the step of each variable to some range or fraction of the Newton step. However, this damping also reduces the methods effectiveness.
Simultaneous correction procedures have shown themselves to be generally fast and reliable, having a locally quadratic convergence rate in the case of Newton's method, and these methods are much less sensitive to difficulties associated with nonideal problems than are tearing methods. They also lend themselves to be easier extended with optimization, parametric sensitivity, or continuation methods.
0 = (1-t) (F(X ) - F(X)) + t F(X)
o
where t varies from 0 (where X=X_o) to 1 (where F(X)=0). Better
continuation methods can be formulated while using a parameter which
has some physical significance. In separations problems the most
appropriate choice would be the degree to which mass transfer (between
the present phases) prevails. In the equilibrium model the stage
efficiency and in the nonequilibrium model the mass transfer rates
represent this degree. Thus, they will be multiplied with a parameter
t which will vary from 0 (no separation at all) to 1 (actual
separation).
K = (P / P)
i vap,i
L V
K = (phi / phi )
i i i
where the fugacity coefficients are calculated from an equation
of state. This model is recommended for separations involving
mixtures of hydrocarbons and light gases (hydrogen,
carbondioxide, nitrogen, etc.) at low and high pressures. It is
not recommended for nonideal chemical mixtures at low
pressures. The EOS must be able to predict vapour as well as
liquid fugacity coefficients.
* * V
K = (gamma phi P PF / phi P)
i i i i i i
This option should be used when dealing with nonideal fluid
mixtures. It should not be selected for separations at high
pressures.
*
K = (gamma P / P)
i i i
This is the form of the K-value model used in the DECHEMA
compilations of equilibrium data (Hence the name given to this
menu option). DECHEMA uses the Antoine equation to compute the
vapour pressures but ChemSep allows you to
choose other vapour pressure models if you wish. This option
should be used when dealing with non-ideal fluid mixtures. It
should not be selected for separations at high pressures.
1/m_i 2 3 4
K = A + B T + C T + D T + E T
i i i i i i
You must supply the coefficients A through E and the exponent
m.
c
theta = V / sum x V
i i k k k
* c 2
ln gamma = (V / R T) [ delta - sum x delta theta ]
i i i j j j j
*
ln gamma = ln gamma + ln theta + 1 - theta
i i i i
where delta_i is called the solubility parameter and V_i the
molal volume of component i (both read from the PCD-file). This
regular model is also incorporated in the Chao-Seader method of
estimating K-values.
2
ln gamma = [ A + 2 ( A - A ) x ] x
i ij ji ij i j
It can only be used for binary mixtures (i=1, j=2 and i=2,
j=1).
2
ln gamma = (A / (1+(A x / A x )) )
i ij ij i ji j
It can only be used for binary mixtures (i=1, j=2 and i=2,
j=1).
Lambda = (V /V ) exp ( -(mu -mu )/RT )
ij j i ij ii
c
S = sum x Lambda
i j=1 j ij
c
ln gamma = - ln (S ) - sum x Lambda / S
i i k=1 k ki k
The two interaction parameters are (mu_ij-mu_ii) and
(mu_ji-mu_ii) per binary pair of components.
tau = (g -g ) / RT
ij ij ii
G = exp (-alpha tau )
ij ij ij
c
S = sum x G
i j=1 j ji
c
C = sum x G tau
i j=1 j ji ji
c
ln gamma = C / S + sum x G ( tau - C / S ) / S
i i i k=1 k ik ik k k k
The interaction parameters are (g_ij-g_ii), (g_ji-g_ii), and
alpha_ij per binary (only one alpha is required as
alpha_ij=alpha_ji).
c
r = sum x r
i=1 i i
c
q = sum x q
i=1 i i
phi = x r / r
i i
theta = x r / r
i i
tau = exp ( -(lambda -lambda )/RT)
ji ji ii
c
S = sum theta tau
i j=1 j ji
c
ln gamma = (1 - (z / 2) q ) ln (( phi / x )) + (z / 2) q ln ((theta / x )) - (r / r) + (z / 2) q ((r / r) - (q / q) )
i i i i i i i i i i
r c
ln gamma = q ( 1 - ln (S ) - sum (theta tau / S ) )
i i i k=1 k ik k
c r
ln gamma = ln gamma + ln gamma
i i i
The interaction parameters are (lambda_ij-lambda_ii) and
(lambda_ji-lambda_ii) per binary. The parameters r_i and q_i
are read from the component database (PCD file).
*
ln P = A - (B / T + C )
i i i i
Note the natural logarithm. This option should be selected if
you are using activity coefficient models with parameters from
the DECHEMA series. Antoine parameters are available in the
ChemSep data files and need not be loaded.
* G
ln P = A + (B / C + T) + D T + E ln T + F T
i i i i i i i i
The parameters A through G must be supplied by the user. A
library of parameters for some common chemicals is provided
with ChemSep in the file EANTOINE.LIB.
* E
ln P = A + (B / T) + D T + C ln T + D T
i i i i i i i
DIPPR parameters A--E are also available in
ChemSep data files.
6
zeta = 36 / T + 96.7 log T - 35 - T
T r r r
6
zeta = 36 / T + 96.7 log T - 35 - T
Tb rb rb rb
phi = 0.118 zeta - 7 log T
T r
psi = 0.0364 zeta - log T
T r
alpha = (0.136 zeta + log P - 5.01 / 0.0364 zeta - log T )
Tb c Tb rb
*
log P = - phi - (alpha-7) psi
r
* (0) (1)
ln P = f + w f
i i
(0) 6
f = 5.92714 - (6.09648 / T ) - 1.28862 ln (T ) + 0.169347 T
r r r
(1) 6
f = 15.2518 - (15.6875 / T ) - 13.4721 ln (T ) + 0.43577 T
r r r
where T_r=T/T_Ci. Both the Riedel and Lee-Kesler models are
recommended for hydrocarbon mixtures in particular.
P = (R T / V)The Virial and cubic EOS are discussed in the sections below.
P = (R T / V) + (B R T / V)
The method of Tsonopoulous for estimating virial coefficients
is recommended for hydrocarbon mixtures at low pressures. It is
based on an earlier correlation due to Pitzer.
c c
B = sum sum y y B
i=1 j=1 i j ij
(0) (0)
B = (R T P ) ( B + w B )
ij c,ij c,ij ij ij ij
(0) 2 3 8
B = 0.1445 - (0.33 / T ) - (0.1385 / T ) - (0.0121 / T ) - (0.000607 / T )
ij r r r r
(1) 2 3 8
B = 0.0637 + (0.331 / T ) - (0.423 / T ) (-0.0008 / T )
ij r r r
w = (w + w / 2)
ij i j
Z = (Z + Z / 2)
c,ij ci cj
1/3 1/3 1/3
V = (V V / 2)
c,ij ci cj
T = (1 - k ) sqrt( T T )
c,ij ij ci cj
P = (Z R T / V )
c,ij c,ij c,ij c,ij
Binary interaction parameters k_ij must be supplied by the
user. For paraffins k_ij can be calculated with:
1/3 1/3 3
k = 1 - (8 sqrt(V V ) / (V +Vc ) )
ij ci cj ci cj
2
P = (R T / V - b) - (a / V )
with
2 2
a = (27 R T / 64 P )
i ci ci
b = (R T / 8 P )
i ci ci
and the mixing rules:
c c
a = sum sum y y a
i=1 j=1 i j ij
a = sqrt( a a )
ij i j
c
b = sum y b
i=1 i i
P = (R T / V - b) - (a / sqrt(T) V (V + b))
with
2 2.5
a = (Omega R T / P )
i a ci ci
Omega = 0.42748
a
b = (Omega R T / P )
i b ci ci
Omega = 0.08664
b
and the mixing rules:
c c
a = sum sum y y a
i=1 j=1 i j ij
a = (1 - k ) sqrt( a a )
ij ij i j
c
b = sum y b
i=1 i i
where k_ij is a binary interaction parameter (original RK:
k_ij=0).
P = (R T / V - b) - (a / V (V + b))
with
a = a (T ) alpha(T ,w )
i i ci ri i
2 2
a (T ) = (Omega R T / P )
i ci a ci ci
Omega = 0.42747
a
2 2
alpha(T ,w ) = [ 1 + (0.480 + 1.574 w - 0.176 w ) (1-sqrt(T )) ]
ri i i i ri
b = (Omega R T / P )
i b ci ci
Omega = 0.08664
b
and the mixing rules:
c c
a = sum sum y y a
i=1 j=1 i j ij
a = (1 - k ) sqrt( a a )
ij ij i j
c
b = sum y b
i=1 i i
2 2
alpha(T ,w ) = [ 1 + (0.48508 + 1.55171 w - 0.15613 w ) (1 - sqrt(T )) ]
ri i i i ri
and specially for hydrogen:
-0.30288 T_{ri}
alpha(T ,w ) = 1.202 e
ri i
P = (R T / V - b) - (a / V (V + b) + b (V - b))
with
a = a (T ) alpha(T ,w )
i i ci ri i
2 2
a (T ) = (Omega R T / P )
i ci a ci ci
Omega = 0.45724
a
2 2
alpha(T ,w ) = [ 1 + (0.37464 + 1.5422 w - 0.26992 w ) (1 - sqrt(T )) ]
ri i i i r
b = (Omega R T / P )
i b ci ci
Omega = 0.07880
b
and the mixing rules:
c c
a = sum sum y y a
i=1 j=1 i j ij
a = (1 - k ) sqrt( a a )
ij ij i j
c
b = sum y b
i=1 i i
2 3
H = A + B T + C T + D T
i i i i i
You must enter the coefficients A through E in the "Load Data"
option of the Properties menu for vapour and liquid enthalpy
for each component.
Table 2.1. Temperature correlations
| Equation | Parameter(s) | Formula |
| number | ||
| 2 | A,B | A + B T |
| 3 | A-C | A + B T + C T^2 |
| 4 | A-D | A + B T + C T^2 + D T^3 |
| 10 | A-C | exp ( A - B / C+T ) |
| 100 | A-E | A + B T + C T^2 + D T^3 + E T^4 |
| 101 | A-E | exp ( A + B / T + C ln T + D T^E ) |
| 102 | A-D | A T^B / 1 + C/T + D/T^2 |
| 103 | A-D | A + B exp ( -C / T^D ) |
| 104 | A-E | A + B / T + C / T^3 + D / T^8 + E / T^9 |
| 105 | A-D | A / B^(1 + (1 - T / C)^D) |
| 106 | A-E | A (1-T_r)^(B + C T_r + D T_r^2 + E T_r^3 ) |
| 107 | A-E | A + B ( C / T/sinh (C / T)^2 + D ( D / T/cosh (D / T)^2 |
Table 2.2. Component properties with the typical correlation number
| Liquid density | 105 |
| Vapour pressure | 101 |
| Heat of vaporisation | 106 |
| Liquid heat capacity | 100 |
| Ideal gas heat capacity | 107 |
| Second virial coefficient | 104 |
| Liquid viscosity | 101 |
| Vapour viscosity | 102 |
| Liquid thermal conductivity | 100 |
| Vapour thermal conductivity | 102 |
| Surface tension | 106 |
| Ideal gas heat capacity (Reid Prausnitz and Poling) | 4 |
| Antoine | 10 |
| Liquid viscosity (Reid, Prausnitz and Sherwood) | 2 |
Table 2.3. Default physical property correlations
| Mixture liquid density | Rackett |
| Component liquid density | Polynomial |
| Vapour density | Cubic EOS |
| Mixture liquid viscosity | Molar averaging |
| Component vapour viscosity | Polynomial/Letsou-Stiel |
| Mixture vapour viscosity | Brokaw |
| Component vapour viscosity | Polynomial |
| Mixture liquid thermal conductivity | Molar average |
| Component liquid thermal conductivity | Polynomial |
| Mixture vapour thermal conductivity | Molar average |
| Component vapour thermal conductivity | Polynomial/9B-3 |
| Liquid diffusivity | Kooijman-Taylor/Wilke-Chang |
| Vapour diffusivity | Fuller et al. |
| Mixture surface tension | Molar average |
| Component surface tension | Polynomial |
| Liquid-liquid interfacial tension | Jufu et al. |
c
T = sum x T
c,m i=1 i c,i
c
V = sum x V
c,m i=1 i c,i
c
Z = sum x Z
c,m i=1 i c,i
P = Z R T / V
c,m c,m c,m c,m
c
M = sum x M
m i=1 i i
which will be referred as the "normal" mixing rules. Reduced
properties will be calculated by:
T = T / T r c P = P / P r c V = V / V r cunless specified otherwise.
L c L
(1 / rho ) = sum (x / rho )
m i=1 i i
where the component liquid densities, rho^L_i, are computed as
discussed below.
c
T = sum x T
c,m i=1 i c,i
c
Z = sum x Z
R,m i=1 i R,i
T = T / T
r c,m
(1+(1-T_r)^{2 / 7})
F = Z
z R,m
c
A = sum (x T / M P )
i=1 i c,i i c,i
L c
rho = (1 / A R F sum x M )
m z i=1 i i
If the reduced temperature, T_r, is greater than unity a
default value of 50 kmol/m^3 is used.
1/3
T = (1-T )
* r
2 3
A = 17.4425 -214.578 Z +989.625 * Z -1522.06 Z
c c c
2 3
Z <= 0.26: B = -3.28257 +13.6377 Z +107.4844 Z -384.211 Z
c c c c
2 3
Z > 0.26: B = 60.20901 -402.063 Z +501 Z + 641 Z
c c c c
L 2 4
rho = (1 + A T + B T + (0.93 - B) T / V )
m * * * c
* c * c * 2/3 c * 1/3
V = (1/ 4) ( sum x V + 3 (sum x (V ) ) (sum x (V ) ) )
m i=1 i i i=1 i i i=1 i i
c c * *
T = (sum sum x x V T / V )
c,m i=i j=i i j ij c,ij m
c
w = sum x w
SRK,m i=1 i SRK,i
Z = 0.291 - 0.08 w
c,m SRK,i
*
P = Z R T / V
c,m c,m c,m m
If the reduced temperature is larger than unity a default value
of 50 kmol/m^3 is used, otherwise the saturated liquid volume
(V_s) is calculated from:
* (0) (delta)
(V / V ) = V (1 - w V )
s m R SRK,m R
(0) 1/3 2/3 4/3
V = 1 + a (1-T ) + b (1-T ) + c (1-T ) + d (1-T )
R r r r r
(delta) 2 3
V = (e + f T + g T + h T / (T - 1.00001))
R r r r r
where
| a=-1.52816 | e=-0.296123 |
| b= 1.43907 | f= 0.386914 |
| c=-0.81446 | g=-0.0427258 |
| d= 0.190454 | h=-0.0480645 |
For the density of compressed liquids the saturated liquid volume is corrected (Thomson et al., AIChE J, 28, 671, 1982):
V = V ( 1 - c ln (beta + P / beta + P ) )
s vpm
1/3 2/3 4/3
beta / P = -1 + a (1-T ) + b (1-T ) + d (1-T ) + e (1-T )
c r r r r
2
e = exp ( f + g w + h w
SRK,m SRK,m
c = j + k w
SRK
where
| a=-9.070217 | g= 0.250047 |
| b= 62.45326 | h= 1.14188 |
| d=-135.1102 | j= 0.0861488 |
| f= 4.79594 | k= 0.0344483 |
P = P P
vpm c,m rm
(0) (1)
log P = P + w P
rm rm SRK,m rm
(0)
P = 5.8031817 log T + 0.07608141 alpha
rm rm
(1)
P = 4.86601 log T + 0.03721754 alpha
rm rm
6
alpha = 35 - 36 / T -96.736 log T + T
rm rm rm
T = T / T
rm c,m
This method should be used for reduced temperatures from 0.25
up to the critical point.
(1+(1-T_r)^{2 / 7})
F = Z
z R
L
rho = (P / R T F )
m c c z
L ig -1 1/3 -1 -1
C - C = 1.45 + 0.45 (1-T ) + 0.25 w [ 17.11 + 25.2 (1-T ) T + 1.742 (1-T ) ]
p,i p r r r r
However, in ChemSep the temperature correlation is
used for all temperatures to prevent problems arrising from using
different liquid heat capacity methods in the same column (which
especially trouble nonequilibrium models). Liquid heat capacities
could also be computed from the selected thermodynamic models to
circumvent this problem.
L c L
ln eta = sum z ln eta
m i=1 i i
where z_i are either the mole fractions (for molar averaging, the
default) or alternatively the weight fractions for mass averaging. A
better method is from Teja and Rice (1981, A479). However, this method
requires interaction parameters. Here a different mixing rule (for
T_cijV_cij) is used which improves the model predictions with unity
interaction coefficients:
c
w = sum x w
m i=1 i i
c
M = sum x M
m i=1 i i
c c
V = sum sum x x V
cm i=1 j=1 i j cij
1/3 1/3 3
V = (( V + V ) / 8)
cij ci cj
c c
T = (sum sum x x T V / V )
cm i=1 j=1 i j cij cij cm
T V = psi (T V + T V / 2)
cij cij ij ci ci cj cj
where psi_ij is set to unity for all components. The liquid viscosity
of the mixture is computed from two reference components
ln (epsilon eta ) = ln (epsilon eta ) + [ ln (epsilon eta ) - ln (epsilon eta ) ] (( w -w / w - w ))
m m 1 1 2 2 1 1 m 1 2 1
with epsilon defined as
2/3
epsilon = (V / sqrt(T M ))
i ci ci i
and the reference component vioscosities are evaluated at T T_ci /
T_cm. Component liquid viscosities are calculated from the liquid
viscosity temperature correlation if the temperature is within the
valid range. Otherwise the component viscosity is computed with DIPPR
procedure 8G, the Letsou-Stiel method (1973, see A471):
1/6 2/3
xi = 2173.424 T / sqrt(M ) P
c,i i c,i
(0) 2 -5
xi = ( 1.5174 - 2.135 T + 0.75 T ) 10
r r
(1) 2 -5
xi = ( 4.2552 - 7.674 T + 3.4 T ) 10
r r
L (0) (1)
eta = ( xi + w xi ) / xi
i
Alternatively the simple temperature correlation given in Reid et
al. (RPS liquid viscosity, see A439) can be used:
log eta = A + B / TA high pressure correction by Lucas (A436) is used to correct the influence of the pressure on the liquid viscosity:
A
eta = (1 + D (Delta P / 2.118) / 1 +C w Delta P ) eta
r i r SL
where eta_SL is the viscosity of the saturated liquid at P_vp, and
Delta P = (P - P ) / P
r vp ci
-4 -0.03877
A = 0.9991 - [4.674 10 /(1.0523 T - 1.0513)]
r
2 3
C = - 0.07921 + 2.1616 T - 13.4040 T + 44.1706 T
r r r
4 5 6 7
- 84.8291 T + 96.1209 T - 59.8127 T + 15.6719 T
r r r r
2.573 0.2906
D = [0.3257/(1.0039-T ) ]-0.2086
r
L c L
eta = sum (x eta / sum x phi )
m i=1 i i i ij
where the interaction parameters phi_ij can be calculated by Wilke's
(1950) method:
1/4 2
phi = ((1 + sqrt( \eta / \eta ) (M / M ) ) / sqrt( 8 (1 + M / M ) ) )
ij i j i j i j
or by Brokaw's method:
phi = S A sqrt( \eta / \eta )
ij i j
1/4
sm = ( 4 / (1+M /M ) (1+M /M ) )
j i i j
0.45 0.45
A = (sm / sqrt(M /M )) ( 1 + ((M /M - (M /M ) ) / 2 (1 + M /M )) + ((1 + (M /M ) ) / sqrt(sm) (1 + M /M )) )
i j i j i j i j i j i j
If the Lennard-Jones energy parameter, epsilon (in Kelvin), and the
Stockmayers polar parameter, delta, are known, S is calculated from:
2 2
S = (1 + sqrt( (T/\epsilon ) (T/\epsilon ) ) + (delta delta /4) / sqrt(1 + T/\epsilon + \delta / 4) sqrt(1 + T/\epsilon + \delta / 4))
i j i j i i j j
otherwise it is approximated by S=1. epsilon and delta can be
estimated from:
3.6
epsilon = 65.3 T Z
c,i c,i
5 2
delta = 1.744 10 9 (mu / V T )
b,i b,i
Where mu is the dipole moment in Debye. Vapour viscosities are a
function of pressure and a correction is normally applied. Mixture
properties are computed with the "normal" mixing rules. DIPPR
procedure 8E can be used to compute the high pressure viscosity:
rho = 1 / V
c c,m
rho = rho / rho
r c
1/6 2/3
xi = 2173.4241 T / sqrt(M ) P
c,m m c,m
1.85
A = exp (1.4439 rho ) - exp (-1.111 rho )
r r
-7
B = 1.08 10 A / xi
eta = eta + B
hp
where rho is the vapour mixture molar density.
Both Wilke's and Brokaw's method require pure component viscosities. These are normally obtained from the vapour viscosity temperature correlations, as long as the temperature is within the valid temperature range. If not, then the viscosity can be computed with the Chapman-Enskog kinetic theory (see Hirschfelder et al. 1954 and A391-393):
*
T = T / epsilon
* -b * *
Omega = a (T ) + c / exp (d T ) + e / exp (f T )
v
V - 2 2 *
eta = 26.69 10 7 M T / sigma (Omega + 0.2 delta / T )
v
where the collision integral constants are a=1.16145, b=0.14874,
c=0.52487, d=0.77320, e=2.16178, and f=2.43787. The viscosity may also
be computed with the Yoon and Thodos method (DIPPR procedure 8B):
1/6 2/3
xi = 2173.4241 T / sqrt(M ) P
i c,i i c,i
V b 8
eta = (1 + a T - c exp (d T-r) + e exp (f T ) / 10 xi)
i r r
where the constants a-f are given in Table 2.4.
Table 2.4. Constants for the Yoon-Thodos method
| Hydrogen | Helium | Others |
| a=47.65 | a=52.57 | a=46.1 |
| b=0.657 | b=0.656 | b=0.618 |
| c=20.0 | c=18.9 | c=20.4 |
| d=-0.858 | d=-1.144 | d=-0.449 |
| e=19.0 | e=17.9 | e=19.4 |
| f=-3.995 | f=-5.182 | f=-4.058 |
-7 0.618
eta = 10 [0.807 T - 0.357 exp (-0.449T ) +
r r
o o
0.340 exp (-4.058 T ) + 0.018] F F / xi
r p q
3 -5 4 1/6
xi = 0.176 (( T / M (10 P ) ))
c c
where F_p^o and F_q^o are polarity and quantum correction factors. The
polarity correction depends on the reduced dipole moment:
-30 2 -5 2
mu = 52.46 ((mu / 3.336 10 ) (10 P ) / T )
r c c
If mu_r is smaller than 0.022 then the correction factor is unity,
else if it is smaller than 0.075 it is given by
o 1.72
F = 1 + 30.55 (0.292 - Z )
p c
else by
o 1.72
F = 1 + 30.55 (0.292 - Z ) [0.96 + 0.1 (T - 0.7)]
p c r
The quantum correction is only used for quantum gases He, H_2, and
D_2,
o 0.15 2 1/M
F = 1.22 Q (1 + 0.00385[(T -12) ] sign(T -12) )
q r r
where Q=1.38 (He), Q=0.76 (H_2), Q=0.52 (D_2). There is also a
specific correction for high pressures (A421) by Lucas.
o
eta = Y F F eta
p q
e f d -1
Y = 1 + (a P / b P + (1+c P ) )
r r r
o -3 o
F = (1 + (F - 1) Y / F )
p p p
o -1 4 o
F = (1 + (F - 1) [Y - 0.007 (ln Y) ] / F )
q q q
where eta^o refers to the low-pressure viscosity (note that the
original Lucas method has a different rule for Y if T_r is below
unity, however, this introduces a discontinuity which is avoided
here). The parameters a through f are evaluated with:
-3 -0.3286
a = (1.245 10 / T ) exp 5.1726 T
r r
b = a (1.6553 T - 1.2723)
r
-37.7332
c = (0.4489 / T ) exp 3.0578 T
r r
-7.6351
d = (1.7368 / T ) exp 2.2310 T
r r
e = 1.3088
0.4489
f = 0.9425 exp -0.1853 T
r
where, in case T_r is below unity, T_r is taken to be unity. For
mixtures the Lucas model uses the following mixing rules:
c
T = sum y T
cm i=1 i ci
c
V = sum y V
cm i=1 i ci
c
Z = sum y Z
cm i=1 i ci
P = R T Z / V
cm cm cm cm
c
M = sum y M
m i=1 i i
o c o
F = sum y F
pm i=1 i pi
o c o
F = A sum y F
qm i=1 i qi
where A is a correction factor depending on the molecular weights of
the components in the mixture. Let H denote the component of highest
molecular weight and L of lowest, then if M_H/M_L > 9:
0.57
A = 1 - 0.01 (M / M )
H L
else A=1. The mixture vapor viscosity is computed with the Lucas
method as for a component which has the mixture properties T_cm, P_cm,
M_m, F_pm^o, and F_qm^o. Therefore, the method is not interpolative in
the same way as the techniques of Wilke and Brokaw (that is, the
method does not necessarily lead to pure component viscosity eta_i
when all y_j=0 except y_i=1).
L c L
lambda = sum x lambda
m i=1 i i
c L
F = x / sum x / rho
v,i i i=1 i i
lambda = 2 / ( 1 / lambda + 1 / lambda )
ij i j
L c c
lambda = sum sum F F lambda
m i=1 j=1 v,i v,j ij
L c L 2
(1 / sqrt(\lambda )) = sum (w / (lambda ) )
m i=1 i i
where w_i is the weight fraction of component i.
1.2 1.4
lambda = (0.63 T P / (30 + P ) + 0.98 + 0.0079 P T ) lambda
hp r r r r r
This is DIPPR procedure 9G-1 where the mixture parameters are computed
by the "normal" mixing rules. Component liquid thermal conductivities
are calculated from one of the following methods:
2/3
f = 3 + 20 (1 - T )
r
2/3
b = 3 + 20 (1 - 273.15 / T )
c,i
-4 x L
lambda = c 10 M rho (f/b)
i i i
for straight chain hydrocarbons c=1.811 and x=1.001 else
c=4.407 and x=0.7717.
L 0.38 1/6
lambda = (A (1-T ) / T )
i r r
* alpha beta gamma
A = (A T / M T )
b i c
where parameters A^*, alpha, beta, and gamma depend on the
class of the component as shown in Table 2.5.
Table 2.5. Parameters for the Latini equation for liquid thermal conductivity
| Family | A^* | alpha | beta | gamma |
| Saturated hydrocarbons | 0.0035 | 1.2 | 0.5 | 0.167 |
| Olefins | 0.0361 | 1.2 | 1.0 | 0.167 |
| Cycloparaffins | 0.0310 | 1.2 | 1.0 | 0.167 |
| Aromatics | 0.0346 | 1.2 | 1.0 | 0.167 |
| Alcohols, phenols | 0.00339 | 1.2 | 0.5 | 0.167 |
| Acids (organic) | 0.00319 | 1.2 | 0.5 | 0.167 |
| Ketones | 0.00383 | 1.2 | 0.5 | 0.167 |
| Esters | 0.0415 | 1.2 | 1.0 | 0.167 |
| Ethers | 0.0385 | 1.2 | 1.0 | 0.167 |
| Refrigerants: | ||||
| R20, R21, R22, R23 | 0.562 | 0.0 | 0.5 | -0.167 |
| Others | 0.494 | 0.0 | 0.5 | -0.167 |
V c V
lambda = sum x lambda
m i=1 i i
V c V c
lambda = sum (x lambda / sum x phi )
m i=1 i i j=1 j ij
where interaction parameters phi_ij are computed from:
3/4 2 2
phi = 0.25 (1+sqrt((\eta \over \eta ) (M \over M ) (T+1.5 T \over T+1.5 T ))) (T+sqrt(1.5 T T ) / T+1.5 T )
ij i j j i b,i b,j b,i b,j b,i
Note that the component viscosities are required for this
evaluation.
rho = (1 / V )
c c,m
rho = (rho / rho )
r c
If the reduced density is below 0.5 then a=2.702, b=0.535, and c=-1;
if the reduced density is witin [0.5,2] then a=2.528, b=0.67, and
c=-1.069; otherwise a=0.574, b=1.155, and c=2.016. The high pressure
thermal conductivity correction is then calculated from:
-8 1/6 2/3 5
Delta lambda = (a 10 (exp (b rho ) + c) / (( sqrt(M ) T / P ) ) Z )
r m c,m c,m c,m
which must be added to the calculated thermal conductivity for low
pressure.
Pure component vapour thermal conductivities are estimated from the following methods:
V V
lambda = (1.15 (C - R) + 16903.36) eta / M
i p i i
V V
lambda = (1.3 (C - R) +14644 -2928.8/T ) eta / M
i p r i i
V V
lambda = 2.5 (C - R) eta / M
i p i i
1/6 2/3
xi = 2173.424 T / sqrt(M ) P
c,i i c,i
-7
lambda = 4.91 10 T C / xi
i r p
1/6 2/3
xi = 2173.424 T / sqrt(M ) P
c,i i c,i
-8 1/6
lambda = 11.05 10 (14.52 T -5.14) C / xi
i r p
k o
DD = D , k=i
ij ij
k o
DD = D , k=j
ij ji
k o o
DD = sqrt(D D ), k != i, k != j
ij ik jk
c k x_k
DD = sum (DD )
ij i=1 ij
Liquid binary infinitive diffusion coefficients (D^o_ij) are normally
computed by the Wilke-Chang method unless selected otherwise. The
following models are available:
o -16 0.6
D = 1.1728 10 ( sqrt(\phi M ) T / eta V )
ab b b b b,a
where phi_b is association factor for the solvent (2.26 for
water, 1.9 for methanol, 1.5 for ethanol and 1.0 for
unassociated solvents).
o -14 -1.14 -0.589
D = 8.62 10 eta V
aw w b,a
o -0.19 -13 (0.00958/V_{b,a} -1.12) 1.52
D = (3.36 V - 3.65) 10 (1000 eta ) T
aw b,a w
o -18 -0.19 0.2 -0.4 1.7
D = 4.3637 10 eta r r T
ab b a b
o -14 -071 (0.0102/V_{b,a} -0791)
D = 9.859 10 V (1000 eta )
ab b,a b
o -16 -0.5473 -1.026
D = 5.6795 10 V eta T
aw b,a w
o -15 -0.45 0.265 -0.907
D = 5.2383 10 V V eta T
ab b,a b,b b
o -12 2/3
D = 5.927 10 (T r / eta r )
ab b b a
o -12 2 1/6 0.6
D = 8.93 10 ( (V / V ) ) ( (P / P ) ) (T / eta )
ab b,a b,b b a b
This method is not yet implemented!
V -2 1.75 2
D = 1.013 10 (T sqrt((1/M +1/M )) / P (sqrt(V ) + sqrt(V )) )
ab a b a b
where V_a and V_b are the Fuller molecular diffusion volumes
which are calculated by summing the atomic contributions from
Table 2.6. This table also lists some special diffusion volumes
for simple molecules.
Table 2.6. Fuller diffusion volumes
| Atomic and structural diffusion volume increments | |||
| C | 15.9 | F | 14.7 |
| H | 2.31 | Cl | 21.0 |
| O | 6.11 | Br | 21.9 |
| N | 4.54 | I | 29.8 |
| Ring | -18.3 | S | 22.9 |
| Diffusion volumes of simple molecules | |||
| He | 2.67 | CO | 18.0 |
| Ne | 5.98 | CO2 | 26.9 |
| Ar | 16.2 | N2O | 35.9 |
| Kr | 24.5 | NH3 | 20.7 |
| Xe | 32.7 | H2O | 13.1 |
| H2 | 6.12 | SF6 | 71.3 |
| D2 | 6.84 | Cl2 | 38.2 |
| N2 | 18.5 | Br2 | 69.0 |
| O2 | 16.3 | SO2 | 41.8 |
| Air | 19.7 | ||
sigma = (sigma + sigma ) / 2
ab a b
epsilon = sqrt( \epsilon + \epsilon )
ab a b
The diffsuion collision integral is
*
T = T / epsilon
ab
* -b * * *
Omega = a (T ) + c / exp (d T ) + e / exp (f T ) + g / exp (h T )
D
where the collision integral constants are a=1.06036, b=0.1561,
c=0.193, d=0.47635, e=1.03587, f=1.52996, g=1.76474, and
h=3.89411. If Stockmayer polar parameters are known the
integral gets corrected with:
*
Omega = Omega + (0.19 delta delta / T )
D,c D a b
and the diffusion coefficient is
V 3/2
D = C (T sqrt(1/M +1/M ) / P sigma w )
ab a b ab D
where constant C=1.883 10^-2.
-2 -3
C = 2.1987 10 - 5.07 10 sqrt(1/M +1/M )
a b
c
sigma = sum x sigma
m i=1 i i
c L 2 c L L c L 2
sigma = (sum ( ((x / rho )) + sum ((x x sqrt(\sigma \sigma ) / rho rho )) ) / (sum ((x / rho )) ) )
m i=1 i i j=1 i j i j i j i=1 i i
c
T = sum x T
c,m i=1 i c,i
c
T = sum x T
b,m i=1 i b,i
c
sigma = sum x sigma
b,m i=1 i b,i
Then it corrects the sigma_b,m with:
*
T = ((1/T -1) / (1/T -1))
r rb
* 1.118091
sigma = 1.002855 (T ) (T / T ) sigma
b r
T = T / T
br b,i c,i
Q = 0.1207 ( 1 + T ( ln (P ) - 11.526 ) / (1 - T ) ) - 0.281
br c,i r
-7 2/3 1/3 11/9
sigma = 4.6 10 P T Q (1 - T )
i c,i c,i r
*
T = ((1/T -1) / (1/T -1))
r rb
* 1.118091
sigma = 1.002855 (T ) (T / T ) sigma
b r
,
sigma = sigma + sigma - 1.1 sqrt( \sigma \sigma )
1 2 1 2
This method generally overpredicts the interfacial tension for aqueous
systems. We use a general method from Jufu et al. (1986):
,, ,
X = - ln ( x + x + x )
1 2 3r
,, ,
sigma , = (K R T X / A exp (X) (x q + x q + x q ))
w0 1 1 2 2 3r 3
with A_w0=2.5 10^5 (m^2/mol), R=8.3144 (J/mol/K), K=0.9414 (-), and
q_i is the UNIQUAC surface area parameters of the components i. The
components are ordered in such a manner that component 1 and 2 are the
dominating components in the two liquid phases. Then the rest of the
components are lumped into one mole fraction, x_3. This lumped mole
fraction is taken for the phase which has the largest x_3 (the
richest). q_3 is the molar averaged q for that phase for all
components except 1 and 2.
a, b Cubic EOS parameters
B Second virial coefficient (m3/kmol)
c Number of components
C_p Mass heat capacity (J/kg.K)
D_ij Binary diffusion coefficient (m^2/s)
DD_ij Binary Maxwell-Stafan diffusion coefficient (m^2/s)
D^o_ij Infinite dilution binary diffusion coefficient (m^2/s)
K_i K-value of component i, equilibrium ratio (K_i=y_i/x_i)
k_ij Binary interaction coefficient (for EOS)
M Molecular mass (kg/kmol)
R Gas constant = 8134 (J/kmol K)
r Radius of gyration (Angstrom)
P Pressure (Pa)
P^*,P_vap Vapour pressure (Pa)
P_i Parachor (m^3 kg^1/4 / s^1/2) of component i
PF Poynting correction
q UNIQUAC surface area parameter
T Temperature (K)
T_r Reduced temperature (T_r=T/T_c)
T_b Normal boiling temperature (K)
V Molar volume (m^3/kmol)
V_b Molar volume at normal boiling point (m^3/kmol)
V_s Saturated molar volume (m^3/kmol)
w Weight fraction (of component)
x Liquid mole fraction (of component)
y Vapour mole fraction (of component)
Z Compressibility
Z_R Racket parameter
Greek:
alpha Attractive parameter in EOS
w Acentric factor
Omega_a, Omega_b EOS parameters
Omega_v Collision integral for viscosity
Omega_D Collision integral for diffusion
gamma Activity coefficient
delta Stockmayer parameter
epsilon Molecular energy parameter (K)
lambda Thermal conductivity (W/m.K)
rho Molar density (kmol/m^3)
eta Viscosity (Pa.s)
phi_i Fugacity coefficient of component i
phi_s Association factor for solvent s (Hayduk-Laudie)
phi_ij Interaction parameter for viscosities
phi_i Fugacity coefficient of component i
phi^* Pure fugacity coefficient at saturation
sigma Surface tension (N/m)
Collision diameter (Angstrom)
sigma_b Surface tension at T_b (N/m)
sigma , Liquid-liquid interfacial tension (N/m)
mu Dipole moment (Debye)
xi Inverse viscosity (defined in text)
Superscripts:
L Liquid
V, G Vapour, gas
* Saturated liquid,
T/epsilon
Subscripts:
b at normal boiling point
c critical
i of component i
j of component j
m mixture
r reduced
s saturated liquid
Abbreviations:
EOS Equation of State
RK Redlich-Kwong
SRK Soave Redlich-Kwong
PR Peng-Robinson
E.N. Fuller, K. Ensley, J.C. Giddings, ``A New Method for Prediction of Binary Gas-Phase Diffusion Coefficients'', Ind. Eng. Chem., Vol. 58, (1966) pp. 19--27.
E.N. Fuller, P.D. Schettler, J.C. Giddings, ``Diffusion of Halogonated Hydrocarbons in Helium. The Effect of Structure on Collision Cross sections'', J. Phys. Chem., Vol. 73 (1969) pp. 3679--3685.
W. Hayduk, H. Laudie, ``Prediction of Diffusion Coefficients for Nonelectrolytes in Dilute Aqueous Solutions'', AIChE J., Vol. 20, (1974) pp. 611--615.
W. Hayduk, B.S. Minhas, ``Correlations for Prediction of Molecular Diffusivities in Liquids'', Can. J. Chem. Eng., 60, 295-299 (1982); Correction, Vol. 61, (1983) pp. 132.
J.O. Hirschfelder, C.F. Curtis, R.B. Bird, Molecular Theory of Gases and Liquids, Wiley, New York (1954).
F. Jufu, L. Buqiang, W. Zihao, ``Estimation of Fluid-Fluid Interfacial Tensions of Multicomponent Mixtures'', Chem. Eng. Sci., Vol. 41, No. 10 (1986) pp. 2673--2679.
R. Taylor, H.A. Kooijman, ``Composition Derivatives of Activity Coefficient Models (For the estimation of Thermodynamic Factors in Diffusion)'', Chem. Eng. Comm., Vol. 102, (1991) pp. 87--106.
A. Letsou, L.I. Stiel, AIChE J., Vol. 19 (1973) pp. 409.
Lielmezs, Herrick, Chem. Eng. J., Vol. 32 (1986) pp. 165.
J.M. Prausnitz, T. Anderson, E. Grens, C. Eckert, R. Hsieh, J. O'Connell, Computer Calculations for Multicomponent Vapor-Liquid and Liquid-Liquid Equilibria, Prentice-Hall (1980).
J.S. Rowlinson, Liquids and Liquid Mixtures, 2nd Ed., Butterworth, London (1969).
R.C. Reid, J.M. Prausnitz, T.K. Sherwood, Properties of Gases and Liquids, 3rd Ed., McGraw-Hill, New York (1977).
R.C. Reid, J.M. Prausnitz and B.E. Poling, The Properties of Gases and Liquids, 4th Ed., McGraw-Hill, New York (1988).
M.A. Siddiqi, K. Lucas, ``Correlations for Prediction of Diffusion in Liquids'', Can. J. Chem. Eng., Vol. 64 (1986) pp. 839--843.
A.S. Teja, P. Rice, ``Generalized Corresponding States Method for the Viscosities of Liquid Mixtures'', Ind. Eng. Chem. Fundam., Vol. 20 (1981) pp. 77-81.
S.M. Walas, Phase Equilibria in Chemical Engineering, Butterworth Publishers, London (1985).
C.R. Wilke, J. Chem. Phys., Vol. 18 (1950) pp. 617.
C.R. Wilke, P. Chang, ``Correlation of Diffusion Coefficients in Dilute Solutions'', AIChE J., Vol. 1 (1955) pp. 264-270.
C.R. Wilke, C.Y. Lee, Ind. Eng. Chem., Vol. 47 (1955) pp. 1253.
P.H. Winterfeld, L.E. Scriven, H.T. Davis, AIChE J., Vol. 24 (1978) pp. 1010.
V + L - F = 0
V y + L x - Fúz = 0
i i i
K x - y = 0
i i i
c
sum (y - x ) = 0
i=1 i i
V L F
V H + L H - F H + Q = 0
T
(X) = (V, y , y ... y , T, x , x ... x , L)
1 2 c 1 2 c
the vector of functions, (F), is:
T
(F) = (TMB ,CMB , CMB ... CMB , H, EQM , EQM ... E , SUM),
1 2 c 1 2 c
The structure of the Jacobian matrix [J] is shown below:
King, C.J., Separation Processes, Second Edition, McGraw Hill (1980).
The equations that model equilibrium stages are called the MESH equations. The MESH equations for the interior stages of a column together with equations for the reboiler and condenser (if they are needed) are solved together with any specification equations to yield, for each stage, the vapour mole fractions; the liquid mole fractions; the stage temperature and the vapour and liquid flowrates.
Since the late 1950's, hardly a year has gone by without the publication of at least one (and usually more than one) new algorithm for solving the equilibrium stage model equations. One of the incentives for the continued activity has always been (and remains) a desire to solve problems with which existing methods have trouble. The evolution of algorithms for solving the MESH equations has been influenced by, among other things: the availability (or lack) of sufficient computer storage and power, the development of mathematical techniques that can be exploited, the complexity of physical property (K-value and enthalpy) correlations and the form of the model equations being solved.
It is not completely clear who first implemented a simultaneous correction method for solving multicomponent distillation and absorption problems. As is so often the case, it would appear that the problem was being tackled by a number of people independently. Simultaneous solution of all the MESH equations was suggested as a method of last resort by Friday and Smith (1964) in a classic paper analysing the reasons why other algorithms fail. They did not, however, implement such a technique. The two best known and most frequently cited papers are those of Goldstein and Stanfield (1970) and Naphtali and Sandholm (1971), the latter providing more details of an application of Newton's method described by Naphtali at an AIChE meeting in May 1965.
To the best of our knowledge, a method to solve all the MESH equations for all stages at once using Newton's method was first implemented by Whitehouse (1964) (see, also, Stainthorp and Whitehouse, 1967). Among other things, Whitehouse's code allowed for specifications of purity, T, V, L or Q on any stage. Interlinked systems of columns and nonideal solutions could be dealt with even though no examples of the latter type were solved by Whitehouse. Since the pioneering work of Whitehouse, Naphtali and Sandholm and Goldstein and Stanfield, many others have employed Newton's method or one of its relatives to solve the MESH equations.
Simultaneous correction procedures have shown themselves to be generally fast and reliable. Extensions to the basic method to include complex column configurations, interlinked columns, nonstandard specifications and applications to column design result in only minor changes in the algorithm. In addition, simultaneous correction procedures can easily incorporate stage efficiencies within the calculations (something that is not always possible with other algorithms). Developments to about 1980 have been described in a number of textbooks (see, for example, Holland, 1963, 1975, 1981; King, 1980; Henley and Seader, 1981) and a recent review by Wang and Wang (1980). Seader (1985) has written an interesting history of equilibrium stage simulation.
Seader (1986) lists a number of things to be taken into consideration when designing a simultaneous correction method; a revised and extended list follows and is discussed in more detail below.
The equations that model equilibrium stages are termed the MESH equations, MESH being an acronym referring to the different types of equations that form the mathematical model. The M equations are the Material balance equations, the E equations are the Equilibrium relations, the S equations are the Summation equations and the H equations are the entHalpy balances:
T
M = (W + V ) + (U + L ) - V - L - F = 0
j j j j j j+1 j-1 j
M = (W + V ) y + (U + L )úx - V y - L x - F z
ij j j ij j j ij j+1 i,j+1 j-1 i,j-1 j ij
E = K x - y = 0
ij ij ij ij
V c
S = sum y - 1 = 0
j i=1 ij
L c
S = sum x - 1 = 0
j i=1 ij
V L V L F
H = (W + V ) H + (U + L ) H - V H - L H - F H + Q = 0
j j j j j j j j+1 j+1 j-1 j-1 j j j
If we count the equations listed, we will find that there are 2c+4 equations per stage. However, only 2c+3 of these equations are independent. These independent equations are generally taken to be the c component mass balance equations, the c equilibrium relations, the enthalpy balance and two more equations. These two equations can be the two summation equations or the total mass balance and one of the summation equations (or an equivalent form). The 2c+3 unknown variables determined by the equations are the c vapour mole fractions, y; the c liquid mole fractions, x; the stage temperature, T and the vapour and liquid flowrates, V and L. For a column of s stages, we must solve s(2c+3) equations. The table below shows how we may easily end up having to solve hundreds or even thousands of equations.
| c | s | s(2c+3) |
| 2 | 10 | 70 |
| 3 | 20 | 180 |
| 5 | 50 | 650 |
| 10 | 30 | 690 |
| 40 | 100 | 8300 |
The first entry in this table corresponds to a simple binary problem that could easily be solved graphically. The second and third are fairly typical of the size of problem encountered in azeotropic and extractive distillation processes. The last two entries are typical of problems encountered in simulating hydrocarbon and petroleum mixture separation operations.
In the case of a total condenser, the vapour phase compositions used in the calculation of the equilibrium relations and the summation equations are those that would be in equilibrium with the liquid stream that actually exists. That is for the total condenser, the vapour composition used in the equilibrium relations is the vapour composition determined during a bubble point calculation based on the actual pressure and liquid compositions found in the condenser. At the same time, these compositions are not used in the component mass balances since there is no vapour stream from a total condenser.
MV MV
E = E K x - y - (1 - E ) y = 0
ij j ij ij ij j i,j+1
where E^MV_j is the Murphree vapour efficiency for stage j. If the
efficiency is unity we obtain the original equilibrium relation from
above. The Murphree efficiency must be specified for all stages except
the condenser and reboiler which are assumed to operate at
equilibrium.
T
(X ) = (V , y , y ... y , x , x ... x , L )
j j 1j 2j cj 1j 2j cj j
and a vector of functions for the j-th stage (F_j):
T T L-V
(F ) = (M , M , M ... M , H , E , E ... E , S )
j j 1j 2j cj j 1j 2j cj j
where
L-V c
S = sum (x - y ) = 0
i=1 ij ij
If the equations and variables are grouped by stage we have
T T T T
(F) = ((F ) , (F ) ... (F ) )
1 2 s
T T T T
(X) = ((X ) , (X ) ... (X ) )
1 2 s
| [B1] | [C1] | |||||
| [A2] | [B2] | [C2] | ||||
| [A3] | [B3] | [C3] | ||||
| . | . | . | ||||
| . | . | . | ||||
| [Am-1] | [Bm-1] | [Cm-1] | ||||
| [Am] | [Bm] |
in which each entry [A], [B], [C] is a matrix in its own right. The [A] submatrices contain partial derivatives of the equations for the j-th stage with respect to the variables for stage j-1. The [B] submatrices contain partial derivatives of the equations for the j-th stage with respect to the variables for the j-th stage. Finally, the [C] submatrices contain partial derivatives of the equations for the j-th stage with respect to the variables for stage j+1. The structure of the submatrices [A], [B] and [C] is shown below.
The elements of [B] include partial derivatives of K-values with respect to temperature and composition. Since it is rather a painful experience to differentiate, for example, the UNIQUAC equations with respect to temperature and composition, in some SC codes this differentiation is done numerically. This can be an extremely time consuming step. However, neglect of these derivatives is not recommended unless one is dealing with nearly ideal solutions, since, to do so, will almost certainly lead to an increase in the required number of iterations or even to failure.
Almost all of the partial derivatives needed in ChemSep are computed from analytical expressions. The exception is the temperature derivatives of the excess enthalpy which requires a second differentiation with respect to temperature of the activity and fugacity coefficient models. Coding just the first derivatives was bad enough.
If the column has pumparounds extra matrices will be present which are not on the diagonal and the use of block tridiagonal methods becomes less straight forward. Similar problems arise with non-standard specifications that are not on the variables of the condenser and reboiler. When we solve multiple interlinked columns (currently not supported by ChemSep) a special ordering is required to maintain the diagonal structure of the jacobian. Therefore, ChemSep now uses a sparse solver to solve the system of equations involved. In the case of non-standard specifications (which allow "wild" specification such as the temperature on a stage in the middle of a column, instead at one of its ends) the specification equation must be calculated using finite differencing.
Still further improvements in the block elimination algorithm for solution of separation process problems can be effected if we take advantage of the special structure of the submatrices [A], [B] and [C]. In fact, [A] and [C] are nearly empty. ChemSep uses the sparse matrix solver NSPIV for solving the sparse linear system effectively.
As ChemSep calculates the Jacobian mostly analytically, the Jacobian is computed for each iteration of the Newton's method. This has shown to be the most effective method until now. Alternatively, the Jacobian could be held constant for a number of iterations or approximated using quasi-Newton methods. ChemSep currently does not implement these options. Another possibility is to selectively retain the derivatives of physical properties in the calculation of the Jacobian. For dynamic simulations this has shown to be counterproductive (Kooijman, 1995) and therefore ChemSep doesn"t support such a technique.
As ChemSep does not scale the equations it is possible that this sum of squares remains too large for such small convergence criterions. As this is not always forseeable, we extended the convergence test in ChemSep with a check on the change in the variables. If the absolute relative value of the change of each variable is smaller than the convergence criterion, ChemSep also considers the problem solved.
Finally, in order to prevent calculations to continue endlessly, there is a maximum number of iterations after which the Newton's method exits. ChemSep will then ask the user whether to continue the calculations, and for how many additional iterations. Typically, we found that 30 Newton iterations suffices to solve most separation problems at hand. Often, when more iterations are needed, there are other problems which prevent the method to converge. For example, specifications can be impossible to attain, or the specifications define a problem with multiple solutions, where the Newton's method oscillates between two or more solutions (such problems have only been found and documented recently).
An exception to this observation of convergence form the simulations of packed columns using the nonequilibrium with plug flow model approximations, where convergence tends to be slow.
The next step is to estimate the compositions and temperatures. This is done iteratively. We start by estimating the K-values assuming ideal solution behavior at the column average pressure and an estimate of the boiling point of the combined feeds. (This eliminates the normal requirement of estimates of end stage temperatures).
Mole fractions of both phases are estimated by solving the material balance and equilibrium equations for one component at a time.
We define a column matrix of discrepancy functions
T c c c
(F) = (ME , ... ME , ... ME )
i,1 i,j i,s
where ME^c_ij is the component material balance equation combined with
the equilibrium equations to eliminate the vapour phase mole
fractions. Each equation depends on only three mole fractions. Thus,
if we define a column matrix of mole fractions (X) by
T
(X) = (x , ... x , x , x , x )
i,1 i,j-1 ij i,j+1 i,s
we may write
(F) = [ABC](X) - (R) = (0)With the equations and variables ordered in this way, the coefficient matrix [ABC] has three adjacent diagonals with coefficients:
A = L j j-1 B = - (V K j + L ) j j i j C = V K j j j iThe right hand side matrix (R) has elements
R = - F z j j ijThis linear system of equations can be solved for the mole fractions very easily using Gaussian elimination. Temperatures and vapour compositions are computed from a bubble point calculation for each stage. The bubble point computation provides K-values for all components on all stages. So we solve the tridiagonal system of equations again using the old flow rates and the new K-values. Temperatures are recomputed as before. This procedure is repeated a third time before proceeding with the main simulation.
For columns without condenser and reboiler a different initialization is used where the compositions of the liquid are set equal to the top liquid feed compositions and the compositions of the vapour equal to the bottom vapour feed compositions. Temperatures in the whole column are set equal to the temperature of the first feed specified.
For columns with either a condenser or a reboiler the compositions are initialized to the overall compositions of all feeds combined, and the temperatures to the bubble point at the column pressure and the overall feed compositions.
Currently there is no special initialization routine that will handle pumparounds if they are present. The result is that columns with large pumparoundflows will require flow initialization by the user, or, to repetitively solve the problem using the old results as the initialization and increasing pumparoundflow (see below).
Several methods have been used to improve the reliability of Newton's method; damping the Newton step, use of steepest descent (ascent) formulations for some of the iterations, and combination with relaxation procedures; none of which has proven to be completely satisfactory. The methods most recently proposed for assisting convergence of Newton's method are continuation methods.
In default mode, ChemSep does NOT use any of these techniques, other than a check to make sure that all quantities remain positive. Mole fractions, for example, are not permitted to take on negative values. The user does have the option of supplying damping factors.
If for some reason your column simulation does not converge, changing the damping factors might help. If you know the iteration history (by either limiting the number of iterations or by printing out the intermediate answers. Both can be done in the Options - Solve options menu) you can adapt the factors so the column simula