A full-plane strain complex variable analytical model considering three-dimensional principal stress rotation and its engineering applications

Abstract

To explore the influence of stress space rotation on roadway stability, this study establishes Cartesian and principal stress axes, obtains principal stress characteristic parameters and direction cosine matrix, and determines the conversion process between coordinate axes. By breaking through the traditional plane strain mechanics model, a full plane strain model of rectangular roadway considering the deflection effect of principal stress axis is established. The analytical solution of stress complex variation of any element around the roadway is derived, and the calculation method of plastic zone is given. The influence of stress deflection on the stress of roadway and the shape distribution of plastic zone is discussed and the control system of roadway under the action of principal stress deflection in spatial region is established. The research results show that: (1) the shear stress and normal stress formed in the remote of the roadway during the stress rotation are in the same order of magnitude, and the plastic quantitative calculation of the roadway cannot simply consider the normal stress. (2) Stress deflection causes the distribution of stress and plastic zone around the roadway to show ‘asymmetric’ failure characteristics, and the transient deflection of stress may cause the rapid expansion of the plastic zone of the roadway. (3) The distribution characteristics of plastic zone calculated by this study are compared with the deformation and failure feature and peep results of the return air roadway in Cuncaota Mine, and the comparison results are good. The supporting system and technology of roadway under the action of principal stress deflection in space area are put forward, and the field application effect is good.

Introduction

The plastic failure mode of the surrounding rock becomes more and more complex due to the diversification of mining depth and geological conditions. The safety evaluation and failure zone prediction of roadway service in the whole cycle have always been the focus of attention in the engineering field (Pu et al. 2024; Zuo et al. 2024; Liu et al. 2024b, c, d, a). The quantification of plastic zone of roadway usually includes the establishment of analytical solution model and finite element model (Zhang et al. 2023; Wang et al. 2023). The analytical solution model has the characteristics of fast calculation speed and convenient use when supplemented by mathematical software. At the same time, it can be used as a verification and supplement of the finite element model. Therefore, the analytical solution model is an important means to preliminarily evaluate the shape of plastic zone. Through the stress loading level, the existing analytical models can be divided into two-way equal pressure model, two-way unequal pressure model (Peng et al. 2024; Ma et al. 2023), three-dimensional hydrostatic pressure model and three-dimensional non-hydrostatic pressure model (Liu et al. 2023b, a). However, these models are generally based on the principal stress axis is coaxial with the Cartesian coordinate axis, and mainly for the circular roadway section, which is difficult to accurately reflect the complex stress state of the roadway in practical engineering.

Engineering practice shows, under the action of mining, the stress principal axis and the earth principal axis will shift from the coaxial state to the non-coaxial state (Cui et al. 2024). Wang et al. (2020a, b) proposed that under the strong mining, the surrounding rock will have the phenomenon of ' dynamic migration of mining-induced stress '. According to the stress rotation characteristics, the advancing direction of the working face is optimized. Kang et al. (2019) and Su et al. (2024) used hollow inclusion strain gauge to monitor the dynamic changing process of principal stress azimuth and dip angle of roadway during excavation and mining, and obtained the stress deflection characteristics of roadway under excavation and mining disturbance. Liu et al. (2024b, c, d, a) and Song et al. (2024) carried out a hollow torsion shear test, and simulated the failure modes of rock mass under different stress paths by applying torque to rock mass. Liu et al. (2023b, a) proposed a calculation process of stress rotation angle on the side of the stope under plane strain conditions by establishing an equivalent hole model. The above studies have found that the principal stress deflection characteristics of roadway are common, and different stress deflection paths lead to different damage of roadway. Therefore, when establishing the analytical solution model, the case of non-coaxial stress should be considered.

The existing analytical solution model has some limitations in analyzing the principal stress deflection problem of roadway far-field stress. Most of these models are based on the plane strain assumption of circular roadways, which simplifies the problem to a two-dimensional stress state, that is, only the maximum and minimum principal stresses in the plane are considered. However, in practical engineering, the principal stress deflection is the result of the combined action of three-dimensional spatial stress field. This simplified treatment has two significant defects. First, the existing model can only reflect the plane rotation of the principal stress around the center point of the roadway, but cannot characterize the true rotation characteristics of the principal stress in three-dimensional space. Secondly, the model completely ignores the rotation effect of the third principal stress and its influence on the mechanical response of surrounding rock. This theoretical limitation makes it difficult for the existing methods to accurately describe the real mechanical response of non-circular roadways such as rectangular roadways in three-dimensional stress fields.

Therefore, through the complex variable theory, the analytical method of plastic zone in rectangular roadway under three-dimensional non-hydrostatic pressure environment is determined. Combined with the stress path equation, the analytical model of rectangular roadway considering spatial stress deflection effect is established. The correctness and adaptability of the theoretical analysis are verified by engineering examples. Based on the principal stress rotation characteristics, the supporting system and supporting supporting technology of roadway surrounding rock under the action of principal stress deflection in spatial region are proposed. This study provides a method for calculating the plastic zone under the effect of principal stress axis deviation and provides some engineering guidance value for the rapid support design of roadway.

Stress conversion under the action of principal stress deflection in spatial region

Stress deflection path calculation

The stress vector includes stress magnitude and stress direction. In order to determine the stress direction, the Cartesian coordinate system is first introduced as a reference. The Cartesian coordinate system is composed of the x, y, and z-axes, and the planes formed by the coordinate axes, namely the xoy plane, yoz plane, and xoz plane. The rectangular roadway is parallel to the y-axis. The roadway mechanical analysis model is shown in Fig. 1.

[See PDF for image]

Fig. 1

Mechanical model of roadway under spatial principal stress deflection

Where: the principal stresses are denoted as σ₁ (maximum), σ₂ (intermediate), and σ₃ (minimum), respectively. β₁, α₁, β₂, α₃, β₃ and α₃ are the azimuth and dip angle of σ₁, σ₂ and σ₃, respectively.

The stress axis composed of three-dimensional principal stress directions is called the principal stress axis. Based on the Cartesian coordinate axis, the principal directions of principal stresses σ₁, σ₂ and σ₃ can be quantified by characteristic parameters. The angle of the projection of the principal stress axis on the Cartesian axis xoy plane deviating from the y-axis is set to the azimuth angle β_i, which is positive clockwise. The angle between the principal stress axis and its projection on the horizontal plane is set as the dip angle α_i, which is positive counterclockwise. Therefore, it can be determined that the direction cosines of β_i and α_i relative to the Cartesian coordinate axis are:

1

$$\text{C}={\left|\begin{array}{ccc}{\text{L}}_{\text{i}}& {\text{M}}_{\text{i}}& {\text{N}}_{\text{i}}\end{array}\right|}_{3\times 3}=\left|\begin{array}{ccc}cos{\alpha }_{1}sin{\beta }_1& cos{\alpha }_{1}cos{\beta }_1& sin{\alpha }_1\\ cos{\alpha }_{2}sin{\beta }_2& cos{\alpha }_{2}cos{\beta }_2& sin{\alpha }_2\\ cos{\alpha }_{3}sin{\beta }_3& cos{\alpha }_{2}cos{\beta }_2& sin{\alpha }_3\end{array}\right|$$where L_i, M_i, and N_i are the cosines of the principal stress axes relative to the Cartesian axes in the x, y, and z axes, respectively, C is the direction cosine matrix.

Liu et al. (2024b, c, d, a) gave the rotation formula of the stress conversion relationship between the new and old coordinate axes. Therefore, the stress component of the roadway under the Cartesian coordinate axis can be obtained. The principal stress axis is expressed as:

2

