Anura3D MPM Research Community

Anura3D MPM Software

Scientific Manual

Version: 2022
21 April 2022
Anura3D MPM Software, Scientific Manual

With the contributions of:

Francesca Ceccato Università degli Studi di Padova, Italy
James Fern Amberg Engineering Ltd / EPFL, Switzerland
Mario Martinelli Deltares Delft, The Netherlands
Alexander Rohe Deltares Delft, The Netherlands
Alba Yerro VirginiaTech, United States

Contributors are listed in alphabetical order. Many people contributed to the creation
of this tutorial. If you feel that you are one of them but your name does not appear
in this list, please contact us.

Published and printed by: Contact and information:

Anura3D MPM Research Community e-mail: [email protected]
c/o Stichting Deltares web: www.Anura3D.com
Boussinesqweg 1
2629 HV Delft (The Netherlands)
Contents v

Contents
1 Introduction 3

2 General Framework 5

3 Mathematical formulation 7
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
3.2 One-phase Single-point formulation for solid analysis . . . . . . . . . . . . . 7
3.2.1 Governing equations . . . . . . . . . . . . . . . . . . . . . . . . . . 7
3.2.1.1 Mass conservation . . . . . . . . . . . . . . . . . . . . . 7
3.2.1.2 Momentum conservation . . . . . . . . . . . . . . . . . . 7
3.2.2 Weak form . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
3.3 One-phase Single-point formulation for liquid analysis . . . . . . . . . . . . 8
3.3.1 Governing equations . . . . . . . . . . . . . . . . . . . . . . . . . . 9
3.3.1.1 Mass conservation . . . . . . . . . . . . . . . . . . . . . 9
3.3.1.2 Momentum conservation . . . . . . . . . . . . . . . . . . 9
3.3.2 Weak form . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
3.4 Two-phase Single-point formulation for saturated soil . . . . . . . . . . . . . 10
3.4.1 Governing equations . . . . . . . . . . . . . . . . . . . . . . . . . . 10
3.4.1.1 Mass conservation . . . . . . . . . . . . . . . . . . . . . 11
3.4.1.2 Momentum conservation . . . . . . . . . . . . . . . . . . 12
3.4.1.3 Boundary conditions . . . . . . . . . . . . . . . . . . . . 12
3.4.2 Weak form . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
3.5 Two-phase Single-point formulation for unsaturated soil . . . . . . . . . . . . 13
3.5.1 Governing equations . . . . . . . . . . . . . . . . . . . . . . . . . . 14
3.5.1.1 Mass conservation . . . . . . . . . . . . . . . . . . . . . 14
3.5.1.2 Hydraulic constitutive equation . . . . . . . . . . . . . . . 15
3.5.1.3 Momentum conservation . . . . . . . . . . . . . . . . . . 15
3.5.1.4 Boundary conditions . . . . . . . . . . . . . . . . . . . . 16
3.5.2 Weak form . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

4 Explicit-Dynamic Formulation 19
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
4.2 Notations and variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
4.3 One-phase Single-point formulation for solid analysis . . . . . . . . . . . . . 19
4.3.1 Discretized momentum balance equation . . . . . . . . . . . . . . . 19
4.3.2 Initialization of material points . . . . . . . . . . . . . . . . . . . . . 21
4.3.3 Calculation of internal forces . . . . . . . . . . . . . . . . . . . . . . 23
4.3.4 Calculation of external forces . . . . . . . . . . . . . . . . . . . . . 24
4.3.5 Mass matrix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
4.3.6 Time discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
4.3.7 Solution algorithm for single time step . . . . . . . . . . . . . . . . . 26
4.3.8 Considerations on the use of one-phase solid analysis with porous
media . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
4.4 One-phase Single-point formulation for liquid analysis . . . . . . . . . . . . . 28
4.4.1 Discretized momentum balance equation of liquid . . . . . . . . . . 28
4.4.2 Calculation of internal forces . . . . . . . . . . . . . . . . . . . . . . 28
4.4.3 Calculation of external forces . . . . . . . . . . . . . . . . . . . . . 29
4.4.4 Mass matrix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
4.4.4.1 Time discretization . . . . . . . . . . . . . . . . . . . . . 30
vi Anura3D MPM Software, Scientific Manual

4.4.5 Solution algorithm for single time step . . . . . . . . . . . . . . . . . 31

4.5 Two-phase Single-point formulation for saturated soil . . . . . . . . . . . . . 31
4.5.1 Discretized momentum balance equation of liquid phase . . . . . . . 32
4.5.2 Discretized momentum balance equation of mixture . . . . . . . . . 32
4.5.3 Calculation of internal forces . . . . . . . . . . . . . . . . . . . . . . 33
4.5.4 Calculation of external forces . . . . . . . . . . . . . . . . . . . . . 33
4.5.5 Lumping procedure: Mass matrix and Drag matrix . . . . . . . . . . 34
4.5.6 Time discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
4.5.7 Solution algorithm for single timestep . . . . . . . . . . . . . . . . . 36
4.6 Two-phase Single-point formulation for unsaturated soil . . . . . . . . . . . . 37
4.6.1 Discretized momentum balance of the liquid . . . . . . . . . . . . . 37
4.6.2 Discretized momentum balance of the mixture . . . . . . . . . . . . 38
4.6.3 Time discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
4.6.4 Solution algorithm for single time step . . . . . . . . . . . . . . . . . 39
4.7 Stability criterion: the Critical timestep . . . . . . . . . . . . . . . . . . . . . 39
4.7.1 Criterion for One-phase-solid . . . . . . . . . . . . . . . . . . . . . 40
4.7.2 Criterion for one-phase-liquid . . . . . . . . . . . . . . . . . . . . . 40
4.7.3 Criterion for two-phase-coupled . . . . . . . . . . . . . . . . . . . . 40
4.8 4-node tetrahedral element . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

5 The boundary conditions 45

5.1 Zero kinematic boundary conditions . . . . . . . . . . . . . . . . . . . . . . 47
5.2 Traction boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . . . 47
5.3 Hydraulic boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . . 49
5.3.1 Total head BC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
5.3.2 Infiltration BC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
5.3.3 Potential seepage face BC . . . . . . . . . . . . . . . . . . . . . . . 52
5.4 Absorbing boundaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52

6 Advanced concepts 55
6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
6.2 Contact problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
6.2.1 Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
6.3 Rigid body . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
6.4 Local damping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
6.4.1 Single-phase problem . . . . . . . . . . . . . . . . . . . . . . . . . 62
6.4.2 Two-phase problem . . . . . . . . . . . . . . . . . . . . . . . . . . 63
6.5 Convergence to static equilibrium . . . . . . . . . . . . . . . . . . . . . . . 63
6.6 Mass scaling for quasi-static problems . . . . . . . . . . . . . . . . . . . . . 64
6.7 Cell-crossing instability and MPM-MIXED . . . . . . . . . . . . . . . . . . . 64
6.8 Pathological locking in low order elements . . . . . . . . . . . . . . . . . . . 65
6.8.1 Strain smoothing procedure . . . . . . . . . . . . . . . . . . . . . . 65
6.8.2 Pore pressure increment smoothing . . . . . . . . . . . . . . . . . . 66
6.9 Liquid free surface detection . . . . . . . . . . . . . . . . . . . . . . . . . . 67
6.10 Moving mesh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68

Bibliography 71
Contents 1

Anura3D - is a software for the numerical modelling and simulation of large deformations and
soil–water–structure interaction using the material point method (MPM). Copyright (C) 2020 Mem-
bers of the Anura3D MPM Research Community.
Anura3D is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser
General Public License as published by the Free Software Foundation <https://www.gnu.org/licenses/>,
either version 3 of the License, or (at your option) any later version.
Anura3D documentation is furnished under the License and may be used only in accordance with
the terms of such license. It is advised to consult the manuals before applying the software.

Anura3D is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without
even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
See the GNU Lesser General Public License for more details.

Anura3D MPM Research Community Contact and information

c/o Stichting Deltares e-mail: [email protected]
Boussinesqweg 1 web: www.Anura3D.com
2629 HV Delft (The Netherlands)
2 Anura3D MPM Software, Scientific Manual
1 Introduction
This manual presents the the theories and numerical methods used in Anura3D. It contains
a chapter on the mathematical formulation and the dynamic-explicit material point formulation
of the 1-phase solid, 1-phase liquid and 2-phase mixture. The boundary conditions and ad-
vanced concepts used in the current version of Anura3D are also described in two separated
sections. Lastly, a dedicated chapter on material models available in this code is also added
for completeness. Additional scientific material can be found in the Publications section of the
website.
4 Anura3D MPM Software, Scientific Manual
2 General Framework
The material point method (MPM) was initially developed by Harlow [1] in Los Alamos National
Laboratory. The basis of that method, called initially Particle-In-Cell method (PIC), was to
represent fluid flow by a set of material points moving through a background fixed mesh. Later,
in 1994 in the New Mexico University, Sulsky and Schreyer [2–4] extended the approach for
modelling problems of solid mechanics, where basically the equation of dynamic momentum
balance is discretized.

MPM is considered as a method in between the particle-based methods and the Finite Ele-
ment Method (FEM). This is because it discretizes the media in two different frames. First,
the continuum is divided into a set of material points (MPs). Each material point represents
a portion of the domain and the mass of such sub-domain, assumed to be fixed during all of
the calculation (to ensure mass conservation). In the classical MPM, it is considered to be
concentrated at the corresponding material point. Other quantities such as velocities, strains
and stresses, are also initialized and carried by the MPs. Each material point moves attached
with the deformations of the body and this provides the Lagrangian description of the media.

The second frame is a computational mesh. It is the same as the one used in the conven-
tional FEM and it is built to cover the full domain of the problem. The discrete governing
equations are solved at the nodes of the computational mesh. The variables required to solve
the equations in the mesh at any step of the analysis are transferred from the material points
to the nodes of the mesh by using mapping functions. For instance, the typical linear shape
functions used in the FEM. Boundary conditions can be imposed at the mesh nodes or at
the material points and the governing equations are solved by using an incremental scheme.
Then the quantities carried by the material points are updated through the interpolation of the
mesh results, using the same mapping functions. The information associated with the mesh
is not required for the next step of the analysis; therefore it can be discarded avoiding any
mesh distortion.

Fig. 2.1 is a typical problem setup. The material domain is represented by a set of material
points. The computational mesh covers all the computational domain.

At the beginning of each time increment, the discretized governing equations are defined
by mapping information from the MPs to the computational nodes of the mesh by means of
the interpolation functions (Fig. 2.2a). The governing equations of motion are solved for the
primary unknown variables, e.g. the nodal accelerations (Fig. 2.2b). This nodal values are
used to update acceleration, velocity and position of MPs, as well as to compute strains and
stresses at the MPs (Fig. 2.2c). No permanent information is stored in the mesh, thus it can
be freely redefined at the end of the time step, but commonly it is kept fixed. The assignment
of MPs to finite elements is updated after mesh adjustment (Fig. 2.2d).
6 Anura3D MPM Software, Scientific Manual

Figure 2.1: Space discretisation. Nodes of the computational mesh and material points,
after [5]

Figure 2.2: Computation scheme of MPM: (a) map information from MPs to grid nodes;
(b) solve governing equations of motion at the nodes; (c) update MP informa-
tion; (c) update MP position and housekeeping
 3 Mathematical formulation

3.1 Introduction
This chapter presents the basic equations for the dynamic equilibrium of the solid body and
the porous media within the framework of continuum mechanics. The variational form, or
weak form, of the corresponding partial differential equations is also developed.

3.2 One-phase Single-point formulation for solid analysis

The MPM formulation for mechanical problems (1-phase MPM formulation) is presented in this
Section. It is based on [3], in which only a single phase is considered. Because the MPM can
be viewed as an extension of the FEM procedure, is important to note that they have much in
common. For instance, the weak form of the governing equations as well as the final system
of equations posed at the nodes of the computational mesh are identical in both schemes.
The most distinguishing feature between them is the calculation of the numerical integrations
over the volume of a finite element. Meanwhile in FEM the integrations are carried out using
Gaussian quadrature, in MPM they are calculated based on the material points.

Three-dimensional domains are considered in this Chapter. Vectors and tensors are identified
with bold type.

3.2.1 Governing equations

From a mathematical point of view, the continuum can be described by a set of differential gov-
erning equations. This includes conservation of mass, conservation of momentum, balance of
energy and the corresponding boundary conditions. Furthermore, constitutive equations such
as stress-strain relationship, which are characteristic of the material forming the continuum,
are needed.

3.2.1.1 Mass conservation

The conservation of mass of the solid phase in Lagrangian form is written with the following
partial differential equation:

d
ρ + ρ∇ · vS = 0 (3.1)
dt
in which v is the velocity vector of the solid phase. One of the features of the MPM formulation
for a 1-phase material is that the mass of each material point mMP remains constant during
all of the calculation. This fact implies that mass conservation is automatically satisfied.

3.2.1.2 Momentum conservation

The momentum conservation of the continuum can be written with the differential Eq. 3.2.
A dynamic formulation is considered, which means that the acceleration term is taken into
account.
dvS
ρ = ∇ · σ S + ρb (3.2)
dt
 8 Anura3D MPM Software, Scientific Manual

dv
The density of the continuum is represented by ρ, dtS is the acceleration, σ S is the Cauchy
stress tensor, and b is the body force vector.

Two kind of boundary conditions are defined: prescribed traction (Eq. 3.3) and prescribed
displacements (Eq. 3.4). Each one is applied on the corresponding domain, ∂Ωσ and ∂Ωu
respectively,

σ S (x, t) · n = btS (t) (3.3)

uS (x, t) = u
bS (t) (3.4)

where σ S (x, t) is the stress tensor, n is the outward unit normal vector of the free surface, and
btS (t) is the surface traction vector. x is the position vector, uS is the displacement vector, and
t is time.

3.2.2 Weak form

The conservation of the momentum (Eq. 3.2) is a key equation in continuum mechanics, rep-
resenting the equation of motion of the continuum. To discretize this equation, its strong form
has to be transformed into the weak form or the well-known virtual work equation. The mo-
mentum equation is now multiplied by a test function or a virtual velocity δ v and is integrated
over the current domain Ω occupied by the continuum, i.e.,
Z Z Z
dv
δ vS ρ S dΩ = δ vS (∇ · σ S )dΩ + δ vS ρgdΩ (3.5)
dt
Ω Ω Ω

with δ v = 0 on ∂Ωu .

The first term on the right-hand side of Eq. (3.5) can be expressed as
Z Z Z
δ vS (∇ · σ )S dΩ = (∇ · δ vS σ S ) dΩ – (∇ · δ vS )σ S dΩ (3.6)
Ω Ω Ω

Using Gauss theorem, also called divergence theorem, the first term on the right-hand side of
Eq. (3.6) can be written as
Z Z
(∇ · δ vS σ S ) dΩ = δ vS τ S d∂Ω (3.7)
Ω ∂Ωσ

Substituting all previous terms in Eq. (3.5), it follows

Z Z Z Z
dv
δ vS ρ S dΩ = δ vS τ S d∂Ω – (∇ · δ vS )σ S dΩ + δ vS ρgdΩ (3.8)
dt
Ω ∂Ωσ Ω Ω

which will be used then in the formulation of the discrete equations.

3.3 One-phase Single-point formulation for liquid analysis

 Mathematical formulation 9

3.3.1 Governing equations

In this section, the conservation equations of mass and linear momentum are written with
respect to the liquid phase.

3.3.1.1 Mass conservation

The conservation of mass of the liquid phase in a lagrangian form is written as follows:
dρL
+ ρL ∇ · vL = 0 (3.9)
dt
in which vL is the velocity vector of the liquid phase. As written for the solid phase, one of
the features of the MPM formulation is that the mass of each material point mMP remains
constant during all of the calculation. This fact implies that mass conservation is automatically
satisfied.

3.3.1.2 Momentum conservation

The momentum conservation of the continuum can be written with the differential Eq. 3.10.
A dynamic formulation is considered, which means that the acceleration term is taken into
account.
dvL
ρL = ∇ · σ L + ρL b (3.10)
dt
The liquid density of the continuum is represented by ρL , vL is the velocity, σL is the Cauchy
stress tensor, and b is the body force vector.

Two kind of boundary conditions are defined: prescribed traction (Eq. 3.11) and prescribed
displacements (Eq. 3.12, Fig. 3.1). Each one is applied on the corresponding domain, ∂Ωσ
and ∂Ωu respectively. The boundary condition along ∂Ωσ is defined as follows:

σ L (x, t) · n = bt(x, t) (3.11)

where σ L (x, t) is the stress tensor, n is the outward unit normal vector of the free surface and
bt(x, t) is the surface traction vector. The boundary condition along ∂Ωu is defined as follows:

u(x, t) = u
b(t) (3.12)

x is the position vector, u

b is the displacement vector and t is time.

The boundary between liquid and gas phase is the free surface and is modelled as a boundary
condition on ∂Ωσ . In case of zero surface tension and setting the atmospheric pressure to
zero, such a boundary condition is set as bt(x, t) = p(x, t) = 0.

3.3.2 Weak form

