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 "()"):
ChemSep supports the following printers and plotters
and desktop publishing file formats:
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:
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:
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:
Table 7.1. Interface Assigned Types
Table 7.2. Assigned Internals and Vapor flow, Liquid flow, and
Entrainment Models
Table 7.3. Assigned MTC and PD Models
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 .
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.
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).
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:
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.
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)
7.6 The SEP-file format
The SEP-files are written in a format which is readible by a human as
well as by the calculation program. However, there are some strict
rules ! The SEP-file is constructed using delimiters in square
brackets: []. The order of the delimeters is not directly of
importance (although ordering delimiters enhances the speed of reading
a SEP-file). The sections under the delimiters are ordered and have a
fixed format. Usually they consist of lines with a selected value and
a comment. ChemSep uses two major sections: the INPUT
and RESULTS sections. The INPUT starts with the delimiter [ChemSep]
and ends with [End of Input]. The Results starts with [Results] and
ends with [End of Results]. Within these sections there are sub
sections. Note that a "*" denotes that the value is not yet set by the
user. In some cases the Interface might create cryptic "* *" lines
where the first star denotes the (not yet known) value of a selector
and the second the description of the (not yet) selected item.
[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 Area
For 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.
7.7 Printing graphs in ChemSep
ChemSep supports dot matrix printers, laser printers,
inktjet printers, or HP plotters to print/plot its graphs. Besides
printing directly to a printer, ChemSep can also write
the output to a file. ChemSep supports nine DeskTop
Publishing (DTP) file formats as well. In order to print to a device
or write to an output file ChemSep needs to know what
type of device you have. Select device, mode, port (file) and work
path (for temporary files written while generating output) in the
printer setup under the graphs or the output setup in the options. Of
course each printer has usually several different modes to print. By
default none of the graphs are printed in color, unless a color device
is selected. ChemSep lets you choose between three
page formats:
However, these page formats vary slightly from printer to printer.
Plotters usually plot only at FULL size. Besides these three page
formats there are (usually) additional modes. Select "Other" and type
the mode number you want (see table below). The lowest mode number is
zero and will always work.
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 7.8 Model Definition and Selection
ChemSep reads a definitions file (CHEMSEP.DEF) at
startup, where models for the mass transfer coefficients, pressure
drop, flow models, entrainment, and holdup are defined. This
alleviates us from adapting the ChemSep interface upon
any addition or modification of a model. The defintion file also
includes the definitions of any non-standard type of internal. In case
no definitions file is found, the nonequilibrium part of
ChemSep is disabled.
[Operation]
ID=1
Name=Vapor-Liquid
Short=VL
The [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).
[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=Single
The [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).
[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.
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,ENTR
The builtin types are listed in Table 7.2. This table also lists
defined vapor and liquid flow models and models for entrainment.
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 tray
Assigned 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.
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) 7.9 Author and program information
The authors can be reached through
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.html
We 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.
Acknowledgements
WordPerfect is a product of the WordPerfect Corporation. Microsoft
Word is a product of Microsoft Incorporated. DOS4GW is a product of
Rational. CauseWay is a product of Devore software Incorporated.
Multi-Edit is a product of American Cybernetics. WATCOM F77 is a
product of WATCOM Incorporated. Turbo Pascal is a product of Borland
International Incorporated.
References
J.M. Prausnitz, T.F. Anderson, E.A. Grens, C.A. Eckert, R. Hsieh, J.P.
O'Connell, Computer Calculations for Multicomponent Vapor-Liquid
and Liquid-Liquid Equilibria, Prentice-Hall, Englewood Cliffs, NJ
(1980).
Chapter 8. FlowSheeting
By combining different ChemSep models you can actually
simulate (small) flowsheets. Several other common unit operations
(reactor, make-up stream, and stream splitter) are available to make
this possible. Simulating flowsheets with the utility program
fs is illustrated with several examples in this chapter.
8.1 Flowsheet Input File
The flowsheet utility uses a text file as input. You will have to make
this file yourself (with, for example, the edit option in
ChemSep's file menu). The ChemSep distribution
contains various examples (in the fs directory). To explain
the format of this file we will use an example where we simulate the
production of diethylether from an ethanol-water (85%-15%) feed, as
shown in Figure 8.1.
Figure 8.1. Diethylether Production
The conversion of the ethanol in the reactor is only 50% and the
unreacted ethanol has to be recycled. Pure ether (99.5%) and water
(with 1% ethanol) are the products. The reaction is:
2 C H OH -> (C H ) O + H O
2 5 2 5 2 2
In 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 comment
The 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).
[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,8
The 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).
[Feed]
Name=1
Temperature=40C
Pressure=1.01325bar
Rate=20mol/s
Z=0.15,0.85,0
The 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,17
where 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.
[Estimate]
Name=7
Temperature=80C
Pressure=1.01325bar
Rate=10mol/s
Z=0.85,0.15,0
The 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).
8.2 Flowsheet execution
If the flowsheet input is correct the flowsheet program starts the
execution of each units program with the specified file written to the
ChemSep.Fil file in the default directory (
ChemSep programs read this file to obtain the Sep-file to
simulate). The order of execution is the order as was specified in the
flowsheet section (so you can manipulate it). If a unit file or type
is missing execution of that unit is skipped.
8.3 Flowsheet Analysis
The flowsheet analysis of our depropanizer example is rather simple.
Flowsheet reports the incidence, adjacency, distance, and cycle
matrices for the specified flowsheet.
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-2
These 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.
8.4 Convergence
This short section displays whether the flowsheet was converged and,
if iteration was required, the attained convergence criterium and
number of iterations.
Flowsheet is solved.
Attained convergence = 0.000826326
Number of iterations = 10
Used relative differences.
8.5 Mass Balances
In the balances section total mass balance is given for each unit, and
for the overall flowsheet. If a flowsheet is converged, each balance
should equal zero or some small number (relative to the in- and outlet
flows of the specific unit). The balances are first written as the
names of the units and associated streams, then in the total molar
flowrates:
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-08
Care 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.
8.6 Stream Report
The stream report lists the streams after the execution has stopped.
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
8.7 Commandline Options
The flowsheet program has several options which can be specified on
the commandline when you start it. These are described below in some
examples followed by explanation (after the equal).
\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.
8.8 Other Unit Operations
If you want to simulate a complex flowsheet you need several other
unit operations models besides flash and column operations. The flash
unit can also be used to model heaters, coolers, or pumps. However, a
nonsense stream has to be added as most streams are either a vapor or
liquid. Here we discuss several other unit operations which are
commonly used in flowsheets. They have there own little programs for
which we supply the pascal code (to be compiled with Borland Pascal,
v7 or later). Most of them are limited to only a couple of hundred
lines which mostly cover the in- and output. If you require other unit
operations you could code them yourself (we welcome your pascal code
for unit operations to enhance the flowsheet capabilities).
8.8.1 Simple Reactor
A unit operation used in almost any flowsheet simulating a chemical
plant is the reactor. We have implemented a simple reactor model which
can handle multiple reactions with specified conversion(s). The
conversion is adapted in case one of the component flows is
constraining. The input file has a similar style as that of the
sep-files and flowsheet input file. The reactor specifications
(conversion, outlet temperature, pressure drop, and the reaction
stoichiometry) are given under the [Reactor] section which is
followed by a [Feeds] section following the sep-file format
(only one feed is allowed). The first line in the [Reactor]
section contains the reactor's name, followed by the outlet
temperature, pressure drop, number of components, and number of
reactions. Then, for each reaction, lines for the base component, the
conversion based on the base component feedflow, and the reaction's
stoichiometry coefficients (coefficients for each component; negative
for reactants, positive for products).
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 diphenyl
After 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
8.8.2 Make-Up Feeds
Often, when a flowsheet contains a recycled solvent or base material,
a (small) make-up feed is required. Since the exact flowrate of the
make-up feed is unknown, we made a unit operation that will add a
make-up feed (the first feed) to another stream (the second feed) to
obtain a specified total flowrate.
[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 flow
When 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 flow
Since 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.
8.8.3 Stream Splitter
Similarly to the Make-up unit, flowsheets containing recycles
sometimes need to purge some part of the recycle to prevent build-up
of various components in the cycle. This means that the recycle stream
must be split into two parts. The splits program implements
the stream split operation. The input file for the splitter consists
of a [Splitter] section, containing the name of the unit and a
splitfactor (which can range from zero to unity):
[Splitter]
Recycle purger
0.1
This 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.
8.9 Examples
Here we discuss several examples which have been run with the
flowsheet program. We will not supply all the details as they can be
found in the various files that are distributed with
ChemSep. Note that all these examples can be solved using the
nonequilibrium column models.
8.9.1 Extractive Distillation (PH)
We need to separate a equimolar mixture of methylcyclohexane (MCH) and
toluene, and do this by extractive distillation with phenol as
solvent. The flowsheet is shown in Figure 8.2. Valve trays are used
for both the columns using the design mode nonequilibrium simulator.
The Phenol recycle is cooled to 100 C. For a high purity of the
products the solvent feed to MCH/Toluene feed ratio as well as the
reflux ratio need to be sufficiently high (for the extractive column).
Figure 8.2. Extractive Distillation
8.9.2 Distillation with a Heterogeneous Azeotrope (BW)
Water and butanol form an azeotrope, so that they cannot be seaparated
by conventional distillation. However, at not too high temperature,
they form two liquid phases; an aqueous phase with little butanol and
a butanol phase with a large mole fraction of water. We can use a
decanter to separate the two liquid phases. A feed which is
predominatly water but contaminated with butanol (1 mole%) can be
separated into two pure products using such a decanter. The flowsheet
is shown in Figure 8.3. The first column produces a water bottom
product, and the vapor is fed to a condenser/cooler and then to the
decanter. There we obtain a water- and a butanol-rich stream which get
recycled to the columns. In a second column we can then produce
butanol as the butanol-rich recycle contains much more butanol than
the butanol-water azeotrope.
8.9.3 Distillation of a Pressure Sensitive Azeotrope (MA)
Methanol and acetone form also an azeotrope. However, the composition
of the azeotrope is sensitive to the pressure. We can make use of this
to separate the two components into pure products by operating two
columns at different pressures to change the azeotrope composition.
The separation of an equimolar feed of methanol and acetone is shown
in Figure 8.4. This type of azeotropic distillation is rare as the
azeotrope composition needs to be quite sensitive to the pressure in
order to obtain a recycle stream which is not unreasonably large, and
that the pressures are such that no special columns or equipment is
required.
8.9.4 Petyluk Columns (PETYLUK)
To lower the energy consumption of separation trains, two columns
separating 3 components can be replaced by one column with a condenser
and reboiler plus one column without condenser and reboiler. The feed
is fed to this (first) column, which receives a vapor to the bottom
and a liquid at the top from the (second) column with the condenser
and reboiler. In turn, the products of the first column are fed into
the second column at the sidestream stages. Figure 8.5 shows such a
configuration for the separations of three alcohols. To get this
flowsheet to run requires that you can solve both the columns
separately, which is not easy. The second column needed some specific
user initialization information to run. Afterwards convergence can be
obtained more easily by using the old results as initialization.
8.9.5 Extraction with Solvent Recovery (BP)
Extraction can be used to separate aromatics from parrafins. This is a
common type of separation in the crude oil refining. a simplified
example is shown here where we separate benzene from n-pentane, using
an extraction with sulfolane as solvent. The extractor is a rotating
disk contactor (RDC) operating at 50 C. The solvent is recovered from
the extract by "ordinary" distillation. UNIQUAC parameters for the
components are given in Table 8.1. The extraction is dependent on the
temperature as can be seen from the interaction parameters, which can
be used to manipulate the separation. See Figure 8.6 for the
flowsheet.
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