Elastic–plastic criterion solution of deep roadway surrounding rock based on intermediate principal stress and Drucker–Prager criterion

Abstract

Excavation-induced rock failure and deformation near an underground opening boundary is closely associated with the intermediate principal stress, strain softening and rock mass dilation. By combining theoretical analysis and numerical simulation to explore the mechanical evolution of the roadway surrounding rock during excavation. The elastic–plastic criterion solutions for surrounding rock stress, displacement, and the plastic zone of a circular roadway were deduced by including the intermediate principal stress, strain softening and rock mass dilation, based on the Drucker–Prager criterion and the nonassociated flow rule. Furthermore, ABAQUS finite element software was utilized for numerical simulations to scrutinize the influence of mechanical parameters on stress redistribution, surface displacement, and plastic range. The results of the numerical simulations verify the reliability of the theoretical analysis and affirm the high consistency of the results with the theoretical solution. The research findings indicate that the strength of surrounding rock is significantly influenced by the intermediate principal stress, and the deformation and failure of rock mass are primarily impacted by rock mass dilation, while strain softening mainly affects the thickness of the fractured zone. Moreover, the study reveals that enhancing the residual strength of the surrounding rock can improve its resistance to deformation and failure. Furthermore, the excavation of the roadway is observed to increase the original rock stress in the surrounding rock, but increasing the ground support strength effectively controls the deformation and failure of surrounding rock. Ultimately, the research outcomes offer valuable references for engineering calculations and ground support design of surrounding rock in deep roadways.

INTRODUCTION

With the demand for coal resources and mining degree, the technical problems of mining stability and support design of deep surrounding rock are gradually emerging.^1–3 There will be strong stress disturbance in rock after excavating the roadway, destressing and stress concentration occur in the surrounding rock, and the rock mass around the roadway shows strong post-peak strain softening and expansion features, resulting in huge plastic deformation and rupture in the surrounding rock mass.⁴ Elastic–plastic criterion solution of the surrounding rock can effectively explore the influence mechanism of stability. According to previous studies that elastic-plastic solutions of roadway surrounding rock could be obtained by models.^5–9

In geotechnical engineering research, most scholars used Mohr–Coulomb strength criteria or Hock-Brown strength criteria for analysis.^10-13 Because the influence of intermediate principal stress on the deformation and failure of roadway surrounding rock is not taken into account, there is a discrepancy between the proposed results and the engineering practical results. This shows that there is a great correlation between rock strength and intermediate principal stress. Therefore, this influence factor must be considered when analyzing the failure and deformation of the roadway surrounding rock.¹⁴ There is currently a widespread application of strength criteria considering the influence of intermediate principal stresses within the study of roadways surrounding rock and related rock engineering.^15-17 These applications mainly include the Drucker−Prager strength criterion, Matsuoka–Nakai strength criterion, and the Unified Strength Theory.^11,18 Based on the Unified Strength Theory, Yu et al.¹⁹ proposed that intermediate principal stress could be represented by a parameter b and the other two principal stresses in the plastic zone. However, this assumption lacks a strong theoretical basis. Zhang et al.²⁰ applied the Unified Strength Theory to study the effect of intermediate principal stress on the range of the fractured zone and deformation of the roadway surrounding rock. However, they did not provide a method for calculating intermediate principal stress. The Matsuoka–Nakai criterion considered the effect of intermediate principal stress but did not consider the internal cohesion of geotechnical materials on the strength of surrounding rock.^21-23 Drucker–Prager criterion introduces the influence of intermediate principal stress, thereby providing a more practical solution to this kind of problem.²⁴ Zhang et al.^25,26 based on the Drucker–Prager criterion, the rock strength increased by about 30% when considering the effect of intermediate principal stress. Atsushi²⁷ verified the importance of intermediate principal stress in the evolution of plastic volumetric strain by comparing the volume strain rates of potential functions with different yield criteria, and also indicated that the Drucker–Prager criterion is more feasible in the study of surrounding rock stability in a deep large chamber. The effect of intermediate principal stress on the dynamic mechanical properties and fracture mode of rock was investigated by a Numerical Servo Triaxial Hopkinson Bar (NSTHB) system.²⁸ The results showed that the intermediate principal stress significantly affects the dynamic failure mode of rock, and the strength and elastic modulus increase with the increase of the intermediate principal stress in a certain range. Haimson and Chang²⁹ obtained a true triaxial strength criterion for the rock by testing granite, and the results showed that the minimum least horizontal stress strength is enhanced by about 50% over that determined in the conventional triaxial test after considering the intermediate principal stress, and the elastic range and crack inclination of the rock were significantly increased. It can be seen that the intermediate principal stress has a significant effect on the strength and mechanical response of rock mass.

Some researchers^30,31 divided the roadway surrounding rock into an elastic zone and plastic softening zone, and obtained the elastic-plastic solution through the Drucker–Prager criterion, but ignored the highly excavation damaged zone (HDZ) with poor mechanical properties of the rock mass. This area is positioned near the edge of the excavation area, where macroscale fracturing, splitting or spalling may occur and the surrounding rock is seriously damaged.^32,33 The rock strength within HDZ is relatively low, making it prone to safety accidents under strong external stress disturbances.³⁴ Therefore, HDZ cannot be ignored in the process of elastic–plastic analysis of the roadway surrounding rock.

In excavating deep caverns, the process is influenced by not only intermediate principal stress, but also by the superimposed effects of post-peak strain softening, dilation angle, and other rock properties.^35-41 Yu et al.⁴² investigated the combined impact of dilation angle and intermediate principal stress on the excavation of deep roadway surrounding rock, and introduced the coefficient b to quantify the effect of intermediate principal stress. The analytical expression of an ideal elastic-plastic material was derived by adjusting the values of the coefficient in the Drucker–Prager criterion to accurately match different strength criteria. The calculation results confirmed the influence of rock mass dilation and intermediate principal stress on rock strength. However, it is important to note that the consideration of strain softening was not accounted for in the analysis. Based on the nonassociated flow rule and the plane strain assumption. Liu and Yu³⁰ accurately matched the Drucker–Prager strength criterion to the Mohr–Coulomb strength criterion for theoretical analysis of deep circular roadways. The ideal elastic–plastic stress and displacement analytical formulas under the impact of intermediate principal stress and dilation angle were derived. However, there are still challenges in exploring and analyzing the stage after the stress of surrounding rock reaches the peak. Researchers have conducted many discussions on the variation trend of friction angle and cohesion of rock after the peak process. Hajiabdolmajid et al.⁴³ propose a model of cohesion weakening and internal friction angle strengthening (CWFS) based on plastic strain, indicating that cohesion strength and friction strength of rock mass cannot be mobilized simultaneously. Some researchers^44,45 studied the simulation effects of various softening models on deformation and failure of roadway surrounding rock, and showed that CWFS models can better reflect the failure characteristics.

