Directionality of butterfly leaves and nonuniform deformation mechanism in gob-side entry driving roadway


 Existing rock pressure theory, with the pressure-plastic zone as the core, cannot reveal the nonuniform large deformation mechanism of the surrounding rock in a deep roadway. To address this problem, the 11030 gob-side entry driving roadway of the Zhaogu No. 2 Mine is used as an example, in conjunction with elastoplastic mechanics theoretical analysis, to study the butterfly plastic region and the direction of the butterfly leaf in the homogeneous circular roadway. In addition, the finite-element method was used to study the environmental characteristics of the stress field of the surrounding rock before the excavation of the gob-side entry driving roadway and to analyze the effect of the principal stress on the butterfly plastic region of the surrounding rock. The results showed that, under the combined action of high ground stress and disturbed stress, the surrounding rock of the deep front edge of the gob in the central area of the main stress could form a large range for the butterfly leaf plastic area. The deflection angle of the four butterfly leaves in the plastic region was the same as that of the main stress direction. The butterfly plate plastic zone that formed in the deep surrounding rock of the gob-side entry driving roadway led to not only extension, degradation and disappearance to different degrees, but also deflection under the action of the main stress. This led to large deformation of the rock strata in the butterfly leaf plastic region, ultimately causing uneven deformation and failure of the 11030 roadway.


INTRODUCTION
In recent years, as deep mining technology with high resource recovery rate has been promoted, a driving roadway along a gob as a retaining method has been widely employed in China's coal mines [1][2][3]. Although driving along the gob at an extreme depth can improve the resource recovery rate of a coal pillar, there are some outstanding problems regarding the aspects of strata properties and the surrounding rock control. Existing related rock pressure theory, with the supporting pressure-plastic zone as the core, only focuses on the influence of the supporting pressure on the morphological characteristics of the surrounding rock plastic zone while ignoring the influence of the stress direction [4]. As a result, it cannot sufficiently reveal the nonuniform large deformation mechanism. The nonuniform deformation and failure mechanism of a roadway-surrounding rock, as well as the deformation and failure characteristics of the roadwaysurrounding rock, are closely related to the distribution characteristics of the plastic zone. Therefore, studying the influence of the stress direction on the distribution characteristics of the plastic zone is essential to better understand the nonuniform deformation and failure mechanism of the surrounding rock in deep roadways.
At present, according to the classical theory of the surrounding rock plastic region of a roadway, the roadway-surrounding rock is primarily round or oval in shape [5,6]. The classical Fenner and Kastner formulas are mainly applicable to the surrounding rock stress field of a shallow roadway. Based on this, many academics have researched surrounding rock plastic regions. For example, Jia et al. [7] found that the magnitude and direction of the main stress in the roadway-surrounding rock changed constantly, indicating that the connection phenomenon was common in the process of plastic region amplification, and the composite roof failure order was shallow plastic failure, deep penetration plastic failure and middle layer failure. Yuan et al. [8] found that the effect of mining stress on a plastic region was more significant in a deep roadway than a shallow roadway. Even though the change of the mining effect was small, it could cause serious plastic failure of the surrounding rock and sudden instability of the roadway. Under different confining pressures, the plastic region of the surrounding rock of a circular roadway presented different shapes in the nonuniform stress field. Guo et al. [9] also provided an index that reflected the form features of the plastic region. Based on the generalized plane strain problem and equal strength theory, Zheng et al. [10] derived the equations of the stress and plastic fracture range in the elastic-plastic zone of a roadway floor for three-dimensional in situ stresses. Based on the study of the evolution law of a plastic zone in a mine roadway, Ma et al. [11] proposed that the stress direction determined the maximum plastic failure depth of the surrounding rock and controlled the range of the potential caving zone. Kang et al. [12] proposed that because of large deep ground stresses, the gate driving along the next gob was more strongly affected by mining, the plastic zone of the coal increased, the coal support, drum out, roof subsidence and cyclotron increased, and the floor heave was serious. However, these results mainly focused on the effect of stress on the plastic region features of the surrounding rock, and the relationship between the plastic zone distribution features and stress direction remains understudied.
Therefore, based on the production and geological conditions of the 11030 gob-side entry driving roadway of the Zhaogu No. 2 Mine in Henan Province in China, the distribution features of the plastic zone of the surrounding rock along the gob were studied. Additionally, the influence of the main stress direction on the distribution features of the plastic region was studied via theoretical analysis and finite-element method (FEM).

