Development of a damage rheological model and its application in the analysis of mechanical properties of jointed rock masses

Fossil fuel resources (eg, coal, oil, and natural gas) exist in underground reservoirs, where discontinuities are widespread and influence the stability of rock structures. The mechanical behavior of jointed rock masses significantly affects the long‐term stability of rock engineering. A major challenge in this area is to link the time‐dependent deformation with damage influence induced by distribution of joints in rock masses. In this paper, a damage mechanical theory is adopted which deals with several groups of joints distributed in rock masses. Based on the geometrical distribution of joints, a damage model considering the influence of normal vector and area density of joints is used to describe the discontinuities. A damage rheological model for jointed rock masses is developed and programmed in the finite difference software FLAC3D, and an example of jointed rock masses is used to verify the validity of the model. The damage rheological model could predict the visco‐elastic strain which corresponds to primary creep property, visco‐plastic strain which corresponds to steady state creep property, and damage strain which reflects the damage effect of joints to rock masses. Finally, we use this model to predict the displacements and damage zones in the rock masses surrounding an underground opening.


| INTRODUCTION
The development and utilization of deep energy resources are often affected by discontinuities in rock masses. Discontinuities (joints, cracks, and weak intercalary strata) are widespread in rock masses, for example: coal mining, metal mining, hot rock reservoir for geothermal energy utilization, and petroleum reservoir. These joints and cracks significantly influence the stability of rock structures. How to evaluate the influence of these discontinuities on the properties of rock masses, especially the long-term behavior, is a big challenge. Figure 1 shows some adverse engineering conditions of jointed rock mass. Figure 1A shows slow time-dependent deformation of jointed rock masses in a deep gold mine in South African, 1 which showed time-dependent squeezing in jointed rock masses can occur, even though these tunnels are developed in hard crystalline rock. Figure 1B shows a pillar that is weakened by two pre-existing discontinuities that contributed to failure of this pillar at a relatively low pillar stress, which have resulted in safety concerns. 2

YANG et Al.
A comprehensive understanding of the mechanical properties of jointed rock mass is of fundamental importance in almost all rock engineering design and construction, including slopes, foundations, or underground excavations in rock. 3,4 Many studies have been conducted on the mechanical behaviors of specimens containing multiple joints through experimental or numerical simulation methods. [5][6][7][8][9][10][11][12][13][14][15] However, these studies focused on the propagation and coalescence of cracks in rock specimens. The stability of engineering projects was not analyzed, because it is very hard to apply these results to predict the mechanical behavior of in situ rock masses.
Field evidence has shown that damage of a jointed rock mass is a time-dependent behavior accompanied by the joint opening, closing, and extension. [16][17][18][19] In most rock masses, multiple sets of joints have been encountered. Depending on their relative scale compared with the structure, the joints can be analyzed in different ways. The rock masses containing relatively small joints can be idealized to be continuum bodies. The joint with relatively large dimensions can be represented by a layer of solid elements in the numerical modeling. However, effects of joints with intermediate size are hard to be estimated in a simple way, leading to the complexity in dealing with the mechanical behavior of such joints. To this end, damage mechanics provides an effective approach to tackle the difficulty. Kachanov 20 first introduced the damage concept to predict the time of creep rupture for metal under tensile stress. The theory used a scalar variable (ratio of voids in a cross-section) as the damage variable which is only applicable to one-dimensional problems. Hayhurst, 21 Dragon and Mróz , 22 and Murakami and Ohno 23 applied this theory to multidimensional problems, which is valid for softening behavior of rocks and concrete.
A stress-strain relation based on the strain equivalence hypothesis was proposed by Lemaitre 24 which supposes that the strain behavior of a damaged material can be represented by the constitutive equation of the material with no damage, by changing the stress to the net stress (effective stress). Kawamoto et al 25 developed a damage model to describe the mechanical behaviors of rock mass containing joints with intermediate-sized. In this model, the many joints with intermediate size (neither sufficiently small nor large compared to the structure) existing in the rock mass are considered as damage to the mass. By separately describing the properties of intact rock and the geometry of the joints, the mechanical responses of rock mass can be treated conveniently. Zhou et al [12][13][14][15] developed a micromechanics-based model to evaluate the complete stress-strain relationship and strength for crack-weakened rock subjected to compressive loading, dynamic tensile loading, and dynamic compressive loading, which is in good agreement with experimental results. Zhou et al 26 proposed a threedimensional long-term strength criterion of rocks based on micromechanical method of penny-shaped microcrack in soft rocks. However, most previous studies focused on damage characterization of elasto-plastic model or micromechanical model, but the influence of damage on the time-dependent behavior of jointed rock mass has rarely been investigated.
In our recent study, we have developed a geometric damage model in FLAC 3D which could reflect the influence of joint density and joint direction. 27 In this study, we extend the model into a damage rheological model for jointed rock masses and further developed this model in FLAC 3D software. We present the method of geometric damage mechanics, and the practical procedure for the long-term behavior of an underground opening.

