On the Coupled Solution of a Combined Population Balance Model

7994

Ind. Eng. Chem. Res. 2009, 48, 7994–8006

On the Coupled Solution of a Combined Population Balance Model Using the Least-Squares Spectral Element Method Zhengjie Zhu,*,† Carlos A. Dorao,*,‡ D. Lucas,*,§ 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, and Forschungszentrum Dresden-Rossendorf e.V. POB 510119, D-01314 Dresden, Germany

In this work, a cross-sectional averaged two-fluid model combined with a population balance model is applied to simulate the flow field and the bubble size distributions in a two-phase bubble column. The MartinezBazan breakage kernel and a modified Prince and Blanch coalescence kernel have been chosen to describe bubble breakage and bubble coalescence, respectively. In the present study, we discuss the use of a higherorder spectral element methodsthe least-squares methodsto compute the system of equations in a coupled manner. The least-squares method is highly accurate and has a number of advantages over the conventional numerical methods like the finite difference and finite volume methods. In contrast to the finite volume method, when desinging an overall solution algorithm, this least-squares method ensures that all the continuity equations are satisfied individually and it deals with both the convective and diffusive terms stably and accurately. The novel iterative algorithm solves the flow and the population balance equation in a coupled manner. The model has been validated against experimental data obtained for two-phase flow in a bubble column. The predicted bubble size distribution and other flow quantities are in good agreement with the experimental data. Introduction Bubble columns are widely used for carrying out gas-liquid and gas-liquid-solid reactions in a variety of industrial applications. The primary advantages of bubble columns are their simple construction due to having no moving parts, high gas-liquid interfacial area, good mass/heat transfer rates between gas and liquid phase, and large liquid holdup which is favorable for slow liquid phase reactions.1 Various methods have been implemented for the simulation of dispersed two-phase flows. They range from direct numerical simulations (DNS)2,3 to continuum-based models requiring a considerable number of closure relationships.4,5 In principle, multiphase flow can be mathematically described by the basic equations of fluid dynamics. However, due to the excessive computational cost, DNS simulations are unsuitable for practical industrial situations. Both the Eulerian-Lagrangian and Eulerian-Eulerian formulations may be applied to model the gas-liquid flows including stirred tanks and bubble columns. The EulerianLagrangian formulation6-8 treats the liquid phase as continuum while it tracks the motion of dispersed phase particles. Due to the large computational cost, the Eulerian-Lagrangian formulation is limited by the number of bubbles and size of the reactors.9 The Eulerian-Eulerian model (multifluid model)9-12 describes the motion of the two-phase mixture in a macroscopic sense. The two-fluid model consists of a set of averaged transport equations for each phase. The gas holdup and interfacial transfer terms couples these two phases. In the Eulerian-Eulerian approach, the computational cost is tractable and does not increase with the number of the bubbles. Finite * To whom correspondence should be addressed. E-mail: zhengjie.zhu@ chemeng.ntnu.no (Z.Z.); [email protected] (C.A.D.); [email protected] (D.L.); [email protected] (H.A.J.). † Department of Chemical Engineering, Norwegian University of Science and Technology. ‡ Department of Energy and Process Engineering, Norwegian University of Science and Technology. § Forschungszentrum Dresden-Rossendorf e.V. POB 510119.

volume discretization and iterative solutions, mostly based on the SIMPLE13 or PISO algorithms14 and their extensions, have been applied for numerical solution to the two-fluid models both in commercial software such as CFX, FLUENT, etc. and with in-house codes such as ESTEEM.15 In bubble columns, the dispersed phase flows can be described on the basis of a statistical Boltzmann-type equation called the population balance equation (PBE) that determines the temporal and spatial rates of change of a suitably defined distribution function (Figure 1).16 Reliable closures are required for the unknown terms. In chemical engineering, Coulaloglou et al.17 introduced the simpler macroscopic formulation, describing the interaction processes in agitated liquid-liquid dispersions. Drifting from the fundamental equations, the closures became an integrated part of the discrete numerical discretization scheme adopted. Glasgow et al.18 and Lehr et al.19 adopted the basic ideas of Coulaloglou17 to formulate the population balance

Figure 1. Bubble column. Bubbles flow upward with different bubble size distributions at different axial positions. Mean bubble volumes are indicated by dashed lines.

10.1021/ie900088q CCC: $40.75  2009 American Chemical Society Published on Web 07/02/2009

Ind. Eng. Chem. Res., Vol. 48, No. 17, 2009

