On the numerical study of bubbly flow created by ventilated cavity in vertical pipe

International Journal of Multiphase Flow 37 (2011) 756–768

Contents lists available at ScienceDirect

International Journal of Multiphase Flow j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / i j m u l fl o w

On the numerical study of bubbly flow created by ventilated cavity in vertical pipe M. Xiang a,b, S.C.P. Cheung b, G.H. Yeoh c,d, W.H. Zhang a, J.Y. Tu b,⇑ a

Institute of Aerospace and Material Engineering, National University of Defense Technology, Changsha 410073, PR China School of Aerospace, Mechanical and Manufacturing Engineering, RMIT University, Victoria 3083, Australia c School of Mechanical and Manufacturing Engineering, University of New South Wales, Sydney, NSW 2052, Australia d Australian Nuclear Science and Technology Organisation (ANSTO), PMB 1, Menai, NSW 2234, Australia b

a r t i c l e

i n f o

Article history: Received 5 September 2010 Received in revised form 25 January 2011 Accepted 26 January 2011 Available online 5 February 2011 Keywords: Ventilated cavity Population balance Bubbly flow CFD

a b s t r a c t The dispersion of bubbles into a down-liquid flow in a vertical pipe is investigated. At low flow rates, the intended design of a swarm of discrete bubbles is achieved. At high flow rates, a ventilated cavity is nonetheless formed, which is attached close to the gas sparger. Behind this ventilated cavity, three different flow regimes characterize the complex bubbly flow field downstream of the down-liquid flow: vortex region with high void fraction, transitional region and pipe flow region. In this study, a numerical model that solved the entire development of the gas–liquid flow including the extended single-phase liquid region upstream to the wall-jet and recirculating-vortex zones in order to allow a more realistic determination of the boundary conditions of the down-liquid flow was adopted. Coupling with the Eulerian– Eulerian two-fluid model to solve the respective gas and liquid phases, a population balance model was also applied to predict the bubble size distribution in the wake right below the cavity base as well as further downstream in the transitional and fully-developed pipe flow regions. The numerical model was evaluated by comparing the numerical results against the data derived from theoretical, numerical and experimental approaches. Prediction of the Sauter mean bubble diameter distributions by the population balance approach at different axial locations confirmed the dominance of breakage due to the high turbulent intensity below the ventilated cavity which led to the generation of small gas bubbles at high void fraction. Further downstream, the coalescence effect dominated leading to merging of the small bubbles to form bigger bubbles. Ó 2011 Elsevier Ltd. All rights reserved.

1. Introduction The flow phenomenon of a ventilated cavity is frequently encountered in industrial devices such as U-tube fermenters and airlift bioreactors (Bacon et al., 1995). Within such systems, gas bubbles being injected through a sparger into the down flowing liquid have been found to rapidly coalesce thereby forming a ventilated cavity instead of the intended design to produce a continuous swarm of dispersed bubbles. Behind this ventilated cavity such as depicted in Fig. 1, three different flow regimes characterize the complex bubbly flow field. The first regime can be identified as the highly turbulent and vigorous vortex flow region. Herein, small bubbles that tear off from the fixed cavity are entrained by the annular wall-jet at the base of the cavity. These bubbles are then shed downstream transiting through the second regime, transitional flow region, and into the third regime known as the fullydeveloped pipe flow region. The presence of a ventilated cavity has a significant influence on both the liquid motion as well as ⇑ Corresponding author. Address: SAMME, RMIT University, Bundoora, Melbourne, Victoria 3083, Australia. Tel.: +61 3 9925 6191; fax: +61 3 9925 6108. E-mail address: [email protected] (J.Y. Tu). 0301-9322/$ - see front matter Ó 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijmultiphaseflow.2011.01.014

the gas dispersion characteristics of the flow field. More importantly from an operational standpoint, the performance of such gas–liquid contacting devices where a ventilated cavity exists is strong governed by a number of important parameters: cavity length, bubble entrainment rate, length of vortex region, bubble size and void fraction distribution. Owing to the complex flow behaviors that can be experienced in this kind of two-phase flow, numerous studies have been undertaken to better understand the basic mechanisms behind the different flow structures. Bacon et al. (1995), Su (1995), Lee et al. (1997, 1999) and Kockx et al. (2005) have established relationships for the cavity length and bubble entrainment rate correlating with various gas and liquid flow rates while Van Hout et al. (1992), Rigby et al. (1997), Thorpe et al. (2001), Delfos et al. (2001), Sotiriadis and Thorpe (2005a) and Sotiriadis et al. (2005b) have concentrated on describing the flow patterns in the bubbly flow region below the ventilated cavity. Significant efforts have also been made to develop suitable numerical models (Evans et al., 2004; Sotiriadis and Thorpe, 2005a) to predict the critical parameters for such kind of two-phase flow in order to provide better control and management of these gas–liquid contacting devices. In particular, Evans et al. (2004) have adopted a population balance model to predict the

757

M. Xiang et al. / International Journal of Multiphase Flow 37 (2011) 756–768

Fig. 1. Flow field below the ventilated cavity (photograph after Delfos et al., 2001).

