Abstract
Universal distinct element code (UDEC) is a simulation software based on the discrete element method, widely used in geotechnical mining. However, in the UDEC, when simulating large-scale excavation, the subsidence of the fractured zone is almost equal to the mining height, which makes the deformation value calculated in the study of gob-side entry retention too large. To solve this problem, in this paper, the double-yield constitutive model is applied to the whole caving zone to analyze the deformation and failure characteristics of surrounding rock along gob-side entry retaining by roof cutting. The results of the simulation are in good agreement with the result of drilling peeking (drilling observation by borehole televiewer) and field condition (observation and measurement in the field). Finally, by using this numerical method, the effects of roadway width, temporary support, and coal side support on the failure of the roof and the arc coal side are studied.
1. Introduction
Longwall mining is the main method of coal mining in Russia, Poland, and China. Generally, two roadways are needed to meet the basic requirements of coal mining and ventilation. There are two sources of roadway. The first is to dig new roadway, and the second is to retain the old roadway. The new and original roadways need to leave a certain width of coal pillar, also known as coal pillar mining, which will lead to a waste of resources. On the contrary, retaining the old roadway (gob-side entry retention belongs to this category) does not need to leave coal pillars, which is called no pillar mining. Gob-side entry is to set roadside support (supporting materials or equipment) along the edge of gob in the original roadway and retain the roadway for the use of the next working face. The roadside support has developed from timber props, pack walls, blocks, and newly developed filling materials, such as the material with high water content [1] and the flexible formwork and pump-concreted support [2], to name a few. This method has gained wide applicability due to safety and economic benefits. However, when the height of the roadside backfill exceeds 3.5 m, the roadside backfill stability decreases [3], which limits the application of this method in thick coal seams. To address this limitation, a new method of gob-side entry retention by roof cutting (gob-side entry retention hereafter) was proposed by He et al. [4]. In this method, the roof of the roadway along the gob was cut through the bilateral cumulative tensile explosion, and the roof along the kerf side was supported by the roof cutting support. The overlying strata in the caving zone automatically form the roadside with periodic weighting to retain the roadway. To date, this method has been successfully implemented in more than 100 coal mines with varied geological conditions in China.
Since the roof cutting and the pressure release technologies advanced, associated theoretical frameworks have been developed rapidly as well. Wang et al. [5] establishes a mechanical model of cutting cantilever beam through material mechanics. He et al. [6] studied the influence of key parameters (height and angle of kerf) of roof cutting on the effect of retaining roadway, the filling characteristics of gob. Yang et al. [7] studied the effects of different cutting positions on the periodic weighting. Little research has been carried out to study the deformation and failure of the surrounding rock of the gob-side entry retention by roof cutting. Due to the complex field conditions under which a lot of assumptions have to be made to simplify the conditions, the conclusions drawn from structural or elastic mechanics perspective have neglected the discontinuity of rock mass, thus unable to effectively guide the field applications. UDEC though considers the discontinuity of the rock mass, when simulating large-scale excavation, the subsidence of fractured rock strata is almost equal to the mining height, which makes the deformation value calculated in the study of surrounding rock deformation and failure along gob-side entry retention too large. However, in the process of coal mining, the rock strata in the caving zone fall and fill the gob. When the rock strata in the fractured zone sink and compact the gangues, the pressure of the gauges increases rapidly and the subsidence velocity of the rock strata in the fractured zone moves from fast to slow until an equilibrium is reached. In the end, the subsidence of fractured rock strata is much less than the mining height.
In order to eliminate the defects of UDEC, the double-yield model (a constitutive model built in the software), which is similar to the broken gangue compression process, is proposed for simulating the compaction process of the whole caving zone. By comparing the values (vertical pressure in the gob and the roof subsidence) of the proposed simulation model, the traditional simulation model, and the field, the rationality of the proposed simulation model is verified. In addition, combining the numerical simulation and field observations, the deformation and failure characteristics of the roof and the curved coal wall of the roof cutting and pressure release technology in a thick coal seam are explored. Results suggest potential applications to guiding the site operation and optimizing the support scheme.
2. Profile of the Coal Mine Studied
Ningtiaota Coal Mine is located in northern Shaanxi Province, China. The main coal seam currently being developed is the 2-2 #coal seam of the Mid-Jurassic Yan’an Formation. The thickness of the coal seam in the working face S1201-II is about 4 m, the depth of burial is about 150 m, and the coal seam is nearly horizontal with an inclination angle of <2°. The thick seam is mined using a fully mechanized coal mining technology with a large mining height. The plane view of the S1201-II working face is shown in Figure 1.

