The Eulerian–Eulerian two-fluid model (EE) is a powerful general model for multiphase flow computations. However, one limitation of the EE model is that it has no ability to estimate the local bubble sizes by itself. In this work, we have combined the discrete phase model (DPM) to estimate the evolution of bubble sizes with the EE model. In the DPM, the change of bubble size distribution is estimated by coalescence, breakup, and volumetric expansion modeling of the bubbles. The time-varying bubble distribution is used to compute the local interface area between gas and liquid phase, which is then used to estimate the momentum interactions such as drag, lift, wall lubrication, and turbulent dispersion forces for the EE model. In this work, this newly developed hybrid model Eulerian–Eulerian discrete-phase model (EEDPM) is applied to compute an upward flowing bubbly flow in a vertical pipe and the results are compared with previous experimental work of Hibiki et al. (2001, “Axial Interfacial Area Transport of Vertical Bubbly Flows,” Int. J. Heat Mass Transfer, 44(10), pp. 1869–1888). The EEDPM model is able to reasonably predict the locally different bubble size distributions and the velocity and gas fraction fields. On the other hand, the standard EE model without the DPM shows good comparison with measurements only when the prescribed constant initial bubble size is accurate and does not change much. Parametric studies are implemented to understand the contributions of bubble interactions and volumetric expansion on the size change of bubbles quantitatively. The results show that coalescence is larger than other effects, and naturally increases in importance with increasing gas fraction.

## Introduction

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 [36]. 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 [1012]. 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) [1318]. 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 [2123].

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 [3133]. 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].

## Model Description: Governing Equations

The current model solves the following equations for the liquid and gas phases:
$∂(αgρg)∂t+∇·αgρgug=0$
(1)

$∂(αlρl)∂t+∇·αlρlul=0$
(2)

$∂αgρgug∂t+∇·αgρgugug=−αg∇p+∇·(μgαg∇ug+∇ugT)+αgρgg+αg∑Fg$
(3)

$∂αlρlul∂t+∇·αlρlulul=−αl∇p+∇·(μlαl∇ul+∇ulT)+αlρlg+αl∑Fl$
(4)

$ρgVpdvidt=∑fp=fDp+fLp+fWp+fTp+fVp+fPp+fBp$
(5)

$dxidt=vi$
(6)

$∑Fg=−∑Fl=FD+FL+FW+FT+FV+FP$
(7)

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.

For the drag force, the Tomiyama drag model [42] is chosen. This model considers deformation of bubble shape by including the Eotvos number $Eo$ in the drag coefficient calculation. Thus, it can be applied to a wide range of bubble shape regimes such as spherical, ellipsoidal, and spherical caps. The drag law is given as
$FD=−34CDdρl|ug−ul|ug−ul$
(8)

$fDp=−34CDdρlVp|vi−ul|vi−ul$
(9)

$CD=maxmin24Rep1+0.15Rep0.687,72Rep,83EoEo+4$
(10)
where
$Eo=gρl−ρgd2σ,ReP=ρlug−uldμlfororEEorReP=ρlvi−uldμl(forDPM)$
(11)
The lift force is important for lateral migration of bubbles. It is known that bubbles migrate differentially depending on their size. Large bubbles have more chance to be deformed due to their smaller surface tension forces and the substantial deformation changes the lift force direction. The Tomiyama lift model [43] captures this sign inversion of lift force at bubble diameter $dcr$ = 5.8 mm based on the bubble shape through the Eotvos number
$FL=CLρlug−ul×∇×ul$
(12)

$fLp=CLρlVpvi−ul×∇×ul$
(13)

