Department of Civil and Environmental Engineering

The Hong Kong University of Science & Technology

CIVL 4750 Numerical Solu5ons to

Geotechnical Problems

Lecture 7: 2D FEM for Stress-Strain Analysis

Jidong Zhao

 Last & This Lecture
•  Last lecture
–  ConGnuum mechanics
–  ConsGtuGve models for soils

•  This lecture
–  RealisGc simulaGon of geotechnical problem
•  Stress-strain, force-seQlement, excavaGon etc
–  Basic procedures and issues
–  Using 2D FEM as an effecGve tool / PLAXIS

CIVL 4750 | Lecture 7 | HKUST 2

 2D Finite Element Method
•  Basic soluGon procedures by FEM
•  Mesh consideraGon
•  Element Issues
•  Convergence of soluGon
•  IniGal and in-situ condiGons
•  Displacement boundary condiGons
•  CreaGng construcGon sequences
•  Undrained/drained analysis and consolidaGon

CIVL 4750 | Lecture 7 | HKUST 3

 Revisi5ng the Spring Example
Line elasGc spring with sGffness k

d
f Basic unknowns in sGffness method of
Downwards structural analysis are nodal displacements.
Load applied
displacement of The final problem wriQen as linear
to the Spring
the Spring equaGons in terms of nodal displacements.
•  CompaGbility

Spring 1: k

1
d1 f1 •  Material behavior
Spring 2: k

2
d2 f2

CIVL 4750 | Lecture 7 | HKUST 4

 Revisi5ng the Spring Example
•  Equilibrium at each nodal point
tb

ta
Node 1 Node 2
tb f1

f2

•  Matrix form
⎡ 2k −k ⎤ ⎡ d1 ⎤ ⎡ f1 ⎤
⎢ −k k ⎥ ⎢ d ⎥ = ⎢ f ⎥
⎣ ⎦ ⎣ 2⎦ ⎣ 2⎦

Global sGffness matrix Unknown displacements Known loads

CIVL 4750 | Lecture 7 | HKUST 5

 Solu5on Procedures by FEM
1.  DefiniGon of the Problem and its Domain
2.  DiscreGzaGon of the Domain
3.  IdenGficaGon of State Variable(s)
4.  Establish Local and Global Coordinate Systems.
5.  Construct ApproximaGon FuncGons for the Elements
and Formulate Finite Element
6.  Determine ConsGtuGve EquaGons and Nature of
Two‑dimensional ApproximaGon.
7.  Obtain the Element Matrices and EquaGons
8.  Coordinate TransformaGons
9.  Assembly of the Element EquaGons.
10.  SoluGon of the Final Set of EquaGons.
11.  CalculaGon of Secondary Variables and InterpretaGon
of the Results

CIVL 4750 | Lecture 7 | HKUST 6

 1. Defining the Problem & Domain
•  Geometrical domain
•  Global coordinate axes
–  Cartesian/cylindrical/spherical coordinates
•  DefiniGon of domain: finite or infinite?
•  Interface
•  Physical domain
•  SimplificaGon of problems
–  Props, beams: rod
–  Pile cap: membrane
–  Soi seams or joints: interface
•  Soil straGficaGon
–  SimplificaGon
–  ConsGtuGve relaGons
CIVL 4750 | Lecture 7 | HKUST 7
 2. Discre5zing the Domain
•  Bounded: finite elements
•  Unbounded: infinite elements or boundary
•  Sub-domains
•  Element types
•  Refinement of elements:
•  finer meshes, beQer results, more Gme
consuming
•  higher order elements allows a coarser mesh
•  AdapGve mesh refinement
•  Successively refining the mesh in areas where
the stresses and displacements are highest.
•  AutomaGc mesh generaGon: pre-processor

CIVL 4750 | Lecture 7 | HKUST 8

 Mesh Genera5on

CIVL 4750 | Lecture 7 | HKUST 9

 Automa5c Mesh Genera5on

CIVL 4750 | Lecture 7 | HKUST 10

 Element Types

CIVL 4750 | Lecture 7 | HKUST 11

 Coupled Analysis

•  In coupled analysis of pore water pressure and stress,

choosing a shape funcGon for nodal head of one degree
lower than the shape funcGon for displacement ensures
that pore pressures are interpolated over the element
in the same manner as the effecGve stress in the soil
skeleton.
CIVL 4750 | Lecture 7 | HKUST 12
 3. Inden5fying State Variables
