Orthogonal Analysis and Numerical Simulation of Rock Mechanics Parameters in Stress Field of Shaft Heading Face

+is paper focuses on improving the blasting effect of the drilling and blasting method in the deep rock mass and solves the problems of blasthole collapse and misfire accident in the process of drilling and blasting construction of heading face. FEM software, ABAQUS, is used to simulate the stress distribution around the blasthole by extending a certain depth in the vertical direction of the shaft heading face. +e sensitivity of different depths, different heading face sizes, and different lithologies on the horizontal stress distribution is analyzed by using a six-factor four-level orthogonal analysis method. +e results show that the change of the radius of the heading face has the most considerable influence on the distance of the distressed zone and the stress concentration zone, followed by the lithology and the excavation depth. Also, the excavation depth has the most significant influence on the peak stress value.+rough the industrial field experiment, the in situ stress of the shaft heading face is tested, and the numerical simulation results are consistent with the field monitoring results. +e results reveal the law of stress distribution near the heading face, which can provide some reference for the design of blasthole depth in the drilling and blasting construction scheme.


Introduction
With the decrement of shallow mineral resources, the development of deep mineral resources cannot be delayed. e excavation of deep shafts has always been a difficult problem to be solved in the event of mineral resources [1]. Currently, the drilling and blasting method is still one of the main techniques of underground engineering excavation in the world. Compared with the use of tunnel boring machines, rock breakers, and other excavation machinery, the drilling and blasting method has low initial investment, low blasting energy consumption, relatively mature technology, and other advantages [2].
Experts and scholars have studied a lot of research on the stress distribution of underground engineering excavation. Jiang et al. [3] based on the established joint strength that can comprehensively consider the tensile and shear characteristics of loess rationally evaluated the stress and displacement calculation of the loess tunnel and obtained the displacement expression of the tube. Luo et al. [4] based on the study of the stress distribution characteristics of the surrounding rock concluded that under the action of high ground stress, the plastic zone of the roadway gradually expands from the periphery to the depth. e pressure of the surrounding rock gradually increases in the plastic area in the deep direction. Qi et al. [5] made some assumptions on the coal in front of the heading face, then analyzed the stress distribution of the coal in front of the heading face from the mechanical point of view, and obtained the formula for calculating the length of the distressed zone to prevent an outburst. Fu et al. [6] pointed out that the plastic zone of the roadway is the main reason for affecting the support of the roadway. e influence law and statistical formula of boundary parameters of the plastic zone at the top and side of the arched roadway are obtained by numerical simulation and the orthogonal test. Yin et al. [7] based on field experiments at the fully mechanized coal mining face studied the relationship between the number of drill cuttings, mine pressure, and gas pressure. It shows that the mine pressure distribution law of the heading face after the pressure relief will form several areas owing to stress redistribution. e width of the pressure relief area is 5-11 m for the plastic area, 11-20 m for the elastic area, and the initial stress area after 20 m. Meng et al. [8] selected a total of 6 different crosssectional roadways including rectangle, trapezoid, straight wall arch, horseshoe, ellipse, and circle for optimization research. Based on FLAC3D simulation, the deformation characteristics and plastic zone distribution of surrounding rock of these six typical roadways after excavation were studied, and the effects of different lateral pressure ratios on them were analyzed. Cardu and Seccatore [9] gave the general relation curve between the tunnel cross section and explosive consumption by analyzing the database of tunnel blasting scheme. Liu et al. [10] used a numerical simulation method to carry out an orthogonal experimental study of the roadway side deep-hole blasting. In addition, many scholars have studied the stability of surrounding rocks in underground engineering through numerical simulation and orthogonal analysis [11][12][13][14][15].
ere have been some new breakthroughs in studying the stability of underground projects in recent years. Xu et al. [16] based on the complex function analyzed the stress distribution law of surrounding rock and the effects of defect location and defect protrusion degree on the stress distribution of the orifice. Wang et al. [17] used the large-scale geomechanical model test system to study the law of stress evolution of the surrounding rock during roadway driving and mining. Shu et al. [18] used FLAC3D to analyze the disturbance stress and the distribution characteristics of the failure zone of the heading face in the coal roadway. ey established the coal roadway heading face of the risk evaluation model. In the process of tunnel excavation, there are apparent mining stress field zoning and permeability field zoning in front of the coal tunnel heading face. e direction of the maximum principal stress controls the development characteristics of the distressed zone and plastic deformation zone. e risk of protrusion of the driving surface increases when the three-dimensional stress is inhomogeneous. Wang et al. [19] believed that roadway excavation would cause rocks near the tunnelling face to flow into the tunnel, which will cause the stress state near the tunnelling face to be perpendicular to the axis of the tunnel which was not a plane strain, but a three-dimensional state. e influence of the internal friction angle on the zonal cracking of the surrounding rock in the gradually excavated arched roadway was studied. Jendryś [20] conducted a numerical simulation study on the surface defects of the shaft wall and analyzed the influence of the width and depth of the defects on the generation of tensile stress in the lining and the change in strength coefficient. Corkum et al. [21] used FLAC3D software to predict the horizontal ground stress on the shaft. Amirata et al. [22] considered the coupling between axial and lateral resistance when modelling large diameter borehole shafts. It is found that the maximum values of shaft bending moment and deflection can be reduced by 12% and 19%, respectively, after coupling analysis.
However, there are few relevant works of literature about the distribution of the stress field in the front heading face of shaft excavation. Because the shaft heading face of drilling blasthole mostly depends on the drilling cutting weight or ground stress test, there is no scientific understanding about the length of blasthole, and the blasting effect will be better. However, it takes time and energy to drill on-site every time, and some areas with complicated geological conditions are not easy to drill. In the construction process of the drilling and blasting method, blasthole collapses and misfire accidents occur from time to time. All these factors are closely related to the stress distribution around the blasthole. erefore, the problem prompted researchers to think about solutions and carry out similar parameter tests and mechanical analysis.

