Scholarly article on topic 'Least Squares Higher Order Method for the Solution of a Combined Multifluid-Population Balance Model: Modeling and Implementation Issues'

Least Squares Higher Order Method for the Solution of a Combined Multifluid-Population Balance Model: Modeling and Implementation Issues Academic research paper on "Chemical engineering"

Share paper
Academic journal
Procedia Engineering
OECD Field of science
{Multifluid-model / "population balance" / "density function" / "breakage rate" / "coalescence rate" / "growth velocity"}

Abstract of research paper on Chemical engineering, author of scientific article — Z. Borka, H.A. Jakobsen

Abstract In this paper we attempt to highlight the main challenges of the implementation of the presented combined Eulerian multifluid-Population balance model for bubbly flows as applied to simulate the behavior of vertical bubble column reactor. The 1–D steady-state model equations for the bubble column with water as the liquid continuous phase and air bubbles of diameter as a dispersed phase, constitutes the liquid continuity and momentum equations and the dispersed phase population balance equation and the dispersed phase momentum equation. In addition, the ideal gas law, two moments of the distribution (void and SMD), and a set of constitutive equations are required.

Academic research paper on topic "Least Squares Higher Order Method for the Solution of a Combined Multifluid-Population Balance Model: Modeling and Implementation Issues"

Available online at

Procedia Engineering


Procedía Engineering 42 (2012) 1121 - 1 132

20th International Congress of Chemical and Process Engineering CHISA 2012 25 - 29 August 2012, Prague, Czech Republic

Least squares higher order method for the solution of a combined multifluid-population balance model: Modeling and

implementation issues

Department of Chemical Engineering, Norwegian University of Science and Technology (NTNU), N — 7491 Trondheim, Norway

In this paper we attempt to highlight the main challenges of the implementation of the presented combined Eulerian multifluid-Population balance model for bubbly flows as applied to simulate the behavior of vertical bubble column reactor. The 1-D steady-state model equations for the bubble column with water as the liquid continuous phase and air bubbles of diameter E, as a dispersed phase, constitutes the liquid continuity and momentum equations and the dispersed phase population balance equation and the dispersed phase momentum equation. In addition, the ideal gas law, two moments of the distribution (void and SMD), and a set of constitutive equations are required.

© 2012 Published by Elsevier Ltd. Selection under responsibility of the Congress Scientific Committee (Petr Kluson)

Keywords: Multifluid-model; population balance; density function; breakage rate; coalescence rate; growth velocity

1. Introduction

Bubble columns [1] are used in many different industries, such as chemical, biochemical, petrochemical and pharmaceutical. Common usage among others includes hydrogenation, chlorination, oxidation, biotechnological applications, cleaning of chemical gas. Bubble column reactors are also employed in the processes of partial oxidation of ethylene to acetaldehyde, oxidation of p-xylene to

* Corresponding author. Tel.: +47 735 941 35; fax: +47 735 940 80. E-mail address:

Z. Borka a*, H. A. Jakobsen


1877-7058 © 2012 Published by Elsevier Ltd. doi:10.1016/j.proeng.2012.07.504

dimethylterephthalate, Fischer-Tropsch synthesis, synthesis of methanol, hydrocarbons, phosgene hydrolysis, ozonization of waste water and biological waste water treatment.

As a gas-liquid contactor, bubble columns are often preferred because of their simplicity of operation and low operating costs. Parameters such as local void fraction, gas-liquid interfacial area, mass transfer coefficients between phases, dispersion coefficients, the distribution of liquid velocity, turbulent kinetic energy etc. have been studied for decades for the design and scale-up of these reactors. During the last two-three decades, the possibility to utilize computational fluid dynamics (CFD) in the modeling of multiphase reactors has been of great interest in the reaction engineering community.

In general, there are two main approaches for computing multiphase flows, the Eulerian-Eulerian model and the Eulerian-Lagrangian model [2]. The Eulerian-Lagrangian model considers the liquid phase to be continuous, while the dispersed phase is represented as discrete particles. The resulting computationally demanding direct numerical simulations of the particle interactions require detailed specification of interface properties to predict particle deformation and explicit rules for coalescence and breakage.