The discretization of the equation of conservation of linear momentum for pure liquid is written
as follows. The strong form of the equation is transformed into the weak form or the well-known
virtual work equation. The momentum equation is now multiplied by a test function or a virtual
velocity δ v and is integrated over the current domain Ω occupied by the continuum, i.e.,
Z Z Z
dv
δ vL ρ L dΩ = δ vL (∇ · σ L )dΩ + δ vL ρgdΩ (3.13)
dt
Ω Ω Ω
 10 Anura3D MPM Software, Scientific Manual

Figure 3.1: Water column collapse. Boundary conditions: displacements on ∂Ωu and
pressure on ∂Ωσ .

with δ v = 0 on ∂Ωu .

The first term on the right-hand side of Eq. (3.13) can be expressed as
Z Z Z
δ vL (∇ · σ L ) dΩ = (∇ · δ vL σ L ) dΩ – (∇ · δ vL )σ L dΩ (3.14)
Ω Ω Ω

Using Gauss theorem, also called divergence theorem, the first term on the right-hand side of
Eq. (3.14) can be written as
Z Z
(∇ · δ vL σ L ) dΩ = δ vL τ L d∂Ω (3.15)
Ω ∂Ωσ

Substituting all previous terms in Eq. (3.5), it follows

Z Z Z Z
dv
δ vL ρ L dΩ = δ vL τ L d∂Ω – (∇ · δ vL )σ L dΩ + δ vL ρgdΩ (3.16)
dt
Ω ∂Ωσ Ω Ω

which will be used then in the formulation of the discrete equations.

3.4 Two-phase Single-point formulation for saturated soil

3.4.1 Governing equations

This section is devoted to the presentation of the governing equations that are implemented
in Anura3Dfor the solution of coupled dynamic, two-phase problems. It’s worth mentioning
that several coupled two-phase formulation have been discussed in the literature since the
pioneering work by Terzaghi and Biot. The reader can refer to the work by Zienkiewicz and
co-workers ([6], [7] and [8]) for a comprehensive discussion on this topic.
 Mathematical formulation 11

With the formulation used in Anura3D, the equilibrium equations are solved for the acceler-
ations of pore liquid and soil skeleton as the primary unknown variables. This formulation
proved to be able to capture the physical response of saturated soil under dynamic as well as
static loading [9]. The governing equations are the mass conservation, and the momentum
conservation of each phase and the mixture.

3.4.1.1 Mass conservation

The mass conservation of the solid phase is expressed as:

d
[(1 – n)ρS ] + (1 – n)ρS ∇ · vS = 0 (3.17)
dt
in which vS is the velocity vector of the solid phase and n is the porosity of the solid skeleton.

On denoting the components of the vector of the (true) velocity of the liquid phase as vL , the
conservation of mass of this phase can be written as:

d(nρL )
+ nρL ∇ · vL = 0 (3.18)
dt

When considering incompressible solid grains and disregarding the spatial variations in den-
sities and porosity, one can reduce the expression for the conservation of mass of the solid
and liquid phases to

dn
– + (1 – n)∇ · vS = 0 (3.19)
dt
and
dn dρL
ρL +n + nρL ∇ · vL = 0 (3.20)
dt dt

respectively. Substituting Equation 3.19 into Equation 3.20 allows to eliminate the term dn
dt .
Hence,

n dρL
(1 – n)∇ · vS + + n∇ · vL = 0 (3.21)
ρL dt

The liquid is considered as weakly compressible material. The effective volumetric strain ε̄L
in the liquid is defined as

dε̄L 1 dρL
=– (3.22)
dt ρL dt

Substituting Equation 3.22 into Equation 3.21 and rearranging terms yields:
 
dε̄L 1
= (1 – n)∇ · vS + n∇ · vL (3.23)
dt n

Equation 3.23 represents the conservation of mass of the saturated soil. It is also known as
storage equation.
 12 Anura3D MPM Software, Scientific Manual

3.4.1.2 Momentum conservation

The conservation of momentum of the solid phase can be expressed as

dvS n2 ρL g
(1 – n)ρS – ∇ · σ 0 – (1 – n)∇pL – (1 – n)ρS g – (vL – vS ) = 0 (3.24)
dt k
where k is the isotropic Darcy permeability. It can be expressed in terms of the intrinsic
permeability κ and the dynamic viscosity of the liquid µd as:
ρL g
k=κ (3.25)
µd

The conservation of momentum of the liquid phase is written as:

dvL n2 ρL g
nρL – n∇pL – nρL g + (vL – vS ) = 0 (3.26)
dt k
the term (vL – vS ) represents the relative velocity of the liquid respect to the solid.

Adding the momentum of the solid phase (Eq. 3.24) to the momentum of the liquid phase
(Eq. 3.26), the momentum conservation for the mixture can be written as:

dvL dvL
(1 – n)ρS + nρL = ∇ · σ + ρsat g (3.27)
dt dt

Summarizing, the two-phase problem is described by two momentum equations, i.e. 3.26 for
the liquid and 3.27 for the mixture, the storage Eq. 3.23, and the constitutive equation for the
soil skeleton. These equations are derived neglecting the spatial variation of densities and
porosity, assuming incompressible soil grains, and assuming the validity of the Darcy’s law.

3.4.1.3 Boundary conditions

The proposed formulation requires that the boundary of the domain is the union of the follow-
ing components

∂Ω = ∂Ωu ∪ ∂Ωσ = ∂ΩvL ∪ ∂Ωp (3.28)

where ∂ΩvL and ∂Ωp are the prescribed velocity and prescribed pressure boundaries of the
liquid phase, respectively, whereas ∂Ωu is the prescribed displacement (velocity) boundary of
the solid phase and ∂Ωσ is the prescribed total stress boundary.

The following conditions should also be satisfied at the boundary

∂Ωu ∩ ∂Ωσ = ∅ and ∂ΩvL ∩ ∂Ωp = ∅ (3.29)

3.4.2 Weak form

Before the discretization, the strong form of the governing equations has to be transformed in
the weak form. This is achieved by multiplying Equations 3.26 and 3.27 by weighting function
t and integrating over the current domain Ω:
Z Z Z Z
dv nρ g
δ tρL L dΩ = δ t∇pL dΩ + δ tρL gdΩ – δ t L (vL – vS )dΩ (3.30)
Ω dt Ω Ω Ω k
 Mathematical formulation 13

Figure 3.2: Displacement and traction boundary conditions for two phase problem (Al-
Kafaji 2013) [10].

Z Z Z Z
dv dv
δ t(1 – n)ρS S dΩ = δ t∇ · σ 0 dΩ + δ tρsat gdΩ – δ tnρL L dΩ (3.31)
Ω dt Ω Ω Ω dt

Applying the divergence theorem and the traction boundary conditions, the final weak forms
are:
Z
dvL
δ tρL dΩ =
Ω dt
Z Z Z Z
nρL g
δ tp̄L dS – ∇δ tpL dΩ + δ tρL gdΩ – δt (vL – vS )dΩ (3.32)
∂Ωp Ω Ω Ω k

Z
dvS
δ t(1 – n)ρS dΩ =
Ω dt
Z Z Z Z
dvL
+ δ tτ dS – ∇δ tσ dΩ + δ tρsat gdΩ – δ tnρL dΩ (3.33)
∂Ωσ Ω Ω Ω dt

The left-hand side terms in Equations 3.32 and 3.33 represent the inertia. In the right-hand
side, the first terms represent the external load applied at the boundary, the second terms
represent the internal force, and the third terms represent the gravity load. The last term in
Equation 3.32 is the drag force. The last term in Equation 3.33 is the liquid inertia, where the
porosity is taken into account.

3.5 Two-phase Single-point formulation for unsaturated soil

The two-phase formulation for unsaturated soil implemented in Anura3D refers to the work
by [11]. It is derived considering dynamic momentum balance of the liquid phase, dynamic
momentum balance of the mixture, and mass balances and constitutive relationships of solid
 14 Anura3D MPM Software, Scientific Manual

and liquid phases, for this reason it can be considered a two-phase formulation and it is
an extension of the formulation presented in Section 3.4. All dynamic terms are taken into
account; acceleration of the solid skeleton aS and acceleration of the pore liquid aL are the
primary unknowns. Unsaturated conditions are accounted for considering that the liquid phase
does not entirely fill the voids.

3.5.1 Governing equations

3.5.1.1 Mass conservation

When considering incompressible solid grains and disregarding the spatial variations in den-
sity and porosity, the expressions for the conservation of mass of the solid and the liquid
reduce to Eqs. 3.34 and 3.35 respectively.

dnS
+ nS div(vS ) = 0 (3.34)
dt

d (ρL nL )
= (vS – vL )∇(nL ρL ) – nL ρL div(vL ) (3.35)
dt

where nL and nS are the liquid and solid volumetric fraction, respectively.

Including Eq. 3.34 into Eq. 3.35, taking into account the definitions of liquid volumetric con-
centration ratios in terms of porosity and degree of saturation, nL = SL n, nS = 1 – n, and
rearranging terms gives Eq. 3.36

d (ρL SL )  
n = div ρL nSL (vS – vL ) – ρL SL div(vS ) (3.36)
dt

Finally, the material derivative in Eq. 3.36 is solved assuming liquid pressure as state variable,
which yields to Eq. 3.37.

 
∂ρL ∂ SL dpL  
n SL + ρL = div ρL nSL (vS – vL ) – ρL SL div(vS ) (3.37)
∂ pL ∂ pL dt

The derivative of liquid density with respect to pressure is given by the state equation of
liquid; while the derivative of the degree of saturation is given by the soil-water retention curve
(SWRC) (see Sec. 3.5.1.2).
 Mathematical formulation 15

3.5.1.2 Hydraulic constitutive equation

The relationship between pore liquid pressure and degree of saturation or liquid content is
essential to model the behaviour of unsaturated soil. This is given by the SWRC which can
assume different analytical expressions. In Anura3D, two alternative relationships are used:
(i) a linear relation (Eq. 3.38) where av is a constant parameter, Smin is the residual degree
of saturation and Smax is the maximum degree of saturation, and (ii) the Van Genuchten
relationship (Eq. 3.39) where pref , λ, are fitting parameters.

SL = Smin + (Smax – Smin )pL av (3.38)

   1 –λ
pL 1–λ
SL = Smin + (Smax – Smin ) 1 + (3.39)
pref

Partial saturation modifies the soil hydraulic conductivity; in general partially saturated soils
are less permeable than fully saturated soils. The ratio between actual hydraulic conductivity
and saturated hydraulic conductivity k/ksat is a function of the degree of saturation as given
by the hydraulic conductivity curve (HCC). A number of relationships have been proposed in
the literature, the functions proposed by Hillel [12] (Eq. 3.40) and Mualem [13] (Eq. 3.41) are
implemented in the applied software, in which r and λ are fitting parameters.

k
= SrL (3.40)
ksat

k p h 1
λ i2
= SL 1 – 1 – SLλ (3.41)
ksat

3.5.1.3 Momentum conservation

The momentum balance of the liquid phase per unit of liquid volume is given in Eq. 3.42,

ρL aL = ∇pL – fdL + ρL g (3.42)

where ρL is the liquid density, fdL is the drag force which accounts for solid-fluid interaction, pL
is the liquid pressure and g is the gravity vector.

The flow is considered laminar and stationary in the slow velocity regime. Hence, the drag
force is governed by Darcy’s law (Eq. 3.43),
nL µL
fLd = (vL – vS ) (3.43)
κL
where µL is the dynamic viscosity of the liquid, κL is the liquid intrinsic permeability, nL is the
liquid volumetric fraction and vL , vS are the absolute liquid and solid velocities respectively.
 16 Anura3D MPM Software, Scientific Manual

The isotropic intrinsic permeability of the liquid κL can also be expressed in terms of Darcy
permeability, or hydraulic conductivity, kL (Eq. 3.44).
µL
κL = kL (3.44)
ρL g

The mixture dynamic momentum conservation can be written as Eq. 3.45,

nS ρS aS + nL ρL aL = div(σ ) + ρm g (3.45)

where ρS is the solid grain density, nS is the volumetric concentration ratio of solid, and ρm =
nS ρS + nL ρL is the density of the mixture. Note that nS = 1 – n and nL = SL n, where n is the
porosity of the solid skeleton and SL is the degree of saturation. σ is the total stress tensor,
which can be computed with the Bishop’s effective stress equation for unsaturated soils and
has the form of Eq. 3.46, where χ is an effective stress parameter, here assumed equal to SL ,
and m is the unit vector, equal to (1 1 1 0 0 0)T in 3D. In this paper, stresses and pressures
are positive for tension, thus suction is equal to pL .

σ = σ 0 + χpL m (3.46)

3.5.1.4 Boundary conditions

The proposed formulation requires the definition of the following boundary conditions (BC):

1 Prescribed liquid displacement or velocity on ∂ΩvL and prescribed pressure on ∂Ωp ,

where ∂Ω = ∂ΩvL ∪ ∂Ωp and ∂ΩvL ∩ ∂Ωp = 0
2 Prescribed solid displacement or velocity on ∂Ωu and prescribed traction on ∂Ωτ , where
∂Ω = ∂ΩvS ∪ ∂Ωτ and ∂ΩvS ∩ ∂Ωτ = 0.

Essential boundary conditions on ∂ΩvL and ∂ΩvS are imposed on the nodes of the computa-
tional grid. Natural boundary conditions on ∂Ωp and ∂Ωτ are included in the weak form of the
momentum balance equations.

In typical problems with partially saturated soil in geomechanics, prescribed liquid velocity
can be applied on infiltration boundaries, and prescribed pressures can be applied either by
defining a pressure load bpL or by assigning a total hydraulic head H.
b Assuming the validity of
Bernoulli’s equation and neglecting the kinematic head, the total hydraulic head can be written
as
bpL
H
b = hg – (3.47)
ρL g
pL /(ρL g) is the pressure head.
where hg stands for the potential head or geometric head and b
The minus sign in Eq. 5.5 is introduced because pressure is assumed negative for compres-
sion.

Sometimes, the boundary condition is part of the problem solution, i.e. it is not known a priori
if the boundary belongs to ∂ΩvL or ∂Ωp . This is typical of free surface flows across porous
media along the so-called potential seepage face.

Figure 3.3 represents schematically how these BCs simulate different hydraulic loading acting
on a levee. The implementation details of these boundary conditions can be found in 5.3.
 Mathematical formulation 17

(a) Rainfall
Infiltration BC

Potential
River level seepage face BC

Total head BC

Impermeable BC 𝑣𝑣𝐿𝐿 = 0

Impermeable layer

(b) Load
Traction BC

Fixities 𝑣𝑣𝑆𝑆 = 0

Figure 3.3: Typical boundary conditions for liquid phase (a) and mixture (b), from [11].

3.5.2 Weak form

The governing equations (liquid and mixture momentum balance) are transformed in the weak
form by multipling by a weighting function w and integrating over the domain Ω
Z Z h
n µ i
ρL aL · δ tdV = ∇pL + ρL g – L L (vL – vS ) · δ tdV (3.48)
κL
V V

The application of the divergence theorem to Eq. 3.48 results in Eq. 3.50

Z Z Z Z
ρL aL · δ tdV = _∂Ωp p bL d∂Ωp – pL ∇ · δ tdV + ρL g · δ tdV– (3.49)
ZV h V V
nL µL i
(vL – vS ) · δ tdV (3.50)
V κL

The mixture momentum balance equation in the weak form corresponds to Eq. 3.51
Z Z h i
[(1 – n)ρS aS ] · udV = – nL ρL aL + ∇σm + ((1 – n)ρS + nL ρL )g · δ tdV (3.51)
V V

The application of the divergence theorem to Eq. 3.51 results in Eq. 3.53
Z Z h Z
[(1 – n)ρS aS ] · δ tdV = – (nL ρL aL ) · δ tdV + τ dδ t∂Ωt (3.52)
V V ∂Ωt
Z Z
– σm : (∇δ t)dV + [((1 – n)ρS + nL ρL )g] · δ tdV (3.53)
V V

Where p bL is the prescribed liquid pressure in the boundary ∂Ωp and τ is the prescribed
traction vector in the boundary ∂Ωτ .
18 Anura3D MPM Software, Scientific Manual
 4 Explicit-Dynamic Formulation

4.1 Introduction
This chapter presents the computation scheme adopted in the dynamic explicit implementa-
tion of Anura3D. It explains how the system of solving equations of the problem is assembled
within the MPM framework and the main computational steps of the MPM solution procedure.

4.2 Notations and variables

The strain and stress tensors can be represented in vector form, taking advantage of the
symmetry of the strain and stress tensors

ε(x, t) = [ε11 ε22 ε33 2ε12 2ε23 2ε31 ]T (4.1)

T
σ (x, t) = [σ11 σ22 σ33 σ12 σ23 σ31 ] (4.2)

where the first three terms in the vectors (σij and εij with i = j) are the normal components
along the axis xi of the coordinate system and the other three represent the shear terms acting
on the xi xj planes. This way to represent the stress and strain is known in the literature as
Voigt notation and is very useful as it reduces the order of the symmetric tensors.

It should be realized that the representation of the strain is different from stress in Voigt nota-
tion. The last three terms of the strain vector are represented by γij = 2εij . The reason for this
is to ensure that energy is preserved and that different expressions of energy using tensors
or vectors are equal,i.e.

εij σi j = εT σ (4.3)

4.3 One-phase Single-point formulation for solid analysis

4.3.1 Discretized momentum balance equation

