# Complex Analytical Study of the Stability of Tunnel-Surrounding Rock in a Layered Jointed Rock Mass

Authors: Su Wu; Kegang Li; Zaihong Yang; Shunchuan Wu; Ting Wang; Qingci Qin; Mingliang Li
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

A barrier function method based on the fmincon optimization function in MATLAB was used to determine the function to map a tunnel boundary to the unit circle in the complex plane, and a structural failure criterion of the mapping convergence was established based on the reliability theory of underground engineering. The joints were approximated as cracks around the tunnel, the anisotropy of the stress intensity factor due to the crack inclination and position changes was studied, and the modified scoring parameter of the layered joints around the tunnel for the rock geomechanics classification (RMR) was proposed and used at points around the tunnel. The calculation method for the supporting stress was required to balance the bias load. The results showed the following: (1) When taking the structural failure criterion as the convergence condition of the mapping function, the mapping function converged when the mapping accuracy wasδ≥31.74%. (2) The crack with an inclination angle of β=45° was the dominant structural plane of the jointed rock mass around the tunnel. The anisotropy of the stress intensity factor KII at the tip of the mode II crack indicated that the corresponding cracks at the various points around the tunnel had inconsistent influences on the tunnel. The position had the greatest influence, followed by the straight wall area, and then by the top and floor areas. (3) When the crack inclination β was equal to the inclination angle β0 of the dominant joint plane, the secondary crack was parallel to the unloading surface of the corresponding tunnel. (4) The bias load formed by the layered joints around the tunnel reduced the stress threshold of the failure of the rib spalling.


---

## Body

## 1. Introduction

A layered jointed rock mass is a sedimentary rock with a layered structure composed of an array of parallel joint planes. The anisotropy of the mechanical behavior of the rock strata and the low-strength characteristics of the joint planes make the stability of the surrounding rock complicated for underground engineering [1]. Many results have been obtained using numerical simulation and model test methods to study the mechanical characteristics of the surrounding rock with different rock inclination angles. Wu et al. [2], Sha et al. [3], Chen et al. [4], and Tian et al. [5] used various numerical simulations to establish a tunnel model with layered joints, and they analyzed the tunnel-surrounding rock with different joint dips, structural deformation, force characteristics, and failure modes. Ma et al. [6] used a large-scale three-dimensional model similarity test system to analyze the deformation and failure behaviors of multilayer jointed rock masses with different inclination angles for high ground stresses. Zhou et al. [7, 8], Li and Zhu [9], and Liu et al. [10] approximated joints as cracks around tunnels, and they used different finite element analysis software programs to simulate the stability of the tunnel-surrounding cracks to the surrounding rock of straight wall arch tunnels. These simulations were combined with model tests to study the damage and destruction behaviors under confining pressure. The results showed that the cracks at the shoulder and foot of the tunnel were weaker locations of the tunnel-surrounding rock damage. When the crack inclination angle was β=45°, the cracks had the greatest impact on the overall stability and strength, and the tunnel shear failure was more significant.Although numerical simulations and model testing are very effective methods, the correctness of the analysis results depends on the rationality of the selected models and materials, which are affected by the subjective judgement of researchers. These individual results that do not reveal the full underlying physics. Using conformal mapping theory of complex analysis can provide a method to study the surrounding rock stress in underground engineering to obtain a closed solution via mathematical mechanics, i.e., using a complex analytical method. Huangfu et al. [11], Zhu et al. [12], and Li and Liu [13] used different algorithms to obtain the conformal mappings of complex cross sections to the unit circle in the complex plane. Huangfu et al. [14], Li et al. [15], Zhu et al. [16], and Li and Chen [17] determined the analytical solutions of the surrounding rock stresses of underground tunnels and chambers by obtaining the conformal mappings. The distribution characteristics of the stress in the surrounding rock of a tunnel were studied. However, the complex analysis method has seldom been used to study the mechanical behaviors of tunnel-surrounding rock in a layered jointed rock mass. Additionally, the convergence condition of the existing mapping function cannot describe the deviation degree of the mapping value of each sampling point from the actual value, and the impact of the random uncertainty of a sampling point on the accuracy has not been considered.Based on the results of previous research and their deficiencies, in this study, the barrier function method based on the fmincon optimization function was used to establish the conformal mapping from a tunnel boundary to the unit circle and develop a structural failure criterion for the convergence of the mapping function. Based on this, the joint was approximated as a crack around the tunnel, and the complex analytical expression of the stress intensity factor at the crack tip under confining pressure was derived. The anisotropy of the stress intensity factor due to crack inclination and position changes was studied. The internal relationship between the cracks at different locations around the tunnel and the corresponding inclination and initiation angles of the dominant joint plane were revealed. The instability propagation criterion of the crack tip of the dominant joint plane at any location around the tunnel was deduced for the maximum circumferential stress. By combining the fragmentation mechanism of the fractured rock mass and considering the relationship between the secondary cracks and the unloading surface of the tunnel excavation, a modified scoring parameter for the inclination of the layered joints around the tunnel was proposed for the Rock Mass Rating (RMR). Based on the complex variable function form of the stability judgement of the layered joint plane at any point around the tunnel, a method was proposed to calculate the supporting stress required to balance the bias load at any point.


## 2. Analysis of Complex Change of the Surrounding Rock Stress in Deep Underground Tunnel

### 2.1. Conformal Mapping from Tunnel Boundary to Unit Circle

When using functions of complex variables to study elastic mechanics plane problems, a given area in the physical plane can be transformed into another area with simple boundary shape in the image plane [18] to study the mechanical behavior of complex boundary engineering. As shown in Figure 1, the counterclockwise direction of the tunnel boundary was specified as the positive direction, and the outer domain of the tunnel boundary was conformally mapped to the outer domain of the unit circle in the complex plane through the mapping function z=ωζ [19]. The mapping function can be expanded in the form of a Laurent series [20] as follows:(1)z=ωζ=Rζ+∑k=0∞Ckζ−k,where ζ is the coordinate of any point outside the boundary of the mapping in the complex plane, R is a positive real number that comprehensively reflects the shape and size of the tunnel section, and Ckk=0,1,2,…,n… represents the kth constant term of the Laurent series, which depends the shape of the tunnel section. These are real constants when discussing axisymmetric problems.Figure 1
Conformal mapping from the exterior region of the tunnel to the exterior region of the unit circle.


### 2.2. Complex Variable Function Solution for Surrounding Rock Stress at the Tunnel Boundary

According to Muskhelishvili’s complex function method [21], the stress analysis of the surrounding rock of a deep underground tunnel can be classified as the problem of excavation in an infinite plane. The stress component around the tunnel after conformal mapping can be expressed in the complex plane as(2)σρ=124ReΦζ−Re2ζ2¯¯ω′ζωζ¯Φ′ζ+ω′ζΨζ,σθ=124ReΦζ+Re2ζ2¯¯ω′ζωζ¯Φ′ζ+ω′ζΨζ,τρθ==2ζ2¯¯ω′ζωζ¯Φ′ζ+ω′ζΨζ,where Φζ=φζ´/ωζ´; Ψζ=ψζ´/ωζ´; ζ is a multivalued function defined as ζ=eiθ=cosθ+isinθ; ωζ¯ is the conjugate complex of ωζ; σρ, σθ, and τρθ are the stress components; and i is the imaginary unit. Re represents the real part of the complex number, and Im represents the imaginary part of the complex number.The analytical functionsφζ and ψζ in equation (2) are defined as follows:(3)φζ=σh∞+σv∞4ωζ+φ0ζ,ψζ=σh∞+σv∞4ωζ+φ0ζ,where σh∞ and σv∞ are the horizontal and vertical stresses at infinity, respectively, and φ0ζ and ψ0ζ are the analytic functions in the neighborhood of infinity, which can be determined using the power series expansion method [15] or the Cauchy integral method [22].The stress boundary condition is(4)φζ+ωζ¯¯ω′ζφ′ζ¯+ψζ¯=i∫abFx+iFyds,where a is the integral starting point taken on the boundary, b is any point on the boundary, and Fx and Fy are the given surface force components along the x-axis and y-axis on the boundary, respectively. When a tunnel is excavated without support and there is no external load around the tunnel, the surface force components Fx and Fy are both zero, and the integral term on the right-hand side of the equation is zero.


### 2.3. Stability Discriminant of Layered Joint Plane around Tunnel