In practice, the Eulerian-Eulerian approach is typically used to study large-scale flow structures and dense dispersed systems. It has lower computational demand compared to the Eulerian-Lagrangian model. A disadvantage of the model is the need for closure laws for the inter-phase transport of mass and momentum. Assuming no chemical reaction or phase changes and considering non soluble gas-liquid flows the model has to account for the mass exchange between bubbles of different sizes by bubble breakage and coalescence and for various interfacial forces which cause momentum transfer between the bubbles and the ambient liquid. In particular, the lift and drag forces have been shown to have the greatest influence on the bubble size distribution along a particular flow.


Latin letters c [1/s]

d32 [m]

Eo [-] fd [kg/m3m] fdrag [N/m3m] fw [-] hi [m] hb [1/m] hf [m] j [kg/m3s] jm [kg.m/s/m3s] Mw [kg/mol] p [Pa]

coalescence rate Sauter mean diameter Eotvos number mass density function drag force

wall friction coefficient initial film thickness daughter size redistribution rapture thickness of the film mass flux momentum flux air molecular weight pressure

R [J/K/mol] universal gas constant vd [m/s] dispersed phase (gas) velocity vl [m/s] liquid velocity vslip [m/s] slip velocity v^ [m/s] growth velocity Greek letters

ad [-] gas volume fraction

Q [m] bubble diameter

e [m2/s3] turbulent energy dissipation

r [kg/m/s]

liquid dynamic viscosity coalescence rate

Pd [kg/ m3] Pi [kg/ m3]

dispersed phase density dispersed phase density

A number of different approaches have been proposed to account for the different bubble sizes in conjunction with the classic two-phase Eulerian flow equations. In an attempt to solve the convergency problems associated with the large number of equations resulting from the direct extension of the two-phase equations to polydispersed flows, the homogeneous multiple size group (MUSIG) [3] introduced in ANSYS CFX discretizes the bubble size distribution profile into several classes based on bubble diameter, i.e. it divides the dispersed phase into M non-overlapping size fractions, while it solves one momentum equation for all bubble classes. A further extension of this approach is the inhomogeneous MUSIG model [4] that uses several dispersed phase velocity groups, each with its own bubble size fractions.

The particular population balance framework applied in this work is derived from a kinetic theory approach averaging the equations over the particle velocity space, but not over the size property space as in the conventional population balance formulation described by Ramkrishna [5]. The result is a set of non-linear integro-differential equations expressed in terms of the size dependent number density function and particle velocities for the bubble phase fields. The dispersed phase model equations are combined with the conventional multi-fluid model equations for the continuous phase. The joint set of model equations constitutes a detailed and computationally feasible approach for the modeling of bubbly flows. Given the associated mathematical complexity and the number of different closure models currently available, a robust yet flexible numerical method is desirable. In this regard, we found the finite element based method of the least-square kind used in this work to be highly adequate [6].

Previous work in our group denotes a combined multifluid-population balance modeling approach for bubble columns. Zhengjie [7] combined the two fluid-Population balance equation model using an average bubble velocity for all bubble sizes. Patruno [8] extended the modeling framework allowing different velocities for different fluid particle sizes. The model was employed for droplets in mist flows in scrubbers. Nayak et al [9] applied the combined multifluid-Population balance model to bubbly flows. Sporleder [10] used a similar model and evaluated the use of different internal coordinates. In the present

work, the model formulation, the sub-models or closures used and model implementation issues are outlined in detail.

Models for both bubble mass and volume based population balances with various closure laws have been proposed. This work presents a mass based approach pointing out model modifications and mathematical constraints required for the closure laws with regard to fulfill the physical constraints of the modeled system (e.g., the mass conservation property) separating the modeling task from the numerical and implementation issues.

2. Model description

2.1. Governing equations

