Analysis for Overburden Destruction on Lateral Boundary of Stope Based on Viscoelastic-Plastic Finite Element Method

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 classiﬁed into four categories, i.e., Hard-Soft, Soft-Hard, Hard-Hard, and Soft-Soft. In order to analyze the eﬀect of the above four structures, we adopt viscoelastic theory to the ﬁnite element method (FEM) and deﬁne the point safety factor to evaluate the rock damage. The accuracy of programming is veriﬁed through example veriﬁcation. A modiﬁed 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.


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 mineinduced 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. e lateral boundary is on either side of the goaf. For overlying rock layers with different lithological structures, the destruction patterns are different. e 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. ough automatically formed entry retaining technology, Manchao He [14][15][16][17][18] developed a nopillar 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.
e overlying rock will fall and fill the goaf after coal extraction. en 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. e 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. e 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.

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. en, it is subsequently verified by the classical solution of a viscoelastic-plastic (V-P) mechanical problem. e iterative method of viscoelasticity and plasticity is based on the approach proposed by Owen and Hinton [31].

Constitutive Equation.
e 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 ε v . erefore, 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. e yield function is introduced to determine whether there is viscoplastic strain: where F 0 is the unidirectional yield stress. e D-P yield criterion is as follows: F 0 � k. It is assumed here that viscoplastic flow only occurs when F � F 0 . 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 c is the flow parameter. According to the plastic flow criteria, Q ≡ F. 〈Φ(x)〉 is a switching function, which can be defined as follows: e general expression of Φ can be written in the form of (F-F 0 )/F 0 , i.e., 2 Advances in Civil Engineering where M and N are constants. Equation (7) is used in this program, and N � 1.

Viscous Strain Increment.
e viscous strains ε v can be decomposed into viscoelastic strain ε ve and viscoplastic strain ε vp . e viscoelastic strain [35,36] is expressed as follows: en e viscoplastic strain increment [37] is as follows: ε vp t can be obtained by the following equation: where c is the flow parameter. e vector {a} T of the solid element is as follows: e yield function of a plane element is as follows: en the vector {a} T is as follows: where σ z′ , and σ z′y′ are the stress components in the local coordinate system of the plane element, respectively. (9) and (10), we could obtain the below equation:

Stress Increment. From the increments in equations
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 [D] is the elasticity matrix. e viscous strain increment [38] is as follows: us the total strain is as follows:

Equilibrium Equation.
At any time t, the following equilibrium equation needs to be satisfied: where

Advances in Civil Engineering
Substitute equation (17) into equation (20) to obtain the following equilibrium equation: where [K] is elastic stiffness matrix [39], which can be expressed as follows: us, the total stress and total displacement at t + ∆t are as follows:

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. e following method is applicable to correct the equilibrium. First, the residual force at t + ∆t can be calculated.
e residual force is added to the external force increment in the next step. is method can not only eliminate the need for an iterative process but also reduce errors.

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 _ ε vp [40] is usually used to judge, which is defined as follows: erefore, the convergence criterion can be expressed as where V a is the allowable value. e 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 steadystate displacement, the size of each time step ∆t should be properly selected. ere are two ways to limit ∆t. e first empirical expression is as follows: where ε is the equivalent strain of total strain, and _ ε vp is the equivalent strain rate of viscoplastic strain rate. e minimum value of _ ε vp 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. e 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. en, based on the derived equations above, the problem can be solved step by step. e calculation steps are as follows.
(1) Assume that the equilibrium state has been reached at time t, and the values of and {F} t and the stress state of element Gaussian point are known, the following equations can be calculated: where c is a flow parameter.

(2) Calculation of displacement increment and stress increment
Δu (3) Calculation of total displacement and total stress (4) Calculation of viscoplastic strain rate Advances in Civil Engineering Substitute the total stress into the equilibrium equation to obtain the residual force.
en add the residual force to the next increment.
(6) Selection of next time step size (7) Verification of calculation convergence e 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. e programming flowchart is shown in Figure 2.

Numerical Example.
e example provided by Owen and Hinton [31] is used for verification. is 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. e study object can be considered as a planar state to compare with the calculation examples in the above literature. e internal pressure is p � 1.4 MPa. e elastic modulus is 2.1 GPa, the Poisson's ratio is 0.3, and the uniaxial yield stress is σ c � 2.4 MPa. e strain strengthening is not considered, i.e., H′ � 0, and the flow parameter is c � 0.001/d. e initial time step is 0.01 d, the step size is 0.01, and k � 1.5. Figure 4 shows the variation of radial displacement of the inner surface when p � 1.4 MPa.

Time-Stepped Loading (Loading Twice).
In the first loading, p � 1.2 MPa, then in the second loading, p � 0.2 MPa. Figure 5 shows the variation of radial displacement of the inner surface over time with the stepped loading. e distribution of circumferential stress is shown in Figure 6. e above calculation results were in good agreement with the results in citations.

Definition of Point Safety Factor.
In the FEM analysis, the M-C criterion is used to judge the destruction of the rock layer. e 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 (x ′ , y ′ , z ′ ) 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. e 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.

A Brief of Simulation
In the analysis, the inclination angle is assumed to be zero. e 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.

Four Kinds of Overburden Structures.
According to the actual rock characteristics, the overburden structures at the Advances in Civil Engineering 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. e 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. e 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. Obtain the residual force by the total stress and equilibrium equation.
And iterate the residual force in the next incremental calculation.
Calculate the time step size of the next iteration step on account of the viscoplastic strain rate.
Judge whether Δt meets with the model setting At t, converge to equilibrium state and obtain all the physical quantities.

End and export
Model setup: geometry, property, initial and boundary conditions.

Yes
Update the strain, stress, variables, and time.
Calculate the stiff matraix.
Save the strain, stress, variables, and time.  Advances in Civil Engineering Figure 8 shows the contour of the point safety factor in the overburden strata.

Point Safety Factor.
(1) Hard-Soft Structure. e 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. e 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. ere 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. is 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. is 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.
e point safety factor of rock mass is less than 1.0 in front of the lateral boundary. us this area can be regarded as damaged. Compared with State    Advances in Civil Engineering 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. Figure 9 is the contour of maximum principal stress above the lateral boundary.

Maximum Principle Stress.
(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. e 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.

Advances in Civil Engineering
(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. e 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. e maximum value is −4.1 kPa, which is close to tensile stress. e 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. e 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. Figure 10 shows the distribution of minimum principal stress.

Minimum Principal Stress.
(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.        Advances in Civil Engineering e analyzing results are essentially consistent with general cases [41][42][43][44][45], which means that the V-P FEM is useful to detect the failure mode.
is study provides a convenient and economical method for determining the safety of overburden strata on the lateral boundary.

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. e 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. e fault zone of overburden rock layers inclines toward the goaf. e 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) e 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. e above results are qualitative conclusions for the specific overburden structures, which are also the four most typical scenarios. erefore, 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
e authors declare that they have no known conflicts of interest.