Application of N-phase algebraic slip model and direct quadrature method of moments to the simulation of air-water flow in vertical risers and bubble column reactor

Computers and Chemical Engineering 90 (2016) 151–160

Contents lists available at ScienceDirect

Computers and Chemical Engineering journal homepage: www.elsevier.com/locate/compchemeng

Application of N-phase algebraic slip model and direct quadrature method of moments to the simulation of air-water flow in vertical risers and bubble column reactor K. Swiderski a,b,∗ , C. Narayanan a , D. Lakehal a a b

ASCOMP AG, Technoparkstrasse 1, 8005 Zurich, Switzerland Institute of Fluid Dynamics, D-MAVT, ETH Zurich, Switzerland

a r t i c l e

Article history: Received 25 September 2015 Received in revised form 14 April 2016 Accepted 16 April 2016 Available online 25 April 2016 Keywords: DQMOM Multiphase flow Algebraic slip model Population balance Bubble column Pipe flow

a b s t r a c t

i n f o

In the present work the Direct Quadrature Method of Moments (DQMOM) has been implemented into the commercial CFD code TransAT. DQMOM has recently become a very attractive approach for solving population balance equation (PBE) due to its capability of representing the most interesting properties of the population, eg. Sauter mean diameter, void fraction, number of particles. The DQMOM technique was coupled with the turbulent N-phase Algebraic Slip Model (ASM) model in order to extend the model to handle dispersed phase populations such that each class has its own velocity field. The results compared to experimental data show that the developed numerical model accurately predicts void fraction profile in a long riser within a bubbly flow regime. Moreover the model is used for the simulation of bubble column, proving that it accurately predicts the gas hold up and the Sauter mean diameter. © 2016 Elsevier Ltd. All rights reserved.

1. Introduction In many industrial applications and in nature as well, bubbles are dispersed in a liquid continuum. The understanding and the modeling of dispersed flows physics is therefore of high importance for many applications including chemical engineering (bubble columns), water treatment (ozonation and purification), nuclear industry (steam generators, accidental depressurization), naval transport (skin drag reduction) and medicine or biotechnology (contrast agent, micro bubbles bursting, virus replication in cell structures), see e.g. (Latorre, 1997; Yadigaroglu et al., 2008; Hsiao et al., 2010; Ishikawa et al., 2013; Jakobsen, 2014; Dürr et al., 2015). The dynamics of bubbly flows is very complex in nature. Typically various complex phenomena are present simultaneously, such as bubble collisions and coalescence or breakage. Almost all interactions between bubbles occur in the presence of turbulence or hydrodynamic mixing, which usually dictate their rate. Under such conditions, the assumption of constant bubble size may lead to incorrect predictions of gas-liquid multiphase flow behavior. The evolution and creation of bubbles of different sizes require the

∗ Corresponding author at: ASCOMP AG, Technoparkstrasse 1, 8005 Zurich, Switzerland. E-mail address: [email protected] (K. Swiderski). http://dx.doi.org/10.1016/j.compchemeng.2016.04.023 0098-1354/© 2016 Elsevier Ltd. All rights reserved.

representation of the particle size distribution. This is typically achieved by the Population Balance Models (PBM) (Ramkrishna and Mahoney, 2002), where the statistical distribution of the dispersed phase can be tracked. The N-phase model with algebraic slip between phases (Manninen et al., 1996) as implemented in TransAT (Lakehal et al., 2013) needs to be extended in this regard. Since the N-phase model solves for the volume fraction of each phase, a natural way to extend it would be direct discretization of the distribution into a number of size bins with birth and death terms accounting for processes such as breakup, coalescence, growth, etc. (Frank et al., 2005; Yeoh et al., 2012; Krepper et al., 2008). Among the other available PBM methods, such as the Method of Moments (MOM), Quadrature Method of Moments (QMOM) and their different variants (Marchisio and Fox, 2013), DQMOM (Marchisio and Fox, 2005) is chosen in this study. The main reason to not use the direct discretization method is the fact that momentbased methods are more efficient in capturing the distribution function. It has been shown that whereas 12–18 bins were required for the direct discretization method (Sanyal et al., 2005; Selma et al., 2010), equal or more accurate results were obtained with only 6 moments using the QMOM method. However, the QMOM method remains difficult to implement because the moment equations are directly solved, whose closure requires the recreation of the quadrature nodes and weights from the moments in order to

152

K. Swiderski et al. / Computers and Chemical Engineering 90 (2016) 151–160

Nomenclature CD CL CWL D Eo N Nc Pr Re R Sc S Y a˜ b˜ d d32 g j k kV  n p t u x yw

Drag force coefficient Lift force coefficient Wall force coefficient Diffusivity [m2 s−1 ] Eötvös number Total number of phases Total number of classes Prandtl number Reynolds number Radius [m] Schmidt number Shear rate tensor [s−1 ] Mass fraction Breakage frequency [m−3 s−1 ] Daughter size distribution function Diameter [m] Sauter mean diameter [m] Gravity acceleration [ms−2 ] Superficial velocity [ms−1 ] Turbulence kinetic energy [m2 s−2 ] Bubble shape factor Normal vector Pressure [Pa] Time [s] Velocity coordinate [m/s] Spatial coordinate [m] Wall distance [m]

Subscripts: b Bubble Class index c d Dispersed f Fluid i,j,k,l Vector index l Liquid Mixture m t Turbulent Superscripts: d Drift m Mixture Slip s td Turbulent dispersion Greek letters: ␣ Volume fraction ˜ ˇ Coalescence kernel [m3 s−1 ] ␰ ␻  ␮ ␳ ␴ ␧

Abscissa, node [m] Weight [m−3 ] Kinematic viscosity [m2 s−1 ] Dynamic viscosity [Pa s] Density [kg m−3 ] Surface tension [N m−1 ] Turbulence dissipation rate [m2 s−3 ]