bubble size distribution just in the wake below the base of the ventilated cavity, which can be divided into the wall-jet and recirculating-vortex zones. Numerical calculations were carried out separately for each zone. Considering that the main feature of such systems is governed predominantly by bubbly flow, particular emphasis is directed towards a better understanding via computational fluid dynamics investigations on the complex interactions between the gas bubbles emanating from the wall-jets of the base of the ventilated cavity and the wake of the fluid below it and, more importantly, the subsequent two-phase structures due to coalescence and breakage downstream from the wake to the transitional and fully-developed pipe flow regions. A numerical model is developed capable of not only predicting the different flow regions that can be vividly encountered but also reproducing the flow field characteristics below the base of the ventilated cavity. The most fundamental difference between the current numerical model in comparison to the model applied by Evans et al. (2004) is the removal of the assumption where the specific energy dissipation rate in the wall-jet zone has been arbitrarily prescribed as three times the average rate in the wake. In this present study, the numerical model solves the entire development of the gas–liquid flow including the extended single-phase liquid region upstream to the wall-jet and recirculating-vortex zones in order to allow a more realistic determination of the boundary conditions of the liquid flow. Coupling with the Eulerian–Eulerian two-fluid model to solve the respective gas and liquid phases, the population balance model based on the class method – MUlti-SIze-Group (MUSIG) model – is also applied to predict the bubble size distribution in the wake as well as further downstream in the transitional and fullydeveloped pipe flow regions. Bubble breakup and coalescence processes are taken into account via appropriate models. The numerical model is verified and validated against data derived from theoretical, numerical and experimental approaches exemplified in Su (1995), Thorpe et al. (2001) and Evans et al. (2004). 2. Mathematical models It is assumed that the ventilated cavity represented by the shape of a Taylor bubble remains unchanged in time. The interface between the gas and liquid phases represented by the cavity boundary of the ventilated cavity as illustrated in Fig. 2 can thus

Fig. 2. Schematic diagram of the simplified physical model.

be taken to be stationary in the flow field, which can be determined according to the theoretical model formulated by Dumitrescu (1943). At the rim of the cavity base, the sheared-off gas bubbles from the Taylor bubble due to surface instability, which are significantly enhanced by the high turbulent intensity and active eddy-bubble interaction in the wake region, assume the entrainment rate to be equivalent to the gas ventilation rate. The mechanism for shearing-off effect is reflected by the fact that the gas bubbles shearing-off from the Taylor bubble contribute the primary source for spherical bubbles being injected into the wake region below the base of the ventilated cavity. The entrained bubble size is taken to be the maximum stable bubble diameter in the mixing zone proposed by Thorpe et al. (2001). Particular attention is directed towards better comprehending the complex physical behaviors of gas bubbles within the liquid flow after shearing off from the cavity base into the wake, transitional and fully-developed pipe flow regions through solving the two-phase gas–liquid flow via the two-fluid and population balance models. 2.1. Theoretical model of ventilated cavity 2.1.1. Cavity boundary It is imperative that the cavity boundary is appropriately defined in order to accurately calculate the liquid film thickness and superficial velocity at the cavity base. Here, the Taylor bubble shape is calculated based on the formula proposed by Dumitrescu (1943):

l=R ¼

" # 8 rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  2 > > < 0:75 1  1  ðr=RÞ l=R 6 0:5 0:75 > > :

0:125 ð1ðr=RÞ2 Þ2

ð1Þ

l=R > 0:5

where l refers to the axial distance from the artificial top of the cavity where the cavity width is zero, r represents the distance from

758

M. Xiang et al. / International Journal of Multiphase Flow 37 (2011) 756–768

the pipe centerline and R is the pipe radius. In the case of a long confined cavity, if the gravity acting on the film is balanced by the shear forces, the film thickness d will reach a constant value which can be calculated from the following correlation given by Henstock and Hanratty (1979):

 2=15  1=3 d ¼ 0:4868Rel2:5 þ 7:8895  108 Re4:5 m2l =g l Rel ¼

u0 D

ð2Þ

*

where a, q and u are the void fraction, density and velocity vector of each phase. Subscripts i = l or g denote the liquid and gas phase respectively. The momentum equation can be expressed as: *

** * @ðqi ai u i Þ þ r  ðqi ai u i u i Þ ¼  ai rP0 þ ðqi  qref Þai g @t h  * i * þ r  ai lei r u i þ ðr u i ÞT þ F i

ð3Þ

ð8Þ

where u0 is the liquid velocity at the inlet, vl is the liquid kinematic viscosity and D represents the pipe diameter.

where ql is adopted as the reference density qref to calculate the buoyancy force. In ANSYS CFX (Ansys, 2007), the hydrostatic pressure gradient contributed by the reference density is absorbed by 0 the modified pressure, P . Fi represents the interfacial forces for * the interfacial momentum transfer and g is the gravity acceleration vector. It is noted that the interfacial forces appearing in the momentum equation strongly govern the distribution of the liquid and gas phases within the flow volume. Details of the different forces acting within the two-phase flow are described in the next section. The void fraction of both phases is determined from the transported volume conservation equation:

ml

2.1.2. Bubble inlet conditions at the rim of cavity base The sheared-off gas bubbles with uniformly distributed bubble size are injected from a velocity inlet boundary set at the rim of the cavity base. At steady state, the bubble entrainment rate should equal to the total rate of gas ventilation rate and bubble reentrained rate. But according to Su (1995), the re-entrained rate can be ignored when compared to the entrainment rate. Therefore, the entrainment rate at the bubble inlet can be assumed to be equivalent to the gas ventilation rate. As information on the size of the inception bubbles is limited, the entrained bubble size is taken to be the maximum stable Sauter mean bubble radius in the mixing zone given by Thorpe et al. (2001):

r max

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ! u u Wec rðR  dc Þ4 1 3 t ¼ 0:9 128ql ðKulf Þ2 ðR  dc  r max =0:6Þ2

ð4Þ

where ql is the liquid density, r refers to the surface tension coefficient and Wec is the critical Weber number for breakup. A value of 4.7 is adopted for Wec and K is a constant taken to have a value of 3.8. dc refers to the thickness of the wall-jet at the section crossing the vortex center below cavity base. It is approximated using the film thickness at the cavity base db as:

dc ¼ db þ 0:068R

ð5Þ

The liquid film velocity ulf at the cavity base is obtained through the conservation of mass according to:

ulf ¼

4Q l

"

1

p D2  ðD  2db Þ2

ð9Þ

This volume conservation equation is solved separately to obtain phasic volume distribution. To enhance the stability of the simulation, a direct coupled solver is adopted to resolve the pressure and velocities coupling issue in one coupled matrix (Ansys, 2007). The fundamental basis of this coupled solver derives from the significant development of Rhie and Chow (1983) in which they have proposed a formulation by expressing the pressure in terms of its neighboring velocity values. Using this formulation, a pressure equation can be derived and directly coupled with neighboring velocities. A system of equations comprising the momentum and * pressure is solved hence forth to yield u i and P. Details of the Rhie and Chow formulation and the pressure equation can also be found in Yeoh and Tu (2009).