$${\left|{\sigma }_{\text{x}}^{\prime },{\sigma }_{\text{y}}^{\prime },{\sigma }_{\text{z}}^{\prime },{\tau }_{\text{xy}}^{\prime },{\tau }_{\text{xz}}^{\prime },{\tau }_{\text{yz}}^{\prime }\right|}^{\text{T}}{\text{= C}}^{\text{T}}{\left|{\sigma }_1,{\sigma }_2,{\sigma }_3\right|}^{\text{T}}\text{C}$$where σ'_x, σ'_y, σ'_z, τ'_xy, τ'_yz and τ'_zx represent the stress components under the Cartesian coordinate axis. C^T is the transpose of matrix C.

Substituting Eqs. (1) into (2), each stress component under the Cartesian coordinate axis can be obtained as follows:

3

$$\begin{array}{c}{\sigma }_{\text{x}}^{\prime }={cos}^2{\alpha }_1·{sin}^2{\beta }_1·{\sigma }_1+{cos}^2{\alpha }_2·{sin}^2{\beta }_2·{\sigma }_2+{cos}^2{\alpha }_3·{sin}^2{\beta }_3·{\sigma }_3\\ {\sigma }_{\text{y}}^{\prime }={cos}^2{\alpha }_1·{cos}^2{\beta }_1·{\sigma }_1+{cos}^2{\alpha }_2·{cos}^2{\beta }_2·{\sigma }_2+{cos}^2{\alpha }_3·{cos}^2{\beta }_3·{\sigma }_3\\ {\sigma }_{\text{z}}^{\prime }={sin}^2{\alpha }_1·{\sigma }_1+{sin}^2{\alpha }_2·{\sigma }_2+{sin}^2{\alpha }_3·{\sigma }_3\\ {\tau }_{\text{xy}}^{\prime }={cos}^2{\alpha }_1·sin{\beta }_1·cos{\beta }_1·{\sigma }_1+{cos}^2{\alpha }_2·sin{\beta }_2·cos{\beta }_2·{\sigma }_2+{cos}^2{\alpha }_3·sin{\beta }_3·cos{\beta }_3·{\sigma }_3\\ {\tau }_{\text{yz}}^{\prime }=cos{\alpha }_1·sin{\beta }_1·sin{\alpha }_1·{\sigma }_1+cos{\alpha }_2·sin{\beta }_2·sin{\alpha }_2·{\sigma }_2+cos{\alpha }_3·sin{\beta }_3·sin{\alpha }_3·{\sigma }_3\\ {\tau }_{\text{zx}}^{\prime }=cos{\alpha }_1·cos{\beta }_1·sin{\alpha }_1·{\sigma }_1+cos{\alpha }_2·cos{\beta }_2·sin{\alpha }_2·{\sigma }_2+cos{\alpha }_3·cos{\beta }_3·sin{\alpha }_3·{\sigma }_3\\ \end{array}$$

Rectangular complex variable analytical solution of full plane strain model

The principal stress function around the roadway is the basis to analysis the roadway's stress and plastic failure. The stress component of any unit is set as follows:

4

$$\text{D}=\left[\begin{array}{ccc}{\sigma }_{\text{x}}& {\tau }_{\text{xy}}& {\tau }_{\text{xz}}\\ {\tau }_{\text{yx}}& {\sigma }_{\text{y}}& {\tau }_{\text{yz}}\\ {\tau }_{\text{zx}}& {\tau }_{\text{zy}}& {\sigma }_{\text{z}}\end{array}\right]$$where D is the matrix expression of stress component at any point around the roadway, σ_x, σ_y, σ_z, τ_xy, τ_yx, τ_xz, τ_zx, τ_yz are the stress components of any element. Shear stress reciprocal theorem shows: τ_xz = τ_zx, τ_yz = τ_zy, τ_xy = τ_yx.

Expanding Eq. (4) by determinant can get 0:

5

$${\sigma }^3-{\text{I}}_1{\sigma }^2-{\text{I}}_2\sigma -{\text{I}}_3=0$$where I₁, I₂ and I₃ are the first, second, and third invariants of the stress tensor, respectively. Their expressions are:

6

$$\left\{\begin{array}{c}{\text{I}}_1=({\sigma }_{\text{x}}+{\sigma }_{\text{y}}+{\sigma }_{\text{z}})\\ {\text{I}}_2=-({\sigma }_{\text{y}}{\sigma }_{\text{z}}+{\sigma }_{\text{z}}{\sigma }_{\text{x}}+{\sigma }_{\text{x}}{\sigma }_{\text{y}}-{\tau }_{\text{xy}}^2-{\tau }_{\text{yz}}^2-{\tau }_{\text{zx}}^2)\\ {\text{I}}_3={\sigma }_{\text{x}}{\sigma }_{\text{y}}{\sigma }_{\text{z}}-{\sigma }_{\text{x}}{\tau }_{\text{yz}}^2-{\sigma }_{\text{y}}{\tau }_{\text{zx}}^2-{\sigma }_{\text{z}}{\tau }_{\text{xy}}^2+2{\tau }_{\text{xy}}{\tau }_{\text{yz}}{\tau }_{\text{zx}}\\ \end{array}\right)$$

Equation (5) is a cubic equation. Based on the Cardano's formula, the σ₍₁₎, σ₍₂₎ and σ₍₃₎ at arbitrary positions around the roadway can be obtained (Lv and Zhang 2007). Equation (6) indicates that the principal stress characteristic equation of roadway is a function of σ_x, σ_y, σ_z, τ_xy, τ_yz and τ_zx. Therefore, It is necessary to derive the stress function at any point around the roadway.

At present, the common rectangular roadway model is the plane strain model. In this paper, considering that the roadway is in true triaxial stress, the roadway is arranged along the y-axis of the Cartesian coordinate axis, and the following assumptions are made: (a) ignore the volume force and gravity in the surrounding rock; (b) The surrounding rock of the roadway is an elastomer with unified lithology; (c) The stress and strain of each section of the roadway along the y-axis direction are the same. Therefore, based on the full plane strain problem, the mechanical model of roadway can be divided into plane strain model (σ_x, σ_z and τ_zx), out-of-plane shear model (τ_yz and τ_xy) and the uniaxial compression model (σ_y).

The dimensions of the roadway are defined by its width, denoted as W, and its height, designated as H, by mapping, the mapping function can be determined to be as follow (Ma et al. 2022):

7

$$\left\{\begin{array}{c}w(\epsilon )=R(\epsilon +C_1{\epsilon }^{-1}+C_3{\epsilon }^{-3}+C_5{\epsilon }^{-5})\\ R=\frac{W}{1+C_1+C_3+C_5},\frac{H}{W}=\frac{1-C_1+C_3-C_5}{1+C_1+C_3+C_5}\\ \end{array}\right)$$where R, C₁, C₃ and C₅ are parameters related to roadway section parameters W and H, and w (ε) is a mapping function.

The stress function of the rectangular full plane strain model first needs to give the complex variable basis function. The expressions of φ₀(ε), φ(ε), ψ₀(ε) and ψ(ε) can be obtained by complex variable calculation:

8

$${\phi }_0(\epsilon )=-S_1{\epsilon }^{-1}-S_3{\epsilon }^{-3}-2BR\left(C_1{\epsilon }^{-1}+C_2{\epsilon }^{-3}+C_5{\epsilon }^{-5}\right)-\left(B^{\prime }-i{C}^{\prime }\right)R{\epsilon }^{-1}$$

9

$$\phi (\epsilon )=-S_1{\epsilon }^{-1}-S_3{\epsilon }^{-3}-BR\left(C_1{\epsilon }^{-1}+C_2{\epsilon }^{-3}+C_5{\epsilon }^{-5}\right)+BR\epsilon -\left(B^{\prime }-i{C}^{\prime }\right)R{\epsilon }^{-1}$$

10