•  Physical phenomena:
•  heat transfer
•  fluid flow
•  solid mechanics
•  consolidaGon
•  Similar ellipGc boundary value problem
•  State variables:
•  hydraulic head/temperature
•  displacement
•  strain, stress, Darcy velocity

CIVL 4750 | Lecture 7 | HKUST 13

 Excava5on Example

Mesh

Domain

CIVL 4750 | Lecture 7 | HKUST 14

 4. Establishing local/global coordinates
•  Reasons for the use of local coordinate system:
•  allows simple interpolaGon funcGons that can
approximate values of state variables across the element
•  provides a simple method of numerical integraGon with
the element.
•  transforms arbitrary‑shaped user‑defined elements to
orthogonal coordinates to facilitate local-global
transformaGons

x
CIVL 4750 | Lecture 7 | HKUST 15
 FEM Shape Func5on
The elements can be rectangle, straight-
sided and curved side quadrilaterals.
In all cases, the local (or natural)
coordinate system is defined such that
the coordinates vary from -1 on one side
to +1 at the other, taking zero over the
quadrilateral medians.

CIVL 4750 | Lecture 7 | HKUST 16

 Iso-parametric Mapping

CIVL 4750 | Lecture 7 | HKUST 17

 Use of Shape Func5on

CIVL 4750 | Lecture 7 | HKUST 18

5. Construc5ng approxima5on func5ons for elements

•  Two enGGes approximated inside each finite

element in terms of values at the node points:
•  the physical (state) variable (the unknown: e.g. displacement,
fluid head)
•  element geometry (shape of the element).
•  Isoparametric element:
•  interpolaGon funcGon and the shape funcGon are the same
•  Approximate funcGons must be compaGble
•  to the state variables to be described
•  linear, quadraGc, cubic, or higher order
•  elements must be convex (interior angle less than 180 degree)
•  permissible elements

CIVL 4750 | Lecture 7 | HKUST 19

 6. Specifying Cons5tu5ve Rela5ons
•  A suitable consGtuGve relaGon has to be specified at
each material (Gauss) point of the element

•  For a domain with different materials, different

consGtuGve relaGon can be used for different zone
of elements

•  For geotechnical problem, consGtuGve relaGons

may involve pressure and water flow, in addiGon to
stress-strain (force-displacement) relaGon.

CIVL 4750 | Lecture 7 | HKUST 20

 7. Nature of 2D approxima5on
⎡ 1 −ν −ν ⎤
⎢ E 0 0 0⎥
E E
⎢− ν 1 −ν ⎥
⎡ εx ⎤ ⎢ 0 0 0⎥ ⎡σ x ⎤
⎢ε ⎥ ⎢ E E E ⎥ ⎢σ ⎥
⎢ y ⎥ ⎢− ν −ν 1 ⎢ y⎥
z
y
⎢ εz ⎥ 0 0 0⎥ ⎢σ z ⎥
⎢ E E E ⎥
⎢ ⎥ = ⎢ ⎥ ⎢ ⎥
γ 1 τ
⎢ xy ⎥ ⎢ 0 0 0 0 0⎥ ⎢ xy ⎥
⎢ γ yz ⎥ ⎢ G ⎥ ⎢τ yz ⎥
⎢ ⎥ ⎢ 0 1 ⎢ ⎥
x
x
⎢⎣ γ zx ⎥⎦ 0 0 0 0⎥ ⎢⎣ τ zx ⎥⎦
yz ⎢ G ⎥
⎢ 1⎥
⎢ 0 0 0 0 0 ⎥
⎣ G⎦
⎡1 − ν ν ν 0 0 0 ⎤
⎡σx ⎤ ⎢ ν 1− ν ν 0 0 0 ⎥⎥ ⎡ εx ⎤
⎢σ ⎥ ⎢ ⎢ε ⎥
⎢ ⎥ y ⎢ ν ν 1− ν 0 0 0 ⎥ ⎢ y⎥
⎢ σz ⎥ ⎢ 1 ⎥ ⎢ εz ⎥
⎢ ⎥ = A ⎢ 0 0 0
2
−ν 0 0 ⎥ ⎢ ⎥
⎢τ xy ⎥ ⎢ 1 ⎥ ⎢ γ xy ⎥
⎢τ yz ⎥ ⎢ 0 0 0 0 −ν 0 ⎥ ⎢ γ yz ⎥
⎢ ⎥ ⎢ 2 ⎥ ⎢ ⎥
τ
⎣⎢ zx ⎦⎥ ⎢ 1 ⎥ ⎣⎢ γ zx ⎦⎥
⎢⎣ 0 0 0 0 0
2
− ν⎥
⎦
CIVL 4750 | Lecture 7 | HKUST 21
 2D or 3D