# ð6Þ

where Ql is the volumetric flow rate at the liquid inlet. The cavity boundary is prescribed as free-slip walls for the liquid phase as well as for the gas phase while the other walls of the physical model are modeled as non-slip for the liquid phase and free-slip for the gas phase. 2.2. Two-fluid model The two-fluid model based on the Eulerian–Eulerian framework solves the ensemble-averaged of mass, momentum and energy whereby the liquid is considered as the continuum phase and the gas is treated as disperse phase. Two sets of governing equations are solved for each phase. Interactions between phases are effected via interfacial transfer terms for heat, mass and momentum exchange. Since there is no interfacial mass or heat transfer between the phases in the present study, the energy equation is not needed to be solved. In the absence of interfacial mass transfer, the continuity equation of the two-phases can be written as: * @ðqi ai Þ þ r  ðqi ai u i Þ ¼ 0 @t

 X 1 @ðq ai Þ * i þ r  ðqi ai u i Þ ¼ 0 @t qi i

ð7Þ

2.3. Interfacial forces The interfacial forces can be categorized into drag, lift, virtual mass, wall lubrication, turbulence dispersion and Basset force. The virtual mass force is proportional to relative phasic acceleration and the Basset force accounts for the deviation of flow pattern from a steady-state condition. In this study, the governing equations for mass, momentum and scalar quantities conservation are solved using the ‘‘false transient’’ approach in order to enhance numerical stability, which marches the system of equation towards steady state. These equations are solved in a fully implicit manner at any given time step until the steady-state condition is reached. The transient terms in the governing equations will eventually become zero or have negligible influence when the simulation achieves steady-state condition. Therefore, the virtual mass and Basset forces are negligible once the simulation is converged to steady state. The total interfacial forces for the interfacial momentum transfer in Eq. (8) are thus: wall lubrication F i ¼ F lg ¼ F drag þ F lift þ F dispersion ¼ F gl lg þ F lg lg lg

ð10Þ

In the above, Flg denotes the momentum transfer from the gas phase to the liquid phase and vice versa for Fgl.

M. Xiang et al. / International Journal of Multiphase Flow 37 (2011) 756–768

2.3.1. Drag force The interfacial momentum transfer between gas and liquid due to the drag force which is resulted from shear and form drag is modeled according to Ishii and Zuber (1979) as:

F drag ¼ lg

* * * * 1 C D aif ql u g  u l ðu g  u l Þ ¼ F drag gl 8

ð11Þ

where aif is the interfacial area density and CD is the drag coefficient which can be evaluated by correlation of several distinct Reynolds number regions for individual bubbles. The model also takes into account dense particle effects, which is suitable for the two-phase flow containing high void fraction. 2.3.2. Lift force Mainly driven by the velocity gradients in radial direction, bubbles rising in a liquid are subjected to a lateral lift force. Inclusion of the lift force is important to obtain the correct spreading of the gas phase. The force can be correlated to the relative velocity and the local liquid vorticity from Drew and Lahey (1979) as: *

*

*

lift F lift lg ¼ C L ql ðu g  u l Þðr  u l Þ ¼ F gl

ð12Þ

759

2.4. Population balance model The population balance method is adopted to predict the size distribution of the poly-dispersed bubbles. The analysis is simplified by considering bubble sizes change due to breakup and coalescence events only, i.e. nucleation, growth and/or dissolution of bubbles due to absorption (desorption) or boiling are not considered. Assuming that each bubble class travel at the same mean algebraic velocity, the number density equation based on Kumar and Ramkrishna (1996) can be expressed as: * @ni þ r  ðu g ni Þ ¼ BC þ BB  DC  DB @t

ð16Þ

where ni is the average bubble number density of the ith bubble class, the source terms BC, BB, DC and DB are, respectively, the birth rates due to coalescence and breakage and the death rate due coalescence and breakage of bubbles. In a form consistent with the variables used in two-fluid model, the transport equation in terms of size fraction becomes:

@ðqg ag fi Þ * þ r  ðqg u g ag fi Þ ¼ BC þ BB  DC  DB @t

ð17Þ

For the lift coefficient, CL, the Eötvos number dependent correlation proposed by Tomiyama (1998) is adopted. More details can be found in Cheung et al. (2007).

For the MUSIG model, which employs multiple discrete bubble size groups to represent the population balance of bubbles, the birth and death rates can be written in terms of the size fraction according to

2.3.3. Wall lubrication force This force constitutes another lateral force due to the surface tension thereby allowing the bubble to concentrate in a region close to the wall but not immediately adjacent to the wall. According to Antal et al. (1991), this force can be modeled as

BC ¼ ðqg ag Þ2

i h * * * * * * 2   a q ðu g  u l Þððu g  u l Þ n w Þn w g l * D s nw F lubrication ¼  C w1 þC w2 lg yw Ds ¼ F lubrication gl

ð13Þ

where Ds refers to the Sauter mean bubble diameter, yw is the dis* tance from the wall boundary and n w is the outward vector normal to the wall. The wall lubrication constants determined through numerical experimentation for a sphere are: Cw1 = 0.01 and Cw2 = 0.05. To avoid the emergence of attraction force, the force is set to zero for large yw, viz.,

F lubrication ¼ F lubrication ¼ 0 if yw > lg gl

C w2 Ds C w1



¼ C TD C D

v t;g ral rag rt;g

al



ag

N X

fi fj

j¼1

BB ¼ qg ag

N X

1 aðM i ; M j Þ Mj

rðM j ; Mi Þfj

DB ¼ qg ag fi

N X

¼

F dispersion gl

ð15Þ

where CTD, vt,g and rt,g represent the turbulent dispersion coefficient, turbulent kinematic viscosity for the gas phase and the turbulent Schmidt number of the gas phase, respectively. Here, the turbulent dispersion coefficient CTD = 0.2 and the turbulent Schmidt number rt,g = 0.963 are used.