$${\psi }_0\left(\epsilon \right)=\frac{-\overline{\omega }\left(\frac{1}{\epsilon }\right)}{{\omega }^{\prime }\left(\frac{1}{\epsilon }\right)}{\phi }_0^{\prime }\left(\epsilon \right)+\overline{S_1}\epsilon +\overline{S_3}{\epsilon }^3+\overline{S_5}{\epsilon }^5-2BR{\epsilon }^{-1}-\left(B^{\prime }+i{C}^{\prime }\right)R\left(C_1{\epsilon }^{-1}+C_3{\epsilon }^{-3}+C_5{\epsilon }^{-5}\right)$$

11

$$\psi \left(\epsilon \right)=\frac{-\overline{\omega }\left(\frac{1}{\epsilon }\right)}{{\omega }^{\prime }\left(\frac{1}{\epsilon }\right)}+\overline{S_1}\epsilon +\overline{S_3}{\epsilon }^3+\overline{S_5}{\epsilon }^5+\left(B^{\prime }+i{C}^{\prime }\right)R\epsilon -2BR{\epsilon }^{-1}$$where B, B' and C' are functions in relation to the far-field boundary stresses σ_x, σ_z and τ_zx. S₁ and S₃ are coefficients. φ₀(ε), φ(ε), ψ₀(ε) and ψ(ε) analytic functions of complex potentials.

(1) Plane strain model

The expressions of σ_x, σ_z and τ_zx are Eqs. (12) and (13) when the plane strain problem is obtained by the complex variable theory:

12

$${\sigma }_x+{\sigma }_z=4\text{Re}\left[\frac{{\phi }^{\prime }(\epsilon )}{w^{\prime }(\epsilon )}\right]$$

13

$${\sigma }_z-{\sigma }_x+2i{\tau }_{\mathit{zx}}=\frac{2}{{\omega }^{\prime }(\epsilon )}\left[\overline{\omega (\epsilon )}\left(\frac{{\phi }^{\prime }(\epsilon )}{w^{\prime }(\epsilon )}\right)+{\psi }^{\prime }(\epsilon )\right]$$

The stress functions of σ_x, σ_z and τ_zx around the roadway are obtained by deriving Eqs. (7), (9) and (11) into Eqs. (12) and (13).

(2) Out-of-plane shear model

Fan (1990) gave the expression of shear stress τ_xy and τ_yz at any point around the roadway under out-of-plane shear:

14

$${\tau }_{\text{xy}}-\text{i}{\tau }_{\text{yz}}=\frac{{\tau }_{\text{xy}}^{\prime }\left(1-\frac{1}{{\epsilon }^2}\right)-\text{i}{\tau }_{\text{yz}}^{\prime }\left(1-\frac{1}{{\epsilon }^2}\right)}{1-\sum _{\text{k}=1}^{\text{n}}{\text{kC}}_{\text{k}}{\epsilon }^{-k-1}}$$

Substituting Eq. (7) into (14), the stress functions of τ_xy and τ_yz are:

15

$${\mathrm{\tau }}_{\text{xy}}=\text{Re}\left[\frac{{\tau }_{\text{xy}}^{\prime }\left(1-\frac{1}{{\epsilon }^2}\right)-{i{\tau }^{\prime }}_{\text{yz}}\left(1-\frac{1}{{\epsilon }^2}\right)}{1-{\text{c}}_1{\epsilon }^{-2}-3{\text{c}}_3{\epsilon }^{-4}-5{\text{c}}_3{\epsilon }^{-6}}\right]$$

16

$${\mathrm{\tau }}_{\mathit{yz}}=-\text{Im}\left[\frac{{\tau }_{\text{xy}}^{\prime }(1-\frac{1}{{\epsilon }^2})-i{\tau }_{\mathit{yz}}^{\prime }(1-\frac{1}{{\epsilon }^2})}{1-c_1{\epsilon }^{-2}-3c_3{\epsilon }^{-4}-5c_3{\epsilon }^{-6}}\right]$$

(3) Uniaxial compression model

From the physical equation, the relationship between σ_y and σ_x, σ_z can be formulated as:

17

$${\mathrm{\sigma }}_{\text{y}}{=E\epsilon }_{\text{y}}\text{+ v}({\mathrm{\sigma }}_{\text{x}}{+\sigma }_{\text{z}})$$where E is the elastic modulus, v is the Poisson's ratio, and ε_y is the axial strain.

After the excavation of the roadway, the stress distribution shows a dynamic change. The stress values of σ_y, σ_x and σ_z in the deep roadway are close to the far field boundary stress. that is, σ_x = σ'_x, σ_y = σ'_y, σ_z = σ'_z. The expression of Eε_y is obtained by substituting σ'_x, σ'_y and σ'_z into the Eq. (17):

18

$$E{\epsilon }_y={{\sigma }^{\prime }}_y-\text{v}({{\sigma }^{\prime }}_x{+{\sigma }^{\prime }}_z)$$

Substituting Eqs. (18) into (17), we get:

19

$${\mathrm{\sigma }}_y={{\sigma }^{\prime }}_y-\text{v}({{\sigma }^{\prime }}_x+{{\sigma }^{\prime }}_z)+v({\sigma }_x+{\sigma }_z)$$

Therefore, based on the Eqs. (12), (13), (15), (16) and (19), the functions σ_x, σ_y, σ_z, τ_xy, τ_yz and τ_zx at any point around the rectangular roadway can be obtained. By substituting the stress function into the Eq. (5), the principal stress of any element can be obtained.

Stress and plastic evolution characteristics of roadways under different stress paths

Stress scheme setting

The location of the three-dimensional principal stresses in space can be described by two characteristic parameters: the azimuth angle and the dip angle. In order to establish an intuitive stress rotation scheme, It is essential to define the relationship between the stress characteristic parameters.

The sum of cosine square vectors in the x, y and z axes in three dimensions is 1. Therefore, the direction cosines in Eq. (1) can be expressed as Eqs. (20), (21), (22) (Liu et al. 1991):

20

$${\text{L}}_1^2+{\text{L}}_2^2+{\text{L}}_3^2=1$$

21

$${\text{M}}_1^2+{\text{M}}_2^2+{\text{M}}_3^2=1$$

22

$${\text{N}}_1^2+{\text{N}}_2^2+{\text{N}}_3^2=1$$

The relationship between the dip angles of σ₁, σ₂ and σ₃ can be satisfied by Eq. (23):

23

$${sin}^2{\alpha }_1+{sin}^2{\alpha }_2+{sin}^2{\alpha }_3=1$$

If the α₁ and α₂ are given, the α₃ satisfies the Eq. (24):

24

$${\alpha }_3=arcsin\sqrt{{cos}^2{\alpha }_1-{sin}^2{\alpha }_2}$$

The direction cosine satisfies the orthogonal relationship pairwise perpendicular and can be expressed as:

25

$$\begin{array}{c}{\text{L}}_1{\text{L}}_2+{\text{M}}_1{\text{M}}_2+{\text{N}}_1{\text{N}}_2=0\\ {\text{L}}_2{\text{L}}_3+{\text{M}}_2{\text{M}}_3+{\text{N}}_2{\text{N}}_3=0\\ {\text{L}}_1{\text{L}}_3+{\text{M}}_1{\text{M}}_3+{\text{N}}_1{\text{N}}_3=0\\ \end{array}$$

The simultaneous Eqs. (1) and (25) obtain:

26

$$\left\{\begin{array}{c}cos{\alpha }_1·sin{\beta }_1·cos{\alpha }_2·sin{\beta }_2+sin{\alpha }_1·sin{\alpha }_2+cos{\alpha }_1·cos{\beta }_1·cos{\alpha }_2·cos{\beta }_2=0\\ cos{\alpha }_2·sin{\beta }_2·cos{\alpha }_3·sin{\beta }_3+sin{\alpha }_2·sin{\alpha }_3+cos{\alpha }_2·cos{\beta }_2·cos{\alpha }_3·cos{\beta }_3=0\\ \end{array}\right)$$