From the above studies, it can be seen that the intermediate principal stress, rock mass dilation, and the change process of strength parameters in the post-peak stage will all have an impact on the stability of the roadway surrounding rock, and only by fully considering these factors can a more reasonable theoretical solution be derived. Therefore, the present research divides three zones of roadway surrounding rock, namely elastic zone, plastic softening zone and fractured zone, analyzing the effect of intermediate principal stress on the roadway surrounding rock using the Drucker–Prager criterion, and introducing the dilation model to characterize the mechanism of dilation angle and CWFS to describe the change process of strength parameters in the post-peak stage. On the basis of the nonassociated flow rule and geotechnical engineering theory, the theoretical equations of stress, radius of plastic softening zone and fractured zone and surface displacement of roadway are derived, also giving a reasonable expression of intermediate principal stress. The research results provide some reference for the engineering calculation and the support design of surrounding rock in deep roadway.

THEORETICAL ANALYSIS MODEL Mechanical model of circular roadway

After excavation of the surrounding rock that failure occurs in the fractured zone, also rock strength is the lowest and close to the residual strength values. However, there were no large deformation characteristics in plastic softening and elastic zones, with the rock strength being between the residual value and peak value, which still has a high supporting capacity.

The mechanical model of the roadway is established in Figure 1. The section of the roadway is set as the center of the circle, and rock is defined as a homogeneous, isotropic, and continuous material. Assume that the roadway depth is deep enough and the working face length is long enough, and the rock mass is in a uniform primary rock stress field. The excavation process of circular roadway can be simplified to plane strain conditions. Assuming that the excavation radius of roadway is R₀, the original rock stress is p₀, support stress is p_i, the radius of the fractured zone is R_s, the radius of the plastic softening zone is R_p, and the radius of the elastic zone is R_e, respectively.

Figure 1. Mechanical model of circular roadway.

σ_θ, σ_z, and σ_r are the tangential, axial, and radial stresses of the roadway surrounding rock respectively, representing the maximum, intermediate, and minimum principal stresses, namely σ₁, σ₂, and σ₃. The relationship between them is satisfied with σ₁ > σ₂ > σ₃. In addition, ε_θ, ε_z, and ε_r are tangential, axial, and radial strains of roadway surrounding rock, respectively.

Drucker–Prager criterion under plane strain condition

Researchers commonly utilize the Mohr-Coulomb yield criterion to present the strength properties of geotechnical materials in engineering. However, the criterion does not account for the influence of the intermediate principal stress. As a result, researchers have developed numerous modifications to address this limitation.^46–49 The most common approach is to use Drucker–Prager criterion to match the Mohr–Coulomb criterion,⁵⁰ resolving the issue. Thus, the expression of Drucker–Prager strength criterion is derived. [Image Omitted. See PDF]where I₁ is the first invariant of the stress tensor, J₂ is the second shear stress invariant tensor, and α, k are constants related to the internal friction angle φ and cohesion c of geotechnical materials. The expressions for I₁ and J₂ are: [Image Omitted. See PDF]

In engineering, the Lode parameter was often introduced to represent the relationship between three principal stresses. The expression of μ_σ is: [Image Omitted. See PDF]

On the basis of the Equation (3) that expression of σ₂ is solved as $<math ><mrow><mrow><msub><mi>\unicode{x003C3}</mi><mn>2</mn></msub><mo>\unicode{x0003D}</mo><mfrac><mrow><msub><mi>\unicode{x003C3}</mi><mn>1</mn></msub><mo>\unicode{x0002B}</mo><msub><mi>\unicode{x003C3}</mi><mn>3</mn></msub></mrow><mn>2</mn></mfrac><mo>\unicode{x0002B}</mo><mfrac><mrow><msub><mi>\unicode{x003BC}</mi><mi>\unicode{x003C3}</mi></msub><mrow><mo stretchy="false">(</mo><mrow><msub><mi>\unicode{x003C3}</mi><mn>1</mn></msub><mo>\unicode{x02212}</mo><msub><mi>\unicode{x003C3}</mi><mn>3</mn></msub></mrow><mo stretchy="false">)</mo></mrow></mrow><mn>2</mn></mfrac></mrow></mrow></math>$, and the relationship between σ₂ and I₁, J₂ is obtained by taking σ₂ into Equation (2). The new expression of yield function of Drucker–Prager criterion is obtained by combining Equation (1). [Image Omitted. See PDF]where M represents the parameter linked to α, N represents the parameter linked to k, then $<math ><mrow><mrow><mi>M</mi><mo>\unicode{x0003D}</mo><mfrac><mrow><msqrt><mrow><mn>9</mn><mo>\unicode{x0002B}</mo><mn>3</mn><msubsup><mi>\unicode{x003BC}</mi><mi>\unicode{x003C3}</mi><mn>2</mn></msubsup></mrow></msqrt><mo>\unicode{x02212}</mo><mn>3</mn><mi>\unicode{x003B1}</mi><msub><mi>\unicode{x003BC}</mi><mi>\unicode{x003C3}</mi></msub><mo>\unicode{x0002B}</mo><mn>9</mn><mi>\unicode{x003B1}</mi></mrow><mrow><msqrt><mrow><mn>9</mn><mo>\unicode{x0002B}</mo><mn>3</mn><msubsup><mi>\unicode{x003BC}</mi><mi>\unicode{x003C3}</mi><mn>2</mn></msubsup></mrow></msqrt><mo>\unicode{x02212}</mo><mn>3</mn><mi>\unicode{x003B1}</mi><msub><mi>\unicode{x003BC}</mi><mi>\unicode{x003C3}</mi></msub><mo>\unicode{x02212}</mo><mn>9</mn><mi>\unicode{x003B1}</mi></mrow></mfrac></mrow></mrow></math>$, $<math ><mrow><mrow><mi>N</mi><mo>\unicode{x0003D}</mo><mfrac><mrow><mn>6</mn><mi>k</mi></mrow><mrow><msqrt><mrow><mn>9</mn><mo>\unicode{x0002B}</mo><mn>3</mn><msubsup><mi>\unicode{x003BC}</mi><mi>\unicode{x003C3}</mi><mn>2</mn></msubsup></mrow></msqrt><mo>\unicode{x02212}</mo><mn>3</mn><mi>\unicode{x003B1}</mi><msub><mi>\unicode{x003BC}</mi><mi>\unicode{x003C3}</mi></msub><mo>\unicode{x02212}</mo><mn>9</mn><mi>\unicode{x003B1}</mi></mrow></mfrac></mrow></mrow></math>$.

