A Hybrid Eulerian–Eulerian Discrete-Phase Model of Turbulent Bubbly Flow

Analysis of multiphase flow systems has received great attention as a grand challenge problem in computational fluid dynamics due to its importance for a wide variety of industrial processes (cooling, energy generation, material processing, chemical reactions, etc.) [1,2]. In spite of numerous studies to date, the modeling of liquid–gas systems is still difficult due to its complexity and lack of fundamental understanding. The Eulerian–Eulerian (EE) [1,2] model has been demonstrated to provide some success in simulating practical multiphase flow problems [3–6]. However, the accuracy of this model is limited by the absence of reliable models for interphase coupling. EE models require additional help from measurements [7] or additional modeling of interfacial coupling such as bubble characteristics (size, shape, and so on) to calculate momentum interactions [8,9].

The present paper specifically addresses multiphase flows in the bubbly and transition regime to slug flow. Here, we have combined and improved a model for the bubbly flow regime within the Eulerian–Eulerian model. One of the uncertainties in bubbly EE models is the bubble size distribution [10–12]. Previous works show that there are several approaches to analyze the evolution of bubble size distribution. One of the most popular approaches, especially in the chemical engineering field, is population balance theory (PBT) [13–18]. In this model, additional transport equations for number density of several bubble sizes are solved together with the EE model, which consists of at least two sets of continuity and momentum equations. Coalescence and breakup effects are modeled as source terms in the number density transport equations as birth and death rates. Numerous works have been published to close these source terms [19,20].

In the discrete approach of PBT, the bubble size distribution is discretized into prefixed finite size groups [14,15,18]. The gas volume is exchanged among the size groups through the coalescence and breakup mechanisms. This approach is attractive since it provides an approximated bubble size distribution directly by the series of discrete size groups in the predefined range. A difficulty of using this approach in practical industrial problems is the high computational cost as the number density transport equation is a complex integro-differential equation and the required number of the transport equations increases with the number of bubble size groups tracked in the simulation. In a complex three-dimensional (3D) multiphase flow system, number of size groups can be as large as ten depending on the user [21–23].

Several different models have been developed to alleviate those issues in the approach of method of momentum (MOM) [13,24]. Instead of the number density of bubbles, this approach solves transport equations for lower-order moments of the bubble size distribution. This approach reduces the computational burden by decreasing the number of transport equations tracking in the calculation to the several moments transport equations, and avoids difficulties associated with the bubble size group selection. A popular method of this approach is the interfacial area concentration model [25]. This model expresses the bubble distribution as interface area of bubbles (second moment) and tracks it by solving a transport equation for the interface area instead of the number density.

The biggest limitation of the MOM approach is that estimation of the bubble size distribution is not available anymore because the obtained bubble size from the MOM is locally averaged (Sauter mean). Also, a closure problem arises due to the mathematical complexity of the source terms for coalescence and breakup [26], and it requires further simplifications [17] or introduction of empirical modeling [27,28].

In contrast with the above Eulerian-based approaches, there are discrete approaches using Lagrangian point-particle tracking, called discrete phase models (DPMs) [2,7]. With DPM, each particle making up the discrete phase is tracked individually as a point mass using Newton's equations of motion. Based on the interaction with the continuous phase, it is classified to be one-way coupled (only the continuous phase affects the discrete phase) [29] or two-way coupled (both phases interact with each other) [30]. Recently, four-way coupled simulations have also been introduced by including particle–particle interactions such as collisions, coalescence, and breakup [31–33]. This approach allows simulation of high gas volume fraction bubbly flows. Limitations of DPM is that a Lagrangian point-particle is not suitable to represent large-size bubbles such as Taylor bubbles or gas pockets, and also large number of bubbles because the computational cost is proportional to the number of bubbles instead of mesh size or number of bubble size groups.

There were several efforts to combine an Eulerian framework model with a Lagrangian framework model to take advantage of the both approaches (combination of EE and discrete element method [34], mixture model and DPM [35], volume of fluid and DPM [36], etc.). The primary benefit of those hybrid models is in their numerical efficiency and simplicity. A realistic bubble size distribution is easily simulated in the Lagrangian framework by solving the simple Newton's equation of motion. Thus, those hybrid models have advantages over the PBT methods when relatively small number of bubbles is simulated.

