# CFD Approaches for Modelling Bubble Entrainment by an Impinging Jet

Authors: Martin Schmidtke; Dirk Lucas
Journal: Science and Technology of Nuclear Installations
 (2009)

---

## Abstract

This contribution
presents different approaches for the modeling
of gas entrainment under water by a plunging
jet. Since the generation of bubbles happens on
a scale which is smaller than the bubbles, this
process cannot be resolved in meso-scale
simulations, which include the full length of
the jet and its environment. This is why the gas
entrainment has to be modeled in meso-scale
simulations. In the frame of a Euler-Euler
simulation, the local morphology of the phases
has to be considered in the drag model. For
example, the gas is a continuous phase above the
water level but bubbly below the water level.
Various drag models are tested and their
influence on the gas void fraction below the
water level is discussed.
The algebraic interface area density (AIAD) model applies a drag
coefficient for bubbles and a different drag coefficient for the
free surface. If the AIAD model is used for the simulation of
impinging jets, the gas entrainment depends on the free parameters
included in this model. The calculated gas entrainment can be
adapted via these parameters. Therefore, an advanced AIAD approach
could be used in future for the implementation of models (e.g.,
correlations) for the gas entrainment.


---

## Body

## 1. Introduction

This work concerns the evaluation of the capabilities
of the CFX-11 software for the numerical predictions of gas entrainment in the
case of a plunging jet configuration. The
configuration of an impinging jet occurs in different scenarios of reactor
safety analysis.In the scenario of an emergency core cooling (ECC), water is injected into
the cold leg. The pipe may only be partially filled with hot water, if a loss
of coolant accident occurs. In this case, the injected cold water impinges as a
jet on the surface of the hot water. Depending on the velocity of the jet,
steam bubbles may be entrained below the surface by the impinging jet. These
bubbles contribute to heat exchanged and mixing of the fluids. Heat transfer
between cold and hot water and mixing in the cold leg play an important role since the
mixed water enters  the reactor pressure vessel and may cause high temperature
gradients at the wall of the vessel. These gradients cause mechanical stress in
the wall due to thermal shock, which can have a negative effect on the
durability of the reactor vessel.An impinging jet may also occur, when an emergency coolant tank is
filled up with water and the initial water level is below the inlet. Here, the
mixing of the injected water and the water in the tank is a point of interest
if the temperatures or the boron concentrations are different.Another scenario for the occurrence of plunging jet phenomena can be
found in the case of a break, when insulation material of components is released
by the break. The fibrous material is transported into the reactor sump and
might there perturb the core cooling system. During this situation, the reactor
sump is partially filled with water. The jet from the break impinges at the
sump water surface and causes a fluid flow in the sump, which influences the
transport of the fibrous insulation material towards the sump strainers. The gas
entrainment and its influence on the fluid flow field and the transport of the
fibrous insulation are of particular interest.Generally for
the CFD modeling of large hydrodynamic configurations with multiphase flow, the
Euler-Euler approach is used. The physical process of bubble generation near a
plunging jet occurs on a very small scale, which cannot be resolved in a mesoscale
simulation. Therefore, the gas entrainment has to be physically modeled in
simulations of plunging jets. The aim of this study is to find an approach for
the simulation of plunging jet, where the gas entrainment can be deliberately
tuned to some extent (e.g., in terms adjusting free parameters), in order for a
physical model or correlation of the entrainment process to be implemented into
future simulations.In the plunging jet configuration, gas has two
different morphologies (see Figure1). The gas above the water level is a
continuous phase, whereas the gas below the water level is bubbly, that is, a
dispersed phase. The water can be regarded as a continuous phase everywhere.
For modelling this with the Euler-Euler method, two approaches are possible.Figure 1
Morphologies of the phases near an impinging jet.One can use two different phases for the two morphologies of gas. Then, water
is treated as a third phase. Gas entrainment near the jet and degassing at the
water surface has to be modelled with sources and sink terms that describe the conversion
of gas from a continuous to a dispersed (bubbly) morphology and vice versa.
This requires algorithms that identify the regions of entrainment and of
degassing.The other approach uses only two phases, one for water
and one for gas. The different morphologies of the gas then have to be
reflected by different coefficients in the closures for the momentum transfer
between the gas and water phases.The first simulations presented here are performed
with water as a continuous phase and gas as a dispersed phase. Thus, the gas is
assumed bubbly everywhere in the domain, and a constant drag coefficient is
applied. The influence of the magnitude of the drag coefficient is
investigated. Then, a more complex drag
model is tested, which take into account the different morphologies of the gas
phase.


## 2. Definition of the Test Case

### 2.1. Geometry and Mesh

A cylindrical tank with a diameter of 100 cm is filled through a nozzle. The
water level is 50 cm below the nozzle and 150 cm above the bottom of the
tank. The nozzle diameter is 19 cm (see Figure2). To reduce the costs
of computation time, only a section of five degrees is used for the simulation
to be performed as 2D axisymmetric calculation (see Figure 3).Figure 2
Geometry of a cylindrical tank filled with water through a nozzle above the free surface.Section (five degrees) of the cylindrical tank; geometry and mesh.
(a)(b)The structured mesh has 125 uniform
cells for the total height of the domain. For the radius of the water inlet,
seven uniform cells are used and 30 uniform cells for the opening (see Figure3).
The tank is quite small compared to the jet diameter and the height of
the nozzle above the water level is small enough for physicality of the result to
be influenced significantly by the walls. This disadvantage is accepted, since
we concentrate here on the gas entrainment, which takes place where the jet
hits the water surface. It can be assumed that effects far away from this area
do not influence the gas entrainment. The limitations of the geometry and the
low mesh resolution are meant to reduce the computational costs. This is
important for parametric studies. The simplicity of the geometry is accepted
here since this investigation is meant to study concepts for modelling the gas
entrainment. For some of the cases, the grid is refined by reducing the cell
size by a factor two in each dimension.


### 2.2. Fluid Properties, Initial Conditions, and Boundary Conditions and Turbulence