According to the Mohr–Coulomb criterion, when the joint plane around the deep underground tunnel is in a limit equilibrium state, the shear stress on the joint plane should be less than the shear strength. The discriminant of the stability of the joint plane is as follows [23]:(5)σ1cosβsinφj−β+σ3sinβcosφj−β+cjcosφj≥0,where σ1 and σ3 represent the large and small principal stresses around the tunnel, β is the dip angle of the joint, cj is the cohesion of the jointed rock mass, and φj is the internal friction angle of the jointed rock mass. When equation (5) is equal to zero, it represents the limiting equilibrium state. If the left-hand side of the equation is less than zero, the joint plane is in an unstable state.Equation (5) can be used to judge the stability of the boundary of an underground tunnel in a jointed rock mass. The excavated tunnel is radially unloaded, and the tangential stress is concentrated, i.e.,(6)σ1=σα,σ3=σr=0,where σα and σr are the tangential and radial stresses at any point of the tunnel boundary, respectively, and α is the polar angle at any point of the tunnel boundary. Equation (5) can thus be written as(7)σαcosβsinφj−β+cjcosφj≥0,where σα can be obtained from the polar stress component σθ in the complex plane using equation (2), and equation (7) can be substituted to obtain(8)2ReΦζ+Reζ2¯¯ω′ζωζ¯Φ′ζ+ω′ζΨζ≥cjcosφjcosβsinβ−φj,where each item on the left-hand side is a function of ζ. Equation (8) can be expressed as follows:(9)Fζ≥cjcosφjcosβsinβ−φj.The coordinate of any point on the boundary of the tunnel is expressed asAjrj,αj, and the coordinate of the mapping point on the corresponding unit circle is Bj1,θj, as shown in Figure 1. The mapping function given by equation (1) can be regarded as the superposition of several fractional linear mappings. According to the shape-preserving property of the fractional linear mapping, the inverse mapping relationship z−1rj,αj=ω−1ζ1,θj can be obtained as follows:(10)rj=Rcosαj−θj+∑k=0∞Ckcosαj−kθj,αj=1+∑k=0∞Ckθj1−∑k=1∞kCk.The stability of the jointed rock mass around the tunnel is not only related to the occurrence of the joint itself, namely, the dip angleβ and the rock mechanical parameters (cohesion cj and internal friction angle φj) but also to the position ζ1,θj and the position and stress magnitude Fζ of the joint endpoint at the tunnel boundary. Equations (8) and (10) together form a complex function for assessing the stability of the layered joint plane at any point around the tunnel.


## 2.1. Conformal Mapping from Tunnel Boundary to Unit Circle

When using functions of complex variables to study elastic mechanics plane problems, a given area in the physical plane can be transformed into another area with simple boundary shape in the image plane [18] to study the mechanical behavior of complex boundary engineering. As shown in Figure 1, the counterclockwise direction of the tunnel boundary was specified as the positive direction, and the outer domain of the tunnel boundary was conformally mapped to the outer domain of the unit circle in the complex plane through the mapping function z=ωζ [19]. The mapping function can be expanded in the form of a Laurent series [20] as follows:(1)z=ωζ=Rζ+∑k=0∞Ckζ−k,where ζ is the coordinate of any point outside the boundary of the mapping in the complex plane, R is a positive real number that comprehensively reflects the shape and size of the tunnel section, and Ckk=0,1,2,…,n… represents the kth constant term of the Laurent series, which depends the shape of the tunnel section. These are real constants when discussing axisymmetric problems.Figure 1
Conformal mapping from the exterior region of the tunnel to the exterior region of the unit circle.


## 2.2. Complex Variable Function Solution for Surrounding Rock Stress at the Tunnel Boundary

According to Muskhelishvili’s complex function method [21], the stress analysis of the surrounding rock of a deep underground tunnel can be classified as the problem of excavation in an infinite plane. The stress component around the tunnel after conformal mapping can be expressed in the complex plane as(2)σρ=124ReΦζ−Re2ζ2¯¯ω′ζωζ¯Φ′ζ+ω′ζΨζ,σθ=124ReΦζ+Re2ζ2¯¯ω′ζωζ¯Φ′ζ+ω′ζΨζ,τρθ==2ζ2¯¯ω′ζωζ¯Φ′ζ+ω′ζΨζ,where Φζ=φζ´/ωζ´; Ψζ=ψζ´/ωζ´; ζ is a multivalued function defined as ζ=eiθ=cosθ+isinθ; ωζ¯ is the conjugate complex of ωζ; σρ, σθ, and τρθ are the stress components; and i is the imaginary unit. Re represents the real part of the complex number, and Im represents the imaginary part of the complex number.The analytical functionsφζ and ψζ in equation (2) are defined as follows:(3)φζ=σh∞+σv∞4ωζ+φ0ζ,ψζ=σh∞+σv∞4ωζ+φ0ζ,where σh∞ and σv∞ are the horizontal and vertical stresses at infinity, respectively, and φ0ζ and ψ0ζ are the analytic functions in the neighborhood of infinity, which can be determined using the power series expansion method [15] or the Cauchy integral method [22].The stress boundary condition is(4)φζ+ωζ¯¯ω′ζφ′ζ¯+ψζ¯=i∫abFx+iFyds,where a is the integral starting point taken on the boundary, b is any point on the boundary, and Fx and Fy are the given surface force components along the x-axis and y-axis on the boundary, respectively. When a tunnel is excavated without support and there is no external load around the tunnel, the surface force components Fx and Fy are both zero, and the integral term on the right-hand side of the equation is zero.


## 2.3. Stability Discriminant of Layered Joint Plane around Tunnel

According to the Mohr–Coulomb criterion, when the joint plane around the deep underground tunnel is in a limit equilibrium state, the shear stress on the joint plane should be less than the shear strength. The discriminant of the stability of the joint plane is as follows [23]:(5)σ1cosβsinφj−β+σ3sinβcosφj−β+cjcosφj≥0,where σ1 and σ3 represent the large and small principal stresses around the tunnel, β is the dip angle of the joint, cj is the cohesion of the jointed rock mass, and φj is the internal friction angle of the jointed rock mass. When equation (5) is equal to zero, it represents the limiting equilibrium state. If the left-hand side of the equation is less than zero, the joint plane is in an unstable state.Equation (5) can be used to judge the stability of the boundary of an underground tunnel in a jointed rock mass. The excavated tunnel is radially unloaded, and the tangential stress is concentrated, i.e.,(6)σ1=σα,σ3=σr=0,where σα and σr are the tangential and radial stresses at any point of the tunnel boundary, respectively, and α is the polar angle at any point of the tunnel boundary. Equation (5) can thus be written as(7)σαcosβsinφj−β+cjcosφj≥0,where σα can be obtained from the polar stress component σθ in the complex plane using equation (2), and equation (7) can be substituted to obtain(8)2ReΦζ+Reζ2¯¯ω′ζωζ¯Φ′ζ+ω′ζΨζ≥cjcosφjcosβsinβ−φj,where each item on the left-hand side is a function of ζ. Equation (8) can be expressed as follows:(9)Fζ≥cjcosφjcosβsinβ−φj.The coordinate of any point on the boundary of the tunnel is expressed asAjrj,αj, and the coordinate of the mapping point on the corresponding unit circle is Bj1,θj, as shown in Figure 1. The mapping function given by equation (1) can be regarded as the superposition of several fractional linear mappings. According to the shape-preserving property of the fractional linear mapping, the inverse mapping relationship z−1rj,αj=ω−1ζ1,θj can be obtained as follows:(10)rj=Rcosαj−θj+∑k=0∞Ckcosαj−kθj,αj=1+∑k=0∞Ckθj1−∑k=1∞kCk.The stability of the jointed rock mass around the tunnel is not only related to the occurrence of the joint itself, namely, the dip angleβ and the rock mechanical parameters (cohesion cj and internal friction angle φj) but also to the position ζ1,θj and the position and stress magnitude Fζ of the joint endpoint at the tunnel boundary. Equations (8) and (10) together form a complex function for assessing the stability of the layered joint plane at any point around the tunnel.


## 3. Solution for Tunnel Boundary Mapping Function

### 3.1. Constrained Optimization Model of Mapping Function

Lettingα=β=0 in equation (10), we can obtain the following expression for R:(11)R=r01+∑k=0∞Ck,where r0=μH, μ∈0,1, and H is the height of the tunnel.When solving an actual problem, the mapping function is composed of a finite number of coefficientsCk. To obtain the mapping relationship between the tunnel boundary and the unit circle, a mathematical optimization model is used to determine the undetermined coefficients Ck.To satisfy equation (10), the following objective function is defined:(12)f=∑j=1mrj−Rcosαj−θj+∑k=0nCkcosαj+kθj2.To ensure thatz=ωζ is a univalent function outside the unit circle, the constant term Ck must satisfy the following expression:(13)∑k=1nkCk<1.Thus, the constrained optimization mathematical model of the mapping function is as follows:(14)minfs.t.sinαj−θj+∑k=0nCksinαj+kθj=0∑k=1nkCk<1j=0,1,2,…,mk=0,1,2,…,n.The objective function is the sum of the squared errors between the actual points on the tunnel boundary and the corresponding mapping points on the unit circle after conformal mapping. When the objective function value is approximately zero, the optimal solution is the undetermined coefficientCkk=0,1,2,…,n.


### 3.2. Barrier Function Method Based on fmincon