Recently, we have [37,38] combined a discrete approach for estimating the evolution of bubble sizes with a Eulerian–Eulerian model. The evolution of bubble size distribution is captured by DPM, and the bubble size information is used to calculate local momentum interactions in the EE model. In this paper, this model is applied in an upward-flowing bubbly flow in a vertical pipe and the results are validated against previous experimental work of Hibiki et al. [39].

As the name of the model suggests, the governing equations are composed of two parts: the EE model has two continuity equations (Eqs. (1) and (2)) and two momentum equations (Eqs. (3) and (4)) for the two phases to calculate velocity and volume fraction fields of each phase and a shared pressure field. Here, summation of the volume fraction of each phase becomes 1 (i.e., $αl+αg=1$⁠). The DPM equations (Eqs. (5) and (6)) track each bubble as a point mass and store its position, velocity, and acceleration. The two models are run together as separate models in the same domain, but are coupled by calculating $∑fp$ for the DPM bubbles using the Eulerian liquid phase flow field. The DPM bubbles do not affect the liquid phase field directly, but the DPM bubble size distribution influences the calculation of the EE model momentum interactions $∑Fl$ and $∑Fg$ by transferring the DPM bubble sizes to the EE model. Since this will change the liquid phase flow field in the end, this model is two-way coupled in one sense. The motivation of Eulerian–Eulerian discrete-phase model (EEDPM) is to improve the EE model results (velocities and gas volume fraction) by supplying a locally and time-varying bubble size distribution using the DPM. The EEDPM uses calculated values of gas velocity and gas volume fraction from the EE model to estimate bubble transport, bubble coalescence, breakup, and volumetric expansion. This provides a bubble size distribution and all the forces acting on the bubble, which are then included in the gas/liquid phase momentum equations. The role of DPM model is solely to calculate the evolution of bubble size distribution. It is important to decide which model is used for modeling momentum interaction terms (Eqs. (5) and (7)$)$⁠. Based on previous works [40,41], each force is modeled as follows: Same models are chosen for momentum interactions in the EE model (⁠$∑Fg$⁠) and DPM (⁠$∑fp$⁠) for the consistency of gas behavior between the EE model and the DPM model.

However, Hibiki et al. [39] observed that this sign inversion occurs when the bubble size becomes $dcr$ = 3.6 mm, not $dcr$ = 5.8 mm in a multibubble situation. In this study, this lift force model is modified to match this observation by shifting the Tomiyama lift coefficient, $CL$⁠, as shown in Fig. 1.

By solving this equation in time for the size of each DPM bubble, the gas volume change from the liquid pressure change is taken into account.

Existing models for bubble breakup and coalescence are mostly developed in the framework of PBT. Since bubbles are expressed with an Eulerian description in PBT, a suitable adjustment is required to transform these theories to a Lagrangian framework for applying them to DPM bubbles. Through this process, the complex integro-differential equations in PBT are simplified to ordinary differential equations and algebraic equations, which are more intuitive and computationally inexpensive.

Coalescence efficiency e determines the probability of coalescence: if two bubbles merge, the coalesced bubble size and velocity are determined by mass and momentum conservation. Otherwise, the two bubbles bounce apart via an elastic collision. This is a reasonable assumption for small bubbles: strong surface tension makes bubbles behave as hard spheres. Our first trial with this approach caused overestimation of coalescence since the contact time becomes too large due to low turbulent dissipation rate except near the wall. We could obtain a reasonable result with e = $10−2$⁠. A more accurate model for coalescence efficiency is required in the future.

Since the daughter bubble must satisfy all three criteria, the smaller upper bound is chosen for the maximum diameter among the energy and mass criteria. If the upper bound is greater than the lower bound, it means there is a daughter bubble size that satisfies all three criteria. Then, a bubble size is randomly picked in the diameter range with the uniform probability density function [54]. The diameter of another daughter bubble is then calculated after subtracting the first bubble volume from the parent bubble volume.

