Microcomputer calculation of the bearing capacity of circular footings on cohesionless soil using the method of characteristics

Microcomputer calculation of the bearing capacity of circular footings on cohesionless soil using the method of characteristics J. F. L Y N E S S a n d W. G I V E N Department o f Ovil Engineering, University o f Ulster at Jordanstown, Co. Antrim, Northern Ireland

J. G. S T U A R T Department o f Ovil Engineering, Queen's University, Belfast, Northern Ireland

A finite difference scheme based on the characteristic network representing the slip surface field under a non-tilting circular footing is used to give stress distributions and bearing capacities. Allowance is made for full friction between the footing and the soil, and account is taken of penetration o f the footing. Output from the microcomputer algorithm is compared with model test results and Terzaghi's formula.

computers are very suitable for the solution of slip surface fields using explicit schemes due to small storage requirements for the solution of field variables. Atkinson et al. s present a microcomputer based algorithm for the calculation of slip surfaces under strip footings.

ASSUMPTIONS Cohesionless rigid-plastic soil obeying the Mohr-Coulomb failure criterion, r f = a n tan¢

INTRODUCTION In engineering practice the bearing capacity of soils under circular footings is generally calculated from Terzaghi's formula. This formula is an empirically modified version of that derived by Terzaghi for long strip footings. Sokolovski 1 has given a numerical method for the calculation of stresses under strip footings. Cox et a l ) present a similar method of characteristics solution to calculate bearing capacities under circular footings. The soil is considered to be rigid-plastic and the equations of equilibrium within the slip field are combined with the Mohr-Coulomb failure criterion to give a set of hyperbolic differential equations. In a later paper, Cox 3 developed the solution further by allowing for weight of the soil. Cox et ai. showed that solution for stress or velocity field within the slip surfaces could be extended. By calculating both stress and velocity fields a kinematically admissible failure state or upper bound was shown to exist. This was also demonstrated to be the lower bound solution, by extending the stress field throughout the entire rigid region without violating the conditions of equilibrium and yield. Plane strain solutions taking account of footing roughness and variation of friction angle ~b, with stress level, have been given by Graham and Stuart. 4 The present algorithm extends that of Cox, allowing for soil/footing friction and penetration of the footing. The equations are written in terms of engineering stress parameters, and the resulting explicit numerical scheme was coded in BASIC for a 32K PET microcomputer. MicroAccepted June 1985. Discussion closes March 1986. 0141-1195/86/010025-07 $2.00 © 1986 Computational Mechanics Publications

(1)

o n = normal stress

¢ = angle of shearing resistance For the axisymmetric case the circumferential principal stress is taken as equal to the minor principal stress in the r - - z plane. This assumption is the same as that used by Meyerhof 6 and is in agreement with the hypothesis of Haar and yon Karman,7 which states that under axially symmetric conditions, the circumferential principal stress is equal to one of the other two principal stresses acting in an axial plane. The assumption neglecting all elastic strain before yield is not in agreement with observed test results in which settlements of the order of 0.1 times the diameter of the model footing were observed before failure. The boundary conditions for the numerical solution scheme have been formulated to account for significant geometry changes preceding failure.

DERIVATION OF THE CHARACTERISTIC EQUATIONS For the axisymmetric case the equations of equilibrium reduce to ~Or + ~Trz ar aTr z

--

Or

az

Or - - 0 0 + - -o r

~Oz

+ -- =7

az

r

+-

Trz

(2)

(3)

~, = unit weight of soil

Adv. Eng. Software, 1986, Vol. 8, No. 1

25

The stress components in the co-ordinate directions can be written as o r = a a (1

+ sin ~ cos 2rt)

(4)

o z = a a (1

-- sin ~bcos 2rt)

(5)

Zrz = a a

sin q~sin 2*/

o1 + 0 3

(6)

mean of principal stresses

oa - - -

2