source term directly on the averaging scales to perform an analysis of bubble breakage and coalescence in turbulent gas-liquid dispersions. Luo and Svendsen20 extended the work of Coulaloglou17 whom formulated 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. A summary of fluid dynamic bubble columns with population balance modeling has been given in refs 21 and 22. Bubble size distributions strongly influence mass transfer and hydrodynamics in a bubble column reactor. Therefore, a combination of flow and PBE is recommended, where the hydrodynamics may be resolved with high resolution of the two phase flow in external coordinate, while the PBE accounts for the bubble breakup and coalescence and predict the size distribution of the dispersed phase. Several methods have been proposed for combining PBE computations with commercial computational fluid dynamics (CFD) codes. Among these are the Monte Carlo method,16,23,24 the methods of the classes (CM),16 the quadrature method of moments (QMOM),25,26 and the direct quadrature method of moments (DQMOM).27 When adapting one of the above methods to solve the PBE together with a commercial CFD code, the particle size distribution is represented by averaged particle sizes like the Sauter mean diameter for fluid dynamics. Two exceptions are the recently developed inhomogeneous MUSIG model28,29 and a similar model by Sha et al.,30 which allows the consideration of size dependent momentum equations. The main focus of the present work is the investigation of a high-order spectral element methodsthe least-squares method (LSM)sto solve the Eulerian-Eulerian two-fluid model coupled with the PBE. 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 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 finite volume method (FVM) requires a large number of grid cells. (ii) The weak formulation of the least-squares method is uniVersal to all types of partial differential equations. Unlike other numerical methods which use central schemes to handle elliptic problems and upwinding schemes to handle hyperbolic equations, the least-squares method solve these within one mathematical/computational framework. (iii) The least-squares formulation has a so-called a posteriori error indicator in the form of the residue. This reliable indicator can be applied to check convergence or be used for adaptive mesh refinement. No other numerical method can provide such information without additional calculation. (iv) A common limitation with the FVM in computational fluid dynamic is the lack of possibilities to construct pressure correction and void schemes which enforce that all mass balances are fulfilled separately.21,22 That is, for each conservation equtation expressed in the flux form, the FVM ensures that all the fluxes are balanced (both locally and globally); hence, the continuity equations that are solved are satisfied. However,

7995

to design proper equations to determine the volume fraction and pressure correction variables particular combinations of the phasic continuity and momentum equations are used. For the convenience of making generic software and stability and accuracy reasons, both the volume fraction and pressure correction variables are often computed from a weighted sum of the phasic equations; hence, we cannot require mass conservation of each individual phase separately. The LSM is not restricted by such limitations, and no pressure correction equation is required. (v) For a coupled problem of the flow and the PBE, in the FVM, operator splitting is normally employed solving the flow equations in a separate loop and the PBE in a second operator.31 This splitting introduces additional errors, and the decoupled approach is not optimal considering computational cost whereas the LSM does not need such an approach and it minimizes the overall residue of the system. The performance of the LSM to solve the hyperbolic equation,32 the Navier-Stokes equation and incompressible flow problems33 have already been investigated by several groups. Recently, Dorao and Jakobsen34 discussed the applicability of the least-squares spectral element method (LSSEM) to solve the PBE. Subsequent works35-37 have been conducted to apply the LSM to solve the more complicated PBE with different terms. Despite the fact that several attempts have been made to explore the capability and numerical properties of the LSM with the help of a constructed analytical solution,32,33,35-38 only a few contributions has been reported applying this method to real physical problems such as multiphase chemical reactor flows. In particular, the use of the least-squares method in multiphase flow calculations combined with the PBE has not yet been investigated in spite of its proven capability. The complexity of this method is a main issue that hinders its broad applications. In addition, the LSM requires more time to evaluate each iteration than the low-order methods. However, this disadvantage associated with the LSM is generally compensated because this method requires considerably fewer discretization points to ensure the same solution accuracy, and the number of interations needed is normally much lower. Hence, considering the overall sequence of solution algorithm operations, the high-order LSM is usually more efficient than low-order methods. A brief review of the previous work and a discussion on the advantages of the LSM is included here to point out the need to study the applicability of using the LSM to solve the multifluid model that is coupled with the PBE. To that end, in the present article, a steady-state two-fluid model coupled with the PBE has been built for an air-water bubble column. Rather than solving a stand-alone PBE, the formulation used in this work consists of a cross-sectional averaged, two-fluid model that accounts for the characteristics of gas-liquid flow, and the population balance equation that predicts the bubble size distribution. A correlation given by Tomiyama39 has been applied for the drag force coefficient for the dispersed phase. Martinez-Bazan’s40-42 breakage kernel and a modified Prince and Blanch’s coalescence kernel43 have been chosen for bubble breakage and coalescence, respectively. In the iteration scheme, the flow calculation and the PBE solution are coupled and optimal solution is found by the LSM via minimizing the overall residue consisting of flow and PBE parts. The fact that no PBE operator splitting is required results in a simpler coupled iteration scheme. The predicted bubble size distribution as well as gas void fraction and bubble mean