close the terms with higher order moments. Also, the velocity field associated with a given moment is difficult to interpret if it depends strongly on the internal coordinate of the population balance equation; as it is typically in case of bubbles of widely different sizes. DQMOM, on the other hand, directly solves for the nodes and weights of the quadrature approximation instead of the chosen moments, and is therefore easier to implement, still offering good

accuracy at acceptable computational effort (Cheung et al., 2013). This is important in the N-phase context wherein the population balance might need to be solved for multiple dispersed phases. It is worth mentioning that DQMOM was originally developed for homogeneous systems (no spatial gradients of particle size of each class) and its application to inhomogeneous systems poses several challenges. One of the critical issues of DQMOM is the problem of zero weights, which makes it not possible to calculate the dispersed phase size when solving for the weighted abscissa. For practical problems, zero weights (or number density) can occur naturally due to segregation of particles due to the flow or due to condensation of bubbles. Another issue is the condition where two particle classes converge to the same size making one of them redundant. This leads to ill-conditioning of the system. This paper presents the implementation of monovariate DQMOM into TransAT with all the necessary improvements for handling inhomogeneous flows, such as dealing with the ill-conditioning of the system due to the above mentioned phenomena. The method is then successfully applied to the simulation of complex flows such as vertical risers and bubbles columns which are of direct industrial relevance. 2. Mathematical model In order to describe the dynamics of N-phase flows an Algebraic Slip Model (ASM) extended to N-phases is used in the present work. Within this framework the DQMOM method is used to capture the most important properties of the dispersed phase like bubble size, void fraction or interfacial area. In this study, only the population of one of the phases is represented using DQMOM. 2.1. N-phase algebraic slip model Multiphase gas-liquid flows can be tackled using either interface tracking methods (ITM) or phase-average models, for both laminar and turbulent flows. Specifically, the Level-Set approach, the phasefield variant and the Volume-of-Fluid methods can be employed as ITM’s. In the context of phase-average models, eg. the mixture (homogeneous) approach applied to gas-liquid systems, transport equations are solved for the mixture quantities rather than for the phase-specific quantities, unlike in the two-fluid model. This implies that one mixture momentum equation is solved for the entire flow system, reducing the number of equations to be solved in comparison to the two-fluid model. In many situations however, the model must be employed with a prescribed closure law for the interphase slip velocity and associated stresses. In this case, the model is known as homogeneous ASM. The N-phase approach is an extension of the homogeneous ASM introduced above, and is invoked in situations involving more than two fluid phases. The N-phase approach could as well be used in the two-fluid flow context. In the homogeneous ASM framework, the N-phase features a modified scheme where mass conservation and energy equations are solved for each phase to better cope with interphase mass transfer, whereas the momentum is solved for the mixture. Further, the model can be used either under its homogeneous form or by adding an algebraic slip velocity to separate the phases. The N-phase ASM represents multiphase flow in an ensemble averaged sense, where the involved phases move at different or equal velocities under the assumption of equilibrium within short spatial length scales (Manninen et al., 1996). The model also requires that the relaxation time for dispersed phase is short in comparison with changes in the flow, therefore Stokes number «1. This condition is generally met when dispersed phase particles are small or carrier fluid viscosity is high. Thus the model is suitable for modeling rising bubbles in a liquid. The model is given in the

K. Swiderski et al. / Computers and Chemical Engineering 90 (2016) 151–160

form of one mass conservation equation for each of the N phases and one momentum conservation equation for the mixture, see Eqs. (1) and (2) below. Momentum conservation equation may contain additional terms to account for velocity differences between the phases. The relative velocities must be computed from the force balance for the dispersed phases. The N-phase ASM for incompressible isothermal flow reads as follows:

∂k ˛k ∂ ∂ (m) (d) (td) + [k ˛k (uj + ujk + u˜ j )] = ∂t ∂xj ∂xj

 

t ∂k ˛k Sct ∂xj



,

∂ ∂ ∂p ∂ (m) (m) (m) (m) (m ui uj ) = − + 2 (m + t ) Sij (m ui ) + ∂t ∂xj ∂xi ∂xj ∂  (d) (d) − ( ˛k k uik ujk ) + m gi + M  . ∂xj N

(1)

 (2)

k=1

density, p =

N

˛ p k=1 k k

is the mixture

N

 ˛ is k=1 k k (m) pressure, Sij

2.2. Closures for drag, lift and wall lubrication forces The interphase forces acting between bubbles and liquid phase are complex due to the fact that bubbles can deform in three dimensions and therefore resulting shape of the interface is not necessarily spherical. This makes difficulties in terms of modeling the interphase momentum exchange. Tomiyama et al. (1998) investigated drag coefficient for a single bubble moving in a stagnant liquid. The result of their work is the expression for drag coefficient for pure, slightly contaminated and fully contaminated systems. Here, we use the latter one which yields: CD = max

 16  Reb

8 Eo 

1 + 0.15Re0.687 , b

3 Eo + 4

(5)

,

where Reb = (l db |ujkl (s) |)/l is the bubble Reynolds number, and

The velocity associated with the phase k is a sum of mixture, drift and turbulent drift contributions coming from components ˜ j (td) respectively. In the above equations, ˛k ujk = uj (m) + ujk (d) + U is the volume fraction of phase k, m =

153

the mixture

2

Eo = [g(l − d )db ]/ is the Eötvös number. The lift force, which is primarily responsible of transverse bubble migration in a pipe flow, is modeled by employing Tomiyama et al. (2002) correlation for the lift force coefficient as:

is the shear