The Ningtiaota Coal Mine has an annual production of 18 million tons of coal and is a super-large mine. In order to meet the requirements of high production and high efficiency, the auxiliary transportation system suggests use of rubber-tired trucks, which requires a larger cross section size of the auxiliary transportation roadway. The height of the auxiliary transportation roadway in S1201-II working face is 3.75 m, and the designed section area is up to 24.3 m2; the roadway is an extralarge section with no other successful precedence to date to retain the roadway.
The gob-side entry retention by roof cutting was tested in the S1201-II auxiliary transportation roadway. The difference of gob-side entry retention by roof cutting and by packing is shown in Table 1 and Figure 2, and the technological process of this new method with roadside cutting is shown in Figure 3.

(a)

(b)

(a)

(b)

(c)

(d)

(e)

(f)
2.1. Initial Expansion Coefficient of the Caving Zone
Figure 4 shows the diagram of roof falling. After coal seam mining, the immediate roof starts to fall gradually, and after a certain distance behind the working face, the height of gangue increases up and touches the roof plate. The rock in the gob is then in a loose state with small vertical stress. With the working face advanced, the roof falls down to compress the gob, and the vertical stress of gangue increases until the roof reaches a new balance.

Huang et al. [8] studied the working face with a large mining height in the Shenfu mining area (the Ningtiaota Coal Mine is in the area) in China. The relationship between the height of the caving zone and the mining height was examined by physical models using a large number of physical models. When the mining height is 4 m, the height of the caving zone is about 2.5 times that of the mining height. Bai et al. [9] statistically analyzed a large number of mines in China and United States and concluded that the height of the caving zone (H) satisfies empirical equation (1), which is widely accepted and applied to practical use:where M indicates the height of the coal seam and c1 and c2 are the correlation coefficients. In this case, the compression strength of the roof rock layer is between 20–40 MPa; therefore, c1 = 4.7 and c2 = 19.
The heights of the caving zone obtained using the two methods above are 10 m and 10.6 m, respectively, and the average value is 10.3 m. The expansion coefficient (the ratio of the height of rock in loose state after caving to the height of rock in whole state before caving) when the gob contacts the roof is considered as the average initial expansion coefficient K. Equation (2) can be obtained by the above definition and be substituted with data to get K = 1.39:
2.2. Stress-Strain Relationship in the Compression Process of the Caving Zone
Researchers studied the stress-strain relationship in the compression process of the various fractured rock mass using special compression devices in the laboratory [10–13], where equation (3) is widely accepted and applied. The type of curve represented by equation (3) has the following characteristics; the vertical stress of the fractured rock mass in the initial state is small. As the strain increases, the vertical stress of the fractured rock mass increases rapidly.where E0 represents the initial tangent modulus of the caving zone, ε stands for strain, and εm (εm = (K − 1)/K) is the limit strain when the expansion coefficient is 1 in the compression.
The initial elastic modulus E0 is mainly affected by the expansion coefficient of rock mass and rock uniaxial compressive strength. Based on a large number of uniaxial compression data obtained from test sands through regression analysis, Yavuz [14] proposed equation (4). The limit strain εm indicates that the broken and expansive gangue recovers to the original state under the roof pressure. It is no more than a value under the ideal state. As per K = 1.39, it leads to εm = 0.28.where σc represents the rock uniaxial compressive strength and K represents the average initial expansion coefficient.
This section may be divided by subheadings. It should provide a concise and precise description of the experimental results, their interpretation as well as the experimental conclusions that can be drawn.
3. Numerical Simulation Method for Load-Bearing Characteristics of Caving Zone
3.1. Denumerical Simulation of the Caving Zone
UDEC is a two-dimensional numerical simulation software based on discrete elements. It is used to stimulate the fracture, caving, and compaction of overlying strata. However, when UDEC simulates large-scale excavation, the subsidence of fractured zone rock strata is almost equal to the mining height, which means the UDEC cannot reflect the expansion characteristics of the caving zone and makes the calculated deformation value of surrounding rock and failure along gob-side entry retention too large. The viewpoint is supported by other researchers [15]. In a study by Lu et al. [16], a customized constitutive model with the consideration of expansion was applied to UDEC, but this way is pretty complicated. When the object to be studied is other than the caving zone, but such as in this paper, the analysis of the movement of surrounding rock of gob-side entry retaining roadway and the caving zone load-bearing characteristic can be stimulated through the double-yield constitutive model.
The double-yield constitutive model was used in a numerical simulation within the compaction of the caving zone to find the mechanical parameters that satisfy equation (5). The parameters obtained for the gob are shown in Table 2, and the results are shown in Figure 5.

