- Email: [email protected]
Three-dimensional elastic stress distribution around the flat end of a cylindrical cavity
Int. J. Rock Mech. Mitt. Sci. & Geomech. ,4bstr. Vol. 13. pp. 331 337. Pergamon Pre~ 1976. Printed in Great Britain
Three-Dimensional Elastic Stress Distribution around the Flat End of a Cylindrical Cavity G. H O C K I N G *
The three-dimensional elastic stress distribution around tile flat end of a cylindrical cat'ity has been preriottsly determined by experimental and numerical methods. Tile stress concentration factors determit~ed by pret~ious workers are not in good agreement, and tile effect of non-coincident principal stress direction with catrity axis has not been considered. Therefore, the stress distributions around the flat eml of a cylindrical carity for a range of Poisson's ratio .fi'om 0 to 0.475 hat,e been determined by the botuuklry integral equation method (B.I.E.M.). Close agreement for the stress concentration fiwtors a and c was Jbtmd between an experimental method and the B.I.E.M. To enable stress distributions to be calculated at particular points for principal stress.fiehls, non-coincident with tile carity axis, f o , r stress fiehls p~ = 1.0, p~ = 1.0, sy~ = 1.0 and s~ = 1.0 were applied at infinity. The cylindrical cat:ity was also soh:ed for fieht stresses p~ = 1.0 and p~ = py = Np~ applied at infinity with N = 0, 0.33, 1 and 2 and a Poisson's ratio of 0.25.
INTRODUCTION The stress and displacement distribution around the end of a cylindrical cavity in an infinite, homogeneous linear elastic medium is of interest to a number of problems:(1) Stress concentrations around the fiat end of a borehole for the determination of field stresses in rock. (2) Estimating the rock mass modulus from displacement monitoring as the face advances. (3) Determination of fracture or yield of the rock surrounding the excavation due to over-stressing. (4) Extent of the zone of influence around the end of a circular tunnel for stress field determination in boreholes driven from the face. A number of workers have investigated the stress concentrations aroung the flat end of a borehole by emperimental and numerical methods 1,1-7,1. The values of the stress concentrations around the fiat end of a borehole obtained by different workers are not in good agreement, and results have not been obtained for the case of borehole axis non-coincident with a principal stress direction. It was therefore considered necessary to determine the complete three-dimensional stress distribution around the fiat end of a cylindrical cavity by an independent method.
The problem posed is the linear elastostatic one. which is the determination of stresses and/or displacements at any point in or on the body due to imposed boundary conditions (whether traction, displacement or mixed). For convenience of explanation, suppose we are interested in the finite body problem of surface area. S, with applied boundary tractions and displacements (t and u, respectively) Fig. l(a). Now the essence of the B.I.E.M. is to formulate a relationship between u and t (i.e. a boundary constraint equation). To obtain this relationship we require a particular solution which involves tractions and displacements. In this case the simplest solution is that of Kelvin, the point load in an infinite mass. Discretize the surface S in Fig. 1 into n elements of which i and j are two. In the same body as in Fig. l(a), apply a unit force at i in the direction of 1, Fig. l(b). The tractions and displacements at all elements around the boundary can be calculated from Kelvin's solution 1-101. At a typical element, j, the effect of the unit force at i is to impose tractions T1. T2 and T3, U3 A T3
u3
! t3
t2
BOUNDARY I N T E G R A L EQUATION M E T H O D The integral equation method used has been described previously 1,8,9,1, however, a simplified explanation is presented here for completeness. ~) *Lecturer in Rock Mechanics. Imperial College of Science and Technology. London. U.K.
331
(b)
Fig. 1. The linear elastic body of surface S.
332
G. Hocking
and displaccmcnts UI, /_2 and /_'3. Fig. l/bi. From Betti's Reciprocal Work theorem [10] we obtain:
½uli + ,_~ (Tliul j + T2~u2j + T3:13j)
j=l
= ~ { U l j t l j + U2jt2i+ U3#3i)
j=l
i~j
(1)
which is simply the tractions of one system times the displacements of a second equal the tractions of the second system times the displacements of the first. The reason for ½uli is because the unit solution is divided equally between the interior and exterior region if S is smooth at i. Therefore by applying unit solutions in all 3 directions at each element in turn the following matrix identity is obtained:
I T ] [,,] = l-U] [t] where [T] and [U] are 3n x 3n
(2}
[u] and [t] are 3n x 1. Once the matrix identity, equation (2), has been formulated, then if either displacements (u) or tractions (t), or a combination of the two (u and t) are known, the system of equations can be solved. Cruse [l 1,12] has used constant and linear variation of u and t over each element and suggested that higher orders of variation would achieve more accurate results. On this suggestion, Lachat and Watson [-8], numerically modelled both geometry and boundary functions (u and t) by parametric representation. The surface is represented by 8 noded quadrilateral and 6 noded triangular elements (Fig. 2). The cartesian co-ordinates of an arbitrary point on an element are expressed in terms of the nodal co-ordinates and shape functions of the intrinsic co-ordinates. x~(~) = N"(~) x~" (i = 1, 2, 3)
(3)
where NI(~) = ¼(~, + 1)(,,."2 + 1)({t + 5, + 1) N,~(~) = ~(~.1 i ,- + 1)(1 -¢~.), -, etc.
{4)
The displacement u~ and Q may be considered to vary linearly, quadratically and cubically with respect to the intrinsic co-ordinates. The value 4) of a function at an arbitrary point on an element is
(5)
4)(~) = m:(~)4):,
where 4): are the nodal values for the function 4), and M:(~), the shape function, which was chosen to be of
g 2 " ~ _ <. / / / ~ ~
'~
8 Fig. 2. Surfaceelements.
1
quadratic variation, scc equation {4l, I\~r the problem presented here. Numerical integration of the boundar~ functions is bv Gaussian quadrative formulae with weight function 1.0. The 2 cases where the load point lies (a) within the element under consideration, and (b) outside the element are considered separately because of the variation of the kernals {U and TI with distance from the load point. If the surface S is not smooth at the load point, e.g. at corners and edges, then the singularity is not equally divided between the interior and exterior regions. Fortunatelx this effect can be calculated implicitly by considering rigid body translation [12]. When the identity (2) is assembled, the coefficients are scaled to obtain a non-dimensional and therefore numerically stable system. The coefficients of [U] are multiplied by the Young's modulus (El and divided by the maximum distance between 2 elements (rm~0; given values of displacement are divided by r,.... and given values of traction are divided by E. The unknowns are arranged to be placed in [u]. the knowns in [t], and the columns of [U] and [T] are interchanged where necessary. The system of equations is reduced bv Gaussian elimination using a block solver [8]. After solution, the calculated values of u and t are multiplied by rm= and E. respectively. Boundary stresses are calculated from surface traction and tangential strain, as mentioned by Cruse [-12]. Displacements at internal points are found by applying equation (1} with i. being the internal point in question. For the calculation of stresses at internal points, the equations for displacement at i are differentiated to give the strain components from which stresses can be deduced. Watson [13] has re-written the three-dimensional B.I.E.M. program of Lachat and Watson [8] and modified it to handle infinite body problems with field stresses applied at infinity. The modified program has been checked for the problems of a spherical cavity in a uniaxial stress field, and a spheroidal cavity in a triaxial stress field [14]. Stresses determined for the spherical and spheroidal cavities are within 1.5% of those given by Neuber's theoretical solutions [15].
C Y L I N D R I C A L CAVITY PROBLEM The B.I.E.M. program has been used to calculate the elastic stress distributions around a cylindrical cavity of unit radius and 10 diameters long (Fig. 3) for a variety of stress fields. The boundary element configuration is shown in Fig. 4; because the problem has symmetry with respect to all 3 axes. only one eighth of the geometry need be represented. To illustrate that stress distributions determined for this problem can be applied to the case of an infinitely long tunnel, the cavity was subjected to field stresses p: = I and Px = P: = Np: (Fig. 3) applied at infinity parallel to the co-ordinate axes with Poisson's ratio, v = 0.25. Boundary stresses in the y-z plane at point
Three Dimensional Elastic Stress Distribution
333
Sz.,~ x . '
1/
. z7" 4 I
~ v O ' S -~o..s.' X PLAN
Fig. 3. Geometry of the cylindrical cavity problem.
0 for field strcsscs of N = 0. 0.33, 1, and 2 are shown in Fig. 5. It is important to note that all 4 of these distributions correspond with those given by plane strain analysis [15]. STRESS C O N C E N T R A T I O N FACTORS
case II: pr = p.- = 0 p:, = 1.0, The factor c is the stress concentration in the direction normal to the uniaxial principal stress due to the presence of the flat borehole cnd. Thus. due to radial symmetry, the stress concentration factor c at point A (Fig. 3) is defined as:
If the axis of the cylindrical cavity is coincident with a principal stress direction, then stress concentration values, a, b and c are defined as follows: case I: p,~
= p y = 0 p : = 1.0.
The factor a is the stress concentration in the direction of the uniaxial applied principal strcss due to the presence of the flat borehole end. Since a_.: at point A (Fig. 3) acts in the direction of p= then a--
O'.~
p__' The factor b is the stress concentration in the direction normal to the uniaxial applied principal stress due to the presence of the flat borehole end. Since ayy at point A (Fig. 3) acts in a direction normal to p_- then b --
O'yy
Gyy -- Gz:
C --
P.~
Px
For the purposes of comparison, the stress concentration factors a, h and c were determined by the B.I.E.M. program for values of Poisson's ratio from 0 to 0.475. A summary of stress concentrations determined by previous workers is presented in Table 1. In physical modelling the difficulties associated with the finite size of the model can lead to significant errors (greater than 2 0 ° ) [15], and the resulting solutions have to be regarded as approximate. The axisymmetric finite element solution [7] is also approximate for this problem, because the determination of stress concentration factors a, b and c is a truly three-dimensional problem.
P:
(2"
¢
z 5.0 40 3.0 Y 20 1.0 0
0 Fig. 4. Boundary element configuration of the cylindrical cavity.
1.0
20
3.0
40
50
-,,'-(7
Fig. 5. Tangential stresses at section 10 radii from the flat end for N = O. 0.33. 1 and 2.
334
G. ttocking T~,BLL 1, SLMMAR'~ ()k SI RI-SS ( - ( I N I ' < T R . \ T I ) x .
Reference
Method
Galle & W i l h i o t [ I ] Lceman [2]
Photocla>ticit 3 Photochtqidt~ Experimental Experimental Experimental Experimental Photoelasticit.',
Hoskins [3] Bonncchere [4] Van Heerden [5]
FA('T()RS Ol I I R M I N I . D I'h PRI \ [()l N \\()RKIR5;
Poisson's ratio
,~
1,
,
O.47 0.4,~ O.29 0.26 0.22 0.38 0.48 0.35 0.30 0,26 0.44 0.37 0.29 0.24 11.4 0.3 0.2 0.0 0.475 0.4 0.3 0.25 0.2 0.1 0.0
1.56 ].55 1.51 1.53 1,56 1.25 1.25 1.28 1.24 1.22 1.42 1.39 1.32 I, 36 1.45 1,41 1.40 1.36 1.39 1.38 1.36 1.36 1.35 1.34 1.33
0 0 () 0 0 0 - 0.071 -0.022 - 0.065 -0.098 - 0.060 -0.204 -0.289 - 0.304 1).0 -11.1)4 - 0.08 -0.12 0.08 0.03 - 0.03 - I).I)5 - 0.07 -0.10 -0.13
- 0.85 -0.74 -0.71 -0.67 - 1, I 0 -0.98 -O.82 - 0.69 -0.91 -0.84 - 0.75 -0.52 - 0.91 -0.82 - 0.70 - 0.64 - 0.58 -0.48 -0.37
Ex pc rimcn tat1 Experirrlerltal Experimcntal
Experimental
Hiramatsu & Oka [6]
Experimental
Experimental Experimental A x i s y m m e t r i c finite
Coates & Yu [ 7 ]
elements
Axisymmetric Axisymmetric Axisymmetric Hocking (present results)
3-D 3-D 3-D 3-D 3-D 3-D 3-D
B.I.E.M. B.I.E.M. B.I.E.M. B.I.E.M. B.I.E.M. B.I.E.M. B.1.E.M
The stress concentration factors a, b and c determined by the B.I.E.M. program are approximately constant over a quarter of the flat end of the cylindrical cavity. Stress concentration factor c shows the greatest variation, approx 10%, within a radius 0.5 from point A. STRESS FIELDS N = 0, 0.33, 1 A N D 2 The B.I.E.M. program has been used to calculate the elastic stresses and displacements around the cylindrical cavity under the influence of 4 field stress conditions. Field stresses of p= = 1.0 and p., = p~. = Np= were applied at infinity in directions parallel to the coordinate axes (Fig. 3). Calculations were carried out for N = 0, 0.33, 1 and 2 and v = 0.25. Boundary stresses along the centre line of the cylindrical cavity are shown in Fig. 6. Along the crown of the cavity, a.,x approaches the plane strain values after
Stress
coi1CCiltr:_!.tion fach~r~
4.0 radii. The boundary stresses at the sharp corner will be in error, due to the lack of concentration of elements at this corner, but the stresses calculated at the mid-side node (located 0.25 radii from the corner) will be correct. It is interesting to note the approximately constant value of stress at the flat end of the cavity for a distance of 0.5 from point A. Tangential stresses are plotted in Fig. 7 [a) and (b) for sections perpendicular to the cavity axis at 3.5 and 0.75 radii from the flat end of the cavity. In Fig. 7(a) the tangential stresses are approximately the same as for the plane strain case. In Fig. 7(b) the tangential stresses are slightly reduced from those in Fig. 7(a), indicating that the effect of the flat end on tangential stresses 0.75 radii from the cavity end is small (the maximum variation from plane strain is 17°o). Boundary stresses at points A, B and C (Fig. 3) located at the flat end of the cavity are presented in Table 2. An interesting point is that the stress fields o4.0
. . . .
TENSION
/ I/
CC~PRESSION -N=20
3.0 2.0 1.0
.
" "~;.0" CT " -
• • 2.0 10
0
80
7.0
6.0
50
40
3.0
2.0
10
X "~
SECTION
AA
Fig. 6. B o u n d a r y stresses along ca'dry centre-line for N = 0, 0.33. 1 a n d 2.
0
Three Dimensional Elastic Stress Distribution O" 5.0
335
(7" --~---"-----. ~..._.
-
-
-
-
-
-
-
-
...... ..'.
TENSIOr'~
5.0
COMFRESS'ON
4.0
4.0
3,0
3.O
20
2.O
1.0
1.0
0
0
0
1.0
2.0
3.0
4.0
5.0
~'~-~-.
.-1..o0
(ca) SECTION BB
COMPIRESS;CN
7' x 1.0
2.0
3.0
4.0
' :'
--~O"
50
(b) SECTION CC
Fig. 7. Tangential stresses at sections (a) 3.5, and (b) 0.75 radii from the flat end for N = 0. 0.33. 1 and 2.
for N = 0 and N = 2 yield approximately the same compressive stresses for all 3 points A, B and C. The uniformity of the. boundary stresses in the region of points A, B and C is worth highlighting with respect to placing a strain rosette on this face. The extent of the zone of influence of the cavity measured along the cavity axis from the flat end into the solid, is presented in Table 3 for stress fields N = 0. 0.33, 1 and 2 with v = 0.25. To define a zone of influence the following rule is applied: a point lies within the zone of influence if the stresses acting at that point depart from the stress fields applied at infinity by 0.05 or more. The displacements in the crown of the cylindrical cavity in the z direction for p= = 100, N = 0, 0.33, 1 and 2 with Young's Modulus E = 21,000 and v = 0.25 are presented in Fig. 8. Eighty % of the vertical displacement difference between the sharp corner D and a point E remote from the fiat end occurs within 0.75 radii from the flat end. CAVITY AXIS N O N - C O I N C I D E N T W I T H P R I N C I P A L STRESS DIRECTION Since the cylindrical cavity is axisymmetric, then by applying the 4 unit stress fields p:,, p.., %.. and s.-x (Fig. 3), it is possible to determine the stress distributions for any orientation of principal stresses with respect
to the geometrical axes Oxyz of the cylindrical cavity. The boundary element configuration (Fig. 4) was altered when non-symmetric loading was imposed. If the principal stress directions are defined by rectangular axes Ox'y'z' which have direction cosines relative to the geometrical axes of the cylindrical cavity Oxyz of (11, ml, nt), (12, m2, t12), (13- m3. 113). respectively. then [16]
[,j :
.s...~ s.,y
= J
"T
n~
| 2n~lt L2/~ml
n3
2n212
211313
212m2
2/3nl 3
Stress
0
A
art a_._. are. a:. ayr a....
- 0.05 1.36 - 0.04 1.38 - 0.06 1.36
B C
Negative values denote tensile stresses.
0 0
"
0 0 0 0
where Px', Py' and p_.' are the principal stresses acting in the directions Ox' Oy' and 0z'. respectively. The stresses due to the four unit stress fields Px, P_-, %.. and s:., for points A, B and C (Fig. 31 are presented in Table 4: however, for the sake of simplicity, the cases of stress field Pr and s.,.r which are the same as p.. and s:x are included. Therefore by the use of the three-dimensional transformation matrix and superposition, the stresses at points A. B and C can be determined for any magnitude and orientation of stress field applied to the cylindrical cavity.
TABLE 2. BOtL~DaRY STRFSSES A t POINTS A, B AND C FOR .'kr = 0. 0.33. 1 AYD 2
Point
300
Stress f i e l d - - v a l u e of N 0.33 l
0.19 1.13 0.20 1.13 O. 16 1.13
0.67 0,67 0.68 0.61 0.60 0.66
2
1.38 - 0.02 1.39 -0.17 1.26 - 0.04
( i flocking
336 IXBIt
3
ZI)NFI1F INFL~IX(I
t{)R
\
=
1.~3
] \Ni} ~
Stres~ lield
Z o n e o1 influence m e a s u r e d a l o n g ca~ it~ axis from Ilat end into solid
(1
2.0
I
2.5
T A B I , E 4, STRESSES AT P{)INI'S A . B AND C FOR U N I T SI'Rt(SS H E I I ) S
APPLIEI) AF INFINITY
w i n 1 v = 0.25
Point A
1',-= I
&. = 1
0.64 0.64 0 0.65 (/.72 0 0.72 0.67 0
- 1.36 -0.05 0 1.36 -0.06 0 1.38 - 0.04 0
o',. &. r~: ay,. ~r._ r.~: a~,~ a:: ry:
B
C
Stress field at infinit) p: = I s,: = I s:., = 1 -0.05 1.36 0 -(I.(14 1.38 0 -0.06 1.36 (1
0 0 1.41 0 0 1.42 0 0 1.42
0 0 0 0.05 -0.37 0 (1 O -0.19
~,-
1
0 I) 0 (1 0 - O. 19 -0.37 0.05 (1
N e g a t i v e values d e n o t e tensile stresses.
i
0
I dz 0
I i
0.01
' " N : 0.33
dz 0 ! t
i
"~"N =2.0
I 0
9.0
8D
Q01
dz o
I
7.0
60
O01
[
I
5.0
/-,.0
I 3.0
2.0
E 1.0
z 0
X,f
SECTION
AA
Fig. 8. B o u n d a r y d i s p l a c e m e n t a l o n g cavity centre-line in the z direction for N = 0. 0.33, 1 a n d 2.
CONCLUSIONS The three-dimensional stress distributions around the flat end of a cylindrical cavity in an infinite, homogeneous, isotropic, linear elastic medium have been determined by the BTE.M. for a variety of stress fields applied at infinity. Stress concentration factors a, b and c due to radial and axial field stresses were determined and found to be in close agreement with previous experimental results. The stress concentration factors were found to be approximately constant over a concentric circular region of the flat end of radius half that of the cavity. The stress concentration factor c showed the maxitnum variation, approximately lff'i, within this region.
The stress distributions due to the stress fields. N = 0, 0.33, l and 2. indicate that departure from plane strain conditions close to the flat end of the cylindrical cavity is slight. For stress field determination using boreholes driven from the face of circular tunnels, the measuring device must be located at least 3.7 radii from the face to ensure that the influence of the tunnel on measurements will be small for most stress fields experienced in situ.
The stresses at points on the flat end of the cylindrical cavity have been determined for 4 stress fields applied at infinity. By applying superposition, the stresses at these points can be found for any orientation of principal stresses with respect to the geometrical axes of the cavity. From the vertical displacement profile. 80". of the vertical elastic displacement of the crown of a tunnel will occur within 0.75 radii from the face. An estimate of the rock mass modulus could be calculated from these displacement distributions, providing displacements can be monitored as the face axlvances. The BTE.M. is well suited to solving problems such as the one mentioned here. Thirty input cards were needed to formulate the cylindrical cavity problem, and computational requirements for solving this problem with 5 load cases was 30.6 s central processor time. 1.7 s input/output time, and 36k (decimal) small core memory of a CDC 7600 computer. AcknoMedgements--The w o r k described in this p a p e r flwms part of a n i n v e s t i g a t i o n into the a p p l i c a t i o n of the B . I . E M . to stress analysis of u n d e r g r o u n d o p e n i n g s being carried out at I m p e r i a l College, London. T h e a u t h o r ~ i s h e s to t h a n k Dr. J. O. W a t s o n for access to the t h r e e - d i m e n s i o n a l p r o g r a m , and for help in m a k i n g ,.he p r o g r a m operative, and Dr. E. T. Brm~n R~r ad',ice and g u i d a n c e in prepa r a t i o n of the m a n u s c r i p t .
Rvccircd 18 .4U~lU~t 1976.
Three Dimensional Elastic Stress Distribution REFERENCES 1. Galle E. M. & Wilhiot J. Stresses around a well bore due to internal pressure and unequal geostatic stresses. J. Soc. Petrol. En 9. (AIME} 2. 145-155 (1962). 2. Leeman E. The measurement of stress in rock--I. II and 1I!. JI. S. Aft. Inst. Min. Metall. 65, 45--114, 254-284 (1964). 3. Hoskins E. An investigation of strain-relief methods of measuring rock stress, int. J. Rock Mech. Min. Sci. 4, 155-164 (1967). 4. Bonnechere F. & Fairhurst C. Determination of the regional stress field from doorstopper measurements. Jl. S..4/?. Inst. Min. Metall. 69. 520-544 (1968). 5. Van Heerden E. W. Stress concentration factors for the flat borehole end for use in rock stress measurements. Enqm3 Geol. 3, 307-323 (1969). 6. Hiramatsu Y. & Oka Y. Determination of the stress in rock unaffected by boreholes or drifts from measured strains or deformations. Int. d. Rock Mech. Min. Sci. 5, 337-353 (1968). 7. Coates D. F. & Yu Y. S. A note on the stress concentrations at the end of a cylindrical hole. Int. J. Rock Mech. Mi~, Sci. 7, 585-588 (1970). 8. Lachat J. C. & Watson J. O. A second generation boundary integral equation program for three dioaensional elastic analysis. Boundary Integral Eq,uttion Method: Computational Applications
R.~l.~.s. 13 12
t
9.
10. 11. 12. 13. 14.
15. 16. 17.
337
in Applied Mechanics, (Edited by Cruse T. A. & Rizzo F. J.) pp. 85 100. ASME, NY (1975). Lachat J. C. A further development of the boundary integral technique for elastostatics. Ph.D. Thesis, University of Southampton, UK., 1975. Love A. E. H. A Treatise on the Mathematical Theory of Elasticity 4th edn., Cambridge University Press, Cambridge (1927}. Cruse T. A. Numerical solutions in three-dimensional elastostatics. Int. J. Solids Structures 5, 1259-1274 (1969). Cruse T. A. An improved boundary integral equation method for three-dimensional elastic stress analysis. Comput. & Struct. 4, 741-754 (1974). Watson J. O. Personal communication. Hocking G., Brown E. T. & Watson J. O. Three-dimensional elastic stress analysis of underground openings by the boundary integral eqnation method. 3rd Symp. Engineering Applications of Solid Mechanics, Toronto (1976)(in press}. Terzaghi K. & Richart F. E. Stresses in rock around cavities. Geotechnique 3, 57 90 0952). Riley W. F. Stresses at tunnel intersections. J. Enyny Mech. Div. Am. Soc. cir. Engrs 90. 167-179 (1964). Jaeger J. C. & Cook N. G. W. Fundamentals of Rock Mechanics. Chapman & Hall, London (1971).
Three-Dimensional Elastic Stress Distribution around the Flat End of a Cylindrical Cavity G. H O C K I N G *
The three-dimensional elastic stress distribution around tile flat end of a cylindrical cat'ity has been preriottsly determined by experimental and numerical methods. Tile stress concentration factors determit~ed by pret~ious workers are not in good agreement, and tile effect of non-coincident principal stress direction with catrity axis has not been considered. Therefore, the stress distributions around the flat eml of a cylindrical carity for a range of Poisson's ratio .fi'om 0 to 0.475 hat,e been determined by the botuuklry integral equation method (B.I.E.M.). Close agreement for the stress concentration fiwtors a and c was Jbtmd between an experimental method and the B.I.E.M. To enable stress distributions to be calculated at particular points for principal stress.fiehls, non-coincident with tile carity axis, f o , r stress fiehls p~ = 1.0, p~ = 1.0, sy~ = 1.0 and s~ = 1.0 were applied at infinity. The cylindrical cat:ity was also soh:ed for fieht stresses p~ = 1.0 and p~ = py = Np~ applied at infinity with N = 0, 0.33, 1 and 2 and a Poisson's ratio of 0.25.
INTRODUCTION The stress and displacement distribution around the end of a cylindrical cavity in an infinite, homogeneous linear elastic medium is of interest to a number of problems:(1) Stress concentrations around the fiat end of a borehole for the determination of field stresses in rock. (2) Estimating the rock mass modulus from displacement monitoring as the face advances. (3) Determination of fracture or yield of the rock surrounding the excavation due to over-stressing. (4) Extent of the zone of influence around the end of a circular tunnel for stress field determination in boreholes driven from the face. A number of workers have investigated the stress concentrations aroung the flat end of a borehole by emperimental and numerical methods 1,1-7,1. The values of the stress concentrations around the fiat end of a borehole obtained by different workers are not in good agreement, and results have not been obtained for the case of borehole axis non-coincident with a principal stress direction. It was therefore considered necessary to determine the complete three-dimensional stress distribution around the fiat end of a cylindrical cavity by an independent method.
The problem posed is the linear elastostatic one. which is the determination of stresses and/or displacements at any point in or on the body due to imposed boundary conditions (whether traction, displacement or mixed). For convenience of explanation, suppose we are interested in the finite body problem of surface area. S, with applied boundary tractions and displacements (t and u, respectively) Fig. l(a). Now the essence of the B.I.E.M. is to formulate a relationship between u and t (i.e. a boundary constraint equation). To obtain this relationship we require a particular solution which involves tractions and displacements. In this case the simplest solution is that of Kelvin, the point load in an infinite mass. Discretize the surface S in Fig. 1 into n elements of which i and j are two. In the same body as in Fig. l(a), apply a unit force at i in the direction of 1, Fig. l(b). The tractions and displacements at all elements around the boundary can be calculated from Kelvin's solution 1-101. At a typical element, j, the effect of the unit force at i is to impose tractions T1. T2 and T3, U3 A T3
u3
! t3
t2
BOUNDARY I N T E G R A L EQUATION M E T H O D The integral equation method used has been described previously 1,8,9,1, however, a simplified explanation is presented here for completeness. ~) *Lecturer in Rock Mechanics. Imperial College of Science and Technology. London. U.K.
331
(b)
Fig. 1. The linear elastic body of surface S.
332
G. Hocking
and displaccmcnts UI, /_2 and /_'3. Fig. l/bi. From Betti's Reciprocal Work theorem [10] we obtain:
½uli + ,_~ (Tliul j + T2~u2j + T3:13j)
j=l
= ~ { U l j t l j + U2jt2i+ U3#3i)
j=l
i~j
(1)
which is simply the tractions of one system times the displacements of a second equal the tractions of the second system times the displacements of the first. The reason for ½uli is because the unit solution is divided equally between the interior and exterior region if S is smooth at i. Therefore by applying unit solutions in all 3 directions at each element in turn the following matrix identity is obtained:
I T ] [,,] = l-U] [t] where [T] and [U] are 3n x 3n
(2}
[u] and [t] are 3n x 1. Once the matrix identity, equation (2), has been formulated, then if either displacements (u) or tractions (t), or a combination of the two (u and t) are known, the system of equations can be solved. Cruse [l 1,12] has used constant and linear variation of u and t over each element and suggested that higher orders of variation would achieve more accurate results. On this suggestion, Lachat and Watson [-8], numerically modelled both geometry and boundary functions (u and t) by parametric representation. The surface is represented by 8 noded quadrilateral and 6 noded triangular elements (Fig. 2). The cartesian co-ordinates of an arbitrary point on an element are expressed in terms of the nodal co-ordinates and shape functions of the intrinsic co-ordinates. x~(~) = N"(~) x~" (i = 1, 2, 3)
(3)
where NI(~) = ¼(~, + 1)(,,."2 + 1)({t + 5, + 1) N,~(~) = ~(~.1 i ,- + 1)(1 -¢~.), -, etc.
{4)
The displacement u~ and Q may be considered to vary linearly, quadratically and cubically with respect to the intrinsic co-ordinates. The value 4) of a function at an arbitrary point on an element is
(5)
4)(~) = m:(~)4):,
where 4): are the nodal values for the function 4), and M:(~), the shape function, which was chosen to be of
g 2 " ~ _ <. / / / ~ ~
'~
8 Fig. 2. Surfaceelements.
1
quadratic variation, scc equation {4l, I\~r the problem presented here. Numerical integration of the boundar~ functions is bv Gaussian quadrative formulae with weight function 1.0. The 2 cases where the load point lies (a) within the element under consideration, and (b) outside the element are considered separately because of the variation of the kernals {U and TI with distance from the load point. If the surface S is not smooth at the load point, e.g. at corners and edges, then the singularity is not equally divided between the interior and exterior regions. Fortunatelx this effect can be calculated implicitly by considering rigid body translation [12]. When the identity (2) is assembled, the coefficients are scaled to obtain a non-dimensional and therefore numerically stable system. The coefficients of [U] are multiplied by the Young's modulus (El and divided by the maximum distance between 2 elements (rm~0; given values of displacement are divided by r,.... and given values of traction are divided by E. The unknowns are arranged to be placed in [u]. the knowns in [t], and the columns of [U] and [T] are interchanged where necessary. The system of equations is reduced bv Gaussian elimination using a block solver [8]. After solution, the calculated values of u and t are multiplied by rm= and E. respectively. Boundary stresses are calculated from surface traction and tangential strain, as mentioned by Cruse [-12]. Displacements at internal points are found by applying equation (1} with i. being the internal point in question. For the calculation of stresses at internal points, the equations for displacement at i are differentiated to give the strain components from which stresses can be deduced. Watson [13] has re-written the three-dimensional B.I.E.M. program of Lachat and Watson [8] and modified it to handle infinite body problems with field stresses applied at infinity. The modified program has been checked for the problems of a spherical cavity in a uniaxial stress field, and a spheroidal cavity in a triaxial stress field [14]. Stresses determined for the spherical and spheroidal cavities are within 1.5% of those given by Neuber's theoretical solutions [15].
C Y L I N D R I C A L CAVITY PROBLEM The B.I.E.M. program has been used to calculate the elastic stress distributions around a cylindrical cavity of unit radius and 10 diameters long (Fig. 3) for a variety of stress fields. The boundary element configuration is shown in Fig. 4; because the problem has symmetry with respect to all 3 axes. only one eighth of the geometry need be represented. To illustrate that stress distributions determined for this problem can be applied to the case of an infinitely long tunnel, the cavity was subjected to field stresses p: = I and Px = P: = Np: (Fig. 3) applied at infinity parallel to the co-ordinate axes with Poisson's ratio, v = 0.25. Boundary stresses in the y-z plane at point
Three Dimensional Elastic Stress Distribution
333
Sz.,~ x . '
1/
. z7" 4 I
~ v O ' S -~o..s.' X PLAN
Fig. 3. Geometry of the cylindrical cavity problem.
0 for field strcsscs of N = 0. 0.33, 1, and 2 are shown in Fig. 5. It is important to note that all 4 of these distributions correspond with those given by plane strain analysis [15]. STRESS C O N C E N T R A T I O N FACTORS
case II: pr = p.- = 0 p:, = 1.0, The factor c is the stress concentration in the direction normal to the uniaxial principal stress due to the presence of the flat borehole cnd. Thus. due to radial symmetry, the stress concentration factor c at point A (Fig. 3) is defined as:
If the axis of the cylindrical cavity is coincident with a principal stress direction, then stress concentration values, a, b and c are defined as follows: case I: p,~
= p y = 0 p : = 1.0.
The factor a is the stress concentration in the direction of the uniaxial applied principal strcss due to the presence of the flat borehole end. Since a_.: at point A (Fig. 3) acts in the direction of p= then a--
O'.~
p__' The factor b is the stress concentration in the direction normal to the uniaxial applied principal stress due to the presence of the flat borehole end. Since ayy at point A (Fig. 3) acts in a direction normal to p_- then b --
O'yy
Gyy -- Gz:
C --
P.~
Px
For the purposes of comparison, the stress concentration factors a, h and c were determined by the B.I.E.M. program for values of Poisson's ratio from 0 to 0.475. A summary of stress concentrations determined by previous workers is presented in Table 1. In physical modelling the difficulties associated with the finite size of the model can lead to significant errors (greater than 2 0 ° ) [15], and the resulting solutions have to be regarded as approximate. The axisymmetric finite element solution [7] is also approximate for this problem, because the determination of stress concentration factors a, b and c is a truly three-dimensional problem.
P:
(2"
¢
z 5.0 40 3.0 Y 20 1.0 0
0 Fig. 4. Boundary element configuration of the cylindrical cavity.
1.0
20
3.0
40
50
-,,'-(7
Fig. 5. Tangential stresses at section 10 radii from the flat end for N = O. 0.33. 1 and 2.
334
G. ttocking T~,BLL 1, SLMMAR'~ ()k SI RI-SS ( - ( I N I ' < T R . \ T I ) x .
Reference
Method
Galle & W i l h i o t [ I ] Lceman [2]
Photocla>ticit 3 Photochtqidt~ Experimental Experimental Experimental Experimental Photoelasticit.',
Hoskins [3] Bonncchere [4] Van Heerden [5]
FA('T()RS Ol I I R M I N I . D I'h PRI \ [()l N \\()RKIR5;
Poisson's ratio
,~
1,
,
O.47 0.4,~ O.29 0.26 0.22 0.38 0.48 0.35 0.30 0,26 0.44 0.37 0.29 0.24 11.4 0.3 0.2 0.0 0.475 0.4 0.3 0.25 0.2 0.1 0.0
1.56 ].55 1.51 1.53 1,56 1.25 1.25 1.28 1.24 1.22 1.42 1.39 1.32 I, 36 1.45 1,41 1.40 1.36 1.39 1.38 1.36 1.36 1.35 1.34 1.33
0 0 () 0 0 0 - 0.071 -0.022 - 0.065 -0.098 - 0.060 -0.204 -0.289 - 0.304 1).0 -11.1)4 - 0.08 -0.12 0.08 0.03 - 0.03 - I).I)5 - 0.07 -0.10 -0.13
- 0.85 -0.74 -0.71 -0.67 - 1, I 0 -0.98 -O.82 - 0.69 -0.91 -0.84 - 0.75 -0.52 - 0.91 -0.82 - 0.70 - 0.64 - 0.58 -0.48 -0.37
Ex pc rimcn tat1 Experirrlerltal Experimcntal
Experimental
Hiramatsu & Oka [6]
Experimental
Experimental Experimental A x i s y m m e t r i c finite
Coates & Yu [ 7 ]
elements
Axisymmetric Axisymmetric Axisymmetric Hocking (present results)
3-D 3-D 3-D 3-D 3-D 3-D 3-D
B.I.E.M. B.I.E.M. B.I.E.M. B.I.E.M. B.I.E.M. B.I.E.M. B.1.E.M
The stress concentration factors a, b and c determined by the B.I.E.M. program are approximately constant over a quarter of the flat end of the cylindrical cavity. Stress concentration factor c shows the greatest variation, approx 10%, within a radius 0.5 from point A. STRESS FIELDS N = 0, 0.33, 1 A N D 2 The B.I.E.M. program has been used to calculate the elastic stresses and displacements around the cylindrical cavity under the influence of 4 field stress conditions. Field stresses of p= = 1.0 and p., = p~. = Np= were applied at infinity in directions parallel to the coordinate axes (Fig. 3). Calculations were carried out for N = 0, 0.33, 1 and 2 and v = 0.25. Boundary stresses along the centre line of the cylindrical cavity are shown in Fig. 6. Along the crown of the cavity, a.,x approaches the plane strain values after
Stress
coi1CCiltr:_!.tion fach~r~
4.0 radii. The boundary stresses at the sharp corner will be in error, due to the lack of concentration of elements at this corner, but the stresses calculated at the mid-side node (located 0.25 radii from the corner) will be correct. It is interesting to note the approximately constant value of stress at the flat end of the cavity for a distance of 0.5 from point A. Tangential stresses are plotted in Fig. 7 [a) and (b) for sections perpendicular to the cavity axis at 3.5 and 0.75 radii from the flat end of the cavity. In Fig. 7(a) the tangential stresses are approximately the same as for the plane strain case. In Fig. 7(b) the tangential stresses are slightly reduced from those in Fig. 7(a), indicating that the effect of the flat end on tangential stresses 0.75 radii from the cavity end is small (the maximum variation from plane strain is 17°o). Boundary stresses at points A, B and C (Fig. 3) located at the flat end of the cavity are presented in Table 2. An interesting point is that the stress fields o4.0
. . . .
TENSION
/ I/
CC~PRESSION -N=20
3.0 2.0 1.0
.
" "~;.0" CT " -
• • 2.0 10
0
80
7.0
6.0
50
40
3.0
2.0
10
X "~
SECTION
AA
Fig. 6. B o u n d a r y stresses along ca'dry centre-line for N = 0, 0.33. 1 a n d 2.
0
Three Dimensional Elastic Stress Distribution O" 5.0
335
(7" --~---"-----. ~..._.
-
-
-
-
-
-
-
-
...... ..'.
TENSIOr'~
5.0
COMFRESS'ON
4.0
4.0
3,0
3.O
20
2.O
1.0
1.0
0
0
0
1.0
2.0
3.0
4.0
5.0
~'~-~-.
.-1..o0
(ca) SECTION BB
COMPIRESS;CN
7' x 1.0
2.0
3.0
4.0
' :'
--~O"
50
(b) SECTION CC
Fig. 7. Tangential stresses at sections (a) 3.5, and (b) 0.75 radii from the flat end for N = 0. 0.33. 1 and 2.
for N = 0 and N = 2 yield approximately the same compressive stresses for all 3 points A, B and C. The uniformity of the. boundary stresses in the region of points A, B and C is worth highlighting with respect to placing a strain rosette on this face. The extent of the zone of influence of the cavity measured along the cavity axis from the flat end into the solid, is presented in Table 3 for stress fields N = 0. 0.33, 1 and 2 with v = 0.25. To define a zone of influence the following rule is applied: a point lies within the zone of influence if the stresses acting at that point depart from the stress fields applied at infinity by 0.05 or more. The displacements in the crown of the cylindrical cavity in the z direction for p= = 100, N = 0, 0.33, 1 and 2 with Young's Modulus E = 21,000 and v = 0.25 are presented in Fig. 8. Eighty % of the vertical displacement difference between the sharp corner D and a point E remote from the fiat end occurs within 0.75 radii from the flat end. CAVITY AXIS N O N - C O I N C I D E N T W I T H P R I N C I P A L STRESS DIRECTION Since the cylindrical cavity is axisymmetric, then by applying the 4 unit stress fields p:,, p.., %.. and s.-x (Fig. 3), it is possible to determine the stress distributions for any orientation of principal stresses with respect
to the geometrical axes Oxyz of the cylindrical cavity. The boundary element configuration (Fig. 4) was altered when non-symmetric loading was imposed. If the principal stress directions are defined by rectangular axes Ox'y'z' which have direction cosines relative to the geometrical axes of the cylindrical cavity Oxyz of (11, ml, nt), (12, m2, t12), (13- m3. 113). respectively. then [16]
[,j :
.s...~ s.,y
= J
"T
n~
| 2n~lt L2/~ml
n3
2n212
211313
212m2
2/3nl 3
Stress
0
A
art a_._. are. a:. ayr a....
- 0.05 1.36 - 0.04 1.38 - 0.06 1.36
B C
Negative values denote tensile stresses.
0 0
"
0 0 0 0
where Px', Py' and p_.' are the principal stresses acting in the directions Ox' Oy' and 0z'. respectively. The stresses due to the four unit stress fields Px, P_-, %.. and s:., for points A, B and C (Fig. 31 are presented in Table 4: however, for the sake of simplicity, the cases of stress field Pr and s.,.r which are the same as p.. and s:x are included. Therefore by the use of the three-dimensional transformation matrix and superposition, the stresses at points A. B and C can be determined for any magnitude and orientation of stress field applied to the cylindrical cavity.
TABLE 2. BOtL~DaRY STRFSSES A t POINTS A, B AND C FOR .'kr = 0. 0.33. 1 AYD 2
Point
300
Stress f i e l d - - v a l u e of N 0.33 l
0.19 1.13 0.20 1.13 O. 16 1.13
0.67 0,67 0.68 0.61 0.60 0.66
2
1.38 - 0.02 1.39 -0.17 1.26 - 0.04
( i flocking
336 IXBIt
3
ZI)NFI1F INFL~IX(I
t{)R
\
=
1.~3
] \Ni} ~
Stres~ lield
Z o n e o1 influence m e a s u r e d a l o n g ca~ it~ axis from Ilat end into solid
(1
2.0
I
2.5
T A B I , E 4, STRESSES AT P{)INI'S A . B AND C FOR U N I T SI'Rt(SS H E I I ) S
APPLIEI) AF INFINITY
w i n 1 v = 0.25
Point A
1',-= I
&. = 1
0.64 0.64 0 0.65 (/.72 0 0.72 0.67 0
- 1.36 -0.05 0 1.36 -0.06 0 1.38 - 0.04 0
o',. &. r~: ay,. ~r._ r.~: a~,~ a:: ry:
B
C
Stress field at infinit) p: = I s,: = I s:., = 1 -0.05 1.36 0 -(I.(14 1.38 0 -0.06 1.36 (1
0 0 1.41 0 0 1.42 0 0 1.42
0 0 0 0.05 -0.37 0 (1 O -0.19
~,-
1
0 I) 0 (1 0 - O. 19 -0.37 0.05 (1
N e g a t i v e values d e n o t e tensile stresses.
i
0
I dz 0
I i
0.01
' " N : 0.33
dz 0 ! t
i
"~"N =2.0
I 0
9.0
8D
Q01
dz o
I
7.0
60
O01
[
I
5.0
/-,.0
I 3.0
2.0
E 1.0
z 0
X,f
SECTION
AA
Fig. 8. B o u n d a r y d i s p l a c e m e n t a l o n g cavity centre-line in the z direction for N = 0. 0.33, 1 a n d 2.
CONCLUSIONS The three-dimensional stress distributions around the flat end of a cylindrical cavity in an infinite, homogeneous, isotropic, linear elastic medium have been determined by the BTE.M. for a variety of stress fields applied at infinity. Stress concentration factors a, b and c due to radial and axial field stresses were determined and found to be in close agreement with previous experimental results. The stress concentration factors were found to be approximately constant over a concentric circular region of the flat end of radius half that of the cavity. The stress concentration factor c showed the maxitnum variation, approximately lff'i, within this region.
The stress distributions due to the stress fields. N = 0, 0.33, l and 2. indicate that departure from plane strain conditions close to the flat end of the cylindrical cavity is slight. For stress field determination using boreholes driven from the face of circular tunnels, the measuring device must be located at least 3.7 radii from the face to ensure that the influence of the tunnel on measurements will be small for most stress fields experienced in situ.
The stresses at points on the flat end of the cylindrical cavity have been determined for 4 stress fields applied at infinity. By applying superposition, the stresses at these points can be found for any orientation of principal stresses with respect to the geometrical axes of the cavity. From the vertical displacement profile. 80". of the vertical elastic displacement of the crown of a tunnel will occur within 0.75 radii from the face. An estimate of the rock mass modulus could be calculated from these displacement distributions, providing displacements can be monitored as the face axlvances. The BTE.M. is well suited to solving problems such as the one mentioned here. Thirty input cards were needed to formulate the cylindrical cavity problem, and computational requirements for solving this problem with 5 load cases was 30.6 s central processor time. 1.7 s input/output time, and 36k (decimal) small core memory of a CDC 7600 computer. AcknoMedgements--The w o r k described in this p a p e r flwms part of a n i n v e s t i g a t i o n into the a p p l i c a t i o n of the B . I . E M . to stress analysis of u n d e r g r o u n d o p e n i n g s being carried out at I m p e r i a l College, London. T h e a u t h o r ~ i s h e s to t h a n k Dr. J. O. W a t s o n for access to the t h r e e - d i m e n s i o n a l p r o g r a m , and for help in m a k i n g ,.he p r o g r a m operative, and Dr. E. T. Brm~n R~r ad',ice and g u i d a n c e in prepa r a t i o n of the m a n u s c r i p t .
Rvccircd 18 .4U~lU~t 1976.
Three Dimensional Elastic Stress Distribution REFERENCES 1. Galle E. M. & Wilhiot J. Stresses around a well bore due to internal pressure and unequal geostatic stresses. J. Soc. Petrol. En 9. (AIME} 2. 145-155 (1962). 2. Leeman E. The measurement of stress in rock--I. II and 1I!. JI. S. Aft. Inst. Min. Metall. 65, 45--114, 254-284 (1964). 3. Hoskins E. An investigation of strain-relief methods of measuring rock stress, int. J. Rock Mech. Min. Sci. 4, 155-164 (1967). 4. Bonnechere F. & Fairhurst C. Determination of the regional stress field from doorstopper measurements. Jl. S..4/?. Inst. Min. Metall. 69. 520-544 (1968). 5. Van Heerden E. W. Stress concentration factors for the flat borehole end for use in rock stress measurements. Enqm3 Geol. 3, 307-323 (1969). 6. Hiramatsu Y. & Oka Y. Determination of the stress in rock unaffected by boreholes or drifts from measured strains or deformations. Int. d. Rock Mech. Min. Sci. 5, 337-353 (1968). 7. Coates D. F. & Yu Y. S. A note on the stress concentrations at the end of a cylindrical hole. Int. J. Rock Mech. Mi~, Sci. 7, 585-588 (1970). 8. Lachat J. C. & Watson J. O. A second generation boundary integral equation program for three dioaensional elastic analysis. Boundary Integral Eq,uttion Method: Computational Applications
R.~l.~.s. 13 12
t
9.
10. 11. 12. 13. 14.
15. 16. 17.
337
in Applied Mechanics, (Edited by Cruse T. A. & Rizzo F. J.) pp. 85 100. ASME, NY (1975). Lachat J. C. A further development of the boundary integral technique for elastostatics. Ph.D. Thesis, University of Southampton, UK., 1975. Love A. E. H. A Treatise on the Mathematical Theory of Elasticity 4th edn., Cambridge University Press, Cambridge (1927}. Cruse T. A. Numerical solutions in three-dimensional elastostatics. Int. J. Solids Structures 5, 1259-1274 (1969). Cruse T. A. An improved boundary integral equation method for three-dimensional elastic stress analysis. Comput. & Struct. 4, 741-754 (1974). Watson J. O. Personal communication. Hocking G., Brown E. T. & Watson J. O. Three-dimensional elastic stress analysis of underground openings by the boundary integral eqnation method. 3rd Symp. Engineering Applications of Solid Mechanics, Toronto (1976)(in press}. Terzaghi K. & Richart F. E. Stresses in rock around cavities. Geotechnique 3, 57 90 0952). Riley W. F. Stresses at tunnel intersections. J. Enyny Mech. Div. Am. Soc. cir. Engrs 90. 167-179 (1964). Jaeger J. C. & Cook N. G. W. Fundamentals of Rock Mechanics. Chapman & Hall, London (1971).