In this study, a barrier function is used to solve the above optimization problem with constraints. Starting from the interior point, a barrier function is defined to keep the search within the feasible region. The constrained optimization problem is converted into an unconstrained problem, and then the optimization iteration process is used to continuously update the barrier function to make the algorithm converge. The algorithm flow is shown in Figure2.Figure 2
Program diagram of the barrier function.The algorithm convergence condition of the barrier function is as follows:(15)ϕx∗rk,rk−ϕx∗rk−1,rk−1ϕx∗rk−1,rk−1≤ε1,(16)x∗rk−x∗rk−1≤ε2.Between two adjacent iterations, equation (15) requires the relative change of the barrier function value to be sufficiently small, and equation (16) indicates that the unconstrained minimum point is sufficiently close. When the unconstrained minimum point x∗rk that satisfies the convergence condition approaches the optimal point of the original problem, the iteration is terminated. The optimal solution of the original problem is as follows:(17)x∗=x∗rk,fx∗=fx∗rk.The above algorithm was implemented using the fmincon function of the MATLAB software, and the constrained optimization given by equation (14) can be approximated by the following problem:(18)minfx−μ∑ilnsis.thx=0,gx+s=0,where s is a slack variable, lnsi guarantees that s>0, gx+s=0 is equivalent to gx≤0, and μ>0. When μ⟶0, equation (18) is equivalent to the constrained optimization model given by equation (14). During the operation, the fmincon function expands the initial point x0 in the feasible region based on the Taylor series in the neighborhood, and its Hessian matrix is as follows:(19)∇2φx,s,λ,μ=∇2f+∑iλhi∇2hi+∑iλgi∇2gi0JhTJgT0μs20IJh000JgI00,where Jh and Jg are the Jacobian matrices, J=∇f1T⋯∇fnT.The algorithm convergence condition of the fmincon function is(20)max∇φx,s,λ,μ<ε.The Hessian matrix is iteratively calculated until the convergence condition is met. The convergence is equivalent to equations (15) and (16), and the global optimal solution obtained at that time is the undetermined coefficient Ckk=0,1,2,…,n of the mapping function.


### 3.3. Tunnel Mapping Example and Analysis of Convergence Conditions

Following the tunnel clearance section of the “Code for Design of Highway Tunnels (JTG 3370.1-2018),” the barrier function method was used to obtain the mapping of several main types of tunnel boundaries, as shown in Figure3. The red circles in the figure represent the mapping points, and the solid black lines represent the actual boundary of the tunnel. Since the tunnel was axisymmetric, taking the positive half axis of the ordinate x-axis as the starting point, the x-axis was divided into 18 equal parts. That is, m=18 sampling points were selected in a counterclockwise order for mapping, and the number of coefficients Ck of the mapping function was n=6. The mapping coefficients are shown in Table 1.Figure 3
Conformal mapping for the tunnel. (a) Fourth-class highway two-lane tunnel (20 km/h). (b) Third-class highway two-lane tunnel (40 km/h). (c) Second-class highway two-lane tunnel (80 km/h). (d) First-class highway two-lane tunnel (100 km/h). (e) First-class highway three-lane tunnel (120 km/h).
(a)(b)(c)(d)(e)Table 1
Temperature and wildlife count in the three areas covered by the study.
TunnelC0C1C2C3C4C5a0.1207−0.09360.0603−0.03160.00940.0014b0.1596−0.11880.0713−0.03300.00570.0047c0.1457−0.16740.0678−0.02650.00130.0083d0.1846−0.13390.0765−0.03200.00200.0074e0.1472−0.24350.0688−0.0166−0.00130.0075Figure3 shows that the mapping points were very close to the actual outline of the tunnel. To further describe the accuracy of the mapping, the objective function value, average absolute error, maximum absolute error, average relative error, and maximum relative error of the sampling points of the tunnel section were calculated. The objective function value was the sum of the squared errors, the absolute error was the difference between the mapped value and the actual value, and the relative error was the percentage of the absolute error compared to the actual value. As shown in Table 2, the objective function values of the mapping results were all less than 0.015 m2, the average absolute error was in the range of 15–25 mm, and the maximum absolute error was in the range of 45–70 mm. The average relative error was less than 0.5%, and the maximum relative error was less than 1.2%. For the excavation section of an underground tunnel, over- or underexcavation can occur. When the maximum absolute error was used as the convergence condition for determining the mapping function, the mapping accuracy met the construction acceptance standard for the allowable deviation (−30 to 150 mm) of the excavation section of underground tunnel engineering. In this sense, when using the mapping function obtained for this condition to calculate the stress of the surrounding rock of the tunnel, the analytical solution obtained could be considered to be accurate [12].Table 2
Error analysis of the sampling points.
TunnelObjective function value (m2)Average absolute error (m)Maximum absolute error (m)Average relative error (%)Maximum relative error (%)Total length of real tunnel boundary (m)Mapping accuracyδ1 (%)Mapping accuracyδ2 (%)a0.006620.01510.04590.391.1126.17240.0644.12b0.010850.02000.05010.481.1529.19120.0742.96c0.015260.02220.06850.441.1732.84410.0745.94d0.013910.02200.06370.441.1132.90700.0744.12e0.014900.02430.05890.380.8242.01920.0641.22The use of the ratio of the average absolute error of a sampling point to the total side length of a real tunnel boundary was previously proposed as the criterion for judging the mapping accuracy [11], that is,(21)δ1=Δs¯L,where Δs¯ is the average absolute error of the sampling point and L is the total length of the real tunnel boundary. When the mapping accuracy δ1≤5%, the accuracy requirement of engineering calculation was achieved. The values of δ1 shown in Table 2 were all less than 0.1%, and the mapping accuracy met the convergence conditions proposed in this research.The convergence condition of the mapping accuracy was established previously [12] based on the fact that the maximum absolute error meets the engineering acceptance standard, that is, the engineering deviation criterion. The maximum absolute error can reflect the range of the absolute error between the mapped and actual values, but the criterion cannot describe the degree of deviation between the mapped value of each sampling point and the actual value. In a previous publication [11], the ratio of the average absolute error to the circumference was taken as the convergence condition of the mapping accuracy, that is, the average error criterion. However, this criterion did not account for the impact of the random uncertainty of the sampling points on the accuracy. Based on the above deficiencies, in this study, the structural failure criterion of the underground structure reliability analysis theory [24] was used to establish the convergence conditions of the mapping accuracy.The probability of a structure completing a predetermined function under specified conditions is called the reliability of the structure. The reliability and failure of a structure are two mutually exclusive events, and the sum of their probabilities is equal to 1, i.e.,(22)pr+pf=1,where pr represents the reliable probability that the structure completes the predetermined function and pf represents the failure probability that the structure does not complete the predetermined function.X=X1,X2,…,XnT is a basic random variable that affects the structure function. In this study, the random function(23)Z=gx=gX1,X2,…,Xnrepresents the random error between the mapped value and the actual value of each sampling point, which obeyed a normal distribution. Because it is a continuous random variable, the failure probability is defined as follows:(24)pf=PZ≤0=∫−∞0fzzdz.By assuming thatZ∼Nμz,σz2, the mean is μz, the standard deviation is σz2, and the probability density function is fzz, then(25)fzz=12πσzexp−z−μz22σz2.By definingY=Z−μz/σz, Z is transformed into a standard normal distribution variable Y∼N0,1, and its probability density function and cumulative distribution function are as follows:(26)φy=12πexp−y22,Φy=∫−∞yφydy.Thus, equation (24) can be written as(27)pf=∫−∞−μz/σzφydy=Φ−μzσz.The reliability index of the structure is defined as(28)βz=μzσz.This expression states that the standard deviation can be used to measure the distance from the origin to the average. Equation (27) can then be expressed as follows:(29)pf=Φ−βz=1−Φβz.Because the probability density functionfzz is symmetric about the axis z=μz, the relationship between the mapping accuracy δ2 and the failure probability pf of the structure proposed in this research is expressed as follows:(30)δ2=2pf=21−Φβz.Based on the deviation between the error of each sampling point and the mean value, to make the random uncertainty of the mapping error as small as possible, the mean value of the error deviation of each sampling point multiplied by the standard deviationσz was taken as the structural failure criterion to establish the mapping accuracy convergence condition. Thus, δ2≥21−Φ1=31.74%. The values of δ2 in Table 2 are all greater than 41%, satisfying the structural failure criterion proposed in this study.In summary, the barrier function method implemented by the fmincon optimization function is efficient and easy to be implemented, and it produced a high-precision mapping of the tunnel section. The convergence condition of the mapping accuracy satisfied the engineering deviation criterion, average error criterion, and structural failure criterion.


## 3.1. Constrained Optimization Model of Mapping Function