If the α₁, α₂ and β₁ are known, the β₂ can be obtained from the first equation of Eq. (26) as follows:

27

$${\beta }_2={\beta }_1-arccos\left(-tan{\alpha }_1·tan{\alpha }_2\right)$$

If the α₂, α₃ and β₂ are known, the β₃ satisfies the Eq. (28):

28

$${\beta }_3={\beta }_2-arccos\left(-tan{\alpha }_2·tan{\alpha }_3\right)$$

Therefore, when the α_i and β_i of any principal stress are determined, the α_i and β_i parameters of all principal stresses can be obtained based on the α_i of any of the other two principal stresses.

Taking the in-situ stress test data as an example (Cai et al. 2013): σ₁, σ₂ and σ₃ are β₁ = 228.1°, β₂ = 60.50° and β₃ = 137.50°, and the dip angles are α₁ = 13.30°, α₂ = 76.20° and α₃ = -−2.50°. The β₁ and α₁ of σ₁ are fixed, and the α₂ of σ₂ is set to −50°, −30°, −10°, 10°, 30° and 50°. The β₁, α₁ and α₂ are inserted into Eqs. (24), (27) and (28) respectively to obtain the characteristic parameters of principal stress rotation under different α₂, as shown in Table 1. Based on the principle of stereographic projection, the characteristic parameters of the principal stress can be expressed in the Wulff net projection network, which can reflect the rotational path of the σ₁, σ₂ and σ₃. The stereographic projection is shown in Fig. 2.

Table 1. The principal stress rotation characteristic parameters

Plain	α1 (°)	β1 (°)	α2 (°)	β2 (°)	α3 (°)	β3 (°)
Plan 1	13.30	228.10	− 50	154.46	36.88	127.88
Plan 2	13.30	228.10	− 30	145.94	56.61	117.09
Plan 3	13.30	228.10	− 10	140.49	73.25	86.35
Plan 4	13.30	228.10	10	135.71	73.25	9.85
Plan 5	13.30	228.10	30	130.26	56.61	− 20.89
Plan 6	13.30	228.10	50	121.74	36.88	− 31.68

[See PDF for image]

Fig. 2

Characteristics of principal stress distribution based on Wulff net projection network

During the rotation of the principal stresses, the magnitudes of σ₁, σ₂ and σ₃ are kept constant at 45, 30, and 20 MPa, respectively. The stress characteristic parameters β_i, α_i (i = 1–3) and the stress value are substituted into Eq. (3) to obtain the calculation results of the remote stress component as shown in Fig. 3.

[See PDF for image]

Fig. 3

The distribution characteristics of far-field normal stress and shear stress under the action of principal stress deflection

From the external load components of the two groups of roadways, the far-field stress will form different degrees of tangential stress components (τ_xy, τ_yz and τ_zx) during the rotation of the principal stress. Taking the numerical values of τ_xy and σ_z as examples, τ_xy accounts for 37.37%, 34.80%, 32.42%, 32.01% and 33.88% of σ_z under the five groups of schemes, respectively. It shows that the tangential stress component generated by spatial stress rotation is close to or even more than 1/3 of the normal stress in value, and its influence cannot be ignored. Ignoring these tangential stress components may lead to insufficient assessment of roadway support and surrounding rock deformation, which in turn affects the safety of the project. Therefore, accurate calculation and analysis of the shear stress component at the far field is very important to improve the accuracy of roadway design and engineering safety.

The size and shape of stress distribution around the roadway

Set W and H as 4 and 3 m, and the lithology parameters of the roadway are cohesion of C = 2 MPa, internal friction angle of φ = 28°, and Poisson's ratio of v = 0.2. The remote field stress components of the roadway in Fig. 3 are substituted into the Eqs. (12), (13), (15), (16) and (19) to obtain the values of σ_x, σ_y, σ_z, τ_xy, τ_yz and τ_zx. Based on the root formula, the σ₁, σ₂ and σ₃ are obtained. The size and distribution of the σ₁ are finally obtained through the transformation of rectangular coordinates and polar coordinates, as shown in Fig. 4.

[See PDF for image]

Fig. 4

Distribution and morphological characteristics of principal stress in roadway under different principal stress rotation conditions

At present, the stress values formed at the four corners of the roadway obtained in the existing research are completely consistent, and the stress characteristics around the roadway are completely symmetrically distributed with the symmetry axes of 90° and 180°(Shan et al. 2024; Zhao et al. 2015). This is because the common roadway model does not consider the rotation effect of the principal stress, and the model does not give the shear stress effect in the far field. Figure 4 shows that the stress distribution presents a continuous undulating ' wave ' feature, showing a maximum value at the roadway's four corners. The stress peaks at the corner points near 36.87° and 216.87° are 252.53, 247.64, 233.47, 217.26, 207.54 and 208.70 MPa, respectively, at 143.13° and 323.13°. The stress peaks near are 140.31, 127.53, 133.25, 154.65, 180.38 and 196.43 MPa, respectively. The peak stress near the corner decreases first and subsequently escalates as the magnitude of α₂ rises. The difference of the peak stress at the two extreme points is 112.22, 120.08, 100.22, 62.61, 27.16 and 12.27 MPa, respectively. That is, the increase of α₂ leads to a general trend of decreasing difference in peak stress at the two extreme points, and the stress distribution characteristics gradually change from ' asymmetric ' to ' symmetric ', indicating that the far-field tangential stress controls the stress rotatio.

The stress form can directly reflect the stress concentration position of the roadway. Under the influence of stress concentration at the corner, the stress distribution pattern of the roadway presents an X-type expansion as a whole. With the increase of α₂ value, the expansion characteristics of X-type single wing decrease, and the overall stress distribution pattern transits to symmetrical distribution. It shows that stress rotation will result in the expansion of specific areas of roadway or the transfer of stress concentration areas. Therefore, the stress rotation has an important influence on the stress distribution of the roadway. It must be fully considered in the design and analysis process to ensure that the roadway control can effectively deal with the stress changes and ensure the long-term stability of the roadway.

Prediction of the failure zone in the roadway

Kastner proposed to substitute the elastic stress function into the classical yield strength criterion, which can roughly predict the plastic failure state of roadway surrounding rock. At present, the commonly used classical strength criteria include Mohr–Coulomb criterion, Hoek–Brown or Drucker-Prager criterion, Zienkiewice-Pande criterion and so on. Liu et al. (2023b, a) explored the influence of different strength criteria on the stability of roadway surrounding rock, and obtained the consistent characteristics of the morphological evolution characteristics of the plastic zone of roadway surrounding rock under different strength criteria. In this paper, the common M-C criterion is taken as an example to explore the failure characteristics of plastic zone of roadway surrounding rock under different principal stress rotation. The common M-C criterion expression is the relationship between σ₍₁₎ and σ₍₃₎, which does not consider the influence of intermediate principal stress. Therefore, this paper uses the three-dimensional M-C criterion to judge the plastic failure state around the roadway.

29

$$\sqrt{\frac{J_2}{3}}sin{\theta }_{\sigma }tan\phi -\sqrt{J_2}\frac{cos{\theta }_{\sigma }}{cos\phi }-\frac{{\sigma }_{(1)}+{\sigma }_{(2)}+{\sigma }_{(3)}}{3}tan\phi -c=0$$where θ_σ, J₂ and J₃ are the Lode angle, the second deviatoric stress invariant, and the third deviatoric stress invariant, respectively. Their expressions are:

30

