7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
A 3D TwoPhase Numerical Model for Sediment Transport
Julien Chauchatt,* Malika OuriemiP,* Pascale Aussillous*,
Marc M6dale* and Elisabeth Guazzelli*
IUSTI CNRS UMR 6595 Polytech'Marseille AixMarseille Universit6 (Ul), France
t Laboratoire des Ecoulements G6ophysiques et Industriels, Universit6 Joseph Fourier, INPG, CNRS, Grenoble, France
SIFPLyon, Solaize, France
julien.chauchat@polytech.univmrs.fr and elisabeth.guazzelli@polytech.univmrs.fr
Keywords: Sediment transport, twophase model, granular rheology, numerical model
Abstract
We have developed a three dimensional numerical model based on the twophase equations to study the bedload
transport ( I.iiiIli.ii and M6dale (2009). We have considered two formulations of the model based on a two fluid or a
single mixedfluid description. The governing equations are discretized by a finite element method and a penalisation
method is introduced to cope with the incompressibility constraint whereas a regularisation technique is used to deal
with the viscoplastic behaviour of the granular phase. The accuracy and efficiency of the numerical models have
been compared with the analytical solution of Ouriemi et al. (2009a). It turns out that one must takes a smaller
regularisation parameter (one order of magnitude) in the mixedmodel than in the twofluid one for a comparable
accuracy. Using an Arc Length Continuation algorithm coupled with this model we have investigated the evolution of
the solution in terms of the height of the flowing granular layer and the particle flux against the longitudinal pressure
gradient. The results of these simulations are in good agreement with the analytical solution in the 2D case and
threedimensional computations have been carried out in a square and circular crosssection ducts showing that the
effect of the geometry is nontrivial.
Introduction
The transport of sediment or more generally the trans
port of particles by a fluid flow is a problem of major
importance in geophysical flows such as coastal or river
morphodynamic or in industrial flows with the hydrate
or sand issues in oil production and granular transport in
food or pharmaceutical industries. This problem has been
extensively studied in the literature since the middle of the
twentieth century but poorly understood actually (Einstein
1942; MeyerPeter and Muller 1948; Bagnold 1956; Yalin
1963).
Recently, Ouriemi et al. 12', a ,., have proposed a two
phase model describing the bedload transport in lami
nar flows that allows to incorporate more physics than
in previous modelling based on particle flux or erosion
deposition approaches. This twophase model is based
on a Newtonian rheology for the fluid phase and a fric
tional rheology for the particulate phase (p(I) Forterre
and Pouliquen (2008) or Coulomb friction) while the fluid
particle interaction is assumed to follow a Darcy law. This
approach allows to predict the threshold of motion for the
particle phase and give a description of the flow of inside
the flowing granular layer. Away from the threshold of
motion, a simpler analytical model for the particle flux
is obtained which gives a quite satisfactory description of
experimental observations of bedload transport in pipe
flows Ouriemi et al. (2009a).
Based on this theoretical model we have developed a
3D numerical model that allows to simulate bedload trans
port in 2D or 3D configurations (( Ii.iiili.ii and M6dale
2009). It is restricted to the cases where the granular bed
does not change its shape in the course of time, conse
quently ripples and dunes formation are beyond the scope
of this paper. We consider a mixedfluid and a twofluid
formulation, the models equations are discretized by a fi
nite element method and a penalisation method is used to
impose the incompressibility constraint. The granular rhe
ology is analoguous to a viscoplastic behaviour especially
by the existence of a yield stress with the particularity that
this yield stress depends on the depth inside the granular
layer. A regularisation technique is used to deal with the
viscoplastic rheology and is shown to give satisfactory re
sults.
Twophase model
The present model is based on Jackson (2000, 1997) aver
aged equations using the closures developed by Ouriemi
et al. (2009a). These equations are summarized hereafter
in dimensionless form using the following scaling: the
length is scaled by H, the channel height (see figure 1),
and the stresses is scaled by ApgH, and therefore the time
is scaled by l/ApgH where Ap pp pf. The problem
is expressed in terms of the solid volume fraction y, the
mixture velocity u" and the particulate velocity up.
Twofluid model
V. (u"
V. up
H3 DD
H3 DuD
Ga RpD
d3 DL
Vp + V 
VP V.
H 2 Up
K +A 1Ap1 7
Vp Vp7f
+V. 7 + v. )
(1 _ )H72 t $
K V II 
(1)
In these equations, k Vu + (V.. F with k = m
or k p, R, pf/pp represents the density ratio and
Ga = d3pfApg/r2 is the Galileo number where d is the
particle diameter. The Galileo number is a Reynolds num
ber based on the settling velocity of particles.
Mixedfluid model
1 R3 )D uL
C., H 1+ R) D
D
Vp P V pp + PM
Ap II ( II
+V. T I +
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
with ps the internal friction coefficient and the permeabil
ity K (1 )3d2
ity K = with kOK w 180 (Happel and Brenner
1973).
Numerical model
The Finite Element Method (FEM) is based on the dis
cretisation of the variational formulation associated with
equations (1) for the twofluid model and on equation (2)
for the mixedfluid model.
The twofluid formulation is based on the solution of
the system (1) in a weakly coupled way. The algorithm
associated with the mixedfluid formulation is based on
the solution of the system (2). The nonlinearities in the
governing equations are solved by a NewtonRaphson al
gorithm by taking the first variation of the variational for
mulations.
The specific issue raised by the previous twophase model
lie in the calculation of the frictional stress. Following
(Jop et al. 2006), the frictional stress can be written as
TP = lP
with I. = i'/I P . The particulate viscosity di
verges as the particulate shear rate tends toward zero (i.e.:
in the static zone) raises obvious numerical problems. The
basic idea to overcome this issue consists in regularizing
the viscosity by adding a small quantity (A) to the denom
inator of the particulate viscosity ,. ' 1'/( +A)
then the divergence is controlled by this parameter and the
viscosity is kept finite. In other words, the static zone in
the frictional rheology is replaced by a very viscous zone.
In the next section of this paper we will discuss the influ
ence of the value of the parameter A on the model solution.
In our implementation, we use piecewise quadratic poly
nomial approximation for the velocity and piecewise lin
ear discontinuous approximation for the pressure. In the
computations, we have employed a 27nodes hexahedra
element (H27) for the velocities. The incompressibility
constraint is solved by a penalisation method. The code
is developed with the PETSc library Balay et al. (2001,
2004, 1997) which provides several parallel iterative and
direct solvers. As we use a penalisation method to cope
with the incompressibility constraint, all the algebraic sys
tems have been solved by the MUMPS direct solver Amestoy
et al. (2000, 2001, 2006) with a penalty parameter set to
109 for all the simulations presented in this paper.
Results
To close these equations we need to prescribe the effec
tive fluid viscosity rie = (1 + 5/2y) (Einstein 1906), the
particulate viscosity = /i 1'/ I Ij  (Jop et al. 2006)
We present the results of the previous twophase numeri
cal model applied to the flow of a Newtonian fluid over a
granular bed. First we compare the results of the numer
ical model with the analytical solution of Ouriemi et al.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
z have performed simulations on the half of the domain in
u= v= = o the transverse direction for obvious symmetry reasons.
H Figure 4 shows the velocity profile in a cross section
,f = o 0 = wf = o of the ducts. The contour colors represent the xvelocity
OfS a s of the mixture (um). The horizontal thick solid line at
ax d P ax z 0.5 represents the position of the granular bed. The
fluid and the mixture are sheared in both z and y direc
Stions inducing an increase in the friction compared with
vip = o vf wP= 0 the twodimensional case. Due to this shear increase the
oup Hp OuP velocity is lower than in the twodimensional case. Figure
ax ax 5 shows the velocity profiles of the fluid phase velocity in
__ blue and the particulate phase velocity in red (an offset of
ad af 103 has been added to make the particulate phase veloc
v1 = w = v = uv = wp = 0 and =0
az ity visible) obtained with the twofluid model. These re
Figure 1: Sketch of the flow of a Newtonian fluid over a sults illustrate the good behaviour of the numerical model
granular bed. for threedimensional flow configurations.
12i i 1.11 for the bedload transport in laminar shearing flows
to validate quantitatively the two formulations of the two
phase flow model. We have also looked at the computa Analytic
tional efficiency of the numerical model associated with UP
both formulations.
The sketch of the problem and boundary conditions are 0.8
given in figure 1. The lower half of the domain is filled
with particles at = 0.55 immersed in a fluid and the 0.6
upper part is filled with pure fluid ( 0). Therefore N
in this problem the values of the dimensionless numbers 0.4
are: Re 2 102, Ga 11, R = 0.4 and d/H
30. The regularisation parameter is set to A 10 6 s 1 0.2
when its value is not mentioned in the figure captions.
We solved by FEM the two formulations of the twophase 0 0.001 0.002 0.003 0.004 0.005
flow model for a 6x1x40 mesh with a requested absolute U
residual lower than 1011 per degree of freedom. (a) Twofluid model
Figures 2a and 2b present the comparison of the hori
zontal velocity profiles obtained for the two formulations 1 Analytic
of the twophase flow model (twofluid and mixedfluid) UM
by numerical simulations compared with the analytical
solution proposed by Ouriemi et al. (2009a). In both fig 0
ures the black solid line represents the analytical solution.
The good agreement between the numerical solution and 0.6
the analytical one gives a first qualitaitve validation of the N
FEM model for the bedload transport. 0.4
The spatial convergence analysis for both formulations
(see figure 3) shows that the regularisation of the granular 0.2
rheology reduces the spatial order of convergence to order
one whereas a second order is theoretically hopped. Also, O ' 0 0 0
0.001 0.002 0.003 0.004 0.005
to reach the same accuracy with both formulations we U
observe that the regularisation for the mixedfluid model (b) Mixedfluid model
must be decreased by one order of magnitude compared
ith te d easo d bone Figure 2: Longitudinal velocity profiles for the flow of
with the twofluid one.
wie ow l model t a Newtonian fluid over a granular bed between two in
We now apply the model to threedimensional configu
finite parallel planes obtained by numerical simulations
rations: a square and a circular crosssection ducts, with
the same values of the dimensionless numbers Re 2 102 for a) the twofluid model and b) the mixedfluid model
Ga = 11, R 0.4, d/H 30 and Bn 2 104. We 'compared with the analytical solution of Ouriemi et al.
(2009a).
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Error e = f(h)
le08
le09
x' O
0
S ..... e_p (le6) 
S(.. e_f (le6) ..... 
e_m(le6) El
e_m(le7) 0
0.01 0.1
Figure 3: RMS Error against analytical solution for the flow of a Newtonian fluid over a moving granular bed between
two infinite parallel planes: ep and ef stands for the particulate phase and the fluid phase error respectively for the
twofluid model whereas em designates the mixture velocity error for the mixedfluid model. The value in brackets is
the value of the regularisation parameter A. The RMS error is defined as: e 7= ( (U U ana)2)1/2 where N is
the number of nodes in the mesh.
025
(a) Square duct
U
S00040
S 00035
0 0030
0 0025
00020
0 0015
S 00010
0 0005
0 0000

0 00025 0005 0 5
(b) Cylindrical duct
Figure 4: Velocity profile obtained by numerical simulations with the mixedfluid model for a) the square crosssection
duct (6x20x40) and b) the circular crosssection duct (6x896).
(a) Square duct
025 
N 0
0 25 F
W F   
0 0 002 0 004 5
u u
(b) Cylindrical duct
Figure 5: Velocity profile obtained by numerical simulations with the twofluid model for a) the square crosssection
duct (6x20x40) and b) the circular crosssection duct (6x896). The fluid phase velocity is in blue and the particulate
phase velocity is in red. An offset of 10 3 has been added to the velocity of the particulate phase (up) to make it
visible.
025
Figure 6 shows the application of the mixedfluid model
coupled with an Arc Length Continuation algorithm that
allows us to explore the parameter space in terms of the
longitudinal pressure gradient dp/Ox. The agreement be
tween the analytical solution of Ouriemi et al. 12'1 1'.1 and
the numerical model for the height of the flowing granular
layer Hp He and the particle flux Qp for a wide range
of longitudinal pressure gradient gives another validation
of our three dimensional numerical model. As for the ve
locity profiles we have explored the parameter space for
three dimensional configurations that allows us to obtain
the height evolution of the flowing granular layer Hp He
and the particle flux Qp in rectangular and circular cross
section ducts. These results are presented in Figure 7 and
illustrate the nontrivial effect of the geometry on the be
haviour of the flowing granular layer. This method ap
plied to complex geometry is a powerful tool to compare
the prediction of the initial twophase model with experi
ments in real geometry such as rectangular or cylindrical
ducts. This is a strong argument for the development of a
three dimensional numerical model.
0.05 0.1
dpdx
0.15 0.2
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Conclusions
In summary, we have developed a three dimensional nu
merical model to simulate incompressible particulate two
phase flows for the flow of a Newtonian fluid over a granu
lar bed. The incompresibility is imposed by a penalisation
method and the viscoplastic characteristic of the granular
rheology is dealt with a regularisation technique. We have
studied the accuracy of our numerical solution by compar
ison with the analytical solution derived by Ouriemi et al.
'2'i "1.,, in function of the spatial discretisation and the
value of the regularisation parameter for the mixedfluid
and the twofluid model. We have concluded that one
must take a regularisation parameter one order of mag
nitude lower for the mixedfluid model than for the two
fluid one to reach the same accuracy.
The Finite Element method allows us to deal with ar
bitrary geometries and we have shown numerical results
0.1 0.2 0.3
dpdx
(a) Height of the flowing layer
(a) Height of the flowing layer
3e04
2e04
le04
Oe+00
0 0.05
dpdx
(b) Particle flux
6e04
4e04
2e04
Oe+00
0.1 0.15
Figure 6: Comparison of the numerical results for the
2D (mesh size 4x1x640) configurations with the analyti
cal solution of Ouriemi et al. 12'1, '.ii in terms of a) the
height of the flowing granular layer Hp He and b) the
particle flux Q, versus the longitudinal pressure gradient.
2D
Rectangular .....
Cylindrical
0 0.1 0.2
dpdx
(b) Particle flux
0.3
Figure 7: Comparison of the numerical results for the 2D
and 3D configurations (rectangular and circular cross sec
tion ducts) in terms of a) the height of the flowing granular
layer Hp He and b) the particle flux Qp versus the lon
gitudinal pressure gradient. The mesh sizes for the 2D
case is 4x1x640, for the rectangular duct (aspect ratio of
W/H 0.2692) is 1x160x320 and for the cylindrical
duct is 1x179998.
 I .............
