Hydrodynamics of Bubble Columns : Turbulence and Population Balance Model

This paper presents an in-depth numerical analysis on the hydrodynamics of a bubble column. As in previous works on the subject, the focus here is on three important parameters characterizing the flow: interfacial forces, turbulence and inlet superficial Gas Velocity (UG). The bubble size distribution is taken into account by the use of the Quadrature Method of Moments (QMOM) model in a two-phase Euler-Euler approach using the open-source Computational Fluid Dynamics (CFD) code OpenFOAM (Open Field Operation and Manipulation). The interfacial forces accounted for in all the simulations presented here are drag, lift and virtual mass. For the turbulence analysis in the water phase, three versions of the Reynolds Averaged Navier-Stokes (RANS) k-ε turbulence model are examined: namely, the standard, modified and mixture variants. The lift force proves to be of major importance for a trustworthy prediction of the gas volume fraction profiles for all the (superficial) gas velocities tested. Concerning the turbulence, the mixture k-ε model is seen to provide higher values of the turbulent kinetic energy dissipation rate in comparison to the other models, and this clearly affects the prediction of the gas volume fraction in the bulk region, and the bubble-size distribution. In general, the modified k-ε model proves to be a good compromise between modeling simplicity and accuracy in the study of bubble columns of the kind undertaken here.


Introduction
Two-phase flow is a major topic of study in diverse research projects, due simply to the fact that it is a phenomenon occurring in many industrial flow situations, and because a rich body of information has been accrued to support advanced fluid dynamic studies of the subject for model validation.
Bubble columns in particular are widely used as test cases for experimental analysis of dispersed two-phase flows [1][2][3], since they entail simple geometric configurations, and are able to provide useful data for the validation of the CFD models [4].Table 1 provides a comprehensive list of some of the most referenced experimental studies in this area.
Both the global and local properties of the flow in bubble columns are related to the existing flow regime, which can be defined as homogeneous or heterogeneous, depending on the prevailing bubble size [5].If the bubbles are in the range of 0.001-0.007m and the superficial gas velocity is less than 0.05 m/s, the flow can be defined as homogeneous [5][6][7].
According to Besagni et al. [8], homogeneous flow (as that studied in the present work) is characterized by the gas phase being uniformly distributed in the liquid phase as discrete bubbles [9].This characterization can be further designated mono-dispersed flow or poly-dispersed flow (also known as "pseudo-homogeneous flow").Mono-dispersed flow, as the name implies, is characterized by a single bubble-size distribution and a flat local void fraction, while poly-dispersed flow involves bubbles of different sizes and a local void fraction profile with a central peak [5,8].
Most of the cited experimental papers provide local flow data (i.e., velocities and gas volume fractions) for different superficial inlet gas velocities (in homogeneous and heterogeneous ranges), but only a few provide data on the turbulence characteristics [10,11].This is unfortunate, since turbulence in the liquid phase is known to play an important role in the local gas distribution in the numerical simulation of dispersed bubbles [10].
Many numerical works have studied the effects of different turbulence models in bubble columns for constant-diameter bubbles (as seen in Table 2).As reported by Liu and Hinrichsen [12], this assumption is justified for bubble columns with low inlet gas velocity (up to 0.04 m/s), and for low void fraction.However, more advanced models are required to predict other phenomena occurring in the flow, such as coalescence and breakup of the bubbles, induced by the turbulent eddies, especially in flows where mass transfer (i.e., boiling) plays an important role independently of the inlet gas velocity [13].Consequently, bubble size distribution has normally been coupled with CFD by means of models capable of solving for the population balance.
Population balance is a continuity condition written on terms of the Number Density Function (NDF), which itself is a variable containing information on how the population of the dispersed phase, particles or bubbles, inside an infinitesimal control volume is distributed over the properties of interest (volume, velocity or length, for instance) of the continuous phase.NDF is defined in terms of an internal coordinate of the elements of the dispersed phase, i.e., particle/bubble size, volume or velocity, and in terms of the external coordinates, i.e., physical position and time [14].
The population balance equation, as defined by Ramkrishna [15], is a partial integro-differential equation that defines the evolution of the NDF (see Section 3.1).In the present work the length-based NDF, n(L; t), is considered, since our primary interest is in the evolution of the diameter of a bubble.
Different methods exist for numerically solving the population balance equation [13,14], such as the class method [16], Monte Carlo [17], the Method of Moments (MOM) [18], as well as the Quadrature Method of Moments (QMOM) and Direct QMOM (DQMOM) [19].The QMOM approach was first proposed by McGrow [18] for aerosol distributions, but it has also been shown to be suitable for modelling different dispersed flow regimes involving aggregation and coagulation issues [14].
QMOM, in contrast to the class method, solves the population balance equation by considering the concept of moments, which are obtained after integration of the NDF transport equation along the particle/bubble length, using a quadrature approximation to estimate this integral (see Section 3.1).Normally, just three quadrature nodes giving six moments, are sufficient to provide a reliable bubble distribution.This method is then generally less expensive in terms of CPU time than the discrete method [13].
In several industrial applications of bubble reactors, especially for those in which mass transfer and chemical reactions occur, bubble coalescence and break-up lead to the formation of thin films, fine drops and threads of much smaller scale than that of the average bubble [20,21].Interactions taking place at these smaller scales could eventually result in a larger range of scales needing to be considered, and treated differently.Aboulhasanzadeh et al. [20] have developed a multi-scale model that allows mass diffusion to be incorporated into Direct Numerical Simulation (DNS) of bubbly flows and predict in more detail the mass transfer between moving bubbles.This multi-scale approach, together with population balance models, can provide very deep understanding of the various phenomena occurring in multi-phase flows [20,21].
In the present work, Euler-Euler simulations are carried out for turbulent "pseudo-homogeneous" two-phase flows in a scaled-up bubble column similar to the experimental test case of Pfleger et al. [22], but with the bubble size distribution determined using QMOM.Three different superficial gas velocities, 0.0013, 0.0073 and 0.02 m/s are studied, for which bubble coalescence rather than bubble breakup is the predominant bubble interaction mechanism for the dispersed phase.It is important to mention here that the use of a population balance model for the considered inlet velocities is merely to provide a deeper understanding of the possible phenomena occurring in the dispersed flow, and not to provide any comparison with constant bubble diameter results, which has already been done by Bannari et al. [23], and Gupta and Roy [13], for instance.
Three different variants of the k-ε turbulence model are used to predict the liquid flow patterns inside the bubble column.The choice of Reynolds Averaged Navier-Stokes (RANS) models is based on the fact that they represent a simple implementation of the turbulence effects, with much lower computational effort, compared to Large Eddy Simulation (LES) and the anisotropic Reynolds Stress Model (RSM), for example, and that they are already widely used for two-phase flow simulations [24].
The main contribution of this work lies in the attempt to better understand the flow by specific investigation of the bubble diameter distribution, the gas-holdup, the axial liquid velocity, as well as the turbulent dissipation rate, which represents the coupling between the turbulence and the bubble-size distribution.This work is motivated by the need for better understanding of the hydrodynamics of gas-liquid flow in the design of bubble reactors, particularly for those who are unfamiliar with the topic of two-phase flow modeling.Very careful selection of the experimental data suitable for the validation of the models is presented, and the underlying physics explained in detail.The ultimate goal is to provide an overview of the main parameters affecting the flow patterns of the bubbles in water columns, and, through comparison with experimental data, to provide recommending options for solving similar problems.Mude and Simonin [28] Rectangular column [29]; V i = 0.20 m/s D-Vm 0.003 m -k-ε [30] and low-k-ε No significant discrepancy on the results were found between the turbulence models investigated.
Li et al. [31] Cylindrical column [32]; VSG = 0.10 m/s; D-L-Vm-WL-TD Class method Break: L-S; Coal: P-B Modified-k-ε [33] The configuration and location of gas distributor affects the overall gas hold-up and bubble sizes.
Ekambara and Dhotre [24] Cylindrical column [34]; RSM and LES provided a better prediction near the sparger, where the flow is more anisotropic.
Zhang et al. [37] Rectangular column [38]; The lift force strongly influences the bubble plume dynamics, and Tomiyana lift coefficient provided better results.