eoretical Model.
According to the theory of elastoplasticity mechanics, the excavated shaft heading face releases the potential energy of the compression deformation of part of the initial rock through the displacement of the rock mass. With the change of the surrounding rock stress distribution and large deformation, the axial stress drop at the heading face of the shaft is 0. Also, with the increase of the research distance, the vertical stress and horizontal stress tend to the initial rock stress.
After several numerical simulation experiments, the stress path of 20 meters downward extension of the simulated model is analyzed. Figure 1 shows the horizontal stress distribution and horizontal stress simulation results. e vertical stress σ v and horizontal stress σ H reach the maximum stable value at a distance of L from the heading face to the second crossing.
rough the numerical simulation results, the vertical stress around the blasthole has been rising and finally tends to the initial rock stress. Because there is no peak stress, and the law is evident, there is not much research here. However, the horizontal stress simulation results show that there are two points of intersection between the horizontal stress and the initial rock stress, which are regarded as the "threshold" of whether the rock mass around the blasthole enters the stress concentration zone or not. Among them, the distressed zone is from the heading face to the first point of intersection, the stress concentration zone is from the first point of intersection to the second point of intersection, and the initial rock stress area is after the second point of intersection. Further, the value range of these two "thresholds" and the magnitude of peak stress need to be also explored by the orthogonal test. e vertical stress distribution and vertical stress simulation results are plotted in Figure 2. e results show that the vertical stress σ v does not converge with the initial rock stress in the stress path extending 20 meters downward.
ere is no stress concentration in the vertical stress, and the development trend is that in the process of extending vertically from the heading face to the underground rock, the pressure increases gradually and finally tends to the initial rock stress. On the other hand, the horizontal stress will climb from the distressed zone at the beginning and then continue to rise to the peak value after the intersection with the initial rock stress and then gradually decrease to the initial rock stress.
For shaft drilling and blasting construction, the in situ stress around the blasthole will increase as the vertical shaft depth increases. During the excavation process, the weight around the working face is redistributed. When the stress distribution on the heading face extends deeper, it can be divided into a distressed zone, stress concentration zone, and initial rock stress area. Peak stress is located in the stress concentration zone, which can be expressed by σ p . When the depth of the blasthole is in the distressed zone, the surrounding rock clips used during blasting are relatively small, and this is the case for shallow hole blasting. When the depth of the blasthole is located in the stress concentration zone, the in situ stress of the blasthole is relatively large, that is, the confining pressure constraint is relatively more significant, and this is the case for deep-hole blasting. When the blasthole is located in the distressed zone, the stress concentration zone, and the initial rock stress area, the blasting effect will be different under the same construction scheme. erefore, when drilling is in the distressed zone, the stress concentration zone, and the initial rock stress area, it is necessary to adjust the blasthole depth, charge structure, charge quantity, and other essential blasting parameters to achieve a better blasting effect.
is paper Advances in Materials Science and Engineering intends to simulate the stress redistribution of the shaft heading face in different lithologies, different section sizes, and different depths by using ABAQUS software and orthogonal analysis. e relationship between stress distribution range and depth, lithology, and excavation section is revealed, which guides blasting design and construction of the shaft and roadway.

