7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
A generalized subgrid air entrainment model for RaNS modeling of bubbly flows
around ship hulls
Jingsen Ma* Assad A. Oberai* Mark C. Hyman! Richard T. Lahey, Jr*
and Donald A. Drew*
Center for Multiphase Research. Rensselaer Polytechnic Institute. 110 8th St. Troy, NY 12180. USA.
t Naval Surface Warfare CenterPanama City, Panama City, Florida. 32407 USA
maj6@rpi.edu, oberaa@rpi.edu, mark.c.hyman@navy.mil, laheyr@rpi.edu and drewd@rpi.edu
Keywords: subgrid air entrainment, bubbly flows, multiphase simulation, ship hulls
Abstract
The quantitative prediction of bubbly flow around a maneuvering Naval surface ship is critical in determining its
hydrodynamic performance and its acoustic and optical signatures. It is a challenging multiscale problem that relies
heavily on subgrid models of turbulence and air entrainment. In this manuscript we analyze this problem using a
novel air entrainment model that predicts the location and rate of air entrainment around a surface ship. This model
was coupled with a twofluid Reynolds averaged Navier Stokes (RaNS) model and used to evaluate the flow field
around a Naval surface ship in straight ahead and turning motions. For straight ahead motion the predicted void
fraction distributions aft of the stern were compared with experiments at three different speeds and good agreement
was found. The qualitative differences in the location of the air entrainment and the resulting bubbly flow between
straight ahead and steady turning motions are also discussed and compared with experimental observations. To our
knowledge this study contains the first quantitative comparisons between predicted and experimental void fraction
distributions for a fullscale surface ship.
Introduction
A maneuvering surface ship at sea typically entrains
large quantities of air due to breaking bow waves, un
steady transom flows and wave/structure interactions
along the hull. This is evidenced by the "white water"
that extends many ship lengths behind the stern wake.
Bubbly flow around a Naval surface ship is important
because it determines its acoustic and optical signatures
and may alter its drag characteristics. Thus the simula
tion and prediction of the air entrainment and the result
ing bubbly flow around a surface ship is of interest in
Naval and marine hydrodynamics.
Over the last few decades igmiik.illm progress (see for
example Carrica et al. 1999; Hyman 2006; Moraga et al.
2008) has been made in predicting bubbly flows around
surface ships by utilizing multifluid models described in
Drew and Passman (1999) and Lahey (2009). However,
to date quantitative comparisons with experimental data
have not been feasible partly due to the lack of a ro
bust and accurate subgrid air entrainment model. Clearly
the location and rate of air entrainment play a critical
role in determining the accuracy of a multiphase simu
lation. Moraga et al. (2008) have previously proposed
a subgrid air entrainment model and applied it to simu
lating the flow around a surface ship with some success.
This model assumed that bubbles are present everywhere
close to the free surface and they are entrained once the
downward liquid velocity exceeds the bubble rise veloc
ity. However, this model does not account for the fact
that if the free surface moves downward with the same
velocity as the liquid, no bubbles are entrained. Further,
it only provides an expression for the location of air en
trainment and not its rate, which must be chosen empir
ically. Several models for determining the air entrain
ment rate have also been proposed for canonical flows
such as plunging liquid jets (see for example Bin 1993;
( li.iii n 1996; Ma et al. 2010b) and hydraulic jumps
(see for example ( II.iMi, in 1996) and they can be used in
select regions around a surface ship, but none that can
capture all the sources simultaneously.
Recently, we have developed a subgrid air entrain
ment model (Ma et al. 2009, 2010a) that assumes that
the turbulence in the liquid produces a rough air/liquid
interface with cavities of a certain size. These cavities
are entrained into the liquid once the downward liquid
velocity exceeds the downward velocity of the interface.
This leads to a simple expression for the location and
rate of air entrainment that is proportional to the product
liquid's turbulent kinetic energy and the gradient of the
liquid velocity in the downward (pointed away from the
interface) direction. This model has been implemented
in the framework for a twofluid model (Drew and Pass
man 1999; Lahey 2009) and applied to canonical prob
lems such as a plunging liquid jet and a hydraulic jump.
In both cases the same model yielded accurate results
(Ma et al. 2009, 2010a). In this manuscript we apply it
to model the twophase flow around the research vessel
(RV) Athena. We consider both straight ahead motion
and a steady turn. For the former we compare our pre
dictions with atsea measurements of void fraction. For
the latter no such data is available but we are able to
make qualitative comparisons between the two motions
and aerial photographs from a sea test. We note that the
air entrainment model contains one parameter that de
pends on material properties of the liquid and the gas
components of the bubbly flow. Since in all the simu
lations in this manuscript the liquid is sea water and the
gas is air, we have frozen this parameter.
The format of the remainder of this manuscript is as
follows. In Section 2 we describe the RaNS twofluid
model and the subgrid air entrainment model that are
used in this study. In Section 3, we present simulations
of bubbly flows around the RV Athena, in both straight
ahead and steady turn motions. We discuss the results
and compare the predicted void fraction profiles with the
available experimental data. We end with Conclusions in
Section 4.
Multiphase Modeling Framework
In this section we summarize the derivation of the sub
grid air entrainment model and describe its implementa
tion within a RaNStype, twofluid, computational mul
tiphase fluid dynamics (CMFD) code. For simplicity,
we restrict our simulations to monodispersed bubbles
and to oneway coupling. We model the fluid as com
posed of a continuous liquid phase and a dispersed gas
phase comprised of bubbles of a given diameter. This
diameter was selected to be the typical bubble diameter
measured in the experiments. Readers interested in the
effects of twoway coupling and polydisperse modeling
are referred to the work of Moraga et al. (2008).
Air entrainment model The subgrid air entrain
ment model we used may be derived as follows. Con
sider a gas/liquid interface, Fi, with roughness and air
cavities of size z a (see Fig. 1). In addition F is an
imaginary surface below the interface at a distance a/2.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
We require that F moves with the same downward veloc
ity as the interface. Thus any air that passes through F is
entrained. The average volume of air passing through F
per unit area and time is then given by:
V (o)
V,(a/2)
'7/,:
Figure 1: Schematic of a rough air/water interface.
q C= C(u,(a/2)
~
SC(u,(a/2)
C2Ou" (\))
w ^wnO~a
Un(0))
Un(0))
where Un is the downward liquid velocity at the inter
face, an overbar indicates an ensemble average and
Sf f >0
/(f) 1 f < 0
A simple approximation for the roughness (Hunt 1984)
is a u'2/g = k/g, where k is the local turbulent ki
netic energy and g is the acceleration due to gravity. So,
C 9un, k
2 On g
Based on this, the equivalent volume sources per unit
volume of water in a liquid layer of thickness Pent is
given by:
q C 0O, k
Q = = (4)
Pent 2 dn ,entg
Thus the source term for the number density of bubbles
per unit volume of water is:
Q
CentK( an~ k
On . 
where vg is the volume of an average bubble, Cent = 2
is the model coefficient that depends on the physical
properties of the liquid and the gas in the multiphase
flow.
Ja
Mass conservation of the dispersed phase
Due to the assumption of momodispersed bubbly flow,
bubble coalescence and air dissolution play no role. In
this case, the conservation equation for the bubble num
ber density, N.", for bubbles of a characteristic diam
eter, D,, moving with a velocity of ug, often referred
to as the population balance equation (see for example
MartinezBazan et al. 1999), is given by:
a N1'
at+ V. ( N[") = (6)
a t
where , is rate at which bubble density increases due
to air ingestion and is given by eq. (5). Note that other
sources of bubbles on the right hand side of eq. (6), such
as those due to bubble breakup and coalescence, can also
be included if appropriate.
Momentum balance of the dispersed phase
The ensembleaveraged balance of momentum equation
for the dispersed phase is given by Moraga et al. (2006),
o(vN "p.ug)
dU) + V. (v9N .'pd, u) (7)
at
Y N"Vc + 7N"'p g + M
where vg is the volume of a bubble with diameter Dg,
Pd is the air density, and pc is the pressure of the liquid.
Note that Mg, the fluctuating interfacial force density,
must be constituted. This term primarily determines mo
mentum transfer between the continuous phase and the
dispersed phase. It may be partitioned to account for
different types of interactions,
MT MD MTD MV ML + M (8)
where the right hand side contains contributions due to
drag (D), turbulent dispersion (TD), virtual mass (VM),
lift (L) and wallinduced forces (W). For the specific ex
pressions for each of theses terms the reader is referred
to Moraga et al. (2008).
Conservation laws for the continuous phase It
is assumed that the continuous liquid phase is incom
pressible and is not affected by the dispersed phase (one
way coupling). As a result the continuity equation for
the liquid phase simplifies to,
V uc = 0. (9)
The ensembleaveraged statement of the balance of lin
ear momentum for the continuous phase is (Moraga et al.
2008),
p + u (10)
at+ V. ue uc = V (T + TR") +pg (10)
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Where T = pl + 2pcD is the Cauchy stress for a
Newtonian fluid, and TR" is the Reynolds stress tensor
modeled as:
In the expressions above I is the identity tensor, /c is
the liquid's viscosity, Dc is the rate of strain of the liq
uid phase, and k and pt are the turbulent kinetic energy
and the turbulent viscosity, respectively. In this work,
a blended (k c)/(k c) turbulence model developed
by Menter (1994) was used to determine k and pt, and
hence to construct TZ". All model coefficients are given
in Menter (1994). For further details on the implementa
tion of this model, the reader is referred to Moraga et al.
(2008).
Modeling the free surface The free surface of the
air/water mixture is represented using a singlephase
level set function p described in Sussman et al. (1994),
which represents the signed distance from the free sur
face, with the level set = 0 representing the free sur
face. Its evolution is governed by:
 + u, VO = 0
at
For more details on the application of the singlephase
level set method to bubbly flow simulations the reader
is referred to the work of Carrica et al. (2006); Moraga
et al. (2008).
Simulation of bubbly flow around the Naval RV
Athena
Using the methodology described in the previous section
we have performed simulations of bubbly flows around
the hull of the research vessel Athena which is operated
by the Carderrock Division of the US Navy Naval War
fare Center. Athena results have been previously used by
Hyman (2006) and Moraga et al. (2008) as a test case for
verifying CMFD models. We first simulate the straight
ahead motion of Athena at three different ship velocities
and compare the void fraction profiles aft of the transom
with the experimental measurements reported by Terrill
(2006). Then we simulate the flow around Athena dur
ing a steady turn. Some interesting differences between
the straight ahead and steady turning motions are ob
served and discussed.
Straight ahead ship motion We model onehalf of
the ship by assuming symmetry about the hull and ignor
ing any asymmetrical instabilities in the flow. The grid
consisted of 3.7 x 106 cells with one and two levels of
Te _(2pk)I +2/tDc,
C
nested Chimera overset grids that were used to resolve
the hull appendages (i.e., the skeg and rudder). This grid
emphasizes a heavily refined overset block in the near
wake region but not in the bow region since our focus
was on the transom flows for which experimental mea
surements have been reported (Terrill 2006). The reader
is referred to Hyman (2006) and Moraga et al. (2008) for
further details about the grid and boundary conditions
employed in our numerical simulations.
Figure 2: Predicted rate of air entrainment close to the
free surface for RV Athena at 9 knots: (a) Overall view
in which the entrainment at the transom is clearly seen;
(b) Magnified view of the masker; (c) Magnified view of
the hullairwater contact line.
We have simulated the straight ahead motion of
Athena at 6, 9 and 10.5 knots, which corresponds to
Reynolds numbers of 1.46 x 108, 2.18 x 108 and 2.55 x
108, respectively, and Froude numbers of 0.14, 0.21 and
0.25, respectively, based on the ship speed and ship
length, L = 47.2m. The air entrainment model was
switched on once a statistical stationary state for the
liquid flow was obtained. A uniform bubble size of
0.08mm was assumed. For the air entrainment model,
the entrainment depth was set to Pent = 0.002L and
Cent = 11.8 was used. Using a time step of 0.005 in
nondimensionalized units, which translates to a dimen
sional time interval of about 0.06s, it took about 300
steps for bubbles to traverse the computational domain.
We started to collect data after this interval and moni
tored the evolution of timeaveraged values for the gas
velocities and void fractions until they attained a statis
tically stationary state.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Fig. 2 is a snapshot of the source term E, close to the
free surface. We observe strong entrainment just aft of
the transom stern and also see some scattered sources in
the nearwake region. In addition to these regions, air
entrainment was also observed along the hullairwater
contact line and at the masker. All these regions of en
trainment are consistent with visual observations at sea
(see Terrill 2006). The only exception is the bow region,
where we did not predict the expected entrainment be
cause we were unable to resolve the breaking bow wave
due to insufficient grid resolution (also see Moraga et al.
2008; Hyman 2006).
Figure 3: (a) Timeaveraged void fraction aft of the
transom stem along the ship's center plane, at a ship
speed of 9 knots; (b) Enlarged view of the region close
to the transom, that is the boxed part of (a).
The timeaveraged plots of void fraction at the cen
ter plane along the ship length are shown in Fig. 3.
The zoomedin view near the transom stern indicates the
presence of a large vortical structure that is responsible
for much of the air entrainment and is reminiscent of the
air entrainment mechanisms in a hydraulic jump (see Ma
et al. 2010a, 2009).
We have computed the timeaveraged void fraction
data along the center plane of the ship, 2m aft of the
transom. We compare our simulations with at sea mea
surements for RV Athena (Terrill 2006). This compar
ison was performed at three different ships speeds and
the results are presented in Fig. 4. The predicted results
are in good agreement with measured data. It should
be noted that we achieve this agreement using the same
value of the air entrainment parameter Cent at all speeds.
The main deviation occurs well below the interphase, at
a depth greater than im, where our simulation appears to
underpredict the void fraction likely because we did not
do a polydispersed prediction and thus did not predict
the smaller bubbles. Anyway the measured void frac
tion in this region is very low (around 10 4), and thus
it is prone to experimental error. Furthermore, some of
these bubbles have been generated due to cavitation at
the propeller, which was not modeled in our simulation.
,,
I~
'i ~~~~
'1.
Model
Exper.
103 102
mean void fraction
Model
 Exper.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
In any event our simulation appears to accurately predict
the void fraction distribution in the transom stem region.
Steady turn We have also computed the flow around
Athena undergoing a steady starboard turn with a turn
ing radius of Rt 4L. For this simulation the flow was
asymmetric and a grid around the entire ship was gen
erated having approximately 6.5 million points in 100
computational blocks. The ship speed during the turn
was U 10.5 knots, resulting in a Reynolds number
of 2.55 x 108 and a Froude number of 0.25. The val
ues used for the bubble diameter and the air entrain
ment coefficient were the same as for straight ahead ship
motion. The background grid extended sufficiently far
away from the ship for the assumption of undisturbed
waters in the farfield boundaries to hold. As a result
the boundary conditions in the far field were the same as
for straight ahead motion. The simulation was initiated
with the ship and the background grid advancing in the
prescribed turn and used a time step of At 0.01L/U.
This condition was continued for 100 time steps in order
to dissipate the startup transient. At time t L/U, the
bubble entrainment model was switched on and it con
tinued on for another 100 time steps.
10' 102
mean void fraction
103 102
mean void fraction
Figure 4: Predicted and measured void fraction distri
butions 2m aft of the transom stern along the ship's cen
ter plane, as a function of depth for straight ahead ship
speeds of of 6 (top), 9 (center) and 10.5 knots (bottom).
Figure 5: Predicted rate of air entrainment close to the
free surface for RV Athena in a steady turn: (a) Overall
view in which nonsymmetric air entrainment at the tran
som is clearly seen; (b) Magnified view of the masker;
(c) Magnified view of the hullairwater contact line.
Fig. 5 shows a snapshot of the predicted bubble
source strength close to the free surface around RV
Model
Experiment
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
straight ahead ship motion. From Fig. 7 we note that in
the turning ship case there is some entrainment ahead of
the transom, toward the left side of the picture, that is
missing from the simulation. We believe that this is due
to the interaction of a wavy sea with the contact line and
is not captured because we are simulating steady seas.
Figure 6: Void fraction close to the free surface for
Athena in straight ahead motion (right) and during a
steady turn (left).
Athena during a steady turn. We observe a wave orig
inating on the starboard (inside turn) corer of the tran
som. By comparing Figs. 2 and 5 we conclude that
in both the straight ahead and steady turn cases air en
trainment is seen at the contact line along the hull, at
the masker and particularly at the transom. This en
trainment leads to a region of "white water" aft of the
transom. Some ,ignilik.iIii differences are observed be
tween the bubbly wake of the straight ahead and steady
turn cases and they are highlighted in Figs. 6 and 7.
First, as expected, the bubbly wake for the steady turn
case is curved toward the path traced by the ship. Sec
ond, the void fraction distribution in the near wake is
markedly different. For the straight ahead case we ob
serve two symmetric streaks of high void fraction re
gion (void fraction over 0.2) perpendicular to the tran
som, originating at the the two ends of the transom. This
was also observed in the sea tests and is believed to be
caused by the diverging waves due to the transom's ge
ometry (see Fig. 7). In the turning case these streaks are
also observed, however they are closer together and they
last longer, implying higher rates of air entrainment. In
addition a third streak was observed at the inside cor
ner of the turn though it's shorter and weaker than the
former two (as seen in the left plot of Fig. 6). We be
lieve that this extra entrainment is caused by the comer
wave induced by the turning motion. Once again this
pattern was also observed in the sea tests (see Fig. 7).
In general it appears that the void fraction distribution
for the steady turn case is more heterogeneous than for
Figure 7: Top: predictions of void fraction close to
the free surface for Athena in straight ahead (right) and
steady turn (left) motion. Bottom: Corresponding sea
test pictures of the bubbly wake taken along a similar
orientation for same ship motions.
Conlusions
Three dimensional simulations of the bubbly flow
around the Naval research vessel, Athena, was per
formed using a RaNSbased twofluid modeling ap
proach and a novel subgrid air entrainment model. Con
sistent with experimental observations, the model pre
dicted entrainment at the masker, at the contact line, and
particularly at the transom of the ship. For straight ahead
motion the timeaveraged values of the predicted void
fraction profiles near the transom were compared with
experimental data at three different speeds and good
agreement was found. A single simulation was also
performed for a steady turn of the ship and qualitative
comparisons were made between the predicted straight
ahead and steady turn cases and the RV Athena visual
ization data taken at sea. It was found that the ship ma
neuvers in a turn leads to greater rates of air entrainment
and to a more heterogeneous bubbly wake.
Acknowledgments
This work was sponsored by the Office of Naval Re
search (ONR), grant N000140310826, under the ad
ministration of Dr. Patrick Purtell, and was supported
in part by a grant of computer time from the DOD High
Performance Computing Modernization Program at the
Maui High Performance Computing Center (MHPCC),
US Army Engineering and Research Development Cen
ter (ERDC) and Arctic Region Supercomputing Center
(ARSC). We owe thanks to Dr. Francisco Moraga at
the General Electric Corporate Research Center for very
helpful discussions about CMFD models and technical
details about ship modeling. We also give thanks to Dr.
Eric Terrill at Scripps Institution of Oceanography for
providing the experimental aboard RV Athena void frac
tion data.
References
A.K. Bin. Gas entrainment by plunging liquid
jets. Chemical Engineering Science, 48(21):35853630,
1993.
PM Carrica, D. Drew, F Bonetto, and RT Lahey Jr. A
polydisperse model for bubbly twophase flow around a
surface ship. International journal of multiphase flow,
25(2):257305, 1999.
PM Carrica, RV Wilson, and F. Stem. An unsteady
singlephase level set method for viscous free surface
flows. International Journal for Numerical Methods in
Fluids, 53(2):229256, 2006.
H. ( li.iii,'n Air Bubble Entrainment in FreeSurface
Turbulent Shear Flows. Academic Press, 1996.
D.A. Drew and S.L. Passman. Theory of multicompo
nentfluids. Springer Verlag, 1999.
J.C.R. Hunt. Turbulence structure and turbulent diffu
sion near gasliquid interfaces. In W. Brutsaert and G.H.
Jirka, editors, Gas Transfer at Water Surfaces, pages 67
 82. Reidel, Dordrecht, The Netherlands, 1984.
M. Hyman. Unsteady bubbly flow modeling near ship
hulls. In ONR Ship WaveBreaking and BubbleWake
Review, March 2006. Organized by Program Manager
Patrick Purtell.
R. T. Lahey, Jr. On the computation of multiphase flow.
J. Nuclear T,. ih. 1.... ., 167:117, 2009.
J. Ma, A.A. Oberai, D.A. Drew, R.T. Lahey Jr, and
M.C. Hyman. A comprehensive subgrid air entrain
ment model for Reynoldsaveraged simulations of free
surface bubbly flows. In APS DFD Meeting Abstracts,
2009.
J. Ma, A.A. Oberai, D.A. Drew, R.T. Lahey Jr, and
M.C. Hyman. A comprehensive subgrid air entrainment
model for rans modeling of bubbly flows near the free
surface, to be submitted, 2010a.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
J. Ma, AA Oberai, DA Drew, RT Lahey Jr, and FJ Mor
aga. A quantitative subgrid air entrainment model for
bubbly flowsplunging jets. Computers & Fluids, 39(1):
7786, 2010b.
C. MartinezBazan, JL Montanes, and JC Lasheras. On
the breakup of an air bubble injected into a fully devel
oped turbulent flow. Part 1. Breakup frequency. Journal
of Fluid Mechanics, 401:157182, 1999.
F.R. Menter. Twoequation eddyviscosity turbulence
models for engineering applications. AIAA journal, 32
(8):15981605, 1994.
FJ Moraga, AE Larreteguy, DA Drew, RT Lahey, et al. A
centeraveraged twofluid model for wallbounded bub
bly flows. Computers & Fluids, 35(4):429461, 2006.
FJ Moraga, PM Carrica, DA Drew, and RT Lahey Jr. A
subgrid air entrainment model for breaking bow waves
and naval surface ships. Computers & Fluids, 37(3):
281298,2008.
M. Sussman, P. Smereka, and S. Osher. A level set ap
proach for computing solutions to incompressible two
phase flow. Journal of computational Physics, 114(1):
146159, 1994.
E. J. Terrill. Measurements of the air entrainment around
the full scale surface ships. In ONR Ship WaveBreaking
and BubbleWake Review, March 2006.