rate tensor for the mixture, m is the mixture viscosity given by a rheological law, e.g. (Ishii and Zuber, 1979), t is the turbulent eddy viscosity, gi is gravitational acceleration, and M  represents the effect of surface tension (ignored in our consideration). If k is the dispersed phase and l is the continuous phase, then the drift velocity N of phase l is given as ujl (d) = Y u (s) , where Yk denotes mass k=1 k jkl fraction of phase k, and ujkl (s) is the slip between the phases l and k given by:

CL =

⎧  ⎨ min [0.288 tanh (0.121Reb ) , f (Eo )] ⎩

forEo < 4

f (Eo )

for 4 ≤ Eo ≤ 10 ,

−0.27

forEo > 10

(6)

where f (Eo ) = 0.00105Eo 3 − 0.0159Eo 2 − 0.0204Eo + 0.474. The sign of CL changes at db ≈ 2.8mm from positive to negative, which favors smaller bubbles migration to the wall and larger to

(3) which represents the balance between the drag, lift, wall lubrication and buoyancy forces. Note that for the basic N-phase ASM model, Rk is an input for each phase and remains constant. The normal vector nW points away from the wall and (ujl − ujk )|| denotes phase relative velocity component tangential to the wall surface. For incompressible flows, the mixture continuity equation is given as,

∂u(m) j ∂xj

+

∂u˜ (td) j ∂xj (td)

∂  (d) ( ˛k ujk ) = 0 ∂xj

+



(4)

k=1

where U˜j = t ∇ xm /(Sct m ) is the turbulent drift correction proposed by Ishii (1975), and t /Sct denotes the mixture turbulent diffusivity, Sct is the turbulent Schmidt number, and t = t m is the mixture eddy viscosity. Eq. (4) represents the volume balance and is used to derive the pressure correction equation. Turbulence is modeled with the standard k − ε model for the mixture (Launder and Spalding, 1974). Bubble induced turbulence was not considered in this work. This effect might be important however is still unresolved issue. Therefore we do not consider it in our model. Bubbly flows are turbulent from their nature. From practical point of view two methods of modeling turbulent flows are established well to date, Large Eddy Simulation (LES) and Reynolds Averaged Navier Stokes (RANS) equations based methodology. Both can be applied for bubbly flow simulation in a pipe and are capable of predicting fundamental turbulent quantities such as dissipation rate and turbulent kinetic energy (Vijiapurapu and Cui, 2010). However RANS models can suffer in predicting flow separations because they cannot resolve complex vertical motions. Therefore if the flow is strongly instationary LES methodology is generally better choice.



0.757 1/3

Eötvös number and dH = db 1 + 0.163Eo is the maximum horizontal bubble dimension (Wellek et al., 1966). Wall lubrication force tends to push bubbles away from the wall. The model developed by Antal et al. (1991) is used here. The wall force coefficient is given as follows:



CWL = max

N

2

the pipe center. In Eq. (6) Eo = [g(l − d )dH ]/l is the modified

CW 1 CW 2 0, + yw db



,

(7)

where the non-dimensional coefficients CW 1 = −0.01 and CW 2 = 0.05 and yw is the distance to the nearest wall. 2.3. Direct quadrature method of moments The DQMOM method is based on direct tracking of the weights ωc and weighted abscissas ωc c appearing in the quadrature approximation for the first 2Nc moments of the discrete phase distribution. In the univariate version, DQMOM is developed for Nc classes of the internal variable (node, abscissa) c , which in our case represents the dispersed phase diameter. The volume fraction of the c-th class is given as ˛c = kV ωc c3 where the shape factor kV = /6 for spherical particles. 2.3.1. Transport equations and moment conservation constraint equation In order to extend the N-phase method to include a population balance of the dispersed phase, the mass balance for the dispersed phase is replaced by the solution of the DQMOM size classes represented by the weights and weighted nodes of the class. In this study,

154

K. Swiderski et al. / Computers and Chemical Engineering 90 (2016) 151–160

only two-phase flows are considered where one of the phases is chosen as the dispersed phase. The following system:

∂ωc ∂ ∂ ∂ωc ujc ωc = (Dx ) + ac , + ∂t ∂xj ∂xj ∂xj (8) ∂ ∂ ∂ ∂ ujc ωc c = (Dx ωc c ) + bc , c = 1, ..., Nc , ωc c + ∂t ∂xj ∂xj ∂xj is solved for the dispersed phase with the use of the implicit solver in TransAT. We use first-order scheme for DQMOM equations in the implicit solver to avoid problems with non-conservation as reported in (Mazzei, 2013). Let us recall that the velocity associated with c-th class is given in the same way as in Eq. (2) for k-th phase ˜ j (td) . as, ujc = uj (m) + ujc (d) + U The algebraic linear system, which can be referred as moment conservation constraint, for the RHS terms of Eq (8), namely ac andbc , emerging from DQMOM framework reads as follows: Nc 

(1 − k)

ck ac

+k

Nc 

c=1

ck−1 bc

c=1

= Sk + Ck ,

k = 0, ..., 2Nc − 1.

(9)

N

c Here CK = k(k − 1)Dx c=1 ( c k−1 ωc ∇ x c · ∇ x c) is the diffusion correction source term for k-th moment of the distribution. Dx is the diffusivity in physical space and represents any diffusion source that may affect the simulation, such as turbulent diffusion, molecular diffusion or numerical diffusion. Sk stands for physical interactions source term, such as coalescence and breakage between particles.

2.3.2. Closures for bubble coalescence and breakage The general form for bubble coalescence and breakage kernels is given as (Marchisio and Fox, 2005): Sk =

Nc Nc    1

2

3 3 c1 + c2



k/3

k ˜ c1,c2 ωc1 ωc2 ˇ − c1

c1=1c2=1

 Nc

+

(k) b¯ c1 a˜ c1 ωc1

−