$${\theta }_{\sigma }=\frac{1}{3}{sin}^{-1}\left(-\frac{3\sqrt{3}}{2}·\frac{J_3}{\sqrt{J_2^3}}\right)$$

31

$$J_2=\frac{1}{6}\left[{({\sigma }_{(1)}-{\sigma }_{(2)})}^2+{({\sigma }_{(2)}-{\sigma }_{(3)})}^2+{({\sigma }_{(1)}-{\sigma }_{(3)})}^2\right]$$

32

$$J_3=\frac{1}{3}\left[(2{\sigma }_{(1)}-{\sigma }_{(2)}-{\sigma }_{(3)})(2{\sigma }_{(2)}-{\sigma }_{(3)}-{\sigma }_{(1)})(2{\sigma }_{(3)}-{\sigma }_{(1)}-{\sigma }_{(2)})\right]$$

Therefore, the implicit formulation representing the plastic region of the roadway can be expressed as:

33

$$\text{T}=\left({\sigma }_1,{\sigma }_2,{\sigma }_3,{\beta }_{\text{i}},{\alpha }_{\text{i}},\text{W},\text{H},\text{c},\phi ,\text{v}\right)$$

The above stress characteristic parameters, stress magnitude, roadway section parameters and surrounding rock mechanical parameters are inserted into Eq. (33) to traverse the elastoplastic state of each point around the roadway, and the boundaries of all failure points are converted from rectangular coordinates to polar coordinates to determine the failure boundaries. At the same time, in order to proof the rationality of the theoretical analysis, the equivalent circle method and the complex variable analytical are compared. The plastic zone solution method of the equivalent circular tunnel model mainly includes: (1) Stress conversion between the principal stress and the Cartesian coordinate axes; (2) The irregular roadway is equivalent to the circumcircle equivalent circle to determine the equivalent radius. (3) Determine the normal stress σ_r, σ_t, σ_v and shear stress τ_rt, τ_rv, τ_tv of any element around the roadway; (4) Based on the polar coordinate stress function, the principal stress of any element around the roadway is obtained, and the damage range is determined based on Eq. (29). The calculation results of plastic zone are shown in Fig. 5.

[See PDF for image]

Fig. 5

Characteristics of plastic zone of surrounding rock under different principal stress rotation angles based on equivalent circle and complex variable analytical solution

Comparing the results of equivalent circle and complex variable analytical solution, from the perspective of the overall deflection of morphology, under the action of principal stress rotation, the deflection degree of plastic zone morphology obtained by the two calculation methods is consistent. Based on the equivalent circle method, the overall failure orientation and damage range of rectangular roadway can be roughly judged. From the area of plastic zone, there is a certain deviation between the two calculation results. This is because the complex variable theory considers the mechanical bearing characteristics of the corner and the boundary of the rectangular roadway, and the equivalent circle simplifies the boundary of the roadway. At the same time, the plastic zone area obtained based on the equivalent circle method is larger than the complex variable analytical solution. The reason for this is that the radius in the equivalent circle method is derived from the radius of the circumscribed circle of the rectangular roadway. When the aspect ratio (irregularity) is larger, the roadway section area of the equivalent circle method is much larger than the actual roadway section area, which will lead to a larger error in the calculation result, and the calculation results are more conservative. Therefore, in the calculation of complex geometric shapes and boundary conditions, the accuracy of the calculation results obtained by the complex variable analytical solution is relatively higher.

The maximum failure polar angles of the roadway are near 80.83°, 85.55°, 88.62°, 86.27°, 85.07 and 86.98° respectively, and the polar angles of the maximum failure points are close to 90°. From the point of view of the shape, the overall plastic zone morphology is deflected from the first quadrant and then the deflection degree gradually decreases with the increasing α₂. The morphology of the roadway's plastic zone steadily changes from ' asymmetric ' to ' symmetric ', which aligns with the evolutionary features of the principal stress morphology above.

According to the traditional plastic zone theory of roadway surrounding rock, it is considered that a straight line of any length is made along the center of the circular roadway. When there are countless intersections, 4 intersections and 8 intersections with the boundary of the plastic zone, the plastic zone of the surrounding rock of the roadway is circular, elliptical and butterfly, respectively. The contour of the plastic zone is a central symmetric graph, and there are at least two symmetric axes. However, this definition is based on the plane two-dimensional stress of the surrounding rock of the roadway, without considering the spatial stress rotation of the roadway. It can be seen from the plastic zone distribution map obtained by the two calculation methods that the stress rotation will lead to the asymmetric change of the plastic zone shape of the roadway. Different from the common asymmetric failure of the roadway, the plastic zone shape shows a one-way expansion trend, and the plastic zone contour is only central symmetry, no symmetry axis. This asymmetry of the single wing not only changes the shape of the plastic zone, but also may affect the stability of the roadway and the failure mode during excavation.

Therefore, incorporating the effects of principal stress rotation within the spatial domain is very important for accurately predicting the plastic failure mode of the roadway, which can provide a more reliable foundation for support design and avoid the potential safety hazards that traditional calculation methods may bring. In the actual roadway support design, the support design should be based on the failure depth of the plastic zone at different locations, combined with the alienation characteristics of the surrounding rock, especially the key support control in the asymmetric area of the plastic zone, to ensure that the support design can deal with the stress distribution and deformation requirements under different surrounding rock conditions, and enhance safety of the roadway.

Roadway support system and technology under the influence of principal stress deflection in space area

Engineering background

The Cuncaota mining is located in Ordos City, China. The main research object is the return airway of 31205 working face of 3₁ coal seam. The cross-section shape of the return airway of the working face is rectangular, and the width and height of the cross-section are 5.4 m and 3.6 m, respectively. The average depth is 300 m. The roof of the working face is sandy mudstone, siltstone, fine-grained sandstone, siltstone, fine-grained sandstone and medium-grained sandstone, and the average layer thickness is 2.5, 11.53, 2.58, 3.12 and 9.86 m, respectively. The bottom plates are sandy mudstone, siltstone and medium grained sandstone, with an average thickness of 3.78, 9.58 and 17.36 m, respectively. The drilling core is processed and based on the rock mechanics experiment, the lithology parameters of surrounding rock are obtained, such as bulk modulus K, shear modulus G, density ρ, cohesion c, internal friction angle φ and tensile strength σ_t, as shown in Fig. 6. The roadway should be used as both the auxiliary transportation roadway in the upper section and the air inlet roadway in the section during the whole service period. Therefore, it is affected by multiple disturbances during the whole service period, and the roadway has a serious mine pressure phenomenon.

[See PDF for image]

Fig. 6

Roadway layout position and surrounding rock parameters

Stress is the basis for analyzing the plastic failure of roadway. The existing in-situ stress measurement technology can accurately measure the original rock stress state, and the stress state under the action of mining is mostly determined by numerical inversion (Liu et al. 2024b, c, d, a). This study establishes a numerical calculation model with x, y and z of 811, 1000 and 100 m, respectively. 150 m and 100 m coal pillars are reserved along the x and y axes to mitigate boundary effects. z-axis positive direction of the model adopts stress constraint, and the rest positions adopt displacement constraint. The stress boundary values are: S_zz = 7.5 MPa, S_xx = S_yy = 8.25 MPa. The constitutive model adopts the M-C constitutive model, and the parameters are shown in Fig. 6. In order to fully consider the compaction characteristics of the goaf, The goaf is backfilled using the double yield model, and the surrounding rock values in the goaf are obtained by the inversion trial and error method as shown in Table 2.

Table 2. Double yield values of goaf

G (GPa)	K (GPa)	ρ (kg/m3)	c (MPa)	ψ (°)
2.8	4.5	1800	0.001	7

Evolution process of principal stress vector

