Available online at www.sciencedirect.

com
Journal
of
Terramechanics
Journal of Terramechanics 44 (2007) 339–346
www.elsevier.com/locate/jterra

Numerical simulation of soil–water interaction using smoothed

particle hydrodynamics (SPH) method
a,*
Ha H. Bui , K. Sako b, R. Fukagawa a,b

a
Department of Civil Engineering, Ritsumeikan University, BKC Campus, Kusatsu 525-8577, Japan
b
Center of Excellence (COE) Promotion Program, Ritsumeikan University, BKC Campus, Kusatsu 525-8577, Japan

Accepted 18 September 2007

Available online 5 November 2007

Abstract

An application of smoothed particle hydrodynamics (SPH) to simulation of soil–water interaction is presented. In this calculation,
water is modeled as a viscous fluid with week compressibility and soil is modeled as an elastic–perfectly plastic material. The Mohr–Cou-
lomb failure criterion is applied to describe the stress states of soil in the plastic flow regime. Dry soil is modeled by one-phase flow while
saturated soil is modeled by separate water and soil phases. Interaction between soil and water is taken into account by means of pore
water pressure and seepage force. Simulation tests of soil excavation by a water jet are calculated as a challenging example to verify the
broad applicability of the SPH method. The excavations are carried out in two different soil models, one is dry soil and the other is fully
saturated soil. Numerical results obtained in this paper have shown that the gross discontinuities of soil failure can be simulated without
any difficulties. This supports the feasibility and attractiveness of this a new approach in geomechanics applications. Advantages of the
method are robustness, conceptual simplicity and relative ease of incorporating new physics.
 2007 ISTVS. Published by Elsevier Ltd. All rights reserved.

Keywords: Smoothed particle hydrodynamics (SPH); Soil–water interaction; Dry soil; Saturated soil; Seepage force

1. Introduction cially it has the potential to handle problems with extre-

mely large deformation whereas the grid-based techniques
Smoothed particle hydrodynamics (SPH) was first devel- require additional computational effort.
oped to simulate astrophysical fluid dynamics [1,2] then it In this paper, such advantages of SPH will be exploited
has been successfully applied to a vast range of problems. for the first time to simulate the soil–water interaction in a
These include elastic–plastic flow [3,4], quasi-incompress- soil. Simulation of soil excavation by a water jet will be per-
ible flow [5,6,7], and fracture of brittle solids [8]), etc. formed as an example in order to test our models for dry
SPH is a fully Lagrangian and mesh free method. A com- soil and for saturated soil. The dry soil is modeled by a
putational domain is represented by computational points layer of a single type of particles generated in a soil
(particles), typically of fixed mass. Each particle carries domain; the soil is assumed to be an elastic–perfectly plas-
all field variable information like density, pressure, veloc- tic material and the stress states of soil in the plastic flow
ity, etc., and moves with the material velocity. This method regime are followed by the Mohr–Coulomb failure crite-
has some advantages over grid-based techniques because rion. The saturated soil is modeled by mixing both soil par-
its concept is simple and it is relatively easy to incorporate ticles and water particles together in the same soil domain.
complicated physical effects into the SPH formalism. Espe- The effect of pore water pressure on soil particles is calcu-
lated by allowing the pressure of water particles to contrib-
ute to the pressure of soil particles through the SPH
*
Corresponding author. Tel.: +81 77 566 1111x8717. summation. Furthermore, the seepage force acting on the
E-mail address: gr022026@se.ritsumei.ac.jp (H.H. Bui). soil particles due to the effect of water flow is introduced

0022-4898/$20.00  2007 ISTVS. Published by Elsevier Ltd. All rights reserved.

doi:10.1016/j.jterra.2007.10.003
340 H.H. Bui et al. / Journal of Terramechanics 44 (2007) 339–346