c1=1

Nc 

(10)

where (k) b¯ c1 =

∞

˜ k b( | c1 )d

(11)

0

˜ is fragmentation kernel and b( | c1 ) is daughter size distribution ˜ function, ˇc1,c2 is bubble coalescence kernel between the bubble of class c1 and c2, a˜ c1 denotes breakage frequency. It is very important for the mass conservation to normalize the integral value of Eq. (11) after solving it numerically. In the present work we employ model used by Laakkonen et al. (2007) for breakage and Coulaloglou and Tavlarides (1977) expression for coalescence frequency and efficiency:



˜ c1,c2 = Cagg ε1/3 c1 + c2 ˇ

 exp −Ccoll l l

ε 2





  Nf =

2 3 c1

Cfrag 4 + 3 3

2/3

c1 c2 c1 + c2

5/3





3 c1

2 

1/2

4  ,

1/2 

l c1 ε2/3 3

2/3

c1 + c2

C1 

a˜ c1 = Cbreak ε1/3 erfc

˜ b( | c1 ) = Nf

2 

C2 l + √ 4/3 c1 ε1/3 l d

1−

3 3 c1

2.3.3. Numerical method The time scales of the dispersed phase interaction process are typically much smaller than the convection and diffusion time scale. This makes the combined DQMOM system very stiff leading to poor convergence of the equations. This problem is overcome by using the fractional step approach. The convection-diffusion equations for the weighted nodes and weights are first solved without considering the effect of the source terms (Sk = 0) to obtain inter∗ mediate values (ωc ∗ and c ) as follows:

∂ ∂ ∂ωc ∂ωc ujc ωc = (Dx ) + a∗ c , + ∂t ∂xj ∂xj ∂xj (13) ∂ ∂ ∂ ∂ ujc ωc c = (Dx ωc c ) + b∗ c , c = 1, ..., Nc . ωc c + ∂t ∂xj ∂xj ∂xj The set of Eq. (9) must be solved in each computational cell. In order to solve Eq. (9) we use DGETRS routine from LAPACK library (Anderson et al., 1999). When there is no inhomogeneity of the spatial distribution of abscissas (particle size c for c-th class is not changing in spacex), the right hand side of Eqs (8) and (13) can be simply set to zero. Otherwise the linear system (9) is solved with respect to diffusion correction term only Ck . In general, it can be shown that in the absence of dispersed phase interactions (DQMOM nodes remain spatially constant), the volume fraction advected via the ASM model and the volume fraction calculated for a given DQMOM class should match exactly (Swiderski et al., 2016). After the converged solution for the convection-diffusion step is obtained, the effect of the source terms is added in a second step.

∂ω∗ c = ac , ∂t

k a ˜ c1 ωc1 , c1

c1=1

The constants of the model that we use are as follows: Cagg = 0.88, Ccoll = 6e9, Cbreak = 6, C1 = 0.04, C2 = 0.01. For binary breakage Cfrag = 2. The employed bubble coalescence model is one of the oldest known, however still offers good accuracy when compared to more recent models, especially for the void fraction profile and bubble size distribution (Deju et al., 2013; Deju et al., 2015).

,

Cfrag

,

3 2 9 + 16.5Cfrag + 9Cfrag + 1.5Cfrag



.

(12)

∂(ω∗ c ∗ c ) = b c , c = 1, ..., Nc , ∂t

(14)

where a c and b c are evaluated from Eq. (9) by considering physical processes present in Sk term. The (*) denotes intermediate state of the solution obtained after the convection-diffusion step. Here the Ck term is not considered since it is already applied in convection-diffusion update. The Double precision Variable-coefficient Ordinary Differential Equation solver, with fixed-leading-coefficient implementation (DVODE) (Brown et al., 1989) is used to solve the set of ordinary differential equations. For the details how to compute kernel term Sk the reader can refer to (Marchisio and Fox, 2005; Marchisio et al., 2003). 2.3.4. Ill-conditioned matrix treatment During the implementation of DQMOM into any CFD code, one can meet issues regarding the aforementioned algebraic linear system ill-conditioning. This especially can happen when two or more abscissas are not sufficiently distinct from each other within the same computational cell (singularity problem), and/or the weights or abscissas values tend to zero values. To overcome this problem we use a condition number based algorithm to detect singularities where we consider the system as poorly conditioned if the relative difference between two nodes is less than 1–3% of their arithmetic mean value (non-distinct nodes), and if they appear, we use the adaptive perturbation method. It means that if the perturbation value is not enough for condition number reduction we increase its value until the system can be considered as well conditioned. Typically, if the condition number is higher than 106 , we consider the matrix as ill-conditioned. We also incorporated a heuristic rule of

K. Swiderski et al. / Computers and Chemical Engineering 90 (2016) 151–160

155

Fig. 1. Flow map regime. Picture adapted from Ref. (Yeoh et al., 2013).

Fox (2006) that is essentially rescaling the nodes in Eq. (9) by local maximal value of abscissa which reduces the condition number. Moreover we use iterative algorithm described in (Press et al., 1992) to minimize round-off error when computing Eq. (9). Three iterations are usually enough in this regard (Fox, 2006). We also found that the methods reported by other researchers (Fan et al., 2004), like local source term averaging or setting source term to zero for problematic cells, lead to non-conservation of the moments and solver failure. The method to deal with very low values of the weights can be enforcing a certain fractional number of particles per unit volume, e.g. 1–50 particles so that their impact on the mixture composition is negligible for practical applications. Because any intervention made to the nodes distribution should be mass conservative, if any fraction is added to i-th class then the same fraction is subtracted from j-th class when possible, or when impossible, then the artificial fraction is subtracted from continuous phase and used as corrective one.

