Failure mechanism and fracture aperture characteristics of hard thick main roof based on voussoir beam structure in longwall coal mining

The mechanism and fracture aperture of thick hard roof failure were theoretically analyzed from the perspective of elastoplastic mechanics. A uniform load cantilever beam was established by considering the voussoir beam structure and stress condition. According to the stress inverse principle in solving the plane strain problem, the internal stress component analytic solution was deduced. Considering the lateral thrust and shear force and by applying the Mohr‐Coulomb failure criterion, the plastic zone boundary equation was obtained. Further analysis of the equivalent stress, equivalent strain, elastic energy density, distortion energy density, and strain energy density was studied with blocks of four different lengths, combined with the rotation angles after failure, to analyze the fissure aperture between straight and curved fracture traces. The theoretical analysis reveals that the deviatoric effect plays the lead role in rupture initiation and propagation; the cracks will deflect during failure, and the fracture trace presents a curved shape. The variation of fracture aperture exhibits a tendency to increase first, then decrease, and then increase along the fracture trace. The analytical approach presented in this paper provides a theoretical basis for determining the main seepage channel flow in underground mining.


| INTRODUCTION
Thick hard main roof instability failure is one of the most hazardous and serious dynamic disasters in underground mining. [1][2][3][4][5] The sudden failure exhibits instantaneous, strong impact, accompanied by energy release. In addition, mining-induced fracture aperture is the main channel for aquifer water inrush and gob gas migration. [6][7][8][9] Understanding the failure mechanisms and fracture aperture characteristics plays an important role for establishing a working face early warning system.
In order to explore the failure mechanisms, much research work has been devoted to study thick and hard roof breaks | 341 ZHANG et Al.
using various methods, including theoretical research, [10][11][12] in situ investigation, [13][14][15][16] and physical and simulation analysis. [17][18][19][20][21][22] On the basis of spatial structure evolution and load conditions, analyses typically simplify the unit thickness on the working face middle to a beam model, [23][24][25][26] such as clamped beam, simply supported beam, and voussoir beam; the voussoir beam approach is a kind of structure widely used and was developed rapidly. These studies on the structural model of the main roof overlying strata often model the rock beam as a continuous medium and then analyze the stability and failure of the rock beams. Shabanimashcool & Li 12 studying the effect of horizontal load on the stability of voussoir beam, the horizontal stress improve the stability of the beam against buckling. Zhao et al 27 analyzed the structure instability process of the roof rock beam and point out the stability of the cracked roof rock beam depend on the vertical load and horizontal thrust interaction. Please C.P et al 28 taken the roof layer as a clamped beam and using simple strut and beam theory investigated the fracture of the roof. Most analysis methods were based on static equilibrium, and the internal stress distribution was neglected. Gu et al 29 and Zhang et al 30 illustrate the distribution of deflection, bending moment, and shear force on the roof beam, and without considering the thrust pressure subjected on the end of beam. This approach is based mainly on the relationship between deflection and load distribution to predict the failure region and position, ignoring the influence of internal stress distribution on the failure. However, the distribution characteristics of the internal stress field, especially the stress deviation effect, is the key factor in revealing the failure mechanism and fracture traces, and it is also the basis for further study of fracture aperture seepage.
Because of the main seepage channels caused by the failure of overburdened strata, further study of the fracture features is required in order to determine the diversion capacity. Many researchers have predicted the height of the water-conducting fissures zone. [31][32][33][34][35][36] Chi et al 37 analyzed the water pressure by considering the mining height and aquifer position, and take the ration of water pressure change and mining height as the failure criterion to describe the development of fracture. Bai and Tu 38 point out the mining-induced fracture as main fluid channels and summarized the current study of mining-induced fracture near face regions. Li and Bai 39 studied the activation and redevelopment process of the fracture and the results revealed the mechanism of the secondary failure of the fracture caused by caving. Chen et al 40 considering the water pressure and backfill support and using the Winkler elastic foundation beam model calculated the height of the water-flow fracture zone, which decreases with the increases of the mining levels. Li et al 41 adopted the geo-electrical response, physical model, and numerical method to study the evolution of mining-induced fractured zone. However, most of the studies focused on the macroscopic fracture zone evolution and the height of the water-conducting fissures, but there is little research on the main seepage channels inside the fracture zone. In order to describe the opening degree of the fracture diversion capacity, Xu et al 42 showed that the key strata have a controlling effect on fracture development. Huang et al 43 proposed the concept of fracture aperture through degree and showed that the degree of fracture aperture is linearly reduced along the crack tip at the critical opening displacement. Although past studies made significant progress, there are still questions that require further interpretation, including what mechanism is responsible for the rupture evolution during crack propagation and the variation of crack opening degree along fracture trace theoretically.
In this study, to reveal the mechanism of thick hard roof breaking from the internal stress components and analyze the fracture aperture combined with the rotation angle. We established a cantilever beam model considering the hinged structure and deduced the analytical solution of the internal stress component equation according to the solution of the plane strain problem in elastic mechanics. Moreover, combined with the Mohr-Coulomb failure criterion, the plastic zone boundary equation also obtained. We further analyzed the crack propagation from the deriatoric effect and then discussed the opening degree of fracture combined with the rotation angle.