7996


diameter have been compared to the measurements carried out by Krepper et al.,29 and good agreement has been found. This study should be regarded as the first attempts to elucidate the application of using the LSM to solve the combined multifluid and PBE model that is based on the physical processes occurring in the bubble column. A three-dimensional and transient extension of the study should be proposed for future investigations in order to capture the 3D and transient structure of bubble columns. Population Balance Formulation of the Population Balance Equation. The population balance equation is a statistical equation that describes the size distribution of the dispersed phase in a multiphase flow. Let us characterize a two-phase flow field by a distribution function f′(Vb,z,t), such that f′(Vb,z,t) dVb represents the number of particles per unit volume at time t, at position z, with size between Vb and Vb + dVb. In the PBE, the convective transport of bubbles and their coalescence and breakup as well as mass transfer and an expansion of the gaseous phase are taken into account. A derivation of the PBE is given by several authors19,44 in a similar form. In general, a volume-based PBE that includes coalescence, breakup, and source terms can be written as follows:

(2) rate of birth of bubbles of volume Vb due to breakup of larger bubbles, (3) rate of death of bubbles of volume Vb due to coalescence with other bubbles, (4) rate of birth of bubbles of volume Vb due to coalescence of smaller bubbles. In which, b′(Vb) is the breakup frequency, h′(Vb,Vb′) is the breakup kernel that represents the probability that a bubble of volume Vb′ splits up creating two bubbles of volumes Vb and Vb′ - Vb. The coalescence kernel c′(Vb,Vb′) represents the probability of coalescence per unit time so that two bubbles with volumes Vb - Vb′ and Vb′ to merge resulting in a bubble with the volume Vb. Vb,min and Vb,max are the lower and upper limits of the bubble size, respectively. Equation 3 can also be expressed in terms of the number density function using bubble length as an internal coordinate. The transient term and all the time dependences of each terms can be removed at steady-state. ∂(Vg f (ξ, z)) ) Bd(ξ, z) + Bb(ξ, z) + Cd(ξ, z) + Cb(ξ, z) ∂z (8) in which, Bd(ξ, z) ) -b(ξ)f(ξ, z)

D f′(Vb, z, t) ) B ′(Vb, z, t) + C ′(Vb, z, t) + S ′(Vb, z, t) Dt

Bb(ξ, z) )

∫

(1) 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 particles occurs instantaneously and the collision product may be treated as a single particle. The derivative of the left-hand side is given by dVb ∂ ∂ D ) + ∇ · Vg + Dt ∂t dt ∂Vb

(2)

in which Vg is the gas phase velocity (m/s). Since our investigation is limited to breakup and coalescence, the 1D population balance can be written as ∂(Vgf ′(Vb, z, t)) ∂f ′(Vb, z, t) + ) Bd′(Vb, z, t) + Bb′(Vb, z, t) + ∂t ∂z Cd′(Vb, z, t) + Cb′(Vb, z, t) (3)

Cb(ξ, z) )

ξ2 2

∫

h(ξ, ζ)b(ζ)f(ζ, z) dζ

(ξmax3 - ξ3)1/3

ξmin

c(ξ, ζ)f(ζ, z) dζ

(10) (11)

c((ξ3 - ζ3)1/3, ζ) f(ζ, z)f((ξ3 - ζ3)1/3, z) dζ 3 3 2/3 ξmin (ξ - ζ ) (12)

∫

ξ

where ξ and ζ represent the bubble diameter. The detail transformation from the volume-based PBE to the length-based PBE is explained by Marchisio et al.25 The moments of the distribution function of different orders can be defined as mk(z) )

∫ ξ f(ξ, z) dξ k

(13)

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ξ

(14)

(ii) the gas void fraction or gas holdup (volume of particles per unit volume)

in which, Bd′(Vb, z, t) ) -b′(Vb)f′(Vb, z, t) Bb′(Vb, z, t) )

∫

Vb,max

Vb

Cd′(Vb, z, t) ) -f′(Vb, z, t)

h′(Vb, V′b)b′(V′b)f′(V′b, z, t) dV′b

∫

Vb,max-Vb

Vb,min

(4) (5)

c′(Vb, Vb′ )f′(Vb′ , z, t) dVb′ (6)

Cb′(Vb, z, t) )

Cd(ξ, z) ) -f(ξ, z)

ξmax

ξ

(9)

1 2

∫

Vb

Vb,min

c′(Vb - V′b, V′b)f′(V′b, z, t)f′(Vb V′b, z, t) dV′b

R(z) )