Numerical Model.
Because the shaft is axisymmetric, to reduce the amount of calculation, the model is built into a cube of 1/4. erefore, the numerical model is a cube with 25 m * 25 m * 40 m. e Mohr-Coulomb constitutive model is used in this model. e boundary conditions of the model are set to constrain the normal displacements of the bottom and the surrounding areas. Also, the pressure is applied at the top. H is the buried depth of the shaft (m) and g is the gravitational constant. e analysis step is set to balance the ground stress first and then realize the excavation by "killing" some model elements. e rock mass mechanical parameters are isotropic. is numerical simulation scheme is carried out under the condition of no support to avoid the influence on the change of the surrounding rock due to the difference of the support method. e geometric appearance of the numerical model is plotted in Figure 3.

Orthogonal Numerical Simulation
Test. e influence of various factors on horizontal stress of the heading face is studied by the orthogonal analysis method. e factors that affect the stress evolution of heading face include lithology, depth, and excavation size. erefore, this experiment intends to adopt the L 25 (6 4 ) six-factor four-level experimental design.
e design and results of the orthogonal test are shown in Table 1, where φ is the internal friction angle, c is cohesion, μ is Poisson's ratio, E is elastic modulus, h is excavation depth, r is excavation radius, L 1 is the distance of distressed zone, L 2 is the distance of stress concentration zone, and σ p is the peak stress.

Range Analysis.
As the method of studying the influence degree of skillful factors, range analysis has the advantages of intuitionistic results and less calculation. e first step of range calculation is to list the test calculation table and calculate the difference between the maximum mean value and the minimum mean value of each influencing factor [23]. Firstly, we should list test tables and calculate the difference between the maximum average and minimum average of each influencing factor. Secondly, we need to estimate the sensitivity of each element by comparing the difference.
e R value is plotted in Figure 4. From the R value of range in Figure 4, we can see that all factors influence distressed zone distance, stress concentration zone distance, and peak stress, but the degree of influence is different. According to the influence degree of each factor on the distressed zone, the order is as follows: According to the effect degree of each factor on the stress concentration zone distance, the order is as follows: According to the sensitivity of each factor on the peak stress, the order is as follows: e main factors affecting distressed zone distance and stress concentration zone distances are r and μ. e main factors affecting peak stress are h and μ. To study the influence degree of each factor more intuitively, the orthogonal results are analyzed by variance analysis.