Based on measurements by Hibiki et al. [39], a vertical acrylic resin pipe with 50.8 mm diameter (D) and 3.061 m length test section is considered in this work as a test problem for the EEDPM model (Fig. 2). This test section is chosen as the model domain (Z = 0 is the domain inlet at the bottom, and Z = 60. 3D is the domain outlet at the top, which is assumed to be close to the outlet of the experimental system itself, so is set to 1 atm pressure). Measurements of velocity of both phases, gas fraction, and bubble size were taken on two measurement planes, Z = 6D and 53.5D. By assuming axisymmetric flow, a 1/6th sector of a pipe is used for the 3D domain, using a 42,000 hexahedral mesh. Because the incompressible EE models cannot account for gas expansion and the corresponding increase in gas fraction and velocity with distance up the pipe, the inlet velocity and gas fraction are specified as constants to match the average experimental data at the outlet measurement plane from Hibiki et al. [39]. Table 1 shows the inlet boundary conditions of the three cases modeled in this study.

where $δ=4$⁠, $d¯32$ is from the measurement on Z = 6D associated with the randomly determined radial injection point r, $R$ is a uniform probability density function that varies from 0 to 1. The pipe wall is assumed to be a smooth wall, and a no-slip boundary condition is used for the liquid while a free-slip boundary condition is used for the gas. A wall function for single phase turbulent flow is applied for the liquid at the pipe wall. For DPM bubbles, a reflection boundary condition is applied at the wall. Elastic collisions are assumed when the distance from the wall becomes smaller than the bubble radius. This boundary condition is important to estimate an accurate bubble size near the wall since the DPM model allows a bubble to approach the wall until its center hits the wall. Side faces of the 1/6th sector of the pipe are prescribed as symmetric boundary conditions. For the EE model and the EE part of the EEDPM model outlet, constant pressure boundary condition (P = 1 atmosphere pressure) is assumed at the outlet of the domain. For turbulence modeling, the shear stress transport $k–ω$ model is used, which transitions from $k–ε$ in the bulk to $k–ω$ near the wall through a blending function, as described elsewhere [40]. A time-step of 0.001 s is chosen for the EE model, based on satisfying the CFL number condition. $10−5$ s is chosen for the time-step of DPM calculation to resolve the bubble motion accurately. Radial Sauter mean bubble size distributions are calculated from vertically integrated DPM bubbles in five equally divided zones in height and the transiently updated radial bubble size distributions (d = d(r)) are used for the EE model momentum interaction calculation of each zone. Velocity, gas fraction, and bubble size distribution on the Z = 53.5D plane are compared to the measurements. The new hybrid model is implemented in the commercial software ANSYSfluent through new user-defined subroutines.

First, the standard Eulerian–Eulerian model is used with a constant bubble size to simulate the three cases. The original Tomiyama lift force model is used to calculate the lift force [13]. Cases 1 and 2, in the bubbly flow regime, have nearly equal gas fractions, but case 2 has twice the Reynolds number of case 1 for both phases. The average bubble sizes measured on the measurement plane Z = 53.5D are 2.6 mm for case 1, and 3.0 mm for case 2. It was observed in the experiments that small bubbles (d < 3.6 mm) moved toward the wall due to the lift force and caused peaks of gas fraction and bubble diameter at ∼2.5 mm from the wall.

Figure 3 shows velocity, gas fraction, and local bubble size distribution profiles at Z = 53.5D for case 1. We observe a reasonable agreement with the experimental data, because the local bubble sizes do not deviate much from the average prescribed input bubble size of 2.6 mm used in the simulation.

Figure 4 shows the same comparisons for case 2. These results underestimate the velocities in the core region for both phases. The gas fraction peak near the wall is overestimated slightly.

Figure 5 shows the same comparisons for case 3. The experiments observed that in case 3 the flow regime transitions from a bubbly flow to a slug flow. The big bubbles created by the coalescence process migrate toward the center of the pipe and cause a high gas fraction and large bubble sizes near the center. Due to the significant change of bubble size in the radial direction, the average bubble size of 4.0 mm cannot adjust the local bubble size properly. This causes a large deviation of the velocity magnitudes from the measurements. The gas fraction profiles are completely different. The simulations show a wall-peak profile, but measurements show a core-peak profile. One reason for this disagreement is the lift force model. The Tomiyama lift force estimates a critical bubble size for the inversion of force direction as 5.8 mm, so the lift force acts toward the wall for the 4.0 mm bubbles.

To improve the lift force calculation, a modified lift force model is used to recompute case 3. Before doing that, the radially varying bubble size distribution from the measurement at Z = 53.5D is used everywhere along the pipe instead of a constant averaged bubble size of 4.0 mm. Figure 6 shows the EE model results, with the original Tomiyama lift model. Even though the bubble size distribution is imposed to match the measurements, there is still great disagreement of the velocity and gas fraction profiles.

