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 owhere 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 oThis 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) owhere 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 iwhere 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 iThis 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 iThis 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 iYou 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 iwhere 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 jIt 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 jIt 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 kThe 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 kThe 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 iThe 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 iNote 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 iThe 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 iDIPPR 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 rwhere 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,ijBinary 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 ciand 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 band 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 iwhere 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 band 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 riand 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 band 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 iYou 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 iwhich 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 iwhere 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 iIf 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 mIf 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 rwhere
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 SRKwhere
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,mThis 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 rHowever, 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 iwhere 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 cjwhere 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 1with epsilon defined as
2/3 epsilon = (V / sqrt(T M )) i ci ci iand 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 iAlternatively 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 SLwhere 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 ijwhere 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 jor 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 jIf 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 jotherwise 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,iWhere 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 hpwhere 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 ) vwhere 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 rwhere 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 cwhere 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 cIf 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 celse by
o 1.72 F = 1 + 30.55 (0.292 - Z ) [0.96 + 0.1 (T - 0.7)] p c rThe 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 rwhere 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 qwhere 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 rwhere, 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 qiwhere 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 Lelse 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 iwhere 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 rThis 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 ifor 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 cwhere 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 ijwhere 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,iNote that the component viscosities are required for this evaluation.
rho = (1 / V ) c c,m rho = (rho / rho ) r cIf 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,mwhich 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 ijLiquid 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,awhere 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 bThis 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 bwhere 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 bThe diffsuion collision integral is
* T = T / epsilon ab * -b * * * Omega = a (T ) + c / exp (d T ) + e / exp (f T ) + g / exp (h T ) Dwhere 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 band the diffusion coefficient is
V 3/2 D = C (T sqrt(1/M +1/M ) / P sigma w ) ab a b ab Dwhere 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,iThen 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 2This 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 3with 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 cthe vector of functions, (F), is:
T (F) = (TMB ,CMB , CMB ... CMB , H, EQM , EQM ... E , SUM), 1 2 c 1 2 cThe 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+1where 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 jand 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 jwhere
L-V c S = sum (x - y ) = 0 i=1 ij ijIf 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,swhere 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,swe 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 simulation might converge. Note that convergence is mostly slower when you start to apply extra damping by making the factors smaller, the Newton method looses its effectiveness when damped. Nor does damping guarantee convergence.
The user can specify either temperature or flow profiles, or both. The only requirement is that values for the first and last stages are provided. Missing temperatures on intermediate stages are computed by linear interpolation, missing flows are computed on a constant flow from stage to stage basis. Therefore, it is better to specify the flows of the first and last two stages in case a condenser and reboiler are present. Composition profiles are computed through the method described above, however, temperatures are not computed using the bubble point calculations. If both user specified temperature and flow profiles are incomplete ChemSep switches to the automatic initialization method.
R.P. Goldstein, R.B. Stanfield, "Flexible Method for the Solution of Distillation Design Problems using the Newton-Raphson Technique", Ind. Eng. Chem. Process Des. Dev., 9, 78 (1970).
E.J. Henley, J.D. Seader, Equilibrium-Stage Separation Operations in Chemical Engineering, Wiley (1981).
C.D. Holland, Multicomponent Distillation, Prentice Hall Inc., NJ (1963).
C.D. Holland, Fundamentals and Modelling of Separation Processes, Prentice Hall Inc., NJ (1975).
C.D. Holland, Fundamentals of Multicomponent Distillation, McGraw-Hill Inc.; New York (1981).
C.J. King, Separation Processes, Second Edition, McGraw Hill (1980).
L.M. Naphtali, "The Distillation Column as a Large System", presented at AIChE 56-th National Mtg., San Francisco, May 16, (1965).
L.M. Naphtali, D.P. Sandholm, "Multicomponent Separation Calculations by Linearization", AIChE J., Vol 17 (1), 148 (1971).
J.D. Seader, "The BC (Before Computers) and AD of Equilibrium Stage Operations", Chem. Eng. Ed., Spring, 88, (1985a).
J.D. Seader, Chapter on Distillation in Chemical Engineers Handbook, (Green D. Editor), 6th Edition, McGraw Hill, New York, (1986)
J.D. Seader, Computer Modelling of Chemical Processes, AIChE Monograph Series, No. 15, 81 (1986).
B.D. Smith, Design of Equilibrium Stage Processes, McGraw-Hill, New York, (1964).
F.P. Stainthorp, P.A. Whitehouse, "General Computer Programs for Multi Stage Counter Current Separation Problems - I: Formulation of the Problem and Method of Solution", I. Chem. E. Symp. Ser., Vol 23, 181 (1967).
J.C. Wang, Y.L. Wang; "A Review on the Modeling and Simulation of Multi-Stage Separation Processes" in Foundations of Computer-Aided Chemical Process Design, vol. II; R.S.H. Mah and W.D. Seider, eds.; Engineering Foundation; 121 (1981).
S.M Walas, Phase Equilibria in Chemical Engineering, Butterworth Publishers, (1985).
P.A. Whitehouse, A General Computer Program Solution of Multicomponent Distillation Problems, Ph.D. Thesis in Chem.Eng., University of Manchester, Institute of Science and Technology, Manchester, England (1964).
Figure 5.1. Schematic diagram of a nonequilibrium stage (after Taylor and Krishna, 1993). Figure 5.1 also serves to introduce the notation used in writing down the equations that model the behavior of this nonequilibrium stage. The flow rates of vapor and liquid phases leaving the j-th stage are denoted by V_j and L_j respectively. The mole fractions in these streams are y_ij and x_ij. The N_ij are the molar fluxes of species i on stage j. When multiplied by the area available for interphase mass transfer we obtain the rates of interphase mass transfer. The temperatures of the vapor and liquid phases are not assumed to be equal and we must allow for heat transfer as well as mass transfer across the interface.
If Figure 5.1 represents a single tray then the term phi_j^L is the fractional liquid entrainment defined as the ratio of the moles of liquid entrained in the vapor phase in stage j to the moles of downflowing liquid from stage j. Similarly, phi_j^V is the ratio of vapor entrained in the liquid leaving stage j (carried down to the tray below under the downcomer) to the interstage vapor flow. For packed columns, this term represents axial dispersion. Weeping in tray columns may be accounted for with a similar term.
The component material balance equations for each phase may be written as follows:
V V V V V n V M == ( 1 + r + phi ) V y - V y - phi V y - f - sum G + N ij j j j ij j+1 i,j+1 j-1 j-1 i,j-1 ij nu =1 ij nu ij = 0 i = 1,2,...,c
L L L L L n L M == ( 1 + r + phi ) L x - L x - phi L x - f - sum G - N ij j j j ij j-1 i,j-1 j+1 j+1 i,j+1 ij nu =1 ij nu ij = 0 i = 1,2,...,cwhere G_ij nu is the interlinked flow rate for component i from stage nu to stage j, and n is the number of total stages (trays or sections of packing). The last terms in Equations () and () are the mass transfer rates (in kmol/s), where mass transfer from the ``V" phase to the ``L" phase is defined as positive. At the V/L interface we have continuity of mass and, thus, the mass transfer rates in both phases must be equal.
The total material balances for the two phases are obtained by summing Equations () and () over the component index i.
V V V V V c n V M == ( 1 + r + phi ) V - V - phi V - F - sum sum G + N t_{j} j j j j+1 j-1 j-1 j i=1 nu=1 ijnu tj = 0
L L L L L c n L M == ( 1+r + phi )L - L - phi L - F - sum sum G - N t_{j} j j j j-1 j+1 j+1 j i=1 nu=1 ijnu tj = 0F_j denotes the total feed flow rate for stage j, F_j=sum _i=1^c f_ij. Here total flow rates and mole fractions are used as independent variables and total as well as component material balances are included in the set of independent model equations. In the nonequilibrium model of Krishnamurthy and Taylor (1985a) component flow rates were treated as variables.
The nonequilibrium model uses two sets of rate equations for each stage:
V V R == N - N = 0 i = 1,2 ...,c-1 ij ij ij L L R == N - N = 0 i = 1,2 ...,c-1 ij ij ijwhere N_ij is the mass transfer rate of component i on stage j. The mass transfer rate in each phase is computed from a diffusive and a convective contribution with
V I V N = a J + y N ij j ij ij tj L I L N = a J + x N ij j ij ij tjwhere a^I_j is the total interfacial area for stage j and N_tj is the total rate on stage j (N_tj=sum ^c_i=1 N_ij). The diffusion fluxes J are given by (in matrix form):
V V V V I (J ) = c [k ]((y - y )) t L L L I L (J ) = c [k ]((x - x )) twhere (y^V - y^I) and (x^I - x^L) are the average mole fraction difference between the bulk and the interface mole fractions (Note that the fluxes are multiplied by the interfacial area to obtain mass transfer rates). How the average mole fraction differences are calculated depends on the selected flow model. The matrices of mass transfer coefficients, [k], are calculated from
P P -1 P [k ] = [R ] [Gamma ]where [Gamma^P] is a matrix of thermodynamic factors for phase P. For systems where an activity coefficient model is used for the phase equilibrium properties the thermodynamic factor matrix Gamma (order c-1) is defined by
Gamma = delta + x ( (d ln gamma / d x ) ) ij ij i i j T,P,x_k,k != j=1 ... c-1If an equation of state is used gamma_i is replaced by phi_i. Expressions for the composition derivatives of ln gamma_i are given by Taylor and Kooijman (1991). The rate matrix [R] (orderc-1) is a matrix of mass transfer resistances calculated from the following formulae:
P P c P R = (z / k ) + sum (z / k ) ii i ic k=1,k != i k ik P P P R = -z ( (1 / k ) - (1 / k ) ) ij i ij icwhere k_ij^P are binary pair mass transfer coefficients for phase P. Mass transfer coefficients, k_ij, are computed from empirical models (Taylor and Krishna, 1993) and multicomponent diffusion coefficients evaluated from an interpolation formula (Kooijman and Taylor, 1991). Equations () and () are suggested by the Maxwell-Stefan equations that describe mass transfer in multicomponent systems (see Taylor and Krishna, 1993). The matrix of thermodynamic factors appears because the fundamental driving force for mass transfer is the chemical potential gradient and not the mole fraction or concentration gradient. This matrix is calculated from an appropriate thermodynamic model. The binary mass transfer coefficients are estimated from empirical correlations as functions of column internal type as well as design, operational parameters, and physical properties including the binary pair Maxwell-Stefan diffusion coefficients. Thus, the mass transfer coefficient models form the basis of the nonequilibrium model and it is possible to change the behavior of a column by selecting a different mass transfer coefficient correlation.
Note that there are c times c binary pair Maxwell-Stefan diffusion coefficients, but only c-1 times c-1 elements in the [R^P] and [k^P] matrices and, therefore, only c-1 equations per set of rate equations. This is the result of the fact that diffusion calculations only yield relative transfer rates. We will need an extra equation that will "bootstrap" the mass transfer rates: the energy balance for the interface. Note also that in this model the flux correction on the mass transfer coefficients has been neglected.
The energy balance equations on stage j are written for each phase as follows:
V V V V V V V V VF n V V E == ( 1 + r + phi ) V H - V H - phi V H - F H - sum (G) H j j j j j j+1 j+1 j-1 j-1 j-1 j j nu =1 jnu jnu V V + Q + e = 0 j j
L L L L L L L L LF n L L E == ( 1 + r + phi ) L H - L H - phi L H - F H - sum (G) H j j j j j j-1 j-1 j+1 j+1 j+1 j j nu=1 jnu j nu L L + Q - e = 0 j jwhere G_jnu is the interlink flow rate from stage nu to stage j. The last term in the left-hand-side of Equations () and (), e_j, represents the energy transfer rates for the vapor and liquid phase which are defined by
V I V V I c V V e = a h (T - T ) + sum N ( H) j j i=1 ij ij L I L I L c L L e = a h (T - T ) + sum N ( H) j j i=1 ij ijwhere H_ij are the partial molar enthalpies of component i for stage j. We also have continuity of the energy fluxes across the V/L interface which gives the interface energy balance:
I V L E == e - e = 0 j j jwhere h^V and h^L are the vapor and liquid heat transfer coefficients respectively, and T^V, T^I, and T^L the vapour, interface, and liquid temperatures. For the calculation of the vapour heat transfer coefficients the Chilton-Colburn analogy between mass and heat transfer is used:
Le = (lambda / D C rho) = (Sc / Pr) p V 2/3 h = k rho C Le pFor the calculation of the liquid heat transfer coefficients a penetration model is used:
L h = k rho C sqrt(Le) pwhere k is the average mass transfer coefficient and D the average diffusion coefficient.
In the nonequilibrium model of Krishnamurthy and Taylor (1985a) the pressure was taken to be specified on all stages. However, column pressure drop is a function of tray (or packing) type and design and column operating conditions, information that is required for or available during the solution of the nonequilibrium model equations. It was, therefore, quite straightforward to add an hydraulic equation to the set of independent equations for each stage and to make the pressure of each stage (tray or packed section) an unknown variable. The stage is assumed to be at mechanical equilibrium so p_j^V = p_j^L = p_j.
In the second generation model, the pressure of the top tray (or top of the packing) is specified along with the pressure of any condenser. The pressure of trays (or packed sections) below the topmost are calculated from the pressure of the stage above and the pressure drop on that tray (or over that packed section). If the column has a condenser (which is numbered as stage 1 here) the hydraulic equations are expressed as follows:
P == p - p = 0 1 c 1
P == p - p = 0 2 spec 2
P == p - p - (Delta p ) = 0 j = 3, 4, ..., n j j j-1 j-1where p_c is the specified condenser pressure, p_spec is the specified pressure of the tray or section of packing at the top of the column, and Delta p_j-1 is the pressure drop per tray or section of packing from section/stage j-1 to section/stage j. If the top stage is not a condenser, the hydraulic equations are expressed as
P == p - p = 0 1 spec 1
P == p - p - (Delta p ) = 0 j = 2, 3, ..., n j j j-1 j-1In general we may consider the pressure drop to be a function of the internal flows, the fluid densities, and equipment design parameters.
V L Delta p = f(V , L , rho , rho , Design) j-1 j-1 j-1 j-1 j-1The pressure drop term, Delta p_j-1, is calculated from liquid heights on the tray (from various correlations, see Lockett, 1986, and Kister, 1992) or specific pressure drop correlations for packings (see the section below on pressur drop models).
For bubble cap trays the procedures described by Bolles (1963) can be adapted for computer based calculation. Kister (1992) also covers methods available for estimating the pressure drop in dumped packed columns. The pressure drop in structured packed columns may estimated using the method of Bravo et al. (1986).
Phase equilibrium is assumed to exist only at the interface with the mole fractions in both phases related by:
I I I Q == K x - y = 0 i = 1,2, ...,c ij ij ij ijwhere K_ij is the equilibrium ratio for component i on stage j. The K_ij are evaluated at the (calculated) temperature, pressure, and mole fractions at the interface.
The mole fractions must sum to unity in each phase:
V c S == sum y - 1 = 0 j i=1 ij
L c S == sum x - 1 = 0 j i=1 ijas well as at the interface:
VI c I S == sum y - 1 = 0 j i=1 ij
LI c I S == sum x - 1 = 0 j i=1 ijTable 5.1 lists the type and number of equations for the nonequilibrium model. To solve the model we have 5c+6 equations and variables, where c is the number of components. They are solved simultaneously using Newton's method.
Table 5.1. Nonequilibrium model equations type and number
Equation | Number |
Material balances | 2c+2 |
Energy balances | 3 |
transfer Rate equations | 2c-2 |
Summations equations | 2 |
Hydraulic equation | 1 |
interface e Quilibrium relations | c |
Total MERSHQ | 5c+6 |
If we solve the nonequilibrium model with Newton's method, we also require initial guesses for all the variables. ChemSep uses the same automatic initial guess routine for the nonequilibrium model as the equilibrium model. The temperatures of the vapour, interface, and liquid are initialized all equal to the temperature from the automatic guess. Mass and energy transfer rates are initialized as zero and the interface mole fractions are set equal to the bulk mole fractions which are also provided by the initial guess. Pressure drops are initially assumed to be zero.
The nonequilibrium model, in comparison with the equilibrium model, requires the evaluation of many more physical properties and of the heat and mass transfer coefficients. In addition, a nonequilibrium simulation cannot proceed without some knowledge of the column type and the internals layout. Tray type and mechanical layout data, for example, is needed in order to calculate the mass transfer coefficients for each tray. For packed columns the packing type, size and material must be known. Libraries with standard tray and packing data are available on-line. Table 5.2 lists the currently supported types of column internals.
Table 5.2. Currently supported column internals for the nonequilibrium model
Bubble-cap trays |
Sieve trays |
Valve trays (including double weight valves) |
Dumped packings |
Structured packings |
Equilibrium stage (with Murphree stage efficiency) |
Rotating Disk Contactor (RDC) compartment (for extraction) |
For the evaluation of the heat and mass transfer coefficients, pressure drop, and the entrainment/weeping flows a nonequilibrium simulation needs the following:
Table 5.3. Available mass transfer coefficient correlations per internals type
Bubble-Cap | Sieve | Valve | Dumped | Structured |
tray | tray | tray | packing | packing |
AIChE | AIChE | AIChE | Onda 68 | Bravo 85 |
Hughmark | Chan-Fair | Bravo 82 | Bravo 92 | |
Zuiderweg | Billet 92 | Billet 92 | ||
Harris | Nawrocki 91 | |||
Bubble-Jet | ||||
V V V k = (N) / t a V L L L k = (N) / t a Lwhere the vapor and liquid areas are calculated with
V a = a / epsilon h d f L a = a / alpha h d fthe interfacial area density is computed according Zuiderweg (1982).
V (N) = (0.776 + 4.57 h - 0.238 F + 104.8 Q / W ) / sqrt(Sc ) w s L l V V F = u sqrt(\rho ) s s t V V V Sc = mu / rho D V t L L (N) = 19700 sqrt(D ) (0.4 F + 0.17) t s L t = h Z W / Q L L l LThe clear liquid height h_L is computed from Bennett et al. (1983):
0.67 h = alpha ( h + C ( Q / alpha W ) ) L e w L e l V L V 0.5 0.91 alpha = exp ( -12.55 (u (rho /(rho -rho )) ) ) e s C = 0.50 + 0.438 exp (-137.8 h ) w
V V (N) = ( 10300 - 8670 F ) F sqrt(D ) t / sqrt(h ) f f V L t = ( 1 - alpha ) h / ( alpha u ) V e L e sFor the liquid number of transfer units the same correlations as given for the AIChE method is used (h_L and alpha_e are also computed with the correlation of Bennett et al.).
V V V 2 k = 0.13 / rho - 0.065 / (rho ) t tin which k^V becomes independent of the diffusion coefficient. The liquid mass transfer coefficient is computed from either:
L -5 L -0.25 k = 2.6 10 (mu )or
L L 0.25 k = 0.024 (D )The interfacial area is computed in the spray regime from:
0.3 2 V 0.37 a h = (40 / phi ) (( U rho h FP / sigma )) d f s t Lor in the froth-emulsion regime:
0.3 2 V 0.37 a h = (43 / phi ) (( U rho h FP / sigma )) d f s t LThe transition from the spray to mixed froth-emulsion flow is described by:
FP > 3 b h Lwhere b is the weir length per unit bubbling area:
b = W / A l band the clear liquid height is given by:
0.5 0.25 h = 0.6 h (p FP / b) L w
V (N) = (0.051 + 0.0105 F ) sqrt(\rho \over F ) s L s L 4 (N) = ( -44 +10.7747 10 Q / W + 127.1457 F ) sqrt(D ) A / Q L l s L bub L
V (N) = (0.3 + 15 t / sqrt(Sc )) G G L (N) = (5 + 10 t ( 1 + 0.17 (0.82 F - 1) (39.3 h + 2) ) / sqrt(Sc )) L s w L
t = (h / u ) V l s F = U sqrt(\rho ) s s V V 0.1 0.14 2 2 1/3 (N) = 11 (1 / eta beta ) (( rho F / sigma )) sqrt(D t ) L L s V Vand for the liquid
t = (rho / rho ) t L L V V L 0.1 0.14 2 2 1/3 (N) = 14 (1 / eta beta ) (( rho F / sigma )) (( V / L )) sqrt(D t ) L L s L L
0.9857 d = 1.36 d s h 2 u = sqrt(((2.14 \sigma + 0.505 \rho g d \over \rho d )) s l s l s 0.8464 0.21 d = 0.8868 d u b h h u = (u / (1-f )(1-h /h )) - (u f / (1-f )) b v s cl f s s s 1.32 1.33 f = 165.65 d beta s h f = 1 - f b s h = h - h b f j beta = (A / A ) h bub -6 h = 2.853 10 Re j h Re = (d u rho / eta ) h h h g g d = 1.1 d + 0.25 h j h cl 2 2 u = (u d / (1-FLC) d ) j h h jwhere parameters as h_f and h_cl can be computed by empirical correlations for the tray (here sieve tray since the correlations were obtained with that particular tray type). The fraction of inactive holes, FLC, can be set to zero or estimated by
-1.602 0.524 0.292 FLC = 1836.97 u Q h h L wCurrently this model is not available in ChemSep.
V 0.7 0.333 V -2 k = A Re Sc (a D ) (a d ) V V p p pwhere A=2 if d_p < 0.012 and A=5.23 otherwise. Vapour and liquid velocities are calculated by
V V u = V M / rho A V t L L u = L M / rho A L tand Reynolds and Schmidt numbers:
V V Re = (rho u / (eta a )) V V p L L Re = (rho u / (eta a )) L L p V V V Sc = (eta / (rho D )) V t L L L Sc = (eta / (rho D )) L tThe liquid phase mass transfer coefficient is
L I 2/3 -0.5 0.4 L L 1/3 k = 0.0051 (Re ) Sc (a d ) (eta g / rho ) L L p pwhere Re^I_L is the liquid Reynolds number based on the interfacial area density
I L L Re = (rho u / (eta a )) L L dThe interfacial area density, a_d (m^2/m^3), is computed from
0.75 0.1 -0.05 0.2 a = a [ 1 - exp ( -1.45(sigma /sigma) Re Fr We ) ] d p c L L Lwhere
2 Fr = (a u / g) L p L 2 We = (rho u / a sigma) L L L p
0.392 -0.4 a = 19.78 (Ca Re ) sqrt(\sigma) H a d L V pwhere H is the height of the packed section and Ca_l is the capillary number
L L Ca = u eta / rho sigma L LSince the interfacial area density is used in the calculation for the liquid Reynolds number the Bravo and Fair method will result in different mass transfer coefficients.
L 1/6 L 1/3 k = C ((g rho / eta )) sqrt(D \over d ) ((u / a )) l l l h L p V V 3/4 1/3 k = C ((1 / sqrt(\epsilon - h ))) sqrt(a \over d ) D (Re ) (Sc ) v t h V Vwhere Reynolds and Schmidt numbers are calculated as in Onda et al.. The hydraulic diameter d_h is
d = 4 epsilon / a h pand the liquid holdup fraction, h_t, is calculated as described below under the pressure drop section. The interfacial area density is given by:
L L -0.2 2 L 0.75 2 -0.45 a = a (1.5 / sqrt(a d )) (u d rho / eta ) (u rho d / sigma) (u / g d ) d p p h L h L h L h
0.8 0.333 Sh = 0.0338 Re Sc V V Vand is defined by
V V Sh = (k d / D ) V eqThe equivalent diameter of a channel is given by
d = B h [ 1 / (B+2S) + 1/ 2 S ] eq cwhere B is the base of the triangle (channel base), S is the corrugation spacing (channel side), and h_c is the height of the triangle (crimp height). The vapour phase Reynolds number is defined by
V V Re = (d rho (u +u ) / eta ) V eq V,eff L,effThe effective velocity of vapour through the channel, u_Ve, is
u = u / (epsilon sin theta) V,eff V(u_V is the superficial vapour velocity, epsilon the void fraction, and theta the angle of the channel with respect to the horizontal). The effective velocity of the liquid is
L L 2 L 1/3 u = (3 Gamma / 2 rho ) (( (rho ) g / 3 eta Gamma )) L,effwhere Gamma is the liquid flow rate per unit of perimeter
L Gamma = rho u / P Lwhere P is the perimeter per unit cross-sectional area, computed from
P = (4S+B) / B h cThe penetration model is used to predict the liquid phase mass transfer coefficients with the exposure time assumed to be the time required for the liquid to flow between corrugations (a distance equal to the channel side):
t = S / u L L,eff L L k = 2 sqrt(D \over \pi t ) L
K = 0.614 + 71.35 S 2The mass transfer calculations are dependent on the pressure drop/holdup calculations. The effective area can be adjusted with the surface enhancement factor Fse, and the liquid resistance with a correction on the surface renewal following the penetration model (parameter C_e). Effective velocities are computed with
u = u / epsilon h sin theta L,eff L t u = u / epsilon (1-h ) sin theta G,eff V twhere h_t is the fractional liquid holdup (see below at the section on pressure drop calculation). Reynolds numbers and liquid mass transfer coefficient is now calculated as in Bravo et al. (1985) but with
t = C S / u L e L,effHowever, the vapour phase mass transfer coefficient is obtained from
V V 0.8 1/3 k = 0.054 ((D / S)) Re Sc V Vwhere the equivalent diameter is replaced with the channel side S and a different coefficient is used. The assumption of a completely wetted packing is dropped, the interfacial area density is given by
a = F F a d t se p 0.15 0 0.2 0.6 0.3 F = (29.12 (We Fr ) S .359 / Re epsilon (sin theta) (1 - 0.93 cos gamma)) t L L Lwhere cos gamma is equal to 0.9 for sigma<0.0453, otherwise it is computed by
-16.835 sigma cos gamma = 5.211 10Note that a different switch point is used than reported by Bravo et al. (1992) to guarantee continuity in cos gamma.
V I V I ((y - y )) = (y - y ) I L I L ((x - x )) = (x - x )this keeps the rate equations (relatively) simple and only a function of the mole fractions leaving the current stage. However, on a tray where the vapour bubbles through a liquid which flows from one downcomer to the opposite downcomer this model is not accurate. Indeed, only for very small diameter columns will the mixed flow model give reasonable results. The mixed model is the most simple flow model and is the easiest to converge. For packed columns the convergence to the true column profiles by using increasing number of stages can be quite slow using the mixed flow model.
V I V V I ((y - y )) = Omega[-(N) ] (y - y ) I L L I L ((x - x )) = Omega[-(N) ] (x - x )where the mole fractions are of the leaving streams and the number of mass transfer units (N) for the vapour and liquid are defined as:
V V V V (N) = c k a h A / V t f b L L L L (N) = c k a h A / L t f band Omega[M] is a matrix function defined as
-1 -1 -1 Omega[M] = [exp [M] - [I]] [M] [exp [M]] = [exp [-M] - [I]] [-M]Using this model predicted efficiencies for tray column experiments can be more accurately described. The plug flow model can also be used for packed columns, providing a much faster convergence to the true column profiles compared to the mixed flow model.
Currently no correction terms is applied to the plug flow model to correct for the change in mole fractions over the integration (as is discussed by Kooijman and taylor, 1994).
I L -1 -1 -1 -1 -1 ((x - x )) = [ [p] [exp [m] - [I]] [m] [exp [m]] - [m] [exp [p] - [I]] [p] [exp [p]] ] [b] ((X ) / 2) outwhere we have defined
L a = Pe/2 = (L Z / D W h c ) e cl L 1/2 [ b ] = a [ 2 [(N) ] / a + [I] ] [ p ] = a [I] + [b] [ m ] = a [I] - [b]Currently only a binary implementation is working for the dispersion model. Eddy dispersion coefficient are computed from Zuiderweg's (1982) correlation (this model is recommended by Korchinsky, 1994).
Results of the dispersion flow model are close to the plug flow model. How close depends on the eddy dispersion coefficient. Expect that the dispersion coefficient is larger for smaller diameter columns and trays with small weirs (or low liquid heights). We intend to extend the number of correlations predicting this coefficient. For now, it is advised to not use this flow model and it is not available.
Table 5.4. Pressure drop correlations per internals type
Bubble-Cap | Sieve | Valve | Dumped | Structured |
tray | tray | tray | packing | packing |
Fixed | Fixed | Fixed | Fixed | Fixed |
Estimated | Estimated | Estimated | Ludwig 79 | Billet 92 |
Leva 92 | Bravo 86 | |||
Billet 92 | Stichlmair 89 | |||
Stichlmair 89 | Bravo 92 | |||
h = h + h wt d lwhere h_d is the dry tray pressure drop liquid height and h_l the liquid height:
h = h + h + (h / 2) l cl r lgThe clear liquid height, h_cl, is calculated with
h = alpha h + h cl w owwhere the liquid fraction alpha is computed with the Barker and Self (1962) correlation:
alpha = (0.37 h + 0.012 F + 1.78 Q / W + 0.024 / 1.06 h + 0.035 F + 4.82 Q / W + 0.035) w s L l w s L lThe choice of correlation for the liquid fraction turns out to be important as certain correlations are dynamicly unstable. The height of liquid over the weir, h_ow, is computed by various correlations for different types of weirs (see Perry) and a weir factor (F_w) correction (see Smith, pp. 487) is employed. For example for a segmental weir:
2/3 h = 0.664 F ((Q / W )) ow w L l w = (W / D ) l c 3 2 2.5 2/3 2 2 F = w / 1 - (F w ((1.68 Q / W )) + sqrt(1-w )) w w L lwhere Q_L is the volumetric flow over the weir per weir length. The residual height, h_r, is only taken into account for sieve trays. Bennett's method (see Lockett, pp. 81) is:
2/3 1/3 h = ((6 / 1.27 rho )) ((sigma / g)) ((rho - rho / d )) r L L V hDry tray pressure, h_d, is calculated with:
2 h = K (rho / rho ) u d G L h K = (xi / 2 g)where the orrifice coefficient xi for sieve trays is computed according to Stichlmair and Mersmann (1978). For valve trays we use the method of Klein (1982) as described in Kister (1992, pp. 309--312) where K is given for the cases with the valves closed or open. It is extended for double weight valve trays as discussed by Lockett (1986, pp. 82--86). The dry tray pressure drop is corrected for liquid fractional entrainment.
The froth density is computed with
h = (h / alpha) f clThe liquid gradient, h_lg, is computed according to Fair (Lockett, 1986, pp. 72):
R = (W h / W + 2 h ) h f f U = (Q / W h ) f L cl Re = (R U rho / eta ) f h f L L 4 -1.06 f = 7 10 h Re w f 2 h = (Z f U / g R ) lg f hwhere W is the average flow-path width for liquid flow, and Z the flow path length. The height of liquid at the tray inlet is:
2 2 h = sqrt( (2 \over g) \left((Q \over W )\right) \left((1 \over h ) - (1 \over h )\right) + (2 \alpha h \over 3) ) i L l cl c fwhere h_c is the height of the clearance under the downcomer. The pressure loss under downcomer (expressed as a liquid height) is
2 h = ((1 / 2 g)) ((Q / C W h )) udc L d l cwhere C_d=0.56 according to Koch design rules. The height of liquid in the downcomer can now be calculated with the summation:
h = h + h + h db wt i udcBubble-cap liquid heights are done according to Perry's (1984) and Smith (1963). Additionally the liquid fraction of the froth is computed according to Kastanek (1970).
L L = L M / A a t V V = V M / A a t L u = L / rho L a V u = V / rho V a
2 V B (0.06243 L_a) (Delta p / Delta z) = 3.281 242 A ((0.2048 V ) / (0.06243 rho )) 10 awhere 3.281 242, 0.2048, and 0.06243 are conversion factors so that we can use A and B parameters from Wankat. Its accuracy is limited since the influence of physical properties as viscosity or surface tension on A and B are not included. Even more, the fitted parameters can be flow regime dependent. The loading regime is not well described with the simple exponent term.
L 0.2 2 0.035 L_a phi V (Delta p / Delta z) = 22.3 F (eta ) phi V (10 / g rho ) p awith phi=rho_water/rho_l=1000/rho_l. This result is similar to the Ludwig (1979) equation with corrections for the influence of the liquid density and viscosity. The only parameter is the packing factor F_p which can be obtained from dry pressure drop experiments (see Leva, 1992) or computed by the specific packing area over the void fraction cubed. Again, the loading regime is not well described with the simple exponent term. This is model is the default pressure drop model for random packings if no model is specified, since it requires only the packing factor.
V 2 4.65 (Delta p / Delta z) = 0.75 f (1 - epsilon ) rho * U / (d epsilon ) 0 p V p pwhere the void fraction of the irrigated bed, equivalent packing diameter, Reynolds number, and friction factor for a single particle are:
epsilon = epsilon - h p t d = 6 (1 - epsilon ) a p p p V V Re = u d rho / eta V V p f = C / Re + C / sqrt(Re ) + C 0 1 V 2 V 3Iteration is started by assuming a dry bed for which epsilon_p=epsilon and the holdup fraction is computed with the liquid Froude number:
2 4.65 Fr = u a / g epsilon L L p 1/3 h = 0.555 Fr t LThe liquid holdup is limited to 0.5 in order to handle flooding.
Packing dimension, hydraulic diameter and F-factor are
d = 6 (1 - epsilon) / a p p d = 4 epsilon / a h p V F = u sqrt(\rho ) s VLiquid Reynolds and Froude number are
L Re = u rho / eta a L L L p 2 Fr = u a / g L L pIf Re_L<5 then
0.15 0.1 aha = C Re Fr h L Lelse
0.25 0.1 aha = 0.85 C Re Fr h L L
L 2 1/3 hl1 = ((12 eta a u / rho g)) p L L 2/3 hl2 = hl1 aha L 0.05 h = 0.3741 epsilon ((eta rho / eta rho )) l,fl w w L
L V L V 0.2 epsilon = ((u / u )) sqrt(\rho \over \rho ) ((eta / eta )) fl L V 2 -0.39 epsilon = g / (C epsilon ) fl fl fl 1.5 L V u = sqrt(2 g / \epsilon ) (epsilon-h ) sqrt(h / a ) sqrt(\rho / \rho ) / sqrt(\epsilon) v,fl fl l,fl l,fl pif u_V>u_V,fl then h_t=h_l,fl else
13 h = hl2 + (h - hl2) (u / u ) t l,fl V V,flThe pressure drop is then
K = 1 + (2/3) (1 / 1-epsilon) (d / D ) 1 p c V V Re = u d rho / (1 - epsilon) eta K V V p 1 0.08 0.3 phi = C (64 / Re + 1.8 / Re ) exp (Re /200) (h / hl1) l1 p V V L t 3 2 (Delta p / Delta z) = phi (a / (epsilon - h ) ) (F / 2) K l1 p t s 1
V 2 5 (Delta p / Delta z) = (0.171 + 92.7 / Re ) (rho u / d ) ((1 / (1 - C sqrt(Fr)))) V V,eff eq 3where
u = u / (epsilon sin theta) V,eff V V V Re = (d rho u / eta ) V eq V,eff 2 Fr = u / d g L L eq
2 L We = u rho S/ sigma L L 2 Fr = u / (S g) L L L L Re = u S rho / eta L LThe effective g (as a function of h_t) is calculated:
L V L g = (1 - (dPdZ / dPdZ ) ) ((rho - rho ) / rho ) g eff floodThen F_t (see above), h_t, and dPdZ are computed
2/3 L V 1/3 h = (( 4 Ft / S) ) ((3 eta u / rho sin theta epsilon g )) t L eff V 2 2 A = (0.177 rho / S epsilon (sin theta) ) V 2 B = (88.774 eta / S epsilon sin theta) 2 5 (Delta p / Delta z) = (A u + B u ) ((1 / 1 - K h )) V V 2 tThe calculation is repeated until pressure drop converges or when it becomes larger than the pressure drop at flood. The iteration can actually have problems in convergencing.
Entrainment is computed from the fractional liquid entrainment which is computed from Hunt's correlation and from figure 5.11 of Lockett (1986) for sieve trays:
L -5 3.2 phi = 7.75 10 ( (0.073 / sigma) ) M ( (U / T - 2.5 h ) ) v v s clThe weeping factor is estimated from a figure from Smith (1963, plot on page 548), which was fitted with the following correlation
WF = (0.135 phi ln (34 (H +H ) + 1) / (H + H )) w ow d rwhere phi is the open area ratio.
Different design methods can be employed:
Tray layout parameters that specify a complete design (for the calculation of mass transfer coefficients and pressure drops) are shown in Table 5.5. For packings only the column diameter and bed height are design parameters, other parameters are fixed with the selection of the type of packing (such as void fraction, nominal packing diameter, etc.). The packed bed height must be specified since it determines the desired separation and the capacity.
Table 5.5. Tray layout data
General (sieve) tray layout data: | |
Column diameter | Active area |
Number of flow passes | Total hole area |
Tray spacing | Downcomer area |
Liquid flow path length | Weir length |
Hole diameter | Weir height |
Hole pitch | Deck thickness |
Downcomer clearance | |
Additional data for bubble caps: | |
Cap diameter | Slot area |
Slot height | Riser area |
Skirt clearance | Annual area |
Additional data for valves: | |
Closed Loss K | Open Loss K |
Eddy Loss C | Ratio Valve Legs |
Valve Density | Valve Thickness |
Fraction Heavy Valves | Heavy Valve Thickness |
The initial free area ratio is taken to be 15 % of the active area. The active area is determined with capacity factor calculation with internals specific methods (for sieve and bubble-cap trays the default is Fair's correlation by Ogboja and Kuye (19), and the Glitsch method for valve trays). The tray spacing is initially set to the default value (of 0.5 m) and the downcomer area is calculated according the Glitsch manual (limited by a minimum time residence check). From the combined areas the column diameter is computed. The number of liquid passes on a tray is initially set by the column diameter; under 5 ft one pass, under 8 ft two, 10 ft three, under 13 ft four, else five passes. With the number of passes and the column diameter the total weir length is computed. Once the weir length is determined the liquid weir load is checked, if too high the number of passes is incremented and a new weir length is evaluated until the weir load is below a specified maximum.
Initial weir height is taken as 2", but limited to a maximum of 15 % of the tray spacing. For notched or serrated weirs the notch depth is a third of the weir height. For serrated weirs the angle of serration is 45 degrees. Circular weirs have diameters 0.9 times the weir length. Hole diameter is set to 3/16" for sieve trays and tray thickness 0.43 times the hole diameter (or 1/10"). The hole pitch is computed from the free area ratio and hole diameter according to a triangular pitch. The default downcomer clearance is 1.5" but is limited by the maximum allowed downcomer velocity according to the Glitch method de-rated with the system factor. The clearance is set to be at least half an inch lower than the weir height to maintain a positive liquid seal but is limited to a minimum of half an inch.
For bubble-cap trays the cap diameter is 3" for column diameters below 4.5 ft and 4" for above. The hole diameter can vary between 60 % to 71 % of the capdiameter, and default taken as 70 %. Default skirt clearance is 1" with minimum of 0.5" and maximum of 1.5". slot height can vary in between 0.5" and 1.5", default 1" for cap diameters below 3.5" and 1.25" for larger cap diameters. The pitch can vary from 1.25" to half the flow path length (minimum number of rows is two), default set to 1.25".
Valve trays are initialized to be Venturi orifice uncaged, carbon steel valves of 3 mm thick with 3 legs (see Kister, 1992, p312). The hole diameter is 1" for column smaller than 4.5 ft, otherwise 2". No double weight valves are present.
The second task in the fraction of flooding method consists of finding the proper free area ratio (beta = A_h / A_b = hole area / active area) so that no weeping occurs. This ratio can vary between a minimum of 5% (for stable operation) and a maximum of 20%. To test whether weeping occurs, we use the correlation by Lockett and Banik (1984): Fr_hole>2/3. The method requires all liquid heights to be evaluated at weep rate conditions. This task is ignored for bubble-cap trays. The weep test is done at weeping conditions, with a weep factor at 60 % (this can be changed). Calculating liquid heights is done by adding various contributions with correlations from Lockett (1986) and Kister (1992), see Appendix A. If weeping occurs at the lower bound for the free area ratio, a flag is set for the final task to adapt the design. The final task consists of evaluating all liquid heights at normal conditions and to do a number of checks:
Table 5.6. Tray design checks and adjustments
Problem | Test | Adjustments |
Bubble cap vapor distribution | h_lg/h_d > 0.5 | p+f_1, |
h_skirt+f_2, | ||
h_slot+f_3, | ||
d_h-f_3 | ||
Weeping | Fr_h / (2/3) < 1-f_a | |
free < 0.05 | ||
A_b < A_bf: | A_b=A_bf | |
W_f+f_1 | ||
else: | A_b-f_1 | |
d_h-f_3 | ||
h_w-f_3 | ||
t_v+f_2 (vt) | ||
Hydrodynamic (downcomer) flooding) | T_s < h_db/FF | T_s+f_1 |
A_d+f_1 | ||
h_w+f_2 | ||
h_c+f_3 | ||
Excessive liquid entrainment | A_b+f_1 | |
T_s+f_1 | ||
d_h-f_2 | ||
h_w-f_3 | ||
Froth height limit | h_f > 0.75 T_s | A_b+f_1 |
T_s+f_2 | ||
h_w-f_3 | ||
Excessive pressure drop | g rho h_wt > Delta p_max | A_b+f_1 |
h_w-f_1 | ||
d_h+f_2 | ||
p+f_1 (bc) | ||
h_skirt+f_2 (bc) | ||
h_slot+f_3 (bc) | ||
Excessive vapor entrainment | A_d+f_1 | |
Here we discuss the most important parameters of the file. The file starts with a comment on the first line. The second line specifies the factors f_1, f_2, and f_3 for adapting the design layout parameters. The third line specifies the fraction of change allowed in the flows before a redesign occurs. It also specifies the fraction of deviation allowed in downcomer and bubbling area between current and design values. Line 11 specifies the volumetric weir load after which the number of passes is incremented. Line 14 specifies the maximum froth height as fraction of the tray spacing as is used in the froth height limit check. Line 15 specifies the criterium to which the free area ratio has to conmverge. Line 16 sets the maximum allowed pressure drop for the excessive pressure drop check. Line 20 specifies the maximum number of loops for the design method. Line 21-23 specify the methods to calculate capacity factors for bubble-cap, sieve, and valve trays. Line 24-25 set the downcomer are method and velocity check. Line 45 sets a flag to generate tray parameter output and line 46 sets a flag for intermediate design messages
To determine the packed column diameter, the diameter that gives rise to the flooding pressure drop (as specified) is computed using the selected pressure drop model. The resulting diameter is corrected for the fraction of flooding and the system factor:
D = (D / sqrt(FF~SF)) c c,floodThis does make the resulting column diameter depend on the selected pressure drop model. If no pressure drop model is selected the Leva (1992) model is selected (which is only a function of the packing factor). If no pressure drop at flood is specified, it is estimated with Kister's correlation (1992) (which is only a function of the packing factor). Thus, as long as the packing factor is known, this method will not fail.
Packing design automatically finds the diameter resulting in the specified pressure drop (with the selected pressure drop model). This is done by using a linear search technique as the different packing pressure drop correlations can behave quite irregularly. The maximum allowed pressure drop is the flooding pressure drop as specified or computed from Kister's correlation and the packing factor. If the pressure drop is specified to be very low the column diameter might converge to unrealistic diameters. A zero or larger than flooding pressure drop specification results in a 70 % fraction of flooding design.
a_d Interfacial area density (m^2/m^3) a^I Interfacial area (m^2) A_h Hole area (m^2) A_b, A_bub Bubbling area (m^2) A_d Downcomer area (m^2) c Number of components, Molar concentration (kmol/m^3) d_h Hole diameter (m) D Binary diffusivity coefficient (m^2/s) D_c Column diameter (m) D_e Eddy dispersion coefficient (m^2/s) e Energy transfer rate (J/s) f_ij Component i feed flow to stage j (kmol/s) f_1, f_2, f_3 Design adjustment factors F_j Total feed flow rate to stage j (kmol/s) F_p Packing factor (1/m) F_s F factor F_s=U_v sqrt(\rho_V) (kg^0.5/m^0.5s) FF Fraction of flooding FP Flow parameter FP=M_L/M_V sqrt(\rho^V_t/\rho^L_t) Fr Froude number g Gravitational constant, 9.81 (m/s^2) G Interlinked flow rate (kmol/s) h Heat transfer coefficient (J/m^2 K s) h_c Clearance height under downcomer (m) h_cl Clear liquid height (m) h_d Dry tray pressure drop height (m) h_db Downcomer backup liquid height (m) h_f Froth height (m) h_i Liquid height at tray inlet (m) h_lg Liquid gradient pressure drop height (m) h_l, h_L Liquid pressure drop height (m) h_ow Height of liquid over weir (m) h_r Residual pressure drop liquid height (m) h_wt Wet tray pressure drop liquid height (m) h_w Weir height (m) h_udc Liquid height pressure loss under downcomer (m) H Molar enthalpy (J/kmol) H_i Partial molar enthalpy of component i (J/kmol) J Molar diffusion flux (kmol/m^2 s) k Binary mass transfer coefficient (m/s) K_i K-value or equilibrium ratio component i: K_i=y_i/x_i L Liquid flow rate (kmol/s) Le Lewis number ( Le= Sc/ Pr) M Mass flow rate (kg/s) N Mass transfer rate (kmol/s) n Number of stages p Hole pitch (m), Pressure (Pa) Delta p Pressure drop (Pa) Delta P_max Maximum design pressure drop (Pa/tray or Pa/m) Pr Prandtl number Q Heat input (J/s) Q_L Volumetric flow over the weir (m^3/s) r Ratio sidestream to internal flow [R] Matrix defined by () and () Sc Schmidt number SF System derating factor t Residence time (s) t_v Valve thickness (m) T Temperature (K) T_s Tray spacing (m) V Vapor flow rate (kmol/s) We Weber number W_l Weir length (m) x Liquid mole fraction y Vapor mole fraction z Mole fraction Greek: alpha Fraction liquid in froth beta Fractional free area beta=A_h/A_b, phi Fractional entrainment rho Density (kg/m^3) sigma Surface tension (N/m) eta Viscosity (Pa s) [Gamma] Thermodynamic matrix lambda Heat conductivity (W/m/K) Superscripts: I Interfacial L Liquid P Phase P V Vapor Subscripts: flood at flooding conditions i component i j stage j, component j spec specified t total nu from interlinking stage nu
S.D. Barnicki, J.F. Davis, ``Designing Sieve-Tray Columns, Part 1: Tray Design'', Chem. Engng., Vol. 96, No. 10, (1989) pp. 140--146.
S.D. Barnicki, J.F. Davis, ``Designing Sieve-Tray Columns, Part 2: Column Design and Verification'', Chem. Engng., November, pp. 202--212 (1989).
D.L. Bennett, R. Agrawa, P.J. Cook, ``New Pressure Drop Correlation fo Sieve Tray Distillation Columns'', AIChE J., Vol. 29, No. 3 (1983) pp. 434.
R. Billet, M. Schultes, "Advantage in correlating packed column perfomance", IChemE. Symp. Ser. No. 128, B129 (1992).
R. Billet, Distillation Engineering?, Heyden (1979?).
W.L. Bolles, in B.D. Smith, Design of Equilibrium Stage Processes, Chap. 14, McGraw-Hill (1963).
J.L Bravo, J.R. Fair, Ind. Eng. Chem. Proc. Dev., 21, 163 (1982).
J.L. Bravo, J.A. Rocha, J.R. Fair, "Mass transfer in gauze packings", Hydrocarbon Processing, January, 91 (1985).
J.L. Bravo, J.A. Rocha, J.R. Fair, "Pressure drop in structured packings", Hydrocarbon Processing, March (1986).
J.L. Bravo, J.A. Rocha, J.R. Fair, "A comprehensive model for the performance of columns containing structured packings", IChemE. Symp. Ser. No. 128, A439 (1992).
Chan, J.R. Fair, ``Prediction of point efficiencies on sieve trays'', Ind. Eng. Proc. Des. Dev., Vol. bf 23, 814 (1984)
Chen and Chuang, Ind. Eng. Chem. Res, Vol. 32, 701--708 (1993).
Gerster et al., AIChE J., (1958).
I.J. Harris, `` Optimum Design of Sieve Trays'', Brit. Chem. Engng, Vol. 10, No. 6 (1965) pp. 377.
G.A. Hughmark, ``Mdels for Vapour Phase and Liquid Phase Mass Transfer on Distillation Trays'', AIChE J., Vol. 17, No. 6 (1971) pp. 1295.
F. Kastanek, ``Efficiencies of Different Types of Distillation Plate'', Coll. Czech. Chem. Comm., Vol. 35 (1970) pp. 1170.
H.Z. Kister, Distillation Design, McGraw-Hill, New York (1992).
G.F. Klein, Chem. Engng, May 3, 81 (1982).
H.A. Kooijman, R. Taylor, ``On the Estimation of Diffusion Coefficients in Multicomponent Liquid Systems'', Ind. Eng. Chem. Res., Vol 30, No. 6, (1991) pp. 1217--1222.
W.J. Korchinsky, ``Liquid Mixing in Distillation Trays: Simultaneous Measurement of the Diffusion Coefficient and Point Efficiency'', Trans. I. Chem. E., Vol. 72, Part A, (1994) 472-478.
R. Krishnamurthy, R. Taylor, ``A Nonequilibrium Stage Model of Multicomponent Separation Processes. Part I: Model Description and Method of Solution'', AIChE J., Vol. 31, No. 3 (1985), pp. 449--455.
Leva, (1953?).
M. Leva, ``Reconsider Packed-Tower Pressure-Drop correlations", Chem. Eng. Prog. January, 65 (1992).
M.J. Lockett, Distillation Tray Fundamentals, Cambridge University Press (1986).
M.J. Lockett, S. Banik, ``Weeping from Sieve Trays'', AIChE Meeting, San Francisco, Nov. (1984).
E.E. Ludwig, Applied Process Design for Chemical and Petrochemical Plants, Vol. 2, 2nd Ed., Gulf Pub. Co., Houston, TX, (1979).
O. Ogboja, A. Kuye, ``A procedure for the design and optimisation of sieve trays'', Trans. I. Chem. E., Vol. 68, Part A, Sept. (1990) pp. 445-452.
K. Onda, H. Takeuchi, Y. Okumoto, "Mass transfer coefficients between gas and liquid phases in packed columns", J. Chem. Eng. Jap., Vol. 1, No.1, 56 (1968).
R.H. Perry and D. Green, Perry's Chemical Engineering Handbook, 6th edition, section 18, Liquid-Gas System, 18-8 -- 18-12 (1984).
M. Prado, The bubble-to-Spray Transition on Sieve Trays: Mechanisms of the Phase Inversion, Ph.D. thesis, University of Texas, Austin (1986).
B.D. Smith, Design of Equilibrium Staged Processes, McGraw-Hill, New York (1963)
J. Stichlmair, A. Mersmann, ``Dimensioning Plate Columns for Absorption and Rectification'', Chem. Ing. Tech., Vol. 45, No. 5 (1978) pp. 242.
J. Stichlmair, J.L. Bravo, J.R. Fair, Gas. Sep. Purif., Vol. 3, 19 (1989).
R. Taylor, H.A. Kooijman, ``Composition derivatives of Activity Models (for the estimation of Thermodynamic Factors in Diffusion)'', Chem. Eng. Comm., Vol. 102 (1991) pp. 87--106.
R. Taylor, H.A. Kooijman, J-S. Hung, ``A second generation nonequilibrium model for computer simulation of multicomponent separation processes'', Comput. Chem. Engng., Vol. 18, No. 3, pp. 205--217 (1994).
R. Taylor, R. Krishna, Multicomponent Mass Transfer, Wiley, New York (1993).
P.C. Wankat, Separations in Chemical Engineering - Equilibrium Staged Separations, Elsevier, 420--428 (1988).
F.J. Zuiderweg, ``Sieve Trays - A View of the State of the Art'', Chem. Eng. Sci., 37, 1441--1461 (1982).
Current limitations of the nonequilibrium extraction model are:
x = sqrt(\sigma \over \Delta \rho g) d = 1.8 x hbut d_h is limited (if supplied) by:
0.5 x < d < pi x hand the practical limits (overriding):
3 mm < d < 8 mm hThe hole velocity is computed with:
2 (Eo) = (Delta rho g d / sigma) h -0.26 (We) = 4.33 (Eo) U = sqrt( We \sigma \over \rho d ) h d hIf the hole velocity is less than 0.15 (m/s) then its design value is kept at 0.15 (m/s). The Froude number is computed from
2 (Fr) = (U / g d ) h hFor Eo is less than 0.4 the Sauter mean droplet diameter is computed by:
-0.4 0.67 d = (Eo) ( 2.13 ( (Delta rho / rho ) ) + exp (-0.13 (Fr)) ) d p d hotherwise
-0.42 0.42 d = (Eo) ( 1.24 + exp (-(Fr) ) ) d p hThe hole area is
A = (Q / U ) h d hThe ratio of the hole area over the active area (free area ratio, f) is limited between 1 and 20 %.
A = A / f a hThe hole pitch can be computed if the hole diameter and free area ratio are known. The downcomer velocity can be computed if a minimum droplet diameter, d_min, is assumed which will not be entrained. The downcomer velocity equals the velocity of the continuous phase, U_c:
2 0.33 U = 0.249 d ( ((g Delta rho) / rho eta ) ) c min c cThis droplet diameter is taken to be 0.5 mm. With U_c known we can compute the downcomer area:
A = Q / U d c cThe total area is equal to two downcomer areas plus the active area and 0.5 % area for support etc.:
A = (A + 2 A ) / 0.995 t a dWith the total tray area known the column diameter can be computed. The net area for the disperse phase, A_n, and the disperse velocity, U_d, are:
A = A + A n A d U = (Q / A ) d d nNext the dispersed phase velocity holdup and slip velocity are computed. The slip velocity (V_s) is guessed at one tenth of the disperse phase velocity, making the disperse phase holdup equal to a tenth since it is defined as
phi = (U / V ) d d sThe slip velocity (which is a function of the dispersed phase holdup and needs to be obtained iteratively) can be calculated from:
0.33 1.834 V = sqrt( 2.725 g d \left( (\Delta \rho \over \rho ) \right) \left( ( (1 - \phi ) \over (1 + \phi ) ) \right) ) s p c d dAfter the dispersed phase holdup is computed (it depends on V_s) it is checked to be within 1 and 20 % for standard operation conditions. If it is too small the free area ratio is increased, if it is too large the free area ratio is decreased (each by 5 %) till it is within the desired range.
The Weber number
2 (We) = rho U d / sigma d h pmust be larger than 2 to ensure all holes produce drops (i.e. to avoid inactive holes, see Seibert and Fair, 1988).
The height of the coalesced layer is (according to Treybal, 1980):
2 2 2 h = ((U - U ) rho / 2 g C Delta rho) + (4.5 U rho / 2 g Delta rho) + (6 sigma / d g Delta rho) c h d d d c c p(with C_d=0.67). The first term is height to overcome flow through the orrifices, the second for friction losses due to contraction/expension on entry/exit (0.5 + 1.0) and change of direction (2 times 1.5 velocity heads), and the third term for the interfacial tension effects at the holes. The height needs to be larger than 4 cm (to ensure safe operation). If not, then the hole diameter is decreased by 5 % and we repeat the procedure from the hole velocity calculation (). This design is for a one pass sieve tray, and flow path length, L_f is computed from geometric relations. The weir length is (segmental downcomer):
W = A / L l a fThe tray thickness is defaulted to a tenth of an inch. To prevent entrainment of droplets, the flow under the downcomer is only allowed to be 50 % higher than the downcomer velocity. If higher, then the downcomer clearance is enlarged until this requirement is met. The tray spacing is adjusted so that the coeleseced layer and coalescence zone divided over the length of the downcomer equals the fraction of flooding (multiplied with the system factor).
A = (6 phi / d ) i d pand the drop rising zone (where mass transfer is assumed to take place):
h = t - h drop s cwhere t_s is th etray spacing. The volume for mass transfer on a tray is
V = A h i n dropThe mass transfer coefficients for transport from the disperse phase are (Handlos and Baron, 1957):
k = (0.00375 V / ( 1 + eta / eta ) ) d s d cand for transport from the continuous phase are (Treybal, 1963):
-0.43 -0.58 k = 0.725 (Re) (1 - phi ) V (Nu) c,ij c d s cwith
(Re) = d V rho / eta c p s c c c (Nu) = eta / rho D c c c ijNote that k_d is not a function of the diffusion coefficient and, thus, is the same for all components.
d = 1.15 x sqrt( \sigma \over \Delta \rho g) pThe slip velocity of a single droplet at zero disperse phase holdup is given by
o V = sqrt( 4 \Delta \rho g d \over 3 \rho C ) s p c dwhere Cd=0.38 (for high values of Reynolds). Static disperse holdup is:
phi = (0.076 a d / xi) ds p pwhere a_p is the packing area and xi the packing void fraction. The static holdup area and total area are:
a = 6 0.076 a s p a = a + a p sThe tortuosity is defined as
zeta = (a d / 2) pThe superficial velocity of the continuous phase at the flood point is
e = cos ((pi zeta / 4)) o 2 U = (0.192 xi * V / (1.08 + (Q / Q ) / e )) cf s d cThis needs to be corrected for the fraction of flooding (and system factor):
U = SF FF U c cfto give the net area
A = (Q / U ) n c cfrom which the packed column diameter (D_c) can be calculated.
log ( E / V d ) = 0.046 (V / V ) + 0.301 d d p c d log ( E / V d ) = 0.161 (V / V ) + 0.347 c c p c dwhere d_p is the packing diameter.
U = (Q / A ) c c n U = (Q / A ) d d nThe drop diameter d_p, slip velocity V^o_s, area a, static holdup area a_s, and tortuosity zeta are calculated as above. Then the disperse phase holdup, phi_d, is determined iteratively (starting at 0.1.) from:
f(phi ) = exp ( (-6 phi / pi)) d d o 2 phi = (U / xi (V f(phi ) - U ) e ) d d s d cThen the slip velocity is
o V = V f(phi ) e + (1 - e) U s s d csince U_d=V^o_s f(phi_d). The mass transfer coefficient for the disperse phase is computed by:
phi = (sqrt(Sc ) / ( 1 + eta / eta ) ) d d c phi > 6 : k = (0.023 V / sqrt((Sc) )) d,ij s d phi < 6 : k = (0.00375 V / ( 1 + eta / eta ) ) d s d c (Sc) = ( eta / rho D ) d d d d,ijIf phi is larger than 6 the Laddha and Degaleesan correlation is used otherwise the Handlos-Baron method. For the mass transfer coefficient in the continuous phase:
0.5 0.4 (Sh) = 0.698 (Re) (Sc) (1-phi ) c c c d k = ( (Sh) D / d ) c,ij c c,ij p (Re) = ( rho V d / eta ) c c s p c (Sc) = ( eta / rho D ) c c c c,ijThe interfacial area per unit volume is:
a = ( 6 xi phi / d ) i d pThe total interfacial area in a stage is the stage height times the net area times the interfacial area per unit volume:
a = a A h i,tot i n stage
alpha = ( Q / Q ) d cThe maximum stable drop diameter is
5/21 6/21 10/21 1/21 u = 0.9 ( ( g Delta rho ) sigma / rho rho ) 0 c d 2 d = (sigma / rho u ) p,max c 0A stable drop diameter is selected as half of the maximum diameter
d = 0.5 d p p,maxand the require power input (P_i=N^3 R^5 / H D^2) is computed
0.6 2.5 e = ( (0.25 ((sigma / rho ) ) / d ) ) c p P = (pi e / 4 C ) i p(C_p=0.03 for Re>10^5). If no column diameter is known, an estimate is made from assuming a cross-sectional area for a combined velocity of 0.05 m/s with:
A = (Q + Q ) / 0.05 c c dThe required rotational speed (using these standard ratios) is then
5 2 0.33 N = ( ( P ((0.1 / (0.6) ) ) / D ) ) i cNow the slip-velocity can be calculated using a correlation from Kung and Beckman (1961):
0.9 2.3 0.9 2.6 2 (V * eta / sigma) = (( Delta rho / rho )) ((S / R)) ((H / R)) ((R / D)) ((g / R N )) s c cThe disperse holdup at flood is determined from
2 phi = (sqrt(\alpha +8 \alpha) - 3 alpha / 4 ( 1 - alpha)) dfrom which the continuous phase velocity at flood can be determined with
2 U = V (1 - phi ) (1 - 2 phi ) c,f s d dCorrection for fraction of flooding (and system factor) gives
U = SF FF U c cffrom which the column area and diamater can be calculated
A = ( Q / U ) c c cThe rotor diameter R, stator diameter S, and the height of the compartment have standard ratios with respect to the column diameter (D_c)
R = 0.6 D c S = 0.7 D c H = 0.1 D cso the size of the column is determined. Below a Renolds number of 10^5 C_p becomes a function of the Renolds number. Normally RDC's are operated in the regime above 10^5 so the Renolds number is computed by
2 Re = (rho N R / eta ) d d dand a smaller diameter is selected (and the calculations repeated) if necessary. On re-design the layout of the stage with the largest diameter is used for the entire section.
5 2 (( 10 / Re )) dStemerding et al. (1963) gave a correlation for the axial dispersion coefficient for the continuous phase
2 (E / V H) = 0.5 + 0.012 N R (S/D) / V c c cThe disperse dispersion coefficient is set to twice this number.
C = 0.03 p e = (4 C P / pi) p i 0.6 0.4 d = 0.25 ((sigma / e)) / rho p cThe dispersed holdup phi_d is calculated iteratively as above and the slip velocity is determined as described above (with ). The mass transfer coefficients are:
(Sh) = 10.0 d k = ( (Sh) D / d ) d,ij d d,ij p 0.62 0.36 (Sh) = 2 + 0.42 (Re) (Sc) c c c k = ( (Sh) D / d ) c,ij c c,ij pwith
1.33 0.33 (Re) = ( rho d e / eta ) c c p c (Sc) = ( eta / rho D ) c c c c,ijThe interfacial area per unit volume is
a = ( 6 phi / d ) i d pAlternatively the "Rose-Kintner-Garner-Tayeban" method can be used:
(Sh) = 0.6 sqrt((Re) ) sqrt((Sc) ) c c c 0.225 b = d / 1.242 p w = (8 sigma b / dp) (n (n+1) (n-1) (n+2) / (n+1) rho + n rho ) d c k = 0.45 sqrt(D \omega) d d,ijwhere n=2, and d_p is in cm for the calculation of b and w.
A = Q / U o d oThe droplet diameter can be calculated from (Vedaiyan et al, 1972):
2 -0.067 d = 1.592 (( U / 2 g d )) sqrt( \sigma \over g \Delta \rho) p o oThe flood velocity of the continuous phase is (Treybal, 1963):
0.28 0.075 0.056 2 U = (0.3894 Delta rho / [ 0.2165 eta sqrt(\rho ) + 0.2670 d sqrt(\rho \alpha) ] ) cf c c p dwhere alpha=Q_d/Q_c. The disperse holdup at flood is
2 phi = (sqrt(\alpha + 8 \alpha) - 3 alpha / 4 ( 1 - alpha )) dfThe velocity of the continuous phase is then
U = FF SF U c cfand the column area
A = Q / U c c cfrom which the column diameter can be calculated (The column area must also be larger then the total hole area, if not, the column area is set to four times the hole area).
FF = (U / SF U ) c cfwhere U_cf is computed as in the spray column design and U_c=A_c/Q_c. No lower operating limit is calculated. The dispersion coefficient for the continuous phase is (Vermeulen et al., 1966):
(E / V H) = 7.2 sqrt(U D ) c c d cSince the dispersion coefficient for the disperse phase is unknown it is set equal to that for the continuous phase.
2 4 4 P = (rho sigma / g eta Delta rho) c c 0.15 d = 7.25 sqrt( \sigma \over g \Delta \rho P ) p,tThe drop terminal velocity is (Satish et al., 1974):
2 -0.082 2 1/4 V = 1.088 (( U / 2 g d )) (( sigma g Delta rho / rho )) t o o cWith the continous operating and flood velocities the fraction of flooding is calculated and then the disperse phase holdup is
phi = FF phi d dfand the slip velocity
V = (1 - phi ) V s d tIf the drops are stagnent (d_p < d_p,t) the disperse MTC is computed from
k = 18.9 D / d d,ij d,ij pelse the Handlos-Baron correlation (1957) is used:
k = (0.00375 V / (1 + eta / eta )) d s d cFor the continuous phase MTC we use (Ruby and Elgin, 1955)
-0.43 -0.58 k = 0.725 Re Sc (1 - phi ) V c c c d swhere
Re = (d V rho / eta ) c p s c c Sc = (eta / rho D ) c c c c,ijThe interfacial area for mass transfer per unit of volume is
A = (6 phi / d ) i d p
alpha = (E / V H) - 0.5 d d d alpha = (E / V H) - 0.5 c c cwhere alpha is the fractional backflow ("entrainment") in the stage, and H is the stage height.
For spray columns (Perry, 198x):
E = 7.2 sqrt( V D ) c d cFor packed columns:
log (E / V d ) = 0.046 (V / V ) + 0.301 d d F c d log (E / V d ) = 0.161 (V / V ) + 0.347 c c F c dFor a RDC:
E = 0.5 H V + 0.012 R N H ( ( S / D ) ) c c c E = F E d cwhere F is calculated by
5 2 3.3 F = (4.2 10 / D ) ( ( V / h ) ) c dand must be larger or equal than one. Krishna uses:
E = (0.5 H U / (1-Phi )) + 0.012 R N H ( ( S / D ) ) c c d c E = (0.5 H U / Phi ) + 0.024 R N H ( ( S / D ) ) d d d c
a_p Packing area per unit volume (m^2/m^3) a_s Static holdup area per unit volume (m^2/m^3) A_a Total tray active area (m^2) A_d Downcomer area (m^2) A_i Interfacial area per unit volume (m^2/m^3) A_h Total tray hole area (m^2) A_n Netto tray area (m^2), A_n=A_a+A_d A_t Total tray area (m^2) C_d Drag coefficient D Binary diffusion coefficient (m^2/s) D_c Column diameter (m) d_e Effective drop diameter (m) ? d_h Hole diameter (m) d_min Minimum droplet diameter (m) d_p Sauter mean drop diameter (m) Eo Eotvos number ( Delta rho g d_h / sigma) f Free area ratio (A_h/A_a) F Molar flow (kmol/s) Fr Froude number (U^2_h / g d_h) FF Fraction of flooding g, g_c Gravitational constant, 9.81 (m/s^2) H RDC compartment height (m) h_c Height of coalesced layer (m) h, h_drop Height of drop rising zone (m) h_stage Stage height for packed column (m) k Binary mass transfer coefficient (m/s) M_w Molecular weight (kg/kmol) N Rotation speed (rad/s) Nu Nusselt number Pe Peclet number P_i Power input (?) Q Volumetric flow (m^3/s) R Rotor diameter (m) Re Reynolds number S Inner stator diameter (m) Sc Schmidt number Sh Sherwood number SF System derating factor t Contact time (s) t_s Tray spacing (m) U_c,U_d Continuous, disperse velocity (m/s) U_cf Continuous phase superficial velocity at flood (m/s) U_h Hole diameter (m/s) V_i Tray volume for interfacial mass transport (m^3) V_s Slip velocity (m/s) V^o_s Slip velocity at zero disperse phase holdup (m/s) We Weber number ( rho_d U^2_h d_p / sigma) W_l Weir length (m) Greek: alpha Phase ratio (Q_d / Q_c) rho Mass density (kg/m^3) phi_d Disperse phase holdup fraction phi_ds Static disperse phase holdup fraction sigma Interfacial tension (N/m) eta Liquid viscosity (Pa.s) mu Kinematic viscosity (eta / rho) zeta Tortuosity kappa delta Subscripts: c Continuous phase d Disperse phase, Downcomer i Interface, Component i j Component j
R.M. Griffith, Chem. Eng. Sci., 12, 198 (1960).
A.E. Handlos, T. Baron, ``Mass and Heat Transfer from Drops in Liquid-Liquid Extraction'', AIChE J., 3 (1957) pp. 127--136.
A.E. Handlos, T. Baron, AIChE J., 6, 145 (1957).
Hughmark, Ind. eng. Chem. Fundam., 6, 408 (1967).
D.G. Jordan, Chemical Process Development, Part 2, John Wiley, New York (1968).
W.J. Korchinsky, "Liquid-Liquid Extraction Column Modelling: Is the Forward Mixing Influence Necessary?", Trans. I. Chem. E., Vol. 70, Part A, 333--345 (1992).
R. Krishna, S.M. Nanoti, A.N. Goswami, "Mass-Transfer Efficiency of Sieve Tray Extraction Columns", Ind. Eng. Chem. Res., Vol. 28 (1989) 642-644.
R. Krishna, Design of Liquid-Liquid Extraction Columns, University of Amsterdam (NL), (1993).
R. Kronig, J.C. Brink, Appl. Sci Res., A2, 142 (1950). A. Kumar, S. Hartland, "Prediction of Axial Mixing Coefficients in Rotating Disc and Asymmetric Rotating Disc Extraction Columns", Can. J. Chem. Eng, Vol. 70, 77--87 (1992).
A. Kumar, S. Hartland, "Prediction of drop size, dispersed-phase holdup, slip velocity, and limiting throughputs in packed extraction columns", Trans. IChemE., 72, Part A, 89--104 (1994).
M. Lao et al., "A Nonequilibrium Stage Model of Multicomponent Separation Processes VI: Simulation of Liquid-Liquid Extraction", Chem. Eng. Comm., 86, p73--89 (1989).
G.S. Laddha, T.E. Degaleesan, Transport Phenomena in Liquid Extraction, McGraw-Hill (1978).
T.C. Lo, M.H.I. Baird, C. Hanson, Handbook of Solvent Extraction, John Wiley, NY (1983).
J.A. Rocha, J.L. Humphrey, J.R. Fair, "Mass transfer Efficiency of Sieve Tray Extractors", Ind. Eng. Chem. Process Des. Dev., Vol. 25, No. 4 (1986) pp. 862--871.
P.N. Rowe, K.T. Claxton, J.B. Lewis, Transact. Inst. Chem. Eng., 43, T14 (1965).
P.M. Rose, R.C. Kinter, ``Mass Transfer from Large Oscillating Drops'', AIChE J., Vol. 12, No. 3 (1966) pp. 530.
E.Y. Kung, R.B. Beckman, ``Dispersed-Phase Holdup in a Rotating Disk Extraction Column'', AIChE J., Vol. 7, No. 2 (1961) pp. 319-324.
A.M. Rozen, A.I. Bezzubova, Theor. found. Chem. Eng., 2, 715 (1968); translated from Teor. Osnovy Khim. Tekh, 2, 850 (1968).
Ruby, Elgin, Chem. Eng. Prog., 51, Sump. Ser. 16, 17 (1955).
L. Satish, T.E. Dagaleesan, G.S. Laddha, Indian Chem. Eng., Vol. 16 (1974) pp. 36.
A.F. Seibert, J.R. Fair, "Hydrodynamics and Mass Transfer in Spray and Packed Liquid-Liquid Extraction Columns", Ind. Eng. Chem. Res., 27, No.3, 470--481 (1988).
A.H.P. Skelland, R.M. Wellek, AIChE J., 10, 491 (1964).
A.H.P. Skelland, W.L. Conger, Ind. Eng. Chem. Process De. Devel., 12, 448 (1973).
A.H.P. Skelland, D.W. Tedder, Handbook of Separation Processes, Ed. R.W.Rousseau, Wiley (1987).
S. Stemerding, E.C. Lumb, J. Lips, ``Axiale Vermischung in einer Drehscheiben-Extraktions Kolonne'', Chem. ing. Tech, Vol. 35 (1963) pp. 844--850.
G. Thorsen, S.G. Terjesen, Chem .Eng. Sci., 17, 137 (1962).
R.E. Treybal, Mass Transfer Operations, 3rd ed., McGraw-Hill, New York (1980)
R.E. Treybal, Liquid Extraction, 2nd ed., McGraw-Hill, New York (1963).
S. Vedaiyan, T.E. Degaleesan, G.S. Laddha, HE. Hoelscher, AIChE J., Vol. 18 (1972) pp. 161.
Vermeulen et al., Chem. Eng. Prog., Vol. 62, No. 9 (1966) pp. 95.
M.E. Weber, Ind. Eng. Chem. Fund., 14, 165 (1975).
0.5 0.33 (Sh) = A + B Re (Sc) c d cwith A=2 and B=0.79.
-0/5 0.64 (Sh) = (2 / pi) sqrt( 1 - (Re) (2.89 + 2.15 \mu )) sqrt(Pe) c d rwhere
(Pe) = (Re) (Sc) d c
0.5 0.42 (Sh) = -126 + 1.8 Re (Sc) c d c
0.5 0.33 (Sh) = -178 + 3.62 Re (Sc) c d c
(Sh) = (k phi d / D ) = 16.7 d d d p dfor Re_d < 50.
k = (0.00375 V / 1 + (mu / mu )) d s d c
0.68 3 2 4 0.10 2 -0.14 (Sh) = 0.32 (Re) ( (sigma rho / g mu Delta rho) ) ( (4 D t / d ) ) d c c d,ij c d
0.63 0.50 -0.5 a: (Sh) = 0.32 (Re) (Sc) ( 1 + (mu / mu ) ) d d d d c -5 2.0 0.56 -0.5 b: (Sh) = 7.5 . 10 (Re) (Sc) ( 1 + (mu / mu ) ) d d d d cfor medium (a) and large (b) droplets.
(Re) = ( d rho V / eta ) p s (Sc) = ( eta / rho D ) ij (Sh) = ( k d / D ) ij ij p ij mu = ( eta / rho )Table 1 of chapter 10 in the Handbook of Solvent Extraction supplies us with three more models for the drop rise zone. One for stagnant drops (Skelland and Conger, 1973):
0.5 0.5 k = - (( d / 6 t )) (( rho / M )) ln ( 1 - ( pi D t / 0.5 d ) ) d e d d av vd e 0.5 0.333 k = 0.74 (( D / d )) (( rho / M )) (( d V rho / mu )) (( mu / rho D )) c vc e c c av e s c c c c vcfor circulating drops (Treybal, 1963):
2 -0.34 -0.125 2 0.37 k = 31.4 (( D / d )) (( rho / M )) (( 4 D t / d )) (( mu / rho D )) (( d V rho / sigma )) d vd e d d av vd e d d vd e s c -0.34 -0.58 k = 0.725 (( rho / M )) (( d V rho / mu )) (( mu / rho D )) V ( 1 - phi ) c c c av e s c c c c vc s dand for oscillating drops (Skelland and Conger, 1973):
2 -0.14 0.68 3 2 4 0.10 k = 0.32 (( D / d )) (( rho / M )) (( 4 D t / d )) (( d V rho / mu )) (( sigma rho / mu g Delta rho )) d vd e d d av vd e e s c c c -0.34 1.0 0.7 k = (( D / d )) (( rho / M )) (( d V rho / mu )) [( 50 + 0.0085 (( d V rho / mu )) (( mu / rho D )) )] c vc e c c av e s c c e s c c c c vcwhere t=h/V_s (with h as the height of the drop rise zone) and V_s=V_t(1-phi_d). Perry's also supplies us with some more correlations. There we find that () is from Ruby and Elgin (1955) and is to be applied for circulating drops. Another correlation for the continuous mass transfer coefficients for circulating drops is by Hughmark (1967):
0 0.339 1/3 1/3 0.072 (k d / D ) = [ 2 + 0.463 (Re) .484 Sc (( d g / D )) ] F c p c c p c 2 3 F = 0.281 + 1.615 kappa + 3.73 kappa -1.874 kappa 1/8 1/4 1/6 kappa = Re (( mu / mu )) (( mu V / sigma g )) c d c s cwhere Re is the droplet Reynolds number. A correlation for the disperse mass transfer coefficient for oscillating droplets by Rose and Kinter (1966) is:
2 k = sqrt( (4 D \omega \over \pi) (1 + \delta + (3 \over 8) \delta ) ) d d 3 w = (1 / 2 pi) sqrt( 192 \sigma g b \over d (3\rho +2\rho ) ) c p d c 0.225 b = 1.052 d pwhere delta can be taken as 0.2 if unknown.
2 N = (( N D rho / mu )) 1 r c c 2 N = (( N D / g )) 2 r N = (( mu / sqrt( \sigma \rho D ) )) 3 c c r N = (( rho / rho )) 4 d c 2 N = (( D rho g / sigma )) 5 r c N = (( H / D )) 6 r 2 2 2 N = (( D H rho g / D sigma )) 7 s c c N = (( Delta rho / rho )) 8 c 0.25 0.25 0.75 N = (( mu g / rho sigma )) 9 c c 4 N = (( V rho / g sigma )) 10 d c N = (( g sigma / rho )) 11 c N = (( D / D )) 12 r c 4 0.25 0.25 0.25 N = (( V rho / g sigma )) 13 d c N = (( N D / V )) 14 r c N = (( N D / V )) 15 r d 2 N = (( D g Delta rho / sigma )) 16 r N = (( V D rho / mu )) 17 c r c c N = (( D / H )) 18 c N = (( D / D )) 19 s cThe Sauter droplet size is computed by the high Reynolds formula from Kumar and Hartland (1986):
0.55 -1.3 0.75 -0.3 0.28 (d / D ) = k N exp ((-0.23 N )) N N N N 32 r 1 2 3 4 5 6with 10^3 k=7.01 for no mass transfer. The Kumar and Hartland disperse holdup:
n_1 n_2 n_3 n_4 0.22 0.35 phi = [ k + k N ] N N N N (( 1 + (V / V ) )) d 1 2 2 7 8 9 10 c dwhere (using all data) k_1=65.73, k_2=74.20, n_1=1.24, n_2=-0.34, n_3=-.049, and n_4=0.53. The slip velocity is computed by
0.52 0.25 -0.45 0.08 1.03 0.51 0.28 V = [ k + k exp ((-1.28 N )) ] N N N N N N N s 6 7 2 8 11 9 5 6 12 13with (for all data) 10^2 k_6=-5.11 and k_7=0.20. the continuous phase dispersion coefficient is given by
-0.08 -0.16 0.1 2 (E / V H) = 0.42 + 0.29 (V / V ) + ( 0.0126 N + ( 13.38 / 3.18 + N ) ) N 7 N N N c c d c 14 14 1 12 18 19and the disperse phase coefficients
-0.64 -0.7 -0.9 (E / V H) = 0.3 (( V + V / V )) + 9.37 N N N phi d d c d d 15 16 12 d
1/4 1/4 3/4 0.19 d = C [( mu g rho / Delta rho sigma rho )] sqrt( \sigma \over g \Delta \rho ) p 1 w w dwhere C_1 is 2.54, 2.24, or 3.13 for no mass transfer, transfer from c to d, and from d to c. The dispersed phase holdup is
-1.11 -0.50 2 2 1/3 -0.72 0.10 1/3 1.03 1/3 phi = C e (( Delta rho / rho )) [( (1 / a ) (( rho g / mu )) )] (( mu / mu )) [ Vd (( rho / g mu )) ] exp [ 0.95 V ( rho / g mu ) ] d 2 c p c c d c c c c c cwhere C_2 is 5.34, 6.16, or 3.76. The slip velocity is
-0.11 0.40 2 2 1/3 0.61 -0.10 -1/3 V = C e (( Delta rho / rho )) [( (1 / a ) (( rho g / mu )) )] (( mu / mu )) (( rho / g mu )) slip 3 c p c c d c c cwhere C_3 is 0.24, 0.21, or 0.31. Or, as function of the dispersed phase holdup:
-0.17 0.41 2 2 1/3 0.59 -0.10 -1/3 V = C e (( Delta rho / rho )) [( (1 / a ) (( rho g / mu )) )] (( mu / mu )) (( rho / g mu )) (1-phi ) slip 4 c p c c d c c c dwhere C_4 is 0.30, 0.27, or 0.38. The flooding velocity is:
2 1.54 0.41 2 2 1/3 0.30 0.15 V (1+sqrt(R)) sqrt(a / g) = alpha C e (( Delta rho / rho )) [( (1 / a ) (( rho g / mu )) )] (( mu / sqrt(\Delta \rho \sigma / a ) )) c,f p 1 d p c c c por
1.54 0.41 2 2 1/3 0.30 0.15 V sqrt(a / g) = alpha C e (( Delta rho / rho )) [( (1 / a ) (( rho g / mu )) )] (( mu / sqrt(\Delta \rho \sigma / a ) )) d,f p 1 d p c c c pwhere alpha is 1.0 for continuous phase packing wetting and 1.29 for dispersed phase packing wetting.
SET CAUSEWAY=[setting_1;] [setting_2;] [setting_n;]Seven options are available:
The "CW-" in front of the calculation programs denotes the DOS extender type that is linked to the executables. CauseWay executables include the DOS extender, the Rational DOS-Extender uses a separate file (DOS4GW.EXE). The executables often need at least two MegaByte of Extended memory in order to run. They also require a 386-based system (minimum) to run. The different DOS extender executables can be selected under Options/DOS-Extender. Check this setting when you get error messages stating that executables with the "CW-" or "4G-" are not found. CauseWay is currently our default.
As we support the simulators on other platforms as well, they must also run without the ChemSep interface. That is why we store all information (in- and output) in the problem files, with the .SEP extension. As the simulator must know what problem file to run it checks for the existance of the CHEMSEP.FIL file. If it exists, it will read this name from this file. If the file doesn't exist or no CHEMSEP.FIL file was found in the current directory, the simulator will prompt you for a SEP file. DOS and Windows executables of the simulators also allow specification of the SEP file on the commandline, making the CHEMSEP.FIL unnecessary.
To build in maximum flexibility we use the CS and CSW drivers. Both these drivers can run the simulators upon cooperation with the CS2.EXE interface. The default configuration makes the interface run the CSW driver which in turn runs the necessary simulator. The CSW driver insures that all the output of the simulator gets written into a window, giving the illusion that the simulator is an integral part of the interface. While solving the problem the interface is swapped out of DOS memory, and swapped back in upon termination of the simulator. In this case the CS driver is not doing anything and the interface can be invoked directly as well.
Alternatively, the interface can write a CHEMSEP.FIL file, exit, and the CS driver reads it. It determines which simulator is requested to run and start it. The simulator will read the CHEMSEP.FIL again and read the name of the SEP file to solve. Upon termination the CS driver restarts the ChemSep interface, telling it which SEP file it was solving. The interface reads the SEP file and deletes the CHEMSEP.FIL file. Again, the CS driver also ensures the output to be written into a window on the screen. If the interface is loaded directly, without the driver, this is not possible. However, the interface will still run the simulator by swapping itself out of memory, and back into memory after the simulator returns control. It will call the simulator directly without writing a CHEMSEP.FIL file. Without any driver the interface can't redirect the ouput to a window and thus the screen is cleared before and restored after the simulator runs.
Switching between these modes of solving SEP files is done by specifying the user program under the solve options. The default setting is to call the CSW driver (the interface will locate it, don't specify its path!). If the user program is left empty, the interface will write a CHEMSEP.FIL and quit if the CS driver was loaded, else it will call the simulator directly. Only in this case the screen will be cleared (it also requires the least amount of DOS memory). If the interface is swapped out of memory, it is written to extended memory or disk. If, for whatever reason, the swapfile on disk is deleted the interface will not be able to recover and abort to DOS.
If you are running the simulator programs by your self (that is not abnormal) you might like the small utility MAKEFIL to generate your CHEMSEP.FIL file. Since the simulator programs do not delete the CHEMSEP.FIL file you have only to create it once for each new problem. The MAKEFIL program takes as first commandline argument the SEP-file and as (optional) second argument the scrap file name:
MAKEFILThe CHEMSEP.FIL file will be written in the current directory. If you run under DOS or Windows, you can also specify the SEP-file pathname as commandline parameter to the simulator, and avoid the use of the CHEMSEP.FIL file all together.[Temporary file]
In order to make ChemSep as versatile as possible, we implemented the User Program entry under the "User program" option in the "Solve Options" menu. If you want to run your own program you can enter its full path and name (with the extension!) there and that program will be run no matter what kind of operation is selected. In order for a user program to have access to the SEP-file, the interface provides the user program with its name as the second commandline parameter. In case it also wants to use the temporary file, that is supplied as the third parameter. Running a user program will clear the screen before it starts executing the user program.
It is logical to suppose that the user just might want to run several programs or his own program(s) before/after calling ChemSep's original calulation programs. To allow this you can use a batch file as the user program (use complete path, name and extension !). Look up in your DOS manual how to make batch files. In order to run the original calculations program we provide it to the batch file as the first commandline parameter and the SEP-file as the second. Even running the user program, the Interface generates a CHEMSEP.FIL file to be read by the calculation program. Here's an example of such a batch file that will type the problem SEP-file first before running (note how the parameters are accessed with %1 and %2):
@Echo off Rem ----------------------------------------------------------- Rem Echo is set off to avoid to show this batch file is running Rem use the "@" to suppress echoing of commands to the screen Rem ----------------------------------------------------------- Rem Type the SEP-file: Rem ------------------ Type %2 Rem Pause for the user to strike a key: Rem ----------------------------------- Pause Rem Run the appropriate ChemSep calculation: Rem --------------------------------------- %1 Rem Done!Of course you can leave out all the Rem(arks) if you type this batch file but it is a good habit to document your code ! You might create a complete set of batch files for different problems. Programs that might be run after the calculation might be cost estimation or design (in case of equilibrium simulation) programs. Using command line parameters (or the CHEMSEP.FIL file) your program might append its results to the SEP-file, so all results will be collected in there. Here is an example of a batch file that will automatically run ChemProp to generate physical property information in the sep-file:
@Echo off echo Running ChemSep with physical property information generation rem Run executable %1 rem Run ChemProp to generate physical property information c:\chemsep\bin\cp /c %2 rem Done!The ultimate freedom is allowed by typing "DODOS" as User Program in the Interface. The driver will automatically locate DOS and run it. All Dos commands will be available to you. The only way to access the SEP-file and other information is to read the CHEMSEP.FIL file. When you are ready to go back to the Interface you type "EXIT" and press Enter.
By allowing you to shell to DOS or run batch files from within ChemSep we have created the maximum flexibility. Although ChemSep takes a lot of effort to prevent your system or ChemSep from crashing, it is possible to do so, using batch files or the DOS shell. Avoid deleting crucial files (such as executable files), or changing system parameters while the ChemSep-Driver is loaded. Do not load any TSR (Terminate and Stay Resident) programs or device drivers, since these programs will be removed from memory after the Driver takes over again. However, the interrupt to trigger these programs usually remains active. Changing directories is allowed, but remember the current CHEMSEP.FIL is written in the current directory ! The driver will change the directory back to what it started running from, when returning to the Interface. If you want to change the current directory use the Directory option in the File menu ! Do not load ChemSep again by invoking the driver, as multiple copies will be loaded into memory.
The Interaction Parameter Data (IPD) files are ASCII text files with the interaction parameters for activity coefficient models and equations of state. Currently we have the following IPD files:
[IPD] Comment=DECHEMA NRTL data @ 1atm. Units=cal/mol # 1101 1921 -189.0469 792.8020 0.2999 Methanol/Water p61 1/1aLines starting with a "#" are comment lines which may appear anywhere. Since the interface will only start reading the file from the [IPD] keyword on, you can start the file with some text describing where the data was obtained and remarks on who/when/how changed the file. Most of the IPD files contain information from the DECHEMA series, a very extensive collection of interaction parameters. The Hayden O'Connell virial parameters are from Prausnitz et al. (1980).
Polynomial K-value and enthalpy correlation coefficients as well as extended Antoine coefficients are stored as component LIBraries (LIB files), which are ASCII files as well. The default LIB files are
[LIB] Comment=Extended Antoine Prausnitz et al. # # ID A B C D E F G 902 3.15799e+01 -3.2848e+2 0. 0. -2.5980e+0 0. 2.0 HydrogenGroup Component Data (GCD) files are binary files containing data for group contribution method such as UNIFAC or ASOG. Currently only UNIFAC files are present (UNIFACRQ.GCD, UNIFACVL.GCD and UNIFACLL.GCD) for Vapour-Liquid and Liquid-Liquid systems. These files need not be changed, unless for example the UNIFAC group tables have changed. Several GCD-files are used for the estimation of pure component data in ChemLib and UserPcd.
Internal layout data (ILD) files are text files which store tray or packing layouts for use by the nonequilibrium model. This way a specific design can be saved and reloaded upon demand. Two formats are used, which resemble eachother very much. The difference lies in the use of an equal sign in the internals type. Without the equal sign we are dealing with one of the fixed internal layouts. If an equal sign is present then the two following numbers define the type of operation (VL=1 or LL=2) and the number of parameters that describe the internalstype. For example the variable internals type for a baffle tray column section is:
[Section] 10 Discrete Baffle = 1 3 1.2 Column diameter (m) 0.5 Baffle area (%) 0.7 Baffle spacing (m)Note that the parameters are -always- written in base units, not in the units that are used for display/input. See also the discussion of the model and internals definitions in ChemSep.
Packing Data (CHEMSEP.*PD) files contain many physical and model parameters for various random and structured packings. They are text files that might be modified by the user with an ASCII editor, though there is one restriction, namely that the first line should not be changed! A shortened version of the structured packing data file CHEMSEP.SPD is shown below.
# CHEMSEP SPD Structured Packing Data # # Type (Name): Specific Equiv. Channel Packing Void ... # packing diam. flow factor: fract: ... # surface: angle: !---------------------- Koch Flexipac 1 M 558/m 0.00897m 45ø 98/m 0.91 ... Koch Flexipac 2 SS 223/m 0.01796m 45ø 43/m 0.95 ... Koch Flexipac 3 M 135/m 0.03592m 45ø 26/m 0.96 ... Koch Flexipac 4 M 69/m 0.07183m 45ø 20/m 0.98 ... @---------------------- Glitsch Gempak 1A M 131/m 0.03592m 45ø 30/m * ... Glitsch Gempak 2A SS 223/m 0.01796m 45ø 52/m 0.95 ... Glitsch Gempak 2AT SS 223/m 0.01796m 45ø * 0.96 ... Glitsch Gempak 3A M 394/m 0.01346m 45ø 69/m * ... Glitsch Gempak 4A M 525/m 0.00897m 45ø 105/m * ...Lines starting with "#" are comment lines and are ignored (except for the first line). A line starting with "!" is used to set the length of the packing type identifiers which is set equal to the length of that line. Lines starting with an "@" will insert a separator in the list with packings and blank lines will be ignored. As you can see units may be added as long as there is no space between the number and the unitstring (otherwise errors will occur in reading this file!).
Miscellaneous ChemSep files include:
Number (in Reference units) = fm * ( Number (in Units) - fo )For example 22 C = 1.0 ( 22 - -273.15) = 295.15 K. You can inspect ASCII text files with ChemSep's file viewer ( F7) or make simple modifications with edit-file. However, we strongly suggest you do not change the original data files that come with ChemSep. We carefully selected and typed the data into these files and other users might use your changed data and obtain erroneous results. We urge you first to copy the file to another name before you change or add anything in these files. Errors in the unit conversion can be very frustrating so it is good to check some results if you have changed or added a unit definition. To encourage you to do so we made CHEMSEP.UDF a read-only file !
The CHEMSEP.SYN is a file containing synonyms for over 1000 compounds. While searching for a special component name you can use synonyms if you have selected to do so in the options interface spreadsheet. You must select this file as your synonyms file. The synonym search does not work while typing in a searchlist for a synonyms name. You will have to issue a search under the synonyms name ! The synonyms file is an ASCII file you can modify to your needs.
The CHEMSEP.HLP file contains the information to provide you with help when you press ( F1) for help. It is a binary file that can not be changed. It is generated with the MAKEHELP utility from help source (HSR) files. This utility and the source files are not part of the ChemSep distribution. If you find errors or shortcomings in the help please notify the authors.
The CHEMSEP.TXT file contains the ASCII text of this book. Normally ChemSep will be configured with function key ( F9) assigned to automatically load this file. There are also other formats available on our ftp-site (see Author information). The ASCII version is somewhat limited in available characters, sub- or superscripts, and equations.
The CHEMSEP.CNF file is used to store the default configuration, such as the macro definitions, directory structure, selected video and printer devices, and solve options. CHEMSEP.CNF is the default configuration file which will be loaded upon startup of the interface. If none is found in the current directory, the file supplied in the original distribution is used. This allows one to have multiple CHEMSEP.CNF files in different directories, which automatically configure the interface to (a) specific problem(s).
Finally, the CHEMSEP.SCR file contains the additional introduction
screen(s) that are shown on startup of the program. These are used to
stipulate conditions of the use of the program, but could can be
adapted to suit the users needs (for example if
ChemSep is installed on a network, the operator can place
important notes here). If the file contains more lines than can be
shown on the screen, the user will have to press
Remember that normally the Interface will not read any comments you have added yourself to the SEP-file and thus, not save them again when saving from within the Interface. To overcome this problem we have added the [User-Data] and [End User-Data] delimiters. When loading a SEP-file it is the last section that is looked for. If found, it is read into a buffer and written back to the file if saved again. Note that the delimiters each have to be on a separate line. You can use this data block to save information about the problem or to store parameters for your own programs that process SEP-files. You can edit User Data within the interface of ChemSep under the Solve Options ( F6). The following keywords are used to switch several hidden features (with explanation in "()"):
[D-Models] Diffusivity model (0=Maxwell-Stefan, 1=Effective) Liquid MS-diffusivity model (0=Kooijman-Taylor, 1=Wesselingh-Krishna) [No user interaction] (if present the user is not asked for more iterations if maximum number has been reached, but the program exits) [Sensitivity] (sensitivity factors) Vapour/Light-liquid Mass Transfer Coefficient Liquid/Heavy-liquid Mass Transfer Coefficient Interfacial AreaFor a column problem the INPUT subsections are:
[CHEMSEP] Version number and SEP file name [Paths] Current directories [Units] Current set of units [Components] Number of components, for each component library offset (in the PCD file), Index, Name, and PCD-library filename [Operation] Operation type and kind, condenser and reboiler types, number of stages, feeds, sidestreams, and pumparounds. [Properties] Property selections: [Thermodynamics] K model, Activity coefficient, Wilson model, UNIQUAC model, Equation of State, Cubic EOS, Virial EOS, Vapour pressure models [Enthalpy] Enthalpy model [Physical Properties] Physical Properties model selections [Property Data] Pure component data or interaction parameters if required [Specifications] Specifications: [Heaters/Coolers] Number and stage with duty, if specified. [Sections] Number and for each section: section number, begin and end stages, model selections, and tray/packing layout data. [Efficiencies] default efficiency, number of exceptions and stage with value, if specified. [Pressures] Type of pressure specification, condenser, top, bottom pressures, pressure drop. [Feeds] Number and for each feed: feed state, stage, temperature, pressure, vapour fraction, number of componentflows, molar component flows. [SideStreams] Number and for each sidestream: stage, phase, specification type and value. [Condenser] Specification type and value(s) if present [Reboiler] Specification type and value(s) if present [Solve options] Initialization type, solving method, damping factor (if present), accuracy, maximum number of iterations, and print options. [Programs] Temporary file, and user program [User-Data] Here the user data text is written. [End User-Data] [End of Input]For an equilibrium column the [Sections] subsection will be absent, for a nonequilibrium model the [Efficiencies] subsection is not needed. For a flash a different set of specifications is present consisting of the [Feeds] subsection and a [Flash] specifcation subsection where flash type and specifactions are made. The RESULTS section for a nonequilibrium problem looks like:
[Results] [Profiles] [Temperatures] [Vapour phase compositions] [Liquid phase compositions] [Interface vapour mole fractions] [Interface liquid mole fractions] [Murphree efficiencies] [Mass transfer rates] [Condenser Heat Duty] [Reboiler Heat Duty] [K-values] [Feed streams] [Top product] [Bottom product] [Sidestreams] [Designed Sections] [Operating Limits] [End of Results]For equilibrium problems the Interface mole fractions, Murphree efficiencies, mass transfer rates, designed sections, and operating limits will be missing from the above list.
ChemSep supports the following printers and plotters and desktop publishing file formats:
Dot matrix printers: | Max.Mode | Desktop Publishing: | Max.Mode |
Epson 9-pin dot matrix | 8 | Zsoft PCX | 1 |
Color Epson 9-pin dot matrix | 5 | Windows 3 BMP | 1 |
Epson 24-pin dot matrix | 8 | Gem IMG | 2 |
Color Epson 24-pin dot matrix | 5 | TIFF compressed | 2 |
IBM Proprinter X24 | 8 | TIFF uncompressed | 2 |
IBM Quietwriter | 8 | ANSI CGM | 1 |
Toshiba 24-pin dot matrix | 2 | AutoCad DXF | 0 |
OkiData ML-92 dot matrix | 2 | Video Show | 0 |
Word Perfect Graphic | 1 | ||
Laser/Inktjet printers: | Max.Mode | Hewlett-Packard plotters: | Max.Mode |
LaserJet II | 8 | HP 7090 | 3 |
LaserJet III | 8 | HP 7470 | 1 |
DeskJet | 8 | HP 7475 | 7 |
Color DeskJet | 8 | HP 7550 | 7 |
PaintJet | 14 | HP 7585 | 9 |
Postscript | 11 | HP 7595 | 9 |
The definitions file must start with "[ChemSep Definitions]" followed by a line which defines the version (the current version is 3). Only after these two lines content may appear. Everything after a "mbox char35" is regarded as comment. Comments may appear on any content line, and may also be appended after any content. There are five different definitions types that are read from the definition file: [InternalType], [Operation], [ModelType], [Internal], and [Model]. Other entries are ignored. The five different types of definitions may be mixed throughout the definitions file. A definition may have no blank lines in it as entries are ended by a blank line!
The [InternalType] defines the nature of an internal, that is whether it is continuous or discrete. The [Operation] defines the nature of the operation, that is whether we are contacting a vapor and a liquid (VL) or two liquids (LL). The [ModelType] defines the type of models, that is whether the model describes mass transfer coefficients, pressure drops, etc. Each of these has the following fields: ID, Name, and Short, for example:
[Operation] ID=1 Name=Vapor-Liquid Short=VLThe [Internal] definitions link these types together: its "Type" must be one of the defined internal types, and its "Operation" one of the defined operation types. It also defines which "Models" need to be specified for the internal as well as which layout "Parameters" the internal has. These parameters and models are separated by semicolons (;) or comma's (,). The parameters might include a units string in parenthesis ("(",")") and a default value (with different units) after a colon (:). Note that multiple "Parameters" entries are allowed but that that its total lentgth is limited to 255 characters (these rules also apply for the "Models" and "Internals" entries).
A parameter entries may also contain a picklist in brackets ("[","]"), separated by vertical lines. Then instead of a numeric entry the interface will ask the user to select the parameter from a picklist of items. The first item in the list will result in the value 1, the second in 2, etc. A default item may also be specified using the colon and the numeric value of the parameter. For the model entries, either ID numbers or short notation may be used, as long as they are defined ModelTypes. For example the definition for a baffle tray may be:
[Internal] ID=10 Name=Baffle tray Type=Discrete Operation=VL Short=Baffle Models=MTC,PD,VF,LF,DSGN # multiple parameter lines: Parameters=Column diameter (m):120cm;Baffle area (%):50;Baffle spacing (m) Parameters=Baffle type[Disk & donut|Single flow|Multiflow]:2 #default=SingleThe [Model] definitions has instead of a "Models" entry the "Internals" entry, where the types of internals for which the model can be used are defined. It also has the "Parameters" where model parameters can be defined. For the definition of a baffle tray above the parameters are a column diameter (with meters as units and a default value of 1.2m), the baffle area (as percentage, default 50%), the baffle spacing (in meters, without any default), and the baffle type (with a pick list between 3 different types where the default is the single flow).
The models that are required for the baffle tray are a mass transfer model (MTC), a pressure drop model (PD), models for the vapor and liquid flow (VF and LF), as well as a design model (DSGN). Note that the design method is also considered as a model. A sample design method definition is:
[Model] ID=1 Name=Fraction of flood Type=DSGN Operation=VL Short=%flood Internals=Bubble cap;Sieve;Valve Parameters=Fraction of flooding:75%;Fraction of weeping:70%Model parameters can be loaded from libraries (*.PAR) where for each internal different parameter values are stored. Note that a model defintion may also be repeated (with the same id, name, type, and operation) so to add another internalstype to the list to which the model applies. Short fields are optional, and have a maximum length of ten characters, used for displaying selected models etc. ID fields associate a unique number to the definition. Only for the internal and model definitions non-unique numbers are allowed. When the interface reads the definitions file it uses the Short descriptions to assign the ID's for the Operations, InternalTypes and ModelTypes, see Table 7.1 (the Short descriptions used by the interface are in parenthesis). Thus, you will have to use these Short descriptions but are free to change the ID numbers or names.
Table 7.1. Interface Assigned Types
Operation | InternalType | ModelType |
Vapor-Liquid (VL) | Discrete (DIT) | Mass transfer coefficient (MTC) |
Liquid-Liquid (LL) | Continuous (CIT) | Pressure drop (PD) |
Vapor flow (VF) | ||
Liquid flow (LF) | ||
Entrainment (ENTR) | ||
Holdup (HOLD) | ||
Light liquid flow (LLF) | ||
Heavy liquid flow (HLF) | ||
Backmixing (BACK) | ||
Design method (DSGN) | ||
[Internal] ID=1 Name=Bubble cap tray Type=Discrete Short=Bubble cap Operation=VL Models=MTC,PD,VF,LF,ENTRThe builtin types are listed in Table 7.2. This table also lists defined vapor and liquid flow models and models for entrainment.
Table 7.2. Assigned Internals and Vapor flow, Liquid flow, and Entrainment Models
Internal | Vapor Flow / | Liquid Flow / | Entrainment / |
Light liquid flow | Heavy liquid flow | Backmixing | |
Bubble cap tray (1) | Mixed (1) | Mixed (1) | None (1) |
Sieve tray (2) | Plug flow (2) | Plug flow (2) | Estimated (2) |
Valve tray (3) | |||
Dumped packing (4) | |||
Structured packing (5) | |||
Equilibrium stage (6) | |||
RDC compartment (7) | |||
Spray column stage (8) | |||
[Model] ID=1 Name=AIChE Type=MTC Operation=VL Internals=Bubble cap,Sieve tray,Valve trayAssigned models for Mass Transfer Coefficient and Pressure Drop models are listed in Table 7.3. Model parameters may be read from parameter libraries (*.PAR) that have a format like the packing data files. The first line of such a ASCII text file must start with "# CHEMSEP xxxx" where xxxx must be replaced by the full name of the model (as defined in CHEMSEP.DEF). Upon choosing the library option in the interface the library will be automatically pre-selected if the name of the library file is the same as the short name of its model.
Table 7.3. Assigned MTC and PD Models
Mass Transfer Coefficient | Pressure Drop |
AIChE (1) | Fixed (1) |
Chan Fair (2) | Estimated (2) |
Zuiderweg (3) | Ludwig 1979 (3) |
Hughmark (4) | Leva GPDC (4) |
Harris (5) | Billet-Schultes 1992 (5) |
Onda et al. 1968 (6) | Bravo-Rocha-Fair 1986 (6) |
Bravo-Fair 1982 (7) | Stichlmair-Bravo-Fair 1989 (7) |
Bravo-Rocha-Fair 1985 (8) | Bravo-Rocha-Fair 1992 (8) |
Bubble-Jet (9) | |
Bravo-Rocha-Fair 1992 (10) | |
Billet-Schultes 1992 (11) | |
Nawrocki et al. 1991 (12) | |
Chen-Chuang (13) | |
Handlos-Baron-Treybal (14) | |
Seibert-Fair (15) | |
Kronig-Brink-Rowe (16) | |
Rose-Kintner-Garner-Tayeban (17) | |
Sherwood (20) | |
Regular mail: Ross Taylor / Harry Kooijman Department of Chemical Engineering Clarkson University Potsdam, NY 13699, USA Telephone: Ross Taylor (315) 268 6652 Harry Kooijman (908) 771 6544 Electronic mail: taylor@clarkson.edu kooijman@clarkson.edu WWW page: http://members.home.net/hkooijman/chemsep.htm http://www.clarkson.edu/~chengweb/faculty/ taylor/chemsep/chemsep.htmlWe typed most of the code with the Multi-Edit text editor. It allows us to switch between the source code of the drivers, interfaces, and column simulation executables, as well as the text files for the help and documentation. Each different type of file has its own commands associated with it (compile source code, run latex on the documentation, etc.). As the project is over several hundred thousands lines of source and text, the editor has proven to be very valuable to us.
The drivers and interfaces are written in Turbo Pascal (version 7.0), the column and flash programs are written in standard Fortran 77, which we compile with WATCOM Fortran and link with the CauseWay DOS extender. The source code for the simulation executables has also been successfully compiled and executed on a range of operating systems and platforms .
2 C H OH -> (C H ) O + H O 2 5 2 5 2 2In the flowsheet input file (ep.fs) you define the units and streams in your flowsheet. The input file consists of four parts. The first part consists of one section where all the components, units, streams, feeds, and stream estimates are declared as well as the output file, executable directory, and the method and accuracy used in solving the flowsheet.
[Flowsheet] Comment=Ethylether Production Components=Water,Ethanol,Diethylether Units=Mixer,Reactor,Sep-1,Sep-2 Streams=1,2,3,4,5,6,7,8 Feeds=1 Estimates=7 Output=ep.fs ExeDir=c:\cs\exes\ MaxIter=20 Method=Direct Accuracy=0.001 ! this is a commentThe section starts with the [Flowsheet] identifier. The next lines are part of this section until an empty line is found. You can add (non empty) comment lines in the input file by, for example, using any punctuation character to start the line. The flowsheet program looks for the specific keywords and if none is found it regards the line as a comment. However, comment lines do not get copied to the output file. In each section the keywords can be different. The keywords are the identifiers left of the equal signs in the above flowsheet section. They are not case-sensitive and you may enter them in any order (some restrictions do apply, though).
Comma's are used to separate the various components, streams, units, feeds, or stream estimates (therefore, they can not be part of a name. This is important as IUPAC component names sometimes use comma's. In that case omit the comma's in the component name. Currently the component names are not checked and it is assumed that the same components in the same order are used in all the units! This is very important. The case of the names you specify is not important but the one used in the flowsheet section is used in the rest of the report.
In our example we have the following streams we have given our streams numbers from 1 to 8 (instead of numbers names may also be used). We have only one feed, 1, and we will make a stream estimate for the recycle ( 7). Stream 2 is a stream that is not used (we need it since the flash program that simulates the mixer produces a vapor product stream which we have to assign). Our unit operations are a mixer, reactor, and two columns ( Sep-1 and Sep-2) where we separate the products. Our output file is ep.fs (the same name as our input file which will be overwritten) and the flowsheet method is Direct substitution with a convergence criterium of 0.001. The maximum number of iterations is specified here as 20.
The next part describes the units and, therefore, can consist of multiple sections. We have four units in our example. Each unit section starts with the [Unit] identifier and has five keywords: Name, File, Type, Inlets, and Outlets.
[Unit] Name=Mixer File=ep-mix.sep Type=4g-flash.exe Inlets=1,7 Outlets=2,3 [Unit] Name=Reactor File=ep-reac.sep Type=reactors.exe Inlets=3 Outlets=4 [Unit] Name=Sep-1 File=ep-sep1.sep Type=4g-neq2.exe Inlets=4 Outlets=5,6 [Unit] Name=Sep-2 File=ep-sep2.sep Type=4g-neq2.exe Inlets=6 Outlets=7,8The specified unit name must match the one given in the flowsheet section. The unit file is the associated data file, for most units it is the Sep-file which was generated by ChemSep. The unit type is the program that must be run to simulate this specific unit. For a ChemSep flash that is 4g-flash.exe, for an equilibrium column 4g-col2.exe, and for a nonequilibrium column 4g-neq2.exe. This allows you to make your own unit simulation program that reads a datafile with feed and product section as in a Sep-file. In the case of the reactor this is done by the reactors.exe program (which is described below).
The units are connected by in- and outlets streams, which were declared in the flowsheet section above. Be sure to use the same names (case is not important) as otherwise the flowsheet will be incorrect or incomplete. If you just want to analyze a flowsheet or not to simulate the unit it, ommit the unit file and type (they will not be executed).
[Feed] Name=1 Temperature=40C Pressure=1.01325bar Rate=20mol/s Z=0.15,0.85,0The next part consists of the [Feed] sections where the feed streams are defined. The specified feed name must match the one given in the flowsheet section (again, case is not important). Furthermore, the feed stream pressure, temperature, flowrate, and composition must be specified. Append units when specifying temperatures, pressures, and flows. The default units (which can be omitted) are temperatures in degrees Celcius, pressure in bar (absolute), and flows in mol per second. The flowrate and compositions can also be set by specifying the component flows, for example for the definition of the feed above we could use as well
[Feed] Name=Feed Temperature=40C Pressure=1atm F=3,17where F is the list of component flows (default units mol/s). If all the component flows are specified, the total flowrate and compositions are computed. A partial list (like Z=*,0.85) can be specified as well, which is useful for supplying stream estimates.
All streams are reset at the start, and when the inlets are written to the unit-file only the specifications supplied in the input file are written to the unit file. If stream values are reset the values already present in the unit-file will be used.
The last part consists of [Estimate] sections where estimates of stream variables (flow, temperature, pressure, or composition) can defined. This part is optional and only required if you specified streams under the estimates keyword in the flowsheet section. Stream estimates have the same input format as feeds have, except for the different identifier, of course. In our case we estimate the recycle stream 7 to start with a better value of the flowrate to the reactor:
[Estimate] Name=7 Temperature=80C Pressure=1.01325bar Rate=10mol/s Z=0.85,0.15,0The files that make up our example all start with EP (ep.fs, ep-mix.sep, ep-reac.sep, ep-sep1.sep, ep-sep2.sep) and are included in the ChemSep distribution (in the fs directory).
A unit evaluation takes the unit inlet streams and puts them (in the same order!) into the file under the [Feeds] section. Then it runs the associated program and reads the [Top product] section as the first outlet stream, the [Bottom product] section as the second outlet stream, and each following [Sidestream section as the third, fourth, etc. outlet streams.
In the case that the flowsheet analysis finds one or more cycles it will inform you and set the recycle flag to start an iterative run. The criterium is the maximum relative (or absolute) differences in stream variables when any stream gets updated. For each iteration this is initially set to zero and computed over the simulation of all the units. Convergence is obtained if the criterium is below the specified accuracy (default is 10^-2) or if the maximum number of iterations (default is 20) is attained (In the case that an absolute criterium is used the temperature difference is divided by 10 and for the pressure by 10^4 to scale the differences).
Sometimes it may be handy to abort the simulation or to stop it and inspect certain streams during the simulation. This can be accomplished by holding the Shift keys or toggleing Caps Lock on. The program will beep and display the next unit to be simulated, the current attained convergence, and the current iteration number. A simple menu allows you to change the maximum number of iterations or the accuracy, display a stream, continue, or to quit the simulation. It will also allow you to swap to DOS to do other work, like loading Sep-files into the ChemSep interface to make changes or to evaluate the intermeadiate results. Typing "exit" will then return you to the flowsheet program (don't forget this).
Once the flowsheet is convergenced the output file is written, consisting of four parts: the input file (generated from the information read in), a flowsheet analysis, the mass balances, and a report of all the streams.
Incidence Matrix ---------------- Unit\Flow: 1 2 3 4 5 6 7 8 Mixer + - - + Reactor + - Sep-1 + - - Sep-2 + - - Adjacency Matrix ---------------- Unit: Mixer Reactor Sep-1 Sep-2 Mixer 3 Reactor 4 Sep-1 6 Sep-2 7 Distance Matrix --------------- Unit: Mixer Reactor Sep-1 Sep-2 Mixer 4 1 2 3 Reactor 3 4 1 2 Sep-1 2 3 4 1 Sep-2 1 2 3 4 Cycle Matrix ------------ Cycle: Rank: Streams: C1 4 3,4,6,7 Stream: Frequency: 3 1 4 1 6 1 7 1 Cycle: Rank: Units: C1 4 Mixer,Reactor,Sep-1,Sep-2 Evaluation order from analysis = Reactor Mixer,Sep-1,Sep-2These matrices can be useful in assessing the structure of the flowsheet and the flowsheet evaluation order of the units. The flowsheet analysis will also give an ordering of the units that might be better than the specified order.
Flowsheet is solved. Attained convergence = 0.000826326 Number of iterations = 10 Used relative differences.
Balances in mol/s Total Molar Balances: Mixer = 1+7-2-3 = 20+19.7097-0-39.694 = 0.0157021 Reactor = 3-4 = 39.694-39.694 = 0 Sep-1 = 4-5-6 = 39.694-8.45862-31.2354 = -0.0000204891 Sep-2 = 6-7-8 = 31.2354-19.7097-11.5257 = 0 Overall = 1-2-5-8 = 20-0-8.45862-11.5257 = 0.0156797 Water Balances: Mixer = 3+2.799822-0-5.797864 = 0.00195764 Reactor = 5.797864-14.21037 = -8.412508 Sep-1 = 14.21037-0.0000697587-14.21033 = -0.0000232831 Sep-2 = 14.21033-2.799822-11.41044 = 0.0000614673 Ethanol Balances: Mixer = 17+16.66376-0-33.65003 = 0.0137314 Reactor = 33.65003-16.82501 = 16.82502 Sep-1 = 16.82501-0.0460432-16.77897 = 0 Sep-2 = 16.77897-16.66376-0.115257 = -0.0000505315 Diethylether Balances: Mixer = 0+0.246113-0-0.246111 = 1.920853E-06 Reactor = 0.246111-8.658619 = -8.412507 Sep-1 = 8.658619-8.412504-0.246113 = 1.949957E-06 Sep-2 = 0.246113-0.246113-3.203753E-08 = -3.203753E-08Care must be taken that the individual component balances are also satisfied, they are written after the total molar balances. Also, a reactor total molar balance will not be zero if the reaction changes the number of moles in the mixture. A reactor component balance will not be zero if that component is involved in one of the reactions that has a nonzero conversion. In the case of our example we see that the reaction rate is 8.4 mol/s dietheylether.
Stream 1 Temperature (C) 40 Pressure (bar) 1.01325 Flowrate (mol/s) 20 Zwater 0.15 Zethanol 0.85 Zdiethylether 0 Stream 2 is zero Stream 3 Temperature (C) 58.68402 Pressure (bar) 1.01325 Flowrate (mol/s) 39.694 Zwater 0.146064 Zethanol 0.847736 Zdiethylether 0.00620021 Stream 4 Temperature (C) 49.85001 Pressure (bar) 1.01325 Flowrate (mol/s) 39.694 Zwater 0.357998 Zethanol 0.423868 Zdiethylether 0.218134 Stream 5 Temperature (C) 34.76102 Pressure (bar) 1.01325 Flowrate (mol/s) 8.45862 Zwater 8.24705E-06 Zethanol 0.00544335 Zdiethylether 0.994548 Stream 6 Temperature (C) 77.918 Pressure (bar) 1.01325 Flowrate (mol/s) 31.2354 Zwater 0.454943 Zethanol 0.537178 Zdiethylether 0.0078793 Stream 7 Temperature (C) 76.11902 Pressure (bar) 1.01325 Flowrate (mol/s) 19.7097 Zwater 0.142053 Zethanol 0.84546 Zdiethylether 0.0124869 Stream 8 Temperature (C) 96.62302 Pressure (bar) 1.01325 Flowrate (mol/s) 11.5257 Zwater 0.99 Zethanol 0.01 Zdiethylether 2.77966E-09
\x = Skip execution \r = Force iteration \utF = Set default temperature units to degrees Fahrenheit \utpsia = Set default pressure units to psia \utkmol/s = Set default flow units to kmol/s \da = Use absolute difference criterium \dr = Use relative difference criterium (default)Although the streams in the input file can be defined by any units, the output file will be written with the default units (including the part with the input file!). Unit definitions are read from the FS.UDF file, or if that file is not found, from the CHEMSEP.UDF file.
Production of benzene by hydrogenation of toluene ------------------------------------------------- Reactions: C5H5CH3 + H2 -> CH4 + C6H6 2 C5H5CH3 + H2 -> 2 CH4 + (C5H5)2 Components: Toluene, Hydrogen, Methane, Benzene, Diphenyl [Reactor] Benzene-Reactor 400 K outlet temperature 10000 Pa pressure drop 5 components 2 reactions 1 base component r1 0.9 conversion r1 -1 -1 1 1 0 stoichiometry coef. r1 1 base component r2 0.2 conversion r2 -2 -1 2 0 1 stoichiometry coef. r2 [Feeds] 1 number * * 400 K temperature 101325 Pa pressure * 5 components 0.1 kmol/s toluene 0.1 kmol/s hydrogen 0 kmol/s methane 0 kmol/s benzene 0 kmol/s diphenylAfter running the reactors program the output is written to the same file and contains a small [Results] section which reports whether the specified convergence(s) was attained followed by a [Top product] section with the reactors outlet stream.
[Results] Conversion(s) on base components was limited by component 1 to 90% of specified conversion(s)If there are multiple reactions and one (or more) of the component feedflows is constraining the reactions, the conversions are all equally affected. For the ether reactor of our example the reactors input file looks like:
[Reactor] Reactor 323 K outlet temperature 0 Pa pressure drop 3 components 1 reactions 2 base component r1 0.5 conversion r1 1 -2 1 stoichiometry coef. r1 [Feeds] 1 Number 1 Feed state T & P * 331.834 Temperature 101325 Pressure * Vapour fraction 3 componentflows 0.00579786 Component 1 flow 0.03365 Component 2 flow 0.000246111 Component 3 flow
[Make-up] Make-up Unit .06 Total flowrate (always in kmol/s) [Feeds] 2 Number 1 Feed state T & P 1 stage 373.15 Temperature 100000 Pressure * Vapour fraction 3 componentflows 0 Component 1 flow 0 Component 2 flow 0.001 Component 3 flow 1 Feed state T & P 1 stage 373.15 Temperature 101325 Pressure * Vapour fraction 3 componentflows 5.041514E-09 Component 1 flow 0.0000207905 Component 2 flow 0.0598143 Component 3 flowWhen the make-up program runs it appends a small [Results] section displays the total flowrate of the make-up feed followed by a [Top product] section containing the resulting stream.
[Results] 0.0001649 = Make-up flowSince the make-up feed has to be specified at the beginning this could cause the mass balances to be incorrect. Therefore, if the flowsheet encounters a make-up unit, and the first inlet stream is declared as a feed, its flowrate is adapted on writing the output.
[Splitter] Recycle purger 0.1This is followed by the regular [Feeds] section (see for example, the reactor description), where only one feed may be specified. After running the splits program top and bottom product streams will be appended to the input file, where the top product stream has a flowrate equal to the feedflowrate multiplied with the splitfactor, and the bottom product contains the rest of the feedflow.
Figure 8.3. Distillation with Decanter
Figure 8.4. Azeotropic Distillation at Two Pressures
Figure 8.5. Petyluk Columns
Table 8.1. UNIQUAC LLE interaction parameters (K)
Components i-j | A_ij | A_ji |
$25^o C$: | ||
pentane-benzene | 295.38 | -171.63 |
pentane-sulfolane | 532.04 | 187.84 |
benzene-sulfolane | 151.02 | -36.08 |
$50^o C$: | ||
pentane-benzene | 179.86 | -95.70 |
pentane-sulfolane | 375.93 | 247.77 |
benzene-sulfolane | 131.51 | -6.28 |