VOUSSOIR BEAM STRUCTURE
It is generally acknowledged that movements of overlying rock can be divided into three distinct zones: [44][45][46][47] the caved zone, fracture zone, and continuous deformation zone. The evolution of the overburdened structure is directly related to fracture of the rock strata; in particular, thick and hard roof formations play a key role in the development and evolution of the three zones. For their unique features of high strength, great thickness, and good integrity, thick and hard roof formations suspended over a large space can bear the weight of the overlying weak rock formations and present a block structure after fracture. 48 To describe this phenomenon, the voussoir beam model was first proposed by Evans 49 to study the stability of a fractured beam and then further developed into the compressive arch and multi-jointed voussoir beam. [50][51][52] Recently, the voussoir beam has been widely applied. 53 In China, Qian et al 54 raised a voussoir beam structure based on field measurements and showed that the hard rock was broken and formed an arrangement of articulated rock blocks under the horizontal thrust force, as shown in Figure 1(a); it can be seen that the voussoir beam model forms a large structure. The thick hard roof breaking and movement evolution are from A to C, and the balance of the key blocks A, B, and C near the coal wall in the separation zone has a direct effect on the control of the overlying stratum, as shown in Figure 1 On this basis, Qian subsequently proposed the restriction conditions for sliding failure and crushing failure. 53 According to the solution of Qian et al, 54 the lateral thrust is a necessary condition for the formation of a complete voussoir beam; the vertical shear force and horizontal thrust acting on block A under uniform load are given by.
where T and Q are horizontal thrust and vertical shear force, respectively.
According to Qian's theory, the sinking amount and rotation angle of block A after fracture can be calculated by.
Where θ is the rotation angle; W is the sinking amount of the rock mass, which is related to the mining height M h , the total thickness of the immediate roof Σh, and the buckling coefficient K p , which is obtained experimentally from the laboratory testing and field statistics; L is the periodic weighting step distance obtained by field mine pressure monitoring; h is the thickness of the thick hard rock layer; and q is the overlying load and its self-weight.
It should be noted that this study mainly analyzed the breaking of block A. Equation (1) is used to determine the load on block A and provides the stress boundary conditions for internal stress analysis; Equation (2) gives the rotation angle of block A after breaking, which provides a basis for the subsequent calculation of fracture aperture.

| Mechanical model
According to the voussoir beam structure evolution described in the previous section, block A was simplified in this study to a cantilever beam model as shown in Figure 2. The model assumes homogeneous materials of thickness h and length L. In the model, there is a uniformly distributed load at the top of the beam, and no restriction at the bottom of the beam. The righthand end is regarded as a fixed boundary condition before the crack occurs, the left-hand end is subjected to the action of lateral thrust T and shear force Q from block B, and the point of force application is a/2 away from the beam bottom. It should be noted that in the stress components calculation, we shift the action point to the origin point according to the translation principle of force.
In Figure 2, q = γH, where γ is the average unit weight of the rock, and H is the total thickness of the weak rock formations above the thick hard roof. If a is the contact length of the extrusion area after rotation, the distance of the hinged point is a/2 away from the bottom edge of the beam; therefore, the bending movement expression is. To describe the failure mechanism according to the stress inverse method from elastic mechanics, three sets of equations (equilibrium equations, compatibility equation, and continuity equations) are combined to obtain the internal stress solution, enabling analysis of the stress deviator effect and energy density. On this basis, the solution of the plastic zone boundary was deduced by applying the Mohr-Coulomb criteria, 55-58 as described in the following discussion.

