SCIENCE CHINA Physics, Mechanics & Astronomy • Article •

February 2012 Vol.55 No.2: 297–304 doi: 10.1007/s11433-011-4614-4

Numerical investigation of turbulent bubbly wakes created by the ventilated partial cavity XIANG Min1*, ZHANG WeiHua1, CHEUNG S. C. P.2 & TU JiYuan2 1

Institute of Aerospace and Material Engineering, National University of Defense Technology, Changsha 410073, China; 2 School of Aerospace, Mechanical and Manufacturing Engineering, RMIT University, Victoria 3083, Australia Received March 6, 2011; accepted June 29, 2011; published online January 9, 2012

This paper presents a numerical study on the turbulent bubbly wakes created by the ventilated partial cavity. A semi-empirical approach is introduced to model the discrete interface of the ventilated cavity and its complex gas leakage rate induced by the local turbulent shear stress. Based on the Eulerian–Eulerian two-fluid modeling framework, a population balance approach based on MUltiple-SIze-Group (MUSIG) model is incorporated to simulate the size evolution of the sheared off microbubbles and its complex interactions with the two-phase flow structure in the wake region. Numerical predictions at various axial locations downstream of the test body were in satisfactory agreement with the experimental measurements. The captured bubbly wake structure illustrates that the bubbles may disperse as a twin-vortex tube driven by gravity effect. The predicted Sauter mean bubble diameter has confirmed the dominance of the coleascense process in the axial direction. As the bubbles develop downstream, the coleascense and breakup rate gradually reach balance, resulting in the stable bubble diameter. A close examination of the flow structures, gas void fraction distributions and the bubble size evolution provides valuable insights into the complex physical phenomenon induced by ventilated cavity. ventilated cavity, population balance, bubbly wake, computational fluid dynamics PACS number(s): 47.55.Ca, 47.55.df, 47.85.Dh, 47.27.wb Citation:

Xiang M, Zhang W H, Cheung S.C.P, et al. Numerical investigation of turbulent bubbly wakes created by the ventilated partial cavity. Sci China-Phys Mech Astron, 2012, 55: 297304, doi: 10.1007/s11433-011-4614-4

Ventilated partial cavity, referring to the ventilated cavity terminating on the vehicle body with a re-entrained jet, is frequently observed in the transient launch stage or for vehicle manoeuvring by exposing part of the body to contact with water [1]. A complex multiphase flow field is produced by the ventilated partial cavity as shown in Figure 1. A clear continuous cavity is formed at the front part of the submerged body; at the cavity tail, due to high turbulent shear stress, the continuous cavity is broken up by the re-entrained jet. Gas pockets are sheared off at the cavity base, giving birth to numerous small bubbles under the action of high turbulence. The generated bubbles are then

*Corresponding author (email: [email protected]) © Science China Press and Springer-Verlag Berlin Heidelberg 2012

dispersed downstream, creating bubbly wakes behind the test body. The interaction between a ventilated cavity and its turbulent bubbly wake has to be fully understood for better design and control of high-speed cavity-running bodies. This paper aims at gaining a better understanding of the complex interactions between the sheared off microbubbles

Figure 1 Schematic diagram of multiphase flow field created by the ventilated partial cavity. phys.scichina.com

www.springerlink.com

298

Xiang M, et al.

Sci China-Phys Mech Astron

and the two-phase flow structure in the wake region created by the ventilated cavity in the re-entrant jet regime using a numerical approach. The theme of this paper therefore focuses on establishing a complete numerical scheme with a particular emphasis on predicting the gas bubble dispersion effect, but incorporating essential consideration of the cavity free-surface interface. Based on our previous studies [2,3], a numerical scheme has been proposed to model the entire development of the ventilated partial cavity situation, starting from the single-phase liquid region upstream to the discrete bubble dispersion region downstream. In conjunction with the Eulerian–Eulerian two-fluid model, which solves the phasic distribution explicitly with fundamental considerations of interfacial momentum transfer, a population balance approach based on the MUltiple-SIze-Group (MUSIG) model [4–7] is adopted to capture the evolution of the bubble size distribution caused by bubble coalescence and breakage events. Numerical predictions are validated against the recent experimental work carried out in the St. Anthony Falls Laboratory at the University of Minnesota [8,9]. The numerical results will provide detailed flow field information for better understanding the underlying phenomena embedded in ventilated partial cavity.

