# CFD Modeling of Gas-Liquid Bubbly Flow in Horizontal Pipes: Influence of Bubble Coalescence and Breakup

Authors: K. Ekambara; R. Sean Sanders; K. Nandakumar; J. H. Masliyah
Journal: International Journal of Chemical Engineering
 (2012)

---

## Abstract

Modelling of gas-liquid bubbly flows is achieved by coupling a population balance equation with the three-dimensional, two-fluid, hydrodynamic model. For gas-liquid bubbly flows, an average bubble number density transport equation has been incorporated in the CFD code CFX 5.7 to describe the temporal and spatial evolution of the gas bubbles population. The coalescence and breakage effects of the gas bubbles are modeled. The coalescence by the random collision driven by turbulence and wake entrainment is considered, while for bubble breakage, the impact of turbulent eddies is considered. Local spatial variations of the gas volume fraction, interfacial area concentration, Sauter mean bubble diameter, and liquid velocity are compared against experimental data in a horizontal pipe, covering a range of gas (0.25 to 1.34 m/s) and liquid (3.74 to 5.1 m/s) superficial velocities and average volume fractions (4% to 21%). The predicted local variations are in good agreement with the experimental measurements reported in the literature. Furthermore, the development of the flow pattern was examined at three different axial locations ofL/D = 25, 148, and 253. The first location is close to the entrance region where the flow is still developing, while the second and the third represent nearly fully developed bubbly flow patterns.


---

## Body

## 1. Introduction

Gas-liquid, two-phase flow in horizontal pipes is encountered often in a number of industrial processes. Common applications include chemical plants, evaporators, oil wells and pipelines, fluidized bed combustors, and evaporators. Horizontal bubbly flows have received less attention in the literature than vertical flows, even though this flow orientation is equally important in industrial applications such as hydrotransport, an important technology in bitumen extraction. Experimental observations are also difficult in this case, as the migration of dispersed bubbles towards the top of the pipe, due to buoyancy, causes a highly nonsymmetric volume distribution in the pipe cross-section. This density stratification is not often accompanied by a strong secondary flow. Gas volume fraction, interfacial area concentration, and mean bubble diameter are the three characterizing field variables that characterize the internal flow structure of two-phase, gas-liquid flows in horizontal pipe. In various industrial processes, the gas volume fraction parameter is required for hydrodynamic and thermal design. The interfacial transport of mass, momentum, and energy is proportional to the interfacial area and the driving forces. This is an important parameter required for a two-fluid model formulation. The mean bubble diameter serves as a link between the gas volume fraction and interfacial area concentration. An accurate knowledge of local distributions of these three parameters is of great importance to the eventual understanding and modelling of the interfacial transfer processes [1].In the past two decades, significant developments in the modeling of two-phase flow processes have occurred since the introduction of the two-fluid model. In the volume averaged, two-fluid model, the interfacial transfer terms are strongly related to the interfacial area concentration and the local transfer mechanisms such as the degree of turbulence near the interfaces. Fundamentally, the interfacial transport of mass, momentum, and energy are proportional to the interfacial area concentration (aij) and driving forces. Since the interfacial area concentration aij represents the key parameter that links the interaction of the phases, significant attention has been paid towards developing a better understanding of the coalescence and breakage effects due to interactions among bubbles and between bubbles and turbulent eddies for gas-liquid bubbly flows [2–5]. The population balance method is a well-known method for tracking the size distribution of the dispersed phase and accounting for the breakage and coalescence effects in bubbly flows (see, e.g. [6–14]). This approach is concerned with maintaining a record of the number of bubbles initially and tracking their evolution in space over time.In this work, an attempt has been made to demonstrate the possibility of combining population balance with computational fluid dynamics (CFD) for the case of gas-liquid bubbly flow in the horizontal pipe. The MUSIG model has been implemented in CFX-5.7 to account for the nonuniform bubble size distribution in a gas-liquid flow [7, 15, 16]. The gas volume fraction, interfacial area concentration, Suater mean diameter, and axial liquid velocity have been predicted for a wide range of gas and liquid flow condition. Further, the development of flow pattern has been studied at three different axial locations. The model predictions are compared with available experimental data from the literature.


## 2. Mathematical Modelling

### 2.1. Population Balance Model

Population balance modelling is used in computing the size distribution of the dispersed phase and in accounting for the breakage and coalescence effects in bubbly flows. A general form of the population balance equation is(1)∂ni∂t+∇⋅(ugni)=BB+BC-DB-DC,
where ug is the gas velocity, ni represents the number density of size group i, and terms on the right hand side BB, BC, DB, and DC are, respectively, the “birth” and “death” due to breakup and coalescence of bubbles. The left hand side tracks the spatial and temporal evolution of a class of bubbles, while the right hand side models the exchange between classes due to breakup and coalescence of bubbles. The bubble number density ni is related to the gas volume fraction αg by(2)αgfi=niVi,
where fi represents the volume fraction of bubbles of group i, and Vi is the corresponding volume of a bubble of group i. It is necessary to provide individual models for each of the breakup and coalescence processes as it depends on the mechanisms and is sensitively dependent on the presence of surfactants, turbulence levels, and so forth. These models are discussed next.


#### 2.1.1. Bubble Breakup Model

The breakup of bubbles in turbulent dispersions employs the model developed by Luo and Svendsen [18]. Binary break-up of the bubbles is assumed, and the model is based on the theories of isotropic turbulence. For binary breakage, a dimensionless variable describing the sizes of daughter drops or bubbles (the breakage volume fraction) can be defined as(3)fBV=ViV=di3d3=di3di3+dj3,
where di and dj are diameters (corresponding to Vi and Vj) of the daughter bubbles in the binary breakage of a parent bubble with diameter d (corresponding to volume V). The value interval of the breakage volume fraction is between 0 and 1. The break-up rate of bubbles of volume Vj into volume sizes of Vi (=VfBV) can be obtained as(4)Ω(Vj:Vi)(1-αg)nj=C(∈dj2)1/3∫ζmin⁡1(1+ζ)2ζ11/3exp⁡(-12cfσβρl∈2/3dj5/3ζ11/3)dζ,
where ∈ is the rate of energy dissipation per unit of liquid mass; ζ=λ/dj is the size ratio between an eddy and a particle in the inertial subrange and consequently ζmin=λmin/dj; C and β are determined, respectively, from fundamental consideration of drops or bubbles breakage in turbulent dispersion systems to be 0.923 and 2.0 in Luo and Svendsen [18]; cf is the increase coefficient of surface area given by(5)cf=[fBV2/3+(1-fBV)2/3-1].
The birth rate of group i bubbles due to break-up of larger bubbles is(6)BB=∑j=i+1NΩ(Vj:Vi)nj.
The death rate of group i bubbles due to break-up into smaller bubbles is(7)DB=ΩiniwithΩi=∑k=1iΩki.


#### 2.1.2. Bubble Coalescence Model

The coalescence of two bubbles is assumed to occur in three steps. The first step where the bubbles collide and trap a layer of liquid between them, a second step where this liquid layer drains until it reaches a critical thickness, and a last step during which this liquid film disappears and the bubbles coalesce. The collisions between bubbles may be caused by turbulence, buoyancy, or laminar shear. Only the first cause of collision (turbulence) is considered in the present model. Indeed collisions caused by buoyancy cannot be taken into account as all the bubbles from each class move at the same speed. The coalescence rate considering turbulent collision taken from Prince and Blanch [19] can be expressed as(8)χ=θijexp⁡⁡(-tijτij),
where tij  is the contact time for two bubbles given by (dij/2)2/3/∈1/3.When bubbles collide, a small amount of liquid is entrapped between them, forming a small circular lens or film of radiusR and thickness h. The forces causing the film or lens to grow thinner in pure systems arise from capillary pressure, augmented by compression from a close range Hamaker force which accounts for the mutual attraction of water molecules on opposite sides of the liquid film [19]. For equal size bubbles, Oolman and Blanch [20] derived the thinning formula(9)-dhdt={8R2ρL[h2(2σrb+A6πh3)]}1/2.
Prince and Blanch [19] solved the above equation numerically and show that tij, the time required for two bubbles, having diameters di and dj to coalesce is estimated to be {(dij/2)3ρl/16σ}1/2ln⁡(h0/hj). The equivalent diameter dij is calculated as suggested by Chesters and Hoffman [21]: dij=(2/di+2/dj)-1. The parameters h0 and hj represent the film thickness when collision begins and critical film thickness at which rupture occurs, respectively. The values of the above parameters depend mainly on the physical properties of the liquid phase and have been experimentally computed for the air-water system. According to Prince and Blanch [19], for air-water systems, experiments have determined h0 and hj to be 1 × 10−4 m [22] and 1 × 10−8 m [23], respectively.The turbulent collision rateθij for two bubbles of diameters di and dj is given by(10)θij=π4[di+dj]2(uti2+utj2)1/2,
where the turbulent velocity ut in the inertial subrange of isotropic turbulence [24] is,(11)ut=1.4∈1/3d1/3.
The birth rate of group i due to coalescence of group k and group l bubbles is:(12)BC=12∑k=1N∑l=1Nχi,klninj.
The death rate of group i due to coalescence with other bubbles is:(13)DC=∑j=1Nχijninj.