3.2. Comparison of Numerical Simulation Methods
The numeric computation model is shown in Figure 6. The model only consists of one part of S1201-II working face at a length of 110 m. In this figure, Trigon-BMP is used to divide random triangular blocks, and ordinary rectangular blocks were used elsewhere.

The traditional UDEC simulation method and the new method considering the double-yield constitutive model proposed in this paper are mentioned above. The numerical experimental results obtained by simulating the coal seam excavation of these two methods are compared below. There is one difference between the two numerical simulation methods. In the first method, the traditional method, after the initial stress is balanced, the coal seam in the working face and roadway is excavated and a weaker joint parameter is applied to the kerf. In the second method, after the initial stress reached a balance, the constitutive model for the falling zone within the red line is set as double-yield and the parameters are defined according to Table 2.
Figure 7 shows the monitored subsidence of the roof near the kerf at 0.5 m and in the middle, in which Figure 7(a) is related to the first method and Figure 7(b) is related to the second method.

(a)

(b)

(c)
When simulating using the first method, the subsidence of point 1 and point 2 shown in Figure 7(c) is 251 mm and 195 mm, respectively, and when the time step of computing is 150000 steps, the roof subsidence still tends to increase rapidly; whereas when simulating using the second method, the subsidence of point 3 and point 4 shown in Figure 7(c) is 116 mm and 107 mm, respectively, and when the time step of computing is 60000 steps, the roof subsidence tends to become a constant. As shown in Figure 7(c), at the initial stage, no significant difference of the numerical simulation presents between the first and the second method and the results of both methods showed a rapid increase. However, in method 2, the supporting pressure of gob rock mass to the roof increases gradually as the time step increased and the roadway surrounding rock reaches the stable stage faster compared with method one.
The field monitoring results shows that, in the normal gob-side entry retention, when the roadway roof reaches a stable stage, the average amount of roof subsidence near the kerf is 102 mm, indicating a better fitting of stimulation results obtained using the second method to the actual data obtained in the field.
3.3. Calibration of Mechanical Parameters Used in Modeling with UDEC
According to the SB01 drilling in the S1201-II working face, the intact rock mechanics parameters and rock quality designation (RQD) values of the coal measure strata were obtained. The elastic modulus, uniaxial compressive strength, and tensile strength of rock mass were obtained using equations reported in paper [17].
In Figure 6, the rectangular region outlined by the dotted line is the key region in this study. Otherwise, traditional rectangular block was used and the intact rock parameters were used as the block parameters. In the key region, the block was generated through Trigon-BMP, the Mohr–Coulomb model was used, and the Coulomb slip model was applied to joints. The blocks and joints parameter were calibrated as shown in Table 3. The final parameters of blocks and joints are shown in Table 4.
4. Numerical Analysis of Deformation and Failure Characteristics for the Surrounding Rock of the Retaining Roadway
4.1. Basic Characteristics of Deformation and Failure for the Surrounding Rock under Normal and Abnormal Conditions: Normal and Abnormal Retaining Roadways
When the working face was advanced to a distance of about 860 m from the open-off cut, lagging working face of about 100 m, the top of the curved coal wall started to display rib spalling, and a crack is presented in the middle of the roof. This might be due to the following reasons. (1) Adjusting working face support offset results to enlarging the width of retained roadway even up to 7 m. (2) The coal side was supported lagged. Because of the interference of the coal mining and roadway retaining and not enough constructing labors, coal side support was much lagged from the working face. (3) The temporary support scheme was changed. Temporary support refers to the support members that can be withdrawn after the gob is stabilized. The original temporary support scheme includes one roof cutting support, two π-shaped steel beams, and four single props. About 760 m off from the open-off cut, the temporary support scheme was changed to one roof cutting support, one π-shaped steel beams, and two single props, as shown in Figure 8. Among the temporary support, the π-shaped steel beams were located at the top of the roof cutting supports and single props. The effects of these factors interfere with each other impacting the retained roadway, and this is referred to as abnormal roadway retaining.