3. Model validation Two problems are considered in order to perform the validation of computational simulation model, namely: vertical riser and bubble column reactor. The vertical riser has been experimentally investigated by Szalinski et al. (2010) and also by Hibiki et al. (2001). The available data are used here to confront simulation results. Five cases are modeled here according to the following flow map regime shown in Fig. 1. The solid line represents flow regime transition boundary according to the model of Taitel et al. (1980). From the problem setup point of view all the considered cases have one inflow and outflow boundary condition, and wall aside. At bottom of each one air is injected (as a pure gas or in a liquid jet). The problem sketch is presented in Fig. 2. For the bubble column reactor work of Cachaza Gianzo (2011) is used as reference data. Both the considered problems are within the bubbly flow regime, according to Yeoh et al. (2013), therefore the same models and their constants given in Eq (12) were used for bubble coalescence and breakage. The boundary conditions are presented in Table 1. All the inflow boundaries assume lognormal distribution with standard deviation of 15%.

Fig. 2. Sketch of simulated cases.

3.1. Pipe flows 3.1.1. Vertical riser of Szalinski The experimental facility consists of 6 m long and 67 mm diameter pipe and the mixing device located at bottom of the pipe. Water and air come into the mixer and then are distributed by the sparger being a perforated plate with a number of 3 mm holes. According to Szalinski et al. (2010) the flow is well developed for x/D = 40 (where z is vertical coordinate), which corresponds to x = 2.7 m. Therefore in the computational model the pipe is truncated at the end so that its total length is 3.1m. The experimental measurements have been made at x/D = 74.6 and are used here to compare simulation results probed at the end of the modeled pipe. The grid used for the validation consisted of 750 cells in vertical and 34 cells in radial direction, corresponding to the minimal cell size of 2.7 and 0.35 mm, respectively. The grid was refined in vertical direction for the distance of 1 m as well as near to the wall through the whole pipe. Turbulence intensity at the inflow was set to 10% and eddy viscosity ratio as 120. Bubble population was represented by two classes. The simulation was unsteady and the final solution was captured when a quasi-steady state was obtained. The instantaneous field of water volume fraction, Sauter mean diameter and DQMOM weights are shown in Fig. 3. As expected the flow pattern is not changing from x = 2.7 m, which is in perfect agreement with experiment. Most of breakage process takes place nearby the wall where turbulent dissipation rate is significant in comparison to the rest of the pipe. On the other hand coalescence dominates in the core of the riser. Small bubbles are pushed towards the wall by the lift force, whereas the larger ones are in the middle of the pipe. The computed and measured air void fraction are presented in Fig. 4. The overall match is very good in the middle part of the riser as well as in the vicinity of the wall.

156

K. Swiderski et al. / Computers and Chemical Engineering 90 (2016) 151–160

Fig. 3. Water volume fraction, Sauter mean diameter, number of bubbles (weights) for two DQMOM classes at the final state of simulation.

Fig. 4. Air void fraction profile for fully developed flow. Table 1 Boundary conditions of considered test cases. SI units, except d32 [mm]. Case Szalinski Hibiki

Cachaza Gianzo

Inflow BC

Outflow BC

Wall BC

1 2 3 4

d32 = 7.5, jd = 0.05, jf = 0.25, ␣d = 0.063 d32 = 2.5, jd = 0.0275, jf = 0.491, ␣d = 0.05 d32 = 2.5, jd = 0.0556, jf = 0.491, ␣d = 0.10 d32 = 2.5, jd = 0.0473, jf = 0.986, ␣d = 0.05 d32 = 2.5, jd = 0.1130, jf = 0.986, ␣d = 0.10

p=0 p=0

no-slip no-slip

1 2 3 4 5

d32 = 6.8, jd = 2.3e-3, jf = 0, ␣d = 0.50 d32 = 6.8, jd = 7.1e-3, jf = 0, ␣d = 0.50 d32 = 6.8, jd = 11.9e-3, jf = 0, ␣d = 0.50 d32 = 6.8, jd = 16.6e-3, jf = 0, ␣d = 0.50 d32 = 6.8, jd = 21.3e-3, jf = 0, ␣d = 0.50

p=0

slip

K. Swiderski et al. / Computers and Chemical Engineering 90 (2016) 151–160

157

Fig. 5. Time averaged gas void fraction at z/D = 6.

Fig. 6. Time averaged gas void fraction at z/D = 53.5.

3.1.2. Vertical riser of Hibiki Similarly to the previous experiment, two-phase bubbly flow was investigated in a vertically positioned pipe by Hibiki et al. (2001). Here the pipe length is 3.06 m and its diameter 50.8 mm. Experimental data are available at z/D = 6 and 53.5, and contain gas volume fraction and gas/liquid phase velocity distribution along radial coordinate. Moreover we confront here also MUSIG model results presented in (Yeoh et al., 2013). The MUSIG model employs coalescence and breakage mechanisms of Prince and Blanch (1990) and Luo and Svendsen (1996). Furthermore models of Wu et al. (1998), Hibiki and Ishii (2002) and Yao and Morel (2004) of the single transport equation of the average number density for the investigated problem are also compared here. All the reference numerical results used here have been achieved with the use of ANSYS-CFX by Yeoh et al. (2013).

The following values for the turbulence intensity and eddy viscosity ratio were used at the inflow: 10% and 120. Three classes of DQMOM variables were used in the simulation. From time averaged void fraction profile shown in Fig. 5 for z/D = 6, it can be concluded that DQMOM performs well and qualitatively better than MUSIG for the four cases investigated here. The highest discrepancies from experimental data are observed for jf = 0.491m/s especially in the vicinity of the wall. For jf = 0.986m/s the numerical solution is almost in line with experiment and a low deficit of the void fraction is close to the pipe center. At the second location z/D = 53.5 (Fig. 6) void fraction profile still matches well. It is a little underpredicted for jg = 0.0556m/s however there is strong variation in experimental data close to the wall, which can be dictated by a strong gradient in transverse velocity caused by the lift force. Nevertheless DQMOM produces acceptable results which

