Group Title: 7th International Conference on Multiphase Flow - ICMF 2010 Proceedings
Title: 11.1.1 - Simulations of the Bubbly Two-Phase Flow Around the Research Vessel Athena
ALL VOLUMES CITATION THUMBNAILS PAGE IMAGE ZOOMABLE
Full Citation
STANDARD VIEW MARC VIEW
Permanent Link: http://ufdc.ufl.edu/UF00102023/00273
 Material Information
Title: 11.1.1 - Simulations of the Bubbly Two-Phase Flow Around the Research Vessel Athena Bubbly Flows
Series Title: 7th International Conference on Multiphase Flow - ICMF 2010 Proceedings
Physical Description: Conference Papers
Creator: Castro, A.M.
Carrica, P.M.
Publisher: International Conference on Multiphase Flow (ICMF)
Publication Date: June 4, 2010
 Subjects
Subject: polydisperse bubbly flow
multigroup approach
ship hydrodynamics
breakup and coalescence
air entrainment
 Notes
Abstract: Simulations of the full-scale unsteady two-phase flow around the research vessel Athena are presented. A three dimensional polydisperse model for bubbly two-phase 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.
General Note: The International Conference on Multiphase Flow (ICMF) first was held in Tsukuba, Japan in 1991 and the second ICMF took place in Kyoto, Japan in 1995. During this conference, it was decided to establish an International Governing Board which oversees the major aspects of the conference and makes decisions about future conference locations. Due to the great importance of the field, it was furthermore decided to hold the conference every three years successively in Asia including Australia, Europe including Africa, Russia and the Near East and America. Hence, ICMF 1998 was held in Lyon, France, ICMF 2001 in New Orleans, USA, ICMF 2004 in Yokohama, Japan, and ICMF 2007 in Leipzig, Germany. ICMF-2010 is devoted to all aspects of Multiphase Flow. Researchers from all over the world gathered in order to introduce their recent advances in the field and thereby promote the exchange of new ideas, results and techniques. The conference is a key event in Multiphase Flow and supports the advancement of science in this very important field. The major research topics relevant for the conference are as follows: Bio-Fluid Dynamics; Boiling; Bubbly Flows; Cavitation; Colloidal and Suspension Dynamics; Collision, Agglomeration and Breakup; Computational Techniques for Multiphase Flows; Droplet Flows; Environmental and Geophysical Flows; Experimental Methods for Multiphase Flows; Fluidized and Circulating Fluidized Beds; Fluid Structure Interactions; Granular Media; Industrial Applications; Instabilities; Interfacial Flows; Micro and Nano-Scale Multiphase Flows; Microgravity in Two-Phase Flow; Multiphase Flows with Heat and Mass Transfer; Non-Newtonian Multiphase Flows; Particle-Laden Flows; Particle, Bubble and Drop Dynamics; Reactive Multiphase Flows
 Record Information
Bibliographic ID: UF00102023
Volume ID: VID00273
Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
Resource Identifier: 1111-Castro-ICMF2010.pdf

Full Text



7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 -June 4, 2010


Simulations of the Bubbly Two-Phase Flow Around the Research Vessel Athena


Alejandro M. Castro and Pablo M. Carrica

IIHR-Hydroscience 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 full-scale unsteady two-phase flow around the research vessel Athena are presented. A three
dimensional polydisperse model for bubbly two-phase 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

Two-phase 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 two-phase mixture in
the wake that can be kilometers long. Recently, bubble
induced gas-reduction 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 two-phase 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.


Two-Fluid Model

The modeling of two-phase flow phenomena followed in
this work is based on a two-fluid formulation developed
by Drew & Lahey (1979) and later applied to the pre-
diction of two-phase 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


Navier-Stokes 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 two-phase 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 two-phase 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 two-way 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 g-group bubble density number is defined as


N,(r,t) M 1/2
T9 g -1/2


dm f(m, r, t)


(11)


And the g-group 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
Young-Laplace 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 look-up 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
1-1Y
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 two-phase flows. Two-equation turbulence models
for two-phase flows have been developed in the past as
an extension to the k-E 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 two-phase 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 two-phase 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)
9-1

Hence, the turbulent viscosity is modeled by


Vt = VSI + VBI


Entrainment Model


One of the main issues when simulating the two-phase
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 sub-grid 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 Ad-hoc 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 CFDShip-Iowa v4.0 (Carrica et al. 2007a,b)
was extended to solve the two-phase model equations. It
uses multiblock structured body-fitted 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 TVD-Superbee 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 TVD-Superbee 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 non-orthogonal 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
TVD-Superbee 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)
Mg-e g-e 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 PG-84 Asheville-class patrol gunboat trans-
formed into high-speed 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 sub-grid air en-
trainment model from Ma (2009) and probably others
will be tested. Being this a sub-grid 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


