Computational Fluid Dynamics Modelling of Two-Phase Bubble Columns: A Comprehensive Review

: Bubble columns are used in many different industrial applications, and their design and characterisation have always been very complex. In recent years, the use of Computational Fluid Dynamics (CFD) has become very popular in the ﬁeld of multiphase ﬂows, with the ﬁnal goal of developing a predictive tool that can track the complex dynamic phenomena occurring in these types of reactors. For this reason, we present a detailed literature review on the numerical simulation of two-phase bubble columns. First, after a brief introduction to bubble column technology and ﬂow regimes, we discuss the state-of-the-art modelling approaches, presenting the models describing the momentum exchange between the phases (i.e., drag, lift, turbulent dispersion, wall lubrication, and virtual mass forces), Bubble-Induced Turbulence (BIT), and bubble coalescence and breakup, along with an overview of the Population Balance Model (PBM). Second, we present different numerical studies from the literature highlighting different model settings, performance levels, and limitations. In addition, we provide the errors between numerical predictions and experimental results concerning global (gas holdup) and local (void fraction and liquid velocity) ﬂow properties. Finally, we outline the major issues to be solved in future studies.


Introduction
Many industrial sectors have to deal with multiphase flows, especially gas-liquid flows.Bubble columns represent a class of gas-liquid multiphase reactors widely used in the chemical, biochemical, and petrochemical industries [1] (Table 1).
They offer several advantages due to the absence of moving parts or mechanical stirring equipment, such a simple and compact structure, excellent heat and mass transfer properties, low maintenance and operating costs, and high durability.
The simplest bubble column layout consists of a vertical vessel without internals, where a gas distributor at the bottom disperses the gas phase into bubbles entering the column (Figure 1).The liquid phase is supplied in batch mode, or may be co-current or counter-current to the rising bubbles.Despite the simple arrangement of bubble columns, their fluid dynamics are extremely complex, and a comprehensive understanding of their properties has not been fully achieved yet despite the significant improvements obtained over the last few years.
The overall behaviour of the flow is affected by either phenomena occurring at the macroscale (reactor-scale) or at the the local scale, as many operating variables tend to influence each other.The main parameters of interest are usually the gas holdup and the axial liquid velocity, along with the size distribution of the dispersed phase.
Historically, the design of bubble columns has relied on empirical or semi-empirical correlations:; unfortunately, this approach has severe drawbacks, as the use of correlations is constrained to specific vessel geometries, fluid properties, and operating conditions similar to those under which they were derived [2].The recent increase in interest in Computational Fluid Dynamics (CFD) has spurred substantial research efforts into determining numerical models that can obtain reasonably accurate predictions with limited computational time, thereby overcoming the limitations of traditional empirical methods.The rest of this paper is organised as follows.In Section 2, we discuss the flow regimes occurring in a vertical pipe, and more in detail in large-diameter bubble columns.In Section 3, we present numerical approaches to simulating two-phase bubble columns.In particular, within the context of the Eulerian-Eulerian approach, an overview of the population balance model used to predict the bubble size distribution follows a description of the different models used to consider the momentum exchange between the phases.Finally, models for considering the influence of the dispersed phase on the continuous phase turbulence are provided.In Section 4, we present different studies concerning the numerical simulation of two-phase bubble columns and provide their model settings and performances.

Flow Regimes
In multiphase reactors, the flow regime provides information about the behaviour of the dispersed phase and its interaction with the continuous phase.In the most general case, four flow regimes are encountered in vertical pipes when increasing the gas flow rate using a fixed system design and phase properties (Figure 2): (1) homogeneous regimes, (2) slug regimes, (3) heterogeneous (churn turbulent) regimes, (4) and annular flow regimes.The slug flow regime is not observed when dealing with industrial applications owing to Rayleigh-Taylor instabilities [5] which arise because of large diameter effects.Quantification of the Rayleigh-Taylor instabilities at the reactor scale is obtained by comparing the dimensionless column diameter (D * H ) with a critical value (D * H,cr ): In Equation (1), D H is the hydraulic diameter of the column, σ is the surface tension, g is the acceleration due to gravity, ρ L is the liquid phase density, and ρ G is the gas phase density.A bubble column is classified as having a large diameter if D * H is greater than D * H,cr = 52 [6], i.e., when d c ≥ 0.13-0.15m at ambient conditions.Under the constraint of Equation ( 1), the fluid dynamics of large-diameter bubble columns is explicated in six flow regimes and interpreted by a function of two global fluid dynamics parameters, namely, the drift flux and the gas holdup [7].
The drift flux (J T ) is defined as the volumetric flux of either component relative to a surface moving at volumetric average velocity, and is expressed as follows: In Equation ( 2), ε G is the gas holdup (defined as the ratio between the volume occupied by the gas phase and the total volume), U G is the superficial gas velocity, and U L is the superficial liquid velocity.The sign on the right-hand side depends on the operation mode of the bubble column, i.e., co-current mode (+) or counter-current mode (−); in batch mode, U L = 0.
First, the mono-dispersed homogeneous and the pseudo-homogeneous flow regimes are progressively observed; a mono-dispersed bubble size distribution characterises the former, while a poly-dispersed one in which large bubbles with a negative lift force coefficient move through the centre of the column characterises the latter.With an increase in the gas flow rate, the number of bubbles having a negative lift force coefficient increases, and these continue to migrate towards the centre of the column (Figure 4).Then, the transition flow regime (i.e., the transition flow regime without coalescence-induced structures) begins to be established by non-stable coalescence-induced structures, which causes oscillations in the gas holdup.Subsequently, the coalescence-induced structures stabilise, characterising the transition flow regime with coalescence-induced structures.Despite the presence of coalescence-induced structures, they are not entirely developed.Indeed, increasing the gas flow rate leads to a decrease in the gas holdup, thereby unveiling the imbalance between a higher gas flow rate and the formation of coalescence-induced structures, which indicates a transition process from the prevailing fluid dynamics.
Finally, stable coalescence-induced structures emerge, leading to the pseudo-heterogeneous and purely heterogeneous flow regimes.The former is characterised by an equilibrium between the increase in the flow rate (in terms of drift flux) and the formation of coalescenceinduced structures.This leads to constant holdup with respect to an increase in the drift flux, indicating that the coalescence-induced structures are well established.On the contrary, in the pure-heterogeneous flow regime increasing the gas flow rate increases coalescenceinduced structures, and the increase in drift flux is followed by an increase in the gas holdup; for more details, see the paper by Besagni (2020) [7].