Lettingα=β=0 in equation (10), we can obtain the following expression for R:(11)R=r01+∑k=0∞Ck,where r0=μH, μ∈0,1, and H is the height of the tunnel.When solving an actual problem, the mapping function is composed of a finite number of coefficientsCk. To obtain the mapping relationship between the tunnel boundary and the unit circle, a mathematical optimization model is used to determine the undetermined coefficients Ck.To satisfy equation (10), the following objective function is defined:(12)f=∑j=1mrj−Rcosαj−θj+∑k=0nCkcosαj+kθj2.To ensure thatz=ωζ is a univalent function outside the unit circle, the constant term Ck must satisfy the following expression:(13)∑k=1nkCk<1.Thus, the constrained optimization mathematical model of the mapping function is as follows:(14)minfs.t.sinαj−θj+∑k=0nCksinαj+kθj=0∑k=1nkCk<1j=0,1,2,…,mk=0,1,2,…,n.The objective function is the sum of the squared errors between the actual points on the tunnel boundary and the corresponding mapping points on the unit circle after conformal mapping. When the objective function value is approximately zero, the optimal solution is the undetermined coefficientCkk=0,1,2,…,n.


## 3.2. Barrier Function Method Based on fmincon

In this study, a barrier function is used to solve the above optimization problem with constraints. Starting from the interior point, a barrier function is defined to keep the search within the feasible region. The constrained optimization problem is converted into an unconstrained problem, and then the optimization iteration process is used to continuously update the barrier function to make the algorithm converge. The algorithm flow is shown in Figure2.Figure 2
Program diagram of the barrier function.The algorithm convergence condition of the barrier function is as follows:(15)ϕx∗rk,rk−ϕx∗rk−1,rk−1ϕx∗rk−1,rk−1≤ε1,(16)x∗rk−x∗rk−1≤ε2.Between two adjacent iterations, equation (15) requires the relative change of the barrier function value to be sufficiently small, and equation (16) indicates that the unconstrained minimum point is sufficiently close. When the unconstrained minimum point x∗rk that satisfies the convergence condition approaches the optimal point of the original problem, the iteration is terminated. The optimal solution of the original problem is as follows:(17)x∗=x∗rk,fx∗=fx∗rk.The above algorithm was implemented using the fmincon function of the MATLAB software, and the constrained optimization given by equation (14) can be approximated by the following problem:(18)minfx−μ∑ilnsis.thx=0,gx+s=0,where s is a slack variable, lnsi guarantees that s>0, gx+s=0 is equivalent to gx≤0, and μ>0. When μ⟶0, equation (18) is equivalent to the constrained optimization model given by equation (14). During the operation, the fmincon function expands the initial point x0 in the feasible region based on the Taylor series in the neighborhood, and its Hessian matrix is as follows:(19)∇2φx,s,λ,μ=∇2f+∑iλhi∇2hi+∑iλgi∇2gi0JhTJgT0μs20IJh000JgI00,where Jh and Jg are the Jacobian matrices, J=∇f1T⋯∇fnT.The algorithm convergence condition of the fmincon function is(20)max∇φx,s,λ,μ<ε.The Hessian matrix is iteratively calculated until the convergence condition is met. The convergence is equivalent to equations (15) and (16), and the global optimal solution obtained at that time is the undetermined coefficient Ckk=0,1,2,…,n of the mapping function.


## 3.3. Tunnel Mapping Example and Analysis of Convergence Conditions

Following the tunnel clearance section of the “Code for Design of Highway Tunnels (JTG 3370.1-2018),” the barrier function method was used to obtain the mapping of several main types of tunnel boundaries, as shown in Figure3. The red circles in the figure represent the mapping points, and the solid black lines represent the actual boundary of the tunnel. Since the tunnel was axisymmetric, taking the positive half axis of the ordinate x-axis as the starting point, the x-axis was divided into 18 equal parts. That is, m=18 sampling points were selected in a counterclockwise order for mapping, and the number of coefficients Ck of the mapping function was n=6. The mapping coefficients are shown in Table 1.Figure 3
Conformal mapping for the tunnel. (a) Fourth-class highway two-lane tunnel (20 km/h). (b) Third-class highway two-lane tunnel (40 km/h). (c) Second-class highway two-lane tunnel (80 km/h). (d) First-class highway two-lane tunnel (100 km/h). (e) First-class highway three-lane tunnel (120 km/h).
(a)(b)(c)(d)(e)Table 1
Temperature and wildlife count in the three areas covered by the study.
TunnelC0C1C2C3C4C5a0.1207−0.09360.0603−0.03160.00940.0014b0.1596−0.11880.0713−0.03300.00570.0047c0.1457−0.16740.0678−0.02650.00130.0083d0.1846−0.13390.0765−0.03200.00200.0074e0.1472−0.24350.0688−0.0166−0.00130.0075Figure3 shows that the mapping points were very close to the actual outline of the tunnel. To further describe the accuracy of the mapping, the objective function value, average absolute error, maximum absolute error, average relative error, and maximum relative error of the sampling points of the tunnel section were calculated. The objective function value was the sum of the squared errors, the absolute error was the difference between the mapped value and the actual value, and the relative error was the percentage of the absolute error compared to the actual value. As shown in Table 2, the objective function values of the mapping results were all less than 0.015 m2, the average absolute error was in the range of 15–25 mm, and the maximum absolute error was in the range of 45–70 mm. The average relative error was less than 0.5%, and the maximum relative error was less than 1.2%. For the excavation section of an underground tunnel, over- or underexcavation can occur. When the maximum absolute error was used as the convergence condition for determining the mapping function, the mapping accuracy met the construction acceptance standard for the allowable deviation (−30 to 150 mm) of the excavation section of underground tunnel engineering. In this sense, when using the mapping function obtained for this condition to calculate the stress of the surrounding rock of the tunnel, the analytical solution obtained could be considered to be accurate [12].Table 2
Error analysis of the sampling points.
TunnelObjective function value (m2)Average absolute error (m)Maximum absolute error (m)Average relative error (%)Maximum relative error (%)Total length of real tunnel boundary (m)Mapping accuracyδ1 (%)Mapping accuracyδ2 (%)a0.006620.01510.04590.391.1126.17240.0644.12b0.010850.02000.05010.481.1529.19120.0742.96c0.015260.02220.06850.441.1732.84410.0745.94d0.013910.02200.06370.441.1132.90700.0744.12e0.014900.02430.05890.380.8242.01920.0641.22The use of the ratio of the average absolute error of a sampling point to the total side length of a real tunnel boundary was previously proposed as the criterion for judging the mapping accuracy [11], that is,(21)δ1=Δs¯L,where Δs¯ is the average absolute error of the sampling point and L is the total length of the real tunnel boundary. When the mapping accuracy δ1≤5%, the accuracy requirement of engineering calculation was achieved. The values of δ1 shown in Table 2 were all less than 0.1%, and the mapping accuracy met the convergence conditions proposed in this research.The convergence condition of the mapping accuracy was established previously [12] based on the fact that the maximum absolute error meets the engineering acceptance standard, that is, the engineering deviation criterion. The maximum absolute error can reflect the range of the absolute error between the mapped and actual values, but the criterion cannot describe the degree of deviation between the mapped value of each sampling point and the actual value. In a previous publication [11], the ratio of the average absolute error to the circumference was taken as the convergence condition of the mapping accuracy, that is, the average error criterion. However, this criterion did not account for the impact of the random uncertainty of the sampling points on the accuracy. Based on the above deficiencies, in this study, the structural failure criterion of the underground structure reliability analysis theory [24] was used to establish the convergence conditions of the mapping accuracy.The probability of a structure completing a predetermined function under specified conditions is called the reliability of the structure. The reliability and failure of a structure are two mutually exclusive events, and the sum of their probabilities is equal to 1, i.e.,(22)pr+pf=1,where pr represents the reliable probability that the structure completes the predetermined function and pf represents the failure probability that the structure does not complete the predetermined function.X=X1,X2,…,XnT is a basic random variable that affects the structure function. In this study, the random function(23)Z=gx=gX1,X2,…,Xnrepresents the random error between the mapped value and the actual value of each sampling point, which obeyed a normal distribution. Because it is a continuous random variable, the failure probability is defined as follows:(24)pf=PZ≤0=∫−∞0fzzdz.By assuming thatZ∼Nμz,σz2, the mean is μz, the standard deviation is σz2, and the probability density function is fzz, then(25)fzz=12πσzexp−z−μz22σz2.By definingY=Z−μz/σz, Z is transformed into a standard normal distribution variable Y∼N0,1, and its probability density function and cumulative distribution function are as follows:(26)φy=12πexp−y22,Φy=∫−∞yφydy.Thus, equation (24) can be written as(27)pf=∫−∞−μz/σzφydy=Φ−μzσz.The reliability index of the structure is defined as(28)βz=μzσz.This expression states that the standard deviation can be used to measure the distance from the origin to the average. Equation (27) can then be expressed as follows:(29)pf=Φ−βz=1−Φβz.Because the probability density functionfzz is symmetric about the axis z=μz, the relationship between the mapping accuracy δ2 and the failure probability pf of the structure proposed in this research is expressed as follows:(30)δ2=2pf=21−Φβz.Based on the deviation between the error of each sampling point and the mean value, to make the random uncertainty of the mapping error as small as possible, the mean value of the error deviation of each sampling point multiplied by the standard deviationσz was taken as the structural failure criterion to establish the mapping accuracy convergence condition. Thus, δ2≥21−Φ1=31.74%. The values of δ2 in Table 2 are all greater than 41%, satisfying the structural failure criterion proposed in this study.In summary, the barrier function method implemented by the fmincon optimization function is efficient and easy to be implemented, and it produced a high-precision mapping of the tunnel section. The convergence condition of the mapping accuracy satisfied the engineering deviation criterion, average error criterion, and structural failure criterion.