ð18bÞ

ð18cÞ

rðM i ; Mj Þ

ð18dÞ

j¼1

where Mj refers to the mass of the jth bubble group. The coalescence mass matrix Xjk is defined as:

X jki ¼

8 ðM þM ÞM j i1 k > > < Mi Mi1

Miþ1 ðMj þMk Þ > > Miþ1 Mi

:

0

if Mi1 < M j þ M k < M i if Mi < M j þ M k < M iþ1

ð18eÞ

otherwise

For the merging of bubbles, the coalescence of two bubbles is assumed to occur in three stages. The first stage involves the bubbles colliding thereby trapping a small amount of liquid between them. This liquid film then drains until it reaches a critical thickness. The last stage features the rupturing of the liquid film subsequently causing the bubbles to coalesce. The collisions between bubbles may be caused by turbulence buoyancy and laminar shear. Only random collisions driven by turbulence are usually considered for bubbly flow conditions. The coalescence kernel considering turbulent collision taken from Prince and Blanch (1990) can be written as the product of the collision rate and coalescence efficiency:

aðM j ; M k Þ ¼ hjk PC ðdj ; dk Þ



ð18aÞ

j¼iþ1

ð14Þ

2.3.4. Turbulent dispersion force The turbulent dispersion force is a result of the turbulence fluctuating components of the liquid phase acting on the bubbles driving it to be transported from high concentration region to low concentration region. Formulation of turbulent dispersion proposed by Burns et al. (2004) is adopted in this work. With the advantage of eliminating higher order constitution terms, the model was re-derived based on double (Farve and ensemble) averaging technique. Details of the derivation can be found in the work by Burns et al. (2004).

F dispersion lg

DC ¼ ðqg ag Þ2

i X i 1X Mj þ Mk fj fk aðMj ; M k ÞX jki 2 j¼1 k¼1 Mj Mk

ð19Þ

where the collision rate can be expressed as the product of the collision cross-sectional area A and characteristic velocity u0 as

hjk ¼

 1=2 ½dj þ dk 2 u2tj þ u2tk 4 |fflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflffl}

p

A

u0

while the coalescence efficiency is

ð20Þ

760

M. Xiang et al. / International Journal of Multiphase Flow 37 (2011) 756–768

  tjk PC ðdj ; dk Þ ¼ exp 

ð21Þ

sjk

ðdjk =2Þ2=3

ð22Þ

ðeÞ1=3

where e refers to the turbulent kinetic energy dissipation rate for the continuous phase, while the time required for two bubbles to coalesce tjk can be estimated to be

ql ðdjk =2Þ3 16r

t jk ¼

!1=2 ln



h0 hf

 ð23Þ

The equivalent diameter djk appearing in Eqs. (22) and (23) is calculated as suggested by Chesters and Hoffman (1982):

djk ¼



2 2 þ dj dk

1 ð24Þ

where the diameters dj and dk are evaluated based on the spherical bubble assumption as

dj ¼

6 1

p qg

!1=3 Mj

6 1

dk ¼

p qg

!1=3 Mk

 exp 

!1=3 Z

e 2 di

1

fmin

r

ð1 þ fÞ2 f11=3 !

12cf 2=3 d5=3 f11=3 l i

bq e

df

ð26Þ

where f ¼ k=dj is the size ratio between an eddy and a particle in the inertial sub-range and consequently fmin ¼ kmin =dmin . The lower limit of the integration is given by

dmin ¼

6 1

p qg

!1=3

M min

lt;l ¼ ql C l

k

2

ð29Þ

e

where the turbulent kinetic energy k and its dissipation rate e in Eq. (29) of the continuous liquid phase are determined from transport equations which are straightforward extensions of the standard k–e model (Lopez de Bertodano et al., 1994). Effect of turbulence in the liquid phase on turbulence in the gas phase may be modeled by setting the viscosity to be proportional to the liquid turbulent viscosity:

leg ¼

qg e l ql l

ð30Þ

The extra bubble induced turbulent viscosity in Eq. (28) is evaluated according to the model by Sato et al. (1981):

*



lt;b ¼ ql C lb ag Ds u g  u l with C lb ¼ 0:6 *

ð31Þ

ð25Þ

For air–water systems, experiments have determined h0, initial film thickness and, hf, critical film thickness at which rupture occurs to be 1  104 m and 1  108 m respectively. The turbulent velocity ut in the inertial sub-range of isotropic turbulence by Rotta (1972) is given by ut = 1.4e1/2d1/3. For the consideration of the breakage of bubbles, the model developed by Luo and Svendsen (1996) is employed for the breakage of bubbles in turbulent dispersions. In this model, binary breakage of the bubbles is assumed and the model is based on surface energy criterion and isotropic turbulence. The breakage frequency for binary breakage of bubbles is:

rðM i ; M j Þ ¼ Cð1  ag Þ

ð28Þ

Two-equation k–e model is applied to determine the liquid turbulent viscosity:

The contact time sjk when two bubbles come together is

sjk ¼

lel ¼ ll þ lt;l þ lt;b

kmin ¼ 11:4g

ð27Þ

 1=4 in which g ¼ m3l =e where ml denotes the kinematic viscosity of the continuous phase. The constants C and b are determined respectively from fundamental consideration of drops or bubbles breakage in turbulent dispersion systems to be 0.923 and 2.0. The variable cf denotes the increase coefficient of surface area: h 2=3 cf ¼ fBV þ ð1  fBV Þ2=3  1 where fBV is the breakage fraction: Mj/Mi. 2.5. Turbulence model Turbulence model is required to close the terms of effective viscosity in the Navier–Stokes equations. In this paper, separate turbulence model is adopted for the gas and liquid phases. The effective viscosity of the liquid phase is considered as being composed of the molecular viscosity ll, the turbulent viscosity lt,l and the bubble induced turbulent viscosity lt,b:

3. Experimental details Experiments by Evans et al. (2004) and Su (1995) have been chosen to verify and validate the numerical model described in the previous section. In Evans et al. (2004), experiments were carried out on a 1.8 m long vertical pipe with an internal diameter of 0.050 m. A flow straightener was fitted at the top of the vertical pipe. A gas injection, passing through the center of the flow straightener, was located along the vertical axis of the pipe. In the experiments, the fluids used were air and water. When a steady flow of water was achieved based on a particular liquid flow rate, gas was gradually introduced and a ventilated cavity was subsequently formed at the injection pipe. Based on the range of liquid and gas flow rates, the equilibrium cavity length was measured. Using the recorded images of the bubbly flow taken at the downstream exit point of the wake, volume-based spherical equivalent bubble diameters were calculated. The void fraction was obtained from pressure transducer readings. Based on a corresponding experimental condition in Evans et al. (2004), numerical predictions by the current model have been compared against the average turbulent dissipation rate in the wall-jet and wake regions and the cumulative size fraction of the bubble diameter at the exit of the vortex region. In Su (1995), experiments were performed on a 3 m long vertical pipe with an internal diameter of 0.058 m. A gas injection pipe located along the vertical axis of the pipe was used to produce a large attached ventilated cavity. An elongated shape of the ventilated cavity was formed as more gas was injected into the downward flow of the liquid. Under different superficial liquid and gas velocities, the cavity length, void fraction, bubble size and velocity distributions were measured using the double-wired conductivity bubble detection probe. More details of the experimental apparatus can be found in Su (1995). Based on the experimental conditions of different liquid flow and gas ventilation rates with equilibrium cavity length measured in Su (1995), numerical predictions by the current model have been compared against the measured axial and radial distributions of the void fraction, bubble size and velocity, which focused on the ability of the model to capture the two-phase flow in the wake below the ventilated cavity, transitional and fully-developed pipe flow regions. The fluids used were also air and water.

761

M. Xiang et al. / International Journal of Multiphase Flow 37 (2011) 756–768

4. Numerical details The generic CFD code ANSYS CFX11 was utilized to investigate the liquid flow around the ventilated cavity and the two-phase flow behavior downstream of the cavity. Numerical calculations were performed on a 45° radial sector of the pipe with symmetry boundary conditions imposed at both vertical sides of the pipe. Fig. 3 illustrates the computational mesh corresponding to the schematic geometry illustrated in Fig. 2. Taking the cavity base as the reference datum, the inlet boundary was located at a distance 0.9 m upstream while the outlet boundary was placed at a distance of 1.7 m downstream, which resulted in a total pipe length of 2.6 m. An outlay mesh of 60,000 hexahedral elements was constructed for the entire flow domain. Mesh sensitivity on the numerical predictions was performed by doubling the number of elements. Comparing the predicted results between the two meshes, small discrepancies were observed. The maximum differences between these two mesh levels were less than 5%. The transport equations of the two-fluid and population balance models were discretized by the finite volume approach. The convection terms were approximated by a high order resolution scheme while the diffusion terms were approximated by the second-order central difference scheme. Convergence was achieved within 2500 iterations when the RMS (root mean square) residual dropped

below 105. A fixed physical time scale of 0.001 s was adopted to achieve steady state solutions. For the MUSIG model, bubbles ranging from 0 mm to 10 mm diameter were equally divided into 10 size groups. 5. Results Two flow cases with the same equilibrium cavity length (i.e. as summarized in Table 1) have been selected as the basis of verifying and validating the numerical model as well as better understanding the different two-phase flow characteristics in the wake, transitional and fully-developed flow regions. In Table 1, Ql and Qg refer to the liquid flow and gas ventilation rates while Lv represents the equilibrium cavity length. 5.1. Flow characteristics in the wake region Fig. 4 clearly indicates the distinctive recirculating-vortex regions just below the base of the ventilated cavity followed by the transitional and fully-developed flow regions within the downward pipe flow. For Cases 1 and 2, the vortex regions were characterized by the high void fraction distributions and the nature of the swirling streamlines just below the base of the ventilated cavity. The gas bubbles being sheared off from the ventilated cavity were seen to follow closely with the streamlines of the vortex within the wake. Hence, the residence time of the gas bubbles in this region increased as the region was more populated by the gas phase. The thickness of the annular wall-jet defined by the outer boundary of the vortex increased with the flow developing downstream in this region. Further downstream, the gas bubbles right after the vortex were seen to migrate into the transitional region which could be identified by the spatial void fraction distribution of a ‘‘bottleneck’’ shape (see Fig. 4b). These bubbles were subsequently swept into the bubbly flow region where the gas void fraction remained almost uniform throughout the pipe section. One simplifying assumption that has been considered in the numerical model of Evans et al. (2004) was the distribution of the energy dissipation rate between the wall-jet and recirculating-vortex regions within the wake. The assumption was that the wall-jet region was defined with a specific volume where the turbulent energy dissipation rate was three times the average of that in the wake region. Table 2 shows the comparison of current numerical model predictions against the calculations by Evans et al. (2004) of similar flow condition, where Lw refers to the equilibrium wake region length, ag,pipe represents the void fraction in the fully-developed pipe flow region, and dw is the average Sauter mean bubble diameter at the exit of the wake region. Good agreement was achieved for the length of the wake region and the gas void fraction in the fully-developed region. Furthermore, the ratio of the average energy dissipation rate of the wall-jet region to the whole wake region represented by ej =ew was found to be of the similar order to the prescribed value suggested in Evans et al. (2004). The corresponding liquid recirculation streamline and different axial distributions of the turbulent energy dissipation rate in the wall-jet region for Case 1 are depicted in Fig. 5. More importantly, the current numerical model, by including the single-phase liquid flow upstream of the ventilated cavity in the calculations,

Table 1 Flow conditions for the liquid and gas phases and equilibrium cavity length.

Case 1 Case 2 Fig. 3. Visualization and mesh distribution of computational model.

Liquid flow rate, Ql (m3/s)

Gas ventilation rate, Qg (m3/s)

Equilibrium cavity length, Lv (m)

1.00e3 0.69e3