The weak form of the linear momentum equation over the domain Ω is given by
Z Z Z Z
dv
δ vρ dΩ = δ v τ d∂Ω – (∇ · δ v)σ dΩ + δ vρgdΩ (4.4)
dt
Ω ∂Ωσ Ω Ω

The domain Ω is decomposed into finite subdomains Ωel called finite elements. The union of
these subdomains comprises the total domain, Ω = ∪nelm el=1 Ωel , where nelm denotes the total
number of finite elements in the mesh. Each element is jointed with its surrounding elements
by points called nodes. The state variable is assumed to have pre-defined interpolation func-
tions within the element and the solution is then obtained at the nodes. Hence, the equilibrium
is satisfied at the nodes. It is usual in the finite element method to use matrix notations in the
discretization of the virtual work equation.

In finite element method, the discrete form is obtained by approximating the displacement and
20 Anura3D MPM Software, Scientific Manual

eventually the velocity and acceleration as

u(x, t) ≈ N̄(x)û(t)
v(x, t) ≈ N̄(x)v̂(t) (4.5)
a(x, t) ≈ N̄(x)â(t)

respectively. The corresponding virtual quantities are approximated the same way, for in-
stance,

δ u(x, t) ≈ N̄(x)δ û(t) (4.6)

δ v(x, t) ≈ N̄(x)δ v̂(t) (4.7)
δ a(x, t) ≈ N̄(x)δ â(t) (4.8)

The interpolation function or the shape function matrix N̄ has the following form
 
N̄(x) = N̄1 (x) N̄2 (x) ... ... N̄nT (x) (4.9)

with
 
N̄i (x) 0 0
N̄i (x) =  0 N̄i (x) 0  (4.10)
0 0 N̄i (x)

with nT denoting the total number of nodes in the mesh. The bar superscript implies that the
shape functions are written in terms of the global coordinate system. The vectors of nodal
displacements, velocities and accelerations are denoted as û, v̂ and â, respectively. These
vectors can be written as follows

û(t) = [û11 û12 û13 ... ... ... ûnT 1 ûnT 2 ûnT 3 ]T (4.11)
T
v̂(t) = [v̂11 v̂12 v̂13 ... ... ... v̂nT 1 v̂nT 2 v̂nT 3 ] (4.12)
â(t) = [â11 â12 â13 ... ... ... ânT 1 ânT 2 ânT 3 ]T (4.13)

where v12 indicates the velocity at node 1 in the direction of coordinate x2 .

The kinematic relation can be written using matrix notation as

ε̇(x, t) = Lv(x, t) (4.14)

with L being a linear differential operator, that has the following form
 
∂ 0 0
∂ x1

 0 ∂ 0 

 ∂ x2 
 0 0 ∂ 
L= ∂
 ∂ x3  (4.15)
 ∂x ∂ 0 

 2 ∂ x1 
 0 ∂ ∂ 
 ∂ x3 ∂ x2 
∂ 0 ∂
∂ x3 ∂ x1

Substituting the second of Equation (4.3.1) into Equation (4.14) yields

ε̇(x, t) = LN(x)v̂(t) = B(x)v̂(t) (4.16)

 Explicit-Dynamic Formulation 21

in which B is the strain-displacement matrix that can be written as

B(x) = [B1 (x) B2 (x) B3 (x) ... ... ... BnT (x)] (4.17)

with
 ∂ N̄ (x) 
i
0 0
 ∂ x1 ∂ N̄i (x) 
 0 0 
 ∂ x2 

 0 ∂ N̄i (x) 
0 ∂ x3 

Bi (x) = 
 ∂ N̄i (x) ∂ N̄i (x) (4.18)
0 

 ∂ x2 ∂ x1
∂ N̄i (x) ∂ N̄i (x) 
 
 0
 ∂ x3 ∂ x2 
∂ N̄i (x) ∂ N̄i (x)
∂ x3 0 ∂ x1

The weak form yields as follow:

Z Z
T
δ v̂T N t · n d∂Ω – δ v̂
T
BT σ dΩ+
∂Ω Z Ω 
Z
T T T T
+δ v̂ N ρg dΩ – δ v̂ N ρN dΩ â = 0 (4.19)
Ω Ω

where n is the unit vector which is normal to the boundary of the domain.

The discretized form yields:

fext – fint = Mâ (4.20)

with

Z Z
ext T T
f = N t · n d∂Ω + ρN g dΩ (4.21)
Z ∂Ω Ω

fint = BT σ 0 dΩ (4.22)
Z Ω
T
M= ρN NdΩ (4.23)
Ω

where fext is the vector of nodal external forces, fint is the vector of nodal internal forces and
M is the nodal mass matrix.

4.3.2 Initialization of material points

The MP carry all the information of the continuum. In this section, we discuss the initialization
of MP within the background mesh. This includes association of mass, body forces, trac-
tions and other properties of the continuum to MP. Elements filled with MP are called active
elements and their nodes contribute to the solving system of equations; on the contrary the
empty elements, i.e. those without any MP, are ignored thus reducing the computational cost.

Let us consider a single tetrahedral element to explain the full procedure of initialization of
MP information. Each MP is initially positioned at a predefined local position inside the parent
22 Anura3D MPM Software, Scientific Manual

Figure 4.1: Initialization of surface traction. (a) tetrahedral element (b) triangular element
(c) traction is mapped to the boundary MP, after [10].

element, and hence the local position vector ξMP of the material point MP is initialized. The
global position vector xMP is then obtained as
noX
nodes
t
xMP (ξMP ) ≈ Ni (ξMP )xi (4.24)
i=1

in which nonodes denotes the number of nodes per element, Ni (ξMP ) is the shape function of
node i evaluated at the local position of material point MP and xi are the nodal coordinates.

Volumes associated with MP are calculated so that all the MP inside the element have initially
the same portion of the element volume, i.e.
Z noq,el
1 1 X
ΩMP = dΩ ≈ wMP |J(ξMP )| (4.25)
noMP,el Ωe noMP,el
q=1

where ΩMP is the volume associated with the material point MP, noMP,el denotes the number
of MP in the element, noq,el is the number of Gauss points in the element, wMP is the local
integration weight associated with Gauss point MP, and J is the Jacobian matrix. This implies
that, at the beginning of the calculation, all the active elements are assumed to be fully filled
by the continuum body. An element is said to be partially filled if the sum of the volumes of
the contained MP is less than the element volume.

The mass mMP is then calculated as

mMP = ΩMP ρMP (4.26)

with ρ being the mass density of the material to which the material point MP belongs.

grav
The gravity force fMP is simply calculated using the mass of the MP and the vector of gravi-
tational acceleration g as
grav
fMP = mMP g (4.27)

The external forces applied at the traction boundary are mapped to the MP located next to the
element border, also called boundary MP (Fig. 4.1).
 Explicit-Dynamic Formulation 23

These MP carry surface traction throughout the computations. Considering a tetrahedral ele-
ment, the traction vector τe applied at the triangular surface is interpolated from the nodes of
this surface to the boundary MP. Hence, the traction at boundary material point MP is
ntri
X
τe (xp ) ≈ Ni (ξq )τe (xi ) (4.28)
i=1

where Ni is the shape function of node i of the triangular surface element and ξq are the
coordinates of the boundary MP p inside the parent triangular element. These coordinates
simply represent the projection of the MP on the triangular surface element. The traction force
vector ftrac
p is then

ntri
Se Se X
ftrac
p = τe (xMP ) = Ni (ξMP )τe (xi ) (4.29)
nebMP nebMP
i=1

in which nebMP denotes the number of boundary MP located next to the loaded surface and
Se is the area of the corresponding loaded surface of the element e.

In the initialization of MP, initial conditions, material parameters and constitutive variables
are assigned to them as well. Furthermore, book-keeping is initialized at this step, including
information such as to which element each particle initially belongs and the initial number of
particles per each active finite element.

4.3.3 Calculation of internal forces

In finite element formulation the numerical integration of all integrals above is not performed
along the global coordinate system. Usually, it is more convenient to convert the global coor-
dinates (x) of each element into a reference (or parent) element system (ξ ), and the mapping
procedure is done using the Jacobian matrix. In this work, a unit 4-node tetrahedron with
linear shape functions is chosen to be the reference element. Details about the mapping
procedure are described in Section 4.8. The internal force becomes

no
X elm Z no
X elm Z
fint = BT σ dΩ = BT σ |J| d (4.30)
el=1 Ωel el=1 

where Ωel is the volume of the element el in the global coordinate system, noelm is the number
of active elements,  is volume of the parent element, d = dξ1 dξ2 dξ3 is the infinitesimal
volume in the parent element system, |J| is the determinant of the Jacobian, and B is the
matrix of the shape function gradients calculated at location ξ with respect to the parent
coordinate system. If numerical integration is applied, the equations yield as follows

no
X elm noNodes,el
X noX int,el

fint = BT
i (ξ k )σ k |Jk |Wk (4.31)
el=1 i=1 k=1

where noNodes,el and noint,el are respectively the number of nodes and integration points
inside the element el, and Wk is the weight factor (or integration weight) of the integration
point k.

In the material point method the MP carry all the information of the continuum, including the
stresses. The internal forces are computed in the following form
 24 Anura3D MPM Software, Scientific Manual

no
X elm noNodes,el
X noX MP,el

fint = BT
i (ξ MP )σ MP ΩMP (4.32)
el=1 i=1 MP=1

where noMP,el is the number of material points within element el.

4.3.4 Calculation of external forces

The external forces can be subdivided in two parts. The body force, due to the gravity accel-
eration, and the external loading due to distributed forces applied to the body.

The body force can be computed in the following way

no elm Z no elm Z
T
X X
fext,grav = ρN g dΩ = ρNT g|J| d (4.33)
el=1 Ωel el=1 

where Ωel is the volume of the element el in the global coordinate system, noelm is the number
of active elements,  is volume of the parent element, d = dξ1 dξ2 dξ3 is the infinitesimal
volume in the parent element system, |J| is the determinant of the Jacobian, and N is the
shape function matrix calculated at location ξ with respect to the parent coordinate system.
Applying the material point integration, the equation yield as follows

no
X elm noNodes,el
X noX MP,el
ext,grav
f = mMP NT
i (ξ MP )gΩMP (4.34)
el=1 i=1 MP=1

where noNodes,el and noMP,el are respectively the number of nodes and material points inside
the element el.

The external forces due to traction are computed in the following way

no elm Z no elm noNodes,el

X noX MP,el
T
X X
ext,trac
f = N t · n d∂Ω = NT (ξ MP ) fext
MP (4.35)
el=1 ∂Ωel el=1 i=1 MP=1

where fext
MP is the force stored in the material point due to distributed external forces applied at
the boundary of the body, as described in Section 4.3.2.

The complete form of the external force is summarized below:

no elm noNodes,el
X noX MP,el no elm noNodes,el
X noX MP,el

(ξ MP ) fext,trac
X X
ext T
f = N MP + mMP NT
i (ξ MP )gΩMP
el=1 i=1 MP=1 el=1 i=1 MP=1
(4.36)

4.3.5 Mass matrix

To solve the system of equations (Eq. 4.20), the mass matrix has to be inverted. In practice,
to simplify computations, a lumped mass matrix may be used instead of the consistent mass
matrix given by Eq. 4.23. The lumped mass matrix is a diagonal matrix in which each entry mi
 Explicit-Dynamic Formulation 25

is obtained by summing over the corresponding row of the consistent mass matrix. Matrix in-
versions become trivial, although the result of using a lumped mass matrix is some dissipation
Nn
P p
of kinetic energy that has been quantified by [14]. Using the property Nj = 1 it becomes
j=1

no
X elm noNodes,el
X noX MP,el
lump
M = NT (ξ MP ) mMP (4.37)
el=1 i=1 MP=1

Hereinafter, the superscript lump is removed from the lumped mass matrix Mlump to simplify
the notation, and mass matrix will always refer to a lumped matrix.

In this study the 4-node tetrahedral element is used and the lumping procedure gives the
following expression for the mass matrix M.

nelm
X
M= Mel (4.38)
el=1
 
m1 0 0 ... 0
0 m2 0 ... 0 
Mel = 
0
 (4.39)
... ... mi 0 
0 0 0 ... mNodeEl
   
mi 0 0 0 0 0
mi =  0 mi 0  ; 0 = 0 0 0 (4.40)
0 0 mi 0 0 0
noMP,el
X
mi = mMP Ni (ξ MP ) (4.41)
MP=1

4.3.6 Time discretization

The time has been discretized into different instants (k). Considering ∆t the time step size,
the time discretization approach is tk+1 = tk + ∆t.

If the general system of equations (Eq. 4.20) is posed at time tk , it can be rewritten as shown
k
in Eq. 4.42, where the acceleration â is the unknown.

k k k
Mk · â = fint + fext (4.42)

An explicit Euler time integration scheme is used to update the velocity (Eq. 4.43). This is a
first-order numerical procedure for solving ordinary differential equations (ODEs) with a given
k
initial value. Being v̂ the velocity at time tk , the velocity at the next time step tk+1 is calculated
using the acceleration at time tk as follows.
k+1 k k
v̂ = v̂ + 4t â (4.43)

On the other hand, the displacements at time tk+1 are calculated using the updated velocity
k+1
v̂ as indicated in Eq. 4.44.
k+1 k k+1
û = û + 4t v̂ (4.44)
 26 Anura3D MPM Software, Scientific Manual

4.3.7 Solution algorithm for single time step

The algorithm presented here is based on the work presented by [3], which was improved
respected the original one [2]. The key point is to work with momentum instead of velocity
as much as possible, thus avoiding divisions by nodal masses. In each time step, the MPM
computational cycle can be listed as follows.

1 The nodal mass is calculated using the shape functions and the lumped mass matrix at
time tk is formed (Eq. 4.37). The internal and external forces are evaluated in the nodes
(Eqs.4.36 and 4.32).
k
2 The momentum balance equation (Eq. 4.42) is solved and the nodal accelerations âi are
determined.
k
âi = [Mki ]–1 (fext,k
i – fint,k
i ) (4.45)

3 The velocity at the material points is updated considering Eq. 4.43 as

noX
Nodes
k
vk+1
MP = vkMP + ∆t Ni (ξ kMP ) âi (4.46)
i=1

4 Update nodal momentum

noel,i noMP,el
X X
Pk+1
i = mMP Nj (ξ kMP ) vk+1
MP (4.47)
el=1 MP=1

5 Update nodal velocities

k+1 Pk+1
i
v̂i = (4.48)
Mki

6 Compute incremental nodal displacement for solid and liquid constituent

∆ûk+1
i
k+1
= ∆t v̂i (4.49)

7 Compute strain increment

k+1
∆εk+1
MP = B(xMP )∆ûi

8 The stresses are updated using a material constitutive model

9 Update volume and density of MP

ρkMP
Ωk+1 k
MP = ΩMP (1 + ∆εvol,MP ) ρk+1
MP = (4.50)
(1 + ∆εvol,MP )

10 Particle positions are updated considering Eq. 4.44 as

noX
Nodes
k
xk+1 k
MP = xMP + Ni (ξ kMP ) ∆ûi (4.51)
i=1

11 The computational grid is initialized for the next step, nodal values are discarded and the
material points carry all the updated information.
 Explicit-Dynamic Formulation 27

4.3.8 Considerations on the use of one-phase solid analysis with porous media
Porous media are, in general, a mixture of three phases (solid, liquid and gas), which inter-
act between each other, thus determining the mechanical response of the system. Taking
rigorously into account these interactions may be in many cases unacceptably complicated,
computationally expensive, and even unnecessary for engineering applications. For example,
in numerical analyses dealing with drained and undrained conditions the presence of the pore
liquid (e.g. water) can be considered in a simplified way. In the first case, the excess pore
pressure is assumed to be zero, thus the presence of the liquid can be neglected and the solid
skeleton (e.g. soil) can be regarded as one-phase even though it is saturated. In the latter,
because of the negligible relative movement between solid and pore liquid, the equilibrium of
the solid-liquid mixture can be considered rather than the equilibrium of solid skeleton and
pore liquid as separate phases.

In undrained conditions, the stress state can be described in terms of total stresses or effective
stresses. In the second case the excess pore pressures can be computed by means of the
so-called Effective Stress Analysis, which is based on the assumption of strain compatibility
between the solid skeleton and the enclosed liquid [15]. With the Effective Stress Analyses
the pore pressure increment is computed as:

∆p = KL ∆εvol (4.52)

where ∆εvol is the volumetric strain and KL is the bulk modulus of the liquid.

Modeling of dry soil, and saturated soil in drained and undrained conditions are the field of
applicability of the one-phase formulation. For partially drained conditions a fully-coupled
two-phase formulation is necessary. This is presented in Section 4.5.
 28 Anura3D MPM Software, Scientific Manual

4.4 One-phase Single-point formulation for liquid analysis

4.4.1 Discretized momentum balance equation of liquid

The weak form can be written as follow:

Z Z
T
δ v̂TL N t · n d∂Ω – δ v̂L
T
BT σ L dΩ+
∂Ω ZΩ 
Z
T T
v̂˙ L = 0
T T
+δ v̂L N ρL g dΩ – δ v̂L N ρL N dΩ (4.53)
Ω Ω

This equation has the same form as Eq. 4.3.1. More info about the derivation of the weak
form of the momentum balance equation, the reader should refer to Sec.4.3.1.

The discretized form can be written in the following yields:

fext int
L – fL = ML â (4.54)

with