In this study, taking 400 m from the open-off cut of 31205 working face as an example, a 50 m monitoring line is arranged along the Y axis in the roadway center, and the σ₁, σ₂ and σ₃ are extracted. Under the influence of multiple disturbances, the σ₁, σ₂ and σ₃ are 20.42 MPa, 12.43 MPa and 11.91 MPa, respectively. As the distance from the working face increases, mining effects on the roadway decline, σ₁ and σ₃ gradually decrease and return to the residual mining stress state in the upper section, while σ₂ increases first and then decreases. Figure 7 shows that during the advancing process of the working face, the roof of the goaf will gradually bend and sink with the increase of the goaf area in the upper section and the working face in this section, and finally break into neatly arranged rock blocks above the goaf. The hinge effect between the rock blocks will lead to the deflection of the three-dimensional stress orientation of the coal and rock mass in a certain range in space.

[See PDF for image]

Fig. 7

Evolution characteristics of principal stress of roadway

Equation (33) indicated that the stress vector contains the principal value σ_i, the stress azimuth β_i and the dip angle α_i (i = 1–3). FLAC^3D cannot directly call the dip angle and azimuth angle of the element, so the stress characteristic parameters are obtained based on FLAC^3D embedded Fish language. The main writing steps include: (1) Call the σ₁, σ₂ and σ₃ at the center line of the roadway and the cosine values L₁, L₂ and L₃ of the X-axis corresponding to the principal stress, the cosine values N₁, N₂ and N₃ of the Y-axis, and the cosine values M₁, M₂ and M₃ of the z-axis. (2) Based on Eq. (1), the third formula is used to obtain the dip angles α₁, α₂ and α₃ through the anti-sine function. (3) The first and second formulas of the simultaneous Eq. (1) are obtained based on the arctangent function to obtain the azimuth angles β₁, β₂ and β₃. Based on the Fish language, the stress characteristic parameters along the center line of the roadway are illustrated in Fig. 8.

[See PDF for image]

Fig. 8

Evolution law of principal stress characteristic parameters of roadway

After the advancement of the upper section mining, under the impact of residual mining stress, the azimuth and dip angle of σ₁, σ₂ and σ₃ of roadway on the 50 m monitoring line change little, and the overall change is linear. Under the influence of mining, an obvious secondary mining influence area is formed, and the stress parameters vary significantly within the affected zone. Under the action of lateral mining, the dip angles of σ₁, σ₂ and σ₃ in front of the working face are 58.35°, 0.10° and -31.65° respectively, and the azimuth angles are -89.33°, 0.84° and -89.09° respectively. Under the action of superimposed mining, the dip angles of σ₁, σ₂ and σ₃ are 14.35°, 5.3° and 19.17° respectively, and the azimuth angles are 21.82°, 36.16° and 39.09° respectively. The principal stress characteristic parameters are drawn in the spatial micro unit, and the stress rotation feature at the center line of the roadway is obtained as shown in Fig. 9.

[See PDF for image]

Fig. 9

Principal stress rotation path at the center of roadway under primary mining and repeated disturbance

Figure 8 and Fig. 9 show that the superimposed mining influence of 31,204 and 31,205 working faces leads to the enhancement of stress rotation effect in front of the working face, and the variation in azimuth and dip angle becomes more significant. As the distance from the working face increases, the principal stress characteristic parameter curve gradually approaches the residual mining influence of the upper section. The stress rotation speed gradually decreases, and the rotation trajectory gradually returns to the upper section mining influence.

Roadway failure characteristics and support optimization

To determine the plastic failure condition of the roadway in the forward section, considering the region 10 m in front of the working face as a case study, σ₁, σ₂ and σ₃ of the center line of the roadway at 10 m ahead of the working face are 18.54 MPa, 12.24 MPa and 11.20 MPa, respectively. The dip angles α₁, α₂ and α₃ are 67.63°, 5.90° and − 21.49°, respectively, and the azimuth angles are − 71.63°, 32.93° and − 54.74°, respectively. The principal stress value and characteristic parameters are substituted into Eq. (3) to obtain the stress values in the Cartesian coordinate system: σ_x = 12.46 MPa. σ_y = 12.03 MPa, σ_z = 17.49 MPa, τ_xy = 0.15 MPa, τ_yz = 0.90 MPa, τ_zx = -2.39 MPa. Based on the above, c = 1.2 MPa, φ = 31°, v = 0.159, W = 5.4 m, H = 3.6 m. The stress parameters, roadway surrounding rock parameters and size parameters are inserted into Eq. (33) to obtain the plastic zone morphology at 10 m ahead of the working face, as shown in Fig. 10. It should be noted that since the advancing direction is the opposite direction of the y axis during modeling, the stress parameters extracted by numerical simulation are based on the positive orientation of the y-axis. Therefore, it is necessary to mirror the obtained plastic zone morphology.

[See PDF for image]

Fig. 10

Comparison of roadway failure characteristics based on theoretical calculation, numerical calculation and field measurement

Due to the impact of superimposed mining, the principal stress in the space area of the roadway is deflected, which leads to the inhomogeneity of the stress vector of the roadway. The inhomogeneity is mainly manifested in the inhomogeneity of stress direction and stress size. The inhomogeneity of stress direction determines the easy failure position and asymmetric characteristics of roadway surrounding rock, and the inhomogeneity of stress size determines the failure depth of roadway. Theoretical calculation results show that the regional stress deflection of the roadway under the influence of mining leads to the asymmetric failure characteristics of the roadway. In order to verify the correctness and adaptability of the theoretical analysis, by applying the slicing function, the plastic failure features at 10 m in front of the working face are analyzed. Compared with the theoretical analysis and numerical calculation results, the plastic zone morphology of the roadway is similar to the ' butterfly ' failure morphology. This plastic zone morphology is different from the circular, elliptical and butterfly plastic zones calculated by ignoring the stress direction, mainly due to the three-dimensional tangential stress formed by the deflection of the principal stress. The plastic failure depth of roadway roof and floor is as follows: (1) Roof: coal wall > coal pillar; (2) Floor: coal pillar > coal wall. However, there are some differences between the results of theoretical calculation and numerical calculation in the upper left corner and the lower right corner of the roadway. This is because the numerical calculation is the result of the progressive calculation of the cell layers around the roadway, and the theoretical calculation result is based on the stress state of the center line of the roadway. Therefore, in the local position of the roadway, the calculation results of the two methods are different, but the two are consistent in the overall plastic zone failure contour. Borehole peeping was carried out on the left, middle and right sides of the roadway roof respectively. The peeping results show that the asymmetric failure characteristics of the roof are consistent with the theoretical calculation results. In addition, the failure of the side and bottom plates photographed on site is also consistent with the calculation results, which verifies the rationality and adaptability of the calculation method of plastic zone considering stress deflection in this paper. Based on the comparison of the plastic zone results obtained by the complex variable analysis method and the equivalent circle method, under the action of considering the regional stress deflection, the plastic zone morphology of the roadway surrounding rock obtained by the equivalent circle method also presents an irregular part, but the overall irregularity is small. The macroscopic distribution of the plastic zone is similar to the micro-ellipse or circular distribution, which is quite different from the numerical calculation, on-site peeping and on-site failure.

Based on the calculation results and the initial support design, the primary causes of roadway failure are identified as: (1) The actual failure depth of the roadway is greater than the plastic zone of the original support design. After the large-scale deformation of the roadway, the anchor cable lacks a stable force point, resulting in the failure of the overall extraction of the on-site support body. (2) The three-dimensional shear stress generated by the stress rotation leads to the formation of torsion and shear stress in the supporting body. The poor shear and torsion resistance of the supporting body leads to the shear failure of the supporting body. (3) The effect of stress rotation leads to asymmetric plastic deformation failure of the surrounding rock, and symmetrical support is adopted in the support design. The original support design lacks the key reinforcement of the severely damaged area. (4) Under the influence of strong mining, the stress concentration degree of the shallow part of the roadway surrounding rock is high, and the supporting stress formed by the supporting body such as bolt and cable can not completely resist the high mining stress.

