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