Stress intensity factors for a tunnel containing a radial crack under compression

Abstract

In this article, based on crack surface boundary conditions, the complex stress functions for a tunnel containing a radial crack under compression was obtained, and boundary collocation method was applied to the determination of the unknown coefficients of the complex stress functions. The stress intensity factors were obtained through the complex stress functions, and the results have been validated by numerical simulation using ABAQUS code and by the experiments of tunnel models under compression. The results show that the numerical results using ABAQUS code and the experiment results agree well with the theoretical results by the complex functions, which indicates that the analytical method for determining the stress function and the stress intensity factors proposed in this article are effective and reliable.

Introduction

With increasing scale of tunnel constructions in mining and transport engineering, a series of tunnel accidents have been exposed, such as rock burst and roof fall. Natural joints or cracks are frequently encountered in rock mass, and they are distributed randomly around tunnels. Meanwhile, in the process of rock excavation, under the action of explosive detonation, blast-induced radial cracks will occur in the surrounding rock,^1–3 which could largely weaken tunnel stability. It is expected that as crack orientation differs, the effect extent on tunnel stability will vary, but these are not clear yet. In order to clearly understand the effect of cracks on tunnel stability, it is essential to implement the corresponding theoretical, experimental, and numerical studies so as to obtain a better understanding of the dominant parameters that control tunnel stability.

Using complex stress functions, the basic equations for solving the plane problems of elasticity were proposed by Muskelishvili ⁴ in his famous book. For the problems of a crack emanating from a hole, researchers usually used analytical method and numerical method to study. For example, Bowie ⁵ used a conformal mapping technique to study the problem of a radial crack situated at the edge of a circular hole. Berezhnitskii ⁶ studied the problem of a crack emanating from an elliptical hole using the integral equation method. Tweed and Rooke ⁷ studied this problem using Mellin transform technique, but this method is unsuitable to treat complex configurations. Rubinstein and Sadegh ⁸ considered a circular hole with an arbitrary positioned edge crack using dislocation density distribution function. Woo ⁹ studied the mixed mode problem of cracks emanating from a circular hole in a finite plate. Wang et al. ¹⁰ proposed a group of stress functions for the problems of a crack emanating from a hole with different shapes (including circular, elliptical, rectangular, or rhombic hole). Wang ¹¹ studied the stress intensity factors (SIFs) for an infinite plate with many collinear radial cracks emanating from many circular holes under biaxial loading. Wang and Chau ¹² derived a new boundary integral equation for plane elastic bodies containing cracks and holes and subjected to mixed displacement/traction boundary conditions, and proposed a new boundary element method based on this formulation. Yan ¹³ studied SIFs for cracks emanating from a triangular or square hole in an infinite plate subjected to internal pressure by means of a boundary element method. Zhao et al. ¹⁴ proposed an analytical method for calculating the SIF of an infinite plate containing multiple hole-edge cracks by complex variable functions. Through photoelastic experiments, Wang et al. ¹⁵ analyzed the SIFs of cracks around a tunnel. Through tunnel model tests and numerical analysis, Zhu et al. ¹⁶ and Li et al.^17,18 studied the effects of principal stress orientation and crack orientation on tunnel stability, and Fan et al. ¹⁹ investigated the effect of the ratio of tunnel depth to its width on tunnel stability. Using multi-domain boundary element method, Wang et al. ²⁰ studied the influence of cracks on the stability of the surrounding rock of a laneway. Although the problems of a crack emanating from a hole have been studied by many researchers, less attention has been paid to the issue of a crack emanating from a horseshoe-shaped tunnel edge and the effect of crack orientation on tunnel stability.

In this article, a horseshoe-shaped tunnel in an infinite rock mass with a radial crack emanating from the tunnel edge is considered which is a through-thickness crack, and for simplicity, a plane strain tunnel model is established to yield the analytical solution of SIF. The tunnel model is under biaxial compression, which is different from those references mentioned above in which the loads were tensile. Under compression, cracks in rock mass are closed in general, and crack surface friction will affect crack tip SIFs significantly; thus, it should be considered in the study of cracks under compression.^21–24 The tunnel is inside an infinite plane, and the corresponding complex stress function can be obtained based on the stress conditions of crack surfaces. However, the stress functions consist of a group of unknown coefficients, and they have to be determined in order to obtain the complex stress functions. In this study, boundary collocation method (BCM) has been applied in the determination of the unknown coefficients, which is a simple and accurate numerical method, and has been well established by many researchers.^25–28 After the unknown coefficients have been obtained, the crack tip SIFs can be obtained according to the definition of mode I and mode II SIF.