$CL=min0.288tanh0.121ReP,fEo′(Eo′≤4)=fEo′(4
(14)
where
$Eo′=gρl−ρgdh2σ,dh=d1+0.163Eo0.75713$
(15)

$fEo′=0.00105Eo′3−0.0159Eo′2−0.0204Eo′+0.474$
(16)

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.

A wall lubrication force model [44] is introduced to account for hydrodynamic forces near the wall. Basically, this force always pushes bubbles away from the wall so that bubbles are kept detached from the wall. For small bubbles, the wall lubrication force acts in opposite direction to the lift force: the balance between these lateral forces, i.e., lift and wall lubrication forces, plays a key role in determining the radial gas fraction profile. The Hosokawa wall lubrication force model is used in this work
$FW=CWLρlug−ul2nW$
(17)

$fWp=CWLρlVpvi−ul2nW$
(18)

$CWL=Cwd21yw1D−yw2$
(19)

$Cw=max7Rep1.9,0.0217Eo$
(20)
A turbulent dispersion force model [45] is also included to consider a force from turbulent fluctuations. Burns et al. [45] derived this force through the Favre average of drag force. Numerically, this force smooths out the gas fraction field
$FT=−34CDd|ug−ul|μt,lσlg(∇αgαg−∇αlαl)$
(21)

$fTp=−34CDdVp|vi−ul|μt,lσlg(∇αgαg−∇αlαl)$
(22)
The effect of bubbles on the turbulence of the liquid phase is modeled as a source term of turbulent viscosity [46]
$μt,l=Cμ,lkl2εl+0.6αgd|ug−ul|$
(23)
For transient forces, the virtual mass force and pressure gradient force are added [46,47]. The virtual mass force (Eqs. (24) and (25)) is an additional force required to accelerate the surrounding fluid when the bubble is accelerated. The pressure gradient force (Eqs. (26) and (27)) arises when there is a nonuniform pressure distribution around the bubble
$FV=0.5ρl(DulDt−DugDt)$
(24)

$fVp=0.5ρlVp(DulDt−dvidt)$
(25)

$FP=ρlDulDt$
(26)

$fPp=ρlVpDulDt$
(27)
The buoyancy force for a DPM bubble is calculated as follows:
$fBp=gVp(ρg−ρl)$
(28)

## Volumetric Expansion

Gas bubbles expand as they rise according to the surrounding liquid pressure field. To calculate the bubble size change due to liquid pressure changes, a cubic equation with respect to $dnew$ (new bubble diameter) is derived from the Young–Laplace equation and the ideal gas law
$pl,newdnew23+2σdnew22−pg,olddold23=0$
(29)

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.

## Bubble Interaction Modeling

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.

In the coalescence model, the collision frequency in PBT is easily handled from the calculation of distances between a pair of DPM bubbles. Here, only the distances between pairs located in the same computational cell are calculated to decrease the computational cost from $n2ton$ ($n$ = number of bubbles). Once the distance between the pair is smaller than the sum of the radii of two bubbles, the pair is counted as a collided pair. And then, the coalescence efficiency e is estimated through calculations of drainage time and contact time. According to Prince and Blanch [48], the drainage time ($tdrainage$) for the liquid film that forms between the collided pair and the contact time ($tcontact$) of the pair are calculated as follows:
$e=exp(−tdrainagetcontact)$
(30)

$tdrainage=deq3ρl128σ0.5ln(hihf)$
(31)

$tcontact=(deq/2)eq23ε13$
(32)

$deq=21/d1+1/d2$
(33)

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.

For the breakup model, the works of Luo and Svendsen [49] and Wang et al. [50] are chosen. In Luo and Svendsen's theory, a bubble breaks up when it meets an eddy that has size smaller than the bubble [51], but enough kinetic energy to create a new surface caused by breakup. To decide which eddy size hits a bubble, an eddy size is randomly determined in the range of $λmin<λ<λmax$ based on the eddy size distribution function. In Luo and Svendsen's work, the eddy size distribution is derived from the number density of eddy $nλ˙dλ=c3(1−αg)/λ4dλ$
$Fλ=numberdensityofeddythathassizeλmin∼λtotalnumberdensityofeddies=∫λminλnλ˙dλ∫λminλmaxnλ˙dλ=λmax3λmin3λmax3−λmin3(1λ3−1λmin3)$
(34)

$λmax=D6,λmin=11.4×η$
(35)
Here, $λmax$ and $λmin$ stand for the maximum and minimum eddy sizes in the inertial subrange since the number density expression is derived under an assumption of isotropic turbulence in inertial subrange. Once an eddy size is determined randomly with the cumulative probability density function $Fλ$, it is checked that the eddy size is smaller than the bubble diameter. If so, the range of diameter for a daughter bubble is calculated based on mass, force, and energy balance criteria. The mass balance criterion is based on criterion that the daughter bubble cannot be larger than the parent bubble. The energy balance criterion is obtained from the balance of eddy kinetic energy and surface creation energy as follows:
$ρlVλu¯λ22≥πd12σ+πd22σ−πd2σ=cfπd2σor$

$d1≤ρlVλu¯λ22cfπσ$
(36)
Breakup happens when the eddy has enough kinetic energy so that the energy is enough to create new interface area for the daughter bubbles. Here, $u¯λ=2.00.5ελ13$ [52,53] is assumed. For the force balance criterion, balance between inertial force (dynamic pressure) and surface tension force (capillary pressure) is considered as follows:
$ρlu¯λ22≥σd1orσρlu¯λ22≤d1$
(37)
Luo and Svendsen's theory had only an energy criterion and it encountered an unphysical result such that tiny bubbles are created extensively compared to experiment since it does not have any restriction for the minimum bubble size. To improve this model, Wang et al. [50] added the minimum bubble size from the force balance criterion previously. By combining all three criteria, we get
$dmin≤d1≤dmaxorσρlu¯λ22≤d1≤min(d,ρlVλu¯λ22cfπσ)$
(38)

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.

## Numerical Model Setup for Benchmark Test Problem

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.

Cases 1 and 2 are typical bubbly flows, and case 3 is in a transition between a bubbly and a slug flow regime. DPM bubbles are injected at Z = 6D, based on the measured bubble sizes at that plane with gas flow rates taken from measurements at Z = 6D, and inlet pressure from the EE model. The injection point ($r,θ$) on the cross section of Z = 6D is randomly determined for each bubble as follows:
$r=D2×R,θ=π3×R$
(39)
where $R$ is a uniform probability density function that varies from 0 to 1. These distribution functions make a uniform distribution on a fan-shape cross-sectional area. Once the injection point is determined, the injected bubble size is determined by the radial bubble size distribution from measurement data ($d¯32=d¯32(r)$). Instead of using the Sauter mean bubble size $d¯32$ from the measurement directly, a Rosin–Rammler size distribution is assumed for the injected bubble sizes
$d=d¯32×ln11−R1δ$
(40)

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.

## Eulerian–Eulerian Model Results With Tomiyama Lift Force

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.

## Eulerian–Eulerian Model Results With Modified Lift Force

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.

## Eulerian–Eulerian Discrete-phase Model Results

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).

