- Email: [email protected]
Two-scale modeling of jointed rock masses
ARTICLE IN PRESS
International Journal of Rock Mechanics & Mining Sciences 43 (2006) 426–436 www.elsevier.com/locate/ijrmms
Two-scale modeling of jointed rock masses Jeen-Shang Lina,, Cheng-Yu Kub a
Department of Civil and Environmental Engineering, University of Pittsburgh, Pittsburgh, PA 15261, USA b Geotechnical Research Center, SinoTech Consulting Engineers, Taipei, Taiwan Accepted 31 July 2005 Available online 28 September 2005
Abstract Numerical modeling of the mechanical behavior of jointed rock masses is a difficult task as the discontinuities not only require special modeling consideration, but also require different treatments depending upon the scale of problem involved. A framework is proposed here that is based upon a meshed based partition of unity method, also known as the numerical manifold method. Specifically, the discontinuities posed by joints or faults are divided into two groups. The primary discontinuity set is the one with a wider spacing and dominates the kinematics of the system. It is modeled explicitly as physical discontinuities. The secondary discontinuity set, characterized by high density and small spacing, is modeled as an equivalent continuum by incorporating joint compliance and strength characteristics. Following the presentation of the formulation, two examples are also provided at the end to highlight the ease with which the proposed approach may be used in mechanical analysis of a complex jointed rock system. r 2005 Elsevier Ltd. All rights reserved. Keywords: Joints; Scale; Manifold method; Partition of unity; Discrete element; Discontinuities; Stability; Slope; Tunnel
1. Background The mechanical behavior of a jointed rock mass is difficult to predict or model and has long been an important focus of research [1–4]. The difficulties arise mostly from the complexity posed by the discontinuities introduced by joints or faults. The complexity comes mainly from two sources. First, there is the problem of the geometrical and mechanical representation of discontinuities. Second, depending upon the scale of interest, a jointed rock mass may exhibit behavior extending from that of an intact rock to that of a near-homogeneous highly fractured medium [4]. Moreover, the density and pattern of joints also play a crucial role. Faced with these difficulties, engineers often attempt to resort to an empirical relationship such as the Hoek–Brown strength criterion [2] in quantifying the jointed rock mass properties. It is well recognized, however, that there are severe limitations on the applicability of empirical relationships because of the Corresponding author. Benedum Hall 937, University of Pittsburgh, Pittsburgh, PA 15261, USA. E-mail address: [email protected] (J.-S. Lin).
1365-1609/$ - see front matter r 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijrmms.2005.07.009
variability encountered and the simplification involved. Neither are there, as of now, versatile empirical means for projecting the jointed rock mass deformation characteristics. An alternative is to obtain a general assessment of jointed rock mass behavior through numerical means. Four approaches are often adopted in rock mechanics, namely the equivalent continuum approach, the continuum with joint interface approach, the discrete element approach, and the discrete–continuum approach. The equivalent continuum approach [5,6] modifies the rock constitutive law to include the mechanical effects of joints that are dense and follow a regular pattern. The continuum with joint interface approach [7–10] introduces discontinuous interfaces to model each joint in a rock mass. This approach is simple to implement into finite elements but its application is limited to cases where the geometry does not undergo substantial changes; particularly no significant gaps or separations are to take place across discontinuities under loading. A discrete element approach [11–13], on the other hand, models the kinematics across each joint explicitly and is a powerful tool if a rock mass is delimited into blocks by joints. When discontinuities do not intersect into blocks or when the density of the joints is high, the
ARTICLE IN PRESS J.-S. Lin, C.-Y. Ku / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 426–436 Primary joints Secondary joints
Fig. 1. Illustration of the two-scale concept.
discrete element approach becomes impractical. The discrete–continuum approach [14–17] encompasses the strength of both the continuum and the discrete approach; it has, however, in the past been used mostly as an extended discrete element method to raise the degrees of freedom of the discrete elements. Recent advances in discontinuity modeling, particularly with the introduction of the manifold method, provide an opportunity for visiting the problem with a new perspective. Specifically, a two-scale modeling concept is proposed here within the framework of the manifold method. The underlying thinking of the proposed approach is that among the discontinuities encountered, there are generally few prevailing sets that dominate the gross kinematics. These dominating discontinuities are often widely spaced and, within the region bounded by them, are the set of the much more narrowly spaced joints. As a result, these contained discontinuities are the sources of only minor movements and have less impact on global kinematics. With these considerations, this study categorizes the joints in a rock mass into two sets according to their impact on overall kinematics. They are denoted as the primary discontinuity set and the secondary discontinuity set. It is further proposed that the primary discontinuity set be modeled explicitly, whereas the secondary set be modeled through an equivalent continuum. This two-scale concept is illustrated in Fig. 1. This paper first presents an overview of the numerical manifold method, which is followed by a discussion of the modeling consideration. Their implications are explored in numerical investigations. Two example applications follow at the end. 2. An overview of the manifold method Over the past decade, many numerical schemes have evolved around the principle of partition of unity. It all
427
started with the diffusion element in 1992 [18] that was immediately followed by the element free Galerkin method [19] and a variety of different approaches such as the reproduction kernel method [20], h–p clouds [21], and the partition of unity element method [22]. All these methods share the common thread of being based upon the ‘‘partition of unity’’. In the field of geotechnical engineering, particularly in rock mechanics, a parallel development can also be traced. In 1992, Shi [14] has proposed the framework of ‘‘the manifold method’’ as an extension to an implicit discrete element method, the discontinuous deformation analysis, or DDA [13], so that each element can have many more degrees of freedom. In essence, the manifold method employs two sets of meshes in modeling a problem: mathematical mesh and physical mesh. A physical mesh defines the integration domain. A mathematical mesh is a simple template pattern often consisting of triangular or rectangular grids. The procedure is to first develop a physical portrait of a given problem. This portrait contains both problem boundaries and the discontinuities that are present. A mathematical mesh is then constructed using a template that often consists of triangles or rectangular grids. The template is only required to be large enough to cover every point of the physical portrait. The grid size in a template is dictated by the analysis precision requirements. The part of the template that does not intersect with the physical portrait is trimmed, and the remaining grid is the mathematical mesh. The overlap region of the mathematical mesh and physical portrait gives the physical mesh. Upon the resulting mathematical mesh, nodes and elements are built. Fig. 2 depicts, on the left, a simple triangular mesh after trimming out the excess nodes. Often those nodes and edges that lie outside the problem domain are not shown to avoid clutter as depicted on the right in Fig. 2. The nodes of the mathematical mesh are used to build generalized elements, while the physical mesh is used in defining the integration bound for each element. In the example depicted, an element denoted as element 3 is defined by three nodes—5, 8, and 9—in the mathematical mesh, while its integration domain is defined by 8-9-e-f-8 in the physical mesh. Because of this unique construct the manifold method is much more flexible than other discrete–continuum approaches when it comes to tackling discontinuities that do not delimit a continuum into discrete blocks. Fig. 3 shows three discontinuities in a single rock specimen and that the manifold method treats it no differently from any other problem: it has a rectangular boundary and three discontinuities in the physical portrait. This further allows the manifold method to serve as a framework for tracking cracks as they grow and even after they break a rock into many blocks [23,24]. Lin [12,25,26] has presented two interpretations of the manifold method: one in terms of finite elements and the other in terms of partition of unity. The partition of unity view is more fundamental as the manifold method is in nature a finite cover approximation. Associated with a
ARTICLE IN PRESS J.-S. Lin, C.-Y. Ku / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 426–436
428
Fig. 2. Two meshes used in the manifold method.
Fig. 3. Manifold meshes in the presence of non-delimiting discontinuities.
node, a, on a mathematical mesh is a finite cover, oa. When a triangular mathematical mesh is used, a cover oa is the union of all triangles sharing this a node. Subordinates to each cover is a partition of unity function, ja(x), that satisfies the following: N X
ja ðxÞ ¼ 1
for x 2 O,
(1)
a¼1
ja ðxÞ 2 C s0 ðoa Þ
for 1 a N; s 0,
(2)
where N is the total number of mathematical nodes, O denotes the problem domain, and C s0 ðoa Þ denotes the class of functions that is s differentiable and whose its support is inside oa. On each cover the unknown is represented by an _ approximation function u a ðxÞ, often as a polynomial as follows: _ u a ðxÞ
¼ pT ðxÞ aa
m X
pk ðxÞaak ,
(3)
k¼0
where m is the order of the polynomial, p(x) a vector of polynomial basis, and aa a vector of the coefficients. In this view, an element is an area overlapped by three covers
when triangular meshes are used. The solution approximation within an element, uh(x), becomes uh ðxÞ ¼
3 X
_
ji ðxÞ u i ðxÞ.
(4)
i¼1
Since the physical domain is not represented by the mathematical mesh, the integration of the discrete governing equation has to be modified. This is done by introducing a weighting function, W(x,y), to track the domain of integration. Using the stiffness matrix evaluation as an example, the interpolation function and their derivatives are expressed in terms of the mathematical mesh nodes just like finite elements, but the integration is evaluated on the physical element as follows: ZZ ½K e ¼ W ðx; yÞ½Bðx; yÞT ½Cðx; yÞ½Bðx; yÞ dA, (5) Oe
where the familiar finite element form is modified by a W(x,y), and W(x,y) is one inside the physical mesh of an element and zero elsewhere. With this modification, the integration domain may not be convex.
ARTICLE IN PRESS J.-S. Lin, C.-Y. Ku / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 426–436
3. Modeling the primary discontinuity set
Y
Primary discontinuity sets are explicitly modeled, namely, each of this discontinuity is traced and introduced into the physical mesh. The strength of the manifold method lies in its ingenious way of tackling discontinuities: a mathematical mesh does not have to conform to discontinuities. A cover becomes disjoint whenever it is cut by a discontinuity into parts where no path can connect them without crossing the discontinuity. If a cover becomes disjoint due to a discontinuity, a cover is partitioned and cover counts are increased. This results in an increase in the number of nodes. Details of this construct can be found in Lin [26]. The mechanical modeling of explicit discontinuities implements two requirements: first, the surfaces of discontinuities can come into contact but cannot overlap; second, each contact must not violate the Mohr–Coulomb failure criterion. Contact detection and modeling algorithms, as that developed for DDA, can be applied with slight modification. When two boundaries come into contact, two pairs of penalty springs, each consisting of one normal and one shear, are inserted at the opposing contact boundaries. The normal spring, which is stiff, keeps the contact from overlapping; the shear spring enforces the Mohr–Coulomb law. If the shear strength is reached, the spring yields, i.e., jks sj kn n tan f þ C
for n40,
429
(6)
where ks, kn are the shear and normal penalty spring stiffness, respectively, s and n are shear and normal deformations at the contact point, and C is the cohesion force derived from the contact area that is attributed to the contact point. A contact may take a form as a vertex against an edge, or as an edge against another edge. The former is a point contact that has zero contact area and thus has no cohesion. Cohesion only enters the picture when a contact is made in the latter form. 4. Modeling of the secondary discontinuity set The secondary discontinuity set is modeled by adopting an equivalent continuum formulation. Several approaches have been used in the modeling of a jointed system as a continuum, among them the compliance equivalence approach [27] and the multi-laminated approach [28]. This study adopts a hybrid of these two approaches: the deformation equivalence is established through the equivalent compliance model, whereas the strength equivalence is established through a simplified version of the multilaminate model. 4.1. Deformation considerations The compliance equivalent approach uses joint compliances as a basis for establishing an equivalent continuum model for jointed rock masses. This study adopts the
Joint set 1
n
s S1
S2
X Joint set 2 Fig. 4. Joint sets modeled by equivalent compliance.
compliance matrix approach of Amadei and Goodman [29]. For the present 2-D problem, the joint pattern under this formulation is limited to two sets of perpendicular joints as depicted in Fig. 4. The compliance matrix of an equivalent continuum is simply the sum of the joint compliance matrix and the intact rock compliance matrix. Selecting a local coordinate system (n,s) such that n and s correspond to the normal and shear direction, respectively, of one joint set, the compliance matrix of an equivalent continuum, [C]n,s, in the local (n,s) coordinate system under plane strain condition can be derived as follows: 21 3 1 n 0 E þ kn1 S 1 E 6 n 7 1 1 7, 0 ½C n;s ¼ 6 (7) E E þ kn2 S 2 4 5 1 1 1 þ þ 0 0 2G ks1 S 1 ks2 S 2 where E is the elastic modulus and n is the Poisson’s ratio of the intact rock, while kn , ks and S are the normal stiffness, the shear stiffness and the spacing , respectively, of a joint set. Subscript 1 denotes those associated with the first joint set and subscript 2 denotes the second set. Using the fact that the strain energy is invariant to coordinate transformation, this compliance matrix in terms of the global (x,y) coordinates, [C]x,y, can be found as ½Cx;y ¼ ½T½Cn;s ½T1 ,
(8)
where y is the angle between x and s axis and 2 3 cos2 y sin2 y 2 cos y sin y 6 7 ½T ¼ 4 sin2 y cos2 y 2 cos y sin y 5. cos y sin y
cos y sin y
2
2
cos y sin y
(9)
ARTICLE IN PRESS J.-S. Lin, C.-Y. Ku / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 426–436
4.2. Strength considerations To account for the shear strength of joints, we require that the shear stress on joints obeys the Mohr–Coulomb failure criterion. If the Mohr–Coulomb criterion is violated, the overshot stresses are redistributed through iteration. Along a joint, the Cauchy formula gives the normal and shear stresses as sn ¼ T n, (10) t where T is the stress tensor and n is the unit normal to the joint. The overshoot stress is obtained as follows: 8 < jtj s tan f þ c for sn 40; n j j Dt ¼ (11) :t for sn o0; where fj and cj represent the friction angle and the cohesion of joints, respectively. Here, we consider that the joint under tension has negligible shear strength. Moreover, fj and cj are binary functions of joint displacement. When the shear displacement of a joint exceeds a certain limit, fj and cj drop to a set of residual values. The model adopted here is similar to the ‘‘ubiquitous joint model’’ [30]. 5. Comparison between discrete and equivalent continuum approaches How does the equivalent continuum model so defined compare with the discrete formulation? This question is
addressed by comparing computational results from the continuum approach with those from DDA. Specifically, we looked at a simple rock specimen that contains one single joint and analyzed the problem in two different ways. In DDA, the joint was modeled explicitly and it delimited the sample into two simply deformable blocks. A simply deformable block is a block with constant strains throughout and thus having only six degrees of freedom. In the continuum formulation, the joint was treated implicitly but had the same parameters as the explicit joint. The sample in the continuum approach was covered by 944 triangular elements. In both approaches, two additional blocks were used to model the top and bottom plates of the test setup. The upper plate was the loading plate for the strain controlled test, and the base plate was fixed. The contact between the loading plates and the sample is smooth. The sample had a nominal width of 7.3 cm and a height of 20 cm, and an elastic modulus of 105 MPa. A Poisson ratio of 0.2 was used. The friction angle on the joint interface, fj, was set at 301, and the joint cohesion was 50 kPa. Presented here is the case where the joint has an inclination of 501. Results of the analysis are presented in Fig. 5, which shows that these two approaches gave almost identical results prior to failure. In the case of DDA, the block system became unstable under the strain-controlled test when its joint strength was reached. The bottom half of the sample started to slide toward the right for lack of lateral constraint, and the top half thus moved down without much resistance. This resulted in a substantial drop in the
0.25 Discrete Continuum
0.20
Analytical solution
Stress (MPa)
Continuum model
Discrete model Loading plate
0.15
7.3 cm
944 elements
0.10
4 blocks 20 cm
430
Joint
50 degrees
0.05
Base plate
0.00 0.00E+0
1.00E-6
2.00E-6 3.00E-6 Strain
4.00E-6
Fig. 5. Stress–strain curves obtained from the numerical tests.
5.00E-6
ARTICLE IN PRESS J.-S. Lin, C.-Y. Ku / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 426–436
3.0 Primary joint set only Primary and secondary joints secondary joints only
2.0 Stress (MPa)
measured loading from the top plate. The continuum model, on the other hand, first yielded at a stress slightly higher than the analytical strength of 0.1969 MPa, but gradually approached that value as the strain was increased. This is due to the fact that strains were not uniform throughout the continuum sample, and even though most elements yielded when the analytical strength was reached, elements that remained in the elastic state pushed its strength higher and as they entered failure state under further straining, the analytical strength was approached. Overall, these two results are rather close initially. The diverging results upon reaching failure reflect the impact of the kinematical constraints on the discrete system. This has an important implication: the equivalent continuum model works only if the discrete system it represents is properly constrained. This study thus uses the equivalent continuum representation only for the secondary joints as they are confined by the primary joints and thus have limited kinematics.
431
1.0
0.0 0.0E+ 0
1.0E -5
2.0E -5
3.0E -5
4.0E -5
5.0E -5
Fig. 7. Stress–strain curves obtained.
6. Example applications As an illustration of the two-scale analysis framework, two example applications are presented here. The first example demonstrates the use of the approach in obtaining mechanical properties of a sample with a fairly complex pattern of joints and discontinuities. The second example is a stability study of an excavation and tunneling through a jointed rock slope. 6.1. Modeling the mechanical behavior of a jointed rock We considered a sample which resembles somewhat the one posed in Fig. 1, except that here we have two explicit discontinuities cut through the sample as depicted in Fig. 6. Three numerical uniaxial tests were conducted. First, other than the two discontinuities present, the rock sample itself is deemed intact. Thus, this is the case with primary joints
Fig. 6. A manifold representation of the test sample.
only. Second, we considered that the underlying rock is also filled with secondary joints, namely, both the primary and secondary joints are present. Third, only secondary joints are present throughout the sample, and no explicit discontinuities are included. For the analysis, the intact rock has an E of 2700 MPa and a Poisson ratio of 0.2. For the explicit discontinuity, a friction angle, f, of 301, and a cohesion of 30 kPa were adopted. The equivalent model considers the secondary joints to have a joint inclination of 601 from the horizontal, and the joint friction and cohesion, fj and cj, are 301 and 400 kPa, respectively. The resulting stress–strain curves of these tests are presented in Fig. 7; the displacement and stress field are depicted in Fig. 8. Under uniaxial loading, the primary joint set, when present, affects only the pattern of displacement, but does not facilitate movement along them. Therefore, the sample in the first case only deviates slightly from the elastic behavior in contrast to an intact rock. In essence, only the global elastic modulus is reduced due to slight movement along the explicit discontinuities. With the additional secondary joints, the second case followed initially the first case in deformation, but yielded as failure took place along the secondary joints. For the last case, the sample behaved as an elasto-plastic material and maintained the elastic modulus of the intact rock before failure was reached. It is well-known that stress trajectories rotate in the presence of discontinuities [31], which form an important basis for a lower bound solution. The stress field obtained shows strong stress rotation across the explicit discontinuity in the first and second cases. In the third case, the stress field reflects a joint pattern in the continuum. This example demonstrates that the stress–strain behavior is a function of the discontinuities embedded in a sample, and its characteristics can be very difficult to
ARTICLE IN PRESS 432
J.-S. Lin, C.-Y. Ku / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 426–436
Fig. 8. Displacement and stress fields of various samples at the end of tests.
generalize without a proper analysis. This is particularly so since even under a uniaxial loading, a simple sample is capable of exhibiting rich behavior. The numerical testing scheme as presented here may potentially be a viable tool to guide engineering judgments when dealing with such a material.
6.2. A stability study of excavation and tunneling through a rock slope The second example analyzes a failure that took place during the construction of Hsintan No. 1 tunnel, which is
ARTICLE IN PRESS J.-S. Lin, C.-Y. Ku / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 426–436
433
Fig. 9. Plane view of the slope–tunnel system as designed.
one of the 24 tunnels constructed for a new freeway at northern Taiwan at the outskirts of the capital city, Taipei. The tunnel passed through a rock slope at the point where the slope reached a height of about 140 m. The tunnel was originally designed to have a net width of 14.2 m and an estimated total excavated cross sectional area of 157 m2. The highway itself consists of both northbound and southbound lanes: only the northbound lanes were designed to pass via tunneling while the southbound lanes were constructed using open cut by benching down the slope. Fig. 9 shows a plan view of the tunnel and the slope. The construction of tunneling and the slope excavation proceeded from opposite ends toward each other. When these two lanes met during the construction, slope instability was observed and a large tunnel deformation occurred. This case has previously been studied using DDA [32] considering bedding planes and joints to be discontinuities that cut the slope into a large number of blocks. Such a complete discrete analysis, however, was found to be difficult to conduct because the results are very sensitive to how blocks were formed. Often false failure mechanisms may be obtained, or actual failure mechanisms missed. With a large number of blocks, fictitious key blocks were often formed that prevented actual failure to develop and caused unrealistic high local stress concentrations. Furthermore, a completely discrete formulation is also computationally very demanding. Much of this can be attributed to the need of determining the contact forces that are transmitted to each block. This crucial step requires, in each time step, a search for each block’s neighbors of potential contacts, and a determination of the contact status. A contact may newly form or becomes disengaged within a time step, and this status can only be identified through iteration as change in the status of one contact may impact another. On the other hand, by substantially
reducing the number of discrete components, the proposed approach is able to avoid much of the aforementioned difficulties. The geological unit making up the site is known locally as the Nan-Kang formation, which is composed of massive sandstone (MSS) in the upper slope and of alternate thinbedded sandstone and shale (SS/SH) in the lower slope as shown in Fig. 10. The formation beds dip into the slope at a relatively high angle: the rock strata strike approximately parallel to the slope and dip 55–651 into the slope. The major weak planes controlling the slope stability appear to be the beddings and stress relief joints. Clay seams and slickensides are prevalent in these planes. The stress relief joints, which run roughly parallel to the slope surface, are well developed in the surface layers as revealed by the drilling cores and slope excavation. The joint spacing is in the range of 10–80 cm. Both of the highway lanes were located within the weak and weathered sandstone/shale alternations, and they run almost parallel to the strike of the bedding. From the geological conditions of the site as described above, this analysis divides the rock mass into four distinct zones as depicted in Fig. 11. They are zone (1) for fresh massive sandstone, zone (2) for well-jointed, highly weathered sandstone or shalestone on the slope surface, zone (3) for sandstone/shale interbeds, and zone (4) for weathered massive sandstone near the toe of the slope. The material to be excavated does not impact on the analysis, and zone (2) properties were used. Each zone is modeled with the equivalent continuum model with two joint sets. The elastic moduli of rock masses and joint spacing for different zones used are summarized in Table 1. The interface between the different zones, and the boundaries of excavation are treated as explicit discontinuities. The interface parameters used for these discontinuities were a friction angle of 261 and a cohesion of 20 kPa.
ARTICLE IN PRESS J.-S. Lin, C.-Y. Ku / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 426–436
434
Weathered MSS Stress relief joint Original Ground Surface
Massive SS
Excavated Slope Surface
Bedding
SS/SH Interbeds
Weathered MSS
Fig. 10. Rock properties and joint orientation at the studied site.
Fig. 11. Geologic zoning used in the analysis. Table 1 Rock properties and joint spacing used in the analysis Zone
Elastic modulus (kPa)
cj (kPa)
fj (1)
Joint spacing (m)
1 2 3 4
700,000 200,000 500,000 500,000
700 20 50 50
45 26 38 38
2.0 0.8 2.0 2.0
Based on the physical mesh, a mathematical mesh with 2858 elements was generated as shown in Fig. 12. The in-situ stresses before the construction started were obtained by gravity turn-on. This is followed by the removal of the excavated material. During the construction as the instability was observed, anchors were introduced. This study modeled each anchor as an elastic bar and the pre-stress as initial tension. A pre-stress of 600 kN was applied to each anchor.
ARTICLE IN PRESS J.-S. Lin, C.-Y. Ku / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 426–436
435
Table 2 Comparison of computed results with field data
This analysis Field observation
Fig. 12. Mathematical and physical meshes adopted for the analysis.
(a)
60 cm Displacement field around the tunnel
800 KPa (b)
Stress field around the tunnel
Fig. 13. Displacement and stress field after tunnel and slope excavation.
The analysis meshes are depicted in Fig. 12, where the physical portrait includes the anchors and material zones boundaries are presented in Fig. 11. Results of the analysis, as summarized in Fig. 13, showed that after the excavation was completed, the tunnel underwent large deformation. The displacement vector around the tunnel and on the excavated surface of the slope lies in the direction of nearly 45–551 downward from the horizontal. These displacement orientations are more or less parallel to the beddings of zone (2) and indicate that large movements also occurred along weak planes between two zones. Furthermore, the
Slope movement (cm)
Tunnel settlement (cm)
Anchor force (kN) No. LC3560
32 (V), 26 (H) 35 (V), 28 (H)
34 (V) 35 (V)
660 600
displacements near the surface are found to be larger than those beneath the highly weathered layer, which is consistent with what was observed in the field. The maximum principal stresses around the tunnel also became parallel to the perimeter of the tunnel. Moreover, the minor principal stress is relatively small, and a detrimental uniaxial stress field around the tunnel perimeter was formed making the rock susceptible to failure as its confinement was lost. Among the seven layers of anchors used for stabilizing the slope, only the force of anchor LC3560, or the second lowest layer, was available for comparison. The field instrumentation data also included the movement of the slope at mid-height, the tunnel settlement at crown. Table 2 provides a summary of the comparison between computation and observation. The agreement between these two sets of data appears to be good in view of limited instrumentation data. 7. Conclusions The recent advance in the modeling of discontinuities provides an impetus to explore the difficult issue of scales when modeling a jointed rock system. In this respect this study proposes to tackle discontinuities at two levels and treat them either explicitly or implicitly, or both. The potential of the proposed approach has been illustrated with case studies that are difficult to perform using traditional methods. This study clearly demonstrates that the proposed approach can serve as a useful framework for applications from estimating rock mass properties to engineering analysis and design. At the present time, the equivalent continuum model adopted only accommodates up to two sets of mutual orthogonal joints, and the constitutive models for joints are very basic. Much improvement is to be carried out in these areas. Moreover, to provide credence to the proposed approach, future research is also needed in calibrating the analysis with field observation and large-scale model tests. References [1] Goodman RE. Introduction to rock mechanics. New York: Wiley; 1980. [2] Brady BHG, Brown ET. Rock mechanics for underground mining. London: Chapman & Hall; 1993. [3] Pande GN, Beer G, Williams JR. Numerical methods in rock mechanics. New York: Wiley; 1990.
ARTICLE IN PRESS 436
J.-S. Lin, C.-Y. Ku / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 426–436
[4] Hoek E, Brown ET. Underground excavation in rock. London: Institute of Mining and Metallurgy; 1980. [5] Gerrard CM. Elastic models of rock masses having one, two and three sets of joints. Int J Rock Mech Min Sci Geomech Abstr 1982;19:15–23. [6] Amadei B. The influence of rock anisotropy on measurement of stresses in situ. Berlin: Springer; 1983. [7] Goodman RE, Taylor RL, Brekke TL. A model for the mechanics of jointed rock. J. Soil Mech Foundations Div ASCE 1968;94: 637–59. [8] Richard EG, Christopher SJ. Finite element analysis for discontinuous rocks. In: Desai CS, Christian JT, editors. Numerical methods in geotechnical engineering. New York: McGraw-Hill; 1977. p. 155–75. [9] Desai CS, Zamman MM, Lightner JG, Siriwardane HJ. Thinlayer element for interfaces and joint. Int J Numer Anal Methods Geomech 1984;8:19–43. [10] Belytschko T, Plesha M, Dowding CH. A computer method for the stability analysis of caverns in jointed rock. Int J Numer Anal Methods Geomech 1984;8:473–92. [11] Cundall PA, Hart RD. Development of generalized 2-D and 3-D distinct element programs for modeling jointed rock. ICG Report, US Army Engineering Waterways Experiment Station, published as Misc. Paper SL-85-1, 1983. [12] Cundall PA, Hart RD. Numerical modeling of discontinua. J Eng Comput 1992;9:101–13. [13] Shi G-H. Block system modeling by discontinuous deformation analysis. Boston: Computational Mechanics Publications; 1992. [14] Shi G-H. Modeling rock joints and blocks by manifold method. In: Proceedings of the 33rd US symposium rock mechanics, Santa Fe, Mexico, 1992 p. 639–48. [15] Munjiza A. The combined finite-discrete element method. New York: Wiley; 2004. [16] Cundall PA. UDEC—a generalized distinct element program for modelling jointed rock. Peter Cundall Associates, Report PCAR-180, European Research Office, US Army, Contract DAJA37-79-C0548, 1980. [17] Lin J-S. A unified framework for discrete and continuum analysis. In: Cook BK, Jensen RP, editors. Discrete element methods—numerical modeling of discontinua. Geotechnical Special Publication No. 117. New York: ASCE; 2002. p. 145–50.
[18] Nayroles B, Touzot G, Villon P. Generalizing the finite element method: diffuse approximation and diffuse elements. Comput Mech 1992;10:307–18. [19] Belytschko T, Krongauz Y, Organ D, Fleming M, Krysl P. Meshless methods : an overview and recent developments. Comput Methods Appl Mech Eng 1996;139:3–47. [20] Liu WK, Jun S, Zhang YF. Reproducing kernel particle method. Int J Numer Methods Eng 1995;20:1081–106. [21] Duarte CAM, Oden JT. An h-p adaptive method using clouds. Comput Methods Appl Mech Eng 1996;139:237–62. [22] Babus˘ ka I, Melenk JM. The partition of unity finite element method: basic theory and applications. Comput Methods Appl Mech Eng 1996;139:289–314. [23] Ku C-Y. Modeling of jointed rock masses based on the numerical manifold method. Doctoral dissertation, Department of Civil and Environmental Engineering, University of Pittsburgh, 2001. [24] Lin J-S, Ku C-Y. Modeling the evolution of slope failure a crack propagation problem. In: Proceedings of the 11th international congress on fracture, 2005, paper 3617. [25] Lin J-S. Continuous and discontinuous analysis using the manifold method. In: Proceedings of the working forum on the manifold of material analysis, waterways experiment station, US Army Corp of Engineers, 1995. p. 1–20. [26] Lin J-S. A mesh based partition of unity method. Comput Methods Appl Mech Eng 2003;192:1515–32. [27] Gerrard CM. Joint compliances as a basis for rock mass properties and the design of supports. Int J Rock Mech Min Sci Geomech Abstr 1982;19:285–305. [28] Zienkiewics OC, Pande GN. Time-dependent multi-laminated model of rocks—a numerical study of deformation and failure of rock masses. Int J Numer Anal Methods Geomech 1977;1:219–47. [29] Amadei B, Goodman RE. Formulation of complete plane strain for regulariy jointed rocks. In: Proceedings of the 22nd symposium on rock mechanics, Cambridge, MA, 1981. p. 245–51. [30] Itasca Consulting Group Inc. FLAC—fast Lagrangian analysis of continuum. User’s Manual, Minneapolis, MN, USA, 1995. [31] Parry RHG. Mohr circles, stress path and geotechniques. E & FN Spon; 1995. [32] Chen S, Chern JC, Koo C. Study on performance of tunnel near slope by DDA. In: Proceedings of the first international conference on the analysis of discontinuous deformation, Chung-li, Taiwan, 1995. p. 109–23.
International Journal of Rock Mechanics & Mining Sciences 43 (2006) 426–436 www.elsevier.com/locate/ijrmms
Two-scale modeling of jointed rock masses Jeen-Shang Lina,, Cheng-Yu Kub a
Department of Civil and Environmental Engineering, University of Pittsburgh, Pittsburgh, PA 15261, USA b Geotechnical Research Center, SinoTech Consulting Engineers, Taipei, Taiwan Accepted 31 July 2005 Available online 28 September 2005
Abstract Numerical modeling of the mechanical behavior of jointed rock masses is a difficult task as the discontinuities not only require special modeling consideration, but also require different treatments depending upon the scale of problem involved. A framework is proposed here that is based upon a meshed based partition of unity method, also known as the numerical manifold method. Specifically, the discontinuities posed by joints or faults are divided into two groups. The primary discontinuity set is the one with a wider spacing and dominates the kinematics of the system. It is modeled explicitly as physical discontinuities. The secondary discontinuity set, characterized by high density and small spacing, is modeled as an equivalent continuum by incorporating joint compliance and strength characteristics. Following the presentation of the formulation, two examples are also provided at the end to highlight the ease with which the proposed approach may be used in mechanical analysis of a complex jointed rock system. r 2005 Elsevier Ltd. All rights reserved. Keywords: Joints; Scale; Manifold method; Partition of unity; Discrete element; Discontinuities; Stability; Slope; Tunnel
1. Background The mechanical behavior of a jointed rock mass is difficult to predict or model and has long been an important focus of research [1–4]. The difficulties arise mostly from the complexity posed by the discontinuities introduced by joints or faults. The complexity comes mainly from two sources. First, there is the problem of the geometrical and mechanical representation of discontinuities. Second, depending upon the scale of interest, a jointed rock mass may exhibit behavior extending from that of an intact rock to that of a near-homogeneous highly fractured medium [4]. Moreover, the density and pattern of joints also play a crucial role. Faced with these difficulties, engineers often attempt to resort to an empirical relationship such as the Hoek–Brown strength criterion [2] in quantifying the jointed rock mass properties. It is well recognized, however, that there are severe limitations on the applicability of empirical relationships because of the Corresponding author. Benedum Hall 937, University of Pittsburgh, Pittsburgh, PA 15261, USA. E-mail address: [email protected] (J.-S. Lin).
1365-1609/$ - see front matter r 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijrmms.2005.07.009
variability encountered and the simplification involved. Neither are there, as of now, versatile empirical means for projecting the jointed rock mass deformation characteristics. An alternative is to obtain a general assessment of jointed rock mass behavior through numerical means. Four approaches are often adopted in rock mechanics, namely the equivalent continuum approach, the continuum with joint interface approach, the discrete element approach, and the discrete–continuum approach. The equivalent continuum approach [5,6] modifies the rock constitutive law to include the mechanical effects of joints that are dense and follow a regular pattern. The continuum with joint interface approach [7–10] introduces discontinuous interfaces to model each joint in a rock mass. This approach is simple to implement into finite elements but its application is limited to cases where the geometry does not undergo substantial changes; particularly no significant gaps or separations are to take place across discontinuities under loading. A discrete element approach [11–13], on the other hand, models the kinematics across each joint explicitly and is a powerful tool if a rock mass is delimited into blocks by joints. When discontinuities do not intersect into blocks or when the density of the joints is high, the
ARTICLE IN PRESS J.-S. Lin, C.-Y. Ku / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 426–436 Primary joints Secondary joints
Fig. 1. Illustration of the two-scale concept.
discrete element approach becomes impractical. The discrete–continuum approach [14–17] encompasses the strength of both the continuum and the discrete approach; it has, however, in the past been used mostly as an extended discrete element method to raise the degrees of freedom of the discrete elements. Recent advances in discontinuity modeling, particularly with the introduction of the manifold method, provide an opportunity for visiting the problem with a new perspective. Specifically, a two-scale modeling concept is proposed here within the framework of the manifold method. The underlying thinking of the proposed approach is that among the discontinuities encountered, there are generally few prevailing sets that dominate the gross kinematics. These dominating discontinuities are often widely spaced and, within the region bounded by them, are the set of the much more narrowly spaced joints. As a result, these contained discontinuities are the sources of only minor movements and have less impact on global kinematics. With these considerations, this study categorizes the joints in a rock mass into two sets according to their impact on overall kinematics. They are denoted as the primary discontinuity set and the secondary discontinuity set. It is further proposed that the primary discontinuity set be modeled explicitly, whereas the secondary set be modeled through an equivalent continuum. This two-scale concept is illustrated in Fig. 1. This paper first presents an overview of the numerical manifold method, which is followed by a discussion of the modeling consideration. Their implications are explored in numerical investigations. Two example applications follow at the end. 2. An overview of the manifold method Over the past decade, many numerical schemes have evolved around the principle of partition of unity. It all
427
started with the diffusion element in 1992 [18] that was immediately followed by the element free Galerkin method [19] and a variety of different approaches such as the reproduction kernel method [20], h–p clouds [21], and the partition of unity element method [22]. All these methods share the common thread of being based upon the ‘‘partition of unity’’. In the field of geotechnical engineering, particularly in rock mechanics, a parallel development can also be traced. In 1992, Shi [14] has proposed the framework of ‘‘the manifold method’’ as an extension to an implicit discrete element method, the discontinuous deformation analysis, or DDA [13], so that each element can have many more degrees of freedom. In essence, the manifold method employs two sets of meshes in modeling a problem: mathematical mesh and physical mesh. A physical mesh defines the integration domain. A mathematical mesh is a simple template pattern often consisting of triangular or rectangular grids. The procedure is to first develop a physical portrait of a given problem. This portrait contains both problem boundaries and the discontinuities that are present. A mathematical mesh is then constructed using a template that often consists of triangles or rectangular grids. The template is only required to be large enough to cover every point of the physical portrait. The grid size in a template is dictated by the analysis precision requirements. The part of the template that does not intersect with the physical portrait is trimmed, and the remaining grid is the mathematical mesh. The overlap region of the mathematical mesh and physical portrait gives the physical mesh. Upon the resulting mathematical mesh, nodes and elements are built. Fig. 2 depicts, on the left, a simple triangular mesh after trimming out the excess nodes. Often those nodes and edges that lie outside the problem domain are not shown to avoid clutter as depicted on the right in Fig. 2. The nodes of the mathematical mesh are used to build generalized elements, while the physical mesh is used in defining the integration bound for each element. In the example depicted, an element denoted as element 3 is defined by three nodes—5, 8, and 9—in the mathematical mesh, while its integration domain is defined by 8-9-e-f-8 in the physical mesh. Because of this unique construct the manifold method is much more flexible than other discrete–continuum approaches when it comes to tackling discontinuities that do not delimit a continuum into discrete blocks. Fig. 3 shows three discontinuities in a single rock specimen and that the manifold method treats it no differently from any other problem: it has a rectangular boundary and three discontinuities in the physical portrait. This further allows the manifold method to serve as a framework for tracking cracks as they grow and even after they break a rock into many blocks [23,24]. Lin [12,25,26] has presented two interpretations of the manifold method: one in terms of finite elements and the other in terms of partition of unity. The partition of unity view is more fundamental as the manifold method is in nature a finite cover approximation. Associated with a
ARTICLE IN PRESS J.-S. Lin, C.-Y. Ku / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 426–436
428
Fig. 2. Two meshes used in the manifold method.
Fig. 3. Manifold meshes in the presence of non-delimiting discontinuities.
node, a, on a mathematical mesh is a finite cover, oa. When a triangular mathematical mesh is used, a cover oa is the union of all triangles sharing this a node. Subordinates to each cover is a partition of unity function, ja(x), that satisfies the following: N X
ja ðxÞ ¼ 1
for x 2 O,
(1)
a¼1
ja ðxÞ 2 C s0 ðoa Þ
for 1 a N; s 0,
(2)
where N is the total number of mathematical nodes, O denotes the problem domain, and C s0 ðoa Þ denotes the class of functions that is s differentiable and whose its support is inside oa. On each cover the unknown is represented by an _ approximation function u a ðxÞ, often as a polynomial as follows: _ u a ðxÞ
¼ pT ðxÞ aa
m X
pk ðxÞaak ,
(3)
k¼0
where m is the order of the polynomial, p(x) a vector of polynomial basis, and aa a vector of the coefficients. In this view, an element is an area overlapped by three covers
when triangular meshes are used. The solution approximation within an element, uh(x), becomes uh ðxÞ ¼
3 X
_
ji ðxÞ u i ðxÞ.
(4)
i¼1
Since the physical domain is not represented by the mathematical mesh, the integration of the discrete governing equation has to be modified. This is done by introducing a weighting function, W(x,y), to track the domain of integration. Using the stiffness matrix evaluation as an example, the interpolation function and their derivatives are expressed in terms of the mathematical mesh nodes just like finite elements, but the integration is evaluated on the physical element as follows: ZZ ½K e ¼ W ðx; yÞ½Bðx; yÞT ½Cðx; yÞ½Bðx; yÞ dA, (5) Oe
where the familiar finite element form is modified by a W(x,y), and W(x,y) is one inside the physical mesh of an element and zero elsewhere. With this modification, the integration domain may not be convex.
ARTICLE IN PRESS J.-S. Lin, C.-Y. Ku / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 426–436
3. Modeling the primary discontinuity set
Y
Primary discontinuity sets are explicitly modeled, namely, each of this discontinuity is traced and introduced into the physical mesh. The strength of the manifold method lies in its ingenious way of tackling discontinuities: a mathematical mesh does not have to conform to discontinuities. A cover becomes disjoint whenever it is cut by a discontinuity into parts where no path can connect them without crossing the discontinuity. If a cover becomes disjoint due to a discontinuity, a cover is partitioned and cover counts are increased. This results in an increase in the number of nodes. Details of this construct can be found in Lin [26]. The mechanical modeling of explicit discontinuities implements two requirements: first, the surfaces of discontinuities can come into contact but cannot overlap; second, each contact must not violate the Mohr–Coulomb failure criterion. Contact detection and modeling algorithms, as that developed for DDA, can be applied with slight modification. When two boundaries come into contact, two pairs of penalty springs, each consisting of one normal and one shear, are inserted at the opposing contact boundaries. The normal spring, which is stiff, keeps the contact from overlapping; the shear spring enforces the Mohr–Coulomb law. If the shear strength is reached, the spring yields, i.e., jks sj kn n tan f þ C
for n40,
429
(6)
where ks, kn are the shear and normal penalty spring stiffness, respectively, s and n are shear and normal deformations at the contact point, and C is the cohesion force derived from the contact area that is attributed to the contact point. A contact may take a form as a vertex against an edge, or as an edge against another edge. The former is a point contact that has zero contact area and thus has no cohesion. Cohesion only enters the picture when a contact is made in the latter form. 4. Modeling of the secondary discontinuity set The secondary discontinuity set is modeled by adopting an equivalent continuum formulation. Several approaches have been used in the modeling of a jointed system as a continuum, among them the compliance equivalence approach [27] and the multi-laminated approach [28]. This study adopts a hybrid of these two approaches: the deformation equivalence is established through the equivalent compliance model, whereas the strength equivalence is established through a simplified version of the multilaminate model. 4.1. Deformation considerations The compliance equivalent approach uses joint compliances as a basis for establishing an equivalent continuum model for jointed rock masses. This study adopts the
Joint set 1
n
s S1
S2
X Joint set 2 Fig. 4. Joint sets modeled by equivalent compliance.
compliance matrix approach of Amadei and Goodman [29]. For the present 2-D problem, the joint pattern under this formulation is limited to two sets of perpendicular joints as depicted in Fig. 4. The compliance matrix of an equivalent continuum is simply the sum of the joint compliance matrix and the intact rock compliance matrix. Selecting a local coordinate system (n,s) such that n and s correspond to the normal and shear direction, respectively, of one joint set, the compliance matrix of an equivalent continuum, [C]n,s, in the local (n,s) coordinate system under plane strain condition can be derived as follows: 21 3 1 n 0 E þ kn1 S 1 E 6 n 7 1 1 7, 0 ½C n;s ¼ 6 (7) E E þ kn2 S 2 4 5 1 1 1 þ þ 0 0 2G ks1 S 1 ks2 S 2 where E is the elastic modulus and n is the Poisson’s ratio of the intact rock, while kn , ks and S are the normal stiffness, the shear stiffness and the spacing , respectively, of a joint set. Subscript 1 denotes those associated with the first joint set and subscript 2 denotes the second set. Using the fact that the strain energy is invariant to coordinate transformation, this compliance matrix in terms of the global (x,y) coordinates, [C]x,y, can be found as ½Cx;y ¼ ½T½Cn;s ½T1 ,
(8)
where y is the angle between x and s axis and 2 3 cos2 y sin2 y 2 cos y sin y 6 7 ½T ¼ 4 sin2 y cos2 y 2 cos y sin y 5. cos y sin y
cos y sin y
2
2
cos y sin y
(9)
ARTICLE IN PRESS J.-S. Lin, C.-Y. Ku / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 426–436
4.2. Strength considerations To account for the shear strength of joints, we require that the shear stress on joints obeys the Mohr–Coulomb failure criterion. If the Mohr–Coulomb criterion is violated, the overshot stresses are redistributed through iteration. Along a joint, the Cauchy formula gives the normal and shear stresses as sn ¼ T n, (10) t where T is the stress tensor and n is the unit normal to the joint. The overshoot stress is obtained as follows: 8 < jtj s tan f þ c for sn 40; n j j Dt ¼ (11) :t for sn o0; where fj and cj represent the friction angle and the cohesion of joints, respectively. Here, we consider that the joint under tension has negligible shear strength. Moreover, fj and cj are binary functions of joint displacement. When the shear displacement of a joint exceeds a certain limit, fj and cj drop to a set of residual values. The model adopted here is similar to the ‘‘ubiquitous joint model’’ [30]. 5. Comparison between discrete and equivalent continuum approaches How does the equivalent continuum model so defined compare with the discrete formulation? This question is
addressed by comparing computational results from the continuum approach with those from DDA. Specifically, we looked at a simple rock specimen that contains one single joint and analyzed the problem in two different ways. In DDA, the joint was modeled explicitly and it delimited the sample into two simply deformable blocks. A simply deformable block is a block with constant strains throughout and thus having only six degrees of freedom. In the continuum formulation, the joint was treated implicitly but had the same parameters as the explicit joint. The sample in the continuum approach was covered by 944 triangular elements. In both approaches, two additional blocks were used to model the top and bottom plates of the test setup. The upper plate was the loading plate for the strain controlled test, and the base plate was fixed. The contact between the loading plates and the sample is smooth. The sample had a nominal width of 7.3 cm and a height of 20 cm, and an elastic modulus of 105 MPa. A Poisson ratio of 0.2 was used. The friction angle on the joint interface, fj, was set at 301, and the joint cohesion was 50 kPa. Presented here is the case where the joint has an inclination of 501. Results of the analysis are presented in Fig. 5, which shows that these two approaches gave almost identical results prior to failure. In the case of DDA, the block system became unstable under the strain-controlled test when its joint strength was reached. The bottom half of the sample started to slide toward the right for lack of lateral constraint, and the top half thus moved down without much resistance. This resulted in a substantial drop in the
0.25 Discrete Continuum
0.20
Analytical solution
Stress (MPa)
Continuum model
Discrete model Loading plate
0.15
7.3 cm
944 elements
0.10
4 blocks 20 cm
430
Joint
50 degrees
0.05
Base plate
0.00 0.00E+0
1.00E-6
2.00E-6 3.00E-6 Strain
4.00E-6
Fig. 5. Stress–strain curves obtained from the numerical tests.
5.00E-6
ARTICLE IN PRESS J.-S. Lin, C.-Y. Ku / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 426–436
3.0 Primary joint set only Primary and secondary joints secondary joints only
2.0 Stress (MPa)
measured loading from the top plate. The continuum model, on the other hand, first yielded at a stress slightly higher than the analytical strength of 0.1969 MPa, but gradually approached that value as the strain was increased. This is due to the fact that strains were not uniform throughout the continuum sample, and even though most elements yielded when the analytical strength was reached, elements that remained in the elastic state pushed its strength higher and as they entered failure state under further straining, the analytical strength was approached. Overall, these two results are rather close initially. The diverging results upon reaching failure reflect the impact of the kinematical constraints on the discrete system. This has an important implication: the equivalent continuum model works only if the discrete system it represents is properly constrained. This study thus uses the equivalent continuum representation only for the secondary joints as they are confined by the primary joints and thus have limited kinematics.
431
1.0
0.0 0.0E+ 0
1.0E -5
2.0E -5
3.0E -5
4.0E -5
5.0E -5
Fig. 7. Stress–strain curves obtained.
6. Example applications As an illustration of the two-scale analysis framework, two example applications are presented here. The first example demonstrates the use of the approach in obtaining mechanical properties of a sample with a fairly complex pattern of joints and discontinuities. The second example is a stability study of an excavation and tunneling through a jointed rock slope. 6.1. Modeling the mechanical behavior of a jointed rock We considered a sample which resembles somewhat the one posed in Fig. 1, except that here we have two explicit discontinuities cut through the sample as depicted in Fig. 6. Three numerical uniaxial tests were conducted. First, other than the two discontinuities present, the rock sample itself is deemed intact. Thus, this is the case with primary joints
Fig. 6. A manifold representation of the test sample.
only. Second, we considered that the underlying rock is also filled with secondary joints, namely, both the primary and secondary joints are present. Third, only secondary joints are present throughout the sample, and no explicit discontinuities are included. For the analysis, the intact rock has an E of 2700 MPa and a Poisson ratio of 0.2. For the explicit discontinuity, a friction angle, f, of 301, and a cohesion of 30 kPa were adopted. The equivalent model considers the secondary joints to have a joint inclination of 601 from the horizontal, and the joint friction and cohesion, fj and cj, are 301 and 400 kPa, respectively. The resulting stress–strain curves of these tests are presented in Fig. 7; the displacement and stress field are depicted in Fig. 8. Under uniaxial loading, the primary joint set, when present, affects only the pattern of displacement, but does not facilitate movement along them. Therefore, the sample in the first case only deviates slightly from the elastic behavior in contrast to an intact rock. In essence, only the global elastic modulus is reduced due to slight movement along the explicit discontinuities. With the additional secondary joints, the second case followed initially the first case in deformation, but yielded as failure took place along the secondary joints. For the last case, the sample behaved as an elasto-plastic material and maintained the elastic modulus of the intact rock before failure was reached. It is well-known that stress trajectories rotate in the presence of discontinuities [31], which form an important basis for a lower bound solution. The stress field obtained shows strong stress rotation across the explicit discontinuity in the first and second cases. In the third case, the stress field reflects a joint pattern in the continuum. This example demonstrates that the stress–strain behavior is a function of the discontinuities embedded in a sample, and its characteristics can be very difficult to
ARTICLE IN PRESS 432
J.-S. Lin, C.-Y. Ku / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 426–436
Fig. 8. Displacement and stress fields of various samples at the end of tests.
generalize without a proper analysis. This is particularly so since even under a uniaxial loading, a simple sample is capable of exhibiting rich behavior. The numerical testing scheme as presented here may potentially be a viable tool to guide engineering judgments when dealing with such a material.
6.2. A stability study of excavation and tunneling through a rock slope The second example analyzes a failure that took place during the construction of Hsintan No. 1 tunnel, which is
ARTICLE IN PRESS J.-S. Lin, C.-Y. Ku / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 426–436
433
Fig. 9. Plane view of the slope–tunnel system as designed.
one of the 24 tunnels constructed for a new freeway at northern Taiwan at the outskirts of the capital city, Taipei. The tunnel passed through a rock slope at the point where the slope reached a height of about 140 m. The tunnel was originally designed to have a net width of 14.2 m and an estimated total excavated cross sectional area of 157 m2. The highway itself consists of both northbound and southbound lanes: only the northbound lanes were designed to pass via tunneling while the southbound lanes were constructed using open cut by benching down the slope. Fig. 9 shows a plan view of the tunnel and the slope. The construction of tunneling and the slope excavation proceeded from opposite ends toward each other. When these two lanes met during the construction, slope instability was observed and a large tunnel deformation occurred. This case has previously been studied using DDA [32] considering bedding planes and joints to be discontinuities that cut the slope into a large number of blocks. Such a complete discrete analysis, however, was found to be difficult to conduct because the results are very sensitive to how blocks were formed. Often false failure mechanisms may be obtained, or actual failure mechanisms missed. With a large number of blocks, fictitious key blocks were often formed that prevented actual failure to develop and caused unrealistic high local stress concentrations. Furthermore, a completely discrete formulation is also computationally very demanding. Much of this can be attributed to the need of determining the contact forces that are transmitted to each block. This crucial step requires, in each time step, a search for each block’s neighbors of potential contacts, and a determination of the contact status. A contact may newly form or becomes disengaged within a time step, and this status can only be identified through iteration as change in the status of one contact may impact another. On the other hand, by substantially
reducing the number of discrete components, the proposed approach is able to avoid much of the aforementioned difficulties. The geological unit making up the site is known locally as the Nan-Kang formation, which is composed of massive sandstone (MSS) in the upper slope and of alternate thinbedded sandstone and shale (SS/SH) in the lower slope as shown in Fig. 10. The formation beds dip into the slope at a relatively high angle: the rock strata strike approximately parallel to the slope and dip 55–651 into the slope. The major weak planes controlling the slope stability appear to be the beddings and stress relief joints. Clay seams and slickensides are prevalent in these planes. The stress relief joints, which run roughly parallel to the slope surface, are well developed in the surface layers as revealed by the drilling cores and slope excavation. The joint spacing is in the range of 10–80 cm. Both of the highway lanes were located within the weak and weathered sandstone/shale alternations, and they run almost parallel to the strike of the bedding. From the geological conditions of the site as described above, this analysis divides the rock mass into four distinct zones as depicted in Fig. 11. They are zone (1) for fresh massive sandstone, zone (2) for well-jointed, highly weathered sandstone or shalestone on the slope surface, zone (3) for sandstone/shale interbeds, and zone (4) for weathered massive sandstone near the toe of the slope. The material to be excavated does not impact on the analysis, and zone (2) properties were used. Each zone is modeled with the equivalent continuum model with two joint sets. The elastic moduli of rock masses and joint spacing for different zones used are summarized in Table 1. The interface between the different zones, and the boundaries of excavation are treated as explicit discontinuities. The interface parameters used for these discontinuities were a friction angle of 261 and a cohesion of 20 kPa.
ARTICLE IN PRESS J.-S. Lin, C.-Y. Ku / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 426–436
434
Weathered MSS Stress relief joint Original Ground Surface
Massive SS
Excavated Slope Surface
Bedding
SS/SH Interbeds
Weathered MSS
Fig. 10. Rock properties and joint orientation at the studied site.
Fig. 11. Geologic zoning used in the analysis. Table 1 Rock properties and joint spacing used in the analysis Zone
Elastic modulus (kPa)
cj (kPa)
fj (1)
Joint spacing (m)
1 2 3 4
700,000 200,000 500,000 500,000
700 20 50 50
45 26 38 38
2.0 0.8 2.0 2.0
Based on the physical mesh, a mathematical mesh with 2858 elements was generated as shown in Fig. 12. The in-situ stresses before the construction started were obtained by gravity turn-on. This is followed by the removal of the excavated material. During the construction as the instability was observed, anchors were introduced. This study modeled each anchor as an elastic bar and the pre-stress as initial tension. A pre-stress of 600 kN was applied to each anchor.
ARTICLE IN PRESS J.-S. Lin, C.-Y. Ku / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 426–436
435
Table 2 Comparison of computed results with field data
This analysis Field observation
Fig. 12. Mathematical and physical meshes adopted for the analysis.
(a)
60 cm Displacement field around the tunnel
800 KPa (b)
Stress field around the tunnel
Fig. 13. Displacement and stress field after tunnel and slope excavation.
The analysis meshes are depicted in Fig. 12, where the physical portrait includes the anchors and material zones boundaries are presented in Fig. 11. Results of the analysis, as summarized in Fig. 13, showed that after the excavation was completed, the tunnel underwent large deformation. The displacement vector around the tunnel and on the excavated surface of the slope lies in the direction of nearly 45–551 downward from the horizontal. These displacement orientations are more or less parallel to the beddings of zone (2) and indicate that large movements also occurred along weak planes between two zones. Furthermore, the
Slope movement (cm)
Tunnel settlement (cm)
Anchor force (kN) No. LC3560
32 (V), 26 (H) 35 (V), 28 (H)
34 (V) 35 (V)
660 600
displacements near the surface are found to be larger than those beneath the highly weathered layer, which is consistent with what was observed in the field. The maximum principal stresses around the tunnel also became parallel to the perimeter of the tunnel. Moreover, the minor principal stress is relatively small, and a detrimental uniaxial stress field around the tunnel perimeter was formed making the rock susceptible to failure as its confinement was lost. Among the seven layers of anchors used for stabilizing the slope, only the force of anchor LC3560, or the second lowest layer, was available for comparison. The field instrumentation data also included the movement of the slope at mid-height, the tunnel settlement at crown. Table 2 provides a summary of the comparison between computation and observation. The agreement between these two sets of data appears to be good in view of limited instrumentation data. 7. Conclusions The recent advance in the modeling of discontinuities provides an impetus to explore the difficult issue of scales when modeling a jointed rock system. In this respect this study proposes to tackle discontinuities at two levels and treat them either explicitly or implicitly, or both. The potential of the proposed approach has been illustrated with case studies that are difficult to perform using traditional methods. This study clearly demonstrates that the proposed approach can serve as a useful framework for applications from estimating rock mass properties to engineering analysis and design. At the present time, the equivalent continuum model adopted only accommodates up to two sets of mutual orthogonal joints, and the constitutive models for joints are very basic. Much improvement is to be carried out in these areas. Moreover, to provide credence to the proposed approach, future research is also needed in calibrating the analysis with field observation and large-scale model tests. References [1] Goodman RE. Introduction to rock mechanics. New York: Wiley; 1980. [2] Brady BHG, Brown ET. Rock mechanics for underground mining. London: Chapman & Hall; 1993. [3] Pande GN, Beer G, Williams JR. Numerical methods in rock mechanics. New York: Wiley; 1990.
ARTICLE IN PRESS 436
J.-S. Lin, C.-Y. Ku / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 426–436
[4] Hoek E, Brown ET. Underground excavation in rock. London: Institute of Mining and Metallurgy; 1980. [5] Gerrard CM. Elastic models of rock masses having one, two and three sets of joints. Int J Rock Mech Min Sci Geomech Abstr 1982;19:15–23. [6] Amadei B. The influence of rock anisotropy on measurement of stresses in situ. Berlin: Springer; 1983. [7] Goodman RE, Taylor RL, Brekke TL. A model for the mechanics of jointed rock. J. Soil Mech Foundations Div ASCE 1968;94: 637–59. [8] Richard EG, Christopher SJ. Finite element analysis for discontinuous rocks. In: Desai CS, Christian JT, editors. Numerical methods in geotechnical engineering. New York: McGraw-Hill; 1977. p. 155–75. [9] Desai CS, Zamman MM, Lightner JG, Siriwardane HJ. Thinlayer element for interfaces and joint. Int J Numer Anal Methods Geomech 1984;8:19–43. [10] Belytschko T, Plesha M, Dowding CH. A computer method for the stability analysis of caverns in jointed rock. Int J Numer Anal Methods Geomech 1984;8:473–92. [11] Cundall PA, Hart RD. Development of generalized 2-D and 3-D distinct element programs for modeling jointed rock. ICG Report, US Army Engineering Waterways Experiment Station, published as Misc. Paper SL-85-1, 1983. [12] Cundall PA, Hart RD. Numerical modeling of discontinua. J Eng Comput 1992;9:101–13. [13] Shi G-H. Block system modeling by discontinuous deformation analysis. Boston: Computational Mechanics Publications; 1992. [14] Shi G-H. Modeling rock joints and blocks by manifold method. In: Proceedings of the 33rd US symposium rock mechanics, Santa Fe, Mexico, 1992 p. 639–48. [15] Munjiza A. The combined finite-discrete element method. New York: Wiley; 2004. [16] Cundall PA. UDEC—a generalized distinct element program for modelling jointed rock. Peter Cundall Associates, Report PCAR-180, European Research Office, US Army, Contract DAJA37-79-C0548, 1980. [17] Lin J-S. A unified framework for discrete and continuum analysis. In: Cook BK, Jensen RP, editors. Discrete element methods—numerical modeling of discontinua. Geotechnical Special Publication No. 117. New York: ASCE; 2002. p. 145–50.
[18] Nayroles B, Touzot G, Villon P. Generalizing the finite element method: diffuse approximation and diffuse elements. Comput Mech 1992;10:307–18. [19] Belytschko T, Krongauz Y, Organ D, Fleming M, Krysl P. Meshless methods : an overview and recent developments. Comput Methods Appl Mech Eng 1996;139:3–47. [20] Liu WK, Jun S, Zhang YF. Reproducing kernel particle method. Int J Numer Methods Eng 1995;20:1081–106. [21] Duarte CAM, Oden JT. An h-p adaptive method using clouds. Comput Methods Appl Mech Eng 1996;139:237–62. [22] Babus˘ ka I, Melenk JM. The partition of unity finite element method: basic theory and applications. Comput Methods Appl Mech Eng 1996;139:289–314. [23] Ku C-Y. Modeling of jointed rock masses based on the numerical manifold method. Doctoral dissertation, Department of Civil and Environmental Engineering, University of Pittsburgh, 2001. [24] Lin J-S, Ku C-Y. Modeling the evolution of slope failure a crack propagation problem. In: Proceedings of the 11th international congress on fracture, 2005, paper 3617. [25] Lin J-S. Continuous and discontinuous analysis using the manifold method. In: Proceedings of the working forum on the manifold of material analysis, waterways experiment station, US Army Corp of Engineers, 1995. p. 1–20. [26] Lin J-S. A mesh based partition of unity method. Comput Methods Appl Mech Eng 2003;192:1515–32. [27] Gerrard CM. Joint compliances as a basis for rock mass properties and the design of supports. Int J Rock Mech Min Sci Geomech Abstr 1982;19:285–305. [28] Zienkiewics OC, Pande GN. Time-dependent multi-laminated model of rocks—a numerical study of deformation and failure of rock masses. Int J Numer Anal Methods Geomech 1977;1:219–47. [29] Amadei B, Goodman RE. Formulation of complete plane strain for regulariy jointed rocks. In: Proceedings of the 22nd symposium on rock mechanics, Cambridge, MA, 1981. p. 245–51. [30] Itasca Consulting Group Inc. FLAC—fast Lagrangian analysis of continuum. User’s Manual, Minneapolis, MN, USA, 1995. [31] Parry RHG. Mohr circles, stress path and geotechniques. E & FN Spon; 1995. [32] Chen S, Chern JC, Koo C. Study on performance of tunnel near slope by DDA. In: Proceedings of the first international conference on the analysis of discontinuous deformation, Chung-li, Taiwan, 1995. p. 109–23.