Analysis of Variance.
e range analysis of the above data has been carried out, and the order of the influence degree of the leading mechanical parameters on the surrounding rock deformation has been preliminarily analyzed. However, range analysis cannot give an accurate quantitative estimate of the significance of the factors. erefore, with the help of variance analysis, the stress distribution on the heading surface of the shaft is further quantitatively analyzed. e variance analysis is mainly through the significance test to judge the significance level of each factor. e results of the variance analysis for the distance of the distressed zone are shown in Table 2. e F value is the test statistic, which is proportional to the degree of influence. e P value is defined as significant, which is inversely proportional to the degree of importance. e F value of the whole model is 265.132, and the P value is 0.0001. According to the F value and P of r, φ, c, h, and E, it is shown that r has a significant influence on the distressed zone. e F value of μ is relatively less; P value is less than 0.05, which means that this factor has little effect on the distressed zone. e F value of φ, c, h, and E is very small. Compared with r and μ, the P value is far higher than 0.05, which signified that the φ, c, h, and E have little influence on the damaged area. Hence, the main factors affecting the distressed zone are r and Poisson's ratio, followed by other factors. e results of variance analysis for the distance of the stress concentration zone are shown in Table 3. e value of F is a test statistic, which is proportional to the degree of influence. e P value is significant, which is inversely proportional to the degree of importance. e F value of the whole model is 7.010, and the P value is 0.012. According to the F value and P value of different factors, which shows that the r has a significant influence on the distressed zone. e F value of μ is relatively less, and the P value is less than 0.05, which means that this factor has little effect on the distressed zone. e F value of φ, c, h, and E is minimal. Compared with r and μ, the P value is far higher than 0.05, which signified that the φ, c, h, and E have little influence on the damaged area. Hence, the main factors affecting the stress concentration zone are r and Poisson's ratio, followed by other factors. e results of variance analysis for the distance of peak stress are shown in Table 4. e F value is a test statistic, which is proportional to the degree of influence. e value of P is significant, which is inversely proportional to the degree of importance. e F value of the whole model is 216.455, and the P value is 0.0001. According to the F value and P of r, φ, c, h and E, it is shown that h has a significant influence on the distressed zone. e F value of μ and r is relatively less, and the P value is less than 0.05, which means that these factors have little effect on the peak stress. e F value of φ, c, and E is minimal. Compared with r and μ, the P value is far higher than 0.05, which signified that the φ, c, h, and E have little influence on the peak stress. Hence, the main factors affecting the peak stress are h and Poisson's ratio, followed by other factors.

Analysis of Change of Single Variable.
To analyze the effect of different factors on the distressed zone, stress concentration zone, and stress peak value more intuitively, the change of the single factor is simulated. e E of the numerical model is 10 GPa, μ is 0.22, c is 10 MPa, φ is 32°, and density is 2500 kg/m 3 .
It can be seen from Figure 5 that when μ and r of the surrounding rock change, distressed zone distance increases with the rise of these two factors. However, when μ, h, and r of the surrounding rock change, stress concentration zone distance decreases with the increase of these three factors. When μ, h, and r of the surrounding rock change, peak stress increases with the increase of μ and h and decreases with the     e shift in h does not affect distressed zone distance. e E, c, and φ are stable and do not affect distressed zone distance, stress concentration zone distance, and peak stress. It can be seen that φ, c, and E have little effect on the horizontal stress of distressed zone distance, stress concentration zone distance, and peak stress, which is consistent with the results of the orthogonal test.

Influence of the Other Factors
In addition to lithology, excavation depth, excavation radius, the water content of rock mass, joints of rock mass, and lateral pressure ratio of surrounding rock are the factors influencing the stress field of the heading face. Among these factors, the lateral pressure ratio is an essential factor that cannot be ignored, so it is studied separately. When the rock mass gradually enters the high stress state due to the increase of mining depth, the value of horizontal stress is often greater than the value of vertical stress. e influence of horizontal in situ stress on the deformation and failure of surrounding rock near the heading face is more prominent. Due to the increase of horizontal stress, the difficulty in drilling and the charge quantity increase. e lateral pressure ratio reflects the value of horizontal stress. erefore, through the study of lateral pressure ratio, the influence of horizontal stress on the stress field of the heading face is studied indirectly, which can provide more accurate guidance for drilling and blasting construction. e buried depth of the shaft is 500 m, and the lateral pressure ratio is defined as λ and it is set every 0.4 from 0.6 to 1.8, a total of four groups. e E of the numerical model is 10 GPa, μ is 0.22, φ is 32°, c is 10 MPa, and density is 2500 kg/ m 3 . In the same way as the above analysis, the changes of the distressed zone, stress concentration zone, and peak stress are studied. As displayed in Figure 6, when the λ increases from 0.6 to 1.8, the distance of the distressed zone decreases from 3.26 m to 1.24 m, a decrease of 2.63 times, and peak stress increases from 8.51 MPa to 25.51 MPa, an increase of 2.99 times. As λ increases, the distressed zone distance drops and the peak stress gradually rise. e lateral pressure coefficient has a significant effect on the peak stress and stress concentration zone distance. e distance of the stress concentration zone fluctuates when the λ goes up. With the increase of the lateral pressure coefficient, the distance of the stress concentration zone decreases first and then increases.

