#### Abstract

In longwall mining, the deformation and destruction of overlying strata always lag behind coal extraction. The overlying strata characteristics at the lateral boundary of the stope can be classified into four categories, i.e., Hard-Soft, Soft-Hard, Hard-Hard, and Soft-Soft. In order to analyze the effect of the above four structures, we adopt viscoelastic theory to the finite element method (FEM) and define the point safety factor to evaluate the rock damage. The accuracy of programming is verified through example verification. A modified viscoelastic-plastic FEM model is applied to analyze the performance of four overburden structures. The numerical computation results show the following: From the rupture of overburden rock to its stabilization, the duration time of four typical structures can be sorted as “Soft-Soft < Hard-Soft < Soft-Hard < Hard-Hard”. The fracture direction and dip angle of each structure vary as well. The fracture zone of the H-S structure is inclined toward the goaf, while that of the S-H structure is inclined to the lateral boundary of the stope. The fracture zone of the H-H structure is also inclined toward the lateral boundary, with a greater angle than the S-H structure, while the fracture zone of the S-S structure is inclined to goaf, with a greater angle than the H-S structure.

#### 1. Introduction

As the coal extracting face pushes forwards, the overlying strata continually crack and expand their damaged range. As the working face is pushed to a sufficient length on distance, it extends the maximized range. Studying the destruction law of overlying rock has great practical significance for the arrangement of coal pillar [1], the determination of mining limit [2, 3], and the selection of mining procedure.

Following the pressure-arch hypothesis proposed by German scholars W. Hack and G. Gilizer in 1928, researchers have proposed many hypothesis and theories for overburden rock structures, i.e., “cantilever beam,” “preformed fracture,” “articulated rock block,” “transfer rock beam,” “masonry beam,” and “key stratum” [4, 5]. On the basis of hypotheses, scientists conducted in-depth researches on specific issues. In order to derive abutment pressure of gobs, Zhang et al. [6] established a calculation model of isolated coal pillars, which is also applicable to assess the burst possibility. To judge whether the key strata can reduce the rotating instability of the main roof, Shi et al. [7] established a criterion on the consideration of stratum rotation. Behera et al. [8] investigated four longwall panels of major Indian coalfields and thus supplemented methodology on longwall mining with various geological situations. Based on field monitoring data, Ju and Xu et al. [9] defined three kinds of the structural model of overburden strata. Kang et al. [10] simulated massive roof collapse by physical modeling. Wang et al. [11] analyzed the mechanical attribute of rock arch in key strata and derived the conditions and evolution of its formation. Wen et al. [12] established a roof structural model for large mining-height stopes and introduced its control criteria. Mondal et al. [13] collected mine-induced microseismic data to calculate the Correlation Fractal Dimension and then analyzed stress distribution.

At the stoppage line of the working face, the stope boundary can be classified as the strike boundary and the lateral boundary. The lateral boundary is on either side of the goaf. For overlying rock layers with different lithological structures, the destruction patterns are different. The destruction law of overlying strata on the lateral boundary is of guide significance for the mining layout. During roadway excavating and supporting, we must take into account the stability of lateral pillar. Though automatically formed entry retaining technology, Manchao He [14–18] developed a no-pillar extracting method. When it comes to the superimposed continuous laminate model, Nong Zhang [19, 20] derived equations about support resistance on conditions of roadway deformation. Tan et al. [21] proposed an artificial composite wall to meet the needs of strong weighting mining conditions. Zhang et al. [22] simplified the roof mechanical model to analyze the entry retaining. By analyzing disasters from overburden separation, Yuan et al. [23] created a physical model and adopted the digital speckle correlation method.