into the momentum equation as an external force that, in There are many possible choices of kernel functions
this study, is computed based on Darcy’s law [9]. Numeri- which satisfy (2)–(4). However, the most popular one is
cal results obtained in this study have shown that SPH the cubic spline interpolation function, which was pro-
could be a valuable method for simulation of complex posed by Monaghan and Lattanzio [10]. If we define R as
problem in soil mechanics. Advantages of the method are the relative distance between points x and x 0 ,
robustness, conceptual simplicity and relative ease of incor- R = jx  x 0 j/h, then the cubic spline function has the fol-
porating new physics. lowing form:
8 2 3
>
< 1:5  R þ 0:5R 0 6 R < 1
2. Simulation method using SPH
W ðR; hÞ ¼ ad  ð2  RÞ3 =6 16R<2 ð5Þ
>
:
2.1. Basic ideas of SPH 0 RP2
In one-, two-, and three-dimensional space, ad = 1/h, 15/
The foundation of SPH is interpolation theory. The con-
(7ph2), 3/(2ph3), respectively (see Fig. 1b).
servation laws of continuum fluid dynamics, in the form of
The ‘‘kernel estimation’’ for the spatial derivative of a
partial differential equation, are transformed into integral
vector quantity $ Æ f(x) is obtained simply by substituting
equations through the use of an interpolation function that
f(x) with $ Æ f(x) into Eq. (1), which gives
gives the ‘‘kernel estimate’’ of the field variable at a point Z
[4]. The term ‘‘kernel estimate’’ refers to a weight function hr  f ðxÞi ¼ ½r  f ðx0 ÞW ðx  x0 ; hÞdx0
and defines how much of each field variable contributes to
ZX
the field variable at a point. For a given function f(x), the
¼ r  ½f ðx0 ÞW ðx  x0 ; hÞdx0
‘‘kernel estimation’’ of this function, denote by Æf(x)æ, is X
defined as Z
Z  f ðx0 Þ  rW ðx  x0 ; hÞdx0 ð6Þ
hf ðxÞi ¼ f ðx0 ÞW ðx  x0 ; hÞdx0 ð1Þ X
X
Using the divergence theorem, the first integral on the
where X is the volume of the integration that contains x and right-hand side of Eq. (6) can be converted into an integral
x 0 ; and h is the smoothing length defining the influence do- over the surface S of the domain inside of integration X,
main of the ‘‘kernel estimate’’. W(x  x 0 , h) is the kernel, or Z
smoothing function, which must satisfy the following three hr  f ðxÞi ¼ f ðx0 ÞW ðx  x0 ; hÞ  ~
ndS
properties: the first one is the normalization condition, S
Z
Z
 f ðx0 Þ  rW ðx  x0 ; hÞdx0 ð7Þ
W ðx  x0 ; hÞdx0 ¼ 1 ð2Þ X
X
the second condition is the Delta function property that is where ~ n is the unit vector normal to the surface S. Since the
observed when smoothing length approaches zero, smoothing function W is defined to have compact support,
the surface integral on the right-hand side of Eq. (7) is zero.
lim W ðx  x0 ; hÞ ¼ dðx  x0 Þ ð3Þ
h!0 Therefore, the ‘‘kernel estimation’’ of the spatial derivative
and the third condition is the compact condition, $ Æ f(x) now becomes,
Z
W ðx  x0 ; hÞ ¼ 0 when jx  x0 j > jh ð4Þ
hr  f ðxÞi ¼  f ðx0 Þ  rW ðx  x0 ; hÞdx0 ð8Þ
X
with j is a constant related to the smoothing function for
point at x and defines the influence domain of the ‘‘kernel In the SPH method, the calculation domain is represented
estimate’’ (see Fig. 1a). by a finite number of points (particles), which carry mass

a Smoothing function b
S
Ω
Function × αd

rij κh
j
Smoothed function
Derivative of the smoothed function

(x-x’)/h