CA SE STUDY 2.1 Geological and mining conditions
The average thickness of the main No. 2 1 coal seam in the Zhaogu No. 2 Mine of the Henan Energy Chemical Group was 6.1 m with an average dip angle of ∼5°. The buried depth of the site was ∼700 m. The roof lithology of the No. 2 1 coal seam was primarily mudstone, sandy mudstone and sandstone. The immediate roof thickness range was between 1 and 6.5 m, and the floor was dominated by sandy mudstone and siltstone. The geological column is shown in Fig. 1.
To improve the resource recovery rate, an 8-m narrow coal pillar was set in the gob of the 11011 working face in the southwest area of the mine to excavate the 11030 roadway (Fig. 2). The roadway was driven along the coal seam roof, and the driving section was rectangular in shape, ∼4.8 m wide and 3.3 m high.
The Zhaogu No. 2 mining area has previously experienced many instances of tectonic stress. The coal seam is deep, and under the joint action of high crustal stress and tectonic stress, the underground roadway-surrounding rock is soft and broken, resulting in large deformation of the surrounding rock. During the driving of the 11030 roadway, the section shrinkage was acute, and the surrounding rock showed obvious heterogeneity and large deformation features (Fig. 3). However, the 11030 roadway-surrounding rock produced a nonuniform shape of the plastic zone, which led to large deformation of the surrounding rock [13]. Therefore, this study takes the 11030 roadway as the research object, studies the influence of the stress direction on the distribution characteristics of the butterfly plastic zone in the surrounding rock of the roadway, explores the formation mechanism of the "butterfly" directionality of the butterfly plastic zone and provides a theoretical basis for revealing the nonuniform deformation and failure mechanism of the 11030 roadway.

Mechanical parameter experiment
To investigate the physical and mechanical parameters of the surrounding rocks of the 11030 roadway, and to provide basic data for theoretical calculations and numerical simulations, cores were taken from the roof and floor of the No. 2 1 coal seam. The uniaxial compression, conventional triaxial compression and Brazilian splitting rock mechanics were tested by the HCT microcomputer-controlled electro-hydraulic servo pressure testing machine (Fig. 4) in the coal mine and in the Key Laboratory of Western Mines of the Ministry of Education in China. The tensile strength, compressive strength, elastic modulus, deformation modulus and Poisson's ratio of 2 1 coal and different rocks from roof and floor were measured. According to the experimental requirements, the number of samples in each state is generally ≥3, a total of 45 rock samples, and the samples are loaded at a loading rate of 1.0 kN/s until they fail. The core, part of the test piece and failure form of the test piece are shown in Fig. 5; the test results are shown in Table 1.
The mechanical parameters of sandstone and floor limestone of the No. 2 coal roof are relatively high, while the mechanical parameters of the 2 1 coal are relatively low, indicating that sandstone and floor limestone have better strength and integrity, and are not easy to deform and fail. Moreover, the gap between the mechanical parameters of the mudstone in the top plate and the sandy mudstone in the bottom plate is small, but their tensile strengths are low, indicating that the original rock layer structure is developed, resulting in poor strength and integrity. Under a high-stress environment, it is easy to produce deformation and damage.

METHODOLOGY
Because the 11030 roadway was driven along the 11011 gob with 8-m coal pillars, the stress redistribution on the coal pillar side caused by the 11011 working face could be regarded as the in situ stress before the excavation of the 11030 roadway. Relevant theoretical research and engineering practices [14,15] have shown that for the stress field environment of a deep roadway, the main stress of the surrounding rock can comprehensively reflect the distribution states of the vertical stress, horizontal stress and shear stress, which can reveal the mechanical behaviors of the deformation and failure of the deep rock mass. Thus, the FEM was employed for the study of the surrounding rock failure on one side of the 11011 gob and its influence on the distribution law of the surrounding rock plastic area. Furthermore, the distribution law of the main stress field of the surrounding rock that was highly relevant to the failure of the surrounding rock on the other side of the gob was obtained. This distribution law provided a mechanical basis for analyzing the distribution features of the plastic zone of the surrounding rock in the 11030 roadway. Based on this, the main stress vector effect on the distribution characteristics of the butterfly-shaped plastic region and the orientation of the butterfly-shaped leaves was studied using elastic theory and FEM.