| Internal stress component solution
On the basis of the plane strain problem in elasticity mechanics, the internal stress component solution can be transformed into the solution of a biharmonic equation. 59 From the governing equation and boundary conditions of the beam model, the trial stress function V can be formulated as.
Substituting Equation (5) into Equation (4), the stress function can be written as.
The corresponding stress components are derived as follows: On the basis of the stress boundary conditions on the beam at y = ±h to determine the unknown coefficients, The following is thus derived: According to Equation (9), the expressions of the obtainable coefficients A, B, C, D, and F can be solved, along with the relation between E and G: Considering the stress boundary conditions at both ends of the cantilever beam, and combining with Equation (11), other unknown coefficients can be obtained. The stress boundary condition of the left end of the cantilever beam (x = 0) is: y 4 + Hy 3 + Ky 2 + x Ey 3 + Fy 2 + Gy Substituting Equation (7) into Equation (12) and combining with Equation (11) yields the other coefficient expressions: Substituting all the coefficients into Equation (7) obtains the analytical formula of the stress components: When substituting the stress components of Equation (14) into the right-end boundary conditions to verify the correctness, if all the integral equations are established, then the solution is correct. For the right stress boundary condition of the cantilever beam (x = l), We substituted the stress component Equation (14) into the right-end boundary condition of Equation (15) for verification and found that the integral equations are all established, so the solution of this model is correct.
To further predict the failure characteristics from the perspective of stress deviator and energy distribution, it is necessary to decompose the element stress into hydrostatic stress and deviator stress. The hydrostatic stress and deviator stress represent a combination of shear effect; from the mesoscopic mechanism, the hydrostatic component causes only volume change, and the stress deviation causes the shape change. 60,61 The strain of the material deformation depends on the stress deviation; therefore, the deviatoric component plays an important role in the yield failure. The second stress tensor invariant and strain tensor invariant are calculated as follows: Stress invariant: Strain invariant: Through these two invariants, a synthesized index representing the partial deviation effect can be obtained, such as equivalent stress and equivalent strain. By applying the same principle, splitting the strain energy density into volume-change energy density and deforming energy density yields the following: Strain energy density: Volume-change energy density: Distortion energy density: Thus, using the above Equations, the equivalent stress, equivalent strain, strain energy density, volume-changing energy density, and distortion energy density can be easily calculated to further analyze the failure mechanism.

| Solution of the plastic zone boundary
To study the plastic zone boundary, on the basis of the stress deviation and energy distribution, a failure criterion that accords with the rock features must be applied. For the positive and negative signs in elasticity and rock mechanics are different, it should be noted that we assumed compression stress as positive and tensile stress as negative. The expression of maximum principal stress σ 1 and minimum principal stress σ 3 can be obtained by: The Mohr-Coulomb shear failure yield criterion can be described as.
where φ and c are the internal friction angle and cohesion strength, respectively. We employ function f (x, y) to predict the yield failure point: (23), the function f (x, y) can be written as.

By substituting Equation (21) into Equation
Consequently, recalling Equation (14), we have. where: Equation (25) is the expression of plastic zone boundary based on the Mohr-Coulomb shear failure criterion; the corresponding solution of f (x,y) = 0 is the rupture line of the main roof, where K 1 , K 2 , and K 3 are coefficients related to loads.

| RESULTS AND DISCUSSION
In this section, we discuss the analysis of the equivalent stress, equivalent strain, distortion energy density, bulk energy density, and strain energy density on the basis of the stress components obtained in the Section 3.2 to predict the possible fracture mode. In the fracture mode analysis, taking the yield failure boundary as the fracture trace and combined with the rotation angle, the crack opening degree along the fracture trace was analyzed and compared with the present existing hypothesis. The geological and mechanical coefficients used for calculation are summarized in Table 1.
In order to draw the main roof rupture line using Equation (25), according to the geological data in Table 1,  the periodic weighting step on the workface was determined as L = 20 m, and the bulking coefficient K p was set as 1.25 on the basis of a large number of field measurements and empirical estimates.

F I G U R E 3
With the main roof fracture as an example, the calculation procedure is as follows:  (25), and by calculating and solving the equations with MATLAB software, the rupture trace corresponding to the right side of the block A can be obtained with respect to periodic fracture.

| Synthesized analysis of the deviation effect
We calculated the deviation effect and energy distribution with four different block lengths (12 m, 20 m, 28 m, and 36 m). With the length of 28 m as an example, the equivalent stress, equivalent strain, elastic energy density, and distortion energy density are shown in Figure 3.
On the whole, the equivalent stress, equivalent strain, elastic energy density, and distortion energy density all show asymmetrical distribution in the beam. From the free end to the fixed end, the variations become more pronounced, and in the upper and lower corner points at the fixed-end section, the parameters all reach maximum values. The maximum value of equivalent stress, equivalent strain, elastic energy density, and distortion energy density are 109 MPa, 3.97, 90 kJ/m 3 , and 218 kJ/m 3 , respectively. It is also clear that the distortion energy density is much larger than the elastic energy density, which indicates that the shape change is more obvious than the volume change.
In addition, the reason for the asymmetric distribution is related more to the stress deviator, because the upper part of the beam is in a state of tensile stress and tensile strain occurs; on the contrary, the corresponding low part is in a state of extrusion and compression strain occurs. So, the direction of stress and shear deformation has a deflection in the cross-section. The curves of equivalent stress, equivalent strain, elastic energy density, and distortion energy density at the fixed end are shown in Figure 4. The strain energy density distribution with four different beam lengths is shown in Figure 5. With increasing beam length, the strain energy density also increases sharply; the energy rapidly transfers to the fixed end and accumulates at the upper and lower corners at the fixed end. The maximum strain energy density appears at the upper right corner of the beam; the strain energy density of beams with length 12 m, 20 m, 28 m, and 36 m are 22 kJ/m 3 , 48 kJ/m 3 , 292 kJ/m 3 , and 804 kJ/m 3 , respectively. In addition, it can be seen that the distribution of strain energy density on each beam presents nonuniformity, and with the increasing of beam length the deviatoric effect more obvious. This is due to the external loading environment and constraint conditions, the maximum bending moment at the fixed end, and the fixed-end section undergo the transformation from tensile stress to compressive stress the minimum strain energy density at the critical point.

