7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
A polydisperse twofluid model for bubble plume under breaking waves
Gangfeng Ma*, Fengyan Shi*and James T. Kirby*
Center for Applied Coastal Research, University of Delaware, Newark, DE 19716, USA
gma@udel.edu, fyshi@udel.edu and kirby@udel.edu
Keywords: polydisperse twofluid model, bubble plume, breaking waves
Abstract
Wave breaking in the surf zone entrains large volumes of air bubbles into the water column, forming a complete
twophase bubbly flow field. Numerical study of this bubbly flow is largely restricted by the lack of robust and
comprehensive bubble entrainment models. In this paper, we propose a new model that connects bubble entrainment
with turbulent dissipation rate at the airwater interface. The initial bubble size distribution follows the bubble
size spectrum observed in laboratory experiments. The locations where bubbles are entrained are bounded by a
threshold of turbulence dissipation rate. The model has two free parameters: bubble entrainment coefficient cb and
critical turbulence dissipation rate c. The bubble entrainment model as well as a polydisperse twofluid model are
incorporated into a 3D VOF code Truchas. The model is employed to study bubbly flow under a laboratory surfzone
breaking wave.
Introduction
Wave breaking in the surf zone generates intense tur
bulence and coherent structures eddy within the un
derlying flow field. As the wave breaks, a rather
twodimensional flow structure rapidly becomes three
dimensional, evolving into obliquely descending eddies
(Nadaoka et al., 1989) and downbursts of turbulent fluid
(Ting, 2008). These processes can entrain large vol
umes of air bubbles into the water column, enhancing
wave energy dissipation and airwater mass transfer. De
pending on their concentrations and size distribution,
the entrained bubbles can significantly change the op
tical properties of water (Terrill et al., 2001), introduc
ing potentially large errors in the opticallybased mea
surements. Bubbles can also affect the dynamics of the
flow field. Large void fractions near the surface create
a stable stratification, which can hinder vertical mixing
(Sullivan and McWilliams, 2010).
Early investigations on bubble entrainment and evolu
tion under surfzone breaking waves are mostly through
laboratory measurements. Deane and Stokes (1999,
2002) have conducted photographic studies on air en
trainment mechanism and bubble size distribution un
der laboratory plunging breaking waves. They revealed
that the bubble creation is driven by two largescale pro
cesses: the jet/waveface interaction and the collapsing
cavity. The first process is primarily responsible for the
formation of small bubbles with radius less than Hinze
scale (w 1mm), while the latter is mainly responsible
for the generation of bubbles larger than Hinze scale.
The bubble size spectrum of their measurements satis
fies a 3/2 power law for small bubbles and a 10/3 power
law for large bubbles.
Compared with laboratory experiments, numerical
studies of wave breaking induced twophase bubbly flow
field are still rare. The main reason is perhaps due to the
lack of robust and comprehensive bubble entrainment
models. Carrica et al. (1999) developed a polydisperse
twofluid model to study the bubbly flow field around
a surface ship, but didn't take into account bubble en
trainment processes. The bubbles were introduced into
the computation through measured data in plunging jet
experiments. Moraga et al. (2008) proposed a subgrid
air entrainment model for breaking bow waves. In their
model, bubble entrainment was modeled through a vol
ume source term in the bubble number density equation.
The locations where bubbles are entrained is determined
by the mean downward liquid velocity, which should be
greater than 0.22 m/s. The bubble source intensity is
specified to obtain a good comparison with measured
data. Their model has no criterion to specify the bub
ble source intensity which has spatial and temporal vari
ations. Additionally, the approach to determine bubble
entrainment locations is questionable, considering that
bubbles can also be entrained in regions where liquid
velocity is not downward, for example, bubble entrain
ment in a hydraulic jump. Shi et al. (2008) presented a
polydisperse twofluid bubbly flow model based on mix
ture theory. They formulated the air entrainment by con
necting it with turbulence production at the airwater in
terface. Simulation results showed that the model can
successfully capture the evolution pattern of void frac
tion. Subsequent work (Shi et al., 2010) revealed that
the model is very sensitive to simulation setup. They ar
gued that it is necessary to develop a more theoretically
justifiable air entrainment formulation.
The main objective of this paper is to develop a
physicallybased wave breaking induced bubbly flow
model. In our model, bubble entrainment under break
ing waves is correlated with turbulence intensity at the
free surface. It is assumed that the total energy required
for bubble formation is linearly proportional to the tur
bulent dissipation rate c. The entrained bubbles strictly
follow the bubble size spectrum as observed by Deane
and Stokes(2002). The locations where bubbles are en
trained are determined by turbulence dissipation rate e
which should be greater than a critical value c. The
model results in two free parameters: air entrainment
coefficient Cb and critical turbulence dissipation rate c%.
These two parameters have to be calibrated in the simu
lation. The developed air entrainment model as well as
a polydisperse twofluid model (Carrica et al., 1999) are
incorporated into a VOF code TRUCHAS and applied to
study the bubbly flow under a laboratory spilling break
ing wave.
Twofluid model
To simulate polydisperse twofluid flow, the dispersed
bubbles are separated into NG classes or groups. Each
class has a characteristic bubble diameter dbi, i
1, 2, .. NG, and a corresponding volume fraction of
a,9i. By definition, the volume fractions of all of the
phases must sum to one:
NG
a0+ Z g,= 1 (1)
al 1
where al is the volume fraction of liquid phase.
The polydisperse bubbly flow model in the current pa
per is based on the analysis of Carrica et al. (1999) who
neglected the inertia and shear stress tensors for the gas
phase due to the relatively small gas volume and den
sity. Following Moraga et al. (2008), we neglect bub
ble coalescence and gas dissolution. Thus the governing
equations become
O(pi) + V (tlplu) 0 (2)
8(alptlul)
( + V (alpluulu)
t
lVp + alpig(
(3)
+V [.. .,,(Vui + VTu,)] + Mgt
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
ON *
t + V (u9,N,Ni) =B,, + S,i (4)
agiVp + aip A + Mlg,i = 0 (5)
where pi is liquid density, ul is liquid velocity, p is pres
sure which is identical in phases, p,,i is the bubble den
sity of groupi, g is gravity, Peff,i is the effective vis
cosity of liquid phase, N,,i is bubble number density
of groupi, u,,i is bubble velocity, B,i is groupi bubble
source due to air entrainment, S,,i is the intergroup mass
transfer which only accounts for bubble breakup. The
bubble breakup model proposed by MartinezBazan et
al. (1999a,b) is used in the present paper. Mg1 and Ml,,i
are the momentum transfer between phases, which sat
isfy the following relationship
NG
MgY + Mli
i 1
Bubble source due to air entrainment
Deane and Stokes (2002) divided the lifetime of wave
generated bubbles into two phases; the acoustic phase,
where bubbles are entrained and fragmented inside the
breaking wave crest, and the quiescent phase, where
bubbles evolve under the influence of turbulent diffu
sion, advection, buoyant degassing and dissolution. The
acoustic phase is short lived and the time scale of bub
ble fragmentation is on the order of milliseconds. Given
these features, direct simulations of the acoustic phase
require high temporal and spatial resolutions in order to
capture the details of the air entrainment process, mak
ing their applications on the surfzonescale domain in
feasible. A practical way to introduce bubbles into the
computation is to prescribe air bubbles in a twofluid
model using a bubble entrainment formulation (Shi et
al., 2010). The model fed with the initially entrained
bubbles basically simulates bubbly flows in the quies
cent phase, thus requires much less spatial and temporal
resolutions.
As mentioned by Moraga et al. (2008), there are two
options to model air entrainment, namely a boundary
condition at the interface or a volumetric source in a
region close to the interface. The first option is prob
lematic because we normally do not resolve small scales
necessary for predicting bubble entrainment. Follow
ing Baldy (1993), we denote EA as the characteristic
turbulent energy available for the formation of bubbles
under breaking waves, Eb(a) as the energy required to
entrain a single bubble with a radius of a, and r as
the Kolmogorov dissipation length scale. In breaking
events, we will only consider a bubble range such that
EA > Eb(a) and a > r. In this range, the statisti
cal state of the smallscale disturbances is selfsimilar.
Therefore, bubble entrainment related to smallscale dis
turbances is also selfsimilar and independent of large
scale conditions (Baldy, 1993). Thus, it is justified to as
sume that bubble creation under breaking waves is solely
determined by e.
Based on the above analysis, we can develop the bub
ble entrainment model quantitatively. The energy re
quired to entrain a single bubble with a radius of a is as
sociated with surface tension, which is given by (Buck
ingham, 1997)
Eb(a) 47ra2c (7)
Where a is the surface tension. If we assume that the en
ergy required for bubble creation per second is linearly
proportional to c, we can get
Eb(a)B(a) CbPlC
where B(a) is the rate of bubble creation per cubic me
ter, cb is air entrainment coefficient which has to be cal
ibrated in the model. Thus, the bubble creation rate can
be evaluated as
B(a) = b (rT l2 (9)
47 pi
Equation (9) is quite similar to the model developed
by Baldy (1993) through dimensional consideration, but
it's only developed for the entrainment of a single size
group bubbles. In reality, bubbles under breaking waves
experience a large span of size distribution from mi
crometers to centermeters. In the following analysis,
the bubble size follows the laboratory measurements of
Deane and Stokes (2002). The minimum bubble radius
is taken as 0.1 mm, and the maximum is 10 mm. Since
our studies will be focused on laboratoryscale break
ing waves, the bubble size change due to the variation
of pressure and temperature is negligible. Therefore, we
separate bubbles into NG classes (or groups) with con
stant bubble radius. The characteristic bubble radius of
each class is al, a2, . aN. The width of each class
is ai+l/2 ai1/2, where i 1, 2, .. NG, ai+l/2 =
(ai + ai+l)/2. Thus we have a11/2 0.1mm and
aNG+/ 2 10mm. The bubble size spectrum is given
by Deane and Stokes (2002), who revealed that bubbles
larger and smaller than Hinze scale ( 1mm) respec
tively vary as 10/3 and 3/2 power law at the acoustic
phase.
S(a) oc a 10/3 if a > 1.0mm (10)
S(a) oc a3/2 if a < 1.Omm (11)
where S(a) is the bubble size spectrum.
Then, the bubble entrainment rate per cubic meter for
groupi can be written as
B(ai) = 5oS(ai)Aai (12)
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
where So is a coefficient, B(ai) is the entrainment rate
for groupi bubbles. We assume that the total energy
required for bubble entrainment is linearly proportional
to turbulent dissipation rate, then
NG
4) =
Plugging equation (12) into (13) gives
=Cb E (14)
47 pi G a S(ai)Aai
The polydisperse bubble entrainment model is then
given by
B(ai) = ) 1 (a)A(15)
4w P1 1 aS(ai)Aai
The air entrainment coefficient is expected to be less
than unity, because only a portion of turbulence energy is
used to entrain bubbles. In the numerical simulation, the
air entrainment coefficient has to be calibrated with mea
sured data. To complete the formulation, we still need to
describe how to select the grid points where bubbles are
entrained. It is straightforward to set a critical dissipa
tion rate c, that no bubbles will be entrained if e < e, at
the free surface. The criterion of choosing ec is to make
sure bubbles are only entrained after wave breaks. Nor
mally c, is taken as the turbulence dissipation rate of the
free surface cells at breaking point.
Momentum transfer
The momentum transfer between phases includes virtual
mass, lift force and drag force which is given as
Mlg,i = m + 'N + l" (16)
The virtual mass force which accounts for the acceler
ation of the liquid in the wake of the bubbles is given
by
Sl" agpiCvM( Dul Du ) (17)
Dt Dt
where CvM is the virtual mass coefficient with a con
stant value of 0.5. The D/Dt operators denote the sub
stantial derivatives.
Bubble rotation with finite relative velocity, or fluid
velocity gradients (shear motion), will induce a trans
verse component in the hydrodynamic force, which is
known as the lift force. The effect of lift force is mod
elled by
AlT ,iPlCL(Ul U9,i) x (V x U1)
where CL is the lift force coefficient, which is set to 0.5.
The drag force is originally due to the resistance ex
perinced by bubbles moving in the liquid. The momen
tum transfer by drag force is written as the following
form (Clift et al, 1978)
3 CD
I ag,iPl 8 R (u u,i) u1 ui (19)
8 R,i
where CD is drag coefficient depending on the flow
regime and liquid properties. For rigid spheres the drag
coefficient is usually approximated by the standard drag
curve (Clift et al, 1978)
24
CD g4(1 + 0.15Re0687) (20)
where Re is bubble Reynolds number
where Re,,i is bubble Reynolds number
Re9,i
aipi Iul u9,,i b
Turbulence model
We use a nonlinear k c model which is modified for
twophase bubbly flow (Troshko and Hassan, 2001) to
calculate turbulent eddy viscosity. The conservation
equations of turbulent kinetic energy k and turbulence
dissipation rate c are formulated as
(apik) + V (alpulk) V(aT Vk)
at a( (22)
+ a (G Pie) + Sbk
(alplc) + V (alplue) V (a Ve)
Ot (a T,
62
C,2Pi ) + Sbe
k
where the standard constants for k c model are ak =
1.0, ar, 1.3, Ci 1.44, C2 = 1.92. The term G is
the production of turbulent kinetic energy and described
by G = T : Vul. Ti is shear stress of liquid phase which
is calculated from the nonlinear Reynolds stress (Lin and
Liu, 1998).
The last two terms Sbk and Sb, are bubble induced
turbulence production. Troshko and Hassan (2001) have
proposed formulations for the bubbly flow in the vertical
duct
NG 3(D
NSbk d ,iPi Ur 3 (24)
i 1
Sb, Ce3bSbk (25)
where u, is relative velocity. C,3 is a new constant
which is found to be 0.45 (Troshko and Hassan, 2001).
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
wb is the characteristic frequency of bubbleinduced tur
bulence destruction originally given by Lopez de Berto
dano (1994)
S2Cvmdb, (26
wb =(3 (26)
where C,, and Cd are virtual mass coefficient and drag
coefficient.
Free surface tracking
In this study, a singlephase VOF method was used
to track the free surface. Unlike the twophase VOF
method, a singlephase VOF method does not solve
NavierStokes equations in the air region. This treat
ment has been found to significantly increase the stabil
ity of simulating highdensity ratio flows. It was proved
that singlephase VOF method is sufficiently accurate in
predicting complex flows such as wave breaking in surf
zone (Wu, 2004). To facilitate the numerical implemen
tation, we define the bubble velocities and void fraction
in the air region as zero. At the surface cells, the air and
dispersed bubbles could coexist. In reality, these cells
involve intensive interactions between air and dispersed
bubbles. For example, air will breakup to form dispersed
bubbles, and bubbles will reversely coalescence to be
come air. These interactions are neglected in the present
model. This ensures that the only bubbles introduced
into the flow are through air entrainment model. No flux
boundary conditions at free surface are set to make sure
bubbles can leave the computational domain freely.
Bubbly flow under breaking waves
Model setup
The polydisperse twofluid model is employed to study
the bubbly flow under a laboratory surfzone breaking
wave. Our attention will be focused on the void fraction
distribution and bubble evolution after wave breaking.
Laboratory measurements by Cox and Shin (2003) are
selected to test the model's capability. The experiment
was conducted in a 36 m long by 0.95 m wide by 1.5
m high glasswalled flume. A beach with constant slope
of 1:35 was installed with the toe 10 m from the wave
maker and intersecting the still water line at x=27.85 m.
The flume was filled with tap water to a depth of h=0.51
m. The selected test case is characterized by a spilling
breaker with incident wave height of 0.11 m and wave
period 2.0 s.
To reduce computational effort, a 2D simulation has
been conducted. The computational domain is taken
as 21 m long by 0.7 m high, with the beach toe 1.0
m from the left boundary (see figure (1)). The mea
sured breaking point is located at Xb = 13.07m.
+ a, (Cl G
k
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
0 1
0 2
Figure 1: Computational domain and measurement lo
cations
Three measurement sections are respectively located at
x=13.81m,13.94m,14.17m. Because the velocity mea
surement is only in streamwise direction, the turbulent
kinetic energy is not evaluated. Therefore, our compar
isons with the measured data will be focused on free sur
face elevation, streamwise velocity and void fraction.
The whole computational domain is discretized by a
uniform grid with Ax 0.025m, Ay 0.01m. The
time step is automatically adjusted during the computa
tion to satisfy stability constraints. Both mean veloci
ties and free surface displacement are specified on the
left boundary. On the top, pressure is set to be zero.
Adjacent to a solid wall, the lawofthewall boundary
conditions for k and c are applied (Wu, 2004). Bub
bles are divided into NG 20 groups with a logarith
mic distribution of bubble sizes (Figure (2)). Moraga
et al. (2008) pointed out that logarithmic distribution is
preferable over an equally spaced distribution as it helps
ensure that the ratio Adbi/dbi is small even for small dbi,
where Adbi is the width of the bin centered at dbi. They
also found no significant differences in the simulation
results for NG 15, 30 or 60, but NG < 15 should be
avoided.
There are two free parameters in the bubble entrain
ment model that have to be determined during the simu
lation: the critical turbulent dissipation rate c, and bub
ble entrainment coefficient Cb. The critical turbulent dis
sipation rate determines when and where the bubbles are
entrained. Its value is equal to the turbulent dissipation
rate on the free surface cell at breaking point. In the cur
rent case, we take c = 0.01m2/s3. As long as the bub
bles are entrained into the water column, the void frac
tion is mostly determined by the bubble entrainment co
efficient Cb. Therefore, the calculated void fraction dis
tribution is not very sensitive to ec. We found no big dif
ference in the numerical results if c, 0.005m2/s3 for
the current case. The bubble entrainment coefficient cb
has to be calibrated in the simulation. We take cb = 0.15
in the present simulation.
Figure 2: Initial bubble size distribution
Results
In this section, comparisons between experimental data
and numerical results are presented for free surface ele
vation, streamwise velocity and void fraction. In the ex
periment, these mean quantities were obtained by phase
averaging after the waves reached quasisteady state. In
our computation, in order to reduce computational ef
fort, we did not simulate the wave propagation for the
entire experimental period which is nearly 100 waves.
Instead, the entire computational period is 30 sec. The
numerical results show that the computed waves in the
surf zone are nearly in steady state after 30 seconds of
simulation.
Figure (3) shows the comparison of simulated and
measured wave height distribution along the beach. The
wave height is estimated from the computed free sur
face elevation between 26.0s and 28.0s. The sampling
frequency is 20Hz. As we can see, the simulated wave
height agrees reasonably well with measurement. The
wave height before wave breaking is a little overesti
mated. The model also predicts the wave breaking a
little earlier. These discrepancies could be partially be
cause the traditional k c turbulence closure model can
not accurately predict the initiation of turbulence in a
rapidly distorted shear flow such as breaking waves (Lin
and Liu, 1998). Fortunately, the distance between the
simulated and measured breaking point is short, which is
extremely important for the bubble simulation because
the location ofthe breaking point determines the initia
tion of bubble entrainment.
Cox and Shin (2003) showed that the temporal varia
tion of void fraction above/below the SWL normalized
by wave period and average void fraction appears to be
selfsimilar and can be modeled simply by linear growth
and exponential decay. These variations are also ob
served in our numerical results which are dispicted in
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
20
15 
F 10
> 5
a
16 18 20 22 24 26 28
Figure 3: Comparison of simulated and measured wave
height distribution
4
3
am 2
3
2a
0 0.02
0 0.02
0.04 0.06 0.08 0.1
t/T
Figure 4: Temporal variation of simulated void frac
tions above (upper panel) and below (lower panel) the
still water level
figure (4). We can clearly see that the exponential de
cay of simulated void fractions which are normalized
by wave period and average void fractions is well simu
lated. This result proves that the twofluid model is able
to successfully simulate bubble transport phenomenon.
The linear growth of the void fraction is captured as
well, but the growth rate is underpredicted comparing to
Cox and Shin's estimate. The reason is because the cur
rent bubble entrainment model is largely simplified by
assuming air entrainment is linearly related to the turbu
lence dissipation rate. In addition, the interactions be
tween large air cavities and dispersed bubbles are not
considered in the current model. But we should point
out that this level of accuracy is a considerable improve
ment over previous works in which the bubble source is
manually specified (Carrica et al., 1999; Moraga et al.,
2008).
As the bubble entrainment in our model is directly re
lated to turbulent dissipation rate, it is instructive to take
a look at turbulent dissipation rate and void fraction dis
0 0'
Figure 5: Temporal variation of simulated turbulent dis
sipation rate with time interval 0.2T
tribution simultaneously. These snapshots are shown in
figure (5) and (6), where the turbulent dissipation rate is
normalized by g g(h + H), g is gravity acceleration, h
is still deep water depth, and H is incident wave height.
The plotted void fractions are bounded by 0. :'. follow
ing Lamarre and Melville (1991, 1994) who used 0. '.
as a threshold to evaluate various moments of the void
fraction field. In these two figures, xb is the measured
breaking point. We can see that the evolution patterns of
turbulent dissipation rate and void fraction are similar.
When wave starts to break, the turbulent dissipation rate
will be greater than the critical turbulent dissipation rate,
thus bubbles start to be entrained into the water column.
As the breaking bore moves forward, the turbulent dis
sipation rate and void fraction all increase rapidly and
attain their maximum value shortly after wave break
ing. The peak void fraction appears around 0.5m on
shore from the breaking point. This result is consistent
with the findings of Mori et al. (2007), who argued that
the peak of void fraction happens at 0.1 to 0.2 wave
length onshore from the breaking point. After attaining
their maximum, both the turbulence dissipation rate and
void fraction start to decay. The high turbulence region
is persistently located at the breaking wave crest with a
moderate time variation in the value of dissipation rate,
10 (ttb)/T=0 1
0 0 078
0002
10
20
10 (ttb)/T=03
0'0 0547
0
10 0002
20
10, (ttb)/r=0 5
0 050 0408
10 0x002
20
10 (ttbYT= 07
0 0i02
10 000
20
10 (ttb)/T=0 9 00193
0002
10
20
0 05 1 15 2 25 3
x Xb (m)
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
while the void fraction varies by an order of magitude in
a wave period. The higher void fractions are restricted
in the roller region, but some bubbles are spread down
stream and form a long tail of the bubble cloud under the
wave surface.
Conclusions
We have developed a polydisperse twofluid model for
simulating bubbly flow under surfzone breaking waves.
An air entrainment model was proposed to account for
the bubble generation after wave breaking. This model
connects bubble entrainment to the turbulent dissipation
rate at the airwater interface. The model was utilized
to study bubbly flow under a laboratory surfzone break
ing wave. Numerical results showed that the model can
reproduce the experiment measurements for a spilling
breaking wave. The predicted temporal variation of void
fraction can be modeled with linear growth and expo
nential decay, which is consistent with the laboratory
measurements. The evolution of simulated void fraction
follows that of turbulent dissipation rate. The peak void
fraction appears around 0.5m onshore from the break
ing point. The high void fraction which is mostly con
tributed by large bubbles is generally located in the roller
region.
In the current model, the interactions between air cav
ity and dispersed bubbles have not been considered be
cause they are negligible in a spilling breaking wave.
In a plunging breaking wave, air cavities can be formed
during the splashup cycle. The air cavities can breakup
to form dispersed bubbles under the action of turbulence,
dispersed bubbles can reversely coalescence to become
air. In order to better simulate bubble entrainment pro
cesses in the surf zone, our future work will be focused
on incorporating these interactions into our model.
Acknowledgments
This study was supported by Office of Naval Research,
Coastal Geoscience Program under grant N00014101
0088. Sincere thanks also go to Dr. Cox and Dr. Shin
for making their experimental data available to us.
References
Baldy S. (1993), A generationdispersion model of am
bient and transient bubbles in the close vicinity
of breaking waves, J. Geophy. Res., 98, 18,277
18,293
Buckingham M.J. (1997), Sound speed and void frac
tion profiles in the sea surface bubble layer, Applied
Acoustics, 51, 225250
20
10
O
10
20
10
10
10
20
10
0
10
?n
0o
0 05 1 15 2 25 3
xx, (m)
Figure 6: Temporal variation of simulated void fraction
with time interval 0.2T
Carrica, P. M., Drew, D., Bonetto, F. and Lahey, R. T
(1999), A polydisperse model for bubbly twophase
flow around a surface ship, Int. J. Multiphase Flow,
25, 257305
Cox D. and Shin S. (2003), Laboratory measurements
of void fraction and turbulence in the bore region
of surf zone waves, J. Eng. Mech., 129, 11971205
Clift R., Viollet K.R. and Weber M.E. (1978), Bubbles,
Drops and Particles, Academic Press, New York,
USA
Deane G.B. and Stokes M.D. (1999), Air entrainment
processes and bubble size distributions in the surf
zone, J. Phys. Oceanogr, 29, 13931403
Deane G.B. and Stokes M.D. (2002), Scale dependence
of bubble creation mechanisms in breaking waves,
Nature, 418, 839844
Lamarre E. and Melville W.K. (1991), Air entrainment
and dissipation in breaking waves, Nature, 351,
469472
Lamarre E. and Melville W.K. (1994), Voidfraction
measurements and soundspeed fields in bubble
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
plumes generated by breaking waves, J. Acoust. Wu T.R. (211114), A numerical study of three
Soc. Am., 95, 13171328 dimensional breaking waves and turbulence effects,
Ph.D dissertation, Cornell University.
Lin P. and P.L.F Liu (1998), A numerical study of
breaking waves in the surf zone, J. Fluid Mech.,
359, 239264
Lopez de Bertodano, Lahey R.T. Jr. and Jones O.C.
(1994), Development of a k c model for bubbly
twophase flow, J. Fluids Eng., 116,128134
MartinezBazan C., Montanes J.L. and Lasheras J.C.
(1999), On the breakup of an air bubble injected
into a fully developed turbulent flow. Part 1.
Breakup frequency, J. Fluid Mech., 401, 157182
MartinezBazan C., Montanes J.L. and Lasheras J.C.
(1999), On the breakup of an air bubble injected
into a fully developed turbulent flow. Part 2. Size
PDF of the resulting daughter bubbles, J. Fluid
Mech., 401, 183207
Mori N., Suzuki T. and Kakuno S. (2007), Experimen
tal study of air bubbles and turbulence characteris
tics in the surf zone, J. Geophy. Res., 112, C'050)14
Moraga F.J., Carrica P.M., Drew D.A. and Lahey R.T.
Jr. (2008), A subgrid air entrainment model for
breaking bow waves and naval surface ships, Com
puters and Fluids, 37, 281298
Nadaoka K., Hino M. and Koyano Y. (1989), Structure
of the turbulent flow field under breaking waves in
the surf zone, J. Fluid Mech., 204, 359387
Shi F., Kirby J.T., Haller M.C. and Catalan P. (2008),
Modeling of surfzone bubbles using a multiphase
VOF model, Proc. 31st Int. Conf Coastal Engi
neering, Hamburg, 157169
Shi F, Kirby J.T. and Ma G. (2010), Modeling qui
escent phase transport of air bubbles induced by
breaking waves, Ocean Modeling, submitted
Sullivan P.P., McWilliams J.C. (2010), Dynamics of
winds and currents coupled to surface waves, Ann.
Rev. FluidMech., 42, 1942
Terrill E.J., Lada G. and Melville W.K. (2001), Surf
zone bubble populations, Pr ...... I.,,i,'l,,,,i,.. of
Acoustics, 23, 212219
Ting F.C.K. (2008), Largescale turbulence under a
solitary wave: Part 2 Forms and evolution of co
herent structures, Coast. Eng., 55, 522536
Troshko A.A. and Hassan Y.A. (2001), A twoequation
turbulence model of turbulent bubbly flows, Int. J.
Multiphase Flow, 27, 19652000