Bannari et al. [23]
Rectangular column [22]; VSG = 0.0014 m/s and 0.0073 m/s D-L-Vm Class Method Break: L-S; Coal: S-T and H k-ε [28] Good agreement with experimental data was found when using PBM for the bubble size.

Physical Model
All the numerical simulations presented here are performed for a scaled-up rectangular bubble column geometry.This configuration is chosen because comprehensive experimental data exists concerning gas hold-up and liquid velocity; those provided by Pfleger et al. [22], Buwa et al. [35] and Upadhyay et al. [11] in particular are used here for model validation purposes.As depicted in Figure 1, the column is of width (W) 0.20 m, height (H) 1.2 m and depth 0.05 m.The level of water corresponds to H w /W = 2.25, a value inspired by all the tests performed by the preceding work of Becker et al. [29].The bubbles are introduced via a sparger at the bottom.
In the experimental facility of Pfleger et al. [22], the sparger consists of a set of 8 holes placed in a rectangular configuration for the gas injection, each hole having a diameter of 0.0008 m in a square pitch of 0.006 m [13].For simplicity, the sparger in our numerical simulations is represented by a simple rectangular area placed at the bottom of the column.The same simplification has been used in other numerical works based on this test [13,39].

Governing Equations
The governing equations of mass continuity and momentum for two-phase flow with no heat transfer, following the two-fluid methodology (e.g., Versteeg and Malalasekera [44]), can be expressed as follows: where i denotes the fluid phase (liquid or gas), U i the phase velocity vector, α i the phase fraction, p the pressure, R i the turbulent Reynolds stress and M i the interfacial momentum transfer between the phases, which in this work are comprised of the drag (F D ), lift (F L ) and virtual mass (F vm ) forces.This choice follows previous studies on this topic, which showed that these three forces are enough to capture the bubble behavior in the velocity range studied here [13,37].The turbulent dispersion force has been neglected, since previous studies have already found that, for the range of superficial gas velocities considered in this work, its influence is insignificant [45,46].The drag force, also known as the resistance experienced by a body moving in the liquid [24], ultimately determines the gas phase residence time and thereby the bubble velocities, and greatly influences the macroscopic flow patterns [47].It is defined as follows: The subscripts G and L in the above equation denote gas and liquid phases respectively, and d S the Sauter mean diameter (explained in Section 3.1).The drag coefficient C D is modelled according to Schiller and Naumann [48], where Re is the relative Reynolds number defined as The lift force, which arises from the interaction between the dispersed phase and shear stress in the continuous phase, is directed perpendicular to the incoming flow, and is described as: One of the lift coefficients C L used here is from Tomiyama et al. [49] (see Equation ( 7)), which has also been used by Gupta and Roy [13], Liu and Hinrichsen [12], and Ekambara and Dhotre [24].For comparison purposes, constant values (C L = 0.14 and C L = 0.5) have also been examined.
in which, and d H given by the Wellek et al. [50] correlation: where σ is the surface tension between air and water (0.07 N/m).Finally, the virtual mass force F vm , which refers to the inertia added to the dispersed phase due to its acceleration with respect to the heavier continuous phase, is defined as Here the coefficient C vm is taken to be 0.5, based on the experimental study of Odar and Hamilton [51].
Illustrations of the three interfacial forces here considered are shown in schematic form in Figure 2, and a summary of the coefficients adopted for these forces is given in Table 3. Table 3. Coefficients of the interfacial forces used in the present work and the respective models employed.