Numerical simulation of the principal stress distribution characteristics
The principal stress distribution characteristics were numerically studied using FLAC 3D numerical simulation in the FEM. A FLAC 3D model was established with the geological conditions of the Zhaogu No. 2 coal mine. The FLAC 3D numerical simulation model was 250 m wide, 350 m long and 32 m high (Fig. 6).
The left and right frontiers of the model were restricted by horizontal displacement, while the bottom frontier was restricted by vertical displacement. In order to approximate the 700 m load of the overburden, the average bulk density of the overburden is 25 000 N/m 3 , and a force of 17.5 MPa is applied to the top of the model. According to the measured geostress data of adjacent mines in the mine area, the Ini command was used to apply the initial geostress field, with a vertical stress of 15.37 MPa (σ 3 ) and horizontal stresses of 29.42 MPa (σ 1 ) and 16.70 MPa (σ 2 ). The Moore-Coulomb criterion was employed for the constitutive relationship of the surrounding rock. The physical and mechanical parameters of each rock stratum are shown in Table 1. When the working face was excavated for 120 m along the x-axis and pushed for 100 m along the x-axis, the main stress distribution characteristics of the gob side of the working face were as shown in Fig. 6. At the same time, based on the rock mechanics test results discussed in Section 2.2, Roclab software was used to calculate the numerical simulation mechanical parameters of the No. 2 coal mine and its roof and floor, Roclab is based on the Hoek-Brown strength criterion, and the physical and mechanical parameters of the test piece can be directly converted into real rock mass parameters through Roclab software, as shown in Table 2.

Theoretical analysis of the butterfly plastic zone and the butterfly leaf direction
Due to the limitations of the existing methods of mathematical mechanics, the plastic zone of a rectangular roadway cannot be accurately obtained by theoretical calculations. Therefore, to determine the feasibility of theoretical research, we performed a theoretical analysis based on the generally accepted theoretical model of a plastic region of the homogeneous surrounding rock for a circular roadway. Although the model is different from the actual situation of the 11030 roadway studied in this research, its qualitative conclusion could be compared with the distribution features of the surrounding rock plastic area of the 11030 roadway, which provided a basis for revealing the large and uniform deformation mechanism of the surrounding rock of the deep gob-side entry.
For the actual stress field of the roadway, the stress field direction of the surrounding rock rotates (no longer being absolutely horizontal or vertical) under the effect of the structural stress or mining stress, especially the gob-side entry driving zone, because various structures that are formed by the surrounding rock of the overlying roof can result in the overlying rock to gob-side entry developing a certain tilt, which can cause the direction of a roadway-surrounding rock stress field to change. This will in turn lead to the redistribution of a plastic region in the deep roadwaysurrounding rock. According to the calculation formula of the surrounding rock stress of a uniform circular roadway under the condition of an inhomogeneous stress field [16] combined with elastic-plastic mechanics theory for the condition of an inhomogeneous stress field, a point stress state in the roadway could be obtained after the clockwise rotation α. The surrounding pressure of the roadway was as follows:   Figure 6 Numerical model of the principal stress distribution characteristics. Note: ϕ is the angle of internal friction, C is the cohesion, γ is the bulk density, G is the bulk modulus, K is the shear modulus and R m is the tensile strength.
where λ is the ratio of the minimum principal stress, P 1 (MPa), to the minimum principal stress, P 3 (MPa), in the roadway area stress field; a (m) represents the road radius, r (m); θ (°) represents the polar coordinates of the surrounding rock; and α (°) represents the rotation angle of the minimum main stress and the minimum main stress (relative to the horizontal-vertical coordinate system, in which clockwise is positive). The Moore-Coulomb criterion could be expressed in the form of the stress in polar coordinates: × sin ϕ cos ϕC − C 2 cos 2 ϕ = 0.
In the formula, ϕ (°) is the internal friction angle of the surrounding rock and C (MPa) is the cohesion of the surrounding rock. After substituting Eq. (1) into Eq. (2), when the roadway confining pressure rotated at α, the boundary equation of the surrounding rock plastic region was When the specific calculation was made according to Eq. (3), the calculation method shown in Fig. 7 could be used.