Although the crack tip SIFs can also be obtained through finite element (FE; ABAQUS code) numerical simulation, the corresponding analytical solution is still necessary and significant because the FE simulation is approximate, and the results may have a large error and uncertainty.

Stress functions $\mathrm{\Phi }(z)$ and $\mathrm{\Omega }(z)$ for a tunnel with a radial crack

For the problem of a tunnel containing a radial crack as shown in Figure 1, we define a basic Plemelj function as

$$X(z)=\frac{1}{\sqrt{z(z-b)}}$$

(1)

where $b=R+a$ , and R is the radius of the tunnel roof and a is the length of the radial crack.

Figure 1.

Sketch of a tunnel with a radial crack under compression.

It can be seen that $z\to \infty ,\mathrm{zX}(z)\to 1$ . Suppose t₀ is a point of crack surfaces, and if X(z) has a unique limit as $z\to t_0$ along any path in the right neighborhood of t₀, we denote this limit by $X^{+}(t_0)$ as shown in Figure 2. Similarly, if a unique limit exists as approaching to t₀ from the left neighborhood, it is denoted by $X^{-}(t_0)$ . In order to investigate the relation between $X^{+}(t)$ and $X^{-}(t)$ as X(z) crosses the cracks, let us trace a path beginning at a point D on a crack surface and leading from the right side of the crack around the crack tip N to approach D from the left as shown in Figure 2. From equation (1), one can have

$$X^{-}(t)=X^{+}(t)·e^{1\times (-0.5)\times (-2\pi )i}=-X^{+}(t)$$

(2)

Figure 2.

Sketch of a simple closed path.

According to Zhu, ²² for a tunnel containing a radial crack under the far-field compressive stress ${\sigma }_V$ and ${\sigma }_H$ , the complex stress functions can be written as

$$[\mathrm{\Phi }(t)+\mathrm{\Omega }(t){]}^{+}+[\mathrm{\Phi }(t)+\mathrm{\Omega }(t){]}^{-}=-2{\sigma }_N+2i{\sigma }_f=Q$$

(3)

$$[\mathrm{\Phi }(t)-\mathrm{\Omega }(t){]}^{+}-[\mathrm{\Phi }(t)-\mathrm{\Omega }(t){]}^{-}=0$$

(4)

where Q is a constant, Φ(t) and Ω(t) are two stress complex functions, and ${\sigma }_N$ and ${\sigma }_f$ are the normal stress acting on the crack surfaces and crack surface friction, respectively, and they can be expressed as

$${\sigma }_N=\frac{{\sigma }_V+{\sigma }_H}{2}+\frac{{\sigma }_V-{\sigma }_H}{2}\cos 2\alpha$$

(5)

$${\sigma }_f={\sigma }_{N}f=(\frac{{\sigma }_V+{\sigma }_H}{2}+\frac{{\sigma }_V-{\sigma }_H}{2}\cos 2\alpha )f$$

(6)

It should be noted that if the shear stress ${\sigma }_S$ ( ${\sigma }_S=({\sigma }_V-{\sigma }_H/2)\sin 2\alpha$ , $\alpha$ is shown in Figure 1), which is parallel to the crack surface, is less than the crack surface friction, that is, ${\sigma }_S<{\sigma }_{N}f$ , then the crack tip SIF is zero due to large crack surface friction. Therefore, only when ${\sigma }_S\ge {\sigma }_{N}f$ , the crack tip SIF exists.

According to the stress single-valued condition in multi-connected zones in an infinite plane, the stress complex functions Φ(z) and Ω(z) can be written as

$$\mathrm{\Phi }(z)=B-\frac{1+\nu }{8\pi }(X+\mathrm{iY})\frac{1}{z}+\frac{a_1}{z^2}+\frac{a_2}{z^3}+···$$

(7)

$$\mathrm{\Omega }(z)=B+C-\mathrm{iD}+\frac{3-\nu }{8\pi }(X+\mathrm{iY})\frac{1}{z}+\frac{b_1}{z^2}+\frac{b_2}{z^3}+···$$

(8)