During the calculations, the fluids
are water for the continuous phase and gas for the dispersed bubbly phase. The
main properties (25°C and atmospheric pressure) for water and gas are
summarized in Table 1.Table 1
Fluid properties.
PropertyWaterGas: airViscosity (kg·m−1·s−1)0.00088991.831e-05Density (kg·m−3)9971.185The domain is partially filled with water up to a level of 150 cm above the
bottom. The distance between the free surface and the water nozzle is then
equal to 50 cm. The initial velocity for both water and gas in the
computational domain is taken equal to 0 m/s in each direction. The hydrostatic pressure is initialized accordingly to the water level in the domain.Liquid Inlet
The jet is injected through the
nozzle with a velocity of 3 m/s. The volume fraction is 1 for water and 0
for gas.Gas Outlet
For the gas outlet, an opening
condition is used. The volume fraction is 1 for gas and 0 for water.
Aconstant relative pressure equal to 0 Pa is assumed. For the fluid,
a velocity normal to the boundary condition is considered.Liquid Outlet
For the liquid outlet, an outlet
condition is used. The volume fraction is 1 for water and 0 for gas. Therefore,
the gas mass flow rate is equal to 0 kg/s at this boundary condition. For
the maintenance of a constant liquid level, the liquid mass flow rate leaving
the domain is defined equal to the liquid mass flow rate introduced by the
injector.Walls
Outer walls are adiabatic walls and
are defined using a no slip boundary condition. For the “inner walls” caused by
limiting the domain to a section, a symmetry boundary condition is applied. In
the case of stratified flows, the buoyancy force causes a separation of gas and
water.Turbulence Model
The homogeneous shear stress
turbulence (SST) model is applied (i.e., no separate calculation of the
turbulence for both phases). In the ANSYS CFX-Solver modelling guide [1], a
homogeneous turbulence model is recommended for separate flow and stratified
flow, whereas for dilute dispersed two-phase flow (e.g., bubbly flow), the
manual recommends using separate turbulence model for each phase. In the plunging jet, separate flow and
bubbly flow coexist in one domain, so none of the turbulence approaches is
suitable everywhere in the domain. The calculations presented below are
calculated with a homogeneous SST model by default.  For comparison, some calculations are repeated
with an inhomogeneous turbulence model, which is the SST model for the liquid
phase and a laminar assumption for the gaseous phase (see Section 3.3).


## 2.1. Geometry and Mesh

A cylindrical tank with a diameter of 100 cm is filled through a nozzle. The
water level is 50 cm below the nozzle and 150 cm above the bottom of the
tank. The nozzle diameter is 19 cm (see Figure2). To reduce the costs
of computation time, only a section of five degrees is used for the simulation
to be performed as 2D axisymmetric calculation (see Figure 3).Figure 2
Geometry of a cylindrical tank filled with water through a nozzle above the free surface.Section (five degrees) of the cylindrical tank; geometry and mesh.
(a)(b)The structured mesh has 125 uniform
cells for the total height of the domain. For the radius of the water inlet,
seven uniform cells are used and 30 uniform cells for the opening (see Figure3).
The tank is quite small compared to the jet diameter and the height of
the nozzle above the water level is small enough for physicality of the result to
be influenced significantly by the walls. This disadvantage is accepted, since
we concentrate here on the gas entrainment, which takes place where the jet
hits the water surface. It can be assumed that effects far away from this area
do not influence the gas entrainment. The limitations of the geometry and the
low mesh resolution are meant to reduce the computational costs. This is
important for parametric studies. The simplicity of the geometry is accepted
here since this investigation is meant to study concepts for modelling the gas
entrainment. For some of the cases, the grid is refined by reducing the cell
size by a factor two in each dimension.


## 2.2. Fluid Properties, Initial Conditions, and Boundary Conditions and Turbulence

During the calculations, the fluids
are water for the continuous phase and gas for the dispersed bubbly phase. The
main properties (25°C and atmospheric pressure) for water and gas are
summarized in Table 1.Table 1
Fluid properties.
PropertyWaterGas: airViscosity (kg·m−1·s−1)0.00088991.831e-05Density (kg·m−3)9971.185The domain is partially filled with water up to a level of 150 cm above the
bottom. The distance between the free surface and the water nozzle is then
equal to 50 cm. The initial velocity for both water and gas in the
computational domain is taken equal to 0 m/s in each direction. The hydrostatic pressure is initialized accordingly to the water level in the domain.Liquid Inlet
The jet is injected through the
nozzle with a velocity of 3 m/s. The volume fraction is 1 for water and 0
for gas.Gas Outlet
For the gas outlet, an opening
condition is used. The volume fraction is 1 for gas and 0 for water.
Aconstant relative pressure equal to 0 Pa is assumed. For the fluid,
a velocity normal to the boundary condition is considered.Liquid Outlet
For the liquid outlet, an outlet
condition is used. The volume fraction is 1 for water and 0 for gas. Therefore,
the gas mass flow rate is equal to 0 kg/s at this boundary condition. For
the maintenance of a constant liquid level, the liquid mass flow rate leaving
the domain is defined equal to the liquid mass flow rate introduced by the
injector.Walls
Outer walls are adiabatic walls and
are defined using a no slip boundary condition. For the “inner walls” caused by
limiting the domain to a section, a symmetry boundary condition is applied. In
the case of stratified flows, the buoyancy force causes a separation of gas and
water.Turbulence Model
The homogeneous shear stress
turbulence (SST) model is applied (i.e., no separate calculation of the
turbulence for both phases). In the ANSYS CFX-Solver modelling guide [1], a
homogeneous turbulence model is recommended for separate flow and stratified
flow, whereas for dilute dispersed two-phase flow (e.g., bubbly flow), the
manual recommends using separate turbulence model for each phase. In the plunging jet, separate flow and
bubbly flow coexist in one domain, so none of the turbulence approaches is
suitable everywhere in the domain. The calculations presented below are
calculated with a homogeneous SST model by default.  For comparison, some calculations are repeated
with an inhomogeneous turbulence model, which is the SST model for the liquid
phase and a laminar assumption for the gaseous phase (see Section 3.3).


## 3. Dispersed Phase Model for Gas