In view of the failure reasons of roadway support, the roadway support system considering the rotation of principal stress is proposed as Fig. 11, which mainly includes (Liu et al. 2024b, c, d, a): (1) Improve the length of the support body. The length of the support body of the roadway roof and side is increased from 2100 to 6500 mm. (2) Improve the diameter and energy absorption characteristics of anchor cable. The φ28.6 mm large-diameter energy-absorbing steel strand anchor cable is used to improve the torsional shear performance of the support by increasing the diameter. (3) Based on the theoretical calculation and on-site asymmetric failure characteristics, two φ28.6 mm × 6500 mm steel strand anchor cables are added to the roof and two sides of the roadway respectively, and the positions where the roadway is easy to fall are emphatically anchored. (4) Five sets of super-strong single hydraulic supports (ZQL2 × 22,500/22/38D) are used in the advance position of the return air roadway, and the secondary expansion of the plastic zone is reduced by the high-strength support resistance. (5) Through directional long borehole hydraulic fracturing, the strength of shallow surrounding rock is weakened, the stress concentration degree of shallow part is reduced, and the high concentrated stress of shallow part is transferred to the deep part of surrounding rock. (Xu et al. 2023).

[See PDF for image]

Fig. 11

Superimposed mining roadway support failure analysis and strengthening support design

After strengthening the support, nine groups of observation points are arranged in the hidden danger area of the roof fall in the return air roadway, and the stress on the support body of the roadway is measured by the CMSW6 mine intrinsically safe non-destructive testing instrument, as shown in Fig. 12. In the whole observation period, the stress state of the anchor cable at each measuring point is greater than the initial working resistance, and the supporting effect of the supporting body is good. In addition, there is no roof fall during the mining period, which shows the reliability of the collaborative support technology of ' anchor cable + roadway drilling pressure relief + single hydraulic support in advance working face '.

[See PDF for image]

Fig. 12

Roadway optimization design field effect verification

Conclusions

1. The principal stress characteristic parameters and direction cosine are determined. The transformation relationship between the Cartesian coordinate axis and the principal stress axis is established, and the stress transformation equation is obtained. The stress deflection effect is introduced into the full plane strain model of roadway, and the normal stress and shear stress functions characterizing the stress state of the surrounding element of rectangular roadway are obtained. The approximate calculation method of plastic zone of roadway considering stress rotation is given.

2. The deflection stress in the spatial area leads to the formation of shear stress components in the roadway, resulting in an ' asymmetric ' failure pattern in the roadway. The transient deflection of the principal stress may cause the rapid expansion of the plastic zone of the roadway surrounding rock.

3. Combined with the engineering case, the theoretical calculation results are verified based on the field failure mode, numerical calculation and peeping results of the roadway, and the comparison results are good. Aiming at the rotation characteristics of principal stress in roadway, the corresponding surrounding rock support system and supporting technology are put forward, which provides a new idea and method for solving the problem of roadway control under complex stress environment.

The full-plane complex variable analytical model considering the spatial stress rotation established in this paper is based on the homogeneous and isotropic surrounding rock, without considering the factors such as stratum structure. In the next study, the factors such as stratum structure (such as dip angle) and heterogeneity will be further combined to improve the model to be closer to the engineering practice. At the same time, the adaptability of different strength criteria is different, and the adaptability of the stress rotation model and the strength criterion will be further explored in the future.

Acknowledgements

This work was supported by the National Natural Science Foundation of China (51774288, 52004289 and U22A20165), the Fundamental Research Funds for the Central Universities ( 2022XJNY01, BBJ2024001).

Author Contributions

ZH Conceptualization, Methodology, Supervision, Writing—original draft. HL Conceptualization, Funding acquisition, Validation. LG Writing—original draft, Methodology. WC Methodology, Software. JL Conceptualization, Methodology. ZC Supervision, Validation. XG Funding acquisition, Resources. HW Software, Validation.

Declarations

Competing interest

The authors declare that they do not have any commercial or associative interest that represents a conflict of interest in connection with the work submitted.

List of symbols

x

Cartesian coordinate system x-axis

y

Cartesian coordinate system y-axis

z

Cartesian coordinate system z-axis

σ₁

Far-field maximum principal stress around the roadway

σ₂

Far-field intermediate principal stress around the roadway

σ₃

Far-field minimum principal stress around the roadway

β_i

Far-field principal stress azimuth

(i = 1–3)

α_i

Far-field principal stress dip angle (i = 1–3)

L_i

Direction cosine of principal stress axes relative to the Cartesian x-axis (i = 1–3)

M_i

Direction cosine of principal stress axes relative to the Cartesian y-axis (i = 1–3)

N_i

Direction cosine of principal stress axes relative to the Cartesian z-axis (i = 1–3)

C

Direction cosine matrix

C^T

Transpose matrix of C

σ´_x

Far-field normal stress in x-direction

σ´_y

Far-field normal stress in y-direction

σ´_z

Far-field normal stress in z-direction

τ´_xy

Far-field shear stress on xy-plane

τ´_yz

Far-field shear stress on yz-plane

τ´_zx

Far-field shear stress on zx-plane

D

Stress component matrix at any point around the roadway

σ_x

Normal stress in x-direction at any point around the roadway

σ_y

Normal stress in y-direction at any point around the roadway

σ_z

Normal stress in z-direction at any point around the roadway

τ_xy = τ_yx

Shear stress on xy-plane at any point around the roadway

τ_xz = τ_zx

Shear stress on zx-plane at any point around the roadway

τ_yz = τ_zy

Shear stress on yz-plane at any point around the roadway

I₁

First invariant of stress tensor

I₂

Second invariant of stress tensor

I₃

Third invariant of stress tensor

σ₍₁₎

First principal stress at any point around the roadway

σ₍₂₎

Second principal stress at any point around the roadway

σ₍₃₎

Third principal stress at roadway periphery

W

Roadway width

H

Roadway height

ε

Complex plane

w(ε)

Mapping function

R

Real number representing roadway cross-section size

C_i

Real constant (i = 1,3,5)

B

Mechanical constant (B = 1/4 (σ´_x + σ´_z))

B´

Mechanical constant (B´ = 1/2 (σ´_z − σ´´_x))

C´

Mechanical constant (C´ = τ´_zx)

S_i

Coefficients (i = 1,3)

φ(ε) ψ(ε)

Complex potential analytic function

φ₀(ε) ψ₀(ε)

Analytic function outside the unit circle in the ε-plane

E

Elastic modulus

v

Poisson’s ratio

ε_y

Axial strain of roadway

c

Cohesion

φ

Internal friction angle

θ_σ

Lode angle of stress

J₂

Second invariant of deviatoric stress

J₃

Third invariant of deviatoric stress

h

Rock/coal seam thickness

K

Bulk modulus

G

Shear modulus

ρ

Density

ψ

Dilatancy angle

σ_t

Tensile strength

S_zz

Numerical simulation z-direction stress boundary value

S_xx

Numerical simulation x-direction stress boundary value

S_yy

Numerical simulation y-direction stress boundary value

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

References