Numerical simulation of the directivity of the butterfly leaves in the butterfly-shaped plastic region
To study the directional characteristics of the butterfly blade in the butterfly plastic region of the surrounding rock, and for comparison with the above theoretical calculation and analysis, the FLAC 3D model in the FEM was used to establish the 11030 roadway. Because this study mainly focused on the directional characteristics of the butterfly blade in the roadway section surrounding rock, the distribution characteristics of the plastic zone along the roadway axial surrounding rock were not considered. Therefore, the FLAC 3D model was built according to the two-dimensional plane model, i.e. 50 m along the x-axis, 40 m along the z-axis and 0.5 m along the y-axis (Fig. 8). Furthermore, to simulate the geological conditions of the coal seam as accurately as possible, the dip angle of the numerical model was set to 5°. The approximate principal stress at the 8-m position on one side of the gob was used as the original stress field (σ 1 = 80 MPa, σ 2 = σ 3 = 30 MPa), which was applied by the model in the simulation. Additionally, the distribution characteristics of the surrounding rock plastic region in the original stress field were simulated when the principal stress direction deflection (clockwise) angle increased from 0°to 90°(each increase of 10°). When the principal stress was deflected in the direction of 0°, σ 1 was in the x-direction, σ 2 was in the y-direction and σ 3 was in the z-direction. Moreover, the front, back, left and right frontiers of the model were horizontal displacement constraints, and the upper and lower boundaries were vertical displacement constraints. According to the formula of the principal stress direction in elastic-plastic mechanics, the principal stresses in different directions were decomposed and applied via the Ini command. For the constitutive relationship of the surrounding rock, we adopted the Moore-Coulomb criterion. The physical and mechanical parameters of each rock are reported in Table 3. Figure 9 illustrates the main stress distribution on the gob side of the 11011 working face, which was obtained using the research method described in Section 3.1. The stress concentration degree of the surrounding rock was the highest within 6-18 m of the gob side, and the minimum main stress reached 70-80 MPa, ∼2.38-2.72 times the minimum main stress of the in situ stress. The minimum main stress concentration also occurred in the surrounding rock within 4-8 m of the gob side. The minimum main stress value was ∼30-35 MPa, which was ∼1.95-2.28 times the minimum principal stress of the ground stress. The main stress ratio of the surrounding rock in the range of 0-2 m at the edge of the gob was the largest, at 4-6. The main stress ratio of the surrounding rock in the 4-10 m range of the gob edge was also high, and the ratio of the maximum stress to the minimum stress was ∼2-4. However, the main stress ratio of the surrounding rock was <2 when the working face was >10 m. The main stress ratio of the surrounding rock within 10 m around the edge of the gob was >2.

The butterfly-shaped plastic zone of the roadway-surrounding rocks
Because there was an 8-m coal pillar left along the 11011 gob, the stress field distribution at the gob side of the 11011 working face was the stress field before the unearthing of the 11030 working face. For the feasibility of the theoretical calculation, the maximum main stress and the minimum main stress at a position 8 m from one side of the 11011 working face were taken as the maximum main stress, P 1 , and minimum main stress, P 3 , in the theoretical calculation model of the homogeneous surrounding rock of the circular roadway, respectively. Then, the maximum main stress, P 1 , was 80 MPa (horizontal direction), and the minimum principal stress, P 3 , was 35 MPa (vertical direction). Taking the radius of the circular roadway as the circumscribed radius of the 11030 roadway, the surrounding rock property was simplified to 2 1 coal. According to Eq. (3) and the calculation method suggested in Fig. 7, the calculation program of the plastic region radius was developed, and the parameter assignment in Table 4 was substituted in order to calculate the plastic region radius of the roadway-surrounding rocks. Figure 10 shows the butterfly-shaped distribution diagram of the plastic region that was drawn from the calculated plastic region radius of the roadway-surrounding rocks. The maximum  radius of the plastic region boundary of the roadwaysurrounding rocks was located in the four quadrants of the coordinate system, and the minimum radius was located in the coordinate axis of the coordinate system. This indicates that the plastic zone frontier is represented by the shape of the four quadrants protruding in the concave shape of the coordinate axis, which is defined as the butterfly shape. The parts of the butterfly wing that protrude in the four quadrants are called butterfly leaves.

