7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Simulations of the Bubbly TwoPhase Flow Around the Research Vessel Athena
Alejandro M. Castro and Pablo M. Carrica
IIHRHydroscience and Engineering, The University of Iowa, Iowa City, IA 52240
amcastro@engineering.uiowa.edu and pcarrica@engineering.uiowa.edu
Keywords: polydisperse bubbly flow, multigroup approach, ship hydrodynamics, breakup and coalescence, air entrainment
Abstract
Simulations of the fullscale unsteady twophase flow around the research vessel Athena are presented. A three
dimensional polydisperse model for bubbly twophase flows is solved. Comparisons are made with and without
breakup and coalescence in order to assess impact in the solution. The model is solved numerically using a multigroup
approach in which bubble sizes are discretized in groups of bubbles with a mean velocity and number density uniform
inside the group, resulting in four scalar equations for each size group. Using available experimental data good
estimations of entrained bubble size distribution and void fraction are imposed as one of the inputs to the model.
Predictions made by the model are then compared with experimental data and further conclusions are drawn.
Introduction
Twophase flows around ships have been studied for
years, mostly in relation with ship acoustic signatures:
ships have been tracked acoustically since before World
War II (Borowski et al. 2008). Bubbles are generated
at the bow and shoulder breaking waves, at the hull/free
surface contact line, the propeller and the highly turbu
lent stem flow. These bubbles are further transported
downstream by the flow forming a twophase mixture in
the wake that can be kilometers long. Recently, bubble
induced gasreduction has attracted increasing interest.
Externally injected bubbles have achieved drag reduc
tions ranging from I' to 2' in ship and flat plates
(Latorre et al. 2003; Takahashi et al. 2003).
One of the purposes of this work is to implement a
three dimensional model for the transport of bubbles
and apply it to solve the twophase flow around the US
Navy research vessel Athena II. The simulations pre
sented in this work are performed using a polydisperse
model (Carrica et al. 1999) that considers several bub
ble sizes and includes intergroup transfer mechanisms
such as breakup and coalescence. The simulations are
performed in full scale and this represents an important
contribution since in previous works the liquid phase
was computed in model scale, while the gas phase was
computed in full scale (Carrica et al. 1998, 1999; Mor
aga et al. 2008). Full scale simulations of bubbly flows
are important since the turbulent quantities need to be
computed in full scale in order to accurately predict
breakup and coalescence rates. Moreover, full scale
simulations provide of the correct turbulent mixing no
present in model scale simulations. With the recent ex
periments performed on the Athena II R/V by Johansen
et al. (2010), experimental data in full scale is available.
This experimental data is used to set the entrainment
model size distribution and strength at the bow. The pre
dicted results from the simulation at the stem are then
compared with the available experimental data at that lo
cation. Comparisons are made between simulations with
and without breakup/coalescence in order to assess their
importance.
TwoFluid Model
The modeling of twophase flow phenomena followed in
this work is based on a twofluid formulation developed
by Drew & Lahey (1979) and later applied to the pre
diction of twophase flows around ships (Carrica et al.
1998, 1999; Moraga et al. 2008).
This model describes one main phase which is
referred herein as the continuous phase and a diluted
phase which is referred herein as the dispersed phase.
For the present case, the continuous phase is water and
the dispersed phase constitutes the gas phase in the from
of bubbles.
Continuous Phase. The presence of one phase
or the other at a certain point in space and time is
described statistically. The equations of motion for
each phase are obtained by ensemble averaging the
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
NavierStokes equations conditioned to be in one phase
or the other.
The ensemble averaged momentum and mass conser
vation equations for the continuous phase can be written
as (Drew & Passman 1999; Carrica et al. 1999)
Du
Dt
1
VP+
V (2acve f VSu) + gad (1
_kP
0 + V (au) 0 (2)
where ac and ad are the volume fractions for the
continuous and dispersed phase respectively and which
are related by ad + ac 1. p is the piezometric pressure
related to the absolute pressure p by p = pcgz
and V"u is the symmetric part of the velocity gradient.
Vf f = v + t is the effective viscosity and vt is the
turbulent viscosity modeled with a suitable turbulence
model.
Polydisperse Model for Bubbly Flows. The sta
tistical description of bubbly twophase flows is
developed based on the Boltzmann theory of gases.
The fundamental quantity in this model is the bubble
mass distribution function f(m, r, t). This field gives
the mean number of bubbles per unit volume located at
position r at time t and which posses a mass within dm
of mass m. The use of the mass m as an independent
variable is advantageous for the simulation of the flow
around a ship since this quantity is preserved under
pressure changes. From the bubble mass distribution
function other important twophase parameters can be
obtained by integration:
Bubble number density
N(r,t)= dm f(m, r, t) (3)
Jo/
Gas volume fraction
a(r,t) dm f(m, r, t) (4)
Jo Pd
where Pd is the dispersed phase density. In the follow
ing, subscripts c denote the continuous phase and sub
scripts d denote the dispersed phase.
The conservation equation for the bubble mass dis
tribution, i.e. the Boltzmann transport equation, can be
written as (Lavalle et al. 1994)
8f (m, r, t)
+f( rt (ud(m, r, t)f(m, r, t))+
at
(m rf(m, r, t))
p(m, r, t) + x(m, r, t) + S(m, r, t) (5)
where 3 is the net source due to breakup, X is the net
source due to coalescence and S is a production term
that models the entrainment of bubbles due to wave
breaking and spilling. The third term on the left of Eq.
(5) represents the mass change of bubbles due to con
densation, evaporation or dissolution. Dissolution is an
important process in the prediction of the bubbly wake
behind a ship specially for small bubble sizes for which
the ratio of interfacial area to volume and the surface ten
sion induced pressure are high. However, this term is not
considered in this work where the main objective is to
implement and evaluate a fully twoway coupled poly
dispersed model for a case with complex wall bounded
geometry. This term must and will be considered in fu
ture work.
If only binary interactions are considered, the breakup
source can be written as
p(m,r,t) +(m,r,t) +P(m,r,t)
3+(m, r, t)
3 (m, r, t)
Sdm' h(m, m')b(m') f(m', r, t)
r(m)f(m, r, t) (6)
where /3 is the birth term due to the breakup of larger
bubbles than m into bubbles of size m and 3 is the
death term due to the breakup of bubbles of size m.
b(m) is the breakup rate of bubbles of size m and
h(m, m') is the daughter bubbles size distribution. Since
two bubbles are assumed to be formed per breakup, the
daughter bubbles size distribution is normalized to two
m
j dm h(m, m') 2
The coalescence sources can be written as
for births
x (m)
S dm' f (m m')f(m')Q(m 
2 0
m', m'X8)
and for deaths
/0
X (m) f(m) din' f(m') Q(m,m ') (9)
where X+ is the birth term due to the coalescence of
bubbles of mass m' and m m' and X is the death
term due to the coalescence of bubbles of size m with
bubbles of any other size. Q(m, m') is the coalescence
kernel.
Multigroup approach
Eq. (5) must be solved for every bubble size m. This
can be accomplished using a multigroup scheme where
the basic idea is that bubble sizes between ,* 1/2 and
S/2 can be represented by a single size ,. Then
the range of bubble masses is discretized into G groups
ranging from the minimum bubble size mi to the maxi
mum bubble size mG. The intermediate masses are de
fined as ,. t/2 ('' 1 + ..'.)/2 and the two ex
tremes are defined as m1/2 = 0 and m +1/ 2 = o.
The integration of Eq. (5) between 1/2 and
/2 gives
8N9
+ V. (,, +)) = 3, + x, + S (10)
where the ggroup bubble density number is defined as
N,(r,t) M 1/2
T9 g 1/2
dm f(m, r, t)
(11)
And the ggroup velocity is defined such that
S(r,t)N9(r,t), ,,
9 1g/2
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
procedure used for the continuous phase. Assuming that
for the dispersed phase inertia and shear stresses can be
neglected the ensemble averaged momentum equation
reduces to
aYVp+ acpcg (1
PC) k +
Where M. accounts for the interfacial momentum trans
fer between the continuous phase and group g. This can
be split into different contributions. The contributions
considered in this work are drag, virtual mass, lift and
turbulent dispersion. Additionally, a packing force that
accounts for the interaction between bubbles is consid
ered. Hence, the interfacial momentum transfer can be
written as
r/, = MD + MM + M + MTD + MP (19)
The drag force may be expressed as
dm ud(m, r, t)f(m, r, t)
39, x, and S, are the breakup, coalescence and so
terms in group g respectively.
In the multigroup approach the bubble distribi
function is assumed to have the form
G
f(m,r,t) = N,6(m
n=a
(12)
'urce
M D _a. 11
Mg C~PcCL 1, I
i "', 9
where u,,, is the relative velocity defined as u, ,
ution . u. The drag coefficient is that given by Tomiyama
(1998) with a void fraction correction based on the work
by Ishii & Zuber (1979) to include the interaction be
tween bubbles. For deformed bubbles in contaminated
(13) water the drag force may be expressed as (Moraga et al.
2008)
Hence, Eq. (3) and (4) reduce to
N(r,t)
ad(r,t)
G
S91 (rt)
9=1
9=1
Sp(r, t) N(t) (15)
1 24 8 Eo
CD  max (1 + 0.168Re075), 8 Eo
ac Rey 3 Eo + 4(
(21)
where the bubble Reynolds number Re, is defined with
the bubble relative velocity u,,, and the bubble diameter
D, as Re ...,,gDg/v6. The Eitvos number Eo is
defined as
S(Pc pd)D
where the group void fraction is defined as
a9(r, t)
d( N(r, t)
Pd7r t)
Using the assumption in Eq. (13) in Eq.
group velocity is
i.e. the dispersed velocity for the group size ..
Dispersed Phase Momentum Equation. In or
der to obtain the group g velocity (r, t) present in
Eq. (10) a momentum equation for every bubble group
size has to be solved. This momentum equation
is obtained by following the same ensemble averaging
where a is the surface tension. The correction for bub
(16) ble deformation in the terminal velocity is important for
large bubble sizes.
(12) the The equivalent radius R, is computed upon the as
sumption of spherical bubbles
3 ,,, 1/3
R Y (23)
R 47 pd(r, t) (23)
where the dispersed phase density is made to be a func
tion of space and time since compressibility effects due
to local pressure are considered by using an ideal gas
law of the form
Pd(r,t) = Pd,O
~ (rt) = Ud(". Tt)
where po is the atmospheric pressure and pd,o is the
dispersed phase density at this pressure. The pressure
in Eq. (24) includes the capillary pressure by using the
YoungLaplace equation.
The virtual mass force accounts for the effect of
accelerating the fluid surrounding the bubble and it can
be modeled as (Drew & Passman 1999)
MVM %=gapCvCVMX
( at v at+..v,)
(25)
where the virtual mass coefficient CVM is taken to be
0.5, the theoretical value for dilute potential flow of
spherical bubbles.
For the lift force, the result obtained from ensem
ble averaging the force over an sphere submerged in a
shear flow is used (Drew & Lahey 1987, 1990)
ML _ypcCLUr X (V X Uc) (26)
Turbulent dispersion is modeled in many works as an ex
tra forcing term in the gas momentum equation (Carrica
et al. 1999; Moraga et al. 2008; Bertodano et al. 2004).
In this work the simpler approach of a diffusion in the
number density equation is taken. Hence, Eq. (10) is
modified to be
at+v., )
V V ) + 3ig + g + S, (27)
where Scb is the Schmidt number here taken to be
Scb = 1.
A packing force is added in order to account for
the repulsive interaction between bubbles that occurs at
high void fractions. Following the developments made
for particle flows this force is modeled as (Gidaspow
1994; Apte et al. 2003)
Mf= aYVup(ad)
where up models a collision pressure between the bub
bles. As no further information is available here we pro
pose for the collision pressure an expression of the form
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
stable. The intention of the model provided by Eq.
(29) is to avoid excessive accumulation of bubbles,
specially near walls regions. According to Eq. (29) the
packing force only acts for void fractions larger than 0.3.
Breakup Modeling. Bubbles can break into smaller
bubbles by interacting with turbulent eddies. The model
developed by Luo & Svendsen (1996) is used
b(m) 
/ )1/3 1 (1 +I f )i r
0.29 2 d dS J Xc(X)
Rb 11/3
(30)
where E is the turbulent kinetic energy dissipation and
Xc (X, ) is the critical dimensionless energy for breakup,
computed as
1.89Z(X)u
xc(X, =)
PX ,2 ( /3R 11/3
Z(X) = 2/3 + (1 X)2/3 1 (31)
The lower limit in Eq. (30) is (,mi = Amin,/Db and
Ai, is taken to be the minimum size of the eddies in
the inertial subrange of isotropic turbulence, Ais. This is
related to the microscale by Ai,/Ad 11.4 31.4. The
value 11.4 is adopted in this work. The microscale can
be estimated by
S12. 3/4
Ad 12.56 3/
1/4
The implementation of a numerical integration for the
double integral in Eq. (30) is very intensive computa
tionally when applied to the solution of large tridimen
sional problems since it has to be computed for every
group at every grid point. The solution adopted to solve
this problem was to numerically compute the breakup in
Eq. (30) and to construct a lookup table in terms of two
dimensionless variables; one was chosen to be Q,, and
the other
1.89(
PE2/3Rb3
A uniform distribution is adopted for the daughter
bubbles size distribution, i.e.
2
h(m, m') =
ad < 0.3
0.3)2 ad > 0.3
where P is an adjustable parameter that we set to be
the largest possible value for which the equations are
Coalescence Modeling. In this work it is assumed
that the dominant mechanism for bubble collisions is
the random motion caused by turbulent fluctuations.
P (ad) P
F ^ J
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
According to the model by Prince & Blanch (1990) the
turbulent collision rate for bubbles of groups j and k is
Tk 1.41(Rj + R E)2 1/3 /3 + /3)1/2 (34)
After collision, the bubbles have to be in contact long
enough to allow the liquid film between them to drain
out. In order to compute the fraction of these colli
sions that result in an actual coalescence it is necessary
to asses the coalescence probability. Prince & Blanch
(1990) propose for the coalescence probability
Cjk = tik/Tjk (35)
where tjk is the time required for the coalescence to oc
cur and Tjk is the contact time between these bubbles.
The time required for coalescence is estimated as
Rk ( ') 1/2 ho
i 16T ) *khf
where ho is the initial thickness and hf is the thickness
for the film to break. These are estimated to be ho
10 4 m and hf 10 7 m. The equivalent radius Rjk
is computed as
11Y
R = 0.5( + ) (37)
\Rj Rk
The contact time depends on the bubble size and in
the turbulent intensity. This is estimated from
R2/3
T r (38)
Turbulence Model
Several turbulence models have been proposed for mod
eling twophase flows. Twoequation turbulence models
for twophase flows have been developed in the past as
an extension to the kE model. This approach is taken in
the works by Lahey et al. (1993), Lahey & Drew (2001)
and Troshko & Hassan (2001) where small variations ex
ist between them mainly in how the effect of the bubble
induced turbulence is taken into account.
In this work a twophase k c blended model is im
plemented. A k c counterpart to the k E model can
be found by following the procedure in Menter (1994)
where the equation for c is obtained by the formal sub
stitution E = 3*ck into the equation for E.
The resulting twophase k c blended model reads
Dk
QCt [,(T+ Vu
V [a(v + ckvt)Vk] ac, (39)
DcU
acD
Dt
Sa C VRe 
k
a3ccJ2 + V V [a.(v + vt)Vcc] +
2a,(1 FI) 2 Vc. Vk (40)
co
where TR, is the Reynolds stress tensor computed as
TRe = 2VtVS'. The turbulent viscosity is modeled
as the superposition of two effects. The first one is the
shear induced turbulence modeled by the k c turbulent
viscosity
k2
vsi = CA (41)
and the second effect is the bubble induced turbulence
here predicted using the multigroup version of the model
proposed by Sato et al. (1981)
G
VBI 1.2 R ., II.,, (42)
91
Hence, the turbulent viscosity is modeled by
Vt = VSI + VBI
Entrainment Model
One of the main issues when simulating the twophase
flow around a ship is the modeling of air entrainment.
Entrainment of air into the liquid phase is observed in
wave breaking/spilling and in the contact line between
the free surface and the hull of the ship. In the past, the
location and strength of the entrainment was set by the
user (Carrica et al. 1998). This required certain experi
ence from the modeler. More recently, some researchers
developed the so called subgrid scale models (Moraga
et al. 2008; Ma 2009). The model by Moraga et al.
(2008) was developed for stationary cases and unphys
ical entrainment occurs for instance in the presence of
heading waves. The main objective in this work is to val
idate the transport and intergroup transfer mechanisms
described herein and an Adhoc entrainment model is
proposed.
The entrainment model proposed here has the form
S(m, r, t) = SoD(m)E(r, t) (44)
where So is the entrainment strength, E(r, t) is a func
tion that varies between zero and one that detects the
entrainment locations and D(m) is the entrainment size
distribution.
The free surface slope is defined to be
0 cos 1(*. k) = cos 1(, ) (45)
an the distance to the free surface is denoted by p. The
idea is to activate the entrainment of bubbles when the
free surface slope is higher than a certain critical value
0,,t and the depth is smaller than a limit value e,,t.
Based on this, the following functional form for
E(r, t) is proposed when 0 > 0,nt and 4 < e,,t
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
of bubbles. In this work the fixed pivot technique
presented in Sanjeev & Ramkrishna (1996) is used.
Using this method the breakup terms appearing in Eq.
(10) are computed as
ad) a,, ad
(46)
E(r, t) = 0 in any other case. Here ac, is a critical
void fraction. The functional dependence in the void
fraction models the fact that there cannot be entrainment
if the whole phase is air. The dependence in the angle
smoothly increases the source strength from the criti
cal value Oent to the extreme case where the free sur
face is upside down. The dependence with depth aims
to smooth the transition between entrainment and no en
trainment.
The value of the entrainment parameters in this work
are a, = 0.7, Oent = 17 and Oent = 0.01.
Numerical Method
The code CFDShipIowa v4.0 (Carrica et al. 2007a,b)
was extended to solve the twophase model equations. It
uses multiblock structured bodyfitted grids with overset
gridding capability to model complex geometries and
perform local refinement where needed.
Number Density Transport. Numerically, the
convection term in Eq. (10) is discretized using a
second order TVDSuperbee scheme as in Ismail et al.
(2010). The use of a TVD scheme for the convection
term reduces the chance of having a nonphysical nega
tive value for N,. The TVD scheme does not guarantee
having a positive N, everywhere in the domain, but
it performs better than linear schemes with the same
order of accuracy. The TVDSuperbee limiter is chosen
since, as shown by Ismail et al. (2010), it outperforms
the usual linear schemes and other TVD schemes when
applied to general nonorthogonal curvilinear grids.
This was verified in bubble density transport tests.
Usual linear schemes present small spurious oscillations
that in many cases compare to the positive value of
the solution, resulting in negative values in some areas
of the domain, specially near entrainment sources and
the walls. This problem is practically removed by the
TVDSuperbee scheme.
Intergroup Transfer. The numerical treatment of
the intergroup transfer terms in Eq. (10) requires special
care. Naive discretization of these terms would result
in a non conservation of the total gas mass and number
G
39 = bN,
where b = b(, .) and h,,,, is computed as
h,,, fn9 1
M 9
dm m
dmn h(m, ,) +)
tl [  tl
121
dm h(mn, ')
The first and second integral terms reduce to zero for
g g' and g = 1 respectively.
For coalescence, the births are computed as
mg/,, ,,g
m= m
s',s"
9 T9O
mr m
(1 1 ,,,"~I/' ,Ng Ng,
2 ( 9 ,9 9
(50)
where = + ', Q9,9" = Q( ') and
{ mg <
r(r) 1 1< < (51)
Mge ge t i .
whereas the deaths are computed according to
G
Xg9 N9 1 Qg,gNg,
In order to obtain a stable time marching scheme, both
the breakup and coalescence terms are solved implicitly.
The death terms for breakup and coalescence are lumped
together to the diagonal term in the resulting discrete
system of equations making the scheme to be very stable
for the time steps used in ship applications.
Simulation of the Research Vessel Athena
The US Navy research vessel Athena II is a decom
missioned PG84 Ashevilleclass patrol gunboat trans
formed into highspeed research vessel in 1976. The wa
ter length of the Athena is L = 47 m and propelled by
a gas turbine it can reach maximum speeds of 18 m/s.
The Athena II R/V is fitted with a skeg, starboard and
port roll stabilizers and a compound masker system to
entrain bubbles and reduce the ship's radiated noise. The
masker is a ring fitted around the hull at approximately
x/L = 0.45. Fig. 1 is a detail of the computational
E(r,t)
(ent ) 18 0ent
Pent /180 ent
geometry used in the simulations. Shown in this figure
are the skeg, stabilizers and masker. The geometry also
includes the propeller shaft and the struts that hold it in
place. The computational domain used in the simula
tions only contains half of the ship and symmetry along
the y axis is assumed. The grid system is shown in Fig.
2 and consists of 21 base grids.
Strut
Shafts
Skeg
Stabilizers
Masker
Figure 1: Computational geometry used in the simula
tions of the Athena II R/V.
Refinement blocks are located at the bow and masker
in order resolve in fine detail the breaking waves that
occur there. There are three levels of refinement in the
bow with 6.3 millions of nodes. The masker has three
levels of refinement as well with a total of 2.9 millions
of nodes. A high level of refinement was used in these
zones since for future simulations the subgrid air en
trainment model from Ma (2009) and probably others
will be tested. Being this a subgrid scale model, a high
level of resolution is needed to capture plunging waves
as required by the model. Additional refinement blocks
are also located at the stem. The mesh system has a total
of 19.6 million grid points.
The conditions for the simulations are set to match
those of the experiments performed by Johansen et al.
(2010), i.e. the simulations are performed at full scale.
It is important to perform the simulations in full scale in
order to correctly predict the turbulent field since this
finally determines the breakup and coalescence rates.
Moreover, the correct prediction of the turbulent field
will reproduce the turbulent mixing more accurately.
From the conditions of the experiment the ship veloc
ity is Uo = 5.4 m/s. From Johansen et al. (2010) it is
known that the ship encountered small amplitude head
ing waves with a wavelength A = 35.2 m. The effect of
the incoming waves will not be considered in the present
simulation and the ship is fixed at even keel condition.
The Froude and Reynolds number for the above velocity
are Fr = 0.252 and Re = 2.53 x 108 respectively. From
now on, unless otherwise specified, all the variables are
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
nondimensionalized with the ship length L and velocity
Uo.
First, a single phase simulation was run for five ship
lengths in order to have a developed solution as a start
ing point for the twophase simulation. The time step
was set to 6t = 0.005 and 1000 time steps were run.
7 wall clock hours were needed to perform the com
putation on a Cray XT5 machine using 174 processors
i.e. 0.42 minutes per time step (The overset solver is not
needed for these computations since the grids are fixed).
The solution for this Froude number is characterized by
a wet transom and a small breaking wave is resolved at
the bow as it was expected, see Fig. 3.
The group sizes were chosen in order to span the
whole range of bubbles present in the experiment. The
entrainment size distribution D(R,) is set to be the nor
malized bubble number density function measured in the
bow by Johansen et al. (2010). This data was first fitted
with a power law of the form f(r) b ra. The exponent
was found to be a = 2.35. f(r) was then integrated
over every group size in order to get the contribution to
that group D,
D, N dr f (r) (53)
m g_ /2
This way
G
ED 1 (54)
9 1
The twophase simulation was set up with a total of 15
group sizes. These are summarized in Table 1 together
with the group distribution D,.
Table 1: Summary of the group
group distribution.
sizes and entrainment
Group Radius [pm] D,
1 20.00 6.04E01
2 29.67 2.39E01
3 44.01 9.46E02
4 65.29 3.75E02
5 96.86 1.48E02
6 143.69 5.87E03
7 213.17 2.32E03
8 316.23 9.19E04
9 469.12 3.64E04
10 695.93 1.44E04
11 1032.40 5.70E05
12 1531.50 2.26E05
13 2272.00 8.93E06
14 3370.50 3.54E06
15 5000.00 2.32E06
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Figure 2: Mesh system used in the computation of the Athena II R/V.
(a) Stern breaking wave and wake.
(b) Masker breaking wave.
(c) Bow breaking waves.
Figure 3: Detail of the breaking waves resolved with the computational mesh used in the simulations.
The strength of the entrainment model was set to
So 3 x 1017 in order to match the void fraction in
the bow found in the experiments. In the simulations the
maximum void fraction in the bow measurement point
(x, y) = (0.284, 0.0779) is 1%. The same time step
used in the single phase simulation was used for the two
phase simulation.
Two cases were run. One case with and another with
out intergroup transfer mechanisms. This allowed to
compare the bubble size distributions and from this com
parison the effect of bubble breakup and coalescence
was established. Five ship lengths were run needing a
total of 20 wall clock hours (3 minutes per time step) for
the full model and 14.2 hours for the same case without
intergroup transfer (2.13 minutes per time step).
Fig. 4a shows the normalized group void fraction ob
tained at the bow where the measurements were taken
when breakup and coalescence are active. The result
shown in this figure is no more than the entrainment dis
tribution since breakup and coalescence are negligible in
this region. In fact, Fig. 4a shows that the normalized
size distribution is independent of depth. Fig. 4b shows
the bubble mass distribution function. This is a func
tion of depth because the void fraction is. The bubble
size distribution at this location in the bow is taken as a
reference and other size distributions are compared with
it.
Fig. 5a shows the normalized group void fraction
computed in the stern at (x, y) (1.0106, 0.021) for
the case in which breakup and coalescence are not con
sidered. This corresponds to the location of the pole in
the experiment by Johansen et al. (2010). Near the free
surface this distribution behaves similarly to that shown
in Fig. 4. This indicates that near the free surface the
bubble size distribution is mainly controlled by the en
trainment of bubbles. As depth is increased the distribu
tion shifts to smaller bubbles. The source of this small
bubbles is analyzed below. Fig. 5b shows the bubble
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
102 103
1
106
S104
102
100
102
102 103
R[p m] R[g m]
(a) Normalized group void fraction. (b) Bubble distribution function
Figure 4: Computed size distribution at the bow. (x, y) = (0.284, 0.0779).
(a) Normalized group void fraction. Intergroup transfer not included. (b) Bubble distribution function. Intergroup transfer not included.
0.25 101
e0.1 m
00.4 m 1010 
0.2 0.518 m
90.6 m 10
0.15 0.7 m 6
10
00 10m
0.05 80.55 m
100 0.6 m
2 0.7 m
0 2 10.2 2
10' 10 10 102 10
R[g m] R[g m]
(c) Normalized group void fraction. Intergroup transfer included. (d) Bubble distribution function. Intergroup transfer included.
Figure 5: Computed size distribution at the stern. (x, y) = (1.0106, 0.021) at several depths.
e0.1 m
0.2 m
0.3 m
0.4 m
0.1
0.05
distribution function for this case.
Figs. 5c and 5d show the normalized group void frac
tion and bubble size distribution computed at the same
location in the stem when bubble breakup and coales
cence are considered. The same shift toward small bub
bles with depth that was observed in the case with no
intergroup transfer is present in this case as well. In
this case however, the effect of breakup near the free
surface is evident due to the drop in the distribution for
large bubbles. For intermediate depths, the distribution
presents a characteristic double dome shape with peaks
in bubble sizes of about 95 /m and 685 /m. Fig. 6
shows the measurements from Johansen et al. (2010) at
this same location. Even when the depths shown in Figs.
5c and 6 are not exactly the same, the qualitative agree
ment with the experiment is excellent. The experiment
shows a double dome shaped distribution as well, with
peaks in bubble sizes 80 pm and 800 im. It is of inter
R[p m]
Figure 6: Normalized group void fraction measured at
the stem. From Johansen et al. (2010)
est to discuss possible reasons for this agreement. The
bubbles shown in these figures could come whether from
the entrainment and then subsequent mixing at the stem
or from below the hull. In order to verify this, the bubble
size distribution is determined at a point well inside the
boundary layer and far away from the entrainment at the
same y location than in Fig. 5. The axial location is set
to be x = 0.978 i.e. at 1.5 m upstream the measurement
location right below the hull. Fig. 7 shows the bubble
size distributions found at this point as a function of the
distance to the hull wall for both cases, with and without
intergroup transfer. Both cases present a shift to small
bubbles with increasing distance to the wall/depth. The
reason for this shift is due to the difference in relative
velocity u,9, for the different bubble sizes. Small bub
bles have a very small relative velocity and therefore are
transported with a velocity that is very close to the liquid
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
phase velocity. Thus large bubbles tend to rise up to the
surface and leave the boundary layer faster than smaller
bubbles. The very interesting feature that is present in
Fig. 7c but is not in Fig. 7a is a second peak at bub
bles of size 685 pm. This peak can only be caused by
the coalescence of smaller bubbles since in Fig. 7a these
large bubbles are not present. The peak finally drops
for bubbles larger than 685 pm evidencing the effect of
breakup. Therefore, the bubble size distribution at this
point in the boundary layer is strongly determined by a
balance between breakup and coalescence. It is inter
esting to observe that the 685 pm peak forms close to
the wall due to the larger turbulent fluctuations in the
boundary layer that control the breakup and coalescence
rates. Moreover, bubbles tend to accumulate near the
wall giving rise to an increase in void fraction that also
increases the coalescence rates. Finally, the compari
son of Figs. 7c and 5c allows to conclude that the bub
bles with the double shaped dome in Fig. 5c are com
ing from below the hull and that this is reflected in the
measurements shown in Fig. 6. These results, however,
are inconclusive since no experimental data is available
inside the boundary layer. Shear induced breakup (Mik
sis 1981; Clift et al. 1978; Acrivos 1983), tipstreaming
(De Buijn 1993) and other effects not considered herein
could, and probably will break the bigger bubbles before
they can effectively accumulate near the wall. Another
interesting feature of the solution obtained is that bub
bles that are suck down below the hull from the bow and
contact line entrainment can be transported downstream
and eventually be trapped by the vortex shedding com
ing from the intersection of the shaft with the hull of the
ship. This can be seen in Fig. 8 where the isosurface
of a 8 x 10 5 is shown. These bubbles travel be
tween the struts, hit the rudder and finally leave a finger
shaped profile in the wake as shown in Fig. 10. It is im
portant to observe that large bubbles are more difficult
to suck down and even in that event, the strong breakup
along the boundary layer breaks them down into small
bubbles. This is the reason why the finger shaped pro
file cannot be observed for the large bubble sizes. A
depleted bubbly wake behind the stabilizers and struts
is also resolved and shown in Fig. 8. The void frac
tion profile all along the domain can be observed in Fig.
9. The several stages in the formation of the finger like
wake profile is clearly shown.
Conclusions
A three dimensional polydisperse twofluid model was
implemented in the code CFDShipIowa v4.0 (Carrica
et al. 2007a,b). The model is then applied to calculate
the twophase flow around the US Navy research vessel
Athena II. The novel features of this simulation is that
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
102 103 104
R[g m]
102 103
R[g m]
(a) Normalized group void fraction. Intergroup transfer not included. (b) Bubble distribution function. Intergroup transfer not included.
R[g m]
(c) Normalized group void fraction. Intergroup transfer included.
R[g m]
(d) Bubble distribution function. Intergroup transfer included.
Figure 7: Computed size distribution at the boundary layer at (x, y) = (0.978, 0.021) at several wall distances.
it includes the detailed complex geometry of the vessel,
it is performed in full scale and uses the latest breakup
and coalescence models available in the literature (Luo
& Svendsen 1996; Prince & Blanch 1990). New avail
able experimental data from Johansen et al. (2010) was
used in order to entrain bubbles with a good estimate for
the bubble size distribution and void fraction. Predic
tions made by the model were then compared with the
available experimental data at the stern of the ship. The
predictions made with the model show excellent quali
tative agreement with the measurements. The numerical
simulations enable analysis of the bubbly flow around
the ship at locations very difficult to reach experimen
tally. One of the points analyzed is the boundary layer
immediately before the transom of the ship. The analy
sis of the bubble size distribution at this point indicates
that a stream of bubbles comes from underneath the hull
of the ship and that these bubbles could be the ones de
tected in the experiments at the stern. The level of res
olution allows analysis of fine details in the transport
of the bubbles. It was found that bubbles are attracted
by vortex shedding generated at the intersection of the
shaft with the hull of the ship. A depleted wake be
hind the stabilizers and struts is also found. The results
of the simulations can be further improved by consider
ing breaking mechanisms that were not included in the
present model such that shear induced breakup (Miksis
1981; Clift et al. 1978; Acrivos 1983) and tipstreaming,
De Buijn (1993), and by using newer entrainment mod
els such as the one by Ma (2009).
Acknowledgements
This research was sponsored by the Office of Naval Re
search, grant N000140811084, under the administra
tion of Dr. Patrick Purtell.
e0.0220 m 
e0.0496 m
0.105 m
v0.270 m
0.435 m
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Figure 8: Isosurface of void fraction a = 8 x 10 5.
Figure 9: Slices at different x locations showing contours of void fraction. Coloring is in logarithmic scale.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
0.02
0.01
0
0.01
N 0.02
0.03
0.04
0.05
0.06
Y Y
(a) Group 5. (b) Group 8.
4.00E03
1.49E03
5.52E04
2.05E04
7.61E05
2.83E05
1.05E05
3.90E06
1.45E06
5.38E07
2.00E07
0.02
0.01
0
0.01
N 0.02
0.03
0.04
0.05
0.06
0
0.02
0.01
0
0.01
N 0.02
0.03
0.04
0.05
0.06
4.00E03
1.49E03
5.52E04
2.05E04
7.61E05
2.83E05
1.05E05
3.90E06
1.45E06
5.38E07
2.00E07
U.Ut U.UU U.UO
Y
(c) Group 10.
U U.U U.U.UU V.VO U.UO U.I IZ
Y
(d) Group 12.
Figure 10: Bubbly wake behind the Athena II R/V at x = 1.01.
0.02
0.01
0
0.01
N 0.02
0.03
0.04
0.05
0.06
References
Acrivos A., The breakup of small drops and bubbles in
shear flows, Ann. New York Acad. Sci., Vol. 404, pp.
111, 1983
Apte S.V., Mahesh K. and Lundgren T., A eulerian
lagrangian model to simulate twophase/particulate
flows, Center for Turbulent Research. Annual Research
Briefs, 2003
Bertodano M.L, Moraga F.J., Drew D.A. and Lahey Jr.
R.T., The modeling of lift and dispersion forces in two
fluid model simulations of a bubbly jet, J. Fluids Eng.,
Vol. 126, pp. 573577, 2004
Borowski B., Sutin A., Roh H.S. and Bunin B., Passive
accoustic threat detection in estuarine environments,
Proc. SPIE 6945, 694513, 2008
Carrica P.M., Bonetto F.J., Drew D.A. and Lahey Jr.
R.T., The interaction of background ocean air bubbles
with a surface ship, Int. J. Numer. Meth. Fluids, Vol. 28,
pp. 571600, 1998
Carrica P.M., Drew D.A., Bonetto F.J. and Lahey Jr.
R.T., A polydisperse model for bubbly twophase flow
around a surface ship, Int. J. Multiphase Flow, Vol. 25,
pp. 257305, 1999
Carrica P.M., Wilson R.V. and Stem F., An unsteady
singlephase level set method for viscous free surface
flows, Int. J. Numer. Meth. Fluids, Vol. 53, pp. 229256,
2007
Carrica P.M., Wilson R.V., Noack R.W. and Stem F.,
Ship motions using singlephase level set with dynamic
overset grids, Comput. Fluids, Vol. 36, pp.14151433,
2007
Clift R., Grace J.R. and Weber M.E., Bubbles, drops and
particles, Academic Press, New York, 1978
De Buijn R.A., Tipstreaming of drops in simple shear
flows, Chem. Eng. Sci., Vol. 48, pp. 277284, 1993
Drew D.A. and Lahey Jr. R.T., Application of general
constitutive principles to the derivation of multidimen
sional twophase flow equations, Int. J. Multiphase Flow,
Vol. 5, pp. 243, 1979
Drew D.A. and Lahey Jr. R.T., The virtual mass and lift
force on a sphere in rotating and straining inviscid flow,
Int. J. Multiphase Flow, Vol. 13, pp. 113121, 1987
Drew D.A. and Lahey Jr. R.T., Some supplement analy
sis on the virtual mass and lift force on a sphere in rotat
ing and straining inviscid flow, Int. J. Multiphase Flow,
Vol. 16, pp. 11271130, 1990
Drew D.A. and Passman S.L., Theroy of Multicompo
nent Fluids, Applied Mathematical Sciences, Vol. 135,
1999
Gidaspow D., Multiphase flow and fluidization: contin
uum and kinetic theory descriptions, Academic Press,
Boston, 1994
Ishii M. and Zuber N., Drag coefficient and relative ve
locity in bubbly, droplet or particulate flows, AIChE J.,
Vol. 25, pp. 843855
Ismail F., Carrica P.M., Xing T. and Stem F., Evalua
tion of linear and nonlinear convection schemes on mul
tidimensional nonorthogonal grids with applications to
KVLCC2 tanker, Int. J. Numer. Meth. Fluids, In press,
DOI:10.1002/fld.2318, 2009
Johansen J.P., Castro A.M. and Carrica P.M., Fullscale
twophase flow measurements on Athena research ves
sel, Submitted to Int. J. Multiphase Flow, 2010
Lahey Jr. R.T., Lopez de Bertodano M., and Jones Jr.
O.C., Phase distribution in complex geometry conduits,
Nucl. Engng. Des., Vol. 141, pp. 177201, 1993
Lahey Jr. R.T. and Drew D.A., The analysis of twophase
flow and heat transfer using a multidimensional, four
field, twofluid model, Nucl. Engng. Des., Vol. 204, pp.
2944, 2001
Latorre R., Miller A. and Philips R., Microbubble resis
tance reduction on model SES catamaran, Ocean Eng.,
Vol. 30, pp. 22972309, 2003
Lavalle G.G., Carrica P.M., Clausse A. and Qasi M.K., A
bubble number density constitutive equation, Nucl. En
gng. Des., Vol. 152, pp. 213224, 1994
Luo H. and Svendsen H.F., Theoretical model for drop
and bubble breakup in turbulent dispersions, AIChE J.,
Vol. 42, pp. 12251233, 1996
Ma J., Oberai A.A., Drew D.A., Lahey Jr. R.T., and Mor
aga F.J., A quantitative subgrid air entrainment model
for bubbly flows plunging jets, Comput. Fluids, In
press, DOI:10.1016/j.compfluid.2009.07.004, 2009
Menter F.R., Twoequation eddyviscosity turbulence
models for engineering applications, AIAA Journal, Vol.
32, pp. 15981605, 1994
Miksis M.J., A bubble in an axially symmetric shear
flow, Phys. Fluids, Vol. 24, pp. 12291231, 1981
Moraga F.J., Carrica P.M., Drew D.A. and Lahey Jr.
R.T., A subgrid air entrainment model for breaking bow
waves and naval surface ships, Computers & Fluids, Vol.
37, pp. 281298, 2008
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Prince M.J. and Blanch H.W., Bubble coalescence and
breakup in airsparged bubble columns, AIChE J., Vol.
36, pp. 14851499, 1990
Sanjeev K. and Ramkrishna D., On the solution of pop
ulation balance equations by discretizationI. A fixed
pivot technique, Chem. Eng. Sci., Vol. 51, pp. 1311
1332, 1996
Sato Y., Sadatomi M. and Sekoguchi K., Momentum and
heat transfer in twophase bubbly flowsI, Int. J. Multi
phase Flow, Vol. 7, pp. 167177, 1981
Takahashi T., Kakugawa A., Makino M. and Kodama Y.,
Experimental study on scale effect of drag reduction by
microbubbles using very large flat plate ships, J. Kanai
Soc. N.A., Vol. 239, pp. 1120, 2003
Tomiyama A., Struggles with computational bubble dy
namics, In: Proc. third int. conf. on multiphase flows,
ICMF98, Lyon, France, 1998
Troshko A.A. and Hassan Y.A.,A two equation turbu
lence model of turbulent bubbly flows, Int. J. Multiphase
Flow, Vol. 27, pp. 19652000, 2001