### 2.2. Flow Equations

The numerical simulations presented are based on the two-fluid, Eulerian-Eulerian model. The Eulerian modelling framework is based on ensemble-averaged mass and momentum transport equations for each phase. Regarding the liquid phase (αl) as the continuum and the gaseous phase (bubbles) as the dispersed phase (αg), these equations without interface mass transfer can be written in standard form as follows.Continuity equation of the liquid phase(14)∂∂t(ρlαl)+∇⋅(ρlαlul)=0.Continuity equation of the gas phase(15)∂∂t(ρgαgfi)+∇⋅(ρgαgugfi)=Si.Momentum equation(16)∂∂t(ρkαkuk)+∇⋅(ρkαkukuk)=-αk∇p+ρkαkgi+αkμk∇2u+Fkm(k,m=l,g).In (15), Si is a source term that captures the coalescence and break-up processes. The right side of  (16) describes the following forces acting on the phase k: the pressure gradient, gravity, and the viscous stress term, and Fkm represents the sum of the interfacial forces that include the drag force FD, the lift force FL, the virtual mass force FVM, the wall lubrication force FWL, and the turbulent dispersion force FTD. Detailed descriptions of each of these forces can be found in Anglart and Nylund [25]; Lahey and Drew [26], and Joshi [27].The origin of the drag force is due to the resistance experienced by a body moving in the liquid. Viscous stress creates skin drag, and pressure distribution around the moving body creates form drag. The formulation of the drag force is a key issue in multiphase flows. Clift et al. [28] and Joshi et al. [29] have given excellent accounts of this subject. The interphase momentum transfer between gas and liquid due to drag force is given by(17)FD=34CDαgρl1dS|ul-ug|(ul-ug).The lift force considers the interaction of the bubble with the shear field of the liquid. It acts perpendicular to the main flow direction and is proportional to the gradient of the liquid velocity field. The lift force in terms of the slip velocity and the curl of the liquid phase velocity can be modelled as [30–33](18)FL=CLαgρl(ug-ul)×∇×ul.The wall lubrication force arises because the liquid flow rate between bubble and the wall is lower than between the bubble and the main flow. This results in a hydrodynamics pressure difference driving bubble away from the wall. This force density is approximated as [34](19)FWL=-αgρl(ur-(ur⋅nw)nw)dSmax⁡⁡[C1+C2dSyw,0].
Here, ur=ul-ug is the relative velocity between phases, dS is the dispersed phase Sauter mean bubble diameter, yw is the distance to the nearest wall, and nw is the unit normal pointing away from the wall.The turbulent dispersion force, derived by Lopez de Bertodano [35], is based on an analogy with molecular movement. The turbulence-induced dispersion is a function of turbulent kinetic energy and gradient of the volume fraction of the liquid:(20)FTD=-CTDρlk∇αl.
The drag coefficient CD in (17) has been modelled using Ishii and Zuber [36] drag model. The lift coefficient CL is theoretically proven to be 0.5 for a spherical bubble in a potential flow [37]. It is also known that (i) CL becomes smaller than 0.5 for a single small bubble in a viscous flow, as show by Lopez de Bertodano et al. [38] and Lance and Lopez de Bertodano, [39] (CL=0.25; 0.1, resp.), and (ii) CL strongly depends on the bubble diameter and decreases with dS [40]. These facts indicate that CL is a function of bubble diameter and fluid properties. From this point of view, it is necessary to consider the lift coefficient also to depend on flow conditions. Recently, Tomiyama et al. [41] have developed an empirical correlation for the lift coefficient as a function of Reynolds number and Eotvos number. We have found that this correlation does not perform well for horizontal flows because of the migration of dispersed bubbles towards the top of the pipe. In view of this, we have developed a correlation for lift coefficient in horizontal flows. An interesting finding and a main contribution in this work is that a wide range of flow behavior of two-phase bubbly flows in horizontal pipes is represented by a unique functional relationship between the lift coefficient and the flow Reynolds number. When such closure relations are tested over a wide range of two-phase flows (not only pipe flows, but also bubble columns, etc.) our confidence in using such models to study practical two-phase problems in process equipment will increase over a period of time. Further explanation about CL is given in the results and discussion. The wall lubrication constants C1 and C2, as suggested by Antal et al., [34], are −0.01 and 0.05, respectively. The coefficient CTD = 0.5 was found to give the good results which is in the recommended range of 0.1 to 1.0 [35]. By definition, the interfacial area concentration aij for bubbly flows can be determined through the relationship(21)aij=6αgdS,
where dS is the bubble Sauter mean diameter. The local bubble Sauter mean diameter based on the calculated values of the scalar fraction fi and discrete bubble sizes di can be deduced from(22)dS=1∑ifi/di.
From the drag and nondrag forces above, it is evident that the interfacial area concentration aij and the bubble Sauter mean bubble diameter in (22) are essential parameters that link the interaction between the liquid and gas (bubbly) phases. In most two-phase flow studies, the common approach of prescribing constant bubble sizes through the mean bubble Sauter diameter is still prevalent. Such an approach does not allow dynamic representation of the changes in the interfacial structure.


### 2.3. Turbulence Equations

For the continuous liquid phase, ak-∈ model is applied with its standard constants: C∈1=1.44, C∈2=1.92, Cμ=0.09, σk=1.0, and σ∈=1.3. No turbulence model is applied on the dispersed gas phase, but the influence of the dispersed phase on the turbulence of the continuous phase is taken into account with Sato’s additional term [42]. The governing equations for the turbulent kinetic energy k and turbulent dissipation ε are
(23)∂∂t(ρlαlk)+∂∂xi(ρlαlulk)=∂∂xi(αl(μl+μl,turσk)∂k∂xi)+αl(G-αlρl∈),∂∂t(ρlαl∈l)+∂∂xi(ρlαlul∈)=∂∂xi(αl(μl+μl,turσ∈)∂∈∂xi)+αl∈k(Cε1G-Cε2αlρl∈),
where G is the turbulence production due to viscous and buoyancy forces, which is modeled using(24)G=μt∇u⋅(∇u+∇uT)-23∇⋅u(3μt∇⋅u+ρlk)-μtρσpρlβg⋅∇T.


## 2.1. Population Balance Model

Population balance modelling is used in computing the size distribution of the dispersed phase and in accounting for the breakage and coalescence effects in bubbly flows. A general form of the population balance equation is(1)∂ni∂t+∇⋅(ugni)=BB+BC-DB-DC,
where ug is the gas velocity, ni represents the number density of size group i, and terms on the right hand side BB, BC, DB, and DC are, respectively, the “birth” and “death” due to breakup and coalescence of bubbles. The left hand side tracks the spatial and temporal evolution of a class of bubbles, while the right hand side models the exchange between classes due to breakup and coalescence of bubbles. The bubble number density ni is related to the gas volume fraction αg by(2)αgfi=niVi,
where fi represents the volume fraction of bubbles of group i, and Vi is the corresponding volume of a bubble of group i. It is necessary to provide individual models for each of the breakup and coalescence processes as it depends on the mechanisms and is sensitively dependent on the presence of surfactants, turbulence levels, and so forth. These models are discussed next.


### 2.1.1. Bubble Breakup Model

The breakup of bubbles in turbulent dispersions employs the model developed by Luo and Svendsen [18]. Binary break-up of the bubbles is assumed, and the model is based on the theories of isotropic turbulence. For binary breakage, a dimensionless variable describing the sizes of daughter drops or bubbles (the breakage volume fraction) can be defined as(3)fBV=ViV=di3d3=di3di3+dj3,
where di and dj are diameters (corresponding to Vi and Vj) of the daughter bubbles in the binary breakage of a parent bubble with diameter d (corresponding to volume V). The value interval of the breakage volume fraction is between 0 and 1. The break-up rate of bubbles of volume Vj into volume sizes of Vi (=VfBV) can be obtained as(4)Ω(Vj:Vi)(1-αg)nj=C(∈dj2)1/3∫ζmin⁡1(1+ζ)2ζ11/3exp⁡(-12cfσβρl∈2/3dj5/3ζ11/3)dζ,
where ∈ is the rate of energy dissipation per unit of liquid mass; ζ=λ/dj is the size ratio between an eddy and a particle in the inertial subrange and consequently ζmin=λmin/dj; C and β are determined, respectively, from fundamental consideration of drops or bubbles breakage in turbulent dispersion systems to be 0.923 and 2.0 in Luo and Svendsen [18]; cf is the increase coefficient of surface area given by(5)cf=[fBV2/3+(1-fBV)2/3-1].
The birth rate of group i bubbles due to break-up of larger bubbles is(6)BB=∑j=i+1NΩ(Vj:Vi)nj.
The death rate of group i bubbles due to break-up into smaller bubbles is(7)DB=ΩiniwithΩi=∑k=1iΩki.


### 2.1.2. Bubble Coalescence Model