∫ ( π6 ξ ) f (ξ, z) dξ 3

(15)

and (iii) the Sauter diameter

∫ ξ f(ξ, z) dξ ∫ ξ f(ξ, z) dξ 3

ξ32(z) )

2

(16)

The mean bubble volume can be computed by (7)

B d′, B b′, Cd′, and Cb′ have the following meaning: (1) rate of death of bubbles of volume Vb due to breakup,

∫ ( π6 ξ ) f(ξ, z) dξ ∫ f(ξ, z) dξ 3

R(z) ) Vb,mean(z) ) N(z)

(17)


The integrals in eqs 14-16 should be evaluated over all possible bubble lengths. Breakup Model. Martinez-Bazan et al.40,41 proposed a bubble breakup frequency model based on the turbulence stress and surface tension force balance which was validated with the authors’ experimental data for a turbulent jet. According to this model, the breakup frequency b(ξ) ([1/s]) is given by

β(εξ)2/3

b(ξ) ) Kg

σl - 12 Fl ξ

ξ

(18)

where the constant β ) 8.2 and Kg ) 0.25 were found experimentally for the case of air bubbles in water. The dimensionless daughter particle redistribution function h*(ξ*) ([-]) defined in terms of the PDF (particle distribution function) of the ratio of diameters ξ* (ξ* ) ξ/ζ, in which ξ is the diameter of one of the daughter bubbles and ζ is the parent bubble) was written as h*(ξ*) )

[ξ*2/3 - Λ5/3][(1 - ξ*3)2/9 - Λ5/3]

∫

* ξmax

* ξmin

[ξ*2/3 - Λ5/3][(1 - ξ*3)2/9 - Λ5/3] d(ξ*) (19)

in which Λ ) (ξmin/ξ)2/5. The redistribution function is thus given by h(ξ,ζ) ) h*(ξ*)/ζ ([1/m]). The peak of the 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: ε)

jggFl(1 - Rg) ) jgg Fl(1 - Rg)

(20)

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.21,22 Coalescence Model. A phenomenological model has been proposed by Prince and Blanch43 for the rates of bubble coalescence in turbulent gas-liquid dispersions. Jakobsen21 (p 884) proposed a modified kernel based on the Prince and Blanch’s work. In this kernel, the volume average coalescence frequency c(ξ,ζ) may be defined as the product of an effective swept volume rate hc(ξ,ζ) and the coalescence probability (efficiency) λc(ξ,ζ). c(ξ, ζ) ) hc(ξ, ζ)λc(ξ, ζ)

(21)

Kamp et al.,45 among others, suggested that microscopic closures can be formulated in line with the macroscopic population balance approach, thus the swept volume rate has been defined as hc,ij vjt,i2

π ) (ξi + ξj)2(v¯t,i2 + v¯t,j2)1/2 4

(22)

in which ) β(εξi) and β ≈ 2.0. A collision efficiency determines what fraction of bubble collisions lead to coalescence events. Only when the two bubbles are sufficiently close to each other so that the liquid film trapped 2/3

7997

between them starts draining, there is a possibility that the bubbles may be taken away from each other by a turbulent eddy. An expression for dimensionless coalescence efficiency proposed by Prince and Blanch43 is given as

[

λc,ij ) exp -

( ) rij3Fl 16σ

1/2

ε1/3 ln

h0 hf

rij2/3

]

(23)

Here h0 ) 10-3 cm is the initial film thickness with hf ) 10-6 cm is the critical film thickness where rupture occurs. The equivalent radius is given as rij )

(

1 1 1 + 2 rbi rbj

)

-1

(24)

where rbi and rbj are the radius of bubbles of class i and j. Finally, assuming that turbulence is the dominant mechanism determining the collision rate, we may define the coalescence frequency function as the collision frequency due to turbulence multiplied by the coalescence efficiency: cij ) hc,ijλc,ij

(25)

Bhole et al.46 suggests that the original Prince and Blanch collision rate model somewhat overestimates the collision rate; hence, a coalescence coefficient c0 (0 < c0 < 1) is added as an adjust factor. Cross-Sectional Averaged Two-Fluid Model Derivations. The cross-sectional averaged two-fluid model is developed by applying the cross-sectional area averaging operator on the two-phase transport equations (Chapter 3 of ref 21). At steady-state, by neglecting the interfacial mass transfer, the mass balance equation for each phase can be expressed as ∂ (R V ) ) 0 ∂z g g

(26)

∂ (R V ) ) 0 ∂z l l

(27)

The momentum equations for each phase are

(

)

∂pg ∂Vg ∂ ∂ (RgFgVgVg) + Rg ) -µg Rg + RgFgg + Mτg ∂z ∂z ∂z ∂z (28)