The directionality of the butterfly leaves in the plastic zone
Under the actual pressure environment of the roadway in the coal mine, the surrounding rock stresses in the directions of the maximum main stress and the minimum main stress were not absolutely horizontal or vertical. In the surrounding rock that was affected by mining, the superposition of the ground stress field and the mining stress field caused the surrounding rock in the uneven stress field, which caused the main stress direction in the surrounding rock stress field to deflect. This showed that when the main stress direction of the uneven stress field changed, the distribution characteristics of the butterfly plastic region of the surrounding rock changed accordingly.

Theoretical analysis results
Using the theoretical calculation method described in Section 4.2, we represented α with the values of 0°, 10°, 20°, 30°, 40°, 60°, 70°, 80°and 90°, and substituted in the calculation parameters in Table 3 in order to calculate and draw the butterfly boundary distribution curve of the surrounding rock plastic region (Fig. 11). The value θ in the polar coordinates of the maximum radius position point on the plastic region frontier represented the direction angle of the butterfly blade in the plastic zone. As Fig. 11 shows, in any direction of stress, the boundary of the surrounding rock plastic region was butterfly-shaped, and the radius of the butterfly plastic region was unchanged. When the main stress direction was clockwise, the four butterfly leaves were all clockwise. The deflection angle of the four butterfly leaves was the same as the corresponding main stress direction. The azimuth difference of the four butterfly leaves after the deflection angle was also a multiple of 90°. This showed that when the direction of the main stress changed, the surrounding rock still formed a butterfly-shaped plastic region. The butterfly plastic area was deflected at the same angle as the main stress direction.

Numerical simulation results
Because the cross section of the 11 030 roadway is rectangular and the surrounding rock is nonuniform, the distribution characteristics of the butterfly plastic region and the uniformity area in the circular roadway could be significantly different. Therefore, the numerical simulation method mentioned above was employed to obtain the plastic distribution of the surrounding  rock for different main stress directions (Fig. 12). The direction angle of the plastic area of the butterfly blade on the roof was defined as the angle between the point line of the maximum damage depth of the plastic area on one side of the roof and the same side of the butterfly blade and the roof surface (Fig. 12). The deflection angle of the plastic area of the butterfly blade referred to the difference between the direction angle of the plastic area of the butterfly blade when the main stress direction did not deviate and the direction angle of the plastic area of the butterfly blade when the main stress direction deviated to a certain extent.
As Fig. 12 shows, when affected by the properties of the surrounding rock, stress direction and section shape, the large stress concentration at the sharp rectangular roadway corner, a butterfly-shaped plastic area (similar to the butterfly leaf shape shown in Fig. 10) was formed in the surrounding rock at the sharp roadway corner. When the main stress direction was not deflected, a large range of butterfly leaf plastic areas formed in the rocks at the sharp corners of both sides of the roadway roof, and the two butterfly leaves overlapped each other in the roof rock. In addition, when the deflection angle of the main stress direction increased from 10°to 70°, the plastic area of the butterfly leaf in the rock stratum on the right side of the roof moved to the right. Additionally, the plastic area of the butterfly leaf in the rock stratum on the roof left side moved to the right, and the plastic area of the butterfly leaf gradually degenerated until it disappeared (the deflection angle of the main stress direction was 30°). When the main stress deflection angle was 80°, the butterfly leaf plastic area in the rock stratum on the roof right side continued to move to the right. The butterfly leaf plastic zone in the rock stratum on the left side of the roof formed again, and Table 4 Parameter assignment for computing the radius of the plastic region.  the butterfly leaf plastic area in the rock stratum on the left side of the floor essentially disappeared. When the deflection angle of the principal stress was 90°, the plastic zone of the butterfly leaves in the rock layer on the right side of the roof continued to deflect to the right, but the butterfly leaves degenerated to some extent. Moreover, the plastic zone of the butterfly leaves in the rock layer to the left of the roof expanded to some extent. In order to further explore the impact of the main stress direction on the butterfly plastic area direction, the direction angle of the butterfly plastic area of the right roof rock was used to draw the change curve of the deflection angle of the butterfly plastic area with the main stress direction, as shown in Fig. 13. As suggested in Fig. 13, the deflection angle of the main stress direction was approximately positive and linear with the deflection angle of the plastic area of the butterfly blade. With the direction of the main stress deflection, the butterfly blade plastic area was deflected. However, in contrast to the circular roadway with a uniform surrounding rock, the deflection angle of the main stress direction of the rectangular roadway was larger than that of the butterfly leaf plastic area.

