11th International Conference of IACMAG, Torino 21 Giugno 2005 Exact bearing capacity calculations using the method of characteristics

1
11th International Conference of IACMAG,
Torino21 Giugno 2005Exact bearing capacity
calculations using the method of characteristics

• Dr C.M. Martin
• Department of Engineering Science
• University of Oxford

2
Outline

• Introduction
• Bearing capacity calculations using the method
  of characteristics
• Exact solution for example problem
• Can we solve the Ng problem this way?
• The fast (but apparently forgotten) way to find
  Ng
• Verification of exactness
• Conclusions

3
Bearing capacity

• Idealised problem (basis of design methods)

Central, purely vertical loading
qu Qu/B
Rigid strip footing
D
B
Semi-infinite soil c, f, g, y f
4
Bearing capacity

• Idealised problem (basis of design methods)

Central, purely vertical loading
qu Qu/B
q gD
q gD
Rigid strip footing
B
Semi-infinite soil c, f, g, y f
5
Classical plasticity theorems

• A unique collapse load exists, and it can be
  bracketed by lower and upper bounds (LB, UB)
• LB solution from a stress field that satisfies
• equilibrium
• stress boundary conditions
• yield criterion
• UB solution from a velocity field that satisfies
• flow rule for strain rates
• velocity boundary conditions
• Theorems only valid for idealised material
• perfect plasticity, associated flow (y f)

Statically admissible
Plastically admissible

Kinematically admissible
6
Method of characteristics

• Technique for solving systems of quasi-linear
  PDEs of hyperbolic type
• Applications in both fluid and solid mechanics
• In soil mechanics, used for plasticity problems
• bearing capacity of shallow foundations
• earth pressure on retaining walls
• trapdoors, penetrometers, slope stability,
• Method can be used to calculate both stress and
  velocity fields (hence lower and upper bounds)
• In practice, often gives LB UB ? exact result
• 2D problems only plane strain, axial symmetry

7
Method of characteristics

• Technique for solving systems of quasi-linear
  PDEs of hyperbolic type
• Applications in both fluid and solid mechanics
• In soil mechanics, used for plasticity problems
• bearing capacity of shallow foundations
• earth pressure on retaining walls
• trapdoors, penetrometers, slope stability,
• Method can be used to calculate both stress and
  velocity fields (hence lower and upper bounds)
• In practice, often gives LB UB ? exact result
• 2D problems only plane strain, axial symmetry

8
Outline

• Introduction
• Bearing capacity calculations using the method
  of characteristics
• Exact solution for example problem
• Can we solve the Ng problem this way?
• The fast (but apparently forgotten) way to find
  Ng
• Verification of exactness
• Conclusions

9
Lower bound stress field

• To define a 2D stress field, e.g. in x-z plane
• normally need 3 variables (sxx, szz, txz)
• if assume soil is at yield, only need 2 variables
  (s, q)

x
t ?
X
q
c
sn
s
s1
s3
f
2q
s3 s R
Z
M-C
s1 s R


z
general
10
Lower bound stress field

• To define a 2D stress field, e.g. in x-z plane
• normally need 3 variables (sxx, szz, txz)
• if assume soil is at yield, only need 2 variables
  (s, q)

x
t ?
b
X
b
q
c
sn
s
s1
s3
f
a
2q
s3 s R
Z
a
a
M-C
s1 s R
b


z
general
11
Lower bound stress field

• To define a 2D stress field, e.g. in x-z plane
• normally need 3 variables (sxx, szz, txz)
• if assume soil is at yield, only need 2 variables
  (s, q)

x
t ?
b
X
b
q
e p/4 f/2
c
2e
sn
s
s1
s3
f
a
2q
s3 s R
2e
Z
a
a
e
e
M-C
s1 s R
b


z
general
12
Lower bound stress field

• Substitute stresses-at-yield (in terms of s, q)
  into equilibrium equations
• Result is a pair of hyperbolic PDEs in s, q
• Characteristic directions turn out to coincide
  with a and b slip lines aligned at q ? e
• Use a and b directions as curvilinear coords ?
  obtain a pair of ODEs in s, q (easier to
  integrate)