According to the theory of rock and soil plasticity, introducing α_ψ, then the expression form of plastic potential function can be deduced as follows. [Image Omitted. See PDF]where ψ is dilation angle at strain softening condition of rock mass.

By the nonassociated flow rule, there is: [Image Omitted. See PDF]where S_ij represents the deviant stress tensor of surrounding rock. Under the plane strain condition that there is dε_z = dε_xz = dε_yz = 0, and according to Equation (6), they are. [Image Omitted. See PDF]

By substituting Equation (7) into Equation (1), the expression of the Drucker–Prager strength criterion is converted to another form, F. The expression for F is as follows. [Image Omitted. See PDF]

Equation (8) is one of the expressions that Drucker–Prager criterion with rock mass dilation under plane strain condition.

The Mohr–Coulomb criterion is often expressed in geotechnical research as: [Image Omitted. See PDF]where c stands for internal cohesion of surrounding rock and φ stands for internal friction angle of surrounding rock.

Assuming that the Drucker–Prager criterion matches the Mohr–Coulomb criterion exactly under the plane strain condition, then Equation (8) should be equivalent to Equation (9), then we can obtain the expression for α and k. [Image Omitted. See PDF]

According to the knowledge of rock and soil plastic mechanics,⁵¹ when ψ = φ, we have α_ψ = α. Substituting it into Equation (10), we can get the general expression of α_ψ. [Image Omitted. See PDF]

Such Equations (10) and (11) are the Drucker–Prager criteria for the exact matching of Mohr–Coulomb criteria under plane strain conditions.

Strength parameters modulus

After excavating the roadway, the rock stress is redistributed. In the stage after rock strength reaches the peak, the mechanical properties of materials deteriorate, and then the rock mass has strain softening behavior. The prepeak strength mobilization is mainly controlled by the cohesive strength of rock, and friction becomes effective only after the yielding and failure process has commenced.⁵² The internal friction angle and internal cohesion are not constant values, and the components of both cannot function simultaneously in the process of rock mass failure. To simplify the model, plastic strain methods can be used to estimate cracks and damage within the rock mass. In the context of a CWFS approach, the internal friction angle increases with the accumulation of microcracks, while the internal cohesion has the opposite effect.^53,54 Figure 2 shows the relationship between strength parameters and tangential strain.

Figure 2. The change of strength parameters with plastic shear strain.

The original internal cohesion c₀ and original internal friction angle φ₀ corresponding to the critical plastic strain are initial values. After entering the residual strength, the internal cohesion of the surrounding rock softens to the residual value, and the internal friction angle hardens to the peak value. For simple calculation, the relationship between the strength parameters and the plastic shear strain is simplified as the linear function, and introduce the strength parameter modulus D_c and D_φ to represent the functional relationships.²⁵ According to Figure 2, the expressions for c′ and φ′ are. [Image Omitted. See PDF]where c′ is the internal cohesion of the plastic softening zone, φ′ is the internal friction angle of the plastic softening zone, and ε_θ is the tangential strain of the plastic softening zone, $<math ><mrow><mrow><msubsup><mi>\unicode{x003B5}</mi><mi mathvariant="italic">\unicode{x003B8}</mi><mi mathvariant="italic">s</mi></msubsup></mrow></mrow></math>$ is the tangential strain of the fractured zone. [Image Omitted. See PDF]where $<math ><mrow><mrow><msubsup><mi>\unicode{x003B5}</mi><mi>e</mi><mi>\unicode{x003B8}</mi></msubsup></mrow></mrow></math>$ is the tangential strain of the elastic zone.

Dilation angle model for rock mass

Lade et al.⁵⁵ shows that some uncorrelated flow potential function that causes localization when the strength of the rock material decreases can be applied to the model strain softening behavior, and thus the nonassociated flow rule is generally used to describe the plastic deformation of rocks. During the excavation of roadway, the plastic volume deformation of geotechnical materials occurs due to shear action, and has different expansion coefficients in each zone. Figure 3 shows the relationship between the dilation coefficient and tangential strain.⁵⁶

Figure 3. The change of dilation coefficient with plastic shear strain.

The expansion coefficients of plastic softening and fractured zone are represented by η₁ and η₂, respectively. The flow rules for each zone are: [Image Omitted. See PDF]

The expression of expansion coefficient of plastic softening zone of deep chamber surrounding rock is deduced, then η₁ = (1 + sinψ)/(1 − sinψ).^14,42

The volume strain increment $<math ><mrow><mrow><mo>\unicode{x02206}</mo><msubsup><mi>\unicode{x003B5}</mi><mi>v</mi><mi>s</mi></msubsup></mrow></mrow></math>$ of surrounding rock is composed of tangential strain increment $<math ><mrow><mrow><mo>\unicode{x02206}</mo><msubsup><mi>\unicode{x003B5}</mi><mi>\unicode{x003B8}</mi><mi>s</mi></msubsup></mrow></mrow></math>$ and radial strain increment $<math ><mrow><mrow><mo>\unicode{x02206}</mo><msubsup><mi>\unicode{x003B5}</mi><mi>r</mi><mi>s</mi></msubsup></mrow></mrow></math>$. According to Fu,⁵⁷ the volume strain increment relation of the surrounding rock in the fractured zone can be expressed as follows. [Image Omitted. See PDF]

By combining Equations (14) and (15), we find that η₂ = 1 + ψ. According to analysis results of Farmer's experimental data⁵⁸ that ψ is generally 0.3–0.5.²⁵

ELASTIC–PLASTIC MECHANICS ANALYSIS OF ROCK Fundamental equation

The surrounding rock's internal stress is redistributed after roadway excavation, satisfying the differential equation for stress distribution and the geometric equation for deformation.

Equilibrium differential equation (without physical force): [Image Omitted. See PDF]