1 Physical and mathematical modeling 1.1 Semi-empirical approach for ventilated partial cavity Theoretically speaking, formation of the ventilated partial cavity is characterized by three main dynamical multiphase physical phenomena: (i) continuous cavity formation; (ii) gas leakage caused by turbulent shear stress at the cavity base and (iii) dispersion of microbubbles and its turbulence modulation to the liquid phase. This paper, to minimize the computational burden, models the complex cavity formation and gas leakage phenomena via the mechanistic model based on fundamental theories and experimental observations. The semi-empirical formula proposed by Savchenko [10] is adopted to determine the cavity profile: cd 2 x  4x 1  R  Dn  2  c 1 1 Dn 2 2  Dn In  In  c c 

2

1/ 2

    

,

(1)

here, c is the cavitation number, cd is the drag coefficient of the cavitator and x represents the axial distance from the artificial original point of the cavity curve where the cavity diameter equals zero. R and Dn refer to the cavity radius and cavitator diameter respectively. At the cavity tail, the gas traveling along the boundary layer of the high velocity cavity wall is entrained by

February (2012) Vol. 55 No. 2

re-entrained jet [11,12], resulting in continuous gas leakage from the cavity. The air entrainment rate is given in accordance with the model proposed by Ma et al. [13]: Q g  2kRlU c eff ,

(2)

where Uc is the flow velocity at the cavity boundary. Rl is the cavity radius at the re-entrainment point and eff is the effective thickness of the cavity boundary layer. k is an empirical constant calibrated against experimental conditions. With application of the velocity profile over a flat plat in the cavity boundary layer [14], the variation of air entrainment rate along with the re-entrianment position can be obtained. Adversely, with the known air entrainment rate which equals the ventilation rate under a steady state, the re-entrained position can be estimated by the air entrainment model. Due to the limited information on the size of the inception bubbles, the entrained bubble size is taken to be the maximum stable bubble diameter in the mixing zone given by Thorpe [14]:  We  R 4  1 rmax  0.9 3  c l 2  , 2 U 128  l c   Rl  rmax / 0.6  

(3)

where Wec is the critical Weber number for breakage with a value of 4.7 adopted.  refers to the surface tension coefficient. 1.2

Two-fluid model

The two-fluid model based on the Eulerian–Eulerian framework solves the ensemble-average of mass, momentum and energy whereby the liquid is considered as the continuum phase and the gas is treated as the dispersed phase. Two sets of governing equations are solved for each phase. Interactions between phases are affected 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 does not need to be solved. In the absence of interfacial mass transfer, the continuity equation of the two-phases can be written as: ( i i )    ( i i ui )  0, t

(4)

where ,  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: ( i i ui )    ( i i ui ui )   i P  ( i   ref ) i g t     i ie ui  (ui ) T   Fi , (5)





Xiang M, et al.

Sci China-Phys Mech Astron

where l is adopted as the reference density ref, Fi represents the interfacial forces for the interfacial momentum transfer and g is the gravity acceleration vector. The interfacial momentum transfer between gas and liquid due to the drag force is modeled according to Ishii and Zuber [16] as: 1 Flg drag  CD aif  l ug  ul  ug  ul    Fgl drag , 8

(6)

where 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 different bubble shape regimes, such as spherical particle regime, distorted particle regime and spherical cap regime. In addition to bubble Reynolds number and deformation, locally high gas void fraction can influence the drag dramatically [17,18]. Here, a “cluster” drag model originated from Johansen and Boysan [19] has been adapted to take into account the high void fraction effect on the bubble drag force: CD  CD 0 (1  [min(0.926,  g )]2 /3 ),

1.3

   g  g fi  t

   (  g ug  g fi )  BC  BB  DC  DB .

1.4

Turbulence model

In this paper, the 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 l , the turbulent viscosity t,l and the bubble induced turbulent viscosity t,b:

le   l  t ,l  t ,b .

ni    (ug ni )  BC  BB  DC  DB , t

t , l   l C

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 to coalescence and breakage of bubbles. In a form consistent with the variables used in the twofluid model, the transport equation in terms of size fraction becomes:

k2



,

(11)

where the turbulent kinetic energy k and its dissipation rate

 in eq. (11) of the continuous liquid phase are determined from transport equations which are straightforward extensions of the standard k- model [23]. 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:

g e 

(8)

(10)

The two-equation k- model is applied to determine the liquid turbulent viscosity:

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 travels at the same mean algebraic velocity, the number density equation based on Kumar and Ramkrishna [20] can be expressed as:

(9)

In this study, the coalescence model considering turbulent collision is adopted from Prince and Blanch [20]. The model developed by Luo and Svendsen [22] is employed considering the breakage of bubbles in turbulent dispersions. For more details regarding the model please refer to ref. [3].

(7)

where CD0 is the original drag coefficient calculated according to Ishii and Zuber [16] and the minimum function is provided to ensure that the corrected drag coefficient maintains minimum 5% of its original value. For details of the other interfacial forces including the lift, lubrication and turbulence dispersion force, please refer to ref. [3].

299

February (2012) Vol. 55 No. 2

g e  . l l

(12)

The extra bubble induced turbulent viscosity in eq. (10) is evaluated according to the model by Sato et al. [24]:

t ,b  l Cb g Ds ug  ul , with Cb  0.6 .

(13)

2 Experimental and numerical details Numerical predictions are validated against the experimental measurements in the high-speed water tunnel with a test section of 0.19 m (W)×0.19 m (H)×1 m (L) at St. Anthony Falls Laboratory of the University of Minnesota [8]. A test body with 10mm cavitator installed at the front was placed at the middle of the test section as shown in Figure 2. The clear cavity boundary in the front and the bubbly flow downstream were observed in the snapshot shown in Figure 3. Measurements were taken in the bubbly wake region downstream of the test body.

300

Xiang M, et al.

Sci China-Phys Mech Astron

3 3.1

Figure 2

Schematic display of the test body and boundary condition.

Figure 3 Snapshot of ventilated cavitation and gas leakage behaviour by Schauer [8].

Three flow conditions as shown in Table 1 are selected for investigation. The estimated re-entrainment position is at the point where a maximum diameter of the cavity occurs, which is in accordance with the experimental observation. As discussed previously, the free surface profile of the continuous ventilated cavity is evaluated in accordance with eq. (1) and is assumed to remain steady and stationary in the computational domain as shown also in Figure 2. As the experimental water tunnel and test body are symmetric, numerical simulations are only performed on a half portion of the water tunnel with symmetry boundary condition imposed at the domain side. The air entrainment rate is specified through an inlet boundary set at the circumference of the cavity base. The cavity interface is set as free-slip boundary while the test body is modeled as no-slip for both the liquid and the gas phases. 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. A directly coupled solver is adopted to solve the pressure and velocities in one coupled matrix. The fundamental basis of this coupled solver derives from the significant development of Rhie and Chow [25] where they have proposed a formulation by expressing the pressure in terms of its neighboring velocity values. A fixed physical time scale of 0.0005 s was adopted to achieve steady state solutions.

February (2012) Vol. 55 No. 2

Results and discussion Gas phase distribution

The predicted gas void fraction and liquid superficial velocity streamline distribution of different cases are shown in Figure 4. As depicted by the swirling streamlines, a toroidal vortex structure was observed right behind the test body for Case 1. This re-circulation motion held up a portion of the leakage bubbles forming a high gas void fraction area within the wake region. Low pressure induced by the vortex will contribute to extra drag force to the body. In the radial sections downstream of the vortex region, a “crescent” shape region with high void fraction was observed around the central area of the section. As the flow developed downstream, the “crescent” region gradually shifted up due to the gravity effect along with the decrease of the peak void fraction. This void fraction distribution indicated that the bubbles in the wake were dispersing as a twin-vortex tube [11]. For Case 2 with a lower cavitation number, the larger cavity delayed the flow reattachment length, which incidentally aligned the main liquid streamline with the profile of the test body. The recirculation behind the test body was therefore avoided. The flow field in the bubbly wake remained almost axi-symmetric and the high void fraction region in the radial section presented itself as a “disk”. This phenomenon provides insights into the transition mechanism between re-entrained jet and twin-vortex for ventilated supercavity. For a re-entrained jet supercavity, vortex stays at the cavity tail. As the cavitation number decreases, gravity effect gradually dominates, resulting in the asymmetry of the vortex region and the dispersion of the floated bubbles as a twin-vortex tube. As the number of the sheared-off bubbles gets larger and the gravity effect becomes strong enough to separate the twin-vortex, two vortex tubes appear

