- Email: [email protected]
Numerical determination of elastic compliance tensor of fractured rock masses by finite element modeling
International Journal of Rock Mechanics & Mining Sciences 70 (2014) 474–482
Contents lists available at ScienceDirect
International Journal of Rock Mechanics & Mining Sciences journal homepage: www.elsevier.com/locate/ijrmms
Technical Note
Numerical determination of elastic compliance tensor of fractured rock masses by finite element modeling Jian-Ping Yang a, Wei-Zhong Chen a,b,n, Yong-hao Dai a, Hong-Dan Yu a a State Key Laboratory of Geomechanics and Geotechnical Engineering, Institute of Rock and Soil Mechanics, Chinese Academy of Sciences, 430071 Wuhan, PR China b Research Center of Geotechnical and Structural Engineering, Shandong University, 250061 Jinan, PR China
art ic l e i nf o Article history: Received 18 March 2013 Received in revised form 5 December 2013 Accepted 22 June 2014
1. Introduction Discontinuities including joints, fissures and interface separation are pervasive in rock masses at or near the Earth’s surface. In the analysis of engineering problems dealing with rock masses, such as design of arch dams, bridge piers, tunnels and underground caverns, the discontinuities are of great importance because of their strong influence on the mechanical and hydraulic properties. The treatment of discontinuities in numerical analysis can be divided in two ways. First, the discontinuities are described explicitly if amount of fractures is not large in the analyzing domain. This approach is represented mostly by distinct element method (DEM) and discontinuous deformation analysis (DDA) [1,2], which considers fractured rock masses as an assemblage of intact rock blocks separated by fracture systems. Second, in areas the rock masses containing extremely large number of fractures compare to engineering domain, it is impractical and impossible to model every fracture, and furthermore, the mechanical behavior of rock masses seems to be continuous or pseudo-continuous. So the equivalent continuum method is often adopted, which assume the global behavior of rock masses can be described by the principles of continuum mechanics. Based on the fact that many engineering facilities are constructed in heavily fractured rock masses, the equivalent continuum method is widely used. In the continuum approach, the equivalent elastic mechanical parameters of fractured rock masses are important input
n Corresponding author at: State Key Laboratory of Geomechanics and Geotechnical Engineering, Institute of Rock and Soil Mechanics, Chinese Academy of Sciences, 430071 Wuhan, PR China. Tel./fax: þ 86 27 87199212. E-mail address: [email protected] (W.-Z. Chen).
http://dx.doi.org/10.1016/j.ijrmms.2014.06.007 1365-1609/& 2014 Elsevier Ltd. All rights reserved.
parameters in any analysis of rock masses deformation behavior. Considerable efforts have been made in the past for determining the equivalent elastic mechanical properties. The determination methods can be divided into field tests [3], empirical relationships [4–8], analytical methods [9–15] and numerical methods [16–18]. Min and Jing summarized the characteristics and shortcomings of these methods, and a methodology determining elastic compliance tensor for fractured rock masses using DEM was established in that study [19]. The main advantage of the developed methodology is that it can consider complex fracture system geometry and various constitutive relations of fractures and rock matrix, and their interactions. However, the DEM modeling model in the study ignored the dead-end of fractures and their effects on equivalent elastic behavior of fractured rock masses. The ignorance would cause change of stress distribution and fracture network, which may induce great deviation in determination of compliance matrix. Within the framework established by Min and Jing [19], equivalent elastic mechanical parameters of fractured rock masses were studied by FEM (Finite Element Method) modeling in this paper. The advantage of the modeling is in developing a new mesh generation program for the assembly of intact rock and fractures, by which dead-end of fractures could be incorporated. Here, we briefly describe the methodology to determine elastic compliance matrix for fracture rock masse, and introduce the loading sets and corresponding numerical data processing method, then develop a new mesh generating technique which consider dead-end of fractures. Finally, we conduct a series of numerical simulations at different sizes and rotate angles to study scale effects and anisotropic characteristics of elastic constants of the studied fractured rock masses, and evaluate the applicability of continuum mechanics.
J.-P. Yang et al. / International Journal of Rock Mechanics & Mining Sciences 70 (2014) 474–482
475
2. Modeling approach using FEM code 2.1. Methodology for the determination of the linear elastic compliance tensor for fractured rock masses by FEM modeling The problem considered here is to estimate the equivalent compliance matrix for fractured rock masses whose long fractures are parallel in the direction of out-of-plane, z-direction. The matrix is assumed isotropic. Elastic properties of this kind of fractured rock masses can be investigated by two dimensional plane strain analysis. The linear elastic relation of anisotropic media can be expressed as:
εmn ¼ Smnpq σ pq 2
For plane strain problems, Eq. (1) can be reduced to 3 3 2σ 3 S12 S13 S16 xx 6σ 7 S22 S23 S26 7 yy 7 7 6 7 7 6 7 S32 S33 S36 5 6 4 σ zz 5 σ xy S62 S63 S66
2 εxx S11 6ε 7 6 S 6 yy 7 6 21 6 7 6 εzz 7 ¼ 6 4 5 4 S31 εxy S61
ð1Þ
ð2Þ
Methodologies for the determination of components of the compliance matrix in Eq. (2) using numerical methods have been studied [19]. For completeness, these methods are briefly summarized here. The elastic modulus and the Poisson’s effects in the z-direction are kept as intact rock since that there is no effect of the out-of-plane fractures. Considering the symmetry conditions, then S13 ¼S31 ¼S23 ¼S32 ¼ ν/E, and S33 ¼1/E. Since shear stress σxy does not affect the z-direction deformation, S36 and S63 can be assumed as zero. Therefore, three independent stress boundary conditions are sufficient to determine all components in the matrix shown in Eq. (2). Three loading sets (set 1 to set 3 in Fig. 1) were applied to determine the nine unknowns in Eq. (3), which are S11, S22, S12, S21, S26, S61, S62 and S66. Theoretically, S12 ¼S21, S16 ¼2S61 and S26 ¼2S62 because of the symmetry conditions, it is necessary to pointed out here that because the tensor shear strain εxy , instead of the engineering shear strain γ xy , is used in Eq. (2), then there exists a factor of 2 between S16, S26 and S61, S62. Components S61, S62 and S16, S26 represent the influence of normal stresses on shear strains and normal stresses on shear strains, respectively. Determination of S16 and S26 requires the calculation of averaged normal strain in x-and y-directions under pure shear stress (loading set 3 in Fig. 1), which is much more accurate and easier compared to the calculation of averaged shear strain under normal stress (loading set 1 and set 2 in Fig. 1). So, seven components-S11, S22, S12, S21, S16, S26 and S66 are determined by numerical experiments, and S61, S62 are deduced from S16, S26. Imposing εzz ¼ 0 and expressing σ zz in terms of σ xx and σ yy , Eq. (2) can be written as
εxx ¼ A11 σ xx þ A12 σ yy þ S16 τxy εyy ¼ A21 σ xx þ A22 σ yy þ S26 τxy εxy ¼ S61 σ xx þ S62 σ yy þ S66 τxy
ð3Þ
where Aij ¼ Sij
Si3 Sj3 ði; j ¼ 1; 2Þ S33
ð4Þ
Through imposing loading set 1, S11 and S21 can be obtained by substituting stress and the calculated strain into Eq. (3). As illustrated in Fig. 2(a), the normal strain εxx and εyy are calculated by dividing sample length L with deformation in corresponding directions, which are ΔL11 þ ΔL12 in the x-direction and ΔL21 þ ΔL22 in y-direction. ΔL11 is the average displacement of all nodes on left boundary. All the other deformations were calculated in the same way. Similarly, S22 and S12 can be obtained through imposing loading set 2. In theory, S21 equals S12 . However,
Fig. 1. Schematic diagram of four loading sets: (1) compression in x direction; (2) compression in y direction; (3) pure shear; (4) biaxial normal stress loading. (The first three loading sets are employed to calculate the compliance matrix of the studied model, while the fourth loading set is employed to calculate the shear modulus in 451 direction based on the current model).
∆L 12
∆L 11
∆L 22 u22 u11 A 11
y ∆L 21
x
u12 u21
L Fig. 2. Schematic diagram of (a) normal strain evaluation and (b) shear strain evaluation.
numerical simulation results would show a little difference between them. By imposing loading set 3, S16, S26 and S66 can be obtained. The shear strain εxy is calculated from the change in angle of adjacent perpendicular boundary as illustrated in Fig. 2(b). Angle change of each boundary is calculated from the area formed by deformed and original boundary. Numerical simulation results showed that angle change of every adjacent boundary differs a little from the others. Averaging method is adopted and εxy is evaluated by ðθ11 þ θ12 þ θ21 þ θ22 Þ=4. To obtain S16 and S26 , the average normal strains in x and y direction are calculated the way Fig. 2(a) illustrated. An additional biaxial loading set 4 (Fig. 1(d)) was applied to calculate equivalent shear modulus in the 451 direction compared to the current direction based on the current model. The biaxial normal stress creates pure shear stress in 451 direction whose value equals p. By calculating the average displacement of all boundary nodes, engineering shear strain γ in the 451 direction is easily obtained from Eq. (5).
π L ðΔL11 þ ΔL12 Þ γ ¼ 2arctan 2 L þ ðΔL21 þ ΔL22 Þ
ð5Þ
476
J.-P. Yang et al. / International Journal of Rock Mechanics & Mining Sciences 70 (2014) 474–482
Force boundary conditions were applied in numerical simulation rather than displacement boundary conditions which would induce additional forces caused by discontinuous deformation along fractures. Besides, the inertia relief technique in FEM program was used to avoid rigid body displacement, and the central point of the sample was chosen as reference fixed rotation point. 2.2. Procedure of directional deformation modulus analysis by FEM modeling While the equivalent continuum approach is employed for mechanical analysis in fractured rock masses, the equivalent elastic compliance tensor of the fractured rock masses is a basic input parameter, and it can be determined by the following steps: Step 1: investigate the distributions of fractures in site, eg. trace lengths, dip directions, fracture spacing, etc. Step 2: generate the fracture network models by Monte Carlo simulation method based on fracture distribution parameters. Step 3: generate the mesh of fractured rock masses based on the fracture network generated in step 2. Step 4: determine the elastic parameters of intact rock and fractures separately by laboratory experiments. Step 5: study the scale effect and anisotropic properties of equivalent elastic constants of fractured rock masses using FEM based on the mesh generated in step 3 and elastic parameters obtained in step4. First, choose a point as the center of study windows. Second, study the equivalent elastic constants at different scales increased step-by-step. Third, rotate the study windows incrementally for each scale and study the elastic constants in different directions. Finally, determine the equivalent elastic compliance matrix and corresponding REV if it exists.
2.3. Techniques for finite element mesh generation of fractured rock masses The generated mesh in step 3 of Section 2.2 is critical for FEM simulation results. In previous studies on elastic moduli and permeability dependence on stress of fractured rock masses, only hydraulic connected network are incorporated in the analysis model, and the dead-end of fractures which have no contribution to fluid transport are ignored [17,19]. The ignorance may alter stress distribution and fracture network as discussed before. In order to avoid these shortages, a new technique generating finite element mesh for fractured rock masses is developed, the detailed processes is described as follows: 1. Mesh the constructed fracture network model with triangular elements. Fracture lines inside the model must be seeded with nodes. 2. Suppose i and j are original nodes on fracture, and triangular elements related to them are E1 : ijk, E2 : ilj. First, copy a pair of nodes i1 , j1 compare to original nodes i, j. Positions of i1 , j1 are calculated according to the positions of i, j and the fracture width. Then, renumber the triangular elements as E1 : i1 j1 k, E2 : ilj, and generate the corresponding fracture element E3 : ijj1 i1 , as shown by Fig. 3(a). 3. On the position of two fracture intersection, another three nodes are needed in addition to the existed node. Positions of the three copied nodes are calculated according to the position of original node and the fracture widths of the two intersecting fractures. So the nodes of elements, which are connected to the intersection node, should be renumbered. In addition, four fractured elements are generated. Modifications of fracture
Fig. 3. Mesh technology of fractures and surrounding matrix.
intersection in Fig. 3(b) are: E1 : ijk-E1 : i2 j1 k; E2 : ikl-E2 : i2 kl1 ; E3 : ilm-E3 : i3 lm1 E4 : imn-E4 : imn; E5 : ino-E5 : ino; E6 : ioj-E6 : i1 o1 j E7 : li3 i2 l1 ; E8 : ioo1 i1 ; E9 : mii3 m1 ; E10 : i1 jj1 i2 4. The cross area of two intersecting fractures (area enclosed by ii1 i2 i3 in Fig. 3(b)) is left blank. Because the mechanical properties are different between normal direction and shear direction of fracture element, a fracture element generated on the cross area would prevent one of the intersecting fractures from sliding along it. However, along which fracture the rock masses would slide is dependent on the stress state. So the sliding direction can’t be preset and the “common region” of intersection fractures is left blank. 5. On the fracture end, one additional node is generated. Position of the new node is calculated by the position of the end node and the fracture width. Modification of the surrounding elements is shown in Fig. 3(c).
3. Verification of the proposed FEM modeling technique In order to verify the proposed modeling technique, a model with two orthogonal fractures sets is used for the comparison between the deformation modulus produced by FEM experiments and the closedform expression proposed by Amadai and Goodman [10]. The basic parameters used are listed as follows. For isotropic intact rock, the elastic modulus and Poisson’s ratio is 50 GPa and 0.25. For both fracture sets, the normal stiffness and shear stiffness are 50 GPa/m and 10 GPa/m, and the fracture spacing is 1.5 m, coupling effects between shear deformation and normal stress, and normal deformation and shear stress are not considered. Considering the geometric symmetry, the computational models are rotated in intervals of 151 from 01 to 451. Three loading sets described in Section 2.1 are applied to obtain the compliance matrix
J.-P. Yang et al. / International Journal of Rock Mechanics & Mining Sciences 70 (2014) 474–482
for 301 model. As shown in Fig. 4, both the numerical results calculated from rotated models in 01, 151, 301 and 451, and the red solid line produced by the compliance matrix derived from the 301 model, are agreed very well with the analytical results. To study the effects of S61 , S62 , S16 and S26 on global elastic properties, a comparison curve (the red dotted line in Fig. 4) is plotted where the four components are assumed to be zero. It can be seen that large
477
deviation against analytical results exists except in the direction of 301. Displacement vector maps (Fig. 5) also show obvious angel change of boundary lines under normal stresses, and elongation in x direction under shear stresses. Quantitative analysis shows that, under compressive stress 1 MPa in x direction, εxx ¼ 5.2 10 5, εxy ¼ 9.5 10 6, under shear stress 1 MPa, εxx ¼ εyy ¼ 2.3 10 5, εxy ¼4.7 10 5. The magnitudes of normal strain and shear strain are almost in the same order. The results prove strong coupling effects between normal stresses and shear strains, and between shear stresses and normal strains. From discussion above, it is reasonable to conclude that the proposed FEM modeling technique and compliance matrix computing method are valid and can be used for more complex random fracture models.
4. Results of scale dependency and tensor characteristics of elastic parameters for fractured rock masses Fig. 4. Verification of numerical modeling results against analytical solution. (The black solid line corresponds to the analytical results. The red solid line is produced from compliance matrix derived from 301 model, while the red dot line is produced from the same matrix where components of S16, S26, S61 and S62 are ignored. Points of symbols corresponds to numerical results of rotated FEM models). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
4.1. Data of fractures and mechanical parameters for the analysis From site investigation on central part of auxiliary tunnel B of JinPing hydropower station II, two sets of fractures are identified and the distribution of mean fracture spacing (distance of adjacent
Fig. 5. Displacement vector maps illustrating coupling behavior between normal stresses and shear strains, and between shear stresses and normal strains (unit: m). (a) Vector of total displacement under normal stress, (b) vector of displacement in x direction under normal stress, (c) vector of total displacement under shear stress and (d) vector of displacement in x direction under shear stress.
478
J.-P. Yang et al. / International Journal of Rock Mechanics & Mining Sciences 70 (2014) 474–482
Table 1 Input parameters for random network of fractures. Fractures set
1 2
Fracture spacing
Trace length
Dip direction
Distribution pattern
Mean/ m
Standard deviation/ m
Distribution pattern
Mean/ m
Standard deviation/ m
Distribution pattern
Mean/ m
Standard deviation/ m
Normal distribution Normal distribution
2.5
0.1
4
1
10
0.1
3
0.7
Normal distribution Normal distribution
150
2.0
Normal distribution Normal distribution
50
7
fracture centroids), fracture strike and mean trace length are listed in Table 1. The intact rock is considered isotropic linear elastic medium. Elastic modulus and Poisson’s ratio of intact rock are 50 GPa and 0.25 by laboratory test results. The constitutive relation for a fracture is expressed as: 9 2 9 8 38 K nn K ns K nt > > = = < Δτn > < Δδn > 7 Δτs ¼ 6 ð6Þ 4 K sn K ss K st 5 Δδs > > > > : Δτ ; K tn K ts K tt : Δδt ; t where the subscript ‘n’ represents the normal direction (i.e. perpendicular) to the fracture plane; ‘s’ and ‘t’represent two orthogonal directions on the fracture plane. Thus, Δτn is the normal stress, Δτs and Δτt are the shear stresses; Δδn is the normal displacement of the fracture, Δδs and Δδt are the shear displacements. K ij (i; j ¼ n; s; t) represents components of the fracture stiffness matrix: K nn represents the traction behavior of a fracture; K ss and K tt represent the shear behavior of a fracture; K ns and K nt represent coupling effect of shear displacement caused by normal stress; K sn and K tn represent coupling effect of normal displacement caused by shear stress; K st and K ts represent coupling effect between two shear stresses. Experiments indicate that the normal behavior and shear behavior of a fracture are often coupled [20], and the fracture stiffness is influenced by many factors [21,22] (particularly at large increment of load). The fracture stiffness may be treated as constants in a small increment of load. For a large increment of load, constant fracture stiffness can only be used as approximate estimation. This study neglects the coupling effects between shear behavior and normal behavior, and treats normal stiffness and shear stiffness of fractures as constants, which are 50 GPa/m and 10 GPa/m determined from small increment load tests. The coupling components and closure-variable stiffness could be incorporated into the FEM modeling easily. Fig. 6 shows the numerical simulation cases in this study. Model sizes range from 2 m to 16 m with 2 m interval, and angles rotated from 01 to 901 with 301 interval. Elastic properties of 01–901, 301–1201 and 601–1501 models are obtained from FEM modeling directly. In addition, the equivalent shear modulus of the fractured rock masses in 151, 451 and 751 directions are evaluated by models of 601, 01 and 301 accordingly.
90° 120°
60° 30°
150°
2m 4m 6m 8m 10 m
0°
16 m
Fig. 6. Study cases of deformation modulus of fractured rock masses.
gradient can be observed for 10 m model due to the collective contributions of large number of fractures (Fig. 7(a)). And this is the physical basis in describing macroscopic mechanical behavior of fractured rock masses by continuum mechanics. Fig. 7(b) shows the displacement contour of fracture network without considering fracture dead-end. A large number of fracture dead-ends are deleted compared with fracture network in Fig. 7(a), and the displacement fields are also totally different. Quantitative evaluation shows that, for model considering dead-end, εxx ¼ 4.15 10 5 and εxx ¼23.4 GPa, while for models without dead-end, εxx ¼ 2.77 10 5 and εxx ¼34.5 GPa. The model with dead-end is much softer because of the existence of more fractures. Fig. 8(a) and (b) show that stresses concentrate around fracture tips. Larger shear stresses along fractures in 901 model imply larger shear strains compared with that of 601 model. Consequently, greater displacement will be induced for 901 model. That is the reason the equivalent deformation modulus is the lowest in 901 and highest in 601 (Fig. 9(a)).
4.2. Results and discussion 4.2.1. Deformation characteristics and stress distributions of FEM modeling results Fig. 7 shows the contours of directional displacements of 01–901 FEM model under 1 MPa compressive stresses (loading set 2). The zero displacement point lies in the center of the model since the balanced stresses were applied on the symmetrical two sides. Displacement jump often exists across fractures in local scale. However, a general trend of roughly uniform deformation
4.2.2. Scale effects and tensor characteristics of fractured rock masses It has been suggested that there are two prerequisites for the appropriateness of continuum approach in mechanical analysis of fractured rock masses [17,23]. 1. A REV must exist, then equivalent parameters can be derived from heterogeneous fractured rock masses. The values of
J.-P. Yang et al. / International Journal of Rock Mechanics & Mining Sciences 70 (2014) 474–482
479
Fig. 7. Contrast of displacement contours under compressive loading with and without fracture dead-ends (unit: m). (a) 10 m model including fracture dead-ends and (b) 10 m model without fracture dead-ends.
Fig. 8. Shear stress contours of 10 m FEM models (unit: Pa). (a) 1 MPa compressive stress in 601 direction (horizontal) and (b) 1 MPa compressive stress in 901 direction (vertical).
equivalent parameters change insignificantly as the sizes of rock masses exceed REV. This problem is referred as size effect. 2. The derived equivalent parameters can be expressed in tensor form for the usage of constitutive equations for continuum method. The problem is related to the anisotropic characteristic of fractured rock masses. Problem 1 is studied by learning the variation trends between equivalent elastic constants and model sizes. If all elastic constants are tending towards stability and the variations are within predefined acceptable ranges, the REV is assumed existing. The minimum size beyond which elastic constants change within predefined ranges is defined as REV. Variation of elastic constants (deformation modulus, shear modulus and Poisson’s ratio) is measured by VEC i ¼
jEC i EC i j EC i
ð7Þ
where EC i is elastic constants evaluated by FEM model with side length i, EC i is averaged elastic constants of FEM models having the same rotation angle and whose side length is large or equal than i. For problem 2, it is expected that a compliance matrix derived from one FEM model can predict elastic constants of other
directions in an acceptable error. Predicting error is evaluated by comparison of compliance matrix components between results from FEM modeling and predicted by compliance matrix of other directions. The predicting error is defined as FEM CM 1 jSij Sij j CE ¼ ∑ N i;j SCM ij
ð8Þ
where SFEM and SCM are components of compliance matrix in the ij ij same reference coordinate system, SFEM are the FEM modeling ij results, SCM are the results transformed from compliance matrix of ij other directions, and N denotes the number of components to be compared. 4.2.3. Results of the scale effects Fig. 9 presents the variation of the equivalent elastic constants (deformation modulus, shear modulus, Poisson’s ration) with increasing side length of FEM model. Fig. 9(a) exhibits that the equivalent deformation modulus increase with increasing side length and become stable after a certain size. Variation of modulus in all six directions is less than 1 GPa when the side length reaches 8 m. For modulus in all directions, the values of 2 m 2 m scale are the smallest in all sizes, because the mean lengths of the two fracture
480
J.-P. Yang et al. / International Journal of Rock Mechanics & Mining Sciences 70 (2014) 474–482
Results of equivalent shear modulus shown by Fig. 9(b) are consistent with the results of deformation modulus. Shear modulus are the lowest in 451 and 601, which are almost parallel with one fracture set. The large shear modulus in 01 and 151 are also related with directional properties of fracture sets. Compared Fig. 9(c) to (a), an opposite relationship between Poisson’s ratio and deformation modulus can be seen at both model scale and directional aspects. Whether it is a special or common relationship for fractured rock masses needs more study. Calculation indicates that the value of νxy =Ex almost equals νyx =Ey . This result proves the symmetry characteristics of compliance matrix from numerical aspect. Quantitative evaluation by Eq. (7) shows that variations of deformation modulus, shear modulus and Poisson’s ratio are below 5%, 10% and 5% in 8 m size. As the size increased to 10 m, the variations of the three studied elastic constants are all decreased below 5%.
Fig. 9. Deformation modulus, shear modulus and Poisson’s ratio with increasing side lengths of FEM models (a) deformation modulus (b) shear modulus (c) Poisson’s ratio.
sets are larger than the computational model size, and the computational model is interpenetrated by several single fractures. Results from Fig. 9(a) also show larger equivalent elastic modulus in direction 601 and 1501 than in direction 01 and 901 as the loading direction of 601 and 1501 is nearly parallel to one fracture set and perpendicular to the other one. It is similar to the verification one with two orthogonal fracture sets. So the fractured rock mass is relatively rigid in the parallel direction and soft in angle bisector direction.
4.2.4. Results of tensor characteristics Figs. 10 and 11 present comparisons of equivalent deformation modulus and shear modulus between FEM modeling and compliance matrix. In the two figures, the black, blue and magenta dash dot line represent directional deformation/shear modulus expected by compliance matrix derived from 01, 301, 601 FEM model respectively. The separate symbols are FEM modeling results. Plotting lines show that the all compliance matrixes express elastic anisotropy characteristics properly even at very small scale (say 2 m). This may be due to the four fractures inside 2 m model, as shown in Fig. 6, have almost the same directions with the global fractures. However, the values of both deformation and shear modulus are not stable in small sizes. As model size increases, the deviation between the three lines decrease, and all the three lines predict FEM modeling results well. Fig. 12 presents mean prediction error by Eq. (8) between results of FEM modeling and matrixes. Because six independent unknowns, S11 , S22 , S12 , S16 , S26 and S66 exist in compliance matrix, and S11 , S22 , S66 are thought as dominating ones [17], tensor predicting errors are evaluated by incorporating six and three components separately for comparison. Fig. 12(a) shows that the predicting error is under 10% in 6 m and 5% in 10 m when three components are compared. However, the predicting error increases greatly while six components are contained. The value is under 20% in 10 m and 12% in 12 m as shown by Fig. 12(b). Results exhibit that the compliance matrix derived from 601 model can predicts FEM modeling results in other directions very well, while the reverse is not. For example, the error between matrix of 601 model and 301 FEM results is 2.1% while the reverse is 12%. Meanwhile, the error between matrix of 601 model and 01 FEM results is 4.3%. This implies that the strategy of averaging compliance matrix [17] may be not the best choice, as perhaps one model could give better predicting results, as shown in the studied case. In this sense, the compliance matrix of 601 model is thought best representing the elastic behavior of the studied fracture rock masses. After translating to the reference coordinate system, the compliance matrix is: 0 1 42:92 17:84 5:00 5:04 B 17:82 46:72 5:00 8:68 C B C 12 S¼B ðpa1 Þ ð9Þ C 10 @ 5:00 5:00 20:00 0:00 A 2:54
4:32
0:00
43:37
5. Conclusions Scale effects and anisotropic characteristics of fractured rock masses are studied by a proposed FEM modeling technique
J.-P. Yang et al. / International Journal of Rock Mechanics & Mining Sciences 70 (2014) 474–482
481
Fig. 10. Comparison of deformation modulus between rotated FEM models and compliance matrixes. (a) 2 m, (b) 6 m, (c) 10 m and (d) 14 m. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Fig. 11. Comparison of shear modulus between rotated FEM models and compliance matrixes. (a) 2 m, (b) 6 m, (c) 10 m and (d) 14 m. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
which consideres dead-end fractures. The main conclusions are summarized below. A mesh generating technique is developed for fractured rock masses by which the dead-end of fractures can be incorporated into FEM modeling. In cases where the fracture spacing is close to or larger than the fracture length, plenty of fractures and fracture sections would become dead-end and ignoring of them may leads to wrong results. In the studied case, the FEM modeling result of 10 m model without dead-end almost overestimate the deformation modulus as much as 50% compared to the model incorporating fracture dead-end. The influence coefficients of normal stresses on shear strains, and shear stresses on normal strains, corresponding to S61 , S62 and S16 , S26 in compliance matrix, plays an important role in describing
elastic mechanical behavior of fractured rock masses properly. They can influence both the deformation modulus and anisotropic characteristic greatly. This implies that carefully measurements and analysis are required in laboratory and field tests of elastic constants of fractured rock masses. Both the variation of elastic constants on scale and tensor prediction error should be evaluated in determining the proper compliance matrix for fractured rock masses. For the studied case, variations of deformation modulus, shear modulus and Poisson’s ratio are below 5% as model size reaches to 10 m, and tensor prediction error is also below 5% when S11 , S22 , S66 are contained in prediction. However, after all the six independent componentsS11 , S22 , S12 , S16 , S26 and S66 in compliance matrix are incorporated to evaluate discrepancies between matrix and FEM modeling
482
J.-P. Yang et al. / International Journal of Rock Mechanics & Mining Sciences 70 (2014) 474–482
Foundation of China (Grant Nos. 51225902, 51274189 and 51109207). We would also like to express our sincere gratitude to the editor and the anonymous reviewers for their valuable comments, which have greatly improved this paper. References
Fig. 12. Mean prediction error of compliance matrix components between FEM results and predicted results from matrixes in other directions (a) mean prediction error of 3 components, S11 , S22 and S66 (b) mean prediction error of 6 components, S11 , S22 , S12 , S16 , S26 and S66 .
results, the predicting errors are 20% in 10 m and 12% after reaching 12 m. Considering the importance of S12 , S16 and S26 , the REV is determined as 12 m, which is about three times as much as mean trace length in the studied case (mean fracture spacing is smaller than mean trace length). Acknowledgements The authors gratefully acknowledge the support of the Chinese Fundamental Research (973) Program through the Grant No. 2013CB036006 and the support of the National Natural Science
[1] Hart RD. An introduction to distinct element modeling for rock engineering. In: Hudson JA, editor. Comprehensive rock engineering. Oxford: Pergamon Press; 1993. p. 245–61. [2] SHI GH. Discontinuous deformation analysis: a new numerical model for the statics and dynamics of block system. Ph.D. thesis. Univ of California, Berkeley; 1988. [3] Bieniawski ZT. Determining rock mass deformability-experience from case histories. Int J Rock Mech Min Sci 1978;15:237–47. [4] Palmstrom A, Singh R. The deformation modulus of rock masses: comparisons between in situ tests and indirect estimates. Tunnelling Underground Space Technol 2001;16:115–31. [5] Barton N. Some new Q-value correlations to assist in site characterisation and tunnel design. Int J Rock Mech Min Sci 2002;39:185–216. [6] Sonmez H, Gokceoglu C, Ulusay R. Indirect determination of the modulus of deformation of rock masses based on the GSI system. Int J Rock Mech Min Sci 2004;41:849–57. [7] Zhang L, Einstein HH. Using RQD to estimate the deformation modulus of rock masses. Int J Rock Mech Min Sci 2004;41:337–41. [8] Hoek E, Diederichs MS. Empirical estimation of rock mass modulus. Int J Rock Mech Min Sci 2006;43:203–15. [9] Singh B. Continuum characterization of jointed rock masses: Part I–The constitutive equations. Int J Rock Mech Min Sci 1973;10:311–35. [10] Amadei B, Goodman RE. A 3-D constitutive relation for fractured rock masses. In: Proceedings of the international symposium on the mechanical behavior of structured media. Ottawa; 18–21 May 1981. p. 249–68. [11] Gerrard CM. Elastic models of rock masses having one, two and three sets of joints. Int J Rock Mech Min Sci 1982;19:15–23. [12] Oda M. Fabric tensor for discontinuous geological materials. Soils Found 1982;22:96–108. [13] Oda M, Suzuki K, et al. Elastic compliance for rock-like materials with random cracks. Soils Found 1984;24:27–40. [14] Hu KX, Huang Y. Estimation of the elastic properties of fractured rock masses. Int J Rock Mech Min Sci 1993;30:381–94. [15] Huang TH, Chang CS, Yang ZY. Elastic moduli for fractured rock mass. Rock Mech Rock Eng 1995;28:135–44. [16] Kulatilake PHSW, Wang S, et al. Effect of finite size joints on the deformability of jointed rock in three dimensions. Int J Rock Mech Min Sci 1993;30:479–501. [17] Min KB, Jing L. Numerical determination of the equivalent elastic compliance tensor for fractured rock masses using the distinct element method. Int J Rock Mech Min Sci 2003;40:795–816. [18] Chen SH, He J, Shahrour I. Estimation of elastic compliance matrix for fractured rock masses by composite element method. Int J Rock Mech Min Sci 2012;49:156–64. [19] Min KB, Rutqvist J, et al. Stress-dependent permeability of fractured rock masses: a numerical study. Int J Rock Mech Min Sci 2004;41:1191–210. [20] Jing L, Nordlund E, Stephansson O. An experimental study on the anisotropy and stress-dependency of the strength and deformability of rock joints. Int J Rock Mech Min Sci 1992;29:535–42. [21] Barton N. The shear-strength of rock and rock joints 1976;13:255–79Int J Rock Mech Min Sci 1976;13:255–79. [22] Bandis SC, Lumsden AC, et al. Fundamentals of rock joint deformation. Int J Rock Mech Min Sci 1983;20:249–68. [23] Long JC, Remer JS, Wilson CR. Porous-media equivalents for networks of discontinuous fractures. Water Resour Res 1982;18:645–58.
Contents lists available at ScienceDirect
International Journal of Rock Mechanics & Mining Sciences journal homepage: www.elsevier.com/locate/ijrmms
Technical Note
Numerical determination of elastic compliance tensor of fractured rock masses by finite element modeling Jian-Ping Yang a, Wei-Zhong Chen a,b,n, Yong-hao Dai a, Hong-Dan Yu a a State Key Laboratory of Geomechanics and Geotechnical Engineering, Institute of Rock and Soil Mechanics, Chinese Academy of Sciences, 430071 Wuhan, PR China b Research Center of Geotechnical and Structural Engineering, Shandong University, 250061 Jinan, PR China
art ic l e i nf o Article history: Received 18 March 2013 Received in revised form 5 December 2013 Accepted 22 June 2014
1. Introduction Discontinuities including joints, fissures and interface separation are pervasive in rock masses at or near the Earth’s surface. In the analysis of engineering problems dealing with rock masses, such as design of arch dams, bridge piers, tunnels and underground caverns, the discontinuities are of great importance because of their strong influence on the mechanical and hydraulic properties. The treatment of discontinuities in numerical analysis can be divided in two ways. First, the discontinuities are described explicitly if amount of fractures is not large in the analyzing domain. This approach is represented mostly by distinct element method (DEM) and discontinuous deformation analysis (DDA) [1,2], which considers fractured rock masses as an assemblage of intact rock blocks separated by fracture systems. Second, in areas the rock masses containing extremely large number of fractures compare to engineering domain, it is impractical and impossible to model every fracture, and furthermore, the mechanical behavior of rock masses seems to be continuous or pseudo-continuous. So the equivalent continuum method is often adopted, which assume the global behavior of rock masses can be described by the principles of continuum mechanics. Based on the fact that many engineering facilities are constructed in heavily fractured rock masses, the equivalent continuum method is widely used. In the continuum approach, the equivalent elastic mechanical parameters of fractured rock masses are important input
n Corresponding author at: State Key Laboratory of Geomechanics and Geotechnical Engineering, Institute of Rock and Soil Mechanics, Chinese Academy of Sciences, 430071 Wuhan, PR China. Tel./fax: þ 86 27 87199212. E-mail address: [email protected] (W.-Z. Chen).
http://dx.doi.org/10.1016/j.ijrmms.2014.06.007 1365-1609/& 2014 Elsevier Ltd. All rights reserved.
parameters in any analysis of rock masses deformation behavior. Considerable efforts have been made in the past for determining the equivalent elastic mechanical properties. The determination methods can be divided into field tests [3], empirical relationships [4–8], analytical methods [9–15] and numerical methods [16–18]. Min and Jing summarized the characteristics and shortcomings of these methods, and a methodology determining elastic compliance tensor for fractured rock masses using DEM was established in that study [19]. The main advantage of the developed methodology is that it can consider complex fracture system geometry and various constitutive relations of fractures and rock matrix, and their interactions. However, the DEM modeling model in the study ignored the dead-end of fractures and their effects on equivalent elastic behavior of fractured rock masses. The ignorance would cause change of stress distribution and fracture network, which may induce great deviation in determination of compliance matrix. Within the framework established by Min and Jing [19], equivalent elastic mechanical parameters of fractured rock masses were studied by FEM (Finite Element Method) modeling in this paper. The advantage of the modeling is in developing a new mesh generation program for the assembly of intact rock and fractures, by which dead-end of fractures could be incorporated. Here, we briefly describe the methodology to determine elastic compliance matrix for fracture rock masse, and introduce the loading sets and corresponding numerical data processing method, then develop a new mesh generating technique which consider dead-end of fractures. Finally, we conduct a series of numerical simulations at different sizes and rotate angles to study scale effects and anisotropic characteristics of elastic constants of the studied fractured rock masses, and evaluate the applicability of continuum mechanics.
J.-P. Yang et al. / International Journal of Rock Mechanics & Mining Sciences 70 (2014) 474–482
475
2. Modeling approach using FEM code 2.1. Methodology for the determination of the linear elastic compliance tensor for fractured rock masses by FEM modeling The problem considered here is to estimate the equivalent compliance matrix for fractured rock masses whose long fractures are parallel in the direction of out-of-plane, z-direction. The matrix is assumed isotropic. Elastic properties of this kind of fractured rock masses can be investigated by two dimensional plane strain analysis. The linear elastic relation of anisotropic media can be expressed as:
εmn ¼ Smnpq σ pq 2
For plane strain problems, Eq. (1) can be reduced to 3 3 2σ 3 S12 S13 S16 xx 6σ 7 S22 S23 S26 7 yy 7 7 6 7 7 6 7 S32 S33 S36 5 6 4 σ zz 5 σ xy S62 S63 S66
2 εxx S11 6ε 7 6 S 6 yy 7 6 21 6 7 6 εzz 7 ¼ 6 4 5 4 S31 εxy S61
ð1Þ
ð2Þ
Methodologies for the determination of components of the compliance matrix in Eq. (2) using numerical methods have been studied [19]. For completeness, these methods are briefly summarized here. The elastic modulus and the Poisson’s effects in the z-direction are kept as intact rock since that there is no effect of the out-of-plane fractures. Considering the symmetry conditions, then S13 ¼S31 ¼S23 ¼S32 ¼ ν/E, and S33 ¼1/E. Since shear stress σxy does not affect the z-direction deformation, S36 and S63 can be assumed as zero. Therefore, three independent stress boundary conditions are sufficient to determine all components in the matrix shown in Eq. (2). Three loading sets (set 1 to set 3 in Fig. 1) were applied to determine the nine unknowns in Eq. (3), which are S11, S22, S12, S21, S26, S61, S62 and S66. Theoretically, S12 ¼S21, S16 ¼2S61 and S26 ¼2S62 because of the symmetry conditions, it is necessary to pointed out here that because the tensor shear strain εxy , instead of the engineering shear strain γ xy , is used in Eq. (2), then there exists a factor of 2 between S16, S26 and S61, S62. Components S61, S62 and S16, S26 represent the influence of normal stresses on shear strains and normal stresses on shear strains, respectively. Determination of S16 and S26 requires the calculation of averaged normal strain in x-and y-directions under pure shear stress (loading set 3 in Fig. 1), which is much more accurate and easier compared to the calculation of averaged shear strain under normal stress (loading set 1 and set 2 in Fig. 1). So, seven components-S11, S22, S12, S21, S16, S26 and S66 are determined by numerical experiments, and S61, S62 are deduced from S16, S26. Imposing εzz ¼ 0 and expressing σ zz in terms of σ xx and σ yy , Eq. (2) can be written as
εxx ¼ A11 σ xx þ A12 σ yy þ S16 τxy εyy ¼ A21 σ xx þ A22 σ yy þ S26 τxy εxy ¼ S61 σ xx þ S62 σ yy þ S66 τxy
ð3Þ
where Aij ¼ Sij
Si3 Sj3 ði; j ¼ 1; 2Þ S33
ð4Þ
Through imposing loading set 1, S11 and S21 can be obtained by substituting stress and the calculated strain into Eq. (3). As illustrated in Fig. 2(a), the normal strain εxx and εyy are calculated by dividing sample length L with deformation in corresponding directions, which are ΔL11 þ ΔL12 in the x-direction and ΔL21 þ ΔL22 in y-direction. ΔL11 is the average displacement of all nodes on left boundary. All the other deformations were calculated in the same way. Similarly, S22 and S12 can be obtained through imposing loading set 2. In theory, S21 equals S12 . However,
Fig. 1. Schematic diagram of four loading sets: (1) compression in x direction; (2) compression in y direction; (3) pure shear; (4) biaxial normal stress loading. (The first three loading sets are employed to calculate the compliance matrix of the studied model, while the fourth loading set is employed to calculate the shear modulus in 451 direction based on the current model).
∆L 12
∆L 11
∆L 22 u22 u11 A 11
y ∆L 21
x
u12 u21
L Fig. 2. Schematic diagram of (a) normal strain evaluation and (b) shear strain evaluation.
numerical simulation results would show a little difference between them. By imposing loading set 3, S16, S26 and S66 can be obtained. The shear strain εxy is calculated from the change in angle of adjacent perpendicular boundary as illustrated in Fig. 2(b). Angle change of each boundary is calculated from the area formed by deformed and original boundary. Numerical simulation results showed that angle change of every adjacent boundary differs a little from the others. Averaging method is adopted and εxy is evaluated by ðθ11 þ θ12 þ θ21 þ θ22 Þ=4. To obtain S16 and S26 , the average normal strains in x and y direction are calculated the way Fig. 2(a) illustrated. An additional biaxial loading set 4 (Fig. 1(d)) was applied to calculate equivalent shear modulus in the 451 direction compared to the current direction based on the current model. The biaxial normal stress creates pure shear stress in 451 direction whose value equals p. By calculating the average displacement of all boundary nodes, engineering shear strain γ in the 451 direction is easily obtained from Eq. (5).
π L ðΔL11 þ ΔL12 Þ γ ¼ 2arctan 2 L þ ðΔL21 þ ΔL22 Þ
ð5Þ
476
J.-P. Yang et al. / International Journal of Rock Mechanics & Mining Sciences 70 (2014) 474–482
Force boundary conditions were applied in numerical simulation rather than displacement boundary conditions which would induce additional forces caused by discontinuous deformation along fractures. Besides, the inertia relief technique in FEM program was used to avoid rigid body displacement, and the central point of the sample was chosen as reference fixed rotation point. 2.2. Procedure of directional deformation modulus analysis by FEM modeling While the equivalent continuum approach is employed for mechanical analysis in fractured rock masses, the equivalent elastic compliance tensor of the fractured rock masses is a basic input parameter, and it can be determined by the following steps: Step 1: investigate the distributions of fractures in site, eg. trace lengths, dip directions, fracture spacing, etc. Step 2: generate the fracture network models by Monte Carlo simulation method based on fracture distribution parameters. Step 3: generate the mesh of fractured rock masses based on the fracture network generated in step 2. Step 4: determine the elastic parameters of intact rock and fractures separately by laboratory experiments. Step 5: study the scale effect and anisotropic properties of equivalent elastic constants of fractured rock masses using FEM based on the mesh generated in step 3 and elastic parameters obtained in step4. First, choose a point as the center of study windows. Second, study the equivalent elastic constants at different scales increased step-by-step. Third, rotate the study windows incrementally for each scale and study the elastic constants in different directions. Finally, determine the equivalent elastic compliance matrix and corresponding REV if it exists.
2.3. Techniques for finite element mesh generation of fractured rock masses The generated mesh in step 3 of Section 2.2 is critical for FEM simulation results. In previous studies on elastic moduli and permeability dependence on stress of fractured rock masses, only hydraulic connected network are incorporated in the analysis model, and the dead-end of fractures which have no contribution to fluid transport are ignored [17,19]. The ignorance may alter stress distribution and fracture network as discussed before. In order to avoid these shortages, a new technique generating finite element mesh for fractured rock masses is developed, the detailed processes is described as follows: 1. Mesh the constructed fracture network model with triangular elements. Fracture lines inside the model must be seeded with nodes. 2. Suppose i and j are original nodes on fracture, and triangular elements related to them are E1 : ijk, E2 : ilj. First, copy a pair of nodes i1 , j1 compare to original nodes i, j. Positions of i1 , j1 are calculated according to the positions of i, j and the fracture width. Then, renumber the triangular elements as E1 : i1 j1 k, E2 : ilj, and generate the corresponding fracture element E3 : ijj1 i1 , as shown by Fig. 3(a). 3. On the position of two fracture intersection, another three nodes are needed in addition to the existed node. Positions of the three copied nodes are calculated according to the position of original node and the fracture widths of the two intersecting fractures. So the nodes of elements, which are connected to the intersection node, should be renumbered. In addition, four fractured elements are generated. Modifications of fracture
Fig. 3. Mesh technology of fractures and surrounding matrix.
intersection in Fig. 3(b) are: E1 : ijk-E1 : i2 j1 k; E2 : ikl-E2 : i2 kl1 ; E3 : ilm-E3 : i3 lm1 E4 : imn-E4 : imn; E5 : ino-E5 : ino; E6 : ioj-E6 : i1 o1 j E7 : li3 i2 l1 ; E8 : ioo1 i1 ; E9 : mii3 m1 ; E10 : i1 jj1 i2 4. The cross area of two intersecting fractures (area enclosed by ii1 i2 i3 in Fig. 3(b)) is left blank. Because the mechanical properties are different between normal direction and shear direction of fracture element, a fracture element generated on the cross area would prevent one of the intersecting fractures from sliding along it. However, along which fracture the rock masses would slide is dependent on the stress state. So the sliding direction can’t be preset and the “common region” of intersection fractures is left blank. 5. On the fracture end, one additional node is generated. Position of the new node is calculated by the position of the end node and the fracture width. Modification of the surrounding elements is shown in Fig. 3(c).
3. Verification of the proposed FEM modeling technique In order to verify the proposed modeling technique, a model with two orthogonal fractures sets is used for the comparison between the deformation modulus produced by FEM experiments and the closedform expression proposed by Amadai and Goodman [10]. The basic parameters used are listed as follows. For isotropic intact rock, the elastic modulus and Poisson’s ratio is 50 GPa and 0.25. For both fracture sets, the normal stiffness and shear stiffness are 50 GPa/m and 10 GPa/m, and the fracture spacing is 1.5 m, coupling effects between shear deformation and normal stress, and normal deformation and shear stress are not considered. Considering the geometric symmetry, the computational models are rotated in intervals of 151 from 01 to 451. Three loading sets described in Section 2.1 are applied to obtain the compliance matrix
J.-P. Yang et al. / International Journal of Rock Mechanics & Mining Sciences 70 (2014) 474–482
for 301 model. As shown in Fig. 4, both the numerical results calculated from rotated models in 01, 151, 301 and 451, and the red solid line produced by the compliance matrix derived from the 301 model, are agreed very well with the analytical results. To study the effects of S61 , S62 , S16 and S26 on global elastic properties, a comparison curve (the red dotted line in Fig. 4) is plotted where the four components are assumed to be zero. It can be seen that large
477
deviation against analytical results exists except in the direction of 301. Displacement vector maps (Fig. 5) also show obvious angel change of boundary lines under normal stresses, and elongation in x direction under shear stresses. Quantitative analysis shows that, under compressive stress 1 MPa in x direction, εxx ¼ 5.2 10 5, εxy ¼ 9.5 10 6, under shear stress 1 MPa, εxx ¼ εyy ¼ 2.3 10 5, εxy ¼4.7 10 5. The magnitudes of normal strain and shear strain are almost in the same order. The results prove strong coupling effects between normal stresses and shear strains, and between shear stresses and normal strains. From discussion above, it is reasonable to conclude that the proposed FEM modeling technique and compliance matrix computing method are valid and can be used for more complex random fracture models.
4. Results of scale dependency and tensor characteristics of elastic parameters for fractured rock masses Fig. 4. Verification of numerical modeling results against analytical solution. (The black solid line corresponds to the analytical results. The red solid line is produced from compliance matrix derived from 301 model, while the red dot line is produced from the same matrix where components of S16, S26, S61 and S62 are ignored. Points of symbols corresponds to numerical results of rotated FEM models). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
4.1. Data of fractures and mechanical parameters for the analysis From site investigation on central part of auxiliary tunnel B of JinPing hydropower station II, two sets of fractures are identified and the distribution of mean fracture spacing (distance of adjacent
Fig. 5. Displacement vector maps illustrating coupling behavior between normal stresses and shear strains, and between shear stresses and normal strains (unit: m). (a) Vector of total displacement under normal stress, (b) vector of displacement in x direction under normal stress, (c) vector of total displacement under shear stress and (d) vector of displacement in x direction under shear stress.
478
J.-P. Yang et al. / International Journal of Rock Mechanics & Mining Sciences 70 (2014) 474–482
Table 1 Input parameters for random network of fractures. Fractures set
1 2
Fracture spacing
Trace length
Dip direction
Distribution pattern
Mean/ m
Standard deviation/ m
Distribution pattern
Mean/ m
Standard deviation/ m
Distribution pattern
Mean/ m
Standard deviation/ m
Normal distribution Normal distribution
2.5
0.1
4
1
10
0.1
3
0.7
Normal distribution Normal distribution
150
2.0
Normal distribution Normal distribution
50
7
fracture centroids), fracture strike and mean trace length are listed in Table 1. The intact rock is considered isotropic linear elastic medium. Elastic modulus and Poisson’s ratio of intact rock are 50 GPa and 0.25 by laboratory test results. The constitutive relation for a fracture is expressed as: 9 2 9 8 38 K nn K ns K nt > > = = < Δτn > < Δδn > 7 Δτs ¼ 6 ð6Þ 4 K sn K ss K st 5 Δδs > > > > : Δτ ; K tn K ts K tt : Δδt ; t where the subscript ‘n’ represents the normal direction (i.e. perpendicular) to the fracture plane; ‘s’ and ‘t’represent two orthogonal directions on the fracture plane. Thus, Δτn is the normal stress, Δτs and Δτt are the shear stresses; Δδn is the normal displacement of the fracture, Δδs and Δδt are the shear displacements. K ij (i; j ¼ n; s; t) represents components of the fracture stiffness matrix: K nn represents the traction behavior of a fracture; K ss and K tt represent the shear behavior of a fracture; K ns and K nt represent coupling effect of shear displacement caused by normal stress; K sn and K tn represent coupling effect of normal displacement caused by shear stress; K st and K ts represent coupling effect between two shear stresses. Experiments indicate that the normal behavior and shear behavior of a fracture are often coupled [20], and the fracture stiffness is influenced by many factors [21,22] (particularly at large increment of load). The fracture stiffness may be treated as constants in a small increment of load. For a large increment of load, constant fracture stiffness can only be used as approximate estimation. This study neglects the coupling effects between shear behavior and normal behavior, and treats normal stiffness and shear stiffness of fractures as constants, which are 50 GPa/m and 10 GPa/m determined from small increment load tests. The coupling components and closure-variable stiffness could be incorporated into the FEM modeling easily. Fig. 6 shows the numerical simulation cases in this study. Model sizes range from 2 m to 16 m with 2 m interval, and angles rotated from 01 to 901 with 301 interval. Elastic properties of 01–901, 301–1201 and 601–1501 models are obtained from FEM modeling directly. In addition, the equivalent shear modulus of the fractured rock masses in 151, 451 and 751 directions are evaluated by models of 601, 01 and 301 accordingly.
90° 120°
60° 30°
150°
2m 4m 6m 8m 10 m
0°
16 m
Fig. 6. Study cases of deformation modulus of fractured rock masses.
gradient can be observed for 10 m model due to the collective contributions of large number of fractures (Fig. 7(a)). And this is the physical basis in describing macroscopic mechanical behavior of fractured rock masses by continuum mechanics. Fig. 7(b) shows the displacement contour of fracture network without considering fracture dead-end. A large number of fracture dead-ends are deleted compared with fracture network in Fig. 7(a), and the displacement fields are also totally different. Quantitative evaluation shows that, for model considering dead-end, εxx ¼ 4.15 10 5 and εxx ¼23.4 GPa, while for models without dead-end, εxx ¼ 2.77 10 5 and εxx ¼34.5 GPa. The model with dead-end is much softer because of the existence of more fractures. Fig. 8(a) and (b) show that stresses concentrate around fracture tips. Larger shear stresses along fractures in 901 model imply larger shear strains compared with that of 601 model. Consequently, greater displacement will be induced for 901 model. That is the reason the equivalent deformation modulus is the lowest in 901 and highest in 601 (Fig. 9(a)).
4.2. Results and discussion 4.2.1. Deformation characteristics and stress distributions of FEM modeling results Fig. 7 shows the contours of directional displacements of 01–901 FEM model under 1 MPa compressive stresses (loading set 2). The zero displacement point lies in the center of the model since the balanced stresses were applied on the symmetrical two sides. Displacement jump often exists across fractures in local scale. However, a general trend of roughly uniform deformation
4.2.2. Scale effects and tensor characteristics of fractured rock masses It has been suggested that there are two prerequisites for the appropriateness of continuum approach in mechanical analysis of fractured rock masses [17,23]. 1. A REV must exist, then equivalent parameters can be derived from heterogeneous fractured rock masses. The values of
J.-P. Yang et al. / International Journal of Rock Mechanics & Mining Sciences 70 (2014) 474–482
479
Fig. 7. Contrast of displacement contours under compressive loading with and without fracture dead-ends (unit: m). (a) 10 m model including fracture dead-ends and (b) 10 m model without fracture dead-ends.
Fig. 8. Shear stress contours of 10 m FEM models (unit: Pa). (a) 1 MPa compressive stress in 601 direction (horizontal) and (b) 1 MPa compressive stress in 901 direction (vertical).
equivalent parameters change insignificantly as the sizes of rock masses exceed REV. This problem is referred as size effect. 2. The derived equivalent parameters can be expressed in tensor form for the usage of constitutive equations for continuum method. The problem is related to the anisotropic characteristic of fractured rock masses. Problem 1 is studied by learning the variation trends between equivalent elastic constants and model sizes. If all elastic constants are tending towards stability and the variations are within predefined acceptable ranges, the REV is assumed existing. The minimum size beyond which elastic constants change within predefined ranges is defined as REV. Variation of elastic constants (deformation modulus, shear modulus and Poisson’s ratio) is measured by VEC i ¼
jEC i EC i j EC i
ð7Þ
where EC i is elastic constants evaluated by FEM model with side length i, EC i is averaged elastic constants of FEM models having the same rotation angle and whose side length is large or equal than i. For problem 2, it is expected that a compliance matrix derived from one FEM model can predict elastic constants of other
directions in an acceptable error. Predicting error is evaluated by comparison of compliance matrix components between results from FEM modeling and predicted by compliance matrix of other directions. The predicting error is defined as FEM CM 1 jSij Sij j CE ¼ ∑ N i;j SCM ij
ð8Þ
where SFEM and SCM are components of compliance matrix in the ij ij same reference coordinate system, SFEM are the FEM modeling ij results, SCM are the results transformed from compliance matrix of ij other directions, and N denotes the number of components to be compared. 4.2.3. Results of the scale effects Fig. 9 presents the variation of the equivalent elastic constants (deformation modulus, shear modulus, Poisson’s ration) with increasing side length of FEM model. Fig. 9(a) exhibits that the equivalent deformation modulus increase with increasing side length and become stable after a certain size. Variation of modulus in all six directions is less than 1 GPa when the side length reaches 8 m. For modulus in all directions, the values of 2 m 2 m scale are the smallest in all sizes, because the mean lengths of the two fracture
480
J.-P. Yang et al. / International Journal of Rock Mechanics & Mining Sciences 70 (2014) 474–482
Results of equivalent shear modulus shown by Fig. 9(b) are consistent with the results of deformation modulus. Shear modulus are the lowest in 451 and 601, which are almost parallel with one fracture set. The large shear modulus in 01 and 151 are also related with directional properties of fracture sets. Compared Fig. 9(c) to (a), an opposite relationship between Poisson’s ratio and deformation modulus can be seen at both model scale and directional aspects. Whether it is a special or common relationship for fractured rock masses needs more study. Calculation indicates that the value of νxy =Ex almost equals νyx =Ey . This result proves the symmetry characteristics of compliance matrix from numerical aspect. Quantitative evaluation by Eq. (7) shows that variations of deformation modulus, shear modulus and Poisson’s ratio are below 5%, 10% and 5% in 8 m size. As the size increased to 10 m, the variations of the three studied elastic constants are all decreased below 5%.
Fig. 9. Deformation modulus, shear modulus and Poisson’s ratio with increasing side lengths of FEM models (a) deformation modulus (b) shear modulus (c) Poisson’s ratio.
sets are larger than the computational model size, and the computational model is interpenetrated by several single fractures. Results from Fig. 9(a) also show larger equivalent elastic modulus in direction 601 and 1501 than in direction 01 and 901 as the loading direction of 601 and 1501 is nearly parallel to one fracture set and perpendicular to the other one. It is similar to the verification one with two orthogonal fracture sets. So the fractured rock mass is relatively rigid in the parallel direction and soft in angle bisector direction.
4.2.4. Results of tensor characteristics Figs. 10 and 11 present comparisons of equivalent deformation modulus and shear modulus between FEM modeling and compliance matrix. In the two figures, the black, blue and magenta dash dot line represent directional deformation/shear modulus expected by compliance matrix derived from 01, 301, 601 FEM model respectively. The separate symbols are FEM modeling results. Plotting lines show that the all compliance matrixes express elastic anisotropy characteristics properly even at very small scale (say 2 m). This may be due to the four fractures inside 2 m model, as shown in Fig. 6, have almost the same directions with the global fractures. However, the values of both deformation and shear modulus are not stable in small sizes. As model size increases, the deviation between the three lines decrease, and all the three lines predict FEM modeling results well. Fig. 12 presents mean prediction error by Eq. (8) between results of FEM modeling and matrixes. Because six independent unknowns, S11 , S22 , S12 , S16 , S26 and S66 exist in compliance matrix, and S11 , S22 , S66 are thought as dominating ones [17], tensor predicting errors are evaluated by incorporating six and three components separately for comparison. Fig. 12(a) shows that the predicting error is under 10% in 6 m and 5% in 10 m when three components are compared. However, the predicting error increases greatly while six components are contained. The value is under 20% in 10 m and 12% in 12 m as shown by Fig. 12(b). Results exhibit that the compliance matrix derived from 601 model can predicts FEM modeling results in other directions very well, while the reverse is not. For example, the error between matrix of 601 model and 301 FEM results is 2.1% while the reverse is 12%. Meanwhile, the error between matrix of 601 model and 01 FEM results is 4.3%. This implies that the strategy of averaging compliance matrix [17] may be not the best choice, as perhaps one model could give better predicting results, as shown in the studied case. In this sense, the compliance matrix of 601 model is thought best representing the elastic behavior of the studied fracture rock masses. After translating to the reference coordinate system, the compliance matrix is: 0 1 42:92 17:84 5:00 5:04 B 17:82 46:72 5:00 8:68 C B C 12 S¼B ðpa1 Þ ð9Þ C 10 @ 5:00 5:00 20:00 0:00 A 2:54
4:32
0:00
43:37
5. Conclusions Scale effects and anisotropic characteristics of fractured rock masses are studied by a proposed FEM modeling technique
J.-P. Yang et al. / International Journal of Rock Mechanics & Mining Sciences 70 (2014) 474–482
481
Fig. 10. Comparison of deformation modulus between rotated FEM models and compliance matrixes. (a) 2 m, (b) 6 m, (c) 10 m and (d) 14 m. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Fig. 11. Comparison of shear modulus between rotated FEM models and compliance matrixes. (a) 2 m, (b) 6 m, (c) 10 m and (d) 14 m. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
which consideres dead-end fractures. The main conclusions are summarized below. A mesh generating technique is developed for fractured rock masses by which the dead-end of fractures can be incorporated into FEM modeling. In cases where the fracture spacing is close to or larger than the fracture length, plenty of fractures and fracture sections would become dead-end and ignoring of them may leads to wrong results. In the studied case, the FEM modeling result of 10 m model without dead-end almost overestimate the deformation modulus as much as 50% compared to the model incorporating fracture dead-end. The influence coefficients of normal stresses on shear strains, and shear stresses on normal strains, corresponding to S61 , S62 and S16 , S26 in compliance matrix, plays an important role in describing
elastic mechanical behavior of fractured rock masses properly. They can influence both the deformation modulus and anisotropic characteristic greatly. This implies that carefully measurements and analysis are required in laboratory and field tests of elastic constants of fractured rock masses. Both the variation of elastic constants on scale and tensor prediction error should be evaluated in determining the proper compliance matrix for fractured rock masses. For the studied case, variations of deformation modulus, shear modulus and Poisson’s ratio are below 5% as model size reaches to 10 m, and tensor prediction error is also below 5% when S11 , S22 , S66 are contained in prediction. However, after all the six independent componentsS11 , S22 , S12 , S16 , S26 and S66 in compliance matrix are incorporated to evaluate discrepancies between matrix and FEM modeling
482
J.-P. Yang et al. / International Journal of Rock Mechanics & Mining Sciences 70 (2014) 474–482
Foundation of China (Grant Nos. 51225902, 51274189 and 51109207). We would also like to express our sincere gratitude to the editor and the anonymous reviewers for their valuable comments, which have greatly improved this paper. References
Fig. 12. Mean prediction error of compliance matrix components between FEM results and predicted results from matrixes in other directions (a) mean prediction error of 3 components, S11 , S22 and S66 (b) mean prediction error of 6 components, S11 , S22 , S12 , S16 , S26 and S66 .
results, the predicting errors are 20% in 10 m and 12% after reaching 12 m. Considering the importance of S12 , S16 and S26 , the REV is determined as 12 m, which is about three times as much as mean trace length in the studied case (mean fracture spacing is smaller than mean trace length). Acknowledgements The authors gratefully acknowledge the support of the Chinese Fundamental Research (973) Program through the Grant No. 2013CB036006 and the support of the National Natural Science
[1] Hart RD. An introduction to distinct element modeling for rock engineering. In: Hudson JA, editor. Comprehensive rock engineering. Oxford: Pergamon Press; 1993. p. 245–61. [2] SHI GH. Discontinuous deformation analysis: a new numerical model for the statics and dynamics of block system. Ph.D. thesis. Univ of California, Berkeley; 1988. [3] Bieniawski ZT. Determining rock mass deformability-experience from case histories. Int J Rock Mech Min Sci 1978;15:237–47. [4] Palmstrom A, Singh R. The deformation modulus of rock masses: comparisons between in situ tests and indirect estimates. Tunnelling Underground Space Technol 2001;16:115–31. [5] Barton N. Some new Q-value correlations to assist in site characterisation and tunnel design. Int J Rock Mech Min Sci 2002;39:185–216. [6] Sonmez H, Gokceoglu C, Ulusay R. Indirect determination of the modulus of deformation of rock masses based on the GSI system. Int J Rock Mech Min Sci 2004;41:849–57. [7] Zhang L, Einstein HH. Using RQD to estimate the deformation modulus of rock masses. Int J Rock Mech Min Sci 2004;41:337–41. [8] Hoek E, Diederichs MS. Empirical estimation of rock mass modulus. Int J Rock Mech Min Sci 2006;43:203–15. [9] Singh B. Continuum characterization of jointed rock masses: Part I–The constitutive equations. Int J Rock Mech Min Sci 1973;10:311–35. [10] Amadei B, Goodman RE. A 3-D constitutive relation for fractured rock masses. In: Proceedings of the international symposium on the mechanical behavior of structured media. Ottawa; 18–21 May 1981. p. 249–68. [11] Gerrard CM. Elastic models of rock masses having one, two and three sets of joints. Int J Rock Mech Min Sci 1982;19:15–23. [12] Oda M. Fabric tensor for discontinuous geological materials. Soils Found 1982;22:96–108. [13] Oda M, Suzuki K, et al. Elastic compliance for rock-like materials with random cracks. Soils Found 1984;24:27–40. [14] Hu KX, Huang Y. Estimation of the elastic properties of fractured rock masses. Int J Rock Mech Min Sci 1993;30:381–94. [15] Huang TH, Chang CS, Yang ZY. Elastic moduli for fractured rock mass. Rock Mech Rock Eng 1995;28:135–44. [16] Kulatilake PHSW, Wang S, et al. Effect of finite size joints on the deformability of jointed rock in three dimensions. Int J Rock Mech Min Sci 1993;30:479–501. [17] Min KB, Jing L. Numerical determination of the equivalent elastic compliance tensor for fractured rock masses using the distinct element method. Int J Rock Mech Min Sci 2003;40:795–816. [18] Chen SH, He J, Shahrour I. Estimation of elastic compliance matrix for fractured rock masses by composite element method. Int J Rock Mech Min Sci 2012;49:156–64. [19] Min KB, Rutqvist J, et al. Stress-dependent permeability of fractured rock masses: a numerical study. Int J Rock Mech Min Sci 2004;41:1191–210. [20] Jing L, Nordlund E, Stephansson O. An experimental study on the anisotropy and stress-dependency of the strength and deformability of rock joints. Int J Rock Mech Min Sci 1992;29:535–42. [21] Barton N. The shear-strength of rock and rock joints 1976;13:255–79Int J Rock Mech Min Sci 1976;13:255–79. [22] Bandis SC, Lumsden AC, et al. Fundamentals of rock joint deformation. Int J Rock Mech Min Sci 1983;20:249–68. [23] Long JC, Remer JS, Wilson CR. Porous-media equivalents for networks of discontinuous fractures. Water Resour Res 1982;18:645–58.