(a)

(b)
After determining the reasonable support parameters, the coal mining process and the roadway work were optimized. The roadway deformed within the reasonable range and was retained well. This is called the normal roadway retaining.
4.2. Comparison of Numerical Simulation Schemes
Considering the three problems mentioned above, we designed different schemes for numerical simulation of normal and abnormal gob-side entry retention. For the normal roadway retention, the width of the roadway was set to 6 m and the coal wall was supported with bolts. For the abnormal roadway retention, the width of the roadway was set to 7 m and the coal wall was not supported. So, the differences in the first two problems are easy to achieve in numerical simulations. For the third problem, approximations were made for the numerical simulation of the temporary support combination of the roof cutting support, π-shaped steel beam, and single prop.
Since the trackless rubber tire vehicle travels in the roadway, and no support should be put in the middle of the roadway. The π-shaped steel beam is to transfer partial supporting force from the roof cutting support and the single props to the roof in the middle of the roadway. Therefore, the combined function of roof cutting support, π-shaped steel beam, and single prop can be considered a common support by 3 support units, which were labeled i, j, and k from left to right. Because the supporting force of the roof cutting support is much greater than that of the single prop, it can be approximated that the supporting force in the middle of the roadway roof comes from the roof cutting support, that is, the sum of the supporting force of support units i and j, which can be determined according to the working resistance of the roof cutting support, while the supporting force of support unit k can be determined according to the single prop. In addition, the decrease in the number of π-shaped steel beam can be approximately considered to be the reduced supporting force to the middle roof.
In summary, in the normal retaining roadway, the maximum supporting force of units i, j, and k was set to 1641 kN, 25 kN, and 250 kN, respectively. In the abnormal retaining roadway, the maximum supporting force of units i, j, and k was set to 1653.5 kN, 12.5 kN, and 125 kN, respectively. In fact, the above values are only qualitative estimates because it is difficult to evaluate the accurate contribution of roof cutting support to the supporting force in the middle of the roadway.
4.3. Analysis of Numerical Simulation Results
Figures 9(a) and 9(b) show the crack growth of the surrounding rock in normal retaining roadway and abnormal retaining roadway, respectively. In the figures, the pink lines represent the open cracks and the blue lines are for the slip cracks. To explore the distribution of cracks, we defined the following zones: zone I is the upper left of the kerf, zone II is the high rock layer on the top of the roadway, zone III is the upper left of the curved coal wall, and zone IV is the curved coal wall, all of which are outlined via a red dotted line in Figure 9.

(a)

(b)
As shown in Figure 9(a), for the normal roadway retaining, in zone I, the open cracks grow in a complex way and the tensile and slip cracks grow in stagger. In zone II, the open cracks grow in a simple pattern and the cracks are long. The open cracks show the fracture location of the high-position strata. In zone III, the slip cracks are mainly distributed near the anchorage area, and the lower-position strata are almost intact. In zone IV, the main cracks are the slip cracks.
Figure 9(b) shows the surrounding rock fissure development of the abnormal roadway retaining. In zone II, the number of open cracks is less than that of the normal roadway retaining. In zone III, the cracks fully grow and the slip cracks near the anchorage area grow longer compared with that of the normal roadway retaining and are at a large inclination. The open cracks appear clearly near the middle of the roof. In general, the cracks in the zone III tend to grow from the lower right to the upper left. The crack growth directly impacts the deformation and failure of the roof in zone III. In zone IV, the density of slip cracks is larger than that of the normal roadway retaining.
To investigate the deformation characteristics of the surrounding rock under roof cutting and roadway retaining, the normal retaining roadway is taken as an example and seven plines in total were set up, three of which are above the roadway. Distances measured from the top of coal seam are 0 m, 5 m, and 10 m for pline 1, pline 2, and pline 3. There were four plines set on the coal wall. Distances measured from the curved coal wall are 0 m, 1 m, 2 m, and 3 m for pline 4, pline 5, pline 6, and pline 7. Among them, pline 1, 2, and 3 are used for monitoring the roof’s vertical displacement. The results are shown in Figure 10. Pline 4, 5, 6, and 7 were used to monitor horizontal displacement of the curved coal wall. The results are shown in Figure 11.


