Ind. Eng. Chem. Res. 2010, 49, 3365–3373
3365
Dynamic Simulation of Continuous-Contact Separation Processes with the Moment Transformation Method Jonas Roininen* and Ville Alopaeus Aalto UniVersity School of Science and Technology, Department of Biotechnology and Chemical Technology, POB 16100, FIN-00076 AALTO, Finland
A model for the simulation of continuous-contact separation processes, based on heat and energy balances, multicomponent mass transfer, and chemical engineering correlations for mass and heat transfer coefficients and liquid holdup is developed and solved with the moment transformation method. Polynomials are used to approximate the concentration and flow rate profiles within the packing. The moment transformation method is used to transform the partial differential equations into a system of differential-algebraic equations that is solved with Newton’s method. The model can be used for both steady-state and dynamic simulation, and it can also be used for systems with axial dispersion. The features of the model are demonstrated with an example of ternary distillation. Introduction For decades, the means of simulating continuous-contact separation processes were limited to the equilibrium stage model or simplifying concepts such as HTU or HETP, which are discussed in great detail in classical chemical engineering textbooks such as King’s.1 In the 1970s and 1980s so-called rate based models emerged that were based on the equations of multicomponent mass transfer, rather than on the concept of the hypothetical equilibrium stage. Those models became quickly popular in both academia and industry, since they have many advantages over the simplified models. Although now considered obsolete, the equilibrium stage model is still widely used in design and simulation since it requires only very little information on the technical details of the column and relies only on an accurate VLE model, in contrast to rate-based models that require large amounts of additional physical and technical data. Mathematically, continuous-contact separation processes are modeled by partial differential equations for mass, momentum, and energy transport in both phases. These equations can be solved using a number of techniques, such as the finite difference method,2,3 the finite volume method, or polynomial approximations such as orthogonal collocation.4,5 The probably most popular rate-based method, the nonequilibrium stage model of Krishnamurthy and Taylor,6,7 is a finite volume model in the sense that the column is divided into a finite number of slices and the mole fractions and temperatures are averaged within those slices. All these models can be used for both steady state and dynamic simulations. Usually, Newton’s method is used to solve the nonlinear system of equations that arises from the temporal and spatial discretization of the partial differential equations. The interfacial mass and heat transfer across the phase boundary is often calculated with the Maxwell-Stefan equations or an appropriate simplified model. The method presented in this paper is related to the orthogonal collocation method that has been used for the simulation of packed-bed separation processes by Joseph and co-workers.4,5 Similar to the orthogonal collocation method, the moment transformation method also uses polynomials for the approximation of the concentration profiles. The difference is that in the * To whom correspondence should be addressed. Tel: +358 9 470 22633. Fax: +358 9 470 22694. E-mail:
[email protected].
moment transformation method the end points of the polynomials do not have to satisfy the boundary conditions exactly; the boundary conditions are automatically satisfied by the moment transformation equation. The moment transformation method in its present form with some relevant applications is presented in papers by Alopaeus and co-workers.8,9 In those papers, the focus is mainly on the dynamic simulation of chemical reactors. In this paper, the methodology is extended to the simulation of continuous-contact separation processes such as distillation, absorption, and stripping. Its features are demonstrated with an example of a distillation of a ternary mixture. Model Equations We start the model development by writing down the component and heat balance equations for the vapor and liquid phases. Component and thermal dispersion are included in the model, but phase dispersion is left unconsidered. A possible reaction, usually taking place in the liquid phase, could also be included in the equations but is not considered here. Further, energy buildup in the solid phase (packing) is neglected. With these assumptions, the partial differential equations are: Total buildup (i.e., the amount of vapor or liquid in moles per unit volume of packing) in either phase (continuity equations): ∂BL ∂L + )∂t ∂z ∂BV ∂V )∂t ∂z
nc
∑N
iVLa
(1)
i)1
nc
∑N
iVLa
(2)
i)1
Individual component buildup in either phase:
(
)
(3)
(
)
(4)
∂biL ∂(xiL) ∂xi ∂ + )D B + NiVLa ∂t ∂z ∂z L L ∂z ∂biV ∂(yiL) ∂yi ∂ + )DVBV - NiVLa ∂t ∂z ∂z ∂z Heat balance (in terms of enthalpy):
10.1021/ie901572a 2010 American Chemical Society Published on Web 03/15/2010
3366
Ind. Eng. Chem. Res., Vol. 49, No. 7, 2010
( (
) )
∂(BLHL) ∂(LHL) ∂HL ∂ + )R B + qVLa ∂t ∂z ∂z L L ∂z
(5)
∂(BVHV) ∂(LHV) ∂HV ∂ + )RVBV - qVLa ∂t ∂z ∂z ∂z
(6)
Equations 1-6 can be used to model almost any kind of continuous-contact separation processes such as absorption, stripping, and distillation. If the operation is countercurrent, then the flow rate in one phase must be defined negative. The axial dispersion coefficients D and R are not to be confused with the physical properties of diffusivity and thermal diffusivity. Instead, they should be understood as parameters in a model for axial dispersion in the packing that are determined experimentally for each packing type. It is also important to note that although the bulk phase mass and energy balances are written with the time derivative, the vapor-liquid interface and the vapor and liquid films are assumed to be at equilibrium, that is, there is no accumulation of mass and energy at the interface. The system of partial differential eqs 1-6 can be closed by using the concept of liquid holdup. Volumetric liquid holdup hL is a common parameter in packed column operation. It is defined as the fraction of the interstitial volume of the packing that is occupied by liquid. Liquid holdup is usually calculated from an empirical correlation for a specific packing type. Although holdup and pressure drop correlations apply strictly in the steady state only, they are frequently used in dynamic simulations.3,5,10,11 It can be assumed that the changes in holdup and pressure drop during the transients are so small that the correlations apply with good accuracy. This applies during transients from one steady state to another but not necessarily during startup and shutdown periods as the column may be initially completely dry. By using correlations for holdup and pressure drop, we can avoid solving the nonlinear equations of motion together with eqs 1-6. BL and BV are related to hL in the following way: hL )
BL cTLε
hV ) 1 - hL )
(7) BV cTVε
(8)
Using a correlation for liquid holdup, this yields two additional algebraic equations that have to be satisfied along with eqs 1-6:
0)1-
∂B*L ∂L* + )∂θ ∂ζ ∂B*V ∂V* )∂θ ∂ζ
nc
∑ (N
iVLa)*
(11)
i)1
nc
∑ (N
iVLa)*
(12)
i)1
∂b*iL ∂(xiL*) ∂x 1 ∂ + )B* i + (NiVLa)* ∂θ ∂ζ Pem,L ∂ζ L ∂ζ
(
)
(13)
∂b*iV ∂(yiV*) ∂y 1 ∂ + )B* i - (NiVLa)* ∂θ ∂ζ Pem,V ∂ζ V ∂ζ
(
)
(14) ∂E*L ∂(L*H*) ∂H* 1 ∂ L + )B*L L + (qVLa)* ∂θ ∂ζ Peh,L ∂ζ ∂ζ
(
)
(15) ∂E*V ∂(V*H*) ∂H* 1 ∂ V + )B*V V - (qVLa)* ∂θ ∂ζ Peh,V ∂ζ ∂ζ
(
)
(16) Equations 11-16 are applicable to both co- and countercurrent flow arrangements. If the flows are countercurrent, then one flow rate (here V) is defined negative. The only difference between the equations for the vapor and liquid phases is the sign of the interfacial mass and heat transfer fluxes. The variables to be solved are total buildup B* in both phases (2 variables), component buildups bi* in both phases (2 × [nc - 1] variables), molar flow rates L* and V* (2 variables), and energy buildup E* in both phases (2 variables). The enthalpy at a given point can be calculated from the energy buildup as E*Href B* The temperature can then be calculated as a function of the enthalpy and specified pressure and the molar composition. The corresponding partial differential and algebraic equations are equations for total buildup, eqs 11 and 12 (2 eqs), equations for individual component buildup, eqs 13 and 14 (2 × [nc - 1] eqs), equations for enthalpy, eqs 15 and 16 (2 eqs), and the holdup requirement, eqs 9 and 10 (2 eqs). As an alternative to using the total buildup, one could choose to use nc equations for individual component buildup and compute total buildup as H)
BL hL,corrcTLε
(9)
BV (1 - hL,corr)cTVε
(10)
0)1-
mol m-3 and uref ) 0.1 m s-1 were used. The enthalpy of an equimolar vapor mixture at column pressure and feed temperature was used as Href. The reference length ∆h is the height of a computational packing element (which in this case are all assumed to be equally high). In nondimensional form, eqs 1-6, then, become
Holdup and pressure drop are usually solved from a single set of correlations. Typical correlations are for example those of Buchanan12 for random packings and Rocha et al.13 for structured packings. In the present model, pressure drop does not need to be modeled (it can be set to zero or a fixed value), but a holdup correlation is strictly required. To apply the moment transformation, we have to cast eqs 1-6 into nondimensional form. For this, the dimensionless variables and quantities shown in Table 1 are defined. Instead of the Pe´clet number, some authors use the Bodenstein number (Bo) that is the inverse of the Pe´clet number. The reference values can be chosen rather arbitrarily; in this work cref ) 1000
nc
Bk )
∑b
ik
i)1
The closure models that are needed to solve the system are correlations for physical properties in either phase, enthalpy correlations for either phase, correlations for liquid holdup (and pressure drop), and model for interfacial mass transfer.
Ind. Eng. Chem. Res., Vol. 49, No. 7, 2010 Table 1. Dimensionless Variables and Quantities in eqs 11-15 dimensionless variables
definition
axial coordinate
ζ)
z ∆h
time
θ)
t τ
total buildup
B* )
B cref
component buildup
bi b* i ) cref
energy buildup
BH E* ) ) B*H* crefHref
enthalpy
H* )
liquid flow rate
L L* ) crefuref
vapor flow rate
V* )
dimensionless quantities volumetric interfacial mass transfer flux
a
(e.g., orthogonal collocation) is that the polynomials used to model the state profiles do not have to meet the boundary conditions exactly; they are automatically satisfied by the moment transformed model equations. The moment transformation technique is explained elsewhere,8,9 and only the most important equations are given here. We assume that the state profiles along the packing length can be approximated with a set of basis functions P(ζ) that are linear in terms of the coefficients. n
definition Niaτ (Nia)* ) cref
a
i
(17)
i
i)0
where f(ζ,θ) may be any of the dimensionless variables B*, b*, L*, V*, E* in either phase. The simplest choice of basis functions P(ζ) are the first n powers of ζ Pn(ζ) ) ζn
(18)
According to the Weierstrass approximation theorem, every continuous function in an interval can be uniformly approximated as closely as desired by a polynomial function and therefore polynomials are a natural choice for the basis functions. The number of transformed moments (the method order) is equal to the number of polynomial coefficients [degree + 1], degree being the maximum degree of the polynomials used as approximating functions. The moments can be calculated from the basis functions as mj(θ) )
qaτ crefHref
∑ w (θ)P (ζ)
f(ζ, θ) )
H Href
V crefuref
3367
∫
1
0
n
f(ζ, θ)ζjdζ )
∑ w (θ) ∫
1
i
i
0
Pi(ζ)ζjdζ, j )
volumetric interfacial heat transfer flux
(qa)* )
mass Pe´clet number
Pem )
uref∆h D
Now it is possible to establish a linear transformation between the moments and the polynomial weights. In a matrix form, it is
heat Pe´clet number
Peh )
uref∆h R
(m) ) [A](w)
reference time (s)
∆h τ) uref
The same reference values are used for vapor and liquid phases.
The pressure profile, if applicable, can be integrated using the pressure drop correlation. Moment Transformation Equations 11-16 constitute a set of partial differential equations (PDE) that have to be solved such that eqs 9 and 10 are satisfied at every time step. The fundamental idea of the moment transformation method is to reduce the system of PDE to a system of ordinary differential equations (ODE) by moment transformation and then follow the time evolution of the moments. In the present model the moment transformation method is applied to eqs 9-16, and we end up with a system of differential-algebraic equations (DAE). The advantage of the moment transformation method over other commonly used methods such as finite difference or finite volume is that the number of variables is greatly reduced, together with improved accuracy of the solution. This holds especially when compared to first order methods such as the nonequilibrium stage model. The advantage over other polynomial approximation methods
0...degree (19)
(20)
and the polynomial weights can be calculated at every time step by simple matrix inversion (w) ) [A]-1(m)
(21)
If the basis functions are ordinary polynomials, then the matrix [A] is the notorious Hilbert matrix, a classic example of ill-conditioned matrices. This means that the degree of the polynomials used is limited to maximally about 10. However, 10th degree polynomials are completely sufficient in most cases, and the degree may be limited to a much lower number for practical reasons (such as computational effort or solver convergence). The numerical behavior of the moment transformation can be improved to some extent by using orthogonal polynomials such as Chebyshev polynomials as basis functions, but as long as the degree of the polynomials is limited to 10, the choice of the basis functions is insignificant.8 Transformation of the Model Equations. Moment transformation of the dimensionless continuity and component and energy buildup equations is carried out in the sense of eq 19. According to the Leibniz rule, the time derivative can be written inside the integral, and since integration is a linear operation, every term can be integrated separately. The transformed equations are total buildup and nc - 1 component buildups in both phases (2 × nc × [degree + 1] × [number of elements] eqs)
3368
Ind. Eng. Chem. Res., Vol. 49, No. 7, 2010
dmBL j ) dθ
∫
1
0
nc
(ζj
∑ (N
iVLa)*
+ jζj-1L*)dζ + δ(0, j)L*in - L*out
i)1
(22) dmbiL j ) dθ
[
(
)]
B*L ∂xi ζ (N a)* + jζ L* dζ + iVL 0 Pem,L ∂ζ B*L ∂xi B*L ∂xi - xiL* (23) δ(0, j) xiL* Pem,L ∂ζ ζ)0 Pem,L ∂ζ ζ)1
dmBV j dθ
∫
1
(
)
j
j-1
) (
0
)
nc
∫ ∑ (-N 1
(ζj
iVLa)*
+ jζj-1V*)dζ + δ(0, j)V*out - V*in
i)1
(24) dmbiV j ) dθ
(
∫
1
0
The enthalpy gradients in eqs 26 and 27 are evaluated similarly using the energy buildup and the total buildup
)] )
( ) (
[
B*V ∂yi dζ + Pem,V ∂ζ B*V ∂yi - yiV* (25) Pem,V ∂ζ ζ)1
ζj(-NiVLa)* + jζj-1 V* -
B*V ∂yi δ(0, j) yiV* Pem,V ∂ζ
ζ)0
( )
*/∂ζ) */∂ζ) ∂H*L B*(∂E - E*(∂B ∂ E*L L L L L ) ) 2 ∂ζ ∂ζ B*L (B*)
Boundary Conditions. If the axial dispersion coefficients are zero (D ) 0, R ) 0, Pe f ∞), then the boundary conditions are straightforward: The inflow boundary value is simply the convective flux either from the separation unit boundary condition or from the previous subinterval. At the feed location, the liquid and vapor fractions of the feed are added to the respective boundary fluxes. The outflow boundary value is calculated from the values of the approximating polynomials at the boundary. If there is axial dispersion present, then Danckwerts boundary conditions apply.9 For the component profile for component i in the liquid phase, the boundary conditions are
(
heat buildup in both phases (2 × [degree + 1] × [number of elements] eqs) dmHL j ) dθ
[
( ) (
)]
B*L ∂H*L ζ (qVLa)* + jζ L*H*L dζ + 0 Peh,L ∂ζ B*L ∂H*L B*L ∂H*L - L*H*L δ(0, j) L*H*L Peh,L ∂ζ ζ)0 Peh,L ∂ζ ζ)1 (26)
dmHV j ) dθ
∫
1
(
∫
(
1
0
j
j-1
[
(
)]
B*V ∂H*V dζ + Peh,V ∂ζ B*V ∂H*V - V*H*V Peh,V ∂ζ ζ)1 (27)
ζj(-qVLa)* + jζj-1 V*H*V -
B*V ∂H*V δ(0, j) V*H*V Peh,V ∂ζ
) ( ζ)0
) )
The holdup discrepancy functions (eqs 9 and 10) are transformed similarly (2 × [degree + 1] × [number of elements] eqs)
∫
FLj ) 0 )
FVj ) 0 )
∫
1
0
1
0
(
(
1-
1-
)
BL ζjdζ hL,corrcTLε
)
BV ζjdζ (1 - hL,corr)cTVε
(28)
(29)
The total number of variables and equations after the moment transformation is then [2 × nc + 2 + 2] × [degree + 1] × [number of elements]. The remaining nonlinear integral in each of eqs 22-29 is easily calculated using an appropriate quadrature rule. In this work, 5-, 10-, and 15-point Gaussian quadratures were tested. With increasing number of variables used the order of the quadrature should be increased, whereas with a low number of variables a low-order quadrature rule suffices. Fifteen-point quadrature was used throughout this work to keep the results comparable. The mole fraction gradients in eqs 23 and 25 are evaluated using total and component buildups (that are calculated directly from their respective approximating polynomials) and the derivation rule for quotients as
( )
*iL /∂ζ) - b*iL(∂B*/∂ζ) ∂xi B*(∂b ∂ b*iL L L ) ) 2 ∂ζ ∂ζ B*L *) (BL
(30)
(31)
L
xiL* -
B*L ∂xi PeL ∂ζ
)
( ) ∂xi ∂ζ
+
) xiL*| ζ)1
(32)
ζ)0
ζ)1
)0
(33)
Physically the Danckwerts boundary conditions mean that the diffusive flux across the physical domain boundaries is zero. Since the moment transformation method allows every part of the boundary flux to be specified separately, the diffusive fluxes at both physical domain boundaries can be simply set to zero for the Danckwerts boundary conditions to be satisfied. There are two ways to increase the accuracy of the solution: by increasing the degree of the polynomials used or increasing the number of subintervals (elements) in the domain. In every element, eqs 22-27 must be solved separately, and some boundary conditions need to be specified at the element boundaries. The specification of these section boundary conditions is also discussed in our previous work.9 The correct boundary condition is continuity of the flux at the interface
(
xiL* -
B*L ∂xi PeL ∂ζ
) ( -
) xiL* -
ζ)1
B*L ∂xi PeL ∂ζ
)
+
(34) ζ)0
Numerically this can be satisfied for example by setting the diffusive flux in both boundaries to the arithmetic average of the values calculated from the approximating polynomials on either side. This boundary condition has the advantage that it is independent of the flux direction at the boundary B*L ∂xi ∂x 1 ≡ B* i PeL ∂ζ 2PeL L ∂ζ
[( )
ζ)1,cacld
(
+ B*L
∂xi ∂ζ
+
)
ζ)0,calcd
]
(35)
This formulation, however, leads to numerical instabilities at the element boundaries. Stable results are obtained when the arithmetic averages of the total buildup and the mole fraction gradients are calculated independently
[
B*L ∂xi 1 + ≡ + B*| (B*|L ζ)0,calcd) × PeL ∂ζ 4PeL L ζ)1,calcd ∂xi ∂xi + ∂ζ ζ)1,calcd ∂ζ
( |
|
+ ζ)0,calcd
)]
(36)
The values for total buildup at the element boundaries are calculated directly from their respective approximating poly-
Ind. Eng. Chem. Res., Vol. 49, No. 7, 2010
nomials and the mole fraction gradients using eq 30. The same holds also for the enthalpy gradients at the element boundaries. The total number of variables in a profile is [number of elements] × [degree + 1]. The method order is independent of the number of elements and is defined by the polynomial degree only. If zeroth degree polynomials on multiple elements are used (constant values of variables within the subintervals), then the model is reduced to the familiar nonequilibrium stage model. Feed Point. In a packed distillation column, the feed point is usually located between two separated sections of packing. In the simulation, the feed point must be correctly isolated, since there is most likely a discontinuity in the composition and temperature profiles at the feed point. For the simulation of a typical distillation column, then, at least two elements are needed, one for the stripping and one for the rectifying section. In the moment method, the treatment of the feed point poses no problem as long as it coincides with an element boundary. (The elements do not need to have same length.) Since the boundary conditions are simply added to the moments, as in eqs 22-27, any specification of the feed is allowed: liquid or vapor only or a specified mixture of liquid and vapor. The feed can also be flashed at column conditions, but this is not necessary. Solution of the System and Time Integration First we define a target vector with length [2 × nc + 2 + 2] × [degree + 1] × [number of elements] containing the functions defined in eqs 22-29
(
)
BL T dmbiL dmBV dmbiV dmHL dmHV j j j j j j ) dmj FLj FVj F dθ dθ dθ dθ dθ dθ (37)
mbiL mBV mbiV mHL mHV mLj mVj )T Sj ) (mBL j j j j j j
(38) The steady-state solution is obtained by solving the system of nonlinear equations
(39)
Time integration is done by adding the discretized time derivative to the residual on the RHS of eqs 22-27. There are three practical methods for the discretization of the time derivative: the implicit Euler method, the three-level implicit method, and the Crank-Nicolson scheme,14 of which the former two are the easiest to implement. In principle, higher-order methods could be used as well but usually second order accuracy is sufficient for most purposes. The three-level implicit method is stable at all time steps and is accurate to O(∆θ2).
(
3 k,n+1 1 k,n-1 m - 2mk,n j + mj dmk,n+1 2 j 2 j 0) dθ ∆θ
)
(40)
The implicit Euler method is also stable at all time steps but is only accurate to O(∆θ). 0)
(
dmk,n+1 mk,n+1 - mk,n j j j dθ ∆θ
)
the three-level implicit scheme should be used. Notice that the discretized time derivatives are added only to the time derivatives of the moments eqs 22-27 but not the transformed discrepancy functions 28 and 29. In this work, the Newton-Raphson method was used to solve the system of nonlinear equations. This works reasonably well for relatively small problems, but convergence and robustness issues arise if the system is large or complicated. This holds especially for polynomials of higher than third degree. Especially if the solution has to start with an initial guess that is far from the actual solution then it might be hard to achieve convergence. Ideally, an initial guess that is similar to the actual solution is used. When such an initial guess is not available, it is best to generate an intermediate solution where the column internal flow rates are close to the actual solution, even though the composition is not. Typically, the solutions at the time steps of a dynamic simulation converge much faster (within a few iterations) than a steady-state solution since the solution from the previous time step can be used as the initial guess. Future work should focus on finding a more reliable method for solving the equations. j with respect to Sj and then One possible idea is to linearize F solve the original system by iteration on the linearized system. Already in conjunction with the orthogonal collocation method of Srivastava and Joseph5 it was suggested that the steady-state solution could be obtained by solving the system of partial differential equations with Newton’s method. For the time integration, however, they used a Runge-Kutta scheme, but for this they needed to evaluate the derivative of the liquid holdup correlation with respect to liquid flow rate. It is clear that the method used here is more flexible in the sense that it does not rely on an analytical derivative of the liquid holdup correlation. Example: Distillation of the Ternary Mixture Benzene-Toluene-m-Xylene
The corresponding solution vector is
j (Sj) ) 0 F
3369
(41)
Both time integration schemes eventually lead to the same steady-state solution, but if exact transients are important, then
The method is demonstrated by performing simulations for the distillation of an equimolar mixture of three aromatic compounds, benzene, toluene, and m-xylene, in a column that consists of two packed sections of 1.0 m, one above and one below the feed point. The column diameter is assumed to be 0.1 m and the packed sections are assumed to be filled with a random packing. The packing characteristics used are those of Raschig rings with nominal packing size 6 mm. The mass and heat transfer fluxes were calculated with the full Maxwell-Stefan matrix method.15 The diffusion coefficients in the liquid phase were estimated using the Wilke-Chang,16 Vignes,17 and Kooijman and Taylor18 correlations. The vapor phase diffusion coefficients were calculated with an expression derived from kinetic gas theory.15 For the mass transfer coefficients and the interfacial area, Onda’s correlation with Bravo and Fair’s modification was used.15 The heat transfer coefficients were estimated from the mass transfer coefficients using the Chilton-Colburn analogy. The Peng-Robinson equation of state was used to calculate the vapor-liquid equilibrium at the interface. The liquid holdup was estimated with the correlation of Buchanan.12 The condenser was modeled as a total condenser without subcooling, and the reboiler was modeled as an equilibrium stage. Condenser and reboiler holdups were not taken into account in this example. The distillation was assumed to take place at atmospheric pressure. Pressure drop was neglected since the relative changes in pressure are usually small in columns operating at moderate pressure.19 (It was also checked, using Buchanan’s correlation,12 that the pressure drop was indeed negligible.) The feed was
3370
Ind. Eng. Chem. Res., Vol. 49, No. 7, 2010
Figure 3. Comparison of steady state liquid composition profiles calculated with 0th degree polynomials and increasing number of elements (equivalent to nonequilibrium stage model) to polynomial approximations. Figure 1. Schematic picture of the distillation column. Table 2. Parameters used in Column Simulation parameter
value
nominal packing size specific packing area void fraction column diameter feed feed temperature bottoms flow rate reflux ratio pressure
0.006 m 984 m-1 0.65 0.1 m 0.03 mol s-1 375.16 K 0.016 mol s-1 8 100 kPa
assumed to be liquid at its boiling point. Figure 1 shows the column configuration. The parameters used in the simulation are summarized in Table 2. The column specifications were fixed reflux ratio and bottom product flow rate. Their values are given in Table 2. Figure 2 shows the steady-state liquid composition profiles in the column with the assumption of no axial dispersion, calculated with 24 variables/profile and polynomials from zeroth to third degree. The respective numbers of elements are 24, 12, 8, and 6. The comparison of the zeroth degree approximation (nonequilibrium stage model) to the higher-order polynomial approximations reveals the magnitude of numerical diffusion
Figure 2. Steady-state liquid composition profiles in the column without axial dispersion. Profiles calculated with 24 variables/profile and increasing polynomial degree. Approximation with 0th degree polynomials is equivalent to the nonequilibrium stage model. Condenser at z ) 0, bottom of the packing at z ) 2 m, feed point at z ) 1 m, reboiler not shown.
in the nonequilibrium stage model. Obviously, higher-order polynomial approximations are capable of modeling the profiles with much lower error than nonequilibrium stage model. Figure 3 shows a comparison between nonequilibrium stage approximations with 24, 40, and 100 elements to a reference solution (see Error Analysis) and a polynomial solution with third degree polynomials and two elements (8 variables). The comparison shows that the nonequilibrium stage model converges only slowly toward the reference solution, whereas with a high-order approximation an accurate solution is obtained already with few variables. Even with 100 elements, the effect of numerical diffusion in the profiles is still clearly visible. On the other hand, the higher-order polynomial approximations are very accurate even with a low number of variables. Axial Dispersion Traditionally, axial dispersion has been neglected in studies of packed-bed separation processes. Obviously, the concept of axial dispersion is only hardly tangible within the framework of the traditional equilibrium and nonequilibrium stage models. In tall industrial distillation and absorption columns, axial dispersion is rarely a problem because of the ratio of packing size to packing height being very small. It can be expected, however that axial dispersion becomes important in small laboratory columns and modern microprocess devices. Axial dispersion is also important in bubble columns and extraction columns. There are a number of studies that address the issue of axial dispersion in column packings.20-24 The effect of axial dispersion on separation efficiency in the context of the classic design methods is also discussed in King’s book.1 With the present model, axial dispersion can be dealt with very easily and intuitively. Since there is no information available on the actual value of the axial dispersion coefficients in this case, we demonstrate a hypothetical case where all axial dispersion coefficients, that is, DL, DV, RL, and RV, were set to the same value. Naturally, the present model allows for independent variation of all these parameters, and it should not be generally assumed that they are the same in both phases. Since the two packing section of the column are physically separated, the Danckwerts boundary conditions were applied separately to both packing sections. Figure 4 shows the steadystate liquid composition profiles with the axial dispersion coefficients set to a plausible value of 0.005 m2 s-1, and Figure 5 shows the composition profiles when the axial dispersion coefficients are set to 0.01 m2 s-1 (a high degree of axial
Ind. Eng. Chem. Res., Vol. 49, No. 7, 2010
3371
Figure 4. Steady-state liquid composition profiles in the column with axial dispersion (DL ) DV ) RL ) RV ) 0.005 m2 s-1). Condenser at z ) 0, bottom of the packing at z ) 2 m, feed point at z ) 1 m, reboiler not shown. Profiles calculated with 24 variables/profile and increasing polynomial degree. Figure 6. Errors in the toluene mole fraction profiles in the base case and the axial dispersion case (DL ) DV ) RL ) RV ) 0.005 m2 s-1) compared to the reference solution.
numbers of elements affects the predicted column performance in the same way as simulated axial dispersion. Error Estimation
Figure 5. Steady-state liquid composition profiles in the column with axial dispersion (DL ) DV ) RL ) RV ) 0.01 m2 s-1). Condenser at z ) 0, bottom of the packing at z ) 2 m, feed point at z ) 1 m, reboiler not shown. Profiles calculated with 24 variables/profile and increasing polynomial degree.
dispersion). As expected, the performance of the column decreases with increasing axial dispersion coefficient. This can be seen in Table 3 that shows the top and bottom product compositions for the different axial dispersion coefficients. Especially the relative changes in the mole fractions of the light component in the bottoms product and the heavy component in the distillate are quite large. Table 3 also shows the product compositions calculated with the zeroth degree approximation (nonequilibrium stage model) for different numbers of elements; It demonstrates again that the numerical diffusion at low
An error estimate for the moment transformation method can be obtained by comparing the solutions calculated with a lower number of variables to a reference solution that can be obtained by the same method with a relatively higher number of variables. In this case, a solution with 10 subintervals and third degree polynomials (40 variables/profile, fourth order method) was found to be a valid reference solution. On the basis of earlier results,8,9 it can be expected that the error of the solution with 40 variables compared to a hypothetical “exact solution” is at least 2 orders of magnitude smaller than that of the solution with 24 variables and method order lower or equal to 4, so that it is meaningful to compare the lower-order solutions to this reference solution. Figure 6 shows the relative errors of the toluene mole fraction profile of different solutions (calculated as a discrete approximation of the root-mean-square) compared to the reference solution. The comparison shows clearly that the error decreases rapidly when the polynomial degree is increased. This tendency is especially clear when zeroth, first, and second order polynomials are compared. Moving on to thirdorder polynomials does not decrease the error drastically anymore. Also the rate at which the error decreases when the number of variables is increased is higher for high-order
Table 3. Comparison of the Distillate and Bottom Product Compositions of the Column with Different Values for the Axial Dispersion Coefficients and Different Number of Elements in the Nonequilibrium Stage Model (Calculated with the Present Model Using 0th Degree Polynomials) Axial dispersion coefficient (m2 s-1)
Nonequilibrium stage model Number of elements
0
0.005
0.01
24
40
100
0.714 0.279 0.007
0.708 0.280 0.012
0.700 0.283 0.017
0.710 0.278 0.013
0.712 0.278 0.010
0.713 0.279 0.009
0.001 0.381 0.619
0.006 0.380 0.615
0.012 0.378 0.610
0.004 0.382 0.614
0.002 0.382 0.616
0.001 0.381 0.618
distillate benzene toluene m-xylene bottom product benzene toluene m-xylene
3372
Ind. Eng. Chem. Res., Vol. 49, No. 7, 2010
mole fraction of toluene in the distillate goes through a maximum and is higher than the mole fraction of benzene between 50 and 175 s. This behavior occurs because the bottoms flow rate is fixed, so that additional liquid is taken out only from the top of the column. The lightest component (benzene) is withdrawn overproportionally before the column stabilizes again. Conclusions This work presents how the moment transformation method developed by Alopaeus and co-workers8,9 can be applied to the simulation of continuous-contact separation processes. The model can be used for both steady-state and dynamic simulations, and axial dispersion that is usually neglected in other models can be easily included. This work focuses on the theoretical framework of the model. Further work will focus on improving the robustness and speed of the numerical solution. In this work, the Newton-Raphson method was used to solve the group of nonlinear equations that is obtained after the moment transformation. Although it works well in cases where the initial guess is not very different from the solution, it is clear that a tailored solver is needed. With these improvements, the moment transformation method can become an interesting alternative to the traditional finite difference and dynamic nonequilibrium stage models currently used for dynamic simulation of continuous-contact separation processes. Nomenclature
Figure 7. Dynamic response of the column to a change of the reflux ratio from 8 to 2 at t ) 0 s. Liquid composition profiles at different column heights. The two profiles at z ) 1 m are at points immediately above and below the feed point, respectively. Solid line: no axial dispersion. Dashed line: DL ) DV ) RL ) RV ) 0.005 m2 s-1. 0: benzene. ): toluene. ∆: m-xylene.
polynomials. Especially the comparison to zeroth degree polynomial approximations (nonequilibrium stage, first order model) to the higher-order approximations is revealing in the sense that even with a very high number of variables, the accuracy of the zeroth degree approximation is far below the accuracy of the higher-order methods. In the axial dispersion case the error does not decrease as rapidly as in the base case but the trend is still the same. Dynamic Simulation To demonstrate this feature of the model, dynamic simulations were made where the reflux ratio was suddenly dropped from 8 to 2. The bottom product flow rate was kept constant at 0.016 mol s-1. This leads to a decrease in the column internal flow rates and decreased performance of the rectifying section of the column. For the reboiler and condenser, the simple total condenser and equilibrium stage models were used. Reboiler and condenser dynamics were not taken into account but could be easily included by defining variables for the reboiler and condenser holdups and solving them along with the column profiles (the same time integration method can be used). All simulations were performed using the three-level implicit method and time step ∆t ) 20 s. Figure 7 shows the dynamic response of the composition profiles at different heights in the column for the cases without and with axial dispersion. The
a ) specific vapor-liquid interfacial area (m-1) [A] ) transformation matrix between moments and polynomial weights biR ) molar buildup of component i in phase R (mol m-3) b*i ) dimensionless molar buildup of component i BR ) total molar buildup of phase R (mol m-3) B* ) dimensionless total buildup cref ) reference concentration (mol m-3) ctR ) total concentration in phase R (mol m-3) degree ) polynomial degree DR ) mass dispersion coefficient in phase R (m2 s-1) ER ) energy buildup in phase R (kJ m-3) E* ) dimensionless energy buildup f ) state profile function j ) vector of target functions F hL ) liquid holdup hV ) vapor holdup HR ) specific enthalpy of phase R (kJ mol-1) Href ) reference enthalpy (kJ mol-1) ∆h ) height of packing element (m) L ) liquid superficial flow rate, based on column cross-section (mol m-2 s-1) L* ) dimensionless liquid superficial flow rate m ) moment of a distribution (m) ) vector of moments NiVL ) vapor-liquid mass transfer flux of component i (mol m-2 s-1) P ) basis function Pem,R ) mass Pe´clet number in phase R Peh,R ) heat Pe´clet number in phase R q ) vapor-liquid heat flux (kW m-2) Sj ) solution vector t ) time (s) uref ) reference velocity (m s-1)
Ind. Eng. Chem. Res., Vol. 49, No. 7, 2010 V ) vapor superficial flow rate, based on column cross-section (mol m-2 s-1) V* ) dimensionless vapor superficial flow rate w ) polynomial coefficient (w) ) vector of polynomial coefficients xi ) mole fraction of component i in liquid phase yi ) mole fraction of component i in vapor phase z ) spatial coordinate (m) Greek Symbols R ) heat axial dispersion coefficient (m2 s-1) δ ) Kronecker delta ε ) void fraction of packing λ ) thermal conductivity (kW m-1 K-1) θ ) dimensionless time τ ) reference time (s) ζ ) dimensionless spatial coordinate Subscripts and Superscripts corr ) calculated from correlation i ) component i j ) index (indicating jth moment) k ) index (indicating kth equation/variable) L ) liquid phase ref ) reference value n ) index (indicating time step, degree of basis function) V ) vapor phase nc ) number of components * ) dimensionless quantity + ) right side of a boundary - ) left side of a boundary
Literature Cited (1) King, C. J. Separation Processes, 2nd ed.; McGraw-Hill: New York, 1980. (2) Hitch, D. M.; Rousseau, R. W.; Ferrell, J. K. Simulation of Continuous-Contact Separation Processes: Multicomponent, Adiabatic Absorption. Ind. Eng. Chem. Process Des. DeV. 1986, 25, 699. (3) Hitch, D. M.; Rousseau, R. W.; Ferrell, J. K. Simulation of Continuous-Contact Separation Processes: Unsteady-State, Multicomponent, Adiabatic Absorption. Ind. Eng. Chem. Res. 1987, 26, 1092. (4) Cho, Y. S.; Joseph, B. Reduced-Order Steady-State and Dynamic Models for Separation Processes. AIChE J. 1983, 29, 261. (5) Srivastava, R. K.; Joseph, B. Simulation of Packed-Bed Separation Processes Using Orthogonal Collocation. Comput. Chem. Eng. 1984, 8, 43.
3373
(6) Krishnamurthy, R.; Taylor, R. A Nonequilibrium Stage Model of Multicomponent Separation Processes. AIChE J. 1985, 31, 449. (7) Krishnamurthy, R.; Taylor, R. Simulation of Packed Distillation and Absorption Columns. Ind. Eng. Chem. Process Des. DeV. 1985, 24, 513. (8) Alopaeus, V.; Laavi, H.; Aittamaa, J. A Dynamic Model for PlugFlow Reactor State Profiles. Comput. Chem. Eng. 2008, 32, 1494. (9) Roininen, J.; Alopaeus, V. The Moment Method for OneDimensional Dynamic Reactor Models with Axial Dispersion. Comput. Chem. Eng., In press. (10) Gunaseelan, P.; Wankat, P. C. Transient Pressure and Flow Predictions for Concentrated Packed Absorbers using a Dynamic Nonequilibrium Model. Ind. Eng. Chem. Res. 2002, 41, 5775. (11) Ruivo, R.; Paiva, A.; Mota, J. P. B.; Simo˜es, P. Dynamic Model of a Countercurrent Packed Column Operating at High Pressure Conditions. J. Supercrit. Fluids 2004, 32, 183. (12) Buchanan, J. E. Pressure Gradient and Liquid Holdup in Irrigated Packed Towers. Ind. Eng. Chem. Fundam. 1969, 8, 502. (13) Rocha, J. A.; Bravo, J. L.; Fair, J. R. Distillation Columns Containing Structured Packings: A Comprehensive Model for Their Performance. 1. Hydraulic Models. Ind. Eng. Chem. Res. 1993, 32, 641. (14) Ferziger, J. H.; Peric´, M. Computational Methods for Fluid Dynamics, 3rd ed.; Springer: Berlin, 2002. (15) Taylor, R.; Krishna, R. Multicomponent Mass Transfer; Wiley: New York, 1993. (16) Wilke, C. R.; Chang, P. Correlation of Diffusion Coefficients in Dilute Solutions. AIChE J. 1955, 1, 264. (17) Vignes, A. Diffusion in Binary Solutions. Ind. Eng. Chem. Fundam. 1966, 5, 189. (18) Kooijman, H. A.; Taylor, R. Estimation of Diffusion Coefficients in Multicomponent Liquid Systems. Ind. Eng. Chem. Res. 1991, 30, 1217. (19) Choe, Y.-S.; Luyben, W. L. Rigorous Dynamic Models of Distillation Columns. Ind. Eng. Chem. Res. 1987, 26, 2158. (20) Dunn, W. E.; Vermeulen, T.; Wilke, C. R.; Word, T. T. Longitudinal Dispersion in Packed Gas-Absorption Columns. Ind. Eng. Chem. Fundam. 1977, 16, 116. (21) Ellenberger, J.; Krishna, R. Counter-current operation of structured catalytically packed distillation columns: pressure drop, holdup and mixing. Chem. Eng. Sci. 1999, 54, 1339. (22) Macı´as-Salinas, R.; Fair, J. R. Axial Mixing in Modern Packings, Gas, and Liquid Phases: II. Two-Phase Flow. AIChE J. 2000, 46, 79. (23) Baten, J. M.; Ellenberger, J.; Krishna, R. Radial and Axial Dispersion of the Liquid Phase within a KATAPAK-S Structure: Experiments vs CFD Simulations. Chem. Eng. Sci. 2001, 56, 813. (24) Zhang, W. W.; Liu, C. J.; Yuan, X. G.; Yu, G. C. Prediction of Axial Mixing for a Structured Packed Column Using a Two-equation Model. Chem. Eng. Technol. 2008, 31, 208.
ReceiVed for reView October 8, 2009 ReVised manuscript receiVed February 1, 2010 Accepted February 24, 2010 IE901572A