The simplest approach for modelling the plunging jet is achieved if the water is treated as continuous  phase and
the gas is a dispersed phase with a constant particle diameter (d=2 mm). This
approach neglects the fact that the particle model is not appropriate for the gas
above the water level.


### 3.1. Drag Model

For the bubbles, a constant drag coefficient is used. The default value used here isCD=0.44, which is the drag coefficient for solid spheres in the Newton
range. To study the effect of the particle drag coefficient on the gas entrainment, the simulations are performed with the drag coefficient CD=0.44 and with a reduced value CD=0.05
for comparison.


### 3.2. Nondrag Forces

Since the focus of this study is on the gas entrainment at the surface, nondrag
forces (lift force and turbulent dispersion force) are neglected here. It is
expected that the turbulent dispersion force causes an increase of the horizontal
extension of the bubble plume. Nevertheless, an application of nondrag forces
above the water level is meaningless. Therefore, nondrag forces are not
modelled here.


### 3.3. Results for the Simulations with Gas as  Dispersed Phase

A few seconds after the jet release from the nozzle the interface becomes stable
and the gas void fraction field also becomes steady (see Figure4). A
reduction of the drag coefficient from CD=0.44 to CD=0.05 has no significant effect on the gas void fraction field. Thus, the drag
coefficient cannot be used as a parameter that influences the gas entrainment in the
simulation. If the SST model is applied for the liquid phase and the turbulence
of the gas is neglected (laminar assumption), the gas void fractions are
similar to those in Figure 4 which have been calculated with a homogenous SST
model. Therefore, the coupling of both phases by sharing the same turbulence
field does not contribute to gas entrainment. In the subsequent simulations,
only the homogenous SST model is used.Gas void fraction fields for a simulation using a dispersed phase model
for gas homogenous SST turbulence for both phases.
(a)(b)


### 3.4. Vertical Gas Fluxes below the Water Level

