Journal of Applied Mathematics

Volume 2013 (2013), Article ID 420536, 19 pages

http://dx.doi.org/10.1155/2013/420536

Research Article

A Model of Anisotropic Property of Seepage and Stress for Jointed Rock Mass

Pei-tao Wang,1,2 Tian-hong Yang,1,2 Tao Xu,1,2 Qing-lei Yu,1,2 and Hong-lei Liu1,2

¹Key Laboratory of Ministry of Education on Safe Mining of Deep Metal Mines, Northeastern University, Shenyang 110819, China

²School of Resources & Civil Engineering, Northeastern University, Shenyang 110819, China

Received 30 March 2013; Revised 18 June 2013; Accepted 19 June 2013

Academic Editor: Pengcheng Fu

Copyright © 2013 Pei-tao Wang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

Joints often have important effects on seepage and elastic properties of jointed rock mass and therefore on the rock slope stability. In the present paper, a model for discrete jointed network is established using contact-free measurement technique and geometrical statistic method. A coupled mathematical model for characterizing anisotropic permeability tensor and stress tensor was presented and finally introduced to a finite element model. A case study of roadway stability at the Heishan Metal Mine in Hebei Province, China, was performed to investigate the influence of joints orientation on the anisotropic properties of seepage and elasticity of the surrounding rock mass around roadways in underground mining. In this work, the influence of the principal direction of the mechanical properties of the rock mass on associated stress field, seepage field, and damage zone of the surrounding rock mass was numerically studied. The numerical simulations indicate that flow velocity, water pressure, and stress field are greatly dependent on the principal direction of joint planes. It is found that the principal direction of joints is the most important factor controlling the failure mode of the surrounding rock mass around roadways.

1. Introduction

Underground mining has been considered a high-risk activity worldwide. Violent roof failure or rock burst induced by mining has always been a serious threat to the safety and efficiency of mines in China. Accurate and detailed characterization for rock masses can control stable excavation spans, support requirements, cavability, and subsidence characteristics, and thus influence the design of mining layouts and safety of mines. Rock mass is a geologic body composing of the discontinuities which have a critical influence on deformational behavior of blocky rock systems [1]. The mechanical behavior of this material depends principally on the state of intact rock whose mechanical properties could be determined by laboratory tests and existing discontinuities containing bedding planes, faults, joints, and other structural features. The distributions and strength of these discontinuities are both the key influencing factors for characterizing the discontinuous and anisotropic materials. Amadei [2] pointed out the importance of anisotropy of jointed rock mass and discussed the interaction existing between rock anisotropy and rock stress. Because of computational complexity and the difficulty of determining the necessary elastic constants, it is usual for only the simplest form of anisotropy, transverse isotropy, to be used in design analysis [3]. The key work is to study the principal direction of elasticity or permeability and then assume the jointed rock mass as transversely isotropic geomaterial.

Extensive efforts have been made to investigate the mechanical response of transversely anisotropic rock material. Zhang and Sanderson [4] used the fractal dimension to describe the connectivity and compactness of fracture network and found that the deformability and the overall permeability of fractured rock masses increase greatly with increasing fracture density. By using the artificial transversely isotropic rock blocks, the mechanical properties with different dip angles were obtained by Tien and Tsao [5]. Brosch et al. [6] evaluated the fabric-dependent anisotropy of a particular gneiss by studying the strength values and elastic parameters along different directions. Exadaktylos and Kaklis [7] presented the explicit representations of stresses and strains at any point of the anisotropic circular disc compressed diametrically. The stress concentration in uniaxial compression for a plane strain model was investigated, and a formulation of failure criteria for elastic-damage and elastic-plastic anisotropic geomaterials is formed by Exadaktylos [8]. A methodology to determine the equivalent elastic properties of fractured rock masses was established by explicit representations of stochastic fracture systems [9], and the conditions for the application of the equivalent continuum approach for representing mechanical behavior of the fractured rock masses were investigated. Such anisotropy has an effect on the interpretation of stress distributions. Amadei’s results indicate that an anisotropy ratio (damaged elastic modulus to intact elastic modulus) of between 1.14 and 1.33 will have a definite effect on the interpreted in situ state of stress [10]. Hakala et al. [11] concluded that the anisotropy ratio of Finland rock specimens is about 1.4 and suggested to be taken into account in the interpretation of stress measurement results.

Moreover, in rock engineering like mining and tunnel engineering, interactions between in situ stress and seepage pressure of groundwater have an important role. Groundwater under pressure in the joints defining rock blocks reduces the normal effective stress between the rock surfaces and therefore reduces the potential shear resistance which can be mobilized by friction. Since rock behavior may be determined by its geohydrological environment, it may be essential in some cases to maintain close control of groundwater conditions in the mine area. Therefore, accurate description of joints is an important topic for estimating and evaluating the deformability and seepage properties of rock masses.

Many researchers have done a lot of studies on the seepage and anisotropic mechanical behaviors of the jointed rock mass. Using homogeneous samples like granite or basalt, Witherspoon et al. [12] investigated and defined the permeability by fracture aperture in a closed fracture. Based on geometrical statistics, Oda [13] has studied and determined the crack tensor of moderately jointed granite by treating statistically the crack orientation data via a stereographic projection. Using Oda’s method, Sun and Zhao [14] determined the anisotropy in permeability using the fracture orientation and the in situ stress information from the field survey. For fractured rock system, attempts have been made by Jing et al. [15–18] and Bao et al. [19] to investigate the permeability. Using UDEC code, Jing et al. studied the permeability of discrete fracture network, such as the existence of REV in [15], relations between fracture length and aperture in [16] or stress effect on permeability in [17], or solute transport in [18]. Bao et al. [19] have discussed the mesh effect on effective permeability for a fractured system using the upscaled permeability field. It would help the workers to spend less computational effort and memory requirement to investigate effective permeability. The key for determining deformation modulus and hydraulic parameters is to study representative elementary volume (REV) and scale effects of fractured rock mass [9, 20–22]. However, the difficulty of testing jointed rock specimens, at scales sufficient to represent the equivalent continuum, indicated that it is necessary to postulate and verify methods of synthesizing rock mass properties from those of the constituent elements like intact rock and fractures.