| Contrast analysis of fracture aperture
For the analysis discussed in this section, we took the plastic zone boundary as the rock block fracture trace combined with the rotation angle, analyzed the crack opening degree with four different lengths, and compared the results with the current hypothesis. The comparative analysis diagram with four beam lengths is shown in Figure 6. The difference between them is obvious; the current general assumption is a 90° straight line, but the result of our calculated rupture trace is a curve with an inflection point, which is consistent with the deviation results of the previous analysis. With the increase in beam length, the curvature increases, which reflects the crack propagation and extension angle; this is more helpful for predicting the fracture angle of the goaf. After the rock block breaks, there will be a rotary movement. Different block lengths correspond to different rotation angles; according to the Equation (2) mentioned above, the rotation angles can be calculated, for rock block lengths of 12 m, 20 m, 28 m, and 36 m, and the corresponding rotation angles are 23°, 14°, 10°, and 8°. With increasing block length, the rotation angle decreases, and when the rock block rotation movements become stable, the fracture will open and separate. We define the relative distance as the crack aperture; the crack aperture varies along the fracture trace. The fracture aperture corresponding to the fracture trace of Figure 6 is shown in Figure 7, which quantitatively compares and analyzes the differences of the curves and straight lines for the four different block lengths in terms of the fracture aperture. As shown, the fracture aperture along the current straight-line hypothesis has an increasing trend, whereas the fracture aperture along the curves calculated in this study exhibit a trend of increasing first, then decreasing, and then increasing again. This trend is more obvious with increasing rock block length.
When the rock block length is 12 m, the difference in the fracture aperture along fracture trace between the curves and straight line is small, only a bit less than the straight line in the range of 4 m < x < 6 m. As the block length increases to 20 m and 28 m, the range of decline gradually increases, and the lowest point is reached at x = 6 m. Meanwhile, the range greater than the straight line appears in the range of 10 m < x < 14 m. It is notable that when the block length reaches 36 m, the fracture aperture declines to zero in the range of 6 m < x < 8 m, which shows that the crack in this range closed because of the fracture rock block rotation. This is more realistic than the current straight-line hypothesis, and further analysis of the difference shows that after a certain angle of curve fracture rotation, the opening width in the inflection point starts to decrease, which means that the flow of seepage changes along the fracture trace. This result is highly significant for the study of breaking crack seepage based on the fracture aperture.
Previously, the fracture trace was generally assumed to be a straight line, which was inconsistent with field observations. The determination of the fracture trace morphology can provide a basis for predicting the main diversion channel of mining-induced fracture fissure and further provides a foundation for the study of water or gas migration in the future.

| CONCLUSIONS
In this study, the failure mechanism and fracture fissure apertures of thick hard roofs were studied from the perspective of elastoplastic mechanics, according to the voussoir beam hinged structure and stress conditions. A model of a uniformly loaded cantilever beam was established, in which the horizontal thrust and vertical shear force are taken into account, and the internal stress component analytic solution was obtained using the stress inverse method. The boundary equation of the plastic zone was obtained by combining the Mohr-Coulomb failure criterion. From the plastic mechanics perspective, the equivalent stress, equivalent strain, elastic energy density, distortion energy density, and strain energy density were analyzed with four blocks of different lengths, and for the rotation angle, the variation of the straight line and curved traces along the fracture trace were compared and analyzed.
The contour distribution of the equivalent stress, equivalent strain, elastic strain density, distortion energy density, and energy strain density are all curved, with the most dramatic changes at the fixed end, which indicates that the cracks will deflect during rupture initiation and propagation. This result is consistent with the plastic boundary shape obtained through the Mohr-Coulomb failure criteria. The failure will occur at the fixed end, and the rupture process will follow the curve from the top to the bottom at the fixed end. Through comparative analysis, it can be seen that with increasing beam length, the rotation angle decreases and the corresponding variation in the fissure aperture along the fracture trace also changes. The variation degree presents a tendency to increase first, then decrease, and then increase along the rupture line; a greater block length exhibits a more obvious trend, which is helpful for determining the main seepage channel of the goaf and plays an important role for further predicting the inrush flow of the broken fracture.