The results obtained for pline 3 show that after the working face and roadway space were excavated, the overlying strata bended. The closer to the gob, the larger the falling slope of the roof is. The results from pline 1 and 2 show a sudden increase in the falling slope at the coal wall. This could be attributed to the fact that the distance near the roadway decreases and the influence of the excavation on the roof subsidence increases. From pline 3 to 1, the roof falling curves develop into a lower convex pattern from a straight line, which shows that the lower the height is, the greater the roof displacement is influenced by the supporting.
Figure 11 shows the horizontal deformation characteristic of the curved coal wall observed from pline 4, 5, 6, and 7 as a function of height from the floor. The horizontal displacements of nearly the entire curved coal wall are to the right, which is quite different from that of the traditional vertical coal wall [18]. The curve of pline 4 showed that the maximum horizontal deformation appeared at the highest point of the curved coal wall. When the height from the floor decreased from 4.0 m to 2.5 m, the horizontal displacement decreased rapidly; as the height decreased from 2.5 m to 1.0 m, the horizontal displacement is almost stable; when the height decreased from 1.0 m to 0 m, the horizontal displacement increased slowly. The results from the pline 5, 6 and 7 showed that as the distance into the coal wall increased, the coal seam height at which the horizontal deformation is minimum decreased. The horizontal deformation of line 7 is approximately linear to its distance to the floor, and the height of coal seam at which the horizontal deformation is minimum is close to 0 m. In summary, the top corner of the curved coal wall is most affected by compression and shear, indicating that this location should be the key area to be reinforced.
4.4. Field Observation
The S1202-II working face has been mined. The normal retaining roadway has reported previously [19, 20]. However, no study addressing the abnormal retaining roadway has been reported in the literature. This section presents the borehole televiewer images taken in the roof and the surrounding rock failure of the abnormal retaining roadway, as shown in Figures 12 and 13.


(a)

(b)

(c)

(d)
At a position 760 m away from the open-off cut, the cracks started to appear in the middle of the roadway roof and the rib spalling began to appear in the curved coal wall as shown in Figures 13(c) and 13(d). At Figure 13(c), the white paint in the middle roof indicates the direction of the development of the cracks above, and Figure 13(d) indicates the rib spalling of the curved coal wall. All pictures were taken after the reinforce support, and after the above phenomenon occurs, the density of the π-shaped steel beams and the single props on the coal side was increased.
To explore the bed separation above the roof, three sections with 790 m, 810 m, and 820 m offset from the open-off cut were selected for drilling peeping. Each section had three drillings holes, and the middle one was located at the center of the roadway. The drilling height is 15 m, the row spacing is 1000 mm, and the drilling diameter is 50 mm. Figures 12(a)–12(c) show the damages of the roof. The purple points represent that cracks are open and the bed separation is greater than 50 mm; the orange points represent that the cracks are open and the bed separation is less than 50 mm, and the green points stand for the closed cracks.
The borehole televiewer images reveal that (1) the distribution of the cracks in the roof of the gob-side entry retaining roadway is mostly inclined, the cracks observed by every hole can be connected as lines, and lines extend from the top left to the bottom right of the roof; (2) there are some differences in the damage characteristics of the roof of the three sections, which might be attributed to a large influence of the periodic weighting on the roof of the roadway; and (3) the cracks in the lower rock of the roof are mostly closed, but the closer to the anchorage area, the greater the degree of crack opening.
The borehole spreading pictures along the same lines above for the drillings about 790 m away from the open-off cut are shown in Figure 13. The principle of the borehole spreading pictures is shown in Figure 12(d), in which h represents the height of the bed separation and ∆h represents the height difference of the cracks, which can be used to calculate the dip corner of the cracks. The final result is shown in Figures 13(a) and 13(b).
In Figures 13(a) and 13(b), the three boreholes are numbered I, II, and III. The borehole televiewer images of each hole were lined up from high to low according to the position. The inclined distribution of the cracks at the upper roadway is obvious, and the crack characteristic varies at different heights. Based on the analysis of the bed separation and the developing direction of the crack, four cracks were selected for further explanation.
The cracking lines 1 and 2 are located near the anchorage area. The borehole spreading pictures show that the cracks in this area had the characteristics of large bed separation and large dip angle. The cracking lines 3 and 4 are located below the anchorage area, and cracks in this area had a small or even zero separation and a relatively small dip angle.
Compared Figure 13 with Figure 9(b), the observation on the abnormal retaining roadway revealed that the results of the open cracks in the center of the roadway, the inclined cracks in the roof, and the compressed-shear deformation at the tip areas of the curved coal sides are in good agreement with the numerical simulation results.
5. Effect of the Roadway Width, the Temporary Support, and the Curved Side Support on Deformation and Failure of Surrounding Rock of Retained Roadway
In the abnormal retaining roadway, the increased width of retaining roadway, the decrease in temporary support strength, and the lag of curved coal wall support contributed to a less stable retaining roadway. However, the effects of these factors on deformation and failure of surrounding rock remain unclear. Further discussion is made combining the field observations and the numerical simulation.
5.1. The Influence of the Roadway Width
When examining the influence of the roadway width on the deformation and failure of the surrounding rock of the roadway retaining, we keep the temporary support, the cable support, and the bolt support the same and the width of the roadway is the only variable. The results of the crack development of the surrounding rock are shown in Figure 14; among them, cables and bolts are not shown. The roadway widths in Figures 14(a)–14(d) are 4 m, 5 m, 6 m, and 7 m, respectively. Here, the pink line represents the open cracks and the blue line stands for the slip cracks.

