CY Ou

Finite Element Analysis of Deep

Excavation

Chang-Yu Ou, Chair Professor

Department of Civil and Construction Engineering
National Taiwan University of Science and Technology

08:01:34
 CY Ou
1. FEA Theory

Linear elastic material

From virtual work principle, we obtain the relationship

between the internal strain energy and external work


vol

{}T [ ] d (vol) = {d }T {F }d (vol ) +
vol

area
{d }T [T ] d (area)

By simplification, we obtain the relationship between internal

force and external force


vol

[ B]T [ ] d (vol ) = [ N ]T {F } d (vol ) +
vol

area
[ N ]T [T ] d (area)

08:01:36
 CY Ou

We therefore establish the element stiffness, global nodal

action and displacement

[K ] =  [ B] [ D][ B]dV
T

vol

[ K ]{d} = {R}

where [D] can be expressed in terms of Young’s modulus

and Poisson’s ratio

08:01:36
 CY Ou
Nonlinear and Plastic material

According to plastic theory, we can derive the following relationship

( )
{} = [ D] − [ D p ] {} = [ Dep ]{}

[ Dep ] = [ De ] − [ D p ]

[D] is the elastic constitutive law (E and m)

[Dp] is the plastic constitutive law (soil model)
[Dep] is the elastoplastic constitutive law.

[K ] =  [ B] [ D
T ep
][ B]dV
vol

08:01:36
 CY Ou
Iteration scheme (e.g. Newton -Raphson Method)