*7 = angle between the major principal stress and the r-axis Substitution of (4), (5) and (6) into (2) and (3) gives two partial differential equations with o a and '7 as dependent variables and r and z as independent variables. Listing these two partial differential equations and the two equations of variation for Oa and r/in matrix form gives (1 + sinq~cos2r/),

sin¢sin2r/

, - - 2 o a sin C sin 27/,

F i g u r e 1. M o h r ' s b e t w e e n oa a n d c *

circle f o r

stress showing

1

OOa 20 a

relationship

sin ~ cos 2,1

--0 a

s i n ¢ - - ( 1 + cos2rt r

sinCsin2r/

,

(1 -- sin q~cos 2r/),

2Oasinq~cos2r/ ,

I 20 a

sin q~sin 2*7

oa

sin ¢) - sin 2r/+ 7

Oz I

dr

dz

,

0

-

-

dr

0

,

0

,

dr

r

I

drl I

dz

--

(7)

d%

i

d~

i

.dzJ The governing simultaneous partial differential equations can be shown to be hyperbolic by setting the determinant of the matrix above, equal to zero and calculating the roots of the resulting quadratic equation in d z / d r . The two real roots for dz/dr give the characteristic directions,

Following Cox, 3 auxiliary variables are introduced to simplify the solution of equations (9) and (10), as follows: h = settlement of footing c* = h " / t a n ¢

-dz -=tan

dr

( 77-- 4 + 2¢) and tan ( 7 / + -4/r - ~)

(11)

(8) A=cotq~ln(

Let

(12)

°asin¢]

\

c*

/

R = radius of footing

rr

a --~R and denote the characteristic directions c~ and 13. Then (Fig. 2), slope of a characteristics = tan qJ slope of/3 characteristics = tan

~ + - + 4~

(13)

C*

The c* term known as 'relative cohesion' will not alter the equations of equilibrium and yield. The relationship between c* and o a is shown in Fig. 1. The ordinary differential equations (9) and (10) can be arranged in the form" 1

If dsa and ds¢ are distances along the a and/3 characteristics, rearrangement of the matrix equation gives the following two ordinary differential equations which apply along the a and 15characteristic directions:

r

(9)

on ~ characteristics, and

(14)

dA -- 2dff = -- - (cos ~ dr -- ( 1 -- sin ¢) dz) r

G - - - - (sin ~bdr -- cos ¢ dz) e x p ( - - A tan q~)

sin ¢

R

(cos(C, + ¢)

+ sin if) dst~ -- 7 cos ~k dsa = 0 on 13characteristics.

26

G + - (sin q~dr + cos¢ dz) e x p ( - - A tan ~) R

1

+ sin ~b) dsa - - 7 sin(~ + ~) dsa = 0

cos$ d o a - - 2 a a sin $ d$

r

on a characteristics,