where ν is Poisson’s ratio; B, C, and D are constants, and they can be determined by the far-field compressive stresses, that is

$$\begin{array}{c}B=\frac{1}{4}({\sigma }_V+{\sigma }_H)\hfill \\ \mathrm{C}=-\frac{{\sigma }_V-{\sigma }_H}{2}\cos 2(\frac{\pi }{2}-\alpha )=\frac{{\sigma }_V-{\sigma }_H}{2}\cos 2\alpha \hfill \\ \mathrm{D}={\tau }_{\mathrm{xy}}=\frac{{\sigma }_V-{\sigma }_H}{2}\sin 2(\frac{\pi }{2}-\alpha )=\frac{{\sigma }_V-{\sigma }_H}{2}\sin 2\alpha \hfill \end{array}$$

From equations (7) and (8), one can have

$$\underset{z\to \infty }{lim}\mathrm{\Phi }(z)=B\mathrm{and}\underset{z\to \infty }{lim}\mathrm{\Omega }(z)=B+C-\mathrm{iD}$$

(9)

Let $E(t)=\mathrm{\Phi }(t)+\mathrm{\Omega }(t)$ and $F(t)=\mathrm{\Phi }(t)-\mathrm{\Omega }(t)$ , equations (3) and (4) can be rewritten as

$$E^{+}(t)+E^{-}(t)=Q$$

(10)

$$F^{+}(t)-F^{-}(t)=0$$

(11)

Equations (10) and (11) are simple Hilbert problem, and dividing by $X^{+}(t)$ and $X^{-}(t)$ , from equations (10) and (2), one can have

$$\frac{E^{+}(t)}{X^{+}(t)}-\frac{E^{-}(t)}{X^{-}(t)}=\frac{Q}{X^{+}(t)}$$

(12)

Let $G(t)=E(t)/X(t)$ and $g(t)=Q/X^{+}(t)$ , then equation (12) can be rewritten as

$$G^{+}(t)-G^{-}(t)=g(t)$$

(13)

According to the theory of Cauchy integral on an arc, the general solution of equation (13) can be written as

$$G(z)=\frac{1}{2\pi i}\underset{L}{\int }\frac{g(t)\mathrm{dt}}{t-z}+P(z)$$

(14)

where $P(z)=\sum _{k=-\infty }^1{C}_k{z}^k=C_{-\infty }{z}^{-\infty }+\cdots +C_{-1}{z}^{-1}+C_0{z}^0+C_{1}z$ .

From equation (14), one can have

$$\mathrm{\Phi }(z)+\mathrm{\Omega }(z)=\frac{X(z)}{2\pi i}\underset{L}{\int }\frac{\mathrm{Qdt}}{X^{+}(t)(t-z)}+X(z)P(z)$$

(15)

Equation (11) can be rewritten as

$$F^{+}(t)=F^{-}(t)$$

(16)

This means that the function F(z) is analytic in the entire plane except possibly at infinite where it can at most have a pole. Hence, by Laurent’s theorem, F(z) is a polynomial, and considering the condition that $z\to \infty ,F(z)$ is limited, then the general solution is

$$F(z)=\mathrm{\Phi }(z)-\mathrm{\Omega }(z)=\mathrm{\Phi }(\infty )-\mathrm{\Omega }(\infty )=-C+\mathrm{iD}$$

(17)

From equations (15) and (17), the two complex stress functions for the case of tunnel with radial crack under compression can be formulated as

$$\mathrm{\Phi }(z)=\frac{X(z)}{4\pi i}\underset{L}{\int }\frac{\mathrm{Qdt}}{X^{+}(t)(t-z)}+X(z)P(z)-\frac{1}{2}(C-\mathrm{iD})$$

(18)

$$\mathrm{\Omega }(z)=\frac{X(z)}{4\pi i}\underset{L}{\int }\frac{\mathrm{Qdt}}{X^{+}(t)(t-z)}+X(z)P(z)+\frac{1}{2}(C-\mathrm{iD})$$

(19a)

where

$$C-\mathrm{iD}=\frac{1}{2}({\sigma }_V-{\sigma }_H)e^{-2i\alpha }$$

(19b)

where $C_k(-\infty <k\le 0)$ is constant, and they will be determined by the BCM later, and the constant $C_1$ can be obtained from equation (9)