## 4. Study of Stability of Tunnel-Surrounding Rock in Layered Jointed Rock Mass

### 4.1. Study of Stress Intensity Factors of Crack Tips of Layered Joints around Tunnels

To study the influence of the layered joints around the tunnel on the stability of the surrounding rock of the tunnel, first, the tunnel was taken as the research object to study the stress distribution of the surrounding rock at the boundary of the tunnel. Taking a certain section of the tunnel project in Yibin, Sichuan, as the research background [8], the straight-walled arched tunnel was 5 m wide and 6 m high, and it had a circular arch radius of 2.5 m. It was subjected to a vertical stress of σv=10MPa and a horizontal stress of σh=2.5MPa. The calculation model is shown in Figure 4. The mapping coefficient and the accuracy of the tunnel are shown in Table 3. The convergence conditions met the requirements of the engineering deviation, average error, and structural failure criteria.Figure 4
Tunnel calculation diagram.Table 3
Coefficients and accuracy of the mapping function for the straight wall arched tunnel.
C0C1C2C3C4C5Maximum absolute error (m)Mapping accuracyδ1 (%)Mapping accuracyδ2 (%)−0.10950.07580.0787−0.10020.03850.00410.08800.1345.94To obtain the mapping function, a program was written in MATLAB to obtain the tangential stressσα of the surrounding rock at each boundary point of the straight wall arch tunnel, and the relationship between the tangential stress σα and the polar angle α around the tunnel was plotted, as shown in Figure 5. The top, spandrel, and foot regions showed stress concentrations, the shear stress around the tunnel reached a maximum at the foot region, and the shear stresses in the straight wall and floor plate areas were lower.Figure 5
Relationship between tangential stress and polar angles.The layered joints were parallel or nearly parallel joints that formed in the same geological period. To facilitate analysis, the layered joints around the tunnel were simplified as a crack with an inclination angle ofβ and a length of 2a, as shown in Figure 6. Taking this crack as the research object, the Westergaard stress function of the crack [25] was(31)ZIz=zz2−a2−1σ,ZIIz=zz2−a2−1τ,where the analytical functions ZIz and ZIz represent the Westergaard stress functions of type I and type II cracks, respectively. The mapping function z=ωζ=aζ+ζ−1/2 was used to transform the crack to the unit circle in the complex plane. z=a+ζ is the coordinate of any point in the complex plane, σ is the normal stress perpendicular to the crack, and τ is the shear stress parallel to the crack. The coordinate transformation formula of the stress components are as follows:(32)σ=σαsin2αcos2β,τ=σαsin2αsinβcosβ,where σα is the tangential stress at any point of the tunnel boundary, which can be obtained from the polar stress component σθ in the complex plane using equation (2), α is the polar angle at any point on the tunnel boundary, and β is the crack inclination angle. The stress intensity factor at the crack tip is(33)KI=limz⟶aZI2πz−a,KII=limz⟶aZII2πz−a,where KI represents the stress intensity factor at the tip of the mode I crack and KII represents the stress intensity factor at the tip of the mode II crack.Figure 6
Stress analysis diagram for the crack on the tunnel boundary.By substituting equations (31) and (32) into equation (33), we obtain the following:(34)KI=πaσαsin2αcos2β,KII=πaσαsin2αsinβcosβ.According to the theory of fracture mechanics, the stress intensity factorK is a physical quantity that characterizes the strength of the stress field at the crack tip. When K≥KC, the crack tip is unstable, and the crack expands. KC represents the fracture toughness, which is the mechanical performance index of the material and is related to the structure and composition of the material. KC can be obtained with a rock sample test [26], empirical formula estimate [27], numerical analysis [28], or theoretical model prediction [29].Equation (34) shows that K is determined by α, β, and σα. This formula is the analytic expression of the complex function of the stress intensity factor of the crack tip around the tunnel under a confining pressure. If α and σα remain unchanged, K is controlled by β. If α and β remain unchanged, K is controlled by σα.In this study, the tunnel crack inclination and position model reported previously [7–9] was substituted into equation (34) to calculate and plot the stress intensity factor of the tunnel surrounding the crack tip. The sign of the compressive stress in the fracture mechanics analysis was defined as negative.Figure7(a) shows that the plots of K versus β at the top of the tunnel (α=0°) followed sine and cosine functions. Since KI was controlled by the normal stress perpendicular to the crack surface, when β=0°, i.e., when the crack was parallel to the horizontal plane, the crack was closed due to the pressure, and KI was the smallest at this time. As β gradually increased, the compressive stress of the crack provided by the tunnel edge tangential stress gradually decreased, so KI gradually increased. At β=90°, the crack was perpendicular to the horizontal plane, the tangential stress around the tunnel no longer provided compressive stress for the crack, and KI was zero. Since K was controlled by the shear stress parallel to the crack surface, when the crack inclination angle was β=45°, KII reached its maximum value. KI and KII were equal in value when the dome position was β=45°, which was the dominant structural plane inclination angle for the crack tip instability growth. The conclusions of this work were consistent with the numerical damage simulations and similar model test results reported previously [7]. Thus, when the RPFA software was used to simulate the numerical damage of a tunnel with a crack inclination of β=45°, at the top of the arch, a shear damage area was formed in the crack tip area at the initial stage of loading. The similar model test results showed that the compressive strength of the tunnel model was the lowest when the preset crack inclination was β=45° at the top of the tunnel arch. When the crack inclination angle was 0° and 90°, the compressive strength of the tunnel model was greater.Figure 7
Stress intensity factorsKI and KII. (a) Values of KI and KII for different crack inclination angles β at the top of the arch. (b) Values of KI and KII for different crack inclination angles β for the spandrel and foot. (c) Values of KI and KII when the crack inclination angle was β=45°.
(a)(b)(c)Figure7(b) shows that at the spandrel (α=60°) and foot (α=140°) regions, KII reached the maximum when the crack inclination angle was β=45°, which was consistent with previously reported numerical simulations and similar model test results [8, 9]. Thus, the Tresca stress value of the tunnel numerical model was the largest when the crack inclination angle was β=45° at the spandrel and foot. KII at the crack tip was also the largest, and the similar model test results showed that the crack inclination angle was β=45° at the spandrel and foot. The stress concentration of the tunnel model was greater, and the strength was lower.Figure7(c) shows that when the crack inclination angle was β=45°, the crack tip K changed with the crack at different positions of the tunnel boundary. K reached a larger value at the spandrel, it reached maximum at the foot, it was lower in the straight wall area, and it was close to zero in the top and floor areas. The research conclusions of this work are consistent with previously reported numerical damage simulations and the similar model test results [7]. Thus, when the cracks were located at different locations in the tunnel, the main form of tunnel failure was the compression and shear failure at the crack tip, and the wall foot and spandrel were the weak parts formed by the tunnel crack.In summary, the mechanical behavior of the crack tip around the tunnel under confining pressure was mainly controlled by the type II stress intensity factorKII, which was anisotropic with the change of the crack inclination and position. The relevant parameters of the tunnel calculation model were substituted into equation (34) to calculate and plot KII for different crack inclination angles and positions around the tunnel, as shown in Figure 8. The inclination and position of the cracks around the tunnel affected the value of KII and showed an evident pattern. When the crack inclination angle was β=45°, it was the dominant structural plane inclination angle for the instability and propagation of the crack tip around the tunnel, and the KII value about this symmetry decreased to zero as the inclination angle approached β=0° and β=90°. The KII value of the crack at the same inclination angle around the tunnel changed with the polar angle α of the tunnel, showing a low and high double peak shape. Its double peaks corresponded to the spandrel and foot, and the stress intensity was lower in the straight wall area and close to zero in the top and floor areas.Figure 8
Stress intensity factorKII versus tunnel polar angle and crack inclination angle.


### 4.2. Scoring Parameters of Crack Propagation around Tunnel and Inclination Angle of Rock Joints