• Solution can be marched out from known BCs

13
Lower bound stress field

• Substitute stresses-at-yield (in terms of s, q)
  into equilibrium equations
• Result is a pair of hyperbolic PDEs in s, q
• Characteristic directions turn out to coincide
  with a and b slip lines aligned at q ? e
• Use a and b directions as curvilinear coords ?
  obtain a pair of ODEs in s, q (easier to
  integrate)
• Solution can be marched out from known BCs

gt 0
14
Lower bound stress field

• Marching from two known points to a new point

(xB, zB, sB, qB)
B
(xA, zA, sA, qA)
A
15
Lower bound stress field

• Marching from two known points to a new point

(xB, zB, sB, qB)
B
(xA, zA, sA, qA)
A
?
?
C
(xC, zC, sC, qC)
16
Lower bound stress field

• Marching from two known points to a new point

(xB, zB, sB, qB)
B
(xA, zA, sA, qA)
A
?
?
C
(xC, zC, sC, qC)
17
Lower bound stress field

• Marching from two known points to a new point
• One-legged variant for marching from a known
  point onto an interface of known roughness

(xB, zB, sB, qB)
B
(xA, zA, sA, qA)
A
?
?
C
(xC, zC, sC, qC)
18
Lower bound stress field

• Marching from two known points to a new point
• One-legged variant for marching from a known
  point onto an interface of known roughness

(xB, zB, sB, qB)
B
(xA, zA, sA, qA)
A
?
?
C
(xC, zC, sC, qC)
FD form
FD form
19
Upper bound velocity field

• Substitute velocities u, v into equations for
• associated flow (strain rates normal to yield
  surface)
• coaxiality (princ. strain dirns princ. stress
  dirns)
• Result is a pair of hyperbolic PDEs in u, v
• Characteristic directions again coincide with the
  a and b slip lines aligned at q ? e
• Use a and b directions as curvilinear coords ?
  obtain a pair of ODEs in u, v (easier to
  integrate)
• Solution can be marched out from known BCs

20
Upper bound velocity field

• Marching from two known points to a new point

x,u
(xB, zB, sB, qB, uB, vB)
B
(xA, zA, sA, qA, uA, vA)
A
z,v
21
Upper bound velocity field

• Marching from two known points to a new point

x,u
(xB, zB, sB, qB, uB, vB)
B
(xA, zA, sA, qA, uA, vA)
A
z,v
?
?
C
(xC, zC, sC, qC, uC, vC)
22
Upper bound velocity field

• Marching from two known points to a new point

x,u
(xB, zB, sB, qB, uB, vB)
B
(xA, zA, sA, qA, uA, vA)
A
z,v
?
?
C
(xC, zC, sC, qC, uC, vC)
23
Upper bound velocity field

• Marching from two known points to a new point
• One-legged variant for marching from a known
  point onto an interface of known roughness

x,u
(xB, zB, sB, qB, uB, vB)
B
(xA, zA, sA, qA, uA, vA)
A
z,v
?
?
C
(xC, zC, sC, qC, uC, vC)
24
Upper bound velocity field

• Marching from two known points to a new point
• One-legged variant for marching from a known
  point onto an interface of known roughness

x,u
(xB, zB, sB, qB, uB, vB)
B
(xA, zA, sA, qA, uA, vA)
A
z,v
?
?
C
FD form
FD form
(xC, zC, sC, qC, uC, vC)
25
Outline

• Introduction
• Bearing capacity calculations using the method
  of characteristics
• Exact solution for example problem
• Can we solve the Ng problem this way?
• The fast (but apparently forgotten) way to find
  Ng
• Verification of exactness
• Conclusions

26
Example problem
Rough base
qu
q 18 kPa
q 18 kPa
B 4 m
c 16 kPa, f 30, g 18 kN/m3
after Salençon Matar (1982)
27
Example problem stress field (partial)
a
s known (passive failure) q p/2
b
28
Example problem stress field (partial)
a
s known (passive failure) q p/2
b
Symmetry q 0 on z axis (iterative construction
reqd)
29
Example problem stress field (partial)
a
s known (passive failure) q p/2
b
Symmetry q 0 on z axis (iterative construction
reqd)