The coalescence of two bubbles is assumed to occur in three steps. The first step where the bubbles collide and trap a layer of liquid between them, a second step where this liquid layer drains until it reaches a critical thickness, and a last step during which this liquid film disappears and the bubbles coalesce. The collisions between bubbles may be caused by turbulence, buoyancy, or laminar shear. Only the first cause of collision (turbulence) is considered in the present model. Indeed collisions caused by buoyancy cannot be taken into account as all the bubbles from each class move at the same speed. The coalescence rate considering turbulent collision taken from Prince and Blanch [19] can be expressed as(8)χ=θijexp⁡⁡(-tijτij),
where tij  is the contact time for two bubbles given by (dij/2)2/3/∈1/3.When bubbles collide, a small amount of liquid is entrapped between them, forming a small circular lens or film of radiusR and thickness h. The forces causing the film or lens to grow thinner in pure systems arise from capillary pressure, augmented by compression from a close range Hamaker force which accounts for the mutual attraction of water molecules on opposite sides of the liquid film [19]. For equal size bubbles, Oolman and Blanch [20] derived the thinning formula(9)-dhdt={8R2ρL[h2(2σrb+A6πh3)]}1/2.
Prince and Blanch [19] solved the above equation numerically and show that tij, the time required for two bubbles, having diameters di and dj to coalesce is estimated to be {(dij/2)3ρl/16σ}1/2ln⁡(h0/hj). The equivalent diameter dij is calculated as suggested by Chesters and Hoffman [21]: dij=(2/di+2/dj)-1. The parameters h0 and hj represent the film thickness when collision begins and critical film thickness at which rupture occurs, respectively. The values of the above parameters depend mainly on the physical properties of the liquid phase and have been experimentally computed for the air-water system. According to Prince and Blanch [19], for air-water systems, experiments have determined h0 and hj to be 1 × 10−4 m [22] and 1 × 10−8 m [23], respectively.The turbulent collision rateθij for two bubbles of diameters di and dj is given by(10)θij=π4[di+dj]2(uti2+utj2)1/2,
where the turbulent velocity ut in the inertial subrange of isotropic turbulence [24] is,(11)ut=1.4∈1/3d1/3.
The birth rate of group i due to coalescence of group k and group l bubbles is:(12)BC=12∑k=1N∑l=1Nχi,klninj.
The death rate of group i due to coalescence with other bubbles is:(13)DC=∑j=1Nχijninj.


## 2.1.1. Bubble Breakup Model

The breakup of bubbles in turbulent dispersions employs the model developed by Luo and Svendsen [18]. Binary break-up of the bubbles is assumed, and the model is based on the theories of isotropic turbulence. For binary breakage, a dimensionless variable describing the sizes of daughter drops or bubbles (the breakage volume fraction) can be defined as(3)fBV=ViV=di3d3=di3di3+dj3,
where di and dj are diameters (corresponding to Vi and Vj) of the daughter bubbles in the binary breakage of a parent bubble with diameter d (corresponding to volume V). The value interval of the breakage volume fraction is between 0 and 1. The break-up rate of bubbles of volume Vj into volume sizes of Vi (=VfBV) can be obtained as(4)Ω(Vj:Vi)(1-αg)nj=C(∈dj2)1/3∫ζmin⁡1(1+ζ)2ζ11/3exp⁡(-12cfσβρl∈2/3dj5/3ζ11/3)dζ,
where ∈ is the rate of energy dissipation per unit of liquid mass; ζ=λ/dj is the size ratio between an eddy and a particle in the inertial subrange and consequently ζmin=λmin/dj; C and β are determined, respectively, from fundamental consideration of drops or bubbles breakage in turbulent dispersion systems to be 0.923 and 2.0 in Luo and Svendsen [18]; cf is the increase coefficient of surface area given by(5)cf=[fBV2/3+(1-fBV)2/3-1].
The birth rate of group i bubbles due to break-up of larger bubbles is(6)BB=∑j=i+1NΩ(Vj:Vi)nj.
The death rate of group i bubbles due to break-up into smaller bubbles is(7)DB=ΩiniwithΩi=∑k=1iΩki.


## 2.1.2. Bubble Coalescence Model

The coalescence of two bubbles is assumed to occur in three steps. The first step where the bubbles collide and trap a layer of liquid between them, a second step where this liquid layer drains until it reaches a critical thickness, and a last step during which this liquid film disappears and the bubbles coalesce. The collisions between bubbles may be caused by turbulence, buoyancy, or laminar shear. Only the first cause of collision (turbulence) is considered in the present model. Indeed collisions caused by buoyancy cannot be taken into account as all the bubbles from each class move at the same speed. The coalescence rate considering turbulent collision taken from Prince and Blanch [19] can be expressed as(8)χ=θijexp⁡⁡(-tijτij),
where tij  is the contact time for two bubbles given by (dij/2)2/3/∈1/3.When bubbles collide, a small amount of liquid is entrapped between them, forming a small circular lens or film of radiusR and thickness h. The forces causing the film or lens to grow thinner in pure systems arise from capillary pressure, augmented by compression from a close range Hamaker force which accounts for the mutual attraction of water molecules on opposite sides of the liquid film [19]. For equal size bubbles, Oolman and Blanch [20] derived the thinning formula(9)-dhdt={8R2ρL[h2(2σrb+A6πh3)]}1/2.
Prince and Blanch [19] solved the above equation numerically and show that tij, the time required for two bubbles, having diameters di and dj to coalesce is estimated to be {(dij/2)3ρl/16σ}1/2ln⁡(h0/hj). The equivalent diameter dij is calculated as suggested by Chesters and Hoffman [21]: dij=(2/di+2/dj)-1. The parameters h0 and hj represent the film thickness when collision begins and critical film thickness at which rupture occurs, respectively. The values of the above parameters depend mainly on the physical properties of the liquid phase and have been experimentally computed for the air-water system. According to Prince and Blanch [19], for air-water systems, experiments have determined h0 and hj to be 1 × 10−4 m [22] and 1 × 10−8 m [23], respectively.The turbulent collision rateθij for two bubbles of diameters di and dj is given by(10)θij=π4[di+dj]2(uti2+utj2)1/2,
where the turbulent velocity ut in the inertial subrange of isotropic turbulence [24] is,(11)ut=1.4∈1/3d1/3.
The birth rate of group i due to coalescence of group k and group l bubbles is:(12)BC=12∑k=1N∑l=1Nχi,klninj.
The death rate of group i due to coalescence with other bubbles is:(13)DC=∑j=1Nχijninj.


## 2.2. Flow Equations

The numerical simulations presented are based on the two-fluid, Eulerian-Eulerian model. The Eulerian modelling framework is based on ensemble-averaged mass and momentum transport equations for each phase. Regarding the liquid phase (αl) as the continuum and the gaseous phase (bubbles) as the dispersed phase (αg), these equations without interface mass transfer can be written in standard form as follows.Continuity equation of the liquid phase(14)∂∂t(ρlαl)+∇⋅(ρlαlul)=0.Continuity equation of the gas phase(15)∂∂t(ρgαgfi)+∇⋅(ρgαgugfi)=Si.Momentum equation(16)∂∂t(ρkαkuk)+∇⋅(ρkαkukuk)=-αk∇p+ρkαkgi+αkμk∇2u+Fkm(k,m=l,g).In (15), Si is a source term that captures the coalescence and break-up processes. The right side of  (16) describes the following forces acting on the phase k: the pressure gradient, gravity, and the viscous stress term, and Fkm represents the sum of the interfacial forces that include the drag force FD, the lift force FL, the virtual mass force FVM, the wall lubrication force FWL, and the turbulent dispersion force FTD. Detailed descriptions of each of these forces can be found in Anglart and Nylund [25]; Lahey and Drew [26], and Joshi [27].The origin of the drag force is due to the resistance experienced by a body moving in the liquid. Viscous stress creates skin drag, and pressure distribution around the moving body creates form drag. The formulation of the drag force is a key issue in multiphase flows. Clift et al. [28] and Joshi et al. [29] have given excellent accounts of this subject. The interphase momentum transfer between gas and liquid due to drag force is given by(17)FD=34CDαgρl1dS|ul-ug|(ul-ug).The lift force considers the interaction of the bubble with the shear field of the liquid. It acts perpendicular to the main flow direction and is proportional to the gradient of the liquid velocity field. The lift force in terms of the slip velocity and the curl of the liquid phase velocity can be modelled as [30–33](18)FL=CLαgρl(ug-ul)×∇×ul.The wall lubrication force arises because the liquid flow rate between bubble and the wall is lower than between the bubble and the main flow. This results in a hydrodynamics pressure difference driving bubble away from the wall. This force density is approximated as [34](19)FWL=-αgρl(ur-(ur⋅nw)nw)dSmax⁡⁡[C1+C2dSyw,0].
Here, ur=ul-ug is the relative velocity between phases, dS is the dispersed phase Sauter mean bubble diameter, yw is the distance to the nearest wall, and nw is the unit normal pointing away from the wall.The turbulent dispersion force, derived by Lopez de Bertodano [35], is based on an analogy with molecular movement. The turbulence-induced dispersion is a function of turbulent kinetic energy and gradient of the volume fraction of the liquid:(20)FTD=-CTDρlk∇αl.
The drag coefficient CD in (17) has been modelled using Ishii and Zuber [36] drag model. The lift coefficient CL is theoretically proven to be 0.5 for a spherical bubble in a potential flow [37]. It is also known that (i) CL becomes smaller than 0.5 for a single small bubble in a viscous flow, as show by Lopez de Bertodano et al. [38] and Lance and Lopez de Bertodano, [39] (CL=0.25; 0.1, resp.), and (ii) CL strongly depends on the bubble diameter and decreases with dS [40]. These facts indicate that CL is a function of bubble diameter and fluid properties. From this point of view, it is necessary to consider the lift coefficient also to depend on flow conditions. Recently, Tomiyama et al. [41] have developed an empirical correlation for the lift coefficient as a function of Reynolds number and Eotvos number. We have found that this correlation does not perform well for horizontal flows because of the migration of dispersed bubbles towards the top of the pipe. In view of this, we have developed a correlation for lift coefficient in horizontal flows. An interesting finding and a main contribution in this work is that a wide range of flow behavior of two-phase bubbly flows in horizontal pipes is represented by a unique functional relationship between the lift coefficient and the flow Reynolds number. When such closure relations are tested over a wide range of two-phase flows (not only pipe flows, but also bubble columns, etc.) our confidence in using such models to study practical two-phase problems in process equipment will increase over a period of time. Further explanation about CL is given in the results and discussion. The wall lubrication constants C1 and C2, as suggested by Antal et al., [34], are −0.01 and 0.05, respectively. The coefficient CTD = 0.5 was found to give the good results which is in the recommended range of 0.1 to 1.0 [35]. By definition, the interfacial area concentration aij for bubbly flows can be determined through the relationship(21)aij=6αgdS,
where dS is the bubble Sauter mean diameter. The local bubble Sauter mean diameter based on the calculated values of the scalar fraction fi and discrete bubble sizes di can be deduced from(22)dS=1∑ifi/di.
From the drag and nondrag forces above, it is evident that the interfacial area concentration aij and the bubble Sauter mean bubble diameter in (22) are essential parameters that link the interaction between the liquid and gas (bubbly) phases. In most two-phase flow studies, the common approach of prescribing constant bubble sizes through the mean bubble Sauter diameter is still prevalent. Such an approach does not allow dynamic representation of the changes in the interfacial structure.