Z Z
T T
fext
L = N tL · n d∂Ω + ρL N g dΩ (4.55)
Z ∂Ω Ω

fint
L = BT σ L dΩ (4.56)
ZΩ
T
ML = ρL N NdΩ (4.57)
Ω

where fext int

L is the vector of nodal external forces, fL is the vector of nodal internal forces and
ML is the nodal mass matrix.

4.4.2 Calculation of internal forces

In finite element formulation the numerical integration of all integrals above is usually per-
formed using a reference (or parent) element system (ξ ), and the mapping procedure is done
using the Jacobian matrix. In this manual, a unit 4-node tetrahedron with linear shape func-
tions is chosen to be the reference element. Details about the calculation of the internal
forces are discussed in Section 4.3.3 and details about the mapping procedure are described
in Section 4.8.

The internal force vector becomes

no
X elm Z no
X elm Z
fint
L = BT σ L dΩ = BT σ L |J| d (4.58)
el=1 Ωel el=1 

where Ωel is the volume of the element el in the global coordinate system, noelm is the number
of active elements,  is volume of the parent element, d = dξ1 dξ2 dξ3 is the infinitesimal
volume in the parent element system, |J| is the determinant of the Jacobian, and B is the ma-
trix of shape function gradients calculated at location ξ with respect to the parent coordinate
system.
 Explicit-Dynamic Formulation 29

In the material point method the MP carry all the information of the continuum, including the
stresses. Thus, the internal forces are computed in the following form

no
X elm noNodes,el
X noX MP,el

fint
L = BT
i (ξ MP )σ L,MP ΩMP (4.59)
el=1 i=1 MP=1

where noMP,el is the number of material points within element el.

4.4.3 Calculation of external forces

As for the one-phase-solid formulation, the external forces can be subdivided in two parts:
the body force, which is generated by gravity acceleration, and the external loading due to
distributed forces applied to the boundary of the liquid body.

The body force can be computed in the following way

no elm Z no elm Z
ext,grav T
X X
fL = ρL N g dΩ = ρNT g|J| d (4.60)
el=1 Ωel el=1 

where Ωel is the volume of the element el in the global coordinate system, noelm is the number
of active elements,  is volume of the parent element, d = dξ1 dξ2 dξ3 is the infinitesimal
volume in the parent element system, |J| is the determinant of the Jacobian, and N is the
shape function matrix calculated at location ξ with respect to the parent coordinate system.
Applying the material point integration, the equation yield as follows

noelm noNodes,el
X noX MP,el
ext,grav
X
fL = mMP NT
i (ξ MP )gΩMP (4.61)
el=1 i=1 MP=1

where noNodes,el and noMP,el are respectively the number of nodes and material points inside
the element el.

The external forces due to traction are computed in the following way

no elm Z no elm noNodes,el

X noX MP,el
T
fext,trac NT (ξ MP ) fext,trac
X X
L = N tL · n d∂Ω = L,MP (4.62)
el=1 ∂Ωel el=1 i=1 MP=1

where fext
L,MP is the force stored in the material point due to distributed external forces applied
at the boundary of the body, as described in Section 4.3.2.

The complete form of the external force is summarized below:

no
X elm noNodes,el
X noX MP,el no
X elm noNodes,el
X noX MP,el

fext
L = N T
(ξ MP ) fext
MP + mMP NT
i (ξ MP )gΩMP
el=1 i=1 MP=1 el=1 i=1 MP=1
(4.63)
 30 Anura3D MPM Software, Scientific Manual

4.4.4 Mass matrix

To solve the system of equations (Eq. 4.54), the mass matrix has to be inverted. In practice,
to simplify computations, a lumped mass matrix may be used instead of the consistent mass
matrix given by Eq. 4.57. The lumped mass matrix is a diagonal matrix in which each entry
mi is obtained by summing over the corresponding row of the consistent mass matrix. Matrix
inversions become trivial, although the result of using a lumped mass matrix is some dissipa-
no
P MP
tion of kinetic energy that has been quantified by [14]. Using the property Ni (ξ MP ) = 1 it
MP=1
becomes
noelm noNodes,el
X noX MP,el
lump
X
ML = NT (ξ MP ) mMP (4.64)
el=1 i=1 MP=1

Hereinafter, the superscript lump is removed from the lumped mass matrix Mlump to simplify
the notation, and mass matrix always will refer to a lumped matrix.

In this study the 4-node tetrahedral element is used and the lumping procedure gives the
following expression for the mass matrix M.
nelm
X
ML = ML,el (4.65)
el=1
 
m1 0 0 ... 0
0 m2 0 ... 0 
ML,el = 
0
 (4.66)
... ... mi 0 
0 0 0 ... mNodeEl
   
mi 0 0 0 0 0
mi = 0
 mi 0 ; 0 = 0 0 0
  (4.67)
0 0 mi 0 0 0
noMP,el
X
mi = mMP Ni (ξ MP ) (4.68)
MP=1

4.4.4.1 Time discretization

The time has been discretized into different instants (k). Considering ∆t the time step size,
the time discretization approach is tk+1 = tk + ∆t.

If the general system of equations (Eq. 4.54) is posed at time tk , it can be rewritten as Eq. 4.42,
k
where the acceleration â is the unknown.
k k ext k
MkL · âL = fint
L + fL (4.69)
An explicit Euler time integration scheme is used to update the velocity (Eq. 4.43). This is a
first-order numerical procedure for solving ordinary differential equations (ODEs) with a given
k
initial value. Being v̂ the velocity at time tk , the velocity at the next time step tk+1 is calculated
using the acceleration at time tk as follows.
k+1 k k
v̂L = v̂L + 4t âL (4.70)
On the other hand, the displacements at time tk+1 are calculated using the updated velocity
k+1
v̂L as indicated in Eq. 4.71.
k+1 k k+1
ûL = ûL + 4t v̂L (4.71)
 Explicit-Dynamic Formulation 31

4.4.5 Solution algorithm for single time step

The algorithm presented here is based on the work presented by [3], which was improved
respected the original one [2]. The key point is to work with momentum instead of velocity
as much as possible, thus avoiding divisions by nodal masses. In each time step, the MPM
computational cycle can be listed as follows.

1 The nodal mass is calculated using the shape functions and the lumped mass matrix at
time tk is formed (Eq. 4.64). The internal and external forces are evaluated in the nodes
(Eqs. 4.63 and 4.59).
k
2 The momentum balance equation (Eq. 4.69) is solved and the nodal accelerations âL,i are
determined.
k
âL,i = [MkL,i ]–1 (fext,k int,k
L,i – fL,i ) (4.72)

3 The velocity at the material points is updated considering Eq. 4.70 as

noX
Nodes
k
vk+1 k
MP = vMP + ∆t Ni (ξ kMP ) âL,i (4.73)
i=1

4 Update nodal momentum

noel,i noMP,el
X X
Pk+1
L,i = mMP Nj (ξ kMP ) vk+1
MP (4.74)
el=1 MP=1

5 Update nodal velocities

k+1 Pk+1
L,i
v̂L,i = (4.75)
MkL,i

6 . Compute incremental nodal displacement

∆ûk+1 k+1
L,i = ∆t v̂L,i (4.76)

7 . Compute strain increment

k+1
∆εk+1
MP = B(xMP )∆ûL,i

8 The stresses are updated using a material constitutive model of water

9 Update volume and density of MP

ρkMP
Ωk+1 k
MP = ΩMP (1 + ∆εvol,MP ) ρk+1
MP = (4.77)
(1 + ∆εvol,MP )

10 Particle positions are updated considering Eq. 4.44 as

noX
Nodes
k
xk+1
MP = xkMP + Ni (ξ kMP ) ∆ûL,i (4.78)
i=1

11 The computational grid is initialized for the next step, nodal values are discarded and the
material points carry all the updated information.

4.5 Two-phase Single-point formulation for saturated soil

 32 Anura3D MPM Software, Scientific Manual

4.5.1 Discretized momentum balance equation of liquid phase

The weak form can be written as follow:
Z Z Z 
T T ρL g
δ v̂TL N tL · n
T
d∂Ω – δ v̂L T
B pL δ
T
dΩ – δ v̂L N n N dΩ (v̂L – v̂S )+
∂Ω Ω Ω k
Z Z 
T T T T
+δ v̂L N ρL g dΩ – δ v̂L N ρL N dΩ âL = 0 (4.79)
Ω Ω
 T
in which δ = 1 1 1 0 0 0 .

The discretized form can be written in the following yields:

fext int d
L – fL – f = ML â (4.80)
with
Z Z
T T
fext
L = N tL · n d∂Ω + ρL N g dΩ (4.81)
Z ∂Ω Ω

fint
L = BT σ L dΩ (4.82)
Z Ω
T ρL g
fd = N n N dΩ (v̂L – v̂S ) (4.83)
k
ZΩ
T
ML = ρL N NdΩ (4.84)
Ω
where fext int d
L is the vector of nodal external forces, fL is the vector of nodal internal forces, f is
the vector of drag forces and ML is the nodal mass matrix.

4.5.2 Discretized momentum balance equation of mixture

The weak form can be written as follow:
Z Z
T
δ v̂TS N t·n
T
d∂Ω – δ v̂S BT σ dΩ +
∂Ω Z Ω  Z 
Z
T T T T T T
+ δ v̂S N ρg dΩ – δ v̂S N n ρL N dΩ âL – δ v̂S N (1 – n) ρS N dΩ âS = 0
Ω Ω Ω
(4.85)

The discretized form can be written in the following yields:

fext – fint – M̄L âL = MS âS (4.86)
with
Z Z
T T
fext = N t · n d∂Ω + ρN g dΩ (4.87)
Z ∂Ω Ω

fint = BT σ dΩ (4.88)
ZΩ
T
M̄L = N n ρL N dΩ (4.89)
ZΩ
T
MS = (1 – n)ρS N NdΩ (4.90)
Ω
 Explicit-Dynamic Formulation 33

where fext is the vector of nodal external forces, fint is the vector of nodal internal forces, M̄L
is the nodal mass matrix (including porosity) of the liquid phase and MS is the nodal mass
matrix of the solid phase.

4.5.3 Calculation of internal forces

In finite element formulation the numerical integration of all integrals above is usually per-
formed using a reference (or parent) element system (ξ ), and the mapping procedure is done
using the Jacobian matrix. In this manual, a unit 4-node tetrahedron with linear shape func-
tions is chosen to be the reference element. Details about the calculation of the internal
forces are discussed in Section 4.3.3 and details about the mapping procedure are described
in Section 4.8.

The internal force vector becomes

no
X elm Z no
X elm Z
fint
L = T
B σ L dΩ = BT σ L |J| d (4.91)
el=1 Ωel el=1 

where Ωel is the volume of the element el in the global coordinate system, noelm is the number
of active elements,  is volume of the parent element, d = dξ1 dξ2 dξ3 is the infinitesimal
volume in the parent element system, |J| is the determinant of the Jacobian, and B is the
matrix of the shape function gradients calculated at location ξ with respect to the parent
coordinate system.

In the material point method the MP carry all the information of the continuum, including the
stresses. Thus, the internal forces are computed in the following form
no
X elm noNodes,el
X noX MP,el

fint
L = BT
i (ξ MP )σ L,MP ΩMP (4.92)
el=1 i=1 MP=1

where noMP,el is the number of material points within element el.

4.5.4 Calculation of external forces

As for the one-phase-solid formulation, the external forces can be subdivided in two parts:
the body force, which is generated by gravity acceleration, and the external loading due to
distributed forces applied to the boundary of the liquid body.

The body force can be computed in the following way

noelm Z no elm Z
ext,grav T
X X
fL = ρL N g dΩ = ρNT g|J| d (4.93)
el=1 Ωel el=1 

where Ωel is the volume of the element el in the global coordinate system, noelm is the number
of active elements,  is volume of the parent element, d = dξ1 dξ2 dξ3 is the infinitesimal
volume in the parent element system, |J| is the determinant of the Jacobian, and N is the
shape function matrix calculated at location ξ with respect to the parent coordinate system.
Applying the material point integration, the equation yield as follows
noelm noNodes,el
X noX MP,el
ext,grav
X
fL = mMP NT
i (ξ MP )gΩMP (4.94)
el=1 i=1 MP=1
 34 Anura3D MPM Software, Scientific Manual

where noNodes,el and noMP,el are respectively the number of nodes and material points inside
the element el.

The external forces due to traction are computed in the following way

no elm Z no elm noNodes,el

X noX MP,el
T
fext,trac NT (ξ MP ) fext,trac
X X
L = N tL · n d∂Ω = L,MP (4.95)
el=1 ∂Ωel el=1 i=1 MP=1

where fext
L,MP is the force stored in the material point due to distributed external forces applied
at the boundary of the body, as described in Section 4.3.2.

The complete form of the external force is summarized below:

no
X X noX
elm noNodes,el MP,el no
X X noX
elm noNodes,el MP,el

fext
L = NT (ξ MP ) fext
MP + mMP NT
i (ξ MP )gΩMP
el=1 i=1 MP=1 el=1 i=1 MP=1
(4.96)

4.5.5 Lumping procedure: Mass matrix and Drag matrix

To solve the system of equations (Eq. 4.105), the mass matrix has to be inverted. In practice,
to simplify computations, a lumped mass matrix may be used instead of the consistent mass
matrix given by Eq. 4.90. The lumped mass matrix is a diagonal matrix in which each entry
mi is obtained by summing over the corresponding row of the consistent mass matrix. Matrix
inversions become trivial, although the result of using a lumped mass matrix is some dissipa-
no
P MP
tion of kinetic energy that has been quantified by [14]. Using the property Ni (ξ MP ) = 1 it
MP=1
becomes
noelm noNodes,el
X noX MP,el
lump
X
M̄L = NT (ξ MP ) mL,MP (4.97)
el=1 i=1 MP=1

lump
Hereinafter, the superscript lump is removed from the lumped mass matrix M̄ to simplify
the notation, and mass matrix always will refer to a lumped matrix.

In this study the 4-node tetrahedral element is used and the lumping procedure gives the
following expression for the mass matrix M̄L .

nelm
X
M̄L = ML,el (4.98)
el=1
 
m1 0 0 ... 0
0 m2 0 ... 0 
M̄L,el = 
0
 (4.99)
... ... mi 0 
0 0 0 ... mNodeEl
   
mi 0 0 0 0 0
mi =  0 mi 0  ; 0 = 0 0 0 (4.100)
0 0 mi 0 0 0
noMP,el
X
mi = mL,MP Ni (ξ MP ) (4.101)
MP=1
 Explicit-Dynamic Formulation 35

The same procedure is used to calculate the lumped mass matrix for the solid phase MS .

The drag force fd can be written as follows

fd = Q (v̂L – v̂S ) (4.102)

with
Z
T ρL g
Q= N n N dΩ
Ω k

The lumping technique is also applied to the drag matrix Q as follows:

no elm noNodes,el
X noX int,el
lump
X ρL g
Q = NT
i nk |Jk | (4.103)
k
el=1 i=1 k=1

Hereinafter, likewise for the mass matrix, the superscript lump is removed from the lumped
drag matrix Qlump to simplify the notation, and drag matrix always will refer to a lumped matrix.
In this study the 4-node tetrahedral element is used and the lumping procedure gives the
following compact expression:
 
q1 0 0 ... 0
nelm
X 0 q2 0 ... 0 
Q= Qel ; Qel = 
0

... ... qi 0 
el=1
0 0 0 ... qNodeEl
   
qi 0 0 0 0 0
qi =  0 qi 0  ; 0 = 0 0 0
0 0 qi 0 0 0
noint,el
X ρL g
qi = nk Ni (ξ k )|Jk |Wk
k
k=1

where noNodes,el and noint,el are respectively the number of nodes and integration points
inside the element el, and Wk is the weight factor (or integration weight) of the integration
point k.

4.5.6 Time discretization

The discrete system of equations can be written as:
k
ML âL = fext,k
L – fint,k
L – fd,k (4.104)
k k k
MtS âS = fext,k – fint,k – M̄L âL (4.105)

Note that in Eqs. 4.105 and 4.105 the subscript L and S denote that the quantity is referred to
the fluid and water phase respectively; no subscript indicates that the quantity belongs to the
mixture.

The Euler-Cromer scheme is used to integrate the equations in time. From Eq. 4.104 the
fluid acceleration at time t is calculated and used to update the fluid velocity vt+∆t . The solid
L
acceleration is calculated solving Eq. 4.105 and used to update the solid velocity vt+ ∆t
S . Incre-
mental strains are calculated at the MP from the updated velocities, after that the constitutive
relations are used to calculate the stresses and pore water pressure.
 36 Anura3D MPM Software, Scientific Manual

In the implementation used for this manual only one set of MP, representing the solid phase
is considered. This means that the MP store all the informations regarding the soild and the
liquid phase, and their positions are updated according to the solid displacement.

4.5.7 Solution algorithm for single timestep

The initialization of MP explained in Section 4.3.2 is easily extended to the two-phase problem,
then no more details are given in this Section.

The solution sequence of a single time step is described in the following:

1 The momentum equations for the fluid and the mixture are initialized by mapping the
significative quantities from the MP to the mesh nodes. The procedure is similar to step 1
of the algorithm presented in Section 4.3.7.
t
2 Eq. 4.104 is solved for âL
t
âL = Mt,–1
 ext,t int,t d,t 