$$C_1=B+\frac{1}{2}(C-\mathrm{iD})=-\frac{1}{4}({\sigma }_V+{\sigma }_H)-\frac{1}{4}({\sigma }_V-{\sigma }_H)e^{-2i\alpha }$$

(19c)

In order to get the solution of $\mathrm{\Phi }(z)$ and $\mathrm{\Omega }(z)$ from equations (18) and (19a), first we have to solve the following equation

$$I(z)=\underset{L}{\int }\frac{\mathrm{Qdt}}{X^{+}(t)(t-z)}$$

(20)

For the radial crack, the integrals along L_k (k = 1) in equation (20) can be expressed in terms of integral along a lacet C_k (k = 1) surrounding it as shown in Figure 3, and the integral over the lacet can be evaluated by the residue theory. Considering the contour integral around the lacet C_k (k = 1), one can have

$$\begin{array}{c}\underset{C_k}{\int }\frac{f(\zeta )d\zeta }{X(\zeta )(\zeta -z)}=\underset{L_k}{\int }\frac{f(t)\mathrm{dt}}{X^{+}(t)(t-z)}+\underset{\rho \to 0}{lim}\underset{|z-b_k|=\rho }{\int }\frac{f(\zeta )d\zeta }{X(\zeta )(\zeta -z)}\hfill \\ -\underset{L_k}{\int }\frac{f(t)\mathrm{dt}}{X^{-}(t)(t-z)}+\underset{\rho \to 0}{lim}\underset{|z-a_k|=\rho }{\int }\frac{f(\zeta )d\zeta }{X(\zeta )(\zeta -z)}\hfill \end{array}$$

(21)

Figure 3.

An arc and a lacet surrounding it.

It can be seen that the second and the fourth integrals in equation (21) tend toward zero as ρ→0, and combining equation (2), one can have

$$\begin{array}{c}\multicolumn{1}{c}{\underset{C_k}{\int }\frac{f(\zeta )d\zeta }{X(\zeta )(\zeta -z)}=\underset{L_k}{\int }\frac{f(t)}{t-z}[\frac{1}{X^{+}(t)}-\frac{1}{X^{-}(t)}]\mathrm{dt}}\\ =2\underset{L_k}{\int }\frac{f(t)\mathrm{dt}}{X^{+}(t)(t-z)}\hfill \end{array}$$

(22)

From equation (22), the relationship between the integrals along a crack L_k and along its lacet C_k can be obtained

$$I(z)=\underset{L}{\int }\frac{f(t)\mathrm{dt}}{X^{+}(t)(t-z)}=\frac{1}{2}\underset{C}{\int }\frac{f(\zeta )d\zeta }{X(\zeta )(\zeta -z)}$$

(23)

Using equation (23), equations (18) and (19a) can be rewritten as

$$\mathrm{\Phi }(z)=\frac{X(z)}{8\pi i}\underset{C}{\int }\frac{\mathrm{Qd}\zeta }{X(\zeta )(\zeta -z)}+X(z)P(z)-\frac{1}{2}(C-\mathrm{iD})$$

(24)

$$\mathrm{\Omega }(z)=\frac{X(z)}{8\pi i}\underset{C}{\int }\frac{\mathrm{Qd}\zeta }{X(\zeta )(\zeta -z)}+X(z)P(z)+\frac{1}{2}(C-\mathrm{iD})$$

(25)

From the residue theory, the integral in equation (24) along the contour C can be expressed in terms of the sum of the residues at point z and at infinity; thus, equation (24) can be rewritten as

$$\begin{array}{c}\multicolumn{1}{c}{\mathrm{\Phi }(z)=\frac{\mathrm{QX}(z)}{4}[\frac{1}{X(z)}+\mathrm{Re}s(\frac{1}{X(t)(t-z)},\infty )]}\\ +X(z)P(z)-\frac{1}{2}(C-\mathrm{iD})\hfill \end{array}$$

(26)

where $\mathrm{Re}s(1/X(t)(t-z),\infty )=-\mathrm{Re}s(1/X(1/t)(1/t-z),0)=-z+(1/2)b$ .

By simplifying, equation (26) can be expressed as

$$\mathrm{\Phi }(z)=\frac{Q}{4}+\frac{\mathrm{QX}(z)}{4}[-z+\frac{1}{2}b]+X(z)P(z)-\frac{1}{2}(C-\mathrm{iD})$$

(27)