Numerical Modelling: The Eulerian-Eulerian Approach
Two main models have been developed to predict the complex fluid dynamics phenomena of a dispersed bubbly flow, namely, the Eulerian-Lagrangian and Eulerian-Eulerian frameworks.The Eulerian-Lagrangian model couples the Eulerian description of the continuous phase with a Lagrangian scheme for tracking the individual particulates.The dynamics of the surrounding fluid (continuous phase) are solved through the governing equations, while the particulates (dispersed phase) are tracked independently through the surrounding fluid by computing their trajectory.While a variety of studies have applied the Eulerian-Lagrangian approach [4,[9][10][11], these have been limited to the simulation of small-scale reactors and low gas holdup scenarios.In the Eulerian-Eulerian model, both phases are treated as interpenetrating continua and an ensemble average is performed on the balance equations, reducing the computational cost of the simulation.The information about the precise location of the bubble interface is lost, and only an average description of the flow is possible; however, this approach is not limited by the number of bubbles to be modelled.Because industrial flows usually have high gas volume fraction, meaning that a very large number of bubbles are present, only the Eulerian-Eulerian method is suitable for CFD implementations.

Governing Equations
In the Euler-Euler method, all phases share a single pressure field, whereas the momentum and continuity equations must be solved for each phase.Considering an isothermal flow without mass transfer, the governing equations for the k-phase read as follows: where τk is the stress-strain tensor for the k-th phase, defined as Here, µ k and λ k express the shear and bulk viscosity of the k phase, respectively.The other terms on the right-hand side of Equation ( 4) represent the pressure gradient and the interfacial momentum exchange between the phases, accounting for all the interaction forces exchanged between the phases (see Section 3.2).

Interfacial Forces
To properly solve the k-th equation, a feasible set of closure relations must be implemented to include all the possible interactions between primary and secondary phases, expressed in the form of momentum transfer per unit volume at the phase interface (Figure 5).Interfacial momentum forces are typically added as source terms in the momentum equation, and can be divided into drag and non-drag components.
The non-drag forces term comprises different physical mechanisms: lift, turbulent dispersion, wall lubrication, and virtual mass.Thus, the total interfacial force responsible for momentum exchange is provided by a linear combination of each individual force [12].Indicating the continuous liquid phase with the subscript "L" and the dispersed gas phase with the subscript "G", the source term in the momentum equation reads as follows:

Drag Force
The drag force reflects the resistance opposing the bubbles' motion relative to the surrounding liquid.Due to rising upward motion through the continuous liquid phase, the dispersed phase particles are accelerated by buoyancy and slowed down by pressure effects and viscous friction at the interface.As the drag force accounts for such decelerating phenomenon, it is always applied in the opposite direction with respect to the bubble motion; therefore, it is flow region-dependent, and should be formulated with the inclusion of local flow parameters such as the relative slip velocity, particle diameter, and void fraction [12]: In Equation ( 7), C D is the drag coefficient, which depends on the bubble size, shape, flow field, and contaminations.Therefore, it is a function of the relative bubble Reynolds number (Re b ), Eötvös number (Eo), and Morton number (Mo) [12], defined as follows: For spherical bubbles without deformation and fully contaminated systems, the functional form of the drag coefficient depends only on the bubble Reynolds number, as determined by Schiller and Naumann (1935) [14].They modified the Stokes model, originally formulated only for creeping flows (Re b << 1), by adding a correction factor to extend its validity up to Re b = 1000.A constant value of 0.44 was considered in order to extend the applicability of the model: This model is suitable only for small bubbles, and is commonly used in bubble columns operated in the homogeneous flow regime.
Similarly to Shiller and Naumann (1935) [14], White (1974) [15] proposed a drag coefficient correlation highlighting the dependency of the drag coefficient only on the bubble Reynolds number: For deformed bubbles, the drag coefficient is affected by the Eötvös and Morton numbers.Grace et al. (1976) [16] were among the first to publish a drag law suitable for deformed bubbles: where In Equation (15), the bubble terminal velocity (U t ) is calculated as follows: where J is provided by the piecewise function and In Equation ( 19), µ ref is set to 0.0009 kg/(m•s).
A model with the same functional form as Equation ( 13) was proposed by Ishii and Zuber (1979) [17], who distinguished different shape regimes and found the following: Except at high values of Eötvös number, this correlation agrees well with experimental results concerning the terminal velocity of bubbles rising in quiescent liquids.Grevskott et al. (1996) [18] experimentally studied the drag coefficient considering the bubble size distribution in pure water by dividing the population of bubbles into small and large bubbles.Fitting the experimental data, they proposed the following drag coefficient correlation: A popular correlation suitable for gas-liquid flows with a range of shapes was derived by Tomiyama et al. (1998) [19].The proposed drag coefficient consists of three equations, respectively corresponding to pure, slightly contaminated, and fully contaminated systems: This drag correlation has a wide range of applicability (10 −3 < R b < 10 6 and 10 −2 < Eo < 10 3 ) and holds for spherical as well as deformed bubbles.Zhang et al. (2006) [20] derived a simple relation that provides results comparable to other models when the bubble diameter is in the range between 4 mm and 10 mm, for which the rising velocity is weakly dependent on the bubble diameter: Recently, Dijkhuizen et al. (2010) [21] used a DNS front-tracking technique to derive a correlation valid for both spherical and deformed bubbles: where The proposed drag closure matches very well with experimental data for ultra-pure liquids (water/glycerine mixture).A comparison between the different drag coefficient correlations is presented in Figure 6.Schiller and Naumann (1933) [14], Grace et al. (1976) [16], Ishii and Zuber (1978) [17], Grevskott et al. (1996) [18] and Tomiyama (1998) [19].

Swarm Factor Definitions
The drag coefficient correlations presented in Section 3.2.1 apply to isolated bubbles rising in initially quiescent liquids, and can only be used for low gas volume fractions.In the case of high gas volume fractions, as is common in industrial bubble columns with high superficial gas velocity, bubbles are very close to each other, resulting in a crowding effect in which bubbles travel through the column in swarms.The bubble boundary layers interact, changing the interphase momentum exchange and, importantly, the drag force.For this reason, a correction factor (h) should be applied to the drag coefficient for isolated bubbles where C D is the actual drag coefficient and C ∞ D is the drag coefficient for an isolated bubble.Several authors have investigated the impact of bubble-bubble interactions on the drag coefficient.The first approach was proposed by Bridge et al. (1964) [22] and Wallis (1969) [23]: for which Bridge et al. (1964) [22] selected n = 1.39 and Wallis (1969) selected [23] Ishii and Zuber (1979) [17] investigated the drag force over a broad range of gas volume fractions and Reynolds numbers, establishing the following relation based on a large sample of experimental data: Rusche and Issa (2000) [24] proposed a swarm factor correlation for oblate bubbles, validated for a local volume fraction α G ≤ 0.45 and 27 ≤ Re b ≤ 960: Roghair et al. (2011) [25], using a DNS front-tracking model, derived a drag-modification term that can be applied for high Reynolds numbers and oblate bubbles in the range Buffo et al. (2016) [26] developed a straightforward power law equation to account for the crowding effect of bubbles based on many investigations into fluidized beds and liquid-liquid systems: The authors concluded from a sensitivity analysis that a value of C A = 1.3 provides the best agreement with experimental data.
The swarm factor correlations discussed above hinder the rising velocity of the bubbles, thereby increasing the drag coefficient of isolated bubbles at high gas volume fractions.They are reported in Figure 7 and compared with the Simmonet et al. ( 2008) [27] model; while reasonable agreement is found at low gas volume fractions, the predictions diverge strongly at high volume fractions.The correction factor proposed by Simmonet et al. (2008) [27] differs from the other models in highlighting the presence of a critical value of the gas volume fraction at 0.15.Below this value, the rising velocity of the bubbles decreases with increasing void fraction; on the contrary, beyond the critical value the wake acceleration effect of large bubbles is predominant, causing an increase in the bubbles' rising velocity and a consequent reduction in the swarm factor, which decays rapidly to zero.For local void fractions less than 0.30 and air-water pure systems, the correlation of Simonnet et al. (2008) [27] reads as follows: where m is an adjustable coefficient.Simonnet et al. (2008) [27] set m = 25, as this showed good agreement with their experimental data.
McClure et al. (2014c) [28] suggested that no hindering effect should be taken into account for bubbles larger than 7 mm, and modified the correlation of Simonnet et al. (2008) [27] accordingly: where h ′ is the original swarm factor of Simonnet et al. (2008) [27].
Performing a new experimental study, McClure et al. (2017) [29] criticized their previous model and observed that the drag correction factor is a function of both the bubble size distribution and the local gas volume fraction.Thus, they proposed a new empirical correlation valid for α G > 0.25: where n and b are model constants calculated with least-square fitting of experimental data.The authors suggested n = 50, while b depends on the gas sparger design (for example, a value of 0.20 was found in the case of a perforated plate distributor with a 0.5 mm hole diameter, while a value of 0.08 was found in the case of a perforated plate with a 3 mm hole diameter).(2008) [27] (validated for α G < 0.3) to obtain a swarm factor feasible for a wide range of operating conditions.In particular, they noted that the decrease in the swarm factor proposed by Simonnet et al. (2008) [27] is too strong for high volume fractions; thus, they added a minimum asymptotic value (h min ) towards which h should tend at high void fractions: The authors suggested h min = 0.5 for industrial scaling-up purposes.Yang et al. (2018) [31] formulated a drag modification relation for bubble swarms that accounts for both the hindering effect of small bubbles and the wake accelerating effect of large bubbles.To model both phenomena, the resulting swarm factor was calculated using the product of two corrective terms: where f b, large is the fraction of large bubbles and α G, small = α G 1 − f b,large .From this model equation, it is possible to observe that an increasing void fraction means that the drag correction factor increases to a maximum and then decreases.Drag correction factors that reduce the drag force in bubble swarms are presented in Figure 8.

Lift Force
A bubble moving through a shear fluid is subjected to a force perpendicular to its motion.This force results from the action of pressure and stress on the surface of the bubble.The corresponding momentum source is where C L is the lift coefficient.Experimental and numerical studies show that the coefficient is positive for small bubbles and negative for large bubbles.Consequently, the lift force direction depends on the bubble size and shape.For small bubbles, the lift force acts in the direction of decreasing liquid velocity, which in the case of batch or co-current mode is toward the pipe wall.Conversely, when large bubbles are considered, a force assimilated to the lift force pushes the bubbles towards the centre of the column.In other words, sign reversal of the lift force is observed for bubbles with a diameter larger than a critical value (d cr ).Smaller bubbles (those with a diameter lower than the critical threshold) tend to move towards the wall, whereas larger bubbles (those for which d b > d cr ) are pushed towards the vessel centre.Tomiyama et al. (2002b) [32] studied the motion of bubbles in a glycerol-water solution, deriving the following correlation for the lift coefficient: where here, Eo ⊥ is the modified Eötvös number calculated considering the maximum horizontal dimension (major axis) of the bubble (d ⊥ ), provided by the empirical correlation for the aspect ratio (E) proposed by Wellek et al. (1966) [33] for contaminated flows: where d eq is the equivalent spherical diameter of the bubble.This correlation applies to larger-scale deformable bubbles in the ellipsoidal and spherical cap regimes.The experimental conditions under which Equation ( 45) was derived were limited to the ranges of −5.5 ≤ log 10 (Mo) ≤ −2.8 and 1.39 ≤ Eo ≤ 5.84.Air-water systems at normal conditions usually have Mo = 2.63 × 10 −11 ( log 10 (Mo) = −10.55),which is quite different from the conditions at which the above correlation was derived; nevertheless, good results have been reported for this case.When coupled with that of Wellek et al. (1966) [33] for the calculation of d ⊥ , the correlation predicts that the change in sign of the lift coefficient occurs at d b = 5.8 mm.
However, Ziegenhein and Lucas (2019) [34] pointed out that the correlation of Wellek et al. (1966) [33] is not suitable for purified air-water systems and leads to an overprediction of the major axis of an elliptic bubble.Besagni et al. (2016) [35] studied the size distribution and shapes of bubbles in an annular gap bubble column and reached the same conclusion.Correct estimation of the critical diameter is essential in CFD modelling, as it allows the bubble population to be divided into a small group and large group, permitting consideration of the poly-disperse nature of the flow in which bubbles of different sizes move in different directions.
To overcome this limitation, Ziegenhein et al. (2018) [36] suggested a quadratic fit correlation for the lift force as a function of the modified Eötvs number and bubble Reynolds number: In their model, the correlation for the bubble major axis is modified based on experimental observation of bubble shapes in six column configurations: The Ziegenhein et al. ( 2018) [36] correlations results in a critical diameter of about 5.3 mm, which is lower than that obtained with the Wellek et al. (1966) [33] relation.Figure 9 shows a comparison between the Tomiyama et al. ( 2002) [32] and the Ziegenhein et al. (2018) [36] lift models.
To conclude, Behzadi et al. ( 2004) [37] proposed a model highlighting that the lift coefficient cannot be independent of the local gas void fraction.After fitting their experimental observations, they proposed the following correlation: It should be pointed out that this simple correlation does not consider the change of sign of the lift coefficient (i.e., C L does not depend on Eo ⊥ ) and, due to the power law formulation, the coefficient tends to infinity for α G −→ 0. This is avoided by limiting Cl by a finite value, which is taken as 0.25 [37].

Turbulent Dispersion Force
The turbulent dispersion force accounts for the fluctuations in the liquid velocity that affect the dispersed phase by scattering it.The turbulent eddies redistribute the bubbles in the lateral direction from regions with high concentration to regions with low concentrations.The effect of the turbulent dispersion force is evident in the gas fraction profiles, as it modulates peaks of small bubbles near the wall pipe and spreads out large bubbles.
The modelling of the turbulent dispersion force is expressed as a force proportional to the void fraction gradient.Simonin and Viollet (1990) [38] proposed the following model for the turbulent dispersion force: In analogy to the molecular diffusion, σ TD is the Schmidt number (a value of σ TD = 0.9 is typically used) and C TD is the turbulent dispersion coefficient, while k LG and D t,LG are the mixture's turbulent kinetic energy and the liquid dispersion scalar, respectively.Burns et al. ( 2004) [39] derived an explicit formulation for the turbulent dispersion force by using Favre averaging on the drag force: This model estimates the dispersion scalar of the Simonin and Viollet (1990) [38] model with the liquid turbulent viscosity µ turb L .An alternative formulation is the one proposed by Lopez De Bertodano et al. (2004) [40]: where St is the Stokes number and u ′ L u ′ L represents the normal Reynolds stresses.

Wall Lubrication Force
The wall lubrication force is caused by the surface tension and prevents the bubbles from touching the walls, ensuring zero void fraction (as found experimentally) near the vertical walls.
The wall lubrication force has the general form where ( u G − u L ) || is the relative velocity component parallel to the wall, ˆ y is the unit normal to the wall pointing into the fluid, and C W L is the wall lubrication coefficient.
Different expressions for C W L have been proposed in the literature.The relation proposed by Antal et al. (1991) [41] states that With C w1 and C w2 being two non-dimensional coefficients.The Antal et al. (1991) [41] model is applicable only if the condition y w ≤ 5d b is respected; thus, it is only active on a sufficiently fine mesh.Tomiyama et al. (1995) [42], studying bubble trajectories in a glycerol-water solution, proposed a different function for the wall lubrication coefficient by considering an additional dependence on the Eötvös number and pipe diameter (D): Although the Tomiyama (1995) [42] model has proven superior to the Antal et al. (1991) [41] model, it is restricted to flows in pipe geometries.Hosokawa et al. (2002) [43] reviewed the experimental results of Tomiyama (1995) [42] and obtained the following functional form for C W L : Here, the coefficient C w depends on the Eötvös number and the phase-relative Reynolds number (Re b ), as indicated below: The Hosokawa et al. ( 2002) [43] model often shows better agreement with experimental data in both vertical and inclined bubbly flow systems.
Later, Frank et al. (2005) [44] proposed a generalized formulation for the wall lubrication coefficient: where C WC and C WD are the cut-off and the damping coefficient, respectively.The Eötvös number coefficient C w is provided by Equation (58).From numerical simulations, the au-thors found that good agreement with the experimental data was obtained for C WC = 10, C WD = 6.8, and p = 1.7.
Figure 10 shows the influence of the above described non-drag forces on the numerical prediction of the local void fraction profile [45].

Virtual Mass Force
The virtual mass force arises from the acceleration of an immersed moving object relative to the surrounding fluid.This effect is significant when the continuous phase density is much higher than that of the dispersed phase.
Here, D j /Dt, D k /Dt denote the total derivative with respect to time for the indicated phase.As reported by Tabib et al. (2008) [46], the contribution of virtual mass is negligible for large-diameter bubble columns (D > 0.15 m).A common practice is to set a constant value C V M = 0.5, which has been obtained experimentally for isolated spherical bubbles, as the virtual mass coefficient.

Population Balance Modelling
Most of the interfacial correlations defined in Section 3.2 require the equivalent diameter of the bubbles as an input; in this regard, three approaches can be used to model the dispersed phase: (1) the mono-dispersed model, (2) the homogeneous Multiple Size Group model (MUSIG), and (3) the inhomogeneous Multiple Size Group model (iMUSIG).The mono-dispersed approximation considers two separate continuity and momentum equations for the gas and liquid phases.The MUSIG model applies to poly-dispersed systems, with the bubble population divided into two groups with different sizes.This model solves a separate continuity equation for each size group and a single momentum equation for all sizes together.Conversely, the iMUSIG model considers several velocity groups and solves a separate momentum equation for each [47].
An uncertainty in bubble column design and scale-up arises from a lack of understanding of the phenomena governing bubble size distribution (BSD), and consequently the column fluid dynamics and interfacial area.Bubble size distributions depend on many phenomena occurring in the bubble column, including the balance between the coalescence and breakup rate, the pressure change, and the gas-liquid mass transfer.In particular, the bubble size distribution mainly results from the coalescence and breakup phenomena occurring in the bubble column.
Consequently, the dispersed phase characteristic properties vary in time and space, meaning that a constant bubble size distribution represents a strong assumption.To account for the changes in the particle population, a balance equation can be included in addition to the mass, momentum, and energy balances within the contest of MUSIG and iMUSIG models.This concept is referred to as the Population Balance Model (PBM).
A population balance model can be described by a transport equation reflecting particles entering or leaving a control volume via several mechanisms.The bubble number density function, called the population balance equation (PBE), reads [48]: where n( x, V b , t) is the bubble density distribution function that specifies the probable number density of bubbles at a given time (t) in the spatial range (d( x)) about a position x for a bubble volume between V b and V b + d(V b ) at time t.Finally, S( x, V b , t) is the following source term: where S b , S c , S ph , S p , S m , and S r are the bubble source/sink terms due to breakup, coalescence, phase change, pressure change, mass transfer, and reaction, respectively.The breakup source term (S b ) reads where the first term is the birth rate of bubbles of volume V b due to the breakup of bubbles with volumes larger than V b , the second term is the death rate of bubbles of volume V b due to breakup, m(V ′ b ) is the mean number of daughter bubbles produced by the breakup of a parent bubble with volume V ′ b , g(V ′ b ) is the breakage frequency (the fraction of particles of volume V ′ b breaking per unit time), and β(V b , V ′ b ) denotes the daughter distribution factor of a particle of volume V b breaking from a parent bubble of volume V ′ b .The coalescence source term (S c ) reads where the first term is the birth rate of bubbles of volume V b due to the coalescence of bubbles of volume V b − V ′ b and V ′ b , the second term is the death rate of bubbles of volume V b due to coalescence with others bubbles, and Γ(V b , V ′ b ) is the coalescence rate between bubbles of volume V b and V ′ b .Neglecting the interfacial mass transfer, the reaction contribution, the phase variations, and the changes due to pressure gradients, the coalescence and breakup phenomena are the only items of interest.Considering a steady state, it follows that Which can be rewritten in terms of bubble diameter as where S c and S b are determined by modelling the breakup and coalescence rates.

Bubble Breakup Phenomena Modelling
The breakup of fluid particles in the dispersed phase depends on the hydrodynamics of the continuous phase and interface interactions.The phenomenon of rupture can be expressed as a balance between two opposite effects, namely, the external stresses of the continuous phase, which attempt to destroy the bubble, and the surface stress of the particle plus the viscous stress of the fluid inside it, which opposes the deformation of the bubble.Therefore, breakage is caused by the hydrodynamic conditions in the surrounding liquid and the properties the bubble itself [49].
Turbulent fluctuations and collisions, in which breakage is mainly caused by turbulent pressure fluctuations along the surface or by particle-eddy collisions.The dominant external force is the dynamic pressure difference around the bubble, meaning that the breakage process can be studied as the balance between the dynamic pressure and surface stresses.

2.
Viscous shear stresses, which cause a velocity gradient around the interface that can deform or break the bubble.In addition, wake effects may be responsible for the formation of shear stresses.Breakage can be modeled as the balance between external viscous stresses and surface tension forces as expressed by means of the capillary number C a .

3.
Shearing-off processes, which occur when small bubbles are sheared off from a large bubble through erosive breakage.4.
Interfacial instabilities, which arise in the absence of a continuous phase net flow.If a significant density difference is present, as in the case of a light liquid accelerated into a heavy fluid, Rayleigh-Taylor instabilities are found; on the other hand, if the density ratio is close to unity Kelvin-Helmholtz instabilities exist.
In general, the most common causes of breakage are associated with turbulent fluctuations and the other mechanisms are often neglected.Modelling of breakup phenomena requires definition of the breakage frequency (g) and daughter size distribution function.Applicable breakup models for bubble columns are reported in Table 2.

Bubble Coalescence Phenomena Modelling
According to Chesters (1991) [56], the coalescence process is more complex than breakup, as it involves interactions of bubbles with the surrounding liquid phase as well as those between bubbles themselves when they are brought into contact by the external flow or body forces [48].Coalescence can be modelled by using empirical or physical models.The former are much easier to implement in CFD simulations, as they typically take the form of power law functions with adjustable coefficients.However, because these models depend on the experimental facility geometrical aspects, the results cannot be extended to other setups.The latter models attempt to describe the underlying physics of the phenomenon by identifying specific physical parameters capable of detailing the bubble collision mechanism.In this approach, the coalescence kernel is determined as the product of the collision frequency (Γ(d b , d ′ b )) and coalescence efficiency λ(d b , d ′ b ): where d b,i and d b,j are the diameters of the two colliding bubbles.The collision frequency represents the number of collisions between bubbles per unit of time.For a turbulent bubbly flow, there are five mechanisms of collision (Figure 12): 1.
Turbulence-induced collisions occur as a result of the random motion of bubbles caused by turbulent fluctuations.

2.
Viscous shear-induced collisions are generated by global liquid velocity gradients, meaning that bubbles in a location of high liquid velocity may collide with those in a region of low liquid velocity.

3.
Eddy-capture, in which collision events are produced by the shear rate of the flow in the turbulent eddy.4.
Bouyancy-induced collisions, in which collisions may occur because bubbles of different size have different rising velocities. 5.
Wake entrainment collisions, which may result when bubbles are accelerated by the wake region behind a large spherical-cap bubble.Thus, the overall rate of collision events is determined by the sum of the effects of all the listed mechanisms.Typically, only turbulence-induced collisions are assumed to be predominant.The others are frequently insignificant except at very high superficial gas velocities (where, for example, wake entrainment effects cannot be ignored due to the formation of large spherical-cap bubbles).
Bubble collisions do not always lead to coalescence, as bubbles may bounce off each other or separate after the collision.Hence, the efficiency term (h(d b , d ′ b )) is introduced in the definition of the coalescence kernel.Three models for calculation of the coalescence efficiency have been introduced: 1.
Film drainage model: first proposed by Shinnar and Church (1960) [57], this is a widely used model of coalescence efficiency.When two bubbles collide, a thin liquid film is trapped between their surfaces and is progressively drained.If the contact time is sufficiently high, the liquid film reaches a minimum thickness, then ruptures, causing coalescence.

2.
Energy model: first proposed by Howarth (1964) [58], this model is based on the concept of collision energy, in which a higher collision energy indicating a higher probability of coalescence.

3.
Critical approach velocity model: in this model, collisions result in coalescence phenomena if the approach velocity of the bubbles exceeds a certain critical value; otherwise, they bump into or bounce off of each other, and do not coalesce.
Several different coalescence frequency and efficiency models have been proposed; the ones most commonly used in bubble column simulations are presented in Table 3.

Solution Methods
Many numerical schemes are available in the literature to solve the population balance in Equation ( 63).These can be classified into four main categories: 1.
Class (or Sectional) Method (CM): the internal coordinate domain is divided into a finite number of intervals (or bins), transforming the population balance equation into a set of balance equations in the physical domain.Any coalescence and/or breakup event is accompanied by the migration of particles from one class to the adjoining classes.The advantages of this method are its robust numerics and that it computes the Particle Size Distribution directly [62].

2.
Monte Carlo Method: this method solves the population balance equation based on the statistical ensemble approach, accurately tracking particulate changes in a multidimensional system.Nevertheless, the method accuracy strongly depends on the number of simulation particles, and requires an extensive computational time to track large numbers of particles.This makes the Monte Carlo method poorly compatible with the conceptual framework of Computational Fluid Dynamics [47].

3.
Standard Method of Moments (SMM): the population balance equation is turned into a set of transport equations for the moment of the particle size distribution.The primary advantage is numerical economy, as it is sufficient to solve a limited number of moment equations.Mathematically, the transformation from the space of particulate size distribution to the space of moments is extremely rigorous, and fractional moments, representing the Sauter mean diameter of the bubbles, present a serious closure problem [47].This closure constraint can be overcome by resorting to a Quadrature Method of Moments (QMOM) approach.4.
Quadrature Method of Moments (QMOM): first suggested by McGraw (1997) [63] for modelling aerosol and coagulation problems, QMOM was later applied by Marchisio et al. (2003) [64] for solving the population balance equation, becoming an attractive alternative.Compared to the SMM method, this approach solves only the transport equations of the low-order moments; however, it is able to overcome the closure problem of the SMM method [64].With the QMOM, the integral terms in the momentum transport equations are approximated by employing an N-node Gaussian quadrature formula.This quadrature approximation requires knowledge of N weights (ω i ) and N nodes of abscissas L i , and determines a sequence of polynomials orthogonal to the unknown number density function.The functional form (for a univariate problem with φ as internal coordinate) reads as follows: where the weights (ω i ) and abscissas (L i ) are determined through the Product-Difference (PD) algorithm from the lower-order moments [64].When the weights and abscissas are known, the source term due to coalescence and breakup can be calculated and the transport equation for moments can be solved.Finally, starting from the moments of the distribution, it is possible to solve the inverse problem of reconstructing the bubble size distribution [62].
When dealing with a multivariate population balance equation, for which the productdifference algorithm can not be applied, other extensions of QMOM are available, such as the Conditional Quadrature Method of Moments (CQMOM) or the Direct Quadrature Methods of Moments (DQMOM).In DQMOM, the transport equations are directly solved for the weights and nodes of the quadrature approximation, whereas CQMOM represents the multivariate extension of QMOM.Moreover, in CQMOM closure is achieved by means of multivariate quadrature approximation, and the transport equations for the moments of the distributions are solved [65].
From the above discussion, it is clear that for univariate problems (which is the case for bubble columns when the mass transfer is not modelled) the class method and quadrature method of moments are the only viable candidates for solving the population balance equation.
In the class method, a high number of bubble classes (i.e., 12-18) are required to describe the evolution of the particle population accurately, which incurs high computational costs.On the contrary, in the quadrature method of moments approach it is sufficient to solve a limited number of moments (i.e., 4-8), yielding a reduction in the dimensionality of the problem compared to the class method, in turn leading to lower computational costs [62].

Turbulence Modelling
Turbulence modelling in bubbly flows has predominant importance due to its significant influence on the local distribution of the dispersed phase and in determining the coalescence and breakage phenomena.The description of the effects of turbulent fluctuations of velocities and other scalar quantities in multiphase simulations is complex, as the number of terms to be modelled in the momentum equations is larger compared to the case of single-phase flows.Moreover, turbulent mixing occurs in bubble columns over a wide range of length and time scales [66]; the largest turbulent scales are of the same order of magnitude as the characteristic scales of the mean flow, and depend on the geometry and boundary conditions, the smaller turbulence scales depend on bubble dynamics and are proportional to the bubble diameter, and the smallest scales are associated with the Kolgomorov scale and are responsible for the dissipation of the turbulent kinetic energy.Considering the importance of turbulence, a variety of numerical methods have been developed to address this issue, which can be grouped into the three following categories: (1) Reynolds-Averaged Navier-Stokes (RANS) models; (2) Large Eddy Simulation (LES); and (3) Direct Numerical Simulation (DNS).
In RANS turbulence models, attention is focused on the mean flow characteristics; the instantaneous flow equations are time-averaged through a Reynolds decomposition approach applied to both pressure and velocity components in the continuity and momentum equations.Six extra terms involving the products of fluctuating velocity components, called Reynolds stresses, appear after performing this time-averaging process.Appropriate closure relations, namely, k − ε and k − ω, are required to solve the system of mean flow equations, making use of the Boussinesq hypothesis based on the assumption of isotropic turbulence.This approximation considers the Reynolds stresses proportional to the mean rates of deformation of the fluid elements.The generic component of the viscous stress tensor τ for a Newtonian fluid results in the following expression: where δ ij is the Kronecker delta and k = 1 2 u ′2 + v ′2 + w ′2 is the turbulent kinetic energy per unit mass.The first term on the right-hand side of Equation ( 71) contains the turbulent eddy viscosity µ k .For k − ε models, the turbulent viscosity is provided by where C µ is a dimensionless constant, often taken as 0.09.With RANS models, the computational resources needed for reasonably accurate flow simulations are modest, making them the mainstay of many engineering flow calculations.Multiphase turbulence RANS models are derived directly from their single phase equivalents.The isotropic turbulence assumption for the core of the flow in RANS models is not strictly satisfied for bubble columns, as the velocity fluctuations in the gravity direction are normally twice those in the other directions [67].As a result, the turbulence generated by the rising bubbles is mainly anisotropic.RANS models fail to exactly represent the directional effects of the Reynolds stress field; the Reynolds Stress Model (RSM) can be used to address these problems adequately, as it incorporates anisotropy by solving the transport equation of each Reynolds stress component.This increases the computer storage and run time, and as a result has found little success in the simulation of bubbly flows.Nonetheless, RSM seems a promising tool to provide more realistic numerical prediction of turbulence quantities that can be used as inputs in bubble coalescence and breakage models [68,69].
A further option is the use of Large Eddy Simulation (LES), an approach that uses an unsteady flow computation to resolve only the turbulent structures that have a length scale larger than a certain cutoff, then describes the smaller eddies and their effects on the larger scales by means of Sub-Grid-Scale (SGS) models.By resolving only large portions of the turbulent motion, LES appears to be more suitable for accounting for anisotropic structures in bubble columns; however, this comes at the cost of a much larger volume of calculations compared to RANS methods.Another issue is selection of the proper mesh size, which when using LES should be bounded within a certain range in order to obtain a correct filter cutoff width; a very dense grid is commonly required.Thus, LES may be unable to consistently represent the correct sub-grid-scale stress for various flow situation [46].
DNS does not include any turbulence modelling, and directly resolves the instantaneous Navier-Stokes equations while considering all turbulence scales of motion and even the fastest fluctuations.Unfortunately, its enormous computational and memory storage requirements make its applicability to industrial-scale reactors unfeasible.

Bubble-Induced Turbulence
Experimental observations have found that the sources of the bulk liquid turbulence can be divided into two categories, namely, (1) Shear-Induced Turbulence (SIT) and (2) Bubble-Induced Turbulence (BIT).The former is independent of the bubble size, whereas the presence of gaseous bubbles generates the latter.When a bubble rises in a pool of liquid in a turbulent flow regime, a portion of its pressure energy is converted into liquid phase turbulence, then into internal energy at the level of the Kolmogorov scale [70].
Within the context of the RANS k − ε and k − ω models, it is common practice to account for the influence of the dispersed phase on the liquid phase turbulence by adding extra source terms (S k,L , S ε,L , or S ω,L ) to the transport equations for the turbulence quantities.A plausible approximation is that all energy lost by bubbles due to the drag force is converted to turbulent kinetic energy (k) in the wake of the bubble [71].Hence, the ksource becomes The ε-equation source term is obtained by dividing the k-equation source term by a characteristic time scale (τ p ): The ε-source term can be transformed into into an equivalent ω-source term as follows: Modelling of the time scale (τ p ) proceeds largely based on dimensional analysis [71].The time scale factors commonly used in numerical simulation of bubble columns are presented in Table 4.  Morel (1997) Unlike the model presented in Table 4, Sato and Sekoguchi (1975) [77] did not explicitly add source terms in the turbulence transport equations, rather, they tried to incorporate the effects of bubble-induced turbulence in the turbulent viscosity, proposing the relation

Literature Survey
Tables 5-8 collect different studies regarding CFD simulations of bubble columns validated with experimental data; Tables 5-7 refer to the operating conditions, physical settings, and numerical settings, respectively.Conversely, Table 8 focuses on model performance.To derive a consistent comparison between different studies, the values in Table 8 have been computed by extracting published values using a consistent procedure.The CFD accuracy is assessed by computing the relative error (δ) between numerical (ω CFD ) and experimental (ω EXP ) data, as follows: When local validation wis performed with N local quantities (i.e., local void fraction profile and local liquid velocity profile), the error (δ) is computed as follows: Pfleger and Becker (2001) [74] performed 3D transient simulations of a 0.288-m inner diameter bubble column operating in the homogeneous flow regime using a mono-dispersed approach.Concerning the interfacial forces, they considered only the drag force with a constant drag coefficient C D = 0.44.The standard k − ε model was used for turbulence in the liquid phase.Their main aim was to study the effect of BIT on the correctness of CFD predictions.The local liquid velocity and local void fraction showed the positive impact of the BIT model, with δ = 30.29%(δ = 51.75% without BIT) and δ = 4.21% (δ = 11.14%without BIT), respectively.The simulation results, including BIT or not, were found to over-predict the global gas holdup, which the authors attributed to the simplified gas inlet modelling.
Chen et al. ( 2005) [78] carried out 2D axisymmetric transient simulations in a 0.19-m inner diameter bubble column operating in the heterogeneous flow regime.The main focus was studying the influence of different coalescence and breakup models.In particular, they compared the coalescence efficiency of Prince and Blanch (1990) [50], Chesters (1991) [60], and Luo (1993) [61] and the breakup kernels of Luo and Svendsen (1996) [51] and Martinez and Bazan (1999) [52,53].The authors found that different bubble breakup and coalescence closures did not significantly impact the simulation results.However, the best agreement with the experimental values was provided by the breakup model of Luo and Svendsen (1996) [51] coupled with the coalescence efficiency of Prince and Blanch (1990) [50] concerning the local void fraction (δ = 22.74%) and the breakup model of Martinez-Bazan [52,53] coupled with the coalescence efficiency of Luo (1993) [61] regarding the local liquid velocity (δ = 51.44%).
Ekambara et al. ( 2010) [79] compared several different turbulence models: standard k − ω, standard k − ε, k − ε RNG, RSM, and LES.They considered a mono-dispersed approach including all the non-drag forces; turbulent dispersion was considered only for the RANS methods.Concerning the local void fraction profile, LES provided the best agreement with the experimental data in the near-sparger region (δ LES = 18%, δ RSM = 29.55%,δ k−εRNG = 31.44%,δ k−ε = 34% and δ k−ω = 36.24%),where the flow is more anisotropic.The same was found for the local liquid velocity (δ LES = 32.08%,δ RSM = 53.59%,δ k−εRNG = 57.06%,δ k−ε = 46.05%,and δ k−ω = 50.65%).No remarkable differences were found in the fully-developed region, indicating that the RANS approaches perform well when the objective is to understand the steady and time-averaged features of the flow.
Ziegenhein et al. (2013) [45] simulated a cylindrical bubble column operating in the homogeneous flow regime considering mono-dispersed and iMUSIG (two bubble classes without coalescence and breakup) approaches.Moreover, they analysed the influence of non-drag forces, with the exception of the virtual mass, which was not included in the simulations.Considering a superficial gas velocity of 0.15 cm/s, the mono-dispersed and poly-dispersed approaches with non-drag forces performed similarly, with δ = 22.02% and δ = 20.89%,respectively.The poly-disperse approach neglecting the non-drag force, despite the error being slightly lower than the other cases (δ = 18.79%), predicted an overly strong centre-picked void fraction profile.When increasing the superficial gas velocity, the poly-disperse models with and without non-drag forces provided similar results (at U G = 0.5 cm/s, δ without NDF = 32.41%, and δ with NDF = 15.56%,respectively).However, the model including the non-drag forces performed better near the walls; conversely, neglecting the non-drag forces resulted in better predicting the void fraction in the column centre.
Ziegenhein et al. ( 2013) [71] carried out transient simulations considering non-drag forces and a swarm factor, fixed poly-dispersity (modelled using the iMUSIG model), and BIT.In particular, they compared the BIT models of Morel (1997) [72], Troshco (2001) [73], Politano (2003) [75], Rezehak (2012) [76], and Sato (1975) [77].Regarding the void fraction profile at U G = 0.3 cm/s, the results provided by the mono-dispersed model without swarm factor matched very well with the experimental data (δ = 4.19%), and the inclusion of BIT did not improve model prediction.When increasing the superficial gas velocity to U G = 1.3 cm/s, the mono-dispersed treatment was not able to reproduce the gradient of the gas volume fraction near the wall, which was better reproduced by the poly-dispersed treatment (δ MONO = 14.42% and δ POLY = 6.48%;BIT of Rzehak (2013) [76]).The same was found for the superficial liquid velocity (δ MONO = 44.27%and δ POLY = 30.20%;BIT of Rzehak (2013) [76]).Regarding the different BIT models, the Rzehak (2013) [76]) model performed better than the others.Xu et al. (2013) [80] simulated a large-diameter bubble column operating in the heterogeneous flow regime using mono-dispersed, MUSIG, and iMUSIG approaches.In the latter two cases, coalescence and breakup were considered by applying the coalescence model of Luo (1993) [61] and the breakup kernel of Luo and Svendesen (1996) [51].Moreover, they studied the influence of the lift force by adopting the correlation of Tomiyama et al. (2002) [32].The Shiller and Naumann (1933) [14] drag coefficient was used in the single-size model, and the Ishii and Zuber (1979) [17] drag coefficient was used for the MUSIG and iMUSIG simulations.None of the mono-dispersed, MUSIG, and iMUSIG models neglecting the lift force were able to reproduce the radial distribution of the local gas holdup, with δ = 21%, δ = 21.44%, and δ = 21.53%,respectively.When the lift force was included, the best agreement with the experimental data was provided by the iMUSIG model with the drag coefficient of Ishii and Zuber (1979) [17] (δ = 9.68%).No remarkable differences between the models were found in predicting the local liquid velocity profile.
McCLure et al. ( 2013) [81], simulating a cylindrical bubble column with a monodisperse approach, investigated the influence of the mean bubble diameter, lift force, and BIT.Considering the void fraction profile at a superficial gas velocity of U G = 8 cm/s and a bubble diameter of d b = 4 mm, the inclusion of the lift force resulted in a better agreement with the experimental data (δ LIFT = 15.37% and δ NO LIFT = 23.31%).The lift force slightly influenced the results when d b = 6 mm (δ LIFT = 20.46% and δ NO LIFT = 22.29%).Regarding the BIT, the CFD prediction of the void fraction profile was in reasonable agreement with experimental data when the BIT was incorporated into the model (δ BIT = 20.80%,δ NO BIT = 75.80%).
Masood at al. (2014) [82] studied the hydrodynamics of a square bubble column and tested different drag closures.In particular, they considered the models of Shiller and Naumann (1933) [14], Grace et al. (1976) [16], and Ishii and Zuber (1979) [17], with a constant drag coefficient of C D = 0.44.Turbulence was modelled using the k − ε RNG model, including BIT of Sato (1975) [77].The constant drag coefficient was the worst in predicting the global gas holdup curve (δ = 22.75%), followed by the simplistic correlation proposed by Schiller and Naumann (1933) [14] (δ = 14.76%).The Ishii and Zuber (1979) [17] model was found to be slightly superior to that of Grace et al. (1976) [16], with δ = 6.18% and δ = 9.45%, respectively.Liu et al. (2014) [83] performed 2D-axisymmetric simulations of a bubble column operating in the heterogeneous flow regime with a MUSIG model considering coalescence and breakup.They included only the drag (Tomiyama et al. (1998) [19] model with the swarm factor of Ishii and Zuber (1979) [17]) and lift forces (Tomiyama et al. (2002) [32] and Behzadi et al. (2004) [37] models); the other forces were neglected.For the turbulence, they compared the standard k − ε and RSM models.The authors performed a sensitivity analysis on the number of bubble classes through 10 and 20 classes for the k − ε and RSM models with BIT.They found that increasing the number of classes did not significantly increase agreement with the experimental data; for example, concerning the void fraction profile and considering the k − ε turbulence model, δ = 8.66% and δ = 5.17% for 10 classes and 20 classes, respectively.However, much more computational time was needed when increasing the number of bubble classes to 20.Regarding the comparison between the lift correlations, a much better agreement with the experimental void fraction profile was achieved using the Tomiyama et al. (2002) [32] lift model.For example, considering the RSM turbulence model, δ TOMIYAMA = 9.60% and δ BEHZADI = 27.70%.Such an increment in the relative error can be explained by the fact that the model of Behzadi et al. (2004) [37] does not predict the change in the sign of the lift coefficient, resulting in a flat profile of the void fraction.No remarkable differences were found between the k − ε and RSM turbulence models.
Syed at al. (2017) [84] carried out 2D-axisymmetric simulations of a bubble column operating in the homogenous flow regime.They performed a sensitivity analysis on the Luo (1993) [61] coalescence kernel by changing the coalescence parameter c λ ; see Table 3 from 0.1 to 1.1.Moreover, they compared the breakup kernels of Luo and Svendsen (1996) [51] and Lehr et al. (2002) [55].The results showed that the radial profiles of void fraction and local liquid velocity were significantly affected by these parameters.Considering the local void fraction profile at U G = 0.38 cm/s and the Luo and Svendsen (1996) [51] breakup kernel, δ c 0 =0.1 = 10.5% and δ c 0 =1.1 = 23.14%;consequently, a lower coalescence parameter (c λ ) led to better prediction of the experimental data.The same was found for the global gas holdup and local liquid velocity profile.Comparing the two breakup models, it was found that the Lehr et al. (2002) [55] model improved prediction of the local void fraction as compared to the Luo and Svendsen (1996) [51] model (δ c λ =0.3,LUO and SVENDSEN = 13.3% and δ c λ =0.3,LEHR = 3.65%).Moreover, when the breakup model of Lehr et al. (2002) [55] was considered, the results regarding the local void fraction profiles were slightly influenced by c λ ; on the contrary, the authors found that a decrease in the coalescence parameter negatively influenced the prediction of the local liquid velocity profiles.
Gemello et al. (2018) [30] proposed an interesting modification of the Simonnet et al. ( 2008) [27] swarm factor that was suitable for very high volume fractions while avoiding stability problems encountered in the original formulation.The results obtained with the proposed swarm model matched very well with the experimental data.Considering the gas holdup, the relative error was δ = 2.83% (δ = 29.82%with the Simonnet et al. ( 2008) [27] swarm factor).The authors obtained good prediction of the local void and liquid velocity profiles, with a relative error at U G = 16 cm/s of δ = 4.07% and δ = 25.18%,respectively.
Zhang et al. (2020) [85] studied the influence of non-drag forces and BIT by simulating a cylindrical bubble column operating in the heterogeneous flow regime with a MUSIG model coupled with the coalescence kernel of Luo (1993) [61] and the breakup model of Luo and Svendsen (1996) [51].Considering the local void fraction profile, good agreement with the experimental data was obtained considering all the non-drag forces and BIT (δ = 8.07%).When the lift force was not considered, the relative error increased to δ = 42.69%,indicating that the lift force could not be neglected.Furthermore, neglecting the turbulent dispersion and wall lubrication forces reduced the model's predictive capacity, increasing the relative error with respect to the experimental data ( δ = 23.26% and δ = 12.69%, respectively).No evident differences were found when neglecting the BIT and virtual mass force.
As can be observed from Table 8, the accuracy of CFD models has been determined by comparing the numerical results against global gas holdup, local void fraction, and in a few cases the local liquid velocity, without providing any information about model performance with respect to turbulent quantities.Indeed, measuring the local liquid velocity and turbulence in large-scale bubble columns with optical methods is complex, and is usually limited to low gas holdups and thin geometries.To overcome this issue, Deen et al. (2000) [86] used Particle Image Velocimetry (PIV) and Laser Doppler Velocimetry (LDA) to study turbulent quantities in a square bubble column.Similarly, Ziegenhein et al. (2019) [87] designed a Particle Shadowgraph Velocimetry (PSV) technique and tested it in a large-diameter bubble column to provide data for CFD model validation.
The main conclusion that emerges from the above literature survey is that several different methods have been applied when dealing with the numerical simulation of twophase bubble columns.As can be noted from Tables 2 and 3, different strategies in terms of physical models, interfacial momentum exchange, turbulence models, and numerical settings have been adopted, with the result that while certain models perform better than others, no general conclusion can be reached.Moreover, available numerical studies are often limited to air-water systems, with the consequence that a model which has outstanding performance on air-water systems may not be suitable in other cases.Consequently, fully predictive CFD models for bubble columns are far from being developed, and more effort must be spent in order to overcome the lack of knowledge regarding numerical simulation of bubble columns.