158

K. Swiderski et al. / Computers and Chemical Engineering 90 (2016) 151–160

Fig. 7. Time averaged gas velocity at z/D = 53.5.

Fig. 8. Time averaged liquid velocity at z/D = 53.5.

look better than the MUSIG reference data where the void fraction is overestimated in the middle part of the pipe and underestimated nearby the wall. For the gas and liquid velocities a good match to experimental data is obtained, as shown in Fig. 7 and Fig. 8. The gas and liquid velocities for jf = 0.491m/s are very close to measurement data at z/D = 53.5. A noticeable difference appears for jf = 0.986m/s, but still DQMOM output looks qualitatively fine and has similar trend as experimental data. Overall DQMOM gives better results in comparison to the other methods presented. More investigation is needed for the better resolution of the void fraction profile nearby the wall, especially for lower superficial velocities (case 1 and 2). The suggestion might be to employ a correction for turbulent viscosity due to presence of

the bubbles. However in this paper we use standard expression for turbulent quantities as dictated by the standard k − ε model.

3.2. Bubble column reactor of Cachaza Gianzo Here a cuboidal polymethyl methacrylate bubble column 0.2 m wide, 1.8 m high, 0.04 m deep was used in experimental setup (Cachaza Gianzo, 2011; Díaz et al., 2008). The experimental inlet pattern was comprised of 8 holes with 1 mm diameter. The computational model is two-dimensional with the inlet section of 2 cm located centrally at bottom of the tank (see Fig. 9). The initial water head is 0.45m. Full slip boundary conditions are used at the walls, and pressure outlet boundary condition with normal conditions is set on top of the computational domain. Three classes of bub-

K. Swiderski et al. / Computers and Chemical Engineering 90 (2016) 151–160

159

diffusion correction, as originally formulated in (Marchisio and Fox, 2005), to suppress non-conservation of moments higher than the first two. We have demonstrated that CFD-PBE based method is a powerful tool for solving complex problems like bubbly flow in pipes and bubble columns. We obtained a good agreement between the numerical results and experimental data for all the investigated pipe flow cases. For bubble columns the model has been verified and gas holdup and Sauter mean diameter results were validated. It indicates that the used modeling strategy with chosen breakage and coalescence mechanisms might be a good compromise among other available models and methods for two-phase bubbly flows. The next step will be to extend the model to a bubble column with mass transfer and validation by experimental data. From an application point of view, extension of the population balance model to handle multiple dispersed phases would be very useful. Acknowledgements

Fig. 9. Instantaneous contours of air void fraction (logscale) and velocity vectors [m/s].

References

Table 2 Gas hold-up and Sauter diameter comparison to experimental data. Case

1: jg = 2.3e-3 2: jg = 7.1e-3 3: jg = 11.9e-3 4: jg = 16.6e-3 5: jg = 21.3e-3

Gas hold-up

This work is partially supported from project TRIUMPH (Treating Urban Micropollutants and Pharmaceuticals in Wastewaters) being a part of ACQUEAU − the Eureka Cluster for Water. Financial support of CTI Switzerland under contract 15301.3 PFNM-NM and contribution of Prof. Patrick Jenny from Institute of Fluid Dynamics, ETH Zurich, are gratefully acknowledged.

d32 at the center of the column

TransAT (DQMOM)

Exp.

Transat (DQMOM)

Exp.

0.66% 1.75% 2.75% 3.31% 4.22%

0.69% 1.81% 2.63% 3.36% 4.10%

7.20 7.45 7.50 7.65 7.70

6.83 7.05 6.50 6.40 7.73

bles represent their population here in context of DQMOM method. The Sauter mean diameter is set as 6.8 mm at the inlet and the air volume fraction is equal 0.5 at inlet as well. The inlet turbulence eddy viscosity ratio is 100 and turbulence intensity 5%. The mesh is uniform and fine with the cell size of 1 mm. Five cases are investigated as shown in Table 1 and Table 2. Fairly good agreement is obtained with experimental data for both gas holdup and Sauter mean diameter. The highest discrepancy is observed for Case 3 and 4 for d32 measured at the center of the column. This is however difficult to explain since the experimental trend is not monotonic function whereas in simulation d32 increases with increasing jg , and for the latter considered Case 5 almost perfect agreement is obtained for d32 . This can suggest that the reference data should be considered with some dose of uncertainty, at least for Case 3 and 4. In spite of aforementioned observation, the gas hold-up value always follows the experimental trend, and the overall match is actually good. 4. Conclusions The paper describes fundamental details of N-phase ASM model and DQMOM. The first strategy is widely used for solving practical multiphase problems while the second one allows incorporating population balance in a feasible way and at a reasonable computational cost. Both methods are applicable for homogeneous and inhomogeneous problems. However special attention must be paid to proper implementation of DQMOM into a CFD framework. Especially for inhomogeneous systems, DQMOM formulation demands