Fig. 1. (a) Particle approximations using particles within the influence domain of the smoothing function W for particle i. The influence domain is circular
with a radius of jh. (b) The cubic spline kernel and its first derivative.
 H.H. Bui et al. / Journal of Terramechanics 44 (2007) 339–346 341

and field variable information such as density, stress, en- is the external force. The Einstein convention for summa-
ergy, etc. Accordingly, the continuous integral representa- tion over repeated indices is applied.
tion for f(x) is approximated in the following form: The stress tensor normally consists of two parts: an iso-
Z tropic pressure p and a viscous shear stress s,
hf ðxÞi ¼ f ðx0 ÞW ðx  x0 ; hÞdx0
X rab ¼ pdab þ sab ð14Þ
X
N
mj For a Newtonian fluid like water, the viscous shear stress is
 f ðxj ÞW ðx  xj ; hÞ ð9Þ
j¼1
qj proportional to the shear strain rate denoted by e through
the viscosity l,
Using Eq. (9), the particle approximation for a function at
particle i can be written as sab ¼ leab ð15Þ

XN where
mj  
hf ðxi Þi ¼ f ðxj ÞW ij ð10Þ
qj ovb ova 2 ovc ab
j¼1 eab ¼ þ  d ð16Þ
oxa oxb 3 oxc
where
Using the SPH approximation discussed in Section 2.1, the
W ij ¼ W ðxi  xj ; hÞ system of partial differential equations (12), (13) and (16)
Eq. (10) states that the value of a function at particle i is can be converted into the SPH formulations which will
approximated using a weighted average of those values of be used to solve the motion of fluid particles. First, the
the function at all other particles in the influence domain right-hand side of the continuity Eq. (12) can be rewritten
of particle i. Following the same argument, the particle as
approximation for the spatial derivative of the function  
ova o a a oq
at particle i is q a ¼  ðqv Þ  v  ð17Þ
ox oxa oxa
XN
mj
hr  f ðxi Þi ¼ f ðxj Þ  ri W ij ð11Þ Substituting Eq. (17) into Eq. (12) and applying the SPH
j¼1
qj particle approximation to the gradients, the continuity
density equation of particle i in the SPH formulation will
where be,
!
xi  xj oW ij xij oW ij Dqi XN XN
ri W ij ¼ ¼ a oW ij a oW ij
rij orij rij orij ¼ m j vj  vi mj a ð18Þ
Dt j¼1
oxai j¼1
oxi
It should be noted here that the negative sign in Eq. (11) or
was already removed because $iWij is taken with respect
to particle i, not to particle j. The more details of the Dqi X N
oW ij
¼ mj ðvai  vaj Þ a ð19Þ
SPH method are given in [11,12]. Dt j¼1
oxi

The derivation of the SPH formulation for the particle

2.2. SPH model for water approximation of momentum evolution is somewhat simi-
lar. Considering the first term on the right-hand side of Eq.
For a fluid like water, with sound speed  103 ms1, the (13),
Mach number is extremely small, and it is customary to  
model it as exactly incompressible. However, the approach 1 orab o rab rab oq
b
¼ b
þ 2 ð20Þ
in SPH is different; the real fluid is approximated by an q ox ox q q oxb
artificial fluid which is more compressible than the real and applying the SPH particle approximation to the gradi-
one, while still possessing a speed of sound which is much ents leads to
larger than the flow speed, and which therefore has very
small density fluctuation. The governing equations for fluid   !
1 orab X
N
rab
i
ab
rj oW ij
flow are the well-known Navier–Stokes equations, which in ¼ mj 2
þ 2 ð21Þ
the Lagrangian description state the conservation of mass q oxb i j¼1
qi qj oxbi
and momentum as follows:
Substituting Eq. (21) into Eq. (13) then the SPH form of
Dq ova momentum equation at particle i is
¼ q a ð12Þ
Dt ox !
Dva 1 orab Dvai XN
rab
i rab
j oW ij
¼ þ fa ð13Þ ¼ mj 2
þ 2 b
þ fia ð22Þ
Dt q oxb Dt j¼1
q i q j ox i