## 2.3. Turbulence Equations

For the continuous liquid phase, ak-∈ model is applied with its standard constants: C∈1=1.44, C∈2=1.92, Cμ=0.09, σk=1.0, and σ∈=1.3. No turbulence model is applied on the dispersed gas phase, but the influence of the dispersed phase on the turbulence of the continuous phase is taken into account with Sato’s additional term [42]. The governing equations for the turbulent kinetic energy k and turbulent dissipation ε are
(23)∂∂t(ρlαlk)+∂∂xi(ρlαlulk)=∂∂xi(αl(μl+μl,turσk)∂k∂xi)+αl(G-αlρl∈),∂∂t(ρlαl∈l)+∂∂xi(ρlαlul∈)=∂∂xi(αl(μl+μl,turσ∈)∂∈∂xi)+αl∈k(Cε1G-Cε2αlρl∈),
where G is the turbulence production due to viscous and buoyancy forces, which is modeled using(24)G=μt∇u⋅(∇u+∇uT)-23∇⋅u(3μt∇⋅u+ρlk)-μtρσpρlβg⋅∇T.


## 3. Method of Solution

The multiple size group (MUSIG) model (CFX 5.7 from ANSYS) used in this study combines the population balance method with the break-up [18] and coalescence [19] models in order to predict the bubble size distribution of the dispersed phase, and it uses the Eulerian-Eulerian two-fluid model. A standard two-phase flow calculation, with equation for continuity, momentum, and turbulence for a continuous and a dispersed phase, can be extended to include mass fraction of bubbles within several size ranges using the MUSIG model. The size range of the bubbles is split into several groups with, for example, bands of equal diameter or equal volume. Equations are solved for the mass fraction in each band. The MUSIG model has been implemented in the CFX-5.7 software to account nonuniform bubble size distribution in a gas-liquid mixture. The MUSIG model has been extensively used for different systems [7, 10, 15, 16, 43, 44]. These size fractions provide a more accurate measure of the interfacial area density and therefore allow better calculation of the heat and mass transfer taking place between the continuous and dispersed phases.In this present study, bubbles ranging from 1 mm to 10 mm diameter are equally divided into 10 classes (see Table1) as the experimental observation of maximum bubble diameter for highest superficial gas velocity is 6 mm. Even if we have considered the range, model predictions picks up the experimental observation bubble size range. The fate of the discrete bubble sizes so prescribed was tracked using the population balance model. Instead of considering 11 different complete phases, it was assumed that each bubble class travels at the same mean algebraic velocity to reduce the computational time and cost. Therefore, it results in 10 continuity equations for the gas phase coupled with a single continuity equation for the liquid phase. Sensitivity of the number of size groups needed to describe a meaningful distribution was examined by dividing the bubble diameters equally into 10, 15, and 20 size groups. The results revealed that no appreciable difference (± 2 %) was found for the predicted maximum bubble Sauter mean diameter using the 10, 15, or 20 bubble size groups. For the subdivision into 10 size groups, the maximum Sauter bubble diameter was under predicted by a maximum difference of 2%. In view of computational resources and times, it was therefore concluded that the subdivision of the bubbles sizes into 10 size groups was sufficient and all subsequent computational results are based on the discretization of the bubble sizes into 10 groups.Table 1
Diameter of each bubble class tracked in the simulation.
Class index12345678910Bubble diameter,di (mm)1.452.353.254.155.055.956.857.758.659.55Solution to the coupled sets of governing equations for the balances of mass and momentum of each phase was obtained using CFX 5.7. The conservation equations were discretized using the control volume technique. Computational grid is based on the unstructured set of blocks each containing structured grid. The structured grid within each block is generated using general curvilinear coordinates ensuring accurate representation of the flow boundaries. In order to select an adequate grid resolution, the effect of changing grid size was investigated. Several simulations were carried out using progressively larger number of grid points of 87156, 152482, 257670, 303245, and 341612. Sample grid sensitivity results are shown in Figure1. It can be seen that there is practically no change in the gas volume fraction and liquid velocity profiles when the grid size increased beyond 257670. In view of the observed effect of grid size, the simulations have been carried out by using 257670 grid points. Initial simulations were carried out with a coarse mesh to obtain an initially converged solution and to obtain an indication of where a high mesh density was needed. However, a dense mesh required additional computational effort. The velocity-pressure linkage was handled through the SIMPLE procedure. Three-dimensional transient simulations were performed. The time stepping strategy used in the transient simulations to reach a steady state was typically a variable step size strategy according to the following scheme: 100 steps at 1 × 10−4 s, followed by 300 steps at 5 × 10−4 s, 400 steps at 1 × 10−3 s, 1400 steps at 5 × 10−3 s, and 8000 steps at 1 × 10−2 s. Underrelaxation factors between 0.6 and 0.7 were adopted for all flow quantities, and pressure was never underrelaxed. The hybrid-upwind discretization scheme was used for the convective terms. At the inlet, gas, liquid, and the average volume fraction have been specified. At the pipe outlet, a relative average static pressure of zero was specified. For initiating the numerical solution, average volume fraction and parabolic liquid velocity profile are specified as initial conditions. The operating conditions are summarized in Table 2. The liquid and gas superficial velocities were varied between 3.74 to 5.1 m/s and 0.25 to 1.34 m/s, respectively. The fluid data are taken at room temperature (25°C) and are treated as isothermal and saturated. Therefore, heat and mass transfer effects are neglected in the simulations.Table 2
Operating conditions.
Geometry50.3 mm IDGas phaseAir at 25°CLiquid phaseWater at 25°CGas superficial velocity0.25–1.34 m/sLiquid superficial velocity3.74–5.1 m/sAverage gas volume fraction0.04–0.205Effect of grid size on gas volume fraction and axial liquid velocity.
(a)(b)


## 4. Results and Discussion

The CFD simulations are carried out for the experimental conditions reported by Kocamustafaogullari and Wang [17] and Kocamustafaogullari and Huang [1]. The local radial profiles of the gas volume fraction, interfacial area concentration, Sauter mean bubble diameter, and liquid velocities are predicted by solving the coupled two-fluid and population balance models. The predicted results are compared with the experimental data at the axial location of L/D=253 and along a vertical and horizontal line passing through the centre of the pipe axis. Here, y/D and x/R are the normalized vertical and horizontal positions in the pipe.


### 4.1. Estimation of Lift Coefficient