L fL – fL – f (4.106)
3 The acceleration vector ats is calculated from Eq. 4.105 as:
t
âtS = Mt,–1
 ext,t int,t
– M̄L atL

S f –f (4.107)
4 The velocities of the MPs are updated using nodal accelerations and shape functions:
nen
∆t
v̂t+ t
)âtL,i
X
t
L,MP = v̂L,MP + ∆tNi (ξMP (4.108)
i=1
nen
∆t
v̂t+ t
)âtS,i
X
t
S,MP = v̂S,MP + ∆tNi (ξMP (4.109)
i=1
∆t ∆t
5 The nodal velocities v̂t+
L and v̂t+
S are then calculated from the updated momentum
solving the following equation
nMP
t ∆t ∆t
M̄L v̂t+ )v̂t+
X
∆t
Pt+
L,i = L ≈ mL,MP NT (ξMP
t
L,MP (4.110)
MP=1
nMP
t t+∆t ∆t
)v̂t+
X
∆t
Pt+
S,i = MS v̂S ≈ mS,MP NT (ξMP
t
S,MP (4.111)
MP=1
6 Nodal velocities are integrated to get nodal incremental displacements
∆t ∆t ∆t t+∆t
∆ût+
L = ∆tv̂t+
L ; ∆ût+
S = ∆tv̂S (4.112)
7 Strains at MPs are calculated as
∆t t+∆t
∆εt+ t
L,MP = B(ξMP )v̂L ∆t (4.113)
∆t
∆εt+ ∆t
S,MP = t
B(ξMP )v̂t+
S ∆t (4.114)
and stresses are updated according to a constitutive relation
8 Water pressures at the material points are updated as:
∆t 1 h ∆t t+∆t
i
∆ε̄t+
L,MP = (1 – ntMP )∆εt+ t
S,MP + nMP ∆εL,MP (4.115)
ntMP
∆t t+∆t
pt+ t
L,MP ≈ pL,MP + KL,MP ∆ε̄L,MP (4.116)
 
where δ = 1 1 1 0 0 0 , ∆εvol,S and ∆εvol,L are the volumetric strain, i.e.
∆εvol = ∆ε11 + ∆ε22 + ∆ε33 , at the MP for the solid and liquid phase respectively.
 Explicit-Dynamic Formulation 37

9 The total stress is calculated as:

t+∆t 0t+∆t ∆t
σMP = σMP + δ pt+
L,MP (4.117)

10 Volumes associated with MP are updated using the solid volumetric strain increment
∆t t+∆t
Ωt+
MP = (1 + ∆vol,S,MP ) (4.118)

11 The positions of MP are updated using the displacements of the solid phase
nnonodes
)∆ût+1
X
∆t
xt+
MP = xtMP + t
Ni (ξMP S,i (4.119)
i=1

12 The book-keeping is updated using the new position of particles

The reader should observe that, similarly to the one-phase solution procedure, MP velocity
are calculated from nodal accelerations and nodal velocities are computed from the nodal
momentum. This is called modified Lagrangian algorithm and allows to overcome the small
mass problem.

4.6 Two-phase Single-point formulation for unsaturated soil

This section summarizes the algorithm for the solution of the two-phase Single-point MPM
formulation for unsaturated soil. Since the procedure is very similar to the two-phase Single-
point formulation for saturated soil, the reader is referred to Section 4.5 for further details.

4.6.1 Discretized momentum balance of the liquid

The discretized momentum balance equation of the liquid is

e L aL = f ext – f int – QL (vL – vS )

M (4.120)
L L

aL , vL , and vS are the nodal acceleration and velocity vectors; M

e L is the intrinsic liquid lumped
ext int
mass matrix (Eq. 4.121); fL , fL are internal and external nodal force vectors of the liquid
phase (Eqs. 4.122-4.123), and QL is the drag force matrix (Eq. 4.124).

nMP
MP
X
eL ≈
M eL N
m (4.121)
MP=1

Z nNP
X
fLext ≈ NT p
bL d∂Ωp + mMP T
L N g (4.122)
∂Ωp MP=1

nMP
X
fLint ≈ BT pMP
L mVMP (4.123)
MP=1
 38 Anura3D MPM Software, Scientific Manual

nMP
X nMP ρL g
QL ≈ NT L MP NVMP (4.124)
kL
MP=1

where p bL is the prescribed liquid pressure in the boundary ∂Ωp , N is the matrix of nodal
shape functions and B is the matrix of the gradients of the nodal shape functions evaluated at
local MP positions. The treatment of boundary conditions for unsaturated soils is described in
details in Sec. 5.3.

The phase mass of the MP is calculated as Eq. 4.125, where VMP is the volume of the MP.
MP
m
eL = ρL VMP (4.125)

4.6.2 Discretized momentum balance of the mixture

The discretized momentum balance equation of the mixture is:

MS aS + ML aL = f ext – f int (4.126)

aS and aS are the nodal acceleration vectors; MS and ML are the solid liquid lumped mass
matrices (Eq. 4.127-4.128); f ext , f int are internal and external nodal force vectors of the
mixture (Eqs. 4.129-4.130).
nMP
X
ML ≈ mMP
L N (4.127)
MP=1

nMP
X
MS ≈ mMP
S N (4.128)
MP=1

Z nMP
X
ext T
f ≈ N τ d∂Ωτ + mMP T
m N g (4.129)
∂Ωτ MP=1

nMP
X
f int ≈ BT σ MP VMP (4.130)
MP=1
Where τ is the prescribed traction vector in the boundary ∂Ωτ .

The phase mass of the MP is calculated as Eqs. 4.131-4.133, where VMP is the volume of
the MP.

mMP MP
S = nS ρS VMP (4.131)
mMP
L = nMP
L ρL VMP (4.132)
mMP
m = ρm VMP (4.133)
 Explicit-Dynamic Formulation 39

4.6.3 Time discretization

Eqs. 4.120 and 4.126 are integrated in time using the Euler-Cromer explicit method. This
time integration scheme is conditionally stable and the critical time step size is computed as
indicated in Sec. 4.7.

k
ML âL = fext,k
L – fint,k
L – fd,k (4.134)

k k k
MtS âS = fext,k – fint,k – M̄L âL (4.135)

The momentum balances (Eqs. 4.134 and 4.135) are solved at the nodes of the mesh. Mass
balances (Eqs. 3.34 and 3.37) and constitutive equations are posed locally at the MPs to
update secondary variables.

4.6.4 Solution algorithm for single time step

The numerical algorithm is based on the modified lagrangian algorithm originally proposed by
Sulsky et al. [3] for one-phase material, and successfully extended to multiphase materials
(see e.g. [11, 16, 17]). The MPM solution scheme for each time step can be summarised as
follows.

1 Liquid momentum balance equation (Eq. 4.120) is assembled and solved for the liquid
nodal acceleration aL .
2 Mixture momentum balance equation (Eq. 4.126) is assembled and solved to obtain the
nodal acceleration of the solid aS .
3 Velocities and momentum of the MPs are updated from nodal accelerations of each phase.
4 Nodal velocities are calculated from nodal momentum and used to compute the strain
increment at the MP location and the terms on the right-end-side of Eq. 3.37.
5 Mass balance equation (Eq. 3.37) and soil stress-strain equation give the increment of
pore pressure and effective stress respectively.
6 State variables at MPs are updated. Degree of saturation and hydraulic conductivity are
updated according to SWRC and HCC respectively.
7 Displacement and position of each MP are updated according to the updated velocity of
the solid phase.
8 Nodal values are discarded, the MPs carry all the updated information, and the computa-
tional grid is initialised for the next time step.

Steps 1 and 2 are commonly referred to as Lagrangian phase, while steps 3 to 8 are called
convective or Eulerian phase of the method to emphasize the fact that MPM combines the
advantages of Lagrangian and Eulerian approaches.

4.7 Stability criterion: the Critical timestep

A numerical algorithm is stable and well posed when the numerical solution is maintained
close to the exact one throughout the calculation, which means that the errors stay bounded.
The stability of a numerical algorithm depends on the integration scheme. Explicit time inte-
gration schemes, such as the one presented above, are conditionally stable.
 40 Anura3D MPM Software, Scientific Manual

The Courant-Friedrichs-Levy condition is considered to obtain the critical time interval in order
to achieve a stable solution.

4.7.1 Criterion for One-phase-solid

The criterion to determine the critical timestep in case of one-phase-solid simulation is the
following
s
L Ec
∆tcr = min ; c = (4.136)
c ρ
where Ec is the constrained modulus of the solid media. The length Lmin is the minimal length
of the element and can be determined as the shortest distance from the side of the maximum
area to the opposite node. The construction is shown in Figure 4.2.

Figure 4.2: Minimum length Lmin for a tetrahedral element.

4.7.2 Criterion for one-phase-liquid

The criterion to determine the critical timestep in case of one-phase-solid simulation is the
following
s
L KL
∆tcr = min ; cL = (4.137)
cL ρL
where KL is the bulk modulus of the water. The length Lmin is the minimal length of the
element and can be determined as the shortest distance from the side of the maximum area
to the opposite node. The construction is shown in Figure 4.2.

4.7.3 Criterion for two-phase-coupled

The stability analysis of the two-phase-coupled formulation is onerous; for this reason, a first
simplified analysis was performed in [18] on a simplified formulation that is derived by assum-
ing an incompressible pore fluid. Based on the eigenvalue analyses of a single linear finite
element, it is observed that in addition to the CFL stability condition, the influence of the per-
meability must be included. The full stability criterion for the two-phase-coupled formulation is
represented by the following equation
Explicit-Dynamic Formulation 41


∆tcr = min ∆tcr,1 ; ∆tcr,2 (4.138)

where
Lmin
∆tcr,1 = m (4.139)
c1
2ρ̃k
∆tcr,2 = (4.140)
ρL g
s
Euc
c1 = (4.141)
ρ
Euc = Ec + KL /n (4.142)
ρ = (1 – n)ρS + nρL (4.143)
1
ρ̃ = ρ + ( – 2)ρL (4.144)
n

Euc is the undrained constrained modulus of the saturated soil, Ec is the laterally confined
modulus of the solid skeleton and KL is the bulk modulus of the liquid. The length Lmin is the
minimal length of the element and can be determined as the shortest distance from the side
of the maximum area to the opposite node. The construction is shown in Figure 4.2.

It is worth noticing that ∆tcr,1 is a function of the size of the mesh, the stiffness, the porosity
and the density of the material, while ∆tcr,2 is a function of the permeability, the porosity and
the density. In many cases, the time step required to ensure the stability in low permeability
soils can be smaller than the one required for high permeability soils, i.e. ∆tcr,2 < ∆tcr,1 .

The condition expressed by Eq. 4.138 can be improved and reformulated in the single equa-
tion [19]:

q p
–2a + 4a2 + 8(b + b2 – 4c)
∆tcr = p (4.145)
b+ b2 – 4c

The coefficients a, b and c are given by

nρg
a= (4.146)
(1 – n)ρS k
4(nρKL + (1 – 2n)ρL KL + nρL Ec )
b= (4.147)
n(1 – n)ρS ρL Lmin
16Ec KL
c= (4.148)
(1 – n)ρS ρL Lmin
 42 Anura3D MPM Software, Scientific Manual

4.8 4-node tetrahedral element

The isoparametric element used in this manual is the 4-node tetrahedral. Figure 4.3 shows
the element in the parent domain, in which the node number varies from 1 to 4, and in the
global domain where the nodes numbers can vary arbitrarily depending on global mesh, and
in this case number 11, 12, 13 and 14 are used as an example. One of the main features
of the isoparametric elements is that the shape functions in both parent and global domain
are the same. Therefore, the coordinates of a point in the global domain are determined as
follows

4
X
x(ξ , t) = Ni (ξ ) xi (t) (4.149)
i=1
where ξ is the vector of the parent coordinate system which is written as follows

ξ = [ξ1 ξ2 ξ3 ]T (4.150)

Figure 4.3: 4-node tetrahedral element. Parent and Global domain.

The shape function is defined in the parent domain, i.e. in the local coordinate system, as
follows:
N1 (ξ ) = 1 – ξ1 – ξ2 – ξ3 (4.151)
N2 (ξ ) = ξ1 (4.152)
N3 (ξ ) = ξ2 (4.153)
N4 (ξ ) = ξ3 (4.154)
where the subscript 1, 2, 3 or 4 represents the node number in the parent domain.

Usually, the derivative of the shape function is also required to compute for instance the B
matrix, which is needed to determine the internal forces of the system. The derivative is
computed using the chain rule as follows
∂ Ni (ξ ) ∂ Ni (ξ ) ∂ξ 1 ∂ Ni (ξ ) ∂ξ 2 ∂ Ni (ξ ) ∂ξ 3
= + +
∂ x1 ∂ξ 1 ∂ x1 ∂ξ 2 ∂ x1 ∂ξ 3 ∂ x1
∂ Ni (ξ ) ∂ Ni (ξ ) ∂ξ 1 ∂ Ni (ξ ) ∂ξ 2 ∂ Ni (ξ ) ∂ξ 3
= + +
∂ x2 ∂ξ 1 ∂ x2 ∂ξ 2 ∂ x2 ∂ξ 3 ∂ x2
∂ Ni (ξ ) ∂ Ni (ξ ) ∂ξ 1 ∂ Ni (ξ ) ∂ξ 2 ∂ Ni (ξ ) ∂ξ 3
= + + (4.155)
∂ x3 ∂ξ 1 ∂ x3 ∂ξ 2 ∂ x3 ∂ξ 3 ∂ x3
Explicit-Dynamic Formulation 43

The compact form can be written as follows

∂ξ3  ∂ Ni (ξ) 
   
∂ξ1 ∂ξ2

∂ N i (ξ ) 
 ∂ξ1 

 ∂ x1 
  ∂ x1 ∂ x1 ∂ x1  
∂ξ3  ∂ N i (ξ )
∂ N i (ξ )  ∂ξ ∂ξ2
= 1 (4.156)
∂ x2   ∂ x2 ∂ x2 ∂ x2   ∂ξ2 
 ∂ N i (ξ )  ∂ξ1 ∂ξ2 ∂ξ3   ∂ N i (ξ ) 

  
∂ x3 ∂ x3 ∂ x3 ∂ x3 ∂ξ3

The derivatives of the natural coordinates ξ to the global coordinates x are not explicitly avail-
able. The (3x3) matrix in Eq.(4.156) is the inverse of the Jacobian matrix J, which is defined
as follows

∂ξ1 ∂ξ2 ∂ξ3 ∂ x1 ∂ x1 ∂ x1 –1

   
∂ x1 ∂ x1 ∂ x1 ∂ξ1 ∂ξ2 ∂ξ3
 ∂ξ
 1 ∂ξ2 ∂ξ3  = J–1 =
 ∂ x2 ∂ x2 ∂ x2 
 ∂ x2 ∂ x2 ∂ x2 

 ∂ξ1 ∂ξ2 ∂ξ3 
 (4.157)
∂ξ1 ∂ξ2 ∂ξ3 ∂ x3 ∂ x3 ∂ x3
∂ x3 ∂ x3 ∂ x3 ∂ξ1 ∂ξ2 ∂ξ3
44 Anura3D MPM Software, Scientific Manual
5 The boundary conditions
From a numerical point of view, two kinds of boundary conditions can be distinguished: natural
and essential. On the one hand, essential boundary conditions are imposed directly on the
solution. In this way, degrees of freedom are directly eliminated from the system of equations.
Dirichlet conditions are an example. On the other hand, natural boundary conditions are
imposed on a secondary variable, such as in Newman or Newton (Robin) conditions.

The boundary conditions required in a mechanical problem are summarized in Table 5.1.
Those are prescribed traction and prescribed velocity, which correspond to natural and essen-
tial conditions respectively. In the numerical implementation, prescribed traction is included in
the weak form of the governing equations, while calculated velocities are directly overwritten
by the prescribed values.

Table 5.1: Types of boundary conditions in a mechanical problem.

Essential Natural
Prescribed velocity Prescribed traction

It is commonly accepted that boundary conditions are specified at the nodes of the mesh
[20, 21]. This idea makes sense considering that the governing equations are solved in the
computational grid, and consequently the MPM formulation persists almost identical to that
used in the FEM.

The problem of applying the boundary conditions at the nodes appears when a moving bound-
ary is considered. In these cases, the boundary condition has to move attached to the body.
However, in the standard MPM formulation, the computational mesh remains constant. There-
fore, in such cases, the boundary condition should be carried by the material points.

This way of handling the boundary conditions to solve the issue of movement leads to some
inaccuracies, especially in natural conditions. If the boundary condition is applied on a ma-
terial point, this condition has to be distributed all over the nodes of the element to solve the
system of equations. If a material point is affected by a moving boundary condition, the ele-
ment containing the point becomes part of the contour. Then, the contour has an equivalent
thickness of the size of the element and the boundary condition is spread affecting all material
points within the cell, even those which are not at the boundary. To minimize numerical errors,
it is recommended to use a fine mesh.

In some cases a boundary condition must be applied at a deforming boundary, i.e. the compu-
tational mesh in MPM does not necessarily align with the boundary of the material. In these
cases, the nodes belonging to the boundary are determined with the procedure proposed
by [11]. For each time step, firstly the active elements (i.e. elements containing MPs) adja-
cent to an empty element (i.e. element without MPs) are detected (boundary elements); then
the nodes belonging to the element side adjacent to an empty element (boundary side) are
marked as boundary nodes. The MPs next to the boundary side are identified as boundary
material points as shown in Fig. 5.1. Finally, if the boundary node lies inside the area where a
specified condition has to be applied, e.g. in the infiltration zone or on the potential seepage
face, the corresponding boundary condition is applied.

Another difficulty with moving boundaries appears when prescribed velocity or traction (pres-
sure) do not have a constant direction, but it is normal to the boundary. This means that if the
shape of the contour changes, the direction of the applied condition has to be updated during
46 Anura3D MPM Software, Scientific Manual

𝒏𝒏

Material point (MP) Boundary

side
Boundary material point (BMP)
Empty
Active node
element
Active
Boundary node (BN)
element
Non-active node

Figure 5.1: Definition of boundary nodes and boundary material points

the calculation. The normal direction at the node is determined by means of the gradient of
mass as Eq. 5.1

PnMP
mMP
m B
n ≈ PMP=1
nMP MP
(5.1)
k MP=1 mm Bk

In Figure 5.2 a very simple problem is presented. It shows the movement of a beam that is
fixed at one end and subjected to a vertical force F at the opposite end. In this example the
fixity can be understood as an essential boundary condition in which the prescribed velocity
is zero. Because it remain motionless, it can be applied at the nodes. On the other hand, the
vertical force is a natural boundary condition which moves attached to the beam, therefore it
should be carried by the material points.

Figure 5.2: Problem of a beam subjected to a vertical force at one end and embedded
at the other one. (a) Initial discretization; (b) final discretization. The bound-
ary conditions are carried in those nodes and material points colored in red.
Yellow nodes are affected by the moving boundary condition although it is not
stored in them [from [5]].
 The boundary conditions 47

5.1 Zero kinematic boundary conditions

Zero kinematic boundary conditions are applied in a way that is identical to that of standard
Lagrangian FEM. Special attention should however be paid when applying those conditions
in the framework of MPM, that is they also have to be applied on those boundary nodes that
might become active during computations. This ensures that, when a previously inactive node
becomes active during computations, zero prescribed kinematic boundary conditions at this
node must be applied (activated). Figure 5.3 illustrates such a situation. Regarding zero
traction boundary conditions, they are automatically enforced to be satisfied by the solution of
equations of motion.

Figure 5.3: Kinematic boundary conditions are prescribed at the active and inactive
nodes (a) inactive fixity (b) active fixity. [from. Al-Kafaij 2013]

5.2 Traction boundary conditions

Zero-traction boundary conditions are automatically enforced to be satisfied by the solution
of equations of motion. Dealing with non-zero traction boundary conditions is more complex.
Traction loads can be applied on either the element boundaries or the material points (MPs).
The first option is applicable only when the boundary of the body remains aligned with loaded
element boundaries throughout the computation. The second option consists in storing the
load on selected MPs that move through the mesh.

When the traction is assigned at the element boundary, the nodal traction force is integrated
like in FEM applying Gauss quadrature and then used in the momentum balance equation
(Fig. 5.4 a). The applied load is thus integrated accurately and the traction nodal force is
non-zero only for the nodes belonging to the loaded surface. When the traction is applied
on the MPs, this is mapped to the material points located next to the element border, also
called boundary material points (Fig. 5.4 b). These MP carry surface traction throughout the
computations.

Considering, for instance, a tetrahedral element, the traction vector τe applied at the triangular
surface is interpolated from the nodes of this surface to the boundary MP. See figure 5.5.

Hence, the traction at boundary material point MP is

ntri
X
τe (xMP ) ≈ Ni (ξMP )τe (xi ) (5.2)
i=1

where Ni is the shape function of node i of the triangular surface element and ξMP are the
coordinates of the boundary material point MP inside the parent triangular element. These
48 Anura3D MPM Software, Scientific Manual

(a) (b)
Loa de d bounda ry

f tra
node
c
f tra
node
c Loa de d MP

f tra
MP
c

Node MP
Figure 5.4: (a) Traction is integrated from element boundary to adjacent nodes. (b) Trac-
tion is mapped from a loaded MP to the nodes of the element.

coordinates simply represent the projection of the MP on the triangular surface element. The
traction force vector ftrac
MP is then
ntri
Se Se X
ftrac
MP = τe (xMP ) = Ni (ξMP )τe (xi ) (5.3)
nebMP nebMP
i=1
in which nebMP denotes the number of boundary MP located next to the loaded surface. This
force vector keep the same direction and intensity throughout the calculation.

Figure 5.5: Traction is mapped from a boundary particle to all nodes of the element. [from.
Al-Kafaij 2013]

During the simulation, the boundary particles move throughout the computational mesh, and a
way to treat such tractions during the MPM simulation is to map them from boundary material
points to all nodes of the element where a boundary particle is located, see Figure 5.4b. Such
mapping procedure is carried out according to the following equation.

noel noX
X MP,el
tract
f = NT (ξMP )ftrac
MP (5.4)
el=1 MP=1
 The boundary conditions 49

The disadvantage of this way of dealing with surface tractions is that the effect of the surface
force is distributed over the layer of elements that borders the boundary. To reduce such
smearing error, the thickness of the elements along the boundary should be small.

5.3 Hydraulic boundary conditions

Typical geotechnical problems with unsaturated soils such as the study of dams, levees, and
natural slopes may need the treatment of the following BCs for the liquid phase (Figure 3.3a):

Impermeable BC: fluid is not permitted to flow across this boundary, thus the fluid velocity
perpendicular to the boundary ∂ΩvL,0 ⊂ ∂ΩvL is zero. This condition is implemented as
explained in Sec. 5.1.
Infiltration BC: imposed flux along a boundary ∂ΩvL,infilt ⊂ ∂ΩvL able to target infiltration and
evaporation processes (see Sec. 5.3.2 );
Pressure BC: prescribed pressure applied on ∂Ωp . This condition is very similar to the traction
boundary condition described in Sec. 5.2.
Hydraulic head BC: this is an alternative definition of pressure BC to simulate water levels in
b is prescribed on the boundary ∂ΩH ⊂ ∂Ωp , which is
reservoirs. A total head value H
related to the applied pressure b
pL by means of Bernoulli’s equation (Eq. 5.5) where hg
is the geometric head

bpL
H
b = hg – (5.5)
ρL g
Seepage face BC: this condition occurs at the interface between soil and the atmosphere,
like the downstream side of a river embankment, where it is not known a priori if the
condition will be Dirichelet or Neumann type. At this boundary, if soil is saturated, then
fluid is free to exit at zero pressure (i.e. b
pL = 0, Neumann BC). Otherwise, the boundary
becomes closed (i.e. prescribed velocity, Dirichelet BC). Implementation details are in
Sec. 5.3.3.

5.3.1 Total head BC

The implementation of the total head boundary condition follows [11, 22]. The total head is
converted to an imposed pressure condition pbL , using Eq. 5.5 solved for the liquid pressure.
This means that a hydrostatic approximation is adopted, neglecting local turbulence induced
by the river flow or by the potential collapse of the internal bank. A series of hydraulic head
values in time, like the ones provided by water gauge readings, is assigned with an input file,
and at each time step the current hydraulic head Hcurr is computed by linear interpolation.

The resulting nodal vector is part of the vector of external forces in the liquid momentum
equation (fext
L in Eq. 4.120) and in the mixture momentum equation (f
ext
in Eq. 4.126) to
account also for water weight.

The potential seepage face BC is automatically assigned at the boundary nodes included in
the hydraulic head BC but lying above.

5.3.2 Infiltration BC
This boundary condition is necessary to simulate rainfall or evaporation boundary condition
and has been elaborated by [11, 23]. It consists of applying a prescribed specific discharge
50 Anura3D MPM Software, Scientific Manual

q = w
b · n along the boundary, where n is the outward normal unit vector (Eq. 5.1) and the
seepage velocity w
b is defined as Eq. 5.6,

w
b = nL (vL – vS ) (5.6)

The application of this boundary condition is based on a predictor-corrector scheme: liquid

and solid velocity are first predicted assuming zero-pressure boundary conditions at the infil-
tration boundary and then (eventually) corrected to ensure the prescribed infiltration rate. The
procedure can be summarized in the following steps, also shown in Fig. 5.6:

1 atL and atS are computed by solving Eqs. 4.120 and 4.126 assuming b pL = 0 at the infiltra-
tion boundary
2 Nodal velocities are predicted as ṽLt+∆t = vtL + aL ∆t and ṽSt+∆t
= vtS + aS ∆t
3 Infiltration condition is checked. If the net infiltration discharge qnet (Eq. 5.7) is positive,
ponding conditions occur and if fluid accumulation above the boundary is not allowed (it
must remain at zero pressure) no correction is necessary. If the net infiltration discharge
is negative, or liquid ponding is allowed above the surface, then liquid and solid velocity
must be corrected to ensure the correct infiltration rate and mixture momentum balance
∆t ∆t
qnet = (nL (ṽt+
L – ṽt+
S )–w b) · n (5.7)

4 If necessary, the liquid nodal velocity is corrected by Eq. 5.8 and 5.9
∆t ∆t
vt+
L = ṽt+
L + ∆vL (5.8)

∆t ∆t
vt+
S = ṽt+
S + ∆vS (5.9)

where
∆t ∆t
(nL (ṽt+
L – ṽt+
S )–w b ) · n)n
∆vL = –   (5.10)
mL
nL 1 + m
S