Conclusions
This paper has presented a detailed review of CFD modelling of two-phase bubble columns.The models most commonly adopted for interfacial momentum exchange, bubbleinduced turbulence, coalescence, and breakup are described.A quantitative comparison between the different modelling approaches is presented in Section 4, considering various studies from the literature and computing the relative errors between the CFD predictions and the experimental data.
On the basis of the literature review presented in this work, the following recommendations are made for future studies: 1.
Concerning the interfacial forces, the momentum transfer between the phases is dominated by the drag force.For a proper description of the drag coefficient, the models of Tomiyama et al. (1998) [19], Grace et al. (1976) [16], and Ishii and Zuber (1978) [17] can be implemented; however, they should be corrected with a swarm factor.When presenting numerical studies, a sensitivity analysis among the different models should be performed.
The lift, wall lubrication, and turbulence dispersion forces should be added to the model to obtain more accurate solutions.In particular, a correlation that predicts the change in the sign of the lift coefficient should be considered.

2.
Concerning the turbulence modelling, RANS models, in particular the k − ε RNG and the k − ω SST models, provide satisfying results in terms of average quantities.The LES turbulence model provides better results in the near-sparger region, where the flow is more anisotropic.However, no remarkable differences compared to the RANS methods have been highlighted in the fully developed region.