Accurate prediction of developing bubbly flows in horizontal pipes cannot be carried out without sufficient knowledge of a transverse lift force acting on the bubbles, the force that governs the transverse migration of a bubble in a shear field. It has been clarified through a number of experiments that the lateral migration strongly depends on bubble size, that is, small bubbles tend to migrate toward the pipe wall which results in a peak in the bubble volume fraction distribution near the wall. Just like the functional form of the drag coefficient for a single particle interaction is extended to multiparticle systems, the functional form of the lift force that captures the lateral migration phenomenon is given by (18) which has CL as the lift coefficient. Just as the drag coefficient is a function of local Reynolds number based on the slip velocity, one can expect the lift coefficient also to vary with local Reynolds number, and in general, it is an unknown function for such a complex flow field. In the literature, it is used as a fitting parameter, but various values have been reported. Further, the most of the correlations available in the literature were for vertical flows. The correlations of Legendre and Magnaudet [45], Tomiyama et al. [41] have been used in the simulations and the simulation results shown in Figure 2(a). It can be observed that the correlation of Legendre and Magnaudet [45], Tomiyama et al. [41] does not perform well for the horizontal flows because of most the dispersed bubbles migrate towards the top of the pipe, due to buoyancy. The negative lift coefficient needed because this force pushes bubbles to the pipe center. For the given simulation, we need a negative lift coefficient to predict near wall peak for the gas volume fraction profile. We need a correlation which gives negative lift coefficient value. In view of this, we have developed a correlation to be a function of Reynolds number.Lift coefficient (CL).
(a)(b)The difference between the model predictions and experimental data on the spatial variation of field quantities such as liquid velocity profiles, volume fraction profiles, and interfacial area measurements is minimized by the tuning of this parameter. Bubbly horizontal pipe flow experiments by Kocamustafaogullari and Wang [17]; Kocamustafaogullari and Huang [1]; Kocamustafaogullari et al. [46] were chosen for tuning the lift coefficient as they have detailed experimental data on the spatial variation of liquid velocity profiles, volume fraction profiles, and interfacial area measurements. In all the experiments, adiabatic, incompressible, air-water bubbly flows at atmospheric pressure and room temperature were used. The main result of tuning this parameter is shown in Figure 2(b). The estimated values of the lift coefficient at different experimental condition of gas and liquid flow rates were not scattered all over, but exhibited a well-defined correlation with the Reynolds number defined as Re=dsVsρL/μL, where dS is the average bubble diameter, and VS is the slip velocity. We capture this relationship by a polynomial expression of the form, CL = aRe3+ bRe2+ cRe + d, where a = −1 × 10−10, b = 2 × 10−7, c = 2 × 10−4, and d = −0.2937. It is worth noting that the correlation is based on the locally measured properties of turbulence as well as the bubble number density, and hence, one can expect it to be valid irrespective of the dimension of the pipe as well as the liquid system. Such a relationship can then be used back in the simulation for predictive purposes at other flow conditions.


### 4.2. Gas Volume Fraction

Figures3(a)–6(a) show the comparison of the predicted gas volume fraction with the experimental data of Kocamustafaogullari and Wang [17] for different superficial gas velocities of 0.25 m/s, 0.50 m/s, 0.80 m/s, and 1.34 m/s at a fixed liquid velocity of 5.1 m/s. Similarly Figures 7(a) and 8(a) show the gas volume fraction for liquid velocities of 3.74 m/s and 4.40 m/s at a fixed gas velocity of 0.51 m/s. The agreement between the predicted and the experimental profiles can be seen to be very good. As the superficial gas velocity increases, the average gas volume fraction also increases. It can be observed from these figures that most of the bubbles tend to migrate towards the top of the pipe wall under the dominating influence of buoyancy force. The balance of buoyancy and lift forces causes the profiles of gas volume fraction to show a distinct peak near the top wall at about y/D = 0.9 to 0.95 for all the flow conditions. A similar observation was made experimentally by Kocamustafaogullari and Wang [17], Kocamustafaogullari and Huang [1], Kocamustafaogullari et al., [46], and Iskandrani and Kojasoy [47]. At a constant gas superficial velocity of 0.51 m/s in Figures 7 and 8, the average and the peak value of the volume fraction decreases with increasing liquid velocity, as expected. The fact that the spatial variation of the gas volume fraction matches well with the experimental data over a wide range of flow conditions gives us confidence that the lift coefficient correlation that we have developed is quite appropriate. The real test of this correlation must of course await testing against similar data in a larger diameter pipe. The challenge in developing multiphase flow models using the volume-averaged framework is to develop such closure relationships for each of the individual mechanisms and test their validity under a wide variety of scales and flow conditions. We will be testing this correlation for bubble columns in the near future. The model prediction of gas volume fraction shows relative mean and maximum errors are ±6% and ±19%, respectively.Comparison of predicted and experimental data of Kocamustafaogullari and Wang [17] for superficial gas velocity of 0.25 m/s, superficial liquid velocity of 5.1 m/s, and volume fraction 0.043: (a) gas volume fraction, (b) interfacial area concentration, (c) Sauter mean bubble diameter, and (d) axial liquid velocity.
(a)(b)(c)(d)Comparison of predicted and experimental data of Kocamustafaogullari and Wang [17] for superficial gas velocity of 0.50 m/s, superficial liquid velocity of 5.1 m/s, and volume fraction of 0.080; (a) gas volume fraction, (b) interfacial area concentration, (c) Sauter mean bubble diameter, and (d) axial liquid velocity.
(a)(b)(c)(d)Comparison of predicted and experimental data of Kocamustafaogullari and Wang [17] for superficial gas velocity of 0.80 m/s, superficial liquid velocity of 5.1 m/s, and volume fraction of 0.139;  (a) gas volume fraction, (b) interfacial area concentration, (c) Sauter mean bubble diameter, and (d) axial liquid velocity.
(a)(b)(c)(d)Comparison of predicted and experimental data of Kocamustafaogullari and Wang [17] for superficial gas velocity of 1.34 m/s, superficial liquid velocity of 5.1 m/s, and volume fraction of 0.204; (a) gas volume fraction, (b) interfacial area concentration, (c) Sauter mean bubble diameter, and (d) axial liquid velocity.
(a)(b)(c)(d)Comparison of predicted and experimental data of Kocamustafaogullari and Wang [17] for superficial gas velocity of 0.51 m/s, superficial liquid velocity of 3.74 m/s, and volume fraction of 0.105: (a) gas volume fraction, (b) interfacial area concentration, (c) Sauter mean bubble diameter, and (d) axial liquid velocity.
(a)(b)(c)(d)Comparison of predicted and experimental data of Kocamustafaogullari and Wang [17] for superficial gas velocity of 0.51 m/s, superficial liquid velocity of 4.4 m/s, and volume fraction of 0.102; (a) gas volume fraction, (b) interfacial area concentration, (c) Sauter mean bubble diameter, and (d) axial liquid velocity.
(a)(b)(c)(d)


### 4.3. Interfacial Area Concentration (IAC)

The current simulation results and the experimental results of Kocamustafaogullari and Wang [17] on the local interfacial area concentration variation along the vertical direction are compared in Figures 3(b)–8(b). The flow conditions remain the same in the previous section. The CFD prediction shows good agreement with experimental data. From these figures, it can be seen that the interfacial area concentration shows characteristics that are similar to the gas volume fraction distribution. But the interfacial area depends not only on the volume fraction of the phase, but also on the bubble size distribution. Since the volume fraction and the interfacial area are independent measurements, the data on the interfacial area variation along the vertical direction provide a valuable test of the model predictions from the population balance models. Thus, the agreement seen with the gas volume fraction variation in the previous section provides a level of confidence in the lift coefficient model, while the agreement seen on the interfacial area measurements in the current section provides a level of comfort that the birth and death processes modeled in the population balance model are adequate to describe the bubble dynamics. Further, it can be seen that the local interfacial area concentration can be as high as 1000 m2/m3 towards the top of the pipe in horizontal two-phase flow. These values are quite high compared to vertical bubbly flows. This will result in increasing the intensity of the interfacial transport of mass, momentum, and heat near the top of the pipe. In addition, it can be observed that increasing the superficial gas velocity or decreasing the superficial liquid velocity would increase the local and overall interfacial area concentration and tend to flatten the interfacial area concentration profile. The model prediction of interfacial area concentration shows relative mean and maximum errors are ±8% and ±22%, respectively.


### 4.4. Sauter Mean Bubble Diameter

The comparison of predicted and experimental data of the local Sauter mean bubble diameter distribution is shown in Figures3(c)–8(c) for various superficial gas and liquid velocities. The Sauter mean bubble diameters are in the range of 1.5–5 mm, depending on the location and flow conditions. It should be noted that the experimental data on Sauter mean diameter is inferred from other measurements, and it is not a directly measured quantity. The scatter in the experimentally derived data is high, particularly in the lower region where the gas volume fraction is low, indicating that perhaps the signals are weaker in that region. Good agreement was achieved against the measured bubble size for all the experimental conditions. From these figures, it can be seen that the bubble size distribution is almost uniform in the pipe cross-section except near the wall region. The Suater mean bubble size tends to reduce close to the top of the pipe wall. This can be attributed to the fact that near the wall a very strong velocity gradient exists, which causes further break-up into smaller bubble sizes. Furthermore, the Suater mean bubble size is seen to increase with increasing the superficial gas velocity (Figures 3(c)–8(c)) and to decrease with increasing superficial liquid velocity (Figures 7(c)–8(c)). The simulation results capture all of these trends faithfully. The model prediction shows relative mean and maximum errors are ±9% and ±24%, respectively.


### 4.5. Axial Liquid Velocity