( )

∂pl ∂Vl ∂ ∂ (R F V V ) + Rl ) -µl R + RlFlg - Mτg ∂z l l l l ∂z ∂z l ∂z (29) where F is density, µ is viscosity, p is pressure, and g is the acceleration of gravity. Closures. The summation of the phase volumes must recover the whole domain, so that the fractions of the phases are coupled by the following relation: R g + Rl ) 1

(30)

It is assumed that the pressure in the liquid phase is equal to that in the gas phase: p g ) pl ) p

(31)

The gas phase interfacial momentum viscous stress term is

7998


3 CD Mτg ) - Fl Rg |Vg - Vl |(Vg - Vl) 4 dS

(32)

The Tomiyama drag coefficient correlation39 is used for obtaining the drag coefficient CD:

{ [

CD ) max min

]

16 48 8 Eo (1 + 0.15Rep0.687), , Rep Rep 3 Eo + 4

}

(33)

algebraic operators, and h1 and h2 are given vector-valued functions on the boundary Γ1 and Γ2, respectively. The definitions of problem operators L 1 and L 2, boundary operators B1 and B2, and source terms g1, g2, h1, and h2, are outlined in the Appendix. It is assumed that the system is well-posed. The least-squares formulation is to minimize the residual of a norm-equivalent functional. The norm-equivalent functional for the system is given by J (f1, f2) ) J1(f1) + J2(f2)

in which the dimensionless groups, the bubble Reynolds number, Rep, and the Eo¨tVo¨s number, Eo, are defined by

(40)

in which, g(Fl - Fg)dB2 Eo ) σl FlVgdB µl

Rep )

in Ω1 ∈ R2 on Γ1

(36) (37)

and a cross-sectional averaged two-fluid model (problem 2): L 2f2 ) g2,

in Ω2 ∈ R

B2f2 ) h2,

on Γ2

(41)

with the norm defined as

The least-squares method is a well established numerical method for solving a wide range of mathematical problems.33,47-49 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. Consider the system of linearized bubble column model that consists of a PBE (problem 1):

B1f1 ) h1,

1 1 | L f - g1 | Y1(Ω1)2 + |B1f1 - h1 | Y1(Γ1)2 2 22 2 1 1 J 2(f2) ) | L 2f2 - g2 | Y2(Ω2)2 + |B2f2 - h2 | Y2(Γ2)2 2 2

J1(f1) )

(35)

Least-Squares Spectral Element Method

L 1f1 ) g1,

{

(34)

1

(38) (39)

with L 1 and L 2 are first-order partial differential operators. Ω1 ∈ R2 and Ω2 ∈ R1 are bounded domains with piecewise smooth boundaries Γ1 and Γ2, respectively (Figure 2). The number of the space dimensions of the PBE problem is 2 as it includes both internal and external coordinates, while that of the twofluid model is 1 as it only has an external coordinate. f1 is the unknown density function n of x1 ) (ξ,z), and f2T ) (Vg,p,Vl,ug,ul) is a vector of five unknown functions of x2 ) z. g1 and g2 are given vector-valued functions. B1 and B2 are boundary

|•| Y(Ω)2 ) 〈•, •〉Y(Ω) )

∫ •• dΩ Ω

(42)

Following the conventional perturbation theory concepts, a system of equations is obtained by substituting f1 ) f1 + ε1v1 and f2 ) f2 + ε2v2 and take the partial derivative of J with respect to ε1 and ε2. The minimization statement is equivalent to the following: find f1 ∈ X1(Ω1) and f2 ∈ X2(Ω2) such that

{

limε1,2f0 limε1,2f0

∂ J (f1 + ε1v1, f2 + ε2v2) ) 0, ∀v1 ∈ X1(Ω1) ∂ε1 ∂ J (f1 + ε1v1, f2 + ε2v2) ) 0, ∀v2 ∈ X2(Ω2) ∂ε2 (43)

where X1(Ω1) and X2(Ω2) are the space of the admissible functions. v1 and v2 are arbitrary trial functions, and ε1 and ε2 are small perturbations. Consequently, the necessary condition can be written as follows: Find f1 ∈ X1(Ω1) and f2 ∈ X2(Ω2) such that

{

A1(f1, v1) ) F1(v1) ∀v1 ∈ X1(Ω1) A2(f2, v2) ) F2(v2) ∀v2 ∈ X2(Ω2)

with

{

A1(f1, v1) F1(v1) A2(f2, v2) F2(v2)

) ) ) )

(44)