Although great progress has been made, it is difficult to study the anisotropic deformation of large-scale rock mass, ground water permeability change, and interaction between stress and seepage due to the geological complexity of discontinuities. Qiao et al. [23] have pointed out that the rock mass which is not highly fractured and has only few sets of joint system usually behaves anisotropically. However, the mechanical properties directionality of highly jointed rock mass is usually ignored. In particular, since the limitation of mesh generation, the numerical method can hardly deal with the highly jointed specimens. The discontinuity like fault could be treated as specific boundary. An effective solution should be found to characterize the influence of joints on the rock mass.

In this paper, based on equivalent continuum theory and theoretical analysis, a mathematical model for anisotropic property of seepage and elasticity of jointed rock mass is described. A permeability analysis code is developed to evaluate the anisotropic permeability for DFN model based on VC++ 6.0 in this paper. The DFN model could be analyzed to evaluate the REV size and anisotropic property of permeability, which would provide important evidence for the finite element model. The outline is as follows. First, 3D images involving detailed geometrical properties of rock mass, such as trace lengths, outcrop areas, joint orientations, and joint spacing, were captured using ShapMetriX3D system. According to the statistical parameters for each set of discontinuities, fracture network is generated using Monte Carlo method, and jointed rock samples in different sections could be captured. Next, permeability tensors and elasticity tensors of rock mass of the sample are calculated by discrete medium seepage method and geometrical damage theory, respectively. Finally, using the finite element method (FEM), an anisotropic mechanical model for rock mass is built. An engineering practice of roadway stability at the Heishan Metal Mine, in Hebei Iron and Steel Group Mining Company, is described. Then the stress and seepage fields surrounding the roadway were numerically simulated. On the basis of the modeled results, the influences of joint planes on stress, seepage, and damage zone were analyzed. It is expected throughout this study to gain an insight into the influences of discontinuities on the mechanical behavior of rock mass and offer some scientific evidence for the design of mining layouts or support requirements.

2. Generation of a Fracture System Model by ShapeMetriX3D

Traditional methods for rock mass structural parameters in mining engineering applications include scanline surveying [24] and drilling core method [25]. The process generally requires physical contact with rock mass exposure and therefore is hazardous. Additionally, taking manual measurements is time consuming and prone to errors due to sampling difficulties or instrument errors. ShapeMetriX3D is a tool for the geological and geotechnical data collection and assessment for rock masses [26]. The equipment is used for the metric acquisition of rock and terrain surfaces and for the contact-free measurement of geological/geotechnical parameters by metric 3D images. The images are captured using Nikon D80 camera with 22.3 megapixel. As shown in Figure 1, two images are acquired to reproduce the 3D rock mass face. Details information related to the measurements is reported in [27] by Gaich et al. The measurements could be taken at any required number and extent, even in regions that are not accessible. During measuring process, two digital images taken by a calibrated camera serve for a 3D reconstruction of the rock face geometry which is represented on the computer by a photorealistic spatial characterization as shown in Figure 1. From it, measurements are taken by marking visible rock mass features, such as spatial orientations of joint surfaces and traces, as well as areas, lengths, or positions. Finally, the probability statistical models of discontinuities are established. It is generally applied in the typical situations such as long rock faces at small height and rock slopes with complex geometries [28]. Almost any rock face can be reconstructed at its optimum resolution by using this equipment and its matching software. By using this system we can increase working safety, reduce mapping time, and improve data quality.

Figure 1: Imaging principle of the 3D surface measurement.

In order to accurately represent rock discontinuities, ShapeMetriX3D (3GSM) in this paper is used for the metric acquisition of rock mass exposure and for the contact-free measurement of geological parameters by metric 3D images. Stereoscopic photogrammetry deals with the measurement of three-dimensional information from two images showing the same object or surface but taken from two different angles, just as shown in Figure 1.

From the determined orientation between the two images and a pair of corresponding image points and , imaging rays (colored in red) are reconstructed whose intersection leads to a 3D surface point . By automatic identification of corresponding points within the image pair, the result of the acquisition is a metric 3D image that covers the geometry of the rock exposure. Once the image of a rock wall is ready, geometric measurements can be taken as shown in Figure 2. There are a total of three groups of discontinuities in this bench face.

Figure 2: Geological mapping and geometric measurements of stereoscopic restructuring model.

Figure 3 is facilitated to show the 3-dimensional distribution of the trace along the section of Profile I-I from the stereoscopic model shown in Figure 2. The height of this rock face is approximately 5.0 m with a distance of 2.0 m perpendicular to the trace. The measured orientation in a hemispherical plot could also be captured as in Figure 4(a). Figure 4(b) shows the results of the distribution of joints.

Figure 3: Length and sketch of Profile I-I by stereoscopic model.

(a)

(b)

Figure 4: Stereographic plot and the distribution of total three sets of joints.

The generated 3D rock face is about 5 m × 5 m with a high resolution enough to distinguish the fractures. According to the survey results, the geologic data like joint density, dip angle, trace length, and spacing can be acquired. Baghbanan and Jing [16] generated the DFN models whose orientations of fracture sets followed the Fisher distribution. In this paper, four different probability statistical models are used to generate the DFN models. The dip angle, dip direction, trace length, and spacing all follow one particular probability statistical model. Table 1 shows the basic information about the fracture system parameters. Type I of the probability statistical model in Table 1 stands for negative exponential distribution, Type II for normal distribution, Type III for logarithmic normal distribution, and Type IV for uniform distribution. Based on the probability models, the particular fracture network could be generated using Monte Carlo method [29, 30].

Table 1: Characteristics of 3 sets of discontinuities for the slope exposure.

Set
Density (m⁻¹)
Fracture characteristics

Dip direction (°)
Dip angle (°)
Trace length (m)
Spacing (m)

Type
Mean value
Standard deviation
Type
Mean value
Standard deviation
Type
Mean value
Standard deviation
Type
Mean value
Standard deviation

1#
1.2
II
90.73
13.88
III
71.24
15.05
II
0.86
0.23
IV
0.83
0.75

2#
1.3
II
289.96
22.89
IV
59.85
8.4
II
0.89
0.28
IV
0.76
0.73

3#
3.3
III
189.05
14.55
III
55.29
13.35
II
0.85
0.25
IV
0.3
0.31

3. Constitutive Relation of Anisotropic Rock Mass

Due to the existence of joints and cracks, the mechanical properties (Young’s modulus, Poisson’s ratio, strength, etc.) of rock masses are generally heterogeneous and anisotropic. Three preconditions should be confirmed in this section: (i) anisotropy of rock mass is mainly caused by IV or V-class structure or rock masses containing a large number of discontinuities; (ii) rock masses according to (i) could be treated as homogeneous and anisotropic elastic material; (iii) seepage tensor and damage tensor of rock mass with multiset of joints could be captured by the scale of representative elementary volume (REV).