Figures3(d)–8(d) show the comparison of predicted and experimental data of axial liquid velocity profiles for different superficial gas and liquid velocities. If only a single liquid phase moves in the pipe, the liquid velocity in the pipe top region will be equal to the velocity in the bottom region, exhibiting a perfect axisymmetry. But these results show that the axial liquid velocity profile has a slight degree of asymmetry due to the presence of gas flow. The degree of asymmetry decreases with increasing liquid flow or decreasing gas flow. For increasingly higher gas velocities (Figures 3(d)–6(d)), the liquid velocity in the upper region of the pipe is slightly lower than in the lower region. This could be attributed to larger volume fraction of gas in the upper region which is the reason for the asymmetric distribution of the liquid velocity. An interesting feature of the velocity profile is that the velocity distribution within the bottom liquid layer resembles closely a fullydeveloped turbulent pipe flow profile irrespective of the liquid and gas superficial velocities. The model prediction of axial liquid velocity shows relative mean and maximum errors are ±5% and ±14%, respectively.


### 4.6. Simulation Results

From the simulation, we can get much more additional information, while some of these quantities are more difficult to measure in an experiment. One such quantity is the slip velocity between the two phases. The variation in the vertical direction of the slip velocity is shown in Figure9 for various combinations of gas and liquid flow rates. The slip velocity is larger in magnitude near the top region of the pipe, while a smaller slip velocity exists in the bottom part of the pipe. The slip velocity is an important characteristic of two-phase flow, particularly because of the large difference in densities between phases. Relatively smaller bubbles and fewer in number are found in the bottom region, and hence, they tend to move with the liquid resulting in a smaller slip velocity, while relatively larger bubbles and more in number are found near the top of the pipe, resulting in a larger slip velocity.Slip velocity at different superficial gas and liquid velocities: (a)VG=0.25 m/s, VL=5.1 m/s; (b) VG=0.50 m/s, VL=5.1 m/s; (c) VG=0.80 m/s, VL=5.1 m/s; (d) VG=1.34 m/s, VL=5.1 m/s; (e) VG=0.51 m/s, VL=3.74 m/s; (f) VG=0.51 m/s, VL=4.4 m/s.
(a)(b)(c)(d)(e)(f)


### 4.7. Average Interfacial Parameters

While we have used the data on the spatial variation of quantities such as volume fraction and interfacial area to tune the model parameters, from a practical view point one is often interested only in a quantity that is averaged over the pipe cross-section. Hence, area averaged gas volume fraction, interfacial area concentration, and mean bubble diameter at the exit plane are shown in Figure10 as a function of superficial gas velocity at various liquid velocities of 5.1 m/s, 4.4 m/s, and 3.74 m/s. The average volume fraction and interfacial area increase significantly with increasing superficial gas velocity, as expected. The influence of superficial liquid velocity on the gas volume fraction and interfacial area concentration are less significant. Figure 10(c) shows that the average bubble diameter increases slightly with increasing superficial gas velocity, all though the influence is not significant. However, the mean bubble size decreases with increasing superficial liquid velocity. This observation supports the fact that the bubble size is determined primarily by liquid flow turbulence in horizontal flows. Figure 10 compares the measured gas volume fraction, interfacial area, and Sauter bubble mean diameter values with those predicted using CFD-PBM model, and the relative mean and maximum errors are ±4% and ±11%, respectively.Effect of superficial gas and liquid velocity on (a) average gas volume fraction; (b) average interfacial area concentration (c) average mean bubble diameter: (1)VL=5.1 m/s; (2) VL=4.4 m/s; (3) VL=3.74 m/s.
(a)(b)(c)


### 4.8. Bubble Size Distribution

The bubble size distribution is determined by bubble coalescence and breakup. In a given system, bubble coalescence and breakup are primarily influenced by the local gas volume fraction and kinetic energy dissipation rate. Because of the nonuniform profiles of the gas volume fraction and dissipation rate, the bubble size distribution varies with the position as well. The spatial evolution of bubble size distribution between the inlet and the outlet of the pipe is shown in Figure11 for a superficial gas velocity of 0.25 m/s and a superficial liquid velocity of 4.67 m/s. While selecting the bubble size, we have considered a range from 1 to 10 mm, and the experimental observation of bubble size for highest superficial velocity is 6 mm. Figures 11(a), 11(c), and 11(e) show the bubble size distribution that was specified at the pipe inlet. These correspond, respectively, to the monosized bubbles of 1.45 mm (Figure 11(a)), 9.55 mm (Figure 11(c)), and a uniform distribution of bubbles in the range of 1 to 10 mm (Figure 11(e)). The corresponding distribution at the pipe exit is shown on the right hand side in Figures 11(b), 11(d), and 11(f), respectively. It can be seen from these figures that the bubble size distribution function reaches an independent state as determined by the balance between birth and death processes that depend on the local flow conditions, and its original state at the inlet has very little impact. Although this kind of distribution function was not measured in the experiments of Kocamustafaogullari and Wang [17], the spatial variation of the bubble sizes was measured as shown in Figure 3. It is comforting to note that the ranges of bubbles sizes measured under similar flow conditions show a range of 2-3 mm, the same range shown in Figure 11, even though extremely small (1.45 mm) and large (9.55 mm) sizes were used at the pipe inlet.The bubble class volume-based p.d.f at inlet and exit of the pipe for superficial gas velocity is 0.25 m/s, superficial liquid velocity is 4.67 m/s, and average volume fraction is 0.043.
(a)(b)(c)(d)(e)(f)


### 4.9. Development of Flow Pattern

To see the development of flow pattern in the axial direction, several three-dimensional simulations were carried out using the coupled two-fluid and population balance models. The flow evolution is shown in Figure12–14 at three different axial locations of L/D = 25, 148, and 253. The first location represents close to the entrance of the pipe region where the internal flow develops, and the second and third locations indicate the extent to which the flow has reached a fully developed state, by the lack of further change in flow profiles. Figures 12, 13, and 14 show, respectively, the development of the local gas volume fraction, interfacial area concentration, and axial liquid velocity in axial direction for the superficial gas velocity of 1.21 m/s and the superficial liquid velocity of 4.67 m/s. Good agreement can be seen between the predicted and experimental data at the axial location of L/D = 25. The gas volume fraction and interfacial area concentration do not show a significant variation in the vertical direction, near the entrance of the pipe (L/D = 25). This is because the bubble residence time was very small, and the transverse phase segregation due to the gravity has not been established yet. However, from first location (L/D = 25) to the second location (L/D = 148), the large differences can be observed. From second location (L/D = 148) to third location (L/D = 253), there is no significant difference was observed, but the fluid segregation due to the buoyancy is still effective. Further, it can be observed from Figure 14 that the axial liquid velocity profile shows nearly the same for all the locations. A slight change in the numerical values of the velocity can be attributed to the expansion of the gas phase associated with the frictional pressure gradient causing a continuous acceleration of the mixture in the axial direction.Gas volume fraction development in axial direction for superficial gas velocity is 1.21 m/s, superficial liquid velocity is 4.67 m/s, and average volume fraction is 0.205. At vertical position, (a)L/D = 25; (b) L/D = 148; (c) L/D = 253, at horizontal position, (d) L/D = 25; (e) L/D = 148; (f) L/D = 253.
(a)(b)(c)(d)(e)(f)Interfacial area concentration (IAC) development in axial direction for superficial gas velocity is 1.21 m/s, superficial liquid velocity is 4.67 m/s, and average volume fraction is 0.205. At vertical position, (a)L/D= 25; (b) L/D= 148; (c) L/D= 253, at horizontal position, (d) L/D= 25; (e) L/D= 148; (f) L/D= 253.
(a)(b)(c)(d)(e)(f)Liquid velocity development in axial direction for superficial gas velocity is 1.21 m/s, superficial liquid velocity is 4.67 m/s, and average volume fraction is 0.205. At vertical position, (a)L/D= 25; (b) L/D= 148; (c) L/D= 253, at horizontal position, (d) L/D= 25; (e) L/D= 148; (f) L/D= 253.
(a)(b)(c)(d)(e)(f)


## 4.1. Estimation of Lift Coefficient