After the excavation of an underground tunnel or cavern, the surrounding rock stress is redistributed. The surrounding rock of the tunnel wall is unloaded in the normal direction, and the tangential stress is concentrated. The instability and propagation of the crack tip of the layered joints around the tunnel cause the surrounding rock to produce radial horizontal tensile cracks approximately parallel to the excavation unloading surface, which extend and expand under the action of the bias load. Due to the cutting action of the structural plane, the tension cracks form a rock slab, which gradually deform inward toward the air. With the continuous adjustment of the stress or the influence of the accumulated blasting disturbance nearby, the rock slab separates and falls off from the parent rock from the surface to the inside, finally forming a slab [30]. Therefore, the existence of layered joints around the tunnel reduces the stress threshold of the rib spalling.Because most of the rock was in a mechanical environment under multidirectional compression, the singular stress field at the crack tip was a I-II compound type [31]. According to the theory of maximum circumferential stress, the crack propagated in the θ direction corresponding to the maximum circumferential stress σθmax, and when the circumferential stress in this direction reached a critical value, the crack started to extend. The expression of the stress component in the polar coordinates of the crack front is as follows:(35)σθθ=122πrcosθ2KI1+cosθ−3KIIsinθ,where θ is the crack initiation angle, r is the distance between the microunit and the crack tip, and KI and KII are the stress intensity factors of the type I and II cracks, respectively.θ in equation (35) is differentiated and set equal to zero to obtain(36)KIsinθ+KII3cosθ−1=0.From this, the initiation angleθ=θ0 of the secondary crack can be determined. At that time, the composite stress intensity factor of the crack tip is as follows:(37)K∗=σθmax2πr=12cosθ02KI1+cosθ0−3KIIsinθ0.WhenK∗≥KC, the crack begins to extend.Considering the anisotropy of the stress intensity factor of the crack tip around the tunnel, equation (34) was substituted into equation (36) to obtain(38)tanβ0=sinθ01−3cosθ0.When the secondary crack expands along the direction of the initiation angleθ0 and when it is parallel to the unloading surface of the tunnel excavation, a radial horizontal tension crack parallel to the unloading surface of the excavation forms under the action of the bias load. Combined with the existing literature on the mechanism of tunnel rib spalling, β0 in equation (38) is the dominant expansion angle of the joint plane of the crack tip due to the instability expansion of the slab. Substituting the geometric relationship θ0+180°=α+β0 into equation (38), the relationship between α, β0, and θ0 is as follows:(39)sinβ01−4cosθ0+sinα=0.Based on this, the corresponding relationship betweenα and θ0 with respect to β0 was plotted, as shown in Figure 9.Figure 9
Relationship betweenα and θ0 with respect to β0 (VA: very advantageous, F: favorable, G: general, UF: unfavorable, and VD: very disadvantageous, as shown in Table 4).Table 4
Modified score parameters of joint obliquity.
Joint inclinationVery advantageous (VA)Favorable (F)General (G)Unfavorable (UF)Very disadvantageous (VD)Parameters suggested in this article0°–15°75°–90°20°–30°60°–75°15°–20°30°–45°45°–60°Score value0−2−5−10−12To verify the corresponding relationship between the angles given in Figure9, the finite element software Abaqus was used to perform a numerical simulation to calculate the Tresca stress cloud diagram of the crack tunnel model, as shown in Figure 10. The elastic modulus of the surrounding rock of the tunnel was E=2.95GPa, Poisson’s ratio was μ=0.25, the vertical stress was σv=10MPa, and the horizontal stress was σh=2.5MPa. The cracks in the model were established using the extended finite element module XEFM. To better simulate the relative relationship between the secondary crack and the unloading surface of the tunnel excavation, the default crack length was 2a=5m, and the crack tip was 1 m from the tunnel periphery. The number of model nodes was 8,796, the number of units was 8,626, and the unit type was CPS4R.Figure 10
Tresca stress nephogram of the tunnel models with different crack inclinations. (a) Top positionβ0=75°. (b) Spandrel position β0=45°. (c) Strait wall position β0=30°.
(a)(b)(c)Figure10 shows that when β=β0, the secondary crack was parallel to the unloading surface of the corresponding tunnel, and a shear stress concentration area formed near the secondary crack. The crack tips at different locations around the tunnel corresponded to different β0. The model simulation results were consistent with the above conclusions.Equation (39) was substituted into equation (37) to obtain(40)K∗=12cosα+β02KI1−cosα+β0−32cosα+β02KIIsinα+β0.If the fracture toughnessKC of the jointed rock mass at a certain position around the tunnel was measured for the maximum circumferential stress, then K∗≥KC could be used to judge the possibility of rib spalling when the rock inclination angle was β=β0 at this position.Considering the anisotropy of the stress intensity factor of the layered joints around the tunnel, a modified scoring parameter was proposed for the inclination of the layered joints around the tunnel in the RMR to further strengthen the adverse effects of the inclination on the surrounding rock of the tunnel. As shown in Table4, the initial value of the RMR was corrected based on the joint inclination range of different parts of the tunnel.


### 4.3. Bias Load and Bolt Supporting Stress around Tunnel

Based on the existing literature and the above research, it was determined that the damage to the surrounding rock of the tunnel caused by the bias load formed by the layered joints around the tunnel is reflected in the following two aspects: microscopically, the crack tip with the inclination angle of the dominant expansion joint plane is expanded, and the secondary crack forms a radial horizontal tension crack parallel to the excavation unloading surface. Macroscopically, the fractured rock mass is sheared and destroyed along the bedding layer, causing interlayer slippage around the tunnel. The two aspects act together to affect the stability of the tunnel-surrounding rock.Due to the influences of the excavation unloading and the blasting disturbance, the cohesive forcecj of the jointed rock mass is reduced. If the tunnel boundary in the layered jointed rock mass does not satisfy the joint plane stability discriminant given by equation (8), the fractured rock mass undergoes shear displacement along the joint plane. At that time, the anchor rod in the surrounding rock support system is subjected to axial tension, exerting a normal stress on the joint plane and restricting its normal separation displacement [32]. At the same time, the anchor rod makes full use of the shear strength of the anchored rock mass, and the anchor rod transmits the tensile force to the surrounding rock through the anchoring solid [33], thereby suppressing the instability and propagation of the crack tip. Thus, the anchor rod provides a support stress to the tunnel-surrounding rock to balance the bias load:(41)pb≥σαcosβsinβ−φj−cjcosφjsinβcosβ−φj,where σα is the tangential stress at any point of the tunnel boundary, which can be obtained from the polar stress component σθ in the complex plane using equation (2), α is the polar angle at any point of the tunnel boundary, β is the inclination angle of the joint, cj is the cohesion of the jointed rock mass, and φj is the internal friction angle of the jointed rock mass. The inequality sign in equation (41) was set to an equal sign to calculate the anchoring stress pb at any point on the side of the tunnel.The Yangjiaping Tunnel of the Chenglan Railway was considered as an example [34]. The tunnel was buried at a depth of 330 m, and it was located in a layered phyllite formation. The strata strike was parallel to the tunnel axis, and the strata dip was in the range of 60°–85°. The tunnel was subjected to the vertical stress σv=9.56MPa and a horizontal stress σh=13.39MPa. The jointed rock mass cohesion was cj=0.10MPa, and the internal friction angle was φj=24°. The tunnel was initially supported by anchor rods and a steel-framed concrete secondary lining. The mapping coefficient and the accuracy of the tunnel section are shown in Table 5. The convergence condition satisfied the engineering deviation, average error, and structural failure criteria at the same time. After substituting different joint inclination angles into equation (41), the distribution diagram of the bolt support stress pb was calculated and plotted, as shown in Figure 11.Table 5
Coefficients and accuracy of the mapping function for the Chenglan tunnel.
C0C1C2C3C4C5Maximum absolute error (m)Mapping accuracyδ1 (%)Mapping accuracyδ2 (%)−0.0669−0.00330.0415−0.03560.00930.00480.06050.0841.22Figure 11
Distribution of the support pressurepb.The supporting stress distribution of the spandrel, wall, and foot shown in Figure11 was consistent with the contact pressure of the surrounding rock and the initial support at each monitoring point in the literature [34]. The support stress was anisotropic under the bias load of the inclined layered surrounding rock, and the rock inclination angle of the specific part of the tunnel must be considered when the support design is carried out.


## 4.1. Study of Stress Intensity Factors of Crack Tips of Layered Joints around Tunnels