mL ∆vL
∆vS = – (5.11)
mS

mL and mS are the liquid and solid nodal masses.

t+∆t – vt
5 The corrected liquid and solid acceleration is computed as atL(S) = (vL(S) L(S) )/∆t
6 The MPM solution scheme can proceed with the convective phase as explained in Sec. 4.6.4.

Note that this boundary condition is applied at nodal level, thus nodal liquid volumetric fraction
is necessary. This can be computed by mapping nL = nSL from the MPs to the nodes of the
mesh.
The boundary conditions 51

Solve momentum balance equations

assuming 𝑝𝑝�𝐿𝐿 = 0 at seepage boundary.

Predict liquid and solid velocity

no
𝑞𝑞𝑛𝑛𝑛𝑛𝑛𝑛 > 0?

yes

Is ponding
allowed?
yes no

Correct velocity and No correction

acceleration required

Go to convective phase

Figure 5.6: Flow chart of the infiltration boundary condition, from [11].
 52 Anura3D MPM Software, Scientific Manual

Absorbing boundary

spring
Kelvin-Voigt
element

dashpot

MP
Boundary node
Internal node

Figure 5.7: Absorbing boundary with Kelvin-Voigt element.

5.3.3 Potential seepage face BC

A potential seepage face can be defined as an interface between soil and atmosphere where
the fluid is free to exit at zero pressure when the soil is saturated, but it cannot enter when the
soil is partially saturated. An example of potential seepage face is the downstream surface
of a river embankment, or it can also arise after a rapid drawdown of the river level. This
boundary condition is necessary where it is unknown if the boundary is an essential or a
natural boundary condition.

The implementation can be considered as a special case of the infiltration boundary con-
dition described in Sec. 5.3.2, where w b = 0. Liquid and solid nodal velocities are pre-
dicted assuming zero-pressure at the potential seepage face (natural boundary condition).
∆t ∆t
If (nL (ṽt+
L – ṽt+
S )) · n > 0 it means that fluid is flowing out of the soil at zero pressure and
∆t ∆t
no correction is required. If (nL (ṽt+
L – ṽt+
S )) · n < 0, then fluid is flowing into the soil at zero
pressure, thus velocity must be corrected with Eqs. 5.8 and 5.9 in which w b = 0 (switch to
essential boundary condition).

5.4 Absorbing boundaries

The use of finite boundaries produces wave reflections that do not characterise a naturally
unbounded domain. Different strategies have been proposed in literature to treat wave re-
flections. One of the most common approaches in dynamic codes is the use of absorbing
boundaries. This method [24] applies viscous damping forces (dashpot) along the artificial
boundary. Al-Kafaji [10] observed that the boundary that is supported by dashpots tends to
creep continuously, particularly under sustained loading. In order to control creep, dashpots
are replaced by Kelvin-Voigt elements, which consider the parallel combination of a spring
and a dashpot as shown in Fig. 5.7.

The traction vector corresponding to the absorbing boundary has three components: a com-
ponent normal to the boundary τnab , and two components in tangential directions τtab
1
and τtab
2
.
The boundary conditions 53

The response of the Kelvin-Voigt element is described by Eq. 5.12.

τnab = –aρcp vn – kp un (5.12a)

τtab
1
= –bρcs vt1 – ks ut1 (5.12b)

τtab
2
= –bρcs vt2 – ks ut2 (5.12c)

where a and b are dimensionless parameters, vn , vt1 and vt2 the velocities, un , ut1 and ut2
the displacements, ρ the unit mass, cp and cs the velocities of the compression and shear
waves, respectively. kp and ks represent the stiffness per unit area associated to the elastic
component.

In an isotropic linear elastic material, the wave velocities (Eq. 5.13) are functions of the con-
strained modulus Ec , the elastic shear modulus G and the mass density ρ.
s s
Ec G
cp = and cs = (5.13)
ρ ρ

The first term in the right-hand side of Eq. 5.12 represents the traction given by the dashpot,
which is proportional to the velocity. The second term represents the traction given by the
spring, which is proportional to the displacement. The matrix form of the discretised equation
of motion, including the force at the absorbing boundary, reads as Eq. 5.14.

Ma = f trac – f ab + f grav – f int (5.14)

where the contribution of the absorbing boundary f ab is computed by integrating τ ab over the
surface ∂Ωab wherever the absorbing boundary is applied. It can then be written in matrix
form as Eqs. 5.15 to 5.17.

fab = Cab v + Kab u (5.15)

Z
ab
C = Nη ab NT dS (5.16)
∂Ωab

Z
ab
K = Nkab NT dS (5.17)
∂Ωab

where the dashpot matrix Cab contains the dashpot coefficients (Eq. 5.16), and the spring
matrix Kab contains the spring coefficients (Eq. 5.17).
54 Anura3D MPM Software, Scientific Manual

The coefficients kp and ks can be expressed as functions of the elastic moduli of the material
adjacent to the absorbing boundary and a virtual thickness δ as Eq. 5.18.

Ec G
kp = and ks = (5.18)
δ δ

The virtual thickness δ can be interpreted as the thickness of a virtual layer, which extends
outside the boundary. Note that for δ → 0 the absorbing boundary reduces to a rigid bound-
ary. It reduces to a dashpot boundary for δ → ∞.

Since MPs move through the mesh, the scheme is implemented in an incremental form such
that displacement and velocity increments of a MP are only accounted for when it enters the
considered element (Eqs. 5.19 and 5.20).

f ab,t = f ab,t–1 + ∆f vb,t (5.19)

∆f vb,t = C t ∆v + K t ∆u (5.20)

For the two-phase formulation, two sets of Kelvin-Voigt elements are defined. For the liquid
phase, the traction of the absorbing boundary has only the normal component and is given by
Eq. 5.21.

pab
L = –aL ρL cL vL,n – kL uL,n (5.21)

where vL,n and uL,n are the normal component of liquid velocity and displacement, respec-
tively. The velocity of the compression wave in the liquid is given by Eq. 5.22.
s
KL
cL = (5.22)
ρL

The stiffness coefficient kL can be expressed as a function of the bulk modulus of the liquid
KL and the virtual thickness δL as Eq. 5.23.

KL
kL = (5.23)
δL

Further details and an example can be found in [25].

 6 Advanced concepts

6.1 Introduction
Some additional features are described in this chapter. The contact algorithm first, then the
local damping for both 1-phase and 2-phase material, and then the mass scaling procedure for
quasi-static problems are initially discussed. Lastly, the strategies for mitigation of pathological
locking for low order elements is also examined.

6.2 Contact problems

Problems of soil-structure interaction are common in geotechnical engineering. MPM is nat-
urally capable of handling non-slip contact between different bodies, indeed inter-penetration
cannot occur because the bodies’ velocities belong to the same vector field. However, when
continuum bodies come into contact, in most cases frictional sliding occurs at the contact
surface. To simulate such a sliding interaction a specific contact algorithm that allows relative
motion at the interface between contacting bodies is required.

In this section the contact algorithm proposed by [26] is presented as well as its extension to
the adhesive contact. The advantage of this algorithm is that it does not require any special
interface element. It proved to be efficient in modeling interaction between solid bodies as
well as shearing in granular materials.

The adhesive type of contact is well suited to simulate soil-structure interaction in case of
cohesive soil under undrained conditions. Indeed, in this case the tangential force cannot
exceed the undrained shear strength.

Figure 6.1: Example of two bodies in contact.

 56 Anura3D MPM Software, Scientific Manual

Figure 6.2: Cases of approaching bodies (left) and separating bodies (right) (Al-Kafaji
2013).

6.2.1 Formulation
The contact algorithm used in Anura3Dcan be considered as a predictor-corrector scheme, in
which the velocity is predicted from the solution of each body separately and then corrected
using the velocity of the coupled bodies following the contact law.

Consider body g (gray in Fig. 6.1) and body b (black in Fig. 6.1), which are in contact at time t.
The procedure starts with the initialization of the equation of motion (Eq. 4.20) for each body
separately, as well as for the combined system.