THEORY FOR DISCONTINUOUS ROCK MASSES
For the one-dimensional case, the damage variable can be expressed as the ratio of effective area over total area. The effective area equals the total area minus the defects. However, in this case, the direction of effective area cannot be considered.
Murakami and Ohno 23 proposed a damaged tensor to define the effects of some groups of planar cracks: where Ω is the area density of a group of cracks, n is the unit vector normal to the cracks, and ⊗ indicates the tensor product.

| Damage tensor for rock masses
In order to apply the damage concept to rock mass behavior, Kawamoto et al 25 improved the original concept proposed by Murakami and Ohno 23 A rock mass is assumed to consist of cell elements of intact rock, and planar joints existing on the boundary surfaces of these elements. Joints could propagate along these boundary surfaces, but cannot penetrate into the intact rock elements. In this model, only the effect of initial damage is considered. Figure 1 shows a cell element and the effective surface area in a rock mass.
The dimensions of a cell element can be determined from average spacing of the joints. In Figure 2, v and V are the volumes of the cell element and the rock mass, respectively. The damage tensor in V can be obtained as 25 where l = v 1/3 , N is the number of joints in V, a k is the area of the kth joint, and n k is the unit normal.

| Normal vector of joints
The direction of each group of joints can be determined based on data observed on the surfaces. If a group of joints have angles θ 1 , θ 2 and θ 3 on coordinate surfaces X 1 , X 2 , and X 3 , respectively (Figure 3), the unit normal n = (n 1 , n 2 , n 3 ) of the group can be obtained as The unit normal n can also be defined as follows or

| Net stress
Since the reduction in the effective area is caused by the distribution of joints, we have where is the Cauchy stress, ̃ is the net stress, and the superscript "−1" denotes the inverse matrix.
The transformation law of Cauchy stress into net stress characterizes the mechanical effects of the damage.
Equation (8) can be further expressed as Equation (9) can be simplified and expressed as where The net stress is then given by 27 Connection between the kinematic field of a damaged body with the net stress can be built by the formation of a constitutive relation where Φ is the function of the constitutive equation, which will be introduced in section 3. ε is the total strain.
As aforementioned, the constitutive law for a rock mass is given in the form of Equation (14), while, if the material has no damage in the body, the net stress ̃ in Equation (14) is regarded simply as the Cauchy stress .
We define the damage zones as follows 27 : A shear damage zone will be defined in the numerical calculation, if the value of shear strain of an element exceeds the threshold of the shear strain.
A tensile damage zone will be defined in the numerical calculation, if the value of tensile strain of an element exceeds the threshold of the tensile strain.