3.1. Stress Analysis

Elasticity represents the most common constitutive behavior of engineering materials, including many rocks, and it forms a useful basis for the description of more complex behavior. The most general statement of linear elastic constitutive behavior is a generalized form of Hooke’ Law, in which any strain component is a linear function of all the stress components; that is,

where is the flexibility matrix, is the strain, and is the stress.

Many underground excavation design analyses involving openings where the length to cross-section dimension ratio is high are facilitated considerably by the relative simplicity of the excavation geometry. In this section, the roadway is uniform cross section along the length and could be properly analyzed by assuming that the stress distribution is the same in all planes perpendicular to the long axis of the excavation (Figure 5). Thus, this problem could be analyzed in terms of plane geometry.

Figure 5: A diagrammatic sketch for underground excavation.

The state of stress at any point can be defined in terms of the plane components of stress (, and ) and the components (, and ). In this research, the direction is assumed to be a principal axis and the antiplane shear stress components would vanish. The plane geometric problem could then be analyzed in terms of the plane components of stress since the component is frequently neglected. Equation (1), in this case, may be recast in the form

where is the flexibility matrix of the material under plane strain conditions. The inverse matrix (or ) could be expressed in the form

where is expressed as follows:

The moduli and and Poisson’s ratios and could be provided by uniaxial strength compression or tension in 1 (or 2) and 3 directions.

The mechanical tests including laboratory and in situ tests for rock masses in large scale can hardly capture the elastic properties directly. Hoek-Brown criterion [31–33] could calculate the mechanical properties of weak rocks masses by introducing the Geological Strength Index (GSI). Nevertheless, the anisotropic properties cannot be captured using this criterion, and thus the method for analyzing anisotropy of jointed rock mass needs a further study.

In this paper, the original joint damage in rock mass is considered as macro damage field. In elastic damage mechanics, the elastic modulus of the jointed material may degrade and the Young’s modulus of the damaged element is defined as follows [34]:

where represents the damage variable and and are the elastic moduli of the damaged and intact rock samples, respectively. In this equation, all parameters are scalar.

With the geometric information of the fracture sample, the damage tensor [35] could be defined as

where is the number of joints, is the minimum spacing between joints, is the volume of rock mass, is the normal vector of the kth joint, and is the trace length of the kth joint (for 2 dimensions).

According to the principal of energy equivalence [36], the flexibility matrix for the jointed rock sample can be obtained as

where is the flexibility matrix for intact rock; and are the principal damage values in and directions, respectively. For plane strain geometric problem, the constitutive relation, where the coordinates and principal damage have the same direction, could be expressed as follows:

where is the Young’s modulus of intact rock and is Poisson’s ratio for intact rock.

The parameters could be easily captured by laboratory tests. Equation (7) gives the principal damage values and . Based on geometrical damage mechanics, all elements in this matrix could be obtained by the method mentioned previously, and the anisotropic constitutive relation of jointed rock sample could finally be confirmed.

3.2. Seepage Analysis

The seepage parameters of rock mass are quantized form of permeability and also are the basis to solve seepage field of equivalent continuous medium. Based on the attribute of fractured rock mass and by taking engineering design into account, fractured rock mass is often considered to be anisotropic continuous medium. In the fracture network shown in Figure 6, a total number of cross points or water heads and line elements are contained. Parameters can be obtained with the model such as water head, the related line elements, equivalent mechanical fissure width, seepage coefficient, and so on. The corresponding coordinates of each point could be acquired. For a fluid flow analysis based on the law of mass conservation, the fluid equations on a certain water head take the form [37]

where is the quantity of flow from line element to water head is the total number of line elements intersect at , and is the fluid source term. In the joint network, each line element would be assigned a length and fissure width to investigate the permeability.

Figure 6: Fracture network schematic diagram in seepage area. and are the constant head boundaries with the hydraulic head of and respectively; and are the impervious boundaries with flow of zero.

According to the hydraulic theory [38], for a single joint seepage, the flow quantity of line element can be expressed as

where is the hydraulic gradient, is the coefficient of flow viscosity, is the density of water, and is the fissure width of joint. On the basis of (9) and (10), the governing equations can be represented as (11) for seepage in fracture network

Based on the discrete fracture network method [39], an equivalent continuum model for seepage has been established [28, 40, 41]. The hydraulic conductivity can be acquired based on Darcy’s Law on the basis of water quantity in the network (see (12)).

Boundary conditions in the model shown in Figure 6 are as follows:

(a) and as the constant head boundary;(b) and as the impervious boundary with flow of 0.

Then the hydraulic conductivity can be defined as

where is the total quantity of water in the model region (m²·s⁻¹), is the water pressure difference between inflow and outflow boundaries (m), is the equivalent hydraulic conductivity coefficient in the direction (m·s⁻¹), and and are the side lengths of the region (m).

Based on Biot equations [42], the steady flow model is given by

where is the hydraulic conductivity and is the hydraulic pressure. For plane problems, the dominating equation of seepage flow is as follows in (14). The direction of joint planes is considered to be the principal direction of hydraulic conductivity

3.3. Coupling Mechanism of Seepage and Stress

3.3.1. Seepage Inducted by Stress

The coupling action between seepage and stress makes the failure mechanism of rock complex. The investigations on this problem have pervasive theoretical meaning and practical value. The principal directions associated with the symmetric crack tensor are coaxial with those of the permeability tensor. The first invariant of the crack tensor is proportional to the mean permeability, while the deviatoric part is related to the anisotropic permeability [13]. Generally, the change of stress which is perpendicular to the joints plane is the main factor leading to the increase or decrease of the ground water permeability. In the numerical model, seepage is coupled to stress describing the permeability change induced by the change of the stress field. The coupling function can be described as follows as given by Louis [43]:

where is the current groundwater hydraulic conductivity, is the initial hydraulic conductivity, is the stress perpendicular to the joints plane, and is the coupling parameter (stress sensitive factor to be measured by experiment) that reflects the influence of stress. The larger is, the greater the range of stress induced the permeability [44].

3.3.2. Stress Inducted by Seepage

On the basis of generalized Terzaghi’s effective stress principle [45], the stress equilibrium equation could be expressed as follows for the water-bearing jointed specimen:

where is the total stress tensor, is the elastic tensor of the solid phase, is the strain tensor, is a positive constant which is equal to 1 when individual grains are much more incompressible than the grain skeleton, is the hydraulic pressure, and is the Kronecker delta function.

3.4. Coordinate Transformation

Generally, planes of joints are inclined at an angle to the major principal stress direction as shown in Figure 7. In establishing these equations, the , and 1, 3 axes are taken to have the same axis, and the angel is measured from the to the 1 axis.

Figure 7: The relationship between planes of joints and the principal stress direction.

Considering the hydraulic pressure needs no coordinate transformation, and (16) will be expressed as

where is the reverse matrix for stress coordinates transformation and is the strain coordinates transformation matrix, and they can be expressed as follows:

Similarly, the coordinate transformation form can be expressed as

where is the hydraulic conductivity coefficient along the distribution direction of joint planes, is perpendicular to the distribution direction, and is measured from the coordinate system to the optimal direction of permeability.

4. A Case Study

4.1. Description of the Area under Study

The model described above is applied to evaluate the seepage field and stress field of jointed rock roadway in the Heishan Metal Mine (Figure 8). The mine is located in Chengde city, Hebei province, in northern China. Heishan Metal Mine has transferred from open pit mining to underground mining since 2009. The deformation and stability of the roadways for mining the hanging-wall ore become the key technique issue. The elevation of the crest of the slope is 920 m. The roadway in the research is in 674 m level of the high northern slope. The rock mass is mainly composed of anorthosite and norite in the northern slope area. The anisotropic properties of seepage and stress field will be discussed in detail.

Figure 8: The rock mass exposure of northern slope in Heishan open pit mine.

4.2. Capture of Joint Network

Depending on the 3D contact-free measuring system, discontinuities on the exposure could be easily captured. The collection process for the geology information of northern slope in Heishan Metal Mine has been discussed previously. A 3D fracture network of rock mass could be obtained using the Monte Carlo method by analyzing the statistical parameters (Figure 9(a)). Finally, the fracture network is generated and expanded to the roadway, and the fracture network perpendicular to the center line of roadway, could be easily captured (Figure 9(b)).

(a)

(b)

Figure 9: Joint network based on the statistical data using Monte Carlo method. (a) 3D joint network in rock mass; (b) profile with 10 m edge length perpendicular to the roadway.

4.3. Permeability Investigation

According to the symmetry of geometry, the method for solving the hydraulic conductivity coefficient in different directions of fracture network is presented to save computational time [41].

Rotated every 15° clockwise and exerted certain water pressure on the boundary shown in Figure 10, the hydraulic conductivity coefficient of verticaldirection , 75° direction , 60° direction , 45° direction , 30° direction , and 15° direction could be acquired, respectively. Then the permeability tensors by fracture parameters can be determined. By enlarging the geometric size of the fracture network, the permeability scale effect can be investigated finally.

Figure 10: Generated fracture network by Monte Carlo method (side length of internal squares is 1 m, 4 m, 7 m, 8 m, and 9 m, resp.).

Based on the algorithm in Section 3.2, the size effect of rock mass in different sizes of statistics window is studied. The region of fracture network is enlarged from 1 m to 9 m with different step lengths until the equivalent parameters in different directions achieve a stable value. With a fixed center of this region, the equivalent permeability is calculated rotating every 15° clockwise as shown in Figure 11.

Figure 11: Fracture network in the study area with edge length of 7 m.

The variation of hydraulic conductivity coefficients with the increase of direction angles and sample size is depicted in Figure 12. The main direction angle is approximated to 15°. It can be seen that the anisotropy of seepage property is apparent with the distribution of joints. Seven hydraulic conductivity coefficients of different sizes (1 m, 3 m, 4 m, 5 m, 6 m, 8 m, and 9 m) are chosen to be compared with those of the sample whose size is 7 m. The coefficient deviation between one particular size and 7 m would be found and presented in Figure 13. Values of permeability decrease with the increase of sample size and tend to stabilize when the sample size comes to 7 m.

Figure 12: Permeability tensor of the rock mass (unit: m/s).

Figure 13: Deviation of the permeability values to that of  7 m sample under different directions of samples in different sizes.

The principal values of permeability in maximum and minimum in a rock sample with joint plane angle being 15° are listed in Table 2. The principal permeability values decrease from 4.23 (×10−6 m/s) to 2.66 (×10−6 m/s) as the sample size increases from 1 m to 9 m along the main direction. According to the results discussed previously, the REV is 7 m × 7 m for this fracture sample, and the related hydraulic conductivity is 2.77 (×10−6 m/s) for maximum and 0.87 (×10−6 m/s) for minimum. The ratio for the maximum to minimum hydraulic conductivity value is 3.18.

Table 2: Principal values of fracture hydraulic conductivity at different sizes.

Sample size (m)
1
4
7
8
9

Principal value of hydraulic conductivity (×10⁻⁶ m/s)
Maximum
4.23
3.28
2.77
2.69
2.66

Minimum
0.80
0.97
0.87
0.88
0.85

Ratio
5.29
3.38
3.18
3.06
3.13

4.4. Damage Tensor

Similar to the algorithm in Section 3.1, the size effect of sample damage in different sizes of statistics window is studied. The region of fracture network is enlarged from 3 m to 12 m with different step lengths until the damage values in different directions stabilize. With a fixed center of this region, the principal damage values are calculated rotating every 15° clockwise. Figure 14 shows the sample damage in different sizes and directions. Similar to the deviation analysis of permeability, the damage values also fluctuatingly tend to stability according to the deviation of the damage tensor in different sizes and directions shown in Figure 15. In this study, the size of representative elementary volume is also 7 m in length, and the initial damage tensor of REV of jointed rock mass could also be obtained. The principal damage values and for fracture sample are 0.17 for minimum and 0.50 for maximum, and the principal direction angle of damage tensor is approximately 105° which is perpendicular to the direction of principal permeability. The related rock parameters, including undamaged Young’s modulus , Poisson’s ratio , and uniaxial tensile strength, et al, are listed in Table 3. It should be noted that the elasticity and strength of rock could be determined in the laboratory tests. Finally, the constitutive relation can be found according to (8). Equations (17) and (19) can be introduced into the FEM simulations.

Table 3: Parameters used in the model to validate the model in simulating the anisotropic properties of stress and seepage.

Material parameters
Values

Undamaged Young’s modulus (GPa)
65