Geometrical equation: [Image Omitted. See PDF]where u is the radial displacement and r is the extreme diameter of the roadway.

Constitutive model of plane strain: [Image Omitted. See PDF]where E is elastic modulus and μ is Poisson ratio of rock mass, ε_θ and ε_r are tangential strain and radial strain of rock mass respectively, σ_θ and σ_r are tangential stress and radial stress of rock mass respectively, and p₀ is initial stress of rock mass.

Theoretical analysis of elastic zone

When tunnel is simplified as a plane axisymmetric problem, the stress is independent of not only z-axis but also θ axis. Assuming that the stress component of elastic zone is represented by the Airy stress function, denoted as Φ. As a result, the normal stress components in cylindrical coordinates can be derived using the polar coordinate stress formula. The expressions for the normal stress components can thus be derived as follows. [Image Omitted. See PDF]where r is the distance from any point in the rock mass to the center of the roadway.

The coordination equation under the condition of plane axisymmetry is solved as follows: [Image Omitted. See PDF]

Its general solution is $<math ><mrow><mrow><mi>\unicode{x003A6}</mi><mo>\unicode{x0003D}</mo><mi>A</mi><mspace width=".25em"/><mi>ln</mi><mspace width=".25em"/><mi>r</mi><mo>\unicode{x0002B}</mo><mi>B</mi><msup><mi>r</mi><mn>2</mn></msup><mspace width=".25em"/><mi>ln</mi><mspace width=".25em"/><mi>r</mi><mo>\unicode{x0002B}</mo><mi>C</mi><msup><mi>r</mi><mn>2</mn></msup><mo>\unicode{x0002B}</mo><mi>D</mi></mrow></mrow></math>$. Where A, B, C, and D are the constant values. For a circular hole in a plane tunnel, the single-valued condition for the displacement in the θ direction requires B = 0. Combining Equation (19), there are: [Image Omitted. See PDF]

Roadway excavation minimally affects the elastic zone, and it can remain stable. This is because the rock mass farther from the excavation face is affected only by ground stress, while the rock mass closer to the excavation face is also affected by excavation stress. According to the boundary conditions of the elastic zone, it can be obtained: when r → ∞, $<math ><mrow><mrow><msubsup><mi>\unicode{x003C3}</mi><mi>e</mi><mi>r</mi></msubsup></mrow></mrow></math>$ = $<math ><mrow><mrow><msub><mi>p</mi><mn>0</mn></msub></mrow></mrow></math>$; when r = R_p, $<math ><mrow><mrow><msub><mi>\unicode{x003C3}</mi><mi>\unicode{x003B8}</mi></msub><mo>\unicode{x0002B}</mo><msub><mi>\unicode{x003C3}</mi><mi>r</mi></msub><mo>\unicode{x0003D}</mo><mn>2</mn><msub><mi>p</mi><mn>0</mn></msub></mrow></mrow></math>$. Incorporate the boundary conditions and yield function into Equation (21) to calculate the values of A and C: [Image Omitted. See PDF]

By substituting Equation (22) into Equation (21), we can derive the stress equations of the elastic zone. [Image Omitted. See PDF]

When combining Equations (23) and (18) then radial strain and tangential strain of the elastic zone are calculated as: [Image Omitted. See PDF]where H₀ represents the parameter linked to M and N, then $<math ><mrow><mrow><msub><mi>H</mi><mn>0</mn></msub><mo>\unicode{x0003D}</mo><mrow><mo stretchy="false">(</mo><mrow><mn>1</mn><mo>\unicode{x0002B}</mo><mi>\unicode{x003BC}</mi></mrow><mo stretchy="false">)</mo></mrow><mrow><mo stretchy="false">[</mo><mrow><msub><mi>p</mi><mn>0</mn></msub><mrow><mo stretchy="false">(</mo><mrow><mi>M</mi><mo>\unicode{x02212}</mo><mn>1</mn></mrow><mo stretchy="false">)</mo></mrow><mo>\unicode{x0002B}</mo><mi>N</mi></mrow><mo stretchy="false">]</mo></mrow><mo>/</mo><mrow><mo stretchy="false">(</mo><mrow><mi>M</mi><mo>\unicode{x0002B}</mo><mn>1</mn></mrow><mo stretchy="false">)</mo></mrow><mi>E</mi></mrow></mrow></math>$.

According to the geometric Equation (17) and the elastic strain Equation (24), the elastic displacement can be obtained as: [Image Omitted. See PDF]

Theoretical analysis of plastic softening zone

Theoretical conditions in the plastic softening zone satisfy both the yield function and the equilibrium differential equation. Substitute Equation (4) into Equation (16) and solve it to obtain: [Image Omitted. See PDF]

M′ and N′ are the physical quantities in the plastic softening zone. C₁ is the integration constant. When r = R_p, there is σ_r + σ_θ = 2p_0. Therefore, C₁ can be calculated as $<math ><mrow><mrow><msub><mi>C</mi><mn>1</mn></msub><mo>\unicode{x0003D}</mo><mfrac><mrow><mn>2</mn><mrow><mo stretchy="false">[</mo><mrow><msub><mi>p</mi><mn>0</mn></msub><mrow><mo stretchy="false">(</mo><mrow><msup><mi>M</mi><mo accent="true">\unicode{x02032}</mo></msup><mo>\unicode{x02212}</mo><mn>1</mn></mrow><mo stretchy="false">)</mo></mrow><mo>\unicode{x0002B}</mo><msup><mi>N</mi><mo accent="true">\unicode{x02032}</mo></msup></mrow><mo stretchy="false">]</mo></mrow></mrow><mrow><msubsup><mi>R</mi><mi>p</mi><mrow><msup><mi>M</mi><mo accent="true">\unicode{x02032}</mo></msup><mo>\unicode{x02212}</mo><mn>1</mn></mrow></msubsup><mrow><mo stretchy="false">(</mo><mrow><msup><mi>M</mi><mo accent="true">\unicode{x02032}</mo></msup><mo>\unicode{x0002B}</mo><mn>1</mn></mrow><mo stretchy="false">)</mo></mrow></mrow></mfrac></mrow></mrow></math>$. Bring C₁ into Equation (26) then the stress in plastic softening zone can be calculated. [Image Omitted. See PDF]