The nodal accelerations for each body and the combined system are calculated solving the
momentum equations and then used to predict the nodal velocities at time t + ∆t as follows:
∆t
vt+
g = vtg + ∆tatg (6.1)
∆t
vt+
b = vtb + ∆tatb (6.2)
t+∆t t t
v = v + ∆ta (6.3)

In this code, the contact surface is predefined by the user during the pre-processing phase.
For each contact node, the algorithm proceeds with checking if the bodies are approaching or
separating. This is done by comparing the normal component of the single body velocity with
the normal component of the combined bodies velocity. Hence, the following two cases are
possible:
∆t
(vt+
g – vt+∆t ) · ntk > 0 ⇒ approaching
∆t
(vt+
g – vt+∆t ) · ntk < 0 ⇒ separating

where ntk is the unit outward normal to body g at node k. The algorithm allows for free
separation, i.e. no correction is required in this case and each body moves with the single
body velocity vt+ ∆t
g(b) . If the bodies are approaching, than we need to check whether sliding
occurs.

The predicted relative normal and tangential velocities at an approaching contact node can
be respectively written as:
∆t ∆t
vt+ vt+ – vt+∆t · ntk ntk


norm = g (6.4)
∆t ∆t t+∆t
vt+ t
vt+ t


tan = nk × g –v × nk (6.5)
Advanced concepts 57

Figure 6.3: Contact forces for the stick case (Al-Kafaji 2013).

These components can be used to predict the contact forces at the node as:

mtk,g
∆t ∆t
ft+
norm = vt+
norm (6.6)
∆t
mtk,g
∆t ∆t
ft+
tan = vt+
tan (6.7)
∆t
where mtk,g is the nodal mass integrated from MP of body g.

The maximum tangential force is:

∆t
fmax,t+
tan = ft+∆t t+∆t
adh + µ ftan (6.8)

where ft+ ∆t
adh is the adhesive force at the contact and µ is the friction coefficient.

Depending on the magnitude of the predicted contact forces we can distinguish between stick
and slip contact:

∆t max,t+∆t
If ft+
tan < ftan ⇒ stick contact
∆t
If ft+
tan
∆t
> fmax,t+
tan ⇒ slip contact

In the first case, i.e. the bodies stick to each others, no correction is required and the velocity
corresponds to vt+∆t . In the second case, i.e. the bodies are sliding one respect to the other,
the velocity needs to be corrected in such a way that no inter-penetration is allowed and the
magnitude of the tangential force respect Eq. 6.8.

The predicted single body velocity vt+∆t is corrected to a new velocity ṽt+∆t such that the
g k,g
normal component coincide with the normal component of the combined bodies velocity, i.e.,
∆t t
ṽt+
g · nk = vt+∆t · ntk (6.9)

which can also be written as

∆t ∆t ∆t
ṽt+
g = vt+
g + ct+
norm (6.10)

where
∆t t+∆t
ct+ – vt+∆t · ntk ntk


norm = – vg (6.11)
58 Anura3D MPM Software, Scientific Manual

is the correction of the normal component of the predicted velocity.

This correction is equivalent to apply the following normal contact force:

t+∆t mtk,g
∆t
f̃norm = ct+
norm (6.12)
∆t

When sliding occurs, the maximum tangential contact force assumes the expression:

t+∆t ∆t
f̃tan = fmax,t+
tan t (6.13)

with t being a unit vector indicating the direction of the tangent.

Substituting Eq. 6.8 into Eq. 6.13 we get:

t+∆t ∆t t+∆t
f̃tan = ft+


adh + µ fnorm t (6.14)

The total contact force is:

t+∆t t+∆t t+∆t
f̃cont = f̃tan + f̃norm (6.15)

which is used to correct the velocity as:

t+∆t
∆t ∆t f̃cont
ṽt+
g = vt+
g + ∆t (6.16)
mtk,g

It can also be written as:

∆t ∆t ∆t t+∆t
ṽt+
g = vt+
g + ct+
norm + ctan (6.17)

where the correction for the tangential component assumes the form:

∆t ∆t ∆t t+∆t
ct+ ft+


tan = adh + µ fnorm t (6.18)
mtk,g

The force introduced by the adhesion a can be expressed as;

∆t
ft+ t
adh = aAk (6.19)

where Atk is the contact area associated with the node k. It is integrated from the contact
elements that share node k.

Substituting Eqs. 6.19 and 6.12 in Eq. 6.18, the corrected velocity can be written as:

∆t ∆t ∆t
ṽt+ = vt+ vt+ – vt+∆t · ntk ntk


g g – g
aAtk ∆t
 
vgt+∆t t+∆t
· ntk


– –v µ+ t
t (6.20)
mk,g

Fig. 6.4 illustrates with a flow chart the main steps of the implemented contact algorithm.
 Advanced concepts 59

Figure 6.4: Flow chart illustrating the contact algorithm.

Having calculated the velocity of the contact node k at time t + ∆t, the corrected acceleration
vector at the node must be recalculated as:
∆t
ṽt+
g – vtg
ãtg = (6.21)
∆t
This corrected acceleration is used to update the MP velocity according to the algorithm pre-
sented previously. It should be remarked that the contact algorithm is applied between the
Lagrangian phase and the convective phase. Indeed the nodal velocities are first predicted in
the Lagrangian phase, then the corrected nodal velocities and accelerations are computed by
the contact algorithm and these new values of nodal accelerations are used to compute the
velocities at the MP in the convective phase. The same procedure explained here for body g
must be applied to body b.

6.3 Rigid body

Anura3D uses explicit analysis, and as a result, it is conditionally stable. The stability condition
(i.e., critical time step) for one-phase analysis (i.e., dry soil) follows the CFL condition [27],
which depends on the bulk modulus of the material into consideration. The larger the elastic
modulus, the lower the critical time step. Therefore, the computational cost of simulating
rigid intrutions is high. For example, the stiffnes of steel is about to 200,000 times higher
than the stiffnes of sand. Simulating rigid intrusors as rigid body save computational time
because i) there is no need for integrating stresses on rigid MPs, and ii) the critical time step
is independent of the rigid body stiffness. In this way the calculation time is minimized.

The rigid body algorithm is embedded into the contact framework and can work together with
the moving mesh algorithm. The proposed algorithm consist of the following steps which are
also illustrated in Figure 6.5 [28]:

1 In the Lagrangian phase, the nodes at the contact between the rigid body and surrounding
material are identified (Figure 6.6), the sum of internal and external forces is balanced
60 Anura3D MPM Software, Scientific Manual

with the rate of momentum, and the acceleration of the rigid body is obtained as the ratio
between the rate of moment and its total mass as follows:

 Pnc 
i
mrb g + trb – i=1 fsoil
arb = (6.22)
mrb

where arb is the acceleration vector of the rigid body, mrb is the total mass of the rigid
body, g is the gravity acceleration vector, trb is the external load vector applied to the rigid
body (equal to 0 for this study), nc is the number of nodes within the contact surface (see
Figure 6.6), i is the node relative index within the contact surface. Finally, fisoil is the force
transmitted from the soil to the FFP at the i–th contact node, which is calculated as the
i–th nodal internal force only considering the contribution from soil material points. Note
that the model analyzed in this work is axisymmetric around the FFP axis; hence, only the
vertical component of the reaction force is required.
2 At the end of the Lagrangian phase, the nodal acceleration of the rigid body (airb ) is up-
dated based on the acceleration of the rigid body calculated in the previous step (arb ):

airb = arb for: 1 ≤ i ≤ nrb (6.23)

where i denotes the node number within a total of nrb nodes (blue and red zone in Figure
6.6)
3 The nodal velocity of the rigid body (virb ) is calculated as:

virb = vi0 + airb ∆t for: 1 ≤ i ≤ nrb (6.24)

n n
with v0i being the velocity at the beginning of each time-step, arbi is the current nodal
acceleration, and ∆t is the time-step increment. Note that in this manner, the rigid motion
is ensured in all the nodes of the rigid body.
4 During the convective phase, the velocity and acceleration of the material points in the
rigid body are updated based on the nodal values. Finally, the position of the material
points is updated accordingly.

In this paper, the rigid body algorithm is used together with the contact algorithm proposed by
[26]. To account for relative displacement and free separation of different bodies, the contact
algorithm requires the evaluation of different acceleration and velocity fields at the nodes: one
for each body (vsoil , asoil , vrb , and arb ) and one for the system (vsys , asys ) (i.e., combined
body fields) [26]. In order to ensure that the contact algorithm is compatible with the rigid
body motion, the acceleration and velocity obtained in steps 2 and 3 must be prescribed to
the system (i.e., vsys = vrb and asys = arb ) before the contact algorithm is applied (Figure
6.5).

Note that the rigid body algorithm proposed here is only valid for rigid body linear motions
since we only consider the linear momentum balance. This could be extended to simulate
rigid body rotation by solving the angular momentum.
Advanced concepts 61

Figure 6.5: Rigid body algorithm flow diagram. The shaded rectangles are processes that
run independently of the rigid body algorithm. MPM phases and the contact
algorithm are included for reference.
 62 Anura3D MPM Software, Scientific Manual

Figure 6.6: Nodes belonging to the rigid intrusor (blue and red zone), and nodes belong-
ing to the contact surface (red zone).

6.4 Local damping

Natural dynamic systems contain some degree of damping of the vibration energy within the
system; otherwise, the system would oscillate indefinitely when subjected to driving forces.
Damping is due, in part, to energy loss as a result of internal friction in the intact material and
slippage along interfaces, if these are present.

6.4.1 Single-phase problem

For a dynamic analysis, the damping in the numerical simulation should reproduce in magni-
tude and form the energy losses in the natural system when subjected to a dynamic loading.
In soil and rock, natural damping is mainly hysteretic, i.e. independent of frequency, see [29]
and [30]. This type of damping is difficult to reproduce numerically. However, if a constitutive
model is found that contains an adequate representation of the hysteresis that occurs in a real
material, then no additional damping would be necessary [31].

Rayleigh damping is commonly used to provide damping that is approximately frequency-

independent over a restricted range of frequencies. However, this kind of damping introduces
body forces that hinder the steady state collapse and might influence the mode of failure [32].
[33] describes a local non-viscous damping to overcome the difficulty associated with the
viscous damping.

An alternative to Rayleigh damping is the so-called local damping. The local damping force
is proportional to the out of balance force f = fext – fint and acts opposite to the direction of
the velocity. For any degree-of-freedom in the considered system, the local damping can be
described as follows

ma = f + fdamp (6.25)
 Advanced concepts 63

where
v
fdamp = –α |f| (6.26)
|v|
is the damping force at the considered degree of freedom. The dimensionless parameter α is
called local damping factor.

6.4.2 Two-phase problem

In the two-phase formulation the fluid and the solid phase are damped separately. At any
degree of freedom the momentum equation for the water assumes the form:
grav drag damp
mw aw = ftrac
w + fw – fint
w – fw + fw (6.27)
grav
if fw = ftrac
w + fw – fint
w is the unbalanced force for the water, the damping force for the water
can be written as:
damp w
fw = –αw |fw | (6.28)
|w|

The momentum equation for the mixture is:

ms as = –m̄w aw + ftrac + fgrav – fint + fdamp (6.29)

where
damp damp
fdamp = fs + fw (6.30)

The out of balance force for the solid phase is f = ftrac + fgrav – fint , hence the damping force
for the solid can be written as:
damp v
fs = –αs |f – fw | (6.31)
|v|
In this manual the local damping factor for the fluid (αw ) and the solid (αs ) phase always
assume the same value.

Local damping was originally designed for static simulations. However, it has some character-
istics that make it attractive for dynamic simulations if proper values of the damping coefficient
are used. In quasi-static problems high value of α, i.e. 0.7-0.8, can be used to accelerate
convergence. In slow-process problems a small value of α, i.e. 0.05-0.15, can simulate natu-
ral energy dissipation of the material, if it is not taken into account by the constitutive model.
Local damping should not be used or used only with high care in highly dynamic problem,
where wave propagation is of great importance.

6.5 Convergence to static equilibrium

The dynamic MPM can also be adopted to model quasi-static problems, therefore, a criteria
for the detection of the static equilibrium is required. The static equilibrium is reached when
both the out-of-balance forces and the kinetic energy of the system vanish. The out-of-balance
convergence criterion is usually adopted in quasi-static formulation to detect equilibrium.

Let us define a dimensionless force ratio F as

kFext – Fint k
F= (6.32)
kFext k
 64 Anura3D MPM Software, Scientific Manual

Together with the out-of-balance criterion, we adopt another criterion to ensure that the kinetic
energy of the system vanishes. Let us define E, a dimensionless energy ratio, as

KE
E= (6.33)
Wext

where KE denotes the kinetic energy of the system and Wext is the work induced by the
external forces. The system is assumed to have reached static equilibrium when F and E
reach a pre-defined tolerance.

6.6 Mass scaling for quasi-static problems

The simulation of quasi-static or slow-process problems with a dynamic explicit code can
require high computational effort because the time step size is bounded by the stability con-
dition. However, if the inertia effect is negligible, then the time step size can be artificially
increased by scaling the density. Introducing the mass scaling factor β , the critical time step
√
size increases by a factor β . Indeed:

β e l
β∆t1crit
p
∆tcrit = q = (6.34)
Ec
βρ

where ∆t1crit is the critical time step for β = 1, i.e. no mass scaling is applied.

Mass scaling is a very useful technique to improve computational efficiency of dynamic codes
in simulating quasi-static ans slow-process motion. However, sensitivity analysis is necessary
to calibrate the mass scaling factor in slow process problems, indeed extremely high values
of β can significantly effect the result.

6.7 Cell-crossing instability and MPM-MIXED

The classical MPM approach suffers from a "cell-crossing instability" for problems involving
large displacements. Whenever material points cross boundaries of any element in the com-
putational background grid, an unphysical unbalance force appears at the nodes that are
shared between previous and new elements of that crossing points. The explanation for that
is the lack of smoothness of the nodal shape functions used in the interpolation of information
between grid and material points. The shape functions are linear C0 , which means that the
gradients of the shape functions within the element are constant but they are discontinuous
at the element borders (they change the sign). These gradients are used in the calculation of
the nodal internal forces as described in Section 4.3.3. Therefore, when one material point
crosses a border, the internal forces of those involved nodes suffer an unphysical instability
due to a jump discontinuity in the gradient of linear shape functions.

In this code there is a procedure to mitigate the cell-crossing instability. This simple procedure
arises from calculating the stress of each element as a constant value that corresponds to the
average of the stresses of the material points that at time tk are located within the element.
Based on that, Gauss integration is adopted to determine the internal forces (as in FEM), in
which a single point with an averaged stress σ av is considered in each element. The stress
averaging is calculated according to the following expression (Eq. 6.35), where noMP,el is
the number of material points within the element, and σ MP and ΩMP are the stress and the
 Advanced concepts 65

volume of a material point, respectively.

noP
MP,el
σ MP ΩMP
MP=1
σ av = noP
(6.35)
MP,el
ΩMP
MP=1

The calculation of the internal forces by means of the Gauss point integration is only consid-
ered in the elements located in the interior of the continuum. The Gauss point integration is
also used for the elements located at the boundaries only if the following condition is fulfilled:
noMP,el
X
ΩMP > Ffill Ωel (6.36)
MP=1

in which the factor Ffill is set to 0.9. It the aforementioned condition (Eq. (6.36)) is not fulfilled,
the internal force in the boundary elements is calculated with the classic MPM procedure.
This is a mixed approach which uses material points (MP) and Gauss points to calculate the
internal force, and in the current code is named MPM-MIXED.

6.8 Pathological locking in low order elements

6.8.1 Strain smoothing procedure

The formulation of the MPM is very close to the FEM in many aspects. At some point, this
similarity is beneficial because most of the existing knowledge can be transferred to the MPM.
However, several shortcomings associated with the FEM have been inherited by the MPM.
Probably one of the most important limitation is the kinematic locking that can occur due to the
linearity of the traditional shape functions constructed on the computational mesh. Difficulties
arise when determining the displacement field for a nearly incompressible solid.

As a simple example of this phenomenon, consider Fig. 6.7. Element e1 , in Fig. 6.7a, is
defined by nodes 1 and 2 on the x axis, and node 3 on the y axis. The area of the triangular
element must remain constant if it is incompressible. If nodes 1 and 2 are fixed, y3 must
remain constant and u3y = 0. Therefore, the remaining degree of freedom is the horizontal
displacement u3x . Similarly, the only remaining degree of freedom in the element e2 (Fig. 6.7b)
defined by nodes 4, 5, and 6, is the vertical displacement u5y . Two triangles may be assembled,
see Fig. 6.7c. Since incompressibility for element e1 requires u3y = 0 and incompressibility for
element e2 requires u3x = 0, node 3 cannot move, and both elements are completely locked
up. With nodes 1 through 4 locked up, the nodes for elements 3, 4, 5, 6, 7 and 8 will also be
locked (Fig. 6.7c). Such locking usually propagates throughout the entire mesh yielding an
unrealistic stiff response and an erroneous velocity field.

Several anti-locking approaches have been presented by different authors to mitigate this nu-
merical problem associated with the linear shape functions [34, 35]. In this work, the Nodal
Mixed Discretization (NMD) technique for linear tetrahedra elements presented by [35] has
been used to mitigate the over-stiff behaviour. This procedure has showed effectiveness in
mitigating the locking associated with (near) incompressible deformations [10]. In this tech-
nique, the element volumetric behaviour is averaged over the elements sharing its nodes via
a least squares smoothing process. The effect of applying the NMD scheme is to increase the
number of degrees of freedom per element. Being ε̇ the strain rate of an element calculated
 66 Anura3D MPM Software, Scientific Manual