Table 1 Flow parameters of the selected flow conditions for validations Fr

Qg

c

Case 1

29

0.08

0.25

Case 2

29

0.16

0.15

Case 3

20

0.18

0.15

Figure 4 Gas void fraction and water velocity streamline in the bubbly wake. (a) Case 1; (b) Case 2.

Xiang M, et al.

Sci China-Phys Mech Astron

at the cavity tail. Figure 5 illustrates the comparison between predicted and measured radial gas void fraction distributions at different axial distances away from the aft of the test body. It was observed that the void fraction gradually reduced along the axial direction. Nonetheless, a steep void fraction gradient was found in the radial direction indicating the rapid dispersion of bubbles due to the combining lift and buoyancy forces. The qualitative agreement between the predicted and measured results indicates that the proposed numerical scheme is capable of modeling the essential physical phenomenon embedded in the partially ventilated cavitation problem. Figure 6 presents the comparison of void fraction distribution under different cavitation numbers. For Case 1, the peak of void fraction profiles gradually shifted upward away from the test body centerline driven by the gravity effect. Compared with Case 2 which contains a

Figure 5 Comparison of the predicted void fraction in the bubbly wake with experimental data (Case 3).

Figure 6 Comparison of the predicted void fraction in the bubbly wake under different cavitation numbers (Case 1 and Case 2).

February (2012) Vol. 55 No. 2

301

higher ventilation rate, a higher void fraction was observed for Case 1 due to the larger vortex behind the test body, which made bubbles concentrate. In comparison with Case 3, the higher velocity for Case 2 contributed to stronger dispersion of bubbles. 3.2

Analysis on the bubbly wake structure

Figure 7 presents water velocity distribution in the bubbly wake. The velocity magnitude within the wake gradually reached its freestream value as the downstream distance increased. A higher water velocity for Case 2 was observed due to the stronger dispersion mentioned above. Figure 8 displays the turbulence intensity distribution in the bubbly wake, which presents similar variation characteristics as the radial void fraction. Both the peak of the velocity and the turbulence intensity profiles were shifted upward under the action of gravity. However, the shifted distance for the tur-

Figure 7 Turbulence intensity distribution in the bubbly wake.

Figure 8

Displacement thickness variation in the bubbly wake.

302

Xiang M, et al.

Sci China-Phys Mech Astron

bulence intensity was higher than the other profiles. This conclusion was in accordance with the experimental observations in the bubbly wake behind a cavitating hydrofoil [9]. Figure 9 gives the variation of the displacement thickness * in the bubbly wake. Here  represents the momentum thickness. The dimensionless parameter in the figure increased linearly along with the dimensionless distance from the test body, indicating the wake for both cavitating and non-cavitating flows grew similarly as the turbulent axisymmetric wake with a high Reynolds number [26]. Figure 10 presents the mean bubble velocity profile scaled with the centerline velocity deficit u0. It is concluded that the self-similarity of the wake shape exists in the bubbly wake proved by the reasonable collapse of the mean velocity profiles at different locations. 3.3

February (2012) Vol. 55 No. 2

tion of bubble size that depended on both void fraction and turbulence dissipation rate. High dispersion of turbulence dominated in the axial direction, leading to the coleascense of bubbles. The same bubble size evolution trend was also observed in the experimental results [9] where comparable void fraction and bubble size existed in the bubbly wake. Figure 12 gives the bubble size distribution profiles in the bubbly wake. The bubble size decreased sharply in the radial direction due to the high dispersion of bubbles. In the axial direction, the maximum bubble size increased and the peak point shifts upward gradually. Smaller bubbles were observed at the same axial position for Case 1, which resulted from the high turbulence and low void fraction shown in Figure 6. For a better understanding of the bubble breakup and coalescence processes, Figure 13 presents the size fraction

Bubble size evolution