The main governing equations of the model outlined in Table 1 are the microscopic mass and momentum conservation laws for both the continuous (liquid) and dispersed (gas) phases. The presented model assumes steady-state operation conditions relevant in industrial processes. It uses one domain variable z in the physical dimension corresponding to the vertical height of the liquid column in the reactor. Although the effect of bubble deformation on the flow [11] has been studied to some extent, most of the work in the literature assumes the collisions of non-deformable spherical bubbles. The bubble diameter § represent the variable of the property space used in the population balance equation.

o cP^) 8 .-o-be^.J

8 0 o° feo ■'ViT)"1

Fig. 1. Vertical bubble column

Overall, the model unknowns are the continuous phase velocity vi(z), pressure p(z), mass distribution function fd(z, §), dispersed phase velocity vd(z, §), dispersed phase volume fraction ad(z), growth velocity v§(z, §), mass flux j(z,§) and momentum flux jm(z, §).

The model variables are normalized and made dimensionless using the values of Dirichlet boundary conditions as summarized in Table 2.

The breakage and coalescence operators in the literature are conventionally defined with |min = 0. However, the existence of bubbles with zero diameters is not physical and will make the dispersed phase momentum balance singular as the balance contains the distribution density function. Hence, the numerical implementation may require special treatment for zero minimum values. Furthermore, some of the proposed breakage rate closures assume a certain critical particle diameter for which the smaller particles can not undergo breakage. In the present model, a non-zero minimum bubble size |min is used. To retain the mass conservation properties of the coalescence operators, it is necessary to modify the integral limits as shown in Table 3. A Heaviside function H(x) was introduced in the coalescence operator terms to restrict the range of physically possible mother particles. The Heaviside function H(x) may be defined as:

H (x) = 0 if x < 0 else 1 (1)

The coalescence source and sink terms can then be specified as: C — Hh- 3/^^TW (h r(hh)1/3 c(z,(?-C3)1/3,C) ft(z,h3-Q3)1/3) f (zO) o (2) cs =H h V2V ))x 2 h (h-c3)2'3 Vh) - no V(C)Pd (z)d0

cd =-H (v V (hx) - V h) -h)x fd (zh)f--|3)1/3 c( zh,C)f[z do (3)

v ! Ji~ V (0)pd (z)

The overwhelming majority of the breakage and coalescence closures were developed with binary collisions in mind. The breakage rate b(h) gives the likelihood of a particle of diameter h breaking up to smaller daughter particles. A particle is assumed to break into exactly two smaller daughter particles, and the daughter sizes are specified by the daughter size redistribution function hb. Of a particular interest are the 0th and 3rd moments of the redistribution function which for a fixed number of daughter particles v and volumes of the parent V(h) and daughter V(ra) particles can be used to check number and mass conservation properties:

\^hh(ah)da=v (4) iim"h(>h)V (a) da = V (h) (5)

»min »min

Z. Borka and H. A. Jakobsen /Procedia Engineering 42 (2012) 1121 - 1132 Table 1. Model equations and constitutive relationships

Continuity equation for continuous(liquid) phase

& pv ) = 0

Dispersed phase volume fraction Jçmin p

Momentum equation for continuous phase

dv, dv, dp (7)