MODEL FOR ROCK MASSES
In this study, the Nishihara model is adopted to model the damage rheological properties of rock masses. As shown in Figure 4, the Nishihara model is an elastic, visco-elastic, and visco-plastic model, composed of a Hookean body, a Kelvin body, and a Bingham body in series. This model is selected as it combines the advantages of the Bingham and generalized Kelvin models.
The Nishihara model can represent the instantaneous deformation, stable creep deformation range for low stress levels, and the unstable deformation after a specified critical stress level (long-term strength) σ s is reached. We adopted the Nishihara model as the basic model to study and interpret the rheological behavior of the rocks in the excavation zones. 28 For the Nishihara rheological model (Figure 4), components in the creep model meet the following formulation where σ is the total stress and ε is the total strain. σ e , σ ve , and σ vp are stresses for the Hookean body, Kelvin body, and Bingham body of the model, respectively; ε e , ε ve , and ε vp are strains for the three parts of the model; E 1 is the elastic modulus in the Hookean body and E 2 are the visco-elastic modulus in the Kelvin body; and η 1 and η 2 are the viscosity coefficients.
The one-dimensional constitutive equation of Nishihara model is: When: σ ≤ σ s : When σ >σ s : According to the principle of strain equivalence, we assume that damage effects are modeled by the transformation of Cauchy stress into net stress. The Cauchy stress σ in Equations (17a-17b) is replaced by the net stress ̃ . Therefore, the one-dimensional Nishihara visco-elasto-plastic damage rheological constitutive equation can be obtained.
In the three-dimensional stress state of the rock, stress tensor can be expressed as the sum of a spherical stress tensor σ m δ ij and a deviatoric stress tensor S ij . Similarly, the strain tensor can be expressed as the sum of a spherical strain tensor ε m δ ij and a deviatoric strain tensor e ij . The expressions are as follows: where δ ij is Kronecker symbol.
The constitutive relationship for an isotropic elastic medium can be written as where K is bulk modulus, G is shear modulus.
The three-dimensional Nishihara visco-elasto-plastic damage rheological constitutive equation can be expressed as: When S ≤ s : where G 1 represents the elastic modulus in the Hookean body, and G 2 represents the visco-elastic shear modulus in the Kelvin body; η 1 and η 2 represent the viscosity coefficients; e ij is the deviatoric strain tensor, and S ij is the net deviatoric stress tensor.

| Difference form of the visco-elasto-plastic damage rheological constitutive equation
The three-dimensional finite difference model is developed and applied in the commercial software FLAC 3D . In order to program this rheological model, the three-dimensional finite difference constitutive equation is derived. 29 Equation (16) is written in the incremental form: The visco-elastic part of the model originating from the Kelvin body can be expressed in the form of central difference: The elastic part of the model, originating from the Hookean body, can also be written in the central difference form as: where the overbar indicates the mean value over the time step Δt: After substituting Equations (24) and (27) into Equation (22), and solving for the new net deviatoric stress component, we get (using the mean value definitions Equations (25) and (26)): where: The volumetric deformation behavior is controlled by For the convenience of programming, the threedimensional Nishihara visco-plastic damage rheological constitutive equation is written in the central difference form: where the superscript N and O represent new and original values in one time step.
The failure criterion used in the FLAC 3D model is a combination of Mohr-Coulomb criterion with tension cutoff, 29,30 which is also used in this study. The details about the viscoplastic correction can be found in references. 29,30 In summary, the stress-strain relationship of the damage visco-elastic-plastic creep model can be expressed using Equations (13) and (32) which are used to program the model.

| Implementation of the damage rheological model
FLAC 3D adopts secondary development by users, and all constitutive models are provided to users as Dynamic Link Library (DLL) files. We programed this damage rheological model in FLAC 3D with the object-oriented language C++. During the computation process, the main program automatically calls the DLL files of the new damage rheological model created by us. Figure 5 shows the program flow diagram of the damage rheological model. Sterpi 31 and Yang et al 30 introduced the methods to determine creep parameters of rock and rock mass. In our study, we determined parameters of rock by triaxial compressive creep tests.

| Parameters of long-term strength
Several triaxial creep tests of specimens were carried out to obtain strain-time relationship (Figure 6A), and isochronous stress-strain curves ( Figure 6B) can be obtained. Figure 6A shows a series of creep strain-time curves under different shear stress levels. The dashed vertical lines were marked at different time and intersect with all curves. The intersection points contain three parameters (stress, strain, and time) for generating the Figure 6B. Each isochronous curve is approximately composed of linear and nonlinear segments, and there is an obvious turning point which decreases with time and tends to reach a certain limit (the red dash line in Figure 6B). The value of this limit is defined as long-term strength of rock (σ s in Figure 6B). Figure 7 shows the variation of long-term triaxial compressive strength as a function of confining stress. Based on Mohr-Coulomb strength criterion, we can get the long-term cohesion c ∞ and long-term friction angle ϕ ∞ . The long-term tensile strength σ t∞ can be obtained by the direct tensile creep test or Brazilian splitting creep test.

