6204
Ind. Eng. Chem. Res. 2010, 49, 6204–6214
Mass Conservative Solution of the Population Balance Equation Using the Least-Squares Spectral Element Method Zhengjie Zhu,*,† Carlos A. Dorao,*,‡ and Hugo A. Jakobsen*,† Department of Chemical Engineering and Department of Energy and Process Engineering, Norwegian UniVersity of Science and Technology, N-7491 Trondheim, Norway
In the standard formulation of the population balance equation that consists of breakage terms, significant loss of mass is observed for the dispersed phase. This mass loss is actually caused by the inexact conservation property reflected by many breakage kernels; hence, incorrect physical interpretations of the model simulations may be drawn. In this work, a constrained method is developed enforcing mass conservation. This numerical property is accomplished by adding an extra restriction to the original population balance equation in terms of the dispersed phase continuity equation through the Lagrange multipliers strategy. The discretized system resulting from applying the method to a two-phase population balance equation problem is symmetric and pseudopositive definite. Numerical experiments are carried out simulating the motion of a two-phase mixture passing through a 2D domain. The results obtained by the modified least-squares spectral element method show that the mass is conserved everywhere in the domain with high accuracy. 1. Introduction In multiphase flow problems, the dispersed phase distribution presents a strong effect on the hydrodynamic properties such as the mass and heat transfer fluxes. For this reason, considerable efforts have been made in order to develop polydispersed multifluid models with an inherent population balance module considering the effects of the variations in the size distribution of the dispersed phase. The dispersed phase size distribution can be described by a statistical Boltzmann-type equation, which is called the population balance equation (PBE). The PBE determines the temporal and spatial rates of change of a suitably defined distribution function.1 The population balance equation is a nonlinear partial integro-differential equation that must be solved numerically. Several methods have been proposed for the PBE computations. Among these are the Monte Carlo method,1-3 the method of the classes (CM),1 the quadrature method of moments (QMOM),4,5 and the direct quadrature method of moments (DQMOM).6 When adopting one of the above methods to solve the PBE together with a commercial computational fluid dynamics (CFD) code for the fluid dynamics of the model, the particle size distribution is represented by an averaged particle size like the Sauter mean diameter. The present status of PBE modeling of dispersed bubbly flows has been examined by Jakobsen et al.7 Recently, Dorao and Jakobsen8 discussed the applicability of the least-squares spectral element method (LSSEM) to solve the PBE. Subsequent works9-11 have been conducted to apply the least-squares method (LSM) to solve the more complicated PBE with different terms. The LSM is based on the minimization of a norm-equivalent functional. This method consists in finding the minimizer of the residual in a certain norm. Compared with other methods, the least-squares method has the following advantages (i) The LSM is highly accurate. The LSM is a higher order spectral element method and therefore possesses all * To whom correspondence should be addressed. E-mail:
[email protected] (Z.Z.);
[email protected] (C.A.D.);
[email protected] (H.A.J.). † Department of Chemical Engineering. ‡ Department of Energy and Process Engineering.
favorable properties of spectral element methods. The use of spectral methods leads to higher order convergence with h-refinement and even exponential convergence with p-refinement if the underlying exact solution is sufficiently smooth. In general, for a given accuracy the high order spectral element method requires only a few collocation points whereas the low order methods such as the finite volume method (FVM), require a large number of grid points. (ii) The weak formulation of the least-squares method is uniVersal to all types of partial differential equations. Unlike other numerical methods like the FDM and FVM which use central schemes to handle elliptic problems and upwind schemes to handle hyperbolic equations, the least-squares method solves these within one computational framework. (iii) The least-squares formulation leads to a symmetric, positive definite operator, which is suitable for both direct and iterative matrix solvers. (iv) The least-squares formulation has a so-called posterior error indicator in form of the residual. This indicator can serve as a convergence criterion for the nonlinear iteration process or be used in an adaptive mesh refinement procedure. None of the other numerical methods can provide such information without additional calculations. Reliable closures are required for the unknown source terms of the PBE. In chemical engineering, Coulaloglou et al.12 introduced a simple macroscopic formulation describing the particle interaction processes in agitated liquid-liquid dispersions. Starting out from the microscopic equations, the closure laws became an integrated part of the discrete numerical discretization scheme adopted. Lee et al.13 and Lehr et al.23 adopted the basic ideas of Coulaloglou et al.12 to formulate the population balance source term directly on the averaging scales to perform analysis of bubble breakage and coalescence in turbulent gas-liquid dispersions. Luo and Svendsen15 extended the work of Coulaloglou et al.12 formulating the population balance directly on the macroscopic scales where the closure laws for the source terms were integrated parts of the discrete numerical scheme used to solve the model equations.
10.1021/ie900710y 2010 American Chemical Society Published on Web 06/15/2010
Ind. Eng. Chem. Res., Vol. 49, No. 13, 2010
For the breakup processes, numerous source term closures have been formulated. Commonly used breakup closures in physical problems include the models proposed by Colaloglou and Tavlarides,12 Diemer and Olson,16 Martinez-Baza´n,17 Konno,18 etc. A shortcoming of the breakup kernels employed in many CFD models is that the conservation of volume/mass is not always fulfilled. Models that do not conserve volume/mass generally lead to an evolution of the predicted volume fraction which is unstable and not meaningful from a physical point of view.19 Inexact conservation of mass/volume has also been observed using the model presented in this work. The numerical verifications performed show that the system tends to lose mass when breakup takes place. Proper investigations of the reason for the poor mass conservation quality of the breakage kernels has not yet been reported. The main problem seems to be related to the properties of the daughter size redistribution function. To the authors’ knowledge, no remedy or modification in attempting to fix the critical problem has been made, either. This reveals the need to find a possible way to improve the mass conservation performance of the existing breakage kernels. The purpose of this article is to pay particular attention to the mass conservation of the PBE and propose a new approach for enforcing mass/volume conservation of the PBE when applied to physical problems. The problem can be considered as an optimality system consisting of a quadratic functional and a linear constraint equation. Different methods can be applied such as Lagrange multiplier methods, sensitivity or adjoint-based gradient methods, quasi-Newton methods, evolutionary algorithms, etc. In the least-squares formulation, the Lagrange multiplier method might be more amenable to the optimality problem due to its efficient uncoupling strategies.20,21 Other alternative is the penalty/least-squares finite element method.22,20 The new approach used in this paper modifies the standard nonconservative least-squares formulation by introducing a Lagrange multiplier associated with the dispersed phase continuity equation. The modified LSM implementation has been tested and compared with the standard one in simulations in which the Martinez-Baza´n breakage model is employed. The results obtained show that the modified LSM implementation can significantly improve the breakage kernel performance predicting a near zero mass conservation residual even when the breakup kernel is nonconservative. However, we shall also show that by introducing the Lagrange multiplier, one loses the positive-definitenesssone of the appealing features of the leastsquares method. A mixed problem for which one needs to establish an inf-sup inequality is obtained, and equal order interpolation is not possible anymore. The outline of this paper is as follows: Section 2 presents the bubble number density equation and the kernels of breakage. Section 3 introduces the basics of the least-squares spectral element method as well as the formulations of the standard and constrained least-squares method for solving the PBE considering breakage processes only. In section 4, we describe a test case and some numerical issues. In section 5, the results obtained from the numerical simulations are discussed and comparisons of the different methods are made with focus on the mass conservation property. Finally, conclusions are drawn in section 6. 2. Population Balance 2.2. Formulation of the Population Balance Equation. The population balance equation is a statistical equation that describes the size distribution of the dispersed phase in
6205
multiphase flows. The particle population may be regarded as being randomly distributed in the particle state space, which includes both the physical space and the space of internal coordinates. It is convenient to distinguish between the external coordinates z ) (zx, zy, zz) which may be employed to denote the position vector of the particle, and the internal coordinates r ) (r1, r2, ..., rd) which represents d different properties associated with the particle. The particle state vector (r, z) accounts for both the internal and external coordinates. We shall further let the domain Ωr represent the internal coordinates and Ωz be the domain of external coordinates, the latter being the set of points in physical space in which the particle is present. Moreover, Ω denotes the total state space that consists of both the internal space (Ωr) and the external coordinates (Ωz), i.e., Ω ) Ωr × Ωz
(1)
Let us characterize a two-phase flow field by a distribution function f′(V, z, t) (V ∈ Ωr and z ∈ Ωz), such that f′(V, z, t) dV represents the number of particles per unit volume at time t, at position z, with size between V and V + dV. In the generalized PBE, the convective transport of bubbles, their coalescence and breakup, and the interfacial mass transfer and the expansion of the gaseous phase may be taken into account. The derivation of the PBE has been outlined by several authors.23,24 A reduced volume-based PBE that includes coalescence, breakup, and source terms can be written as follows: D f' (V, z, t) ) B′(V, z, t) + C′(V, z, t) + S′(V, z, t) Dt
(2)
where B′ is the net bubble source due to breakup, C′ is the net bubble source due to coalescence, and S′ is the source production term. Here, and throughout the paper, “coalescence” means a process of growth by collision in which the merging of the colliding particle centers is instantaneous and the collision product may be treated as a single particle. The derivative of the left-hand side is given by D dV ∂ ∂ +∇·v+ ) Dt ∂t dt ∂V
(3)
Since our investigation is limited to breakup, the population balance can be written as ∂f' (V, z, t) + ∇ · (Vf' (V, z, t)) ) -b'(V)f' (V, z, t) + ∂t Vmax h'(V, W)b'(W)f' (W, z, t) dW (4) V
∫
On the right-hand side of eq 4, from left to right, the terms are recognized as the rate of death of bubbles of volume V due to breakup and the rate of birth of bubbles of volume V due to breakup of larger bubbles, respectively. b′(V) is the breakup frequency, h′(V, W) is the breakup kernel that represents the probability that a bubble of volume W splits up creating two bubbles of volumes V and W - V. Equation 4 can also be expressed in terms of the number density function using bubble length as the internal coordinate (ξ ∈ Ωr) by observing that f′ dV ) f dξ. The transient term and all the time dependences of all the terms can be removed at steady-state; hence ∇ · (Vf (ξ, z)) ) -b(ξ)f (ξ, z) +
∫
ξmax
ξ
h(ξ, ζ)b(ζ)f (ζ, z) dζ (5)
6206
Ind. Eng. Chem. Res., Vol. 49, No. 13, 2010
where ξ and ζ represent the bubble diameter. The details of the transformation of the PBE from the volume-based to the lengthbased is explained by Marchisio et al.4 The moments of the distribution function of different orders can be defined as
∫ ξ f (ξ, z) dξ
mk(z) )
k
(6)
In this work, the moments which have physical significance are (i) the number density (number of particles per unit volume) N(z) )
∫ f (ξ, z) dξ
(7)
ε)
∫ vol(ξ)f (ξ, z) dξ
(8)
vol(ξ) ) (π/6)ξ3 is the characteristic volume of a spherical particle with the diameter of ξ. The mean bubble volume can be computed by
∫ ( π6 ξ )f (ξ, z) dξ ∫ f (ξ, z) dξ 3
R(z) ) Vmean(z) ) N(z)
(9)
The integrals in eqs 7-9 should be evaluated over all possible bubble lengths. 2.2. Breakup Kernels. Martinez-Baza´n et al.17 proposed a bubble breakup frequency model based on a force balance considering turbulence stress and surface tension forces only. The novel breakup kernel was validated with the authors’ own experimental data for a turbulent jet. According to this model, the breakup frequency b(ξ, z, t) (s-1) is given by
β(εξ)2/3 - 12
b(ξ, z, t) ) Kg
σl Flξ
ξ
(10)
where the constant β ) 8.2 and the empirical parameter Kg ) 0.25 was determined experimentally for the case of air bubbles in water. Moreover, system properties like the liquid density (Fl) and liquid-gas interfacial tension (σl) must be provided. The flow property (ε) represents the averaged turbulent kinetic energy dissipation rate per unit mass. The daughter particle probability density function h*(ξ*, ζ, z, t) (m-1) of ξ* (i.e., ξ* ) ξ/ζ, the dimensionless diameter of the daughter bubble with respect to the parent bubble) was written as h*(ξ*, z, t) )
∫ξ
ζ
ξmax
ξmin
vol(ξ)h(ξ, ζ) dξ ) vol(ζ)
holds. Considering the daughter particles generated in all possible breakage events of a parent particle of size ζ, eq 13 implies that the total volume of all the daughter particles having sizes within the size interval [ξmin, ζ) must be equal to the volume of the parent particle. Figure 1 shows the performance of the Martinez-Baza´n breakage kernel. A loss of volume/mass is observed, and there is a significant reduction in the volume compared with the original parent particle volume as the dissipation rate increases, i.e.,
∫ξ
ζ min
(14)
3. Least-Squares Spectral Element Method The least-squares method is a well-established numerical method for solving a wide range of mathematical problems.26-28 The basic idea in the LSM is to minimize the integral of the square of the residual over the computational domain. For the cases in which the exact solution is sufficiently smooth, the convergence rate is exponential by increasing the polynomial degree of the approximation. 3.1. Standard LSM for the PBE. Consider a system of one linearized PBE that contains the convection and breakage terms only:
[ξ*2/3 - Λ5/3][(1 - ξ*3)2/9 - Λ5/3] d(ξ*)
where Λ ) ξc/ξ0, and ξc ) (12σ/(βF))3/5ε-2/5 is a critical diameter. The maximum peak of the particle distribution function (PDF) is located at ξ* ) 0.8. This value corresponds to the formation of two daughter particles of the same volume. The size PDF becomes wider as either ε or ζ are increased. The average dissipation rate is calculated by the ratio of the specific kinetic energy introduced by the gas and the mass of the liquid in the bubble column
vol(ξ)h(ξ, ζ) dξ < vol(ζ)
We can see that the loss of mass during the breakage event is due to the fact that the breakage kernel (daughter particle distribution function) does not fulfill the mass conservation constraint. This unphysical feature is highly undesirable in applications where the composition of the system in terms of disperse phase fraction is a key process variable. This problem will be tackled by the constrained LSM proposed in the next section.
L f ) g,
(11)
(13)
min
[ξ*2/3 - Λ5/3][(1 - ξ*3)2/9 - Λ5/3]
∫**
(12)
in which g represents the gravitational acceleration and jg is the volumetric flux of gas. This estimation of the averaged dissipation rate, ε, may be better than the local dissipation rate predicted by the k-ε turbulent model in 2D/3D simulation because the latter is merely a tuning parameter for the turbulent length scale representing a best fit of the velocity field in pipe flows.7,25 2.3. Mass Nonconservative Breakage Kernels. If the dispersed phase density is assumed to be a constant and there is no loss of mass during the breakage, then
(ii) the gas void fraction or gas holdup (volume of particles per unit volume) R(z) )
jggFl(1 - R) ) jgg Fl(1 - R)
in Ω ) Ωr × Ωz ∈ R3 B f ) gΓ,
on Γ
(15) (16)
with L f ) ∇ · (vf (ξ, z)) + b(ξ)f (ξ, z) -
∫
ξmax
ξ
b(ζ)h(ξ, ζ)f (ζ, z) dζ(17)
where Ω ) (ξmin, ξmax) × (0, L) × (0, H) ∈ R3 are bounded domains with a piecewise smooth boundary Γ, respectively. Ωz ) (0, L) × (0, H) is the external coordinate (Figure 2), while Ωr ) (ξmin, ξmax) is the internal coordinate. ξmin and ξmax are
Ind. Eng. Chem. Res., Vol. 49, No. 13, 2010
6207
Figure 1. Loss of mass of a parent particle of size ζ ) 0.02 m during a breakage event where the Martinez-Baza´n17 breakage kernel is used (tested for the air-water flow with the dissipation rate ε varies from 5 × 10-5 to 100 m2 · s-3).
and Φ is an arbitrary trial function (Φ ∈ Y). The minimization statement is then equivalent to the following: find f ∈ X(Ω) lim µf0
∂ J (f + µΦ) ) 0, ∂µ
∀Φ ∈ X(Ω)
(21)
where X(Ω) is the space of the admissible functions. Φ represents a set of arbitrary trial functions, and µ is a small perturbation. Consequently, the necessary condition can be written as follows: Find f ∈ X(Ω) such that Figure 2. External domain of the test case problem. The gas-liquid mixture is fed from the inlet at a constant velocity. Bubble breakup is taking place within the domain.
the minimum and maximum diameter of the particle, respectively. Note that Γ is the boundary of Ω. f is the unknown density function of x ) (ξ, zx, zy). g is a given set of vectorvalued functions. B is a given set of algebraic boundary operators, and gΓ is a given set of vector-valued functions on the boundary Γ, respectively. The least-squares method amounts to seeking a minimizer of the distance |L f - g|0 + |Bf - gΓ|0 in Y(Ω). We write the quadratic functional as J (f) )
1 1 |L f - g|02 + |B f - gΓ |02 2 2
|•| 02 ) 〈•, •〉0 )
∫ •• dΩ
(19)
|•| 02 ) 〈•, •〉0 )
∫ •• dΓ
(20)
Γ
A system of equations is obtained by substituting f ) f + µΦ and differentiating with respect to µ. Here µ is a real number
∀Φ ∈ X(Ω)
(22)
with A(f , Φ) ) 〈L f, L Φ〉Y(Ω) + 〈B f, B Φ〉Y(Γ)
(23)
F(Φ) ) 〈g, L Φ〉Y(Ω) + 〈gΓ, B Φ〉Y(Γ)
(24)
where A:X × X f R are symmetric, continuous bilinear forms, and F:X f R are continuous linear forms. 3.2. Constrained LSM for the PBE. The mass conservative least-squares solution is based on the minimization of the leastsquares functional (18) subjected to the constraint that mass is conserved ∇ · (Rv) ) 0
(18)
with the norms defined like Ω
A(f , Φ) ) F(Φ)
(25)
Equation 25 is obtained by integrating eq 5 over the entire internal coordinate domain and using the definition eq 8 of the void fraction (R). Equation 25 is also known as the continuity equation of the dispersed phase. Therefore, the constrained minimization problem can be expressed in the form: min J (f) subject to (25) f∈Y
where the constraint eq 25 holds in function space S.
(26)
6208
Ind. Eng. Chem. Res., Vol. 49, No. 13, 2010
To enforce mass conservation, the Lagrange multiplier method is used. The constrained least-squares solution results from the equilibrium point of the following extended Lagrangian functional I (f , λ) ) J (f) +
∫ λ[∇ · ( ∫ Ω
ξmax
ξmin
f (ξ, z)vol(ξ) dξV)] dΩ
(27) where J (f) represents the least-squares functional defined in eq 18. A system of equations is obtained by substituting f ) f + µΦ and λ ) λ + τΨ and taking the partial derivatives of I with respect to µ and τ: lim
µ,τf0
lim
µ,τf0
∂ I (f + µΦ, λ + τΨ) ) 0, ∂µ ∀Φ ∈ X(Ω)
∀Ψ ∈ M(Ω) (28)
∂ I (f + µΦ, λ + τΨ) ) 0, ∂τ ∀Φ ∈ X(Ω)
∀Ψ ∈ M(Ω) (29)
This result yields the following mixed formulations: Find (f, λ) ∈ X(Ω) × M(Ω) such that:
{
A(f, Φ) + B(Φ, λ) ) F(Φ) ∀Φ ∈ X(Ω) B(f, Ψ) ) 0 ∀Ψ ∈ M(Ω)
(30)
where the space X(Ω) consists of the function f and the space M(Ω) consists of the function λ. Note that in general X and M can be different function spaces. The bilinear form A and its corresponding right-hand side linear form F have been defined in eq 23 and eq 24. The bilinear form B(Φ, λ) is given by B(Φ, λ) ) 〈λ, ∇ · [(
∫
Ωr
Φvol(ξ) dξ)v]〉S(Ω)
(31)
Then, the constrained minimization problem (26) is equivalent to the unconstrained optimization problem of finding saddle points (f, λ) ∈ X × M of the Lagrangian functional, eq 31. Let us assume that the system (26) has a unique solution (f, λ) ∈ X × M. Moreover, |f| X + |λ| M e C(|f| Y + |λ| S)
(32)
in which, C is a positive constant and f ∈ X is the unique solution of the constrained minimization problem (26).20 3.3. Spectral Discretization of the Standard LSM. The discretization statement consists of searching the solution in a ˆ ), i.e. reduced subspace XN(Ω ˆf N(ξˆ , zˆx, zˆy) ∈ XN(Ω ˆ ) ) X ∩ PN(Ω ˆ ) ⊂ X(Ω ˆ)
(ξˆ 2 - 1) dLNξ(ξˆ )/dξˆ
φk(ξˆ ) )
The Nξ + 1 GLL points are the roots of the first derivative of the Legendre polynomial (LNξ) of degree Nξ, extended with the boundary points. The 1D basis functions for the other coordinates can be formulated in an analogous manner. The ˆfN(ξˆ , zˆx, zˆy) can then be expressed as a linear combination of the 1D basis functions: Nξ+1 Nx+1 Ny+1
ˆfN(ξˆ , zˆx, zˆy) )
∑ ∑ ∑ k)1
l)1
ˆ )) ) (Nξ + 1) × (Nx + 1) × (Ny + 1) ) Nt in R3 dim(PN(Ω (34) in which, Nt is the total number of degrees of freedom. As an example, the Nξ + 1 one-dimensional basis functions consist of Lagrange polynomials through the Gauss-LobattoLegendre (GLL) collocation points in ξˆ and can be expressed as
ˆfklmφk(ξˆ )φl(zˆx)φm(zˆy) (36)
m)1
where ˆfklm is the basis coefficient associated with each of the basis function φk(ξˆ )φl(zˆx)φm(zˆy), it follows that the basis coefficient ˆfklm is equal to the value of the discrete solution ˆfN(ξˆ , zˆx, zˆy) at the GLL points (or nodes) (ξˆ p, [zˆx]q, [zˆy]r), i.e. ˆf(ξˆ p, [zˆx]q, [zˆy]r) ) ˆfN(ξˆ p, [zˆx]q, [zˆy]r). Due to this property, we also say that we are using a nodal basis. The integration of the function is performed by use of a Gaussian quadrature rule:
∫ ∫ ∫ 1
1
-1
-1
Nξ+1 Nx+1 Ny+1
1
ˆf(ξˆ , zˆx, zˆy) dξˆ dzˆx dzˆy ≈ -1
∑ ∑ ∑ k)1
l)1
wkwlwmˆfklm
m)1
(37) where wk, wl, wm, ξˆ k, [zˆx]l, and [zˆy]m are the weights and points of the quadrature rule.29 If the number of the nodes is equal to the order of the basis functions, the derivatives of the two-dimensional basis function with respect to each reference variables evaluated at the GLL points are given by ∂(φk(ξˆ )φl(zˆx)φm(zˆy)) ∂ξ
|
) φ′k(ξˆ p)φl([zˆx]q)φm([zˆy]r)
ξp,[zˆx]q,[zˆy]r
) [dξ]pk[bx]ql[by]rm
(38) ∂(φk(ξˆ )φl(zˆx)φm(zˆy)) ∂zx
|
) φk(ξˆ p)φ′l([zˆx]q)φm([zˆy]r)
ξp,[zˆx]q,[zˆy]r
) [bξ]pk[dx]ql[by]rm
(39) ∂(φk(ξˆ )φl(zˆx)φm(zˆy)) ∂zy
|
) φk(ξˆ p)φl([zˆx]q)φm′ ([zˆy]r) ξp,[zˆx]q,[zˆy]r
) [bξ]pk[bx]ql[dy]rm
(33)
in which, the accent indicates that we discretized the solution on a three-dimensional reference domain. We can define the three-dimensional high-order polynomial ˆ ) by using the one-dimensional nodal basis space PN(Ω (Lagrangian interpolant) based on a tensor product Gauss-Lobatto (Legendre) grid. Note that
(35)
Nξ(Nξ + 1)LNξ(ξˆ )(ξˆ - ξˆ k)
(40) in which [bξˆ ], [bzˆx], and [bzˆy] are the one-dimensional mass matrices in ξˆ , zˆx, and zˆy, respectively. The [dξˆ ], [dzˆx], and [dzˆy] are the one-dimensional derivatives matrices in ξˆ , zˆx, and zˆy, respectively. Definitions of the one-dimensional mass and derivative matrices can be found in ref 27. Using the the approximations in (36)-(40), the following matrix systems are obtained for the standard PBE problem. Af ) F where the matrix A ∈ R
Nt×Nt
, vectors f, F ∈ R
(41) Nt×1
.
[A]ij ) A(Φj, Φi) ) 〈L Φj, L Φi〉Y(Ω) + 〈BΦj, BΦi〉Y(Ω)
(42)
Ind. Eng. Chem. Res., Vol. 49, No. 13, 2010
[F]i ) F(Φi) ) 〈[g]i, L Φi〉Y(Γ) + 〈[gΓ]i, BΦi〉Y(Γ)
(43)
[f]i ) f([x]i)
(44)
for 1 e i, j e Nt. 3.4. Spectral Discretization of the Constrained LSM. The system of the mass conservative PBE problem is formulated by introducing a similar polynomial expansion for the variable λ.
T
∑ ∑ ∑ k)1
l)1
t
6209
(51)
t
I :RN × RM f R I (w, η) ) J (w) + λT(Bw - G)
(52)
m)1
Inserting for the nonconstraint quadratic functional for the J (w),
ˆ )) ) (Mξ + 1) × (Mx + 1) × (My + 1) ) Mt in R3 dim(PM(Ω
1 J (w) ) wTAw - wTF 2
(53)
1 I (w, η) ) wTAw - wTF + λT(Bw - G) 2
(54)
t
in which M is the total number of degree of freedom. Using the approximation 45, together with (36)-(40) in the weak formulation (30), the following linear algebraic matrix system is obtained for the mass conservative PBE problem.
[ ][ ] [ ]
(46)
Aefe ) Fe
(47)
T
A B B 0
f F ) λ G
we can write
The stationary point, (w, η), of the Lagrangian, I (w, η), corresponds to a saddle point of I, hence
or
I (f, λ) ) min max I (w, λ) ) max min I (w, λ) (55) t
t
w∈RN λ∈RM
where,
[ ]
A BT , B 0
fe )
[]
f , λ
and Fe )
[] F G
in which A, F are defined in (42) and (43), respectively. The t t matrix B ∈ R(M ×N ), [B]ij ) B(Φj, Ψi) ) 〈Ψi, [(
∫
Ωr
Ae ) AeT
[]
yields
[
]
-1 T fe ) A B q * 0 -q
sup
B(fN, λN) g Kb |fN | X |λN | M
(56)
in which Kb is a positive constant which must be satisfied. Due to the similarity with the mixed formulation of the Stokes problem,30 the discrete space for the Lagrange multiplier, λ, needs to be particularly chosen to render a stable discretization. If the variable f is approximated by polynomials of degree Nt, the discrete inf-sup condition can be satisfied by approximating the variable λ by a polynomials of order Mt ) Nt - 2.
In order to assess the overall numerical properties of the standard and mass constrained LSM, the evolution of the bubble size distribution of an air-water two-phase mixture moving through a 2D rectangular channel (external coordinate Ωz), depicted in Figure 2, has been simulated.
p fe ) *0 0
Next consider
t
4. Test Case
therefore Ae is symmetric and Ae has real eigenvalues. Then, consider an arbitrary nontrivial vector
feTAefe ) pTAp > 0
inf
λN∈MN,λN*0 fN∈XN,fN*0
(49)
t
λ∈RM w∈RN
In this saddle point problem, the LadyzhenskajaBabusˇka-Brezzi (LBB) condition, commonly referred to as the (discrete) inf-sup condition owing to the equivalent form yields
Φjvol(ξ) dξ)v]〉Y(Ω) (48)
for 1 e i e Nt and 1 e j e Mt. t t t t t The vectors λ, G ∈ R(M ×1), the matrix Ae ∈ R(N +M )×(N +M ), t t and the vectors fe, Fe ∈ R(N +M )×1. The constrained LSSEM then has Mt more unknowns than the standard LSSEM for the same discretization. The extended matrix Ae resulting from the modified system is symmetric and pseudo-positive definite, which we explain in the following. We can see
and
-1
The discrete form of the Lagrangian, J (w, η), is defined by
λˆ klmψk(ξˆ )ψl(zˆx)ψm(zˆy) (45)
ˆ ) and on PM(Ω
Ae )
T
In conclusion, Ae must have both positive and negative eigenvalues. The algebraic problem does not correspond to either a minimization problem or a maximization problem.
Mξ+1 Mx+1 My+1
λˆ (ξˆ , zˆx, zˆy) )
-1 T
fe Aefe ) -q BA B q ) -p A p < 0 T
(50)
The rectangular channel of height H and length L has been set such that the origin of coordinate coincides with its left bottom corner. The air-water mixture is moving at a constant velocity in the zx-direction from left (inlet) to right (outlet). The dispersed flow is considered to be incompressible, hence the mass/volume conservation of the dispersed phase must hold.
6210
Ind. Eng. Chem. Res., Vol. 49, No. 13, 2010
Figure 3. GLL-GLL grid with eight spectral elements of order (a) 4, (b) 5, and (c) 6 is used in the 2D external coordinate (Ωz) of the PBE model.
Using Gauss’ theorem and the definition of the void fraction eq 8, eq 25 can be cast into
∫Ω ∇ · (Rv) dΩ ) ∫ (Rv) · n dΓ ) ∫ (n RV + n RV ) dΓ ) ∫ [n ( ∫ f(ξ, z)vol(ξ) dΩ )V + n ( ∫ f(ξ, z)vol(ξ) dΩ )V ] dΓ z
z
z
Γz
Γz
x
x
Γz
x
Ωr
y
Ωr
y
y
z
r
r
y
(57)
x
kernel has been adopted to close the PBE. For these simulations, the mixture moves with a velocity of 0.3 m/s along the centerline of a narrow channel of length L ) 1.0 m and height H ) 0.2 m. The breakage coefficient (Kg) varies from 3.0 × 10-3 to 2.0 × 10-2. An h refinement is employed by dividing the external physical domain into several subelements and assigning these subelements with polynomials of a proper degree. Figure 3 depicts the spectral element grids for the external coordinates used in the simulations. The grid consists of eight spectral subelements with polynomial order (Nx and Ny) varying from 4 to 6 for each of them. For the current test case the values of Nx ) Ny ) 5 are used, as illustrated in Figure 3b. For the internal coordinate, a polynomial order of 15 is used in the simulations. The velocity components (v ) [0.3, 0]) are prescribed on the boundary of the column at the in- and outlet. The bubble size distribution at the inlet is approximated as f0(ξ, zx, zy)
|
) zx)0
(
3N0 2 ξ3 ξ exp Vmean,0 Vmean,0
)
(58)
in which N0 ) 1 × 106 m-3 and Vmean,0 ) 3 × 10-7 m3 are the number density and the averaged bubble volume at the inlet, respectively. It can been seen from eq 58 that f0 has no dependence on zx or zy. The boundary integral eq 57 is approximated using a Gauss-Lobatto-Legendre quadrature rule to assess the conservation properties of the different formulations.
z
where Γz represents the boundary of the channel Ωz and where n ) [nx, ny] represents the outward directed unit vector. Here v ) [Vx, Vy] denotes the two-dimensional velocity field. We assume that only breakage is taking place as the mixture is moving through the channel, and the Martinez-Baza´n breakage
5. Results and Discussion For comparison, the absolute value of the boundary integrals eq 57 of the standard LSM eq 22 and the constrained LSM eq 30 has been plotted as a function of the breakage kernel parameter (Kg) in Figure 4.
Figure 4. Boundary integral values predicted by the standard (i.e., mass nonconservative) LSM and the constrained (i.e., mass conservative) LSM.
Ind. Eng. Chem. Res., Vol. 49, No. 13, 2010 5
Figure 5. Value of boundary integral is plotted against the polynomials degree of the internal coordinate at different breakage rates. Table 1. Loss of Total Mass (%) Calculated at Outlet Kg ) 3.00 × 10-3 Kg ) 6.00 × 10-3 Kg ) 1.00 × 10-2 nonconservative LSM conservative LSM
9.83
18.68
29.14
0.45 × 10-4
-0.094 × 10-4
0.26 × 10-4
The breakage parameter (Kg) has been varied from 3.0 × 10-3 to 2.0 × 10-2. The higher the Kg value, the stronger the breakage rate. From the figure, it can be seen that the value of boundary integral of the mass nonconservative LSM grows with an increasing breakage rate; hence, the conservation properties are increasingly affected as the breakage rate grows. In contrast, the result from the constrained LSM shows the divergence of the velocity times the void fraction is close to zero at any breakage rate, indicating that mass conservation has been strictly satisfied to machine accuracy. The absolute values of the boundary integral obtained with different breakage rates have also been plotted versus the polynomial degree of the internal coordinate ξ, Figure 5. The use of the different polynomial degrees allows one to investigate whether the polynomial degree plays a role in the extent of mass lost. The given results reveal that an increase in the polynomial degree gives almost no improvement in the quality of the mass conservation property although a higher polynomial degree is supposed to increase the accuracies of the integration of eqs 22 and 30. The integral values are significantly more affected by the breakage rate than by the polynomial degree. This trend is illustrated in Figure 4. Inferring from this fact, one can conclude that the loss of mass is not an inherent accuracy limitation of the LSM method. The percentages of mass loss are listed in Table 1, where the comparisons between the standard and constrained LSM at different breakage rates are shown. These results show that the mass loss increases from 9.83% at the weakest breakage rate to nearly 30% at the strongest breakage rate for the standard mass nonconservative LSM. The constrained LSM with the Lagrange multipliers, on the other hand, leads to a significant improvement to the mass conservation property. For the constrained LSM, the loss of
6211
mass is approximately 10 times less than for the standard LSM, thus the loss of mass between the inlet and the outlet of the channel can merely be neglected using the constrained LSM. In Figure 6, plots of the different evolutions of the density function f along zx at the centerline of the channel as predicted by the standard LSM and constrained LSM are given. It appears to be no significant deviation between the density functions from these two formulations. However, careful inspection of the plots reveals that notable deviations of the number densities exist for larger bubbles. Outtakes make these differences more obvious as the profile of f of the larger bubble part has been zoomed in. The number of large bubbles predicted by the mass conservative method is higher than that of the nonconservative method, whereas the number of the smaller bubbles remains almost unchanged. A comparison of these number density profiles indicates that the improved mass conservation of the LSM is obtained at the expense of emergence of large bubbles. This behavior can be understood if one realizes that the LSM tries to find a solution of f by minimizing the residuals of both the PBE (eq 5) and the dispersed phase continuity equation (eq 25) which are discretized on the domain of interest, simultaneously. Increasing the number of large bubbles is the numerically optimal way of making up the loss of mass, and keeping the density function satisfying the breakage dominant PBE as much as possible. In other words, the constrained LSM attempts to make “least” modifications on the original nonconservative curve (continuous line in Figure 6) by increasing the number of large bubbles. Figure 7 shows comparisons of the axial profiles of the void fraction R, the number density N, and the average bubble diameter ξmean as predicted by the standard LSM and the constrained LSM. Under the assumption of incompressible flow, for which the velocity v ) [0.3, 0] is constant in the zx direction, the void fraction, which is proportional to the mass of the gas phase, should also be constant along the channel length according to eq 25. Therefore, the measure of the mass conservation quality can, for example, be obtained by the values of the axial changes of the void fraction R. The constrained LSM fulfills the conservation property quite well as shown in Figure 7a. Again, it has been observed that for the standard LSM the calculated mass flux is substantially lower than the mass flux entering the channel. The number density computed by the constrained LSM is slightly higher than the corresponding value calculated with the standard LSM (Figure 7b) due to the fact that the constrained LSM increases the number density mainly in favor of the large bubbles, resulting in a small increase of the total bubble number in comparison to the nonconservative profile. Based on the observations in Figure 7a and b, larger average bubble diameters are expected from the mass conservative formulation due to the larger number of the large bubbles (Figure 7c). Therefore, we may conclude that the constrained LSM tends to reduce the intensity of the bubble breakage. In practical multiphase flow problems where the fluid dynamics part of the model is coupled with the PBE, there are two ways of treating the dispersed phase velocity. The first one is to assume that the velocity is a function of both the external and internal coordinates. The second one, which is more common, assumes that the velocity is only a function of the external coordinates. In the second approach, the PBE is basi-
6212
Ind. Eng. Chem. Res., Vol. 49, No. 13, 2010
Figure 6. Comparisons of the axial evolution of the bubble size distribution at zx ) (a) 0.00, (b) 0.25, (c) 0.50, and (d) 1.00 m at the centerline (zy ) 0.1 m) with Kg ) 6.0 × 10-3.
Figure 7. Comparisons of (a) the void fraction, (b) the number density, and (c) the average bubble diameter axial profiles at the centerline (zy ) 0.1 m) with Kg ) 6.0 × 10-3.
cally applied to compute a mean bubble diameter (ξmean), which is passed to the fluid dynamics part of the calculation and particularly used in the drag force correlation. On the basis of the results shown in Figures 6 and 7, we may conclude that if the constrained LSM is being used in the first approach, an enforcement of the mass conservation might cause some fluid dynamic inaccuracies due to the emergence of extra larger bubbles. In contrast, in the second approach, in which the PBE only contributes to the moments of the density function with a relatively coarse requirement on the density function profile, the constrained LSM guarantees the mass conservation property perfectly without any side effects. Figure 8 shows the axial evolution of the number density function at different breakage rates as obtained using the constrained LSM. As can be observed from all the three
figures, the small bubble number grows in the axial direction. The extent of growth becomes increasingly stronger with the gradually increasing breakage rate, leading to a smaller and smaller average bubble diameter. The results of the constrained LSM are highly desirable as the mass is conserved. Moreover, the shapes of the density functions clearly reveal that the flow is determined by the phenomena of bubble breakage. Although the enforcement of the mass conservation property is achieved at the expense of the emergence of a few large bubbles, the overall quality of the density function is quite satisfactory. Considering that in many engineering problems the moments of the density function usually get more attention than the density function itself, the constrained LSM can be an attractive fix to these mass nonconservative breakage kernels.
Ind. Eng. Chem. Res., Vol. 49, No. 13, 2010
6213
Figure 8. Mass conservative axial evolutions of bubble size distribution at various zx positions at centerline (zy ) 0.1 m) with (a) Kg ) 3.00 × 10-3, (b) Kg ) 6.00 × 10-3, and (c) Kg ) 1.00 × 10-2.
6. Conclusions The mass loss problem of the dispersed gas phase is due to the fact that the breakup kernel does not precisely fulfill the mass conservation property. The problem is merely physical, not numerical, and does not disappear through increasing the numerical accuracy. The practical applicability of the mass conservative formulation eq 30 has been examined in this work by performing numerical simulations of the PBE with the Martinez-Baza´n kernel. The modified mass conservative LSM is defined by introducing a new variable λ associated with the dispersed phase continuity equation which acts as a constraint to the original PBE. This constrained LSM attempts to find an optimal solution by finding the saddle point of the system that consists of the PBE and the dispersed phase continuity equation. The constrained LSM significantly improves the performance of the nonconservative breakup kernels and the given results show that by adding the variable λ we can guarantee mass conservation at the expense of emergence of larger bubbles. Numerically, the use of the constrained LSM is also questionable for three reasons. First of all, the use of the Lagrange multipliers undermines one of the main advantages of the LSM. Indeed, one sacrifices the positive definiteness of the leastsquares formulation which may led to indefinite problems similar to those arising in the standard Galerkin formulations of the Stokes equations.30 Second, in addition to the variables needed to obtain the first order least-squares formulations, one needs an additional Lagrange multiplier, λ. And finally, an additional compatibility requirement (e.g., the LBB condition), not present in the standard least-squares implementations, enters the formulation to guarantee a unique solution.
Nevertheless, this work offers insight into a fundamental, unsolved problem in the modeling of multiphase flows (i.e., the conservation of mass in population balance equation modeling), and the results are potentially useful in a variety of engineering applications. Acknowledgment The Ph.D. fellowship (Z.Z.) financed by the gas technology center at NTNU, SINTEF, and the Research Council of Norway through a strategic University Program (CARPET) is gratefully appreciated. Literature Cited (1) Ramkrishna, D. Population Balances, Theory and Applications to Particulate Systems in Engineering; Academic Press: San Diego, 2000. (2) Robert, C.; Casella, G. Monte Carlo Statistical Methods; Springer: New York, 2004. (3) Smith, M.; Matsoukas, T. Constant-number Monte Carlo simulation of population balances. Chem. Eng. Sci. 1998, 53, 1777–1786. (4) Marchisio, D.; Vigil, R.; Fox, R. Quadrature method of moments for aggregation-breakage processes. J. Colloid Interface Sci. 2003, 258, 322–334. (5) McGraw, R. Description of Aerosol Dynamics by the Quadrature Method of Moments. Aerosol Sci. Technol. 1997, 27, 255–265. (6) Marchisio, D.; Fox, R. Solution of population balance equations using the direct quadrature method of moments. J. Aerosol Sci. 2005, 36, 43–73. (7) Jakobsen, H. A.; Lindborg, H.; Dorao, C. A. Modeling of Bubble Column Reactors: Progress and Limitations. Ind. Eng. Chem. Res. 2005, 44, 5107–5151. (8) Dorao, C. A.; Jakobsen, H. A. A least squares method for the solution of population balance problems. Comput. Chem. Eng. 2006, 30, 535–547. (9) Dorao, C. A.; Jakobsen, H. A. Least-squares spectral method for solving advective population balance problems. J. Comput. Appl. Math. 2007, 201, 247–257.
6214
Ind. Eng. Chem. Res., Vol. 49, No. 13, 2010
(10) Dorao, C. A.; Jakobsen, H. A. Time-space-property least squares spectral method for population balance problems. Chem. Eng. Sci. 2007, 62, 1323–1333. (11) Zhu, Z.; Dorao, C. A.; Jakobsen, H. A. A least-squares method with direct minimization for the solution of the breakage-coalescence population balance equation. Math. Comput. Simul. 2008, 79, 716–727. (12) Coulaloglou, C.; Tavlarides, L. Description of interaction processes in agitated liquid-liquid dispersions. Chem. Eng. Sci. 1977, 32, 1289–1297. (13) Lee, C.; Erickson, L.; Glasgow, L. Bubble breakup and coalescence in turbulent gas-liquid dispersions. Chem. Eng. Commun. 1987, 59, 65–84. (14) Prince, M.; Blanch, H. Bubble coalescence and break-up in airsparged bubble columns. AIChE J. 1990, 36, 1485–1499. (15) Luo, H.; Svendsen, H. Theoretical model for drop and bubble breakup in turbulent dispersions. AIChE J. 1996, 42, 1225–1233. (16) Diemer, R.; Olson, J. A moment methodology for coagulation and breakage problems: Part 3-generalized daughter distribution functions. Chem. Eng. Sci. 2002, 57, 4187–4198. (17) Lasheras, J.; Eastwood, C.; Martinez-Baza´n, C.; Montane´s, J. A review of statistical models for the break-up of an immiscible fluid immersed into a fully developed turbulent flow. Int. J. Multiphase Flow 2002, 28, 247–278. (18) Konno, M.; Aoki, M.; Saito, S. Scale effect on breakup process in liquid-liquid agitated tanks. J. Chem. Eng. Jpn. 1983, 16, 312–319. (19) Zaccone, A.; Ga¨bler, A.; Maaβ, S.; Marchisio, D.; Kraume, M. Drop breakage in liquid-liquid stirred dispersions: Modelling of single drop breakage. Chem. Eng. Sci. 2007, 62, 6297–6307. (20) Bochev, P.; Gunzburger, M. Least-Squares Finite Element Methods for Optimality Systems Arising in Optimization and Control Problems. SIAM J. Numer. Anal. 2006, 43, 2517.
(21) Parrussini, L.; Pediroda, V.; Obayashi, S. Design under uncertainties of wings in transonic field. JSME Int. J. Ser. B 2005, 48, 218–223. (22) Bochev, P.; Gunzburger, M. Least-squares finite-element methods for optimization and control problems for the stokes equations. Comput. Math. Appl. 2004, 48, 1035–1057. (23) Lehr, F.; Millies, M.; Mewes, D. Bubble-Size distributions and flow fields in bubble columns. AIChE J. 2002, 48, 2426–2443. (24) Valentas, K.; Amundson, N. Breakage and Coalescence in Dispersed Phase Systems. Ind. Eng. Chem. Fundam. 1966, 5, 533–542. (25) Jakobsen, H. A. Chemical Reactor Modeling: Multiphase ReactiVe Flows; Springer Verlag: New York, 2008. (26) Jiang, B. The Least-Squares Finite Element Method: Theory and Applications in Computational Fluid Dynamics and Electromagnetics; 1998. (27) Pontaza, J.; Reddy, J. Spectral/hp least-squares finite element formulation for the Navier-Stokes equations. J. Comput. Phys. 2003, 190, 523–549. (28) Proot, M.; Gerritsma, M. Least-Squares Spectral Elements Applied to the Stokes Problem. J. Comput. Phys. 2002, 181, 454–477. (29) Karniadakis, G.; Sherwin, S. Spectral/hp element methods for CFD; Oxford University Press: New York, 1999. (30) Proot, M.; Gerritsma, M. Mass-and Momentum Conservation of the Least-Squares Spectral Element Method for the Stokes Problem. J. Sci. Comput. 2006, 27, 389–401.
ReceiVed for reView May 2, 2009 ReVised manuscript receiVed March 10, 2010 Accepted May 14, 2010 IE900710Y