Figure 11 displays the bubble size distribution in the bubbly wake. As the void fraction distribution shown in Figure 4, large bubbles existed around the central “disk” and “crescent” area in the radial sections. However, as the axial distance increased, bubbles got larger, behaving in a way that was contrary to the void fraction. This was due to the evolu-

Figure 11 Bubble size distribution contours in the bubbly wake. (a) Case 1; (b) Case 2.

Figure 9

Figure 10

Bubble velocity profile in the bubbly wake (Case 2).

Bubble size distribution profile in the bubbly wake (Case 2).

Figure 12 Size fraction of each bubble class along the centerline. (a) Case 1; (b) Case 2.

Xiang M, et al.

Sci China-Phys Mech Astron

Figure 13 (Color online) Size fraction of each bubble class along the centerline. (a) Case 1; (b) Case 2.

of each bubble class along the centerline. As depicted, for both cases, in the region close to the test body (x/Dn<5), most size fraction was taken up by low size bubble groups (Ds<400m). As the bubbles dispersed downstream (x/Dn< 15), coalescense dominated due to the dramatically weakening of turbulence, leading to the merging of small bubble groups, and resulting in a higher fraction of the larger bubble classes. Subsequently, the variation magnitude of the group size fraction became smaller. In the region further downstream (x/Dn>20), the coalescense and breakup rate almost reached an equilibrium whereby the bubble size fraction changed slightly. The aforementioned predicted results clearly demonstrated the dynamical changes of the bubble size due to breakage and coalescence processes in the bubbly wake created by the ventilated partial cavity. It also revealed that the population balance approach based on the MUSIG model could be an investigation tool to gain more insight into and understanding of the complex bubbly flow structure in the bubbly wake.

4

Conclusions

A complete numerical scheme for investigating the bubbly

February (2012) Vol. 55 No. 2

303

wake induced by the ventilated partial cavity was presented. A semi-empirical approach was proposed to model the discrete interface profile and the gas leakage rate of the ventilated cavity. The Eulerian-Eulerian two-fluid model, coupled with the population balance approach, was introduced to handle the two-phase flow field created by the ventilated partial cavity. Numerical predictions at various axial locations downstream of the test body were in satisfactory agreement with the experimental measurements. For the lower cavitation number condition, the flow re-attachment point was delayed by the large cavity and successfully minimized the vortex volume behind the test body, resulting in an axisymmetric bubbly wake region. In contrast, for the higher cavitation number condition, a vortex region with a high void fraction was created by the earlier flow re-attachment. The bubbles were dispersed as a “twin-vortex” tube, creating a “crescent” shape region with a high void fraction in the radial section which gradually shifted upward driven by the gravity effect. Through the velocity analysis, it was found that the bubbly wake grew similarly as the turbulent axisymmetric wake with a high Reynolds number and the self-similarity of the wake shape existed. The predicted Sauter mean bubble diameter confirmed the dominance of the breakup process in the radial direction. However, coalescence dominated in the axial direction near the test body, which generated more large bubbles. As the bubbles developed downstream, the coleascense and breakup rate gradually reached balance, resulting in a stable bubble diameter. This study provides detailed flow field information for complex multiphase ventilated partial cavity situations, which is assessable from the existing measuring technique. Numerical predictions and analysis of this study help to condense more in-depth understanding of and invaluable insights into the problem, which is likely to contribute to higher efficiency for cavitation vehicles in the future. This work was supported by the Chinese Council Scholarship (Grant No. 2009611040) and the Australian Research Council (Grant No. DP0877743). 1

2

3

4

5

6

Abraham N V. High-speed bodies in partially cavitating axisymmetric flow. In: Fifth International Symposium on Cavitation. Osaka, Japan, 2005 Mohanarangam K, Cheung, S C P, et al. Numerical simulation of microbubble drag reduction using population balance model. Ocean Eng, 2009, 36(11): 863–872 Xiang M, Cheung S C P, Yeoh G H, et al. On the numerical study of bubbly flow created by ventilated cavity in vertical pipe. Int J Multiphase Flow, 2011, 37(7): 756–768 Olmos E, Gentric C, Vial C H, et al. Numerical simulation of multiphase flow in bubble column. Influence of bubble coalescence and break-up. Chem Eng Sci, 2001, 56(21-22): 6359–6365 Yeoh G H, Tu J Y. Thermal-hydrodynamic modelling of bubbly flows with heat and mass transfer. Am Inst Chem Eng J, 2005, 51(1): 8–27 Frank T, Zwart P J, Shi J, et al. Inhomogeneous MUSIG model––a population balance approach for polydispersed bubbly flow. In:

304

7

8

9

10

11

12 13

14 15

16

Xiang M, et al.

Sci China-Phys Mech Astron

Proceeding of International Conference for Nuclear Energy for New Europe. Bled, Slovenia, 2005 Cheung S C P, Yeoh G H, Tu J Y. On the numerical study of isothermal vertical bubbly flow using two population balance approaches. Chem Eng Sci, 2007, 62(3): 4659–4674 Schauer T J. An experimental study of a ventilated supercavitating vehicle. Dissertation for the Master Degree. Minneapolis, MN: University of Minnesota, 2003 Qin Q, Wang H, Arndt R E A. An integrated experimental/numerical study of the bubbly wake behind a cavitating hydrofoil. In: 2005 ASME Fluids Engineering Division Summer Meeting and Exhibition. Houston, TX, USA, 2005 Savchenko Y N. Investigation of high-speed supercavitating underwater motion of bodies. In: High-speed motion in water. AGARD Report 827, 20-1-20-12, 1998 Kinzel M P. Computational techniques and analysis of cavitatingfluid flows. Dissertation for the Doctoral Degree. Pennsylvanis: The Pennsylvanis State University, 2008 Fernandes R C, Semiat R. Hydrodynamic model for gas-liquid slug flow in vertical tubes. Am Inst Chem Eng J , 1983, 29(6): 981–989 Ma J, Oberai A A, Drew D A, et al. A quantitative sub-grid air entrainment model for bubbly flows––plunging jets. Comput Fluids, 2010, 39(1): 77–86 Spurk J H. On the gas loss from ventilated supercavities. Acta Mech, 2002, 155(3-4): 125–135 Thorpe R B, Evans G M, Zhang K. Liquid recirculation and bubble breakup beneath ventilated gas cavities. Chem Eng Sci, 2001, 56(21-22): 6399–6409 Ishii M, Zuber N. Drag coefficient and relative velocity in bubbly, droplet or particulate flows. Am Inst Chem Eng J, 1979, 25(5): 843–

17

18

19 20

21

22

23

24

25

26

February (2012) Vol. 55 No. 2

855 Simonnet M, Gentric C, Olmos E, et.al. Experimental determination of the drag coefficient in a swarm of bubbles. Chem Eng Sci, 2007, 62(3): 858–866 Kunz R F, Gibeling H J, Maxey M R. Validation of two-fluid Eulerian CFD Modeling for microbubble drag reduction across a wide range of reynolds numbers. J Fluids Eng, 2007, 129(1): 66–79 Johansen S T, Boysan F. Fluid dynamics in bubble stirred ladles: Part II. mathematical modeling. Metall Trans B, 1988, 19(5): 755–764 Kumar S, Ramkrishna D. On the solution of population balance equations by discretisation—I. A fixed pivot technique. Chem Eng Sci, 1996, 51(8): 1311–1332 Prince M J, Blanch H W. Bubble coalescence and break-up in air sparged bubble columns. Am Inst Chem Eng J, 1990, 36(10): 1485– 1499 Luo H, Svendsen H. Theoretical model for drop and bubble break-up in turbulent dispersions. Am Inst Chem Eng J, 1996, 42(5): 1225– 1233 Lopez de Bertodano M, Lahey Jr R T, Jones O C. Development of a k- model for bubbly two-phase flow. J Fluids Eng, 1994, 116: 128–134 Sato Y, Sadatomi M, Sekoguchi K. Momentum and heat transfer in two-phase bubbly flow—I. Int J Multiphase Flow, 1981, 7(2): 167–178 Rhie C M, Chow W L. Numerical study of the turbulent-flow past an airfoil with trailing edge separation. AIAA J, 1983, 21(11): 1525– 1532 Johanssona P B V, George W K. Equilibrium similarity, effects of initial conditions and local Reynolds. Phys Fluids, 2003, 15(3): 603–617