The 12th International Conference of
International Association for Computer Methods and Advances in Geomechanics (IACMAG)
1-6 October, 2008
Goa, India
A Quasi-Static Method for Large Deformation Problems in
Geomechanics
P. A. Vermeer, L. Beuth, T. Benz
Institute of Geotechnical Engineering, University of Stuttgart, Germany
Keywords: large deformations, meshfree methods, Material Point Method, slope failure
ABSTRACT: The Finite Element Method (FEM) has become the standard tool for the analysis of a wide range of
mechanical problems. However, the classical FEM is not well suited for the treatment of large deformation
problems since excessive mesh distortions require remeshing. The Material Point Method (MPM) represents an
approach in which material points moving through a fixed finite element grid are used to simulate large
deformations. As the method makes use of moving material points, it can be classified as a meshfree method.
With no mesh distortions, it is an ideal tool for the analysis of large deformation problems. All existing MPM codes
found in literature are dynamic codes with explicit time integration and only recently implicit time integration. In
this study, a quasi-static MPM is presented. The paper starts with the description of the quasi-static governing
equations, the numerical discretisation and an explanation of calculation procedures. Afterwards, geotechnical
boundary-value problems are considered.
1 Introduction
When problems involving large deformations are modelled with an Updated Lagrangian FEM, considerable mesh
distortions occur, which require computationally time-consuming remeshing. To overcome the difficulties of the
FEM, meshfree methods have been developed, e.g. the Element-Free Galerkin Method and Smoothed Particle
Hydrodynamics (Li et al., 2002). The Material Point Method might be classified as a meshfree method or an
Arbitrary Lagrangian-Eulerian method (Więckowski, 1999).
The early beginnings of the MPM can be traced back to the work of Harlow (Harlow, 1964), who studied fluid flow
by material points moving through a fixed grid. Sulsky et al. (1996) later extended the approach for the modelling
of solid mechanic problems and called it the Material Point Method. Bardenhagen et al. (2000) extended the
method further to include frictional contact between deformable solid bodies.
The potential of the MPM for simulating granular flow was first recognised by Więckowski (1998). Several papers
on the MPM modelling of granular flow were published (Więckowski et al., 1999; Więckowski, 1998 and 2003).
Coetzee (2004) and Coetzee et al. (2005) extended the method to include a micro-polar Cosserat continuum and
described applications to anchor pull-out and excavator bucket filling.
MPM uses two discretisations of the continuum, one based on a computational mesh and the other based on a
collection of material points or “particles”. All the properties and state parameters of the continuum as well as
external loads are carried by the material points, while the grid carries no permanent information. The
computational grid is used to determine incremental displacements of material points by solving the governing
equations as with the standard FEM. Large deformations are modelled by memorizing incremental deformations
through material points that move through the mesh. By this approach, MPM combines the advantages of both
Eulerian and Lagrangian formulations.
Most MPM implementations developed so far are dynamic codes, which employ an explicit time integration
scheme. Although it is possible to use these programs also for the analysis of quasi-static problems, this is
computationally inefficient, as explicit integration requires very small time steps. For these reasons, it was
decided to develop a quasi-static MPM implementation, which uses an implicit integration scheme, thereby
broadening the possibilities of large deformation analyses for complex, large-scale geotechnical problems.
Meanwhile, the new quasi-static code has been validated (Beuth et al., 2007) by considering some boundaryvalue
problems with known analytical solutions. In fact, excellent agreement was found between results obtained
with the MPM and reference solutions. The method is now being adapted and applied to solving geotechnical
problems like slope failures or, at a later stage, pile penetrations. For the calculation of the presented problems,
high-order elements with quadratic interpolation of displacements have been employed as they reproduce stress
55
and strain fields more accurately than low-order elements which use linear interpolation. Furthermore, high-order
elements are less prone to locking effects often observed when applying low-order elements to elastoplastic
problems.
2 Quasi-static MPM formulation
In this section the field equations of quasi-static large deformations are presented. It is followed by a description
of the numerical scheme.
2.1 Field equations of quasi-static deformation
The virtual work equation obtained from the internal and external static equilibrium of a continuum reads:
∫ = ∫ ⋅ + ∫ ⋅
∂
∂
⋅
S
i i
V
i i
V j
i
ij dV u dV u dS
x
u γ δ τ δ
δ
σ (1)
where xj are Cartesian coordinates, σij denotes the Cauchy stress tensor, γi represents body forces, τi denotes
the boundary traction components and δui a kinematically admissible variation of displacements, i.e. a virtual
displacement. The symbols V and S stand for the volume of the body, the surface of the body with prescribed
traction respectively.
The development of the stress state σij is regarded as an incremental process:
ij ij ij σ = σ 0 + Δσ (2)
In this relation σij represents the actual state of stress at the end of a loading step and σij
0 represents the
previous state of stress at the beginning of this particular loading step. Equation 1 can now be written as:
∫ ∫ ∫ ∫ ∂
∂
= ⋅ + ⋅ − ⋅
∂
∂
Δ ⋅
V j
i
ij
S
i i
V
i i
V j
i
ij dV
x
u
dV u dV u dS
x
u δ
γ δ τ δ σ
δ
σ 0 (3)
So far, the final geometry of the body is taken as a reference configuration. However, at the beginning of an
increment this configuration is not known. Therefore, Equation 3 is reformulated so that the variables of the initial
geometry, i.e. 0
i x and V0, of a loading step are used as the reference configuration:
∫ ∫ ∫ ∫ +
∂
∂
= ⋅ + ⋅ − ⋅
∂
∂
Δ Σ ⋅
0 0 0 0
Higher Order Terms 0 0
0
0 0 0 0
V j
i
ij
S
i i i
V
i
j
i
V
ji dV
x
u
dV u dV u dS
x
u δ
γ δ τ δ σ
δ (4)
where ji Σ is the first Piola-Kirchoff stress tensor and Higher Order Terms express virtual work due to geometric
changes of the surface S. The increment of the first Piola-Kirchoff stress tensor can be written as:
0
0
0
0
k
k
ji
k
j
ji ji ki x
u
x
u
∂
∂Δ
+ ⋅
∂
∂Δ
Δ Σ = Δσ −σ ⋅ σ (5)
Due to the second term, this stress rate is non-symmetric. The Cauchy stress increment Δσij can be written as:
ik kj jk ki
J
ij ij Δσ = Δσ +σ 0 ⋅ Δω +σ 0 ⋅ Δω (6)
where J
ij Δσ is the Jaumann stress rate. The incremental strains ij Δε and rotations ij Δω are given by:
56
⎟ ⎟
⎠
⎞
⎜ ⎜
⎝
⎛
∂
∂Δ
+
∂
∂Δ
Δ = 2 0 0
1
i
j
j
i
ij x
u
x
u ε , ⎟
⎟
⎠
⎞
⎜ ⎜
⎝
⎛
∂
∂Δ
−
∂
∂Δ
Δ = 2 0 0
1
j
i
i
j
ij x
u
x
u
ω (7)
The constitutive relation between the Jaumann stress increment and the strain increment is:
ijkl kl ij kk
J
ij Δσ = D ⋅ Δε −σ 0 ⋅ Δε (8)
where ijkl D is the material stiffness matrix. The left hand side of Equation 4 now becomes:
[ 2 ] Right Hand Side
0 0
0
0
, ,
0
0 ∫ ⋅ Δ ⋅ +∫ ⋅ Δ ⋅ − ⋅ ⋅ Δ ⋅ =
V V
ijkl kl ij ki j k j i ki jk ij D ε δε dV σ u δu σ ε δε dV (9)
where the first term is the usual (small strain) internal virtual work and the second term contains the large
deformation contribution (McMeeking and Rice, 1975).
2.2 Numerical scheme
In its discretised form, the equilibrium equation, given by Equation 9, reads as with the classical FEM:
[ ] Higher Order Terms load nterna
+ ⋅ Δ = − +
i l
K K G v f f (10)
where the stiffness matrices are given by:
= ∫ ⋅ ⋅
0
0
V
K BT D B dV and = ∫ [ ⋅ ⋅ − ⋅ ⋅ ⋅ ]
0
0
0 2 0
V
jpk ki jqi jpk ki iqj
G
qp K L σ L B σ B dV (11)
with:
ijk ij k L N , = and [ ] ijk ijk kji B = ⋅ L + L
2
1
(12)
where Nij denote the shape function for degree of freedom k. The geometric stiffness matrix KG
qp is determined
from the stresses at the beginning of the load step. The external load vector is assembled from body forces and
surface tractions:
=∫ ⋅ +∫ ⋅
0 0
load 0 0
V S
f NT γ dV NT τ dS (13)
while the so-called internal force vector is given by:
0
0
nternal
0
f B dV
V
T
i
= ∫ ⋅σ (14)
In standard FEM Gauss-Legendre integration is used to integrate the stiffness terms over the element volume
with a fixed number of integration points. With the MPM, material points can move through a calculation grid. The
number of material points in a specific element can therefore vary so that integration is performed with a changing
number of material points over an element. Details of MPM integration will be presented in an upcoming paper.
57
3 MPM calculation process
The calculation process of the quasi-static MPM can be divided into three steps: Initialisation Phase, Lagrangian
Phase and Convective Phase (Sulsky, 1996) as explained in the following with reference to the equations given in
the previous section.
Figure 1. Initialisation Phase.
3.1 The Initialisation Phase
This phase starts by generating a finite element mesh. In contrast to the FEM, the mesh is not only generated
where material exists, but over the complete domain where material is expected to move. Material points are
placed inside elements to form (define) the solid body. This is shown in Figure 1 for the simple case of a
cantilever. Elements containing material points are called activated elements and elements containing no material
points are called deactivated elements. As material points move through the grid, elements become activated and
deactivated, the state of an element being determined using a housekeeping algorithm. Similar to the FEM,
constraints are handled by fixing nodes, i.e. zero displacements are enforced at these nodes (Figure 1). During
the Initialisation Phase, loads are assigned to material points, which they carry throughout the deformation
process. The material points are initially evenly distributed inside an element, and thus each given the same
volume Ωp.
3.2 The Lagrangian Phase
With all state variables initialised at material points during the Initialisation Phase, the loading steps can be
executed, during which external loads are incrementally applied. During each load step, the information carried by
the material points is projected onto a background finite element mesh where equations of equilibrium are solved
in an Updated Lagrangian frame as with the FEM. The result of the Lagrangian Phase are nodal displacement
increments Δv as shown in Figure 2(b).
3.3 The Convective Phase
After the Lagrangian Phase, the mesh is reset to its initial configuration as shown in Figure 2(c), while the
material points keep their global positions. The housekeeping algorithm now determines which elements should
be activated and deactivated.
The material point volumes Ωp are updated using:
i
p
p
i
i p
p
i
p J
J
J
Ω = ⋅ ⋅
−
−
0
1
ω 1 (15)
where Jp is the Jacobian matrix evaluated at the material point positions and ω
p
i-1 represents the material point
volumes in the element’s reference coordinate system at the end of the last loading step.
The mass of each material point, used for calculating the body forces, is calculated by:
Material points
Activated elements Deactivated elements
Fixed nodes
Applied force
58
p p m = ρ ⋅ Ω (16)
where ρ is the material density assigned to the material point.
Figure 2. The 3 phases of the MPM calculation process.
4 Benchmark problems
4.1 Slope problem
The first boundary-value problem, as shown in Figure 3, consists of a slope with an inclination of 45° and a height
of 1 m as might be used for a model test. A free-slip boundary condition is applied at the vertical cut through the
slope; at the lower mesh boundary no slip is allowed. The weight of the material is increased under gravity
loading up to 96 kN/m³. Such loading is typical for experiments on model slopes in a geocentrifuge. The problem
was analysed using an elastic-plastic model with Mohr-Coulomb failure envelope, a cohesion c = 1 kPa, friction
angle ϕ = 25°, dilatancy angle ψ = 0°, Poisson’s ratio ν = 0.33 and a Young’s modulus of 50 kPa. As the currently
developed MPM code is a full 3D code, the plane strain problem was analysed in a 3D slice as shown
Figure 3. Slope geometry (left) and initial configuration showing the activated elements (right).
c = 1 kPa
ϕ = 25°
ψ = 0°
E = 50 kPa
ν = 0.333
(a) Initialisation Phase (b) Lagrangian
Phase
(c) Convective Phase
59
Figure 4. Computed load-displacement curves for Point A (left) and deformed mesh at the end of the FEM
calculation (right).
in the right hand side of Figure 3. The mesh with the active and inactive elements at the inset of loading is shown
on the right hand side of Figure 3. 15-noded wedge elements with initially 18 material points per element, each
with an initial weight of ω
p = 1/18, have been used for this calculation.
Figure 4 shows the unit weight of the soil γ plotted as a function of the total displacement for the point located at
the top of the slope, marked as “A” in Figure 3. The load-displacement curve calculated with the MPM is
compared with the result from an Updated Lagrangian FEM simulation. When reaching the range of large
deformations, an increasing stiffening effect due to the changing geometry of the slope can be observed from
both FEM and MPM. However, the MPM calculation is able to proceed much further than the FEM. In fact, the
FEM gives good results up to a displacement of 25 cm, but then it produces very poor results. Figure 4 shows on
the right hand side the deformed mesh at the very end of the FEM calculation. Obviously, large mesh distortions
occur around the toe point, meaning that FEM calculation results cannot be considered reliable beyond a load of
γ ≈ 25 kN/m³; at least not without remeshing.
Figure 5. Deformation of the slope and incremental shears.
Although FEM results may be somewhat improved by remeshing techniques, this would not solve the “flow”
problem at the toe of the slope. A “hard” horizontal boundary in front of the slope is to limit vertical movement of
the slope material as indicated by the dashed line in Figure 4. Preventing material to deform beyond this
boundary in case of the FEM would bring along the additional complexity of a contact algorithm.
Figure 5 shows the computed deformation of the slope represented by particles at different loading stages. As
can be seen, material is deposited within newly activated elements in front of the original slope toe with the MPM,
giving additional support to the slope due to contact with the lower mesh boundary. Compared to the FEM, the
additional computational costs of the MPM are small.
Figure 5 also shows the development of shear bands within the soil. Up to γ = 30 kN/m³, deformation occurs in a
single shear band. However, starting at γ = 30 kN/m³, a new shallower shear band starts to develop. Hence, from
30 to 40 kN/m³, two shear bands are active. At γ = 80 kN/m³, the first one is no longer active but a third shear
band is initiated.
4.2 Retaining wall problem
The retaining wall problem consists of a block being pushed into the soil and pulled back. Two different materials
FEM
γ = 30 kN/m³ γ = 40 kN/m³ γ = 80 kN/m³
60
are involved in this problem: the soil is modelled with an elastic-plastic Mohr-Coulomb constitutive law and the
block is composed of a linear elastic material. Figure 6 shows the geometry and discretisation of the problem.
Like the slope from the previous section, it is a plane-strain problem with a width of 0.1 m perpendicular to the
plane of deformation. The total domain has dimensions 5.3 m x 1.2 m. Full fixities are applied along the lower
edge, and a free-slip condition is applied on the left and right sides. The domain is filled with 15-noded wedge
elements, each element containing 18 particles initially. The block, lightly shaded in Figure 6, has dimensions 0.6
x 0.4 m, and possesses a material weight of γ = 20 kN/m³, a Young’s modulus E = 100 MPa and a Poisson’s ratio
ν = 0.0. The remaining active elements are soil elements with a material weight γ = 20 kN/m³, Young’s modulus E
= 10 MPa and Poisson’s ration = 0.3, cohesion c = 10 kN/m², a friction angle of ϕ = 30° and a dilatancy angle ψ =
0°. The calculation is performed in two phases: firstly, gravity loading is applied, secondly, prescribed block
displacements of 600 mm are applied on the soil-structure contact surface. In the inverse case of an active earth
pressure problem, the block is displaced backwards away from the soil by an (active) displacement of 400 mm.
In this problem, block elements are updated as in conventional Updated Lagrangian FEM, such that the contact
boundary between soil and structure coincides with element boundaries at all stages of the analysis. For the soil
itself, however, the meshfree approach is maintained, but not in the sense that the mesh is fixed. Instead of fixing
the soil mesh at its initial position, the soil mesh moves with the (prescribed) displacement of the block structure
mesh. Hence, instead of resetting the mesh, we set the mesh to a new position.
Figure 7 shows the computational results for the passive earth pressure problem. A completely rough wall is
considered with full adhesion, so that the soil just in front is sticking to the block, but a passive wedge further in
front is pressed out. Some of the soil is falling on top of the block. On the right hand side of Figure 7, shear bands
in front of the block structure as well as the heavy shear deformations occurring at the bottom of the wall can be
observed. A third shear band develops. On the left hand side of Figure 7, the computed load-displacement curve
is displayed. The force consists of the horizontal earth pressure in front of the wall and the horizontal shear force
below the block (for a wall length of 0.1 m). Large displacements cause soil heave in front of the block and
consequently an increase of the earth pressure. For small deformations the force amounts to 5.6 kN, which
complies well with the analytical solution from classical earth pressure theory and the sliding force below the
block.
Figure 6. Geometry of the retaining wall problem (left) and initial configuration of activated elements within the
finite element grid (right).
Figure 7. Passive case: computed load-displacement curve (left) and incremental shear strains inside the
deformed soil after a wall movement of 600 mm (right).
c = 10 kPa
ϕ = 30°
ψ = 0°
E = 10 MPa
ν = 0.3
γ = 20 kN/m³
E = 100 MPa
ν = 0.0
γ = 20 kN/m³
initial location of the block
Initially active soil elements
61
Figure 8. Active case: computed load-displacement curve (left) and shear strains inside the deformed soil after a
wall movement of 400 mm (right).
Figure 8 shows computational results for the active wall problem. The block has lost its contact to the soil, as a
free slope has been formed. A few soil particles are still sticking to the wall, as a cohesive soil with full adhesion
to the wall is considered. On the other side of the block some soil is creeping upwards. An extremely small area
of shear strains can be seen indicating a slope failure zone and a strong horizontal shear band below the wall. On
the left hand side, the computed load-displacement curve is shown, where Fx is a pulling force on the wall due to
the fact that a very cohesive soil with full adhesion at the soil-wall contact is considered. Moreover, all
calculations were done allowing tension. For small deformations, a peak load of about 1.5 kN is found, being the
sum of the active earth pressure and the tangential shear force at the bottom of the block. On continuing
displacement, geometric softening reduces the force to a residual force of 1.0 kN being the shearing resistance of
the block along the bottom.
No doubt, such problems are best modelled by introducing interface elements at the soil-structure contact. As yet
this has not been tried out, but there would seem to be no obstacle for doing so.
5 Conclusions
The concept of the quasi-static Material Point Method and results from a slope and retaining wall problem have
been presented to illustrate the possibilities of the MPM, which arise from the fact that large deformations can be
modelled without the high computational costs of remeshing.
It is shown that the MPM can model the failure of a slope. The calculated load-displacement curves correspond
well with the results from the FEM as long as deformations remain relatively small. Thereafter, the MPM predicts
well the geometric stiffening of the deforming slope taking into consideration the deposition of material in front of
the slope toe, while the FEM encounters problems due to heavy distortions of elements.
The presented active and passive earth pressure problems show, how soil-structure interaction can be
considered with the MPM by modification of the finite element mesh. The presented approach poses no problem
to the MPM, as all permanent information of the continuum is stored with the material points.
In the current MPM, 15-noded wedge elements with quadratic interpolation functions are used, as this yields
more accurate results than with 6-noded elements using linear interpolation of the displacement field. In this
paper, the material was assumed to be elastic-plastic, but in the future, hardening plasticity will be added.
Interface elements are currently being developed in order to improve the simulation of large deformations in
problems that involve two materials, structural elements respectively.
It is concluded that the MPM is well suited for the modelling of large deformation problems. Geotechnical
applications such as pile driving, installation of bucket foundations will be investigated in the near future.
6 Acknowledgements
The authors wish to thank Professor Więckowski from Technical University of Łódź, Poland, and Dr. Bonnier from
Plaxis B.V. in Delft, The Netherlands, for their valuable advise on the development of the quasi-static MPM. The
development of the MPM code is sponsored by Plaxis B.V. and Deltares, both situated in Delft, The Netherlands.
It should also be stated, that the Plaxis 3D Foundation FEM code was used as a starting point for the
development of the MPM.
initial location of the block
62
7 References
Bardenhagen S.G., Brackbill J.U., Sulsky D. 2000. The material-point method for granular materials. Computational Methods in
Applied Mechanics and Engineering, 187, 529-541.
Beuth L., Benz T., Vermeer P.A., Coetzee C.J., Bonnier P., Van Den Berg P. 2007. Formulation and Validation of a Quasi-
Static Material Point Method, 10th International Symposium on Numerical Methods in Geomechanics (NUMOG), Rhodes,
Greece.
Beuth L., Benz T., Vermeer P.A. (2007/08). Large deformation analysis using a quasi-static Material Point Method. Journal of
Theoretical and Applied Mechanics, Bulgarian Academy of Sciences, in print.
Coetzee C.J. 2004. The modelling of granular flow using the particle-in-cell method. PhD Thesis, Department of Mechanical
Engineering, University of Stellenbosch, South Africa.
Coetzee C.J., Vermeer P.A., Basson A.H. 2005. The modelling of anchors using the material point method. International
Journal for Numerical and Analytical Methods in Geomechanics, 29, 879-895.
Harlow F.H. 1964. The particle-in-cell computing method in fluid dynamics. Methods in Computational Physics, 3, 319-343.
Li S., Liu W.K. 2002. Meshfree and particle methods and their applications. Appl. Mech. Rev., 55, No. 1.
McMeeking R.M., Rice J.R. 1975. Finite-element formulations for problems of large elastc-plastic deformation. Int. J. Solids
Struct., 11, 601-616.
Sulsky D., Schreyer H.L. 1996. Axisymmetric form of the material point method with applications to upsetting and Taylor impact
problems. Comput. Methods Appl. Mech. Engng., 139, 409-429.
Więckowski Z. 1998. A particle-in-cell method in analysis of motion of a granular material in a silo. Computational Mechanics:
New Trends and Applications, CIMNE, Barcelona.
Więckowski Z., Youn S.K., Yeon Y.H. 1999. A particle-in-cell solution to the silo discharging problem. Int. J. Numer. Meth.
Engng., 45, 1203-1225.
Więckowski Z. 2003. Modelling of silo discharge and filling problems by the material point method. Task Quarterly, 4, 701-721.
63
A Quasi-Static Method for Large Deformation Problems in Geomechanics
Subscribe to:
Post Comments (Atom)
0 comments:
Post a Comment