When calculating the total strain in the plastic softening zone, the total strain in the plastic zone and the strain on the elastic–plastic boundary should be taken into account, namely $<math ><mrow><mrow><msub><mi>\unicode{x003B5}</mi><mi>r</mi></msub><mo>\unicode{x0003D}</mo><msubsup><mi>\unicode{x003B5}</mi><mi>r</mi><mi>p</mi></msubsup><mo>\unicode{x0002B}</mo><msubsup><mi>\unicode{x003B5}</mi><mi>r</mi><mi>e</mi></msubsup><mo mathvariant="normal">\unicode{x0007C}</mo><mi>r</mi><mo>\unicode{x0003D}</mo><msub><mi>R</mi><mi>p</mi></msub></mrow></mrow></math>$, $<math ><mrow><mrow><msub><mi>\unicode{x003B5}</mi><mi>\unicode{x003B8}</mi></msub><mo>\unicode{x0003D}</mo><msubsup><mi>\unicode{x003B5}</mi><mi>\unicode{x003B8}</mi><mi>p</mi></msubsup><mo>\unicode{x0002B}</mo><msubsup><mi>\unicode{x003B5}</mi><mi>\unicode{x003B8}</mi><mi>e</mi></msubsup><mo mathvariant="normal">\unicode{x0007C}</mo><mi>r</mi><mo>\unicode{x0003D}</mo><msub><mi>R</mi><mi>p</mi></msub></mrow></mrow></math>$. Combined with expansion coefficient of plastic softening zone and Equations (17) and (24), the displacement and strain equations of the plastic softening zone can be obtained as: [Image Omitted. See PDF]

Analysis of the fractured zone

In the fractured zone, the internal friction angle and cohesion are reduced to residual values, while still satisfy the yield function and equilibrium equation. Based on the boundary condition, when r = R₀, σ_r = p_i Combined with Equation (4), the stresses of the fractured zone can be calculated as: [Image Omitted. See PDF]where M* and N* are the physical quantities in the fractured zone, R₀ is excavation radius of roadway surrounding rock. Similarly, on the boundary between the plastic softening zone and fractured zone, there are $<math ><mrow><mrow><msubsup><mi>\unicode{x003B5}</mi><mi>r</mi><mi>s</mi></msubsup><mo>\unicode{x0002B}</mo><msubsup><mi>\unicode{x003B5}</mi><mi>r</mi><mi>p</mi></msubsup><mo mathvariant="normal">\unicode{x0007C}</mo><mi>r</mi><mo>\unicode{x0003D}</mo><msub><mi>R</mi><mi>s</mi></msub></mrow></mrow></math>$, $<math ><mrow><mrow><msub><mi>\unicode{x003B5}</mi><mi>\unicode{x003B8}</mi></msub><mo>\unicode{x0003D}</mo><msubsup><mi>\unicode{x003B5}</mi><mi>\unicode{x003B8}</mi><mi>s</mi></msubsup><mo>\unicode{x0002B}</mo><msubsup><mi>\unicode{x003B5}</mi><mi>\unicode{x003B8}</mi><mi>p</mi></msubsup><mi mathvariant="italic">\unicode{x0007C}r</mi><mo>\unicode{x0003D}</mo><msub><mi>R</mi><mi>s</mi></msub></mrow></mrow></math>$. Combined with the expansion coefficient in the fractured zone and Equation (17) to obtain the displacement of the fractured zone. [Image Omitted. See PDF]

When r = R₀, the tunnel displacement u can be obtained.

According to Equations (17) and (30), the strain equations of fractured zone can be obtained as: [Image Omitted. See PDF]

Analysis of radius of plastic softening and fractured zone

When r = R_s, c′ = c*, φ′ = φ*. From Equations (24) and (28) and the strain condition at the elastic–plastic interface, it can be obtained as: [Image Omitted. See PDF]

Based on Equation (28) and strength parameter modulus D_φ and D_c, the formula of residual internal cohesion and residual internal friction angle in the fractured zone can be derived as: [Image Omitted. See PDF]

At the boundary between the plastic softening and the fractured zone, there is $<math ><mrow><mrow><msubsup><mi>\unicode{x003C3}</mi><mi>r</mi><mi>b</mi></msubsup></mrow></mrow></math>$ = $<math ><mrow><mrow><msubsup><mi>\unicode{x003C3}</mi><mi>r</mi><mi>p</mi></msubsup></mrow></mrow></math>$. Then the radius of plastic softening zone can be calculated by combining Equations (27) and (29). [Image Omitted. See PDF]

Ultimately, the radius of fractured zone can be calculated by combining Equations (33) and (34).

INFLUENCE OF MECHANICAL PARAMETERS ON SURROUNDING ROCK

Based on the above theoretical analysis, the stability of the surrounding rock and the evolution characteristics of partition structure are affected by the intermediate principal stress, strain softening, rock mass dilation, as well as the original rock stress and ground support strength. The excavation radius of the circular roadway R₀ is 3 m with uniform original rock stress p₀ being 15 MPa and support stress p_i being simplified as 1.5 MPa. Table 1 lists the calculation parameters for mechanical model of the roadway surrounding rock.

Table 1 Calculation parameters for mechanical model.

Parameters	Value	Parameters	Value
Original rock stress p0/MPa	15	Support stress pi/MPa	1.5
Elastic modulus E/GPa	2.5	Residual internal cohesion c*/MPa	2
Poisson's ratio μ	0.25	Residual internal friction angle φ*/°	21
Density ρ/kg·m−3	2000	Dilation angle ψ/°	15
Original internal cohesion c0/MPa	5	Expansion coefficient in plastic softening zone η1	1.7
Original internal friction angle φ0/°	10	Expansion coefficient in fractured zone η2	1.4
Lord parameter μσ	0.36

Effect of σ₂ on R_p, R_s, u and stresses of the roadway surrounding rock

The lord parameter μ_σ is used to characterize the magnitude of intermediate principal stress, and its value range is [−1, 1]. Figure 4A shows the effect of the Lode parameter on R_p, R_s, and u. It can be seen that the images are characterized by high ends and low middle, the radius of the plastic softening zone, the fractured zone and the roadway displacement all show the characteristic of decreasing first and then increasing under the action of Lode parameters (i.e., intermediate principal stress). According to Figure 4A, the minimum values of the radius of the plastic softening zone and the fractured zone, and the roadway displacement are obtained when μ_σ is 0.36. While the value of μ_σ is greater than or less than 0.36, the radius of the plastic softening zone and the fractured zone, and the roadway displacement will increase. From Figure 4B, When μ_σ is 0.36, the calculated values of tangential stress and radial stress of the surrounding rock are smaller, and the tangential stress of the surrounding rock is closer to the roadway surface. At the same time, the stress disturbance in the surrounding rock is lower. It shows that the Lode parameter effect has a great impact on surrounding rock strength. When surrounding rock strength can be improved, the radius of the fractured zone, plastic softening zone, and roadway displacement can be reduced. It should be noted, however, that excessively small or large values will reduce the rock strength. In other words, when σ₂ approaches σ₃ or σ₁, the strength of the rock mass will decrease. Consequently, it is crucial to select an appropriate σ₂ in practical engineering applications to mitigate disturbance from excavation on the roadway surrounding rock, thereby allowing the surrounding rock to effectively regulate itself and maintain stability.