• Shape of false head region emerges naturally
• qu from integration of tractions
• Solution not strict LB until stress field
  extended

30
Example problem stress field (complete)
Minor principal stress trajectory
31
Example problem stress field (complete)

• Extension strategy by Cox et al. (1961)
• Here generalised for g gt 0
• Utilisation factor at start of each spoke must
  be ? 1

Minor principal stress trajectory
32
Extension technique
q
z0
s1
z
s3
s1
gz0 q
s1 g(z ? z0)
gz q
33
Extension technique
q
z0
s1
z
s3
s1
gz0 q
Critical utilisation is here
s1 g(z ? z0)
gz q
34
Example problem velocity field
Rigid
Rigid
Rigid
Rigid
Rigid
35
Example problem velocity field
Rigid
Rigid
Rigid
Rigid
Rigid

• Discontinuities are easy to handle treat as
  degenerate quadrilateral cells (zero area)

36
Some cautionary remarks

• Velocity field from method of characteristics
  does not guarantee kinematic admissibility!
• principal strain rates may become
  mismatched with principal stresses s1, s3
• this is OK if f 0 (though expect UB ? LB)
• but not OK if f gt 0 flow rule violated ? no UB
  at all
• If f gt 0, as here, must check each cell of mesh
• condition is sufficient
• Only then are calculations for UB meaningful
• internal dissipation, e.g. using
• external work against gravity and surcharge

37
Example problem velocity field
Rigid
Rigid
Rigid
Rigid
Rigid

• qu from integration of internal and external work
  rates for each cell (4-node ?, 3-node ?)
• Discontinuities do not need special treatment

38
Convergence of qu (kPa) in example
39
Convergence of qu (kPa) in example
LB
40
Convergence of qu (kPa) in example
UB
LB
41
Outline

• Introduction
• Bearing capacity calculations using the method
  of characteristics
• Exact solution for example problem
• Can we solve the Ng problem this way?
• The fast (but apparently forgotten) way to find
  Ng
• Verification of exactness
• Conclusions