〈L1f1, L1v1〉Y1(Ω1) + 〈B1f1, B1v1〉Y1(Γ1) 〈g1, L1v1〉Y1(Ω1) + 〈h1, B1v1〉Y1(Γ1) 〈L2f2, L2v2〉Y2(Ω2) + 〈B2f2, B2v2〉Y2(Γ2) 〈g2, L2v2〉Y2(Ω2) + 〈h2, B2v2〉Y2(Γ2) (45)

where A1: X1 × X1 f R and A2: X2 × X2 f R are symmetric, continuous bilinear forms, and F1: X1 f R and F2: X2 f R are continuous linear forms. Spectral Discretization in R2. The discretization statement ˆ ), consists in searching the solution in a reduced subspace XN(Ω i.e. ˆ ) ) X ∩ PN(Ω ˆ ) ⊂ X(Ω ˆ) ˆfN(xˆ) ∈ XN(Ω

Figure 2. Physical domain of the steady-state, one-dimensional bubble column model with the (a) domain of the PBE and (b) domain of the twofluid equations.

(46)

in which, the accent indicates that we discretize the solution on a two-dimensional reference domain. We can define the two-dimensional high-order polynomial ˆ ) by using the one-dimensional nodal basis space PN(Ω


(Lagrangian interpolant) based on a tensor product Gauss-Lobatto (Legendre) grid; see Figure 3. Note that ˆ )) ) (Nξ + 1) × (Nz + 1) in dim(PN(Ω

R2

[A1]ij ) A1(Φj, Φi) ) 〈L1Φj, L1Φi〉Y(Ω1) + 〈B1Φj, B1Φi〉Y(Ω1) (54) [F1]i ) F(Φi) ) 〈[g1]i, L1Φi〉Y(Γ1) + 〈[h1]i, B1Φi〉Y(Γ1)

(47)

(55)

Therefore, ˆfN(ξˆ ,zˆ) can be expressed in terms of linear combinations of basis functions: Nξ+1 Nz+1

ˆfN(ξˆ , zˆ) )

∑ ∑ r)1

1

(48)

s)1

-1

K+1 L+1

1

ˆf(ξˆ , zˆ) dξˆ dzˆ ≈ -1

∑ ∑ w w ˆf

k l kl

(49)

k)1 l)1

where wk, wl and ξˆ k, zˆl are the weights and points of quadrature.50 If the number of 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 ∂(φr(ξˆ )φs(zˆ)) ∂ξˆ

|

∂(φr(ξˆ )φs(zˆ)) ∂zˆ

|

ξp,zq

ξp,zq

) φr′(ξˆ p)φs(zˆq) ) [dξˆ ]pr[bzˆ]qs

(50)

) φr(ξˆ p)φs′(zˆq) ) [bξˆ ]pr[dzˆ]qs

(51)

in which [bzˆ] and [bξˆ ] are the one-dimensional mass matrices in ξˆ and zˆ, respectively, while [dξˆ ] and [dzˆ] are the onedimensional derivatives matrices in ξˆ and zˆ. Using these approximations in eq 45, the following matrix systems are obtained for the PBE problem and the flow problem, respectively. A1f1 ) F1

(52)

A2f2 ) F2

(53)

where the matrix A1 ∈ R(Nξ×Nz )×(Nξ×Nz ), vectors f1,F1 ∈ R(Nξ×Nz ), t t t and matrix A2 ∈ R(Nz ×5)×(Nz ×5), vectors f2,F2 ∈ R(Nz ×5), are defined like t

t

t

[f1]i ) f([x1]i)

(56)

and ˆfrsφr(ξˆ )φs(zˆ)

where ˆfrs is the basis coefficient associated with the basis function φr(ξˆ )φs(zˆ); it follows that the basis coefficient ˆfrs is equal to the value of the discrete solution ˆfN(ξ,z) at the GaussLegendre-Lobatto (GLL) points (or node) (ξˆ p,zˆq), i.e. ˆf(ξˆ p,zˆq) ) ˆfN(ξˆ p,zˆq). Due to this property, we also say that we are using a nodal basis. The integration of the function is evaluated by the Gaussian quadrature rule. If order of the Gauss quadrature points are set to be equal to the order of the approximating polynomials, then

∫ ∫

7999

t

t

t

[A2]ij ) A2(Φj, Φi) ) 〈L2Φj, L2Φi〉Y(Ω2) + 〈B2Φj, B2Φi〉Y(Ω2) (57) [F2]i ) F(Φi) ) 〈[g2]i, L2Φi〉Y(Γ2) + 〈[h2]i, B2Φi〉Y(Γ2) (58) [f2]i ) f2([x2]i)

(59)