(a)

(b)

(c)

(d)
The roadway width dependent simulation results showed the following. (1) When the width of the roadway retaining changed from 4 m to 7 m, the number of slip and open cracks in zone I decreased gradually while the open cracks in zone II gradually increased. Similarly, the slip and open cracks increased. The results showed that when the roadway width was less than 6 m, the roof above the roadway had good integrity, and the place where the fracture highly developed was located at the upper left of the kerf. As the roadway width was more than 6 m, the place where the fracture highly developed shifted above the roadway. (2) The cracks above the roadway were initially parallel to the direction of the kerf, starting from the top corner of the curved coal wall having an inclined distribution. Then, the cracks expanded on both sides. (3) When the width of the retained roadway changed from 6 m to 7 m, the growth of height of the open crack in the central roof of the roadway enhanced abruptly. The cracks became denser at the top of the curved coal wall.
The subsidence of the roadway roof surface was monitored, and the results are shown in Figure 15. Curves 1 through 4 are the results for the width of the roadway being 4 m, 5 m, 6 m, and 7 m, respectively. Since the widths are different, the x-axis was set to be the ratio of the distance to the kerf to the width of the roadway. When the width of the retained roadway changed from 6 m to 7 m, the roof subsidence on the kerf side increased by 47.4%. Furthermore, when the width of the retained roadway is 7 m, the roof subsidence plot showed a prominent “convex,” and a sharp corner appeared in the roadway center, indicating that the roadway roof generated an open crack.

In general, when the width of retained roadway changes from 6 m to 7 m, the developed height of the open cracks and slip cracks increases in the center of the roadway. The roof bends heavier and subsides larger, and the pressure on the top of the curved coal wall increases, which aggravates the compression-shear damage.
5.2. The Influence of Temporary Supports
To examine the influence of the temporary support on the deformation and failure of the retained roadway, the width of the roadway was maintained at 7 m with the cable and bolt support present. Only the temporary support varied in the form of the “roof cutting support, π-shaped steel beam, and single prop” support arranged. This can be viewed as the change of the central support strength. The simulated crack development of the surrounding rock is shown in Figure 16, in which the cables and bolts are not shown.

(a)

(b)

(c)