For a characterization of the gas
entrainment, it is advantageous to use integral quantities for the intensity
and the geometry of the gas plume. By performing an extensive survey of
experimental studies Bin [2] obtained correlations for the penetration depth
and the entrainment rate. The penetration depth hp is the vertical extension of the gas plume below the
water level. Bin's correlation [2] for the penetration depth in meter is(1)hp=2.1wj0.775d00.67, where d0
is the nozzle diameter
in meter and wj
is the
vertical jet velocity at the water level (in m/s). Due to gravitational
acceleration, the velocity of a free falling jet is increasing until it hits the surface. If hj
is the height of the nozzle above the surface and w0
is the liquid velocity at the nozzle, wj can be calculated as(2)wj=w02+2ghj. For
the height of hj=0.5 and w0=3.0 m/s, one obtains wj=4.3 m/s for the jet
velocity at the water level and a penetration depth of 215 cm according to (1). The predicted
value for the penetration depth is larger then the depth of the water in the
tank. Therefore, the length of the gas plume might be restricted artificially
by the geometry. In fact, according to Figure 4 the gas plumes almost reach
the bottom of the tank. For a better quantification of the vertical
distribution of the gas, the gas void fraction αG is integrated on horizontal planes:(3)AG=∫AαGdA. Since the gas void fraction is dimensionless, the integral (3) yields the dimension
of an area for AG. This
can be interpreted as the area occupied by gas on the horizontal plane A. In
Figure 5, AG is plotted versus the depth
below the water level. Here, AG is normalized by the area of the jet cross-section at the inlet. There is only
a little difference between the values for CD=0.44 and CD=0.05
(see also Figure 4).Figure 5
Normalized gas area as a function of the depth below the water surface.The depth at which the normalized gas area is zero can be used to define a
penetration depth for jets. According to Figure5, the penetration depth is 150 m which means that the plumes reach the tank bottom. This is in accordance with
the prediction of (1).The entrainment rate is the ratio of the gas fluxQG entrained below the
water by the impinging jet and the water flux QL of the jet. The correlation for the entrainment rate
suggested by Bin [2] is(4)QGQL=0.04Fr0.28(hjd0)0.4,withFr=wj2gd0, where hj is the jet height above the water level and g is the gravity.
For the boundary conditions used in the simulations, this correlation yields an
entrainment rate of 0.08. Another correlation was obtained by Ohkawa et al. [3]:(5)QGQL=0.016[Fr0.28(hjd0)0.4]1.17. This correlation yields an entrainment rate of 0.036 for the boundary
conditions used here.To compare the simulation results in terms of entrainment rate with the predictions given
by correlations (4) and (5), the gas fluxes below the water have to be
investigated more closely. The product of the gas void fraction αG and the vertical velocity of the gas wG
defines a vertical gas flux density qG:(6)qG=αGwG. The upward and the downward fluxes can be distinguished by the definition of(7)qG+={qG,ifqG>0,0,else,qG−={qG,ifqG<0,0,else. So the total downward flux at a certain level below the surface is(8)QG−=∫AqG−dA, where A is the horizontal cross-section of the
domain at a certain level below the surface. The total upward flux QG+ is calculated in the same way. Figures 6–8 show the gas
void faction, the vertical velocities, and the vertical gas flux density at a
depth of 30 cm below the water level.Figure 6
Gas void  fraction as a function of the horizontal distance to the jet axis, 30 cm below the
water level. The values forCD=0.44 (blue) and CD=0.95 (mangenta)
differ hardly.Figure 7
Vertical velocities as a function of the horizontal distance to the jet axis,
30 cm below the water level (CD=0.44).Figure 8
Vertical gas flux density as a function of the horizontal distance to the jet
axis, 30 cm below the water level (CD=0.44).In Figure9,
the total upward and downward gas fluxes are shown for the two jets modelled
with the drag coefficients CD=0.44 and CD=0.05. The gas fluxes
below the water level are normalized by the water flux of the jet QL,0
at the nozzle and
plotted as function of the depth below the water level. For CD=0.44,
the upward and downward gas fluxes are similar which means that the solution is
steady. At each level, the same amount is transported upwards and downwards.
For the different drag coefficients, the gas downward fluxes are also similar. Thus,
the drag hardly contributes to gas entrainment. The carry-under of gas therefore seems to be mainly caused
by numerical effects within the solver. The curves for the normalized gas
fluxes show a local minimum ca. 10 cm below the water level. This is the depth
of the deformed water surface (“trumpet”) near the jet. At a depth of 60 cm,
all the curves have a local maximum. This can be explained by the
re-entrainment of bubbles, which are trapped in the vortex caused by the jet.
At the depth of 20 cm, the normalized gas fluxes are about 0.06 which is just
between the predictions for the entrainment rate by Bin (see (4)) and by
Ohkawa (see (5)). Since the entrainment is mainly caused by numerical effects
in this setup, we can expect the value to be sensitive to the geometry and
resolution of the grid. This must be studied in future. However, it seems to be
a coincidence that the simulated entrainment rate in this simulation is in the
range predicted by empirical correlations.Figure 9
Normalized gas
fluxes as function of the depth below the water level. Blue symbols: upward
gas flux (CD=0.44).
Red symbols: downward gas flux (CD=0.44).
Yellow: downward gas flux (CD=0.05).Figure 10
The continuous phase density as a function of the gas void
fraction.Figure 11
The area density as a
function of the gas void fraction.


## 3.1. Drag Model

For the bubbles, a constant drag coefficient is used. The default value used here isCD=0.44, which is the drag coefficient for solid spheres in the Newton
range. To study the effect of the particle drag coefficient on the gas entrainment, the simulations are performed with the drag coefficient CD=0.44 and with a reduced value CD=0.05
for comparison.


## 3.2. Nondrag Forces

Since the focus of this study is on the gas entrainment at the surface, nondrag
forces (lift force and turbulent dispersion force) are neglected here. It is
expected that the turbulent dispersion force causes an increase of the horizontal
extension of the bubble plume. Nevertheless, an application of nondrag forces
above the water level is meaningless. Therefore, nondrag forces are not
modelled here.


## 3.3. Results for the Simulations with Gas as  Dispersed Phase

A few seconds after the jet release from the nozzle the interface becomes stable
and the gas void fraction field also becomes steady (see Figure4). A
reduction of the drag coefficient from CD=0.44 to CD=0.05 has no significant effect on the gas void fraction field. Thus, the drag
coefficient cannot be used as a parameter that influences the gas entrainment in the
simulation. If the SST model is applied for the liquid phase and the turbulence
of the gas is neglected (laminar assumption), the gas void fractions are
similar to those in Figure 4 which have been calculated with a homogenous SST
model. Therefore, the coupling of both phases by sharing the same turbulence
field does not contribute to gas entrainment. In the subsequent simulations,
only the homogenous SST model is used.Gas void fraction fields for a simulation using a dispersed phase model
for gas homogenous SST turbulence for both phases.
(a)(b)


## 3.4. Vertical Gas Fluxes below the Water Level

For a characterization of the gas
entrainment, it is advantageous to use integral quantities for the intensity
and the geometry of the gas plume. By performing an extensive survey of
experimental studies Bin [2] obtained correlations for the penetration depth
and the entrainment rate. The penetration depth hp is the vertical extension of the gas plume below the
water level. Bin's correlation [2] for the penetration depth in meter is(1)hp=2.1wj0.775d00.67, where d0
is the nozzle diameter
in meter and wj
is the
vertical jet velocity at the water level (in m/s). Due to gravitational
acceleration, the velocity of a free falling jet is increasing until it hits the surface. If hj
is the height of the nozzle above the surface and w0
is the liquid velocity at the nozzle, wj can be calculated as(2)wj=w02+2ghj. For
the height of hj=0.5 and w0=3.0 m/s, one obtains wj=4.3 m/s for the jet
velocity at the water level and a penetration depth of 215 cm according to (1). The predicted
value for the penetration depth is larger then the depth of the water in the
tank. Therefore, the length of the gas plume might be restricted artificially
by the geometry. In fact, according to Figure 4 the gas plumes almost reach
the bottom of the tank. For a better quantification of the vertical
distribution of the gas, the gas void fraction αG is integrated on horizontal planes:(3)AG=∫AαGdA. Since the gas void fraction is dimensionless, the integral (3) yields the dimension
of an area for AG. This
can be interpreted as the area occupied by gas on the horizontal plane A. In
Figure 5, AG is plotted versus the depth
below the water level. Here, AG is normalized by the area of the jet cross-section at the inlet. There is only
a little difference between the values for CD=0.44 and CD=0.05
(see also Figure 4).Figure 5
Normalized gas area as a function of the depth below the water surface.The depth at which the normalized gas area is zero can be used to define a
penetration depth for jets. According to Figure5, the penetration depth is 150 m which means that the plumes reach the tank bottom. This is in accordance with
the prediction of (1).The entrainment rate is the ratio of the gas fluxQG entrained below the
water by the impinging jet and the water flux QL of the jet. The correlation for the entrainment rate
suggested by Bin [2] is(4)QGQL=0.04Fr0.28(hjd0)0.4,withFr=wj2gd0, where hj is the jet height above the water level and g is the gravity.
For the boundary conditions used in the simulations, this correlation yields an
entrainment rate of 0.08. Another correlation was obtained by Ohkawa et al. [3]:(5)QGQL=0.016[Fr0.28(hjd0)0.4]1.17. This correlation yields an entrainment rate of 0.036 for the boundary
conditions used here.To compare the simulation results in terms of entrainment rate with the predictions given
by correlations (4) and (5), the gas fluxes below the water have to be
investigated more closely. The product of the gas void fraction αG and the vertical velocity of the gas wG
defines a vertical gas flux density qG:(6)qG=αGwG. The upward and the downward fluxes can be distinguished by the definition of(7)qG+={qG,ifqG>0,0,else,qG−={qG,ifqG<0,0,else. So the total downward flux at a certain level below the surface is(8)QG−=∫AqG−dA, where A is the horizontal cross-section of the
domain at a certain level below the surface. The total upward flux QG+ is calculated in the same way. Figures 6–8 show the gas
void faction, the vertical velocities, and the vertical gas flux density at a
depth of 30 cm below the water level.Figure 6
Gas void  fraction as a function of the horizontal distance to the jet axis, 30 cm below the
water level. The values forCD=0.44 (blue) and CD=0.95 (mangenta)
differ hardly.Figure 7
Vertical velocities as a function of the horizontal distance to the jet axis,
30 cm below the water level (CD=0.44).Figure 8
Vertical gas flux density as a function of the horizontal distance to the jet
axis, 30 cm below the water level (CD=0.44).In Figure9,
the total upward and downward gas fluxes are shown for the two jets modelled
with the drag coefficients CD=0.44 and CD=0.05. The gas fluxes
below the water level are normalized by the water flux of the jet QL,0
at the nozzle and
plotted as function of the depth below the water level. For CD=0.44,
the upward and downward gas fluxes are similar which means that the solution is
steady. At each level, the same amount is transported upwards and downwards.
For the different drag coefficients, the gas downward fluxes are also similar. Thus,
the drag hardly contributes to gas entrainment. The carry-under of gas therefore seems to be mainly caused
by numerical effects within the solver. The curves for the normalized gas
fluxes show a local minimum ca. 10 cm below the water level. This is the depth
of the deformed water surface (“trumpet”) near the jet. At a depth of 60 cm,
all the curves have a local maximum. This can be explained by the
re-entrainment of bubbles, which are trapped in the vortex caused by the jet.
At the depth of 20 cm, the normalized gas fluxes are about 0.06 which is just
between the predictions for the entrainment rate by Bin (see (4)) and by
Ohkawa (see (5)). Since the entrainment is mainly caused by numerical effects
in this setup, we can expect the value to be sensitive to the geometry and
resolution of the grid. This must be studied in future. However, it seems to be
a coincidence that the simulated entrainment rate in this simulation is in the
range predicted by empirical correlations.Figure 9
Normalized gas
fluxes as function of the depth below the water level. Blue symbols: upward
gas flux (CD=0.44).
Red symbols: downward gas flux (CD=0.44).
Yellow: downward gas flux (CD=0.05).Figure 10
The continuous phase density as a function of the gas void
fraction.Figure 11
The area density as a
function of the gas void fraction.


## 4. The Simmer-Drag Model

In the previous simulations, the gas was treated as a
dispersed phase everywhere in the domain. However, the SIMMER model, first introduced into the SIMMER-code
[4], takes into account the distinction in morphology that phases can have in the domain. The
morphology of the phases has to be reflected by appropriate parameters in the
drag force. The magnitude of the force density for the drag is(9)|D|=CDa12ρ|VR|2, where CD
is the drag
coefficient, a is the interfacial area density, and ρ is the density of the continuous phase (if the
other phase is a dispersed phase). VR is the relative velocity between the two phases.In the SIMMER model, the drag force depends on the gas void fractionαG.
The gas is assumed to have the morphology of bubbles where the gas void
fraction is low, that is, αG<0.3. Where the gas void fraction is high (αG>0.7), droplets are supposed to be present in gas. In the intermediate
range (0.3<αG<0.7), a linear interpolation between bubble drag and droplet drag is
performed. This means that the interfacial area density a and the continuous phase density ρ depend on the gas void fraction, which is used
as indicator for the morphology of both phases.


### 4.1. The Continuous Phase Density

If the gas void fraction is low, the liquid phase is the continuous phase (ρ=ρL). For high gas void fractions, the gas
is the continuous phase (ρ=ρG). In the intermediate range, the density ρ is interpolated:(10)ρ={ρL,forαG≤0.3,ρG,forαG≥0.7,ρL+(ρG−ρL)αG−0.30.7−0.3,else.


### 4.2. The Area Density

The total area density for a spherical particle is(11)aP=6αPdP, where αP
is the particle void fraction. The drag coefficients for particles are related
to projected areas. The projected area of a sphere is 1/4th of its total area. Therefore,
the so-called projected area
densities for spherical bubbles and droplets are calculated as(12)aB=32αGdB,aD=32αLdD, where dB
und dD are the bubble
diameter and the droplet diameter, respectively. In the simulations for
simplicity, the same particle diameters are applied for droplets (d=dB=dD). Similar to the continuous
phase density, the global area density a is defined as(13)a={aB,forαG≤0.3,aD,forαG≥0.7,aB*+(aD*−aD*)αG−0.30.7−0.3,else,whereaB*=aB(αG=0.3),aD*=aD(αG=0.7).


### 4.3. The Drag Coefficient

Bubbles are assumed spherical, where the drag coefficient that has a constant value of 0.44
is applied. As an alternative, the Schiller-Naumann drag correlation is used
which reads(14)CD=24Re(1+0.15Re0.687). Since the
material properties of the continuous phase are included in the Reynolds
number, this equation yields two drag coefficients: CD,B for bubbles and CD,D for droplets. The drag coefficient CD
at a position in the domain is calculated according to the local gas void
fraction in the same way as it is done for the continuous phase density and the
area density:(15)CD={CD,B,forαG≤0.3,CD,D,forαG≥0.7,CD,B+(CD,D−CD,B)αG−0.30.7−0.3,else. If a constant drag is used, a case differentiation is not necessary. Then CD=0.44 is applied everywhere.


### 4.4. Results

The development of the gas void fraction near the jet is studied using transient
calculations. The drag is modified by applying either a constant drag
coefficient or the Schiller-Naumann drag correlation. The calculations show a
steady behaviour after a few seconds (see Figure12).Gas void fraction for a constant drag (CD=0.44) and particle diameter d=2 mm.
(a)(b)(c)(d)The influence of the drag model (constant drag versus Schiller-Naumann drag) and of
the particle diameter is very low. The gas entrainment seems to be always overestimated,
since gas void fractions higher than 60% occur below the surface in all
simulations.There is no free parameter inside the SIMMER drag model, which could be used to
adjust the entrainment according to an empirical correlation or another
physical entrainment model. The effect of modified drag coefficients has not
been studied yet. However, using arbitrary drag coefficients causes unphysical
velocities for buoyant particles (e.g., bubbles) and it is therefore meaningless.


## 4.1. The Continuous Phase Density

If the gas void fraction is low, the liquid phase is the continuous phase (ρ=ρL). For high gas void fractions, the gas
is the continuous phase (ρ=ρG). In the intermediate range, the density ρ is interpolated:(10)ρ={ρL,forαG≤0.3,ρG,forαG≥0.7,ρL+(ρG−ρL)αG−0.30.7−0.3,else.


## 4.2. The Area Density

The total area density for a spherical particle is(11)aP=6αPdP, where αP
is the particle void fraction. The drag coefficients for particles are related
to projected areas. The projected area of a sphere is 1/4th of its total area. Therefore,
the so-called projected area
densities for spherical bubbles and droplets are calculated as(12)aB=32αGdB,aD=32αLdD, where dB
und dD are the bubble
diameter and the droplet diameter, respectively. In the simulations for
simplicity, the same particle diameters are applied for droplets (d=dB=dD). Similar to the continuous
phase density, the global area density a is defined as(13)a={aB,forαG≤0.3,aD,forαG≥0.7,aB*+(aD*−aD*)αG−0.30.7−0.3,else,whereaB*=aB(αG=0.3),aD*=aD(αG=0.7).


## 4.3. The Drag Coefficient

Bubbles are assumed spherical, where the drag coefficient that has a constant value of 0.44
is applied. As an alternative, the Schiller-Naumann drag correlation is used
which reads(14)CD=24Re(1+0.15Re0.687). Since the
material properties of the continuous phase are included in the Reynolds
number, this equation yields two drag coefficients: CD,B for bubbles and CD,D for droplets. The drag coefficient CD
at a position in the domain is calculated according to the local gas void
fraction in the same way as it is done for the continuous phase density and the
area density:(15)CD={CD,B,forαG≤0.3,CD,D,forαG≥0.7,CD,B+(CD,D−CD,B)αG−0.30.7−0.3,else. If a constant drag is used, a case differentiation is not necessary. Then CD=0.44 is applied everywhere.


## 4.4. Results

The development of the gas void fraction near the jet is studied using transient
calculations. The drag is modified by applying either a constant drag
coefficient or the Schiller-Naumann drag correlation. The calculations show a
steady behaviour after a few seconds (see Figure12).Gas void fraction for a constant drag (CD=0.44) and particle diameter d=2 mm.
(a)(b)(c)(d)The influence of the drag model (constant drag versus Schiller-Naumann drag) and of
the particle diameter is very low. The gas entrainment seems to be always overestimated,
since gas void fractions higher than 60% occur below the surface in all
simulations.There is no free parameter inside the SIMMER drag model, which could be used to
adjust the entrainment according to an empirical correlation or another
physical entrainment model. The effect of modified drag coefficients has not
been studied yet. However, using arbitrary drag coefficients causes unphysical
velocities for buoyant particles (e.g., bubbles) and it is therefore meaningless.


## 5. The Algebraic Interfacial Area  Density (AIAD) Model

### 5.1. Drag Model

The algebraic interfacial area
density model applies two different drag coefficients,CD,B for
bubbles and CD,S for free surface. The interfacial area density a also depends on the morphology of the phases. For bubbles, the projected
interfacial area density is(16)aB=146dBαG, where dB is the bubble diameter
and αG is the gas void fraction. For a free surface, the interfacial area
density is(17)aS=|∇αG|. Since the concept of a
continuous phase is not meaningful in the range of medium gas void fractions, instead
of a continuous phase density, an average density is applied in  (9). The average density is defined as(18)ρ=ρGαG+ρL(1−αG), where ρL and ρG are the liquid and the gas phase densities, respectively.
In the bubbly regime, where αG is low, the average density according to (18) is close to the liquid phase density ρL, which is the continuous phase
density in this case. According to the flow regime (bubbly flow or stratified
flow with a free surface), the corresponding drag coefficients and interfacial
area densities have to be applied. This can be done by introducing a blending
function f which is 1 for bubbly flow
and 0 for stratified flow. Then, the area density and the drag coefficient are
well defined everywhere in the domain by(19)a=faB+(1−f)aS,(20)CD=fCD,B+(1−f)CD,S. It is not easy to find an algorithm
that recognizes the flow regime of course. A very simple approach identifies
the flow regime by using a gas void fraction limit α0. Bubbly flow is assumed, where αG<α0, and stratified flow everywhere
else. This would mean that blending function f is a step function. To
avoid numerical problems, a continuous blending function is preferred (see Figure 13):(21)f=11+exp(−100(αG−α0)). For a first judgement, the gas
entrainment is quantified by the gas void fractions just below the liquid interface.
These are investigated for various values of the free surface drag coefficient CD,S and gas void fraction
limits α0. For the bubble drag coefficient, a constant value of CD,B=0.44 is taken, based
on the drag of rigid spheres at the medium to high Reynolds number regime. As
bubble diameter dB=2 mm is
chosen.Figure 13
Blending functionf according to (21) for α0=0.1 (black curve) and α0=0.3 (red).


### 5.2. Variation of the Surface Drag Coefficient

It is not clear which surface drag coefficient is appropriate for the situation of the
impinging jet. The value ofCD,S has to include subgrid information of the free surface structure (“rough” or
“smooth”), and this certainly depends on the grid resolution, since with a
finer mesh more details of the surface structure are resolved. Therefore, the
free surface drag coefficient CD,S is varied over several orders of magnitude. Its influence on the gas void
fraction below the water surface is studied while keeping the gas void fraction
limit constant at α0=0.1. Note that the vertical water velocity at the nozzle is kept
constant at w0=3 m/s. The simulation is performed in the transient mode, but the result is almost in
a steady state 10 seconds after the start when the jet is released from the
nozzle.Figure14
shows the gas void fraction for CD,S=10. The gas entrainment seems to be overestimated here, since even at a depth
of 50 cm below the water surface gas void fractions of 0.5 appear. In Figure 15(a),
the corresponding bubble area density aB is displayed. Note that aB is proportional to the gas void fraction in (16), where it is greatest at αG=1. Of course bubbles are not assumed to be present where αG=1. According to (19) and (20), the blending
function switches to free surface area density aS
at high gas void fractions. The free surface area
density aS for this case
is shown in Figure 15(b). In Figure 16, the total area density a according to (19), the total drag
coefficient CD according
to (20), and the product of a and CD are shown.Figure 14
Gas void fractionαG for the free surface drag
coefficient CD=10.(a) Bubble projected area densityaB. (b) Free surface area density aS.
(a)(b)(a) Total area densitya. (b) Total drag coefficient CD.  (c) Product aCD.
(a)(b)(c)As one can see from Figure17 that the gas entrainment below
the surface decreases if the surface drag coefficient is reduced. Note that the
solver does not converge when CD,S=0. Since the maximal gas void fraction below the water surface is similar for CD,S=0.1 and CD,S=0.01, it is obvious
that entrainment cannot be suppressed by further reducing CD,S. Even for a very small bubble drag coefficient (CD,B=0.044), the gas
entrainment is not negligible (see Figure 18). This indicates that numerical
diffusion contributes to entrainment.Gas void fraction for various surface drag
coefficients. Representative plots at timet>10 seconds.
(a)(b)(c)(d)Figure 18
Gas void fraction for a very small bubble drag
coefficient.


### 5.3. Variation of the Blending Function

By changing the gas void fraction limitα0 the blending function
can be modified. The value of α0 has a significant
influence on the gas entrainment (see Figure 19). It is not clear which value
is appropriate since α0 has no physical
meaning. The definition of the blending function in general and of α0 in particular is
arbitrary to some extent. Note that in this case, the gas void fraction αG is used as criterion
to identify the location of the surface. This is a quite simple approach of
course and a more sophisticated blending function could use the gradient of αG to identify the
surface since this gradient is high near the surface.Gas void fraction forCD,S=0.1 and various α0.
(a)(b)


### 5.4. Grid Resolution

To check the influence of the grid resolution on the
numerical solution, one calculation (CD=0.1) is repeated with a
doubled spatial resolution. The gas void fraction fields are similar (see
Figure 20).Gas void fraction for two different grid
resolutions,CD=0.1, representative plots.
(a)(b)Figure20 shows that the gas plume is narrower in the
calculation with the higher resolution. By an integration of the vertical flux
density at this level according to (8) and by normalizing the result with the
water flux at the nozzle QL,0,
the dimensionless entrainment rate is obtained. For the coarse mesh, the
entrainment rate is 3.5% and for the fine mesh it is about 6.4%. Of course the
results of a CFD model should not depend on the grid resolution (see [5]). However, in a simulation of an impinging jet with increasing resolution more details of the complex surface geometry at the impinging zone is resolved. In the borderline case of an infinite resolution, the real
bubble generation process could be captured. With a decreasing resolution, the
geometry of the impinging zone is further simplified. This is the reason why
the resolution has an effect on the simulated entrainment.


## 5.1. Drag Model

The algebraic interfacial area
density model applies two different drag coefficients,CD,B for
bubbles and CD,S for free surface. The interfacial area density a also depends on the morphology of the phases. For bubbles, the projected
interfacial area density is(16)aB=146dBαG, where dB is the bubble diameter
and αG is the gas void fraction. For a free surface, the interfacial area
density is(17)aS=|∇αG|. Since the concept of a
continuous phase is not meaningful in the range of medium gas void fractions, instead
of a continuous phase density, an average density is applied in  (9). The average density is defined as(18)ρ=ρGαG+ρL(1−αG), where ρL and ρG are the liquid and the gas phase densities, respectively.
In the bubbly regime, where αG is low, the average density according to (18) is close to the liquid phase density ρL, which is the continuous phase
density in this case. According to the flow regime (bubbly flow or stratified
flow with a free surface), the corresponding drag coefficients and interfacial
area densities have to be applied. This can be done by introducing a blending
function f which is 1 for bubbly flow
and 0 for stratified flow. Then, the area density and the drag coefficient are
well defined everywhere in the domain by(19)a=faB+(1−f)aS,(20)CD=fCD,B+(1−f)CD,S. It is not easy to find an algorithm
that recognizes the flow regime of course. A very simple approach identifies
the flow regime by using a gas void fraction limit α0. Bubbly flow is assumed, where αG<α0, and stratified flow everywhere
else. This would mean that blending function f is a step function. To
avoid numerical problems, a continuous blending function is preferred (see Figure 13):(21)f=11+exp(−100(αG−α0)). For a first judgement, the gas
entrainment is quantified by the gas void fractions just below the liquid interface.
These are investigated for various values of the free surface drag coefficient CD,S and gas void fraction
limits α0. For the bubble drag coefficient, a constant value of CD,B=0.44 is taken, based
on the drag of rigid spheres at the medium to high Reynolds number regime. As
bubble diameter dB=2 mm is
chosen.Figure 13
Blending functionf according to (21) for α0=0.1 (black curve) and α0=0.3 (red).


## 5.2. Variation of the Surface Drag Coefficient

It is not clear which surface drag coefficient is appropriate for the situation of the
impinging jet. The value ofCD,S has to include subgrid information of the free surface structure (“rough” or
“smooth”), and this certainly depends on the grid resolution, since with a
finer mesh more details of the surface structure are resolved. Therefore, the
free surface drag coefficient CD,S is varied over several orders of magnitude. Its influence on the gas void
fraction below the water surface is studied while keeping the gas void fraction
limit constant at α0=0.1. Note that the vertical water velocity at the nozzle is kept
constant at w0=3 m/s. The simulation is performed in the transient mode, but the result is almost in
a steady state 10 seconds after the start when the jet is released from the
nozzle.Figure14
shows the gas void fraction for CD,S=10. The gas entrainment seems to be overestimated here, since even at a depth
of 50 cm below the water surface gas void fractions of 0.5 appear. In Figure 15(a),
the corresponding bubble area density aB is displayed. Note that aB is proportional to the gas void fraction in (16), where it is greatest at αG=1. Of course bubbles are not assumed to be present where αG=1. According to (19) and (20), the blending
function switches to free surface area density aS
at high gas void fractions. The free surface area
density aS for this case
is shown in Figure 15(b). In Figure 16, the total area density a according to (19), the total drag
coefficient CD according
to (20), and the product of a and CD are shown.Figure 14
Gas void fractionαG for the free surface drag
coefficient CD=10.(a) Bubble projected area densityaB. (b) Free surface area density aS.
(a)(b)(a) Total area densitya. (b) Total drag coefficient CD.  (c) Product aCD.
(a)(b)(c)As one can see from Figure17 that the gas entrainment below
the surface decreases if the surface drag coefficient is reduced. Note that the
solver does not converge when CD,S=0. Since the maximal gas void fraction below the water surface is similar for CD,S=0.1 and CD,S=0.01, it is obvious
that entrainment cannot be suppressed by further reducing CD,S. Even for a very small bubble drag coefficient (CD,B=0.044), the gas
entrainment is not negligible (see Figure 18). This indicates that numerical
diffusion contributes to entrainment.Gas void fraction for various surface drag
coefficients. Representative plots at timet>10 seconds.
(a)(b)(c)(d)Figure 18
Gas void fraction for a very small bubble drag
coefficient.


## 5.3. Variation of the Blending Function

By changing the gas void fraction limitα0 the blending function
can be modified. The value of α0 has a significant
influence on the gas entrainment (see Figure 19). It is not clear which value
is appropriate since α0 has no physical
meaning. The definition of the blending function in general and of α0 in particular is
arbitrary to some extent. Note that in this case, the gas void fraction αG is used as criterion
to identify the location of the surface. This is a quite simple approach of
course and a more sophisticated blending function could use the gradient of αG to identify the
surface since this gradient is high near the surface.Gas void fraction forCD,S=0.1 and various α0.
(a)(b)


## 5.4. Grid Resolution

To check the influence of the grid resolution on the
numerical solution, one calculation (CD=0.1) is repeated with a
doubled spatial resolution. The gas void fraction fields are similar (see
Figure 20).Gas void fraction for two different grid
resolutions,CD=0.1, representative plots.
(a)(b)Figure20 shows that the gas plume is narrower in the
calculation with the higher resolution. By an integration of the vertical flux
density at this level according to (8) and by normalizing the result with the
water flux at the nozzle QL,0,
the dimensionless entrainment rate is obtained. For the coarse mesh, the
entrainment rate is 3.5% and for the fine mesh it is about 6.4%. Of course the
results of a CFD model should not depend on the grid resolution (see [5]). However, in a simulation of an impinging jet with increasing resolution more details of the complex surface geometry at the impinging zone is resolved. In the borderline case of an infinite resolution, the real
bubble generation process could be captured. With a decreasing resolution, the
geometry of the impinging zone is further simplified. This is the reason why
the resolution has an effect on the simulated entrainment.


## 6. Conclusion

Generally for the CFD modelling of large hydrodynamic configurations with multiphase flow,
the Euler-Euler approach is used. This is the reason why the capabilities of
the Euler-Euler approach for the modelling of the impinging jet are
investigated in this contribution. The physical process of bubble generation
near the jet occurs on a very small scale, which cannot be resolved in a large-scale
simulation. Therefore, the gas entrainment has to be described by a model,
which represents the physics of the entrainment (e.g., a correlation). For the
implementation of such a physical model in the frame of a CFD code, a mechanism
is required that allows the adjustment of the gas entrained in the simulation
according to the correlation. Since the gas and liquid phases tend to separate due to the buoyancy force, it is an obvious choice to use its counterpart—the drag force—to obtain gas entrainment below the surface.
However, if the gas is modelled as dispersed phase in an Euler-Euler
simulation, the entrainment barely depends on the magnitude of drag
coefficient, and it is obviously caused by numerical effects. Thus, this
approach is not suitable for the implementation of a physical model for gas
entrainment.The SIMMER model assumes bubbly flow, where the gas void fraction is low and it assumes droplet flow,
and droplets in gas, where the gas void fraction is high. A variation of the
drag force—either by modification of the assumed particle
diameter for bubbles and drops or by using different correlations for the drag
coefficients of spherical particles—does not have a significant effect on the gas
entrainment in the simulations performed with the SIMMER model. The gas
entrainment is overestimated in all simulations, and there is no free parameter
inside the SIMMER model, which allows the modification of the amount of gas
entrainment in the simulation.The algebraic interfacial area density (AIAD) model was found to be a suitable
approach to adjust the entrainment. There are two free parameters inside the
AIAD model that have a strong influence on the suction of gas across the liquid
interface. These parameters are a drag coefficient for the free surface and the
shape of the blending function. The blending function is used to identify
regions of stratified flow (free surface flow) and regions of dispersed phase
flow (bubbly flow) in order to apply the appropriate drag model.Nevertheless, the gas entrainment calculated with the AIAD model is arbitrary, as the model
does not realistically reflect physics of the bubble entrainment. Further
investigations will be performed to improve the parameterization in terms of
the AIAD model.(i)
In the AIAD model, only the magnitude of the gas
void fraction is evaluated by the blending function to identify the location of
the free surface. In a more sophisticated approach, more criteria could be
evaluated such as the gradient of the gas void fraction.(ii)
It is not clear which drag coefficients for the free surface are appropriate. The
drag coefficient should also reflect the roughness of the jet surface, for
example.(iii)
Up to now, the blending function is meant to identify the location of the free
surface. With a more complex algorithm, it might be possible to identify the
region where the jet entrains gas. This would allow applying special closure
models (e.g., drag forces) to obtain a more controlled gas entrainment there.(iv)
The parameters for a realistic entrainment probably depend on the grid resolution.(v)
Up to now, the literature about gas entrainment near impinging jets is rather
fragmentary. More experimental data are necessary to adjust the CFD models and
obtain realistic entrainment in simulations. The Forschungszentrum
Dresden-Rossendorf (FZD) is planning to perform new experiments with impinging
jets. New sensors for multiphase flow measurements, which have been developed
at the FZD, will also be used.



---

*Source: 148436-2008-09-22.xml*