Load
T
2.[ K ] =  [ B] [ D ep [ B]d v

8. update stress { k +1} = { k } + { } and

10. = {Fext } − {Fint } hardening parameter of the yield function
1.{Fext }

T
9.{Fint } =  [ B] { }dv

iter Displacement
12.{d } = {d 0 } +  {d k }
k =1

3.{d 0 } = [ K k =1 ]−1{Fext } 11.{d k } = [ K k ]−1 k Stress integration scheme(e.g. Forward Euler method with subincrements)

[Dep]
σ

4.{ } = [ B]{d } { }

Δσ 5.{ } =
m

6.{ i } = [ D ep ]{ i }

m
7.{ } =  { i }
i =1

ε
Δε
m increment
Δε/m

08:01:36
 CY Ou
Detailed computation procedure (by KH Yang)

08:01:36
 CY Ou
The FEM computation procedure is as follows
(1) Divide the external load into many small increments, {Fext}
(2) Establish the stiffness matrix for each element, [K]
(3) Compute the displacement increment, {d0}, using the standard FEM procedure
(4) Compute the strain increment, {}
(5) Divide the strain increment into m sub-strain increment, {}
(6) Compute the sub-stress increment, {} using Forward Euler method
(7) Add the sub-stress increment together and obtain the stress increment, {}
(8) Update the stress increment and hardening parameters of the yield functions
(9) Obtain the equivalent internal force, {Fint}
(10) Obtain error, ={Fext}-{Fint }(check convergence)
(11) Use  as external force and calculate residual displacement increment, {dk}
(12) Update the displacement increment by adding , {dk}, together

08:01:36
 CY Ou

Simplified computation procedure (by KH Yang)

U   F
internal

displacement strain stress internal force

Finternal-Fexternal ≤ tolerance
Check convergence

08:01:36
 CY Ou
Check convergence

Convergence criteria

Fint − Fext
≤1tolerance
%
Fext

Percentage error 1% by default

3%~5% is more practical

Check convergence.
If not converged, change(increase) displacement u

08:01:36
 CY Ou
Detailed computation procedure (by KH Yang)
Global/Nodal
(Force vs. Displacement) START at
a Given Boundary Conditions
e.g. Prescribed force or displacement

If Yes
If NOT CACULATE Next loading increment
External Force Increment Fext

Iteration Scheme
e.g. Newton-Raphson CACULATE
Displacement Increment u
CHECK
Convergence
CONVERT to Fext –Fint <tolerance
Compatibility Strain Increment 

Constitutive Model
INTEGRATE
Stress Integration Algorithm
Stress Increment 
e.g. forward Euler scheme Local/Gaussian
(Stress vs. Strain)

CACULATE
Internal Force Increment Fint
08:01:36
 CY Ou
2. Types of analysis

Drained analysis: No excess porewater pressure generation

Undrained analysis: With excess porewater pressure generation

Effective stress analysis: soil and water are analyzed

separately
Total stress analysis: soil and water was analyzed as a
single material

Effective stress undrained analysis: c’, f’, g’, with water table

Total stress undrained analysis: c, f=0, gsat, no water table

08:01:36
 CY Ou
Effecive stress undrained analysis
{ } = { } + { w }
{ } = [ D]{ } { w} = [ Dw ]{ }
{ } = [ D]{ } = [ D]{ } + [ Dw ]{ } [ D] = [ D] + [ Dw ]
FE computation:

1. Input soil effective stress parameters（c’,f’, E’,m’）

2. Input water stiffness [Dw]
3. Compute [ D] = [ D] + [ Dw ]
4. Establish [K]
5. Use FEM procedure,[ K ][d ] = [ R] , obtain nodal displacement [d ]
6. Compute[ ] = [ B][d ], obtain strain{}
7. Use { } = [ D]{ } , to obtain{ }
8. Use { w} = [ Dw ]{ }，obtain{ w}

2018/10/31 12
 CY Ou
Total stress undrained analysis

{s}=[D]{e}

FE computation:

1. Input soil total stress parameters（c, f, E, m）

2. Compute [D]
3. Establish [K]
4. Use FEM procedure, [ K ][d ] = [ R] , obtain nodal displacement [d ]
5. Compute [ ] = [ B][d ], obtain strain{}
6. Use [ ] = [E][ ] , to obtain [s]

2018/10/31 13
 CY Ou

uw  a a

 , c , f 
z

wz Ka z K a z +  w z

(a) (b) (c)

Drained analysis

08:01:36
 CY Ou

uw 

 , c , f 

(a) (b)

Effective stress undrained analysis

08:01:36
 CY Ou

 sat , su , f = 0
z

K a sat z − 2 su

Total stress undrained analysis

08:01 16
 CY Ou
3. FEA procedure

1. Define Geometry and Boundary

2. Select Constitutive model (soil model)

3. Generate Mesh

4. Generate Initial Condition

5. Model Construction Procedure

6. Check Convergence

7. Validate Analysis Results

08:01:36
 CY Ou
Validation (KH Yang)
Validate numerical results using

(1) Analytical or Empirical Solutions

(Equations or design charts)

(2) Laboratory Tests

(3) Field Observations

(4) Case History

Always validate your numerical results!!

08:01:36
 CY Ou
4. Selection of soil model 1

3

 3 (large)
Δ z  0

1 − 3
Δ x  0
 z − x Δ z = 0
 3 (intermediate)
Δ z  0
2 Δ x  0 Δ x = 0
C  3(small)

B A
unloading

Δ z  0 Δ z  0 1
Δ x  0 Δ x  0

D E  z + x
F Δ z = 0
2
Δ z  0
Δ x  0
Δ x = 0
Δ z  0
Δ x  0

08:01:36
 CY Ou

5. Mohr-Coulomb (MC) model

2018/10/31 20
 CY Ou

Mohr-Coulomb yield surface

2018/10/31 21
 CY Ou

Two drawbacks:

(1)The magnitude of the plastic volumetric strain (i.e.,

the dilation) is much larger than that observed in real
soils.
(2) Once the soil yields, it will dilate for ever. Real soil,
which may dilate initially on meeting the failure
surface, will often reach a constant volume condition
(i.e., zero incremental plastic volumetric strains) at
large strain

2018/10/31 22
 CY Ou
Use of Plastic potential function

Yield function, F
t, p/2
f f

 
Plastic potential function, P

, p

a
y
 CY Ou

If =f’, F＝P, Associated conditions arise.

If f’ , Non-associated conditions and as  reduces less

dilation is generated.

If =0, zero plastic dilation (i.e., no plastic volume change) occurs

{ } = {[ D] − [ D p ]}{  } = [ D ep ]{ }

Required parameters:

c, f , E, m , 

2018/10/31 24
 CY Ou

1 −3

1 −3
A B  1

E50  3
1
O
1  1
(a) (b)

Required parameters: c , f, E, m, 

Loading stress path: E, m=0.3

Unloading stress path: : Eur, m=0.3
08:01:36
 CY Ou

Undrained Analysis
q
(kPa)

u =10 kPa
30 1

2=3
ESP
TSP

100 110 p or p’ (kPa)

Loading 1 2 3 p p q
stage (kPa) (kPa) (kPa) (kPa) (kPa) (kPa)
Consolidation 100 100 100 100 100 0
Shearing 100 100 100 110 100 30

08:01:36
 CY Ou

q Kf compression line q Kf compression line q

2su Uf Kf compression line

Uf over
estimates Uf
2su 2su
ESP
ESP TSP TSP ESP
TSP
p p p
(a) Real Soil (b) Undrained A (c) Undrained B

MC Undrained A:

MC Undrained B:

08:01:36
 CY Ou

GL.+0.0 1FL GL.+0.0

CL f =30º H300x300x10x15 GL.-2.3
GL.-2.0  =33%-38% GL.-2.8 Stage 1
Sungshan VI LL=33-36 B1F GL.-3.5 GL.-4.9
PI=13-16 Stage 2
GL.-5.6 SM N=4-11

• A good case history,

Sungshan V f  =31º B2F GL.-7.1 GL.-8.6
GL.-8.0
Stage 3
B3F GL.-10.3
GL.-11.8 Stage 4
TNEC is used in B4F GL.-13.7
GL.-15.2
Stage 5
evaluation of the Sungshan IV CL
 =32%-40%
H400x400x13x21 GL.-16.5

B5F GL.-17.1
GL.-17.3 Stage 6
GL.-19.7
LL=29-39 Stage 7
performance of the PI=9-23
f =29º
unit: meter

above models
GL.-33.0
Sungshan III SM N=22-24 f =31º
GL.-35.0 GL.-35.0
Sungshan II CL N=9-11
GL.-37.5 f =29º

SM
Sungshan I N=14-37
f =32º
GL.-46.0
Gravel
GP N>100
2018/10/31 28
 CY Ou

Input parameters for the lateral support

Stage Type A(m2) t(m) s(m) E(MPa) 

1 steel 0.012 - 8 210000 0.2

2 slab - 0.15 - 21000 0.15

3 slab - 0.15 - 21000 0.15

4 slab - 0.15 - 21000 0.15

5 slab - 0.15 - 21000 0.15

6 steel 0.0219 - 3.4 210000 0.2

Note : s= spacing distance between struts ; t = thickness of slab

2018/10/31 29
 CY Ou

2018/10/31 30
 CY Ou

2018/10/31 31
 CY Ou

2018/10/31 32
 CY Ou

Displacement (cm) Displacement (cm) Displacement (cm)

0 1 2 3 4 0 2 4 6 8 0 2 4 6 8 10 12

0 0 0

10 10 10

20 20 20
Depth (m)

30 30 30

40 40 40

(a) Initial stage (b) Middle stage (c) Final stage

50 50 50

Field measurements HSS Model MC Undrained B Model

HS Model MC Undrained A Model

2018/10/31 33
 CY Ou

Distance from the wall (m) Distance from the wall (m) Distance from the wall (m)
0 20 40 60 80 100 0 20 40 60 80 100 0 20 40 60 80 100

0 0 0
Ground settlement (cm)

-1 -2
-0.2

-2 -4
-0.4
-3 -6

-0.6
-4 -8

-0.8 (a) Initial stage -5 (b) Middle stage -10 (c) Final stage
Field measurements HSS Model MC Undrained B Model
HS Model MC Undrained A Model

2018/10/31 34
 CY Ou
6. Modified Cam-clay (MCC) model

Eur is called pure elastic Young’s modulus or unloading-reloading Young’s modulus

mur is called pure elastic Poisson’s ratio or unloading-reloading Poisson’s ratio

2018/10/31 35
 CY Ou

2018/10/31 36
 CY Ou

Assuming in the laboratory triaxial test condition:

1
p  = ( 1 +  2 +  3 ) q =  1 −  3
3 ,

where 1  2 3 are the major, intermediate

and minor principal stresses.

2018/10/31 37
The critical state of soil refers to the state where CY Ou

volumetric strain is not further produced with the

increase of shear strain at large shear strain. The state
is usually either the failure state or the ultimate state.

900
Critical state
 1 −  3 (kPa)

450

0.00
v

0.05

0.10
0 0.1 0.2
1
2018/10/31 38
 CY Ou
The critical state can be expressed by the following two
equations:

q = M p
e = e cs −  ln p 
e= void ratio
M = void ratio when p =1.0 on the critical state line.
ecs = slope of the projection of the critical state line on the
p - q plane.
 = slope of the projection of the critical state line on the
e- ln p  plane.

2018/10/31 39
 CY Ou
Critical state line (CSL)
q e
One dimensional consolidation line
ea Isotropic virgin consolidation line(IVCL)


M e cs
Unloading/reloading
e


p ln p
p = 1

(a)
(b)

2018/10/31 40
 CY Ou
Assume the soil consolidation is in the isotropic consolidation
state ( 1 =  2 = 3 ).
The the e-p’ curve obtained from the isotropic consolidation
is called the virgin isotropic consolidation line, as shown in
Figure 8.17. It can be expressed by the following equation:

e = e a −  ln p 
Where ea is the void ratio when p’ = one unit.
 = the slope. In the e-p’ coordinate, the critical state line, the
virgin one-dimensional consolidation line, and the isotropic
virgin consolidation line parallel one another.

2018/10/31 41
 CY Ou

If we designate the slope as k, the unloading/reloading

equation can be expressed as follows:

e = e  −  ln p 
Where e = void ratio when p’ = one unit.
 is the slope of the unloading/reloading line.

2018/10/31 42
 CY Ou
q

q = Mp 
X
Critical state line (CSL)

X
X
Elastic wall e

F
E Unloading/reloading
Y Y
E
Yielding surface

Isotropic virgin consolidation line

p

FIGURE 8.18 State boundary surface

2018/10/31 43
 CY Ou

Critical state soil mechanics assumes there exists a single state

boundary surface for normally consolidated or lightly
oveconsolidated soils in the e - p - q space,
as shown in Figure 8.18. The state boundary surface intersects
the plane q=0 along line YY, which is called the virgin
isotropic consolidated line. The border line of the state
boundary surface is line XX. The state boundary surface
represents the stress states in the various ultimate conditions
(the coordinates of the stress states are e , p , q).
That is, the stress states of soil exist only on or below the state
boundary surface.

2018/10/31 44
 CY Ou

Assume the stress state of soil is on the swelling line (EF), only
elastic deformation is produced, i.e., no plastic deformation is
produced. Furthermore, all paths that remain on the curved
vertical plane above the swelling line EF, but below the state
boundary surface, will only induce elastic deformation.
Thus, this curve plane is called the elastic wall. Line EX, the
state boundary surface intersecting with the elastic wall, is
then projected onto the plane of e = 0 , which forms line E'X'.
E'X' is generally defined as a yielding surface. As discussed
above, if the state boundary surface is known, we can then
derive the yielding function.

2018/10/31 45
 CY Ou

According to Roscoe and Burland (1968), the equations for

the state boundary surface with the modified Cam-clay
model are

p   2 (1 −  /  )
 ea − e 
=  2
M

2 
p e = exp  
pe  M + q / p    
2

 M2 
The yielding equation p = p 0  
2 
M
2
+q /p 
2

where p0 is the p -value when q = 0

2018/10/31 46
 CY Ou

Yield surface

2018/10/31 47
 CY Ou

With the yielding equation given above, we can then

obtain the constitutive law for the Cam-clay
elastoplastic model

{ } = {[ D] − [ D p ]}{  } = [ D ep ]{ } (6.15)

The parameters related to the state boundary surface or the

yielding surface are M , , and  .
Thus, if the Cam-clay model is adopted as the yielding equation,
in addition to Eur and mur , the extra required input parameters
are M, , and  .

2018/10/31 48
 CY Ou

Input parameters of drained material (sand ) for the

Mohr-Coulomb model

E=2×β×G(1+ υ) b=0.5

G= ρVs

2018/10/31 49
 CY Ou

Input parameters of undrained material (clay) for the

Modified Cam-Clay model

2018/10/31
50
 CY Ou

Comparison
Displacement (cm) Distance from the wall (m)
• MCC model 0
141210 8 6 4 2 0 0 10 20 30 40 50 60 70 80 90 100
0

Settlement (cm)
1 23
5 -2
10 2 4
15
3 5 -4
5 4

Depth (m)
6
20 7 6 excavation stage
-6
25 7
-8
30
35 Field measurement
40 MCC model
45 real soil parameters
parameter
50
Displacement (cm) Distance from the wall (m)
141210 8 6 4 2 0 0 10 20 30 40 50 60 70 80 90 100
0 0

Settlement (cm)
2
5 3 -2
1 4
10
15 2 -4
3 5
Depth (m)

4 6 -6
20 7 6 5 excavation stage
7
25 -8
30
35 Field measurement
40 MCC model
45 2018/10/31
51
50 adjusted parameters
 CY Ou

Adjusted input parameters of undrained material (clay)

for the Modified Cam-clay model

2018/10/31
52
 CY Ou
q

Critical state line

Real soil yield surface

K0 line

B
p’
E C
D

Modified Cam-Clay yield surface

Relationship of stress path in modified Cam-Clay yield surface

and real soil yield surface
 CY Ou
7. HS model q

qa
E50 failure
qf 1

c
b 1 q
a −ε1 =
Ei 1 − q /qa
Eur

1
Eur

p e
ε 1 ε 1
-ε1
ε1

Mohr- Coulomb failure line

Deviatoric stress, q=|σ1-σ3|

2018/10/31 Mean effective stress, p' 54

 CY Ou

2018/10/31 55
 CY Ou

Alternative stiffness parameters

2018/10/31 56
 CY Ou

2018/10/31 57
 CY Ou

Summary of parameter estimation for clay

c f =0
Eurref = f (C s )

E50ref = Eurref /(3 ~ 5) = Eurref / 3.5 (NC Clay)

E50ref = Eurref /( 2 ~ 3) = Eurref / 2.5 (OC Clay)

E ref
50 =E ref
ur /3
ref
Eoed = 0.7 E50ref

m=0.5~1 for various soil; m=1 for soft clay

08:01 58
 CY Ou

Summary of parameter estimation for sand

c f =f−30
E50ref = f ( SPT − N ) 𝑟𝑒𝑓
𝐸𝑢𝑟 = 𝑓(𝑆𝑃𝑇 − 𝑁൯
Eurref = f ( SPT - N )
Eurref = (3 ~ 5) E50ref = 3.5 E50ref (Loose sand)

Eurref = (2 ~ 3) E50ref = 2.5 E50ref (dense sand)

Eurref = 3E50ref (default)

m=0.5~1 for various soil; m=0.5 for sand

08:01 59
 CY Ou

, Comparison
Calvello and Finno (2004) suggested:
Eref50= 1/3 Erefur , Erefoed = 0.7 Eref50

m=1  ur = 0.2
• HS model
Displacement (cm) Distance from the wall (m)
141210 8 6 4 2 0 0 10 20 30 40 50 60 70 80 90 100
0 0

Settlement (cm)
2
5 3 -2
1
10 4
15 2 -4
3 5
Depth (m)

4 -6
20 7 6 6
excavation stage
5 7
25 -8
30
35 Field measurement
40 HS model
45
50
08:01
60
 CY Ou
8 .Hardening soil model with small strain stiffness (HSsmall)

2018/10/31 61
 CY Ou
Additional Soil Parameter :

0.7

2018/10/31 62
 CY Ou

The Additional Stiffness Parameters of HS-Small

2018/10/31 63
 CY Ou

Other Correlation of Threshold Shear Strain

2018/10/31 64
 CY Ou

1
PLAXIS suggested:  0.7 = 9G [2c(1 + cos 2f ) −  1 (1 + K 0 ) sin 2f ]
0

Displacement (cm) Distance from the wall (m)

• HS small model 141210 8 6 4 2 0 0 10 20 30 40 50 60 70 80 90 100
0 0

Settlement (cm)
5 2
1 3 -2
10 2 4
15 54 3 5 -4

Depth (m)
20 6 6 -6
25 7 7 excavation stage -8
30
35 Field measurement
40 HS small model
45 real soil parameter
50
Displacement (cm) Distance from the wall (m)
141210 8 6 4 2 0 0 10 20 30 40 50 60 70 80 90 100
0 0

Settlement (cm)
2
5 1 -2
10 3
2
15 3 4 -4
5
Depth (m)

5 4
20 6 -6
6
25 7 excavation stage -8
7
30
35 Field measurement
40 HS small model
45 2018/10/31
65
50 adjusted parameter,  0.7 = 10 -5
 CY Ou
9. MC model (total stress undrained analysis)

Displacement (cm) Distance from the wall (m)

141210 8 6 4 2 0 0 10 20 30 40 50 60 70 80 90 100
0 0

Settlement (cm)
2
5 1 -2
2 3
10 4
15
3 -4
4
Depth (m)

5 5
20 7 6 -6
6
25 7 excavation stage -8
30
35 Field measurement
40 Mohr-Coulomb model
45
50

2018/10/31 66
 CY Ou
10. USC model (Total stress undrained analysis)

Displacement (cm) Distance from the wall (m)

141210 8 6 4 2 0 0 10 20 30 40 50 60 70 80 90 100

Settlement (cm)
0 12 0
5 1 3
-2
10 2 4
3 5 -4
15 5 4
Depth (m)

6 6 -6
20 7
25 7 excavation stage -8
30
35
40 Field measurement
USC model
45
50
real soil parameter

Stress path dependent total stress model

67 2018/10/31
 CY Ou
11. Stability analysis

FEM with c, f reduction

Strength reduction ratio (SR)

Divergence of
numerical solutions
𝑐𝑜𝑟𝑔𝑖𝑛𝑎𝑙
𝑐𝑖𝑛𝑝𝑢𝑡 = SR1
𝑆𝑅
SR3
∅𝑜𝑟𝑖𝑔𝑖𝑛𝑎𝑙 Upper part
∅𝑖𝑛𝑝𝑢𝑡 =
𝑆𝑅

Lower part

Nodal displacement
Convergence criterion: FS = SR1
68
 CY Ou

Taipei Rebar
Broadway case

69
 CY Ou

Taipei Rebar Broadway Case

B = 25.8 m
5.1 m 5.1 m 5.4 m 5.1 m 5.1 m
GL 0.0 m GL -1.0 m
GL -1.8 m H300x300x10x15 SM  t = 20.3 kN/m3
GL -3.55 m
GL -4.35 m H350x350x12x19

GL -7.65 m
GL -6.85m ML-CL  t = 19.2 kN/m3
H400x400x13x21 Su (kPa)
GL -10.15m SM-SP  t = 19.7 kN/m3 -10 0 50 100
GL -10.95 m H400x400x13x21
GL -13.45 m CL  t = 18.8 kN/m3
-15

Wall
(0.7 m CKoU-AC -20
thickness) CKoU-AE
GL -24.0 m CKoU-DSS
Center post For analysis -25

-30
70
 CY Ou

N, kN/m
-14000 1st strut layer
1
2 2nd strut layer
3 -10000 3rd strut layer
4 -1000 5 4th strut layer
5 5 wall
4
4
-500
Excavation 3 5
stage
5 2 3 1 23 4
0 M kNm/m
-1000 -600 -200 -100 0 4 100 200 600 1000
5
BL of 1st strut layer
500
BL of 2nd strut layer

1000 BL of 3rd and 4th strut layers

10000
BL of wall
14000
71
 CY Ou

32 m

Plastic point plot at SRmax as using EP support system

72