where a, b are the Greek superscripts used to denote the Similarly, the SPH approximation of shear strain rate eab
coordinate directions, rab is the total stress tensor and fa for particle i is
342 H.H. Bui et al. / Journal of Terramechanics 44 (2007) 339–346

XN
mj b oW ij X N
mj a oW ij K for calculation sometime gives large fluctuation in pres-
eab
i ¼ ðvj  vbi Þ a þ ðvj  vai Þ b sure especially for small scale simulation. Therefore, K
j¼1
qj oxi j¼1
qj oxi
! should be chosen carefully in order to insure the near
2X N
mj c oW ij incompressibility of soil while minimizing pressure fluctua-
 ðv  vcj Þ c dab ð23Þ tions. This paper chooses K to be ten times the maximum
3 j¼1 qj j oxi
initial pressure (K = 10Pomax).
The final equation needed to solve the above Navier– The deviatoric shear stress s appearing in Eq. (26) can be
Stokes equations for water is the ‘‘equation of state’’ which estimated using the following equation:
is used to estimate the pressure change of water. Monaghan  
1
[5] applied the following form when he first simulated the s_ ab ¼ 2Ge_ ab ¼ 2G e_ ab  dab e_ cc ð28Þ
3
free surface flows of water using SPH,
"  # where G is the shear modulus, s_ is the stress rate tensor and
k
q e_ is the strain rate tensor. However, for finite rotation this
p¼B 1 ð24Þ
qo equation does not satisfy material frame indifference, or
‘‘MFI’’, i.e. the material response will depend in an
where k is a constant, and set equal to seven in most cir-
unphysical way on rotation of the material. In order to sat-
cumstances; qo is the reference density; B is a problem
isfy MFI, the Jaumann rate is adopted here with the fol-
dependent parameter, which sets a limit for the maximum
lowing constitutive equation:
change of the density. If we choose B such that the density
fluctuation
pffiffiffiffiffiffiffiffiffiffiffiffiof
ffi water is 1%ffi and the sound speed of water is
pffiffiffiffiffiffiffiffiffiffiffiffi s_ ab  sac -bc  scb -ak ¼ 2Ge_ ab ð29Þ
c ¼ op=oq  kB=qo then B will be calculated as
100V 2typ qo where the strain rate tensor e_ and the rotation rate tensor -
B¼ ð25Þ are defined as
k  
where Vtyp is the typical speed of water. It should be noted 1 ova ovb
e_ ab ¼ þ ð30Þ
that Eq. (25) is different from that quoted in [5], 2 oxb oxa
 
B ¼ 100V 2typ =ðqo kÞ; we suspect a typographical error. 1 ova ovb
-ab ¼  ð31Þ
2 oxb oxa
2.3. SPH model for soil
Using the concept of SPH approximation, Eqs. (30) and
(31) can be converted into the approximate SPH form,
Modeling the behavior of soil using the SPH method is
yielding for particle i,
similar to that of water. The SPH form of conservation
equations, Eqs. (19) and (22), are still used to estimate 1X N
mj a oW ij
the density and motion of soil particles. The key difference e_ ab
i ¼ ðvj  vai Þ b
2 j¼1 qj oxi
between these two models is the calculation of the stress !
tensor appearing in Eq. (22) in which the pressure and mj oW ij
stress–strain relationship of soil are calculated differently þ ðvbj  vbi Þ a ð32Þ
qj oxi
from those of water; soil is assumed herein to be an elas-
tic–plastic material. The stress tensor of soil is made up 1X N
mj a oW ij
-ab
i ¼ ðvj  vai Þ b
of two parts: isotropic pressure p and deviatoric shear 2 j¼1 qj oxi
stress s, !
mj oW ij
rab ¼ pdab þ sab ð26Þ  ðvbj  vbi Þ a ð33Þ
qj oxi
The pressure in SPH is normally calculated using an ‘‘equa-
tion of state’’, which has a function of density change. Therefore, combining Eqs. (28), (29), (32) and (33) the
Since soil is assumed to have elastic behavior, the pressure deviatoric shear stress components can be calculated. These
equation of soil will obey Hooke’s law, shear stress components should be restricted to the failure
  surface when plastic flow takes place. In this paper, the
