- Email: [email protected]
Numerical investigation on the bubble size distribution around NACA0015 hydrofoil
Ocean Engineering 172 (2019) 59–71
Contents lists available at ScienceDirect
Ocean Engineering journal homepage: www.elsevier.com/locate/oceaneng
Numerical investigation on the bubble size distribution around NACA0015 hydrofoil
T
Sara Vahajia, Jiang Hanb, Sherman C.P. Cheungc, Guan H. Yeohd, Jiyuan Tuc,* a
School of Engineering, Deakin University Geelong, Victoria, 3217, Australia Naval Architecture and Ocean Engineering College, Dalian Maritime University, Dalian, 116026, China c School of Engineering, RMIT University Plenty Rd, Bundoora, Victoria, 3083, Australia d UNSW, Sydney, Australia b
A R T I C LE I N FO
A B S T R A C T
Keywords: Gas-liquid flow Bubble size distribution Hydrofoil Cavitation
Two-phase bubbly flows occur widely in nature and are extensively applied in industry. The aeration processes underwater is one type of two-phase bubbly flow that directly impacts on the downstream water quality by reducing the oxygen content in the water. The most important influencing factors for optimization design of Auto-venting turbines (AVTs), for solving the low level of dissolved oxygen (DO) in the discharged downstream water, are the quantity of entrained air, the bubble size distribution resulting from coalescence and breakage processes, and the rate of oxygen transfer from the bubbles. In order to better understand the influencing flow conditions on the bubble size distribution, in this paper a numerical investigation for flow around NACA0015 hydrofoil is carried out. The numerical simulations require the consideration of the dynamic behaviors of twophase flow and bubbles undergoing coalescence and breakup. For this purpose, the ensemble-averaged mass and momentum transport equations for continuous and dispersed phases are modeled within the two-fluid modeling framework. These equations are coupled with population balance equations (PBE) to aptly account for the coalescence and break-up of the bubbles. The resultant bubble size, normalised velocity and void fraction distributions for different flow conditions including angle of attack (AOA), air-entrainment coefficient, and Reynolds number are presented and discussed. The results show that varying AOA has the most significant impact on the distribution of the bubbles in the wake.
1. Introduction Two-phase flows are prevalent in nature and in industry. Cavitation is one type of two-phase flows that occurs in many hydrodynamic applications (e.g. hydrofoils, hydrodynamic pumps, propellers, etc.). It happens when the pressure is dramatically reduced - typically due to high-speed flow or existence of a sharp object in the flow - where the liquid medium breaks down leading to creation of some cavities in the liquid or vaporous bubbles. Cavitation could lead to bubbly flows during the underwater ventilation - where the air is entrained inside water – for example in turbine aeration. Clearly, waters in lakes, rivers, reservoirs and underground aquifers are of most importance for people and to the world's ecological systems (Daskiran et al., 2017). However, hydroelectric facilities discharge water with low dissolved oxygen, which impacts on the downstream water quality, especially at greater depths. Auto-venting turbine (AVT) (March et al., 1992) is one of the
*
techniques proposed for solving the low Dissolved Oxygen (DO) level in the discharged downstream water. AVT selects the locations of the turbine where the aeration could be satisfactorily performed (by increasing the DO of the discharged water while keeping the efficiency of the system in the acceptable range) by computing the pressures, velocities, and turbulence levels throughout a turbine's flow passageway. However, there is limited research on optimization design of AVTs in which the most important influencing factors are the quantity of entrained air, the bubble size distribution (resulted from coalescence and breakage processes), and the rate of oxygen transfer from the bubbles. The interactions of the bubbles including the bubble coalescence and breakup as well as the bubble deformation process determine the bubble size distribution in the wake of the hydrofoil. These processes are governed by the flow conditions and system configurations including the Reynolds number, air ventilation flow rate, and angle of attack of hydrofoil. The impact of the operating conditions such as the air ventilation
Corresponding author. E-mail address: [email protected] (J. Tu).
https://doi.org/10.1016/j.oceaneng.2018.11.045 Received 17 April 2018; Received in revised form 2 November 2018; Accepted 25 November 2018 0029-8018/ © 2018 Elsevier Ltd. All rights reserved.
Ocean Engineering 172 (2019) 59–71
S. Vahaji et al.
Continuity equation of vapour phase
rate, liquid velocity, etc. on the resultant bubble size distribution in the wake, which will subsequently influence on the dissolved oxygen amount in the downstream water through the bubbles interfacial area (Thompson and Gulliver, 1997), needs to be investigated. One of the methods to investigate the bubble size distribution, which would help to understand the underlying physics and tackle the problem of optimization design of AVTs, is through computational fluid dynamics. Recently Daskiran (Daskiran et al., 2017) investigated the dissolved oxygen characteristics in water through computational fluid dynamics. They conducted a parametric study for a two-dimensional geometry by adopting Eulerian multiphase model. Their results could help the design and optimization of energy harvesting devices like translating blades. Zhou et al. (2013) employed direct numerical simulations (DNS) coupled with Lagrangian particle dynamics (LPD) to numerically investigate the bubble dispersion over a hydrofoil. They found that bubble stream deflection increases with increasing angle of attack. Hsiao et al. (2013) performed a numerical and experimental investigation on bubble entrainment and the resultant bubbly wake due to horizontal plunging jet. They employed an Eulerian/Lagrangian oneway coupled two-phase flow model, which could capture the flow structures, velocity field, and overall bubble spreading region near the plunging region. Krepper et al. (2009) investigated the polydispersed bubbly two-phase flow around an obstacle. They considered the coalescence and breakup phenomena and momentum interphase transfers related to drag and lift forces. Also, Kerdouss et al. (2006) carried out simulations for a double impeller stirred tank, since the stirred tanks could improve the dissolution performance through the momentum mixing of one or more impellers. They could obtain the bubble diameter in the tank and the gas volume fraction. They observed that the impeller induces small eddies that break the bubbles; hence, smaller bubbles are located around the impeller. Min et al. (2008) conducted CFD simulations for three-impeller agitators. They concluded that bubble coalescence and breakup needs to be accounted for accurately capturing a local gas holdup in the stirred tank. Consequently, a detailed numerical simulation considering coalescence and breakup of the bubbles is necessary to obtain the bubble size distribution. In this paper, the numerical investigation on the bubble size distribution around NACA0015 hydrofoil is conducted. To ensure that our model on population balance modelling (Cheung et al., 2010), is capturing the coalescence and breakage processes (Cheung et al., 2013; Deju et al., 2013), and the proposed mechanistic model for predicting the bubble size distribution in bubbly flows (Cheung et al., 2014; Yeoh et al., 2014; Vahaji et al., 2017) performs correctly, the results are validated against the existing experimental data on a ventilated hydrofoil in a closed loop water channel performed by Karn et al. (Karn et al., 2015a,b,c; Karn et al., 2015a,b,c; Karn et al., 2015a,b,c; Karn et al., 2016). Accordingly, a systematic numerical investigation on bubbly wake of NACA0015 hydrofoil is carried out under different conditions including angle of attack (AOA), air-entrainment coefficient, and Reynolds number to allow the study of their influence on the bubble size distribution in the bubbly wake.
∂ρg αg fi ∂t
(2)
Momentum equation of liquid phase
∂ρl αl ⇀ ul ul ul ) = −αl ∇P + αl ρl⇀ g + ∇ [αl μle (∇⇀ ul + (∇⇀ ul )T )] + ∇ . (ρl αl ⇀⇀ ∂t ug − Γgl ⇀ ul ) + Flg + (Γlg ⇀ (3) Momentum equation of vapour phase
∂ρg αg ⇀ ug ∂t
+ ∇ . (ρg αg ⇀ ug ⇀ ug ) = −αg ∇P + αg ρg ⇀ g + ∇ [αg μge (∇⇀ ug + (∇⇀ ug )T )] + (Γgl ⇀ ul − Γlg ⇀ ug ) + Fgl
(4)
In equation (2) Si denotes the additional source terms due to coalescence and break-up for the range of bubble classes that are present in the population balance model for the vapour phase. The interfacial force Flg appearing in equation (3) is formulated through appropriate consideration of different sub-forces affecting the interface between each phase. For the liquid phase, the interfacial force comprises the sum of the sub-forces such as drag, lift, wall lubrication and turbulent dispersion respectively. Note that for the gas phase, Fgl& #x202F;= − Flg. Detail descriptions of these forces can be found in (Anglart and Nylund, 1996) and (Lahey Jr and Drew, 2001). Briefly, Interphase momentum transfer between gas and liquid due to drag force is given by
Flgdrag =
1 ⎯u − → ⎯u − → CD aif ρl | ⎯→ ul|( ⎯→ ul ) g g 8
(5)
Lift force in terms of the slip velocity and the curl of the liquid phase velocity is described by
⎯u − → Flglift = αg ρl CL ( ⎯→ ul ) × (∇ × → ul ) g
(6)
Wall lubrication force, is the force that appears when the rising bubble approaches the wall. At the close proximity of the wall, drainage of the fluid around the bubble - which is usually uniform - changes significantly. The drainage rate on the wall side decreases (due to noslip boundary condition), and on the opposite side increases, which results in creating a hydrodynamic force that acts in the normal direction away from the wall and decays with distance from the wall. The wall lubrication force is expressed by
Flglubrication = −
⎯u − → αg ρg ( ⎯→ ul ) g Ds
max (0, Cw1 + Cw2
Ds → )n yw
(7)
Turbulent dispersion force acts as a turbulent diffusion in dispersed flows. For example for a boiling flow in a heated vertical pipe, vapour is generated on heated wall surfaces. The turbulent dispersion force plays a crucial role in driving the vapour away from the vicinity of the wall towards the centre of the pipe. Turbulence dispersion taken as a function of turbulent kinetic energy and gradient of the void fraction of the liquid yields in the form of
2. Mathematical modelling 2.1. Two-fluid model Ensemble-averaged of mass, momentum and energy transport equations are considered for each phase in the Eulerian-Eulerian modelling framework. The liquid phase is the continuum, whose volume fraction is shown as αl, and the vapour phase (bubbles) is the disperse phase, whose volume fraction is shown as αg. These equations can be written as: Continuity equation of liquid phase
∂ρl αl ul ) = Γlg + ∇ . (ρl αl ⇀ ∂t
+ ∇ . (ρg αg ⇀ ug ) = Si − fi Γlg
Flgturbulent dispersion = −CTD ρl κ∇αl
(8)
The drag coefficient CD in equation (5) has been correlated for several distinct Reynolds number regions for individual bubbles according to Ishii and Zuber (1979). The constant CL takes a value of 0.01 (Wang et al., 1987). The wall lubrication constants Cw1 and Cw2 as suggested by Antal et al. (1991) are −0.01 and 0.05 respectively. Recommended value for CTD according to Kurul and Podowski (1990) of 0.1 is used for the turbulent dispersion force.
(1) 60
Ocean Engineering 172 (2019) 59–71
S. Vahaji et al.
2.2. Population balance modelling
a (Mi , Mj ) =
Population balance equations (PBEs) have been applied in many diverse applications that involve particulate systems. Theoretically, the changes of bubble size distribution and its evolution can be captured by the population balance equation (PBE). PBEs are developed in integrodifferential form that are accompanied by coalescence and breakup kernels. In order to solve the PBEs for practical systems, many numerical methods have been proposed. The class method (CM) is one of the existing methods for solving PBE that is widely adopted. In CM, the internal coordinate, for example the particle volume, is discretised into a finite series of classes (bins). In literature, encouraging results can be found for using CM in the form of MUSIG (Multi-SIze Group) model for simulating bubbly flows (Olmos et al., 2001; Pohorecki et al., 2001; Frank et al., 2004; Yeoh and Tu, 2004; Yeoh and Tu, 2005; Cheung et al., 2007a,b; Cheung et al., 2007a,b; Cheung et al., 2008). In the MUSIG model, a number of transport equations and mechanistic bubble coalescence and breakup kernels are solved for each discrete bubble class - in the CFD framework - to predict the evolution of the bubble size distribution. In PBE a desired quantity (e.g. volume fraction or the number density) is conserved for each discrete bubble class, whose conservation is evaluated through coalescence and breakup mechanisms. The solution of these equations results in the population balance of bubbles. The particle (bubble) size distribution is calculated according to the population balance equation that is generally expressed in an integro-differential form:
∂f (x , ξ , t ) + ∇ . (V (x , ξ , t ) f (x , ξ , t )) = S (x , ξ , t ) ∂t
∂t
+ ∇ . (ug ρjg αjg fi ) = Si
1 3 /(ϵ l) 3
where τij is the contact time for two bubbles given by (dij /2) and tij is the time required for two bubbles to coalesce having diameters di and dj which is estimated to be [ϑ]0.5 ln(h 0 / hf ) . The equivalent diameter dij is calculated based on the proposal by (Chesters and Hofman, 1982) dij = (2/ di + 2/ dj )−1. The initial film thickness h 0 = 1 × 10−4 m and critical film thickness hf = 1 × 10−8 m at which rupture occurs for air–water systems have been employed. The turbulent velocity ut in the inertial sub-range of isotropic turbulence is given by (Rotta): 1
1
ut = 2 (ϵl) 3d 3 . For breakup of fluid particles, the partial breakage frequency r (Mi , Mj ) is a function of total breakage frequency, r (Mi ) and the daughter size distribution, β (Mi , Mj ) . β (Mi , Mj ) =
r (Mi , Mj ) (13)
r (Mi )
The bubble binary break-up under isotropic turbulence situation is considered for the bubble breakage proposed by Luo and Svendsen (1996). In this model, a stochastic break-up volume fraction fBV is used for the daughter size distribution. The break-up rate in terms of mass can be obtained as:
r (Mi , Mj ) = C (1 − α g )(
ϵ 13 ) dj2 ξ
1
∫
(1 + ξ )2 ξ
min
11
3
× exp(−
12cf σ ) dξ 2 5 11 l βρ (ϵl) 3d 3ξ 3 (14)
2
2
Where cf = [f BV3 + (1 − fBV ) 3 − 1] is the increase coefficient of surface area, ξ = λ / dj is the size ratio between an eddy and a particle in the inertial sub-range and consequently ξmin = λmin / dj and C = 0.923 and β = 2.0 are determined from fundamental consideration of the break-up of drops or bubbles in turbulent dispersion systems.
(9)
3. Experimental details Experimental data of Karn et al. (Karn et al., 2015a,b,c) is used in this study for investigating the bubble size distribution at the wake of the hydrofoil at different flow conditions. The experiments were conducted in the SAFL high-speed water tunnel at University of Minnesota. The test section's dimensions in the water tunnel are 1 m (length) × 0.19 m (width) ×& #x202F;0.19 m (height). The NACA0015 hydrofoil used in the experiments has a span of 190 mm and a chord length of 81& #x202F;mm. A narrow span-wise slot was carved near the leading edge of the hydrofoil so that air could be injected for creating a dense bubbly wake. They used Particle Shadow Velocimetry (PSV) technique to capture the bubble images at different downstream locations in the wake. They carried out three different sets of experiments. In the first set of experiments, the air entrainment coefficient – CQ = Q/ UcS , where Q,U,c and S are air ventilation flow-rate, liquid velocity, hydrofoil chord and span, respectively – and Reynolds number, Re = Uc / ν , were fixed and the angle of attack was changed. The other two sets of experiments were performed at zero angle of attack. In the second set of experiments, the air entrainment coefficient was fixed and experiments were repeated at different Reynolds numbers; whereas, in the third set of experiments, the Reynolds number was fixed and the air entrainment coefficient was changing.
(10)
In the above equation, Si represents the net change in the number density distribution due to coalescence and break-up processes. This entails the use of a fixed non-uniform volume distribution along a grid, which allows a range of large sizes to be covered with a small number of bins and yet still offers good resolution. Such discretisation of the population balance equation has been found to allow for accurate determination of the desired characteristics of the number density distribution. The interaction term Si = (BC + BB + DC + DB ) contains the source rates of BC , BB , DC and DB , which are the birth rates due to coalescence (BC) and break-up (BD) and the death rates to coalescence (DC) and break-up (BB) of bubbles respectively. For coalescence between fluid particles, the coalescence efficiency a (Mi , Mj ) could be calculated as a product of collision frequency, h (Mi , Mj ) and coalescence efficiency, λ (Mi , Mj ) .
a (Mi , Mj ) = h (Mi , Mj ) λ (Mi , Mj )
(12) 2
where f (x , ξ , t ) is the particle (bubble) number density distribution per unit mixture and particle (bubble) volume, V (x , ξ , t ) is velocity vector in external space dependent on the external variables x for a given time t and the internal space ξ whose components could be characteristic dimensions such as volume, mass etc. On the right hand side, the term S (x , ξ , t ) contains the particle (bubble) source/sink rates per unit mixture volume due to the particle (bubble) interactions such as coalescence, break-up and phase change. Homogeneous MUSIG represents the most commonly used technique for solving PBE. The discrete form of the number density equation, expressed in terms of size fraction fi of M bubble size groups, can be written as:
∂ρjg αjg fi
tij π [di + dj]2 (uti2 + utj2)0.5exp(− ) 4 τij
4. Numerical details
(11) The conservation equations for mass, momentum and energy of each phase are discretised using the finite volume technique for a 2D geometry in ANSYS CFX 17.2. The geometry was borrowed from Karn et al. (Karn et al., 2015a,b,c) where a NACA0015 hydrofoil with the span of 0.19 m and a chord length of 0.081 m was
In this paper, the turbulent random collision is considered for the bubble coalescence proposed by Prince & Blanch (Prince and Blanch, 1990). The coalescence rate in terms of mass can be expressed as (Liao and Lucas, 2010): 61
Ocean Engineering 172 (2019) 59–71
S. Vahaji et al.
Fig. 1. Computational domain.
Fig. 3 demonstrates the mesh around the hydrofoil in the computational domain. In this study, sixteen bubble classes with minimum diameter of 0.000125 m and the maximum diameter of 0.002 m are initialised for the dispersed phases to represent the characteristics of the dispersed phase in the MUSIG model. This means sixteen more transport equations which are progressively solved and coupled with the flow equations during the simulations. The turbulent flows are simulated based on Reynolds Averaged Navier-Stokes (RANS) equations for both carrier and dispersed phases. Fig. 2. Definition of injection angle.
5. Results and discussion Two-phase gas-liquid flows in the wake of NACA0015 at different conditions were simulated to study the effects of air discharging configurations on the interactions between the carrier fluid flow and bubbles. The investigated flow conditions are summarized in Tables 1–3. In these tables, Q is air volume flow-rate, U is liquid velocity, c is the hydrofoil chord length (0.081 m), S is hydrofoil span length (0.19 m), ν is kinematic viscosity, μ is Dynamic viscosity, and AOA is angle of attack. Due to the different interactions between the turbulent flow and bubbles, the final distributions of bubbles are varying, depending on the discharging configurations. 5.1. Validation The simulation results are first validated against the existing experimental data of Karn et al. (Karn et al., 2015a,b,c). For this purpose, the bubble size distribution profiles in the form of probability density function (PDF) for two different air-discharging conditions are compared against the existing experimental data in Fig. 4. In this figure, the angle of attack (AOA) is zero and Reynolds number (Re) is equal to 4.1e5. The air-entrainment coefficient (CQ) is 1.6e-4 for Fig. 4a and CQ& #x202F;= 1.1e-4 for Fig. 4b. The length scale used in defining the Reynolds number is hydrofoil chord: Re = Uc / ν and the airentrainment coefficient was borrowed from Karn et al. (Karn et al., Table 1 Flow conditions at different air-entrainment coefficients (AOA =& #x202F;0, U = 5 m/s, Re =& #x202F;4.1e5).
Fig. 3. Computational grid around hydrofoil.
used. Air was discharged from the injection slot at an angle of 45° to the hydrofoil chord. The computational domain, air injection slot and specific definition of injection angle are illustrated in Figs. 1 and 2, respectively. Appropriate meshing techniques were employed to mesh the computational domain with structured grids (quadrilateral grid), resulting in 174200 grids in total. The near wall regions are refined for more accurate predictions of the flow behaviour in the boundary layer.
Air entrainment Coefficient (CQ = 1.1e-4 2.2e-4 4.4e-4
62
Q ) UcS
Volume flow rate (m3/s)
Mass flow rate (kg/s)
8.3e-6 1.69e-5 3.39e-5
9.9e-6 2.01e-5 4.01e-5
Ocean Engineering 172 (2019) 59–71
S. Vahaji et al.
Table 2 Flow conditions at different AOAs (CQ = 1.6e-4, U& #x202F;= 5 m/s, Re = 4.1e5). Volume flow rate (m3/s)
Mass flow rate (kg/s)
AOA (degree)
1.23e-5 1.23e-5 1.23e-5
1.46e-5 1.46e-5 1.46e-5
0 6 10
Table 3 Flow conditions at different Reynolds numbers (AOA = 0, CQ = 1.1e-4). Volume flow rate (m3/s)
Mass flow rate (kg/s)
Re(Re =
4.9e-6 8.3e-6 1.26e-5 1.67e-5
5.9e-6 9.9e-6 1.49e-5 1.98e-5
2.4e5 4.1e5 6.1e5 8.1e5
Uc ) ν
U (m/s) 2.93 5.00 7.44 9.88
Fig. 5. Bubble diameter distribution at 4.65c from the hydrofoil centre of rotation under different AOA.
Fig. 4. Comparison of numerical and experimental PDFs at 4.65c from the hydrofoil centre of rotation.
2015a,b,c) studies, which is defined similar to other ventilated cavitation studies (Laali and Michel, 1984). In these figures the bubble size distribution is bimodal, which is similar to the findings of other researchers (Tse et al., 2003; Karn et al., 2015a,b,c). Two modes are at 0.3 and 0.62 mm, with more bubbles in the smaller size (first mode) that shows many small bubbles exist in the bubbly wake of the hydrofoil. 5.2. Impact of different angle of attacks (AOA) The impact of changing the angle of attack of hydrofoil on bubble size distribution is investigated in this section. Fig. 5 shows the bubble diameter distribution at 4.65c downstream the hydrofoil centre of rotation for different AOAs. By increasing the angle of attack, the bubble size decreases; it is evident from the bubble size distribution profile where the smallest size group significantly increases from 10% at zero AOA to 17% at 6° AOA, to 30% at 10° AOA. The reason goes back to the fact that as the air is discharged from the ventilation slot, it gets disturbed by the liquid turbulent velocity fluctuations, which leads to breakup of the air jet into small bubbles. On one hand, the turbulence in the wake downstream of a hydrofoil directly impacts on the bubble size and the velocity in the wake. On the other hand, the bubble velocity fluctuations impact on the local liquid velocities. Fig. 6 depicts the wake velocity distribution at 4.65c downstream of the hydrofoil centre of rotation for different AOA. This figure shows that the bubble velocity in the wake varies within 15%, with the lowest velocities near the centreline of wake. It shows that the velocities decrease vertically from both sides; however, the wake is not symmetrical about the line where the minimum velocity occurs. This could be caused by the dispersion of the bubbles with different sizes at different locations. For example, larger bubbles tend to rise within the wake due to buoyancy force. Karn et al. (Karn et al., 2015a,b,c) observed that the larger bubbles have slightly lower velocities compared to tiny bubbles that follow the liquid streamlines more closely (Karn et al., 2015a,b,c); however, since in
Fig. 6. Wake velocity distribution at 4.65c from the hydrofoil centre of rotation under different AOA.
these simulations we are adopting homogeneous MUSIG population balance modelling, this phenomenon could not be reflected in our simulation results. In our results, the wake velocity distribution is almost symmetrical around the minimum velocity line. Another interesting observation is that within the same location downstream of the hydrofoil, increasing the AOA causes the minimum velocity to drift upwards and experience more fluctuations in the wake velocity from 9% variation at zero AOA, to 13% at 6° AOA, to 15% at 10° AOA. The void fraction contours for varying AOA are presented in Fig. 7. This figure clearly demonstrates the vortex generation and vortex shedding phenomena by increasing the AOA to 10°. In order to have a closer examination on the void fraction distribution downstream the discharging slot, the local void fractions at 1.34c, 3c, and 4,65c for varying AOA are presented in Fig. 8. In this figure, the maximum local void fraction for AOA = 0 (Fig. 8a) can get as high as 0.951 at 1.34c downstream hydrofoil and is reached below the centreline, which is caused based on the fact that air is discharged 63
Ocean Engineering 172 (2019) 59–71
S. Vahaji et al.
(a) 00
(b) 60
(c) 100 Fig. 7. Void fraction under different AOA.
the void fraction is reduced farther away from the hydrofoil (less void fraction is observed at 3c and much less at 4.65c) as shown if Fig. 8a–c. The greater void fraction close to the discharging slot could be
from the slot that is located at 6 mm below the centreline. The void fraction distribution along transverse (vertical) direction is almost symmetrical about the location of discharging slot. Also, as expected, 64
Ocean Engineering 172 (2019) 59–71
S. Vahaji et al.
(a) 00
(b) 60
(c) 100 Fig. 8. Local void fraction at 1.34c, 3c and 4.65c under different angles of attack.
average is evaluated on the vertical lines (y-dimension) where Y is limited to −0.01 m < Y < & #x202F;0.015 m to capture the location where the most variation of the void fraction is happening. In this figure, it is clearly shown that by increasing the angle of attack, less void fraction is observed at farthest downstream. Also, the oscillation of vortices impacts on the distribution of the void fraction, which is clearly observed in the case of angle of attack equal to 10°.
attributed to both larger number of bubbles and greater size of bubbles. Moreover, the maximum void fraction value decreases by increasing AOA: from 0.951 at AOA = 0, to 0.437& #x202F;at AOA = 6°, to 0.125 at AOA& #x202F;= 10°, which indicates that the bubbles are dispersed more effectively at larger AOAs. The variation of AOA of hydrofoil is expected to also cause a drift in the bubbly wake downstream. In order to examine the dispersion of the gas phase more closely, the area average of gas phase void fraction vs. x-direction for different angles of attack is demonstrated in Fig. 9. In this figure, the area 65
Ocean Engineering 172 (2019) 59–71
S. Vahaji et al.
Fig. 9. Area average void fraction along the x-direction for different angles of attack.
5.3. Impact of different air-entrainment coefficients In this section, the influence of air-entrainment coefficient (CQ) on the bubble size distribution is studied. Fig. 10 presents the bubble size distribution at 4.65c downstream the hydrofoil at the centre of rotation with air-entrainment coefficients changing from 1.1e-4 to 4.4e-4. As the CQ is increased, bubbles' interactions and collisions increase that lead to greater breakup rate of bubbles resulting in smaller size of the bubbles. Therefore, at higher CQ, there is a greater fraction of small size bubbles. Interestingly, there is also a greater fraction of medium size bubbles for higher CQ. The experiments of Karn et al. (Karn et al., 2015a,b,c) suggest that bubbles are more closely spaced at a higher CQ. These bubbles located at the vicinity of each other would interact effectively resulting in the coalescence in the bubbly wake, which leads to increase in the number of bubbles in intermediate size range as well (Tse et al., 2003). It is also expected that since by increasing CQ the net amount of air ventilated per unit time is also increased, larger values of CQ result in a larger number of bubbles. In Fig. 11, the wake velocity distribution for varying CQ at 4.65c downstream hydrofoil is demonstrated. The bubble velocity in the wake for all CQ conditions varies within 10% of the stream velocity. Also, as mentioned earlier, the minimum velocity occurs below the centreline of the hydrofoil due to the air-discharge at 6&
Fig. 11. Wake velocity distribution at 4.65c from the hydrofoil centre of rotation under different air ventilation flow rate.
#x202F;mm below the hydrofoil's centreline. In addition, the largest value of CQ experiences the highest value of velocity amongst the cases studied. In Fig. 12, the void fraction contour plots for different CQ conditions are presented. Clearly, by increasing the CQ, larger areas with high void fraction have resulted. Consequently, when the air is discharged at higher CQ, the bubbles are dispersed at higher concentrations with larger values of void fraction in similar areas compared to lower values of CQ. The wake is spread more in larger CQ value, which leads to greater breakup events and results in having smaller size bubbles as was earlier shown in Fig. 10. Fig. 13 provides the prospect for a closer examination of local void fraction plot. In this figure, the local void fraction contour plots for varying CQ at 1.34c, 3c, and 4.65c downstream of the hydrofoil are presented. The maximum value of void fraction varies by changing the air ventilation rate: 0.53 at CQ equal to 1.1e-4, to 0.998 at CQ equal to 2.2e-4, to 1.0& #x202F;at CQ equal to 4.4e-4. For all of the cases, the maximum void fraction occurs below the hydrofoil centreline due to discharging air from the slot that is located below the centreline. The bubbles disperse in an almost symmetrical distribution in the transverse (vertical) direction resulting in a similar void fraction distribution for all cases.
5.4. Impact of different Reynolds numbers Changing the flow condition by varying the Reynolds number has direct impact on the bubble size distribution at the wake of the hydrofoil; however, it has insignificant impact on both the average bubble size and the void fraction distribution in the wake of the hydrofoil. In Fig. 14 the bubble size distribution at 4.65c downstream of the hydrofoil at the centre of rotation for different Reynolds numbers changing from 2.4e5 to 8.1e5 is presented. This figure shows that by increasing the Reynolds number, the fraction of smallest bubble size increases from 7% to 27% for the range of Reynolds numbers investigated in this study. However, opposite trend is observed for the second smallest bubble size (0.25 mm): by increasing the Reynolds number, the fraction of the bubble decreases from 30% to 13%. Consequently, changing the Reynolds number does not have a significant impact on the average bubble size in the wake of the hydrofoil. This finding is aligned with experimental findings of Karn et al. (Karn et al., 2015a,b,c). In Fig. 15 the wake velocity distribution at different Reynolds numbers is presented. The velocity in this graph is
Fig. 10. Bubble diameter distribution at 4.65c from the hydrofoil centre of rotation under different air ventilation flow rate. 66
Ocean Engineering 172 (2019) 59–71
S. Vahaji et al.
(a) CQ
1.1e-4
(b) CQ
2.2e-4
(c) CQ
4.4e-4
Fig. 12. Void fraction under different air ventilation flow rate.
67
Ocean Engineering 172 (2019) 59–71
S. Vahaji et al.
(a) CQ
1.1e-4
(b) CQ
2.2e-4
(c) CQ
4.4e-4
Fig. 13. Local void fraction at 1.34c, 3c and 4.65c under different air ventilation flow rate.
68
Ocean Engineering 172 (2019) 59–71
S. Vahaji et al.
(a) Re=2.4e5
Fig. 14. Bubble diameter distribution at 4.65c from the hydrofoil centre of rotation under different Re
(b) Re=4.1e5
(c) Re=6.1e5
Fig. 15. Wake velocity distribution at 4.65c from the hydrofoil centre of rotation under different Re
normalised versus the free stream velocity. Therefore, all three conditions exhibit similar behaviour in the wake in terms of non-dimensional velocity distribution. The void fraction contour plots for different Reynolds numbers are presented in Fig. 16. In this figure, by increasing the Reynolds number no significant variation is observed in the void fraction distribution in the wake of the hydrofoil. In order to have a closer examination of the void fraction distribution for varying Reynolds numbers, in Fig. 17 the void fraction at 1.34c, 3c and 4.65c downstream the hydrofoil for 2.1e5, 4.1e5, 6.1e5 and 8.1e5 is displayed. The maximum void fraction in this figure slightly changes from 0.521 to 0.537 by increasing the Reynolds number from 2.1e5 to 8.1e5.
(d) Re=8.1e5 Fig. 16. Void fraction under different Reynolds numbers.
wake of NACA0015 hydrofoil is investigated and validated against the experimental data of Karn et al. (Karn et al., 2015a,b,c). A systematic numerical investigation on bubbly wake of NACA0015 hydrofoil is carried out under different conditions including angle of attack (AOA), air-entrainment coefficient, and Reynolds number to allow the study of the influence of these parameters on the bubble size distribution in the bubbly wake. The results show that AOA has significant impact on the bubble size distribution, velocity distribution and void fraction distribution in the wake. The bubble size is reduced by increasing the
6. Conclusion In this paper, numerical simulation of two-phase bubbly flow in the 69
Ocean Engineering 172 (2019) 59–71
S. Vahaji et al.
(a) Re=2.4e5
(b) Re=4.1e5
(c) Re=6.1e5
(d) Re=8.1e5 Fig. 17. Local void fraction at 1.34c, 3c and 4.65c under different Reynolds numbers.
70
Ocean Engineering 172 (2019) 59–71
S. Vahaji et al.
angle of attack due to the higher turbulence imposed by larger AOA. Change in AOA also creates a drift in the location where the minimum velocity occurs in the wake, which is an upward drift by increasing the AOA. Also, by increasing the AOA, the bubbles are dispersed more effectively within the wake. Increasing the air-entrainment coefficient results in production of smaller size bubbles due to higher collision and break-up rate of bubbles. Also, with higher air-entrainment coefficient, larger areas with higher void fraction are observed downstream the hydrofoil. Changing the Reynolds number does not have significant influence on the average bubble size in the wake of the hydrofoil. Similarly, the void fraction distribution is not drastically changed by varying the Reynolds number. However, Reynolds number has direct impact on the bubble size distribution at the wake of the hydrofoil.
bubble entrainment due to a horizontal plunging jet. Int. Shipbuild. Prog. 60 (1–4), 435–469. Ishii, M., Zuber, N., 1979. Drag coefficient and relative velocity in bubbly, droplet or particulate flows. AIChE J. 25 (5), 843–855. https://doi.org/10.1002/aic. 690250513. Karn, A., Ellis, C., Hong, J., Arndt, R.E., 2015a. Investigations into the turbulent bubbly wake of a ventilated hydrofoil: moving toward improved turbine aeration techniques. Exp. Therm. Fluid Sci. 64, 186–195. Karn, A., Ellis, C.R., Milliren, C., Hong, J., Scott, D., Arndt, R.E., Gulliver, J.S., 2015b. Bubble size characteristics in the wake of ventilated hydrofoils with two aeration configurations. International Journal of Fluid Machinery and Systems 8 (2), 73–83. Karn, A., Monson, G.M., Ellis, C.R., Hong, J., Arndt, R.E., Gulliver, J.S., 2015c. Mass transfer studies across ventilated hydrofoils: a step towards hydroturbine aeration. Int. J. Heat Mass Tran. 87, 512–520. Karn, A., Shao, S., Arndt, R.E., Hong, J., 2016. Bubble coalescence and breakup in turbulent bubbly wake of a ventilated hydrofoil. Exp. Therm. Fluid Sci. 70, 397–407. Kerdouss, F., Bannari, A., Proulx, P., 2006. CFD modeling of gas dispersion and bubble size in a double turbine stirred tank. Chem. Eng. Sci. 61 (10), 3313–3322. Krepper, E., Beyer, M., Frank, T., Lucas, D., Prasser, H.-M., 2009. CFD modelling of polydispersed bubbly two-phase flow around an obstacle. Nucl. Eng. Des. 239 (11), 2372–2381. https://doi.org/10.1016/j.nucengdes.2009.06.015. Kurul, N., Podowski, M.Z., 1990. Multi-dimensional effects in forced convection subcooled boiling. In: Ninth Heat Transfer Conference. Hemisphere Publishing Corporation, Jerusalem, Israel. Laali, A., Michel, J., 1984. Air entrainment in ventilated cavities: case of the fully developed “half-cavity”. J. Fluid Eng. 106 (3), 327–335. Lahey Jr., R.T., Drew, D.A., 2001. The analysis of two-phase flow and heat transfer using a multidimensional, four field, two-fluid model. Nucl. Eng. Des. 204 (1), 29–44. Liao, Y.X., Lucas, D., 2010. A literature review on mechanisms and models for the coalescence process of fluid particles. Chem. Eng. Sci. 65 (10), 2851–2864. https://doi. org/10.1016/j.ces.2010.02.020. Luo, H., Svendsen, H.F., 1996. Theoretical model for drop and bubble breakup in turbulent dispersions. AIChE J. 42 (5), 1225–1233. March, P.A., Brice, T.A., Mobley, M.H., Cybularz, J.M., 1992. Turbines for solving the DO dilemma. Hydro Rev. 11 (1), 30–36. Min, J., Bao, Y., Chen, L., Gao, Z., Smith, J.M., 2008. Numerical simulation of gas dispersion in an aerated stirred reactor with multiple impellers. Ind. Eng. Chem. Res. 47 (18), 7112–7117. Olmos, E., Gentric, C., Vial, C., Wild, G. and Midoux, N., 2001. Numerical simulation of multiphase flow in bubble column reactors. Influence of bubble coalescence and break-up. Chem. Eng. Sci.. 56(21–22), 6359-6365. https://doi.org/10.1016/S00092509(01)00204-4. Pohorecki, R., Moniuk, W., Bielski, P. and Zdrójkowski, A., 2001. Modelling of the coalescence/redispersion processes in bubble columns. Chem. Eng. Sci.. 56(21–22), 6157-6164. https://doi.org/10.1016/S0009-2509(01)00214-7. Prince, M.J., Blanch, H.W., 1990. Bubble coalescence and break-up in air-sparged bubble columns. AIChE J. 36 (10), 1485–1499. https://doi.org/10.1002/aic.690361004. Rotta, J., Turbulente Strömungen. 1972. Teubner, Stuttgart. Thompson, E.J., Gulliver, J.S., 1997. Oxygen transfer similitude for vented hydroturbine. J. Hydraul. Eng. 123 (6), 529–538. Tse, K., Martin, T., McFarlane, C., Nienow, A., 2003. Small bubble formation via a coalescence dependent break-up mechanism. Chem. Eng. Sci. 58 (2), 275–286. Vahaji, S., Cheung, S.C.P., Yeoh, G.H., Tu, J., 2017. Investigation of the influence of elevated pressure on subcooled boiling flow—model evaluation toward generic approach. J. Heat Tran. 139 (7), 074501. Wang, S., Lee, S., Jones, O., Lahey, R., 1987. 3-D turbulence structure and phase distribution measurements in bubbly two-phase flows. Int. J. Multiphas. Flow 13 (3), 327–343. Yeoh, G., Vahaji, S., Cheung, S., Tu, J., 2014. Modeling subcooled flow boiling in vertical channels at low pressures–Part 2: evaluation of mechanistic approach. Int. J. Heat Mass Tran. 75, 754–768. Yeoh, G.H., Tu, J.Y., 2004. Population balance modelling for bubbly flows with heat and mass transfer. Chem. Eng. Sci. 59 (15), 3125–3139. https://doi.org/10.1016/j.ces. 2004.04.023. Yeoh, G.H., Tu, J.Y., 2005. Thermal-hydrodynamic modeling of bubbly flows with heat and mass transfer. AIChE J. 51 (1), 8–27. https://doi.org/10.1002/aic.10297. Zhou, Z., Yuan, Q., Jia, X., Feng, W., Wen, J., 2013. Experimental study and CFD simulation ofMass TransferCharacteristics of a gas-induced pulsating flow BubbleColumn. Chem. Biochem. Eng. Q. 27 (2), 167–175.
Acknowledgement The authors are grateful for the support provided by Australia-China Centre for Maritime Engineering (ACCME) (Project No. ACSRF48199), National Natural Science Foundation of China (No.51436002) and the CSC Scholarship (No.201506575024). References Anglart, H. and Nylund, O., 1996. CFD application to prediction of void distribution in two-phase bubbly flows in rod bundles. Nucl. Eng. Des.. 163(1–2), 81-98. https://doi. org/10.1016/0029-5493(95)01160-9. Antal, S., Lahey Jr., R., Flaherty, J., 1991. Analysis of phase distribution and turbulence in dispersed particle/liquid flows. Chem. Eng. Commun. 174, 85–113. Chesters, A., Hofman, G., 1982. Bubble Coalescence in Pure Liquids. Mechanics and Physics of Bubbles in Liquids. Springer, pp. 353–361. Cheung, S., Vahaji, S., Yeoh, G., Tu, J., 2014. Modeling subcooled flow boiling in vertical channels at low pressures–Part 1: assessment of empirical correlations. Int. J. Heat Mass Tran. 75, 736–753. Cheung, S. C. P., Deju, L., Yeoh, G. H. and Tu, J. Y., 2013. Modeling of bubble size distribution in isothermal gas–liquid flows: numerical assessment of population balance approaches. Nucl. Eng. Des.. 265(0), 120-136. https://doi.org/10.1016/j. nucengdes.2013.08.049. Cheung, S. C. P., Duan, X., Yeoh, G. H., Tu, J., Krepper, E. and Lucas, D., 2010. Modelling of Polydispersed Flows Using Two Population Balance Approaches. 6th International Symposium on Multiphase Flow, Heat Mass Transfer and Energy Conversion. L. J. Guo, D. D. Joseph, Y. Matsumoto, M. Sommerfeld and Y. S. Wang. vol. 1207: 809814. Cheung, S. C. P., Yeoh, G. H. and Tu, J. Y., 2007a. On the modelling of population balance in isothermal vertical bubbly flows—average bubble number density approach. Chem. Eng. Process: Process Intensification. 46(8), 742-756. https://doi.org/10. 1016/j.cep.2006.10.004. Cheung, S.C.P., Yeoh, G.H., Tu, J.Y., 2007b. On the numerical study of isothermal vertical bubbly flow using two population balance approaches. Chem. Eng. Sci. 62 (17), 4659–4674. https://doi.org/10.1016/j.ces.2007.05.030. Cheung, S.C.P., Yeoh, G.H., Tu, J.Y., 2008. Population balance modeling of bubbly flows considering the hydrodynamics and thermomechanical processes. AIChE J. 54 (7), 1689–1710. https://doi.org/10.1002/aic.11503. Daskiran, C., Liu, I.-H., Oztekin, A., 2017. Computational study of multiphase flows over ventilated translating blades. Int. J. Heat Mass Tran. 110, 262–275. Deju, L., Cheung, S. C. P., Yeoh, G. H. and Tu, J. Y., 2013. Capturing coalescence and break-up processes in vertical gas–liquid flows: assessment of population balance methods. Appl. Math. Model.. 37(18–19), 8557-8577. https://doi.org/10.1016/j. apm.2013.03.063. Frank, T., Shi, J., Burns, A.D., 2004. Validation of Eulerian multiphase flow models for nuclear safety application. In: Proceeding of the Third International Symposium on Two-phase Modelling and Experimentation, Pisa, Italy. Hsiao, C.-T., Wu, X., Ma, J., Chahine, G., 2013. Numerical and experimental study of
71
Contents lists available at ScienceDirect
Ocean Engineering journal homepage: www.elsevier.com/locate/oceaneng
Numerical investigation on the bubble size distribution around NACA0015 hydrofoil
T
Sara Vahajia, Jiang Hanb, Sherman C.P. Cheungc, Guan H. Yeohd, Jiyuan Tuc,* a
School of Engineering, Deakin University Geelong, Victoria, 3217, Australia Naval Architecture and Ocean Engineering College, Dalian Maritime University, Dalian, 116026, China c School of Engineering, RMIT University Plenty Rd, Bundoora, Victoria, 3083, Australia d UNSW, Sydney, Australia b
A R T I C LE I N FO
A B S T R A C T
Keywords: Gas-liquid flow Bubble size distribution Hydrofoil Cavitation
Two-phase bubbly flows occur widely in nature and are extensively applied in industry. The aeration processes underwater is one type of two-phase bubbly flow that directly impacts on the downstream water quality by reducing the oxygen content in the water. The most important influencing factors for optimization design of Auto-venting turbines (AVTs), for solving the low level of dissolved oxygen (DO) in the discharged downstream water, are the quantity of entrained air, the bubble size distribution resulting from coalescence and breakage processes, and the rate of oxygen transfer from the bubbles. In order to better understand the influencing flow conditions on the bubble size distribution, in this paper a numerical investigation for flow around NACA0015 hydrofoil is carried out. The numerical simulations require the consideration of the dynamic behaviors of twophase flow and bubbles undergoing coalescence and breakup. For this purpose, the ensemble-averaged mass and momentum transport equations for continuous and dispersed phases are modeled within the two-fluid modeling framework. These equations are coupled with population balance equations (PBE) to aptly account for the coalescence and break-up of the bubbles. The resultant bubble size, normalised velocity and void fraction distributions for different flow conditions including angle of attack (AOA), air-entrainment coefficient, and Reynolds number are presented and discussed. The results show that varying AOA has the most significant impact on the distribution of the bubbles in the wake.
1. Introduction Two-phase flows are prevalent in nature and in industry. Cavitation is one type of two-phase flows that occurs in many hydrodynamic applications (e.g. hydrofoils, hydrodynamic pumps, propellers, etc.). It happens when the pressure is dramatically reduced - typically due to high-speed flow or existence of a sharp object in the flow - where the liquid medium breaks down leading to creation of some cavities in the liquid or vaporous bubbles. Cavitation could lead to bubbly flows during the underwater ventilation - where the air is entrained inside water – for example in turbine aeration. Clearly, waters in lakes, rivers, reservoirs and underground aquifers are of most importance for people and to the world's ecological systems (Daskiran et al., 2017). However, hydroelectric facilities discharge water with low dissolved oxygen, which impacts on the downstream water quality, especially at greater depths. Auto-venting turbine (AVT) (March et al., 1992) is one of the
*
techniques proposed for solving the low Dissolved Oxygen (DO) level in the discharged downstream water. AVT selects the locations of the turbine where the aeration could be satisfactorily performed (by increasing the DO of the discharged water while keeping the efficiency of the system in the acceptable range) by computing the pressures, velocities, and turbulence levels throughout a turbine's flow passageway. However, there is limited research on optimization design of AVTs in which the most important influencing factors are the quantity of entrained air, the bubble size distribution (resulted from coalescence and breakage processes), and the rate of oxygen transfer from the bubbles. The interactions of the bubbles including the bubble coalescence and breakup as well as the bubble deformation process determine the bubble size distribution in the wake of the hydrofoil. These processes are governed by the flow conditions and system configurations including the Reynolds number, air ventilation flow rate, and angle of attack of hydrofoil. The impact of the operating conditions such as the air ventilation
Corresponding author. E-mail address: [email protected] (J. Tu).
https://doi.org/10.1016/j.oceaneng.2018.11.045 Received 17 April 2018; Received in revised form 2 November 2018; Accepted 25 November 2018 0029-8018/ © 2018 Elsevier Ltd. All rights reserved.
Ocean Engineering 172 (2019) 59–71
S. Vahaji et al.
Continuity equation of vapour phase
rate, liquid velocity, etc. on the resultant bubble size distribution in the wake, which will subsequently influence on the dissolved oxygen amount in the downstream water through the bubbles interfacial area (Thompson and Gulliver, 1997), needs to be investigated. One of the methods to investigate the bubble size distribution, which would help to understand the underlying physics and tackle the problem of optimization design of AVTs, is through computational fluid dynamics. Recently Daskiran (Daskiran et al., 2017) investigated the dissolved oxygen characteristics in water through computational fluid dynamics. They conducted a parametric study for a two-dimensional geometry by adopting Eulerian multiphase model. Their results could help the design and optimization of energy harvesting devices like translating blades. Zhou et al. (2013) employed direct numerical simulations (DNS) coupled with Lagrangian particle dynamics (LPD) to numerically investigate the bubble dispersion over a hydrofoil. They found that bubble stream deflection increases with increasing angle of attack. Hsiao et al. (2013) performed a numerical and experimental investigation on bubble entrainment and the resultant bubbly wake due to horizontal plunging jet. They employed an Eulerian/Lagrangian oneway coupled two-phase flow model, which could capture the flow structures, velocity field, and overall bubble spreading region near the plunging region. Krepper et al. (2009) investigated the polydispersed bubbly two-phase flow around an obstacle. They considered the coalescence and breakup phenomena and momentum interphase transfers related to drag and lift forces. Also, Kerdouss et al. (2006) carried out simulations for a double impeller stirred tank, since the stirred tanks could improve the dissolution performance through the momentum mixing of one or more impellers. They could obtain the bubble diameter in the tank and the gas volume fraction. They observed that the impeller induces small eddies that break the bubbles; hence, smaller bubbles are located around the impeller. Min et al. (2008) conducted CFD simulations for three-impeller agitators. They concluded that bubble coalescence and breakup needs to be accounted for accurately capturing a local gas holdup in the stirred tank. Consequently, a detailed numerical simulation considering coalescence and breakup of the bubbles is necessary to obtain the bubble size distribution. In this paper, the numerical investigation on the bubble size distribution around NACA0015 hydrofoil is conducted. To ensure that our model on population balance modelling (Cheung et al., 2010), is capturing the coalescence and breakage processes (Cheung et al., 2013; Deju et al., 2013), and the proposed mechanistic model for predicting the bubble size distribution in bubbly flows (Cheung et al., 2014; Yeoh et al., 2014; Vahaji et al., 2017) performs correctly, the results are validated against the existing experimental data on a ventilated hydrofoil in a closed loop water channel performed by Karn et al. (Karn et al., 2015a,b,c; Karn et al., 2015a,b,c; Karn et al., 2015a,b,c; Karn et al., 2016). Accordingly, a systematic numerical investigation on bubbly wake of NACA0015 hydrofoil is carried out under different conditions including angle of attack (AOA), air-entrainment coefficient, and Reynolds number to allow the study of their influence on the bubble size distribution in the bubbly wake.
∂ρg αg fi ∂t
(2)
Momentum equation of liquid phase
∂ρl αl ⇀ ul ul ul ) = −αl ∇P + αl ρl⇀ g + ∇ [αl μle (∇⇀ ul + (∇⇀ ul )T )] + ∇ . (ρl αl ⇀⇀ ∂t ug − Γgl ⇀ ul ) + Flg + (Γlg ⇀ (3) Momentum equation of vapour phase
∂ρg αg ⇀ ug ∂t
+ ∇ . (ρg αg ⇀ ug ⇀ ug ) = −αg ∇P + αg ρg ⇀ g + ∇ [αg μge (∇⇀ ug + (∇⇀ ug )T )] + (Γgl ⇀ ul − Γlg ⇀ ug ) + Fgl
(4)
In equation (2) Si denotes the additional source terms due to coalescence and break-up for the range of bubble classes that are present in the population balance model for the vapour phase. The interfacial force Flg appearing in equation (3) is formulated through appropriate consideration of different sub-forces affecting the interface between each phase. For the liquid phase, the interfacial force comprises the sum of the sub-forces such as drag, lift, wall lubrication and turbulent dispersion respectively. Note that for the gas phase, Fgl& #x202F;= − Flg. Detail descriptions of these forces can be found in (Anglart and Nylund, 1996) and (Lahey Jr and Drew, 2001). Briefly, Interphase momentum transfer between gas and liquid due to drag force is given by
Flgdrag =
1 ⎯u − → ⎯u − → CD aif ρl | ⎯→ ul|( ⎯→ ul ) g g 8
(5)
Lift force in terms of the slip velocity and the curl of the liquid phase velocity is described by
⎯u − → Flglift = αg ρl CL ( ⎯→ ul ) × (∇ × → ul ) g
(6)
Wall lubrication force, is the force that appears when the rising bubble approaches the wall. At the close proximity of the wall, drainage of the fluid around the bubble - which is usually uniform - changes significantly. The drainage rate on the wall side decreases (due to noslip boundary condition), and on the opposite side increases, which results in creating a hydrodynamic force that acts in the normal direction away from the wall and decays with distance from the wall. The wall lubrication force is expressed by
Flglubrication = −
⎯u − → αg ρg ( ⎯→ ul ) g Ds
max (0, Cw1 + Cw2
Ds → )n yw
(7)
Turbulent dispersion force acts as a turbulent diffusion in dispersed flows. For example for a boiling flow in a heated vertical pipe, vapour is generated on heated wall surfaces. The turbulent dispersion force plays a crucial role in driving the vapour away from the vicinity of the wall towards the centre of the pipe. Turbulence dispersion taken as a function of turbulent kinetic energy and gradient of the void fraction of the liquid yields in the form of
2. Mathematical modelling 2.1. Two-fluid model Ensemble-averaged of mass, momentum and energy transport equations are considered for each phase in the Eulerian-Eulerian modelling framework. The liquid phase is the continuum, whose volume fraction is shown as αl, and the vapour phase (bubbles) is the disperse phase, whose volume fraction is shown as αg. These equations can be written as: Continuity equation of liquid phase
∂ρl αl ul ) = Γlg + ∇ . (ρl αl ⇀ ∂t
+ ∇ . (ρg αg ⇀ ug ) = Si − fi Γlg
Flgturbulent dispersion = −CTD ρl κ∇αl
(8)
The drag coefficient CD in equation (5) has been correlated for several distinct Reynolds number regions for individual bubbles according to Ishii and Zuber (1979). The constant CL takes a value of 0.01 (Wang et al., 1987). The wall lubrication constants Cw1 and Cw2 as suggested by Antal et al. (1991) are −0.01 and 0.05 respectively. Recommended value for CTD according to Kurul and Podowski (1990) of 0.1 is used for the turbulent dispersion force.
(1) 60
Ocean Engineering 172 (2019) 59–71
S. Vahaji et al.
2.2. Population balance modelling
a (Mi , Mj ) =
Population balance equations (PBEs) have been applied in many diverse applications that involve particulate systems. Theoretically, the changes of bubble size distribution and its evolution can be captured by the population balance equation (PBE). PBEs are developed in integrodifferential form that are accompanied by coalescence and breakup kernels. In order to solve the PBEs for practical systems, many numerical methods have been proposed. The class method (CM) is one of the existing methods for solving PBE that is widely adopted. In CM, the internal coordinate, for example the particle volume, is discretised into a finite series of classes (bins). In literature, encouraging results can be found for using CM in the form of MUSIG (Multi-SIze Group) model for simulating bubbly flows (Olmos et al., 2001; Pohorecki et al., 2001; Frank et al., 2004; Yeoh and Tu, 2004; Yeoh and Tu, 2005; Cheung et al., 2007a,b; Cheung et al., 2007a,b; Cheung et al., 2008). In the MUSIG model, a number of transport equations and mechanistic bubble coalescence and breakup kernels are solved for each discrete bubble class - in the CFD framework - to predict the evolution of the bubble size distribution. In PBE a desired quantity (e.g. volume fraction or the number density) is conserved for each discrete bubble class, whose conservation is evaluated through coalescence and breakup mechanisms. The solution of these equations results in the population balance of bubbles. The particle (bubble) size distribution is calculated according to the population balance equation that is generally expressed in an integro-differential form:
∂f (x , ξ , t ) + ∇ . (V (x , ξ , t ) f (x , ξ , t )) = S (x , ξ , t ) ∂t
∂t
+ ∇ . (ug ρjg αjg fi ) = Si
1 3 /(ϵ l) 3
where τij is the contact time for two bubbles given by (dij /2) and tij is the time required for two bubbles to coalesce having diameters di and dj which is estimated to be [ϑ]0.5 ln(h 0 / hf ) . The equivalent diameter dij is calculated based on the proposal by (Chesters and Hofman, 1982) dij = (2/ di + 2/ dj )−1. The initial film thickness h 0 = 1 × 10−4 m and critical film thickness hf = 1 × 10−8 m at which rupture occurs for air–water systems have been employed. The turbulent velocity ut in the inertial sub-range of isotropic turbulence is given by (Rotta): 1
1
ut = 2 (ϵl) 3d 3 . For breakup of fluid particles, the partial breakage frequency r (Mi , Mj ) is a function of total breakage frequency, r (Mi ) and the daughter size distribution, β (Mi , Mj ) . β (Mi , Mj ) =
r (Mi , Mj ) (13)
r (Mi )
The bubble binary break-up under isotropic turbulence situation is considered for the bubble breakage proposed by Luo and Svendsen (1996). In this model, a stochastic break-up volume fraction fBV is used for the daughter size distribution. The break-up rate in terms of mass can be obtained as:
r (Mi , Mj ) = C (1 − α g )(
ϵ 13 ) dj2 ξ
1
∫
(1 + ξ )2 ξ
min
11
3
× exp(−
12cf σ ) dξ 2 5 11 l βρ (ϵl) 3d 3ξ 3 (14)
2
2
Where cf = [f BV3 + (1 − fBV ) 3 − 1] is the increase coefficient of surface area, ξ = λ / dj is the size ratio between an eddy and a particle in the inertial sub-range and consequently ξmin = λmin / dj and C = 0.923 and β = 2.0 are determined from fundamental consideration of the break-up of drops or bubbles in turbulent dispersion systems.
(9)
3. Experimental details Experimental data of Karn et al. (Karn et al., 2015a,b,c) is used in this study for investigating the bubble size distribution at the wake of the hydrofoil at different flow conditions. The experiments were conducted in the SAFL high-speed water tunnel at University of Minnesota. The test section's dimensions in the water tunnel are 1 m (length) × 0.19 m (width) ×& #x202F;0.19 m (height). The NACA0015 hydrofoil used in the experiments has a span of 190 mm and a chord length of 81& #x202F;mm. A narrow span-wise slot was carved near the leading edge of the hydrofoil so that air could be injected for creating a dense bubbly wake. They used Particle Shadow Velocimetry (PSV) technique to capture the bubble images at different downstream locations in the wake. They carried out three different sets of experiments. In the first set of experiments, the air entrainment coefficient – CQ = Q/ UcS , where Q,U,c and S are air ventilation flow-rate, liquid velocity, hydrofoil chord and span, respectively – and Reynolds number, Re = Uc / ν , were fixed and the angle of attack was changed. The other two sets of experiments were performed at zero angle of attack. In the second set of experiments, the air entrainment coefficient was fixed and experiments were repeated at different Reynolds numbers; whereas, in the third set of experiments, the Reynolds number was fixed and the air entrainment coefficient was changing.
(10)
In the above equation, Si represents the net change in the number density distribution due to coalescence and break-up processes. This entails the use of a fixed non-uniform volume distribution along a grid, which allows a range of large sizes to be covered with a small number of bins and yet still offers good resolution. Such discretisation of the population balance equation has been found to allow for accurate determination of the desired characteristics of the number density distribution. The interaction term Si = (BC + BB + DC + DB ) contains the source rates of BC , BB , DC and DB , which are the birth rates due to coalescence (BC) and break-up (BD) and the death rates to coalescence (DC) and break-up (BB) of bubbles respectively. For coalescence between fluid particles, the coalescence efficiency a (Mi , Mj ) could be calculated as a product of collision frequency, h (Mi , Mj ) and coalescence efficiency, λ (Mi , Mj ) .
a (Mi , Mj ) = h (Mi , Mj ) λ (Mi , Mj )
(12) 2
where f (x , ξ , t ) is the particle (bubble) number density distribution per unit mixture and particle (bubble) volume, V (x , ξ , t ) is velocity vector in external space dependent on the external variables x for a given time t and the internal space ξ whose components could be characteristic dimensions such as volume, mass etc. On the right hand side, the term S (x , ξ , t ) contains the particle (bubble) source/sink rates per unit mixture volume due to the particle (bubble) interactions such as coalescence, break-up and phase change. Homogeneous MUSIG represents the most commonly used technique for solving PBE. The discrete form of the number density equation, expressed in terms of size fraction fi of M bubble size groups, can be written as:
∂ρjg αjg fi
tij π [di + dj]2 (uti2 + utj2)0.5exp(− ) 4 τij
4. Numerical details
(11) The conservation equations for mass, momentum and energy of each phase are discretised using the finite volume technique for a 2D geometry in ANSYS CFX 17.2. The geometry was borrowed from Karn et al. (Karn et al., 2015a,b,c) where a NACA0015 hydrofoil with the span of 0.19 m and a chord length of 0.081 m was
In this paper, the turbulent random collision is considered for the bubble coalescence proposed by Prince & Blanch (Prince and Blanch, 1990). The coalescence rate in terms of mass can be expressed as (Liao and Lucas, 2010): 61
Ocean Engineering 172 (2019) 59–71
S. Vahaji et al.
Fig. 1. Computational domain.
Fig. 3 demonstrates the mesh around the hydrofoil in the computational domain. In this study, sixteen bubble classes with minimum diameter of 0.000125 m and the maximum diameter of 0.002 m are initialised for the dispersed phases to represent the characteristics of the dispersed phase in the MUSIG model. This means sixteen more transport equations which are progressively solved and coupled with the flow equations during the simulations. The turbulent flows are simulated based on Reynolds Averaged Navier-Stokes (RANS) equations for both carrier and dispersed phases. Fig. 2. Definition of injection angle.
5. Results and discussion Two-phase gas-liquid flows in the wake of NACA0015 at different conditions were simulated to study the effects of air discharging configurations on the interactions between the carrier fluid flow and bubbles. The investigated flow conditions are summarized in Tables 1–3. In these tables, Q is air volume flow-rate, U is liquid velocity, c is the hydrofoil chord length (0.081 m), S is hydrofoil span length (0.19 m), ν is kinematic viscosity, μ is Dynamic viscosity, and AOA is angle of attack. Due to the different interactions between the turbulent flow and bubbles, the final distributions of bubbles are varying, depending on the discharging configurations. 5.1. Validation The simulation results are first validated against the existing experimental data of Karn et al. (Karn et al., 2015a,b,c). For this purpose, the bubble size distribution profiles in the form of probability density function (PDF) for two different air-discharging conditions are compared against the existing experimental data in Fig. 4. In this figure, the angle of attack (AOA) is zero and Reynolds number (Re) is equal to 4.1e5. The air-entrainment coefficient (CQ) is 1.6e-4 for Fig. 4a and CQ& #x202F;= 1.1e-4 for Fig. 4b. The length scale used in defining the Reynolds number is hydrofoil chord: Re = Uc / ν and the airentrainment coefficient was borrowed from Karn et al. (Karn et al., Table 1 Flow conditions at different air-entrainment coefficients (AOA =& #x202F;0, U = 5 m/s, Re =& #x202F;4.1e5).
Fig. 3. Computational grid around hydrofoil.
used. Air was discharged from the injection slot at an angle of 45° to the hydrofoil chord. The computational domain, air injection slot and specific definition of injection angle are illustrated in Figs. 1 and 2, respectively. Appropriate meshing techniques were employed to mesh the computational domain with structured grids (quadrilateral grid), resulting in 174200 grids in total. The near wall regions are refined for more accurate predictions of the flow behaviour in the boundary layer.
Air entrainment Coefficient (CQ = 1.1e-4 2.2e-4 4.4e-4
62
Q ) UcS
Volume flow rate (m3/s)
Mass flow rate (kg/s)
8.3e-6 1.69e-5 3.39e-5
9.9e-6 2.01e-5 4.01e-5
Ocean Engineering 172 (2019) 59–71
S. Vahaji et al.
Table 2 Flow conditions at different AOAs (CQ = 1.6e-4, U& #x202F;= 5 m/s, Re = 4.1e5). Volume flow rate (m3/s)
Mass flow rate (kg/s)
AOA (degree)
1.23e-5 1.23e-5 1.23e-5
1.46e-5 1.46e-5 1.46e-5
0 6 10
Table 3 Flow conditions at different Reynolds numbers (AOA = 0, CQ = 1.1e-4). Volume flow rate (m3/s)
Mass flow rate (kg/s)
Re(Re =
4.9e-6 8.3e-6 1.26e-5 1.67e-5
5.9e-6 9.9e-6 1.49e-5 1.98e-5
2.4e5 4.1e5 6.1e5 8.1e5
Uc ) ν
U (m/s) 2.93 5.00 7.44 9.88
Fig. 5. Bubble diameter distribution at 4.65c from the hydrofoil centre of rotation under different AOA.
Fig. 4. Comparison of numerical and experimental PDFs at 4.65c from the hydrofoil centre of rotation.
2015a,b,c) studies, which is defined similar to other ventilated cavitation studies (Laali and Michel, 1984). In these figures the bubble size distribution is bimodal, which is similar to the findings of other researchers (Tse et al., 2003; Karn et al., 2015a,b,c). Two modes are at 0.3 and 0.62 mm, with more bubbles in the smaller size (first mode) that shows many small bubbles exist in the bubbly wake of the hydrofoil. 5.2. Impact of different angle of attacks (AOA) The impact of changing the angle of attack of hydrofoil on bubble size distribution is investigated in this section. Fig. 5 shows the bubble diameter distribution at 4.65c downstream the hydrofoil centre of rotation for different AOAs. By increasing the angle of attack, the bubble size decreases; it is evident from the bubble size distribution profile where the smallest size group significantly increases from 10% at zero AOA to 17% at 6° AOA, to 30% at 10° AOA. The reason goes back to the fact that as the air is discharged from the ventilation slot, it gets disturbed by the liquid turbulent velocity fluctuations, which leads to breakup of the air jet into small bubbles. On one hand, the turbulence in the wake downstream of a hydrofoil directly impacts on the bubble size and the velocity in the wake. On the other hand, the bubble velocity fluctuations impact on the local liquid velocities. Fig. 6 depicts the wake velocity distribution at 4.65c downstream of the hydrofoil centre of rotation for different AOA. This figure shows that the bubble velocity in the wake varies within 15%, with the lowest velocities near the centreline of wake. It shows that the velocities decrease vertically from both sides; however, the wake is not symmetrical about the line where the minimum velocity occurs. This could be caused by the dispersion of the bubbles with different sizes at different locations. For example, larger bubbles tend to rise within the wake due to buoyancy force. Karn et al. (Karn et al., 2015a,b,c) observed that the larger bubbles have slightly lower velocities compared to tiny bubbles that follow the liquid streamlines more closely (Karn et al., 2015a,b,c); however, since in
Fig. 6. Wake velocity distribution at 4.65c from the hydrofoil centre of rotation under different AOA.
these simulations we are adopting homogeneous MUSIG population balance modelling, this phenomenon could not be reflected in our simulation results. In our results, the wake velocity distribution is almost symmetrical around the minimum velocity line. Another interesting observation is that within the same location downstream of the hydrofoil, increasing the AOA causes the minimum velocity to drift upwards and experience more fluctuations in the wake velocity from 9% variation at zero AOA, to 13% at 6° AOA, to 15% at 10° AOA. The void fraction contours for varying AOA are presented in Fig. 7. This figure clearly demonstrates the vortex generation and vortex shedding phenomena by increasing the AOA to 10°. In order to have a closer examination on the void fraction distribution downstream the discharging slot, the local void fractions at 1.34c, 3c, and 4,65c for varying AOA are presented in Fig. 8. In this figure, the maximum local void fraction for AOA = 0 (Fig. 8a) can get as high as 0.951 at 1.34c downstream hydrofoil and is reached below the centreline, which is caused based on the fact that air is discharged 63
Ocean Engineering 172 (2019) 59–71
S. Vahaji et al.
(a) 00
(b) 60
(c) 100 Fig. 7. Void fraction under different AOA.
the void fraction is reduced farther away from the hydrofoil (less void fraction is observed at 3c and much less at 4.65c) as shown if Fig. 8a–c. The greater void fraction close to the discharging slot could be
from the slot that is located at 6 mm below the centreline. The void fraction distribution along transverse (vertical) direction is almost symmetrical about the location of discharging slot. Also, as expected, 64
Ocean Engineering 172 (2019) 59–71
S. Vahaji et al.
(a) 00
(b) 60
(c) 100 Fig. 8. Local void fraction at 1.34c, 3c and 4.65c under different angles of attack.
average is evaluated on the vertical lines (y-dimension) where Y is limited to −0.01 m < Y < & #x202F;0.015 m to capture the location where the most variation of the void fraction is happening. In this figure, it is clearly shown that by increasing the angle of attack, less void fraction is observed at farthest downstream. Also, the oscillation of vortices impacts on the distribution of the void fraction, which is clearly observed in the case of angle of attack equal to 10°.
attributed to both larger number of bubbles and greater size of bubbles. Moreover, the maximum void fraction value decreases by increasing AOA: from 0.951 at AOA = 0, to 0.437& #x202F;at AOA = 6°, to 0.125 at AOA& #x202F;= 10°, which indicates that the bubbles are dispersed more effectively at larger AOAs. The variation of AOA of hydrofoil is expected to also cause a drift in the bubbly wake downstream. In order to examine the dispersion of the gas phase more closely, the area average of gas phase void fraction vs. x-direction for different angles of attack is demonstrated in Fig. 9. In this figure, the area 65
Ocean Engineering 172 (2019) 59–71
S. Vahaji et al.
Fig. 9. Area average void fraction along the x-direction for different angles of attack.
5.3. Impact of different air-entrainment coefficients In this section, the influence of air-entrainment coefficient (CQ) on the bubble size distribution is studied. Fig. 10 presents the bubble size distribution at 4.65c downstream the hydrofoil at the centre of rotation with air-entrainment coefficients changing from 1.1e-4 to 4.4e-4. As the CQ is increased, bubbles' interactions and collisions increase that lead to greater breakup rate of bubbles resulting in smaller size of the bubbles. Therefore, at higher CQ, there is a greater fraction of small size bubbles. Interestingly, there is also a greater fraction of medium size bubbles for higher CQ. The experiments of Karn et al. (Karn et al., 2015a,b,c) suggest that bubbles are more closely spaced at a higher CQ. These bubbles located at the vicinity of each other would interact effectively resulting in the coalescence in the bubbly wake, which leads to increase in the number of bubbles in intermediate size range as well (Tse et al., 2003). It is also expected that since by increasing CQ the net amount of air ventilated per unit time is also increased, larger values of CQ result in a larger number of bubbles. In Fig. 11, the wake velocity distribution for varying CQ at 4.65c downstream hydrofoil is demonstrated. The bubble velocity in the wake for all CQ conditions varies within 10% of the stream velocity. Also, as mentioned earlier, the minimum velocity occurs below the centreline of the hydrofoil due to the air-discharge at 6&
Fig. 11. Wake velocity distribution at 4.65c from the hydrofoil centre of rotation under different air ventilation flow rate.
#x202F;mm below the hydrofoil's centreline. In addition, the largest value of CQ experiences the highest value of velocity amongst the cases studied. In Fig. 12, the void fraction contour plots for different CQ conditions are presented. Clearly, by increasing the CQ, larger areas with high void fraction have resulted. Consequently, when the air is discharged at higher CQ, the bubbles are dispersed at higher concentrations with larger values of void fraction in similar areas compared to lower values of CQ. The wake is spread more in larger CQ value, which leads to greater breakup events and results in having smaller size bubbles as was earlier shown in Fig. 10. Fig. 13 provides the prospect for a closer examination of local void fraction plot. In this figure, the local void fraction contour plots for varying CQ at 1.34c, 3c, and 4.65c downstream of the hydrofoil are presented. The maximum value of void fraction varies by changing the air ventilation rate: 0.53 at CQ equal to 1.1e-4, to 0.998 at CQ equal to 2.2e-4, to 1.0& #x202F;at CQ equal to 4.4e-4. For all of the cases, the maximum void fraction occurs below the hydrofoil centreline due to discharging air from the slot that is located below the centreline. The bubbles disperse in an almost symmetrical distribution in the transverse (vertical) direction resulting in a similar void fraction distribution for all cases.
5.4. Impact of different Reynolds numbers Changing the flow condition by varying the Reynolds number has direct impact on the bubble size distribution at the wake of the hydrofoil; however, it has insignificant impact on both the average bubble size and the void fraction distribution in the wake of the hydrofoil. In Fig. 14 the bubble size distribution at 4.65c downstream of the hydrofoil at the centre of rotation for different Reynolds numbers changing from 2.4e5 to 8.1e5 is presented. This figure shows that by increasing the Reynolds number, the fraction of smallest bubble size increases from 7% to 27% for the range of Reynolds numbers investigated in this study. However, opposite trend is observed for the second smallest bubble size (0.25 mm): by increasing the Reynolds number, the fraction of the bubble decreases from 30% to 13%. Consequently, changing the Reynolds number does not have a significant impact on the average bubble size in the wake of the hydrofoil. This finding is aligned with experimental findings of Karn et al. (Karn et al., 2015a,b,c). In Fig. 15 the wake velocity distribution at different Reynolds numbers is presented. The velocity in this graph is
Fig. 10. Bubble diameter distribution at 4.65c from the hydrofoil centre of rotation under different air ventilation flow rate. 66
Ocean Engineering 172 (2019) 59–71
S. Vahaji et al.
(a) CQ
1.1e-4
(b) CQ
2.2e-4
(c) CQ
4.4e-4
Fig. 12. Void fraction under different air ventilation flow rate.
67
Ocean Engineering 172 (2019) 59–71
S. Vahaji et al.
(a) CQ
1.1e-4
(b) CQ
2.2e-4
(c) CQ
4.4e-4
Fig. 13. Local void fraction at 1.34c, 3c and 4.65c under different air ventilation flow rate.
68
Ocean Engineering 172 (2019) 59–71
S. Vahaji et al.
(a) Re=2.4e5
Fig. 14. Bubble diameter distribution at 4.65c from the hydrofoil centre of rotation under different Re
(b) Re=4.1e5
(c) Re=6.1e5
Fig. 15. Wake velocity distribution at 4.65c from the hydrofoil centre of rotation under different Re
normalised versus the free stream velocity. Therefore, all three conditions exhibit similar behaviour in the wake in terms of non-dimensional velocity distribution. The void fraction contour plots for different Reynolds numbers are presented in Fig. 16. In this figure, by increasing the Reynolds number no significant variation is observed in the void fraction distribution in the wake of the hydrofoil. In order to have a closer examination of the void fraction distribution for varying Reynolds numbers, in Fig. 17 the void fraction at 1.34c, 3c and 4.65c downstream the hydrofoil for 2.1e5, 4.1e5, 6.1e5 and 8.1e5 is displayed. The maximum void fraction in this figure slightly changes from 0.521 to 0.537 by increasing the Reynolds number from 2.1e5 to 8.1e5.
(d) Re=8.1e5 Fig. 16. Void fraction under different Reynolds numbers.
wake of NACA0015 hydrofoil is investigated and validated against the experimental data of Karn et al. (Karn et al., 2015a,b,c). A systematic numerical investigation on bubbly wake of NACA0015 hydrofoil is carried out under different conditions including angle of attack (AOA), air-entrainment coefficient, and Reynolds number to allow the study of the influence of these parameters on the bubble size distribution in the bubbly wake. The results show that AOA has significant impact on the bubble size distribution, velocity distribution and void fraction distribution in the wake. The bubble size is reduced by increasing the
6. Conclusion In this paper, numerical simulation of two-phase bubbly flow in the 69
Ocean Engineering 172 (2019) 59–71
S. Vahaji et al.
(a) Re=2.4e5
(b) Re=4.1e5
(c) Re=6.1e5
(d) Re=8.1e5 Fig. 17. Local void fraction at 1.34c, 3c and 4.65c under different Reynolds numbers.
70
Ocean Engineering 172 (2019) 59–71
S. Vahaji et al.
angle of attack due to the higher turbulence imposed by larger AOA. Change in AOA also creates a drift in the location where the minimum velocity occurs in the wake, which is an upward drift by increasing the AOA. Also, by increasing the AOA, the bubbles are dispersed more effectively within the wake. Increasing the air-entrainment coefficient results in production of smaller size bubbles due to higher collision and break-up rate of bubbles. Also, with higher air-entrainment coefficient, larger areas with higher void fraction are observed downstream the hydrofoil. Changing the Reynolds number does not have significant influence on the average bubble size in the wake of the hydrofoil. Similarly, the void fraction distribution is not drastically changed by varying the Reynolds number. However, Reynolds number has direct impact on the bubble size distribution at the wake of the hydrofoil.
bubble entrainment due to a horizontal plunging jet. Int. Shipbuild. Prog. 60 (1–4), 435–469. Ishii, M., Zuber, N., 1979. Drag coefficient and relative velocity in bubbly, droplet or particulate flows. AIChE J. 25 (5), 843–855. https://doi.org/10.1002/aic. 690250513. Karn, A., Ellis, C., Hong, J., Arndt, R.E., 2015a. Investigations into the turbulent bubbly wake of a ventilated hydrofoil: moving toward improved turbine aeration techniques. Exp. Therm. Fluid Sci. 64, 186–195. Karn, A., Ellis, C.R., Milliren, C., Hong, J., Scott, D., Arndt, R.E., Gulliver, J.S., 2015b. Bubble size characteristics in the wake of ventilated hydrofoils with two aeration configurations. International Journal of Fluid Machinery and Systems 8 (2), 73–83. Karn, A., Monson, G.M., Ellis, C.R., Hong, J., Arndt, R.E., Gulliver, J.S., 2015c. Mass transfer studies across ventilated hydrofoils: a step towards hydroturbine aeration. Int. J. Heat Mass Tran. 87, 512–520. Karn, A., Shao, S., Arndt, R.E., Hong, J., 2016. Bubble coalescence and breakup in turbulent bubbly wake of a ventilated hydrofoil. Exp. Therm. Fluid Sci. 70, 397–407. Kerdouss, F., Bannari, A., Proulx, P., 2006. CFD modeling of gas dispersion and bubble size in a double turbine stirred tank. Chem. Eng. Sci. 61 (10), 3313–3322. Krepper, E., Beyer, M., Frank, T., Lucas, D., Prasser, H.-M., 2009. CFD modelling of polydispersed bubbly two-phase flow around an obstacle. Nucl. Eng. Des. 239 (11), 2372–2381. https://doi.org/10.1016/j.nucengdes.2009.06.015. Kurul, N., Podowski, M.Z., 1990. Multi-dimensional effects in forced convection subcooled boiling. In: Ninth Heat Transfer Conference. Hemisphere Publishing Corporation, Jerusalem, Israel. Laali, A., Michel, J., 1984. Air entrainment in ventilated cavities: case of the fully developed “half-cavity”. J. Fluid Eng. 106 (3), 327–335. Lahey Jr., R.T., Drew, D.A., 2001. The analysis of two-phase flow and heat transfer using a multidimensional, four field, two-fluid model. Nucl. Eng. Des. 204 (1), 29–44. Liao, Y.X., Lucas, D., 2010. A literature review on mechanisms and models for the coalescence process of fluid particles. Chem. Eng. Sci. 65 (10), 2851–2864. https://doi. org/10.1016/j.ces.2010.02.020. Luo, H., Svendsen, H.F., 1996. Theoretical model for drop and bubble breakup in turbulent dispersions. AIChE J. 42 (5), 1225–1233. March, P.A., Brice, T.A., Mobley, M.H., Cybularz, J.M., 1992. Turbines for solving the DO dilemma. Hydro Rev. 11 (1), 30–36. Min, J., Bao, Y., Chen, L., Gao, Z., Smith, J.M., 2008. Numerical simulation of gas dispersion in an aerated stirred reactor with multiple impellers. Ind. Eng. Chem. Res. 47 (18), 7112–7117. Olmos, E., Gentric, C., Vial, C., Wild, G. and Midoux, N., 2001. Numerical simulation of multiphase flow in bubble column reactors. Influence of bubble coalescence and break-up. Chem. Eng. Sci.. 56(21–22), 6359-6365. https://doi.org/10.1016/S00092509(01)00204-4. Pohorecki, R., Moniuk, W., Bielski, P. and Zdrójkowski, A., 2001. Modelling of the coalescence/redispersion processes in bubble columns. Chem. Eng. Sci.. 56(21–22), 6157-6164. https://doi.org/10.1016/S0009-2509(01)00214-7. Prince, M.J., Blanch, H.W., 1990. Bubble coalescence and break-up in air-sparged bubble columns. AIChE J. 36 (10), 1485–1499. https://doi.org/10.1002/aic.690361004. Rotta, J., Turbulente Strömungen. 1972. Teubner, Stuttgart. Thompson, E.J., Gulliver, J.S., 1997. Oxygen transfer similitude for vented hydroturbine. J. Hydraul. Eng. 123 (6), 529–538. Tse, K., Martin, T., McFarlane, C., Nienow, A., 2003. Small bubble formation via a coalescence dependent break-up mechanism. Chem. Eng. Sci. 58 (2), 275–286. Vahaji, S., Cheung, S.C.P., Yeoh, G.H., Tu, J., 2017. Investigation of the influence of elevated pressure on subcooled boiling flow—model evaluation toward generic approach. J. Heat Tran. 139 (7), 074501. Wang, S., Lee, S., Jones, O., Lahey, R., 1987. 3-D turbulence structure and phase distribution measurements in bubbly two-phase flows. Int. J. Multiphas. Flow 13 (3), 327–343. Yeoh, G., Vahaji, S., Cheung, S., Tu, J., 2014. Modeling subcooled flow boiling in vertical channels at low pressures–Part 2: evaluation of mechanistic approach. Int. J. Heat Mass Tran. 75, 754–768. Yeoh, G.H., Tu, J.Y., 2004. Population balance modelling for bubbly flows with heat and mass transfer. Chem. Eng. Sci. 59 (15), 3125–3139. https://doi.org/10.1016/j.ces. 2004.04.023. Yeoh, G.H., Tu, J.Y., 2005. Thermal-hydrodynamic modeling of bubbly flows with heat and mass transfer. AIChE J. 51 (1), 8–27. https://doi.org/10.1002/aic.10297. Zhou, Z., Yuan, Q., Jia, X., Feng, W., Wen, J., 2013. Experimental study and CFD simulation ofMass TransferCharacteristics of a gas-induced pulsating flow BubbleColumn. Chem. Biochem. Eng. Q. 27 (2), 167–175.
Acknowledgement The authors are grateful for the support provided by Australia-China Centre for Maritime Engineering (ACCME) (Project No. ACSRF48199), National Natural Science Foundation of China (No.51436002) and the CSC Scholarship (No.201506575024). References Anglart, H. and Nylund, O., 1996. CFD application to prediction of void distribution in two-phase bubbly flows in rod bundles. Nucl. Eng. Des.. 163(1–2), 81-98. https://doi. org/10.1016/0029-5493(95)01160-9. Antal, S., Lahey Jr., R., Flaherty, J., 1991. Analysis of phase distribution and turbulence in dispersed particle/liquid flows. Chem. Eng. Commun. 174, 85–113. Chesters, A., Hofman, G., 1982. Bubble Coalescence in Pure Liquids. Mechanics and Physics of Bubbles in Liquids. Springer, pp. 353–361. Cheung, S., Vahaji, S., Yeoh, G., Tu, J., 2014. Modeling subcooled flow boiling in vertical channels at low pressures–Part 1: assessment of empirical correlations. Int. J. Heat Mass Tran. 75, 736–753. Cheung, S. C. P., Deju, L., Yeoh, G. H. and Tu, J. Y., 2013. Modeling of bubble size distribution in isothermal gas–liquid flows: numerical assessment of population balance approaches. Nucl. Eng. Des.. 265(0), 120-136. https://doi.org/10.1016/j. nucengdes.2013.08.049. Cheung, S. C. P., Duan, X., Yeoh, G. H., Tu, J., Krepper, E. and Lucas, D., 2010. Modelling of Polydispersed Flows Using Two Population Balance Approaches. 6th International Symposium on Multiphase Flow, Heat Mass Transfer and Energy Conversion. L. J. Guo, D. D. Joseph, Y. Matsumoto, M. Sommerfeld and Y. S. Wang. vol. 1207: 809814. Cheung, S. C. P., Yeoh, G. H. and Tu, J. Y., 2007a. On the modelling of population balance in isothermal vertical bubbly flows—average bubble number density approach. Chem. Eng. Process: Process Intensification. 46(8), 742-756. https://doi.org/10. 1016/j.cep.2006.10.004. Cheung, S.C.P., Yeoh, G.H., Tu, J.Y., 2007b. On the numerical study of isothermal vertical bubbly flow using two population balance approaches. Chem. Eng. Sci. 62 (17), 4659–4674. https://doi.org/10.1016/j.ces.2007.05.030. Cheung, S.C.P., Yeoh, G.H., Tu, J.Y., 2008. Population balance modeling of bubbly flows considering the hydrodynamics and thermomechanical processes. AIChE J. 54 (7), 1689–1710. https://doi.org/10.1002/aic.11503. Daskiran, C., Liu, I.-H., Oztekin, A., 2017. Computational study of multiphase flows over ventilated translating blades. Int. J. Heat Mass Tran. 110, 262–275. Deju, L., Cheung, S. C. P., Yeoh, G. H. and Tu, J. Y., 2013. Capturing coalescence and break-up processes in vertical gas–liquid flows: assessment of population balance methods. Appl. Math. Model.. 37(18–19), 8557-8577. https://doi.org/10.1016/j. apm.2013.03.063. Frank, T., Shi, J., Burns, A.D., 2004. Validation of Eulerian multiphase flow models for nuclear safety application. In: Proceeding of the Third International Symposium on Two-phase Modelling and Experimentation, Pisa, Italy. Hsiao, C.-T., Wu, X., Ma, J., Chahine, G., 2013. Numerical and experimental study of
71