1.08e4 0.59e4

0.6 0.6

762

M. Xiang et al. / International Journal of Multiphase Flow 37 (2011) 756–768

Table 2 Comparison of current numerical model predictions against the calculations by Evans et al. (2004) of similar flow condition.

Case 1 Evans et al. (2004)

Ql (m3/ s)

Qg (m3/ s)

Lv (m)

Lw (m)

ag,pipe

dw (mm)

ej =ew

1.00e3 1.00e3

1.08e4 1.18e4

0.6 0.65

0.14 0.12

0.19 0.18

1.5 1.9

2.5 3

Fig. 5. Distribution of turbulent energy dissipation rate in the wall-jet region for Case 1.

Fig. 4. Streamlines of liquid and contours of gas void fraction distributions: (a) Case 1 and (b) Case 2.

allowed the turbulence to be implicitly determined from the development of the flow as it squeezed through the wall-jet region rather than an explicit specification which can be rather arbitrary and not applicable to a wide range of flow conditions.

The velocity fields in the wake region have a great influence on the rate of bubble entrainment and dispersion. Fig. 6 illustrates the radial velocity distributions for the axial velocity component below the ventilated cavity for both cases at the three locations (i.e. S1, S2 and S3), corresponding to the axial section across the vortex eye, the end of the wake and the fully-developed pipe flow region respectively. Positions of the three locations are also indicated in Fig. 4. Here, the experimental data by Sotiriadis and Thorpe (2005a) have been employed to assess the model predictions. In the figure, the axial velocity component uz was normalized as uz/ ulf, due to the fact that the speed of the toroidal vortex could exhibit a direct proportional relationship with the liquid superficial velocity at the base of the cavity, ulf. Good agreement was achieved between the numerical predictions and experimental measurements at all locations for both cases. As can be seen in the figure, the velocity profiles present different flow characteristics according to different flow regions. In the vortex eye section, the axial velocity component changed sign, being directed upward in the center of the pipe and downward closer to the wall. The corresponding maximum axial recirculation velocity, across the vortex eye, was ascertained to be 0.4 times of the velocity of the falling film at the base of the ventilated cavity. This is in close agreement with the value of 0.38 obtained by Thorpe et al. (2001) and

M. Xiang et al. / International Journal of Multiphase Flow 37 (2011) 756–768

Fig. 6. Axial velocity component distributions at various locations below the ventilated cavity: (a) Case 1 and (b) Case 2.

Sotiriadis and Thorpe (2005a). At the entrance to the transitional region, the axial velocity increased slightly from the pipe centerline then began to fall near the pipe wall. In the pipe flow region, the velocity was almost uniformly distributed across the whole section. In Fig. 7, it can however be seen that the radial velocity components were of significantly lower value than the corresponding axial velocity component. In the vortex eye region, the velocity initially increased from the pipe centerline until 0.5R then rapidly decreased as it crossed the x-axis at the vortex eye and became negative with a maximum absolute value at 0.9R. As the flow developed downstream, the radial velocity component attained a more uniform distribution and finally reached a zero velocity in the pipe flow region.

5.2. Location of the vortex center Fig. 8 shows the centre of the recirculating-vortex regions as well as its relative sizes for Cases 1 and 2 beneath the base of

763

Fig. 7. Radial velocity component distributions at various locations below the ventilated cavity: (a) Case 1 and Case 2.

the ventilated cavity. Although no detail measurement is available for the vortex centre, Thorpe et al. (2001) presented a numerical study attempting to simulate the vortex centre within the wake region. By assuming the annular film of 2.5 mm with 3.0 m/s liquid velocity, they reported that the vortex center is approximately located at radial position of rc = 16.3 mm and axial position of hc = 2.5 mm with the original point set at the intersection point of the cavity base and the pipe centerline. Aiming to verify our predictions with their prediction, location of the vortex center is compared based on the dimensionless ratio of rc/R, hc/D and Lw/D as tabulated in Table 3. From the table, it can be observed that the radial centre location of the recirculating-vortex cores predicted by the current model agreed well with Thorpe et al. (2001) while differences were found for the axial centre location of the recirculating-vortex cores which were probably due to different velocities of the annular film at the cavity base. For Case 1, the vortex ceased at approximately 2.8D below the base which was 0.8D further downstream than Case 2 because of the higher liquid film velocity. The predicted lengths of the vortex regions were found to be higher than those of Thorpe et al. (2001). It should be

764

M. Xiang et al. / International Journal of Multiphase Flow 37 (2011) 756–768 Table 3 Comparison of vortex center location against the numerical prediction by Thorpe et al. (2001) of similar flow condition.

Case 1 Case 2 Thorpe et al. (2001)

rc/R

hc/D

Lw/D

0.64 0.64 0.65

0.84 0.68 0.5

2.8 2.0 1.88

current numerical model prediction especially for the higher gas ventilation rate in Case 1. 5.3. Gas phase distributions In comparison to the experimental data of Su (1995), the radial void fraction distribution at different axial locations corresponding to wake, transitional and fully-developed pipe flow regions are presented in Fig. 9 for Case 1, where z refers to the axial distance from the cavity base. At 0.4D below the cavity in the vortex region,

Fig. 8. Size of recirculating-vortex regions: (a) Case 1 and (b) Case 2.

noted that the result from Thorpe et al. (2001) was obtained based on the simulation of a single-phase flow model. In the present study, the presence of gas bubbles would give rise to longer vortex region. The trend has been consistently replicated through the

Fig. 9. Radial void fraction distribution corresponding at different axial positions for Case 1. (a) Predicted void fraction distribution and (b) comparison with experimental data.

M. Xiang et al. / International Journal of Multiphase Flow 37 (2011) 756–768