Figure 4. Effects of μσ on Rp, Rs, u, and stresses of the roadway surrounding rock. (A) μσ—Rp, Rs and u. (B) μσ—stresses.

Effect of residual strength parameters on R_p, R_s, u, and stresses of the roadway surrounding rock

The relationship between the increase in residual strength values and the corresponding decrease in R_p, R_s, and u can be observed in Figure 5. It is evident that when the residual internal friction angle is constant, a 0.5 MPa increase in residual internal cohesion results in a decrease of 13%–35% in R_p, 17%–38% in R_s, and 18%–20% in u, respectively. Similarly, when the residual internal cohesion is fixed, an increase in residual internal friction angle by 5° leads to a decrease of 15%–32% in R_p, 10%–30% in R_s, and 10%–12% in u, respectively. Figure 5D demonstrates a progressive decrease in R_p, R_s, and u with the increase in residual strength values. Further, the impact of residual strength values on the radius of the softening zone, the radius of the fractured zone, and the surface displacement of the surrounding rock is evident. With the excavation of the roadway, the tangential stress of the surrounding rock increases first and then decreases, while the radial stress increases continuously. From Figure 6, when the continuous reduction of the residual strength values occurs, the peak tangential stress of the surrounding rock is farther from the roadway surface and when there is a greater difference between the magnitudes of tangential stress and radial stress, and the stability and maintenance of the roadway are more adversely affected.

Figure 5. Effects of φ* and c* on Rp, Rs, and u of the roadway surrounding rock. (A) φ* and c*—Rp. (B) φ* and c*—Rs. (C) φ* and c*—u. (D) φ* and c*—Rp, Rs, and u.

Figure 6. Effects of φ* or c* on stresses of the roadway surrounding rock. (A) φ*—stresses. (B) c*—stresses.

Effect of rock mass dilation on R_p, R_s and u of the roadway surrounding rock

According to elastic–plastic theory, the expansion coefficient η₁ is linked solely to R_s and u, while the expansion coefficient η₂ is only related to u. Moreover, both η₁ and η₂ are independent of R_p. As shown in Figure 7A, when η₂ is constant, the dilation angle increases by 5° (i.e., η₁ every 19% increase), which results in R_s increasing by 4% and u increasing by 5%. As shown in Figure 7B, when η₂ is larger, there is a greater variation of u with dilation angle. Moreover, when η₁ is constant, for every 0.5 increase in η₂, u increases by about 2%, and the change of u with η₂ becomes larger as η₁ increases. The rock mass dilation has a greater impact on the displacement of the roadway compared to the fractured zone. The findings indicate that rock mass dilation significantly impacts roadway surface deformation and caving, but has minimal effect on surrounding rock stability.

Figure 7. Effects of expansion coefficient on Rs and u of the roadway surrounding rock. (A) η1 and η2—Rs. (B) η1 and η2—u.

Effect of original rock stress and ground support strength on R_p, R_s, u and stresses of the roadway surrounding rock

From Figure 8A–C, under the action of p₀ or p_i, R_p, R_s and u exhibit the same trend changes. As the original rock stress increases, R_p, R_s, and u also increase, as illustrated in Figure 8D. The stress disturbance of surrounding rock in a higher stress environment is more obvious. As ground support strength increases, R_p, R_s, and u rapidly decrease. It shows that ground support strength plays an important role in controlling the deformation and failure of roadway surrounding rock. The distribution of stress field of surrounding rock is significantly affected by the stress of primary rock, as can be seen in Figure 9A. When the original rock stress increases by 2 MPa each time, the peak tangential stress increases by about 2 MPa, but the peak value farther from the roadway wall increases the range of stress disturbance in the internal rock mass. Moreover, increasing the support effect reduces the difference between tangential and radial stress in the surrounding rock and enhances the stress of the roadway wall, as depicted in Figure 9B.

Figure 8. Effects of p0 and pi on Rp, Rs, and u the roadway surrounding rock. (A) p0 and pi—Rp. (B) p0 and pi—Rs. (C) p0 and pi—u. (D) p0—Rp, Rs, and u.

Figure 9. Effects of p0 and pi on stresses the roadway surrounding rock. (A) p0—stresses. (B) pi—stresses.

NUMERICAL ANALYSIS OF ABAQUS FINITE ELEMENT SOFTWARE Establishment of two-dimensional model

The ABAQUS finite element software offers a range of material constitutive relationships and yield criterion models, such as the Mohr–Coulomb model and the Drucker–Prager model.⁵⁹ These models are particularly effective in accurately representing the internal environment of rock mass and the nonlinear behavior of materials, making the software particularly valuable for addressing nonlinear problems in geotechnical engineering. Therefore, ABAQUS demonstrates strong applicability in the field of geotechnical engineering.⁶⁰

Using ABAQUS software, the process of roadway excavation can be simulated based on the elastic-plastic mechanical conditions of the surrounding rock. This involves simplifying the excavated rock mass into a two-dimensional plane strain model to investigate the space-time evolution characteristics of stress, radius, and displacement of the surrounding rock.⁶¹ Before conducting a numerical simulation, the model is set with a size of 50 m × 40 m and a radius of 3 m. The model is divided into a total of 4872 cells, each set to a quadrilateral structure using a 4-node bilinear plane strain quadrilateral (CPE4I) as the cell type. The model is shown in Figure 10. The simulation parameters depend on the mechanical properties of the rock as listed in Table 1. Furthermore, it is essential to convert the mechanical parameters of the Mohr–Coulomb criterion to those applicable to the Drucker–Prager criterion for the simulation to proceed accurately.⁶²

Figure 10. Plane model of roadway surrounding rock.

Mechanical evolution characteristics of the roadway surrounding rock

