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

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 modiﬁed scoring parameter of the layered joints around the tunnel for the rock geomechanics classiﬁcation (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 K II at the tip of the mode II crack indicated that the corresponding cracks at the various points around the tunnel had inconsistent inﬂuences on the tunnel. The position had the greatest inﬂuence, followed by the straight wall area, and then by the top and ﬂoor 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


Introduction
A layered jointed rock mass is a sedimentary rock with a layered structure composed of an array of parallel joint planes. e 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 tunnelsurrounding cracks to the surrounding rock of straight wall arch tunnels. ese simulations were combined with model tests to study the damage and destruction behaviors under confining pressure. e 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. ese 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. e 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. e anisotropy of the stress intensity factor due to crack inclination and position changes was studied. e 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. e 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. 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]. e mapping function can be expanded in the form of a Laurent series [20] as follows:

Analysis of Complex
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 C k (k � 0, 1, 2, . . . , n . . .) represents the kth constant term of the Laurent series, which depends the shape of the tunnel section. ese are real constants when discussing axisymmetric problems.

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. e stress component around the tunnel after conformal mapping can be expressed in the complex plane as ; ζ is a multivalued function defined as ζ � e iθ � cos θ + i sin θ; ω(ζ) 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. e analytical functions φ(ζ) and ψ(ζ) in equation (2) are defined as follows: 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]. e stress boundary condition is where a is the integral starting point taken on the boundary, b is any point on the boundary, and F x and F y 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 F x and F y are both zero, and the integral term on the right-hand side of the equation is zero.

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. e discriminant of the stability of the joint plane is as follows [23]: where σ 1 and σ 3 represent the large and small principal stresses around the tunnel, β is the dip angle of the joint, c j 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. e excavated tunnel is radially unloaded, and the tangential stress is concentrated, i.e., Mathematical Problems in Engineering 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 where σ α can be obtained from the polar stress component σ θ in the complex plane using equation (2), and equation (7) can be substituted to obtain where each item on the left-hand side is a function of ζ.
Equation (8) can be expressed as follows: e coordinate of any point on the boundary of the tunnel is expressed as A j (r j , α j ), and the coordinate of the mapping point on the corresponding unit circle is B j (1, θ j ), as shown in Figure 1. e 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 − 1 (r j , α j ) � ω − 1 [ζ(1, θ j )] can be obtained as follows: e 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 c j 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.

Solution for Tunnel Boundary
Mapping Function