Undamaged Poisson’s ratio
0.23

Density (kg/m⁻³)
2700

Principal direction angle of joint plane (°)
15

Coupling parameter (Pa⁻¹)
0.5

Figure 14: Damage tensor (×10⁻¹).

Figure 15: Deviation of the damage tensor to 7 m sample under different directions of samples in different sizes.

4.5. Geometry and Boundary Conditions

A numerical model of Heishan Metal Mine is established in order to simulate the mechanism of stress and seepage in jointed rock mass, taking into account the anisotropic property, as shown in Figure 16. Table 3 shows the basic mechanical properties. Parameters listed in Table 4 are the hydraulic conductivity coefficients , damage tensor and Young’s modulus , and shear modulus of the fracture sample along joint plane direction used in the finite element code. The model contains two roadways with a 3.5 m × 3 m three-centered arch section within a 50 m × 50 m domain. The bottom boundary of the domain is fixed in all directions, and the left and the right boundaries are fixed in the horizontal direction. In this regard, a pressure () of 5.97 MPa is applied on the top boundary of the model to represent the 221 m deep overburden strata. Under the steady-state groundwater flow condition, a hydrostatic pressure, , of 2.21 MPa is applied upon top boundary. No-flow conditions are imposed on the three boundaries of the rectangular domain. The initial water pressure on the roadways boundaries is 0 MPa. All the governing equations described previously are implemented into COMSOL Multiphysics, a powerful PDE-based multiphysics modeling environment. The model is assumed to be in a state of plane strain (with no change in elastic strain in the vertical direction) and static mechanical equilibrium.

Table 4: The hydraulic conductivity coefficients , damage tensor and Young’s modulus , and shear modulus of the fracture sample along joint plane direction.

Subscripts
( m/s)

(GPa)
(GPa)

11
2.77
0.50
18.83
—

22
0.87
0.17
51.91
—

12
—
—
9.34
10.97

(a)

(b)

Figure 16: Plane strain roadway model: (a) boundary conditions; (b) the finite element mesh.

4.6. Results and Discussion

4.6.1. Stress Distribution

Adverse performance of the rock mass in the postexcavation stress field may be caused by either failure of the anisotropic medium or slip on the weakness planes [3]. The elastic stress distribution around the roadways directly influences the deformation of rock and thus determines the design process.

Figure 17 shows the contour of first principal stress coupled with the seepage process. The orientation and magnitude of maximum principal stress controlled the distribution of the stress concentration in the heterogeneous media. The simulation result shows that the principal stress concentration zones appear mainly in rock surrounding the roadways. There exist maximum stress concentration areas in the arch foot and floor. Measures should be taken to control the deformation and assure the construction safety.

Figure 17: Contours of the first principal stress for ° case.

To characterize the response of the stress to the hydraulic mechanics, a comparison of two scenarios is also presented as shown in Figure 18. The first principal stress in the stress-seepage coupled model along the horizontal section A-A′ where is compared with a decoupled models. The result shows that the first principal stress increases when the seepage process is considered. Figure 19 shows a plot of normal stress in joint plane direction and hydraulic conductivity along the horizontal section A-A′ where as is shown in Figure 16. The coupled and decoupled model could be analyzed using Comsol Multiphysics code. In the coupled case, the governing equations for solid and fluid phase are solved in weakly coupled sense. For the sake of convenient contrast, the normal stress distribution as well as hydraulic conductivity is plotted when no seepage-coupled process is considered. The result shows that the normal stress will be underestimated when no coupling of seepage process is included. Moreover, the compressive stress leads to the decrease of permeability, and the hydraulic conductivity increases when seepage process is included.

Figure 18: Contrast of first principal stress between coupled and decoupled model.

Figure 19: Observed changes of normal stress in joint plane direction and hydraulic conductivity between coupled and decoupled models.

4.6.2. Seepage Distribution

Flow velocity with the angle of joint plane being 15° is shown in Figure 20. All the process is not considered time dependent. Therefore, the velocity change is induced by the permeability variation caused by the compressive stress upon the joint plane. In this case study, water flows along the principal direction of joint plane where the permeability is largest in the model. An asymmetric seepage field is observed, and maximum of flow velocity is distributed on right top of the roadways roof. Moreover, the seepage pressure is also asymmetrically distributed as shown in Figure 21.

Figure 20: Darcy’s velocity magnitude (m/s).

Figure 21: The seepage pressure distribution (Pa).

4.6.3. Damage Zone

For underground mining, two principal engineering properties of the joint planes should be considered. They are low tensile strength in the direction perpendicular to the joint plane and the relatively low shear strength of the surfaces. The anisotropic strength parameters as tensile and compressive strength parameters for jointed rock mass are discussed by Chen et al. [46], Claesson and Bohloli [47], Nasseri et al. [48], Gonzaga et al. [49], and Cho et al. [50]. As discussed previously, the fluid pressure and velocity are sensitive to the joint plane angles . In this section, the Hoffman anisotropic strength criterion is used to assess the damage zone in this numerical model as shown in

The properties and represent the tensile strength along joint plane direction and perpendicular to joint plane direction, respectively. and represent the compressive strength along joint plane direction and perpendicular to joint plane direction, respectively. is the shear strength of the material along the joint direction. represents the normal stress along the principle direction of elasticity and represents the normal stress perpendicular to the principle direction of elasticity. represents the shear stress.

The shear strength of the joints can be described by the simple Coulomb law

where is cohesive strength and is the effective angle of friction of the joint surfaces.

It should be noted that the mechanical parameters be acquired from laboratory or in situ tests. However, it is difficult to directly employ the strength parameters for jointed rock mass due to the inaccessibility of the tests for huge rock mass. Based on the in situ and laboratory tests of Heishan Metal Mine, the tensile, compressive, and shear strength parameters are listed in Table 5.

Table 5: Mechanical properties of rock mass in the direction of joint plane.

Tensile strength (MPa)

16.0

4.0

Compressive strength (MPa)

120

96

Shear strength (MPa)

1.5

50

The direction of joint plane is 15°, and thus the tensile strength in the direction perpendicular to joint plane is relatively low, compared with that in other directions. Figure 22 shows the damage zone in this case study and the failure area mainly distributes in the direction perpendicular to the weakness plane, which gives an illustration of the roadway failure mode in tabular orebodies. Moreover, the failure of covered rock mass and the rock pillar in and between the two roadways does not influence each other in this case study, and thus the choice of roadway’s interval is proper from the aspect of the mechanics analysis.