The overlying rock will fall and fill the goaf after coal extraction. Then the layers above the fallout form the hydraulic fracture zone and bending deformation zone. Research methods for the destruction of overlying layers mainly include the structural mechanics method, material mechanics method, and elastic mechanics method. In addition, researchers also used various numerical methods to study the destruction of overlying layers, such as FEM, FLAC, Discrete Element Code. Wang et al. [24] adopted a physical model test to analyze fracture development in overlying strata. Tien et al. [25] presented a discontinuum modeling approach for top coal caving behaviour. By analyzing data in the top coal caving and studying nonassociated elastoplastic strain softening material behaviour, Alehossein and Brett [26] developed a yield and caveability criterion. Alshkane et al. [27] presented a methodology in predicting rock strength and deformability. Li et al. [28] used UDEC to investigate crack evolution characteristics and illustrated the parameter calibration process for various strata structures. Upon the currently used methods, Li et al. [29] explored a physical method to investigate the failure mode of weak rocks. Shabanimashcool and Li [30] presented a novel numerical approach to investigate the loading process of rock bolts.

The FEM and FDM have unique characteristics and problem-solving functions, while they are feasible for continuous media and thus are inadequate in analyzing rock fracture. The discrete element method can be used to investigate block structure, while the analysis results are related to various factors such as joint property. In addition, due to the relative simplicity of the assumed mechanical model, only a limited number of issues can be solved. Many scholars have been working on the retention or the excavation roadway along the goaf, while the theories and background of overburden failure at the lateral boundary are not consummate. In this scope, the viscoelastic-plastic method is applied to solve problems in material and geotechnical engineering. Based on a developed V-P FEM, we aim to discuss the crack form of overburden strata at the lateral boundary.

#### 2. Viscoelastic-Plastic Finite Element Method

Upon the existing three-dimensional elastic-plastic FEM program, the calculation functions for viscoelasticity and plasticity were added. First, the equations and corresponding formulas were derived to update the program. Then, it is subsequently verified by the classical solution of a viscoelastic-plastic (V-P) mechanical problem. The iterative method of viscoelasticity and plasticity is based on the approach proposed by Owen and Hinton [31].

##### 2.1. Constitutive Equation

The one-dimensional V-P model is presented in Figure 1. In the nonlinear continuum problems, the total strain ** ε** can be decomposed into two parts [32, 33], i.e., elastic strain

*ε*_{e}and viscous strain . Therefore, the total strain is as follows:

Assume that the total stress rate and elastic strain rate obey the generalized Hooke’s law [34]:where [**D**] is the elastic matrix. The yield function is introduced to determine whether there is viscoplastic strain:where **F**_{0} is the unidirectional yield stress. The D-P yield criterion is as follows:

. It is assumed here that viscoplastic flow only occurs when . Assume that the change of viscoplastic strain rate obeys the orthogonal flow rule, then the following equation can be obtained:where *F*_{0} is yield equation, ** Q** is the plastic potential, and

*γ*is the flow parameter. According to the plastic flow criteria,

*Q*≡

*F*. is a switching function, which can be defined as follows:

The general expression of Φ can be written in the form of (*F*–*F*_{0})/*F*_{0}, i.e.,andwhere **M** and **N** are constants. Equation (7) is used in this program, and ** N** = 1.

##### 2.2. Viscous Strain Increment

The viscous strains can be decomposed into viscoelastic strain and viscoplastic strain . The viscoelastic strain [35, 36] is expressed as follows:Then

The viscoplastic strain increment [37] is as follows: can be obtained by the following equation:where *γ* is the flow parameter. The vector {*a*}^{T} of the solid element is as follows:where . In the above equation, and letwhere is the partial stress component. Thus, equation (12) could be rewritten as follows:

The yield function of a plane element is as follows:

Then the vector {*a*}^{T} is as follows:where , and are the stress components in the local coordinate system of the plane element, respectively.

##### 2.3. Stress Increment

From the increments in equations (9) and (10), we could obtain the below equation:

By choosing the appropriate morphological function, the following equation can be obtained:where [** B**] is the geometric matrix and

**{∆**

*u*}_{t}is the displacement increments. Substitute equations (9), (10), and (16) to equation (15),where [

**] is the elasticity matrix. The viscous strain increment [38] is as follows:**

*D*Thus the total strain is as follows:

##### 2.4. Equilibrium Equation

At any time ** t**, the following equilibrium equation needs to be satisfied:where [∆

*F*]

_{t}= 0 indicates the change in load over the time interval

**∆**. Usually, [∆

*t**F*]

_{t}= 0 is satisfied for all time steps except for the first incremental change in the time interval. Substitute equation (17) into equation (20) to obtain the following equilibrium equation:where [

**] is elastic stiffness matrix [39], which can be expressed as follows:**

*K*Thus, the total stress and total displacement at **t** **+** **∆t** are as follows:

##### 2.5. Equilibrium Modification

Since the stress increments are obtained according to the linearized form of incremental equilibrium equation, the total stress **{ σ}**

_{t+∆t}obtained by accumulating all the above incremental stresses is inaccurate and does not accurately satisfy the general equilibrium equation. The following method is applicable to correct the equilibrium. First, the residual force at

**t**

**+**

**∆t**can be calculated.

The residual force is added to the external force increment in the next step. This method can not only eliminate the need for an iterative process but also reduce errors.

##### 2.6. Convergence Criterion and Time Step

At the end of **∆t**, if the viscoplastic strain rate is very small, it can generally be approximated that the displacement has been stabilized. In the program, the equivalent viscoplastic strain rate [40] is usually used to judge, which is defined as follows:

Therefore, the convergence criterion can be expressed aswhere **V**_{a} is the allowable value. The summation in the above equation is performed on all yield Gauss integration points, and **V**_{a} provided by Owen is 0.01. In order to find the steady-state displacement, the size of each time step **∆t** should be properly selected. There are two ways to limit **∆t**. The first empirical expression is as follows:where is the equivalent strain of total strain, and is the equivalent strain rate of viscoplastic strain rate. The minimum value of is calculated by summing the Gaussian integral of all yields. In the explicit method, *τ* can be set to any value within the range of 0.01–0.15 to obtain accurate results.

The second method is suitable for variable step size method, i.e.,

To limit the time step change between any two intervals, usually, ** k** = 1.5. In this program, the two methods, i.e., equations (31) and (32), were used to control the time step.

##### 2.7. Calculation Steps

To solve the viscoelastic-plastic problem, the elastic static definite problem should be calculated from the time point of **t** = 0. In this step, **{u}**_{0}, **{F}**_{0}, **{ ε}**

_{0}, and

**{**

*σ*}_{0}are obtained. Then, based on the derived equations above, the problem can be solved step by step. The calculation steps are as follows.(1)Assume that the equilibrium state has been reached at time

**, and the values of**

*t***{u}**

_{t},

**{**

*σ*}_{t},

**{**

*ε*}_{t},

**{**

*ε*_{ve}

**}**

_{t},

**{**

*ε*_{vp}

**}**

_{t}, and

**{F}**

_{t}and the stress state of element Gaussian point are known, the following equations can be calculated: where

*γ*is a flow parameter.(2)Calculation of displacement increment and stress increment(3)Calculation of total displacement and total stress(4)Calculation of viscoplastic strain rate(5)Equilibrium modification Substitute the total stress into the equilibrium equation to obtain the residual force. Then add the residual force to the next increment.(6)Selection of next time step size(7)Verification of calculation convergence

The convergence of structure is determined. If a satisfactory result is achieved, we can either terminate the solution-seeking process or apply the next load increment; otherwise, the program returns to the first step and repeats the entire process. The programming flowchart is shown in Figure 2.

##### 2.8. Numerical Example

The example provided by Owen and Hinton [31] is used for verification. This example is a thick-walled cylinder, with 100 mm as inner radius and 200 mm as outer radius. A quarter of the cylinder is taken for study (Figure 3). A total of 200 mm in the axial direction is selected with constraints at both ends. The study object can be considered as a planar state to compare with the calculation examples in the above literature. The internal pressure is = 1.4 MPa. The elastic modulus is 2.1 GPa, the Poisson’s ratio is 0.3, and the uniaxial yield stress is *σ*_{γ} = 2.4 MPa. The strain strengthening is not considered, i.e., *H′* = 0, and the flow parameter is *γ* = 0.001/d. The initial time step is 0.01 d, the step size is 0.01, and ** k** = 1.5.

###### 2.8.1. Instantaneous Loading ( = 1.4 MPa)

Figure 4 shows the variation of radial displacement of the inner surface when = 1.4 MPa.

###### 2.8.2. Time-Stepped Loading (Loading Twice)

In the first loading, = 1.2 MPa, then in the second loading, = 0.2 MPa. Figure 5 shows the variation of radial displacement of the inner surface over time with the stepped loading. The distribution of circumferential stress is shown in Figure 6.

The above calculation results were in good agreement with the results in citations.

##### 2.9. Definition of Point Safety Factor

In the FEM analysis, the M-C criterion is used to judge the destruction of the rock layer. The point safety factor of the complete rock layer is defined as follows:where *σ*_{1} and *σ*_{3} are the maximum and minimum principal stress, respectively. For the faults and weak structural surfaces, the point safety factor can be defined as follows:where is local coordinates of faults and weak structural surfaces.

At a greater safety factor, the point is further away from the destruction state. *F*_{p} = 1 indicates that the point is at the limit, and *F*_{p} < 1 indicates that the rock layer has been destroyed.

#### 3. Simulation on the Destruction Law of Overburden Rock at Lateral Boundary

##### 3.1. A Brief of Simulation Model

Based on the viscoelastic-plastic FEM, the destruction law of overburden rock is analyzed. The mining condition is assumed as follows: The coal seam is 6 m in mining height and buried 600 m underground.

According to fully mechanize caving practice, the height of the hydraulic fracture zone is approximately 10 times as large as the thickness of the coal seam. For a coal seam with a thickness of 6 m, the height of the hydraulic fracture range is 60–70 m. The step size of the first weighting is 30–40 m, and that of periodic pressure is 10–15 m. Thus as shown in Figure 7, a profile along the lateral direction of the stope is analyzed. The thickness of the floor in the finite element model is 20 m; the roof is 100 m in height, and 60 m on each side. The front-back distance is 50 m. The top boundary is loaded with gravity stress and other boundaries are constrained.

The overlying rock is composed of 8 layers, and the thickness of each layer is 4 m. In order to reflect the movement between different rock beams, the interface between layers is simulated with a joint of 0.1 m. After destruction occurred, the strength of the rock decreased accordingly. During the analysis, the fallout zones and hydraulic fracture zones in the disruption state can be achieved by lowering the mechanical parameters of rock formations.

In the analysis, the inclination angle is assumed to be zero. The plane strain model is adopted using the 8-node isoparametric element. As shown in Table 1, the rock parameters were obtained from mechanical and creep tests.

##### 3.2. Destruction Law of Overlying Rock Structure

###### 3.2.1. Four Kinds of Overburden Structures

According to the actual rock characteristics, the overburden structures at the lateral boundary of stope can be described by the following four lithological states:

*(1) Hard-Soft Structure*. A hard rock layer is located within 4 m above the coal seam, and then a relatively soft layer is located upward. The duration of deformation and destruction after extraction lasts for 40 days.

*(2) Soft-Hard Structure*. A soft rock layer is located within 4 m above the coal seam, and then a hard rock layer is located above the soft layer. The destruction of the rock mass can last for 45 days. Compared with the first type of structure, the duration of deformation is prolonged, which is related to the fact that most overburden strata are hard rock layers.

*(3) Soft-Soft Structure*. When the overburden strata were composed of soft rock layers, the destruction of rock mass lasts for 37 days. Compared with the first two states, its duration is shortened, which is associated with the overburden being a soft rock layer.

*(4) Hard-Hard Structure*. When overburden strata were composed of hard rock layers, the deformation lasts for 46 days. Compared with the first three states, the duration of failure is extended, which is related to the overburden being hard rock layers.

###### 3.2.2. Point Safety Factor

Figure 8 shows the contour of the point safety factor in the overburden strata.

**(a)**

**(b)**

**(c)**

**(d)**

*(1) Hard-Soft Structure*. The point safety factor of coal seam laterally near the stope is less than or equal to 1.0, indicating that the coal body is already in a yield state. The fracture of the hard rock layer in the upper part of the coal extracting layer is located about 15 m ahead of the lateral boundary. It can be found that the coal body within 25 m in front of the lateral boundary has been destroyed, which is due to that the coal body is pressed by an overlying rock beam. From Figure 8, we conclude that the coal body close to the floor layer has a larger damage range. In addition, as the position of fracture moves upwards away from where the coal seam locates, the fracture gradually slopes toward the goaf. Obviously, this change in the fracture position is related to the characteristics of the overburden structure.

*(2) Soft-Hard Structure*. There is a band of overburden with a point safety equal to and less than 1.0, which indicates that the rock mass in the area has been damaged. It is obvious that the fault zone is deflected upward and develops in front of the coal body. This phenomenon is exactly opposite to the law of fault development in State 1. In addition, the point safety factor of the coal body in front of the lateral boundary is less than 1.0, which indicates a damaged state of this area.

*(3) Soft-Soft Structure*. In a band area of overburden strata in front of stope’s lateral boundary, the point safety is 0.9. In addition, an isoline of 1.0 can be observed further nearby. This result indicates that the overburden rock layer is fractured along this area. Compared with State 1, in which the overburden rock layer is composed of H-S rock layers, the fracture position had an increased inclination angle toward the goaf. Obviously, this is related to the lithology of overburden rock. In addition, in front of the lateral boundary of the stope, the point safety factor is less than 1.0, which indicates that the coal body is in a destructed state.

*(4) Hard-Hard Structure*. The point safety factor of rock mass is less than 1.0 in front of the lateral boundary. Thus this area can be regarded as damaged. Compared with State 2, in which the overburden is composed of S-H rock layers, the fracture position has a larger inclination angle toward the front of the lateral boundary. In addition, the safety factor is observed to be less than 1.0 in some areas of the coal body, indicating that this part of the coal body is destroyed.

###### 3.2.3. Maximum Principle Stress

Figure 9 is the contour of maximum principal stress above the lateral boundary.

**(a)**

**(b)**

**(c)**

**(d)**

*(1) Hard-Soft Structure*. In the contour of maximum principal stress, the tensile stress is positive and the negative value represents compressive stress. Tensile stress occurs in the rock layer above goaf, and it extends throughout the upper rock mass. The maximum value of tensile stress is 32.7 kPa. So, the surrounding rock is easily destroyed in this case. Along with the fracture position of the rock beam, the compressive stress value is in the range of −42.0 ∼ −27.0 kPa.

*(2) Soft-Hard Structure*. Same with that of Case 1, the tensile stress appeared in the rock layer above the goaf, while the maximum value is 5.5 kPa. The maximum principal stress increased from goaf to the front of the lateral boundary. Meanwhile, the compressive stress reaches its maximum value, i.e., −40.2 kPa, near the fault zone of the overburden rock beam.

*(3) Soft-Soft Structure*. The maximum value is −4.1 kPa, which is close to tensile stress. The maximum principal stress at the fracture position of the overburden rock layer is −52.3 kPa. In addition, in the coal body in front of the lateral boundary, the maximum principal stress is −4.1 kPa, and the compressive stress is very small, which is close to tensile stress.

*(4) Hard-Hard Structure*. Tensile stress appears in the overburden strata at the lateral boundary, with a maximum value of 5.4 kPa. At the fracture position of the rock beam, the compressive stress value of maximum principal stress is −59.1 kPa. In addition, in front of the lateral boundary, tensile stress appears at a position close to the coal floor with a maximum value of 53.8 kPa. The appearance of this value may be related to boundary conditions; however, from the perspective of stress distribution, a tensile stress is observed in the coal body at this position.

###### 3.2.4. Minimum Principal Stress

Figure 10 shows the distribution of minimum principal stress.

**(a)**

**(b)**

**(c)**

**(d)**

*(1) Hard-Soft Structure*. Tensile stress appears at the left boundary of upper rock mass, which contributes to the damage of rock mass. In some areas above the goaf, the maximum and minimum stresses are both tensile stresses. In the fracture area of the rock beam, the minimum principal stress is in the range of −202.9 ∼ −167.7 kPa.

*(2) Soft-Hard Structure*. The minimum principal stress kept increasing from the lateral boundary to the fracture position of the rock beam, reaching the maximum value of −134.2 kPa^{.}

*(3) Soft-Soft Structure*. The maximum value of minimum principal stress, −230 kPa, is observed at the fracture position.

*(4) Hard-Hard Structure*. The maximum value of minimum principal stress is observed at the fracture position of the overburden rock layer, which is −228.2 kPa. The minimum principal stress gradually decreases from the fracture position of overlying rock strata to both sides. The minimum principal stress decreases faster towards the lateral boundary of the stope.

The analyzing results are essentially consistent with general cases [41–45], which means that the V-P FEM is useful to detect the failure mode. This study provides a convenient and economical method for determining the safety of overburden strata on the lateral boundary.

#### 4. Conclusions

In this research, a visco-elastoplastic model for analyzing rock strata movement is established on the basis of previous studies. Considering four kinds of overburden rock structure forms, we adopt the developed method to study the lateral boundary rock failure mode of the stope. In addition, we defined the point of safety to judge whether the rock layer is damaged.

According to the characteristics of stope overburden, the mentioned structural models represent different combinations of overburden rock layers. The numerical computation indicates that the fracture modes of rock beam in the coal wall are various at different combinations.(1)For the overburden rock structure composed of Hard-Soft layers, the hard rock beams first break in front of the lateral boundary of the stope, while the fracture position in the upper soft rock beams gradually moves towards the goaf. The fault zone of overburden rock layers inclines toward the goaf. The maximum principal stress above the goaf is tensile stress.(2)For the overburden rock structure composed of Soft-Hard layers, the soft rock layer first breaks in front of the lateral boundary, and then the fracture position in the lateral upper hard rock gradually develops in the rock body in front of lateral boundary, forming a fracture zone.(3)For Hard-Hard rock layers, the location of the fracture zone develops toward the front of the stope’s lateral boundary, but its inclination is greater than that of the second state.(4)For Soft-Soft rock layers, the fracture zone of the overburden layer develops towards the goaf, but the inclination angle is greater than that in the first state.(5)The time from the occurrence of fracture to stabilization varies for different overburden structures, which can be ranked according to the length of time as follows: S-S < H-S < S-H < H-H.

In practice, the location of fractures in overburden rock can vary widely due to the complexity of rock mass. The above results are qualitative conclusions for the specific overburden structures, which are also the four most typical scenarios. Therefore, the calculation conclusions are representative.

#### Data Availability

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

#### Conflicts of Interest

The authors declare that they have no known conflicts of interest.

#### Acknowledgments

This work was supported by the Key Laboratory of Deep Coal Mine Excavation Response and Disaster Prevention and Control, Anhui Province (Anhui University of Science and Technology), Huainan, China, 232001, (KLDCMERDPC16101 and 01CK05002), Jinxiao Liu ([email protected]).