Anderson, E., Bai, Z., Bischof, C., Blackford, L.S., Demmel, J., Dongarra, J.J., et al., 1999. LAPACK Users’ Guide, 3rd ed.). Society for Industrial and Applied Mathematics, Philadelphia, PA, USA. Antal, S.P., Lahey Jr., R.T., Flaherty, J.E., 1991. Analysis of phase distribution in fully developed laminar bubbly two-phase flow. Int. J. Multiph. Flow. 17, 635–652, http://dx.doi.org/10.1016/0301-9322(91)90029-3. Brown, P., Byrne, G., Hindmarsh, A., 1989. VODE: a variable-Coefficient ODE solver. SIAM J. Sci. Stat. Comput. 10, 1038–1051, http://dx.doi.org/10.1137/0910062. Cachaza Gianzo, E.M., 2011. Hydrodynamics and Mass Transfer Effects in Bubble Columns. Universidad de Salamanca. Cheung, S.C.P., Deju, L., Yeoh, G.H., Tu, J.Y., 2013. Modeling of bubble size distribution in isothermal gas–liquid flows: numerical assessment of population balance approaches. Nucl. Eng. Des. 265, 120–136, http://dx.doi. org/10.1016/j.nucengdes.2013.08.049. Coulaloglou, C.A., Tavlarides, L.L., 1977. Description of interaction processes in agitated liquid–liquid dispersions. Chem. Eng. Sci. 32, 1289–1297, http://dx. doi.org/10.1016/0009-2509(77)85023-9. Díaz, M.E., Iranzo, A., Cuadra, D., Barbero, R., Montes, F.J., Galán, M.A., 2008. Numerical simulation of the gas–liquid flow in a laboratory scale bubble column: influence of bubble size distribution and non-drag forces. Chem. Eng. J. 139, 363–379, http://dx.doi.org/10.1016/j.cej.2007.08.015. Dürr, R., Müller, T., Kienle, A., 2015. Efficient DQMOM for multivariate population balance equations and application to virus replication in cell cultures. IFAC-Pap. 48, 29–34, http://dx.doi.org/10.1016/j.ifacol.2015.05.045. Deju, L., Cheung, S.C.P., Yeoh, G.H., Tu, J.Y., 2013. Capturing coalescence and break-up processes in vertical gas–liquid flows: assessment of population balance methods. Appl. Math. Model. 37, 8557–8577, http://dx.doi.org/10. 1016/j.apm.2013.03.063. Deju, L., Cheung, S.C.P., Yeoh, G.H., Qi, F., Tu, J., 2015. Comparative analysis of coalescence and Breakage kernels in vertical gas-liquid flow. Can. J. Chem. Eng. 93, 1295–1310, http://dx.doi.org/10.1002/cjce.22196. Fan, R., Marchisio, D.L., Fox, R.O., 2004. Application of the direct quadrature method of moments to polydisperse gas–solid fluidized beds. Powder Technol. 139, 7–20, http://dx.doi.org/10.1016/j.powtec.2003.10.005. Fox, R.O., 2006. Bivariate direct quadrature method of moments for coagulation and sintering of particle populations. J. Aerosol. Sci. 37, 1562–1580, http://dx. doi.org/10.1016/j.jaerosci.2006.03.005. Frank T., Zwart, P.J., Shi, J.-M., Krepper, E., Lucas, D., Rohde, U., 2005 Inhomogeneous MUSIG model −a population balance approach for polydispersed bubbly flows, Bled, Slovenia. Hibiki, T., Ishii, M., 2002. Development of one-group interfacial area transport equation in bubbly flow systems. Int. J. Heat Mass Transf. 45, 2351–2372, http://dx.doi.org/10.1016/S0017-9310(01)00327-1. Hibiki, T., Ishii, M., Xiao, Z., 2001. Axial interfacial area transport of vertical bubbly flows. Int. J. Heat Mass Transf. 44, 1869–1888, http://dx.doi.org/10.1016/ S0017-9310(00)00232-5. Hsiao, C.-T., Lu, X., Chahine, G., 2010. Three-dimensional modeling of the dynamics of therapeutic ultrasound contrast agents. Ultrasound Med. Biol. 36, 2065–2079, http://dx.doi.org/10.1016/j.ultrasmedbio.2010.08.022.

160

K. Swiderski et al. / Computers and Chemical Engineering 90 (2016) 151–160

Ishii, M., Zuber, N., 1979. Drag coefficient and relative velocity in bubbly, droplet or particulate flows. AIChE J. 25, 843–855, http://dx.doi.org/10.1002/aic. 690250513. Ishii, M., 1975. Thermo-fluid Dynamic Theory of Two-phase Flow. Eyrolles, Paris. Ishikawa, A., Imai, R., Tanaka, T., 2013. Experimental Study on Two-phase Flow in Horizontal Tube Bundle Using SF6-water (accessed June 20, 2015) http://inis. iaea.org/Search/search.aspx?orig q=RN:45069334. Jakobsen, H.A., 2014. Chemical Reactor Modeling. Springer International Publishing Cham (accessed June 20, 2015) http://link.springer.com/10.1007/ 978-3-319-05092-8. Krepper, E., Lucas, D., Frank, T., Prasser, H.-M., Zwart, P.J., 2008. The inhomogeneous MUSIG model for the simulation of polydispersed flows. Nucl. Eng. Des. 238, 1690–1702, http://dx.doi.org/10.1016/j.nucengdes.2008.01.004. Laakkonen, M., Moilanen, P., Alopaeus, V., Aittamaa, J., 2007. Modelling local bubble size distributions in agitated vessels. Chem. Eng. Sci. 62, 721–740, http://dx.doi.org/10.1016/j.ces.2006.10.006. Lakehal, D., Narayanan, C., Caviezel, D., von Rickenbach, J., Reboux, S., 2013. Progress in computational microfluidics using TransAT. Microfluid. Nanofluidics. 15, 415–429, http://dx.doi.org/10.1007/s10404-013-1163-3. Latorre, R., 1997. Ship hull drag reduction using bottom air injection. Ocean Eng. 24, 161–175, http://dx.doi.org/10.1016/0029-8018(96)00005-4. Launder, B.E., Spalding, D.B., 1974. The numerical computation of turbulent flows. Comput. Methods Appl. Mech. Eng. 3, 269–289, http://dx.doi.org/10.1016/ 0045-7825(74)90029-2. Luo, H., Svendsen, H.F., 1996. Theoretical model for drop and bubble breakup in turbulent dispersions. AIChE J. 42, 1225–1233, http://dx.doi.org/10.1002/aic. 690420505. Manninen, M., Taivassalo, V., Kallio, S., 1996. On the Mixture Model for Multiphase Flow, 288. VTT Publ. Marchisio, D.L., Fox, R.O., 2005. Solution of population balance equations using the direct quadrature method of moments. J. Aerosol Sci. 36, 43–73, http://dx.doi. org/10.1016/j.jaerosci.2004.07.009. Marchisio, D.L., Fox, R.O., 2013. Computational Models for Polydisperse Particulate and Multiphase Systems. University Press, Cambridge. Marchisio, D.L., Vigil, R.D., Fox, R.O., 2003. Quadrature method of moments for aggregation–breakage processes. J. Colloid Interface Sci. 258, 322–334, http:// dx.doi.org/10.1016/S0021-9797(02)00054-1. Mazzei, L., 2013. Segregation dynamics of dense polydisperse fluidized suspensions modeled using a novel formulation of the direct quadrature method of moments. Chem. Eng. Sci. 101, 565–576, http://dx.doi.org/10.1016/ j.ces.2013.07.006. Press, W.H., Flannery, B.P., Teukolsky, S.A., Vetterling, W.T., 1992. Numerical Recipes in FORTRAN, 77., 2nd edn. University Press, Cambridge. Prince, M.J., Blanch, H.W., 1990. Bubble coalescence and break-up in air-sparged bubble columns. AIChE J. 36, 1485–1499, http://dx.doi.org/10.1002/aic. 690361004.