Figure 22: Damage zone of the simulating model (the direction of damage zone is approximately perpendicular to that of the joint planes whose angle is 15°).

4.7. Further Discussion

These simulations were performed to develop an understanding of the mechanics of joints and influence on stress and seepage fields and to gauge the ability of the proposed transversely anisotropic model to capture the response of jointed rock mass. For this purpose, a total of six scenarios with the joint plane angles ranging from 0° to 150° with an interval of 30° are simulated in order to examine the effect of joint plane directions; see Figure 7 for the definition of . And the anisotropic properties of seepage field and damage zones are examined in these simulations.

4.7.1. Seepage Field

Figure 23 presents the fluid pressure distribution and flow field arrows with different joint plane directions. It can be seen that the maximum pressure, which is located on the top boundary, equals the initial fluid pressure (2.21 MPa). The fluid pressure distribution on top of the roadways differs with the increase of directions.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 23: Contour of the fluid pressure and flow velocity vectors in different joint plane directions.

When the joint plane is parallel or perpendicular to the floors of roadways, the fluid pressure distributions and flow arrows are all axially symmetric, which agrees with expectations. The fluid pressure in the 30° or 60° case is distributed unsymmetrically as shown in Figures 23(b) and 23(c).

The flow velocity in direction along the horizontal section A-A′ where  m as shown in Figure 16(a) curves are plotted in Figure 24 for °, 30°, 60°, and 90°, respectively. In all cases the absolute value of velocity increases approximately exponentially above the roadways. As shown in Figure 24, the flow velocity in scenario where ° is in the lowest level. The reason is that the joint plane in this scenario is horizontally distributed and the compressive stress leads to the decrease of permeability. When angle of joint plane increases to 30°, flow velocities above both roadways increase and Roadway I has a more higher velocity than Roadway II. The scenario where ° has similar performance. However, when the joint plane is vertically distributed, flow velocity above Roadway I is lower than that where °. The flow velocity in the case where joints are vertically or horizontally distributed is symmetric.

Figure 24: The flow velocity in direction along the horizontal section A-A′ for different numerical models, at θ = 0°, 30°, 60°, and 90°.

4.7.2. Damage Zone

The effect of joint plane angle on the damage zone is illustrated by using the Hoffman anisotropic strength criterion as shown in (23). The damage zones in different angles of joint planes are shown in Figure 25 and could be an index to visualize the potential failure mode of the roadway.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 25: The damaged zone under different principal elastic directions (red arrows represent the flow velocity vector).

It is clear from Figures 25(a) to 25(f) that an increase of the joint plane angle has significantly influenced the shape and size of damage zone. When the joints are horizontally distributed (°), the tensile strength in the direction perpendicular to joint planes is much lower and thus the damage zone mainly concentrates within roof and bottom of roadways. This failure mainly manifests as roof falling and floor heave. Similarly, when the angle increases to 30°, the damage zones mainly concentrate in the rock mass of left roof and right floor within the roadways which is also in the direction perpendicular to joint planes. Compared with the result in Figure 25(b), the direction of damage zone with ° shown in Figure 25(c) rotated significantly and the area also increases. When the joints are vertically distributed (°), the principal direction of damage zone is nearly horizontal. Lateral rock mass surrounding the roadways stabilizes with the increase of joint plane angle in this study, and the pillar between two roadways is also stable. The scenarios where ° or ° have the similar response to the joint plane direction with that of 60° or 30°. Qualitatively, the simulation results are in good agreement with the results by Jia et al. [51] and Huang et al. [52]. The model in this study could to some extent be effectively used to analyze the anisotropic property for jointed rock mass.

5. Conclusion

The main purpose is to introduce the anisotropic model of seepage and stress and apply the model to the jointed rock mass. The coupled finite element analysis was performed, and the effects of joint planes direction on stress field, flow velocity, and Darcy’s velocity were verified. A more reliable model for the stability analysis of rock engineering and risk evaluation of water inrush was provided. Based on the results of a series of numerical simulations under different scenarios, the following conclusions are drawn.

(1)A linkage of digital information of fractures and mechanical analysis is realized. Relations between digital images involving detailed geometrical properties of rock mass and the quantitative determination of hydraulic parameters as well as elastic properties were realized. The results show that the scale for both damage tensor and permeability tensor of the rock mass’s REV in northern slope of Heishan Metal Mine is 7 m 7 m. (2)We examine the model for seepage-stress coupled analysis on anisotropic properties. The numerical simulations in this study have indicated that the existence of joint planes greatly affects the seepage properties and stress field. In anisotropic rock, water flows mainly along the joint planes, and water pressure is asymmetrically distributed where the angle of joint plane is 15° in the northern slope of Heishan Metal Mine.(3)The influence of fractures cannot be neglected in stability analysis of rock mass. The numerical results visualized the damage zones in different directions of joint planes. The direction of damage zones was found perpendicular to that of joint planes which agrees well with the field observations and theoretical analysis.

The work reported in this paper is an initial effort on the influence of joint orientation on the anisotropic property of seepage and stress in jointed rock masses. A model for more complex jointed rock mass remains to be quantified.

Nomenclature

: | Damage variable (dimensionless)
: | Damaged Young’s modulus (Pa)
: | Undamaged Young’s modulus (Pa)
: | Shear modulus (Pa)
: | Hydraulic conductivity coefficient (m/s)
: | Hydraulic conductivity coefficient along joint plane direction (m/s)
: | Fluid pressure (Pa)
: | Flexibility coefficient (Pa⁻¹)
: | The number of joints
: | Minimum spacing between joints (m)
: | Volume of rock mass (m³)
: | The hydraulic gradient;
: | The normal vector of the kth joint (dimensionless)
: | The trace length of the kth joint (m)
: | The fluid source term ()
: | The reverse matrix for stress coordinates transformation
: | The strain coordinates transformation matrix
: | The tensile strength in joint direction (Pa)
: | The tensile strength perpendicular to direction (Pa)
: | The compressive strength in joint direction (Pa)
: | The compressive strength perpendicular to direction (Pa)
: | The shear strength of the material (Pa).

Greek Symbols