Engineering Case
Wanfu Coal Mine is located in the southern end of the Juye coalfield. It is a sizeable modern mine constructed by Yanmei Heze Energy Chemical Co., Ltd. in the Juye mining area in Southwest Shandong Province. e buried depth of the coal seam is about 1100 m, and the mining area is located in the Yellow River alluvial plain, with flat terrain, slightly higher in the West and lower in the East. e thickness of mine alluvium is more than 800 meters, which is the largest mine with the largest alluvium thickness among the existing and under construction mines in China. Vertical shaft development is adopted, and the main shaft, auxiliary shaft, and air shaft are arranged in the industrial square. e shaft diameter of the auxiliary shaft is 7.0 m, the depth is 894.368 m, and the surface soil thickness is 751.8 m. ere are siltstone, medium sandstone, mudstone, and fine sandstone passing through the main roadway in Wanfu Coal Mine, and the rock strength is generally small. e test results show that the uniaxial compressive strength of limestone is the largest, about     Advances in Materials Science and Engineering of the overlying rock mass. e in situ stress can be calculated by the following equations: In equations (1) and (2) where the symbols designate the parameters that are P s is the shut-in pressure, P r is the fracture reopening pressure, and P 0 is the pore water pressure. In equation (3), c is the volume-weight and h is the buried depth of measuring points.
rough the field monitoring of auxiliary shaft heading face and comparing the monitoring data of numerical simulation results, the reliability of numerical simulation and orthogonal test results is verified. e hydraulic fracturing method is used to test a set of data for each meter of the auxiliary shaft, and a total of 20 groups are examined.
ree measuring points are selected on the heading surface of the auxiliary shaft starting from −865 meters. e λ is taken as 2 to simulate the Wanfu auxiliary shaft. In practical problems, the maximum horizontal and vertical in situ stresses have an enormous influence on the blasthole. So only the maximum horizontal stress and vertical stress are compared. e stress program and field measured results are compared with the numerical simulation results, as shown in Figure 7. e difference between the measured data of the three groups of maximum horizontal stress may be due to the difference in tectonic stress. e horizontal stress value σ v of the numerical simulation and the data comparison of the three measuring points show that the overall trend is the same, which is first rising and then decreasing. e difference between the three groups of vertical stress measured data and numerical simulation data is less than 10%, which is within the error range. e range of the distressed zone and pressure concentration zone decreases with the increase of Advances in Materials Science and Engineering 9 maximum horizontal stress. is result is consistent with the previous orthogonal analysis result.

Conclusions
is article takes the shaft heading face as the engineering background. It combines the methods of theoretical analysis, field monitoring, numerical simulation, and orthogonal analysis to get the stress distribution of the shaft heading face. e central summary is as follows: (1) e main factors affecting the distressed zone and stress concentration zone are excavation radius and Poisson ratio. e most crucial factor affecting peak stress is excavation depth, followed by the Poisson ratio and excavation radius. However, the other elements, e.g., the internal friction angle, cohesion, and elastic modulus, have little influence on peak stress. (2) With the increment of the lateral pressure ratio, the distressed zone gradually decreases, peak stress increases slowly, and the stress concentration zone first decreases and then increases. e stress concentration zone is the minimum when the lateral pressure ratio is 1.
(3) By comparing the in situ stress monitoring data and numerical simulation data of the Wanfu auxiliary shaft, the horizontal stress value σ H and the vertical stress value σ v have the same overall trend as the data of the three measuring points. e horizontal stress value σ H initially rises and then decreases. e distance of the stress concentration zone is inversely proportional to the horizontal stress. e vertical stress value σ v has been rising since the heading face and finally tends to the initial rock stress, and the influence of the horizontal stress value on the vertical stress is not apparent.

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

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.   10 Advances in Materials Science and Engineering