| Inversion of the creep parameters
The stepwise iterative optimization inversion method was used to numerically back analyzed the creep parameters from the creep strain data of the rock specimens. In the creep model, there are 5 other creep parameters to be solved except the parameters of long-term strength. We used the following procedures to conduct parameter inversion:

The values of creep strain under a certain level of stress at
time t i were calculated by numerical simulation taking into account the 5 creep parameters assumed in step 1. 3. The objective function for the inversion of creep parameters was established. In this study, the sum of residual square between the calculated and measured creep strain at time t i was taken as the objective function, that is where N is the measured strain; u i (X, t i ) is the strain at time t i calculated from numerical analysis results; and ui is the measured strain at time t i . 4. The objective function (33) was judged whether it reached the minimum value. If yes, the selected creep parameter X is the result sought by the inversion analysis. Otherwise, the design variable will be continuously adjusted and step 2~3 will be iterated, until the objective function reaches the minimum value.

| Parameters of joints
The parameters of joints can be determined as follows: We can observe dominant directions and scales of plane joints in a rock mass. We choose the volumes of the rock mass V and the length of the side of the cell element l according to engineering geology. The damage tensor can be obtained for each set of joints by in situ observation, and we need to collect the data of joints on three independent surfaces of the rock mass ( Figure 8). (X 1 , X 2 , X 3 ) is a local coordinate system. Suppose some sets of joints appear on the coordinate surfaces X 1 , X 2 , and X 3 . The angle θ 1 of a joint on the X 1 -surface is measured from the X 2 -axis in the counterclockwise direction. Similarly, θ 2 and θ 3 can be obtained.
If a rock mass involves multiple sets of joints, we can separately determine the unit normals and the damage tensors for the set. More details are given by Kawamoto 25 and Yang. 27