3.
The modelling approach of the dispersed phase (i.e., mono-dispersed, MUSIG, iMUSIG, PBM) should always be related to the simulated flow regime.A mono-dispersed approximation applies at very low superficial gas velocities.Conversely, multiple-size models that include coalescence and breakup should be considered.4.
Regarding the numerical settings, high-order resolution discretization schemes should be used in order to prevent or reduce numerical discretization errors as much as possible.
Finally, it can be concluded that detailed comprehension of the phenomena governing multiphase flows remains needed, and additional effort should be spend on developing a fully predictive a priori CFD model.

Figure 2 .
Figure 2. Representation of flow regimes in vertical pipes.Reprinted with permission from [4].

Figure 3 .
Figure 3. Flow regimes and regime transitions in a large-diameter bubble column.Reprinted with permission from [7].

Figure 4 .
Figure 4. Flow pattern evolution in two-phase bubble columns.Reprinted with permission from [8].

Figure 5 .
Figure 5. Schematic representation of the interphase forces acting on a bubble.Reprinted with permission from [13].

Figure 10 .
Figure 10.Influence of Non-Drag-Forces (NDF), (with the exception of the virtual mass force) on the local void fraction profile.Reprinted with permission from [45].
i and d b,j refer to the parent and daughter bubble diameters, respectively [m]; λ ε is the turbulent eddy size [m].

Table 2 .
Applicable breakup models for simulated bubble column.

Table 3 .
Applicable coalescence models for bubble column simulations.

Table 4 .
Time scale factors commonly used in numerical simulation of bubble columns.

Table 5 .
Literature survey: column characteristics and operating conditions.
* In the case of a non-cylindrical column, depth and width are provided.** Modelled as a uniform inlet.*** Modelled with mass source points
-: information not provided.