42
Why not?
The solutions obtained from the method of
characteristics are generally not exact collapse
loads, since it is not always possible to
integrate the stress-strain rate relations to
obtain a kinematically admissible velocity ?eld,
or to extend the stress ?eld over the entire
half-space of the soil domain.
Hjiaj M., Lyamin A.V. Sloan S.W. (2005).
Numerical limit analysis solutions for the
bearing capacity factor Ng. Int. J. Sol. Struct.
42, 1681-1704.
43
Ng problem as a limiting case
qu
d
q
q
B
c 0, f gt 0, g gt 0, y f
44
Ng problem as a limiting case
qu
d
q
q
B
c 0, f gt 0, g gt 0, y f
45
Ng problem as a limiting case
qu
d
q
q
B
c 0, f gt 0, g gt 0, y f
46
Ng problem as a limiting case
qu
d
q
q
B
c 0, f gt 0, g gt 0, y f
47
Stress field as gB/q ? ?
c 0, f 30, Rough (d f)
48
Stress field as gB/q ? ?
c 0, f 30, Rough (d f)
49
Stress field as gB/q ? ?
c 0, f 30, Rough (d f)
50
Stress field as gB/q ? ?
c 0, f 30, Rough (d f)
51
Stress field as gB/q ? ?
c 0, f 30, Rough (d f)
52
Stress field as gB/q ? ?
c 0, f 30, Rough (d f)
53
Stress field as gB/q ? ?
c 0, f 30, Rough (d f)
54
Stress field as gB/q ? ?
c 0, f 30, Rough (d f)
55
Stress field as gB/q ? ?
c 0, f 30, Rough (d f)
56
Stress field as gB/q ? ?
c 0, f 30, Rough (d f)
57
Stress field as gB/q ? ?
c 0, f 30, Rough (d f)
58
Stress field as gB/q ? ?
c 0, f 30, Rough (d f)
59
Stress field as gB/q ? ?
c 0, f 30, Rough (d f)
60
Stress field as gB/q ? ?
c 0, f 30, Rough (d f)
61
Stress field as gB/q ? ?
c 0, f 30, Rough (d f)
62
Stress field as gB/q ? ?
c 0, f 30, Rough (d f)
63
Stress field as gB/q ? ?
c 0, f 30, Rough (d f)
64
Stress field as gB/q ? ?
c 0, f 30, Rough (d f)
65
Stress field as gB/q ? ?
c 0, f 30, Rough (d f)
66
Stress field as gB/q ? ?
c 0, f 30, Rough (d f)
67
Stress field as gB/q ? ?
Take as Ng
Fan (almost) degenerate
c 0, f 30, Rough (d f)
68
Velocity field as gB/q ? ?
c 0, f 30, Rough (d f)
69
Velocity field as gB/q ? ?
c 0, f 30, Rough (d f)
70
Velocity field as gB/q ? ?
c 0, f 30, Rough (d f)
71
Velocity field as gB/q ? ?
c 0, f 30, Rough (d f)
72
Velocity field as gB/q ? ?
c 0, f 30, Rough (d f)
73
Velocity field as gB/q ? ?
c 0, f 30, Rough (d f)
74
Velocity field as gB/q ? ?
c 0, f 30, Rough (d f)
75
Velocity field as gB/q ? ?
c 0, f 30, Rough (d f)
76
Velocity field as gB/q ? ?
c 0, f 30, Rough (d f)
77
Velocity field as gB/q ? ?
c 0, f 30, Rough (d f)
78
Velocity field as gB/q ? ?
c 0, f 30, Rough (d f)
79
Velocity field as gB/q ? ?
c 0, f 30, Rough (d f)
80
Velocity field as gB/q ? ?
c 0, f 30, Rough (d f)
81
Velocity field as gB/q ? ?
c 0, f 30, Rough (d f)
82
Velocity field as gB/q ? ?
c 0, f 30, Rough (d f)
83
Velocity field as gB/q ? ?
c 0, f 30, Rough (d f)
84
Velocity field as gB/q ? ?
c 0, f 30, Rough (d f)
85
Velocity field as gB/q ? ?
c 0, f 30, Rough (d f)
86
Velocity field as gB/q ? ?
c 0, f 30, Rough (d f)
87
Velocity field as gB/q ? ?
c 0, f 30, Rough (d f)
88
Velocity field as gB/q ? ?
Take as Ng
Fan (almost) degenerate
c 0, f 30, Rough (d f)
89
Convergence of 2qu/gB when gB/q 109
90
Convergence of 2qu/gB when gB/q 109
LB
91
Convergence of 2qu/gB when gB/q 109
UB
LB
92
Completion of stress field (coarse)
c 0 f 30 gB/q 109 Rough (d f) Ng
14.7543
93
Completion of stress field (fine)
c 0 f 30 gB/q 109 Rough (d f) Ng
14.7543
94
Completion of stress field (fine)
c 0 f 30 gB/q 109 Rough (d f) Ng
14.7543
EXACT
95
It also works for smooth footings
c 0 f 30 gB/q 109 Smooth (d 0) Ng
7.65300
96
and other friction angles
c 0 f 20 gB/q 109 Rough (d f) Ng
2.83894
97
Outline

• Introduction
• Bearing capacity calculations using the method
  of characteristics
• Exact solution for example problem
• Can we solve the Ng problem this way?
• The fast (but apparently forgotten) way to find
  Ng
• Verification of exactness
• Conclusions

98
Notice anything?
c 0 f 30 gB/q 109 Smooth (d 0) Ng
7.65300

• Tractions ? distance from singular point
• Characteristics self-similar w.r.t. singular point

99
Recall Ng problem definition
q 0
Semi-infinite soil c 0, f gt 0, g gt 0
100
Recall Ng problem definition
q 0
y
r
Semi-infinite soil c 0, f gt 0, g gt 0
101
Governing equations

• No fundamental length ? can solve in terms of
  polar angle y and radius r
• Along a radius, stress state varies only in
  scale
• mean stress s ? r
• major principal stress orientation q const
• Combine with yield criterion and equilibrium
  equations to get a pair of ODEs

von Kármán (1926)
102
Direct solution of ODEs
Underside of footing (d 0)
y
r
Edge of passive zone
solve (iteratively)
103
Direct solution of ODEs