(1 -&d ) Pi^r + (1 -&d ) PVi-T + (1 -&=

-(1 -&d) Pig + & Ptg + Fi

Growth velocity

v,( z,4) = -

3Pd ( z )

v ( z,4)

Momentum equation for dispersed(gas) phase

J- ( fd ( z, £)vr (z, £)vr (z,D ) = j ( z, £)

fd (z,4) dp(z)

+ fd (z,4)g + fdrav (z,4)

Pd ( z ) dz

Population balance equation

I" (vrfd ) + ^ Of ) = ( BB + BD + CB + CD ) fd

dz dç

Dispersed phase density

PjL = p 0 0 Pd P

Mass flux

j = fdvg Momentum flux

jm = fdV4Vd Volume fractions

a, + a„ = 1

Table 2. Dimensionless variables and initial and boundary conditions

z = — L

P= 0 P

( Po +Pg Az )Mad

4 (17)

v, (18)

v, = —

(1 -&0) J

v, = v, ——

d d .0

v, = v,

l\z=0 1

P|z=L = P 0

(23) LU,„= 0 (27)

f I = f

Jd |z=0 J d

v,L=0 = ^

vd\4=4m= vi + vsüp (25)

v4\ = 0

■mit (26)

2.2. Initial and boundary conditions

The type of gas spargers used is known to have a strong influence on the hydrodynamic properties of the flow. It induces a hydrodynamically active and highly transient flow pattern close to the inlet of the reactor which is not accounted for in the presented Eulerian equations. Instead, the influence of various gas spargers on the flow is simulated by the population balance using a measured or even an assumed initial mass based bubble size distribution profile fdimt(|).

Table 3. Closure and constitutive relations

Slip velocity [12]

vsi,P =( + v2 )'

Drag force

- 3 P&^^TT- V (Z't) - V (z)| (h (z,4) - vi(z))(31) 4 4 Pd (z)

Drag coefficient [13]

{ ■ [ 16 nxii HP 0.68^ 48 ] 8 Eo } (32) CDT = max\min[-(1 + 0.15Rep ),-],--0—}


Re/ 3 E0 + 4

2ct, gi

Vb2 i T

Reynolds number

Rel ( z ) = PVD

Dissipation rate = jdg

Wall friction factor[14]

f ( z ) = (0.79 ln[ Rei ( z )]-1.64) -

Bubble swarm effect [15,16]

Cd = (1 -a,)-nCD

Gas flux

Jd = Vdad

Mass density function

fd = fpdV (£)

Volume conservation

a +ad= 1

Coalescence source term

Eötvös number

g ( Pi-Px

Sauter mean diameter

d32 = VfMI

Breakage source term

bb = V&f-'h (£,mz,o

Breakage sink term

BD =-b( z,Z)fd ( z ,£)

f ( z V (Ç)

Coalescence sink term

Mm*-i3)"3 Cd =-fd (z,#)f c(z,#, Z)

f (z,Q d ' V (C)Pd ( z )

#2V(#) r«3-iL)"3 c(z,(4 -Z3)1/3,Z) fd(z,(4 -Z3)1/3) f(z,0

(43 - ç3)2

V(i) - v (Z ) v (O Pd ( z )

The liquid velocity vimit at the inlet is usually not known, but the liquid flux or superficial velocity is generally specified. The bubbles at the inlet are normally considered to have reached their terminal velocities For the water-gas systems, a slip velocity model [12] can be used to obtain the initial dispersed phase velocity vd(z,§) profile at the inlet. Further along the reactor, the dispersed phase velocity is calculated from the momentum equation. The dispersed phase velocity for the smallest bubbles with diameter |min is set to the slip velocity. The velocity of the smallest bubbles is then equal to the ambient liquid velocity. The distribution function fd for bubbles with |min ^ 0 is put equal to zero inside the reactor. At the top outlet of the reactor, the pressure p0 equals the atmospheric pressure of 1 atm.

For the hyperbolic equations, initial conditions are required. To preserve mass conservation, the bubble size growth velocity v^(z,§) must be zero at the inlet boundary of the property domain. In addition, the bubbles cannot grow outside the domain as they rise up in the column due to the pressure drop. This put a requirement to the constitutive relation 12 at the property outlet, the growth velocity property must approach zero as ^^ 0 thus v(z, §max) = 0.

3. Numerical solution

3.1. Least-Square method

The resulting set of integro-differential equations is solved using the least square spectral element method [17]. It is a high-order finite element method which uses the roots of Gauss-Lobatto-Legendre (GLL) polynomials as collocation points to minimize the spectral error. The minimization problem is expressed as

3.2. Implementation

The implementation follows a general finite element code pattern. The actual model definition is based on a data driven approach with emphasis on simplicity. Loose coupling between the individual components and unit testing of numerical results proved to be highly useful during the iterative development. The element based code allowed for trivial parallelization.

4. Results and Discussion

The parameters for the model were chosen with regard to available experimental data. Yao [18] measured radial profiles at different locations along the column. Jakobsen [19] provides data on bubble distributions, velocities and gas holdup. The available radial profiles were cross-sectional averaged for comparison with model simulations. No data was found regarding the experimental errors.

The simulations were carried out on an air-water vertical bubble column. The vertical height of the column was H = 4.25 m, with inner diameter D = 29 cm. The computational grid used for the simulations consisted of 16 points in z using 5 spectral Gauss-Lobatto-Legendre elements of order 4. In we used 1 such element of order 20 so that the integral operators for breakage and coalescence were confined to one

e=0 q=0

The system of linear equations for the approximate solution is given by

Af = F A = LWWL F = LTWg

element for all bubble diameters The breakage and coalescence term present a Fredholm and Volterra integral operator of a 2nd kind which can be solved for iteratively [20].

The breakage rate and breakage redistribution model of Coulaloglou [21] was used in the simulations with parameter values k1 = 0.336 and k2 = 0.106. For the coalescence we used the model presented by Prince [22]. The coalescence rate kernel c in this model is assumed to be a product of collision rate hc and coalescence efficiency Xc.

Normal ambient pressure p0 = 1 atm and temperature T = 293K are assumed at the column outlet. At the inlet, the liquid volume flux ji0 is 1 cm/s, the gas volume flux jd0 is 4 cm/s. The inlet void fraction provided by Jakobsen [19] was ad0 = 0.133. The initial bubble mass distribution density was assumed to have a log normal distribution with mean value | = -0.165 and variance ct = 0.364 corresponding to an average bubble diameter value |mean= 5.03 mm as observed by Kim et al [23]. The low end convergence limit of the conjugate gradient solver was 10-14. The relative tolerance of Picard iterations was set to 10-3.

Two sets of results were simulated. In the first calculations we used the drag coefficient provided by Tomiyama (Eq. 34) for slightly-contaminated system correlated to tap water. In the second calculations we considered the swarm effect (Eq. 36) on the drag coefficient as depicted by Roghair et al [15].

Figure 2. Initial volume distribution

Figure 3. Pressure

Figure 2 shows the initial profile for the bubble distribution function fd. The dispersed phase density pd and the pressure p depicted on figure 3 follow a similar descending profile along the vertical column given that the two quantities are related according to the ideal gas law.

Shown on Figure 6 is the average gas velocity with no swarm effect. Simulations considering the bubble swarm effect yielded slightly higher values for the gas velocity profile which didn't seem to agree with the measurements.

Figure 4. Liquid phase velocity

Figure 5. Pressure

The calculated values of the Sauter mean diameter on Figure 7 correspond to the initial bubble volume distribution profile from Figure 2 where the mean diameter |mean = 5.03 mm is in agreement with the calculated values. The four point experimental values for mean bubble diameter seem to average around 3 mm. No data on the distribution profile from the measurements was found, therefore a lognormal distribution from the measurements of Kim et al [23] was assumed without further modifications.

Figure 6. Average dispersed phase velocity 5. Conclusion

Figure 7. Sauter mean diameter

A gas-liquid model for the steady-state bubbly flow in a vertical bubble column has been presented. By solving for the distribution profile directly, the population balance approach with both breakage and coalescence closures resulted in more detailed results from the 1-D model than the methods considering only the moments of the distribution.

Numerical simulations were carried out by using the least squares spectral element method for air-water system. Predictions from the model were in satisfactory agreement with the available experimental data. Data regarding the shape of the initial bubble distribution at the inlet of the reactor were taken from a separate source, thereby resulting in a higher mean bubble diameter compared to the data from the first source.

The effect of breakage and coalescence on the bubble distribution profile along the reactor was observed. The breakage operators shift the distribution towards lower bubble sizes, coalescence has the

opposite effect. Depending on the breakage closure used the shift could occur up to the smallest diameters §mm in the domain or only up to some critical diameter given by the particular model. Further work concerning the validity of various breakage and coalescence models is needed.

A mass conservative formulation of the breakage and coalescence terms was presented. In our implementation we found that the element based approach for the solution of integral terms for breakage and coalescence significantly increased the complexity of implementation and computational effort without improving the convergence of the iterations. While the breakage terms were solved very efficiently by precalculations, the implementation of the nonlinear coalescence terms proved to be troublesome and require further attention.


The Ph.D. fellowship (Z. Borka) financed by the Research Council of Norway through the GASSMAKS program is gratefully appreciated.


[1] Delnoij E. Fluid dynamics of gas-liquid bubble columns. A theoretical and experimental study. Ph.D. thesis, University of Twente, Enschede, Netherlands 1998.

[2] Jakobsen HA. Chemical Reactor Modeling. Multiphase Reactive Flows, 1st Edition, Springer; 2008.

[3] Lo S. Application of the MUSIG model to bubbly flows. AEAT-1096, AEA Technology 1996.

[4] Krepper E, Lucas D, Frank T, Prasser HM, Zwart PJ. The inhomogeneous MUSIG model for the simulation of polydispersed flows. NuclEngDes 2008;238:1690-1702.

[5] Ramkrishna D. Population Balances, Theory and Application to Particulate Systems in Engineering. Academic Press, 2000.

[6] Dorao CA. High order methods for the solution of the population balance equation with applications to bubbly flows. Ph.D. thesis, Norwegian University of Science and Technology 2006.

[7] Zhu Z. The least-squares spectral element method solution of the gas-liquid multi-fluid model coupled with the population balance equation. Ph.D. thesis, Norwegian University of Science and Technology 2009.

[8] Patruno LE. Experimental and numerical investigations of liquid fragmentation and droplet generation for gas processing at high pressures. Ph.D. thesis, Norwegian University of science and technology 2010.

[9] Nayak AK, Borka Z, Patruno LE, Sporleder F, Dorao CA, Jakobsen HA. A combined multifluid-population balance model for vertical gas-liquid bubble-driven flows considering bubble column operating conditions. IndEng Chem Res 2011;50:1786-1798.

[10] Sporleder F, Dorao CA, Jakobsen HA. Model based on population balance for the simulation of bubble columns using methods of the least-square type. Chem Eng Sci 2011;66:3133-3144.

[11] Leon-Becerril E. Effect of bubble deformation on stability and mixing in bubble columns. Chemical Engineering Science 2002;57:3283-3297.

[12] Fan LS, Tsuchiya K. Bubble wake dynamics in liquids and liquid-solid suspensions. Butterworth-Heinemann, 1990.

[13] Tomiyama A. Struggle with computational bubble dynamics. Multiphase Science and Technology 1998;10:369-405.

[14] F. M. White. Fluid mechanics. McGraw, 1979.

[15] Roghair I, Lau YM, Deen NG, Slagter HM, Baltussen MW, Sint Annaland M et al. On the drag force of bubbles in bubble swarms at intermediate and high Reynolds numbers. Chem Eng Sci 2011;66:32043211.

[16] Richardson J, Zaki W. Sedimentation and fluidisation: Part I. Transactions of the Institution of chemical engineers 1954;32:35-53.

[17] Bochev PB, Gunzburger MD. Least-Squares finite element method. Springer Science+Business Media 2009.

[18] Yao BP, Zheng C. Bubble behavior and flow structure of bubble columns. Chem Eng Process 1991;29:65-75.

[19] Jakobsen HA, Hjarbo KW, Svendsen HF. Bubble size distributions, bubble velocities and local gas fraction from conductivity measurements. Tech. report, Sintef 1994.

[20] Wan Z, Chen Y, Huang Y. Legendre spectral galerkin method for second-kind volterra integral equations. Front Math China 2009;4:181-193.

[21] Coulaloglou CA, Tavlarides LL. Description of interaction processes in agitated liquid-liquid dispersions. Chem Eng Sci 1977;32:1289-1297.

[22] Prince MJ, Blanch HW. Bubble coalescence and break-up in air-sparged bubble columns. AIChE Journal 1990;36:1485-1499.

[23] Kim JW, Chang JH, Lee WK. Inhibition of bubble coalescence by the electrolytes. Korean J. of Chem. Eng. 1990;7:100-108.