(d)
The varied temporary support can be considered as varied combinations of these three support units in the numerical simulation, and the support units are numbered as i, j, and k from the left to right in the roadway. Table 5 presents the maximum supporting force combinations for four groups of tests. Since the supporting force of the roof cutting support is much greater than that of the single prop, it is assumed that the supporting force of the support unit j was mainly provided by the roof cutting support, while the supporting force of the support unit k was kept the same.
The simulation results show the following. (1) Without the support force considered at the center of the roadway, the slip cracks and open cracks above the roadway were fully developed, and the height of the cracks in the middle of the roof was large; after a support unit with small support force added to the center of the roadway, the growth of the cracks was significantly reduced in the roof. Meanwhile, the highly developed area of the cracks transferred to the upper left of the kerf from right above the roadway. (2) As the central supporting force increased, the development of the cracks in the roof decreased and then stayed at a fixed level.
The subsidence of the roadway roof was monitored, and the results are shown in Figure 17. With the increase in the supporting force in the middle of the roadway, the subsidence of the roadway roof is obviously reduced, but the decreasing range is gradually reduced. Meanwhile, the falling of coal sides became light and weakened its compression and shear to the coal side top corner. Additionally, the roof crack opened less in the roadway center.

To conclude, increasing the supporting force in the middle of the roof helps to avoid the generation of open cracks in the middle of the roadway. Meanwhile, it reduced the amount of roof subsidence on the side of the curved coal wall and reduced the compression-shear damage on the curved coal wall.
5.3. The Influence of Coal Side Supports
To study the influence of coal supports on the deformation and failure of the surrounding rock in the roadway, the width of the roadway was kept 6 m, and the cable support was kept exactly the same as the previous setting. Only the timing of the coal side support changed. The results of the simulations of the crack development of the surrounding rock are shown in Figure 18, in which the cables and bolts are not shown. In the figure, from (a) to (d) is the coal side support without lag, with a 10000 steps lag, 20000 steps lag, and a 30000 steps lag, respectively.

(a)

(b)

(c)

(d)
The simulation results show that (1) with the lag of coal side support, the supporting capacity of coal wall decreased, and the inclined open cracks started to appear in zone III, which extended continuously from the top of the coal wall to the upper left; (2) the integrity of the rock mass in zone III damaged, which induced the rock strata in zone II to rotate further to the gob, and its open cracks moved to the right; (3) compared with the prompt supporting, the lagged supporting made the open crack in the roof center initiate at a slightly higher place.
The subsidence of the roadway roof was monitored, and the results are shown in Figure 19. Curves 1, 2, 3, and 4 are similar without obvious “lower convex,” and the change in roof subsidence on the kerf side is less than the previous two factors.

The lagging support has insignificant effect on the occurrence of the open cracks at the center of the roof. However, the curved coal wall is affected by compression and shear function. If the coal side support is not installed in time, it is difficult to strengthen the damaged coal wall, and the rib spalling cannot be controlled effectively.
6. Conclusions
The UDEC numerical simulation method cannot characterize the expansion characteristic of the rock layers caving. That is, the subsidence height of model’s top is equal to the excavation height of the coal seam. For the UDEC, the main roof has an excessive rotated deformation and drives the retained roadway roof to deform much larger than the practical conditions. It also takes a long time for the model to reach equilibrium. The proposed method presented here using a double-yielding constitutive model successfully solved these, and the proposed method takes less than 2/5 of the former’s time to reach equilibrium. Our results are consistent with the filed measurements. Especially when the deformation of retaining roadway was studied, the numerical simulation results are in a good agreement with the field measurements. The new method is also highly efficient and reliable.
The results of the numerical simulation using our new method and the field observation indicate the following. (1) Close relationship exists between the surrounding rock deformation, the width of retaining roadway, and the coal side support, which can affect roof subsidence in varying degrees, specifically as follows. The large width is the major cause of roof’s larger falling and the central open cracks; increasing the supporting force in the middle of the roof significantly reduces the generation of the open cracks in the middle of the roof; the top of the curved coal wall is the major place of compression and shear, and it is necessary to reinforce the top of the curved coal wall in time to prevent big deformation of the roof and the rib spalling. (2) Under the condition of gob-side entry retention, most of the roof cracks distribute inclined. In the abnormal retaining roadway, the bed separation and dip angle of cracking are large in the anchorage area. Below the anchorage area, the roof has a little or even no bed separation and a relatively small crack inclination as well.
Data Availability
The 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 conflicts of interest.
Acknowledgments
This research was funded by the National Key Research and Development Program (grant nos. 2018YFC0603705 and 2016YFC0600901), the National Natural Science Foundation of China (grant no. 51674265), and the Yue Qi outstanding scholar award program of China University of Mining and Technology, Beijing.