Constrained Optimization Model of Mapping Function.
Letting α � β � 0 in equation (10), we can obtain the following expression for R: where r 0 � μ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 coefficients C k . To obtain the mapping relationship between the tunnel boundary and the unit circle, a mathematical optimization model is used to determine the undetermined coefficients C k .
To satisfy equation (10), the following objective function is defined: To ensure that z � ω(ζ) is a univalent function outside the unit circle, the constant term C k must satisfy the following expression: us, the constrained optimization mathematical model of the mapping function is as follows: e 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 coefficient C k (k � 0, 1, 2, . . . , n).

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. e 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. e algorithm flow is shown in Figure 2.
e algorithm convergence condition of the barrier function is as follows: 4 Mathematical Problems in Engineering 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 * (r k ) that satisfies the convergence condition approaches the optimal point of the original problem, the iteration is terminated. e optimal solution of the original problem is as follows: e 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: where s is a slack variable, ln(s i ) guarantees that s > 0, g(x) + s � 0 is equivalent to g(x) ≤ 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 x 0 in the feasible region based on the Taylor series in the neighborhood, and its Hessian matrix is as follows: where J h and J g are the Jacobian matrices, . e algorithm convergence condition of the fmincon function is e Hessian matrix is iteratively calculated until the convergence condition is met. e convergence is equivalent to equations (15) and (16), and the global optimal solution obtained at that time is the undetermined coefficient C k (k � 0, 1, 2, . . . , n) of the mapping function.

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 Figure 3. e 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.
at is, m � 18 sampling points were selected in a counterclockwise order for mapping, and the number of coefficients C k of the mapping function was n � 6. e mapping coefficients are shown in Table 1. Figure 3 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. e 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 m 2 , 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. e average relative error was less than 0.5%, and the maximum relative error was less than 1.2%. For the excavation section Begin Generate a feasible initial point x 0 , a penalty factor r 0 , a reduction coefficient c, and a convergence accuracy ε.
Let the number of iterations k← 0 Is the convergence condition satisfied? Mathematical Problems in Engineering 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]. e 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, 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. e 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. e 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. e 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. e probability of a structure completing a predetermined function under specified conditions is called the reliability of the structure. e reliability and failure of a structure are two mutually exclusive events, and the sum of their probabilities is equal to 1, i.e., where p r represents the reliable probability that the structure completes the predetermined function and p f represents the failure probability that the structure does not complete the predetermined function. X � [X 1 , X 2 , . . . , X n ] T is a basic random variable that affects the structure function. In this study, the random function Z � g(x) � g X 1 , X 2 , . . . , X n (23) represents 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: By assuming that Z ∼ N(μ z , σ 2 z ), the mean is μ z , the standard deviation is σ 2 z , and the probability density function is f z (z), then By defining Y � (Z − μz)/σ z , Z is transformed into a standard normal distribution variable Y ∼ N(0, 1), and its probability density function and cumulative distribution function are as follows: us, equation (24) can be written as Table 1: Temperature and wildlife count in the three areas covered by the study.
e reliability index of the structure is defined as is 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: Because the probability density function f z (z) is symmetric about the axis z � μ z , the relationship between the mapping accuracy δ 2 and the failure probability p f of the structure proposed in this research is expressed as follows: 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.
us, 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. e convergence condition of the mapping accuracy satisfied the engineering deviation criterion, average error criterion, and structural failure criterion.

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 � 10 MPa and a horizontal stress of σ h � 2.5 MPa. e calculation model is shown in Figure 4. e mapping coefficient and the accuracy of the tunnel are shown in Table 3. e convergence conditions met the requirements of the engineering deviation, average error, and structural failure criteria.
To 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. e 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. e 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 where the analytical functions Z I (z) and Z I (z) represent the Westergaard stress functions of type I and type II cracks, respectively. e 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. e coordinate transformation formula of the stress components are as follows: σ � σ α sin 2 α cos 2 β, τ � σ α sin 2 α 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. e stress intensity factor at the crack tip is where K I represents the stress intensity factor at the tip of the mode I crack and K II represents the stress intensity factor at the tip of the mode II crack. By substituting equations (31) and (32) into equation (33), we obtain the following: According to the theory of fracture mechanics, the stress intensity factor K is a physical quantity that characterizes the strength of the stress field at the crack tip. When K ≥ K C , the crack tip is unstable, and the crack expands. K C represents the fracture toughness, which is the mechanical performance index of the material and is related to the structure and composition of the material. K C can be obtained with a rock sample test [26], empirical formula estimate [27], numerical analysis [28], or theoretical model prediction [29]. 8 Mathematical Problems in Engineering Equation (34) shows that K is determined by α, β, and σ α . is 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][8][9] was substituted into equation (34) to calculate and plot the stress intensity factor of the tunnel surrounding the crack tip. e sign of the compressive stress in the fracture mechanics analysis was defined as negative. Figure 7(a) shows that the plots of K versus β at the top of the tunnel (α � 0°) followed sine and cosine functions. Since K I 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 K I 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 K I 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 K I was zero. Since K was controlled by the shear stress parallel to the crack surface, when the crack inclination angle was β � 45°, K II reached its maximum value. K I and K II were equal in value when the dome     position was β � 45°, which was the dominant structural plane inclination angle for the crack tip instability growth. e conclusions of this work were consistent with the numerical damage simulations and similar model test results reported previously [7]. us, 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. e 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(b) shows that at the spandrel (α � 60°) and foot (α � 140°) regions, K II 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]. us, the Tresca stress value of the tunnel numerical model was the largest when the crack inclination angle was β � 45°at the spandrel and foot. K II 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. e stress concentration of the tunnel model was greater, and the strength was lower. Figure 7(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. e research conclusions of this work are consistent with previously reported numerical damage simulations and the similar model test results [7]. us, 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 factor K II , which was anisotropic with the change of the crack inclination and position. e relevant parameters of the tunnel calculation model were substituted into equation (34) to calculate and plot K II for different crack inclination angles and positions around the tunnel, as shown in Figure 8. e inclination and position of the cracks around the tunnel affected the value of K II 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 K II value about this symmetry decreased to zero as the inclination angle approached β � 0°and β � 90°. e K II 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.

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. e surrounding rock of the tunnel wall is unloaded in the normal direction, and the tangential stress is concentrated. e 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]. erefore, 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. e expression of the stress component in the polar coordinates of the crack front is as follows: where θ is the crack initiation angle, r is the distance between the microunit and the crack tip, and K I and K II 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 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: When K * ≥ K C , 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 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: sin β 0 1 − 4 cos θ 0 + sin α � 0.