Figure 6.7: A pathological case of volumetric locking in triangular elements. Degrees of

freedom are indicated.

from nodal velocities, it can be divided into deviatoric ė and volumetric ε̇vol components as

ε̇ = ė + ε̇vol δ (6.37)

where δ = (1, 1, 1, 0, 0, 0)T .

The volumetric strain rate for a node i, can be defined as weighted average of the surrounding
element values ε̇vol,e with the following equation:
no
P el,i
ε̇vol,el Ωel
el=1
ε̇vol,i = no
(6.38)
P el,i
Ωel
el=1

where noel,i are the elements surrounding the node i, and Ωel is the volume of the element el.

Then, a mean value for the element, ε̇vol is calculated by taking the average of nodal quantities
as
nonode,el
1 X
ε̇vol = ε̇vol,i (6.39)
nonode,el
i=1

where nonode,el is the number of nodes in an element.

Finally, the element strain rates is redefined by superposition of the deviatoric part and the
volumetric average as

ε̇ = ė + ε̇vol δ (6.40)

6.8.2 Pore pressure increment smoothing

In unsaturated soils, the pore pressure is calculated incrementally according to 3.37, which
can be rewritten as:

dpL e  
= KL div ρL nSL (vS – vL ) – ρL SL div(vS ) (6.41)
dt
 Advanced concepts 67

Figure 6.8: Bulk modulus variation with degree of saturation.

Where

∂ SL –1
  
∂ρL
KL = n SL
e + ρL (6.42)
∂ pL ∂ pL

is an operative bulk modulus depending on the derivative of the SWRC, bulk modulus of the
liquid, porosity, and degree of saturation. When the degree of saturation is close to 1, the oper-
ative bulk modulus increases very rapidly, see Fig. 6.8, generating large pressure oscillations.
To mitigate these oscillations, a pressure increment smoothing procedure is implemented.

With this feature, the liquid pressure increments are initially computed at the MP and then
averaged at the nodal level with a similar procedure used in Sec. 6.8. An average nodal liquid
pressure increment is computed:
 
P P
Elem MP ∆pLp ∗ Ωp
∆pN = P P (6.43)
Elem MP Ωp
P P
Where Elem runs over the elements that share node N and MP loops over the integration
points of the element. For each element, an average pressure increment is computed aver-
aging the nodal contributions using Eq. 6.44 where NEleNodes is the number of nodes of the
element.

P
NEleNodes ∆pN
∆pEl = (6.44)
NElementNodes

This average pressure increment ∆pEl is assigned to each MP inside the element.

6.9 Liquid free surface detection

The free surface problem is an important topic in coastal applications and many engineering
problems. The expression for free surface is due to the contact of two fluids having different
densities where free surface is developed at the interface. It is important to detect this surface
 68 Anura3D MPM Software, Scientific Manual

properly. In most applications where single fluid flow is involved, a fluid experiences a free
surface as exposed to the atmospheric pressure, where the air pressure corresponds to gauge
pressure. The important feature of this type of problems that the shape of the surface is
unknown, a priori, as it depends on the developed flow.

In MPM, mass density field can be evaluated at nodes each time increment. In this approach,
the density field is represented at the grid nodes using the following formula
no
P el,i noP
MP,el
Ni (ξ MP ) mMP
el=1 MP=1
ρL,i = no
(6.45)
el,i
P Ωel
nonodes,el
el=1

It is important to highlight here that the denominator of Eq. 6.45 involves only the non–empty
elements, so that the density is approximated properly.

Consequently, the density at any location x inside elements can be interpolated using
nonode,el
X
ρ̃L,MP = Ni (ξ MP ) ρL,i (6.46)
i=1

In which the ρ̃L,MP is the interpolated liquid density field. Using this procedure, the interpolated
density ρ̃L,MP is evaluated for all material points, which is used in the current MPM formulation
to capture the free surface; i.e.

ρ̃L,MP ≤ FFreeSurf ρL,0 (6.47)

where ρL,0 is the reference value of the liquid density and FFreeSurf is the factor that controls
the continuity of the free surface. The bigger this value is, but always smaller than 1, the
greater the number of particles detected. This factor is mesh dependent and needs to be set
accordingly. The value suggested for this factor is 0.7.

6.10 Moving mesh

The moving mesh procedure exploits the fact that the computational grid in MPM does not
store any permanent information. Therefore, it can be kept fixed or be freely modified at the
end of each time step. This approach is useful when applying non-zero traction and kinematic
boundary conditions because it ensures that the computational mesh aligns with the surface
where tractions (kinematics) are prescribed [36]. It is particularly convenient when simulating
the movement of a body inside the soil such as in pile driving or anchor installation where it is
necessary to keep a refined mesh always near the structure [37–41].

The positions of the grid nodes are modified according to predefined rules after having up-
dated stress, strain, velocity and position of MPs at the end of the convective phase. The MPs
are assigned to the elements and the computation proceeds with a new time step.

The discretised domain is divided as follows with the moving mesh procedure.

Moving mesh portion of the grid whose elements do not deform along the simulation but move
rigidly with certain predefined rules.
Deforming mesh portion of the grid located between the moving mesh and the fix boundary
of the domain that modifies its dimensions with time. Thus, its elements deform accord-
ingly.
Advanced concepts 69

All the elements of the discretised domain must be associated to only one of the aforemen-
tioned parts. In general, the moving mesh portion can translate in one or more directions and
rotate rigidly following predefined rules, but in Anura3D 2021it is limited to translation in one
direction.

The moving mesh usually contains a rigid body and it moves with the same displacement of
this body, thus keeping the boundary of the elements aligned to the boundary of this structure.
This is convenient when the contact algorithm is applied because it prevents the need for
identifying the contact surface and computing the normal directions at each time increment.
It is also possible to attach the moving mesh to a surface over which a traction is applied and
that moves with a prescribed velocity.

Because of the movement of the moving mesh zone, the elements of the deforming mesh
zones change their aspect ratio with time. It is important to choose a discretisation that en-
sures a reasonable aspect ratio of all the elements during the entire simulation.

Further detail and an example can be found in [42].

70 Anura3D MPM Software, Scientific Manual
References
[1] Francis H Harlow. The particle-in-cell computing method for fluid dynamics. Methods in
computational physics, 3(3):319–343, 1964.

[2] Deborah Sulsky, Z. Chen, and H.L. Schreyer. A particle method for hystory-dependent
materials. Computer Methods in Applied Mechanics and Engineering, 118(1-2):179–
196, 1994.

[3] Deborah Sulsky, Shi-Jian Zhou, and Howard L. Schreyer. Application of a particle-
in-cell method to solid mechanics. Computer Physics Communications, 87(1-2):236–
252, May 1995. ISSN 00104655. doi: 10.1016/0010-4655(94)00170-7. URL http:
//linkinghub.elsevier.com/retrieve/pii/0010465594001707.
[4] D Sulsky and HL Schreyer. Axisymmetric form of the material point method with applica-
tions to upsetting and Taylor impact problems. Computer Methods in Applied Mechan-
ics and Engineering, 0457825(96), 1996. URL http://www.sciencedirect.com/
science/article/pii/S0045782596010912.
[5] A. Yerro. MPM modelling of landslides in brittle and unsaturated soils. PhD Thesis,
Universitat Politècnica de Catalunya (UPC Barcelona), Spain, 2015.

[6] OC Zienkiewicz and T Shiomi. Dynamic behaviour of saturated porous media; the gen-
eralized Biot formulation and its numerical solution. International Journal for Numerical
and Analytical Methods in Geomechanicsnumerical and analytical methods . . . , 8:71–93,
1984. URL http://onlinelibrary.wiley.com/doi/10.1002/nag.1610080106/
abstract.
[7] O. C. Zienkiewicz, a. H. C. Chan, M. Pastor, D. K. Paul, and T. Shiomi. Static and
Dynamic Behaviour of Soils: A Rational Approach to Quantitative Solutions. I. Fully
Saturated Problems. Proceedings of the Royal Society A: Mathematical, Physical
and Engineering Sciences, 429(1877):285–309, June 1990. ISSN 1364-5021. doi:
10.1098/rspa.1990.0061. URL http://rspa.royalsocietypublishing.org/cgi/
doi/10.1098/rspa.1990.0061.
[8] Olgierd C Zienkiewicz, AHC Chan, M Pastor, BA Schrefler, and T Shiomi. Computational
geomechanics. Wiley Chichester, 1999.

[9] J Van Esch, Dieter Stolle, Issam Jassim, and J M van Esch. Finite element
method for coupled dynamic flow-deformation simulation. 2nd International Sym-
posium on . . . , (1), 2011. URL http://www.uni-stuttgart.de/igs/content/
publications/218_Issam_Dubrovnik.pdfhttp://www.researchgate.net/
publication/259579508_FINITE_ELEMENT_METHOD_FOR_COUPLED_DYNAMIC_
FLOW-DEFORMATION_SIMULATION/file/e0b4952cbd16168e83.pdf.
[10] I.K.J. Al-Kafaji. Formulation of a dynamic material point method (MPM) for geomechan-
ical problems. PhD Thesis, Institut für Geotechnik der Universität Stuttgart, Germany,
2013.

[11] Francesca Ceccato, Alba Yerro, Veronica Girardi, and Paolo Simonini. Two-phase dy-
namic mpm formulation for unsaturated soil. Computers and Geotechnics, 129:103876,
2021.

[12] D. Hillel. Soil and water physical principles and processes. Academic press, London
(UK), 1971.
72 Anura3D MPM Software, Scientific Manual

[13] Yechezkel Mualem. Hysteretical models for prediction of the hydraulic conductivity
of unsaturated porous media. Water Resources Research, 12(6):1248–1254, dec
1976. doi: 10.1029/WR012i006p01248. URL http://doi.wiley.com/10.1029/
WR012i006p01248.
[14] Jeremiah U Brackbill, Douglas B Kothe, and Hans M Ruppel. Flip: A low-dissipation,
particle-in-cell method for fluid flow. Computer Physics Communications, 48(1):25–38,
1988.

[15] P. A. Vermeer. PLAXIS 2D Reference Manual Version 5. Balkema, Rotterdam / Brook-

field, 1993.

[16] I. Jassim, D. Stolle, and P.A. Vermeer. Two-phase dynamic analysis by material point
method. International Journal for Numerical and Analytical Methods in Geomechanics,
37(15):2502–2522, 2013. doi: 10.1002/nag.2146.

[17] A. Yerro, E.E. Alonso, and N.M. Pinyol. The material point method for unsaturated soils.
Géotechnique, 65(3):201–217, 2015. doi: 10.1680/geot.14.P.163.

[18] M.M.J. Mieremet. Numerical stability for velocity-based 2-phase formulation for geotech-
nical dynamic analysis. Report 15-03, ISSN 1389-6520, Reports of the Delft Institute of
Applied Mathematics, Delft University of Technology, Delft, The Netherlands, 2015.

[19] M.M.J. Mieremet. Numerical stability for velocity-based 2-phase formulation for geotech-
nical dynamic analysis. MSc Thesis, Applied Mathematics, Delft University of Technol-
ogy, Delft, The Netherlands, 2015.

[20] E. Alonso and F. Zabala. Progressive failure of Aznalcóllar dam using the material point
method. Géotechnique, 61(9):795–808, 2011. doi: 10.1680/geot.9.P.134.

[21] E. Love and D.L. Sulsky. An unconditionally stable, energy–momentum consistent

implementation of the material-point method. Computer Methods in Applied Me-
chanics and Engineering, 195(33-36):3903–3925, July 2006. ISSN 00457825. doi:
10.1016/j.cma.2005.06.027. URL http://linkinghub.elsevier.com/retrieve/
pii/S0045782505003592.
[22] Veronica Girardi, Alba Yerro, Francesca Ceccato, and Paolo Simonini. Modelling de-
formations in water retention structures with unsaturated material point method. Pro-
ceedings of the Institution of Civil Engineers-Geotechnical Engineering, 174(5):577–592,
2021.

[23] Mario Martinelli, Wei-Lin Lee, Chjeng-Lun Shieh, and Sabatino Cuomo. Rainfall Bound-
ary Condition in a Multiphase Material Point Method. In Fifth World Landslide Forum,
pages 303–309. 2021. doi: 10.1007/978-3-030-60706-7_29.

[24] John Lysmer and R. L. Kuhlemeyer. Finite dynamic model for infinite media. In Proc. of
ASCE, pages 859–877, 1969.

[25] F. Ceccato and P. Simonini. Numerical Features used in Simulations. In E.J. Fern,
A. Rohe, K. Soga, and E.E. Alonso, editors, Mater. Point Method Geotech. Eng. A Pract.
Guid., chapter 6, pages 101–123. CRC Press, London, 2019. ISBN 1138323314.

[26] S.G. Bardenhagen, J.E. Guilkey, K.M. Roessig, J.U. Brackbill, W.M. Witzel, and J.C. Fos-
ter. An improved contact algorithm for the material point method and application to stress
propagation in granular material. Computer Modeling in Engineering and Sciences, 2:
509–522, 2001.
REFERENCES 73

[27] R. Courant, K. Friedrichs, and H. Lewy. On the partial difference equations of mathemat-
ical physics. IBM Journal of Research and Development, 11:215–234, 1967.

[28] Luis Zambrano-Cruzatty and Alba Yerro. Numerical simulation of a free fall penetrometer
deployment using the material point method. Soils and Foundations, 60(3):668–682,
June 2020. ISSN 0038-0806. doi: 10.1016/j.sandf.2020.04.002. URL http://www.
sciencedirect.com/science/article/pii/S0038080620336179.
[29] Andrew Gemant and Willis Jackson. Xciii. the measurement of internal friction in some
solid dielectric materials. The London, Edinburgh, and Dublin Philosophical Magazine
and Journal of Science, 23(157):960–983, 1937.

[30] RL Wegel and H Walther. Internal dissipation in solids for small cyclic strains. Journal of
Applied Physics, 6(4):141–157, 1935.

[31] PA Cundall. FLAC Manual: A Computer Program for Fast Lagrangian Analysis of Con-
tinua, 2001.

[32] R Hart, PA Cundall, and J Lemos. Formulation of a three-dimensional distinct element

model: Part ii. mechanical calculations for motion and interaction of a system composed
of many polyhedral blocks. In International Journal of Rock Mechanics and Mining Sci-
ences & Geomechanics Abstracts, volume 25, pages 117–125. Elsevier, 1988.

[33] PA Cundall. Distinct element models of rock and soil structure. Analytical and computa-
tional methods in engineering rock mechanics, 4:129–163, 1987.

[34] C.M. Mast, P. Mackenzie-Helnwein, P. Arduino, G.R. Miller, and W. Shin. Mitigating
kinematic locking in the material point method. Journal of Computational Physics, 231
(16):5351–5373, June 2012. ISSN 00219991. doi: 10.1016/j.jcp.2012.04.032. URL
http://linkinghub.elsevier.com/retrieve/pii/S0021999112002148.
[35] C Detournay and E Dzik. Nodal mixed discretization for tetrahedral elements. In 4th
international FLAC symposium, numerical modeling in geomechanics. Minnesota Itasca
Consulting Group, Inc. Paper, number 07-02, 2006.

[36] L. Beuth. Formulation and application of a quasi-static material point method. PhD
Thesis, Institut für Geotechnik der Universität Stuttgart, Germany, 2012.

[37] F. Ceccato, L. Beuth, P.A. Vermeer, and P. Simonini. Two-phase material point method
applied to the study of cone penetration. Computers and Geotechnics, 80:440–452.,
2016. doi: 10.1016/j.compgeo.2016.03.003.

[38] N.T.V. Phuong, A.F. Van Tol, A.S.K. Elkadi, and A. Rohe. Numerical investigation of pile
installation effects in sand using material point method. Computers and Geotechnics,
73:58–71, 2016. doi: 10.1016/j.compgeo.2015.11.012.

[39] Francesca Ceccato, Alberto Bisson, and Simonetta Cola. Large displacement numerical
study of 3D plate anchors Large displacement numerical study of 3D plate anchors.
European Journal of Environmental and Civil Engineering, 2017. ISSN 1964-8189.
doi: 10.1080/19648189.2017.1408498. URL https://doi.org/10.1080/19648189.
2017.1408498.
[40] Mario Martinelli and Vahid Galavi. Investigation of the material point method in the sim-
ulation of cone penetration tests in dry sand. Computers and Geotechnics, 130:103923,
2021.
74 Anura3D MPM Software, Scientific Manual

[41] Kaleigh M Yost, Alba Yerro, Russell A Green, Eileen Martin, and Jon Cooper. Mpm mod-
eling of cone penetrometer testing for multiple thin-layer effects in complex soil stratig-
raphy. Journal of Geotechnical and Geoenvironmental Engineering, 148(2):04021189,
2022.

[42] F. Ceccato and P. Simonini. Cone Penetration Tests. In E.J. Fern, A. Rohe, K. Soga,
and E.E. Alonso, editors, Mater. Point Method Geotech. Eng. A Pract. Guid., chapter 18,
pages 311–324. CRC Press, London, 2019. ISBN 1138323314.