- Email: [email protected]
Efficient mesh generation procedures for finite element analysis of underground structures
Int. J. Rock Mech. Min. Sci. & Geomech. Abstr. Vol. 30, No. 6, pp. 591-600, 1993
0148-9062/93 $6.00 + 0.00 Copyright © 1993 Pergamon Press Lid
Printed in Great Britain. All rights reserved
Efficient Mesh Generation Procedures for Finite Element Analysis of Underground Structures M. W. F. GRABINSKYt J. H. CURRANt
The efficiency and accuracy of the finite element method & dependent on the volume mesh developed for analysis. This can be particularly significant for 3-D analyses of unbounded problems such as underground structures. Several developments in automatic mesh generation have been made in recent years and are increasingly found in commercially available software, although these programs target mainly mechanical engineering applications. The work described in this paper considers application of advanced mesh generation techniques to the analysis of underground structures. The functionality of the software used is described, and some general considerations for mesh development are given. Application of the software is then demonstrated using two mining and two civil engineering structures. It is shown that the newer approaches to mesh generation can provide significant advantages in terms of computational efficiency and modelling functionality. However, further software development to make these modelling techniques more accessible to the geotechnical designer is desirable.
INTRODUCTION Many commercially available finite element programs produce meshes using a technique known as "mapped meshing" in which a topologically regular set of elements is distorted to fit the problem geometry (Fig. 1). While this approach is advantageous from the programmer's point of view, it becomes awkward to apply where the problem geometry is irregular and/or complex. Further, meshes produced using mapping techniques are notoriously inefficient for modelling 3-D problems with "infinite" domains. Issues of computational efficiency have been addressed by developing special elements [1-3] and hybrid computational schemes where a finite element mesh (or similar domain mesh) is mated to a boundary element mesh [4-8]. (It should be noted that hybrid schemes, while theoretically attractive, do not necessarily provide clear computational advantages for all kinds of problems since the specialized solvers developed for finite element matrices are incompatible with matrices produced by boundary elements [28].) Nevertheless, the practitioner is still left with the problem of trying to generate a mesh within and about complex geometries, even if specialized elements or hybrid computational schemes are used. An alternative approach to enhancing the efficiency of finite element analyses lies in a different method of
generating the mesh. "Free meshing" is a term used to describe approaches in which nodes are inserted into the problem domain and then elements are connected to these nodes (Fig. 2). In the node insertion part of the process, it is not necessary for the nodes to be placed in any preordained topology. Various free meshing algorithms have been developed and are increasingly found in many software programs intended primarily for mechanical engineering applications [9-12]. The main advantages of the free meshing technique are that a mesh can
tDepartment of Civil Engineering,Universityof Toronto, Toronto, Ontario, Canada MSS IA4. 591
Fig. 1. Exampleof mapped mesh.
592
GRABINSKY and CURRAN:
MESH GENERATION FOR FEA MESH GENERATION USING FREE MESHING
Fig. 2. Example of free mesh.
be produced in volumes with very general geometries, and greater control is gained over the densities of the nodes in the model (and therefore the model's computational efficiency). The enhanced modelling functionality comes at the price of increased programming effort, though this is not of direct concern to the end user. A commercially available software package called I-DEAS [13] has been used to investigate how recently developed geometric modelling and mesh generation functionality can best be applied to underground geomechanical structures. The methods used to generate finite element meshes using I-DEAS and application of these methods to four case histories are described in this paper. The consequences of using different approaches (free meshing vs mapped meshing with special enhancements) to finite element analysis of underground structures is then considered. The finite element meshes generated using I-DEAS were analyzed using the ABAQUS program [14]. In addition to the geomechanical constitutive models available in ABAQUS (extended Drucker-Prager and Critical State models), custom constitutive models (Hoek-Brown [15], Elastic-Brittle-Damage [16] and ubiquitous joint [17] models) were developed and used. All of the work was carried out on a Silicon Graphics Inc. (SGI) Iris Indigo R3000. This machine has processing capabilities similar to an Intel 486/66 MHz-based personal computer, but is equipped with special video capabilities for real-time animation of 3-D models. The video capabilities greatly enhance the development of model geometries and meshes, and significantly reduce the amount of time required to display 3-D analysis results. tin the context of this discussion, geometric modelling includes all activities used to definethe geometryof the model. It is more usual to considergeometricmodellingas the mathematicaldescriptionof curves and surfaces, and solid modelling as the representation of physical solid objects (i.e. the ordering of curves and surfaces into bounded, sensible solids).
Recent approaches to mesh generation tend to consider goemetric modellingt and element discretization as separate activities. In this way, the discretization is developed on the basis of the underlying goemetric definition, and therefore many different discretizations of the same geometry may be developed independently of each other. This is different from the mesh generation techniques commonly employed in programs based on mapped meshing, where the geometry must be defined in terms of modelling primitives (cubes, wedges and tetrahedrons) which are subsequently subdivided to produce the elements. Hence, geometric modelling and mesh generation are considered separately below. A further issue requiring consideration is the model representation of far-field boundary conditions: where the mesh is truncated and what combination of loads and displacements are prescribed at these outer limits. It will be shown that this issue can be handled easily using free meshing techniques. Geometric modelling Geometric information may be associated with a variety of parameters relevant to rock mechanics modelling, for example: the surfaces of excavations or regions yet to be excavated; the contacts of different geomechanical domains; the locations of faults, shear zones and other significant discontinuities; the attitudes of joint sets, bedding planes and schistocity; the locations of instrumentation cells and anchor points; and the locations of microseismic events. Various CAD and mine planning software tools are widely available and can be used to great advantage in generating, displaying and interpreting such information. However, the geometric modelling functionality of these tools can vary significantly in sophistication, and care must be taken when attempting to use these tools to produce geometry as a basis for generating finite element meshes. Rudimentary geometric modelling tools typically allow the user to define a series of profiles, each profile being a polyline or a string of low-order curves, and to (semi)-automatically connect these profiles using triangular surfaces. The resulting surface which envelopes the profiles is then faceted, not smooth. Functionality may be added to enable the modelling of holes in an object, or various geometric bifurcations (branches, Ys, etc.). The program may be written to calculate volumetric information in an ad hoe manner as the user generates the geometry. Although the foregoing approach to geometric modelling is comparatively unsophisticated, it can adequately fulfill the functional requirements of many geomechanical design projects. For example, the surface of a geologic contact may be interpolated on the basis of a few point locations determined from a drilling program, and fitting anything more sophisticated than a piecewise linear surface through these points is generally unjustified. Similarly, higher-order curves are often used to define the desired excavation profiles of many civil
GRABINSKY and CURRAN: MESH GENERATION FOR FEA engineering structures, but the as-built profiles may vary from the design profile to the same degree as a piecewise linear approximation varies from the higher-order curve. The currently favored approach to geometric modelling involves the description of curved lines and surfaces using non-uniform rational B-splines (NURBS) and the building of objects using constructive solid modelling and/or boundary representation modelling [18, 19]. Various functionality is used to define NURBS-based objects from geometric primitives, extrusions or sweeps of 2-D profiles, skinning of a series of profiles, and so on. Boolean operations may then be used to combine the basic objects, thereby producing more complex objects. Throughout these operations the precise geometric history of each object is maintained. The foregoing approach to geometric modelling has particular advantages for many civil engineering structures, where excavation geometries are often defined using higher-order geometric primitives; for example, the layout of tunnels, surge chambers and caverns in a hydroelectric scheme, or the branching of transportation and conveyance tunnels. Maintaining smooth geometry also makes it easier to generate smooth, graded distributions of nodes on surfaces; in contrast, surfaces defined in a piecewise manner generally impose the additional constraint that nodes must lie on the vertices and edges of the geometry. Flexibility of node placement can be important when trying to optimize the mesh by grading node distributions, as will be illustrated in the following section involving case studies. The I-DEAS program uses a NURBS-based solid modeller, though it is also possible to use the program to define geometries using the more rudimentary techniques typically found in less sophisticated programs. In this way, several approaches to building geometries could be explored for the study described here. Particular modelling techniques are described in the next section. Far-field boundary conditions The effect of mesh truncation and approximated boundary conditions has been investigated analytically for some special cases [20]. For example, consider the case of a circular hole in an infinite homogeneous, isotropic and linear elastic medium under plane strain conditions. Let the approximation to this system be a thick-walled tube, such that the outer radii of the tube represents the truncation of the "mesh" (Fig. 3). Either zero-displacement or in situ stress boundary conditions may be applied to this outer surface. The constitutive, equilibrium and strain~lisplacement equations can be solved for simultaneously subject to these boundary conditions and the results compared with results from the Kirsch solution to assess the influence of the proximity of the outer boundary on the stresses and displacements everywhere in the medium. For stress boundary conditions it is convenient to separate the in situ stress tensor into hydrostatic and pure shear components in order to facilitate solutions. For hydrostatic boundary (and initial in situ) stresses, it turns out that stresses and displacements everywhere in
593
b
Fig. 3. Geometry of the thick-walled tube model used to study mesh truncation effects.
the thick-walled tube are greater than those in an infinite medium (for corresponding locations) by a factor of:
32 f12_ 1'
(1)
where fl = 2, b is the distance to the outer boundary and a is the distance to the inner boundary. This result is plotted in Fig. 4. A heuristic commonly used in practical modelling is fl ,~ 5-7. A similar exercise for pure shear boundary conditions results in more complicated expressions. However, for fl ~ 7 the stresses are overpredicted by no more than about 10%, and displacements are overpredicted by no more than about 25% in proximity of the excavation (~ < ~ 2, where r is the distance into the medium). For zero-displacement boundary conditions, stresses and displacements are underpredicted. For /3 ~ 7 the predicted stresses are no less than 90% of those based on the Kirsch solution (for Poisson's ratios less than 0.4) and the predicted displacements are no less than 80% of those based on the Kirsch solution (for Poisson's ratios less than 0.4 and ~ < ~ 2). For 3-D problems, the solutions for a spherical cavity in an infinite medium and a thick-walled spherical shell can be compared. Under hydrostatic stress boundary conditions stresses everywhere in the thick-walled shell are greater than those in an infinite medium by a factor of:
/33
(2)
/33-- 1 '
3.0
All stresses and displacements at a distance r are greater
by a factor of ~ than the corresponding
stresses and displacement obtained using 2.5
132 t~ = [32_ 1
the Kitsch solution.
2.o
1.5
1.C
o
i
i
i
2
3
4
5
i
6
7
8
9
i0
Fig. 4. The influence of boundary proximity on stresses and displacements for hydrostatic stress boundary conditions.
594
GRABINSKY and CURRAN: MESH GENERATION FOR FEA
Note that this is a higher-order function than that given in equation (1), indicating a smaller influence of the proximity of the outer boundary for the 3-D problem as compared to the 2-D one. The displacements are overpredicted to a larger extent, though this may be neglected for fl ~ 7 and ~ < ~ 2. For zero-displacement boundary conditions and hydrostatic initial stress, stresses and displacements are underpredicted everywhere in the medium. For fl ~ 7 the predicted stresses are no less than about 98% of those based on the spherical cavity solution (for Poisson's ratios less than 0.4) and the predicted displacements are no less than about 90% of those based on the spherical cavity solution (for Poisson's ratios less than 0.4 and ~< ~2). A complete derivation of the equations for the foregoing examples and various resulting graphs can be found in Ref. [20]. In summary, for fl ~, 7, Poisson's ratios less than 0.4 and ~ < ~ 2 errors in predicted stresses and displacements due to "mesh" truncation will be no greater than about -I- 10%. Given the uncertainties normally associated with input values for geomechanical modelling parameters (strength, stiffness, in situ stresses, dependency on direction of loading and heterogeneity), errors of + 10% due to mesh truncation may generally be regarded as insignificant. Indeed, it is more often the case in rock mechanics modelling that a series of models are used to examine trends in behaviour, or the influence of various construction alternatives or changes in modelling assumptions on predicted rock mass response [21,22]. The above may be extended, at least conceptually, to other problem geometries by examining the predicted elastic stress distributions around different shaped openings and noting that, in general, the perturbations to the initial stress field at a distance fl ~ 7 are more extreme for a circular opening than for irregularly shaped openings (e.g. Appendix 3 in Ref. [23]). This assumes that the maximum "radius" of the excavation is used for a in fl = _b In the case of multiple openings it is usually a" conservative to assume that the openings interact, and the "excavation radius" then spans several adjacent openings. It is considered in this work, then, that 3-D meshes truncated at a distance of approximately three "excavation diameters" (i.e. fl ,~, 7) from the underground openings will not produce results that are unduly influenced by the proximity of the outer extents of the model. It is reassuring to be able to demonstrate some theoretical basis for the well-used heuristic, fl ~ 5-7. Subsequent sections will demonstrate that the limit fl ~ 7 is easily and economically achieved using free-meshing techniques, and that mesh extents of fl >>7 can be obtained with very little additional computational expense. Mesh generation
Geometry created in I-DEAS using the program's geometric modeller can be transferred directly to the mesh generation module, thereby maintaining appro-
priate descriptions of curves bounding surfaces, surfaces bounding volumes, surfaces common to adjacent volumes, etc. If the geometry is defined manually. then surface and volume descriptions and connectivity must also be defined manually before the mesh can be generated. In practice it was often found that some user intervention was required even if full solid modelling was used, in order to gain better control over the distribution of elements. The geometry of the excavated regions is typically defined to be located at the center of a bounding cube. The side length of the cube is about seven times the "diameter" of the contained excavation or group of excavations, as described in the previous section. Zerodisplacement boundary conditions are prescribed on the faces of the cube, and internal stresses equivalent to the modelled in situ stress field are defined for each element generated. Local element sizes are given at each of the points used to define the geometry. Desired element sizes are smallest in areas of interest (also typically areas of relative high strain gradient), and largest at the exterior of the bounding cube. The specific choice of local element size depends on several factors, and experience must be used to attempt reasonable mesh sensitivity studies. Some examples are given in the following section. Experience indicates that element sizes should be chosen to maximize the possible rate of growth of element size away from the modelled excavated regions. This is consistent with: (i) the theoretical notion that modelling degrees of freedom (nodes) should concentrate in areas of higher strain gradient; and (ii) the trends exhibited by the Kirsch solution and hollow spherical inclusion solution which indicate high strain gradients (of order r -2 and r 3, respectively, for hydrostatic initial stresses) in proximity to the excavation boundary. Empirically, then, it would appear that these strain gradients decline faster than the node density required to accurately model the gradients can increase. This is considered to be the primary reason why meshes generated using mapped meshing techniques are so inefficient: the nodes of a mapped mesh cannot be appropriately distributed between areas of high and low strain gradient. Given that mapped meshing techniques predominantly have been used to generate finite element meshes, it may then be appreciated why such considerable efforts have been devoted to specialized elements and hybrid computational schemes which reduce the extent of the mapped finite element mesh. The I-DEAS program then attempts to automatically fill specified mesh volumes with tetrahedral elements. Robust algorithms for filling volumes of arbitrary shape have been developed only for tetrahedral elements thus far---efficient methods for filling volumes with brick elements is still an area of active research. Ten-noded tetrahedrons may be preferred over four-noded tetrahedrons, depending on the problem being analyzed (e.g. under conditions of plastic flow [24]). Examples of different volumes of arbitrary shape that have been meshed using the I-DEAS program are illustrated in the
GRABINSKY and CURRAN: MESH GENERATION FOR FEA next section. Such volumes are far more tedious to mesh using mapped meshing techniques, as it is necessary to subdivide the volume into geometric primitives. Further, the examples given in the next section illustrate how free meshing techniques provide easy control of element density in different parts of a given mesh volume. Using the foregoing approach to generate finite element meshes, it has been found that node concentrations are much higher in the center of the model than at the edges of the model. Therefore comparatively few elements are used to model the far-field rock mass, and the mesh may be extended to distances of/~ ~ 7 or more without significantly compromising the efficiency of the solution. As will be explained in the next section, solution accuracy was established by comparing finite element analysis results with 3-D boundary element results. As a result, the meshing techniques described in this section can be used to generate efficient and effective finite element meshes for the analysis of complex underground structures. CASE STUDIES Two civil engineering structures and two Canadian shield mines are used to illustrate the application of geometric modelling and free meshing techniques in practical geomechanical analysis. For each case finite element mesh sensitivity studies were carried out as follows: a boundary element mesh was generated based on the same model geometry used to generate the finite element mesh; successive boundary element meshes were then generated by recursive element subdivision, and solutions from these various meshes were compared to ensure the predicted elastic stress distribution was not influenced by boundary element mesh density; the elastic stress distribution from the boundary element analysis was then compared with the solution obtained using finite elements to ensure that the finite element mesh density was suitable. Boundary element results were used to establish the "benchmark" elastic stress distribution because, using the tools available for this study, bound-
595
ary element sensitivity studies were more easily conducted. In addition, it is reassuring when two different analysis methods produce equivalent analysis results. All boundary element modelling was performed using Examine 3D and Compute 3D [25-28]. Further mesh sensitivity studies were also performed for plasticity analyses, although a detailed description of such studies is beyond the scope of this paper.
Sudbury Neutrino Observatory The Sudbury Neutrino Observatory, at Creighton Mine, Sudbury, Ontario, is a 2 0 m dia. × 30m high cylindrical vault at a depth of 2000 m, constructed as part of an international scientific effort to study solar neutrinos. The main cavity is intersected by two service tunnels. The objective of the analysis was to investigate the distribution of potential rock mass failure around the cavity, particularly where the cavity is intersected by the access tunnels. Figure 5 shows the model geometry developed for the analysis. Since the cavern is axially symmetric, a profile of revolution was defined and swept through 360 ° to create the solid representing the cavern. Solid tunnels were constructed by extruding tunnel profiles. The final model geometry was developed by appropriately orienting the tunnels and using Boolean union operations together with the cavern. The program automatically determines curves of intersection, trims intersecting surfaces, deletes extraneous surfaces, and so on. Precise geometric definitions are maintained so that meshes of arbitrary geometric accuracy can subsequently be generated. The mesh generated on the basis of the model geometry is shown in Fig. 6. The element density at the surface of the excavation is shown in Fig. 6a, while Fig. 6b indicates how it is possible to use the free meshing technique to quickly grade the element size with distance from the excavation. The cavern in the model was extracted in several benches in order to investigate changes in the cavern stability with excavation stage. The largest analysis stage contained about 30,000
(b) (a)
D Fig. 5. Geometry of the Sudbury Neutrino Observatory:(a) profilesof the cavern (rotated) and intersectingtunnels (extruded); and (b) perspectiveview of the solid model built from geometry in (a), above.
596
(IRABINSK'~ and ( ' U R R A N :
(a)
MESH G E N E R A T I O N
FOR t t A
(b)
Fig. 6. Finite element model of the Sudbur~ Neutrino Observatory: (a) Perspective view of element surfaces forming the cavity and tunnel excavations: and (b) perspective view of solid mesh used for analysis. (Half the mesh shown: geometry from (a) is at the center of the model.)
elements and 5000 nodes, and required about 2 hr CPU time (on the SGI Indigo R3000) for an elastic analysis; the smallest stage contained about 10,000 elements and 2000 nodes and required about 12 rain CPU time for an elastic analysis. P o w e r cavern
The potential usefulness of geometric modelling is also demonstrated using a common civil engineering structure - - a n underground power cavern with intersecting penstocks. Figure 7 shows the geometry of the example cavern with penstocks. Advantage may be taken of the repeated symmetry of a long cavern with several parallel penstocks, where profile sections through the centerlines of the penstocks and halfway between parallel penstocks are treated as planes of symmetry. To construct this geometry, the cavern was modelled as an elliptical cylinder which was trimmed by rectangular surfaces to form the stepped base. Penstocks were then formed using circular cylinders which were intersected with the cavern. Again, the solid modelling functionality of the program computes all intersections
and new surfaces in a geometrically precise manner. Constructing the geometry for this model takes only a few minutes. The cavern and penstocks (which are now represented as a solid) are then "subtracted" from a larger rectangular block representing the host rock. The resulting solid is the host rock, and the hole in this solid has the shape of the power cavern and penstocks. The model geometry is brought into the meshing module of the program, local element sizes are defined and the mesh automatically generated. Figure 8 shows the finite element surfaces representing the wall of the power cavern and penstocks (the model in this figure has been reflected about one of the planes of symmetry, in order to show the whole penstocks). Figure 9 shows a solid view of the finite element model, and illustrates how elements are graded in size from smallest in the immediate area of the cavern-penstock intersection to largest at distances farthest from the cavern. The lateral boundaries of this model would be extended for analysis purposes, to ensure appropriate simulation of an infinite rock mass.
/
O
,
I
::"
0
Fig. 7. Power cavern model: geometry of the cavern and intersecting tunnels,
Fig. 8. Power cavern model: perspective view of element surfaces forming the excavated geometry,
GRABINSKY and C U R R A N :
Fig. 9. Power cavern model: perspective view of solid finite element mesh. (Front and back vertical planes are planes of symmetry.)
Ansil mine
Minnova Inc.'s Ansil mine, in northern Qu+bec, is a copper/zinc sulfide deposit located between depths of about 1200 and 1500m. The shape of the deposit is irregular, the top third forming a massive skewed tetrahedron which then tapers down into a more uniform tabular shape dipping about 50 °. The upper stopes were mined as vertical blocks, the largest being about 15 x 20 m wide × 60 m high. Longhole stopes were used in the lower, tabular portion, these stopes being about 15 m along the strike by 5 x 30 m high.
Fig. 10. Ansil mine model: perspective view of element surfaces forming the stopes. (Different stopes are shown in alternating shades of gray.)
MESH G E N E R A T I O N FOR FEA
597
The finite element model generated for the whole mine is shown in Fig. 10. Individual stopes can be identified by alternating shades of gray. Solid modelling was not used for this mine as the modelling tools available in I-DEAS Level V were not particularly well-suited to modelling certain stope geometries. Instead, geometry was created by entering coordinates for the stope edges by hand, defining curves through these points, defining areas (stope surfaces) by connecting curves, and then defining which surfaces enclosed individual stope volumes. The mine geometry was then enveloped with a surrounding box. Constructing this detailed model required about 3 days to idealize the geometry from the mine plans and 5 days to enter coordinate data and build the geometry; however, a given mesh based on this geometry could then be generated in less than 4hr, expediting many different studies. The mesh actually used for analysis consisted of approx. 8000 four-noded tetrahedrons and 1500 nodes. Each excavation step requires about 10 min of CPU time and 40 MB disk space. The element density used in this model was deemed appropriate for investigating mine-wide stress redistributions of interest when considering overall mine stability [22]. Interpretation of the results of the foregoing analysis indicated the potential for significant stability problems in the upper portion of the mine, and subsequently it was decided to perform a more detailed analysis of the area. The portion of the refined mesh corresponding to stope surfaces is shown in Fig. 11, and was based on a selected portion of the geometry constructed for the previous model. Figure 11 also illustrates the control over mesh density possible using the free meshing technique (note that mesh density is much higher at the top of the stopes, where results are of primary interest). This mesh consists of approx. 18,000 four-noded tetrahedrons and 3300 nodes. Each excavation step requires about 20 min of CPU time (on the SGI Indigo R3000) and 80 MB disk space. The element density used in this model was
Fig. 11. Ansil mine model: perspective view of element surfaces forming the excavated stopes (upper horizon analysis).
598
GRABINSKY and CURRAN: MESH GENERATION FOR FEA
Fig. 12. Strathcona mine model: perspectiveview of element surfaces forming the stopes.
deemed appropriate for investigating the stress distribution in the immediate stope hangingwall, which is of interest when conducting local stope stability analyses [22]. Such an analysis is documented elsewhere [29]. Non-linear analyses have also been conducted using the whole mine model, to investigate how the stress distribution in proximity of the mine may have been influenced by strength degradation of the orebody resulting from alteration and mineralization processes [20, 30]. This non-linear analysis was run overnight and required about 120 MB disk space. It should be noted that time required for non-linear analyses are dependent on many factors, and solution times can change by orders of magnitude depending on the selection of values for parameters such as material constants or equilibrium tolerances. Nevertheless, it would appear that the meshes developed using free meshing techniques have similar economic advantages for elastic-plastic analyses as for elastic analyses.
plans and develop an idealized geometry, and about half a day to generate the geometry by hand. Generating individual meshes based on this geometry takes only a few minutes. Because overall pillar stability is of concern here, the element density for this model was selected to provide accurate representation of stresses through the bulk of the pillar; higher element densities at the edges of the pillar would have to be used to investigate the local stability of stope walls adjacent to the pillar. The model used for analysis consisted of approx. 9400 four-noded tetrahedrons and 1700 nodes. Elastic analyses require a little more than l0 min of CPU time; plasticity analyses can require several orders of magnitude more time, depending on the particulars of the analyses, but even the longest of these analyses execute as overnight jobs (on the SGI Indigo R3000).
Strathcona mine
Many researchers and designers have concluded that the finite element method is impractical for 3-D analysis of underground structures due to the need to construct a volume mesh within and between multiple arbitrarily shaped openings and to extend the mesh to a region remote from the openings (for example, see the conclusions presented in Ref. [4]). However, these conclusions should be reconsidered in light of recent advances in engineering software and computing hardware. In developing a 3-D model the user must idealize the geometry of the structure of interest and subsequently generate a computer model of this geometry. The computer model defines solid volumes, either explicitly (i.e. the program maintains a datastructure of relevant geometric information) or implicitly (i.e. the user must intervene to ensure a physically meaningful representation has been constructed). Some underground structures can be modelled appropriately using relatively rudimentary descriptions of surfaces. For example, many mine geometries may be satisfactorily idealized using piecewise linear surfaces. The geometries of other underground structures (notably many civil engineering structures) may be more conveniently constructed using advanced
Falconbridge's Strathcona mine, near Sudbury, Ontario, is a nickel deposit located between depths of about 500 and 1000m. The shape of the deposit is essentially tabular with widths of 15-50 m, dipping at about 45 °. However, a footwall portion to the east of the deposit enlarges forming a more massive zone. A combination of blasthole stope and cut and fill mining techniques have been used. The analysis conducted for this mine was concerned with the response of a regional pillar during late mine life. The pillar measured approx. 200 m along the strike x 60 × 50 m vertically. The finite element mesh used for the pillar analysis is shown in Fig. 12 (note that excavated stopes are shown as the solid objects; the "pillar" is therefore the unmined portion near the center of the model). The geometry of this mine was amenable to the solid modelling tools available in I-DEAS Level V. However, when the geometry was transferred to the meshing module of the program some user intervention was required in order to properly mesh the region between the two bodies shown in Fig. 12. Construction of the model required about 1 day to review mine
SUMMARY AND CONCLUSIONS
GRABINSKY and CURRAN:
MESH GENERATION FOR FEA
solid modelling techniques, as demonstrated by the examples used in this paper. Generalizing the results of analytic models of the effect of finite far-field boundaries indicates that mesh truncation effects will generally not be significant when the diametrical length of the finite element model is about seven times greater than the typical diametrical length of the modelled excavated regions. In this context, "significant" is relative to the other assumptions which must be made in geomechanics modelling, such as assumed values for in situ stresses, rock mass strength, etc. Although the geometric limits indicated above may appear severe, it has been found that these limits can be easily met and exceeded using free meshing techniques. Models may be generated using these techniques such that element density is much greater in the area of interest to design that it is at the model's periphery. Based on comparisons of results obtained using different analysis techniques, it appears that the gradation of element density does not adversely effect analysis results. Free meshing algorithms, however, are non-trivial and generally must be integrated with some form of solid modelling functionality. This makes current engineering software with free meshing capabilities expensive (although typically no more expensive than many of the current mine planning software tools). In addition, it may take several weeks or even months to learn how to fully use the modelling functionality provided by sophisticated programs. Therefore, use of these current tools could probably only be economically justified by specialist organizations such as large consulting firms or industrial technical centers. On the other hand, impressively sophisticated modelling capabilities are possible once investment has been made in learning these techniques. Several practical numerical stress analysis studies have been undertaken using solid modelling and free meshing tools. A few of these studies are mentioned in this paper. Most of these 3-D models consist of several tens-ofthousand elements, and solution times for linear analyses are in the range of 10-30min using current desktop workstations. It bears mentioning that the kind of workstation used in this study is increasingly commonplace at many mines and consulting engineers' offices. Non-linear analyses take longer due to the iterative nature of the solution techniques, but many of these analyses may be run overnight or on weekends. With respect to the economics of labor in the analysis process, it is the authors' experience that the amount of time spent idealizing geometry and generating meshes is negligible compared to the amount of time spent viewing and interpreting analysis results and subsequently preparing diagrams for reports and presentations to clients. It is concluded, therefore, that 3-D finite element analyses of underground structures should be regarded as technically feasible, although further software development to make these modelling techniques more accessible to the geotechnical designer is desirable. Older and more rudimentary mesh generation approaches imposed re-
599
strictions on representation of geometry and on solution efficiency, but these restrictions have been ameliorated by advances in engineering analysis software. The enhanced analysis capability afforded by newer modelling approaches then allows the designer to concentrate on the more salient aspects of rock mechanics in the design of underground structures. Acknowledgement--The financial support of the Natural Sciences and Engineering Research Council of Canada (NSERC) is gratefully acknowledged.
Accepted for publication 5 April 1993.
REFERENCES 1. Beer (3. 'Infinite domain' elements in finite element analysis of underground excavations. Int. J. Numer. Anal. Meth. Geomech. 7, 1-7 (1983). 2. Beer G. and Meek J. L. 'Infinite domain' elements. Int. J. Nermer. Meth. Engng. 17, 43-52 (1981). 3. Zienkiewicz O. C., Emson C. and Betess P. A novel boundary infinite element. Int. J. Numer. Meth. Engng. 19, 393-404 (1983). 4. Beer G. Implementation of combined boundary element-finite element analysis with applications in geomechanics. Developments in Boundary Element Methods (Banerjee P. K. and Watson J. O., Eds), Vol. 4, pp. 191-226 (1986). 5. Beer (3. Finite element, boundary element and coupled analysis of unbounded problems in elastostatics. Int. J. Numer. Meth. Engng. 19, 567-580 (1983). 6. Brady B. H. G. and Wassyng A. A coupled finite elementboundary element method of stress analysis. Int. J. Rock Mech. Min. Sci. & Geomech. Abstr. lg, 475-485 (1981). 7. Lorig L., Brady B. H. G. and Cundall P. A. Hybrid distinct element-boundary element analysis of jointed rock. Int. J. Rock Meth. Min. Sci. & Geomech. Abstr. 23, 303-312 (1986). 8. Pande (3. N., Beer G. and Williams J. R. Numerical Methods in Rock Mechanics, p.327. Wiley, Chichester, U.K. (1990). 9. Cavendish J. C., Field D. A. and Frey W. H. An approach to automatic three-dimensional finite element mesh generation. Int. J. Numer. Meth. Engng. 21, 329-347 (1985). 10. Schroeder W. J. and Shephard M. S. Geometry-based fully automatic mesh generation and the Delaunay triangulation. Int. J. Numer. Meth. Engng. 26, 2503 2515 (1988). 11. Shephard M. S. Approaches to the automatic generation and control of finite element meshes. Appl. Mech. Rev. 41, 169-185 (1988). 12. Shephard M. S., Gallagher R. H. and Abel J. F. The synthesis of near-optimum finite element meshes with interactive computer graphics. Int. J. Numer. Meth. Engng. 15, 1021-1039 (1980). 13. I-DEAS VI Users" Manual, Structural Dynamics Research Corporation, Milford, Ohio (1990). 14. ABAQUS Users' Manual, Version 4.9. Hibbitt, Karlsson & Sorensen, Inc., Pawtucket, RI (1989). 15. Shah S. A study of the behaviour of jointed rock masses. Ph.D. Thesis, p. 252. Dept Civil Engineering, University of Toronto
(1992). 16. OfoegbuG. I. and Curran J. H. Yieldingand damage of intact rock. Can. Geotech. J. 28, 503-516 (1991). 17. Ofoegbu G. I. Development and numerical implementationof a ubiquitous joints model. Dept Civil Engineering,Universityof Toronto (1993). (In preparation.) 18. Farin G. Curves and Surfaces for Computer Aided Geometric Design. A Practical Guide, p. 334. Academic Press, San Diego (1988). 19. Mfintyl/i M. An Introduction to Solid Modelling, p. 401. Computer Science Press, Rockville, MD (1988). 20. Grabinsky M. W. Quantifying rock stress and evaluating its significance in the geomechanical design of underground mines. Ph.D. Thesis, p. 376. Dept. Civil Engineering, University of Toronto (1992). 21. Starfield A. M. and Cundall P. A. Towards a methodology for rock mechanics modelling. Int. J. Rock. Mech. Min. Sci. & Geomech. Abstr. 25, 99-106 (1988).
600
GRABINSKY and CURRAN:
22. Hoek E., Grabinsky M. W. F. and Diederichs M. S. Numerical modelling for underground excavation design. Trans. Instn Min. Metall. (Sect. A: Min. Industry) 100, A22 A30 (1991). 23. Hoek E. and Brown E. T. Underground Excavations in Rock, p. 527. Insitution of Mining and Metallurgy, London (1980). 24. Nagtegaal J. C., Parks D. M. and Rice J. R. On numerically accurate finite element solution in the fully plastic range. Comput. Meth. Appl. Mech. Engng 4, 153-177 (1974). 25. Examine3D Users' Manual, Version 2.0. Dept Civil Engineering, University of Toronto (1993). 26. Shah S. and Curran J. H. Accuracy of the regularizing transformation method for l/r singularities. Int. J. Numerical Meth. Engng (1993). (In preparation.)
MESH GENERATION FOR FEA 27. Shah S. and Curran J. H. Hybrid integration schemes for evaluating non-singular integrals over planar triangular elements. Int J. Rock. Mech. Min. Sci. & Geomech. Abstr. (1993). (In preparation.) 28. Shah S. and Curran J. H. Performance of direct and iterative solvers for BEM systems of equations. Int. J. Numer. Anal. Meth. Geomech. (1993). (In preparation.) 29. Hutchinson D. J. and Grabinsky M. W. Back analysis of stope stability at Ansil mine using instrumentation data and numerical modelling. Rock Support in Mining and Underground Construction (Kaiser P. K. and McCreath D. R., Eds), pp. 167-176 (1992). 30. Grabinsky M. W., Curran J. H. and Bawden W. F. Stress determinations for underground geomechanical design. Can. Geotech. J. (submitted for publication).
0148-9062/93 $6.00 + 0.00 Copyright © 1993 Pergamon Press Lid
Printed in Great Britain. All rights reserved
Efficient Mesh Generation Procedures for Finite Element Analysis of Underground Structures M. W. F. GRABINSKYt J. H. CURRANt
The efficiency and accuracy of the finite element method & dependent on the volume mesh developed for analysis. This can be particularly significant for 3-D analyses of unbounded problems such as underground structures. Several developments in automatic mesh generation have been made in recent years and are increasingly found in commercially available software, although these programs target mainly mechanical engineering applications. The work described in this paper considers application of advanced mesh generation techniques to the analysis of underground structures. The functionality of the software used is described, and some general considerations for mesh development are given. Application of the software is then demonstrated using two mining and two civil engineering structures. It is shown that the newer approaches to mesh generation can provide significant advantages in terms of computational efficiency and modelling functionality. However, further software development to make these modelling techniques more accessible to the geotechnical designer is desirable.
INTRODUCTION Many commercially available finite element programs produce meshes using a technique known as "mapped meshing" in which a topologically regular set of elements is distorted to fit the problem geometry (Fig. 1). While this approach is advantageous from the programmer's point of view, it becomes awkward to apply where the problem geometry is irregular and/or complex. Further, meshes produced using mapping techniques are notoriously inefficient for modelling 3-D problems with "infinite" domains. Issues of computational efficiency have been addressed by developing special elements [1-3] and hybrid computational schemes where a finite element mesh (or similar domain mesh) is mated to a boundary element mesh [4-8]. (It should be noted that hybrid schemes, while theoretically attractive, do not necessarily provide clear computational advantages for all kinds of problems since the specialized solvers developed for finite element matrices are incompatible with matrices produced by boundary elements [28].) Nevertheless, the practitioner is still left with the problem of trying to generate a mesh within and about complex geometries, even if specialized elements or hybrid computational schemes are used. An alternative approach to enhancing the efficiency of finite element analyses lies in a different method of
generating the mesh. "Free meshing" is a term used to describe approaches in which nodes are inserted into the problem domain and then elements are connected to these nodes (Fig. 2). In the node insertion part of the process, it is not necessary for the nodes to be placed in any preordained topology. Various free meshing algorithms have been developed and are increasingly found in many software programs intended primarily for mechanical engineering applications [9-12]. The main advantages of the free meshing technique are that a mesh can
tDepartment of Civil Engineering,Universityof Toronto, Toronto, Ontario, Canada MSS IA4. 591
Fig. 1. Exampleof mapped mesh.
592
GRABINSKY and CURRAN:
MESH GENERATION FOR FEA MESH GENERATION USING FREE MESHING
Fig. 2. Example of free mesh.
be produced in volumes with very general geometries, and greater control is gained over the densities of the nodes in the model (and therefore the model's computational efficiency). The enhanced modelling functionality comes at the price of increased programming effort, though this is not of direct concern to the end user. A commercially available software package called I-DEAS [13] has been used to investigate how recently developed geometric modelling and mesh generation functionality can best be applied to underground geomechanical structures. The methods used to generate finite element meshes using I-DEAS and application of these methods to four case histories are described in this paper. The consequences of using different approaches (free meshing vs mapped meshing with special enhancements) to finite element analysis of underground structures is then considered. The finite element meshes generated using I-DEAS were analyzed using the ABAQUS program [14]. In addition to the geomechanical constitutive models available in ABAQUS (extended Drucker-Prager and Critical State models), custom constitutive models (Hoek-Brown [15], Elastic-Brittle-Damage [16] and ubiquitous joint [17] models) were developed and used. All of the work was carried out on a Silicon Graphics Inc. (SGI) Iris Indigo R3000. This machine has processing capabilities similar to an Intel 486/66 MHz-based personal computer, but is equipped with special video capabilities for real-time animation of 3-D models. The video capabilities greatly enhance the development of model geometries and meshes, and significantly reduce the amount of time required to display 3-D analysis results. tin the context of this discussion, geometric modelling includes all activities used to definethe geometryof the model. It is more usual to considergeometricmodellingas the mathematicaldescriptionof curves and surfaces, and solid modelling as the representation of physical solid objects (i.e. the ordering of curves and surfaces into bounded, sensible solids).
Recent approaches to mesh generation tend to consider goemetric modellingt and element discretization as separate activities. In this way, the discretization is developed on the basis of the underlying goemetric definition, and therefore many different discretizations of the same geometry may be developed independently of each other. This is different from the mesh generation techniques commonly employed in programs based on mapped meshing, where the geometry must be defined in terms of modelling primitives (cubes, wedges and tetrahedrons) which are subsequently subdivided to produce the elements. Hence, geometric modelling and mesh generation are considered separately below. A further issue requiring consideration is the model representation of far-field boundary conditions: where the mesh is truncated and what combination of loads and displacements are prescribed at these outer limits. It will be shown that this issue can be handled easily using free meshing techniques. Geometric modelling Geometric information may be associated with a variety of parameters relevant to rock mechanics modelling, for example: the surfaces of excavations or regions yet to be excavated; the contacts of different geomechanical domains; the locations of faults, shear zones and other significant discontinuities; the attitudes of joint sets, bedding planes and schistocity; the locations of instrumentation cells and anchor points; and the locations of microseismic events. Various CAD and mine planning software tools are widely available and can be used to great advantage in generating, displaying and interpreting such information. However, the geometric modelling functionality of these tools can vary significantly in sophistication, and care must be taken when attempting to use these tools to produce geometry as a basis for generating finite element meshes. Rudimentary geometric modelling tools typically allow the user to define a series of profiles, each profile being a polyline or a string of low-order curves, and to (semi)-automatically connect these profiles using triangular surfaces. The resulting surface which envelopes the profiles is then faceted, not smooth. Functionality may be added to enable the modelling of holes in an object, or various geometric bifurcations (branches, Ys, etc.). The program may be written to calculate volumetric information in an ad hoe manner as the user generates the geometry. Although the foregoing approach to geometric modelling is comparatively unsophisticated, it can adequately fulfill the functional requirements of many geomechanical design projects. For example, the surface of a geologic contact may be interpolated on the basis of a few point locations determined from a drilling program, and fitting anything more sophisticated than a piecewise linear surface through these points is generally unjustified. Similarly, higher-order curves are often used to define the desired excavation profiles of many civil
GRABINSKY and CURRAN: MESH GENERATION FOR FEA engineering structures, but the as-built profiles may vary from the design profile to the same degree as a piecewise linear approximation varies from the higher-order curve. The currently favored approach to geometric modelling involves the description of curved lines and surfaces using non-uniform rational B-splines (NURBS) and the building of objects using constructive solid modelling and/or boundary representation modelling [18, 19]. Various functionality is used to define NURBS-based objects from geometric primitives, extrusions or sweeps of 2-D profiles, skinning of a series of profiles, and so on. Boolean operations may then be used to combine the basic objects, thereby producing more complex objects. Throughout these operations the precise geometric history of each object is maintained. The foregoing approach to geometric modelling has particular advantages for many civil engineering structures, where excavation geometries are often defined using higher-order geometric primitives; for example, the layout of tunnels, surge chambers and caverns in a hydroelectric scheme, or the branching of transportation and conveyance tunnels. Maintaining smooth geometry also makes it easier to generate smooth, graded distributions of nodes on surfaces; in contrast, surfaces defined in a piecewise manner generally impose the additional constraint that nodes must lie on the vertices and edges of the geometry. Flexibility of node placement can be important when trying to optimize the mesh by grading node distributions, as will be illustrated in the following section involving case studies. The I-DEAS program uses a NURBS-based solid modeller, though it is also possible to use the program to define geometries using the more rudimentary techniques typically found in less sophisticated programs. In this way, several approaches to building geometries could be explored for the study described here. Particular modelling techniques are described in the next section. Far-field boundary conditions The effect of mesh truncation and approximated boundary conditions has been investigated analytically for some special cases [20]. For example, consider the case of a circular hole in an infinite homogeneous, isotropic and linear elastic medium under plane strain conditions. Let the approximation to this system be a thick-walled tube, such that the outer radii of the tube represents the truncation of the "mesh" (Fig. 3). Either zero-displacement or in situ stress boundary conditions may be applied to this outer surface. The constitutive, equilibrium and strain~lisplacement equations can be solved for simultaneously subject to these boundary conditions and the results compared with results from the Kirsch solution to assess the influence of the proximity of the outer boundary on the stresses and displacements everywhere in the medium. For stress boundary conditions it is convenient to separate the in situ stress tensor into hydrostatic and pure shear components in order to facilitate solutions. For hydrostatic boundary (and initial in situ) stresses, it turns out that stresses and displacements everywhere in
593
b
Fig. 3. Geometry of the thick-walled tube model used to study mesh truncation effects.
the thick-walled tube are greater than those in an infinite medium (for corresponding locations) by a factor of:
32 f12_ 1'
(1)
where fl = 2, b is the distance to the outer boundary and a is the distance to the inner boundary. This result is plotted in Fig. 4. A heuristic commonly used in practical modelling is fl ,~ 5-7. A similar exercise for pure shear boundary conditions results in more complicated expressions. However, for fl ~ 7 the stresses are overpredicted by no more than about 10%, and displacements are overpredicted by no more than about 25% in proximity of the excavation (~ < ~ 2, where r is the distance into the medium). For zero-displacement boundary conditions, stresses and displacements are underpredicted. For /3 ~ 7 the predicted stresses are no less than 90% of those based on the Kirsch solution (for Poisson's ratios less than 0.4) and the predicted displacements are no less than 80% of those based on the Kirsch solution (for Poisson's ratios less than 0.4 and ~ < ~ 2). For 3-D problems, the solutions for a spherical cavity in an infinite medium and a thick-walled spherical shell can be compared. Under hydrostatic stress boundary conditions stresses everywhere in the thick-walled shell are greater than those in an infinite medium by a factor of:
/33
(2)
/33-- 1 '
3.0
All stresses and displacements at a distance r are greater
by a factor of ~ than the corresponding
stresses and displacement obtained using 2.5
132 t~ = [32_ 1
the Kitsch solution.
2.o
1.5
1.C
o
i
i
i
2
3
4
5
i
6
7
8
9
i0
Fig. 4. The influence of boundary proximity on stresses and displacements for hydrostatic stress boundary conditions.
594
GRABINSKY and CURRAN: MESH GENERATION FOR FEA
Note that this is a higher-order function than that given in equation (1), indicating a smaller influence of the proximity of the outer boundary for the 3-D problem as compared to the 2-D one. The displacements are overpredicted to a larger extent, though this may be neglected for fl ~ 7 and ~ < ~ 2. For zero-displacement boundary conditions and hydrostatic initial stress, stresses and displacements are underpredicted everywhere in the medium. For fl ~ 7 the predicted stresses are no less than about 98% of those based on the spherical cavity solution (for Poisson's ratios less than 0.4) and the predicted displacements are no less than about 90% of those based on the spherical cavity solution (for Poisson's ratios less than 0.4 and ~< ~2). A complete derivation of the equations for the foregoing examples and various resulting graphs can be found in Ref. [20]. In summary, for fl ~, 7, Poisson's ratios less than 0.4 and ~ < ~ 2 errors in predicted stresses and displacements due to "mesh" truncation will be no greater than about -I- 10%. Given the uncertainties normally associated with input values for geomechanical modelling parameters (strength, stiffness, in situ stresses, dependency on direction of loading and heterogeneity), errors of + 10% due to mesh truncation may generally be regarded as insignificant. Indeed, it is more often the case in rock mechanics modelling that a series of models are used to examine trends in behaviour, or the influence of various construction alternatives or changes in modelling assumptions on predicted rock mass response [21,22]. The above may be extended, at least conceptually, to other problem geometries by examining the predicted elastic stress distributions around different shaped openings and noting that, in general, the perturbations to the initial stress field at a distance fl ~ 7 are more extreme for a circular opening than for irregularly shaped openings (e.g. Appendix 3 in Ref. [23]). This assumes that the maximum "radius" of the excavation is used for a in fl = _b In the case of multiple openings it is usually a" conservative to assume that the openings interact, and the "excavation radius" then spans several adjacent openings. It is considered in this work, then, that 3-D meshes truncated at a distance of approximately three "excavation diameters" (i.e. fl ,~, 7) from the underground openings will not produce results that are unduly influenced by the proximity of the outer extents of the model. It is reassuring to be able to demonstrate some theoretical basis for the well-used heuristic, fl ~ 5-7. Subsequent sections will demonstrate that the limit fl ~ 7 is easily and economically achieved using free-meshing techniques, and that mesh extents of fl >>7 can be obtained with very little additional computational expense. Mesh generation
Geometry created in I-DEAS using the program's geometric modeller can be transferred directly to the mesh generation module, thereby maintaining appro-
priate descriptions of curves bounding surfaces, surfaces bounding volumes, surfaces common to adjacent volumes, etc. If the geometry is defined manually. then surface and volume descriptions and connectivity must also be defined manually before the mesh can be generated. In practice it was often found that some user intervention was required even if full solid modelling was used, in order to gain better control over the distribution of elements. The geometry of the excavated regions is typically defined to be located at the center of a bounding cube. The side length of the cube is about seven times the "diameter" of the contained excavation or group of excavations, as described in the previous section. Zerodisplacement boundary conditions are prescribed on the faces of the cube, and internal stresses equivalent to the modelled in situ stress field are defined for each element generated. Local element sizes are given at each of the points used to define the geometry. Desired element sizes are smallest in areas of interest (also typically areas of relative high strain gradient), and largest at the exterior of the bounding cube. The specific choice of local element size depends on several factors, and experience must be used to attempt reasonable mesh sensitivity studies. Some examples are given in the following section. Experience indicates that element sizes should be chosen to maximize the possible rate of growth of element size away from the modelled excavated regions. This is consistent with: (i) the theoretical notion that modelling degrees of freedom (nodes) should concentrate in areas of higher strain gradient; and (ii) the trends exhibited by the Kirsch solution and hollow spherical inclusion solution which indicate high strain gradients (of order r -2 and r 3, respectively, for hydrostatic initial stresses) in proximity to the excavation boundary. Empirically, then, it would appear that these strain gradients decline faster than the node density required to accurately model the gradients can increase. This is considered to be the primary reason why meshes generated using mapped meshing techniques are so inefficient: the nodes of a mapped mesh cannot be appropriately distributed between areas of high and low strain gradient. Given that mapped meshing techniques predominantly have been used to generate finite element meshes, it may then be appreciated why such considerable efforts have been devoted to specialized elements and hybrid computational schemes which reduce the extent of the mapped finite element mesh. The I-DEAS program then attempts to automatically fill specified mesh volumes with tetrahedral elements. Robust algorithms for filling volumes of arbitrary shape have been developed only for tetrahedral elements thus far---efficient methods for filling volumes with brick elements is still an area of active research. Ten-noded tetrahedrons may be preferred over four-noded tetrahedrons, depending on the problem being analyzed (e.g. under conditions of plastic flow [24]). Examples of different volumes of arbitrary shape that have been meshed using the I-DEAS program are illustrated in the
GRABINSKY and CURRAN: MESH GENERATION FOR FEA next section. Such volumes are far more tedious to mesh using mapped meshing techniques, as it is necessary to subdivide the volume into geometric primitives. Further, the examples given in the next section illustrate how free meshing techniques provide easy control of element density in different parts of a given mesh volume. Using the foregoing approach to generate finite element meshes, it has been found that node concentrations are much higher in the center of the model than at the edges of the model. Therefore comparatively few elements are used to model the far-field rock mass, and the mesh may be extended to distances of/~ ~ 7 or more without significantly compromising the efficiency of the solution. As will be explained in the next section, solution accuracy was established by comparing finite element analysis results with 3-D boundary element results. As a result, the meshing techniques described in this section can be used to generate efficient and effective finite element meshes for the analysis of complex underground structures. CASE STUDIES Two civil engineering structures and two Canadian shield mines are used to illustrate the application of geometric modelling and free meshing techniques in practical geomechanical analysis. For each case finite element mesh sensitivity studies were carried out as follows: a boundary element mesh was generated based on the same model geometry used to generate the finite element mesh; successive boundary element meshes were then generated by recursive element subdivision, and solutions from these various meshes were compared to ensure the predicted elastic stress distribution was not influenced by boundary element mesh density; the elastic stress distribution from the boundary element analysis was then compared with the solution obtained using finite elements to ensure that the finite element mesh density was suitable. Boundary element results were used to establish the "benchmark" elastic stress distribution because, using the tools available for this study, bound-
595
ary element sensitivity studies were more easily conducted. In addition, it is reassuring when two different analysis methods produce equivalent analysis results. All boundary element modelling was performed using Examine 3D and Compute 3D [25-28]. Further mesh sensitivity studies were also performed for plasticity analyses, although a detailed description of such studies is beyond the scope of this paper.
Sudbury Neutrino Observatory The Sudbury Neutrino Observatory, at Creighton Mine, Sudbury, Ontario, is a 2 0 m dia. × 30m high cylindrical vault at a depth of 2000 m, constructed as part of an international scientific effort to study solar neutrinos. The main cavity is intersected by two service tunnels. The objective of the analysis was to investigate the distribution of potential rock mass failure around the cavity, particularly where the cavity is intersected by the access tunnels. Figure 5 shows the model geometry developed for the analysis. Since the cavern is axially symmetric, a profile of revolution was defined and swept through 360 ° to create the solid representing the cavern. Solid tunnels were constructed by extruding tunnel profiles. The final model geometry was developed by appropriately orienting the tunnels and using Boolean union operations together with the cavern. The program automatically determines curves of intersection, trims intersecting surfaces, deletes extraneous surfaces, and so on. Precise geometric definitions are maintained so that meshes of arbitrary geometric accuracy can subsequently be generated. The mesh generated on the basis of the model geometry is shown in Fig. 6. The element density at the surface of the excavation is shown in Fig. 6a, while Fig. 6b indicates how it is possible to use the free meshing technique to quickly grade the element size with distance from the excavation. The cavern in the model was extracted in several benches in order to investigate changes in the cavern stability with excavation stage. The largest analysis stage contained about 30,000
(b) (a)
D Fig. 5. Geometry of the Sudbury Neutrino Observatory:(a) profilesof the cavern (rotated) and intersectingtunnels (extruded); and (b) perspectiveview of the solid model built from geometry in (a), above.
596
(IRABINSK'~ and ( ' U R R A N :
(a)
MESH G E N E R A T I O N
FOR t t A
(b)
Fig. 6. Finite element model of the Sudbur~ Neutrino Observatory: (a) Perspective view of element surfaces forming the cavity and tunnel excavations: and (b) perspective view of solid mesh used for analysis. (Half the mesh shown: geometry from (a) is at the center of the model.)
elements and 5000 nodes, and required about 2 hr CPU time (on the SGI Indigo R3000) for an elastic analysis; the smallest stage contained about 10,000 elements and 2000 nodes and required about 12 rain CPU time for an elastic analysis. P o w e r cavern
The potential usefulness of geometric modelling is also demonstrated using a common civil engineering structure - - a n underground power cavern with intersecting penstocks. Figure 7 shows the geometry of the example cavern with penstocks. Advantage may be taken of the repeated symmetry of a long cavern with several parallel penstocks, where profile sections through the centerlines of the penstocks and halfway between parallel penstocks are treated as planes of symmetry. To construct this geometry, the cavern was modelled as an elliptical cylinder which was trimmed by rectangular surfaces to form the stepped base. Penstocks were then formed using circular cylinders which were intersected with the cavern. Again, the solid modelling functionality of the program computes all intersections
and new surfaces in a geometrically precise manner. Constructing the geometry for this model takes only a few minutes. The cavern and penstocks (which are now represented as a solid) are then "subtracted" from a larger rectangular block representing the host rock. The resulting solid is the host rock, and the hole in this solid has the shape of the power cavern and penstocks. The model geometry is brought into the meshing module of the program, local element sizes are defined and the mesh automatically generated. Figure 8 shows the finite element surfaces representing the wall of the power cavern and penstocks (the model in this figure has been reflected about one of the planes of symmetry, in order to show the whole penstocks). Figure 9 shows a solid view of the finite element model, and illustrates how elements are graded in size from smallest in the immediate area of the cavern-penstock intersection to largest at distances farthest from the cavern. The lateral boundaries of this model would be extended for analysis purposes, to ensure appropriate simulation of an infinite rock mass.
/
O
,
I
::"
0
Fig. 7. Power cavern model: geometry of the cavern and intersecting tunnels,
Fig. 8. Power cavern model: perspective view of element surfaces forming the excavated geometry,
GRABINSKY and C U R R A N :
Fig. 9. Power cavern model: perspective view of solid finite element mesh. (Front and back vertical planes are planes of symmetry.)
Ansil mine
Minnova Inc.'s Ansil mine, in northern Qu+bec, is a copper/zinc sulfide deposit located between depths of about 1200 and 1500m. The shape of the deposit is irregular, the top third forming a massive skewed tetrahedron which then tapers down into a more uniform tabular shape dipping about 50 °. The upper stopes were mined as vertical blocks, the largest being about 15 x 20 m wide × 60 m high. Longhole stopes were used in the lower, tabular portion, these stopes being about 15 m along the strike by 5 x 30 m high.
Fig. 10. Ansil mine model: perspective view of element surfaces forming the stopes. (Different stopes are shown in alternating shades of gray.)
MESH G E N E R A T I O N FOR FEA
597
The finite element model generated for the whole mine is shown in Fig. 10. Individual stopes can be identified by alternating shades of gray. Solid modelling was not used for this mine as the modelling tools available in I-DEAS Level V were not particularly well-suited to modelling certain stope geometries. Instead, geometry was created by entering coordinates for the stope edges by hand, defining curves through these points, defining areas (stope surfaces) by connecting curves, and then defining which surfaces enclosed individual stope volumes. The mine geometry was then enveloped with a surrounding box. Constructing this detailed model required about 3 days to idealize the geometry from the mine plans and 5 days to enter coordinate data and build the geometry; however, a given mesh based on this geometry could then be generated in less than 4hr, expediting many different studies. The mesh actually used for analysis consisted of approx. 8000 four-noded tetrahedrons and 1500 nodes. Each excavation step requires about 10 min of CPU time and 40 MB disk space. The element density used in this model was deemed appropriate for investigating mine-wide stress redistributions of interest when considering overall mine stability [22]. Interpretation of the results of the foregoing analysis indicated the potential for significant stability problems in the upper portion of the mine, and subsequently it was decided to perform a more detailed analysis of the area. The portion of the refined mesh corresponding to stope surfaces is shown in Fig. 11, and was based on a selected portion of the geometry constructed for the previous model. Figure 11 also illustrates the control over mesh density possible using the free meshing technique (note that mesh density is much higher at the top of the stopes, where results are of primary interest). This mesh consists of approx. 18,000 four-noded tetrahedrons and 3300 nodes. Each excavation step requires about 20 min of CPU time (on the SGI Indigo R3000) and 80 MB disk space. The element density used in this model was
Fig. 11. Ansil mine model: perspective view of element surfaces forming the excavated stopes (upper horizon analysis).
598
GRABINSKY and CURRAN: MESH GENERATION FOR FEA
Fig. 12. Strathcona mine model: perspectiveview of element surfaces forming the stopes.
deemed appropriate for investigating the stress distribution in the immediate stope hangingwall, which is of interest when conducting local stope stability analyses [22]. Such an analysis is documented elsewhere [29]. Non-linear analyses have also been conducted using the whole mine model, to investigate how the stress distribution in proximity of the mine may have been influenced by strength degradation of the orebody resulting from alteration and mineralization processes [20, 30]. This non-linear analysis was run overnight and required about 120 MB disk space. It should be noted that time required for non-linear analyses are dependent on many factors, and solution times can change by orders of magnitude depending on the selection of values for parameters such as material constants or equilibrium tolerances. Nevertheless, it would appear that the meshes developed using free meshing techniques have similar economic advantages for elastic-plastic analyses as for elastic analyses.
plans and develop an idealized geometry, and about half a day to generate the geometry by hand. Generating individual meshes based on this geometry takes only a few minutes. Because overall pillar stability is of concern here, the element density for this model was selected to provide accurate representation of stresses through the bulk of the pillar; higher element densities at the edges of the pillar would have to be used to investigate the local stability of stope walls adjacent to the pillar. The model used for analysis consisted of approx. 9400 four-noded tetrahedrons and 1700 nodes. Elastic analyses require a little more than l0 min of CPU time; plasticity analyses can require several orders of magnitude more time, depending on the particulars of the analyses, but even the longest of these analyses execute as overnight jobs (on the SGI Indigo R3000).
Strathcona mine
Many researchers and designers have concluded that the finite element method is impractical for 3-D analysis of underground structures due to the need to construct a volume mesh within and between multiple arbitrarily shaped openings and to extend the mesh to a region remote from the openings (for example, see the conclusions presented in Ref. [4]). However, these conclusions should be reconsidered in light of recent advances in engineering software and computing hardware. In developing a 3-D model the user must idealize the geometry of the structure of interest and subsequently generate a computer model of this geometry. The computer model defines solid volumes, either explicitly (i.e. the program maintains a datastructure of relevant geometric information) or implicitly (i.e. the user must intervene to ensure a physically meaningful representation has been constructed). Some underground structures can be modelled appropriately using relatively rudimentary descriptions of surfaces. For example, many mine geometries may be satisfactorily idealized using piecewise linear surfaces. The geometries of other underground structures (notably many civil engineering structures) may be more conveniently constructed using advanced
Falconbridge's Strathcona mine, near Sudbury, Ontario, is a nickel deposit located between depths of about 500 and 1000m. The shape of the deposit is essentially tabular with widths of 15-50 m, dipping at about 45 °. However, a footwall portion to the east of the deposit enlarges forming a more massive zone. A combination of blasthole stope and cut and fill mining techniques have been used. The analysis conducted for this mine was concerned with the response of a regional pillar during late mine life. The pillar measured approx. 200 m along the strike x 60 × 50 m vertically. The finite element mesh used for the pillar analysis is shown in Fig. 12 (note that excavated stopes are shown as the solid objects; the "pillar" is therefore the unmined portion near the center of the model). The geometry of this mine was amenable to the solid modelling tools available in I-DEAS Level V. However, when the geometry was transferred to the meshing module of the program some user intervention was required in order to properly mesh the region between the two bodies shown in Fig. 12. Construction of the model required about 1 day to review mine
SUMMARY AND CONCLUSIONS
GRABINSKY and CURRAN:
MESH GENERATION FOR FEA
solid modelling techniques, as demonstrated by the examples used in this paper. Generalizing the results of analytic models of the effect of finite far-field boundaries indicates that mesh truncation effects will generally not be significant when the diametrical length of the finite element model is about seven times greater than the typical diametrical length of the modelled excavated regions. In this context, "significant" is relative to the other assumptions which must be made in geomechanics modelling, such as assumed values for in situ stresses, rock mass strength, etc. Although the geometric limits indicated above may appear severe, it has been found that these limits can be easily met and exceeded using free meshing techniques. Models may be generated using these techniques such that element density is much greater in the area of interest to design that it is at the model's periphery. Based on comparisons of results obtained using different analysis techniques, it appears that the gradation of element density does not adversely effect analysis results. Free meshing algorithms, however, are non-trivial and generally must be integrated with some form of solid modelling functionality. This makes current engineering software with free meshing capabilities expensive (although typically no more expensive than many of the current mine planning software tools). In addition, it may take several weeks or even months to learn how to fully use the modelling functionality provided by sophisticated programs. Therefore, use of these current tools could probably only be economically justified by specialist organizations such as large consulting firms or industrial technical centers. On the other hand, impressively sophisticated modelling capabilities are possible once investment has been made in learning these techniques. Several practical numerical stress analysis studies have been undertaken using solid modelling and free meshing tools. A few of these studies are mentioned in this paper. Most of these 3-D models consist of several tens-ofthousand elements, and solution times for linear analyses are in the range of 10-30min using current desktop workstations. It bears mentioning that the kind of workstation used in this study is increasingly commonplace at many mines and consulting engineers' offices. Non-linear analyses take longer due to the iterative nature of the solution techniques, but many of these analyses may be run overnight or on weekends. With respect to the economics of labor in the analysis process, it is the authors' experience that the amount of time spent idealizing geometry and generating meshes is negligible compared to the amount of time spent viewing and interpreting analysis results and subsequently preparing diagrams for reports and presentations to clients. It is concluded, therefore, that 3-D finite element analyses of underground structures should be regarded as technically feasible, although further software development to make these modelling techniques more accessible to the geotechnical designer is desirable. Older and more rudimentary mesh generation approaches imposed re-
599
strictions on representation of geometry and on solution efficiency, but these restrictions have been ameliorated by advances in engineering analysis software. The enhanced analysis capability afforded by newer modelling approaches then allows the designer to concentrate on the more salient aspects of rock mechanics in the design of underground structures. Acknowledgement--The financial support of the Natural Sciences and Engineering Research Council of Canada (NSERC) is gratefully acknowledged.
Accepted for publication 5 April 1993.
REFERENCES 1. Beer (3. 'Infinite domain' elements in finite element analysis of underground excavations. Int. J. Numer. Anal. Meth. Geomech. 7, 1-7 (1983). 2. Beer G. and Meek J. L. 'Infinite domain' elements. Int. J. Nermer. Meth. Engng. 17, 43-52 (1981). 3. Zienkiewicz O. C., Emson C. and Betess P. A novel boundary infinite element. Int. J. Numer. Meth. Engng. 19, 393-404 (1983). 4. Beer G. Implementation of combined boundary element-finite element analysis with applications in geomechanics. Developments in Boundary Element Methods (Banerjee P. K. and Watson J. O., Eds), Vol. 4, pp. 191-226 (1986). 5. Beer (3. Finite element, boundary element and coupled analysis of unbounded problems in elastostatics. Int. J. Numer. Meth. Engng. 19, 567-580 (1983). 6. Brady B. H. G. and Wassyng A. A coupled finite elementboundary element method of stress analysis. Int. J. Rock Mech. Min. Sci. & Geomech. Abstr. lg, 475-485 (1981). 7. Lorig L., Brady B. H. G. and Cundall P. A. Hybrid distinct element-boundary element analysis of jointed rock. Int. J. Rock Meth. Min. Sci. & Geomech. Abstr. 23, 303-312 (1986). 8. Pande (3. N., Beer G. and Williams J. R. Numerical Methods in Rock Mechanics, p.327. Wiley, Chichester, U.K. (1990). 9. Cavendish J. C., Field D. A. and Frey W. H. An approach to automatic three-dimensional finite element mesh generation. Int. J. Numer. Meth. Engng. 21, 329-347 (1985). 10. Schroeder W. J. and Shephard M. S. Geometry-based fully automatic mesh generation and the Delaunay triangulation. Int. J. Numer. Meth. Engng. 26, 2503 2515 (1988). 11. Shephard M. S. Approaches to the automatic generation and control of finite element meshes. Appl. Mech. Rev. 41, 169-185 (1988). 12. Shephard M. S., Gallagher R. H. and Abel J. F. The synthesis of near-optimum finite element meshes with interactive computer graphics. Int. J. Numer. Meth. Engng. 15, 1021-1039 (1980). 13. I-DEAS VI Users" Manual, Structural Dynamics Research Corporation, Milford, Ohio (1990). 14. ABAQUS Users' Manual, Version 4.9. Hibbitt, Karlsson & Sorensen, Inc., Pawtucket, RI (1989). 15. Shah S. A study of the behaviour of jointed rock masses. Ph.D. Thesis, p. 252. Dept Civil Engineering, University of Toronto
(1992). 16. OfoegbuG. I. and Curran J. H. Yieldingand damage of intact rock. Can. Geotech. J. 28, 503-516 (1991). 17. Ofoegbu G. I. Development and numerical implementationof a ubiquitous joints model. Dept Civil Engineering,Universityof Toronto (1993). (In preparation.) 18. Farin G. Curves and Surfaces for Computer Aided Geometric Design. A Practical Guide, p. 334. Academic Press, San Diego (1988). 19. Mfintyl/i M. An Introduction to Solid Modelling, p. 401. Computer Science Press, Rockville, MD (1988). 20. Grabinsky M. W. Quantifying rock stress and evaluating its significance in the geomechanical design of underground mines. Ph.D. Thesis, p. 376. Dept. Civil Engineering, University of Toronto (1992). 21. Starfield A. M. and Cundall P. A. Towards a methodology for rock mechanics modelling. Int. J. Rock. Mech. Min. Sci. & Geomech. Abstr. 25, 99-106 (1988).
600
GRABINSKY and CURRAN:
22. Hoek E., Grabinsky M. W. F. and Diederichs M. S. Numerical modelling for underground excavation design. Trans. Instn Min. Metall. (Sect. A: Min. Industry) 100, A22 A30 (1991). 23. Hoek E. and Brown E. T. Underground Excavations in Rock, p. 527. Insitution of Mining and Metallurgy, London (1980). 24. Nagtegaal J. C., Parks D. M. and Rice J. R. On numerically accurate finite element solution in the fully plastic range. Comput. Meth. Appl. Mech. Engng 4, 153-177 (1974). 25. Examine3D Users' Manual, Version 2.0. Dept Civil Engineering, University of Toronto (1993). 26. Shah S. and Curran J. H. Accuracy of the regularizing transformation method for l/r singularities. Int. J. Numerical Meth. Engng (1993). (In preparation.)
MESH GENERATION FOR FEA 27. Shah S. and Curran J. H. Hybrid integration schemes for evaluating non-singular integrals over planar triangular elements. Int J. Rock. Mech. Min. Sci. & Geomech. Abstr. (1993). (In preparation.) 28. Shah S. and Curran J. H. Performance of direct and iterative solvers for BEM systems of equations. Int. J. Numer. Anal. Meth. Geomech. (1993). (In preparation.) 29. Hutchinson D. J. and Grabinsky M. W. Back analysis of stope stability at Ansil mine using instrumentation data and numerical modelling. Rock Support in Mining and Underground Construction (Kaiser P. K. and McCreath D. R., Eds), pp. 167-176 (1992). 30. Grabinsky M. W., Curran J. H. and Bawden W. F. Stress determinations for underground geomechanical design. Can. Geotech. J. (submitted for publication).