Equation (27) can be written in a simple form

$$\begin{array}{c}\multicolumn{1}{c}{\mathrm{\Phi }(z)=\frac{1}{\sqrt{z(z-b)}}\{\sum _{k=-\infty }^1{C}_k{z}^k-\frac{Q}{4}z+\frac{Q}{8}b\}}\\ +\frac{Q}{4}-\frac{1}{2}(C-\mathrm{iD})\hfill \end{array}$$

(28)

Similarly, the function $\mathrm{\Omega }(z)$ can be obtained, and it can be expressed as

$$\begin{array}{c}\multicolumn{1}{c}{\mathrm{\Omega }(z)=\frac{1}{\sqrt{z(z-b)}}\{\sum _{k=-\infty }^1{C}_k{z}^k-\frac{Q}{4}z+\frac{Q}{8}b\}}\\ +\frac{Q}{4}+\frac{1}{2}(C-\mathrm{iD})\hfill \end{array}$$

(29)

where $C_k(-\infty <k\le 0)$ are unknown coefficients. If the coefficient $C_k$ is solved, the stress functions, $\mathrm{\Phi }(z)$ and $\mathrm{\Omega }(z)$ , can be obtained from equations (28) and (29), respectively. In order to acquire the unknown coefficient $C_k$ , the stress boundary conditions will be applied.

Stress boundary condition

Generally, stress functions can be written as

$$U(z)=\frac{1}{2}[\overline{z}\varphi (z)+z\overline{\varphi (z)}+\theta (z)+\overline{\theta (z)}]$$

(30)

where $\phi \prime (z)=\mathrm{\Phi }(z)$ and $\theta ″(z)=\overline{\mathrm{\Omega }(\overline{z})}-\mathrm{\Phi }(z)-z\mathrm{\Phi }\prime (z)$ . According to the theory of elasticity, the relation of the stresses ${\sigma }_x$ , ${\sigma }_y$ , and ${\tau }_{\mathrm{xy}}$ versus the stress function $U(z)$ can be expressed as

$${\sigma }_x=\frac{{\partial }^{2}U}{\partial y^2};{\sigma }_y=\frac{{\partial }^{2}U}{\partial x^2};{\tau }_{\mathrm{xy}}=-\frac{{\partial }^{2}U}{\partial x\partial y}$$

(31)

where $(\partial /\partial x)=(\partial /\partial z)+(\partial /\partial z)$ and $\partial /\partial y=i(\partial /\partial z)-i(\partial /\partial \overline{z})$ . Substituting equation (30) into equation (31), one can obtain the stress^4,22,24 ${\sigma }_x$ , ${\sigma }_y$ , and ${\tau }_{\mathrm{xy}}$ as

$${\sigma }_x=3\mathrm{Re}[\mathrm{\Phi }(z)]-\mathrm{Re}[\mathrm{\Omega }(\overline{z})]-(\overline{z}-z)i\mathrm{Im}[\mathrm{\Phi }\prime (z)]$$

(32)

$${\sigma }_y=\mathrm{Re}[\mathrm{\Phi }(z)]+\mathrm{Re}[\mathrm{\Omega }(\overline{z})]+(\overline{z}-z)i\mathrm{Im}[\mathrm{\Phi }\prime (z)]$$

(33)

$${\tau }_{\mathrm{xy}}=-\mathrm{Im}[\mathrm{\Phi }(z)]-\mathrm{Im}[\mathrm{\Omega }(\overline{z})]-i(\overline{z}-z)\mathrm{Re}[\mathrm{\Phi }\prime (z)]$$

(34)

where $\mathrm{\Phi }(z)$ and $\mathrm{\Omega }(z)$ can be acquired from equations (28) and (29), respectively.

According to Figure 4, the stress boundary condition can be expressed as

$$l{\sigma }_x+m{\tau }_{\mathrm{xy}}=\overline{X}$$

(35)

$$l{\tau }_{\mathrm{xy}}+m{\sigma }_y=\overline{Y}$$

(36)

where $l=\cos (N,x)=\cos \gamma$ , $m=\cos (N,y)=\sin \gamma$ , and N is normal vector of the point P. Substituting equations (32), (33), and (34) into equations (35) and (36), and combining with equations (28) and (29), the boundary conditions $\overline{X}$ and $\overline{Y}$ can be expressed in terms of the coefficients C_k and complex variable z. The parameters of l and m on different tunnel edges in Figure 1 are presented in Table 1.