Interfacial Force Model
Drag (C D ) Schiller and Naumann [48] Lift (C L ) Tomiyama et al. [49] and Constant (0.14 and 0.5) Virtual mass (C vm ) Constant (0.5) The Reynolds stress for phase i, needed in Equation ( 2), is given as: where I is the identity tensor, k i the turbulent kinetic energy for phase i, and ν eff ,i is the effective viscosity, composed of the following contributions: where the subscriptions l, t and L refer, respectively, to laminar, turbulent and liquid phase.The three k-ε turbulence models considered in this work calculate the effective viscosity in different ways.The standard k-ε model calculates ν eff ,L exactly as written in Equation (10).The modified k-ε includes an additional term on the right-hand-side of the equation, accounting for the turbulence induced by the movement of the bubbles through the liquid.Here, this term follows the model proposed by Sato et al. [33]: where d S is the bubble diameter (see Equation ( 35) below), U G the gas velocity, and U L the liquid velocity.
For both the standard and modified k-ε models, ν t,L is defined as: in which C µ = 0.09, k is the turbulent kinetic energy, and ε the rate of dissipation of the turbulent kinetic energy, calculated just for the liquid phase, as described by Equations ( 13) and ( 14) respectively: where 3, and P k stands for the production of turbulent kinetic energy defined, for example, as in Rusche [52]: For the mixture k-ε equations, the turbulence variables are defined as mixture quantities of the two phases, and are calculated here for both the gas and liquid phases in the manner proposed by Behzadi et al. [53]: The coefficients present in the mixture k-ε model are the same as for the standard k-ε model, while the mixture properties are defined as: where