## Evaluation of Bubble Interaction and Expansion Effects

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.

## Conclusions

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.

## Funding Data

• Continuous Casting Consortium.

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

## Nomenclature

• $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

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 =

•
• 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

## References

References
1.
Harlow
,
F. H.
, and
Amsden
,
A. A.
,
1975
, “
Numerical-Calculation of Multiphase Fluid-Flow
,”
J. Comput. Phys.
,
17
(
1
), pp.
19
52
.
2.
Riley
,
J. J.
,
1974
, “
Diffusion Experiments With Numerically Integrated Isotropic Turbulence
,”
Phys. Fluids
,
17
(
2
), p.
292
.
3.
Monahan
,
S. M.
,
Vitankar
,
V. S.
, and
Fox
,
R. O.
,
2005
, “
CFD Predictions for Flow-Regime Transitions in Bubble Columns
,”
AlChE J.
,
51
(
7
), pp.
1897
1923
.
4.
Sanyal
,
J.
,
Vasquez
,
S.
,
Roy
,
S.
, and
Dudukovic
,
M. P.
,
1999
, “
Numerical Simulation of Gas-Liquid Dynamics in Cylindrical Bubble Column Reactors
,”
Chem. Eng. Sci.
,
54
(
21
), pp.
5071
5083
.
5.
Sokolichin
,
A.
, and
Eigenberger
,
G.
,
1994
, “
Gas-Liquid Flow in Bubble Columns and Loop Reactors—Part I: Detailed Modelling and Numerical Simulation
,”
Chem. Eng. Sci.
,
49
(
24
), pp.
5735
5746
.
6.
Law
,
D.
,
Battaglia
,
F.
, and
Heindel
,
T. J.
,
2006
, “
Numerical Simulations of Gas-Liquid Flow Dynamics in Bubble Columns
,”
ASME
Paper No. IMECE2006-13544.
7.
Law
,
D.
,
Jones
,
S. T.
,
Heindel
,
T. J.
, and
Battaglia
,
F.
,
2011
, “
A Combined Numerical and Experimental Study of Hydrodynamics for an Air-Water External Loop Airlift Reactor
,”
ASME J. Fluids Eng.
,
133
(
2
), p.
021301
.
8.
Ma
,
J.
,
Oberai
,
A. A.
,
Drew
,
D. A.
, and
Lahey
,
R. T.
,
2012
, “
A Two-Way Coupled Polydispersed Two-Fluid Model for the Simulation of Air Entrainment Beneath a Plunging Liquid Jet
,”
ASME J. Fluids Eng.
,
134
(
10
), p.
101304
.
9.
Ishii
,
M.
, and
Mishima
,
K.
,
1984
, “
Two-Fluid Model and Hydrodynamic Constitutive Relations
,”
Nucl. Eng. Des.
,
82
(
2–3
), pp.
107
126
.
10.
Grevskott
,
S.
,
Sannaes
,
B. H.
,
Dudukovic
,
M. P.
,
Hjarbo
,
K. W.
, and
Svendsen
,
H. F.
,
1996
, “
Liquid Circulation, Bubble Size Distributions, and Solids Movement in Two- and Three-Phase Bubble Columns
,”
Chem. Eng. Sci.
,
51
(
10
), pp.
1703
1713
.
11.
Krishna
,
R.
,
van Baten
,
J. M.
, and
Urseanu
,
M. I.
,
2000
, “
Three-Phase Eulerian Simulations of Bubble Column Reactors Operating in the Churn-Turbulent Regime: A Scale Up Strategy
,”
Chem. Eng. Sci.
,
55
(
16
), pp.
3275
3286
.
12.
Lehr
,
F.
,
Millies
,
M.
, and
Mewes
,
D.
,
2002
, “
Bubble-Size Distributions and Flow Fields in Bubble Columns
,”
AlChE J.
,
48
(
11
), pp.
2426
2443
.
13.
Hulburt
,
H. M.
, and
Katz
,
S.
,
1964
, “
Some Problems in Particle Technology
,”
Chem. Eng. Sci.
,
19
(
8
), pp.
555
574
.
14.
Krepper
,
E.
,
Lucas
,
D.
,
Frank
,
T.
,
Prasser
,
H. M.
, and
Zwart
,
P. J.
,
2008
, “
The Inhomogeneous MUSIG Model for the Simulation of Polydispersed Flows
,”
Nucl. Eng. Des.
,
238
(
7
), pp.
1690
1702
.
15.
Lo
,
S.
,
1996
, AEA Technology.
16.
Maisels
,
A.
,
Kruis
,
F. E.
, and
Fissan
,
H.
,
2004
, “
Direct Simulation Monte Carlo for Simultaneous Nucleation, Coagulation, and Surface Growth in Dispersed Systems
,”
Chem. Eng. Sci.
,
59
(
11
), pp.
2231
2239
.
17.
McGraw
,
R.
,
1997
, “
Description of Aerosol Dynamics by the Quadrature Method of Moments
,”
Aerosol Sci. Technol.
,
27
(
2
), pp.
255
265
.
18.
Ramkrishna
,
D.
,
2000
,
Population Balances: Theory and Applications to Particulate Systems in Engineering
,
,
London
.
19.
Liao
,
Y. X.
, and
Lucas
,
D.
,
2009
, “
A Literature Review of Theoretical Models for Drop and Bubble Breakup in Turbulent Dispersions
,”
Chem. Eng. Sci.
,
64
(
15
), pp.
3389
3406
.
20.
Liao
,
Y. X.
, and
Lucas
,
D.
,
2010
, “
A Literature Review on Mechanisms and Models for the Coalescence Process of Fluid Particles
,”
Chem. Eng. Sci.
,
65
(
10
), pp.
2851
2864
.
21.
Liao
,
Y. X.
,
Lucas
,
D.
,
Krepper
,
E.
, and
Schmidtke
,
M.
,
2011
, “
Development of a Generalized Coalescence and Breakup Closure for the Inhomogeneous MUSIG Model
,”
Nucl. Eng. Des.
,
241
(
4
), pp.
1024
1033
.
22.
Liu
,
Z. Q.
,
Qi
,
F. S.
,
Li
,
B. K.
, and
Cheung
,
S. C. P.
,
2016
, “
Modeling of Bubble Behaviors and Size Distribution in a Slab Continuous Casting Mold
,”
Int. J. Multiphase Flow
,
79
, pp.
190
201
.
23.
Law
,
D.
, and
Battaglia
,
F.
,
2013
, “
Numerical Simulations for Hydrodynamics of Air-Water External Loop Airlift Reactor Flows With Bubble Break-Up and Coalescence Effects
,”
ASME J. Fluids Eng.
,
135
(
8
), p.
081302
.
24.
Randolph
,
A. D.
, and
Larson
,
M. A.
,
1988
,
Theory of Particulate Process: Analysis and Techniques of Continuous Crystallization
,
,
London
.
25.
Ishii
,
M.
, and
Hibiki
,
T.
,
2006
,
Thermo-Fluid Dynamics of Two-Phase Flow
,
,
New York
.
26.
Frenklach
,
M.
, and
Harris
,
S. J.
,
1987
, “
Aerosol Dynamics Modeling Using the Method of Moments
,”
J. Colloid Interface Sci.
,
118
(
1
), pp.
252
261
.
27.
Hibiki
,
T.
, and
Ishii
,
M.
,
2000
, “
One-Group Interfacial Area Transport of Bubbly Flows in Vertical Round Tubes
,”
Int. J. Heat Mass Transfer
,
43
(
15
), pp.
2711
2726
.
28.
Wu
,
Q.
,
Kim
,
S.
,
Ishii
,
M.
, and
Beus
,
S. G.
,
1998
, “
One-Group Interfacial Area Transport in Vertical Bubbly Flow
,”
Int. J. Heat Mass Transfer
,
41
(
8–9
), pp.
1103
1112
.
29.
Elghobashi
,
S.
,
1991
, “
Particle-Laden Turbulent Flows—Direct Simulation and Closure Models
,”
Appl. Sci. Res.
,
48
(
3–4
), pp.
301
314
.
30.
Elghobashi
,
S.
, and
Truesdell
,
G. C.
,
1992
, “
Direct Simulation of Particle Dispersion in a Decaying Isotropic Turbulence
,”
J. Fluid Mech.
,
242
(
1
), pp.
655
700
.
31.
Vance
,
M. W.
,
Squires
,
K. D.
, and
Simonin
,
O.
,
2006
, “
Properties of the Particle Velocity Field in Gas-Solid Turbulent Channel Flow
,”
Phys. Fluids
,
18
(
6
), p.
063302
.
32.
Yamamoto
,
Y.
,
Potthoff
,
M.
,
Tanaka
,
T.
,
Kajishima
,
T.
, and
Tsuji
,
Y.
,
2001
, “
Large-Eddy Simulation of Turbulent Gas-Particle Flow in a Vertical Channel: Effect of Considering Inter-Particle Collisions
,”
J. Fluid Mech.
,
442
, pp.
303
334
.
33.
Rani
,
S. L.
,
Winkler
,
C.
, and
Vanka
,
S.
,
2004
, “
Numerical Simulations of Turbulence Modulation by Dense Particles in a Fully Developed Pipe Flow
,”
Powder Technol.
,
141
(
1–2
), pp.
80
99
.
34.
Sun
,
J.
,
Battaglia
,
F.
, and
Subramaniam
,
S.
,
2007
, “
Hybrid Two-Fluid DEM Simulation of Gas-Solid Fluidized Beds
,”
ASME J. Fluids Eng.
,
129
(
11
), pp.
1394
1403
.
35.
Ma
,
J. S.
,
Hsiao
,
C. T.
, and
Chahine
,
G. L.
,
2015
, “
Euler-Lagrange Simulations of Bubble Cloud Dynamics Near a Wall
,”
ASME J. Fluids Eng.
,
137
(
4
), p.
041301
.
36.
Pirker
,
S.
,
Kahrimanovic
,
D.
, and
Schneiderbauer
,
S.
,
2015
, “
Secondary Vortex Formation in Bifurcated Submerged Entry Nozzles: Numerical Simulation of Gas Bubble Entrapment
,”
Metall. Trans. B
,
46
(
2
), pp.
953
960
.
37.
Yang
,
H.
, and
Thomas
,
B. G.
,
2017
, “
Multiphase Flow and Bubble Size Distribution in Continuous Casters Using a Hybrid EEDPM Model
,” SteelSIM2017, Qingdao, China, Aug. 16–18.
38.
Yang
,
H.
,
Vanka
,
S. P.
, and
Thomas
,
B. G.
,
2018
, “
Modeling of Argon Gas Behavior in Continuous Casting of Steel
,”
TMS Annual Meeting and Exhibition
, Phoenix, AZ, Mar. 11–15, pp.
119
131
.
39.
Hibiki
,
T.
,
Ishii
,
M.
, and
Xiao
,
Z.
,
2001
, “
Axial Interfacial Area Transport of Vertical Bubbly Flows
,”
Int. J. Heat Mass Transfer
,
44
(
10
), pp.
1869
1888
.
40.
Rzehak
,
R.
, and
Kriebitzsch
,
S.
,
2015
, “
Multiphase CFD-Simulation of Bubbly Pipe Flow: A Code Comparison
,”
Int. J. Multiphase Flow
,
68
, pp.
135
152
.
41.
Frank
,
T.
, and
Menter
,
F. R.
,
2003
, “
Bubble Flow in Vertical Pipes
,”
ANSYS
,
Inc., Canonsburg, PA
.
42.
Tomiyama
,
A.
,
Kataoka
,
I.
,
Zun
,
I.
, and
Sakaguchi
,
T.
,
1998
, “
Drag Coefficients of Single Bubbles Under Normal and Micro Gravity Conditions
,”
JSME Int. J. Ser. B
,
41
(
2
), pp.
472
479
.
43.
Tomiyama
,
A.
,
Tamai
,
H.
,
Zun
,
I.
, and
Hosokawa
,
S.
,
2002
, “
Transverse Migration of Single Bubbles in Simple Shear Flows
,”
Chem. Eng. Sci.
,
57
(
11
), pp.
1849
1858
.
44.
Hosokawa
,
S.
,
Tomiyama
,
A.
,
Misaki
,
S.
, and
,
T.
, 2002, “
Lateral Migration of Single Bubbles Due to the Presence of Wall
,”
ASME
Paper No. FEDSM2002-31148.
45.
Burns
,
A. D.
,
Frank
,
T.
,
Hamill
,
I.
, and
Shi
,
J.-M.
,
2004
, “
The Favre Averaged Drag Model for Turbulent Dispersion in Eulerian Multi-Phase Flows
,”
Fifth International Conference on Multiphase Flow (ICMF)
, pp.
1
17
.
46.
ANSYS
,
2012
, “
ANSYS Fluent Theory Manual
,” ANSYS, Inc., Canonsburg, PA.
47.
Auton
,
T.
,
1987
, “
The Lift Force on a Spherical Body in a Rotational Flow
,”
J. Fluid Mech.
,
183
(
1
), pp.
199
218
.
48.
Prince
,
M. J.
, and
Blanch
,
H. W.
,
1990
, “
Bubble Coalescence and Break-Up in Air-Sparged Bubble-Columns
,”
AlChE J.
,
36
(
10
), pp.
1485
1499
.
49.
Luo
,
H.
, and
Svendsen
,
H. F.
,
1996
, “
Theoretical Model for Drop and Bubble Breakup in Turbulent Dispersions
,”
AlChE J.
,
42
(
5
), pp.
1225
1233
.
50.
Wang
,
T. F.
,
Wang
,
J. F.
, and
Jin
,
Y.
,
2003
, “
A Novel Theoretical Breakup Kernel Function for Bubbles/Droplets in a Turbulent Flow
,”
Chem. Eng. Sci.
,
58
(
20
), pp.
4629
4637
.
51.
Nambiar
,
D. K. R.
,
Kumar
,
R.
,
Das
,
T. R.
, and
Gandhi
,
K. S.
,
1992
, “
A New Model for the Breakage Frequency of Drops in Turbulent Stirred Dispersions
,”
Chem. Eng. Sci.
,
47
(
12
), pp.
2989
3002
.
52.
Kuboi
,
R.
,
Komasawa
,
I.
, and
Otake
,
T.
,
1972
, “
Collision and Coalescence of Dispersed Drops in Turbulent Liquid Flow
,”
J. Chem. Eng. Jpn.
,
5
(
4
), pp.
423
424
.
53.
Kuboi
,
R.
,
Komasawa
,
I.
, and
Otake
,
T.
,
1972
, “
Behavior of Dispersed Particles in Turbulent Liquid Flow
,”
J. Chem. Eng. Jpn.
,
5
(
4
), pp.
349
355
.
54.
Hesketh
,
R. P.
,
Etchells
,
A. W.
, and
Russell
,
T. W. F.
,
1991
, “
Experimental-Observations of Bubble Breakage in Turbulent-Flow
,”
Ind. Eng. Chem. Res.
,
30
(
5
), pp.
835
841
.