the void fraction increased gradually from the centerline of the pipe reaching a maximum near the boundary of the vortex. Then the void fraction decreased sharply to zero adjacent to the wall. At 1D away from the cavity base, the void fraction distribution remained similar to the upstream location, with slightly lower values throughout the cross section. For the transitional region at the bottom of the wake displayed by 2.8D, a substantially lower void fraction compared with that close to the cavity base was predicted. The void fraction first approached maximum at the center, and then decreased gently approaching the wall. This radial distribution indicated that the bubbles were released from the vortex region mainly near the centerline, resulted from the ‘‘bottleneck’’ shape such as seen from the void fraction distribution. In the pipe flow region represented by 8D, the void fraction remained nearly uniform and then dropped off near the wall. At all positions, the void fraction at the immediate vicinity of the wall was nearly zero, indicating a narrow region free of gas bubbles. For the locations 0.4D and 8D, good agreement was achieved between the numerical predictions and experimental measurements. Note that the void fraction was under-predicted near the pipeline centre at 0.4D which was caused by the assumption of the flat boundary assumed for the base of the cavity in the current numerical model. Fig. 10 demonstrates the average void fraction at different sections along the axial direction for both cases. The axial void fraction distribution decreased gradually from the base of ventilated cavity into the vortex region. When the gas bubbles entered the transitional region, the average void fraction decreased linearly as the bubbles were shed quickly downstream into the fully-developed pipe flow region where the void fraction reached a constant downstream. The under-prediction of the void fraction for Case 2 beneath the base of ventilated cavity was due to assumption of the flat boundary imposed during the numerical calculations, which showed a more pronounced effect for a shorter wake length. Nevertheless, as the flow developed further downstream, the difference diminished with good agreement being achieved the numerical predictions and experimental measurements for both cases. It should be noted that with the same cavity length, the void fraction in the vortex region increased with the liquid flow rate due to the stronger vortex recirculation flow while in the pipe flow region the void fraction was almost the same which depended primarily on the ratio of liquid and gas flow rates.

Fig. 10. Area-averaged void fraction distribution along the different axial locations.

765

5.4. Phenomenological bubble size evolution due to breakage and coalescence The cumulative size distribution of the bubble diameter at the exit of the recirculating-vortex region for Case 1 is presented in Fig. 11. Close agreement attained between the numerical predictions and experimental measurements further reaffirmed the application of the current numerical model in aptly predicting the bubble size distribution in the wake region. The predicted radial distributions of the Sauter mean bubble diameter at various axial locations are presented in Fig. 12. In the vortex region, the bubble size decreased from the pipe centerline due to high turbulence intensity near the vortex boundary. Although high void fraction existed in this region, the dominance of the breakage process led to the generation of smaller bubbles. In the transitional region, the coalescence process became more prominent due to the

Fig. 11. Cumulative size fraction of bubble diameter at the exit of the vortex region for Case 1.

Fig. 12. Radial distribution of Sauter mean bubble diameter along the different axial sections.

766

M. Xiang et al. / International Journal of Multiphase Flow 37 (2011) 756–768

decreasing turbulent eddy dissipation rate. The resultant bubble size shed from the end of the vortex was slightly bigger than that in the vortex eye section. In the fully-developed pipe flow region, the bubble coalescence and breakage rates almost reached an equilibrium whereby the bubble diameter distribution was uniform across most of the pipe section, dropping near the pipe wall due to low void fraction. In the experiment, only the Sauter mean bubble diameter in the fully-developed pipe flow region was measured; good agreement was also found between the numerical predictions and experimental measurements. The evolution of the area-averaged size fraction for each bubble group along the axial sections is shown in Fig. 13. As depicted, in both cases, shear-off bubbles with 1.5 mm diameter were mostly disintegrated into the smaller bubble classes (i.e. bubble class 1

with 0.5 mm diameter) due to the high turbulence in the vortex region (i.e. z/D 6 2.8 for Case 1 and z/D 6 2 for Case 2). Subsequently, the coalescence effect became dominant in transition and pipe flow regions leading to merging of the small bubbles to form bigger bubbles, resulting in a higher fraction of the larger bubble classes. This observation can be also affirmed by the trend of area-averaged Sauter mean bubble diameter along the axial sections shown in Fig. 14. Comparing Cases 1 and 2, higher liquid flow rate in Case 1 brought about a stronger turbulence which enhanced the breakage of bubbles in the vortex region and gave birth to smaller bubbles in the pipe flow region. The aforementioned predicted results clearly demonstrated the dynamical changes of the bubble size due to breakage and coalescence processes in the different flow regions below the ventilated cavity. It also revealed

Fig. 13. Area-averaged size fraction of each bubble class along the different axial locations: (a) Case 1 and (b) Case 2.

M. Xiang et al. / International Journal of Multiphase Flow 37 (2011) 756–768

767

cavity imposed in the numerical model which gave an underprediction of the void fraction in the near the pipe centerline just beneath the ventilated cavity, remarkable agreement was nonetheless achieved further downstream of the pipe flow. The velocity distributions were also found to be in good agreement with experimental measurements indicating that qualitative details of the two-phase flow structure were successfully ascertained. The Sauter mean bubble diameter distributions predicted by the population balance approach based on the MUSIG model at different axial locations confirmed the dominance of breakage due to the high turbulent intensity below the ventilated cavity which led to the generation of small gas bubbles at high void fraction. This was in contrast to typical normal bubbly flow where bubbles tended to coalesce to form elongated bubbles at high void fraction. The current proposed numerical model will be able to provide a useful methodology in predicting multiphase and multi-region flow field below the ventilated cavity. The detailed flow field information obtained can lead to invaluable insights in designing and controlling of the two-phase flow systems with an existing ventilated cavity.

Acknowledgments The financial support provided by the Chinese Scholarship Council (CSC 2009611040) and the Australian Research Council (ARC Project ID DP0877743) is gratefully acknowledged.

References

Fig. 14. Area-averaged Sauter mean bubble diameter along the different axial sections: (a) Case 1 and (b) Case 2.