On the other hand, the simulation results are improved when the modified lift model is used, as shown in Fig. 7. Velocity profiles show reasonable agreement with the measurements, and the gas fraction correctly shows a core-peak profile. This result implies that modification of the lift force is necessary when there is a transition of flow regime from bubbly to slug flow in multibubble situations. Also, the EE model has potential to accurately simulate the transition regime when a proper bubble size distribution is provided.

We next recomputed cases 2 and 3 using the new EEDPM model including changes to the bubble size due to coalescence, break up, and volumetric expansion. These bubble sizes are then used to compute the interactions between the bubble and the liquid phase. Figures 8 and 9 show two snapshots of bubble distributions, illustrating two bubbles coalescing (Fig. 8) and another bubble breaking up into two bubbles (Fig. 9).

Figure 10(a) (case 2) shows the DPM bubble sizes at Z = 53.5D for case 2. In spite of the low gas fraction, collisions between bubbles happened and the coalescence effect on the bubble size distribution was not negligible. Modeling of elastic bounce after noncoalescing collisions is crucial to get a realistic bubble motion and spatial distribution of bubbles since the DPM model does not otherwise have a mechanism to avoid overlapping of bubbles. Unrealistic accumulation of DPM bubbles at the wall was observed due to the lift force if the elastic bounce effect is not included. A few breakups are observed near the wall, caused by the high turbulent dissipation rate, which decreases the Sauter mean diameter near the wall.

Figure 11 shows Sauter mean diameters of DPM bubbles at several locations for case 2 and compares them with the measurements. The injected DPM bubbles at Z = 6D plane increase in size due to both volumetric expansion and coalescence as they float upwards in the duct. It is seen that the Sauter mean diameter of the DPM bubbles passing through Z = 53.5D plane matches the measurements well. The Sauter mean diameter profile obtained from vertically integrating over the DPM bubbles in the zone including Z/D = 53.5 shows a similar trend to the results at Z = 53.5D but produces a slightly rough bubble size distribution because it is a spatial average of bubbles in the zone at simulation time t = 18 s. This radial bubble size distribution is transferred to the EE model for the calculation of local momentum interactions at each time-step.

Figure 12 shows the velocity and gas fraction from the EEDPM model and compares them with the measurements. The newly computed velocity fields agree better with the experiment and are higher near the center, but the gas fraction is overestimated near the center, and underestimated near the wall compared to the measurements and to the pure EE model results shown in Fig. 4. This error in gas fraction is due to an underestimation of the lift coefficient by the modified lift force model. This suggests that a more sophisticated lift model is needed to simulate the effect of neighbor bubbles in the case of multiple bubbles.

Figure 13 shows the distribution of Sauter mean diameters of the EEDPM bubbles for case 3 at different axial locations compared with the measurements. As shown in Fig. 10(b) (case 3), the bubble size increases significantly near the center of the pipe due to the coalescence effect. Large bubbles (d > 3.6 mm) migrate toward the center by the lift force and create even larger bubbles through coalescence. On the other hand, breakup of bubbles near the wall creates smaller bubbles and causes a decrease in Sauter mean bubble diameter. The Sauter mean diameters obtained by the DPM model at Z = 53.5D match well with the measurements. Figure 14 compares the velocity and gas fraction distributions with measurements. Compared to the results of the EE model presented in Fig. 5, the EEDPM model shows a better agreement. High gas volume fraction is predicted near the wall because small bubbles injected from the inlet and created by breakup in the vicinity of near-wall region migrate toward the wall by the shear lift force if the bubble size is less than the critical value (⁠$d<dcr$⁠).

To understand the importance of different contributions to the size change of bubbles quantitatively, parametric studies of cases 2 and 3 are conducted by numerically activating only some of the effects, which include bubble collisions, volumetric expansion, breakup, and coalescence. Table 2 lists the activated effects in each run. The collision effect is always activated to avoid the unphysical overlapping of DPM bubbles.