To study the influence of the layered joints around the tunnel on the stability of the surrounding rock of the tunnel, first, the tunnel was taken as the research object to study the stress distribution of the surrounding rock at the boundary of the tunnel. Taking a certain section of the tunnel project in Yibin, Sichuan, as the research background [8], the straight-walled arched tunnel was 5 m wide and 6 m high, and it had a circular arch radius of 2.5 m. It was subjected to a vertical stress of σv=10MPa and a horizontal stress of σh=2.5MPa. The calculation model is shown in Figure 4. The mapping coefficient and the accuracy of the tunnel are shown in Table 3. The convergence conditions met the requirements of the engineering deviation, average error, and structural failure criteria.Figure 4
Tunnel calculation diagram.Table 3
Coefficients and accuracy of the mapping function for the straight wall arched tunnel.
C0C1C2C3C4C5Maximum absolute error (m)Mapping accuracyδ1 (%)Mapping accuracyδ2 (%)−0.10950.07580.0787−0.10020.03850.00410.08800.1345.94To obtain the mapping function, a program was written in MATLAB to obtain the tangential stressσα of the surrounding rock at each boundary point of the straight wall arch tunnel, and the relationship between the tangential stress σα and the polar angle α around the tunnel was plotted, as shown in Figure 5. The top, spandrel, and foot regions showed stress concentrations, the shear stress around the tunnel reached a maximum at the foot region, and the shear stresses in the straight wall and floor plate areas were lower.Figure 5
Relationship between tangential stress and polar angles.The layered joints were parallel or nearly parallel joints that formed in the same geological period. To facilitate analysis, the layered joints around the tunnel were simplified as a crack with an inclination angle ofβ and a length of 2a, as shown in Figure 6. Taking this crack as the research object, the Westergaard stress function of the crack [25] was(31)ZIz=zz2−a2−1σ,ZIIz=zz2−a2−1τ,where the analytical functions ZIz and ZIz represent the Westergaard stress functions of type I and type II cracks, respectively. The mapping function z=ωζ=aζ+ζ−1/2 was used to transform the crack to the unit circle in the complex plane. z=a+ζ is the coordinate of any point in the complex plane, σ is the normal stress perpendicular to the crack, and τ is the shear stress parallel to the crack. The coordinate transformation formula of the stress components are as follows:(32)σ=σαsin2αcos2β,τ=σαsin2αsinβcosβ,where σα is the tangential stress at any point of the tunnel boundary, which can be obtained from the polar stress component σθ in the complex plane using equation (2), α is the polar angle at any point on the tunnel boundary, and β is the crack inclination angle. The stress intensity factor at the crack tip is(33)KI=limz⟶aZI2πz−a,KII=limz⟶aZII2πz−a,where KI represents the stress intensity factor at the tip of the mode I crack and KII represents the stress intensity factor at the tip of the mode II crack.Figure 6
Stress analysis diagram for the crack on the tunnel boundary.By substituting equations (31) and (32) into equation (33), we obtain the following:(34)KI=πaσαsin2αcos2β,KII=πaσαsin2αsinβcosβ.According to the theory of fracture mechanics, the stress intensity factorK is a physical quantity that characterizes the strength of the stress field at the crack tip. When K≥KC, the crack tip is unstable, and the crack expands. KC represents the fracture toughness, which is the mechanical performance index of the material and is related to the structure and composition of the material. KC can be obtained with a rock sample test [26], empirical formula estimate [27], numerical analysis [28], or theoretical model prediction [29].Equation (34) shows that K is determined by α, β, and σα. This formula is the analytic expression of the complex function of the stress intensity factor of the crack tip around the tunnel under a confining pressure. If α and σα remain unchanged, K is controlled by β. If α and β remain unchanged, K is controlled by σα.In this study, the tunnel crack inclination and position model reported previously [7–9] was substituted into equation (34) to calculate and plot the stress intensity factor of the tunnel surrounding the crack tip. The sign of the compressive stress in the fracture mechanics analysis was defined as negative.Figure7(a) shows that the plots of K versus β at the top of the tunnel (α=0°) followed sine and cosine functions. Since KI was controlled by the normal stress perpendicular to the crack surface, when β=0°, i.e., when the crack was parallel to the horizontal plane, the crack was closed due to the pressure, and KI was the smallest at this time. As β gradually increased, the compressive stress of the crack provided by the tunnel edge tangential stress gradually decreased, so KI gradually increased. At β=90°, the crack was perpendicular to the horizontal plane, the tangential stress around the tunnel no longer provided compressive stress for the crack, and KI was zero. Since K was controlled by the shear stress parallel to the crack surface, when the crack inclination angle was β=45°, KII reached its maximum value. KI and KII were equal in value when the dome position was β=45°, which was the dominant structural plane inclination angle for the crack tip instability growth. The conclusions of this work were consistent with the numerical damage simulations and similar model test results reported previously [7]. Thus, when the RPFA software was used to simulate the numerical damage of a tunnel with a crack inclination of β=45°, at the top of the arch, a shear damage area was formed in the crack tip area at the initial stage of loading. The similar model test results showed that the compressive strength of the tunnel model was the lowest when the preset crack inclination was β=45° at the top of the tunnel arch. When the crack inclination angle was 0° and 90°, the compressive strength of the tunnel model was greater.Figure 7
Stress intensity factorsKI and KII. (a) Values of KI and KII for different crack inclination angles β at the top of the arch. (b) Values of KI and KII for different crack inclination angles β for the spandrel and foot. (c) Values of KI and KII when the crack inclination angle was β=45°.
(a)(b)(c)Figure7(b) shows that at the spandrel (α=60°) and foot (α=140°) regions, KII reached the maximum when the crack inclination angle was β=45°, which was consistent with previously reported numerical simulations and similar model test results [8, 9]. Thus, the Tresca stress value of the tunnel numerical model was the largest when the crack inclination angle was β=45° at the spandrel and foot. KII at the crack tip was also the largest, and the similar model test results showed that the crack inclination angle was β=45° at the spandrel and foot. The stress concentration of the tunnel model was greater, and the strength was lower.Figure7(c) shows that when the crack inclination angle was β=45°, the crack tip K changed with the crack at different positions of the tunnel boundary. K reached a larger value at the spandrel, it reached maximum at the foot, it was lower in the straight wall area, and it was close to zero in the top and floor areas. The research conclusions of this work are consistent with previously reported numerical damage simulations and the similar model test results [7]. Thus, when the cracks were located at different locations in the tunnel, the main form of tunnel failure was the compression and shear failure at the crack tip, and the wall foot and spandrel were the weak parts formed by the tunnel crack.In summary, the mechanical behavior of the crack tip around the tunnel under confining pressure was mainly controlled by the type II stress intensity factorKII, which was anisotropic with the change of the crack inclination and position. The relevant parameters of the tunnel calculation model were substituted into equation (34) to calculate and plot KII for different crack inclination angles and positions around the tunnel, as shown in Figure 8. The inclination and position of the cracks around the tunnel affected the value of KII and showed an evident pattern. When the crack inclination angle was β=45°, it was the dominant structural plane inclination angle for the instability and propagation of the crack tip around the tunnel, and the KII value about this symmetry decreased to zero as the inclination angle approached β=0° and β=90°. The KII value of the crack at the same inclination angle around the tunnel changed with the polar angle α of the tunnel, showing a low and high double peak shape. Its double peaks corresponded to the spandrel and foot, and the stress intensity was lower in the straight wall area and close to zero in the top and floor areas.Figure 8
Stress intensity factorKII versus tunnel polar angle and crack inclination angle.


## 4.2. Scoring Parameters of Crack Propagation around Tunnel and Inclination Angle of Rock Joints