Model Implementation There are several possibilities of combining the PBE computation with CFD code. As mentioned in the introduction, various methods are available to solve the PBE numerically. Among these methods, the method of classes is not practical in coupling with flow calculation as it needs at least 12-15 classes to achieve good accuracy. QMOM has been reported as providing efficient and accurate results with the CFD when the number of moments is low. The system size is proportional to the number of moments that needs to be solved. The DQMOM is based on the idea of tracking directly the variables appearing in the quadrature approximation instead of the moments of bubble size distribution. The parallel parent and daughter classes (PPDC) method that is based on modified CM is shown to perform as well as QMOM and DQMOM, using only a few classes (2-3).31 The operator splitting of the PBE is applied such that its convective part is solved together with the continuity and momentum equations, while the breakage/coalescence parts are solved after the flow iteration loop is converged (see Figure 4). Such an algorithm usually needs multiple iterative loopssone loop for the calculation of flow equations combined with the convective term of the PBE (LOOP 2), which contains a nested loop for the convergence of the pressure correction equation (LOOP 1) and another loop for calculation of the nonlinear PBE with the breakage and coalescence terms (LOOP 3). Figure 5 illustrates the algorithm for solving the PBE with multifluid model that is used in the bubble column problem of this work. It starts with two separated guessed values of unknowns f *1 (f *) and f *2 (v*,v c *,P*). d Since the information of the density function is known, moments of any order as well as the void fraction (R*) and the Sauter diameter (ξ32) can be calculated. Together with the guessed values given, the PBE and multifluid linearized opera-

Figure 3. Four of the forty nine two-dimensional basis functions Φr,s(ξˆ ,zˆ) of order Nξ ) Nz ) 6, (a) Φ1,1(ξˆ ,zˆ), (b) Φ3,3(ξˆ ,zˆ), (c) Φ5,5(ξˆ ,zˆ), (d) Φ7,7(ξˆ ,zˆ).

8000


Figure 4. Conventional algorithm for solving the PBE coupled with the two-fluid model.

tors L 1 and L 2 with their corresponding sources terms g1 and g2 are calculated. Computation of matrices A*1 and F*1 , as well as A*2 and F*2 are carried out by the least-squares method matrix constructor based on each individual problem operator. The next step is the assembly procedure. The overall stiffness matrix A* and source vector F* have the following structures: A* )

[

A*1 0 0 A*2

]

F* )

[ ] F*1 F*2

(60)

in which A* is the matrix constructed by diagonalizing the stiffness matrices of the PBE problem A*1 and the multifluid problem A*2 . Note that the prerequest of the above assembly is that both unknowns f *1 and f *2 can be mapped to the L 1(Ω1) and L 2(Ω2) function spaces via problem operators L 1 and L 2, respectively. In comparison to the convectional approach of coupling the PBE and the flow equations, we found that solving these equations by using the least-squares method is simpler and more systematic. Rather than conducting a decoupled iterative process, the least-squares method solves the flow equations and the population balance equation within the same framework. The least-squares approach does not split the PBE operator into a convective part, as well as breakage and coalescence parts so that they are solved sequentially. Therefore, no multiple iterative loops are required. Instead of solving the moments, the leastsquares method solves the bubble size distribution function directly. Moments of any order can be retrieved in postcalculations by eq 13. The number of the moments does not influence the size of the system, which is the case in QMOM. The way that the least-squares method trying to minimize both sets of

Figure 5. Iteration algorithm for solving the PBE coupled with two-fluid model by using the LSM.