•  3D ‑ modelling is recommended if:
–  Strikes of geological structures (faults, joints, etc.) do not
coincide with the long axis of excavaGon.
–  Axes of material anisotropy do not coincide with the axes
of the excavaGon.
–  DirecGons of principal stresses are neither parallel nor
perpendicular to the long axis of excavaGon.
–  The invesGgated geotechnical component (pillar,
opening, etc.) is nearly equi-dimensional, such that end
effects are not negligible.
–  MulGple openings cross each other, such that
the interacGon between adjacent openings is
significant near the interacGon.
–  AcGng forces / displacements leave components in all 3
dimensions (seismic waves).
CIVL 4750 | Lecture 7 | HKUST 22
 8. Obtaining Element Matrices & Eqs.
•  Obtain sGffness matrix
•  Transforming coordinates and assembling element
equaGons.
•  Number of element equaGons set by the number of
degrees of freedom for each state variables
•  Most based on minimizaGon of total potenGal
energy
•  Use approximaGon funcGons to the variaGonal
method to form the equaGon system for state
variable

CIVL 4750 | Lecture 7 | HKUST 23

 9. Specify Boundary Condi5ons

Pore pressure fixiGes applied

Fix displacement boundary to top surface

CIVL 4750 | Lecture 7 | HKUST 24

 10. Solving the Final Eqs Set
•  Using either ⎡ 2k −k ⎤ ⎡ d1 ⎤ ⎡ f1 ⎤
⎢ −k k ⎥ ⎢ d ⎥ = ⎢ f ⎥
–  Direct method ⎣ ⎦ ⎣ 2⎦ ⎣ 2⎦
•  Gauss EliminaGon/Cholesky’s factorisaGon (LU decomposiGon)
•  Accuracy depending on the condiGon of equaGons

–  IteraGve method
•  Gauss-Seidel or Jacobi iteraGons
•  Accuracy depending on number of iteraGons
•  Efficient, but convergence not guaranteed

–  Nonlinear method
•  Newton-Raphson iteraGon
•  For nonlinear consGtuGve relaGons (elastoplasGc)
•  Non-convergence problem
•  Improved by reducing step size of load, refining mesh and others
CIVL 4750 | Lecture 7 | HKUST 25
 Sample FE Formula5on

CIVL 4750 | Lecture 7 | HKUST 26

 Assembling the Matrix
T T T
PI = PE or ∫ B σ dV = ∫ N b dV + ∫ τ dS
N
V V S

Nodal forces corresponding Nodal forces due Nodal forces due to

to internal total stresses to body forces external tracGons

σ
b At each node, there should be
equilibrium of nodal forces

CIVL 4750 | Lecture 7 | HKUST 27

 FEM Mesh

1 3 1 5
2 2 3 4
(a) (b)

2 2
1 3 5 5 6 7
3
4 1 4
(c) (d)
•  There are only enough nodes to crudely approximate the
end nodes for a simply supported beam
•  Five nodes are required to capture the second mode as
shown in (c)
•  Seven nodes are required to capture the third mode in (d)
•  more elements? Or elements of more nodes?
CIVL 4750 | Lecture 7 | HKUST 28
 Guidelines for Meshing
•  Use as uniform a mesh paQern as possible.
•  When making transiGons from coarse‑ to fine‑mesh regions,
do not change the dimensions of adjacent elements by more
than a factor of 2. If necessary, make the transiGon over a
series of elements to maintain this factor.
•  Use quadrilaterals wherever possible. Use triangular
elements only for transiGons or when required by the
geometry.
•  The aspect raGo (length to width raGo) should be kept as
close to unity as possible. Aspect raGos of up to 5.0 are
permissible but should be kept lower than 3.0, if possible.
•  Do not use elements of extremely acute or obtuse angles.
•  Poisson's raGo < 0.5. For materials approaching Poisson's
raGo of 0.5 special elements are required to prevent locking.
CIVL 4750 | Lecture 7 | HKUST 29
 Convergence of Solu5ons