Figure 4.

Sketch of the relation of stresses versus the forces on a boundary.

Table 1.

The parameters of l and m on different tunnel edges according to Figure 1.

Position	l	m
Circle part	$\sin \theta$	$\cos \theta$
Left wall	$\cos \alpha$	$-\sin \alpha$
bottom	$\sin \alpha$	$\cos \alpha$
Right wall	$-\cos \alpha$	$\sin \alpha$

Determining the coefficients $C_k$ by BCM

In order to determine the unknown coefficients $C_k$ involved in equations (28) and (29), one can place a certain number of points along the tunnel edges as shown in Figure 1. The more points, the more accurate the calculation will be. However, when the number of points exceeds a certain number, the influence of the number on the calculation results diminishes. In this study, the total number of points placed is 1560, and the numbers of the points on the tunnel top, bottom, and two walls are 360, 400, 400, and 400, respectively.

Substituting the stress functions $\mathrm{\Phi }(z)$ and $\mathrm{\Omega }(z)$ in equations (28) and (29) into equations (32), (33), and (34), the stress components ${\sigma }_x$ , ${\sigma }_y$ , and ${\tau }_{\mathrm{xy}}$ can be expressed in terms of the coefficients $C_k$ and complex variable z. Substituting equations (32), (33), and (34) into equations (35) and (36), the boundary conditions $\overline{X}$ and $\overline{Y}$ can be expressed in terms of the coefficients $C_k$ and complex variable z.

For each collocated point on the tunnel edge, one can obtain the coordinate z and the direction cosine l and m, as well as the boundary stress conditions $\overline{X}$ and $\overline{Y}$ . If the series coefficients $C_k$ are truncated to a limited items, then equations (35) and (36) will become the equations with respect to the unknowns $C_k$ . Because there are 1560 points placed along the tunnel edge, one can yield 3120 (2 × 1560) such equations, and finally by means of the least square method, the coefficients $C_k$ can be calculated.

In this study, the number of summation terms k involved in equations (28) and (29) is truncated to 17, that is, k = –15, –14, …, –1, 0, 1. Based on the above procedure and parameters, the 17 coefficients $C_k$ can be solved, and therefore, the stress functions $\mathrm{\Phi }(z)$ and $\mathrm{\Omega }(z)$ for tunnel with radial crack can be obtained completely.

It should be noted that one can use resultant stress boundary conditions to determine the unknown coefficients C_k, which has been applied in the BCM by many researchers.^10,21,25,28 However, using general stress boundary condition is simple and easy, whereas using resultant boundary condition sometimes is easily confused with the real and imaginary parts.

SIFs for a tunnel containing a radial crack

The SIFs $K_I$ and $K_{\mathrm{II}}$ can be expressed in terms of stress function^4,22,24 $\mathrm{\Phi }(z)$ as

$$K_I-i{K}_{\mathrm{II}}=\underset{z\to a}{lim}2\sqrt{2\pi (z-b)}\mathrm{\Phi }\left(z\right)$$

(37)

Substituting equation (28) into equation (37), the SIFs can be rewritten as

$$K_I-i{K}_{\mathrm{II}}=\frac{2\sqrt{2\pi }}{\sqrt{b}}\{\sum _{k=-\infty }^1{C}_k{b}^k-0.125\mathrm{Qb}\}$$

(38)

The dimensionless SIFs can be expressed as

$$Y_I-i{Y}_{\mathrm{II}}=\frac{K_I-i{K}_{\mathrm{II}}}{K_0}=\frac{2\sqrt{2}}{\sigma }\{\sum _{k=-\infty }^1{C}_k{b}^{k-1}-0.125Q\}$$

(39)

where $K_0=\sigma \sqrt{\pi b}$ , and $\sigma =max\{{\sigma }_H,{\sigma }_V\}$ , $Y_I$ , and $Y_{\mathrm{II}}$ are the dimensionless $K_I$ and $K_{\mathrm{II}}$ , respectively. $C_k$ and Q are complex, and they can be expressed in terms of real and imaginary parts as

$$C_k=R_k+i{I}_k\mathrm{and}Q=\mathrm{QR}+\mathrm{iQI}$$

(40)

Substituting equation (40) into equation (39), one can have

