Chemical Engineering Research and Design On the prediction of the phase distribution of bubbly ﬂow in a horizontal pipe

Horizontal bubbly ﬂow is widely encountered in various industrial systems because of its ability to provide large interfacial areas for heat and mass transfer. Nonetheless, this particular ﬂow orientation has received less attention when compared to vertical bubbly ﬂow. Owing to the strong inﬂuence due to buoyancy, the migration of dispersed bubbles towards the top wall of the horizontal pipe generally causes a highly asymmetrical internal phase distributions, which are not experienced in vertical bubbly ﬂow. In this study, the internal phase distribution of air-water bubbly ﬂow in a long horizontal pipe with an inner diameter of 50.3 mm has been predicted using the population balance model based on direct quadrature method of moments (DQMOM) and multiple-size group (MUSIG) model. The predicted local radial distributions of gas void fraction, liquid velocity and interfacial area concentration have been validated against the experimental data of Kocamustafaogullari and Huang (1994). In general, satisfactory agree-ments between predicted and measured results were achieved. The numerical results indicated that the gas void fraction and interfacial area concentration have a unique internal structure with a prevailing maximum peak near the top wall of the pipe due to buoyancy effect.


Introduction
In many technological systems, gas-liquid flows are immensely important in the chemical and process industries. Among all the flow regimes that are prevalent in gas-liquid flows, bubbly flow conditions are of the greatest interest because of the capacity of handling processes requiring large scale interfacial areas for heat and mass transfer and efficient mixing processes. Most computational studies of bubbly flows have thus far concentrated on vertical pipe configurations. In vertical flows, gravity mainly affects the axial gas-to-liquid relative velocity, but does not induce any lateral symmetry in either velocity or phase distribution. Nevertheless, in the case of horizontal (or even inclined) flows, the acceleration of gravity not only causes a significant flow asymmetry but also imposes an addi-tional strong radial force. Thereby, under the combination of radial and axial forces, bubbles can travel neither vertically nor horizontally, which increases the difficulty in modelling horizontal bubbly flow in comparison to vertical bubbly flow. Also, the density stratification is often accompanied by a strong secondary flow.
Several measurement studies have been performed to describe the flow patterns in horizontal pipe flow. Govier and Aziz (1972) have classified the flow patterns into five groups, namely, bubbly, plug, slug, wave and annular. Taitel and Dukler (1976) and Weisman et al. (1979) have mapped these flow regimes in a two-dimensional coordinate system and predicted their transition for numerous fluid properties and pipe sizes. Nonetheless, many researchers have focused on the internal structure of horizontal bubbly flow. Some examples are Kocamustafaogullari and Wang (1991), , , Andreussi et al. (1999), Iskandrani  k , B C k mass birth rate due to break-up and coalescence C break-up model constant  Sanders et al. (2004) and Yang et al. (2004). Recently, interfacial structure of horizontal bubbly flow has been observed in 45-degree and 90-degree elbow by Kim et al. (2007Kim et al. ( , 2009). Haoues et al. (2009) and Talley and Kim (2010) have developed a drift flux model to predict the integral flow characteristics of horizontal bubbly flow. Tselishcheva et al. (2010) applied the two-fluid model to simulate the void fraction and velocity profiles in a long straight horizontal pipe and a similar pipe with a 90-degree elbow. The local spatial two-phase geometrical internal structure (bubble diameter or interfacial area concentration) in such flow was found to be affected by the coalescence and break-up through the interactions among bubbles as well as between bubbles and turbulent eddies in turbulent flows. In order to predict the bubble size distribution, the population balance approach has been applied in order to accommodate the complicated bubble interaction mechanisms coupled with the two-fluid model. Ekambara et al. (2008) have applied the MUSIG model while Li et al. (2010) have adopted the Average Bubble Number Density (ABND) approach to investigate the internal phase distribution of a horizontal bubbly flow. For bubble-bubble and bubble-turbulence interactions, the major phenomenological mechanisms in bubbly flow conditions have been identified: (a) coalescence through random collision driven by turbulent eddies, (b) coalescence due to the acceleration of the following bubble in the wake of the preceding bubble, and (c) break-up due to the impact of turbulent eddies. Appropriate mechanistic models by Wu et al. (1998), Hibiki and Ishii (2002) and Yao and Morel (2004) for ABND and Prince and Blanch (1990), Chesters (1991), Luo and Svendsen (1996), Martinez-Bazan et al. (1999a,b), Lehr et al. (2002) and Lo and Zhang (2009) for multi-size bubble consideration have been established.
In the present of work, the method of moments (MOM) is adopted to predict the phase distribution of a horizontal bubbly flow. In contrast to the MUSIG model (Cheung et al., 2007a(Cheung et al., ,b, 2008 which belongs to the class method, the bubble size distribution is tracked through its moments. The main advantage of this approach is essentially the computational economy that it entails. In contract to the MUSIG model where multiple classes (i.e. 10 classes or more) are normally required, the method condenses the problem by only tracking the evolution of a small number of moments (normally 4-6). This becomes ever crucial in the modelling of complex bubbly flows when the bubble dynamics is strongly coupled with already time-consuming calculations of turbulence multiphase flows. Owing to the difficulties experienced in expressing the transport equations in terms of the moments themselves, the DQMOM proposed by Marchisio and Fox (2005) has been applied. From the specific aspect in the modelling of bubbly flow, each node of the quadrature approximation is treated as a distinct gas phase. DQMOM, similar to MUSIG, can offer a powerful approach in describing bubbly flow undergoing coalescence and break-up processes in the context of computational fluid dynamic (CFD) simulations.

Two-fluid model
For isothermal horizontal bubbly flow, the equations for the ensemble-averaged of continuity and momentum governing each phase are solved simultaneously. By denoting the liquid as the continuous phase (˛l) and the gas (i.e. bubbles) as disperse phase (˛g), these equations can be written as:

Momentum equation of liquid phase
∂ ∂t ( l˛l u l ) + ∇ · ( l˛l u l u l ) = −˛l∇p + ˛l l g + ∇ · [˛l l e (∇u l + (∇u l )

Momentum equation of gas phase
∂ ∂t ( g˛g u g ) + ∇ · ( g˛g u g u g ) = −˛g∇p + ˛g g g + ∇ · [˛g g e (∇u g + (∇u g ) The total interfacial force F lg appearing in the liquid phase momentum equation (3) is formulated according to the appropriate consideration of different sub-forces affecting the interface between each phase. As demonstrated by Frank et al. (2004), the total interfacial force for the liquid phase is given by the drag, lift, wall lubrication and turbulent dispersion: Note that the total interfacial force in the gas phase momentum equation (4) is given by F gl = −F lg . The inter-phase momentum transfer between gas and liquid due to the drag force can be obtained via the drag model of Ishii and Zuber (1979) as: where a if is the interfacial area concentration (=6˛g/D s ). The drag coefficient C D in Eq. (6) has been correlated for sev-eral distinct Reynolds number regions for individual bubbles according to Ishii and Zuber (1979). The non-drag forces are those of lift, wall lubrication and turbulent dispersion which these forces are directed perpendicular to the flow direction. Based on the formulation of Drew and Lahey (1979), the lift force in terms of the slip velocity and the curl of the liquid phase velocity can be described according to: Wall lubrication force as proposed by Antal et al. (1991), which acts in the normal direction away from the wall and decays with distance from the wall, can be expressed by: The turbulence induced dispersion developed by Lopez de Bertodano (1998) can be written in the form of: Alternatively, the turbulence induced dispersion based on the consistency of Favre-averaging developed by Burns et al. (2004) has also been applied, viz., In handling bubble induced turbulent flow, unlike single phase fluid flow problem, no standard turbulence model is tailored for bubbly flow. The standard k − ε model based on Launder and Spalding (1972) is considered for the continuous liquid phase with additional buoyancy turbulence source term being accounted within the transport equations. The governing equations for the turbulent kinetic energy k and turbulent dissipation ε are: where P is the shear production defined by: The turbulent viscosity of the liquid phase is given by: where the shear-induced turbulence is given by while the bubble-induced turbulence which accounts for the effect of bubbles on liquid turbulence, based on Sato's bubbleinduced turbulent viscosity model (Sato et al., 1981), can be expressed as The constants for the aforementioned k − ε model are: C = 0.09, k = 1.0, ε = 1.3, C ε1 = 1.44, C ε2 = 1.92, and C p = 1.2. For the gas phase, the dispersed phase zero equation model is utilized; the turbulent viscosity of gas phase can otherwise be obtained as: In Eqs. (3) and (4), the effective viscosity for the liquid and gas phases are subsequently given by:

Population balance models
For the purpose of predicting the Particle Size Distribution (PSD) within the bubbly flow, the population balance models based on the multiple-size group and quadrature approximation of the moment approaches are described below.

DQMOM for bubbly flow
The moments of the PSD can be defined as In general, the first few moments provide important statistical descriptions on the population of particles which can be related directly to some physical quantities. One approach for computing the moment is to approximate the integrals in Eq.
(20) using the numerical quadrature scheme -the quadrature method of moment (QMOM) as suggested by McGraw (1997). Instead of space transformation, the Gaussian quadrature closure is henceforth adopted to approximate the PSD according to a finite set of Dirac's delta functions. For bubbly flows, taking the bubble mass M as the preferred internal coordinate, the PSD takes the form: where N i represents the number density or weight of the ith class and consists of all particles (bubbles) per unit volume with a pivot size or abscissa M i . With the aim of solving multi-dimensional problems, Marchisio and Fox (2005) extended the MOM by developing the DQMOM where the quadrature abscissas and weights are formulated as transport equations. The main idea of the method is to keep track of the primitive variables appearing in the quadrature approximation, instead of the actual moments of the PSD. As a result, the evaluation of the abscissas and weights are obtained using matrix operations; the transport equations for weights and abscissas be derived after some mathematical manipulation according to: where i = N i M i is the weighted abscissas and the terms a i and b i are related to the birth and death rate of population which forms 2n linear equations where the unknowns can be evaluated via matrix inversion according to The 2n × 2n coefficient matrix A = [A 1 A 2 ] in Eq. (24) takes the form: where the 2n vector of unknowns comprises essentially the In order to maintain consistency with the variables employed in the two-fluid model, the weights and abscissas can be related to the size fraction of the dispersed phase (f i ) and an additional variable defined as i = f i /M i . Considering a homogeneous system in the present study in which the bubbles are assumed to travel with a common gas velocity (V g i ≡ u g ), the size fraction of f i is thus related to the weights and abscissas by Using the above expression and the variable i , the transport equations become The moment transform of the coalescence and break-up of the term S k can be expressed as where the terms B and D represent the birth and death rates of the coalescence and break-up of bubbles. From above, the weights N i and N j can be determined according to the definition given in Eq. (27).

MUSIG model for bubbly flow
Instead of inferring the PSD to derivative variables (i.e. moments), the class method directly simulates its main characteristics using primitive variables (i.e. particle number density). Through this consideration, the continuous size range of particles is discretised into a series number of discrete size classes. For each class, a scalar (number density of particles) equation is solved to accommodate the population changes caused by intra/inter-group particle coalescence and breakage. The PSD is thereby approximated as: which incidentally is the same expression as proposed for the QMOM in Eq. (21). However, the groups (or abscissas) of class methods are now fixed and aligned continuously in the state space.
Assuming that each bubble class travel at the common velocity as in DQMOM, the number density equation Kumar and Ramkrishna (1996) can be expressed as In a form consistent with the variables used in two-fluid model, the transport equation in terms of size fraction becomes: The birth and death rates can subsequently be written in terms of the size fraction according to:

Coalescence and break-up kernels for bubbly flow
According to Luo and Svendsen (1996), the bubble break-up rate can be modelled based on the assumption of bubble binary break-up under isotropic turbulence situation. Denoting the increase coefficient of surface area as c where f BV is stochastic break-up volume fraction of the daughter size distribution, the break-up rate in terms of mass can be obtained as where = /d j is the size ratio between an eddy and a particle in the inertial sub-range and consequently min = min /d j and C and ˇ are determined from fundamental consideration of drops or bubbles break-up in turbulent dispersion systems to be 0.923 and 2.0. Bubble coalescence occurs via collision of two bubbles which may be caused by wake entrainment, turbulence random collision and buoyancy. In this present study, the turbulence random collision is only considered; all bubbles are assumed to be of spherical shape (i.e. wake entrainment becomes therefore negligible). The coalescence rate considering the turbulent collision taken from Prince and Blanch (1990) in terms of mass can be expressed as: where ij is the contact time for two bubbles given by (d ij /2) 2/3 /(ε l ) 1/3 and t ij is the time required for two bubbles to coalesce having diameters d i and d j estimated to be [(d ij /2) 3 l /16 ] 0.5 ln(h 0 /h f ). The equivalent diameter d ij is calculated as suggested by Chesters and Hoffman (1982): d ij = (2/d i + 2/d j ) −1 . According to Prince and Blanch (1990), the initial film thickness h o = 1 × 10 −4 m and critical film thickness h f = 1 × 10 −8 m have been assumed. The turbulent velocity u t in the inertial sub-range of isotropic turbulence (Rotta, 1972) which is given by: u t = √ 2(ε l ) 1/3 d i 1/3 . In Eqs. (42) and (43), F B and F C denote the calibration factors for break-up and coalescence.

Experimental details
The internal phase distributions of co-current airwater bubbly flow experiment have been carried out by    still developing while the second and third measuring planes depicted near fully developed flow pattern. Radial profiles were measured at 3 angles: Â = 0 • , 45 • and 90 • where the angle Â was measured from the upper wall of the horizontal pipe orientation. A range of superficial liquid velocities j f and superficial gas velocities j g were performed. More details regarding the experimental set-up can be referred in . The primary aim in this present study has been to assess the performance of DQMOM and compared against those of MUSIG in simulating the bubbly flow region at the nearly fully developed state of L/D = 253. Details of the flow conditions within the bubbly flow regime are summarized in Table 1.

Numerical details
For the two-fluid model, solutions to equations governing conservation of mass and momentum for each phase were obtained via the ANSYS Inc, CFX-5.7.1 computer code. For DQMOM, transport equations governing weights and abscissas were solved to predict the bubble size distribution of which the evaluation of the source terms a i and b i has been achieved through the utilization of a user subroutine incor-porated within the CFD computer code. With regards to the coalescence and break-up kernels, calibration factors F B and F C have been adjusted to 0.15 and 0.05 for DQMOM as well as MUSIG model in this present study to match the experimental data. For comparison purpose, these dimensionless factors were kept constant for all flow conditions during the entire numerical calculations. For DQMOM, four moments were adopted to explicitly track the distribution of bubble sizes ranging from 0 mm to 10 mm. For MUSIG model, 10 bubble classes were, however, equally divided for bubble sizes between 0 mm and 10 mm. For simplicity, all moments or bubbles classes in DQMOM and MUSIG model were assumed to travel in the same gas velocity which has been solved explicitly from the gas-phase of the two-fluid model.
Inlet conditions were assumed to be homogeneous with regards to the superficial liquid and gas velocities, void fractions for both phases and uniformly distributed bubble size in accordance with the flow conditions described in Table 1. At the pipe outlet, a relative average static pressure of zero was specified. A total number of 86,315 grid points was generated for the horizontal cylindrical pipe geometry. An O-grid layout was adopted for the distribution of grid points at the cross-sectional plane of the cylinder with more grid points being concentrated near the wall (see Fig. 1). A denser mesh of 128,765 grid points was tested. Nevertheless, the predicted cross-sectional averaged volume fractions between the two meshes yielded only differences of 2%. The mesh of 86,315 grid points was thereby employed to obtain the predicted results presented in the next section. Reliable convergence was achieved within 1500 iterations for a convergence criterion based on the RMS (Root Mean Square) residuals of 10 −4 and a physical time scale of the fully implicit solution of 0.001 s

Sensitivity of non-drag forces on time-averaged gas void fraction
In order to better understand the effect of different non-drag forces in horizontal bubbly flow, numerical calculations were carried out for three different cases using the DQMOM. The predicted void fraction distributions were compared against experimental data of  at the dimensionless axial position L/D = 253 for j g = 0.419 m/s. For the reference case, the wall lubrication constants C w1 and C w2 were taken to have values of −0.01 and 0.05 as suggested by Antal et al. (1991). According to Ekambara et al. (2008), the lift and wall turbulent induced dispersion constants C L and C TD for horizontal pipe flow took on values of −0.2 and 0.5 respectively.
In the first case, simulations were performed by varying the lift coefficient while other non-drag forces were kept according to the reference case. The lift coefficient based on the formulation by Tomiyama (1998) has been prevalently utilized for vertical bubbly flow such has been demonstrated in Cheung et al. (2007a,b). The lift coefficient relationship takes the form: Eo < 4 where Eo is the Eotvos number while the modified Eotvos number Eo d is defined by Eo d = (g( l − g )D 2 H )/ in which the maximum bubble horizontal dimension can be evaluated through the empirical correlation of Wellek et al. (1966): D H = D s (1 + 0.163Eo 0.757 ) 1/3 . The above relationship allows positive and negative lift coefficients depending on the bubble size and also accounts for the effects of bubble deformation and asymmetric wake of the bubble. The predicted volume fraction profile for a constant C L = −0.2 showed a maximum peak of 0.4 in the vicinity of the top wall and subsequently dropped below this peak at the top wall, which compared rather well with the measured profile. In contrast, the predicted volume fraction profile via the Tomiyama's correlation did not show the same volume fraction profile near the top wall and yielded a maximum value of about 0.9 at the top wall.
In the second case, simulations were carried out by varying the wall lubrication constants while other non-drag forces were kept according to the reference case. According to Krepper et al. (2005), the wall lubrication constants C w1 and C w2 have been adjusted to-0.0064 and 0.016 for vertical bubbly flow. This result showed that the volume fraction dropping only marginally below the maximum peak when compared to the reference case. Here again, better agreement was achieved through adopting the wall lubrication constants C w1 and C w2 of −0.01 and 0.05.
In the third case, simulations were carried out by applying different turbulent induced dispersion forces while other non-drag forces were kept according to the reference case. For the turbulence induced dispersion force based on Burns et al. (2004), the constant C TD was assmued to be unity. It can be observed from Fig. 2(c) that the gas bubbles were found to be well dispersed for the turbulent induced dispersion force based on Lopez de Bertodano (1998). The turbulence induced dispersion force based on Burns et al. (2004) over predicted the maximum peak of the volume fraction by a considerable margin.
From these simulations, the lift coefficient based on Tomiyama's correlation, wall lubrication constants according to Krepper et al. (2005) and Favre-averaged turbulent induced force by Burns et al. (2004), typically applied for simulating vertical bubbly flow, did not show satisfactory agreement with the measured profile. However, the predicted profile by the reference case was found to give good match with the measured profile.

Time-averaged gas void fraction
The predicted void fraction distribution of horizontal bubbly flow comparing against experimental data of  for j g = 0.213 m/s, j g = 0.419 m/s and j g = 0.788 m/s at the dimensionless axial position L/D = 253 are shown in Figs. 3-5. Satisfactory agreement was achieved between the measured and predicted radial profiles at Â = 0 • . A peak persisted for the gas void fraction in the vicinity of the upper wall of the pipe for all the different superficial gas velocities. Physically, this can be explained by the upward migration of gas bubbles due to the upward buoyant force balancing with the downward wall lubrication force to prevent the gas bubbles from collapsing at the upper wall. In contrast to vertical bubbly flow, the movement of bubble towards the wall is generally caused by the balance between the opposing lift and wall lubrication forces. It can also be seen that the peak gas void fraction value increased from the superficial gas velocities from j g = 0.213 m/s to j g = 0.788 m/s to a level of about 0.6, which corresponded to the maximum packing condition of spherical solid particles. At this level, the maximum allowable void fraction has been exceeded indicating that gas bubbles have attained their saturation limit at the upper wall. At Â = 90 • , the experimental radial profile of void fraction was shown to develop gradually from a saddle-type profile to a more parabolic profile with a single maximum at the pipe centre as the gas superficial velocity increased. Although the numerical prediction revealed similar increasing trend, the nature of a parabolic profile at j g = 0.788 m/s was not successfully captured by the current model. This discrepancy could be probably due to the insufficiency of the turbulence-induced force to dramatically push the bubbles away from the pipe wall or the specific requirement to further add a wall reaction force such as suggested by Tselishcheva et al. (2010) to counterbal-ance the buoyant force within horizontal bubbly flow; this wall reaction force remains to be fully tested and deserves separate thorough investigation in the future. In Fig. 4, the effect of buoyancy provoked the migration of gas bubbles towards the  upper wall of the pipe was further evident with the void fraction distribution becoming highly asymmetrical in the pipe cross-section. Both MUSIG and DQMOM gave similar predictions albeit lower peak void fraction was obtained for the latter. The internal flow structure has a general similarity irrespective of the superficial gas velocities.

Time-averaged liquid velocity
The  Reasonably good agreement was achieved between the measured and predicted radial profiles at Â = 0 • and Â = 90 • except for j g = 0.788 m/s. For this particular flow condition, the prevailing asymmetrical distribution was evidently observed which was attributed to the larger void fraction of gas in the upper wall of the pipe such as seen in Fig. 5. Based on the crosssection liquid velocity contours of the MUSIG and DQMOM predictions illustrated in Fig. 7, the increasing asymmetry of the liquid velocity distribution was also clearly observed with increasing gas superficial velocities. Of course, such observa-  tion would not be present in single liquid phase flow where the liquid velocity will move equally in the upper and bottom walls of the pipe, exhibiting a perfect axi-symmetry. On another note, the liquid velocities in bubbly flow system are, in general, greater than those in single phase flow system under the same flow conditions due to the inertial force acting between the gas and liquid phases. Since the liquid phase flow occupied a dominant position in the bottom section of the pipe, an interesting feature of the radial velocity profiles has been the close resemblance to a fully developed turbulent pipe flow irrespective of the superficial gas velocities. Two significant observations can be made between the void fraction and liquid velocity distributions.
Firstly, although the experimentally measured void fraction profiles demonstrated large changes along the radial direction, the changes in the velocity profile shape were comparatively very small for different Â. No apparent peaks were observed in the liquid velocity profiles when compared to the void fraction profiles which showed a peak near the upper wall. According to Beattie (1972), there was no evidence to suggest a strong proportionate correspondence between the void fraction and velocity profiles. Similar phenomena have also been observed in the vertical bubbly flow performed by Hibiki and Ishii (2001) where the "wall peak" void fraction profile essentially has no influence on the velocity distribution.
Secondly, gas bubbles in vertical bubbly flow have a tendency to accelerate along the axial direction, which they are driven by the strong buoyant force. However, the buoyant force in horizontal bubbly flow, which now acts normal to the fluid flow, exerts a lesser contribution in pushing the gas bubbles to move along the axial direction than that in vertical bubbly flow. According to  liquid velocities have been found to be only slightly greater than the velocities of the gas bubbles. The gas bubbles were accelerated by liquid inertia in a very short distance after injection but downstream of the bubbly flow, the local gas phase velocities followed closely to the local liquid phase velocities.

Time-averaged interfacial area concentration (IAC)
The predicted interfacial area concentration (IAC) distribution of horizontal bubbly flow comparing against experimental data at L/D = 253 are shown in Figs. 9-11.
As observed in , the Sauter mean bubble diameter distribution was nearly uniform for any given flow condition. It was therefore not surprising that the IAC followed very closely to the void fraction distributions. Here again, reasonably good agreement was achieved between the measured and predicted radial profiles at Â = 0 • and Â = 90 • except for j g = 0.788 m/s where significant overprediction of the peak values were obtained, more so for MUSIG as compared to DQMOM. This was probably due to the lack of robustness of the population balance approach in predicting the bubble size being closed to the maximum packing limit which may require additional consideration of different bubble mechanistic models. Owing to the smaller Sauter mean diameter, higher than expected IAC was subsequently attained. In Fig. 10, similar to the void fraction distribution, the internal flow structure also showed a general similarity irrespective of the superficial gas velocities as seen from the cross-section contours of IAC.

Conclusion
A two-fluid model coupled with MUSIG and DQMOM has been presented in this paper to handle isothermal bubbly flow in a horizontal pipe. The evaluation of the source terms and the incorporation of appropriate set of equations for the abscissas and weights of DQMOM were implemented within ANSYS-CFX code to determine the temporal and spatial geometrical changes of the gas bubbles. In comparison to the local development of the void fraction, liquid velocity and interfacial area concentration, the DQMOM along with the two-fluid model was found to compare reasonably well with the measured values as well as those of the MUSIG model with the two-fluid model. This demonstrated the strong prospect for the possible application of the DQMOM in handling practical bubbly flow in a horizontal pipe.