7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
A 3D EulerianLagrangian Numerical Model for Sediment Transport
Berend van Wachem* Xiao Yut and TianJian Hsut
Department of Mechanical Engineering, Imperial College London, Exhibition Road, London, SW7 2AZ
t Civil and Environmental Engineering, University of Delaware, Newark, DE 19716, USA
b.vanwachem @imperial.ac.uk and thsu @udel.edu
Keywords: Lagrangian simulations, Sedimentation Transport
Abstract
In this study, a new approach for sediment transport utilizing 3D turbulence resolving EulerianLagrangian framework
is presented. The model employs slightly modified NavierStokes equations to determine the motion of the fluid phase.
The motion of the sediment phase is elucidated by a Lagrangian or Discrete Element Method (DEM), implying that
the individual trajectory of each particle is determined by approximating Newtons second law of motion. The forces
acting on each particle are gravity, the traction force of the fluid phase, and the force resulting from the interaction
with other particles. The traction force of the fluid phase is determined by an empirical equation which has been
validated for low particle Reynolds numbers. The total effect of the traction force summed over all particles is exactly
opposite to the interphase momentum transfer added to the fluidphase momentum equations. The force accounting
for interparticle collisions is based upon a socalled "softsphere" approach, following Tsuji (1994), which estimates
the force based upon the deformation of each particle resulting from the interaction with other particles. We have
performed simulation of sediment transport a 3 dimensional domain of approximately 50x30x50 particle diameters.
The two directions perpendicular to gravity are treated as periodic; i.e. the domain is fictitiously extended to infinity
in those directions. The flow is driven by an external pressure drop that corresponds to Shields parameter around 1.1.
Next to the EulerianLagrangian simulations, turbulenceaveraged 1DV EulerianEulerian simulations based upon
the kinetic theory of granular flow are performed as well. The 3D model results are compared and the results from
the EulerianLagrangian simulations are analyzed to obtain flow statistics, which can be used as a starting point for
EulerianEulerian closure model validation and development.
Introduction
The phenomena of sediment transport is defined as the
movement of solid particles due to a combination of
gravity acting on the particles (i.e. the sediment), the
intraparticle forces, and the movement of the fluid in
which the sediment is entrained. Sediment transport
due to fluid motion occurs in rivers, the oceans, lakes,
due to currents, waves, and tides. Sediment transport
in water involves complicated fluidparticle interactions
and intergranular interactions. A comprehensive
description of various mode of sediment transport,
namely, bedload and suspended load, requires a multi
phase flow approach. Existing modeling approach for
sediment transport utilize Reynoldsaveraged approach
in the EulerianEulerian framework (Dong and Zhang
2002; Hsu et al. 2004; Amoudry and Liu 2009) or
EulerianLagrangian framework (Drake and Calantoni
2001).
Modelling of sediment transport in rivers and coastal
areas is important in identifying erosion, changes
in fluvial and coastal morphology and the health of
ecosystem. Models have been successfully used to
study waveinduced sand transport under energetic
waves, i.e., sheet flow condition (Dong and Zhang
2002; Hsu and Hanes 2004; Yu et al. 2010). Through
Reynoldsaveraging, these models parameterized all the
scales of turbulence and related turbulencesediment
interaction. The coefficients involved in the turbulence
closure hence become highly empirical. Moreover, in
those models where sediment phase is modeled in an
Eulerian framework, kinetic theory of granular flow
is usually adopted for the closure of particle stresses.
Strictly speaking, kinetic theory becomes inappropriate
when sediment volume concentration becomes greater
than randomlosspacking ( .' and hence addi
tional bed treatment in the regime of enduring contact
in sediment transport is also necessary.
Simulations
Fluidphase governing equations
The equations for the fluidphase are derived analo
gously to the NavierStokes equations. However, they
include a local volume fraction (Anderson and Jackson
1967) and additional source terms accounting for the
presence of particles and the driving force. The govern
ing equations are the mass continuity and conservation
of momentum, which are given by equations (1) and (2)
respectively,
D(at f) d(aofpPfv )
at x1
D( 9fpfjV3) 9aftpfV3^) T9(af')
at + axi axi
Op
af O^
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
the Langrangian framework, where particle collisions
are modelled via the softsphere model proposed
by Tsuji et al. (1991b), which accounts for some
nonreversible deformation.
The interactions of particles with other particles
and walls are dynamic of nature. This is because
the particle movements are essentially defined by the
particleparticle interactions; particlewall interactions;
particlefluid interactions and/or body forces (Kuang
et al. 2008). The trajectories of individual particles are
considered (i.e. described in a Langrangian framework)
and Newton's 2nd law is solved for each individual par
ticle, accounting for the fluidparticle, particleparticle,
and particlewall interactions, and approximating the
integral with the Verlet algorithm.
Newton's 2nd law for the particles is written out
dvp P3 (Vf
dt ap 
Vp) + mpg
VpVPf
where af is the local volume fraction of the fluid phase,
,3 represents the local averaged reciprocal of the parti
cle timescale, or drag coefficient arising from the local
behaviour of the particles, Tf represents the source term
linear in the velocity field; S} represent the additional
source terms, which are used for driving the flow, and
T' the stress tensor of the fluid, given by
2 ) Ov'
3 O(3)
f (t vi vj \ (
a zi X3 01(h
where pf is the shear and Af represent the bulk viscos
ity of the fluid. As a first approach, the subgrid scale
stresses are not considered.
The simulations performed in this paper use a fi
nite volume scheme, a second order backward Euler
time discretisation and a second order accurate central
scheme to approximate the advective terms in the conti
nuity and momentum equations. The solving procedure
is made parallel with the MPI libraries. During the
simulations, the CFL number is kept constant at 0.4,
leading to a slight variation in timestep over time.
Particlephase governing equations
A discrete element model (DEM) proposed by Cundall
and Strack (1979) is used to model the particles. The
individual trajectories of the particles are determined in
N
+VpSdrive + S [Fpw + Fpp]
and for the rotational momentum
dw 
dt
where mp is the mass of the particle; Ip is the momen
tum of rotational inertia; Tp the torque of the particle;
w, is rotational velocity; vp is the translational veloc
ity; and 3 is the drag function as proposed by Wen and
Yu (1966), where the reciprocal of the Eulerian fluid
particle timescale is given by
3 apafPflvf
P d
4 dp
Vp 2.65
af
and CD represents the coefficient of drag for an individ
ual particle and af represents the fluid volume fraction.
The right hand side terms of equation (4) are outlined in
table 1.
Implementation of particle collisions
The particleparticle and particlewall interactions as
taken into account in this work are assumed to be
local; i.e. no longrange forces are included. For
establishing the nature of the interaction, each particle
pair could be interrogated. However, this would lead
to a scaling of the computational effort with N2, N
being the number of particles. Instead, a particlemesh
algorithm is adopted, in which each of the particles
is assigned a cell in the particle mesh based upon its
phases J f
+Tffv3+S + 1 p)p1]
Table 1: Terms of right hand side of equation(4)
Term Force type
a3Vf Vp. Drag force
mpg Body force due to gravity
VpVPf Force due to the pressure gradient
VpSdrive External driving force
Fp, Particlewall contact force
Fpp Particleparticle contact force
location. Using this particlemesh, each particle is
only tested for interaction against particles located in
the same or directly neighboring particle mesh cells.
Although there is some additional computational effort
and a slightly more complex algorithm required for
this approach, it leads to a scaling of N logN with the
number of particles, making it a lot more favourable for
large numbers of particles.
The particle collisions are modelled by the soft
sphere model as described by Cundall and Strack
(1979). In brief, this model uses a springdashpotslider
arrangement to describe the particle behaviour before,
during and after a collision. The damping coefficient
as introduced by Tsuji et al. (1991a) accounts for the
nonideal behaviour during the collision of the particles
due to irreversible plastic or viscoelastic deformation.
The normal and tangential contact forces are given
by the sum of forces due to the spring and dashpot.
From Hertzian contact theory the normal and tangential
contact forces are,
Fij = (kn6 i. G n)n (7)
Ftij = kt6t G (8)
where G is the velocity vector of particle i relative to
particle j, Get is the slit velocity vector at the contact
point. The subscripts n and t represent the normal and
tangential components respectively and 6 is the displace
ment, or overlap, of particles i and j during collision.
Three parameters are required by the softsphere model;
stiffness (k), damping coefficient (rI) and the coefficient
of friction (p). p is a well known empirical quantity,
stiffness and the damping coefficient must be estimated
using equations
/n = ac Mk6 (9)
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
mp. The relationship between a and the coefficient of
restitution is well defined by Tsuji et al. (1991b). The
spring constants are, based upon elastic deformation,
4 1 1
S 3 E +
h6g \ A
,, (2 i 2,j ri rj)_2 I
kt 8 G + rir 6 (12)
Gi G \ rirj
where r is the radius of the particles, a the Poisson's ra
tio; E the Young's modulus, G the shear modulus (given
by G = (iE ) and 6, the magnitude of the deforma
tion in the normal direction. Note as above, when a par
ticle reaches a wall, r, + oc, hence, r' + 1.
Simulation Parameters
The fluid properties used in the simulations are those
of water, with density p = 998 and viscosity v
1.00 x 10 3Pa s. The driving force of the fluid is deter
mined by the shields parameter. The shields parameter is
a dimensionless parameter that characterizes the type of
sediment transport regime. The Shields parameter, i.e.,
the nondimensionalized bed shear stress, is given by
T
O (13)
(ps pf) gds
where T is the bed shear stress, g is the gravitational ac
celereation, d, represents the particle diameter, and p,
and pf the solids and fluid density respectively. For sta
tistical steady flow considered in this study, the given
bed shear stress can be further related to pressure drop
in the streamwise direction:
9p T
(14)
Ox h
where h is the height of the domain. The driving force
is made equal to the pressure drop, i.e.
Drive
And this driving force is distributed over the local
fluid and particle phase weighted by the local volume
fraction.
For the simulations, sand (silicon dioxide) particles
were taken and assumed to have constant properties, see
Table 2
Results and Discussion
a/t = a 'kt (10)
where M = m and m is the mass of the particle.
Note that for wall collisions m,  oc, hence M
EulerianLagrangian Simulation Results
The simulations carried out with the Shields parameter
of 1.1. The resulting magnitude of the shear stress is
E r
ri ( rirj
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Table 2: Summary of the sand particle properties used
in computations.
Variable Value Units
d, 200 .106 m
ps 2650 kg/m3
Youngs Modulus 75 .107 Pa
Poisson ratio 0.17
Friction Coefficient 0.45
Restitution Coefficient 0.95
T = 3.89 Pa and the pressure gradient = 55.6 Pa.
According to field observation and laboratory ex
periments, Shields parameter greater than 1.0 represents
typical "sheet flow" condition where intense sediment
transport occurs, the bed is statistically flat (no bedforms
or ripples) with a concentrated layer of sediment moving
near the bed. A snapshot of the numerical simulation is
shown in Figure 1, where the particles are coloured by
their velocity.
The simulations predict the basic characteristics of sheet
flow. Sediment concentration decreases rapidly away
from the bed and all the transport occurs within 2 mm
near the bed. In this moving layer, the fluid velocity
also increases rapidly from 2 cm/s to about 20 cm/s.
Numerical model predicts the statisticallyaveraged
fluid velocity is almost identical to that of particle
velocity for relatively fine sediment in water considered
here.
EulerianEulerian Model Results
A EulerianEulerian twophase model for sediment
transport is adopted in this study. Turbulence
averaged twophase equations are simplified into one
dimensionalvertial (1DV) to model fullydeveloped
sheet flow sand transport driven by a steady or oscil
latory flow forcing (Hsu et al. 2004; Yu et al. 2010).
The eddy viscosity and ke closure for twophase flow
are adopted for fluid Reynolds stress and turbulent
suspension of sediment. Kinetic theory of granular
flow of Jenkins and Savage (1983) and empirically
modified closure for stress due to enduring contact are
utilized for the closure of particle stress. More detailed
model formulations and model applications to sediment
transport can be found in Hsu et al. (2"' 14); Yu et al.
(2010).
The EulerianEulerian model results compared to
the the EulerianLagrangian model results for the same
case are shown in Figure 2 for vertical profiles of
fluid velocity (solid curve in the left panel), sediment
Figure 1: A snapshot of the EulerianLagrangian sim
ulations, showing the individual particles. The particles
are coloured by their velocity.
concentration (solid curve in the middle panel) and fluid
turbulent intensity (2TKE), where TKE represents
the fluid turbulent kinetic energy (right panel). The
fluid stress consists of turbulent Reynolds stress and
viscous stress with the the viscous stress contribution
negligible in this case. EulerianEulerian model predicts
sediment suspended much higher in the water column
due to strong flow turbulence. Strong turbulence can
be observed from large turbulent intensity near the
bed. Some discrepancies with the EulerianLagrangian
model results can be observed. Most evidently,
EulerianEulerian model predicts more intense sediment
suspension throughout the water column and higher
mobility of sediment. Although it is likely that the
flow turbulence and granular temperature may be
overpredicted in the EulerianEulerian model due to
uncertainties in the closures, it is also possible that in
the EulerianLagrangian model, the flow turbulence
is underpredicted due to low resolution used in this
simulation. The underprediction of flow turbulence
can be seen from the almost linear fluid velocity profile
0 50 100 150 0 20 40 60 0 5 10 15
u (cm/s) vol cone (/) sqrt(2*TKE) (cm/s)
Figure 2: EulerianEulerian model results(line)
compared to the the EulerianLagrangian model re
sults(symbols) for the same case, for vertical profiles of
fluid velocity (solid curve in the left panel), sediment
concentration (solid curve in the middle panel) and fluid
stress (right panel, solid curve) and particle stress (right
panel, dotted curve).
away from the bed. On the other hand, Eulerian
Eulerian model also predicts almost identical fluid and
sediment velocities (not shown here), consistent with
EulerianLagrangian model results.
Conclusions
A new approach for sediment transport utilizing 3D
turbulence resolving EulerianLagrangian framework is
presented. Preliminary simulation results suggest the
main characteristics of sheet flow in the concentrated
regime of sediment transport is captured. However,
higher numerical resolution for fluid phase is required,
or subgrid turbulence closure need to be implemented in
order to capture flow turbulence and the resulting turbu
lent suspension of sediment.
Acknowledgements
Financial support for Dr. Berend van Wachem's travel
to ICMF2010 is partly provided by the Royal Academy
of Engineering. Financial support to Hsu and Yu is
provided by U.S. National Science Fundation (OCE
0644497).
References
L.O. Amoudry and P.L.F. Liu. Twodimensional, two
phase granular sediment transport model with applica
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
tions to scouring downstream of an apron. Costal Engi
neering, 56:693702, 2009.
T.B. Anderson and R. Jackson. A fluid mechanical de
scription of fluidized beds. Ind. Engng. Chem. Fundam.,
6(4):527539, 1967.
P.A. Cundall and O.D.L. Strack. A discrete numerical
model for granular assemblies. Geotechnique, 29:47
65, 1979.
P. Dong and K. Zhang. Intense nearbed sediment mo
tions in waves and currents. Coastal Engineering, 45:
7587,2002.
T.G. Drake and J. Calantoni. Discrete particle model
for sheet flow sediment transport in the nearshore. J.
Geophys. Res., 106:859868, 2001.
T.J. Hsu and D.M. Hanes. The effects of wave shape
on sheet flow sediment transport. J. Geophys. Res., 109:
C05025,2004.
TJ. Hsu, J.T. Jenkins, and P. LF Liu. On twophase
sediment transport: sheet flow of massive particles. Proc.
R. Soc. Lond. A, 460:22232250, 2004.
J.T. Jenkins and S.B. Savage. A theory for the rapid flow
of identical, smooth, nearly elastic, spherical particles. J.
Fluid Mech., 130:187202, 1983.
S.B. Kuang, K.W. Chu, A.B. Yu, Z.S. Zou, and Y.Q.
Feng. Computational investigation of horizontal slug
flow in pneumatic conveying. Ind. Eng. Chem. Res., 47:
470480,2008.
Y. Tsuji, T. Kawaguchi, and T Tanaka. Discrete particle
simulation of twodimensional fluidized bed. Powder
Technology, 77:7987, 1991a.
Y. Tsuji, T. Tanaka, and T. Ishida. Lagrangian numerical
simulation of plug flow of cohesionless particles in a hor
izontal pipe. Powder Technology, 71:239250, 1991b.
C. Y. Wen and Y. H. Yu. Mechanics of fluidization.
Chemical Engineering Progress Symposium Series, 62:
100111, 1966.
X. Yu, T.J. Hsu, and D.M. Hanes. Sediment transport
under wave groups: Relative importance between non
linear waveshape and nonlinear boundary layer stream
ing. J. GeophysicalResearch, 115:C02013, 2010.
