- Email: [email protected]
On approximate calculation of multicomponent mass transfer fluxes
Computers and Chemical Engineering 26 (2002) 461– 466 www.elsevier.com/locate/compchemeng
On approximate calculation of multicomponent mass transfer fluxes Ville Alopaeus * Neste Engineering Oy, P.O. Box 310, 06101 Por6oo, Finland Received 29 January 2001; received in revised form 21 November 2001; accepted 21 November 2001
Abstract The two methods for calculating mass transfer with the Maxwell– Stefan approach in unit operation or reactor models are considered. In the first case, mass transfer needs to be calculated separately from the unit operation model in an inner iteration loop. In which case the set of equations to be solved for the multicomponent mass transfer fluxes can be formulated, with suitable simplifications, in a compact way that appears to be noticeably simpler than those presented in the literature so far. In this case, no complicated matrix or scalar function calculations are needed during the iterative solution of the interphasial mass transfer fluxes, even though the multicomponent nature of the Maxwell– Stefan model is retained. Computationally, the solution of mass transfer fluxes is then equal or easier than the solution of an equilibrium flash model. In the second case, mass transfer is calculated simultaneously with other unit operation model equations. Then the mass transfer coefficient correlation and the high flux correction can be combined with the mass transfer flux equation, to be solved with the approximate matrix function method presented earlier by the author. This results in a very simple but still accurate formulation of the mass transfer fluxes. © 2002 Elsevier Science Ltd. All rights reserved. Keywords: Maxwell– Stefan diffusion model; Multicomponent mass transfer; Matrix approximations; Rate-based models
Nomenclature General symbols a a [A] [B] ct D d D – ij f( ) fk ( ) g (J) [K] [k]
linearization parameter for the high flux correction ( ) parameter ( ) matrix (general expression) ( ) matrix function of inverted binary diffusion coefficients (s/m2) total concentration (mol/m3) diffusion coefficient (m2/s) diameter (m) binary diffusion coefficients (m2/s) a function generally (various) mass transfer coefficient correlation function (various) gravitational constant (m/s2) vector of diffusion fluxes (mol/m2 s) matrix of distribution coefficients (n − 1 first as diagonal elements) ( ) mass transfer coefficient matrix (m/s)
* Fax: +358-10-45-27221. E-mail address: [email protected] (V. Alopaeus). 0098-1354/02/$ - see front matter © 2002 Elsevier Science Ltd. All rights reserved. PII: S 0 0 9 8 - 1 3 5 4 ( 0 1 ) 0 0 7 7 0 - 0
V. Alopaeus / Computers and Chemical Engineering 26 (2002) 461–466
462
l m, n n Nt (N) Re Sc Sh 6t (x) (z)
thickness of mass transfer film (m) parameters ( ) number of components ( ) total flux (mol/m2 s) vector of mass transfer fluxes (mol/m2 s) Reynolds number ( ) Schmidt number ( ) Sherwood number ( ) terminal velocity (m/s) vector of component mole fractions. Subscripts: V, dispersed phase; L, continuous phase; I interface (continuous phase side) ( ) vector of mass transfer ratios ( )
Greek letters m v z []] [c] ave D
error value ( ) viscosity (kg/ms) density (kg/m3) high flux correction matrix ( ) mass transfer rate factor matrix ( ) average mass transfer rate factor ( ) Driving force for mass transfer
1. Introduction In many instances, the Maxwell– Stefan approach has been adopted to model the multicomponent mass transfer. This approach is most commonly used with the linearized film model for mass transfer, i.e. constant diffusion coefficient matrix and total concentration in the mass transfer region (Taylor & Krishna, 1993). This approach, however, leads to complicated matrix function calculations both in estimating mass transfer coefficients and in the mass transfer flux expression. Some possibilities to simplify this situation in two different mass transfer calculation situations are discussed in this work. With proper simplifications, the solution of the Maxwell–Stefan equations are, in terms of computational work, quite close to the solution of pseudo-binary mass transfer equations, where the diffusional interactions are neglected. Computationally efficient mass transfer flux calculation methods are needed in several instances. For example, when a rate-based distillation column model is connected to a flowsheet with several other unit operation models, the mass transfer fluxes need to be solved countless times. Another case is a CFD model of two-phase flow, with mass transfer between the phases. Since an accurate solution of the flow fields requires from thousands to millions of control volumes, the mass transfer model needs to be reasonably fast in order to keep the calculation time of the overall model
tolerable. Presently an attempt is made to formulate an accurate, uniform mass transfer calculation approach, which can be used in all situations where mass transfer resistance plays a role, without the need to switch between models of various levels of rigor. Multicomponent mass transfer fluxes can be calculated with the following equation (N)= (J)+ Nt(x)= ct[kave][]I](Dx)+(xI)Nt
(1)
The high flux correction for the film model is []I]= [][exp[] − [I]] − 1
(2)
where the mass transfer rate factor matrix for the linearized mass transfer theory is []= Nt·l/ct[Dave] − 1 = Nt/ct[kave] − 1
(3)
The subscript ‘I’ in the high flux correction matrix refer to the form of the high flux correction equation, which in turn is related to the mass transfer direction and the compositions used in the convective part of the flux equation. The subscript ‘ave’ stands for the average values. Since [k] is assumed constant in the mass transfer region, average compositions should generally be used. Mass transfer correlations are often of the form Sh = a·Re nSc m. This can be generalized for the multicomponent mass transfer coefficients as (Stewart & Prober, 1964; Toor, 1964; Young & Stewart, 1986)
V. Alopaeus / Computers and Chemical Engineering 26 (2002) 461–466
a b·Re n·v m [k] = [D]+ [D]1 − m d·z m d
(4)
According to the Maxwell– Stefan diffusion model, the diffusion coefficient matrix can be calculated from the matrix of inverted binary diffusion coefficients (Krishna & Standart, 1976; Taylor & Krishna, 1993) [D]= [B] − 1
(5)
where Bii =
n xi xk + % ; D – in k = 1 D – jk k"j
Bij = − xi
1 1 − D – ij D – in
(6)
Combining Eqs. (4) and (5) we can write the mass transfer coefficient matrix in a general functional form as b·Re n·v m a [B]m − 1 =fk ([B]) [k]= [B] − 1 + d·z m d
(7)
Note that this formulation is not limited to Eq. (4), but any general mass transfer coefficient correlation can be cast in the same form, provided the matrix of mass transfer coefficients is a function of the diffusion coefficient matrix and scalar properties only. The matrix function in Eq. (7), as well as the other functions of the diffusion coefficient matrix, can be calculated approximately with the method described by Alopaeus and Norde´ n (1999). The approximation is because the diffusion coefficient matrix (and the matrix of inverted binary diffusion coefficients) usually has larger diagonal than off-diagonal elements. Another justification for the approximate method is due to the thermodynamic stability condition, as the diffusion coefficient matrix is positive definite for a stable phase (Taylor & Krishna, 1993). By keeping the diagonal dominating character of the diffusion coefficient matrix in mind, an approximate formulation can be derived. The approximation is based on an assumption that a matrix multiplication [F][M]2 can be neglected compared to [F]2[M], where [F] is the diagonal part of the diffusion coefficient matrix and [M] is the off-diagonal part. (Alopaeus & Norde´ n 1999)
2. Two calculation methods for the interphase mass transfer
2.1. Case 1: iteration of mass transfer fluxes separately from the unit operation model If the interphasial mass transfer fluxes need to be iterated separately from a unit operation or reactor model in an inner loop, the fluxes can be solved quite easily. For example in batch or tube reactor models, mass transfer needs to be calculated at each time or length step when integrating the reactor model. In these
463
cases, the interface compositions and the total flux across the interface need to be iterated in order to get the mass transfer fluxes. In case of non-isothermal system, the heat flux across the interface needs to be calculated as well. In these cases, the mass transfer coefficient matrix can be calculated with the bulk compositions with negligible error (Alopaeus & Aittamaa, 2000; Krishna, 1977; Young & Stewart, 1986). Then the mass transfer coefficient matrix can be calculated before the interface iteration, possibly with a suitable matrix function approximation. Furthermore, the equilibrium coefficients (the K values) and other physical properties can be calculated by using the bulk phase conditions. Then each iteration step requires only a simple square matrix– column matrix multiplication, and calculation of the linearization parameter for the high flux correction, as described by Alopaeus, Aittamaa, and Norde´ n (1999). The following calculation sequence is as follows: first, calculate all physical properties using bulk conditions (which are assumed constant during a mass transfer model solution); secondly, use the matrix function approximation for the mass transfer coefficient correlations; and thirdly, at each step of the iteration, calculate the high flux correction approximately in the following way. First, calculate the high flux linearization parameter for the film theory a=
1 1 − Nt exp(Nt)− 1
(8)
where =
(n− 1) n−1
(9)
ct % [k]ii i=1
is constant during the iteration and Nt is an iterated scalar variable. The linearization parameter is calculated similarly for both the phases. The driving forces and the fluxes are calculated with the following equations: (DxV)= (xV)− [K](xI)
(10)
(DxL)= (xI − xL)
(11)
(NV)= ctV[kV](DxV)− aVNt(DxV)+ Nt(xV)
(12)
(NL)= ctL[kL](DxL)− aLNt(DxL)+ Nt(xI)
(13)
Then the residuals for the equation solver are obtained by equating the mass transfer fluxes at both sides of the interface (n− 1 in number), and from the summation equation for the molar fractions at both sides of the interface (2 in number). Hence, we have n+1 variables and residual equations. This system of equations is quite linear, and with reasonable initial estimates, only one Newton step is
464
V. Alopaeus / Computers and Chemical Engineering 26 (2002) 461–466
usually needed for convergence. The only non-linearity comes from the term Nt(xI). Since the mass transfer fluxes are calculated frequently during the unit operation model calculation, good initial estimates are always available from the previous iteration (except possibly when the mass transfer is calculated for the first time). This approach can be compared to the equilibrium flash calculation needed in the equilibrium models. The number of iterated variables is the same (n + 1, namely the K factors and the vapor fraction), but the equations here are simpler. Hence, it can be assumed that the rate-based approach with proper simplifications is computationally equal or even easier than the equilibrium model. If the distribution coefficients are included in the iterated variables, appropriate residual functions should be added to the system. This residual is obtained by calculating the distribution coefficients using the interface conditions, and comparing these values with those of the guessed ones. Also, non-isothermal systems need an additional energy flux equation. This can, however, be calculated explicitly after the iteration of the mass transfer fluxes, if the high mass transfer rate effect (the Ackermann correction factor) to the heat transfer fluxes is neglected. Even though the formulation of the system of equations is different, the approximation errors resulting from this approach are similar than those considered by Alopaeus and Aittamaa (2000) and Alopaeus et al. (1999). Hence, numerical results when this method is applied are not considered further.
2.2. Case 2: simultaneous solution of mass transfer fluxes and unit operation models If the unit operation model consists of algebraic equations, it is usually preferable to calculate mass transfer as a part of that model. For example, in a multiphase continuous stirred tank reactor with mass transfer resistance, the material and energy balances of the reactor model can be solved simultaneously with the model for mass transfer fluxes. Then nested loops are avoided, and the same calculation tolerance can be used for all modelling equations. For the nested loops, the inner model needs to be calculated into significantly tighter tolerance than the outer loop in order to avoid convergence problems in the outer loop. In the simultaneous solution approach, the bulk properties change during the iterative solution of the unit operation or reactor model. Then the mass transfer coefficient correlations need to be recalculated at each step of the model solution. Here the approximate matrix function approach of Alopaeus and Norde´ n (1999) is extended to the whole multicomponent mass transfer flux equation, resulting in a considerable simplification in the computational
work. After observing that the equation for mass transfer fluxes is a function of diffusion coefficient matrix only, or further, a function of the matrix of inverted binary diffusion coefficients, [B], we can develop a straightforward simplification for the mass transfer flux equation. Combining Eqs. (1)– (3) and (7), we get
(N)= Nt (xI)+ exp
n
Nt f ([B]) − 1 − [I] ct k
−1
(Dx)
(14) Now the approximate method can be used as Aii = f(Bii )=
Aij = Bij
1 Nt exp −1 ct fk (Bii )
f(Bii )− f(Bjj ) Bii − Bjj
(15)
(16)
Note that f(*) states for a function in general, but fk (*) refers specifically to the mass transfer coefficient correlation. Also note also that fk ([B]) is a function of a matrix of size (n− 1)·(n − 1), but f(Bii ) is a scalar function. When the [A] matrix is composed of equations Eqs. (15) and (16), the flux equation (14) can be simply written as (N)= Nt((xI)+ [A](Dx))
(17)
Note that with this simplification no matrix multiplications, matrix inversions or other matrix function calculations are needed. Even then the multicomponent nature of the Maxwell–Stefan equations are retained, with diffusion cross-terms in the mass transfer coefficient matrix and the high flux correction. Eq. (14) and corresponding approximations (15) and (16) are singular at Nt = 0. In that case mass transfer is equimolar and no high flux correction is needed. Then the equations simplify to the following form (N)= ct fk ([B]) − 1 (Dx)
(18)
The matrix approximation is applied to the mass transfer coefficient correlation only, resulting in the following approximation equations: Aii = f(Bii )= 1/fk (Bii )
(19)
f(Bii )− f(Bjj ) Bii − Bjj
(20)
Aij = Bij
The case Nt = 0 is only a single singular point for the total flux, and in practice it is only encountered in cases where the mass transfer model is assumed to be exactly equimolar. This assumption is in most cases unnecessary. However, if the total flux is close to the CPU calculation tolerance during an iterative solution of the fluxes, the zero-flux approximation can be used instead of the complete function at that step.
V. Alopaeus / Computers and Chemical Engineering 26 (2002) 461–466
The second criterion is the error in the mass transfer ratios, i.e.
3. Numerical results Numerical comparisons related to Case 1 are available in the literature and there is no need to reproduce them here. A Monte Carlo type approach is adopted to compare the numerical schemes related to Case 2. Even though the flux calculations presented here are only for one phase, it is expected that the effect of various calculation methods to the overall characteristics of a unit operation model is revealed. Binary diffusion coefficients, density, viscosity, the mass transfer coefficient correlation parameters and the total flux divided by the total concentration were set as random variables within reasonable limits, in addition with the mole fractions at the two ends of a film. A large number of flux calculations were then carried out in order to get representative statistical results from the different methods. The following values were used for the parameters a=0 or 2 (an integer variable) b=0.25–0.75 d =10 − 5 –10 − 2 m v=0.1 –10 mPa s zD = 1.0–1200 kg zC = 600–1200 kg m= 0.25– 0.75 n = 0.25–0.75 Dij =0.5× 10 − 9 –5 × 10 − 9 m2/s xL,i, xI,i =0–1.0 Mole fractions were calculated so that first, a random value between 0 and 1 was given for each component, and then the sum was scaled to unity. In order to get consistent and physically reasonable values for the Reynolds number, Stokes law was used to calculate the terminal velocity. This was used even though the mass transfer model is not restricted to any special geometry, or more specifically to spherical particles moving slowly in continuous fluid. The Stokes law was used only to give reasonable values for the Reynolds number. Then we have gd 2 zD −zC 6t = (21) 18v and z6 d gz zD −zC d 3 Re = t = (22) v 18v 2 The following two criteria were used to compare the models: first, a relative difference between the exact and the approximate [A] matrix; and secondly the difference between the exact and the approximate mass transfer ratios. The first criterion is then m1 =
D
n−1 n−1
% % (A
ij,approx. i=1 j=1 n−1 n−1
−Aij,correct)2
% % (Aij,correct)2
i=1 j=1
465
(23)
(z)= (N)/Nt
D
(24)
n−1
m2 =
% (zi,approx. − zi,correct)2
i=1
(25)
n−1
2 % zi,correct
i=1
For comparisons, the mass transfer calculations were also carried out with Wilke’s effective diffusivity method, where a single diffusion coefficient is calculated for each component according to the following equation Di =
1− xi n x % j D – j=1 ij
(26)
j"i
These diffusion coefficients were then used in the mass transfer correlations for each component. This obviously results in only scalar equations, but otherwise the mass transfer formulation was the same than in the Maxwell –Stefan approach. The effect of high flux correction was also examined by calculating the mass transfer coefficient matrix and the corresponding mass transfer ratios without the high flux correction. Results from these four cases are shown in Table 1. The number of components was kept as an independent variable, since it was believed that the approximation error might be a function of the number of components. The average high flux correction is shown for illustrations as well. This parameter was calculated from the average diagonal values of the exact high flux correction matrix as 1− ]ii . Two sets of numerical comparisons were made, first with a moderate convective flux and the second with quite high convective flux. Both criteria resulted in quite similar trends in error values. The errors in all the methods tested were quite similar for both moderate total flux and for the high total flux, except when the high flux correction was neglected, which obviously resulted in significantly higher errors when the total flux was high. The errors when applying the matrix function approximation to the inverted binary diffusion coefficient matrix were of the order of 2%, and when the matrix approximation was applied to the diffusion coefficient matrix, the errors were about 0.5%. The Wilke method resulted in almost ten times higher errors than the simpler one of the two methods presented here. The same calculation sequence was repeated with a wider range of diffusion coefficients (0.1–10×10 − 9 m2/s). This was expected to increase the errors to some
V. Alopaeus / Computers and Chemical Engineering 26 (2002) 461–466
466
Table 1 Numerical comparison of various approximate methods to the exact eigenvalue solution n
]eff (%)
Criterion 1
Criterion 2
m1 (%)
m2 (%)
m3 (%)
Moderate total flux 3 2.339 4 2.270 5 2.235 8 2.222 10 2.194 15 2.241 30 2.246 50 2.320 Average 2.258
0.9443 1.837 2.330 2.876 2.873 2.695 2.146 1.754 2.182
0.2445 0.5005 0.6348 0.7977 0.8132 0.7689 0.6274 0.5113 0.6123
11.20 14.87 16.65 18.76 19.37 20.05 21.05 21.65 17.95
High total flux 3 46.69 4 45.15 5 44.63 8 44.11 10 43.81 15 44.86 30 45.87 50 45.62 Average 45.09
1.031 1.979 2.574 3.181 3.150 2.915 2.373 1.929 2.392
0.1979 0.4189 0.5411 0.6792 0.6815 0.647 0.5227 0.422 0.5138
11.51 15.30 17.32 19.64 20.09 20.62 21.86 22.45 18.60
m4 (%)
m1 (%)
m2 (%)
m3 (%)
2.094 1.975 1.921 1.899 1.878 1.93 1.957 2.03 1.960
0.8815 1.225 1.315 1.210 1.103 0.8609 0.4964 0.3148 0.926
0.2695 0.3746 0.3895 0.3559 0.3235 0.2500 0.1437 0.0902 0.2746
7.084 7.664 7.736 7.087 6.708 5.898 4.433 3.551 6.270
0.5932 0.7881 0.8445 0.8150 0.7372 0.5628 0.3218 0.2069 0.6087
0.1423 0.1965 0.2059 0.1973 0.1779 0.1354 0.0765 0.0490 0.1476
4.487 4.705 4.678 4.515 4.274 3.688 2.760 2.234 3.918
47.45 44.27 43.54 42.37 42.02 42.98 44.97 44.56 44.02
m4 (%)
2.354 2.243 2.192 2.191 2.162 2.217 2.232 2.310 2.238 23.36 21.59 21.48 21.97 22.05 22.31 22.69 22.65 22.26
m1: average error in the matrix function approximation; m2: average error in the matrix function approximation applied with the [D] matrix instead of [B]; m3: average error in the Wilke effective diffusivity method; m4: average error when the high flux correction is neglected; ]eff: average effect of the high flux correction.
extent, since with a wider range of diffusion coefficients, the non-diagonal elements of the diffusion coefficient matrix were higher. Indeed, the errors in the matrix approximations were roughly doubled, while the errors in the Wilke method increased by about 40%. Still the errors with the matrix approximations were less than 6 and 2% for the two methods, which can be considered acceptable in practical calculations.
4. Conclusion A compact formulation of mass transfer fluxes can be developed when the Maxwell– Stefan approach is adopted. Two cases can be distinguished. In the first case, interphasial mass transfer is calculated separately from the unit operation model. Then a simple iteration scheme can be used, where no complicated matrix calculations are needed. In the second case, mass transfer is calculated simultaneously with the other modelling equations. Then the residuals for the mass transfer fluxes can be calculated with a simple matrix formulation, where the mass transfer coefficient correlation and the high flux correction are combined in the flux expression. The errors in these methods can be considered negligible in practical calculations.
References Alopaeus, V., & Aittamaa, J. (2000). Appropriate simplifications in calculation of mass transfer in a multicomponent rate-based distillation tray model. Industrial and Engineering Chemistry Research, 39, 4336 – 4345. Alopaeus, V., Aittamaa, J., & Norde´ n, H. V. (1999). Approximate high flux corrections for multicomponent mass transfer models and some explicit methods. Chemical Engineering Science, 54, 4267 – 4271. Alopaeus, V., & Norde´ n, H. V. (1999). A calculation method for multicomponent mass transfer coefficient correlations. Computers and Chemical Engineering, 23, 1177 – 1182. Krishna, R. (1977). A film model analysis of non-equimolar distillation of multicomponent mixtures. Chemical Engineering Science, 32, 1197 – 1203. Krishna, R., & Standart, G. L. (1976). A multicomponent film model incorporating a general matrix method of solution to the Maxwell– Stefan equations. AIChE Journal, 22, 383 – 389. Stewart, W. E., & Prober, R. (1964). Matrix calculation of multicomponent mass transfer in isothermal systems. Industrial Engineering and Chemistry Fundamentals, 3, 224 – 235. Taylor, R., & Krishna, R. (1993). Multicomponent mass transfer. New York: Wiley. Toor, H. L. (1964). Solution of the linearized equations of multicomponent mass transfer. AIChE Journal, 10, 460 – 465. Young, T. C., & Stewart, W. E. (1986). Comparison of matrix approximations for multicomponent transfer calculations. Industry and Engineering Chemical Fundamental, 25, 476 – 482.
On approximate calculation of multicomponent mass transfer fluxes Ville Alopaeus * Neste Engineering Oy, P.O. Box 310, 06101 Por6oo, Finland Received 29 January 2001; received in revised form 21 November 2001; accepted 21 November 2001
Abstract The two methods for calculating mass transfer with the Maxwell– Stefan approach in unit operation or reactor models are considered. In the first case, mass transfer needs to be calculated separately from the unit operation model in an inner iteration loop. In which case the set of equations to be solved for the multicomponent mass transfer fluxes can be formulated, with suitable simplifications, in a compact way that appears to be noticeably simpler than those presented in the literature so far. In this case, no complicated matrix or scalar function calculations are needed during the iterative solution of the interphasial mass transfer fluxes, even though the multicomponent nature of the Maxwell– Stefan model is retained. Computationally, the solution of mass transfer fluxes is then equal or easier than the solution of an equilibrium flash model. In the second case, mass transfer is calculated simultaneously with other unit operation model equations. Then the mass transfer coefficient correlation and the high flux correction can be combined with the mass transfer flux equation, to be solved with the approximate matrix function method presented earlier by the author. This results in a very simple but still accurate formulation of the mass transfer fluxes. © 2002 Elsevier Science Ltd. All rights reserved. Keywords: Maxwell– Stefan diffusion model; Multicomponent mass transfer; Matrix approximations; Rate-based models
Nomenclature General symbols a a [A] [B] ct D d D – ij f( ) fk ( ) g (J) [K] [k]
linearization parameter for the high flux correction ( ) parameter ( ) matrix (general expression) ( ) matrix function of inverted binary diffusion coefficients (s/m2) total concentration (mol/m3) diffusion coefficient (m2/s) diameter (m) binary diffusion coefficients (m2/s) a function generally (various) mass transfer coefficient correlation function (various) gravitational constant (m/s2) vector of diffusion fluxes (mol/m2 s) matrix of distribution coefficients (n − 1 first as diagonal elements) ( ) mass transfer coefficient matrix (m/s)
* Fax: +358-10-45-27221. E-mail address: [email protected] (V. Alopaeus). 0098-1354/02/$ - see front matter © 2002 Elsevier Science Ltd. All rights reserved. PII: S 0 0 9 8 - 1 3 5 4 ( 0 1 ) 0 0 7 7 0 - 0
V. Alopaeus / Computers and Chemical Engineering 26 (2002) 461–466
462
l m, n n Nt (N) Re Sc Sh 6t (x) (z)
thickness of mass transfer film (m) parameters ( ) number of components ( ) total flux (mol/m2 s) vector of mass transfer fluxes (mol/m2 s) Reynolds number ( ) Schmidt number ( ) Sherwood number ( ) terminal velocity (m/s) vector of component mole fractions. Subscripts: V, dispersed phase; L, continuous phase; I interface (continuous phase side) ( ) vector of mass transfer ratios ( )
Greek letters m v z []] [c] ave D
error value ( ) viscosity (kg/ms) density (kg/m3) high flux correction matrix ( ) mass transfer rate factor matrix ( ) average mass transfer rate factor ( ) Driving force for mass transfer
1. Introduction In many instances, the Maxwell– Stefan approach has been adopted to model the multicomponent mass transfer. This approach is most commonly used with the linearized film model for mass transfer, i.e. constant diffusion coefficient matrix and total concentration in the mass transfer region (Taylor & Krishna, 1993). This approach, however, leads to complicated matrix function calculations both in estimating mass transfer coefficients and in the mass transfer flux expression. Some possibilities to simplify this situation in two different mass transfer calculation situations are discussed in this work. With proper simplifications, the solution of the Maxwell–Stefan equations are, in terms of computational work, quite close to the solution of pseudo-binary mass transfer equations, where the diffusional interactions are neglected. Computationally efficient mass transfer flux calculation methods are needed in several instances. For example, when a rate-based distillation column model is connected to a flowsheet with several other unit operation models, the mass transfer fluxes need to be solved countless times. Another case is a CFD model of two-phase flow, with mass transfer between the phases. Since an accurate solution of the flow fields requires from thousands to millions of control volumes, the mass transfer model needs to be reasonably fast in order to keep the calculation time of the overall model
tolerable. Presently an attempt is made to formulate an accurate, uniform mass transfer calculation approach, which can be used in all situations where mass transfer resistance plays a role, without the need to switch between models of various levels of rigor. Multicomponent mass transfer fluxes can be calculated with the following equation (N)= (J)+ Nt(x)= ct[kave][]I](Dx)+(xI)Nt
(1)
The high flux correction for the film model is []I]= [][exp[] − [I]] − 1
(2)
where the mass transfer rate factor matrix for the linearized mass transfer theory is []= Nt·l/ct[Dave] − 1 = Nt/ct[kave] − 1
(3)
The subscript ‘I’ in the high flux correction matrix refer to the form of the high flux correction equation, which in turn is related to the mass transfer direction and the compositions used in the convective part of the flux equation. The subscript ‘ave’ stands for the average values. Since [k] is assumed constant in the mass transfer region, average compositions should generally be used. Mass transfer correlations are often of the form Sh = a·Re nSc m. This can be generalized for the multicomponent mass transfer coefficients as (Stewart & Prober, 1964; Toor, 1964; Young & Stewart, 1986)
V. Alopaeus / Computers and Chemical Engineering 26 (2002) 461–466
a b·Re n·v m [k] = [D]+ [D]1 − m d·z m d
(4)
According to the Maxwell– Stefan diffusion model, the diffusion coefficient matrix can be calculated from the matrix of inverted binary diffusion coefficients (Krishna & Standart, 1976; Taylor & Krishna, 1993) [D]= [B] − 1
(5)
where Bii =
n xi xk + % ; D – in k = 1 D – jk k"j
Bij = − xi
1 1 − D – ij D – in
(6)
Combining Eqs. (4) and (5) we can write the mass transfer coefficient matrix in a general functional form as b·Re n·v m a [B]m − 1 =fk ([B]) [k]= [B] − 1 + d·z m d
(7)
Note that this formulation is not limited to Eq. (4), but any general mass transfer coefficient correlation can be cast in the same form, provided the matrix of mass transfer coefficients is a function of the diffusion coefficient matrix and scalar properties only. The matrix function in Eq. (7), as well as the other functions of the diffusion coefficient matrix, can be calculated approximately with the method described by Alopaeus and Norde´ n (1999). The approximation is because the diffusion coefficient matrix (and the matrix of inverted binary diffusion coefficients) usually has larger diagonal than off-diagonal elements. Another justification for the approximate method is due to the thermodynamic stability condition, as the diffusion coefficient matrix is positive definite for a stable phase (Taylor & Krishna, 1993). By keeping the diagonal dominating character of the diffusion coefficient matrix in mind, an approximate formulation can be derived. The approximation is based on an assumption that a matrix multiplication [F][M]2 can be neglected compared to [F]2[M], where [F] is the diagonal part of the diffusion coefficient matrix and [M] is the off-diagonal part. (Alopaeus & Norde´ n 1999)
2. Two calculation methods for the interphase mass transfer
2.1. Case 1: iteration of mass transfer fluxes separately from the unit operation model If the interphasial mass transfer fluxes need to be iterated separately from a unit operation or reactor model in an inner loop, the fluxes can be solved quite easily. For example in batch or tube reactor models, mass transfer needs to be calculated at each time or length step when integrating the reactor model. In these
463
cases, the interface compositions and the total flux across the interface need to be iterated in order to get the mass transfer fluxes. In case of non-isothermal system, the heat flux across the interface needs to be calculated as well. In these cases, the mass transfer coefficient matrix can be calculated with the bulk compositions with negligible error (Alopaeus & Aittamaa, 2000; Krishna, 1977; Young & Stewart, 1986). Then the mass transfer coefficient matrix can be calculated before the interface iteration, possibly with a suitable matrix function approximation. Furthermore, the equilibrium coefficients (the K values) and other physical properties can be calculated by using the bulk phase conditions. Then each iteration step requires only a simple square matrix– column matrix multiplication, and calculation of the linearization parameter for the high flux correction, as described by Alopaeus, Aittamaa, and Norde´ n (1999). The following calculation sequence is as follows: first, calculate all physical properties using bulk conditions (which are assumed constant during a mass transfer model solution); secondly, use the matrix function approximation for the mass transfer coefficient correlations; and thirdly, at each step of the iteration, calculate the high flux correction approximately in the following way. First, calculate the high flux linearization parameter for the film theory a=
1 1 − Nt exp(Nt)− 1
(8)
where =
(n− 1) n−1
(9)
ct % [k]ii i=1
is constant during the iteration and Nt is an iterated scalar variable. The linearization parameter is calculated similarly for both the phases. The driving forces and the fluxes are calculated with the following equations: (DxV)= (xV)− [K](xI)
(10)
(DxL)= (xI − xL)
(11)
(NV)= ctV[kV](DxV)− aVNt(DxV)+ Nt(xV)
(12)
(NL)= ctL[kL](DxL)− aLNt(DxL)+ Nt(xI)
(13)
Then the residuals for the equation solver are obtained by equating the mass transfer fluxes at both sides of the interface (n− 1 in number), and from the summation equation for the molar fractions at both sides of the interface (2 in number). Hence, we have n+1 variables and residual equations. This system of equations is quite linear, and with reasonable initial estimates, only one Newton step is
464
V. Alopaeus / Computers and Chemical Engineering 26 (2002) 461–466
usually needed for convergence. The only non-linearity comes from the term Nt(xI). Since the mass transfer fluxes are calculated frequently during the unit operation model calculation, good initial estimates are always available from the previous iteration (except possibly when the mass transfer is calculated for the first time). This approach can be compared to the equilibrium flash calculation needed in the equilibrium models. The number of iterated variables is the same (n + 1, namely the K factors and the vapor fraction), but the equations here are simpler. Hence, it can be assumed that the rate-based approach with proper simplifications is computationally equal or even easier than the equilibrium model. If the distribution coefficients are included in the iterated variables, appropriate residual functions should be added to the system. This residual is obtained by calculating the distribution coefficients using the interface conditions, and comparing these values with those of the guessed ones. Also, non-isothermal systems need an additional energy flux equation. This can, however, be calculated explicitly after the iteration of the mass transfer fluxes, if the high mass transfer rate effect (the Ackermann correction factor) to the heat transfer fluxes is neglected. Even though the formulation of the system of equations is different, the approximation errors resulting from this approach are similar than those considered by Alopaeus and Aittamaa (2000) and Alopaeus et al. (1999). Hence, numerical results when this method is applied are not considered further.
2.2. Case 2: simultaneous solution of mass transfer fluxes and unit operation models If the unit operation model consists of algebraic equations, it is usually preferable to calculate mass transfer as a part of that model. For example, in a multiphase continuous stirred tank reactor with mass transfer resistance, the material and energy balances of the reactor model can be solved simultaneously with the model for mass transfer fluxes. Then nested loops are avoided, and the same calculation tolerance can be used for all modelling equations. For the nested loops, the inner model needs to be calculated into significantly tighter tolerance than the outer loop in order to avoid convergence problems in the outer loop. In the simultaneous solution approach, the bulk properties change during the iterative solution of the unit operation or reactor model. Then the mass transfer coefficient correlations need to be recalculated at each step of the model solution. Here the approximate matrix function approach of Alopaeus and Norde´ n (1999) is extended to the whole multicomponent mass transfer flux equation, resulting in a considerable simplification in the computational
work. After observing that the equation for mass transfer fluxes is a function of diffusion coefficient matrix only, or further, a function of the matrix of inverted binary diffusion coefficients, [B], we can develop a straightforward simplification for the mass transfer flux equation. Combining Eqs. (1)– (3) and (7), we get
(N)= Nt (xI)+ exp
n
Nt f ([B]) − 1 − [I] ct k
−1
(Dx)
(14) Now the approximate method can be used as Aii = f(Bii )=
Aij = Bij
1 Nt exp −1 ct fk (Bii )
f(Bii )− f(Bjj ) Bii − Bjj
(15)
(16)
Note that f(*) states for a function in general, but fk (*) refers specifically to the mass transfer coefficient correlation. Also note also that fk ([B]) is a function of a matrix of size (n− 1)·(n − 1), but f(Bii ) is a scalar function. When the [A] matrix is composed of equations Eqs. (15) and (16), the flux equation (14) can be simply written as (N)= Nt((xI)+ [A](Dx))
(17)
Note that with this simplification no matrix multiplications, matrix inversions or other matrix function calculations are needed. Even then the multicomponent nature of the Maxwell–Stefan equations are retained, with diffusion cross-terms in the mass transfer coefficient matrix and the high flux correction. Eq. (14) and corresponding approximations (15) and (16) are singular at Nt = 0. In that case mass transfer is equimolar and no high flux correction is needed. Then the equations simplify to the following form (N)= ct fk ([B]) − 1 (Dx)
(18)
The matrix approximation is applied to the mass transfer coefficient correlation only, resulting in the following approximation equations: Aii = f(Bii )= 1/fk (Bii )
(19)
f(Bii )− f(Bjj ) Bii − Bjj
(20)
Aij = Bij
The case Nt = 0 is only a single singular point for the total flux, and in practice it is only encountered in cases where the mass transfer model is assumed to be exactly equimolar. This assumption is in most cases unnecessary. However, if the total flux is close to the CPU calculation tolerance during an iterative solution of the fluxes, the zero-flux approximation can be used instead of the complete function at that step.
V. Alopaeus / Computers and Chemical Engineering 26 (2002) 461–466
The second criterion is the error in the mass transfer ratios, i.e.
3. Numerical results Numerical comparisons related to Case 1 are available in the literature and there is no need to reproduce them here. A Monte Carlo type approach is adopted to compare the numerical schemes related to Case 2. Even though the flux calculations presented here are only for one phase, it is expected that the effect of various calculation methods to the overall characteristics of a unit operation model is revealed. Binary diffusion coefficients, density, viscosity, the mass transfer coefficient correlation parameters and the total flux divided by the total concentration were set as random variables within reasonable limits, in addition with the mole fractions at the two ends of a film. A large number of flux calculations were then carried out in order to get representative statistical results from the different methods. The following values were used for the parameters a=0 or 2 (an integer variable) b=0.25–0.75 d =10 − 5 –10 − 2 m v=0.1 –10 mPa s zD = 1.0–1200 kg zC = 600–1200 kg m= 0.25– 0.75 n = 0.25–0.75 Dij =0.5× 10 − 9 –5 × 10 − 9 m2/s xL,i, xI,i =0–1.0 Mole fractions were calculated so that first, a random value between 0 and 1 was given for each component, and then the sum was scaled to unity. In order to get consistent and physically reasonable values for the Reynolds number, Stokes law was used to calculate the terminal velocity. This was used even though the mass transfer model is not restricted to any special geometry, or more specifically to spherical particles moving slowly in continuous fluid. The Stokes law was used only to give reasonable values for the Reynolds number. Then we have gd 2 zD −zC 6t = (21) 18v and z6 d gz zD −zC d 3 Re = t = (22) v 18v 2 The following two criteria were used to compare the models: first, a relative difference between the exact and the approximate [A] matrix; and secondly the difference between the exact and the approximate mass transfer ratios. The first criterion is then m1 =
D
n−1 n−1
% % (A
ij,approx. i=1 j=1 n−1 n−1
−Aij,correct)2
% % (Aij,correct)2
i=1 j=1
465
(23)
(z)= (N)/Nt
D
(24)
n−1
m2 =
% (zi,approx. − zi,correct)2
i=1
(25)
n−1
2 % zi,correct
i=1
For comparisons, the mass transfer calculations were also carried out with Wilke’s effective diffusivity method, where a single diffusion coefficient is calculated for each component according to the following equation Di =
1− xi n x % j D – j=1 ij
(26)
j"i
These diffusion coefficients were then used in the mass transfer correlations for each component. This obviously results in only scalar equations, but otherwise the mass transfer formulation was the same than in the Maxwell –Stefan approach. The effect of high flux correction was also examined by calculating the mass transfer coefficient matrix and the corresponding mass transfer ratios without the high flux correction. Results from these four cases are shown in Table 1. The number of components was kept as an independent variable, since it was believed that the approximation error might be a function of the number of components. The average high flux correction is shown for illustrations as well. This parameter was calculated from the average diagonal values of the exact high flux correction matrix as 1− ]ii . Two sets of numerical comparisons were made, first with a moderate convective flux and the second with quite high convective flux. Both criteria resulted in quite similar trends in error values. The errors in all the methods tested were quite similar for both moderate total flux and for the high total flux, except when the high flux correction was neglected, which obviously resulted in significantly higher errors when the total flux was high. The errors when applying the matrix function approximation to the inverted binary diffusion coefficient matrix were of the order of 2%, and when the matrix approximation was applied to the diffusion coefficient matrix, the errors were about 0.5%. The Wilke method resulted in almost ten times higher errors than the simpler one of the two methods presented here. The same calculation sequence was repeated with a wider range of diffusion coefficients (0.1–10×10 − 9 m2/s). This was expected to increase the errors to some
V. Alopaeus / Computers and Chemical Engineering 26 (2002) 461–466
466
Table 1 Numerical comparison of various approximate methods to the exact eigenvalue solution n
]eff (%)
Criterion 1
Criterion 2
m1 (%)
m2 (%)
m3 (%)
Moderate total flux 3 2.339 4 2.270 5 2.235 8 2.222 10 2.194 15 2.241 30 2.246 50 2.320 Average 2.258
0.9443 1.837 2.330 2.876 2.873 2.695 2.146 1.754 2.182
0.2445 0.5005 0.6348 0.7977 0.8132 0.7689 0.6274 0.5113 0.6123
11.20 14.87 16.65 18.76 19.37 20.05 21.05 21.65 17.95
High total flux 3 46.69 4 45.15 5 44.63 8 44.11 10 43.81 15 44.86 30 45.87 50 45.62 Average 45.09
1.031 1.979 2.574 3.181 3.150 2.915 2.373 1.929 2.392
0.1979 0.4189 0.5411 0.6792 0.6815 0.647 0.5227 0.422 0.5138
11.51 15.30 17.32 19.64 20.09 20.62 21.86 22.45 18.60
m4 (%)
m1 (%)
m2 (%)
m3 (%)
2.094 1.975 1.921 1.899 1.878 1.93 1.957 2.03 1.960
0.8815 1.225 1.315 1.210 1.103 0.8609 0.4964 0.3148 0.926
0.2695 0.3746 0.3895 0.3559 0.3235 0.2500 0.1437 0.0902 0.2746
7.084 7.664 7.736 7.087 6.708 5.898 4.433 3.551 6.270
0.5932 0.7881 0.8445 0.8150 0.7372 0.5628 0.3218 0.2069 0.6087
0.1423 0.1965 0.2059 0.1973 0.1779 0.1354 0.0765 0.0490 0.1476
4.487 4.705 4.678 4.515 4.274 3.688 2.760 2.234 3.918
47.45 44.27 43.54 42.37 42.02 42.98 44.97 44.56 44.02
m4 (%)
2.354 2.243 2.192 2.191 2.162 2.217 2.232 2.310 2.238 23.36 21.59 21.48 21.97 22.05 22.31 22.69 22.65 22.26
m1: average error in the matrix function approximation; m2: average error in the matrix function approximation applied with the [D] matrix instead of [B]; m3: average error in the Wilke effective diffusivity method; m4: average error when the high flux correction is neglected; ]eff: average effect of the high flux correction.
extent, since with a wider range of diffusion coefficients, the non-diagonal elements of the diffusion coefficient matrix were higher. Indeed, the errors in the matrix approximations were roughly doubled, while the errors in the Wilke method increased by about 40%. Still the errors with the matrix approximations were less than 6 and 2% for the two methods, which can be considered acceptable in practical calculations.
4. Conclusion A compact formulation of mass transfer fluxes can be developed when the Maxwell– Stefan approach is adopted. Two cases can be distinguished. In the first case, interphasial mass transfer is calculated separately from the unit operation model. Then a simple iteration scheme can be used, where no complicated matrix calculations are needed. In the second case, mass transfer is calculated simultaneously with the other modelling equations. Then the residuals for the mass transfer fluxes can be calculated with a simple matrix formulation, where the mass transfer coefficient correlation and the high flux correction are combined in the flux expression. The errors in these methods can be considered negligible in practical calculations.
References Alopaeus, V., & Aittamaa, J. (2000). Appropriate simplifications in calculation of mass transfer in a multicomponent rate-based distillation tray model. Industrial and Engineering Chemistry Research, 39, 4336 – 4345. Alopaeus, V., Aittamaa, J., & Norde´ n, H. V. (1999). Approximate high flux corrections for multicomponent mass transfer models and some explicit methods. Chemical Engineering Science, 54, 4267 – 4271. Alopaeus, V., & Norde´ n, H. V. (1999). A calculation method for multicomponent mass transfer coefficient correlations. Computers and Chemical Engineering, 23, 1177 – 1182. Krishna, R. (1977). A film model analysis of non-equimolar distillation of multicomponent mixtures. Chemical Engineering Science, 32, 1197 – 1203. Krishna, R., & Standart, G. L. (1976). A multicomponent film model incorporating a general matrix method of solution to the Maxwell– Stefan equations. AIChE Journal, 22, 383 – 389. Stewart, W. E., & Prober, R. (1964). Matrix calculation of multicomponent mass transfer in isothermal systems. Industrial Engineering and Chemistry Fundamentals, 3, 224 – 235. Taylor, R., & Krishna, R. (1993). Multicomponent mass transfer. New York: Wiley. Toor, H. L. (1964). Solution of the linearized equations of multicomponent mass transfer. AIChE Journal, 10, 460 – 465. Young, T. C., & Stewart, W. E. (1986). Comparison of matrix approximations for multicomponent transfer calculations. Industry and Engineering Chemical Fundamental, 25, 476 – 482.