that the population balance approach based on the MUSIG model could be an investigation tool to gain more insight and understanding of the complex bubbly flow structure behind the ventilated cavity. 6. Conclusions Bubbly flow below the ventilated cavity can be aptly characterized by three different regions: (i) vortex region with high void fraction, (ii) transitional flow and (iii) fully-developed pipe flow. The numerical model, which comprised the theoretical formulation of the ventilated cavity coupled with the Eulerian–Eulerian twofluid and population balance models, has been found to be able to predict the turbulent two-phase flow structure as well as the bubble size distribution within the three different regions. Comparison of the model predictions were assessed against the data derived from theoretical, numerical and experimental approaches in Su (1995), Thorpe et al. (2001) and Evans et al. (2004). Despite the assumption of the flat boundary for the base of the ventilated

ANSYS, 2007. CFX-11 User Manual. ANSYS CFX. Antal, S.P., Lahey Jr., R.T., Flaherty, J.E., 1991. Analysis of phase distribution in fully developed laminar bubbly two-phase flow. Int. J. Multiph. Flow 17, 635–652. Bacon, R.P., Scott, D.M., Thorpe, R.B., 1995. Large bubbles attached to spargers in downwards two-phase flow. Int. J. Multiph. Flow 21, 949–959. Burns, A.D., Frank, T., Hamill, I., Shi, J.-M., 2004. The favre averaged drag model for turbulence dispersion in Eulerian multi-phase flows. In: Proceedings of the Fifth International Conference on Multiphase Flow, ICMF, Yokohama, Japan. Chesters, A.K., Hoffman, G., 1982. Bubble coalescence in pure liquids. Appl. Sci. Res. 38, 353–361. Cheung, S.C.P., Yeoh, G.H., Tu, J.Y., 2007. On the numerical study of isothermal vertical bubbly flow using two population balance approaches. Chem. Eng. Sci. 62, 4659–4674. Delfos, R., Wisse, C.J., Oliemans, R.V.A., 2001. Measurement of air-entrainment from a stationary Taylor bubble in a vertical tube. Int. J. Multiph. Flow 27, 1769– 1787. Drew, D.A., Lahey Jr., R.T., 1979. Application of general constitutive principles to the derivation of multidimensional two-phase flow equation. Int. J. Multiph. Flow 5, 243–264. Dumitrescu, D.T., 1943. Stromung an einer Luftblase im Senkrechten Rohr. Z. angew. Math. Mech. 23, 139. Evans, G.M., Machniewski, P.M., Bin, A.K., 2004. Bubble size distribution and void fraction in the wake region below a ventilated gas cavity in downward pipe flow. Chem. Eng. Res. Des. 82 (A9), 1095–1104. Henstock, W.H., Hanratty, T.J., 1979. Gas absorption by a liquid layer flowing on the wall of a pipe. AIChE J. 25, 122–131. Ishii, M., Zuber, N., 1979. Drag coefficient and relative velocity in bubbly, droplet or particulate flows. AIChE J. 25, 843–855. Kockx, J.P., Nieuwstadt, F.T.M., Oliemans, R.V.A., 2005. Gas entrainment by a liquid film falling around a stationary Taylor bubble in a vertical tube. Int. J. Multiph. Flow 31, 1–24. Kumar, S., Ramkrishna, D., 1996. On the solution of population balance equations by discretisation – I. A fixed pivot technique. Chem. Eng. Sci. 51, 1311–1332. Lee, Y.H., Scott, D.M., Thorpe, R.B., 1997. The scale-up of large bubbles attached to spargers in downward two-phase flow. Chem. Eng. Sci. 52, 3797–3809. Lee, Y.H., Scott, D.M., Thorpe, R.B., 1999. The scale-up of large bubbles attached to spargers in downward two-phase flow. Chem. Eng. Sci. 54, 3839–3898. Lopez de Bertodano, M., Lahey Jr., R.T., Jones, O.C., 1994. Development of a k–e model for bubbly two-phase flow. J. Fluids Eng. 116, 128–134. Luo, H., Svendsen, H., 1996. Theoretical model for drop and bubble breakage in turbulent dispersions. AIChE J. 42, 1225–1233. Prince, M.J., Blanch, H.W., 1990. Bubble coalescence and breakage in air sparged bubble columns. AIChE J. 36, 1485–1499. Rhie, C.M., Chow, W.L., 1983. Numerical study of the turbulent-flow past an airfoil with trailing edge separation. AIAA J. 21, 1525–1532. Rigby, G.D., Evans, G.M., Jameson, G.J., 1997. Bubble breakup from ventilated cavities in multiphase reactors. Chem. Eng. Sci. 52, 3677–3684. Rotta, J.C., 1972. Turbulente Stromungen. B.G. Teubner, Stuttgart.

768

M. Xiang et al. / International Journal of Multiphase Flow 37 (2011) 756–768

Sato, Y., Sadatomi, M., Sekoguchi, K., 1981. Momentum and heat transfer in twophase bubbly flow – I. Int. J. Multiph. Flow 7, 167–178. Sotiriadis, A.A., Thorpe, R.B., 2005a. Liquid re-circulation in turbulent vertical pipe flow behind a cylindrical bluffbody and a ventilated cavity attached to a sparger. Chem. Eng. Sci. 60, 981–994. Sotiriadis, A.A., Thorpe, R.B., Smith, J.M., 2005b. Bubble size and mass transfer characteristics of sparged downwards two-phase flow. Chem. Eng. Sci. 60, 5917–5929. Su, C., 1995. Gas Entrainment at the Bottom of a Taylor Bubble in Vertical Gas– Liquid Slug Flows. Ph.D Thesis, University of Houston.

Thorpe, R.B., Evans, G.M., Zhang, K., 2001. Liquid recirculation and bubble breakup beneath ventilated gas cavities. Chem. Eng. Sci. 56, 6399–6409. Tomiyama, A., 1998. Struggle with computational bubble dynamics. In: Proceeding of the Third International Conference on Multiphase Flow. Lyon, France. Van Hout, R., Shemer, L., Barnea, D., 1992. Spatial distribution of void fraction within a liquid slug and some other related slug parameters. Int. J. Multiph. Flow 18, 831–845. Yeoh, G.H., Tu, J.Y., 2009. Computational Techniques for Multiphase Flows. Butterworth-Heinemann, New York.