non-dimensionalized 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 two-phase 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 two-phase 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.04E-01
2 29.67 2.39E-01
3 44.01 9.46E-02
4 65.29 3.75E-02
5 96.86 1.48E-02
6 143.69 5.87E-03
7 213.17 2.32E-03
8 316.23 9.19E-04
9 469.12 3.64E-04
10 695.93 1.44E-04
11 1032.40 5.70E-05
12 1531.50 2.26E-05
13 2272.00 8.93E-06
14 3370.50 3.54E-06
15 5000.00 2.32E-06







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

10-2


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
-e--0.1 m
0--0.4 m 1010 -
0.2- -0.518 m
-9--0.6 m 10
0.15 --0.7 m 6
10


00 10m

0.05 -8--0.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.


-e--0.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 two-fluid model was
implemented in the code CFDShip-Iowa v4.0 (Carrica
et al. 2007a,b). The model is then applied to calculate
the two-phase 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 N00014-08-1-1084, under the administra-
tion of Dr. Patrick Purtell.


-e-0.0220 m -
-e-0.0496 m
-0.105 m
-v-0.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.00E-03
1.49E-03
5.52E-04
2.05E-04
7.61E-05
2.83E-05
1.05E-05
3.90E-06
1.45E-06
5.38E-07
2.00E-07


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.00E-03
1.49E-03
5.52E-04
2.05E-04
7.61E-05
2.83E-05
1.05E-05
3.90E-06
1.45E-06
5.38E-07
2.00E-07


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.
1-11, 1983

Apte S.V., Mahesh K. and Lundgren T., A eulerian-
lagrangian model to simulate two-phase/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. 573-577, 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. 571-600, 1998

Carrica P.M., Drew D.A., Bonetto F.J. and Lahey Jr.
R.T., A polydisperse model for bubbly two-phase flow
around a surface ship, Int. J. Multiphase Flow, Vol. 25,
pp. 257-305, 1999

Carrica P.M., Wilson R.V. and Stem F., An unsteady
single-phase level set method for viscous free surface
flows, Int. J. Numer. Meth. Fluids, Vol. 53, pp. 229-256,
2007

Carrica P.M., Wilson R.V., Noack R.W. and Stem F.,
Ship motions using single-phase level set with dynamic
overset grids, Comput. Fluids, Vol. 36, pp.1415-1433,
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. 277-284, 1993

Drew D.A. and Lahey Jr. R.T., Application of general
constitutive principles to the derivation of multidimen-
sional two-phase 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. 113-121, 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. 1127-1130, 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. 843-855

Ismail F., Carrica P.M., Xing T. and Stem F., Evalua-
tion of linear and nonlinear convection schemes on mul-
tidimensional non-orthogonal 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., Full-scale
two-phase 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. 177-201, 1993

Lahey Jr. R.T. and Drew D.A., The analysis of two-phase
flow and heat transfer using a multidimensional, four
field, two-fluid model, Nucl. Engng. Des., Vol. 204, pp.
29-44, 2001

Latorre R., Miller A. and Philips R., Micro-bubble resis-
tance reduction on model SES catamaran, Ocean Eng.,
Vol. 30, pp. 2297-2309, 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. 213-224, 1994

Luo H. and Svendsen H.F., Theoretical model for drop
and bubble breakup in turbulent dispersions, AIChE J.,
Vol. 42, pp. 1225-1233, 1996

Ma J., Oberai A.A., Drew D.A., Lahey Jr. R.T., and Mor-
aga F.J., A quantitative sub-grid air entrainment model
for bubbly flows plunging jets, Comput. Fluids, In
press, DOI:10.1016/j.compfluid.2009.07.004, 2009

Menter F.R., Two-equation eddy-viscosity turbulence
models for engineering applications, AIAA Journal, Vol.
32, pp. 1598-1605, 1994

Miksis M.J., A bubble in an axially symmetric shear
flow, Phys. Fluids, Vol. 24, pp. 1229-1231, 1981

Moraga F.J., Carrica P.M., Drew D.A. and Lahey Jr.
R.T., A sub-grid air entrainment model for breaking bow
waves and naval surface ships, Computers & Fluids, Vol.
37, pp. 281-298, 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
break-up in air-sparged bubble columns, AIChE J., Vol.
36, pp. 1485-1499, 1990

Sanjeev K. and Ramkrishna D., On the solution of pop-
ulation balance equations by discretization-I. 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 two-phase bubbly flows-I, Int. J. Multi-
phase Flow, Vol. 7, pp. 167-177, 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. 11-20, 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. 1965-2000, 2001




University of Florida Home Page
© 2004 - 2010 University of Florida George A. Smathers Libraries.
All rights reserved.

Acceptable Use, Copyright, and Disclaimer Statement
Last updated October 10, 2010 - Version 2.9.7 - mvs