Ramkrishna, D., Mahoney, A.W., 2002. Population balance modeling. Promise for the future. Chem. Eng. Sci. 57, 595–606, http://dx.doi.org/10.1016/S00092509(01)00386-4. Sanyal, J., Marchisio, D.L., Fox, R.O., Dhanasekharan, K., 2005. On the comparison between population balance models for CFD simulation of bubble columns. Ind. Eng. Chem. Res. 44, 5063–5072, http://dx.doi.org/10.1021/ie049555j. Selma, B., Bannari, R., Proulx, P., 2010. Simulation of bubbly flows: comparison between direct quadrature method of moments (DQMOM) and method of classes (CM). Chem. Eng. Sci. 65, 1925–1941, http://dx.doi.org/10.1016/j.ces. 2009.11.018. Swiderski, K., Caviezel, D., Labois, M., Lakehal, D., Narayanan, C., 2016. Computational modelling of gas–liquid multiphase flows with DQMOM and the N-phase Algebraic Slip Model. Comput. Methods Multiph. Flow VIII, 2015. WIT Press, Valencia, pp. 299–310, http://dx.doi.org/10.2495/MPF150261. Szalinski, L., Abdulkareem, L.A., Da Silva, M.J., Thiele, S., Beyer, M., Lucas, D., et al., 2010. Comparative study of gas–oil and gas–water two-phase flow in a vertical pipe. Chem. Eng. Sci. 65, 3836–3848, http://dx.doi.org/10.1016/j.ces.2010.03. 024. Taitel, Y., Bornea, D., Dukler, A.E., 1980. Modelling flow pattern transitions for steady upward gas-liquid flow in vertical tubes. AIChE J. 26, 345–354, http:// dx.doi.org/10.1002/aic.690260304. Tomiyama, A., Kataoka, I., Zun, I., Sakaguchi, T., 1998. Drag coefficients of single bubbles under normal and micro gravity conditions. JSME Int. J. Ser. B. 41, 472–479, http://dx.doi.org/10.1299/jsmeb.41.472. Tomiyama, A., Tamai, H., Zun, I., Hosokawa, S., 2002. Transverse migration of single bubbles in simple shear flows. Chem. Eng. Sci. 57, 1849–1858, http://dx.doi. org/10.1016/S0009-2509(02)00085-4. Vijiapurapu, S., Cui, J., 2010. Performance of turbulence models for flows through rough pipes. Appl. Math. Model. 34, 1458–1466, http://dx.doi.org/10.1016/j. apm.2009.08.029. Wellek, R.M., Agrawal, A.K., Skelland, A.H.P., 1966. Shape of liquid drops moving in liquid media. AIChE J. 12, 854–862, http://dx.doi.org/10.1002/aic.690120506. Wu, Q., Kim, S., Ishii, M., Beus, S.G., 1998. One-group interfacial area transport in vertical bubbly flow. Int. J. Heat Mass Transf. 41, 1103–1112, http://dx.doi.org/ 10.1016/S0017-9310(97)00167-1. Yadigaroglu, G., Simiano, M., Milenkovic, R., Kubasch, J., Milelli, M., Zboray, R., et al., 2008. CFD4NRS with a focus on experimental and CMFD investigations of bubbly flows. Nucl. Eng. Des. 238, 771–785, http://dx.doi.org/10.1016/j. nucengdes.2007.02.047. Yao, W., Morel, C., 2004. Volumetric interfacial area prediction in upward bubbly two-phase flow. Int. J. Heat Mass Transf. 47, 307–328, http://dx.doi.org/10. 1016/j.ijheatmasstransfer.2003.06.004. Yeoh, G.H., Cheung, S.C.P., Tu, J.Y., 2012. On the prediction of the phase distribution of bubbly flow in a horizontal pipe. Chem. Eng. Res. Des. 90, 40–51, http://dx. doi.org/10.1016/j.cherd.2011.08.004. Yeoh, G.H., Cheung, D.C.P., Tu, J., 2013. Multiphase Flow Analysis Using Population Balance Modeling: Bubbles, Drops and Particles, 1st edn. Butterworth-Heinemann, Amsterdam.