# Micro-crack propagation and coalescence during time-dependent deformation of granite based on numerical manifold method

Xian-Yang Yu, Tao Xu, Michael Heap, Patrick Baud, Zhen Heng

► To cite this version:

Xian-Yang Yu, Tao Xu, Michael Heap, Patrick Baud, Zhen Heng. Micro-crack propagation and coalescence during time-dependent deformation of granite based on numerical manifold method. IOP Conference Series: Earth and Environmental Science, 2020, 570, pp.022064. 10.1088/1755-1315/570/2/022064. hal-03446357

HAL Id: hal-03446357

https://hal.science/hal-03446357

Submitted on 28 Nov 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
PAPER • OPEN ACCESS

Micro-crack propagation and coalescence during
time-dependent deformation of granite based on
numerical manifold method

To cite this article: Xian-yang Yu et al 2020 IOP Conf. Ser.: Earth Environ. Sci. 570 022064

View the [article online](#) for updates and enhancements.

You may also like

- [An optimized scheme of dispersion suppression for elastic-wave variable-order rotated staggered-grid forward modeling](#)
Weizhong Wang, Tianyue Hu, Jianyong Song et al.

- [NMMGenerator: an automatic neural mass model generator from population graphs](#)
Maxime Yochum and Julien Modolo

- [Application of GMLS-based numerical manifold method in buckling analysis of thin plates](#)
Hongwei Guo and Hong Zheng

<figure><img src="image_2.png" /></figure>

The Electrochemical Society
Advancing solid state & electrochemical science & technology

241st ECS Meeting

May 29 – June 2, 2022 Vancouver • BC • Canada

Abstract submission deadline: Dec 3, 2021

Connect. Engage. Champion. Empower. Acclimate.
We move science forward

<figure><img src="image_3.png" /></figure>

Submit your abstract
# Micro-crack propagation and coalescence during time-dependent deformation of granite based on numerical manifold method

Xian-yang Yu¹, ², Tao Xu¹*, Michael J. Heap², Patrick Baud² and Zhen Heng¹

¹Center for Rock Instability and Seismicity Research, Northeastern University, Shenyang 110819, China

²Laboratoire de Géophysique Expérimentale, Institut de Physique de Globe de Strasbourg (UMR 7516 CNRS, Université de Strasbourg/EOST), 5 rue René Descartes, 67084 Strasbourg cedex, France

*Corresponding author: xutao@mail.neu.edu.cn

**Abstract.** An understanding of influence of micro-cracks on time-dependent deformation in granite is of fundamental importance in many situations, such as nuclear waste storage facilities or the exploitation of geothermal resources. Time-dependent cracking is considered to be the main mechanism of brittle creep in granite. Creep strain rates are strongly influenced by the density of pre-existing defects that also exert a significant influence on rock physical properties. We introduce a new model that combines the subcritical crack growth (SCG) theory and the numerical manifold method (NMM) to link the local damage (using an exponential material softening law) caused by micro-crack propagation and the macroscopic creep deformation typically observed in granite specimens. In this model, each element contains a virtual micro-crack that can propagate sub-critically following Charles' law. Once the virtual micro-crack length reaches a given value, it will convert to a real micro-crack, which can cut through adjacent elements, open, and slide according to the principle of NMM. We also investigated the influence of virtual micro-crack length, confining pressure and differential stress on creep behaviour. The fact that numerical simulations are in good agreement with experimental results shows that the NMM combined with the SCG theory is a suitable method for modelling the creep behaviour of rocks.

## 1. Introduction

The long-term stability of engineering rock mass structures have been of practical concern for many years[1]. The time-dependent deformation (creep) of rock has a significant influence on the stability of rock slopes and underground structures such as mines and tunnels, and the Earth's upper crust[2]. Creep strain and strain rates are greatly sensitive to differential stress, confining pressure, temperature, etc. Even small changes in any of these parameters will produce order of magnitude changes in creep strain rates[3]. Many researchers have found that micro-crack propagation in a rock mass is one of main parameters that influences the creep strain rate[4]. When slow crack growth (subcritical crack growth (SCG)) occurs, the phenomenon that the extent of slow crack growth affects the fracture stress was first examined by Charles[5]. The SCG theory was initially attributed to stress corrosion. During crack propagation, stress corrosion can explain a certain proportion of relationship between stress intensity

<figure><img src="image_4.png" /></figure>

Content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
factor (K) and crack velocity (V)[6]. SCG theory in rock masses has been also used to explain the growth and development of joints, volcanic eruptions, and underground excavations[7]. The double torsion experimental method was adopted to investigate the SCG theory. There are also many researchers that have used numerical simulations to study the rock creep and SCG. Brantut[8] used a micromechanical model to describe the brittle rock creep, based on the sliding wing-crack model. Konietzky[9] used FLAC to simulate the time-dependent crack growth of granite. These authors used the SCG theory of a virtual crack within an element to describe rock damage. Discontinuous deformation analysis (DDA), analyzing the force-displacement interactions of block systems and proposed by Shi[10], has been also used to investigate creep problems[11]. Xu[12] proposed a thermomechanical time-dependent deformation model based on tests on Beishan granite under constant elevated temperatures. The numerical manifold method (NMM) is a new simulation method which provides a unified framework for solving problems dealing with continuous media, discontinuous media, or both[13]. The NMM has gained a wide attention and application in rock mechanics and engineering due to its efficient treatment of problems involving continuous and discontinuous deformations in a unified way[14]. Using the time step–initial strain method, the creep equation was coupled with the NMM to simulate the time-dependent deformation of rocks[15]. Wu[16] also used the NMM to analyze viscoelastic material creep crack problems by incorporating a generalized Kelvin-Voigt model into the NMM.

## 2. Constitutive model

### 2.1. Crack Initiation and Growth Criteria

Before using the NMM to consider crack problems, the criteria that have an important influence on calculating the crack initiation and propagation angle must be considered. In this paper, there are two different crack growth criteria introduced, which are the maximum circumferential stress criterion (MCSC)[17] and Mohr–Coulomb’s shear strength criterion (MCSSC)[18]. For MCSC, the crack propagates in the direction of the maximum circumferential. The stress field at the crack tip can be formulated in terms of stress intensity factors (SIFs):

$$K_1 \sin \theta_0 + K_{II} (-1 + 3 \cos \theta_0) = 0 \quad (1)$$

and then, the direction of maximum circumferential stress is:

$$\theta_0 = 2 \arctan \frac{1}{4(K_{II}/K_1)} \left\{ 1 \pm \left[ 1 + 8(K_{II}/K_1)^2 \right]^{1/2} \right\} \quad (2)$$

The MCSC can well predict the fracture direction in two mixed modes if there are several cracks in the specimen. In nature, however, there are many cracks, joints, and holes in rock masses, and every flaw will influence the local stress field, which can influence the direction of maximum circumferential stress greatly. We used the MCSSC to deal with this problem.

In classical engineering science, the SIFs are greatly sensitive to the local stress. Coulomb postulated that the crack initiation depends on the local stress relative to the strength of the material rather than on the SIFs[19]. In this framework, the material strength parameters are the cohesion and the angle of internal friction, as shown in Figure 1. The MCSSC can be used to calculate shear and tensile cracks. As shown in Figure 1, $\sigma_1$ and $\sigma_3$ are the maximum and minimum principle stress, respectively, $\sigma_t$ is the tensile stress, and R and r are expressed as:

$$\left\{ \begin{array}{l} R = c \cos \phi + \frac{\sigma_1 + \sigma_3}{2} \sin \phi \\ r = \frac{\sigma_1 + \sigma_3}{2} \end{array} \right. \quad (3)$$

where c and $\phi$ are cohesion and the internal friction angle, respectively.
<figure><img src="image_5.png" /><figcaption>Figure 1 The Mohr–Coulomb criterion for the elements</figcaption></figure>

When the crack satisfies the shear condition:

$$ \begin{cases} R \le r \\ \sigma_3 > -\sigma_t \end{cases} \qquad (4) $$

the crack initiation angle will be:

$$ \theta = \pi/4 + \varphi/2 \qquad (5) $$

where $\varphi$ is the crack initiation angle.

When the crack satisfies the tensile condition:

$$ \begin{cases} R > r \\ \sigma_3 < -\sigma_t \end{cases} \qquad (6) $$

the crack initiation angle will be the direction of the minimum principal stress.

## 2.2. Subcritical Crack Growth (SCG)

In classical fracture mechanics theory, the crack will not propagate if the SIFs at the crack tip are lower than the fracture toughness $K_c$. However, some researchers also find that cracks still propagate if the SIFs are lower than $K_c$, but at a certain lower velocity. This is called the subcritical crack growth (SCG) theory[20]. Figure 2 shows the SCG theory curve between stress intensity factor (K) and crack velocity (v). Figure 2 clearly shows that there are three stages and that $K_{scc}$ is a critical value, which is related to the material. When $K < K_{scc}$, the crack cannot propagate. However, when $K_{scc} < K < K_c$, there are three stages of crack velocity (Figure 2).

<figure><img src="image_6.png" /><figcaption>Figure 2 Schematic stress intensity factor (K) against crack velocity (v) curves</figcaption></figure>

The crack velocity is governed by Charles’ equation:

$$ \begin{cases} v = CK^n \\ K = \sqrt{K_1^2 + K_{II}^2} \end{cases} \qquad (7) $$

where C is a material constant, K is the total SIF, which consists of $K_1$ and $K_{II}$, n is the stress corrosion index, v is the crack velocity, and $\Delta t$ is the time step. The crack length at every step is $v*\Delta t$.

The SIFs $K_1$ and $K_{II}$ are calculated as:

$$ \begin{cases} K_1 = \sigma_n \sqrt{\pi a} \\ K_{II} = \tau_n \sqrt{\pi a} \end{cases} \qquad (8) $$

When the crack at the tension-shear condition, where $\sigma_n$ and $\tau_n$ are the maximum tensile normal stress and maximum shear stress, respectively, a is the half crack length, and $\sigma_n$ and $\tau_n$ are expressed as:
$$ \left\{ \begin{aligned} \sigma_n &= \frac{1}{2}[(\sigma_1 + \sigma_3) + (\sigma_1 - \sigma_3)\cos 2\varphi] \\ \tau_n &= \frac{1}{2}(\sigma_1 - \sigma_3)\sin 2\varphi \end{aligned} \right. \qquad (9) $$

where $\varphi$ is the crack original angle.

When the crack at the compress-shear condition, where $\tau_{eff}$ is effective shear stress, f is the friction coefficient.

$$ \tau_{eff} = |\tau_n| - f|\sigma_n| \qquad (10) \\ K_{ll} = |\tau_{eff}|\sqrt{\pi a} $$

## 2.3. Numerical Manifold Method (NMM)

The NMM consists two of important components: the cover system and the block dynamics. The cover system is a dual cover system, which includes a mathematical cover (MC) of three-node triangle elements and a physical cover (PC) that includes the boundary, the material interface, and the crack. The block dynamics is used to solve the mechanical behaviour of block systems under loading and block contact. The NMM program runs on the MATLAB platform.

Figure 3 shows the dual cover system containing MC and PC. We can clearly see in Figure 3(a) that there are three MCs that are coloured yellow, blue and red, that every MC is a hexagon consisting of 6 triangles, and that the node of every triangle is the centre of the MC. If there a crack that cuts through the triangle ABC (such as D1-D5 in Figure 3), the three MCs are changed. It is clearly shown in Figure 3(b) and (d) that the crack does not cut the MC (A and C) into two parts. Instead, the hexagon and a part of crack consist of a PC, which are A1-A2-A3-A4-D3-D4-D5-D4'-D3'-A5-A6-A1 and C1-C2-C3-C4-C5-D2-D3-D4-D5-D4'-D3'-D2'-C6-C1, respectively. In Figure 3(c), the crack cut the MC (B) into two parts (B and B'), so, there are two PCs, which are B1-B2-D4'-D3'-D2'-D1'-B5-B6-B1 and D4-B3-B4-D1-D2-D3-D4, respectively.

<figure><img src="image_7.png" /><figcaption>Figure 3 Dual cover system and crack cutting theory of NMM</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 4 Three contact conditions of NMM</figcaption></figure>

When the crack is opened in the program, contact theory plays an important role to prevent elements embedding with two springs, which are tension and shear springs. In the NMM, there are three contact conditions, as shown in Figure 4: (a) angle to angle, (b) angle to edge, and (c) edge to edge. Both angle to angle and edge to edge contacts can be converted into angle to edge contacts. As shown in Figure 3, $d_0$ is the one step maximum displacement (safe displacement) and $d_1$ is the general displacement between angle to edge. If $d_1>d_0$, the program does not need to adjust the model. However, if $0<d_1<d_0$ (or $d_1<0$, which means the angle is embedded into the edge), the program uses the open-close iterations method[18] and time steps to push the angle back and keep the displacement larger than $d_0$.

## 2.4. Multi-crack NMM model

A pre-cracked granite specimen model of 50 mm in diameter and 100 mm in length is shown in Figure 5(a). The half real crack length is 2.5mm, and the angle of the crack $\varphi$ is 45°. The meshed specimen that
contains multi-cracks is shown in Figure 5(b). The three-node triangular elements in Figure 5(b) consist of mathematical covers, and the solid black line in Figure 5(a) represents the physical boundary and discontinuities in the specimen. The red solid and dotted lines in Figure 5(c) and (d) represent real and virtual cracks, respectively[9]. The model is consisted of 1108 elements, which means that this model has more than 1000 virtual cracks, every crack has two crack tips, and every crack tip can calculate the SIF at every simulation step. The half initial virtual crack length is 0.1mm, it's smaller enough comparing with real crack. The real crack is a part of the physical cover, and it can therefore cut the element and open and slide. There are no virtual cracks in the area of the real crack, because these elements are damaged. Virtual cracks are only present in the three-node triangular elements, and will not cut the elements. The angle and length of the virtual cracks is arbitrary, and their lengths are equal.

<figure><img src="image_9.png" /><figcaption>Figure 5 NMM model of size and multi-crack distributions</figcaption></figure>

When the SIFs at the crack tips are lower than the fracture toughness Kc, both the real and virtual cracks will propagate according to the SCG theory. The elastic modulus of an element in which the virtual and real crack tips are located will be damaged according to Equation (11)

$$ \begin{cases} E' = D * E \\ D = (1 - L'/L)^n \\ L' = v * \Delta t \end{cases} \qquad (11) $$

where E is the last step elastic modulus, E' is the damage elastic modulus, D is the damage factor, L is the original crack length, L' is the crack propagation length, and n is the damage index.

Because the real crack length is larger than that of a virtual crack, the SIFs of the real crack will be equal to Kc firstly and then propagate according to the MCSSC. The virtual crack is very small and, because the crack length at every step is very small, crack growth will occur along the original direction and will take a long time. Before the virtual crack can change into a real crack, the main role of the virtual crack is to reduce the elastic modulus. The virtual crack will not cut the element. Once the virtual crack length propagates to a pre-defined value (related to model size and mesh density), it will be converted into a real crack according to the principles of NMM. During crack propagation, the boundary and node of the element cannot influence the crack direction. The crack can grow across the element boundary, and crack tip can be located on either the boundary, node, or in the element, which is more freedom to simulate crack propagation.

### 3. Simulation

The improved NMM program was used to simulate the time-dependent deformation of a rock specimen under constant differential stress. The macroscopic mechanical properties of rock are shown in Table 1[21]. An axial stress and a confining pressure of 173.9 and 2 MPa, respectively, were applied to the numerical specimen. The bottom of specimen was fixed.

<table><caption>Table 1 Physico-mechanical parameters of numerical model</caption><thead><tr><td>Items</td><td>Value</td></tr></thead><tbody><tr><td>Elastic modulus: E(GPa)</td><td>72</td></tr><tr><td>Poisson's ratio: υ</td><td>0.25</td></tr><tr><td>Fracture growth constant: C</td><td>1*10<sup>-25</sup></td></tr><tr><td>Stress corrosion index: n</td><td>19.8</td></tr><tr><td>Time step: Δt (h)</td><td>0.25</td></tr></tbody></table>
<figure><img src="image_10.png" /></figure>

**Figure 6** Comparison between experimental and numerical creep strain and strain rate curve of granite The classical creep simulation methods are the empirical method and the rheological component method, which use differential equations to realise the stages of rock creep. In this research, when the SIFs of crack tips are smaller than Kc, both virtual and real cracks propagate according to the SCG theory. During crack growth, the SIFs increase with the increasing of crack length, which can also increase the crack velocity according to Equation (7). The elastic modulus is decreased according to Equation (10) at the same time. When the SIFs are smaller than Kc, the axial strain decreases as an increasing function of time (i.e. the decelerating creep phase). When the SIFs are larger than Kc, the real crack propagates according to the MCSSC. The strain rate also shows decelerating and accelerating creep (Figure 6(b)). The comparison curve of experimental (blue squares) and numerical (red circles) creep strain and strain rate are shown in Figure 6. The experiment is creep test performed on Beishan granite under biaxial compression and room temperature, from Chen (2015), and the numerical model was performed under the same conditions. The experimental results clearly show the two stages of creep curve, decelerating and accelerating creep, and the numerical results are in good agreement with the experiment. In Figure 6(a), the initial axial strain is about 0.24%. In the first 2 h, the initial crack length is very small and, according to Equation (7), the cracks are growing at a slow rate. As a result, the strain increases slowly, and the strain rate decreases at the same time (Figure 6(b)). From 2 to 4 hours, the axial strain increases (Figure 6(a)), and the strain rate decreases at a very slow rate (Figure 6(b)). After 4 h, the axial strain and strain rate increase quickly, the cracks propagate fast and connect together to generate a macro crack zone, which, eventually, results in macroscopic sample failure.

Figure 7 shows the influence of virtual crack length (V-L) and differential stress (D_S) on the time-dependent deformation of granite. As shown in Figure 7(a), the creep strain increases with increasing V-L. The larger the V_L, the sooner the specimen reaches the accelerating creep stage. According to Equation (7) and (8), V-L and SIF are positively correlated: the greater the V_L, the larger the SIF. Further, with an increase of SIF, the velocity of SCG also increases. So, the propagation of micro-cracks can lower the strength of the specimen. Figure 7(b) shows the influence of confining pressure (C_P) on the axial creep strain. C_P can increase the rock strength: the greater the C_P, the smaller the axial strain. C_P can also make the specimen more stable: the larger the C_P, the more difficult it is to break the specimen. In Figure 7(c) and (d), D_S can also influence the axial strain. Even a small increase in D_S can change the axial strain of the specimen, and the failure time of specimen is also decreased. It is easy to explain the strong effect of differential stress on time-to-failure and creep strain. An increase in the D_S can induce an increase of the SIFs at micro-crack tips. Even a modest increase in SIF results in a large increase in the velocity of SCG according to Equation (7).

<figure><img src="image_11.png" /></figure>

**Figure 7** Influence of virtual crack to time-dependent deformation: (a) virtual crack length (V-L), (b) confining pressure (C_P), (c) differential stress (D_S), (d) the influence of D_S to time to failure According to the SCG theory, when K is smaller than Kc, both virtual and real cracks propagate at a very low speed. It is therefore not easy to observe the crack propagation process. On the other hand, if K is larger than Kc, the real crack will propagate according to the MCSSC. Figure 8 shows the real crack
propagation path under time-dependent deformation progress. Figure 8(a) shows the initial real crack location. As SCG progresses, the real crack propagates according to the MCSSC, because $\sigma_1$ is larger than that of $\sigma_3$. As a result, the crack extends in the direction of $\sigma_1$. Because of NMM theory, the crack can propagate by cutting through the elements (Figure 8(c)), which is different from FEM and DEM in which the crack can only propagate along the boundary of the elements. Crack tips in this simulation can not only stay on the node or boundary, but can also stay within the element. After the crack propagates, the crack face is opened too. The program can calculate the contact state at every node of the crack by adding and subtracting the spring to make sure the crack surface is not embedded according to contact theory.

<figure><img src="image_12.png" /><figcaption>**Figure 8** Crack propagation path under time-dependent deformation when K is larger than that of Kc: (a) initial real crack location, (b) real crack propagation path, (c) zoom in of crack open inflection</figcaption></figure>

**4. Results**

We propose that SCG theory can combine with NMM to simulate time-dependent brittle deformation under different loading conditions and different confining pressures. In order to reflect the influence of crack propagation on creep deformation, we use Charles' equation to calculate the crack velocity. The evolution of the crack velocity has three stages. After the validation of SCG theory, we performed a parametric sensitivity analysis to verify the correctness and characteristics of the creep model and find that the virtual crack length (V-L), differential stress (D_S), and confining pressure (C_P) all play a major role in governing elastic deformation. With an increase of V_L, the specimen can reach the accelerating creep stage quickly, and C_P can increase rock strength and slow down the creep deformation progress. Even a small change of D_S can change the axial strain and failure time. After the sensitivity analysis of the parameters, we use this model to realize crack propagation under creep conditions. We find that the cracks propagate along the maximum principle stress direction.

**5. Conclusions**

In this paper, the time-dependent deformation of granite is modelled by combining SCG and NMM on the MATLAB platform. This combined method replaces the traditional creep component model method. In the model, cracks propagate according to the SCG theory. When K is smaller than the fracture toughness Kc, the crack growth length is linked to the local damage (using an exponential material softening law) to reduce the granite strength. When K is larger than Kc, the cracks propagate according to the MCSSC. Parameter sensitivity analysis for the improved NMM program was performed and validated against experimental results. Crack length, confining pressure and differential stress were then investigated. The results show that the microstructure of granite has a great influence on time-dependent deformation. Last but not least, numerical simulations are in good agreement with the experimental results and therefore show that the improved NMM program is suitable for modeling the time-dependent creep deformation of rocks.

**Acknowledgments**

The work was jointly supported by NSFC (51974062, 41672301, 51950410595), National Key Research and Development Program of China (2017YFC1503100), Fundamental Research Funds for the Central
Universities of China (N180101028) and the financial support from the China Scholarship Council (File No. 201906080054). The comments of two anonymous reviewers helped clarify aspects of this manuscript.

## References

[1] Dawson P.R. and Munson D.E. 1983. *International Journal of Rock Mechanics & Mining Sciences & Geomechanics Abstracts.* **20** 33-42

[2] Tsai L.S., Hsieh Y.M., Weng M.C., Huang T.H., et al. 2008. *International Journal of Rock Mechanics and Mining Sciences.* **45** 144-154

[3] Brantut N., Heap M.J., Meredith P.G. and Baud P. 2013. *Journal of Structural Geology.* **52** 17-43

[4] Mongi K. 1962. *Bulletin of the Earthquake Research Institute.* **40** 125-173

[5] Charles R.J. 1958a. *Journal of Applied Physics.* **29** 1549-1553

[6] Evans A.G. 1972. *Journal of Materials Science.* **7** 1137-1146

[7] Olson J.E. 1993. *Journal of Geophysical Research: Solid Earth.* **98** 12251-12265

[8] Brantut N., Baud P., Heap M.J. and Meredith P.G. 2012. *Journal of Geophysical Research Solid Earth.* **117** 1133-1172

[9] Li X. and Konietzky H. 2015. *Acta Geotechnica.* **10** 513-531

[10] Shi G.H., *Discontinuous Deformation Analysis A New Numerical Model for the Static and Dynamics of Block Systems*, in *PhD Dissertation, Dept. of Civil Engineering*. 1989.

[11] Gao Y., Gao F. and Ronald Y.M. 2013. *International Journal of Mining Science and Technology.* **23** 757-761

[12] Xu T., Zhou G.L., Heap M.J., Zhu W.C., et al. 2017. *Rock Mechanics and Rock Engineering.* **50** 2345-2364

[13] Shi G.-H., *Manifold method of material analysis*. 1991, Army Research Office Research Triangle Park NC. p. 51~76.

[14] Ning Y.J., An X.M. and Ma G.W. 2011. *International Journal of Rock Mechanics and Mining Sciences.* **48** 964-975

[15] Yu X.Y., Xu T., Heap M., Zhou G.L., et al. 2018. *International Journal of Geomechanics.* **18** 04018153

[16] He J., Liu Q.S. and Wu Z.J. 2017. *Engineering Analysis with Boundary Elements.* **80** 72-86

[17] Erdogan F. and Sih G.C. 1963. *Journal of Basic Engineering.* **85** 519-525

[18] Wu Z. and Wong L.N.Y. 2012. *Computers & Geotechnics.* **39** 38-53

[19] Brady B.H.G. and Brown E.T., Rock Mechanics for underground mining [M]. 1993.

[20] Atkinson B.K. 1984. *Journal of Geophysical Research Solid Earth.* **89** 4077-4114

[21] Chen L., Liu J.F., Wang C.P., Liu J., et al. 2015. *European Journal of Environmental and Civil Engineering.* **19** S43-S53