: | Positive constant
: | Coupling coefficient (Pa⁻¹)
: | Stress tensor (Pa)
: | The first, second, and third principal stresses (Pa)
: | Total strain tensor (dimensionless)
: | Angle of joint plane (°)
: | Damaged Poisson’s ratio (dimensionless)
: | Undamaged Poisson’s ratio (dimensionless)
: | Density of water (kg/m³)
: | Coefficient of flow viscosity (Pa·s)
: | Effective angle of friction of the joint surfaces.

Conflict of Interests

The authors declare that there is no conflict of interests.

Acknowledgments

This study was performed at Center for Rock Instability & Seismicity Research (CRISR), Northeastern University, Shenyang, China. The work presented in this paper was financially supported by the National Basic Research Program of China (No. 2013CB227900). The Natural Science Foundation of China (Grant No. 51174045, 50904013, 51034001, 50934006, 41172265), the basic scientific research funds of Ministry of Education of China (No. N090101001, N120601002) and the project supported by the Specialized Research Fund for the Doctoral Program of Higher Education of China (Grant No. 20120042120053) also give financial support to this research.

References

W. G. Pariseau, S. Puri, and S. C. Schmelter, “A new model for effects of impersistent joint sets on rock slope stability,” International Journal of Rock Mechanics and Mining Sciences, vol. 45, no. 2, pp. 122–131, 2008.

B. Amadei, “Importance of anisotropy when estimating and measuring in situ stresses in rock,” International Journal of Rock Mechanics and Mining Sciences and Geomechanics, vol. 33, no. 3, pp. 293–325, 1996.

B. H. G. Brady and E. T. Brown, Rock Mechanics For Underground Mining, Springer, Amsterdam, The Netherlands, 2004.

X. Zhang and D. J. Sanderson, “Numerical study of critical behaviour of deformation and permeability of fractured rock masses,” Marine and Petroleum Geology, vol. 15, no. 6, pp. 535–548, 1998.

Y. M. Tien and P. F. Tsao, “Preparation and mechanical properties of artificial transversely isotropic rock,” International Journal of Rock Mechanics and Mining Sciences, vol. 37, no. 6, pp. 1001–1012, 2000.

F. J. Brosch, K. Schachner, M. Blümel, A. Fasching, and H. Fritz, “Preliminary investigation results on fabrics and related physical properties of an anisotropic gneiss,” Journal of Structural Geology, vol. 22, no. 11-12, pp. 1773–1787, 2000.

G. E. Exadaktylos and K. N. Kaklis, “Applications of an explicit solution for the transversely isotropic circular disc compressed diametrically,” International Journal of Rock Mechanics and Mining Sciences, vol. 38, no. 2, pp. 227–243, 2001.

G. E. Exadaktylos, “On the constraints and relations of elastic constants of transversely isotropic geomaterials,” International Journal of Rock Mechanics and Mining Sciences, vol. 38, no. 7, pp. 941–956, 2001.

K.-B. Min and L. Jing, “Numerical determination of the equivalent elastic compliance tensor for fractured rock masses using the distinct element method,” International Journal of Rock Mechanics and Mining Sciences, vol. 40, no. 6, pp. 795–816, 2003.

B. Amadei, Rock Anisotropy and the Theory of Stress Measurements, Springer, 1983.

M. Hakala, H. Kuula, and J. A. Hudson, “Estimating the transversely isotropic elastic intact rock properties for in situ stress measurement data reduction: a case study of the Olkiluoto mica gneiss, Finland,” International Journal of Rock Mechanics and Mining Sciences, vol. 44, no. 1, pp. 14–46, 2007.

P. A. Witherspoon, J. S. Y. Wang, K. Iwai, and J. E. Gale, “Validity of cubic law for fluid flow in a deformable rock fracture,” Water Resources Research, vol. 16, no. 6, pp. 1016–1024, 1980.

M. Oda, “Permeability tensor for discontinuous rock masses,” Geotechnique, vol. 35, no. 4, pp. 483–495, 1985.

J. Sun and Z. Zhao, “Effects of anisotropic permeability of fractured rock masses on underground oil storage caverns,” Tunnelling and Underground Space Technology, vol. 25, no. 5, pp. 629–637, 2010.

K.-B. Min, L. Jing, and O. Stephansson, “Determining the equivalent permeability tensor for fractured rock masses using a stochastic REV approach: method and application to the field data from Sellafield, UK,” Hydrogeology Journal, vol. 12, no. 5, pp. 497–510, 2004.

A. Baghbanan and L. Jing, “Hydraulic properties of fractured rock masses with correlated fracture length and aperture,” International Journal of Rock Mechanics and Mining Sciences, vol. 44, no. 5, pp. 704–719, 2007.

A. Baghbanan and L. Jing, “Stress effects on permeability in a fractured rock mass with correlated fracture length and aperture,” International Journal of Rock Mechanics and Mining Sciences, vol. 45, no. 8, pp. 1320–1334, 2008.

Z. Zhao, L. Jing, I. Neretnieks, and L. Moreno, “Numerical modeling of stress effects on solute transport in fractured rocks,” Computers and Geotechnics, vol. 38, no. 2, pp. 113–126, 2011.

K. Bao, A. Salama, and S. Sun, “Upscaling of permeability field of fractured rock system: numerical examples,” Journal of Applied Mathematics, Article ID 546203, 20 pages, 2012.

W. Z. Chen, J. P. Yang, J. L. Yang, X. B. Qiu, and C. C. Cao, “Hydromechanical coupled model of jointed rock mass and its application to pressure tunnels,” Chinese Journal of Rock Mechanics and Engineering, vol. 27, no. 9, pp. 1569–1571, 2006.

W. Chen, J. Yang, X. Zou, and C. Zhou, “Research on macromechanical parameters of fractured rock masses,” Chinese Journal of Rock Mechanics and Engineering, vol. 27, no. 8, pp. 1569–1575, 2008.

J. P. Yang, Study of macro mechanical parameters of fractured rock mass [Ph.D. thesis], Wuhan institute of the Chinese academy of sciences geotechnical engineering, Wuhan, China, 2009.

C. Qiao, Z. Zhang, and X. Li, “Anisotropic strength characteristics of jointed rock mass,” in Proceedings of the ISRM-Sponsored International Symposium on Rock Mechanics, Rock Characterisation Modelling and Engineering Design Methods, p. 140, 2009.

A. Robertson, “The interpretation of geological factors for use in slope theory,” in Proceedings of the Symposium on Theoretical Background to Planning of Open Pit Mines, pp. 55–71, Johannesburg, South Africa, 1970.