cos~ d o a + 2 o a sin~ d~ + o a sin ~b(cos(ff + ¢)

oa

dA + 2dqJ = - - (cos q~dr + (1 - sin ~) dz)

A d v . E n g . S o l ? w a r e , 1 9 8 6 , Vol. 8, N o . 1

(10)

(is)

on/3 characteristics. These ordinary differential equations (14) and (15), which apply along the characteristic directions, form the basis of a finite difference scheme.

CHARACTERISTICS NETWORK AND BOUNDARY CONDITIONS Figure 2 shows the general form of the characteristics network with three zones delineated, and the a and/3 characteristic directions being marked. Zone I directly under the footing provides boundary conditions for ~O. Sliplines for the family of characteristics are tangential to the underside of the footing so that the angle these characteristics make with the r axis is rr; therefore, ff = rr for 0 ~< r < R . On the top boundary of Zone III under the surcharge, h% the direction of the major principal stress is horizontal; therefore, r / = 0 and ~ = ¢ r / 4 - ¢/2 for r > R . Knowing oa in terms of sin ~ and c*, from Fig. 1, the boundary condition for A under the surcharge can be given as

FOOTING 1 / / / ,,/.,.,

for r > R . The radial zone of c~ characteristics, Zone II, emanates from the discontinuity between the footing and surcharge. This discontinuity is represented in the characteristic network as a degenerate ~ characteristic at point A. Boundary conditions so far have been defined in the regions 0 R leaving conditions at r = R to be defined. Figure 3 shows the radial zone of the network, it can be seen that ~k varies between 7r/4 -- ~b/2 in Zone III and lr in Zone I. Since point A is a degenerate/3 characteristic, equation (15) applies with increments dr and dz, equal to zero, giving (16)

If an angular increment co is considered between Zone I and a radial c~ characteristic then, dA = A2 -- A1

(17)

A2 = A1 + 2d~0 = A1 4- 260

(18)

From Zone III at the surface AI = cot4~ In cot

so.c,A.o~. ~

o(

A = cot ~ In cot (-n -- ¢--~ \4 2/

dA -- 2d~b = 0

SURFACE

r a t \ \ \ \ \ \ \ \ \ \ "

-

Figure 3. Radialzone of ~ characteristics at discontinuity between footing and surcharge

Therefore A2 = cot ~bIn cot

(4

--

+ 260

(19)

For any radial a characteristic at A (r = R ) , therefore, values of ~0 and A can be found, for each increment co, using equation (18) and (19). FINITE DIFFERENCE SCHEME The ordinary differential equations which apply along the characteristic directions can be written in finite difference form, and a two-step explicit iterative scheme can be devised. Referring to Fig. 4, the first step formula for calculation of the position o f p o i n t R , assuming the characteristics to be straight is,

¢_

I j

A

J

J

J

~5 ~ r

TI'T

D

/ C

Figure 2.

Zones of characteristic network representing slipline field

Adv. Eng. Software, 1986, Vol. 8, No. 1

27

iii

p(

./~ kINE ~

~

Equations (22) and (23) give the first step, while equations (24) and (25) provide the second step of the explicit scheme. Subsequent calculations o f r R and z R are improved by using (~1 + ~p)/2, (~1 + ~O)/2 in place of cp and ffQ. Improvement of calculation of A R and fir results fronl the use of ( A p + A I ) / 2 , (A O + A , ) / 2 in place or Ap and AQ ill the exponential terms of equations (24) and (25). These iterations are repeated until acceptable convergence is obtained. Referring to Fig. 5, along the underside of the footing, the boundary condition z = 0 and ~ = n for 0 < r < R simplifies the finite difference form of (24) and (25) so that, rn --

~I~,.O,(L.INE

~R,(r,,

z 0 + r O tan(rr/2 - - ~

~Q)

--

tan0r/2 -- ¢ -- CO)

Z,)

2 - (cos~(r R - r o ) (rR + r O)

AR = A o + 2 ( r t - - ~ O )

G + (1 -- sin ¢) Zo) - - R {sin @(rn - - r o) + z O cos ~}

Figure 4. Intersection o f characteristics showing the first approximation R1 to the true point o f intersection R

x exp(--AQ tan #) Z1

-

-

2p

tan ~ e - - rl - - r e

(20)

on a characteristic, 2 tan

-- ¢--

CQ) =

zl - - z q ro - r 1

(21) 27r

zO

-

-

Q-

zp + rp tan ~p + rQ tan(zr/2 -- ¢ -- ~Q)

(22)

tan ~kp + tan(rr/2 -- ¢ -- CQ) Zl = Zp + (rl --rp) tan ~p

(23)

Writing equations (9) and (10) in finite difference form along the a and t3 characteristic directions and re-arranging terms gives:

COMPUTER

A p + AQ + 2 ~ p - - 2qq2--

SOILC

STRESS

(r, + rp)

× {cos ¢P(rl -- rp) + (1 -- sin ¢) (Zl -- Ze)} 2 -

-

rl + rQ

o z r dr

rrR2

R

PROGRAMS

explicit solution scheme calculating rR, ZR, AR and fiR using equations (22), (23), (24) and (25). explicit solution for r R and A R directly under footing, using equations (26) and (27). evaluation of stress components Or, Oz and rrz from values of A and ~b.

FOOTINOB A S E

I

G -- {sin q~(rl --rQ) -- cos q}(zl --ZQ)} R

I"

~ R =T'f

"}

~r

O( x

~l =

exp(--AQ tan 4 0 ] / 2

A p - - A 1 + 2~bp

AR

Y

:tKNowN Y

(r, + rp)

G + R {sin ¢(rl -- rp) + cos ¢(Zl -- ze) }

Q

0(

x e x p ( - - A p tan @ ) ] / 2

Adv. Eng. Software, 1986, Vol. 8, No. 1

UNKNOWN

(24)

x {cos q~(rl --rp) + (1 -- sin ¢)(Zl --2,o)}

28

(28)

(cos ¢(rl -- ro) -- (1 -- sin ~) (zl -- ZQ)}

G + -- {sin ¢(rl -- rp) + cos #p(zl --Zp)} exp(--Ap tan ¢)

--

L

The finite difference scheme was implemented on a 32K byte PET microcomputer. Three programs were written,

BASEC A~ =

(27)

Once A and ff have been determined under the footing for the slip surface field, oa and rl can be determined from equations (8) and (12) and o z from equation (5). The bearing capacity, Q, of the soil under the circular footing is given by,

on/3 characteristic. Initial calculation for r and z gives rl =

(26)

(25)

^Q

Figure 5. Intersection o f / 3 characteristic with a characteristic under footing base

INPUT

J

t & ~t~,~rlA

IAT P P-Q'

STRESS28

(~ --33"

21,

I

22 20

I

CALCULATE INITIAL VALUES FOR ~"q , ~'¢ I

f

60% RELATIVE S0[L DENSITY

kN/m2 26

IB

n

16 I/*

"~--'~CALCULATE

AI~,/~

12 10

OUTPUT CURRENT

I

B

VALUES AT R FOR

E

.-f", "~c, A , ~

J

+

2

TEST CONVERGENCE I

O 0 005

0 010

0015

0"020

O 025

METRES ~

TERMINATE ITERATIONSI FOR THIS NODE I

- - I "~RECALCU"~-LATE&

,

SOIL

I

Figure 8 Slip line fieM for fully rough solution relative density 60%

Flow chart for program soil C

Figure 6. Table 1.

•

n FOOTJNG

Solution starts from boundary conditions for ~ and A along the upper boundary of Zone III. Proceeding through Zone III to Zone II the extreme/3 characteristic is taken as tangential to the z axis at r = 0. Figure 6 shows the flow chart for program SOILC.

Soil data from model tests

Relative density

"r(kN/m 3)

Settlement at failure (ram)

60%

7.26

33°

3.95

80%

9.68

35.4 °

4.75

90%

10.89

36.9 °

4.90

@

I

L OF FOUNDATION

20 k N/m s

1 • 01

'02

t

03

.

.

.

.

.

O~

-09

.10 ( m e t e r s)

'01

'02

03

(meters)

.04

Figure 7.

Slip distribution under model circular footing - method o f characteristics solution

Adv. Eng. Software, 1986, Vol. 8, No. 1

29

COMPARISON OF NUMERICALLY OBTAINED BEARING CAPACITIES WITH THOSE FROM EXPERIMENTAL RESULTS AND TERZAGHI'S FORMULA

Terzaghi's formula for a circular footing and factors used are given below:

qu = 7'DNq + 0.3B~/'N~

The ultimate bearing capacity values of a sand, at 60%, 80% and 90% relative densities were compared with experimental results. The experimental results come from model tests using a Rowe type cell containing submerged sands with a rough non-tilting circular footing of 25 mm radius. Data obtained from the model tests are given below for full roughness developed between soil and footing.

Nq

exp(37r/2

~') tan~'

:

2 cos20r/4 + 4;/2)

N7 = (Nq - 1) tan(1.44~')

(29)

(30) (31)

D = depth of footing B = diameter of footing

OF FOUNDATION 4

kN/rn z 0

-C -01

-0,2

I

"03

-OZ,

"05

'06

.07

"0

-0,9

.I0

-11 ( m e t e r s )

"02

,03

(meters)

,04

Figure 9.

Slip line field for fully rough solution relative density 80%

¢_ 6F FOUNDATION 60 50 ~z

k N~ m "Z 3O 20 D I S T R I lo / / / / /. ',meters)

•01 d

'02

'03

(meters)

,041

Figure 10.

30

Slip line field for fully rough solution relative density 90%

Adv. Eng. Software, 1986, Vol. 8, No. 1

80

density for experimental results and for both the numerical solution and Terzaghi's formula. It can be seen that Terzachi's formula and the numerical solution produce conservative values for bearing capacity. This may be due to the facts that (a) ~' in Terzaghi's formula was obtained from triaxial tests rather than plane strain tests, as suggested by Bishop,8 and (b) constant soil density was assumed in the plastic zones in the numerical solution.

E z

>~

50

~z

4o

~y

CONCLUDING REMARKS Finite difference schemes based on characteristic networks can be readily put into a form suitable for microcomputer computation. Other slip surface fields can also be accommodated by adjustment o f boundary conditions (for example the axisymmetric equivalent to the plane strain trapped wedge proposed by Terzaghi). Although not described in this paper, the graphical production of the slip surface field during computation is a feature that can be readily incorporated with the appropriate coding.

EXPERIMENTAL

g al .uJ

30

_~

L

i---

REFERENCES ,.._._.~ ~

o

60

~

r

70 RELATIVE

~

T H E OR ET I CA L CURVE [ T,E RZA 6HI )

80

I

90

100

DENSITY %

Figure 11. Comparison o f numerical solution with experimental and theoretical results

Figure 7 shows the vertical stress distribution at failure under the footing obtained with the aid of the computer programs for a sand at 60% relative density. Integration of the stress distribution over the area of the footing gives the ultimate bearing capacity for the soil. Figures 8-10 show the slip surface fields given by ot and ~ characteristics and base stress distributions, obtained for the three relative densities, using approximately 240 nodes. Figure 11 shows an increase in bearing capacity with increasing relative

1 Sokolovski, V. V. Statics of Soil Media, Butterworth Scientific Publications, London, 1960 2 Cox, A. D., Eason, G. L. and Hopkins, H. G. Axially-symmetric plastic deformation in soils, Phil. Trans. Royal Soc. Series A 1961, 254 3 Cox, A. D. Axially-symmetric plastic deformation in soils - II, Indentation of ponderable soils, Int. Journal of Mech. Sci. 1962, 4,371 4 Graham, J. and Stuart, J. G. Scale and boundary effects in foundation analysis, J. of Soil Mech. and Found., Div. ASCE 1971, 97 (SM11), 000, Proc. Paper 8510 5 Atkinson, J. H., Delpak, R. and Jones, G. J. Use of a microcomputer to determine slipqine fields for strip footings, Proc. 3rd Int. Conf. on Eng. Software, Imperial College, London, April 1983 6 Meyerhof, G. G. The ultimate bearing capacity of foundations, Gdotechnique 1951, 2 7 Haar, A. and yon Karman, J. Zur Theorie der Spannungszust~inde in Plastichen and Sandartigen Medien, Goettinger Nachr 1909 math-phys, kl 8 Bishop, A. W. The strengths of softs as engineering materials, Gdotechnique 1966, 16

Adv. Eng. Software, 1986, Vol. 8, No. 1

31