(39)
Based on this, the corresponding relationship between α and θ 0 with respect to β 0 was plotted, as shown in Figure 9.
To verify the corresponding relationship between the angles given in Figure 9, 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. e elastic modulus of the surrounding rock of the tunnel was E � 2.95 GPa, Poisson's ratio was μ � 0.25, the vertical stress was σ v � 10 MPa, and the horizontal stress was σ h � 2.5 MPa.
e 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 � 5 m, and the crack tip was 1 m from the tunnel periphery. e number of model nodes was 8,796, the number of units was 8,626, and the unit type was CPS4R. Figure 10 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. e crack tips at different locations around the tunnel corresponded to different β 0 . e model simulation results were consistent with the above conclusions. Equation (39) was substituted into equation (37) to obtain If the fracture toughness K C of the jointed rock mass at a certain position around the tunnel was measured for the maximum circumferential stress, then K * ≥ K C 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 Table 4, the initial value of the RMR was corrected based on the joint inclination range of different parts of the tunnel.

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. e two aspects act together to affect the stability of the tunnelsurrounding rock.
Due to the influences of the excavation unloading and the blasting disturbance, the cohesive force c j 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. us, the anchor rod provides a support stress to the tunnel-surrounding rock to balance the bias load: 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, c j is the cohesion of the jointed rock mass, and φ j is the internal friction angle of the jointed rock mass. e inequality sign in equation (41) was set to an equal sign to calculate the anchoring stress p b at any point on the side of the tunnel. e Yangjiaping Tunnel of the Chenglan Railway was considered as an example [34]. e tunnel was buried at a depth of 330 m, and it was located in a layered phyllite formation. e strata strike was parallel to the tunnel axis, and the strata dip was in the range of 60°-85°. e tunnel was subjected to the vertical stress σ v � 9.56 MPa and a horizontal stress σ h � 13.39 MPa. e jointed rock mass cohesion was c j � 0.10 MPa, and the internal friction angle was φ j � 24°. e tunnel was initially supported by anchor rods and a steel-framed concrete secondary lining. e mapping coefficient and the accuracy of the tunnel section are shown in Table 5. e 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 p b was calculated and plotted, as shown in Figure 11. e supporting stress distribution of the spandrel, wall, and foot shown in Figure 11 was consistent with the contact pressure of the surrounding rock and the initial support at each monitoring point in the literature [34]. e 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.  Table 4).

Conclusions
(1) e 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) e mechanical behavior of the crack tip of the tunnel edge under confining pressure was mainly controlled by the type II stress intensity factor K II . e value of K II 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 K II value about this symmetry line decreased to zero as the inclination approached β � 0°and β � 90°. e K II 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. e 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.
e crack tips at different parts of the tunnel side corresponded to different inclination angles β 0 of the superior expansion joint plane. e 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) e 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.

Data Availability
e data used to support the findings of this research are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that they have no conflicts of interest.