equations at the same time and searches for the optimal solution to the overall system corresponds to its intrinsic formulation, eq 40. Experimental Setup Figure 6 shows the geometrical construction of the variable gas injection with six almost logarithmic injection modules distributed over the total height. Each module (Figure 7) consists of three chambers. Gas is injected through orifices in the pipe wall. This gas injection offers the advantage that the two-phase flow can rise smoothly to the measurement plane, without being influenced by the feeder within the tube in other height positions. Two of the three chambers (the uppermost and the lowest) are equipped with 72 × 1 mm orifices. The middle chambers have 32 × 4 mm orifices. For rotation-symmetric gas injection, all orifices per chamber are equally distributed over the circumference of the pipe. The chambers are connected with a gas injection pipe, and the compressed air system can be operated separately. The supply of the liquid phase is done from the bottom of the test section by means of an isolating valve and a 90° bend. For these experiments, the measurement plane was always situated at the upper end of the test section (Figure 6). A twolevel low temperature wire-mesh sensor with 64 × 64 wires was used, which corresponds to a lateral distance of wires of 3 mm. Special attention was paid to the pressure and temperature boundary conditions during the experiments. The pressure was kept constant at 0.25 MPa at the location of the activated gas injection to represent the situation of an evaluational flow by using different gas injections. The water temperature was 30 °C ( 1 K.


8001

Table 1. Absolute and Relative Heights at the Test Section Variable Gas Injection injection chamber

position of height

injection length [mm]

Ltube/dtube ratio

1 1 2 2 3

A C D F G

221 335 494 608 1438

1.1 1.7 2.5 3.1 7.4

A detailed explanation on the experimental setup and test procedures is found in refs 51 and 52. Results and Discussion

Figure 6. Vertical test section of the TOPFLOW facility with variable gas injection system. Reprinted with permission from the work of Prasser et al.51 Copyright 2007 Elsevier.

Figure 7. Injection module of the variable gas injection. Reprinted with permission from the work of Prasser et al.51 Copyright 2007 Elsevier.

During each experiment, water is fed through the vertical test section. The control of the water mass flow was done using three parallel switched regulating valves and the associated automatic control loops. For the experiments, water mass flow rates between 1.2 and 48 kg/s were applied. Compressed air was injected into the test section circuit in order to approximately achieve the necessary pressure at downstream the wire-mesh sensor. Afterward the volumetric airflow can be reduced to its set value. In the following, the pressure value can be adjusted to the set value by varying the blow down flow rate. After setting all technical boundary conditions, the measurement can be conducted. For this, the wire-mesh sensor electronics have to be prepared for the measurement and the flow structure has to be checked by means of online visualization. After fulfilling these conditions, the measurement is done offline. In order to be able to determine characteristic flow parameters, such as void fraction profile, velocity profile, and bubble size distribution from the measured data, the time for each single measurement was fixed for all experiments of this test series at 10 s. The measurement frequency is kept at 2.5 kHz. Directly after the measurement, the quality of the data can be checked by visualization on the measurement computer.

In the simulations, instead of the whole column, the first part of it has been considered as the axial domain in order to reduce the influence of the pressure on the bubble sizes. We’ve chosen the tube length from the injection point (z ) 0 m) to point G (z ) 1.438 m). Hence experimental data on the bubble size distribution that is evolved from the injection by the 1 mm orifice at six axial positions (A, C, D, F, and G) are available for comparison. Table 1 lists the vertical distances between the individual gas injections and the first measurement plane of the wire-mesh sensor located in direction of flow. The spectral elements have been used to make the computation more effectively. Relatively intensive changes of the flow parameters and the bubble concentration are observed at the beginning section of the tube; the bubble size distribution function also presents stronger variations at the size range from ca. 4 to 8 mm. Therefore h-refinements are conducted by dividing the physical domain into several subelements in both space and property dimensions and assigning these subelements with the proper number of Gauss quadrature points (Figure 8). The superficial gas and liquid velocities of the experimental set are 0.0151 and 0.405 m/s, respectively. These operating conditions, the system geometry, and the physical data given in Tables 1 and 2 are associated with a pipe Reynolds number on the order of Re ∼ 105; hence, the flow is turbulent. As the data reveals that the breakage is dominant in these experiments, coefficients associated with breakage (Kg) and coalescence (c0) are set to be 6 × 10-3 and 0.0, respectively. The physical properties of the materials used in the simulations are listed in Table 2. Figure 9 shows 3D plot and the contour plot of numerical results of the bubble size distribution n as the function of z and d. It can be seen that the main peak of the bubble number is decreasing along the tube indicating the breakage of bubbles of this size lowered its concentration and results in emergence of a higher number of smaller bubbles. The experimental bubble size distribution for the specified superficial gas and liquid velocities is given in Figure 10 with comparisons from the LSM simulation. The numerical results at the GLL points in axial direction z have been interpolated to the axial experimental measuring points (A, C, D, F, and G). It is seen from the data in Figure 10 that the experimental bubble number density maintains almost the same shape, while it shifts toward the smaller bubbles; there is no accumulations of bubbles of certain sizes. The numerical result presents a relatively steep increase and decrease in number density that are obtained for smaller and larger bubbles, respectively. We can see that the model predicted values correspond to the experimental results at first several points while deviations grow along with the axial coordinate. The numerical results predict a more flattened distribution function than the experiments. Higher numbers of the small

8002


Figure 8. GLL-GL grid with eight spectral elements of order 3 (a), 4 (b), and 5 (c) in the 2D physical domain of the bubble column model. The grid scheme (part b) has been chosen, and the total number of the GLL points in z is Nzt ) 4 × 2 - 1 ) 7. The total number of the GL points in ξ is Ntξ ) 4 × 4 ) 16. Table 2. Physical Properties of the Materials physical property

value

unit

air viscosity water viscosity air density water density surface tension

1.73 × 10-5 1.00 × 10-3 1.2041 1.00 × 103 0.0728

Pa · m Pa · m kg/m3 kg/m3 N/m

bubbles (

On the Coupled Solution of a Combined Population Balance Model

Recommend Documents