• Use any standard adaptive Runge-Kutta solver
• ode45 in MATLAB, NDSolve in Mathematica
• Easy to get Ng factors to any desired precision
• Much faster than method of characteristics
• Definitive tables of Ng have been compiled for
• f 1, 2, , 60
• d/f 0, 1/3, 1/2, 2/3, 1
• Values are identical to those obtained from the
  method of characteristics, letting gB/q ? ?

lt 10 s to generate
104
Selected values of Ng

• Exactness checked by method of characteristics
  LB UB, stress field extensible, match

105
Selected values of Ng

• Exactness checked by method of characteristics
  LB UB, stress field extensible, match

106
Influence of roughness on Ng
d/f 2/3
d/f 1/2
d/f 1/3
Smooth
0.504719
0.500722
0.500043
107
Outline

• Introduction
• Bearing capacity calculations using the method
  of characteristics
• Exact solution for example problem
• Can we solve the Ng problem this way?
• The fast (but apparently forgotten) way to find
  Ng
• Verification of exactness
• Conclusions

108
Ng by various methods
f 30, d f
FELA
Limit Eqm
Characteristics
Upper Bd
FE/FD
ODEs
Formulae
109
Ng by various methods
f 30, d f
FELA
Limit Eqm
Characteristics
Upper Bd
FE/FD
ODEs
Formulae
110
Ng by FE limit analysis
Ukritchon et al. (2003)
Rough
UPPER BOUND
Smooth
LOWER BOUND
Smooth
Rough
111
Ng by FE limit analysis
Hjiaj et al. (2005)
UPPER BOUND
Smooth
Rough
Smooth
LOWER BOUND
Rough
112
Ng by FE limit analysis
Hjiaj et al. (2005)
UPPER BOUND
Smooth
Rough
Smooth
LOWER BOUND
Rough
113
Ng by FE limit analysis
Hjiaj et al. (2005)
UPPER BOUND
Smooth
Rough
Smooth
LOWER BOUND
Rough
114
Ng by FE limit analysis
Hjiaj et al. (2005)
UPPER BOUND
Smooth
Rough
Smooth
LOWER BOUND
Rough

• Structured meshes (different for each f)

115
Ng by FE limit analysis
Makrodimopoulos Martin (2005)
UPPER BOUND
Rough
Smooth
Smooth
Rough
LOWER BOUND
116
Ng by FE limit analysis
Makrodimopoulos Martin (2005)
UPPER BOUND
Rough
Smooth
Smooth
Rough
LOWER BOUND

• Single unstructured mesh (same for each f)

117
Ng by various methods
f 30, d f
FELA
Limit Eqm
Characteristics
Upper Bd
FE/FD
ODEs
Formulae
118
Ng (d f) by common formulae error
119
Bearing capacity factors for design

• If we use Nc and Nq that are exact for y f
• then we should, if we want to be consistent,
  also use Ng factors that are exact for y f
• Then start worrying about corrections for
• non-association (y lt f)
• stochastic variation of properties
• intermediate principal stress
• progressive failure, etc.

120
Bearing capacity factors for design

• If we use Nc and Nq that are exact for y f
• then we should, if we want to be consistent,
  also use Ng factors that are exact for y f.
• Then start worrying about corrections for
• non-association (y lt f)
• stochastic variation of properties
• intermediate principal stress
• progressive failure, etc.

? less capacity!
121
Conclusions

• Shallow foundation bearing capacity is a
  long-standing problem in theoretical soil
  mechanics
• The method of characteristics, carefully applied,
  can be used to solve it ? c, f, g (with y f)
• In all cases, find strict lower and upper bounds
  that coincide, so the solutions are formally
  exact
• If just values of Ng are required (and not proof
  of exactness) it is much quicker to integrate the
  governing ODEs using a Runge-Kutta solver
• Exact solutions provide a useful benchmark for
  validating other numerical methods (e.g. FE)

122
Downloads

• Program ABC Analysis of Bearing Capacity
• Tabulated exact values of b.c. factor Ng
• Copy of these slides

www-civil.eng.ox.ac.uk