Deformation and fracture of the surrounding rock after excavating mostly comes from the roof, bottom and side walls of roadway. Thus, the analysis of rock mass evolution focuses primarily on studying stress distribution, surface displacement, and deformation failure in both vertical and horizontal directions. According to the tangential stress level, the yield failure stage of the surrounding rock can be divided into different zones.

Figures 11, 12, and 13 shows the simulated results of the stress field, surface displacement, and plastic zone of the roadway surrounding rock under residual strength values. The numerical simulation in Figure 11A resulted in a peak tangential stress value of 22.70 MPa, while the elastic-plastic theoretical analysis yielded a value of 22.49 MPa. Similarly, in Figure 12, the magnitude of total displacement (U, Magnitude) of the roadway was found to be 59 mm through numerical simulation, as opposed to the value of 52 mm obtained through elastic–plastic theoretical analysis. The simulation results show that the roof of the roadway has subsided by 58 mm, the two sides have converged by 52 mm, and the bottom floor has bulged by 43 mm. When the lateral pressure coefficient is 1 and the dilation angle is small, the roof of the roadway displacement more than the two sides and the floor. Although the numerical simulation slightly exceeded the theoretical results for the values of tangential stress and displacement, the stress field and displacement field showed consistent developmental characteristics in both sets of data. After excavating the roadway, a significant “layer effect” in space is observed in the distribution of tangential and radial stress, which was shown in Figure 11. The tangential stress value around the roadway starts as the smallest and gradually increases to the maximum value from the inside to the outside before decreasing to the original stress level. Similarly, the radial stress value around the roadway starts as the smallest and gradually increases to the original stress from the inside to the outside. This observation indicates two distinct processes of unloading and loading in stress redistribution.

Figure 11. Distribution of stresses on the roadway surrounding rock. (A) Tangential stress. (B) Radial stress.

Figure 12. The magnitude of total displacement of the roadway surrounding rock.

Figure 13. Equivalent plastic strain of the roadway surrounding rock.

The excavation of a roadway usually results in significant deformation and destruction of the surrounding rock, particularly in the fractured and plastic softening zones. In the study, the equivalent plastic strain (PEEQ) was chosen to indicate the yield state of the surrounding rock, as depicted in Figure 13. In ABAQUS software, PEEQ greater than 0 implies rock mass yielding, and it quantifies the accumulation of plastic strain during excavation and also indicating the range of plastic expansion. In the fracture zone and plastic zone of rock mass, the plastic strain exists, with the highest plastic strain value observed around the roadway, signaling that deformation and failure are most pronounced in this area. Specifically, when the lateral pressure coefficient is 1 under a nonuniform stress field, the plastic strain predominantly accumulates on the two sides of the roadway.

Impact of residual strength parameters on tangential stress distribution, the magnitude of total displacement and plastic distribution ranges of the roadway surrounding rock

Figures 14–17 display the tangential stress of surrounding rock, roadway displacement, and plastic strain distribution under varying residual strength values. When c*=1.5 MPa, φ* = 18°, the stress concentration in the roadway demonstrates an inconspicuous phenomenon, and the numerical simulation results reveal that the maximum tangential stress in the plastic softening zone is 21.67 and 11.90 MPa in the fractured zone. Additionally, the roof of the roadway has subsided by 101 mm, the two sides have converged by 99 mm, and the bottom floor has bulged by 76 mm. While c* = 2.5 MPa, φ* = 26°, the numerical simulation results reveal that the maximum tangential stress in the plastic softening zone is 23.20 and 14.22 MPa in the fractured zone, with the roof of the roadway has subsided by 40 mm, the two sides have converged by 33 mm, and the bottom floor has bulged by 28 mm. The results show that when the residual strength values increase, the stress around the surrounding rock increases continuously, which can improve the bearing capacity. In addition, the surface displacement of surrounding rock decreases continuously, which can effectively improve the deformation and failure effects of surrounding rock. Figure 17 shows that residual strength values significantly affect the plastic strain of surrounding rock. Higher residual strength reduces the likelihood of plastic failure in the roadway and also decreases the disturbances caused by excavation.

Figure 14. Distribution of tangential stresses on the roadway surrounding rock. (A) c* = 1.5 MPa, φ* = 18°. (B) c* = 2 MPa, φ* = 22°. (C) c* = 2.5 MPa, φ* = 26°.

Figure 15. Horizontal displacement of the roadway surrounding rock. (A) c* = 1.5 MPa, φ* = 18°. (B) c* = 2 MPa, φ* = 22°. (C) c* = 2.5 MPa, φ* = 26°.

Figure 16. Vertical displacement of the roadway surrounding rock. (A) c* = 1.5 MPa, φ* = 18°. (B) c* = 2 MPa, φ* = 22°. (C) c* = 2.5 MPa, φ* = 26°.

Figure 17. Equivalent plastic strain of the roadway surrounding rock. (A) c* = 1.5 MPa, φ* = 18°. (B) c* = 2 MPa, φ* = 22°. (C) c* = 2.5 MPa, φ* = 26°.

Influence of rock mass dilation on tangential stress and the magnitude of total displacement of the roadway surrounding rock

The influence of rock mass dilation on the deformation and failure of the roadway surrounding rock was studied for fractured roadway conditions with complex stress. Rock mass typically exhibits rock mass dilation after the peak strain, leading to significant surrounding rock deformation. To investigate this influence, simulation experiments were conducted with dilation angles of 10°, 20°, and 30°. The resulting tangential stresses of the roadway surrounding rock under these different dilation angles are depicted in Figure 18, while the horizontal displacement and vertical displacement of the roadway is presented in Figures 19 and 20. Interestingly, the distribution of the stress field in the surrounding rock did not exhibit significant changes as the dilation angle increased from 10° to 30°. Moreover, the peak value of tangential stress remained consistent across the simulation results. However, as the dilation angle increases, the horizontal and vertical displacements of the roadway also increase significantly. The increase in the dilation angle results in a significant increase in the horizontal and vertical displacements of the roadway. The convergence of the two sides of the roadway has increased from 92 to 130 mm, the subsidence of the roof has increased from 95 to 127 mm, and the bulge of the floor has increased from 71 to 101 mm. As a result, the overall displacement around the roadway has increased by more than 30%. This observation underscores the substantial influence of dilation characteristics on the deformation and failure of surrounding rock in roadways.