After the excavation of an underground tunnel or cavern, the surrounding rock stress is redistributed. The surrounding rock of the tunnel wall is unloaded in the normal direction, and the tangential stress is concentrated. The instability and propagation of the crack tip of the layered joints around the tunnel cause the surrounding rock to produce radial horizontal tensile cracks approximately parallel to the excavation unloading surface, which extend and expand under the action of the bias load. Due to the cutting action of the structural plane, the tension cracks form a rock slab, which gradually deform inward toward the air. With the continuous adjustment of the stress or the influence of the accumulated blasting disturbance nearby, the rock slab separates and falls off from the parent rock from the surface to the inside, finally forming a slab [30]. Therefore, the existence of layered joints around the tunnel reduces the stress threshold of the rib spalling.Because most of the rock was in a mechanical environment under multidirectional compression, the singular stress field at the crack tip was a I-II compound type [31]. According to the theory of maximum circumferential stress, the crack propagated in the θ direction corresponding to the maximum circumferential stress σθmax, and when the circumferential stress in this direction reached a critical value, the crack started to extend. The expression of the stress component in the polar coordinates of the crack front is as follows:(35)σθθ=122πrcosθ2KI1+cosθ−3KIIsinθ,where θ is the crack initiation angle, r is the distance between the microunit and the crack tip, and KI and KII are the stress intensity factors of the type I and II cracks, respectively.θ in equation (35) is differentiated and set equal to zero to obtain(36)KIsinθ+KII3cosθ−1=0.From this, the initiation angleθ=θ0 of the secondary crack can be determined. At that time, the composite stress intensity factor of the crack tip is as follows:(37)K∗=σθmax2πr=12cosθ02KI1+cosθ0−3KIIsinθ0.WhenK∗≥KC, the crack begins to extend.Considering the anisotropy of the stress intensity factor of the crack tip around the tunnel, equation (34) was substituted into equation (36) to obtain(38)tanβ0=sinθ01−3cosθ0.When the secondary crack expands along the direction of the initiation angleθ0 and when it is parallel to the unloading surface of the tunnel excavation, a radial horizontal tension crack parallel to the unloading surface of the excavation forms under the action of the bias load. Combined with the existing literature on the mechanism of tunnel rib spalling, β0 in equation (38) is the dominant expansion angle of the joint plane of the crack tip due to the instability expansion of the slab. Substituting the geometric relationship θ0+180°=α+β0 into equation (38), the relationship between α, β0, and θ0 is as follows:(39)sinβ01−4cosθ0+sinα=0.Based on this, the corresponding relationship betweenα and θ0 with respect to β0 was plotted, as shown in Figure 9.Figure 9
Relationship betweenα and θ0 with respect to β0 (VA: very advantageous, F: favorable, G: general, UF: unfavorable, and VD: very disadvantageous, as shown in Table 4).Table 4
Modified score parameters of joint obliquity.
Joint inclinationVery advantageous (VA)Favorable (F)General (G)Unfavorable (UF)Very disadvantageous (VD)Parameters suggested in this article0°–15°75°–90°20°–30°60°–75°15°–20°30°–45°45°–60°Score value0−2−5−10−12To verify the corresponding relationship between the angles given in Figure9, the finite element software Abaqus was used to perform a numerical simulation to calculate the Tresca stress cloud diagram of the crack tunnel model, as shown in Figure 10. The elastic modulus of the surrounding rock of the tunnel was E=2.95GPa, Poisson’s ratio was μ=0.25, the vertical stress was σv=10MPa, and the horizontal stress was σh=2.5MPa. The cracks in the model were established using the extended finite element module XEFM. To better simulate the relative relationship between the secondary crack and the unloading surface of the tunnel excavation, the default crack length was 2a=5m, and the crack tip was 1 m from the tunnel periphery. The number of model nodes was 8,796, the number of units was 8,626, and the unit type was CPS4R.Figure 10
Tresca stress nephogram of the tunnel models with different crack inclinations. (a) Top positionβ0=75°. (b) Spandrel position β0=45°. (c) Strait wall position β0=30°.
(a)(b)(c)Figure10 shows that when β=β0, the secondary crack was parallel to the unloading surface of the corresponding tunnel, and a shear stress concentration area formed near the secondary crack. The crack tips at different locations around the tunnel corresponded to different β0. The model simulation results were consistent with the above conclusions.Equation (39) was substituted into equation (37) to obtain(40)K∗=12cosα+β02KI1−cosα+β0−32cosα+β02KIIsinα+β0.If the fracture toughnessKC of the jointed rock mass at a certain position around the tunnel was measured for the maximum circumferential stress, then K∗≥KC could be used to judge the possibility of rib spalling when the rock inclination angle was β=β0 at this position.Considering the anisotropy of the stress intensity factor of the layered joints around the tunnel, a modified scoring parameter was proposed for the inclination of the layered joints around the tunnel in the RMR to further strengthen the adverse effects of the inclination on the surrounding rock of the tunnel. As shown in Table4, the initial value of the RMR was corrected based on the joint inclination range of different parts of the tunnel.


## 4.3. Bias Load and Bolt Supporting Stress around Tunnel

Based on the existing literature and the above research, it was determined that the damage to the surrounding rock of the tunnel caused by the bias load formed by the layered joints around the tunnel is reflected in the following two aspects: microscopically, the crack tip with the inclination angle of the dominant expansion joint plane is expanded, and the secondary crack forms a radial horizontal tension crack parallel to the excavation unloading surface. Macroscopically, the fractured rock mass is sheared and destroyed along the bedding layer, causing interlayer slippage around the tunnel. The two aspects act together to affect the stability of the tunnel-surrounding rock.Due to the influences of the excavation unloading and the blasting disturbance, the cohesive forcecj of the jointed rock mass is reduced. If the tunnel boundary in the layered jointed rock mass does not satisfy the joint plane stability discriminant given by equation (8), the fractured rock mass undergoes shear displacement along the joint plane. At that time, the anchor rod in the surrounding rock support system is subjected to axial tension, exerting a normal stress on the joint plane and restricting its normal separation displacement [32]. At the same time, the anchor rod makes full use of the shear strength of the anchored rock mass, and the anchor rod transmits the tensile force to the surrounding rock through the anchoring solid [33], thereby suppressing the instability and propagation of the crack tip. Thus, the anchor rod provides a support stress to the tunnel-surrounding rock to balance the bias load:(41)pb≥σαcosβsinβ−φj−cjcosφjsinβcosβ−φj,where σα is the tangential stress at any point of the tunnel boundary, which can be obtained from the polar stress component σθ in the complex plane using equation (2), α is the polar angle at any point of the tunnel boundary, β is the inclination angle of the joint, cj is the cohesion of the jointed rock mass, and φj is the internal friction angle of the jointed rock mass. The inequality sign in equation (41) was set to an equal sign to calculate the anchoring stress pb at any point on the side of the tunnel.The Yangjiaping Tunnel of the Chenglan Railway was considered as an example [34]. The tunnel was buried at a depth of 330 m, and it was located in a layered phyllite formation. The strata strike was parallel to the tunnel axis, and the strata dip was in the range of 60°–85°. The tunnel was subjected to the vertical stress σv=9.56MPa and a horizontal stress σh=13.39MPa. The jointed rock mass cohesion was cj=0.10MPa, and the internal friction angle was φj=24°. The tunnel was initially supported by anchor rods and a steel-framed concrete secondary lining. The mapping coefficient and the accuracy of the tunnel section are shown in Table 5. The convergence condition satisfied the engineering deviation, average error, and structural failure criteria at the same time. After substituting different joint inclination angles into equation (41), the distribution diagram of the bolt support stress pb was calculated and plotted, as shown in Figure 11.Table 5
Coefficients and accuracy of the mapping function for the Chenglan tunnel.
C0C1C2C3C4C5Maximum absolute error (m)Mapping accuracyδ1 (%)Mapping accuracyδ2 (%)−0.0669−0.00330.0415−0.03560.00930.00480.06050.0841.22Figure 11
Distribution of the support pressurepb.The supporting stress distribution of the spandrel, wall, and foot shown in Figure11 was consistent with the contact pressure of the surrounding rock and the initial support at each monitoring point in the literature [34]. The support stress was anisotropic under the bias load of the inclined layered surrounding rock, and the rock inclination angle of the specific part of the tunnel must be considered when the support design is carried out.


## 5. Conclusions

(1)
The barrier function method based on the fmincon optimization function in MATLAB could achieve a high-precision mapping of the tunnel section, and the algorithm was efficient and easy to implement. Considering the shortcomings of the existing convergence conditions, in this study, the mean value of the error offset of each sampling point multiplied by the standard deviation was taken as the structural failure criterion, and the convergence condition of the mapping accuracy was established, that is, the convergence whenδ≥31.74%.(2)
The mechanical behavior of the crack tip of the tunnel edge under confining pressure was mainly controlled by the type II stress intensity factorKII. The value of KII was anisotropic with the change in the crack inclination and position. When the crack inclination angle was β=45°, it was the dominant structural plane inclination angle for the instability and propagation of the crack tip around the tunnel, and the KII value about this symmetry line decreased to zero as the inclination approached β=0° and β=90°. The KII value of the crack at the same inclination angle around the tunnel changed with the polar angle α of the tunnel, showing a low and high double peak shape. The double peaks corresponded to the spandrel and foot, and the stress intensity was lower in the straight wall area and close to zero in the top and floor areas.(3)
When the crack inclination wasβ=β0, the secondary crack was parallel to the unloading surface of the corresponding tunnel, and a shear stress concentration area was formed near the secondary crack. The crack tips at different parts of the tunnel side corresponded to different inclination angles β0 of the superior expansion joint plane. The existence of layered joints at the edge of the tunnel reduced the stress threshold of the failure of the rib spalling. Based on this, the revised scoring parameters for the inclination of the layered joints around the tunnel in the RMR classification were proposed (Table 4) to further strengthen the adverse effects of the inclination on the surrounding rock of the tunnel.(4)
The damage to the surrounding rock of the tunnel caused by the bias load formed by the layered joints around the tunnel was reflected in the following two aspects: microscopically, the crack tip with the inclination angle of the dominant expansion joint plane was expanded, and the secondary crack formed a radial horizontal tension crack parallel to the excavation unloading surface. Macroscopically, the fractured rock mass was sheared and destroyed along the bedding layer, causing interlayer slippage around the tunnel. Both aspects affected the stability of the tunnel-surrounding rock. Using the complex function of the stability of the layered joint plane at any point of the tunnel boundary, the supporting stress required for balancing the bias load at each point on the tunnel boundary could be calculated.



---

*Source: 5473059-2020-10-16.xml*