Accurate prediction of developing bubbly flows in horizontal pipes cannot be carried out without sufficient knowledge of a transverse lift force acting on the bubbles, the force that governs the transverse migration of a bubble in a shear field. It has been clarified through a number of experiments that the lateral migration strongly depends on bubble size, that is, small bubbles tend to migrate toward the pipe wall which results in a peak in the bubble volume fraction distribution near the wall. Just like the functional form of the drag coefficient for a single particle interaction is extended to multiparticle systems, the functional form of the lift force that captures the lateral migration phenomenon is given by (18) which has CL as the lift coefficient. Just as the drag coefficient is a function of local Reynolds number based on the slip velocity, one can expect the lift coefficient also to vary with local Reynolds number, and in general, it is an unknown function for such a complex flow field. In the literature, it is used as a fitting parameter, but various values have been reported. Further, the most of the correlations available in the literature were for vertical flows. The correlations of Legendre and Magnaudet [45], Tomiyama et al. [41] have been used in the simulations and the simulation results shown in Figure 2(a). It can be observed that the correlation of Legendre and Magnaudet [45], Tomiyama et al. [41] does not perform well for the horizontal flows because of most the dispersed bubbles migrate towards the top of the pipe, due to buoyancy. The negative lift coefficient needed because this force pushes bubbles to the pipe center. For the given simulation, we need a negative lift coefficient to predict near wall peak for the gas volume fraction profile. We need a correlation which gives negative lift coefficient value. In view of this, we have developed a correlation to be a function of Reynolds number.Lift coefficient (CL).
(a)(b)The difference between the model predictions and experimental data on the spatial variation of field quantities such as liquid velocity profiles, volume fraction profiles, and interfacial area measurements is minimized by the tuning of this parameter. Bubbly horizontal pipe flow experiments by Kocamustafaogullari and Wang [17]; Kocamustafaogullari and Huang [1]; Kocamustafaogullari et al. [46] were chosen for tuning the lift coefficient as they have detailed experimental data on the spatial variation of liquid velocity profiles, volume fraction profiles, and interfacial area measurements. In all the experiments, adiabatic, incompressible, air-water bubbly flows at atmospheric pressure and room temperature were used. The main result of tuning this parameter is shown in Figure 2(b). The estimated values of the lift coefficient at different experimental condition of gas and liquid flow rates were not scattered all over, but exhibited a well-defined correlation with the Reynolds number defined as Re=dsVsρL/μL, where dS is the average bubble diameter, and VS is the slip velocity. We capture this relationship by a polynomial expression of the form, CL = aRe3+ bRe2+ cRe + d, where a = −1 × 10−10, b = 2 × 10−7, c = 2 × 10−4, and d = −0.2937. It is worth noting that the correlation is based on the locally measured properties of turbulence as well as the bubble number density, and hence, one can expect it to be valid irrespective of the dimension of the pipe as well as the liquid system. Such a relationship can then be used back in the simulation for predictive purposes at other flow conditions.


## 4.2. Gas Volume Fraction

Figures3(a)–6(a) show the comparison of the predicted gas volume fraction with the experimental data of Kocamustafaogullari and Wang [17] for different superficial gas velocities of 0.25 m/s, 0.50 m/s, 0.80 m/s, and 1.34 m/s at a fixed liquid velocity of 5.1 m/s. Similarly Figures 7(a) and 8(a) show the gas volume fraction for liquid velocities of 3.74 m/s and 4.40 m/s at a fixed gas velocity of 0.51 m/s. The agreement between the predicted and the experimental profiles can be seen to be very good. As the superficial gas velocity increases, the average gas volume fraction also increases. It can be observed from these figures that most of the bubbles tend to migrate towards the top of the pipe wall under the dominating influence of buoyancy force. The balance of buoyancy and lift forces causes the profiles of gas volume fraction to show a distinct peak near the top wall at about y/D = 0.9 to 0.95 for all the flow conditions. A similar observation was made experimentally by Kocamustafaogullari and Wang [17], Kocamustafaogullari and Huang [1], Kocamustafaogullari et al., [46], and Iskandrani and Kojasoy [47]. At a constant gas superficial velocity of 0.51 m/s in Figures 7 and 8, the average and the peak value of the volume fraction decreases with increasing liquid velocity, as expected. The fact that the spatial variation of the gas volume fraction matches well with the experimental data over a wide range of flow conditions gives us confidence that the lift coefficient correlation that we have developed is quite appropriate. The real test of this correlation must of course await testing against similar data in a larger diameter pipe. The challenge in developing multiphase flow models using the volume-averaged framework is to develop such closure relationships for each of the individual mechanisms and test their validity under a wide variety of scales and flow conditions. We will be testing this correlation for bubble columns in the near future. The model prediction of gas volume fraction shows relative mean and maximum errors are ±6% and ±19%, respectively.Comparison of predicted and experimental data of Kocamustafaogullari and Wang [17] for superficial gas velocity of 0.25 m/s, superficial liquid velocity of 5.1 m/s, and volume fraction 0.043: (a) gas volume fraction, (b) interfacial area concentration, (c) Sauter mean bubble diameter, and (d) axial liquid velocity.
(a)(b)(c)(d)Comparison of predicted and experimental data of Kocamustafaogullari and Wang [17] for superficial gas velocity of 0.50 m/s, superficial liquid velocity of 5.1 m/s, and volume fraction of 0.080; (a) gas volume fraction, (b) interfacial area concentration, (c) Sauter mean bubble diameter, and (d) axial liquid velocity.
(a)(b)(c)(d)Comparison of predicted and experimental data of Kocamustafaogullari and Wang [17] for superficial gas velocity of 0.80 m/s, superficial liquid velocity of 5.1 m/s, and volume fraction of 0.139;  (a) gas volume fraction, (b) interfacial area concentration, (c) Sauter mean bubble diameter, and (d) axial liquid velocity.
(a)(b)(c)(d)Comparison of predicted and experimental data of Kocamustafaogullari and Wang [17] for superficial gas velocity of 1.34 m/s, superficial liquid velocity of 5.1 m/s, and volume fraction of 0.204; (a) gas volume fraction, (b) interfacial area concentration, (c) Sauter mean bubble diameter, and (d) axial liquid velocity.
(a)(b)(c)(d)Comparison of predicted and experimental data of Kocamustafaogullari and Wang [17] for superficial gas velocity of 0.51 m/s, superficial liquid velocity of 3.74 m/s, and volume fraction of 0.105: (a) gas volume fraction, (b) interfacial area concentration, (c) Sauter mean bubble diameter, and (d) axial liquid velocity.
(a)(b)(c)(d)Comparison of predicted and experimental data of Kocamustafaogullari and Wang [17] for superficial gas velocity of 0.51 m/s, superficial liquid velocity of 4.4 m/s, and volume fraction of 0.102; (a) gas volume fraction, (b) interfacial area concentration, (c) Sauter mean bubble diameter, and (d) axial liquid velocity.
(a)(b)(c)(d)


## 4.3. Interfacial Area Concentration (IAC)

The current simulation results and the experimental results of Kocamustafaogullari and Wang [17] on the local interfacial area concentration variation along the vertical direction are compared in Figures 3(b)–8(b). The flow conditions remain the same in the previous section. The CFD prediction shows good agreement with experimental data. From these figures, it can be seen that the interfacial area concentration shows characteristics that are similar to the gas volume fraction distribution. But the interfacial area depends not only on the volume fraction of the phase, but also on the bubble size distribution. Since the volume fraction and the interfacial area are independent measurements, the data on the interfacial area variation along the vertical direction provide a valuable test of the model predictions from the population balance models. Thus, the agreement seen with the gas volume fraction variation in the previous section provides a level of confidence in the lift coefficient model, while the agreement seen on the interfacial area measurements in the current section provides a level of comfort that the birth and death processes modeled in the population balance model are adequate to describe the bubble dynamics. Further, it can be seen that the local interfacial area concentration can be as high as 1000 m2/m3 towards the top of the pipe in horizontal two-phase flow. These values are quite high compared to vertical bubbly flows. This will result in increasing the intensity of the interfacial transport of mass, momentum, and heat near the top of the pipe. In addition, it can be observed that increasing the superficial gas velocity or decreasing the superficial liquid velocity would increase the local and overall interfacial area concentration and tend to flatten the interfacial area concentration profile. The model prediction of interfacial area concentration shows relative mean and maximum errors are ±8% and ±22%, respectively.


## 4.4. Sauter Mean Bubble Diameter

The comparison of predicted and experimental data of the local Sauter mean bubble diameter distribution is shown in Figures3(c)–8(c) for various superficial gas and liquid velocities. The Sauter mean bubble diameters are in the range of 1.5–5 mm, depending on the location and flow conditions. It should be noted that the experimental data on Sauter mean diameter is inferred from other measurements, and it is not a directly measured quantity. The scatter in the experimentally derived data is high, particularly in the lower region where the gas volume fraction is low, indicating that perhaps the signals are weaker in that region. Good agreement was achieved against the measured bubble size for all the experimental conditions. From these figures, it can be seen that the bubble size distribution is almost uniform in the pipe cross-section except near the wall region. The Suater mean bubble size tends to reduce close to the top of the pipe wall. This can be attributed to the fact that near the wall a very strong velocity gradient exists, which causes further break-up into smaller bubble sizes. Furthermore, the Suater mean bubble size is seen to increase with increasing the superficial gas velocity (Figures 3(c)–8(c)) and to decrease with increasing superficial liquid velocity (Figures 7(c)–8(c)). The simulation results capture all of these trends faithfully. The model prediction shows relative mean and maximum errors are ±9% and ±24%, respectively.


## 4.5. Axial Liquid Velocity

Figures3(d)–8(d) show the comparison of predicted and experimental data of axial liquid velocity profiles for different superficial gas and liquid velocities. If only a single liquid phase moves in the pipe, the liquid velocity in the pipe top region will be equal to the velocity in the bottom region, exhibiting a perfect axisymmetry. But these results show that the axial liquid velocity profile has a slight degree of asymmetry due to the presence of gas flow. The degree of asymmetry decreases with increasing liquid flow or decreasing gas flow. For increasingly higher gas velocities (Figures 3(d)–6(d)), the liquid velocity in the upper region of the pipe is slightly lower than in the lower region. This could be attributed to larger volume fraction of gas in the upper region which is the reason for the asymmetric distribution of the liquid velocity. An interesting feature of the velocity profile is that the velocity distribution within the bottom liquid layer resembles closely a fullydeveloped turbulent pipe flow profile irrespective of the liquid and gas superficial velocities. The model prediction of axial liquid velocity shows relative mean and maximum errors are ±5% and ±14%, respectively.