Figure 18. Distribution of tangential stresses on the roadway surrounding rock. (A) ψ = 10°. (B) ψ = 20°. (C) ψ = 30°.

Figure 19. Horizontal displacement of the roadway surrounding rock. (A) ψ = 10°. (B) ψ = 20°. (C) ψ = 30°.

Figure 20. Vertical displacement of the roadway surrounding rock. (A) ψ = 10°. (B) ψ = 20°. (C) ψ = 30°.

Influence of original stress and ground support strength on tangential stress and the magnitude of total displacement of the roadway surrounding rock

The distribution of the tangential stress of the surrounding rock and the displacements of roadway under different original rock stresses and ground support strengths are illustrated in Figures 21, 22, and 23. From the numerical simulation results shown in Figures 21A, 22A, and 23A, it is evident that with p₀ = 15 MPa and p_i = 0, the plastic softening zone and fractured zone exhibit maximum tangential stresses of 23.26 and 12.51 MPa, respectively. The roadway's maximum horizontal displacement is 95 mm, while its vertical displacement is 99 mm. The tangential stress is only concentrated at the two sides of the surrounding rock, resulting in a fuzzy demarcation of the boundary of the plastic softening zone. Furthermore, it is observed that under the condition of low original rock stress, the surrounding rock remains relatively unaffected by excavation disturbance, leading to a relatively stable internal zoning structure which possesses the ability to adjust itself to a certain extent. As shown in Figures 21B, 22B, and 23B, when p₀ = 15 MPa and p_i = 2.5 MPa, the range of the plastic region is relatively reduced in space and the surface displacement of the roadway will decrease significantly. Currently, the roadway's maximum horizontal displacement is 37 mm, and its vertical displacement is 43 mm. It shows that the strength of ground support significantly reduces damage around roadway surrounding rock. As shown in Figures 21C, 22C, and 23C, when p₀ = 25 MPa and p_i = 2.5 MPa, numerical simulation yields maximum tangential stresses of 37.09 and 16.12 MPa in the plastic softening and fractured zones, respectively. The roadway's maximum horizontal displacement is 140 mm, while its vertical displacement is 145 mm. These results signify an increase in tangential stress within the rock, along with the gradual expansion of the plastic zone, indicating that the stability of the surrounding rock is significantly impacted by the original rock stress. As the original stress rise, the impact of excavation disturbance on the surrounding rock becomes increasingly pronounced. This leads to a more noticeable phenomenon of tangential stress concentration, along with a continuous expansion of the plastic range. However, enhancing the support effect significantly bolsters the strength of the surrounding rock, effectively controlling the displacement of the roadway surface and thereby ensuring the stability of the surrounding rock.

Figure 21. Distribution of tangential stresses on the roadway surrounding rock. (A) p0 = 15 MPa, pi = 0. (B) p0 = 15 MPa, pi = 2.5 MPa. (C) p0 = 25 MPa, p i = 2.5 MPa.

Figure 22. Horizontal displacement of the roadway surrounding rock. (A) p0 = 15 MPa, pi = 0. (B) p0 = 15 MPa, pi = 2.5 MPa. (C) p0 = 25 MPa, pi = 2.5 MPa.

Figure 23. Vertical displacement of the roadway surrounding rock. (A) p0 = 15 MPa, pi = 0. (B) p0 = 15 MPa, pi = 2.5 MPa. (C) p0 = 25 MPa, pi = 2.5 MPa.

CONCLUSION

This investigation employs the stress–strain model of the three-line segment and is based on the Drucker–Prager criterion and the nonassociated flow rule. By considering the influence of the intermediate principal stress, strain softening, and rock mass dilation, the elastic–plastic criterion solutions for the surrounding rock are determined. The method proposed in this study can be used to study the behavior of completely elastic–plastic materials and rock with cohesion-weakening and frictional-strengthening characteristics. Based on an example, the effects of mechanical parameters on roadway surrounding rock were studied by using the single factor analysis method. Meanwhile, ABAQUS was employed to simulate the spatial evolution of the plastic range of the surrounding rock and to analyze its deformation failure. It is worth noting that the numerical simulation results align with the theoretical analysis results. The main conclusions are as follows:

The intermediate principal stress has a great influence on the range of plastic failure and displacement of the roadway surrounding rock. In a certain range, the greater the intermediate principal stress, the smaller the plastic failure range and displacement, but beyond the critical value, the plastic failure range and displacement increase with the increase of the intermediate principal stress.

Rock mass dilation has a minor impact on the plastic softening zone, but it significantly influences the failure of the roadway. A larger dilation angle is associated with a greater destabilizing effect on roadway stability. Thus enhancing the residual strength values can effectively strengthen the strength of the surrounding rock.

As original rock stress rises, the impact of roadway deformation and failure becomes more apparent. However, it is observed that adequate ground support strength has the potential to significantly improve the instability and deformation of the surrounding rock, while also effectively controlling the surface displacement of the rock.

In the application of the elastic-plastic fractured three-stage model, it is important to note that the model does not encompass the prepeak deformation behavior. It is crucial to acknowledge the impact of mesh size in ABAQUS. Furthermore, there is a need for more engineering examples to refine the selection of strength parameters. Consequently, our forthcoming research phase will involve the utilization of a four-stage stress–strain model to scrutinize the prepeak behavior characteristics of the surrounding rock. It is planned to conduct examples to analyze the effect of intermediate principal stress, strength parameters, and dilation angle on rock mass strength during surrounding rock excavation. This will involve an exploration of the influence of the intermediate principal stress effect on ABAQUS simulation results.

AUTHOR CONTRIBUTIONS

All authors contributed to the study's conception and design. Xinfeng Wang plays a guiding role in conceptualization, method and resources; Youyu Wei mainly takes charge of paper writing, and she made a great contribution to the revision; Tian Jiang and Fuxu Hao main contributions are in data management, resources, writing review and editing; Haofu Xu provided guidance on essay writing. All authors have read and approved the final manuscript.

ACKNOWLEDGMENTS

This work was supported by National Natural Science Foundation of China (grant number 51904266), Excellent youth project of Hunan Provincial Department of Education (grant number 21B0144), Hunan Natural Science Foundation Youth Program (grant number 2023JJ40634) and the Open Fund Project of State Key Laboratory of Mining Response and Disaster Prevention and Control in Deep Coal Mines (grant number SKLMRDPC23KF04).

CONFLICT OF INTEREST STATEMENT

The authors declare no conflict of interest.