F I G U R E 8 Observation of cracks on rock mass surfaces (after
Ref. [25] F I G U R E 9 Mesh sketch of the model  Table 1, c ∞ is the long-term cohesion, ϕ ∞ is the long-term friction angle, and σ t∞ is the long-term tensile strength. Properties of joints are listed in Table 2.

T A B L E 1 Material properties of rock masses in numerical calculation
The complete creep curve of the damage rheological model is shown in Figure 10. The creep curve shows different forms at different levels of σ 1 stress. We can see transient creep (visco-elastic displacement) appears at σ 1 =30, 40, 50, 60, 70 MPa. Steady state creep appears at σ 1 =80 MPa.
Because the damage effect of joints is considered in damage rheological model, the magnitudes of displacements at the applied compressive stresses of 30 MPa to 70 MPa are greater than that of Cvisc model. The displacement rates from the damage rheological model and the Cvisc model in FLAC 3D at the applied compressive stresses of 30 MPa to 70 MPa are exactly the same, confirming the validity of the visco-elastic part of the damage rheological model. However, the Cvisc model does not have any steady state creep part; that is, the rock specimen fractures immediately when the rock comes to plastic yield. Once the applied stress reaches 80 MPa, the displacement of the Cvisc model occurs in the plastic yield stage and increases rapidly. Therefore, the Cvisc model cannot explain the steady state creep caused by visco-plastic properties. However, the damage rheological model can represent the primary creep caused by viscoelastic displacement, steady state creep caused by visco-plastic displacement and also can represent the damage displacement caused by joints. Overall, the damage rheological model can better replicate the creep behavior than the models in FLAC 3D .

| Case 2
We also performed three creep experiments on rock-like material under uniaxial compressive stress states to verify the damage rheological model. Rectangular prismatic blocks with the scale of 300 mm × 150 mm × 150 mm (height × width × thickness) were generated to conduct uniaxial compressive tests. A cement mortar consisting of C32.5 Portland cement, fine sand, and water following the ratio 1: 2.34:0.55 was selected as the rock-like material. The joints were created by inserting steel shims (1.5 mm in thickness, 200 mm in length, and the width depending on the designed width of joints) into the fresh cement mortar paste at the desired location of the joints; the shims were removed after the early condensate of mortar (about 24 hours). The joint geometry was defined by four parameters: joint angle (α), spacing (S), joint length (L), and rock bridge length (B) (Figure 11). Three creep experiments (J1, J2, and J3) with different joints distribution were performed. The joint geometry of specimens is shown in Table 3.
Based on the joint geometry, the damage tensor can be obtained as follows: Then, we calculated the creep strains of samples and compared them with experimental results. Creep properties of samples in numerical calculation are shown in Table 4. Comparison  Figure 12.
From Figure 12, we can conclude that the damage rheological model is capable of replicating the long-term mechanical properties of jointed rock specimens subjected to uniaxial compressive stress states.

| ENGINEERING APPLICATION
We applied this damage rheological model to predict the displacements and fractured zones of a circular underground opening in jointed rock masses after excavation. Figure 13 shows the schematic diagram illustrating an underground opening in jointed rock masses. The scale of the model is 100 m × 100 m (length × height), and the diameter of the circular opening is 10 m. The displacements of the side walls in the X direction and bottom wall in the Y direction are constrained, and the top of the model is a free boundary. The Cvisc model in FLAC 3D was also used for the purpose of comparison.
The lateral coefficients of the initial stress are ϑ z = 0.5 along the opening axis direction (the Z direction) and ϑ x = 0.5 along the direction perpendicular to the opening axis (the X direction). Table 5 lists the material properties of rock masses surrounding the underground opening.
First, we carried out a conventional creep numerical simulation using Cvisc model to analyze the excavation of an underground opening without considering the influence of joints in rock masses. Figure 14 shows the result of creep analysis using Cvisc model. We can see that the opening vault and bottom have large displacements. This is because the horizontal stress is smaller than overburden stress (σ x = 0.5σ y ) before the excavation of the opening.
The figures show that the position of maximum displacement and failure zone shapes change with the distribution of joints, which indicates the damage effect of joints on rock masses.
The displacements (elastic and creep displacements) and fractured elements are not only influenced by the joint angle, which is shown in Figures 15A1-C1 and 15A2-C2, but also by the joint area density, which is shown in Figures 15D1  and D2. A comparison of the results of condition 3 and condition 4 indicates that the greater the joint area density, the greater the displacements and fractured zones.
Displacements and depths of fractured zones of monitoring points A, B, and C are listed in Table 7, and the underline represents the maximum values in A, B, and C. For condition 1 with joint angle θ k = 0°, the maximum displacement and depth of fractured zones appear at point A, whereas the maximum values appear at point C for condition 2 with θ k = 90°. Similarity, the maximum values occur at point B when θ k = 135° for condition 3, which shows the damage rheological model could reflect the influence of joint angle. In addition, compared with condition 3, condition 4 has larger joint area density causing greater displacement and fractured zones, which shows the influence of joint number in this model.

| CONCLUSIONS
Joints play a major role in determining the long-term behavior of structures constructed in rock mass. Based on the basic theory of damage mechanics, we develop a damage rheological model and apply it in describing the long-term behavior of rock masses containing distributed joints. The following conclusions can be drawn: The underline means the maximum value for each condition.
T A B L E 7 Displacements and depths of fractured zones of monitoring points A, B, and C acts on the net surface excluding the damage. The damage tensor is able to describe the joint distribution states in rock mass. By introducing the net stress acting on the net surface excluding the damage, the mechanical effects of joints can be captured.

2.
A new damage rheological model was proposed by combing the damage tensor with Nishihara rheological model (a visco-elastic-plastic model). In the current model, only the effect of initial damage is considered, and the damage evolution with time will be discussed in future work.

3.
The new damage rheological model was applied in FLAC 3D , providing a reference case for the development of other models. The model was validated by conducting a numerical analysis which involved some groups of patterned joints. We illustrated the validity of the damage rheological model through numerical simulations.