DV q plastic flow regime of soil is determined by the Mohr–Cou-
p ¼ K ¼K 1 ð27Þ
V qo lomb failure criterion and, according to an elastic–perfectly
plastic soil model, the deviatoric shear stress components in
where K is the bulk modulus; DV/V is the volumetric this region are scaled back to the maximum shear stress de-
strain; and qo is the initial density. Eq. (27) represents the fined by
change in pressure resulting from the change in density,
sf ¼ c þ p tan / ð34Þ
i.e. equation of state for soil. If the value of K is large, small
variation in density will result in extremely large pressure where c is the cohesion; p is the pressure; and / is the angle
fluctuation, and this requires an extremely small time step of internal friction. The two soil parameters c and / can be
for calculation. Taking the actual value of bulk modulus experimentally determined through direct shear test or
 H.H. Bui et al. / Journal of Terramechanics 44 (2007) 339–346 343

triaxial test. A more accurate constitutive model for soil, Water particle Soil particle
not yet implemented in this research, would take into ac-
count the plastic strain in the plastic flow regime.
Fully saturated
2.4. Soil–water interaction modeling soil domain

Modeling of soil and water in the framework of

smoothed particle hydrodynamics has been presented in Saturated soil is divided
the previous sections however, these models only allow to Water layer into two layers Soil layer
simulate one-phase flow. In order to simulate the soil–
water interaction, which consists of two-phase flows, it is
necessary to derive an equation to describe the interaction Seepage force
between soil and water phases, as follows. When ground-
water is seeping through the pores of a soil, viscous friction
will produce drag on soil particles in the direction of water Pore water pressure
flow, so-called seepage force. This seepage force acts on the
soil particles in addition to the gravitational force, and will
be introduced into the momentum equations for soil and Fig. 2. Modeling of soil–water interaction in SPH simulation.
water as an external force according to the following model
equation based on the Darcy’s law [9], Momentum equation for water phase
XN  
ðvwater  vsoil Þ Dvaa pa pb oW ab
fseepage ¼ cw n ð35Þ ¼ mb 2 þ 2
k Dt b¼1
qa qb oxaa
!
XN
sab sab oW ab XN
where cw is the unit weight of water; n is the porosity; and k þ mb a
þ b
þ f a
 mi