C.-Y. Wang and K. T. Law, “Review of borehole camera technology,” Chinese Journal of Rock Mechanics and Engineering, vol. 24, no. 19, pp. 3440–3448, 2005.

Austrian Startup Company, ShapeMetriX3D Model Merger User Manual, Earth Products China, Shenyang, China, 2008.

A. Gaich, W. Schubert, and M. Poetsch, “Three-dimensional rock mass documentation in conventional tunnelling using Joint-MetriX3D,” in Proceedings of the 31st ITA-AITES World Tunnel Congress, E. Yücel and T. Solak, Eds., pp. 59–64, A.A. Balkema Publishers, Istanbul, Turkey, 2005, Underground Space Use: Analysis of the Past and Lessons for the Future.

T. Yang, Q. Yu, S. Chen, H. Liu, T. Yu, and L. Yang, “Rock mass structure digital recognition and hydro-mechanical parameters characterization of sandstone in Fangezhuang coal mine,” Chinese Journal of Rock Mechanics and Engineering, vol. 28, no. 12, pp. 2482–2489, 2009.

P. T. Wang, T. H. Yang, Q. L. Yu, H. L. Liu, and P. H. Zhang, “Characterization on jointed rock masses based on PFC2D,” Frontiers of Structural and Civil Engineering, vol. 7, no. 1, pp. 32–38, 2013.

Y. Y. Jiao, X. L. Zhang, and T. C. Li, DDARF Method For Simulating the Whole Process of Rock Failure, Science Publishing House, Beijing. China, 2010.

E. Hoek and E. T. Brown, “Practical estimates of rock mass strength,” International Journal of Rock Mechanics and Mining Sciences, vol. 34, no. 8, pp. 1165–1186, 1997.

P. Marinos and E. Hoek, “Estimating the geotechnical properties of heterogeneous rock masses such as flysch,” Bulletin of Engineering Geology and the Environment, vol. 60, no. 2, pp. 85–92, 2001.

E. Hoek, C. Carranza-Torres, and B. Corkum, “Hoek-Brown failure criterion-2002 edition,” in Proceedings of the 5th North American Rock Mechanics Symposium and 17th Tunneling Association of Canada Conference (NARMS-TAC '02), pp. 267–271, 2002.

J. Lemaitre and R. Desmorat, Engineering Damage Mechanics, Springer, Berlin, Germany, 2005.

T. Kawamoto, Y. Ichikawa, and T. Kyoya, “Deformation and fracturing behaviour of discontinuous rock mass and damage mechanics theory,” International Journal for Numerical and Analytical Methods in Geomechanics, vol. 12, no. 1, pp. 1–30, 1988.

F. Sidoroff, “Description of anisotropic damage application to elasticity,” in Proceedings of the IUTAM Colloquium on Physical Nonlinearities in structural analysis, pp. 237–244, 1981.

Y. Wu and Z. Zhang, An Introduction to Rock Mass Hydraulics, Southwest Jiaotong University Press, Chengdu, China, 2005.

D. T. Snow, “Anisotropie Permeability of Fractured Media,” Water Resources Research, vol. 5, no. 6, pp. 1273–1289, 1969.

K. Hestir and J. C. S. Long, “Analytical expressions for the permeability of random two-dimensional Poisson fracture networks based on regular lattice percolation and equivalent media theories,” Journal of Geophysical Research, vol. 95, no. 13, pp. 21565–21581, 1990.

T. Yang and Y. Xiao, “Structural plane survey and analysis of permeability characteristics of open pit mine slope rock mass,” Site Investigation Science and Technology, vol. 16, no. 3, pp. 27–30, 1998.

Y. X. Xiao, C. F. Lee, and S. J. Wang, “Assessment of an equivalent porous medium for coupled stress and fluid flow in fractured rock,” International Journal of Rock Mechanics and Mining Sciences, vol. 36, no. 7, pp. 871–881, 1999.

M. A. Biot, “General theory of three-dimensional consolidation,” Journal of Applied Physics, vol. 12, no. 2, pp. 155–164, 1941.

C. Louis, “Rock hydraulics,” in Rock Mechanics, L. Muller, Ed., pp. 287–299, Springer, New York, NY, USA, 1974.

T. H. Yang, T. Xu, H. Y. Liu, C. A. Tang, B. M. Shi, and Q. X. Yu, “Stress-damage-flow coupling model and its application to pressure relief coal bed methane in deep coal seam,” International Journal of Coal Geology, vol. 86, no. 4, pp. 357–366, 2011.

K. Terzaghi, R. B. Peck, and G. Mesri, Soil Mechanics in Engineering Practice, Wiley-Interscience, New York, NY, USA, 1996.

C.-S. Chen, E. Pan, and B. Amadei, “Determination of deformability and tensile strength of anisotropic rock using Brazilian tests,” International Journal of Rock Mechanics and Mining Sciences, vol. 35, no. 1, pp. 43–61, 1998.

J. Claesson and B. Bohloli, “Brazilian test: stress field and tensile strength of anisotropic rocks using an analytical solution,” International Journal of Rock Mechanics and Mining Sciences, vol. 39, no. 8, pp. 991–1004, 2002.

M. H. B. Nasseri, K. S. Rao, and T. Ramamurthy, “Anisotropic strength and deformation behavior of Himalayan schists,” International Journal of Rock Mechanics and Mining Sciences, vol. 40, no. 1, pp. 3–23, 2003.

G. G. Gonzaga, M. H. Leite, and R. Corthésy, “Determination of anisotropic deformability parameters from a single standard rock specimen,” International Journal of Rock Mechanics and Mining Sciences, vol. 45, no. 8, pp. 1420–1438, 2008.

J.-W. Cho, H. Kim, S. Jeon, and K.-B. Min, “Deformation and strength anisotropy of Asan gneiss, Boryeong shale, and Yeoncheon schist,” International Journal of Rock Mechanics and Mining Sciences, vol. 50, pp. 158–169, 2012.

P. Jia, C.-A. Tang, T.-H. Yang, and S.-H. Wang, “Numerical stability analysis of surrounding rock mass layered by structural planes with different obliquities,” Dongbei Daxue Xuebao/Journal of Northeastern University, vol. 27, no. 11, pp. 1275–1278, 2006.

S. Huang, J. Xu, X. Ding, and A. Wu, “Study of layered rock mass composite model based on characteristics of structural plane and its application,” Chinese Journal of Rock Mechanics and Engineering, vol. 29, no. 4, pp. 743–756, 2010.