Cai, MF; Guo, QF; Li, Y; Du, ZF; Liu, JH. In situ stress measurement and its application in the 10th Mine of Pingdingshan coal group. Chin J Eng; 2013; 35, 11 pp. 1399-1406. [DOI: https://dx.doi.org/10.13374/j.issn1001-053x.2013.11.003]

Cui, XZ; Li, XY; Du, YF; Bao, ZH; Zhang, XN; Hao, JW; Hu, YY. Macro-micro numerical analysis of granular materials considering principal stress rotation based on DEM simulation of dynamic hollow cylinder test. Constr Build Mater; 2024; 412, [DOI: https://dx.doi.org/10.1016/j.conbuildmat.2023.134818] 134814.

Fan, QY. Mechanical analysis for choosing optimal orientation of underground opening. J China Coal Soc; 1990; 15, 3 pp. 62-71.

Kang, HP; Wu, L; Gao, FQ; Lv, HW; Li, JZ. Field study on the load transfer mechanics associated with longwall coal retreat mining. Int J Rock Mech Min Sci; 2019; 124, [DOI: https://dx.doi.org/10.1016/j.ijrmms.2019.104141] 104141.

Liu, YF. The derivation of the direction of the third principal stress and the azimuth (or dip angle) of the second principal stress. Mech Eng; 1991; 6, pp. 60-61.

Liu, HT; Han, ZJ; Guo, XF; Zhou, GD; Wei, SJ; Han, Z; Chen, XG; Cheng, WC. Approximate solution of plastic zone boundary of surrounding rock of circular roadway considering axial stress. Coal Sci Technol; 2023; 51, 10 pp. 12-23. [DOI: https://dx.doi.org/10.13199/j.cnki.cst.2022-1405]

Liu, HT; Chen, ZH; Guo, XF; Han, ZJ; Liu, QY; Han, Z; Zhang, HK. Research on stope equivalent hole model and rotation law of principal stress. J China Coal Soc; 2023; 48, 10 pp. 3646-3658. [DOI: https://dx.doi.org/10.13225/j.cnki.jccs.2022.1594]

Liu, HT; Han, ZJ; Han, Z; Guo, LF; Han, XY; Liu, QY; Chen, ZH; Zhang, RG. Study on borehole stability under different spatial angles in three-dimensional stress field. J China Univ MinTechnol; 2024; 53, 5 pp. 925-942. [DOI: https://dx.doi.org/10.13247/j.cnki.jcumt.20240072]

Liu, HD; Wang, GB; Yang, CH; Zhang, JY; Chen, SW. Rock strength weakening subject to principal stress rotation: experimental and numerical investigations. J Rock Mech Geotech Eng; 2024; 16, 9 pp. 3544-3557. [DOI: https://dx.doi.org/10.1016/j.jrmge.2024.04.004]

Liu, HT; Han, ZJ; Liu, QY; Chen, ZH; Han, Z; Zhang, HK. Low sensitivity research and engineering application of roadway butterfly failure strength criterion. Rock Soil Mech; 2024; 45, 01 pp. 117-130. [DOI: https://dx.doi.org/10.16285/j.rsm.2023.0107]

Liu, HT; Han, ZJ; Chen, ZH; Guo, LF; Liang, JL; Han, Z; Liu, QY; Ji, Y. Study on the relationship between deviatoric stress distribution and failure characteristics of roadway considering axial stress. J China Coal Soc; 2024; 2024, pp. 1-18. [DOI: https://dx.doi.org/10.13225/j.cnki.jccs.2024.0275]

Lv, AZ; Zhang, LQ. Complex variable function method for mechanical analysis of underground tunnel; 2007; Beijing, Science Press:

Ma, YC; Lu, AZ; Cai, H; Zeng, XT. Analytical solution for determining the plastic zones around two unequal circular tunnels. Tunn Undergr Space Technol; 2022; 120, [DOI: https://dx.doi.org/10.1016/j.tust.2021.104267] 104267.

Ma, ZY; Zuo, JP; Zhu, F; Liu, HY; Xu, CY. Non-orthogonal failure behavior of roadway surrounding rock subjected to deep complicated stress. Rock Mech Rock Eng; 2023; 56, pp. 6261-6283. [DOI: https://dx.doi.org/10.1007/s00603-023-03397-x]

Peng, ML; He, MC; Xiao, YM; Cheng, T; Qiao, YF. Asymmetric failure behavior of surrounding rock in the deep roadway: a semi-analytical solution. Eng Fail Anal; 2024; 160, [DOI: https://dx.doi.org/10.1016/j.engfailanal.2024.108075] 108075.

Pu, JY; Yu, QL; Zhao, Y; Li, ZF; Cao, YS; Le, ZH; Yang, ZM; Liu, X. Deformation analysis of a roadway tunnel in soft swelling rock mass based on 3D mobile laser scanning. Rock Mech Rock Eng; 2024; 57, pp. 5177-5192. [DOI: https://dx.doi.org/10.1007/s00603-024-03772-2]

Shan, RL; Bai, HB; Li, YZ; Wu, HT; Sun, P; Zhao, XP; Dou, HY; Wang, HL. Analytical solution of surrounding rock stress and plastic zone of rectangular roadway under non-uniform stress field. Sci Rep; 2024; 14, 27482. [DOI: https://dx.doi.org/10.1038/s41598-024-79221-5]

Song, SX; Wang, P; Yin, ZY; Cheng, Y. Micromechanical modeling of hollow cylinder torsional shear test on sand using discrete element method. J Rock Mech Geotech Eng; 2024; 16, 12 pp. 5193-5208. [DOI: https://dx.doi.org/10.1016/j.jrmge.2024.02.010]

Su, C; Kang, HP; Jiang, PF; Liu, C; Liu, YD; Yi, K. Analysis on mining-induced stress evolution and surrounding rock failure mode of roadway during heading-mining period based on continuous measurement. Chin J Rock Mech Eng; 2024; 43, 9 pp. 2201-2213. [DOI: https://dx.doi.org/10.13722/j.cnki.jrme.2023.1201]

Wang, JC; Wang, ZH; Yang, J; Tang, YS; Meng, LBB; QB, ,. Mining-induced stress rotation and its application in longwall face with large length in kilometer deep coal mine. J China Coal Soc; 2020; 45, 3 pp. 876-888. [DOI: https://dx.doi.org/10.13225/j.cnki.jccs.SJ20.0147]

Wang, JC; Wang, ZH; Li, Y. Longwall top coal caving mechanisms in the fractured thick coal seam. Int J Geomech; 2020; 20, 8 06020017. [DOI: https://dx.doi.org/10.1061/(ASCE)GM.1943-5622.0001722]

Wang, HN; Song, F; Zeng, GS. Research progress on analytical approaches and models for tunnel engineering. Chin Quarterly Mech; 2023; 44, 2 pp. 227-257. [DOI: https://dx.doi.org/10.15959/j.cnki.0254-0053.2023.02.001]

Xu, HS; Li, SJ; Xu, DP; Huang, X; Zheng, MZ; He, JH; Zhao, K. Numerical back analysis method of three-dimensional in situ stress fields considering complex surface topography and variable collinearity. Int J Rock Mech Min Sci; 2023; 170, [DOI: https://dx.doi.org/10.1016/j.ijrmms.2023.105474] 105474.

Zhang, YJ; Xu, M; Liu, SJ; Liu, F; Wang, QS. Rate-dependent constitutive modelling blasting crack initiation and propagation in rock masses. Int J Coal Sci Technol; 2023; 10, 1 83. [DOI: https://dx.doi.org/10.1007/s40789-023-00633-1]

Zhao, GP; Yang, SL. Analytical solutions for rock stress around square tunnels using complex variable theory. Int J Rock Mech Min Sci; 2015; 80, pp. 302-307. [DOI: https://dx.doi.org/10.1016/j.ijrmms.2015.09.018]

Zuo, JP; Ma, ZY; Xu, CY; Zhan, SF; Liu, HY. Mechanism of principal stress rotation and deformation failure behavior induced by excavation in roadways. J Rock Mech Geotech Eng; 2024; 16, 11 pp. 4605-4624. [DOI: https://dx.doi.org/10.1016/J.JRMGE.2024.05.016]