$SSHQGL[&
62/87,216729(5,),&$7,21352%/(06,1
&+$37(5
C.1. INTRODUCTION
7KLVDSSHQGL[FRQWDLQVWKHDQDO\WLFDOVROXWLRQVIRUWKHWKUHHYHULILFDWLRQSUREOHPV
GLVFXVVHGLQ&KDSWHU
C.2. VERIFICATION PROBLEM ONE – TERZAGHI ONEDIMENSIONAL CONSOLIDATION
7HU]DJKLGHULYHGWKHIROORZLQJSDUWLDOGLIIHUHQWLDOHTXDWLRQWRGHVFULEHWKH
VHWWOHPHQWRIDVDWXUDWHGFOD\OD\HUDQGGLVVLSDWLRQRIH[FHVVSRUHSUHVVXUHVZLWK
WLPH
cv
:KHUH
(TXDWLRQ&
∂ 2 u ∂u
=
∂ 2 z ∂t
cv LVWKHFRHIILFLHQWRIFRQVROLGDWLRQ
u LVWKHH[FHVVSRUHZDWHUSUHVVXUHLQWKHVRLO
t LVWLPH
z LVGHSWK
7KHERXQGDU\DQGLQLWLDOFRQGLWLRQVDUH
u ( z ,0 ) = q
u( 0, t ) = 0
0 ≤ z ≤ 2H
IRU t > 0
u( 2 H , t ) = 0
+ LV WKH OHQJWK RI WKH ORQJHVW GUDLQDJH SDWK 7KH VHSDUDWLRQ RI YDULDEOHV
WHFKQLTXHLVXVHGWRREWDLQDVROXWLRQWRWKHSDUWLDOGLIIHUHQWLDOHTXDWLRQIRUWKH
$SSHQGL[&
'DYLG-%HQWOHU
SUHVFULEHGERXQGDU\DQGLQLWLDOFRQGLWLRQV7KHUHVXOWLQJVROXWLRQIRU u( z , t ) LVD
FRQYHUJHQWLQILQLWHVHULHV
∞
u
4
(2 n + 2)π z
(2 n + 1) 2 π 2
=∑
T
sin
exp −
q n = 0 (2 n + 1)π
2
H
4
(TXDWLRQ&
:KHUH
T LVWKHGLPHQVLRQOHVVWLPHIDFWRU =
cv t
H2
7KHGHJUHHRIFRQVROLGDWLRQDWDQ\GHSWK]LVJLYHQE\ U z = 1 −
u
7KHDYHUDJH
q
GHJUHHRIFRQVROLGDWLRQ8LVREWDLQHGE\LQWHJUDWLQJWKHGHJUHHRIFRQVROLGDWLRQ
U z RYHUWKHHQWLUHVRLOFROXPQDQGGLYLGLQJE\WKHKHLJKWRIWKHVRLOFROXPQ
∞
U = 1− ∑
n =0
12n + 16 π
7KHH[SUHVVLRQVIRU
(2n + 1) π T "#
! 4 $
2
8
2
2
exp
2
(TXDWLRQ&
u
DQG U DUHHDV\WRHYDOXDWHXVLQJHLWKHUDSURJUDPPDEOH
q
FDOFXODWRURUDVSUHDGVKHHW
C.3. VERIFICATION PROBLEM TWO – CONSOLIDATION OF
CONTIGUOUS CLAY LAYERS WITH DIFFERENT
PERMEABILITY
*UD\ GHYHORSHG DQDO\WLFDO VROXWLRQV WR JHQHUDO RQHGLPHQVLRQDO SUREOHP
RI WZR FRQVROLGDWLQJ DGMDFHQW FRPSUHVVLEOH VWUDWD *UD\·V VROXWLRQ DQG LWV
QXPHULFDOHYDOXDWLRQIRUYHULILFDWLRQSUREOHPDUHGHPRQVWUDWHGLQWKHIROORZLQJ
0DWKFDGZRUNVKHHW'FRQWLJXRXVOD\HUVPFG
Verification Problem Two - consolidation of contiguous clay layers
with different permeability
The purpose of this Mathcad worksheet is to evaluate the expressions for
the excess pore pressure in two contiguous clay of unlike compressibility
developed by Gray, H. (1944), "Simultaneous consolidation of contiguous
layers of unlike compressible soils", ASCE Transactions, No. 2258, pp.
1327-1356.
I
a1
k1
c1
e1
h1
Case a. Free drainage at
top and bottom
Case b. Free drainage
at top only
Terminology and definitions
stress-strain relations:
the soil is assumed to be linearly elastic
k = coefficient of permeability
av = coefficient of compressibility
e = void ratio
e0 = initial void ratio
c = coefficient of consolidation
p = surface traction applied at t = 0
Hs
H
1
e
z
h
1
e
ks
µ = poisson's ratio
E = Young's modulus
D = constrained modulus
mv = coefficient of compressibility
k
1
D
e
k
c
a v .γ .( 1
z
e2
Hs
II
a2
k2
c2
h2 = νh1
The general cases of two adjacent compressible strata:
E .( 1 µ )
E
G
.
.
.
(1 µ ) (1 2 µ )
2 (1 µ )
mv
e)
1
D
av
mv
1
e0
Define four dimensionless numbers (µ, σ, ν, T):
µ
2
c1
c2
σ
1 . k s1
µ k s2
ν
h2
h1
T
c1
h1
u1 = excess pore pressure in layer I
u2 = excess pore pressure in layer II
U1 = average degree of consolidation in layer I
U2 = average degree of consolidation in layer II
2
.t
Note: The time factor, T, is based only on
the properties of layer I and the time, t.
327
Appendix C
328
Gray (1944) developed the following analytical solutions for the two cases:
Case a
∞
∑
u1 =
C n e − TA n sin( A n
2
n =1
∞
∑
u2 =
z
) sin( µν A n )
h1
C n e − TA n sin( A n ) sin( 1 + ν −
2
n =1
∞
U 1 = 1 − 2∑
sin(µνAn )(σ sin(µνAn ) + sin An )
z
)µ An
h1
(1 − cos An )e −TA
2
n
(An ) (σ sin (µνAn ) + µν sin ( An ))
2 ∞ sin An (σ sin µνAn + sin An )
(1 − cos µνAn )e −TA
U2 = 1−
∑
2
2
2
µν n=1 (An ) (σ sin µνAn + µν sin An )
2
n =1
2
2
2
n
in which An must be a root of
F ( A ) = σ cos A sin µν A + sin A cos µν A
and
Cn = 2 p
σ sin µν A n + sin µν A n
σ A n sin 2 µν A n + µν sin 2 A n
Case b
u1 =
u2 =
∞
∑
n =1
∞
∑
n =1
2
C n e − TA
2
n
∞
U
U
1
2
z
) sin( µν A n )
h1
C n e − TA n cos( A n
= 1 − 2∑
n =1
= 1 −
2
µν
(A n )
cos( A n ) sin( 1 + ν −
z
)µ An
h1
sin( A n )(cos( A n ) sin µν A n )
2
( σ sin
∞
2
cos
∑ (A ) (σ
2
n =1
n
( µν A n ) + µν cos
2
sin
2
2
G ( A ) = σ sin A sin µν A + cos A cos µν A
and
σ A n sin
2
)
µν A n + µν cos
in which An must be a root of
Cn =
2
( A n ))
A n (1 − cos µν A n
2 p cos A n
µν A n + µν A n cos
2
An
e − TA n
2
An
)
e
− TA
2
n
Verification Problem Two - consolidation of contiguous clay layers
with different permeability
329
Develop solution for Case b using Gray's solution. First, set up procedure for finding roots of G(A).
G( A , σ , µ , ν )
σ .sin( A ) .sin( µ .ν .A )
cos( A ) .cos( µ .ν .A )
Define the function zbrak which will bracket the first n roots of G(A,σ,µ,ν). The brackets for each root
are returned as a row in an array. The values of the left and right brackets are in column zero and
one, respectively. The variable inc controls the step size of A that zbrak uses when incrementing A.
zbrak( n , σ , µ , ν , inc )
W
0
A
0
B
A
fa
G( A , σ , µ , ν )
for i ∈ 1 .. n
check
0
while check 0
B
B
inc
G( B , σ , µ , ν )
if fa .fb < 0
fb
check
1
Wi
1,0
A
Wi
1,1
B
A
B
fa
fb
if fa .fb > 0
fa
fb
A
B
W
Define the function brent, which uses Brent's method for finding the root of a function known to lie
between a and b. The root is refined until its accuracy is tol or the maximum no. of iterations is reached.
First, define the functions test and test1 which will perform logical tests for brent.
test( fb , fc )
if fc > 0
if fb > 0
1
break
if fc < 0 if fb < 0 otherwise
1
break
1
test1( e , toli , fa , fb )
1 if
fa > fb
1 otherwise
if
e
toli
Machine ε:
ε mach
1 .10
15
Appendix C
Reference for Brent's method: Brent, R.P. 1973, Algorithms for Minimization without Derivatives
(Englewood Cliffs, NJ: Prentice-Hall), Chapters 3, 4.
brent( a , b , σ , µ , ν , tol )
itmax
fa
100
G( a , σ , µ , ν )
G( b , σ , µ , ν )
break if fa .fb > 0
fb
c
a
fc
fa
d
b
e
d
a
for i ∈ 1 .. itmax
if test( fb , fc ) > 0
c
if
a
fc
fa
d
b
e
d
a
fc < fb
a
b
b
c
c
a
fa
fb
fb
fc
fc
fa
2 .ε mach . b
0.5 .( c b )
toli
xm
0.5 .tol
b
break if xm
toli
break if fb 0
if test1( e , toli , fa , fb ) > 0
s
fb
fa
if a c
p
q2
2 .xm .s
1.0
otherwise
fa
q2
fc
r
fb
s
330
Verification Problem Two - consolidation of contiguous clay layers
with different permeability
r
fc
p s .( 2 .xm .q2 .( r q2 ) ( b a ) .( 1.0
q2 ( q2 1.0 ) .( r 1.0 ) .( s 1.0 )
q2 if p > 0
q2
p
p
if 2 .p < min
e
d
d
p
q2
3 .xm .q2
toli .q2
e .q2
otherwise
d xm
e
d
otherwise
d xm
e
a
b
b
fa
fb
b
b
b
b
fb
d
d if
d > toli
xm
toli .
xm
G( b , σ , µ , ν )
otherwise
r) )
331
Appendix C
332
Define remaining functions necessary for solving an example problem for Case b.
INC
0.1
i
the variable i will determine how many terms are evaluated when
approximating the infinite series.
12
u 1 p , σ , µ , ν , z, h 1 , T
sum
:excess pore pressure
in layer I
0
zbrak( i , σ , µ , ν , INC )
W
for n ∈ 1 .. i
An
brent Wn
1,0
, Wn
1,1
, σ , µ , ν , TOL
2 .p .cos( An )
Cn
σ .An .sin( µ .ν .An )
µ .ν .An .cos( An ) 2
2
term
Cn .exp( T .An ) .cos An .
sum
sum
2
z .
sin( µ .ν .An )
h1
term
sum
u 2 p , σ , µ , ν , z, h 1 , T
sum
0
:excess pore
pressure in
layer II
zbrak( i , σ , µ , ν , INC )
W
for n ∈ 1 .. i
An
brent Wn
1,0
, Wn
1,1
, σ , µ , ν , TOL
2 .p .cos( An )
Cn
σ .An .sin( µ .ν .An )
2
µ .ν .An .cos( An ) 2
term
2
Cn .exp( T .An ) .cos( An ) .sin µ .An . 1
sum
sum
ν
z
h1
term
sum
u p, σ , µ , ν , z, h 1, T
Data p , σ , µ , ν , h 1 , T , inc
u 1 p, σ , µ , ν , z, h 1, T
if z h 1
u 2 p, σ , µ , ν , z, h 1, T
otherwise
a0 , 0
0
a0 , 1
u p , σ , µ , ν , a0 , 0 , h 1 , T
step
ν .h 1
h1
inc
for j ∈ 1 .. inc
a
aj , 0
aj
aj , 1
u p , σ , µ , ν , aj , 0 , h 1 , T
1, 0
step
:general expression for excess pore
pressure at elevation z.
Data is a function for generating an
excess pore pressure isochrones.
Verification Problem Two - consolidation of contiguous clay layers
with different permeability
333
Example Problem (after Gray, H. 1944)
p
H
σ
100
2
µ
ν
5
Data p , σ , µ , ν , h 1 , 100 , 100
4
h1
2
e1 e2
Data p , σ , µ , ν , h 1 , 200 , 100
I
1
0.8
The solution for
this example
agrees with Gray's
solution.
z/H
0.6
0.4
0.2
0
0
0.2
0.4
T = 100
T = 200
0.6
0.8
1
Values of u/p
Example of how the first i roots of G(A) are found
A
0 , .02 .. .7
i 5
First, zbrak is called to bracket the roots
2
G( A , σ , µ , ν )
I
zbrak( i , σ , µ , ν , .05 )
j
1 .. i
Next, brent is called to find the roots.
Brent uses the brackets found by zbrak.
0
rj
,I
brent Ij
1
1, 0 j
2
0
0.5
A
1
brackets
solj
1
G rj
roots
1
,σ,µ,ν
0.071
1.874 10
0.2
0.215
0
I = 0.35 0.4
0.5
0.55
0.65 0.7
r = 0.36
, σ , µ , ν , 1 .10
sol =
2.331 10
0.508
0
0.657
0
13
15
13
check the roots
Value of G at roots
0.05 0.1
0.25
1, 1
Appendix C
334
Define the expressions for the degrees of consolidation of layers I and II for case b.
i
12
INC
.1
Ucon 1 ( σ , µ , ν , T )
sum
0
zbrak( i , σ , µ , ν , INC )
W
for n ∈ 1 .. i
An
brent Wn
An
Ucon 2 ( σ , µ , ν , T )
W
2.
sum sum
2 .sum
sum
, Wn
1, 1
, σ , µ , ν , TOL
sin( An ) .cos( An ) .sin( µ .ν .An )
term
1
1,0
σ .sin( µ .ν .An )
2
µ .ν .cos( An )
.exp( T .An2 )
2
term
0
zbrak( i , σ , µ , ν , INC )
for n ∈ 1 .. i
An
term
sum
1
brent Wn
1,0
, Wn
2
cos( An ) .( 1
1, 1
cos( µ .ν .An ) )
2
An . σ .sin( µ .ν .An )
sum
, σ , µ , ν , TOL
2
µ .ν .cos( An ) 2
.exp( T .An2 )
term
2 .
sum
µ .ν
Define expression for "resultant" degree of consolidation for both layers.
U res( σ , µ , ν , T )
U1
Ucon 1 ( σ , µ , ν , T )
U2
Ucon 2 ( σ , µ , ν , T )
ν .U2
1
U1
ν
Verification Problem Two - consolidation of contiguous clay layers
with different permeability
335
Verification problem 2 for SAGE
p
100
σ
1
2
1
2
µ
ν
1
h1
5
e1 e2
First, plot G(A) to find out nature of function for the parameters of the given problem.
A
0 , .1 .. 10
1
G( A , σ , µ , ν )
0
1
0
5
A
10
INC = 0.1 is a reasonable increment for the function zbrak to bracket the roots of G(A) for this problem.
Note that INC must be small enough to "capture" each root of G(A).
Data p , σ , µ , ν , h 1 , .1 , 100
I
Data p , σ , µ , ν , h 1 , .5 , 100
1
0.8
0.6
z/H
H
0.4
0.2
0
0
0.2
T = 0.10
T = 0.50
0.4
0.6
Values of u/p
0.8
1
Appendix C
T trial( beg , end , n )
T val
k
beg
Define the function Ttrial to generate the time
factors, T, for a plot of resultant U versus the
logarithm of T.
10
0
1
for j ∈ beg .. end
num
10
1
beg is the beginning log cycle (i.e. 0 for 100 )
j
for jj ∈ 1 .. n
T val
k
k
336
num .
k
end is the ending log cycle (i.e. 1 for 101 )
n is the number of increments for T in each
log cycle
10 .
jj
n
1
Example:
1
T val
3.333
6.667
T trial( 0 , 2 , 3 ) = 10
33.333
66.667
l
0 , 1 .. 40
100
:range variable used for plot
Resultant U versus log T for
Verification Problem 2
1
0.8
U
0.6
0.4
0.2
0
1 10 3
0.01
0.1
T
1
10
9HULILFDWLRQSUREOHPWKUHH²VXUIDFHVHWWOHPHQWRIFOD\OD\HUFRQVROLGDWLQJ XQGHUDVWULSIRRWLQJORDG
C.4. VERIFICATION PROBLEM THREE – SURFACE
SETTLEMENT OF CLAY LAYER CONSOLIDATING
UNDER A STRIP FOOTING LOAD
*LEVRQHWDOGHULYHGDQH[SUHVVLRQIRUWKHVHWWOHPHQWRIWKHVXUIDFHRID
FOD\OD\HUORDGHGZLWKDVWULSIRRWLQJ7KHH[SUHVVLRQGHYHORSHGE\*LEVRQHWDO
DQG LWV QXPHULFDO HYDOXDWLRQ DUH GHPRQVWUDWHG LQ WKH IROORZLQJ 0DWKFDG
ZRUNVKHHW*LEVRQHWDOVROXWLRQPFG
Appendix C
338
The purpose of this MathCAD worksheet is to evaluate the expression for the surface settlement
of a clay layer loaded with a strip footing. The original derivations were done by Gibson, R.E.,
Schiffman, R.L.,and Pu, S.L. "Plane Strain and Axially Symmetric Consolidation of a Clay Layer
on a Smooth Impervious Base," Quartly Journal of Mechanics and Applied Mathematics, Vol.
XXIII, Pt. 4, 1970, pp. 505-520. The expression that they developed is difficult to evaluate,
because it involves an semi-infinite integral and an infinite series. Furthermore, it also involves
finding the roots of a characteristic equation.
b
b
f
x,r
z
h
Gibson et al found the following expression for the settlement of the surface of the clay layer.
∞
η. f .
2. G
w o ( x, t )
2
Γ ( x, λ ) .
tanh ( λ ) . 2 . η. M ( λ )
2
2. η 1
λ
P( λ , t )
dλ
0
where:
∞
L( λ )
2
1
λ . η .( 1
λ .csch( λ ) .sech( λ ) )
1
λ
P( λ , t )
2
n=1
1
M ( λ ) λ .η .coth( λ ) .( 1
F α n, λ
1
M( λ )
1.
1
2
λ .csch( λ ) .sech( λ ) )
tan α n
2
1
1
α n .tan α n
αn are the roots of the characteristic equation:
α
2
2 ct
λ .
2
h
F α n, λ
b
x
2 .h .
sin λ . .cos λ . plane strain
.
h
h
πλ
b
r
Γ ( r , λ ) b . J 1 λ . .J 0 λ .
h
h
axisymmetric
material parameters:
ν = poisson's ratio
L( λ ) α .M ( λ ) .tan( α )
2
Γ ( x, λ )
αn
exp
G = shear modulus
k = permeability of the clay
η
1
1
ν
2 .ν
G
E
.
2 (1 ν)
c 2 .G .η .
k
γw
γw = unit weight of pore fluid (water)
c = coefficient of consolidation
Verification Problem Three - surface settlement of clay layer consolidating
under strip footing load
339
For the purposes of an example problem:
h
4
ν
0.30
k
b
4
E
33333
γw
η
1
(1
ν
2 .ν )
G
0.0028
η .( 1
2 .( 1
c
ν)
M( λ )
1.
1
2
1 .10
ε
k
2 .G .η .
γw
1 .10
12
1 .10
14
12
plane strain
case
λ .csch( λ ) .sech( λ ) )
1
F( α , λ )
TOL
ε mach
E
1
L( λ )
1000
62.4
b
x
2 .h .
sin λ . .cos λ .
h
h
π .λ
Γ ( x, λ )
f
tan( α )
2
1 .λ
M( λ )
η .( 1
λ .coth ( λ )
λ .csch( λ ) .sech( λ ) )
tan( α )
2
α
Determine the roots ,αn , of the characteristic equation α2 - L(λ) - α∗M(λ)*tan(α)
q( α , λ )
α
2
L( λ )
α .M ( λ ) .tan( α )
Define the function bisect which returns a vector containing a and b, where a and b define an
interval (a,b) for which f(a)*f(b)<0, using the Bisection Method. For this case f(x) is q(a,l)
bisect ( n , λ , φ , itmax )
a
(n
1 ) .π
b
n
1 .
π
2
iter
The initial estimate of
(a,b) is (nπ,n+1/2π)
0
while iter < itmax
xnew
iter
b
a
a
2
iter
Note: α1 occurs in the interval
(0,π/2), α2 occurs in the interval
(π,3π/2) and so on.
1
a
b
break if
fa
(a
(b
a
b
b
a
π
<φ
q( a , λ )
fnew
The function, bisect, will be used to
determine a very small interval (a,b)
which contains αn .
q ( xnew , λ )
xnew ) if ( fa .fnew > 0 )
xnew ) if ( fa .fnew 0 )
Appendix C
340
Define the function Muller which uses Muller's method to refine the estimate of αn , the root of q(α,λ).
Mulle r 0x, x 1 , x 2 , maxit , λ
for iter ∈ 1.. maxit
h0
x1
x0
h1
x2
x1
P0
q x 0, λ
P1
q x 1, λ
P2
q x 2, λ
P1
δ0
P0
h0
P2
δ1
P1
h1
a0
P2
a2
δ1
δ0
h1
h0
a1
a 2 .h 1
δ1
D
a1
2
4 .a 2 .a 0
E
a1
D if
a1
E
a1
D otherwise
h
2.
D > a1
D
a0
E
xstar
2x
x0
x1
x1
x2
x2
xstar
h
xstar
Reference: Asaithambi, N.S. Numerical Analysis Theory and Practice. Saunders College Publishing.
Verification Problem Three - surface settlement of clay layer consolidating
under strip footing load
341
Define the root finding function alpha. alpha finds the nth root of the characteristic equation: q(α,λ).
Alpha uses the bisection method to narrow the interval around the nth root and then calls Muller to
"polish" the root. Note: Muller's method finds real and complex roots.
alpha( n , λ )
x
bisect ( n , λ , ε , 42 )
Re Muller x0 ,
Example:
an
x0
x1
2
alpha( 2 , 5 )
, x1 , 1 , λ
an = 4.314963
term( n , λ , t )
Define the function,
term, which evaluates
the nth term of P(l,t).
αn
q ( an , 5 ) = 1.776357 10
14
alpha( n , λ )
num
denom
exp
αn
2
2 c .t
λ .
2
h
F α n, λ
num
denom
Define the function, H(λ,t), which estimates P(λ,t) by summing P(λ,t) until the terms computed
become insignificant.
H( λ , t )
H
0
i
350
si
0
10 .TOL
check
while check > TOL
i
i
1
term( i , λ , t )
si
check
H
H
si
H
if
H >0
si
break if i 1
H
Note that roundoff and truncation errors pose serious problems to the evaluation of P(λ,t) using H(λ,t).
Next, Euler's transformation of alternating series and van Wijngaarden's implementation will be
examined as an alternative way of approximating P(λ,t).
Appendix C
342
Define the function Euler, which evaluates Euler's transformation of alternating series with van
Wijngaarden's implementation. One term of the original alternating series are incorporated into the
estimates of the partial differences.
, wk1, wk2 , su m)
sub 1( nter m
if
wk2
wk1
newsum
n
nterm
otherwise
newsum
n
0.5. wk2
sum
Sub1 is a subroutine used by Euler.
Sub1 determines whether a new
partial difference should be
evaluated (increase nterm by 1) or
just revise the estimate of S.
1
sum
wk2
nterm
n
newsum
Eule (r ter ,m jter , wks
m p)
if jter m 1
n
Euler estimates
the summation of
a convergent
infinite series
whose terms
alternate in sign.
1
0.5. term
sum
wksp2
term
otherwise
n wksp0
tmp
wksp2
wksp2
term
for j ∈ 1.. n
1
dum wksp j
wkspj
tmp
wkspn
2
if n 2
2
0.5 . wkspj
1
tmp
dum
2
0.5 . wkspn
a
sub 1 n, wksp
n
n
a0
sum
Wijngaarden's
implementation
adapts Euler's
transformation to
positive or negative
convergent series.
1
1
tmp
, wkspn
2
, wksp1
a1
wksp0
n
wksp1
sum
wksp
Reference: Press, W.H., Teukolosky, S. A., Vetterling, W.T., and Flannery, B.P. Numerical Recipes in
FORTRAN. 2nd ed. Cambridge University Press.
Verification Problem Three - surface settlement of clay layer consolidating
under strip footing load
∞
Set up a trial series to evaluate the function Euler. Solve the problem:
A
( k)
n
n=1
trial( k )
j
for k<1:
0
s
k
A
1
0
W
k
0
check
10 .TOL
W old
0
The function trial uses the Euler
function to estimate A.
while check > TOL
j
j
1
s
0
for t ∈ 8 , 7 .. 0
s
( 1)
s
t
2 .k
s
1.
j
t
2 .j
Use van Wijngaarden's procedure
for evaluating a positive series with
Euler's technique (i.e. convert the
series into an alternating series)
s
Euler( s , j , W )
W
W1
check
W old
W old
W old
if
W old > 0
W1
W1
j
Asum( k , n )
s
Asum evaluates A through n terms.
0
for j ∈ 1 .. n
s
Asum
s
j
k
1
, 22 = 0.333333333333314
4
trial
0.33333333333335
1
=
4
22
A
1
3
for k = 1/4
This verifies that Euler works, but it does not offer significant savings in computations in
this particular case.
343
Appendix C
344
Define J(λ,t), which will estimate P(λ,t) using the Euler's transformation. Van Wijngaarden's
transformation is used to convert P(λ,t) into an alternating series (terms in the sum alternate in
sign). Euler's transformation only works for alternating series.
J( λ , t )
W 0
check
10. TOL
jterm
0
Σ
0
r
0
while check > TOL
r
r
1
jterm
wr
jterm
1
0
for k ∈ 12, 11 .. 0
wr
( 1)
wr
k
k
2 .ter m 2 .r , λ , t
wr
jterm
1.
wr
,W
W Eule r wr , jter m
check
Σ
Σ
W1
Σ
if
Σ >0
W1
Σ
r
Examples of estimating P(λ,t) with J(λ,t) and H(λ,t):
J( 10, 0) =
2.284385
H( 10 , 0 ) = 2.266755
J( .5, 0) =
23
J( 10, 1) =
2.36028 10
0.13361
H( .5 , 0 ) = 0.133549
24
7
H( 10 , 1 ) = 2.36028 10
7
7
J( .5, 1) =
0.087142
H( .5 , 1 ) = 0.087142
6
It is apparent that the functions J and H give slightly different results at times early in the solution,
but there is almost no difference in the values later in the solution. The source of this difference is
probably round-off error. The round-off is probably larger for H since it is a simple summation.
Verification Problem Three - surface settlement of clay layer consolidating
under strip footing load
Define the function Gauleg which will return a matrix whose 1st column contains the Gauss points
and whose 2nd column contains the Gauss weights for the m-point quadrature rule for the interval
(a,b)
Gauleg( a , b , m )
j
m
2
3 .10
0.5 .( a
14
eps
xm
xl
1
b)
0.5 .( b
a)
for i ∈ 1 .. j
z
cos π .
i .25
m .5
2 .eps
check
while check > eps
p1
1
p2
0
for k ∈ 1 .. m
p3
p2
p2
p1
( 2 .k
p1
pp
m.
z1
z
z
z1
1, 0
Am
Ai
i,0
1, 1
Am
i,1
(k
1 ) .p3
z .p1 p2
z .z 1
p1
pp
check
Ai
1 ) .z .p2
k
z
z1
xm
xl .z
xm
xl .z
2.
((1
Ai
xl
.
z z ) .pp .pp )
1,1
A
Examples:
0.354062724
Gauleg( 1 , 1 , 2 ) =
0.5773502692 1
0.5773502692 1
0.872664626
Gauleg( 0 , π , 3 ) = 1.5707963268 1.3962634016
2.7875299296 0.872664626
345
Appendix C
346
Now, define a MathCAD function, w, to evaluate the settlement, w, at any (x,t). The variable
quad control how the integration is performed. Gauss-Legendre quadrature with a quad-point
rule is used. The integral is broken up into subintegrals each of length π. This is done because
the period of the function, w, is π.
w( x , t , qua d)
w
i
0
0
GAUS Gaule g( 0, π , qua d)
check
1
while check > 1 .10
i
i
a
(i
b
i .π
3
1
sum
1) .π
0
for j ∈ 1.. qua d
weight
λ
GAUSj
sum
w
w
GAUSj
1,0
1,1
a
su m Γ ( x , λ ) .
w
sum
check
sum
w
η .f .
w
2 .G
w( 0 , 0 , 5 ) = 0.079367
tan h( λ ) . 2 .η .M ( λ )
2
2 .η 1
λ
2
J( λ , t ) .weight
0
Verification Problem Three - surface settlement of clay layer consolidating
under strip footing load
Define two functions:
wfinal
- evaluates the final settlement at x
wimmed - evaluates the immediate settlement at x
w final( x , n , quad )
w
0
check
1
for i ∈ 1 .. 2 .n
a
(i
b
i.
1).
π
2
π
2
Gauleg( a , b , quad )
GAUS
sum
0
for j ∈ 1 .. quad
weight
λ
GAUSj
sum
w
w
w immed( x , n , quad )
GAUSj
sum
w
sum
check
sum
w
η .f
( 2 .η
1 ) .G
w final( x , n , quad ) .
1, 1
1, 0
Γ ( x , λ ) .tanh( λ )
.weight
λ .( 1 λ .csch( λ ) .sech( λ ) )
.w
2 .η
1
2 .η
For our example problem:
w final( 0 , 20 , 5 ) = 0.111079
G .w final( 0 , 20 , 5 )
w immed( 0 , 5 , 5 ) = 0.07976
G .w immed( 0 , 20 , 5 )
b .f
b .f
= 0.356018
= 0.254298
These results agree with Figure 5 from Gibson et al, 1970.
The following pages contain plots that investigate the nature of the functions involved in Gibson et al's
solution.
347
Appendix C
jj
..1 30
ll
348
1
0
Plot of the first 30 terms of the infinite series P(λ,t).
for λ = 1and t = 0 and 0.02
ter m( j ,j ,l0.02
l )
0.1
ter m( j ,j ,l0l)
0.2
0
10
20
30
jj
The effects of noise and roundoff error appear in
the higher values of αn (~n>10). λ = 1
0
q( alph (a ,j j) , l l)
5 10 9
1 10 8
0
ll
10
jj
20
30
Now set λ = 20 and replot the data.
20
0
J( l ,l )0 =
ter m( j ,j ,l0.02
l )
ter m( j ,j ,l0l)
4.566113
J( l ,l .02
)=
25
0.2
29
H( l l, 0) = 4.4956
0.4
0
10
20
30
jj
H( l l, .02) = 0.428654
Note the "pulse" type shape of the plot of the infinite
series P(λ,t). This "pulse" shape appears when λ is
greater than π.
5 10 10
0
ter m( j ,jπ , 0.02 )
q( alph (a ,j j) , l l)
0.428654
0
0.1
ter m( j ,jπ , 0 )
ter m( j ,jπ
.2 , 0 ) 0.2
5 10 10
0
10
jj
20
1
30
0.3
0
1
jj
2
1
Verification Problem Three - surface settlement of clay layer consolidating
under a strip footing load
ll
349
0 , .2 .. 20
Plot of integrand of wfinal vs. lambda
0.08
This plot clearly shows the damped
oscillating nature of the integrand of
Gibson et al's expression for
settlement.
0.06
0.04
0.02
0
0.02
0
ll
5
10
15
20
20 , 20.2 .. 40
Plot of integrand of wfinal vs. lambda
4 10 4
2 10
4
0
2 10
4
4 10 4
20
25
30
35
40
SAGE results are compared with this solution in Chapter 3.