a
is the soil permeability. This paper choose k = 0.05 cm/s b¼1
q2a q2b oxaa i¼1
for numerical test of saturated soil.
fiaseepage
As saturated soil consists of soil and water mixed  W ia ð37Þ
together while standard SPH models only handle one- qi qa
phase problem, it is necessary to develop a saturated soil where the subscripts i and j represent for soil particles while
model using in SPH simulation. This saturated soil model a and b are used for water particles. It is clear that these
will be described as follows. Assume that the saturated soil momentum equations of saturated soil are different from
domain in SPH can be divided into two separate phases, that of single phase, Eq. (22). The presence of the seepage
which are water phase and soil phase as shown in Fig. 2. force in Eqs. (36), (37), and the contribution of pressure
The motion of SPH particles on each phase is solved sepa- from water to soil in Eq. (36) make them possible to simu-
rately using its own SPH governing equations, which are late the effect of seepage force and pore water pressure in
SPH for soil and SPH for water. These two-phases are then the saturated soil model, as a result the interaction between
superimposed and the interaction between two-phases will soil and water could be simulated through SPH.
be taken into account through the seepage force, which is
introduced into the momentum equation as mentioned 2.5. Artificial viscosity
above. In addition, the water pressure is also allowed to
contribute to the soil pressure during the overlapping pro- In SPH simulation, in order to handle shock waves and
cedure. This allows us to simulate the pore water pressure, to prevent the unphysical penetration of particles, an arti-
which always exists in natural saturated soil. Accordingly, ficial viscosity has been introduced to the pressure term
the momentum equations for saturated soil can be summa- in the momentum equation. Of the many types of artificial
rized as follows: viscosity developed so far, the most widely used is specified
Momentum equation for soil phase as follows [12]:
! ( ac / þb/2
ij ij ij
Dvai XN
p pj oW ij vij  xij < 0
¼ mj 2i þ 2 Pij ¼ q
ij ð38Þ
Dt j¼1
qi qj oxai 0 vij  xij P 0
! where
X
N
sab
i sab
j oW ij
þ mj þ 2 þ fia hij vij  xij 1
j¼1
qi 2
qj oxai /ij ¼ ; cij ¼ ðci þ cj Þ ð39Þ
2
jxij j þ eh2ij 2
X
N
pa oW ia X
N
fiaseepage
 ma
qi qa oxai
þ ma
qi qa
W ia ð36Þ hij ¼ 1 ðhi þ hj Þ; 1
ij ¼ ðqi þ qj Þ
q ð40Þ
a¼1 a¼1 2 2
344 H.H. Bui et al. / Journal of Terramechanics 44 (2007) 339–346

In the above equations, c and v represent the speed of the other hand, if it is too small, some particles will pene-
sound and particle velocity vector; x is the position vector trate the boundary.
of particle; and e is usually taken as 0.01 ([5]). The param-
eter a produces a bulk viscosity while b suppresses particle 2.7. Time integration
interpenetration. Monaghan chose a = 0.01, b = 0 to simu-
late the free surface flow, Libersky [4] chose a = 2.5, Eqs. (19), (29), (41) and (42) are integrated using a stan-
b = 2.5 for material strength simulation. However, neither dard Leap–Frog (LF) algorithm with the following integra-
is suitable for our simulations since we simulate the water tion schemes:
jet flow with high velocity and the soil behavior is different
snþ1 n
ab ¼ sab þ Dt  Dsi ðtÞ ð44Þ
from that of solid. This paper we choose a = 0.01, b = 1 for
nþ1=2 n1=2
water and a = 1, b = 1 for soil. The momentum equations qi ¼ qi þ Dt  Dqi ðtÞ ð45Þ
for saturated soil after introducing the artificial viscosity to nþ1=2 n1=2
vi ¼ vi þ Dt  Dvi ðtÞ ð46Þ
the pressure term are:
nþ1=2
Soil phase xnþ1
i ¼ xni þ Dt  vi ð47Þ
!
Dvai X N
pi pj oW ij where Dt is the time step; n indicates the current time t; and
¼ mj 2 þ 2 þ Pij dij (n + 1) indicates the advanced time (t + D t). The stability
Dt j¼1
qi qj oxai
! of the LF scheme is guaranteed by several criteria for time
X N
sab sab
j oW ij step. The first one is the Courant–Friedrichs–Levy (CFL)
i
þ mj 2 þ 2 a
þ fia condition,
j¼1
q i q j oxi
Dt 6 0:25 minðhi =cÞ ð48Þ
X
N
p oW ia X N
fiaseepage
 ma a þ m a W ia ð41Þ and additional constraints due to the magnitude of particle
qi qa oxai qi qa
a¼1 a¼1 acceleration fi ([5])
Water phase pffiffiffiffiffiffiffiffiffiffi
Dt 6 0:25 minð hi =fi Þ ð49Þ
XN  
Dvaa p p oW ab and viscous diffusion ([7])
¼ mb a2 þ b2 þ Pab dab
Dt q q oxaa
b¼1 a b
! Dt 6 0:125ðh2i =tÞ ð50Þ
XN
sab sab oW ab
þ mb a2 þ b2 a
þ faa Accordingly, the time step Dt for a simulation should be
b¼1
q a q b ox a chosen to be the minimum value of Eqs. (48)–(50).
X
N
fiaseepage
 mi W ia ð42Þ 3. Calculations and discussion
i¼1
qi qa
In order to evaluate our SPH model for soil–water inter-
2.6. Boundary condition action, simulations of excavation by a water jet have been
carried out for dry soil and saturated soil. The dry soil is
To model the solid boundary condition, we simply fol- modeled by one type of particles with uniform material
low the method proposed by Monaghan [5]. This method properties filled in a rectangular reservoir 50 cm in length
assumed the presence of a ghost line of particles located
right on the solid boundary to produce a highly repulsive
force to the real particles near the boundary, and thus to
prevent these particles from penetrating the solid bound-
ary. The repulsive force is calculated similarly to the Len- Nozzle vjet = 25m/s
nard-Jones molecular force as,
8 h n1  n2 i   2cm
Solid boundary
> xij Water jet
< D rrijo  rrijo 2
rij
ro
rij
61
F ¼   ð43Þ
>
:0 ro
rij
> 1

where parameters n1 and n2 are usually taken as 12 and 6, 20cm Soil stratum
respectively; rij is the distance between ghost particle and
real particle near the boundary. D is a problem parameter
chosen to be of the same order as the square of the largest
velocity. The cutoff distance ro must be selected close to the
50cm
initial particle spacing; if ro is too large, a high repulsive
force will push real particles from the solid boundary; on Fig. 3. Initial arrangement of simulation test.
 H.H. Bui et al. / Journal of Terramechanics 44 (2007) 339–346 345

and 20 cm in width. These particles have the following smoothing length of 0.5 cm. The water jet is generated
properties: density q = 2.7 g/cm3, Young’s modulus above the water phase using the same type of water parti-
E = 150 MPa and Poisson’s ratio v = 0.3. The saturated cles with the same smoothing length as above. Excavation
soil model is generated by mixing two types of particles is simulated by applying the jet flow onto the soil domain
that model water and soil phases into the same soil domain. at a speed of 25 m/s. The initial arrangement of the simu-
The soil particles in the saturated soil have the same mate- lation is depicted in Fig. 3.
rial properties with the particles in the above dry soil, and Fig. 4 shows the erosion process of soil excavation by
the water particles have density of q = 1.0 g/cm3 and vis- water jet on the dry soil model at representative times. In
cosity of l = 103 Ns/m2. There are a total of 4000 parti- this simulation, both water pressure and seepage force,
cles generated regularly on each phase with an initial which results from the relative velocity between water par-

Fig. 4. Erosion during excavation by a water jet on dry soil model by SPH simulation. Water jet has speed of vjet = 25 m/s.

Fig. 5. Erosion during excavation by a water jet on saturated soil model by SPH simulation. Water jet has speed of vjet = 25 m/s.
346 H.H. Bui et al. / Journal of Terramechanics 44 (2007) 339–346

ticles in the jet and soil particles, are taken into account. By cous fluid with week compressibility while soil was modeled
this way, transmission of energy from the water jet to the as an elastic–perfectly plastic material. Interaction between
soil particles could be simulated and the erosion process soil and water was modeled by means of pore water pres-
could be reproduced. The soil was easily eroded by the sure and seepage force. Simulation tests of soil excavation
jet flow and soil particles were splashed together with the by a water jet on dry soil and saturated soil have been pre-
water particles. Furthermore, the erosion was also ampli- sented. The results have shown that the extremely large
fied by the overflow of water. The results of this simulation deformation and failure of soil can be handled in SPH
seem reasonable for both soil and water behavior. The soil without any difficulties. The effect of pore water pressure
particles near the boundary are stable during the calcula- and seepage force could be displayed well through SPH
tion process and only particles around the jet flow are simulation. Although numerical results have not been
affected. We have performed a number of simulation tests quantitatively compared with experimental data due to
using different velocities of water jet, from vjet = 1 m/s to the difficulty of experiment, the calculations are stable
vjet = 100 m/s, the same phenomena as the current simula- and the results appear acceptable throughout. We are
tion are obtained over a wide range of jet velocity. encouraged by these results, but recognize the need for
Fig. 5 shows the erosion process of soil excavation on improvement to the constitutive modeling of soil, including
the saturated soil model at representative times. As men- saturated soil. Advantages of the method are its robust-
tioned above, the saturated soil consisted of two-phases, ness, conceptual simplicity, relative ease of incorporating
water and soil, and each phase has to handle separately new physics, and especially its potential to handle large
using its own SPH model in order to obtain the field vari- deformation and failure.
ables of each particle. After the field variables of particles
in each phase are computed, the interaction between soil
References
and water can performed. The water pressure is permitted
to contribute to the pressure of soil particles through the [1] Lucy LB. A numerical approach to the testing of the fission
SPH summation, and thus the pore water pressure can be hypothesis. Astron J 1977;82:1013–24.
simulated. The seepage force is also computed as a function [2] Gingold RA, Monaghan JJ. Smoothed particle hydrodynamics:
of the relative velocity between soil and water particles. theory and application to non-spherical stars. Mon Not Royal
Astron Soc 1977;180:375–89.
The results of excavation on the saturated soil model are
[3] Libersky LD, Petschek AG. Smoothed particle hydrodynamics with
somewhat similar to those on the dry soil, except for the strength of material. In: Proceedings of the next free lagrange
behavior of soil particles in the soil box. In saturated soil conference vol. 395;1991. p. 248–57
simulation, the soil particles tend to run over the soil [4] Libersky LD, Petschek AG, Carney TC, Hipp JR, Allahdahi FA.
box, as seen in the last two images, while this phenomenon High strain Lagrangian hydrodynamics a three-dimensional
SPH code for dynamic material response. J Comput Phys
was not observed in the dry soil simulation. Computation-
1993;109:67–75.
ally, this difference can be explained by the effect of water [5] Monaghan JJ. Simulating free surface flows with SPH. J Comput
in the saturated soil: water is an incompressible material Phys 1994;110:399–406.
(slightly compressible in SPH simulation); therefore, under [6] Takeda H, Miyama M, Sekiya M. Numerical simulation of viscous
the high pressure of the water jet, water particles will tend flow by smoothed particle hydrodynamics. Prog Theor Phys
1994;92(5):939–60.
to run over the soil box; because of the presence of the
[7] Morris JP, Fox PJ, Zhu Y. Modeling low Reynolds
seepage force in the momentum equations (41), (42), these number incompressible flows using SPH. J Comput Phys
water particles will force soil particles to move together, so 1997;136:214–26.
the saturated soil mixture is effectively incompressible. Fur- [8] Benz W, Asphaug E. Simulation of brittle solid using
thermore, comparing the simulation results between two smoothed particle hydrodynamics. Comput Phys Commun
1995:253–65.
cases, it could be seen that the excavation by a water jet
[9] Biot MA. Theory of propagation of elastic waves in a fluid saturated
on dry soil is faster than on a saturated soil. porous solid. I. low frequency range. J Acoust Soc Am
1956;28:168–78.
4. Conclusion [10] Monaghan JJ, Lattanzio JC. A refined particle method for astro-
physical problems. Astron Astrophys 1985;149:135.
[11] Liu GR, Liu MB. Smoothed particle hydrodynamics: a meshfree
The development of smoothed particle hydrodynamics
particle method. World Scientific; 2003.
to simulate the behavior of soil–water interaction has been [12] Monaghan JJ. Smoothed particle hydrodynamics. Annual Rev
described through this paper. Water was modeled as a vis- Astron Astrophys 1992;30:543–74.