of threedimensional simulation of the bedload transport
in a square and a circular crosssection ducts. This model
had also been coupled with an Arc Length Continuation
algorithm to describe the evolution of the height of the
flowing granular layer and the particle flux with the lon
gitudinal pressure gradient in both 2D and 3D configura
tions. This method is particularly interesting in the per
spective of comparing the prediction of this model with
experiments that are always carried out in 3D configura
tions.
Ouriemi et al. (2009b, 2010) have performed a sim
ple linear stability analysis which provides realistic pre
dictions for the formation of small dunes. The aim of
the present numerical model is to perform a full three
dimensional stability analysis that accounts for the two
phase nature of the problem.
Acknowledgement
Funding from the Institut Franqais du P6trole and Agence
National de la Recherche (Project Dunes ANR073_18
3892) are gratefully acknowledged.
References
P. R. Amestoy, I. S. Duff, and J.Y. L'Excellent. Mul
tifrontal parallel distributed symmetricand unsymmetric
solvers. Comput. Methods Appl. Mech. Eng., 184:501
520, 2000.
P. R. Amestoy, I. S. Duff, J. Koster, and J.Y. L'Excellent.
A fully asynchronous multifrontal solver using distributed
dynamic scheduling. SIAM Journal on Matrix Analysis
and Applications, 23(1): 1541, 2001.
P. R. Amestoy, A. Guermouche, J.Y. L'Excellent, and
S. Pralet. Hybrid scheduling for the parallel solution
of linear systems. Parallel Computing, 32(2):136156,
2006.
R. A. Bagnold. The flow of cohesionless grains in fluids.
Phil. Trans. R. Soc. Lond., 249:235297, 1956.
Satish Balay, William D. Gropp, Lois Curfman McInnes,
and Barry F. Smith. Efficient management of parallelism
in object oriented numerical software libraries. In E. Arge,
A. M. Bruaset, and H. P. Langtangen, editors, Modern
Software Tools in Scientific Computing, pages 163202.
Birkhiuser Press, 1997.
Satish Balay, Kris Buschelman, William D. Gropp,
Dinesh Kaushik, Matthew G. Knepley, Lois Curfman
McInnes, Barry F. Smith, and Hong Zhang. PETSc Web
page, 2001. http://www.mcs.anl.gov/petsc.
Satish Balay, Kris Buschelman, Victor Eijkhout,
William D. Gropp, Dinesh Kaushik, Matthew G. Knep
ley, Lois Curfman McInnes, Barry F. Smith, and Hong
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Zhang. PETSc users manual. Technical Report ANL
95/11 Revision 2.1.5, Argonne National Laboratory,
2004.
Julien ( .lii hi.ii .1 d Marc M6dale. A 3D numerical model
for incompressible twophase flow of a granular bed sub
mitted to a laminar shearing flow. Computer Methods in
Applied Mechanics and Engineering (in press), 2009.
A. Einstein. Eine neue bestimmung der molekuldimen
sionen. An. Phys., 19:289 306, 1906.
H.A. Einstein. Formulas for the transportation of bed load.
Transactions of the American Society of Civil Engineers,
2140:561597, 1942.
Y. Forterre and 0. Pouliquen. Flows of dense granular me
dia. Annual Review of Fluid Mechanics, 40:124, 2008.
doi: 10.1146/annurev.fluid.40.111406.102142. URL
http://link.aip.org/link/?PHF/17/103301
J. Happel and H. Brenner. Low Reynolds number hydro
dynamics. Martinus Nijhof, The Hague, 1973.
R. Jackson. Locally averaged equations of motion for
a mixture of identical spherical particles and a newto
nian fluid. Chemical Engineering Science, 52:24572469,
1997.
R. Jackson. The dynamics of fluidized particles. Cam
bridge University Press, Cambridge, 2000.
Pierre Jop, Yoel Forterre, and Olivier Pouliquen. A
constitutive law for dense granular flows. Nature, 441:
727 730, 2006. doi: 10.1038/nature04801. URL
http://dx.doi.org/0101038/nature04801.
E. MeyerPeter and R. Muller. Formulas for bedload
transport. In 2nd Meeting of the International Associa
tion of Hydraulic and Structural Research, pages 3464,
1948.
M. Ouriemi, J. Chauchat, P. Aussillous, M. M6dale, and
E. Guazzelli. Sediment transport and dunes in pipe
flow. In 7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, 2010.
Malika Ouriemi, Pascale Aussillous, and Elisabeth
Guazzelli. Sediment dynamics. Part I: Bedload transport
by shearing flows. Journal of Fluid Mechanics, 636:295
319, 2009a.
Malika Ouriemi, Pascale Aussillous, and Elisabeth
Guazzelli. Sediment dynamics. Part II: Dune formation
in pipe flow. Journal of Fluid Mechanics, 636:321336,
2009b.
S Yalin. An expression for bedload transportation. J.
Hydraul. Division, HY3:221250, 1963.