•  Equilibrium Errors
–  Source 1: yielding and the subsequent stress correcGon always gives rise to
out‑of‑balance loads.
–  Source 2: Due to sharp changes in sGffness or inappropriate Gme.
–  These errors should be <5% or preferably 1%.

•  Numerical Conversion
–  The numerical answer is not dependent on the increment size, both mesh and loading.
–  Refine the mesh and see if the answer changes.
–  Increase the number of iteraGons unGl the final soluGon doesn't change much.
–  Divide load into smaller steps.
–  Replace stress‑control loading with displacement controlled loading.
–  Higher order elements can also be used to provide a beQer soluGon.

•  Check Ques5ons
–  Does the deformaGon paQern look reasonable?
–  Do the stress paths (you need to calculate and plot) behave in an acceptable manner?
–  Are principal stress/strain direcGons consistent with applied loads?
–  Is the distribuGon of pore pressures sensible?
–  Is the distribuGon of the strains credible?

CIVL 4750 | Lecture 7 | HKUST 31

 Ini5al, In-situ and Boundary Condi5ons
•  In-situ stresses: 2D or 3D stress components and pore pressure
values at the start of analysis, such as Boundary condiGons:
boundary fixiGes or linear variaGon of stresses with depth. Use
stress converter to calculate the in-situ condiGons
•  Geometric fixiGes: nodes and edges in the mesh can be forced to
experience a zero or prescribed
displacement in x-, y-, and
z-direcGons.

•  Pore water pressure or boundary flux fixiGes: allows specifying a

pore water pressure or hydraulic gradient fixity along each of the
selected edges.
CIVL 4750 | Lecture 7 | HKUST 32
 Crea5ng Construc5on Sequences
•  Advanced construcGon sequences can be created that
include construcGon of elements, the addiGon and removal
of props (line elements) and the swapping of elements (ea.
replacement of a soil element with a concrete element).
•  Construc5ng Elements
–  Element can be added to the current increment block in one of the
two ways:
–  CreaGon
–  ConstrucGon (of an exisGng element that is currently removed)
•  Removing Elements
•  Reinstalling Elements
–  Elements may start off removed and be constructed later, or they
may start off constructed and remove later. It is also possible to
construct back in an element that has already taken part in the
analysis as a constructed element.

CIVL 4750 | Lecture 7 | HKUST 33

 Excava5on Example

Stage-1:In-situ FE Mesh Stage-2:Install tube tunnels

Stage-3:Install top slab and Stage-4: Excavate beneath top

simulate overhead road load slab to base of second slab
CIVL 4750 | Lecture 7 | HKUST 34
 Excava5on Example

Stage-5:Install second slab Stage-6: Excavate beneath second

slab to base of third slab

Stage-8: Excavate to boQom of base

Stage-7:Install third slab slab and install it
CIVL 4750 | Lecture 7 | HKUST 35
 Run the Analysis

CIVL 4750 | Lecture 7 | HKUST 36

 Undrained Analysis
•  Undrained Analysis Using Kw
–  The effecGve stress law in matrix notaGon:

–  Suppose an element of soil undergoes an incremental total stress

change Δσ which results in a change of pore pressure Δuw

–  CRISP uses this equaGon in the following way:

a)  CRISP uses effecGve stresses.
b)  In calculaGng the element sGffness matrices, CRISP considers the
volumetric sGffness of the pore water.
c)  Following the soluGon of the finite element equaGons CRISP calculates
the changes in effecGve stresses and pore pressure separately
d)  Poisson's raGo close to 0.5.
e)  Kw ≈ 100 K', say
CIVL 4750 | Lecture 7 | HKUST 37
 Drained Analysis & Consolida5on
•  Drained Analysis

•  Consolida5on Analysis

•  Total stress approach

CIVL 4750 | Lecture 7 | HKUST 38

 Summary
•  Numerical simulaGon of geotechnical
problem
•  2D Finite Element Method as an effecGve tool
•  Fundamental soluGon procedures
•  Meshing and Boundary condiGons
•  IniGal and in-situ condiGons
•  ExcavaGon process
•  Undrained/drained analysis
•  Next lecture
•  2D FDM & FLAC

CIVL 4750 | Lecture 7 | HKUST 39