Population Balance Model (PBM)
In this work, we specifically consider the presence of different bubble diameters.The sizes may vary due to the coalescence and break-up processes occurring, caused by the liquid turbulence.We have adopted here the QMOM model of McGrow [18] to determine the bubble-size distribution.
Assuming that both coalescence and breakup processes may occur in the bubble column, the governing equation for the population balance (length-based) can be written, according to Marchisio et al. [14], as: where n(L; t) is the length-based (L) number density function.
Applying the moment transformation and a quadrature approximation for the integral [18], results in the following equation: where the subscript φ refers to the order of each moment, varying from 0 up to (2N − 1), where N is the number of nodes considered in obtaining the weights and abscissae.The coefficients β ij and a i are, respectively, the coalescence and breakage kernels, and b (φ) i denotes the daughter size distribution.Note that specific details of the step-by-step derivation leading to Equation ( 26) can be found in the reference paper of Marchisio et al. [14].
The weights (w i ) and abscissae (L i ) in Equation ( 26) can be obtained through the product-difference approach [54] or Wheeler algorithm [55]; the first approach is adopted in the present work, since it results in a simpler implementation, and is said to be stable for N < 10 and m 1 > 0 [56], which is our case here.
The total number of nodes used is N = 3, which then provides six moments [57].Each moment m φ (t) has different physical meaning depending on the order φ.For φ = 0, m 0 is the total number of particles/bubbles; for φ = 1, m 1 is the length of the particles/bubbles; for φ = 2, m 2 is the surface area of the particles/bubbles; while for φ = 3, m 3 gives the total particle/bubble volume [57].
The coalescence between two colliding bubbles, represented by the term containing the coefficient (also referred as the kernel) β ij in Equation (26), can be expressed as the product of the collision frequency w c (d i , d j ) and the collision efficiency P c (d i , d j ), in which d i and d j are the diameters of the two colliding bubbles, as follows: The coalescence kernel used here is that proposed by Luo [41], which is based on the assumption of isotropic turbulence in the liquid phase [58].This leads to a collision frequency of the following form: where u ij is the mean turbulent velocity of the bubbles (which should be related to the size of the turbulent eddies): where β = 2.0, according to [41].
The collision efficiency P c (d i , d j ) given by Luo [41] is defined as the ratio between the coalescence time t C needed to annihilate the liquid film between the colliding bubbles to a critical rupture thickness, and to the interaction time t I between the two bubbles, as expressed below: Luo [41] considered t C to be a function of Weber number (We ij ), and t I to be a function of the virtual mass coefficient C V M , the physical properties, and the bubble size ratio.The final collision efficiency is then expressed as: where ξ = d i /d j , We ij = ρ L d i u ij /σ and σ is the surface tension.For the breakage kernel a i , we use that suggested by Luo and Svendesson [40], which has been widely applied before [13,23,[59][60][61], and which calculates simultaneously the daughter distribution (b i ) according to the following prescription: where assuming isotropic turbulence, ξ min = λ/d i is the size ratio between an eddy and a bubble in the inertial sub-range, and c f is defined as the increase coefficient of the surface area given by in which f BV is the breakage volume fraction, assumed here to be given by f BV = 0.5, to represent a symmetrical breakage.The integration in Equation ( 32) can be performed by applying an incomplete gamma function approach, as shown in the work of Bannari et al. [23].
After calculating the kernels, the Sauter mean diameter (d S ), which is needed in the momentum and turbulence equations (Equations ( 2) and ( 11)), can be obtained by dividing the third moment (bubble volume) and the second moment (bubble surface area), as suggested by Marchisio et al. [14]: The initial values of the moments may be given by a Gaussian distribution of the density function, as found by Marchisio et al. [56], which results in the following values for the moments: m 0 = 1, m 1 = 5, m 2 = 26, m 3 = 140, m 4 = 778 and m 5 = 4450.

Validation and Verification
Verification and Validation are distinct assessment processes; both have been utilized in the present work to give a measure of confidence in the numerical results obtained.The AIAA [62] has produced a guideline document on verification and validation, in which one can easily appreciate the difference between the two.
In brief, verification is defined as the process to determine whether the mathematical equations have been solved correctly.Assurance of this is obtained by comparing numerical results against analytical data or against numerical data of higher accuracy.Validation is the process of determining the degree to which a model is an accurate representation of the physics, and this can only be achieved by comparison with actual experimental data; i.e., by a simulation featuring the same set-up, boundary conditions, materials and flow conditions of the experiment.The distinction between the two was succinctly stated by Roache [63] as: "verification is solving the equations right" (meaning correctly), "while validation is solving the right equations".
Therefore, before using the models for validation purposes, the QMOM equations implemented in the OpenFOAM code were first checked by comparing the moments obtained (Equation ( 26)) with the analytical solutions of Marchisio et al. [14].As shown in Figures 3 and 4, very good agreement was obtained, so the model was accepted as being correctly implemented for the further validation studies.
Part of the verification process in this paper includes the calculation error where a mesh sensitivity analysis is performed on the test-cases with same setup as the experiment (see Section 2).The time-step here selected is 10 −3 s and remains constant over the simulation time; it respects the CFL condition of Courant number C o < 1, provides low residuals, and small computational time.

Boundary Conditions
The geometry presented in Figure 1 features one inlet at the bottom (sparger), one outlet (top of the water column), and non-slip walls.The superficial gas velocities 0.0013, 0.0073 and 0.020 m/s represent the respective inlet velocities of the air, after multiplying each value by the cross-sectional area of the column and dividing by the sparger area.In the experiment, no liquid was injected with the gas (i.e., α L = 0 and α G = 1 at the inlet), so the inlet liquid velocity is set to zero in the calculation.
The inlet conditions for k and ε are calculated, respectively, as follows: where U = α G U G + α L U L , and I the turbulence intensity, assumed here to be standard 5% (it was not measured), and in which l e is the assumed initial eddy size, calculated as l e = 0.07L c , where L c is a characteristic length, assumed to be the sparger width, and the constant 0.07 is based on the maximum value of the mixing length in a fully-developed turbulent pipe flow [64,65].
For the outlet region, a switch between a Neumann condition (zero-gradient) for outflow and a Dirichlet condition for inflow is used for k, ε, α G , U G and U L .In OpenFOAM, that is imposed through the inletOutlet condition [66].
For pressure, a boundary condition that adjusts the gradient according to the flow is used for both the inlet and the walls (given by the name f ixedFluxPressure in OpenFOAM).A Dirichlet boundary condition is assigned for pressure at the outlet boundary, using a total pressure (atmospheric) of 1.01325×10 5 Pa.
The kernels (Equations ( 27) and ( 32)) are evaluated in each cell of the domain: so, for the inlet, the bubbles are assumed not to have any initial coalescence and breakage rates associated with them.For the outlet, and the walls, a zero-gradient condition is defined for both the coalescence and breakage kernels.The same is applied to the Sauter mean diameter, d S , on the walls and outlet, but at the inlet a fixed value (0.005 m) is assumed, which is the average bubble diameter found in the experiment of Pfleger et al. [22] being simulated here; the bubbles at inlet are assumed to be of spherical shape.
Wall functions are used to estimate k and ε in the viscous layer region near the wall, allowing y+ > 30 [64,67,68].
Table 4 presents all the boundary conditions used in the simulations here discussed.It is important to point out that the simulations with inlet air velocity greater than 0.0013 m/s are initialized using the values obtained from the last time-step of the simulation with the next lower inlet air velocity, to reduce the computational time to reach steady-state conditions.

Numerical Methods
All the simulations are performed using the the open-source CFD code OpenFOAM version 2.3.0.The PIMPLE algorithm, merged PISO [69] and SIMPLE [70], is used for the coupling between the velocity and pressure fields, as recommended by Rusche [52] for bubble-driven flows.Inside the PIMPLE loop, the momentum equation is solved using an initial guess for the pressure (Step 1).Then, at each time iteration, the pressure equation is solved twice with velocity obtained from Step 1, and finally a corrector is applied to provide an updated value of the velocity based on the updated value of the pressure: this is continued to convergence.The procedure is repeated for each time step.The interested reader can obtain a more detailed explanation of the algorithm in Holzmann [71] and the OpenFOAM documentation [72].
A second-order-accurate, central-difference scheme (CDS) is used for the gradient and Laplacian terms in the governing equations, while a second-order-accurate, upwind scheme has been selected for the advection terms: this to ensure stability of the iteration.Time integration is performed using a second-order implicit backward method.The convergence criterion for all the cases is that the average residuals (for mass, pressure, velocity and turbulence variables) are less than 10 −8 .
The Euler-Euler two-phase solver is structured in such a way that first the momentum equations are solved inside the PIMPLE loop, then the turbulence quantities are calculated, and finally the PBM equation is solved to provide the Sauter mean diameter for the bubbles, which is then used in the next time-step to estimate the interfacial forces in the momentum (Equations ( 4) and ( 7)) and turbulence (Equations ( 11) and ( 23)) transport equations.
All the simulations are carried out with a constant time-step of 0.001 s, which ensures a Courant number (C o ) less than 1.For this C o criterion to be achieved, a mesh independence analysis has been performed, as explained in the next section.Under-relaxation factors ( f r ) are used for pressure ( f r = 0.3), velocity ( f r = 0.7), k L ( f r = 0.7) and ε L ( f r = 0.7), in order to promote smooth convergence, especially for the test cases with higher superficial gas velocities (0.0073 m/s and 0.020 m/s).
Time-averaging in each simulation was taken over a period of 200 s, which was sufficient for the steady-oscillation regime to be achieved.Table 5 summarizes all the test cases featured in this work.

Mesh Sensitivity Study
A mesh sensitivity analysis has been performed for three mesh sizes, the details of which are specified in Table 6.The analysis was performed assuming a constant bubble size of 0.005 m, and with the standard k-ε turbulence model.As shown in Table 6, the mesh refinement was undertaken in the transversal direction to the flow, plane xz (see Figure 1), in order to better capture the major velocity gradients, which are higher in this direction than axially, in accordance with the conclusions of Ziegenhein et al. [73] and Guédon et al. [5] in their mesh sensitivity studies.
Another important piece of information concerning the meshes it that, as discussed comprehensively by Milelli [74], an optimum ratio of the grid to the bubble diameter should be maintained of around 1.5, since larger values than this could lead to transfer of a large portion of the resolved, energy-containing scales into sub-grid-scale motion; hence in the present mesh sensitivity study a ratio less than 1.5 was specified in all cases.
The order of grid convergence (p) is here estimated based on the vertical water velocity component integrated over the column width (s i ), as provided by each i-th mesh index (see Table 6), namely by: where r = A 1 /A 2 = A 2 /A 3 = 2 is the grid refinement ratio.
In order to calculate the theoretical asymptotic solution (s 0 ) for zero grid size, also known as the Richardson extrapolation method [75], the value of p calculated from Equation ( 38) is used in the following equation: Another parameter useful in indicating the degree of mesh independence is the Grid Convergence Index (GCI), which is an estimate of the discretization error, and how far the solution is from its asymptotic value.For RANS computations, a GCI of less than 5% is considered acceptable [63,76].According to Roache [63]: Using GCI ij , one can calculate the asymptotic range (ca) of convergence: which indicates if the solution is in the error band, ca ≈1 [77].
Figure 5 shows that Mesh 2 is already in the asymptotic range, representing a good compromise between computational effort and acceptable accuracy.All the details from this mesh sensitivity study are summarized in Table 7.For all the further simulations, Mesh 2 has been the one used.

Influence of the Interfacial Forces on the Fluid Dynamics
For a better and more robust prediction of the flow patterns in a bubble column, and in two-phase flow simulations in general, an accurate description of the relevant interfacial forces is an absolute necessity.Therefore, before investigating the turbulence models in detail, the influence of the drag and lift forces on the fluid flow behaviour has first to be scrutinized.As concluded by Gupta and Roy [13], the virtual mass force does not greatly influence the fluid flow patterns for this kind of bubble column once the flow is established, so this effect has not been studied in detail in this paper.
Pourtousi et al. [47] have provided a comprehensive literature review on the topic of interfacial forces in bubble columns from which they concluded that, for the drag and lift coefficients, the most-used models are respectively those of Schiller and Naumann [48] and Tomiyama et al. [49], respectively.Consequently, this paper concentrates on these models for the prediction of the velocity and gas hold-up profiles.
It is important to point out that for the bubble-size range studied in this paper (≈0.005-0.0054m), the Tomiyama et al. [49] coefficient in the lift force is positive, and acts to drive the bubbles towards the wall.For bubbles larger than a critical diameter (0.0058 m at atmospheric pressure), the opposite occurs, and the bubbles would move towards the center of the water column [78][79][80].
As shown in Figure 6a, a near-zero velocity is found at a distance less than 0.020 m from the wall.This situation occurs for downward flow direction close to the walls, which is not the case for the interior of the column.Later, Figure 11a will confirm this behavior.
An interesting observation is that, for the velocity profiles, the combination of both models (Schiller and Naumann [48] for the drag coefficient and Tomiyama et al. [49] for the lift one) works well, providing good comparisons against data from different experiments, and from associated LES data: see Figure 6a.However, the same is not true for the gas hold-up: see Figure 6b.It is also important to note that the magnitude of the lift coefficient provided by the Tomiyama et al. [49] model is C L = 0.2 (similar observation was obtained by Kulkarni [81]).To test the sensitivity to this value, two further simulations were carried out using smaller (C L = 0.14) and higher (C L = 0.5) lift coefficients.In addition, one extra calculation is presented for which the lift force has been ignored completely, i.e., for which C L = 0.
As can be seen in Figure 6, without the lift force there is no lateral force to encourage the bubbles to migrate towards the walls (in the absence of the turbulent diffusion force), hence both velocity and gas hold-up profiles are too peaked.Including a non-zero lift coefficient, the bubbles are seen to indeed migrate to the walls, and flatten the profiles, particularly so for the gas hold-up profile.A fixed lift coefficient of C L = 0.14 appears to be the optimum choice for this particular experiment.
Previous studies on the topic [8] presented some discrepancies between Tomiyama et al. lift force [49] and the experimental data, which were explained to be originated from the Wellek et al. correlation [50].The latter was originally obtained for droplets in liquids [82] and it might not be suitable for the present poli-disperse bubbly flow [8,82,83].
In Figure 6b, the reason why their experimental data is shifted off-center, according to Buwa et al. [35], is due to the fact that, for the considered superficial gas velocity (0.0013 m/s), the bubble plume is narrow, and its oscillation period large, making it difficult to collect data over a sufficiently long time to obtain a symmetric gas hold-up profile.In the present work, the errors are within the range of the experimental uncertainties [35].
Therefore, we conclude that the lift force has an essential role to play on the prediction of the gas hold-up, which, in comparison to that of the vertical water velocity, is much more sensitive to the choice of lift coefficient.

Influence of the Turbulence Model on the Hydrodynamics
In addition to the interfacial forces, turbulence modelling is an important component in the accurate prediction of the flow patterns in bubble columns.The most used turbulence model, in two-phase flows is still the k-ε model, due to its simplicity and the low computational effort [46,47].For this reason, in this paper, we have studied the turbulence in the context of three variants of the k-ε model: standard, modified and mixture.The purpose here is to reach some conclusion on which choice provides better agreement with experimental data, at least in the context of the present application.
It transpires that both the standard and modified k-ε models produce similar profiles for axial velocity, gas hold-up and Sauter mean diameter.However, the same is not true for the mixture k-ε model, as shown in Figure 7a.
One can also observe that for all the turbulence models, the profile of the bubble size distribution can be explained by the effect of the lift force, which enhances the movement of larger bubbles towards the center of the bubble column, and the smaller bubbles towards the walls, resulting in a parabolic shape of the size distribution for the considered superficial gas velocity (0.0013 m/s).A similar conclusion was reached in the experimental work performed by Besagni and Inzoli [84].As it will be shown in the Section 4.4, as the superficial gas velocity increases, the recirculation of the flow becomes more intense, causing a spread of the bubbles all over the domain (even far from the sparger), and consequently a more homogeneous bubble size distribution (see Figures 10d and 12b).Turbulent boundary conditions affect the results mostly on the bottom wall close to the sparger, where as discussed by Ekambara et al. [24], the flow is more anisotropic and requires different modeling such as Reynolds Stress Model (RSM) and Large Eddy Simulation (LES).Away from the sparger the numerical solutions given by isotropic turbulence models, such as the k-ε model, proves to have a good agreement with the experimental data [24,29,35].
Figure 8 shows the profiles of the time-averaged turbulent dissipation rate (ε) in two different regions: one close to the sparger (y = 0.13 m), and the other close to the liquid/gas interface (y = 0.37 m).As can be observed, near the sparger (Figure 8a,c) the higher values of ε occur in the central region, for all the turbulence models here studied, as a direct consequence of the turbulence created by the insertion of bubbles into the water column.
However, in quantitative terms, the mixture k-ε model predicts ε to be one order of magnitude higher than those of the other two models in the region near the sparger.This can be explained by the high influence the model has in respect to the gas volume fraction through the ratio of the dispersed phase velocity fluctuations to the continuous phase fluctuations, as represented by the coefficient C t : see the variable definitions below Equation (23).This ratio is used in the model to obtain ε for each phase.
According to Behzadi et al. [53], for bubbly flows for which ρ L >> ρ G (such as for water-air flows), the mixture turbulence equations tend to those of the continuous phase alone, except in regions where α G ≈ 1, as happens very near the sparger, explaining the behavior of ε L in Figure 8a, but not elsewhere.
In the region close to the upper water level, where flow recirculation is a predominant mechanism, higher values of ε occur near the wall, and the difference between the mixture and the other k-ε models is even more evident (Figure 8b,d), suggesting this model is not the most adequate for a deeper analyses of this kind of bubble column.
The high values of ε provided by the mixture k-ε model near the sparger, which as a consequence dissipates more effectively the bubbles over the height of the water column, is also reflected in under-prediction of the gas hold-up, as depicted in Figure 9.The results shown in Figure 9a,b are from different test-cases, and taken at different positions, for validation purposes, once the experiments (depicted by green stars in the graphs) were obtained by different superficial gas velocities, at different locations in the bubble column.
It is worth noting that if one is only interested in obtaining velocity profiles of the continuous phase, the mixture k-ε model could be used with some confidence (as can be seen in Figure 7a).The same does not hold for the gas hold-up, which is the variable highly affected by the over-prediction of the turbulent kinetic energy dissipation rate, as shown in Figure 9.

Influence of the Air Superficial Gas Velocity at Inlet on the Fluid Behavior
As well as the specification of the interfacial forces and the turbulence modeling, the superficial gas velocity at the inlet plays an important role on the predicted hydrodynamics of the bubbles in the column.To study this influence, three different air injection velocities have been investigated: UG = 0.0013, 0.0073 and 0.02 m/s.It is to be noted that none of these generate turbulent eddies with enough energy to induce bubble break-up, as observed for UG = 0.03 m/s [85].Consequently, we focus here on the fluid dynamics for conditions under which coalescence is the dominant bubble interaction phenomenon; discussion of results for the UG = 0.03 m/s case is reserved for future work.
Figure 10 shows the gas hold-up and Sauter mean diameter profiles for two different regions of the water column, for three superficial gas velocities.As can be observed, near the sparger (y = 0.02 m) there is a higher concentration of bubbles, albeit in a narrow region, resulting in similar behavior for the bubble size distribution for all the three superficial air inlet velocities.Figure 10d shows that, for UG = 0.02 m/s, the bubbles are predicted to be of uniform size across the width of the column, and reach the top with a higher probability of interaction with each other in the central region due to the higher population density there (Figure 10c).Figure 11a shows downward flow along both walls and upward flow in the bulk region.A rotational pattern can also be observed in the middle of the column, which is stronger in the lower half of the column.This flow pattern can help explaining the fact of 0.0054-m bubbles staying in the centre region of the column, as depicted by Figure 12a.A downward flow along the left wall and on the lower portion of the right column wall can be discerned in Figure 11b, as well as a strong upward flow in the middle of the lower portion of the column, which tends to drift towards the upper-left wall.This cross-mixing behaviour may explain why the 0.0054-m bubbles spread throughout the upper portion of the column as shown by Figure 12b.
Figure 13a,b show, in shaded contour form, the behavior of the bubble volume fraction inside the water column, for UG = 0.0013 and 0.02 m/s.The greater axial velocity and stronger recirculation provided by the UG = 0.02 m/s case (Figure 11b), results in a wider spread of the bubbles across the column (Figure 13b).
For the superficial gas velocity highlighted here (UG = 0.02 m/s), the coalescence kernel (that of Luo [41]) is driven by the collision frequency (w c in Equation ( 28)), this is, in regions where ε is non-negligible the frequency of bubble coalescence is more predominant than the collision efficiency (P c ), which is proportional to the mean turbulent velocity as seen in Equation (31).A similar observation was made by Deju et al. [59], who studied different combination of coalescence and breakup kernels in a circular rather than square bubble column.In their study, they found that the collision frequency diminishes when the eddy dissipation rate becomes insignificant at the centre of the bubble column, while the opposite occurs in respect to the collision efficiency: both calculations incorporated the model proposed by Prince and Blanch [85] for the coalescence kernel.
For the given superficial gas velocities adopted in this work, a bi-dispersed model, as that proposed by Guédon et al. [5], proved to be a very good alternative to the use of a population balance model.However, if one is interested in capturing the size distribution of the bubbles, and understanding the phenomena that take place in the evolution of the bubble sizes, specially for higher inlet velocities (>0.04 m/s), use of population balance approach is recommended [85].

Conclusions
A hydrodynamic investigation of two-phase flow, with different spherical bubble size distributions, in a scaled-up rectangular bubble column geometry using CFD has been presented.Three parameters are identified as playing important roles on the fluid dynamics for this kind of flow: the gas/liquid interfacial forces, the liquid turbulence, and the superficial gas velocity at inlet.
From our examination of the interfacial forces, it is concluded that, though the use of the bubble lift force suggested by Tomiyama et al. [49] has resulted in reasonable comparisons of vertical water velocity compared with the experimental data, the same is not true for the measured gas hold-up, at least for the cases investigated here.In fact, a constant lift coefficient (C L = 0.14) proved to be the most suitable choice for better axial velocity and gas hold-up comparisons with experimental data.The lift force formulation is seen to have had a strong effect on the prediction of the gas hold-up, and consequently on the bubble size profiles, at different locations in the bubble column.
Concerning turbulence modeling, the mixture k-ε model predicted much higher turbulent kinetic dissipation rates, combined with lower values of gas volume; other variants of the k-ε model perform better.From the k-ε turbulence models here investigated, the modified k-ε model proved to be a good compromise between modeling simplicity and accuracy of predictions.
It is seen that the superficial gas velocity at inlet (UG) has an important influence on the characteristics of the fluid dynamics in the bubble column, and ultimately on the bubble size distribution within it.By increasing UG, there is an increase in the gas hold-up, even near the upper interface of the liquid, resulting in a more uniform bubble-size distribution throughout the column.
More work is needed for cases with higher injection velocities than those considered here, for which few experimental data are currently available from the point of view of deciding on the turbulence modeling approaches.

Figure 1 .
Figure 1.Sketch of the bubble column configuration: the sparger is located at the centre of the bottom face, as in the experimental facility of Pfleger et al. [22].

Figure 3 .Figure 4 .
Figure 3. Verification of QMOM equations by comparison with the analytical solution provided by Marchisio et al. [14]: (a) time history of moment 0 and (b) its respective time-history error.

Figure 5 .
Figure 5. Asymptotic averaged-axial-water-velocity behavior for a normalized grid spacing tending to zero.

Figure 6 .Figure 7 .
Figure 6.(a) Axial mean water velocity profile at y = 0.25 m; and (b) gas hold-up profiles at y = 0.37 m for the case UG = 0.0013 m/s from various sources.All the results have been obtained using the modified k-ε model, and are represented by the line with circles (No lift), the joined triangles (C L = 0.14), diamonds (C L = Tomiyama et al. [49] ≈ 0.2), and the straight line (C L =0.5).

Figure 8 .Figure 9 .
Figure 8. Time-averaged profiles of ε across the column width at (a); (c) y = 0.13 m and (b); (d) at y = 0.37 m, for a superficial gas velocity at inlet of 0.0013 m/s.

Figure 10 .
Figure 10.Time-averaged of gas hold-up (left) and Sauter mean diameter (right) provided by Modified k-ε model at (a,b) y = 0.02 m and (c,d) y = 0.37 m.

Figure 11 .
Figure 11.Time-averaged water velocity field with instantaneous water velocity vectors provided by the Modified k-ε model after 100 s of simulation at (a) 0.0013 m/s; and (b) 0.02 m/s.

Table 1 .
Summary of experimental works of turbulent flows in bubble columns

Table 2 .
Summary of numerical works of turbulent flows in bubble columns with and without bubble-size-distribution modelling.

Table 4 .
Boundary conditions for all the variables involved in the simulations of the bubble column here studied.

Table 5 .
Summary of all the cases investigated by this work.

Table 6 .
Details of the meshes used in the sensitivity analysis.
i (10 −6 m 2 ) * * V (m/s) Normalized Grid Space * * A i /A f inest ; **A i = area of the cells along the side walls in the xz plane.

Table 7 .
Convergence criteria used for choosing the appropriate mesh for the further simulations.