Paper No 7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Fine Particle Laden Flow Simulations with Simplified EulerianEulerian Approach
Celalettin E. Ozdemir*, TianJian Hsu** and S. Balachandar***
*University of Florida, College of Engineering, Civil & Coastal Eng.
University of Florida, Gainesville, FL, 32611, U.S.A.
Email: eozdemir@ufl.edu
**University of Delaware, College of Engineering, Civil & Environmental Eng.,
Center for Applied Coastal Research, DE, 19716, U.S.A.
Email: thsu@udel.edu
***University of Florida, College of Engineering, Mech. & Aerosp. Eng.
University of Florida, Gainesville, FL, 32611, U.S.A.
Email: sbalal @ufl.edu
Keywords: Particle laden flow, waveinduced sediment transport, oscillatory boundary layer
Abstract
Studying particleladen oscillatory channel flow constitutes an important step towards understanding practical applications
such as sediment transport in coastal environments and aerosol transport in respiratory airways, etc. This study aims to take a
step forward in our understanding of the role of turbulence on fine particle transport in an oscillatory channel and the back
effect of particles on turbulence modulation using an EulerianEulerian framework. The governing equations are solved by a
3D pseudospectral NavierStokes solver previously used for direct numerical simulation (DNS) of turbulent flows. The
EulerianEulerian modeling framework makes the twoway coupled problem computationally feasible without compromising
accuracy as long as particles are of small response time (Shotorban & Balachandar, 2006). As a first step, the instantaneous
particle velocity is calculated as the superposition of local fluid velocity and particle settling velocity while higher order
particle inertia effect neglected. Correspondingly, the only modulation of carrier flow is due to particleinduced density
stratification quantified by the bulk Richardson number, Ri. In the present investigation we fixed the Reynolds number whose
length scale is the Stokes boundary layer thickness to be ReA = 1000 and varied the bulk Richardson number over a range (Ri =
0.0, 104, 3x104, 6x104). The simulation results reveal critical mechanisms due to different degrees of particleturbulence
interaction. Essentially, four different regimes of particle transport for the given ReA are observed: i) virtually no turbulence
modulation in the case of very dilute condition, i.e., Ri 0; ii) slightly modified regime where slight turbulence attenuation is
observed near the top of oscillatory boundary layer. However, in this regime significant change can be observed in the
concentration profile with the formation of a lutocline; iii) regime where flow laminarization occurs during peak flow,
followed by shear instability during flow reversal. Significant reduction of oscillatory boundary layer thickness is also
observed; iv) complete laminarization due to strong particleinduced stable density stratification.
Introduction
We are motivated to identify main physical mechanisms in
waveinduced fine sediment transport as it is one of the
major processes shaping the coastal morphodynamics and
has a considerable impact on the fate of terrestrial sediment
in the coastal ocean (e.g., Traykovski et al. 2000; Harris,
Traykovski, & Geyer, 2005). Therefore, to understand the
role of turbulencesediment interactions in the resulting
sediment transport at the wave boundary layer scale,
carrying out highly resolved particle laden flow in
oscillatory channel constitutes a major step.
Waveresolving Reynoldsaveraged twophase models are
available to study the effect of waveinduced sediment
transport. The effects of turbulenceparticle interaction (i.e.,
twoway coupling, Conley et al., 2008; Hsu, Ozdemir, &
Traykovski, 2009) and particleparticle interaction
(fourway coupling, Hsu & Hanes, 2004; Yu, Hsu, &
Hanes, 2010) have been investigated. However,
Reynoldsaveraged models suffer from uncertainties in
turbulence closure such as the use of semiempirical
parameters in the transport equations and critical
assumptions, such as eddy viscosity, isotropy and
similarity relationships are used in the standard
twoequation models. On the other hand, existing
fullyresolved Direct Numerical Simulations (DNS)
studies (Spalart & Baldwin, 1989; Vittori & Verzicco,
1998; Salon, Armenio, & Crise, 2007) investigate
oscillatory boundary layer without particles. More recent
work of Chang & Scotti (2006) utilizes Large Eddy
Simulations (LES) for particleladen oscillatory boundary
layer with particle phase modeled to be completely passive
to carrier flow. According to Noh & Fernando (1991),
assuming particles to be completely passive is valid only
under extremely dilute condition. As we shall discuss later,
even at low concentration where particleparticle
interaction can be ignored, the back effect of suspended
particles on the carrier phase can be significant. In
particular, the magnitude of carrier phase turbulence is
modulated and the nature of twophase flow turbulence,
such as the degree of isotropy and kinetic energy budget,
are altered. This change in carrier phase turbulence in turn
Paper No
influence the resulting particle suspension/deposition
processes in an oscillatory flow. The unsteady nature of the
oscillatory particulate flow presents a challenge compared
to steady flow or tidal variations. Under ideal sinusoidal
forcing, the flow velocity varies between positive and
negative peaks with flow reversal in between. If the
instantaneous flow velocity becomes sufficiently large,
then perturbations in the flow start to grow during the
acceleration phase and depending upon the extent of the
acceleration phase the flow might become fully turbulent.
During deceleration, the grown perturbations may fully
decay and the flow might laminarize. Or depending on the
maximum velocity that flow reaches, the perturbations
may only partially decay and the flow can remain turbulent
throughout the entire oscillatory cycle. This unsteady
nature of the flow has motivated prior studies to address
the onset of turbulence (Hino, Sawamoto, & Takasu, 1976)
and evolution of vortex structures (Sarpkaya, 1993). These
studies have classified four flow regimes based on the
Stokes Reynolds number: (1) laminar regime (ReA < 200);
(2) disturbed laminar regime (200< ReA < 400); (3)
intermittently turbulent regime (400< ReA < 1200); and
(4) fully turbulent regime (ReA >1200). The Stokes
Reynolds number ReA is defined as:
ReA = (1)
Vf
where U is the dimensional maximum free stream
velocity, A is the dimensional Stokes boundary layer
thickness defined as = 2v with oa being the
dimensional angular frequency of oscillatory forcing, and
Vs is the kinematic viscosity of the fluid. In the laminar
regime, the flow stays laminar over the entire wave cycle
and in the fully turbulent regime it stays turbulent
throughout the wave cycle. In the perturbed laminar regime,
the flow is perturbed during acceleration phase and
laminarizes during deceleration without ever becoming
fully turbulent. Finally, in the intermittently turbulent
regime the flow becomes turbulent during acceleration and
tends to decay during deceleration without complete
laminarization.
Direct numerical simulations of oscillatory channel flow in
a clear fluid was initiated by Spalart & Baldwin (1989) and
continued in the works of Akhavan, Kamm & Shapiro
(1991), Vittori & Verzicco (1998) and Costamagna, Vittori
& Blondeaux (2003). Recently, Ozdemir, Hsu, &
Balachandar (2010) performed numerical simulations in
the intermittent turbulent regime for very dilute particle
concentration, where particles can be assumed to be
passive and not influence the carrier flow. They use a
higherorder accurate pseudospectral flow solver and
employed the same domain size and resolution as Spalart
& Baldwin (1989). These simulations are of high fidelity
since they use a large computational domain, a fine grid
resolution and spectrally accurate numerical methodology.
Ozdemir, Hsu, & Balachandar (2010) observe randomly
distributed energetic vortices present at the maximum free
stream velocity. These vortices start to lose their energy
during flow reversal and are uplifted towards the upper
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
portions of the channel, where they manifest themselves as
random spikes shown in concentration and vorticity
isocontours. These vortices then transform into organized
long streaks at the early stages of acceleration.
Correspondingly, the peak turbulent fluctuation occurs
close to the bottom boundary when free stream velocity
reaches its maximum. Then, turbulent intensity reduces
and the peak shifts away from the bottom boundary.
The problem of turbulence modulation by particles in an
oscillatory flow is far more challenging. A discussed above,
even in the limit of clear fluid, depending on the Reynolds
number, turbulence may be limited to only part of the wave
cycle. Consequently, the interaction between suspended
particles and carrier flow turbulence can be complex and
highly variable over the wave cycle. In this study, we shall
investigate the following issues of fine particle transport in
an oscillatory boundary layer using twoway coupled
numerical simulations and discuss their relevance to
existing field and laboratory observations: (1) The
formation of lutocline in an oscillatory boundary layer, (2)
the observation of instability and bursting events during
flow reversal in stratified turbulence, and (3) complete
laminarization due to particleinduced stratification.
Numerical Scheme
We idealize the case of waveinduced turbulence at the
bottom boundary layer in a typical coastal settling. Surface
waves at continental shelf have a wave period of about 10
sec and the corresponding Stokes boundary layer thickness
is =1.8 mm. Here, we consider a Reynolds number (see
equation (1)) of Re, =1000, which corresponds to an
oscillatory flow velocity amplitude of 0.56 m/s. This value
is similar in magnitude to that observed at continental shelf
where fluid mud is suspended in the wave boundary layer
(e.g., Traykovski et al. 2000; Traykovski, Wieberg, &
Geyer 2007). We consider a suspension of fine silt in
marine environment. Typical silt size of 24 gm with a
specific gravity of 2.65, results in a still fluid settling
velocity of about 0.5 mm/s. The particle time scale,
defined as ppd2 /18,uf can be computed to be 8.5x105 s,
where p is particle density, d is particle diameter and
uf is dynamic viscosity of water.
For the current case of, Re =1000 as will be shown
below, the flow is turbulent provided the back effect from
the suspended particles is small. The flow thus presents a
range of time scales from the integral to the Kolmogorov
scale. The time scale of the Stokes layer (A/U0) is more
than an order of magnitude larger than the particle time
scale. From the numerical results we have verified that the
particle time scale is also smaller than the Kolmogorov
time scale. It can thus be established that the particles are
sufficiently small that their Stokes number, defined as the
ratio of particle to fluid time scale, is less than 1. As a
result, here we ignore the inertial effect of particles (Ferry
& Balachandar, 2001) and simply define the
nondimensional particle velocity Up as the sum of fluid
velocity, uf and particle settling velocity, V as:
Up = U +V
Paper No
where superscripts f and s stand for the fluid and particle
phases, respectively. Substituting the above algebraic
relationship in the EulerianEulerian twophase equations
and applying Boussinesq approximation valid for relatively
small concentrations, the resulting governing equations can
be greatly simplified (e.g., Cantero, Balachandar, & Garcia,
2008). The continuity (equation (3)) and momentum
equations (equation (4)) of the fluid phase can be written
as:
V.U= 0 (3)
DUf 2 cos( 2 t)e,VP RiCe2+ V2(Uf) (4)
Dt Re, Re, Re,
These equations are nondimensionalized by length,
velocity, and time scales selected to be the Stokes
boundary layer thickness ( A ), maximum free stream
velocity (C o), and A / respectively.
In the present problem of an oscillatory channel flow, the
time variation of the dimensional free stream velocity is
given as follow
U=0 ..... )eI (5)
Where e and e2 correspond to unit vectors along the
streamwise and wallnormal directions, respectively. Far
from the bottom boundary, well outside the oscillatory
boundary layer, the momentum equation can be simplified
for inviscid flow to yield the mean streamwise pressure
gradient to be
oP 0 (6)
S n .... (..6)
The above equation upon nondimensionalization, along with
the definition of Stokes layer thickness and Reynolds
number, reduces to the first term on the righthandside of
equation (4). The second term corresponds to the fluctuating
pressure gradient. The third term on the righthandside
represents the back effect of the particle on the fluid phase
via particleinduced density stratification, which is only
effective in the vertical direction. Here the bulk Richardson
number Ri is defined as:
Ri ('p)g (7)
pil2o
in which, p is fluid density and
volumeaveraged concentration over
computational domain, i.e.
c JC(x,y,z)dV
V
C is
the e
the
entire
(8)
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Here, V is the volume of the domain. In the present
formulation, C is normalized by C. Consistent with
equation (2), in the present formulation the momentum
equation does not include frictional or drag forces arising
from the inertial effect of particles. Hence, preferential
accumulation of particles is not expected to be observed in
our simulations, which is appropriate for particles of very
small response time.
Finally, the particle concentration is calculated by the
conservation of particle mass:
(C)+V[(US+Ve2)C] 1 V2C
t s Re. Sc
The term on the righthandside represents an effective
diffusive flux. Here, the Schmidt number is Sc = V /lK,
where K is the effective diffusivity of the particles. The
governing equation of particle motion, when applied in the
Lagrangian framework, contains no diffusive effect, i.e.
Sc oo However, in the present simulations, a diffusive
term is required not only for numerical stability but also to
account for subscale particle motion through subgrid
interactions. Birman, Martin, & Meiburg (2005) show that
the effect of Sc within the range of (0.2, 5.0) to be quite
insensitive in moderately energetic gravity flows. Similarly
Bonometti & Balachandar (2008) have investigated the
effect of Sc in the range of [1,oo) for gravity currents with
similar governing equations and observed Schmidt number
effect to be small in high Reynolds number flows. It is
observed that in vigorous flows, mixing is dominated by
advection and the effect of molecular mixing represented
by diffusion to be quite small. We have selected Sc to be
0.5 in the present simulations.
The velocity boundary conditions at both the top and
bottom boundaries are noslip wall boundary conditions.
Periodic boundary conditions are imposed at the lateral
boundaries of the computational domain. For particle
concentration, the following relation is imposed at both top
and bottom boundaries of the computational domain:
ac 1 2C (10)
V, (10)
s y Re, Sc ay2
This boundary condition imposes no net transport of
particles across the top and bottom boundaries. Periodic
boundary conditions for particle concentration are applied
along the streamwise and spanwise directions. These
boundary conditions in turn imply that the total volume of
particles remains constant in the domain throughout the
computation. This integral property has been confirmed as
a verification test in the present computations.
The size of the computational domain is selected to be the
60A x 60A x 30A along the streamwise, vertical and
spanwise directions, respectively. Thus, the distance from
the bottom boundary to the center plane of the channel is
chosen to be 30 times the Stokes layer thickness. The
selected domain size is exactly the same as that of Spalart
Paper No
& Baldwin (1989) and their domain is the largest and the
most conservative among the different numerical studies.
By examining the two point correlation function Ozdemir,
Hsu, & Balachandar (2010) showed that this domain size is
sufficiently large to capture the main turbulent
characteristics at ReA = 1000 in the case of clear fluid.
The flow field is timeadvanced using CrankNicolson
scheme for the diffusion terms. The nonlinear advection
terms are dealised with the Arakawa method and advanced
with a thirdorder implicit RungeKutta scheme due to its
low storage requirement. More details on the
implementation of this numerical scheme can be found in
Cortese & Balachandar (1995). The time step chosen for
the simulations is 1/96000 of the wave period, which
maintains the CFL number to be below 0.5. Thus, each
wave is resolved with almost 105 time steps. The
simulation was typically evolved about 6 cycles or more in
order for the initial transients to decay and the turbulence
statistics to reach the stationary state.
Results and Discussion
Simulation results for the four different particle
concentrations (Ri = 0.0, 104, 3x104, 6x104) present
rather different outcomes of particleturbulence interaction,
which shall be qualitatively identified and described in
terms of particle concentration (Figure 1 and 2). Turbulent
structures shown here do not include the case of Ri =
6x104 as the turbulence is fully suppressed by the
suspended particles. Hence, the flow is completely
laminarized and there are no vortex structures.
The effect of increasing particleinduced density
stratification is observed in Figure 1 and Figure 2 where
isosurface of constant concentration is plotted at three
different phases spanning half a wave cycle. These three
instances are selected from the accelerating phase () = 7/3),
one from the flow reversal () = 7/2), and from the
decelerating phase () = 21/3). For the case of Ri = 0.0,
particles do not affect the carrier phase and thus turbulence
corresponds to that of clear fluid discussed by Spalart &
Baldwin (1989). Let us first discuss this clear fluid limit.
At ) = 7/3 the far field (in the present case the
midchannel) a dense pack of nearwall vortical structures
can be observed close to the bottom boundary. The
concentration isocontours observed at ) = 7, when the far
field reverse flow reaches its peak, are qualitatively similar
although not shown here owing to temporal periodicity. As
the free stream decelerates () = 7/3) the vortex structures
thin out and move away from the bottom boundary. This
process of outward migration continues during flow
reversal () = 7/2) and even past flow reversal () = 21/3).
At ) = 5jr/6 when the reverse flow is accelerating in
strength, only few streaky vortical structures can be
observed. These streaky vortices are close to the bottom
boundary and it appears that between ) = 21/3 and ) =
57/6 the turbulence that was observed away from the
boundary at ) = 21/3 decays rapidly and the near wall
structures observed later at the accelerating phase grow
rapidly into mature turbulence when the reverse flow peaks
at ) = 7t.
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
The isocontours of concentration for the case of Ri = 104
shown in Figure 1 are statistically similar to those for the
passive case. Similar observations can be made for the
isosurface of particle concentration (Figure 1) with the
exception that the isosurface level for Ri = 104 is chosen
to be two times the volumeaveraged concentration, which
is slightly higher than that for the passive case. As will be
discussed in detail when we present the mean
concentration profiles, increasing Richardson number from
0 to 104 has a strong influence on mean particle
concentration. Since particle concentration in the bottom
half of the channel increases substantially, a
correspondingly higher contour level must be chosen to
extract similar turbulent structures.
For the case of Ri = 3x104, concentration isocontours
show a drastically different behavior than those observed
in the passive and Ri = 104 cases. First of all, the flow
shows a very quiescent behaviour at peak flow and during
the entire deceleration phase. This can be explained by
turbulence suppression due to particleinduced stable
density stratification. At flow reversal () = 7c/2 in Figure 2)
the effect of incipient instability in terms of small
amplitude waves can be observed in the isosurface of
particle concentration in Figure 2. These instabilities
tend to grow from an earlier nearlaminar condition () =
7/6, 7/3) and reach peak amplitudes in the early
acceleration phase around ) = 2/13. The instability appears
as well organized spanwise vortex rollers formed as a
result of shear instability. Signature of secondary
instability can be seen as spanwise waviness (see ) = 27/3
in Figure 2). While relatively clean, small amplitude, and
nearly twodimensional waves can observed at ) = 7/2, a
short time later at ) = 27/3 a more complex wave train
with spanwise disturbances resembling wavebreaking
formation is observed. At later stages of acceleration and
the isosurface of concentration show decay of the
instabilities and only a remnant of the vortex structure and
its impact on particle concentration can be observed. Note
that the concentration contour shown in Figure 2 is ten
times the volumeaverage concentration, which is much
higher than those of previous cases of lower Ri. For Ri =
3x104 the flow largely laminarizes and as a result, mean
particle concentration close to the bottom wall
significantly increases. In order to extract the effect of
instability, which is located close to the bottom wall, it is
necessary to choose a larger concentration value for the
isosurface of concentration field.
These instability features for the case of Ri = 3x104
observed in Figure 1 resemble the shear instability noted
by Strang & Fernando (2001). They investigate the
behavior of these interfacial shear instabilities with respect
to changes in the bulk Richardson number. While for small
bulk Richardson number they observe pure
KelvinHelmholtz (KH) billows, with the increase in Ri
they observe a combination KH billows plus H61mb6e
waves. In the present problem it is difficult to separate out
the exact characteristic of these structures and the
associated instabilities as they occur in a transient manner
over only a short duration near flow reversal. This is unlike
the case considered in Strang & Fernando (2001) where
the instabilities occur in a quasisteady manner and the
effect of stratification can be considered almost constant.
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Interaction between carrier flow and particles for Ri = 0,
104 and 3x104 are further analyzed in terms of
statisticallyaveraged particle concentration and velocity
statistics. Statistics are obtained by averaging over
horizontal planes and also over time instants of identical
phase. These averages will be referred to as
"phaseaverage" and denoted by < >.
The averaged concentration and streamwise velocity
profiles at six different distinct phases during the wave
cycle for the cases of Ri = 0, 104 and 3x104. Note that the
concentration scale for the Ri =3x104 case is different
from that of the other two cases in order to accommodate
the rapid increase in concentration towards the bottom wall.
The averaged streamwise velocity profiles for Ri = 0 and
104 are almost identical. The corresponding
rootmeansquare (RMS) turbulent velocities in the
streamwise direction at the six different phases are
presented in Figure 1 and Figure 2. Again, it can be
observed that the turbulent intensities for Ri = 0 and 104
cases are of similar magnitude over much of the channel.
In particular, there is little difference in turbulent
intensities close to the bed between y = 0 and y = 10.
However, closer to the midchannel (from y = 12 to y = 30),
although the overall magnitude of turbulent intensities
become small, the turbulence intensity for Ri = 104 ranges
only between 10 to 50 percent of that in the Ri = 0 case.
The effect of reduced turbulence above y = 15 can be
clearly seen in the mean concentration profiles shown in
Figure 1. For the passive case, the concentration profiles
are more or less linear and particles are wellmixed over
the entire channel, except for a thin layer of large
concentration gradient very close to the bottom. On the
other hand, for Ri = 104 case, it can be observed that the
particles are entrapped in the lower half of the domain and
the concentration profile shows a clear feature of
"lutocline" (sharp negative concentration gradient). This
shouldershaped concentration profile, called lutocline, is
predicted by Noh & Fernando (1991) for fine particle
suspended in oscillating grid flow due to the damping of
turbulence kinetic energy by particleinduced stable
density stratification. It is also observed in the laboratory
experiment with similar oscillating grid setup (Huppert,
Turner, & Hallworth, 1995). In the present simulation, the
formation of lutocline can be clearly seen for fine particle
transport in the oscillatory boundary layer. Closer to the
midchannel, turbulence and mixing are suppressed due to
the stabilizing effect of the large concentration gradient.
Suppressed turbulence further prevents upward migration
of particles, which would otherwise occur driven by
turbulence. In this case, since suppression of turbulence
occurs near the mid channel where the overall magnitude
of turbulent kinetic energy is already small, turbulence
modulation does not make significant difference in the
mean velocity profile.
L ..3 0
10g i
_______________ _, __ _o_ __ __ n_ I _________________
Figure 1: Comparison simulation results for passive and Ri
= 1 x 104 case at column (A) ) = 7/3, at column (B) ) = 7/2,
and at column (C) ) = 2j/3. The phase information is given
at the first row, the isocontours of concentration are given
for passive case and Ri = 1 x 104 at the second and the third
rows, respectively, averaged concentration, velocity, and
RMS of streamwise fluctuating velocities are given at rows
3, 4, and 5, respectively.
20L 203
()6. /2 C
S so = 1
40 5 c =0g\l 0
10 o0
02 04 06 004 0 (2 02 00 04 0
) = r/2, and at column (C) ) = 2r13. The phase information
is given at the first row, the isocontours of concentration
are given for passive case and Ri = 1 x 104 at the second
and the third rows, respectively, averaged concentration,
velocity, and RMS of streamwise fluctuating velocities are
given at rows 3, 4, and 5, respectively.
Numerical simulations presented here qualitatively identify
Paper No
Jo
20  C' 10 I =l ai
0
07 02 04g 06
a ,
D204 D 02
Paper No
several regimes of fine particle transport and their
corresponding dynamics in oscillatory channel flow. These
regimes are consistent with prior laboratory and field studies
on fine sediment transport in wave boundary layer (e.g.,
Traykovski et al. 2000, 2007; Lamb et al. 2004; Foster et al.
2006). Essentially, four different regimes of particle
transport for the given Rea are observed: i) virtually no
turbulence modulation in the case of very dilute condition,
i.e., Ri ~ 0; ii) slightly modified regime where slight
turbulence attenuation is observed near the top of oscillatory
boundary layer. However, in this regime significant change
can be observed in the concentration profile with the
formation of a lutocline; iii) regime where flow
laminarization occurs during peak flow, followed by shear
instability during flow reversal. Significant reduction of
oscillatory boundary layer thickness is also observed; iv)
complete laminarization due to strong particleinduced
stable density stratification.
Conclusions
Highly accurate numerical simulations have been carried
out for fine particle transport in oscillatory boundary layer.
Twoway coupled simulations adopting
EquilibriumEulerian method is used to simplify particle
phase formulation appropriate for small particle response
time. In this study, we further neglect higher order particle
inertia terms and the resulting turbulence modulation is
only due to particleinduced density stratification. Flow
turbulence is fully resolved at a Stokes Reynolds number
of ReA= 1000. We present four cases with nondimensional
particle settling velocity of V, = 9x104 at different bulk
Richardson number, i.e. Ri = 0, 104, 3x104 and 6x104,
representing various degree of particleinduced stable
density stratification.
At Ri = 104, which corresponds to a near bed sediment
concentration of 10 g/1 for the Reynolds number and
settling velocity selected, sediment induced stable density
stratification attenuates flow turbulence and reduces
mixing of sediment near the top of the wave boundary
layer, which gives the formation of a sharp concentration
gradient, i.e., lutocline. This test case also gives a
qualitative estimate on the minimum concentration for
which particle cannot be considered as passive to the
carrier flow turbulence. On the other hand, turbulence near
the bed is not affected by the sediment and mean flow
velocity is almost identical to clear fluid condition. At Ri =
3x104 (near bed concentration about 50 g/l), particle
induced density stratification is strong enough to attenuate
turbulence in the entire wave boundary layer. Flow tends to
be laminarized and mean flow velocity and concentration
profiles become similar to laminar solution. However, flow
instability, which can be seen clearly both in coherent
vortical structure and instantaneous particle concentration,
occurs during flow reversal which trigger large fluctuation
production and particle suspension that last about onethird
of the wave period. Finally at Ri = 6x104 (near bed
concentration more than 100 g/l), oscillatory boundary
layer flow become completely laminarized at all times due
to intense sedimentinduced stable density stratification.
At Ri = 3x104 or greater, turbulence is attenuated across
the entire oscillatory boundary layer and hence the
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
oscillatory boundary layer thickness and friction are
reduced due to particleinduced stable density stratification.
In most cohesive sediment transport studies, when
sediment concentration is greater than 0(100) g/l,
theological stresses caused by interactions among the floc
aggregates and interstitial water need to be considered.
Rheological stress gives enhanced effective viscosity that
can further increase energy dissipation and friction factor.
These nonNewtonian flow features shall be studied as
future work. This study identifies several distinct regimes
of particleladen oscillatory boundary layer according to
different magnitudes of bulk Richardson number (or
sediment concentration). Although the exact value of Ri for
these regimes must depends on Reynolds number and
settling velocity, the characteristics of these regimes are
useful to explain wavedriven cohesive sediment transport
process in coastal environments. Future work shall be
devoted to explore these characteristics for different
settling velocity and Reynolds number. Finally, more
detailed laboratory and field measurements shall be carried
out to verify several interesting observations revealed by
the present numerical simulations, such as flow instability
that occurs during flow reversal.
Acknowledgements
This study is supported by the U.S. Office of Naval
Research (N000140910134) and National Science
Foundation (OCE0913283) to University of Delaware and
by U.S. Office of Naval Research grant
(N000140710494) to University of Florida. This work is
partially supported by the National Center for
Supercomputing Applications under OCE70005N and
OCE80005P utilized the NCSA Cobalt and PSC Pople.
Also, the authors acknowledge the University of Florida
HighPerformance Computing Center for providing
computational resources and support that have contributed
to the research results reported in this paper.
References
Akhavan, R., Kamm, R. D. & Shapiro, A. H. 1991 An
investigation of transition to turbulence in bounded
oscillatory Stokes flows. Part 2. Numerical
simulations. J. Fluid Mech. 225, 423.
Birman, V., Martin, J.E., Meiburg, E. 2005 The
nonBoussinesq lockexchange problem. Part 2.
Highresolution simulations. J. Fluid Mech. 537,
125144.
Bonometti, T. & S. Balachandar 2008 Effect of Schmidt
number on the structure and propagation of density
currents Theor. Comput. Fluid Dyn. ,22, 341361 DOI
10.1007/s0016200800852
Chang Y. S. & A. Scotti 2006, Turbulent convection of
suspended sediments due to flow reversal, J. Geophys.
Res., 111, C07001, doi:10.1029/2005JC003240.
Costamagna, P., Vittori, G. & Blondeaux, P. 2003 Coherent
structures in oscillatory boundary layers. J. Fluid
Mech. 474, 1.
Conley, D. C., Falchetti, S., Lohmann, I. P., Brocchini, M.,
2008 The effects of flow stratification by
Paper No
noncohesive sediment on transport in highenergy
wavedriven flows, J. Fluid Mech. ,610, 4367.
Cortese, T. & Balachandar, S., 1995 High performance
spectral simulation of turbulent flows in massively
parallel machines with distributed memory.
International Journal of Supercomputer Applications
9 (3), 187204.
Ferry J. & S. Balachandar, 2001 A fast Eulerian method for
disperse two phase flow, Int. J. Multiphase Flow 27,
1199.
Harris C.K., Traykovski, PA., Geyer, W.R., 2005 Flood
dispersal and deposition by nearbed gravitational
sediment flows and oceanographic transport: a
numerical modeling study of the Eel River shelf,
northern California. J. Geophys. Res., 110(C9),
(2512516).
Hino, M., Sawamoto, M. & Takasu, S., 1976 Experiments
on transition to turbulence in an oscillatory pipe flow,
J. Fluid Mech. 75, 193.
Hsu, TJ & Hanes, D.M., 2004 The effects of wave shape
on sheet flow sediment transport, J. Geophys. Res.,
109(C5), C05025, doi:10.1029/2003JC002075.
Hsu TJ., C. E. Ozdemir, & P. A. Traykovski 2009 High
resolution numerical modeling of wavesupported
sediment gravitydriven mudflows J. Geophys. Res.,
114, C05014, doi: 10.1029/2008JC005006.
Huppert, H. E., J. S. Turner, & M. A. Hallworth, 1995
Sedimentation and entrainment in dense layers of
suspended particles stirred by an oscillatinggrid. J.
Fluid Mech., 289, 263 293.
Lamb, M. P., D'Asaro, E. & Parsons, J. D. 2004 Turbulent
structure of highdensity suspensions formed under
waves. J. Geophys. Res. 109, C12026.
Noh, Y., & Fernando, J. S., 1991 Dispersion of suspended
particles in turbulent flow, Physics of Fluids A, 3(7),
17301740.
Ozdemir, C. E., T.J. Hsu, S. Balachandar, 2010 Simulation
of fine sediment transport in oscillatory boundary
layer, Journal of Hydroenvironment Research,3(4),
247259.
Salon, S., V. Armenio, & A. Crise 2007 A numerical
investigation of the Stokesboundary layer in the
turbulent regime, J. Fluid Mech. 570, 253.
Sarpkaya, T. 1993 Coherent structures in oscillatory
boundary layers. J. Fluid Mech. 253, 105.
Shotorban and Balachandar, 2006 Particle concentration in
homogeneous shear turbulence simulated via
Lagrangian and equilibrium Eulerian approaches,
Phys. of Fluids, 18, 065105.
Spalart, P. R. & Baldwin, B. S. 1988 Direct simulation of a
turbulent oscillating boundary layer. Turbulent Shear
Flows 6 Springer.
Strang, E.J. & Fernando, H.J.S., 2001 Entrainment and
mixing in stratified shear flows. Journal of Fluid
Mechanics 428, 349386
Styles, R. & S. M. Glenn, 2000: Modeling stratified wave
and current bottom boundary layers in the continental
shelf, J. Geophys. Res., 105, 2411924139.
Traykovski, P., W. R. Geyer, J. D. Irish, & J. F Lynch,
2000 The role of waveinduced fluid mud flows for
crossshelf transport on the Eel River continental
shelf, Cont. Shelf Res., 20, 2113 2140,
doi:10.1016/S02784343(00)000716.
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Traykovski, P., P.L. Wieberg, & W. R. Geyer, 2007
Observations and modeling of wavesupported
sediment gravity flows on the Po prodelta and
comparison to prior observations from the Eel shelf.
Cont. Shelf Res., 27, 375399.
Vittori G. & Verzicco R., 1998 Direct simulation of
transition in an oscillatory boundary layer. J. Fluid
Mech. 371, 207232.
Yu, X, Hsu, T.J., & Hanes, D. M., 2010 Sediment
transport under wave groups the relative importance
between wave shape and boundary layer streaming, J.
Geophys. Res., 115, C02013,
doi:10.1029/2009JC005348.