Large deformation mechanism of the surrounding rock
Because of the 8-m pillar driving along the transportation roadway, the roadway-surrounding rock exhibited a high stress concentration in the stress field under the action of the high ground stress and the stress superposition effect of the mining in the adjacent working face. The proportion of the maximum principal stress to the minimum principal stress was >2, which led to the formation of the butterfly plastic zone. Existing research [17] shows that the main stress deflection direction of the 8-m pillar was ∼70°and the butterfly leaf in the plastic area has directivity, leading the 11030 roadway to form similarly to the plastic area shown in Fig. 12h. A large range of the butterfly leaf plastic areas was formed in the rock stratum on the roadway roof right side, and the butterfly leaf plastic area was formed in the surrounding rock on the upper right side. Due to the large plastic failure range and depth in the butterfly leaf plastic region, the surround- ing rock deformation in the butterfly leaf plastic region was larger than that in the other parts of roadway, resulting in the uneven deformation and failure features of the surrounding rock in the 11030 working surface, as illustrated in Fig. 14.

CONCLUSIONS
To reveal the nonuniform large deformation mechanism of surrounding rocks in the 11030 roadway, the 11030 roadway of the Zhaogu No. 2 Mine was used as a case study, and theoretical calculations and numerical simulations were carried out to study the directivity of the butterfly leaves in the plastic region of the deep roadway-surrounding rock. We reached the following conclusions: (1) Under the superimposition of high ground stress and mining stress, obvious stress concentrations occurred in the surrounding rock within 6-10 m of the gob side of the 11011 working surface, and the proportion of the maximum main stress to the minimum main stress was ∼2-4. When the 11030 haulage roadway was excavated 8 m away from the gob, the roadway was in the stress concentration region, and the surrounding rock could easily form a large-scale plastic area. (2) When the main stress direction of the homogeneous surrounding rock of the circular roadway changed, the deflection angle of the four butterfly leaves in the butterfly plastic region was the same as that in the main stress direction, and the azimuth difference of the four butterfly leaves after the deflection was a multiple of 90°. (3) In the 11030 roadway-surrounding rock, a butterflyshaped plastic zone was formed that was similar to the butterfly leaf shape in the homogeneous surrounding rock of a circular roadway. However, the butterfly leaf in the butterfly plastic area expanded, degenerated or disappeared to different degrees. As the main stress direction of the surrounding rock was deflected, the butterfly leaves in the plastic zone were deflected, but the deflection angle of the main stress direction was larger than that of the butterfly leaves in the plastic region. (4) Because of the directivity of butterfly leaves in the butterfly leaf plastic area of the surrounding rock, a large area of the butterfly leaf plastic area was formed in the rock strata on the roof right side and the surrounding rock on the upper right side. This caused the deformed butterfly leaf plastic region in the surrounding rock of this part to be larger than the regions of other parts. This eventually led to the nonuniform large deformation failure of the 11030 roadway-surrounding rock.
According to existing rock pressure theory, the influence of the stress direction on the distribution characteristics of the surrounding rock plastic zone is not generally considered. Thus, the influence of the principal stress direction on the distribution characteristics of the surrounding rock plastic zone can be better illustrated using a theoretical basis for revealing the mechanism of nonuniformity deformation and failure. However, because this article takes actual engineering as a research case and is limited by the geological conditions of the study area, the conclusions obtained have certain limitations. Therefore, future work will include conducting an in-depth study on the influence of stress direction on the distribution characteristics of surrounding rock plastic zones under various geological conditions.