$$Y_I=\frac{2\sqrt{2}}{\sigma }\{\sum _{k=-\infty }^1{R}_k{b}^{k-1}-0.125\mathrm{QR}\}$$

(41)

$$Y_{\mathrm{II}}=-\frac{2\sqrt{2}}{\sigma }\{\sum _{k=-\infty }^1{I}_k{b}^{k-1}-0.125\mathrm{QI}\}$$

(42)

Equations (41) and (42) are the formulas of the dimensionless SIFs for tunnel containing a radial crack under compression, and the SIFs $K_I$ and $K_{\mathrm{II}}$ can be obtained from equation (38).

The calculation results show that the mode I SIF K_I is negative, which means that the cracks under compression will close, and the corresponding crack tips do not have stress concentration. Actually, the compressive load will cause compressive stress at crack tips, which will have a negative effect on crack propagation, but the effect extent of the negative K_I is not clear yet. ²² Therefore, in this study, the negative K_I is not considered, and only mode II K_II is calculated.

Figure 5 shows the analytical results obtained by equation (42) for an infinite plane containing a tunnel with a radial crack under biaxial compression with different ratios of ${\sigma }_H/{\sigma }_V$ . In the calculating, the tunnel’s height H and width W are 5 and 3.5 m, respectively, and the length of the radial crack equals to W/2, that is, a = 1.75 m, R = 1.75 m, and b = a + r = 3.5 m. It can be seen that the dimensionless SIF Y_II varies with crack inclination angle α, and varies with the ratio of ${\sigma }_H/{\sigma }_V$ . As ${\sigma }_H/{\sigma }_V=0.0$ , the corresponding Y_II–α curve is the upmost, and as ${\sigma }_H/{\sigma }_V$ increases from 0.0 to 1.0, the corresponding Y_II–α curve tends toward flat, and the corresponding Y_II decreases; and as ${\sigma }_H/{\sigma }_V$ continues to increase from 1.0 to ∝, the corresponding Y_II–α curves change from convex curves to concave curves. Here, it should be noted that sign of Y_II values only express crack propagation direction, and the positive sign expresses that the wing crack initiated at a mode II crack tip in Figure 1 will curve to left side, and the negative sign expresses that the wing crack will curve to right side.

Figure 5.

Analytical results of dimensionless SIF Y_II versus crack inclination angle α from equation (42) for an infinite plane containing a tunnel with a radial crack under biaxial compression.

Validation of the theoretical results

In this chapter, a numerical simulation using FE code ABAQUS is implemented, and the results are compared with the theoretical results. In addition, experimental study is also implemented to validate the theoretical solutions.

Validation through FE simulation by ABAQUS code

In order to accurately simulate a tunnel with a crack, the dimensions of the numerical model should be large enough because the theoretical results were obtained based on an infinite plane containing a tunnel. In this study, a numerical model of a square plate measuring 120 m in length and 120 m in width is established, as shown in Figure 6. The tunnel’s height H and width W are 5 and 3.5 m, respectively, and the length a of the radial crack equals to W/2, that is, a = 1.75 m. In this FE simulation, the vertical stress ${\sigma }_V$ is kept as a constant 10 MPa, and ${\sigma }_H$ varies which are 0, 2, 4, and 6 MPa, respectively, and accordingly, the corresponding ratios of horizontal stress to vertical stress ${\sigma }_V/{\sigma }_H$ are 0, 0.2, 0.4, and 0.6, respectively. Crack surface friction coefficient f is 0.0, and crack inclination angle $\alpha$ varies in the range from 0° to 90°.

Figure 6.

A finite element model used in ABAQUS code.

Figure 7 shows the analytical solution obtained from equation (42) and from the FE simulation results using ABAQUS code. The dimensionless Y_II by ABAQUS code basically agrees with those analytical results from equation (42), but there are still some slight errors between them in Figure 7. The reason is that the theoretical results of SIFs is based on an infinite rock plane, whereas the FE simulation results by ABAQUS code is based on a finite square plate with a centralized “tunnel.”

Figure 7.

Comparison of the dimensionless Y_II between the analytic results from equation (42) and the simulation results from ABAQUS code: (a) σ_H/σ_V = 0.00, (b) σ_H/σ_V = 0.20, (c) σ_H/σ_V = 0.40, and (d) σ_H/σ_V = 0.60.

Tests of tunnel models with a radial crack

