- Email: [email protected]
A numerical modelling criterion for the analysis of underground openings using infinite elements
A numerical modelling criterion for the analysis of underground openings using infinite elements Prabhat Kumar Structural Engineering Research Centre, Roorkee 247667, India (Received April 1985; revised October 1985)
Underground openings subject to an initial stress field are usually analysed by the finite element method (FEM) by truncating the analysis domain at a large but finite distance from the zone of disturbance. In this paper, the problem is analysed by the method of finite (FE) and infinite elements (IE) to represent near- and far-fields, respectively. Deep openings of circular and elliptical shape and shallow circular openings are analysed to establish the approach. The numerical results of stress concentration factors are compared with analytical results and a satisfactory agreement is found. The mesh configurations and the finite/infinite element (FIE) interface locations are systematically varied to establish a numerical modelling criterion for the most effective use of this approach. It is asserted that the truncation approach may now be abandoned in view of the proposed approach, its applicability and advantages. Keywords: infinite elements, finite element analysis, geotechnical engineering, numerical modelling, rock mechanics, stress concentration, tunnels, underground openings Underground excavations redistribute the initial stresses of the earthmass such that the modified stress pattern is characterized by a high concentration of compressive stresses in some locations around the edge of the opening, while undesirable tensile stresses occur in others. The purpose of the analysis presented is to locate, for all shapes, sizes and depths of opening, the magnitude of these stresses and the extent to which the effect of underground excavation is felt. The displacement at the edges of the opening may also be of interest. The inputs to the analysis are the initial stress pattern, geological features and mechanical characteristics of the underground medium. When the opening shape is simple and the geology of the earthmass regular, a closed-form analytical solution for the modified stresses and displacements can be derived by the theory of elasticity. L2 In most situations of practical interest, such ideal conditions are not found and a numerical method, such as FEM has to be employed. One disadvantage of this method of analysis in the context of underground openings is that it can 0307-904X/86/05357-10/$03.00 © 1986 Butterworth & Co. (Publishers) Ltd
not take into account the unbounded (or large) extent of the analysis domain, which is invariably present. Several approaches are available which enable the application of the finite element method to the above class of problems. The FIE method used in this paper effectively overcomes this problem. The method and the infinite elements used in this study are subsequently described and several problems are solved to evolve guidelines for the most effective use of this approach in solving the underground opening problem.
State of existing research Besides FEM, the boundary element method (BEM) has also been used in the analysis of underground openings. 3.4 A combination of the two has also been used successfully, s.6 The application of FEM with truncation of the analysis domain at a sufficiently large distance from the zone of disturbance continues to be the simplest and most widely used method. A numerical
AppI. Math. Modolling, 1986,Vol. 10, Octobor
357
Analysis of underground openings: P. Kumar modelling criterion for its most effective use is also available. 7.8. This criterion is evolved from the study of deep circular openings located in a homogeneous and isotropic medium, subject to a homogeneous stress field. According to this criterion, the truncation boundary should be located at least six radii away from the opening centre. The region of analysis so obtained must be divided into 125-150 FEs. The FE mesh should preferably be non-uniform so that the element size increases with the distance from the opening. The results of the analysis will then be of acceptable accuracy. The IEs were introduced by Bettess and Zienkiewicz9.m to analyse the hydrodynamic problems involving wave propagation. Their application was subsequently extended to the problems of geomechanics and elastostatics. ~ These IEs need a special numerical integration scheme extended over a semi-infinite range in the evaluation of their mass and stiffness properties. Recently, a new formulation of the IEs was proposed. ~2.~3 The most attractive feature of these IEs is that no special numerical integration scheme is needed in the evaluation of stiffness and mass properties. The drawback, on the other hand, is that their applicability is limited to the elastostatics problems only. The IEs derived from this formulation were onedimensional and an extension of these to a two-dimensional situation was subsequently achieved. ~4-~6 These IEs have since been used to solve several important elastostatic problems; more of these are reported herein. The other IEs, the formulation of which falls in this class are due to Bear and Meek. ~7The mapping function used in the derivation of these is relatively complicated. Nevertheless, these have been used in the stress analysis of underground excavation.~8 This study does not attempt to establish any numerical modelling criterion to enable their most efficient use.
Object and scope of this study The development of IEs has been relatively recent. In the near future, the use of these elements is likely to become as popular and routine as that of FEM itself. The following study is conducted with two objectives. Firstly, to demonstrate the applicability of the IEs to this class of problem and, secondly, to establish guidelines for their most efficient use. In this study, deep as well as shallow openings of circular and elliptical shape are analysed. For deep openings a purely vertical initial stress field is applied, while for shallow openings the stress filed is purely horizontal. Such stress fields are known to produce high stress concentrations. The medium containing these openings is taken as homogenous and isotropic. The modulus of elasticity and the Poisson ratio of this medium are 1000 MPa and 0.2, respectively. To establish a numerical modelling criterion, values of the stress concentration factors as obtained from the FIE analysis are compared with their exact values. The displacements are not considered because their analytical values for all cases of this study are not available. The merits of the FIE method can be established by comparison with the already available numerical modelling criterion for all-FE analysis. 7
358
Appl. Math. Modelling, 1986,Vo1.10, October
5
3
A W
y
5"
2
r
w
I
4
Figure 1 Infinite element geometry
Description of the infinite element and its implementation The coordinate transformation and the unknown function transformation for the IE are given by equations (1) and (2), respectively: /, X = NIx I + N2x2 +
• • • = "~
N~ri
( 1)
t=l
P
u = N,,, +
+.
= \7
(2)
i=1
Ni are shape functions given by: Nl = rs(1 --s)/(1 - r) N, = 2r(1 - s2)/(1 - r) N 3= rs(1 + s ) / ( 1 - r) N4= 0.5(1 + r)(1 - s ) / ( 1 - r) N.~ = 0.5(1 + r)(1 + s)/(1 - r) Ni are the shape functions for an eight-noded FE, x, and ug are the x-coordinate and u-displacements, respectively of node i, and p is the number of nodes in the element. Substituting the shape function Ni in equation (1) and solving for r gives: r = 1 - [ ( - x, - x 3 + x 4 + xs) -- s(x 4 - -
X I -~" X 3 - -
X5)]/
X4
Let x5 = Z~3 and x~ = 2x l, such that equation (3) reduces to:
r = 1 - [(x, +x3) + s ( - x l
+x3)]/x
(4)
Equation (4) shows that the IE shown in Figure 1 displays a reciprocal decay in the trial displacement function. Rearranging, this equation can be written as: x = [(x I + x3) + s ( - x| + x3)]/(1 - r)
(5)
The mapping by equation (5), described in Table 1, is valid when x 5 = 2x 3 and x4 = 2xl, which fixes the location of mid-side nodes. An examination of the mapping as described in Table I shows that equation (1) transforms the IE in the physical plane into a convenient regular shape in the natural plane. To evaluate its mass and stiffness properties, an integration over this finite regular shape can easily be performed by the Gauss-Legendre numerical integration scheme, which is usually available in FE analysis
Table 1
M a p p i n g b y e q u a t i o n (5)
f
S
Analysis of underground openings: P. Kumar Curvilinear mesh configuration
X
+1
All S
=
0
- 1 + 1
x = 2 x l ( = x 4) x = 2x3( = x 5)
-1
-1
x=x
0 +1
1 x l + x3
x=--
2
(=x2)
x=x 3
computer programs. Thus, an implementation of these elements into an existing FE analysis computer program presents no difficulty. A fundamental difference between the FE and IF subroutines lies in the different sets of shape functions which are to be used in the coordinate and function transformations of the IE. Out of these, one set of shape functions is common between them. Once this aspect is sorted out, an implementation of the proposed IE in an existing FEM computer program is straightforward, j~ The five-node IE shown in Figure 1 can be coupled with any two-dimensional quadratic FE. These two may require different orders of numerical integration as has been found in this study.
Method of analysis The analysis of an underground opening must take account of the fact that the entire opening is not excavated at once. Secondly, the excavated surface cannot have any stress on it except in the circumferential direction. A good account of the methods of analysis to simulate these effects is given by Kulhawy. 7 It has been proved that in the case of a homogeneous isotropic medium, stresses appearing on the excavated surface are independent of the cutting sequence.19 The medium considered in this study possesses these properties and, hence, a sequential construction procedure is not followed. This construction procedure corresponds to full-face tunneling. The balancing forces to make the excavated surface stress-free are computed from equations (6) and (7): F, = Shcos'-0 + S,.sin20
(6)
Ft = - (Sv - S,)sin0 cos0
(7)
where F, and Ft are the normal and tangential balancing forces, respectively, Sv and Sh are the intensities of the vertical and horizontal initial stress fields, respectively, and 0 is the inclination from horizontal of the outward normal at the excavated surface.
Study of deep openings In this part of the study, the analysis domain is divided into eight-node quadratic FEs and five-noded IEs as shown in Figure 1. A 3 × 3 Gaussian numerical integration is employed over FEs, the basis of which is available in standard text books. A 2 × 2 integration is used for the IEs. It is found that this order of integration is sufficiently accurate as well as economical. A higher order integration tends to make the IEs too stiff. The study variables are the number of segments into which the opening boundary is subdivided and the distance of the FIE interface from the origin.
To study the influence of the first variable, the FIE interface is fixed at 3.5 R where R is the radius of a circular opening or the largest dimension of a non-circular opening. This value is approximately half the distance of the truncation boundary according to Kulhawy's FE modelling criterion. 7 The number of subdivisions of the opening surface is then progressively increased from three to eight. The resulting mesh configurations are shown in Figures 2-4 for circular, tall elliptical and long elliptical openings, respectively. The number of radial lines in a region is chosen according to the magnitude of stress concentration to be expected there. This is necessary for an adequate representation of the stress gradient in that region.
Resldts and discussion. The analysis results and computer times for the solutions are given in Tables 2 and 3, respectively. It can be seen that while a stress concentration factor of up to 3 can be satisfactorily obtained by the simplest of FIE mesh A, a stress concentration factor in excess of 3 requires further effort. In the latter situation, a minimum of eight subdivisions of the excavated surface is found to be adequate. At least half of these should be located in the region of high stress concentration. When the stress concentration factor is below 3, a finer subdivision of the excavated surface offers no additional advantage, in addition, it imposes the penalty of higher analysis cost. The stress concentration factors at location I of the tall elliptical opening and location II of the long elliptical opening are of concern. Their values progressively improve as the number of radial lines increase. While the former continues to be underestimated, a trend of continuous improvement is upset when the number of IEs increase from six to seven. Interface location To study the influence of this variable, the number of subdivisions of the excavated surface is fixed at eight or nine and the FIE interface is moved towards the origin in discrete jumps by progressively stripping the outer layers of FEs from mesh configuration E. In mesh configurations F, G and J, the interface is placed at 2.5 R, 1.75R and 1.25 R, respectively, which corresponds to stripping one, two and four outer layers of FEs. All other aspects of the analysis are the same as those described in the preceding section.
Results and discussion. The results and computer times for the solutions are given in Tables 4 and 5, respectively. It can be seen that in the case of circular and tall elliptical openings, where the stress concentration factor extends to 3, a change in the interface location has no significant effect on the accuracy of results. For the case of long elliptical openings for which the stress concentration factor exceeds 3, the best results are obtained for the interface location at 1.75 R. The results of analysis with mesh J indicate that if the object of analysis is to estimate the stress concentration factors, one layer of FEs around the excavation surface followed by the IEs is adequate. The important question is to fix the optimum location of the FIE interface. The choice is not straightforward
Appl. Math. Modelling, 1986, Vol. 10, October
359
Analysis of underground openings: P. Kumar
..p
a
r~
I 4 ....
d
C
I
e Figure
2
Finite element mesh configurations for circular opening: (a) mesh A; (b) mesh B; (c) mesh D; (d) mesh E; (e) mesh K
because, in this study, the interface is placed at discrete selected locations. The best values of stress concentration factor at location I of the tall opening and location II of the long opening are obtained for the interface locations at 2.5 R and 1.75 R, respectively. The cost of analysis is smaller for the latter location. To seek a compromise between the accuracy and cost of analysis, the optimum interface location may lie within this range. There are other considerations which may force the choice of interface location. Theoretically, the FEs are employed to simulate the near-field while the IEs simulate the far-field in which the solution is governed by a decay characterizing the phenomenon. To utilize this criterion, the extent of the near-field has to be known,
360 Appl. Math. Modelling, 1986, Vol. 10, October
but is more often than not unavailable. Alternatively, it may become necessary to determine, as accurately as possible, the stress distribution in a certain region. Such is the case in a materially nonlinear analysis where the extent of plasticized zone needs to be determined. This can best be achieved by covering the entire anticipated plastic zone by FEs for the reasons described in a subsequent section. A similar situation arises when important geological features of the earthmass must be included in the analysis domain. Here also it can best be carried out by extending the region of FEs past this feature location. If none of these criteria are relevant, then, the results of the preceding paragraph may be used to fix the interface location.
Analysis of underground openings: P. Kumar
I
I
I I I - rI -
b
a
I I I I -{-__
I
c
d
I
I I
-Ie
Figure 3 Finite element mesh configurations for tall elliptical opening: (a) mesh A; (b) mesh B; (c) mesh D; (d) mesh E; (e) mesh K
Rectangular mesh configuration The mesh configurations of the preceding section are adequate when the earthmass containing the openings is either known to be isotropic, elastic and homogeneous, or can be assumed to be so. An analyst is not always so fortunate, in that the geology of the earth mass is far too complex to fit into the above simplified description. Its idealization as a stratified mass may, however, be appropriate for most practical purposes. In such situations the element boundaries are made to coincide with the plane of stratification. The resulting mesh configuration is rectangular. To cater for this situation, the three opening shapes of this study are also analysed in a rectangular mesh configuration while
retaining the earthmass to be isotropic, elastic and homogeneous. The excavation surface is split into nine subdivisions and the FIE interface is placed at 2.5 R. This mesh contains 165 nodes, 40 FEs and 13 IEs. All other details remain unchanged.
Results and discussion. An essential difference between the rectangular mesh K and curvilinear mesh E is that the former has an additional element boundary at 45 ° inclination and a large number of FEs are located around this line. The results for the rectangular mesh are given in the last column of Table 4. It can be seen that, for a maximum stress concentration factor of up to 3, the results of curvilinear
Appl. Math. M o d e l l i n g , 1986,Vol. 10, O c t o b e r
361
Analysis of underground openings: P. Kumar
I "-F . . . .
--L
L) .ulrll
a
. . . .
b
I
I
f
C
d
4
--
I
.
I
. . . . . I i
f Figure 4 Finite element mesh configurations for long elliptical opening: (a) mesh A; (b) mesh B; (c) mesh C; (d) mesh D; (e) mesh E; (f) mesh K
Table 2
Numerical and analytical values of stress concentration factors for finite/infinite interface fixed at 3.5 times opening radius Numerical value for curvilinear mesh A
B
C
D
E
85 nodes 20 FE 5 IE
100 nodes 24 FE 6 IE
138 nodes 35 FE 7 IE
156 nodes 40 FE 8 IE
Opening shape
Location
Theory
44 nodes 9 FE 3 IE
Circular
I Crown II Side wall
1 3
0.9202 2.9388
0.9993 3.0062
---
1.0082 3.0157
0.9932 3.0047
Tall elliptical
I Crown II Side wall
1 2
0.5128 2.0735
0.7338 2.0067
---
0.8263 2.0073
0.8158 2.0031
Long elliptical
I Crown II Side wall
1 5
1.0844 3.7441
1.0082 4.3337
1.0118 4.6806
1.0108 4.5954
1.0024 4.7594
362
Appl. Math. Modelling, 1986, Vol. 10, October
Analysis o f u n d e r g r o u n d openings: P. K u m a r Table 3 Computer time (s) for solutions in Table2 Shape of opening
Mesh A
Mesh B
Mesh C
Mesh D
Mesh E
Circular
10.59
21.06
--
33.68
39.67
9.88
20.80
--
34.05
39.22
10.12
20.01
23.70
32.37
36.46
Tall elliptical Long elliptical
Table 4 Numerical and analytical values of stress concentration factors for variable location of finite/infinite element interface Numerical values Curved mesh E
Curved mesh F
Curved mesh G
Curved mesh J
Rectangular mesh K
130 nodes 32 FIE 8 IE
104 nodes 24 FE 8 IE
58 nodes 9 FE 9 IE
165 nodes 40 Fe 15 IE
Opening shape
Location
Theory
156 nodes 40 FE 8 IE
Circular
I Crown II Side wall
1 3
0.9932 3.0047
1.0074 3.0173
1.0038 3.0136
1.0047 3.0044
1.0234 3.0325
Tall elliptical
I Crown II Side wall
1 2
0.8158 2.0031
0.8648 2.0486
0.8210 2.0079
0.7583 2.0133
0.8662 2.0519
Long elliptical
I Crown II Side wall
1 5
1.0024 4.7594
1.0106 4.7798
1.0129 4.7884
1.0393 4.5836
1.2216 4.9171
3.5 R
2.5 R
1.75 R
1.25 R
2.5 R
Location of finite/infinite element interface
Table 5 Computer time (s) for solutions in Table 4 Opening shape
Mesh E
Mesh F
Mesh G
Mesh J
Mesh K
Circular
39.67
32.97
24.05
13.36
52.80
Tall elliptical
39.22
32.70
23.35
13.33
54.25
Long elliptical
36.46
31.95
23.76
13.26
52.02
and rectangular mesh configurations are of the same order of accuracy. In the case of maximum stress concentration beyond 3, there is an improvement of about 2% at location II, offset by the stress concentration factor at location I which degrades by about 20%. In addition, this analysis can be seen to be the most expensive
(Table5).
Study of shallow openings The problem of shallow circular openings in a gravitational field was analysed by Mindlin. 2° Subsequently, this analysis was extended to shallow openings subject to a horizontal initial stress field. 21 Also, this theory was verified through photoelastic experiments. Two of these openings are analysed by the method of FIEs. The ratios of the depth of the opening from the free surface h to the opening radius R are 1.1 and 1.5. All other details of the medium containing these openings and those of the FE and IEs remain unchanged.
Shallow opening with h / R = 1.5 The FIE mesh used in this study contains 135 nodes, 32 FEs and 9 IEs and is shown in Figure 5. The FIE interface is placed at a distance of 3 R from the centre of the opening. The variation of stresses along the edge of the opening and along the free surface as obtained from the analysis are shown in Figures 6 and 7, respectively. The analytical curve in these figures corresponds
Figure 5 Finite element mesh for shallow circular opening (h/R = 1.5)
to h/R = 1.54, for which the analytical values are tabulated in reference 21. The experimental values are obtained from Figures 11 and 12 of this reference. Some salient stress concentration factors are given in Table 6.
Shallow opening with h/R =1.I In this part of the study, three FIE meshes are employed. The complete mesh geometry is shown in Figure 8(a). The details of the box portion is shown, on an enlarged scale, in Figures 8(b) and (c).
Appl. Math. Modelling,
1 9 8 6 , V o l . 10, O c t o b e r
363
Analysis of underground openings: P. Kumar I0
T
7 6 5 4q
3 2 I
0 -I
0
I
I
I
I
40
80
120
160
180
Angle of i n c l i n a t i o n , O ° Variation of stress along edge of opening: (--) theory; (O) photoelastic experiments21; (O) finite and infinite element analysis Figure 6
for shallow openings with h/R = 1.5, the stress concentration factors are in better agreement with the analytical than the experimental values. A similar mesh when used for a shallow opening with h/R = 1.1, is satisfactory except for the stress concentration factor at location m. Mesh A is slightly modified to obtain mesh B and the stress concentration factor at location n improves by about 1% while that at location m improves considerably. It is difficult to establish if this improvement results due to a reduction in the FE size in the region of high stress above the opening or due to a shift of FIE interface, or both. The performance of mesh C is highly satisfactory. In this study, all the stress concentration factors are satisfactorily obtained without having to change the FIE interface. There are more difficult opening shapes which may need more effort than this. These can easily be considered in this study, but then, the expected results from an alternative source are not available. In the absence of this comparison, this kind of exercise appears meaningless. More research in this direction is recommended.
Finite/infinite element approach against truncation approach
I... "-
3 2
I1 0
~-o
0
I
I
I
I
04
08
I-2
I-6
1"8
x/R Figure 7 Variation of stress along free surface: (--) theory; (O) photoelastic experiments21; (O) finite and infinite e l e m e n t
analysis
To obtain mesh configuration A, the element boundaries shown by dotted lines are chosen, ignoring the boundaries shown by the dot-dash line. In mesh configuration B, the dot-dash boundaries are used and dotted boundaries are ignored. In the latter configuration, the FIE interface location varies from 3 to 1.75 R. Other details of the mesh are the same as described in the preceding paragraph. In mesh configuration C, the detail of the box portion is as shown in Figure 8(c). The total mesh now contains 163 nodes, 40 FEs and 9 IEs. The variation of stresses along the opening edge and along the free surface are shown in Figures 9 and 10, respectively. The analytical curve shown in these figures correspond to h/R = 1.08 for which the tabulated values are available. The experimental values were read directly from Figure 11 of reference 21. Salient stress concentration factors are given in Table 6. Discussion A study of Figures 6, 7, 9 and 10 confirm the utility of the FIE approach in the analysis of shallow circular openings. It can be seen from a study of Table 6 that
364
Appl. Math. Modelling, 1986, Vol. 10, October
One major criticism of the truncation approach is that it introduces artificial boundaries, the location of which must be determined before an analysis can proceed. In the case of the FIE approach, no artificial boundaries are introduced but the location of the interface needs to be fixed. This study demonstrates that for the analysis of underground openings, the choice of interface location is not exceedingly critical and that the size of the problem to be solved is substantially smaller than in the other approach. Exceptions to this are those situations in which the plastic zone or some important geological feature of the earthmass is involved. The IEs do not permit a satisfactory modelling of variation in either material properties or geological features. In the case of IEs the Gauss points are usually far apart so that the resulting variation modelling is approximate. On the other hand, in a FE with a 3 × 3 integration, nine closely-spaced Gauss points are available, thereby, permitting a more favourable modelling of variation. So, if variation modelling is involved, it can be performed better by FEs, after which come either the FIE interface or the truncation boundary.
Computer times The solutions of this study were obtained on a D E C 2050 computer at Roorkee University Regional Computer Centre ( R U R C C ) running under the TOPS 20 Version 3A operating system. The software was developed by the present author and the reported computer times are all for computations in double precision.
Extensions and further applications The problems chosen for analysis in this paper have known solutions (either analytically or experimentally determined). Although the problems are somewhat artificial the analysis provides invaluable experience which can be successfully utilized in solving real world problems for which a solution is not known apriori. Without
Analysis of underground openings: P. Kumar Table 6 Stress concentration factors in shallow circular openings h/R
Mesh
Location
Theory*
m n o p
0.5 4.6 -3.3
0.53 5.08
m n
0.3 9.2
0.38 11.28
B
p m n
3.55 0.3 9.2
4.30 0.38 11.28
C
p m n
3.55 0.3 9.2
4.30 0.38 11.28
p
3.55
4.3
1.5
1.1
A
[
",
Experiment
3.54
Numerical value
Computer time (s)
0.4424 4.6287 1.1255 3.3164
37.71
1.1443 8.3974 0.8291 3.5668 0.9173 8.4661 0.6996 3.5816 0.3082 9.1269 0.7089 3.5838
33.92
35.12
41.76
*These values have been read from Figure lOof reference 21 and are approximate
ik II
i
,o
.,
/-
4 3
•
0
o
-I
~
•
-2-~ 0
I. . . . .
I
• I
• I
I
I
the background of this experience, the accuracy of the numerical solution cannot be established. Such a problem would be an (effectively) infinite plate containing rivet holes.
\ ~
40
•
80 120 160 180 Angle of inclination,8 ° Figure 9 Variation of stress along edge of opening: (--) theory; (O) photoelastic experiments21; (©) finite and infinite element analysis
a
I I
q
"\~,
Conclusions
'X,
"~ . . . . .
b
C Figure 8 Finite element mesh for shallow circular opening (h ~R =1.1): (a) complete mesh; (b) enlargement of boxed area (meshes A a n d B); (c) enlargement of boxed area (mesh C)
IEs with an inverse type of decay are suitable in the analysis of underground openings. To evaluate their stiffness properties, a 2 × 2 Gauss integration is adequate. It is known that accurate results can be obtained by judiciously selecting the FE mesh. This study attempts to quantify this statement with reference to the underground opening problem. In the analysis, the following guidelines, based on limited evidence, may be used. Many more opening shapes need to be analysed either to confirm or to improve upon these guidelines. • Certain advantages can be gained in the FIE analysis if the maximum value of the stress concentration factor can be estimated a priori. An example within the scope of this study is as follows. • When the maximum stress concentration factor is estimated to be below 3, three subdivisions of the
Appl. Math. Modelling, 1986, Vol. 10, October
365
Analysis of underground openings: P. Kumar 8
6
F.
0 0.0
I 0.4
I 0-8
I 1.2
I 1.6
1.8
x/,q Figure 10 Variation of stress along free surface: (--) theory; (Q) photoelastic experiments21; (©) finite and infinite element analysis
•
•
•
•
excavation surface are adequate. A further subdivision does not improve the results and raises the cost of analysis. When the highest stress concentration factor is estimated to be above 3 (or no estimate is available) up to nine subdivisions of the excavated surface are needed, at least half of which should be located in the region of expected high stress concentration. The location of the FIE interface depends on the objective of the analysis. If a rough estimate of stress concentration factors is required, the interface may be located close to the excavation surface with as little as one layer of FEs in between. If, however, an accurate analysis is required, the interface may be placed at a distance of 1.75 to 2.5 R from the origin, with four or more layers of FEs of gradually increasing size in between. There are, however, situations in which the interface may have to be placed beyond this suggested location. A change from a curvilinear mesh to a rectangular mesh configuration is not recommended unless forced by the local geological features. This numerical modelling criterion is based on stress factors. The displacements are not considered because results are not available for all cases for an alternative method.
The analysis of shallow underground openings further confirms the applicability of these IEs. The overall agreement of the numerical with analytical results is more satisfactory than that of the experimental results. Results of stresses and displacements for other opening shapes from alternative methods are not available and the present study of shallow openings suffers because of this omission. Further research to fill this gap is recommended.
Nomenclature E F.
modulus of elasticity normal balancing force on excavation surface
366 Appl. Math. Modelling, 1986,Vo1.10, October
Ft h k Ni ,~/ p r R s Sn S,. T x x~ u ui 0 o-¢ o-~
tangential balancing force on excavation surface depth of centre of opening from free surface ratio of horizontal to vertical initial stress field shape function at node i for infinite element shape function at node i for finite element number of node in element natural coordinate opening radius or largest dimension of opening natural coordinate intensity of horizontal initial stress field intensity of vertical initial stress field horizontal tensile stress field Cartesian coordinate or distance along free surface Cartesian coordinate of node i displacement displacement at node i angle of inclination circumferential stress in shallow opening stress at free surface
References 1 Savin, G. N. 'Stress concentration around holes', Pergamon Press, New York. 1961 2 Terzaghi, K. and Richart, F. E. 'Stresses in rock around cavities', Geotechn., 1952, 3 (2). 57-90 3 Crouch. S. L. and Starfield, A. M. 'Boundary element method in solid mechanics', George Allen and Unwin, London, 1983 4 Hoek, E. and Brown, E. T. 'Underground excavation in rock', Inst. Mining Metal., London, 1980 5 Brady, H. H. G. and Wassyng, A. "A coupled finite element boundary element method of stress analysis', Int. J. Rock Mech. Mining Sci., 1981, 18 (6), 475--486 6 Zienkiewicz, O. C., Kelly, D. W. and Bettess, P. "The coupling of finite element and boundary element procedures', Int. J. Numer. Methods Eng., 1977, I 1 (2), 355--375 7 Kulhawy, F. H. 'Finite element modelling criterion for underground openings in rock'. Int. J. Rock Mech. Mining Sci., 1974, 11 (12), 455-472 8 Christian, J. T. and Wong, T. H. "Errors in simulating excavation in elastic media by finite elements', Soils Found., 1973, 13 (1), 1-10 9 Bettess, P. "Infinite elements', Int. J. Numer. Methods Eng., 1977, I1 (1), 53--64 10 Beness, P. and Zienkiewicz, O. C. 'Diffraction and refraction of surface waves using finite and infinite elements', Int. J. Numer. Methods Eng., 1977, 11 (8), 1271-1290 11 Bettess, P. 'More on infinite elements', Int. J. Nurner. Methods Eng., 1980, 15 (11), 1613-1626 12 Curnier, A. 'Static infinite element', Int. J. Numer. Methods Eng., 1983, 19 (I0), 1479-1488 13 Zienkiewicz, O. C., Emson, C. and Bettess, P. "A novel boundary infinite element', Int. J. Numer. Methods Eng., 1983, 19 (3), 393--404 14 Kumar, P. "Novel infinite boundary element', Letter to Editor Int. J. Numer. Methods Eng., 1984, 20 (6), 1173-1174 15 Marques, J. M. M. C. and Owen, D. R. J. 'Infinite elements in quasi-static materially nonlinear problems', Computers and structures, 1984, 18 (4), 739-751 16 Kumar, P. 'On the static infinite element formulations', Am. Soc. Civ. Eng., J. Struct. Eng. 1986. 111 (STI 1), 2355-2372 17 Beer, G. and Meek, J. L. 'Infinite domain elments', Int. J. Numer. Methods Eng., 1981, 17 (I), 43-52 18 Beer, G. 'Infinite domain elements in finite element analysis of underground excavations', Int. J. Numer. Anal. Methods Geomech., 1983, 7 (1), 1-8 19 Ishihara, K. 'Relations between process ofcutting and uniqueness of solutions', Soils and Found., 1970, 10 (3), 50-65 20 Mindlin, R. D. 'Stress distribution around tunnel', Trans. Am. Soc. Cir. Eng., 1940, 105, 1117-1153 21 Mindlin, R. D. 'Stress distribution around hole near the edge of a plate under tension', Proc. Soc. Exp. Stress Anal. 1948,
s (2), 56-68
Underground openings subject to an initial stress field are usually analysed by the finite element method (FEM) by truncating the analysis domain at a large but finite distance from the zone of disturbance. In this paper, the problem is analysed by the method of finite (FE) and infinite elements (IE) to represent near- and far-fields, respectively. Deep openings of circular and elliptical shape and shallow circular openings are analysed to establish the approach. The numerical results of stress concentration factors are compared with analytical results and a satisfactory agreement is found. The mesh configurations and the finite/infinite element (FIE) interface locations are systematically varied to establish a numerical modelling criterion for the most effective use of this approach. It is asserted that the truncation approach may now be abandoned in view of the proposed approach, its applicability and advantages. Keywords: infinite elements, finite element analysis, geotechnical engineering, numerical modelling, rock mechanics, stress concentration, tunnels, underground openings Underground excavations redistribute the initial stresses of the earthmass such that the modified stress pattern is characterized by a high concentration of compressive stresses in some locations around the edge of the opening, while undesirable tensile stresses occur in others. The purpose of the analysis presented is to locate, for all shapes, sizes and depths of opening, the magnitude of these stresses and the extent to which the effect of underground excavation is felt. The displacement at the edges of the opening may also be of interest. The inputs to the analysis are the initial stress pattern, geological features and mechanical characteristics of the underground medium. When the opening shape is simple and the geology of the earthmass regular, a closed-form analytical solution for the modified stresses and displacements can be derived by the theory of elasticity. L2 In most situations of practical interest, such ideal conditions are not found and a numerical method, such as FEM has to be employed. One disadvantage of this method of analysis in the context of underground openings is that it can 0307-904X/86/05357-10/$03.00 © 1986 Butterworth & Co. (Publishers) Ltd
not take into account the unbounded (or large) extent of the analysis domain, which is invariably present. Several approaches are available which enable the application of the finite element method to the above class of problems. The FIE method used in this paper effectively overcomes this problem. The method and the infinite elements used in this study are subsequently described and several problems are solved to evolve guidelines for the most effective use of this approach in solving the underground opening problem.
State of existing research Besides FEM, the boundary element method (BEM) has also been used in the analysis of underground openings. 3.4 A combination of the two has also been used successfully, s.6 The application of FEM with truncation of the analysis domain at a sufficiently large distance from the zone of disturbance continues to be the simplest and most widely used method. A numerical
AppI. Math. Modolling, 1986,Vol. 10, Octobor
357
Analysis of underground openings: P. Kumar modelling criterion for its most effective use is also available. 7.8. This criterion is evolved from the study of deep circular openings located in a homogeneous and isotropic medium, subject to a homogeneous stress field. According to this criterion, the truncation boundary should be located at least six radii away from the opening centre. The region of analysis so obtained must be divided into 125-150 FEs. The FE mesh should preferably be non-uniform so that the element size increases with the distance from the opening. The results of the analysis will then be of acceptable accuracy. The IEs were introduced by Bettess and Zienkiewicz9.m to analyse the hydrodynamic problems involving wave propagation. Their application was subsequently extended to the problems of geomechanics and elastostatics. ~ These IEs need a special numerical integration scheme extended over a semi-infinite range in the evaluation of their mass and stiffness properties. Recently, a new formulation of the IEs was proposed. ~2.~3 The most attractive feature of these IEs is that no special numerical integration scheme is needed in the evaluation of stiffness and mass properties. The drawback, on the other hand, is that their applicability is limited to the elastostatics problems only. The IEs derived from this formulation were onedimensional and an extension of these to a two-dimensional situation was subsequently achieved. ~4-~6 These IEs have since been used to solve several important elastostatic problems; more of these are reported herein. The other IEs, the formulation of which falls in this class are due to Bear and Meek. ~7The mapping function used in the derivation of these is relatively complicated. Nevertheless, these have been used in the stress analysis of underground excavation.~8 This study does not attempt to establish any numerical modelling criterion to enable their most efficient use.
Object and scope of this study The development of IEs has been relatively recent. In the near future, the use of these elements is likely to become as popular and routine as that of FEM itself. The following study is conducted with two objectives. Firstly, to demonstrate the applicability of the IEs to this class of problem and, secondly, to establish guidelines for their most efficient use. In this study, deep as well as shallow openings of circular and elliptical shape are analysed. For deep openings a purely vertical initial stress field is applied, while for shallow openings the stress filed is purely horizontal. Such stress fields are known to produce high stress concentrations. The medium containing these openings is taken as homogenous and isotropic. The modulus of elasticity and the Poisson ratio of this medium are 1000 MPa and 0.2, respectively. To establish a numerical modelling criterion, values of the stress concentration factors as obtained from the FIE analysis are compared with their exact values. The displacements are not considered because their analytical values for all cases of this study are not available. The merits of the FIE method can be established by comparison with the already available numerical modelling criterion for all-FE analysis. 7
358
Appl. Math. Modelling, 1986,Vo1.10, October
5
3
A W
y
5"
2
r
w
I
4
Figure 1 Infinite element geometry
Description of the infinite element and its implementation The coordinate transformation and the unknown function transformation for the IE are given by equations (1) and (2), respectively: /, X = NIx I + N2x2 +
• • • = "~
N~ri
( 1)
t=l
P
u = N,,, +
+.
= \7
(2)
i=1
Ni are shape functions given by: Nl = rs(1 --s)/(1 - r) N, = 2r(1 - s2)/(1 - r) N 3= rs(1 + s ) / ( 1 - r) N4= 0.5(1 + r)(1 - s ) / ( 1 - r) N.~ = 0.5(1 + r)(1 + s)/(1 - r) Ni are the shape functions for an eight-noded FE, x, and ug are the x-coordinate and u-displacements, respectively of node i, and p is the number of nodes in the element. Substituting the shape function Ni in equation (1) and solving for r gives: r = 1 - [ ( - x, - x 3 + x 4 + xs) -- s(x 4 - -
X I -~" X 3 - -
X5)]/
X4
Let x5 = Z~3 and x~ = 2x l, such that equation (3) reduces to:
r = 1 - [(x, +x3) + s ( - x l
+x3)]/x
(4)
Equation (4) shows that the IE shown in Figure 1 displays a reciprocal decay in the trial displacement function. Rearranging, this equation can be written as: x = [(x I + x3) + s ( - x| + x3)]/(1 - r)
(5)
The mapping by equation (5), described in Table 1, is valid when x 5 = 2x 3 and x4 = 2xl, which fixes the location of mid-side nodes. An examination of the mapping as described in Table I shows that equation (1) transforms the IE in the physical plane into a convenient regular shape in the natural plane. To evaluate its mass and stiffness properties, an integration over this finite regular shape can easily be performed by the Gauss-Legendre numerical integration scheme, which is usually available in FE analysis
Table 1
M a p p i n g b y e q u a t i o n (5)
f
S
Analysis of underground openings: P. Kumar Curvilinear mesh configuration
X
+1
All S
=
0
- 1 + 1
x = 2 x l ( = x 4) x = 2x3( = x 5)
-1
-1
x=x
0 +1
1 x l + x3
x=--
2
(=x2)
x=x 3
computer programs. Thus, an implementation of these elements into an existing FE analysis computer program presents no difficulty. A fundamental difference between the FE and IF subroutines lies in the different sets of shape functions which are to be used in the coordinate and function transformations of the IE. Out of these, one set of shape functions is common between them. Once this aspect is sorted out, an implementation of the proposed IE in an existing FEM computer program is straightforward, j~ The five-node IE shown in Figure 1 can be coupled with any two-dimensional quadratic FE. These two may require different orders of numerical integration as has been found in this study.
Method of analysis The analysis of an underground opening must take account of the fact that the entire opening is not excavated at once. Secondly, the excavated surface cannot have any stress on it except in the circumferential direction. A good account of the methods of analysis to simulate these effects is given by Kulhawy. 7 It has been proved that in the case of a homogeneous isotropic medium, stresses appearing on the excavated surface are independent of the cutting sequence.19 The medium considered in this study possesses these properties and, hence, a sequential construction procedure is not followed. This construction procedure corresponds to full-face tunneling. The balancing forces to make the excavated surface stress-free are computed from equations (6) and (7): F, = Shcos'-0 + S,.sin20
(6)
Ft = - (Sv - S,)sin0 cos0
(7)
where F, and Ft are the normal and tangential balancing forces, respectively, Sv and Sh are the intensities of the vertical and horizontal initial stress fields, respectively, and 0 is the inclination from horizontal of the outward normal at the excavated surface.
Study of deep openings In this part of the study, the analysis domain is divided into eight-node quadratic FEs and five-noded IEs as shown in Figure 1. A 3 × 3 Gaussian numerical integration is employed over FEs, the basis of which is available in standard text books. A 2 × 2 integration is used for the IEs. It is found that this order of integration is sufficiently accurate as well as economical. A higher order integration tends to make the IEs too stiff. The study variables are the number of segments into which the opening boundary is subdivided and the distance of the FIE interface from the origin.
To study the influence of the first variable, the FIE interface is fixed at 3.5 R where R is the radius of a circular opening or the largest dimension of a non-circular opening. This value is approximately half the distance of the truncation boundary according to Kulhawy's FE modelling criterion. 7 The number of subdivisions of the opening surface is then progressively increased from three to eight. The resulting mesh configurations are shown in Figures 2-4 for circular, tall elliptical and long elliptical openings, respectively. The number of radial lines in a region is chosen according to the magnitude of stress concentration to be expected there. This is necessary for an adequate representation of the stress gradient in that region.
Resldts and discussion. The analysis results and computer times for the solutions are given in Tables 2 and 3, respectively. It can be seen that while a stress concentration factor of up to 3 can be satisfactorily obtained by the simplest of FIE mesh A, a stress concentration factor in excess of 3 requires further effort. In the latter situation, a minimum of eight subdivisions of the excavated surface is found to be adequate. At least half of these should be located in the region of high stress concentration. When the stress concentration factor is below 3, a finer subdivision of the excavated surface offers no additional advantage, in addition, it imposes the penalty of higher analysis cost. The stress concentration factors at location I of the tall elliptical opening and location II of the long elliptical opening are of concern. Their values progressively improve as the number of radial lines increase. While the former continues to be underestimated, a trend of continuous improvement is upset when the number of IEs increase from six to seven. Interface location To study the influence of this variable, the number of subdivisions of the excavated surface is fixed at eight or nine and the FIE interface is moved towards the origin in discrete jumps by progressively stripping the outer layers of FEs from mesh configuration E. In mesh configurations F, G and J, the interface is placed at 2.5 R, 1.75R and 1.25 R, respectively, which corresponds to stripping one, two and four outer layers of FEs. All other aspects of the analysis are the same as those described in the preceding section.
Results and discussion. The results and computer times for the solutions are given in Tables 4 and 5, respectively. It can be seen that in the case of circular and tall elliptical openings, where the stress concentration factor extends to 3, a change in the interface location has no significant effect on the accuracy of results. For the case of long elliptical openings for which the stress concentration factor exceeds 3, the best results are obtained for the interface location at 1.75 R. The results of analysis with mesh J indicate that if the object of analysis is to estimate the stress concentration factors, one layer of FEs around the excavation surface followed by the IEs is adequate. The important question is to fix the optimum location of the FIE interface. The choice is not straightforward
Appl. Math. Modelling, 1986, Vol. 10, October
359
Analysis of underground openings: P. Kumar
..p
a
r~
I 4 ....
d
C
I
e Figure
2
Finite element mesh configurations for circular opening: (a) mesh A; (b) mesh B; (c) mesh D; (d) mesh E; (e) mesh K
because, in this study, the interface is placed at discrete selected locations. The best values of stress concentration factor at location I of the tall opening and location II of the long opening are obtained for the interface locations at 2.5 R and 1.75 R, respectively. The cost of analysis is smaller for the latter location. To seek a compromise between the accuracy and cost of analysis, the optimum interface location may lie within this range. There are other considerations which may force the choice of interface location. Theoretically, the FEs are employed to simulate the near-field while the IEs simulate the far-field in which the solution is governed by a decay characterizing the phenomenon. To utilize this criterion, the extent of the near-field has to be known,
360 Appl. Math. Modelling, 1986, Vol. 10, October
but is more often than not unavailable. Alternatively, it may become necessary to determine, as accurately as possible, the stress distribution in a certain region. Such is the case in a materially nonlinear analysis where the extent of plasticized zone needs to be determined. This can best be achieved by covering the entire anticipated plastic zone by FEs for the reasons described in a subsequent section. A similar situation arises when important geological features of the earthmass must be included in the analysis domain. Here also it can best be carried out by extending the region of FEs past this feature location. If none of these criteria are relevant, then, the results of the preceding paragraph may be used to fix the interface location.
Analysis of underground openings: P. Kumar
I
I
I I I - rI -
b
a
I I I I -{-__
I
c
d
I
I I
-Ie
Figure 3 Finite element mesh configurations for tall elliptical opening: (a) mesh A; (b) mesh B; (c) mesh D; (d) mesh E; (e) mesh K
Rectangular mesh configuration The mesh configurations of the preceding section are adequate when the earthmass containing the openings is either known to be isotropic, elastic and homogeneous, or can be assumed to be so. An analyst is not always so fortunate, in that the geology of the earth mass is far too complex to fit into the above simplified description. Its idealization as a stratified mass may, however, be appropriate for most practical purposes. In such situations the element boundaries are made to coincide with the plane of stratification. The resulting mesh configuration is rectangular. To cater for this situation, the three opening shapes of this study are also analysed in a rectangular mesh configuration while
retaining the earthmass to be isotropic, elastic and homogeneous. The excavation surface is split into nine subdivisions and the FIE interface is placed at 2.5 R. This mesh contains 165 nodes, 40 FEs and 13 IEs. All other details remain unchanged.
Results and discussion. An essential difference between the rectangular mesh K and curvilinear mesh E is that the former has an additional element boundary at 45 ° inclination and a large number of FEs are located around this line. The results for the rectangular mesh are given in the last column of Table 4. It can be seen that, for a maximum stress concentration factor of up to 3, the results of curvilinear
Appl. Math. M o d e l l i n g , 1986,Vol. 10, O c t o b e r
361
Analysis of underground openings: P. Kumar
I "-F . . . .
--L
L) .ulrll
a
. . . .
b
I
I
f
C
d
4
--
I
.
I
. . . . . I i
f Figure 4 Finite element mesh configurations for long elliptical opening: (a) mesh A; (b) mesh B; (c) mesh C; (d) mesh D; (e) mesh E; (f) mesh K
Table 2
Numerical and analytical values of stress concentration factors for finite/infinite interface fixed at 3.5 times opening radius Numerical value for curvilinear mesh A
B
C
D
E
85 nodes 20 FE 5 IE
100 nodes 24 FE 6 IE
138 nodes 35 FE 7 IE
156 nodes 40 FE 8 IE
Opening shape
Location
Theory
44 nodes 9 FE 3 IE
Circular
I Crown II Side wall
1 3
0.9202 2.9388
0.9993 3.0062
---
1.0082 3.0157
0.9932 3.0047
Tall elliptical
I Crown II Side wall
1 2
0.5128 2.0735
0.7338 2.0067
---
0.8263 2.0073
0.8158 2.0031
Long elliptical
I Crown II Side wall
1 5
1.0844 3.7441
1.0082 4.3337
1.0118 4.6806
1.0108 4.5954
1.0024 4.7594
362
Appl. Math. Modelling, 1986, Vol. 10, October
Analysis o f u n d e r g r o u n d openings: P. K u m a r Table 3 Computer time (s) for solutions in Table2 Shape of opening
Mesh A
Mesh B
Mesh C
Mesh D
Mesh E
Circular
10.59
21.06
--
33.68
39.67
9.88
20.80
--
34.05
39.22
10.12
20.01
23.70
32.37
36.46
Tall elliptical Long elliptical
Table 4 Numerical and analytical values of stress concentration factors for variable location of finite/infinite element interface Numerical values Curved mesh E
Curved mesh F
Curved mesh G
Curved mesh J
Rectangular mesh K
130 nodes 32 FIE 8 IE
104 nodes 24 FE 8 IE
58 nodes 9 FE 9 IE
165 nodes 40 Fe 15 IE
Opening shape
Location
Theory
156 nodes 40 FE 8 IE
Circular
I Crown II Side wall
1 3
0.9932 3.0047
1.0074 3.0173
1.0038 3.0136
1.0047 3.0044
1.0234 3.0325
Tall elliptical
I Crown II Side wall
1 2
0.8158 2.0031
0.8648 2.0486
0.8210 2.0079
0.7583 2.0133
0.8662 2.0519
Long elliptical
I Crown II Side wall
1 5
1.0024 4.7594
1.0106 4.7798
1.0129 4.7884
1.0393 4.5836
1.2216 4.9171
3.5 R
2.5 R
1.75 R
1.25 R
2.5 R
Location of finite/infinite element interface
Table 5 Computer time (s) for solutions in Table 4 Opening shape
Mesh E
Mesh F
Mesh G
Mesh J
Mesh K
Circular
39.67
32.97
24.05
13.36
52.80
Tall elliptical
39.22
32.70
23.35
13.33
54.25
Long elliptical
36.46
31.95
23.76
13.26
52.02
and rectangular mesh configurations are of the same order of accuracy. In the case of maximum stress concentration beyond 3, there is an improvement of about 2% at location II, offset by the stress concentration factor at location I which degrades by about 20%. In addition, this analysis can be seen to be the most expensive
(Table5).
Study of shallow openings The problem of shallow circular openings in a gravitational field was analysed by Mindlin. 2° Subsequently, this analysis was extended to shallow openings subject to a horizontal initial stress field. 21 Also, this theory was verified through photoelastic experiments. Two of these openings are analysed by the method of FIEs. The ratios of the depth of the opening from the free surface h to the opening radius R are 1.1 and 1.5. All other details of the medium containing these openings and those of the FE and IEs remain unchanged.
Shallow opening with h / R = 1.5 The FIE mesh used in this study contains 135 nodes, 32 FEs and 9 IEs and is shown in Figure 5. The FIE interface is placed at a distance of 3 R from the centre of the opening. The variation of stresses along the edge of the opening and along the free surface as obtained from the analysis are shown in Figures 6 and 7, respectively. The analytical curve in these figures corresponds
Figure 5 Finite element mesh for shallow circular opening (h/R = 1.5)
to h/R = 1.54, for which the analytical values are tabulated in reference 21. The experimental values are obtained from Figures 11 and 12 of this reference. Some salient stress concentration factors are given in Table 6.
Shallow opening with h/R =1.I In this part of the study, three FIE meshes are employed. The complete mesh geometry is shown in Figure 8(a). The details of the box portion is shown, on an enlarged scale, in Figures 8(b) and (c).
Appl. Math. Modelling,
1 9 8 6 , V o l . 10, O c t o b e r
363
Analysis of underground openings: P. Kumar I0
T
7 6 5 4q
3 2 I
0 -I
0
I
I
I
I
40
80
120
160
180
Angle of i n c l i n a t i o n , O ° Variation of stress along edge of opening: (--) theory; (O) photoelastic experiments21; (O) finite and infinite element analysis Figure 6
for shallow openings with h/R = 1.5, the stress concentration factors are in better agreement with the analytical than the experimental values. A similar mesh when used for a shallow opening with h/R = 1.1, is satisfactory except for the stress concentration factor at location m. Mesh A is slightly modified to obtain mesh B and the stress concentration factor at location n improves by about 1% while that at location m improves considerably. It is difficult to establish if this improvement results due to a reduction in the FE size in the region of high stress above the opening or due to a shift of FIE interface, or both. The performance of mesh C is highly satisfactory. In this study, all the stress concentration factors are satisfactorily obtained without having to change the FIE interface. There are more difficult opening shapes which may need more effort than this. These can easily be considered in this study, but then, the expected results from an alternative source are not available. In the absence of this comparison, this kind of exercise appears meaningless. More research in this direction is recommended.
Finite/infinite element approach against truncation approach
I... "-
3 2
I1 0
~-o
0
I
I
I
I
04
08
I-2
I-6
1"8
x/R Figure 7 Variation of stress along free surface: (--) theory; (O) photoelastic experiments21; (O) finite and infinite e l e m e n t
analysis
To obtain mesh configuration A, the element boundaries shown by dotted lines are chosen, ignoring the boundaries shown by the dot-dash line. In mesh configuration B, the dot-dash boundaries are used and dotted boundaries are ignored. In the latter configuration, the FIE interface location varies from 3 to 1.75 R. Other details of the mesh are the same as described in the preceding paragraph. In mesh configuration C, the detail of the box portion is as shown in Figure 8(c). The total mesh now contains 163 nodes, 40 FEs and 9 IEs. The variation of stresses along the opening edge and along the free surface are shown in Figures 9 and 10, respectively. The analytical curve shown in these figures correspond to h/R = 1.08 for which the tabulated values are available. The experimental values were read directly from Figure 11 of reference 21. Salient stress concentration factors are given in Table 6. Discussion A study of Figures 6, 7, 9 and 10 confirm the utility of the FIE approach in the analysis of shallow circular openings. It can be seen from a study of Table 6 that
364
Appl. Math. Modelling, 1986, Vol. 10, October
One major criticism of the truncation approach is that it introduces artificial boundaries, the location of which must be determined before an analysis can proceed. In the case of the FIE approach, no artificial boundaries are introduced but the location of the interface needs to be fixed. This study demonstrates that for the analysis of underground openings, the choice of interface location is not exceedingly critical and that the size of the problem to be solved is substantially smaller than in the other approach. Exceptions to this are those situations in which the plastic zone or some important geological feature of the earthmass is involved. The IEs do not permit a satisfactory modelling of variation in either material properties or geological features. In the case of IEs the Gauss points are usually far apart so that the resulting variation modelling is approximate. On the other hand, in a FE with a 3 × 3 integration, nine closely-spaced Gauss points are available, thereby, permitting a more favourable modelling of variation. So, if variation modelling is involved, it can be performed better by FEs, after which come either the FIE interface or the truncation boundary.
Computer times The solutions of this study were obtained on a D E C 2050 computer at Roorkee University Regional Computer Centre ( R U R C C ) running under the TOPS 20 Version 3A operating system. The software was developed by the present author and the reported computer times are all for computations in double precision.
Extensions and further applications The problems chosen for analysis in this paper have known solutions (either analytically or experimentally determined). Although the problems are somewhat artificial the analysis provides invaluable experience which can be successfully utilized in solving real world problems for which a solution is not known apriori. Without
Analysis of underground openings: P. Kumar Table 6 Stress concentration factors in shallow circular openings h/R
Mesh
Location
Theory*
m n o p
0.5 4.6 -3.3
0.53 5.08
m n
0.3 9.2
0.38 11.28
B
p m n
3.55 0.3 9.2
4.30 0.38 11.28
C
p m n
3.55 0.3 9.2
4.30 0.38 11.28
p
3.55
4.3
1.5
1.1
A
[
",
Experiment
3.54
Numerical value
Computer time (s)
0.4424 4.6287 1.1255 3.3164
37.71
1.1443 8.3974 0.8291 3.5668 0.9173 8.4661 0.6996 3.5816 0.3082 9.1269 0.7089 3.5838
33.92
35.12
41.76
*These values have been read from Figure lOof reference 21 and are approximate
ik II
i
,o
.,
/-
4 3
•
0
o
-I
~
•
-2-~ 0
I. . . . .
I
• I
• I
I
I
the background of this experience, the accuracy of the numerical solution cannot be established. Such a problem would be an (effectively) infinite plate containing rivet holes.
\ ~
40
•
80 120 160 180 Angle of inclination,8 ° Figure 9 Variation of stress along edge of opening: (--) theory; (O) photoelastic experiments21; (©) finite and infinite element analysis
a
I I
q
"\~,
Conclusions
'X,
"~ . . . . .
b
C Figure 8 Finite element mesh for shallow circular opening (h ~R =1.1): (a) complete mesh; (b) enlargement of boxed area (meshes A a n d B); (c) enlargement of boxed area (mesh C)
IEs with an inverse type of decay are suitable in the analysis of underground openings. To evaluate their stiffness properties, a 2 × 2 Gauss integration is adequate. It is known that accurate results can be obtained by judiciously selecting the FE mesh. This study attempts to quantify this statement with reference to the underground opening problem. In the analysis, the following guidelines, based on limited evidence, may be used. Many more opening shapes need to be analysed either to confirm or to improve upon these guidelines. • Certain advantages can be gained in the FIE analysis if the maximum value of the stress concentration factor can be estimated a priori. An example within the scope of this study is as follows. • When the maximum stress concentration factor is estimated to be below 3, three subdivisions of the
Appl. Math. Modelling, 1986, Vol. 10, October
365
Analysis of underground openings: P. Kumar 8
6
F.
0 0.0
I 0.4
I 0-8
I 1.2
I 1.6
1.8
x/,q Figure 10 Variation of stress along free surface: (--) theory; (Q) photoelastic experiments21; (©) finite and infinite element analysis
•
•
•
•
excavation surface are adequate. A further subdivision does not improve the results and raises the cost of analysis. When the highest stress concentration factor is estimated to be above 3 (or no estimate is available) up to nine subdivisions of the excavated surface are needed, at least half of which should be located in the region of expected high stress concentration. The location of the FIE interface depends on the objective of the analysis. If a rough estimate of stress concentration factors is required, the interface may be located close to the excavation surface with as little as one layer of FEs in between. If, however, an accurate analysis is required, the interface may be placed at a distance of 1.75 to 2.5 R from the origin, with four or more layers of FEs of gradually increasing size in between. There are, however, situations in which the interface may have to be placed beyond this suggested location. A change from a curvilinear mesh to a rectangular mesh configuration is not recommended unless forced by the local geological features. This numerical modelling criterion is based on stress factors. The displacements are not considered because results are not available for all cases for an alternative method.
The analysis of shallow underground openings further confirms the applicability of these IEs. The overall agreement of the numerical with analytical results is more satisfactory than that of the experimental results. Results of stresses and displacements for other opening shapes from alternative methods are not available and the present study of shallow openings suffers because of this omission. Further research to fill this gap is recommended.
Nomenclature E F.
modulus of elasticity normal balancing force on excavation surface
366 Appl. Math. Modelling, 1986,Vo1.10, October
Ft h k Ni ,~/ p r R s Sn S,. T x x~ u ui 0 o-¢ o-~
tangential balancing force on excavation surface depth of centre of opening from free surface ratio of horizontal to vertical initial stress field shape function at node i for infinite element shape function at node i for finite element number of node in element natural coordinate opening radius or largest dimension of opening natural coordinate intensity of horizontal initial stress field intensity of vertical initial stress field horizontal tensile stress field Cartesian coordinate or distance along free surface Cartesian coordinate of node i displacement displacement at node i angle of inclination circumferential stress in shallow opening stress at free surface
References 1 Savin, G. N. 'Stress concentration around holes', Pergamon Press, New York. 1961 2 Terzaghi, K. and Richart, F. E. 'Stresses in rock around cavities', Geotechn., 1952, 3 (2). 57-90 3 Crouch. S. L. and Starfield, A. M. 'Boundary element method in solid mechanics', George Allen and Unwin, London, 1983 4 Hoek, E. and Brown, E. T. 'Underground excavation in rock', Inst. Mining Metal., London, 1980 5 Brady, H. H. G. and Wassyng, A. "A coupled finite element boundary element method of stress analysis', Int. J. Rock Mech. Mining Sci., 1981, 18 (6), 475--486 6 Zienkiewicz, O. C., Kelly, D. W. and Bettess, P. "The coupling of finite element and boundary element procedures', Int. J. Numer. Methods Eng., 1977, I 1 (2), 355--375 7 Kulhawy, F. H. 'Finite element modelling criterion for underground openings in rock'. Int. J. Rock Mech. Mining Sci., 1974, 11 (12), 455-472 8 Christian, J. T. and Wong, T. H. "Errors in simulating excavation in elastic media by finite elements', Soils Found., 1973, 13 (1), 1-10 9 Bettess, P. "Infinite elements', Int. J. Numer. Methods Eng., 1977, I1 (1), 53--64 10 Beness, P. and Zienkiewicz, O. C. 'Diffraction and refraction of surface waves using finite and infinite elements', Int. J. Numer. Methods Eng., 1977, 11 (8), 1271-1290 11 Bettess, P. 'More on infinite elements', Int. J. Nurner. Methods Eng., 1980, 15 (11), 1613-1626 12 Curnier, A. 'Static infinite element', Int. J. Numer. Methods Eng., 1983, 19 (I0), 1479-1488 13 Zienkiewicz, O. C., Emson, C. and Bettess, P. "A novel boundary infinite element', Int. J. Numer. Methods Eng., 1983, 19 (3), 393--404 14 Kumar, P. "Novel infinite boundary element', Letter to Editor Int. J. Numer. Methods Eng., 1984, 20 (6), 1173-1174 15 Marques, J. M. M. C. and Owen, D. R. J. 'Infinite elements in quasi-static materially nonlinear problems', Computers and structures, 1984, 18 (4), 739-751 16 Kumar, P. 'On the static infinite element formulations', Am. Soc. Civ. Eng., J. Struct. Eng. 1986. 111 (STI 1), 2355-2372 17 Beer, G. and Meek, J. L. 'Infinite domain elments', Int. J. Numer. Methods Eng., 1981, 17 (I), 43-52 18 Beer, G. 'Infinite domain elements in finite element analysis of underground excavations', Int. J. Numer. Anal. Methods Geomech., 1983, 7 (1), 1-8 19 Ishihara, K. 'Relations between process ofcutting and uniqueness of solutions', Soils and Found., 1970, 10 (3), 50-65 20 Mindlin, R. D. 'Stress distribution around tunnel', Trans. Am. Soc. Cir. Eng., 1940, 105, 1117-1153 21 Mindlin, R. D. 'Stress distribution around hole near the edge of a plate under tension', Proc. Soc. Exp. Stress Anal. 1948,
s (2), 56-68