Figures 15 and 16 show the comparison of these effects for case 2. It turned out that the coalescence effect is larger than the expansion effect in spite of the low gas fraction (⁠$α≅6%$⁠), and both effects are not negligible. Increase of the Sauter mean diameter from Z = 6D to Z = 53.5D by coalescence was 13.1%, compared to 7.0% by expansion. The total average increase in diameter from Z = 6D to Z = 53.5D was 18.9%. Especially, the coalescence effect is important to capture the peak near r/R$≅0.95$⁠. Decrease of the Sauter mean diameter by breakup averaged 2.8%. The breakup effect is small and decreases average bubble sizes mostly near the wall.

Figures 17 and 18 display the comparison of the effects in case 3. The coalescence effect is dominant compared to other effects due to the higher gas fraction (⁠$α≅27%$⁠). The computed results with collision and coalescence roughly match the measurements at Z = 53.5D. The total increase of Sauter mean diameters from Z = 6D to Z = 53.5D of 30.0% consisted of 20.6% from coalescence and 10.9% by expansion. Decrease of the Sauter mean diameter by breakup was negligible. From the comparison with case 2, while the volumetric expansion effect is determined mostly by height (both absolute height and height difference) and initial bubble size, the coalescence effect depends greatly on local flow conditions.

A new hybrid EEDPM model for gas–liquid multiphase flow in gas–liquid systems has been developed and is tested for three cases of upward bubbly flow. The EEDPM model gives improved results compared to the EE model with a constant bubble size in both the bubbly flow regime (cases 1 and 2) and the transition regime (case 3). This is due to improved calculation of the local bubble size distribution, which evolves dynamically space and time by coalescence, breakup, and volumetric expansion, using a modified lift force relation and coalescence efficiency of 0.01. Further work is needed to improve the internal models for lift force and coalescence efficiency. The parametric studies of bubble interactions and volumetric expansion in cases 2 and 3 show that the coalescence effect is larger than other effects, and the importance is dependent on the flow conditions such as the gas fraction.

• Continuous Casting Consortium.

• Division of Civil, Mechanical and Manufacturing Innovation (15-63553).

$CD$ =

drag coefficient

$CL$ =

lift coefficient

$CTD$ =

turbulent dispersion coefficient

$CWL$ =

wall lubrication coefficient

$Cμ,l$ =

empirical constant of turbulence model

$d$ =

bubble diameter, m

D =

pipe diameter, m

$d¯32$ =

Sauter mean diameter, m

Eo =

Eotvos number

F =

cumulative probability density function

$F$ =

force, $N/m3$

$fp$ =

DPM force, $N$

$g$ =

gravity acceleration, $m/s2$

$hf$ =

final film thickness, m

$hi$ =

initial film thickness, m

$k$ =

turbulent kinetic energy, $m2/s2$

$Kkq$ =

momentum transfer coefficient from phase k to q

$mb$ =

mass of bubble, kg

$nW$ =

unit normal wall vector

$nλ˙$ =

number density of turbulent eddy, 1/$m3$

$p$ =

pressure, Pa

r =

radial position in a pipe, m

$R$ =

a uniform probability density function that varies from 0 to 1

$ReP$ =

particle Reynolds number

$t$ =

time, s

$u$ =

velocity magnitude, m/s

$u$ =

velocity field, m/s

$u¯λ$ =

turbulent eddy velocity, m/s

$vi$ =

ith DPM bubble velocity, m/s

$Vp$ =

DPM bubble volume, $m3$

$x$ =

DPM bubble position, m

$yw$ =

distance from a wall, m

$α$ =

volume fraction

$Γ$ =

gamma function

$δ$ =

parameter of Rosin–Rammler distribution

$ε$ =

turbulent dissipation rate, $m2/s3$

$η$ =

Kolmogorov eddy scale, m

$θ$ =

angular position in a pipe

$λ$ =

turbulent eddy size, m

$μ$ =

viscosity, Pa·s

$ρ$ =

density, kg/$m3$

$σ$ =

surface tension coefficient, N/m

Subscripts

b =

bubble

$B$ =

buoyancy

c =

computational cell

cr =

critical bubble size that sign inversion happens in the lift force model

D =

drag

di =

ith detached bubble

eq =

equivalent

g =

gas phase

$G$ =

gravity

i =

ith DPM bubble

l =

liquid phase

L =

lift

max =

maximum

min =

minimum

new =

new position

old =

old position

P =

pressure gradient

t =

turbulence

T =

turbulent dispersion

V =

virtual mass

W =

wall lubrication

1 =

a smaller bubble in a pair

2 =

a larger bubble in a pair

$λ$ =

turbulent eddy