In order to validate the theoretical results, experimental study using tunnel models is implemented. The tunnel model is a plate specimen with a centralized horseshoe-shaped “tunnel” with a single radial crack, as shown in Figure 8, and the plate is 350 mm in length and 350 mm in width, and 100 mm in thickness, and the “tunnel” is 60 mm in height and 50 mm in width. The tunnel’s arch consists of a semi-circle, and the radius of the circular arch is 25 mm. The crack inclination angle $\alpha$ in tunnel models as shown in Figure 8 was designed as 0°, 15°, 30°, 45°, 60°, 75°, and 90°, and for each inclination angle, three specimens were made by cement mortar, and the ratio of cement: sand: gypsum: water is 1:2.5:0.6:0.45 by weight. The radial crack length is 50 mm, and it was made using a very thin plastic film, which was placed inside the tunnel models during the process of casting in a mold until they were loaded.

Figure 8.

A tunnel mode and the loading device.

The tunnel models were loaded vertically by electrohydraulic servo control device (500T), and the lateral deformation was confined by the hydraulic jack and the steel plates as shown in Figure 8. Based on the lateral deformation is zero, the relationship between the horizontal stress σ₃ and the vertical stress σ₁, according to the elastic theory, is σ₃ = vσ₁, where ν is Poisson’s ratio.

In order to decrease the influence of uneven surface, the specimen surfaces were polished using grinder before loading. The specimen surfaces were also smeared with a kind of lubricant, which can avoid the effect of the friction between the specimen and the loading device. These measures can ensure a better loading condition and improve the accuracy of testing results.

The fracture patterns of the tunnel models containing a crack emanating from the tunnel roof with different inclination angles α under biaxial compressions are shown in Figure 9. It can be seen that crack inclination angles affect tunnel model failure largely.

Figure 9.

Fracture patterns of tunnel models with a radial crack under compression: (a) α = 0°, (b) α = 30°, (c) α = 45°, (d) α = 60°, (e) α = 75°, and (f) α = 90°.

In order to compare the critical compressive stresses predicted by the analytical solution equation (42) with those testing results, based on the relation of mode II SIF K_II versus the major compressive stress σ₁, that is, $K_{\mathrm{II}}=Y_{\mathrm{II}}{\sigma }_1\sqrt{\pi b}$ , the critical compressive stress σ_C can be expressed in terms of mode II fracture toughness K_IIC as

$${\sigma }_C=\frac{K_{\mathrm{IIC}}}{Y_{\mathrm{II}}\sqrt{\pi b}}$$

(43)

where the fracture toughness K_IIC of the cement mortar material measured is $2.8\mathrm{MPa}\sqrt{m}$ , and the dimensionless SIF Y_II is presented in Figure 5 under different ratios of ${\sigma }_H/{\sigma }_V$ .

Using equations (43) and (42), the critical compressive stresses for the tunnel models with an inclination crack can be predicted, and the results are compared with the testing results in Figure 10. It can be seen that basically, the critical compressive stresses predicted by equations (42) and (43) agree with the tunnel model testing results, and as the crack inclination angle is 45°, the corresponding critical compressive stress is the lowest under the conditions that ${\sigma }_V>{\sigma }_H$ .

Figure 10.

Comparison of the tunnel model testing results with the prediction results by equations (42) and (43).

Conclusion

For a tunnel with a radial crack emanating from the tunnel edge, the complex stress functions were obtained according to the crack surface stress condition. BCM was applied in the determination of the coefficients of the stress functions, and the SIFs can be determined using the stress functions. In order to validate the theoretical result, numerical simulation using ABAQUS code and tunnel model compressive tests were implemented, and from this study, the following conclusion can be obtained:

1. For a tunnel with a radial crack emanating from the tunnel roof, the stress functions can be obtained based on crack surface boundary conditions, and the coefficients of the stress functions can be effectively determined by BCM.

2. The theoretical result agrees well with the numerical simulation results and the tunnel model test results, which indicates that the method determining the stress functions and the SIFs for a tunnel with a radial crack proposed in this article is effective and reliable.

3. Figure 5 shows that crack inclination angles and the ratios of horizontal to vertical stress affect the SIFs significantly.

4. For a tunnel with a 45° inclination crack emanating from the tunnel roof, the tunnel model compressive strength is low, and the corresponding tunnel stability is low.