## 4.6. Simulation Results

From the simulation, we can get much more additional information, while some of these quantities are more difficult to measure in an experiment. One such quantity is the slip velocity between the two phases. The variation in the vertical direction of the slip velocity is shown in Figure9 for various combinations of gas and liquid flow rates. The slip velocity is larger in magnitude near the top region of the pipe, while a smaller slip velocity exists in the bottom part of the pipe. The slip velocity is an important characteristic of two-phase flow, particularly because of the large difference in densities between phases. Relatively smaller bubbles and fewer in number are found in the bottom region, and hence, they tend to move with the liquid resulting in a smaller slip velocity, while relatively larger bubbles and more in number are found near the top of the pipe, resulting in a larger slip velocity.Slip velocity at different superficial gas and liquid velocities: (a)VG=0.25 m/s, VL=5.1 m/s; (b) VG=0.50 m/s, VL=5.1 m/s; (c) VG=0.80 m/s, VL=5.1 m/s; (d) VG=1.34 m/s, VL=5.1 m/s; (e) VG=0.51 m/s, VL=3.74 m/s; (f) VG=0.51 m/s, VL=4.4 m/s.
(a)(b)(c)(d)(e)(f)


## 4.7. Average Interfacial Parameters

While we have used the data on the spatial variation of quantities such as volume fraction and interfacial area to tune the model parameters, from a practical view point one is often interested only in a quantity that is averaged over the pipe cross-section. Hence, area averaged gas volume fraction, interfacial area concentration, and mean bubble diameter at the exit plane are shown in Figure10 as a function of superficial gas velocity at various liquid velocities of 5.1 m/s, 4.4 m/s, and 3.74 m/s. The average volume fraction and interfacial area increase significantly with increasing superficial gas velocity, as expected. The influence of superficial liquid velocity on the gas volume fraction and interfacial area concentration are less significant. Figure 10(c) shows that the average bubble diameter increases slightly with increasing superficial gas velocity, all though the influence is not significant. However, the mean bubble size decreases with increasing superficial liquid velocity. This observation supports the fact that the bubble size is determined primarily by liquid flow turbulence in horizontal flows. Figure 10 compares the measured gas volume fraction, interfacial area, and Sauter bubble mean diameter values with those predicted using CFD-PBM model, and the relative mean and maximum errors are ±4% and ±11%, respectively.Effect of superficial gas and liquid velocity on (a) average gas volume fraction; (b) average interfacial area concentration (c) average mean bubble diameter: (1)VL=5.1 m/s; (2) VL=4.4 m/s; (3) VL=3.74 m/s.
(a)(b)(c)


## 4.8. Bubble Size Distribution

The bubble size distribution is determined by bubble coalescence and breakup. In a given system, bubble coalescence and breakup are primarily influenced by the local gas volume fraction and kinetic energy dissipation rate. Because of the nonuniform profiles of the gas volume fraction and dissipation rate, the bubble size distribution varies with the position as well. The spatial evolution of bubble size distribution between the inlet and the outlet of the pipe is shown in Figure11 for a superficial gas velocity of 0.25 m/s and a superficial liquid velocity of 4.67 m/s. While selecting the bubble size, we have considered a range from 1 to 10 mm, and the experimental observation of bubble size for highest superficial velocity is 6 mm. Figures 11(a), 11(c), and 11(e) show the bubble size distribution that was specified at the pipe inlet. These correspond, respectively, to the monosized bubbles of 1.45 mm (Figure 11(a)), 9.55 mm (Figure 11(c)), and a uniform distribution of bubbles in the range of 1 to 10 mm (Figure 11(e)). The corresponding distribution at the pipe exit is shown on the right hand side in Figures 11(b), 11(d), and 11(f), respectively. It can be seen from these figures that the bubble size distribution function reaches an independent state as determined by the balance between birth and death processes that depend on the local flow conditions, and its original state at the inlet has very little impact. Although this kind of distribution function was not measured in the experiments of Kocamustafaogullari and Wang [17], the spatial variation of the bubble sizes was measured as shown in Figure 3. It is comforting to note that the ranges of bubbles sizes measured under similar flow conditions show a range of 2-3 mm, the same range shown in Figure 11, even though extremely small (1.45 mm) and large (9.55 mm) sizes were used at the pipe inlet.The bubble class volume-based p.d.f at inlet and exit of the pipe for superficial gas velocity is 0.25 m/s, superficial liquid velocity is 4.67 m/s, and average volume fraction is 0.043.
(a)(b)(c)(d)(e)(f)


## 4.9. Development of Flow Pattern

To see the development of flow pattern in the axial direction, several three-dimensional simulations were carried out using the coupled two-fluid and population balance models. The flow evolution is shown in Figure12–14 at three different axial locations of L/D = 25, 148, and 253. The first location represents close to the entrance of the pipe region where the internal flow develops, and the second and third locations indicate the extent to which the flow has reached a fully developed state, by the lack of further change in flow profiles. Figures 12, 13, and 14 show, respectively, the development of the local gas volume fraction, interfacial area concentration, and axial liquid velocity in axial direction for the superficial gas velocity of 1.21 m/s and the superficial liquid velocity of 4.67 m/s. Good agreement can be seen between the predicted and experimental data at the axial location of L/D = 25. The gas volume fraction and interfacial area concentration do not show a significant variation in the vertical direction, near the entrance of the pipe (L/D = 25). This is because the bubble residence time was very small, and the transverse phase segregation due to the gravity has not been established yet. However, from first location (L/D = 25) to the second location (L/D = 148), the large differences can be observed. From second location (L/D = 148) to third location (L/D = 253), there is no significant difference was observed, but the fluid segregation due to the buoyancy is still effective. Further, it can be observed from Figure 14 that the axial liquid velocity profile shows nearly the same for all the locations. A slight change in the numerical values of the velocity can be attributed to the expansion of the gas phase associated with the frictional pressure gradient causing a continuous acceleration of the mixture in the axial direction.Gas volume fraction development in axial direction for superficial gas velocity is 1.21 m/s, superficial liquid velocity is 4.67 m/s, and average volume fraction is 0.205. At vertical position, (a)L/D = 25; (b) L/D = 148; (c) L/D = 253, at horizontal position, (d) L/D = 25; (e) L/D = 148; (f) L/D = 253.
(a)(b)(c)(d)(e)(f)Interfacial area concentration (IAC) development in axial direction for superficial gas velocity is 1.21 m/s, superficial liquid velocity is 4.67 m/s, and average volume fraction is 0.205. At vertical position, (a)L/D= 25; (b) L/D= 148; (c) L/D= 253, at horizontal position, (d) L/D= 25; (e) L/D= 148; (f) L/D= 253.
(a)(b)(c)(d)(e)(f)Liquid velocity development in axial direction for superficial gas velocity is 1.21 m/s, superficial liquid velocity is 4.67 m/s, and average volume fraction is 0.205. At vertical position, (a)L/D= 25; (b) L/D= 148; (c) L/D= 253, at horizontal position, (d) L/D= 25; (e) L/D= 148; (f) L/D= 253.
(a)(b)(c)(d)(e)(f)


## 5. Conclusions

A two-fluid model coupled with population balance approach is presented in this paper to handle gas-liquid bubbly flows in horizontal pipe. To demonstrate the application of the population balance approach, the average bubble number density transport equation was formulated and implemented for gas-liquid bubbly flows in the CFD code CFX 5.7 to determine the temporal and spatial geometrical changes of the gas bubbles. Population balance combined with coalescence and break-up models were taken into consideration. A detailed comparison has been presented between the CFD simulation and the experimental data reported by Kocamustafaogullari and Wang [17] and Kocamustafaogullari and Huang [1]. Good agreement was seen between the predicted and the experimental data of the volume fraction, interfacial area concentration, Sauter mean bubble diameter, and liquid velocity for a range of superficial gas (0.25 to 1.34 m/s) and liquid (3.74 to 5.1 m/s) velocities and volume fraction (4 to 21%). The experimental and simulated results indicate that the volume fraction and interfacial area concentration have local maxima near the upper pipe wall, and the profiles tend to flatten with increasing liquid flow rate. It was observed that the mean bubble diameter ranged from 1.5 to 5 mm, depending on the location and flow conditions. Further, it was found that increasing the gas flow rate at fixed liquid flow rate would increase the local volume fraction and interfacial area concentration. The simulation results were consistent with experimental observed from the literature. Further, the development of flow pattern was examined at three axial locations L/D = 25, 148, and 253. It was found that the prediction shows good agreement with experimental data. The axial liquid mean velocity showed a relatively uniform distribution except near the upper pipe wall. The flow in the bottom part of the pipe exhibits a fully developed turbulent pipe flow profile, whereas in the top of the pipe a different flow exists.



---

*Source: 620463-2012-04-02.xml*