There are many investigations that focuses on solid mechanics of fractures in different mining cases (Zou et al.2022; Hannah et al. 2022; Li et al. 2022) Various mathematical methods have been proposed to simulate the height of the pressure relief zone based on the geometric form of mining disc-shaped roof strata, including five novel methods introduced by Abbas et al. (2012). Guo et al. (2009) employed a combination of similar and numerical simulations to determine that the failure modes of overburden strata in coal mining are primarily bending, fracture, and separation. Similarly, Ye et al. (2018) conducted similar simulation experiments to discover that cracks frequently occur at both ends of the goaf and are distributed in a trapezoidal pattern under the influence of coal-rock mining effects. Chen et al. (2022) indicated that the dynamic load can connect the goaf and the separation layer by studying the dynamic load caused by the separation layer. Ma et al. (2019) reported that during the mining process, the overlying strata exhibited characteristics of periodic collapse, and the separation cracks also underwent a cycle of 'generation-development-expansion-continuous-disappearance-compaction'.

Guo et al. (2021) established two corresponding mechanical models for overlying strata deformation and obtained the maximum settlement equation of overlying strata bending basin. Jia et al. (2022) found through similar simulation experiments that the subsidence of overburden rock is nonlinear, and the form of subsidence is asymmetric. Zhang et al. (2010) discovered that the "arch-beam" structure formed by surrounding rock cracks is the main factor causing damage and instability, leading to roof collapse. Ren et al. (2013) found through similar simulation experiments that in the mining process of shallow-buried longwall working faces, the fracture of the main roof can easily extend the fracture zone to the surface. Liu (2011) observed that in the coal seam mining process, the overlying strata collapse begins from the immediate roof and gradually expands to the loose layer.

Zhang et al. (2014) demonstrated that the numerical simulation method can be effectively applied to study the rock fracture migration law. Xiao et al. (2022) determined that after a limited number of periodic weightings, the height of the fault zone no longer changes vertically but only horizontally with mining. Liu et al. (2022) analyzed the subsidence of overburden rock, noting that the overall subsidence trend features significant subsidence in the middle and upper parts and minor subsidence in the lower part, with an asymmetric arch collapse shape.

Based on similarity simulation theory and displacement sensing technology, Zhao et al. (2018) captured the entire process from the gradual development of the abscission layer to catastrophic evolution. Li et al. (2005) determined that the limit span of the roof is 110 m by combining numerical and similar simulations. Zhang et al. (2014) proved the reliability and complementarity between similar simulation experiments and field detection studies. Li et al. (2018) derived the fracture evolution law of floor strata under the influence of coal seam mining through similar simulation experiments. Yang et al. (2009) analyzed the mechanical characteristics of rock strata with different mining thicknesses using similar simulations and determined the displacement change of overburden rock and the distribution law of collapsed rock.

Through the verification of field monitoring data, Coggan John et al. (2012) demonstrated that numerical simulation can provide valuable insights into the behavior of coal mine working face roofs and related reinforcement designs. They also suggested that simulation can confirm the adverse effects of unfavorable geostress direction. Ma et al. (2021) identified the failure characteristics and fracture evolution patterns of overburden rock during the mining process by combining numerical simulation with similar experimental methods. Cai et al. (2007) established the effectiveness of PFC and FLAC software and their coupling method in engineering analysis. Yang et al. (2006) discovered that in PFC2D, the stress–strain relationship of the model is primarily influenced by the elastic modulus and is minimally affected by particle size. Yasitli et al. (2005) reported that coal above the fracture zone cannot be completely shattered or broken into large blocks through numerical simulation and field observation.

This study analyzes the failure mechanisms of strata in the working face by examining the development trends of overburden fractures and the displacement of survey lines and validating the experimental results through numerical simulation. It also investigates the subsequent mining process, which is challenging to achieve in actual experiments. In addition, this research emphasizes conducting specific investigations based on the local conditions of different geological occurrences. It explores the overburden failure of the 40204 working face in a Chinese coal mine, using similar simulation experiments and the discrete element simulation technique. It focuses on the rock mass's failure characteristics by analyzing the roof's stress distribution and fracture development characteristics during deformation and failure. The findings provide insights for identifying air leakage positions in the goaf and predicting subsequent pressure build-ups.

The fully mechanized caving face, 40204, of a coal mine in China, examined in this study, is part of the fourth coal seam. The 40204 working face is the second fully mechanized caving face situated in the southern wing of the 402 mining area. It boasts a strike length of 1871 m, a dip length of 220 m, a recoverable length of 1721 m, and is buried at a depth of approximately 590 m. The average thickness of the coal seam is 11.6 m, with a recoverable thickness of 10.0 m. This working face is configured with one access point and three slopes for transportation, return air, and drainage. To the east of the 40204 working face lies the 40202 goaf, with the goaf of the Jiangjiahe mine field to the south, solid coal to the west, and the protective coal pillar of the fourth coal roadway to the north. The geological cross-section of the 40204 working face is depicted in Fig. 1.

40204 working faces primary condition

The study considers the lithology, thickness, and mechanical characteristics within 50 m above the fourth coal seam through the physical and mechanical properties of the stratum analogous to that of the 40204 working faces found in the mine's geological data, as listed in Table 1.

3.1 Laying experimental model

Based on the principles of similarity theory, accurate experimental results reflecting the actual conditions can only be achieved when the scale model and the prototype share geometric and mechanical parameter similarities (Yang 2005). Given the proportionality between the coal seam thickness and the overlying rock thickness in this experiment, a geometric similarity ratio \(C_{\text{L}}\) of 1:200 was adopted for the model. Based on this ratio, the model measures 190 cm in length, 22 cm in thickness, and 110 cm in height, as illustrated in Fig. 2. During the simulation test, the overburden load was regulated using a load control system.

Actual model diagram

Similar materials primarily consist of fillers (or aggregates) and cements. In this experiment, river sand was selected as the filler for the model and mixed with gypsum in proportion. Calcium carbonate was added to enhance its strength. Layers of rock with varying strengths were created using different material ratios. The ratio for each rock layer was determined based on the mechanical parameters of the rock layer in Table 1. The quantities of sand, calcium carbonate, and gypsum in each layer were calculated according to the layer's bulk density. The water content was approximately half that of the gypsum. These layers were constructed in multiple segments to simulate the stratification phenomenon in thicker rock layers accurately. Table 2 details the specific layout and material usage for each layer.

During the model construction, the model's front and rear restraint guard plates were first attached to the experimental frame sequentially with the layering work. Then, the materials required for each layer were measured according to Table 2, mixed with water, and stirred until homogeneous. The mixture was then poured into the model, leveled, and compacted with a scraper. After completing a layer, a thin layer of mica powder was spread to separate it from the subsequent rock layer, continuing upwards until reaching the designed height. Due to spatial constraints, external loads were applied at the model's top to simulate the pressure of the unreached rock formations.

After the construction, to precisely monitor the overburden rock movement and surface deformation, grid lines were drawn on the model's front using ink, with a grid size of 10 × 10 cm. Reflective sheets, sized 2 × 2 cm, were placed at grid intersections to serve as markers, facilitating displacement observations. A total of ten displacement observation lines were established above the coal seam, monitoring ten data points from bottom to top: A-10.0 cm survey line (10.0 cm above the coal seam), B-20.0 cm survey line, C-30.0 cm survey line, D-40.0 cm survey line, E-50.0 cm survey line, F-60.0 cm survey line (with subsequent lines following in sequence), and each line has 18 measurement points.

3.2 Analysis of overlying rock failure law

The test model was dried approximately 15 days after formation and then excavated. The open-off cut was positioned 15 cm from the model's left boundary to mitigate boundary effects, establishing a boundary-protective coal pillar of 15 cm. The model extended 190 cm in length, with 160 cm being recoverable.

Kong et al. (2019) discovered through rock failure experiments and numerical simulations that the uniaxial compressive strength of mudstone is low, making it susceptible to destruction during the mining process at the working face, leading to roof collapse. The influence range of the abutment pressure on the working face extends approximately 40 m. Therefore, during coal seam mining, excavation proceeded in 20 cm increments, allowing the bedrock to stabilize for 2 h after each excavation. In order to compare the results against the field monitoring measurements, the data used in the subsequent analysis are scaled field data at a 1:200 ratio.

The process of the overlying strata's movement and the roof's subsidence are depicted in Fig. 3.

Roof failure diagram under different excavation lengths

When the working face advances to 40 m, the direct roof does not collapse, and the bending settlement of the overlying strata is hardly noticeable. As the advance reaches 80 m, the direct roof begins to collapse, but the surrounding rock's integrity remains unchanged. Monitoring the subsidence curve of the roof survey line indicates that roof subsidence is almost negligible during this phase.

As the working face progresses to 120 m, the rock beam on both sides of the model roof fractures and begins to show cracks. A distinct delamination fracture occurs approximately 20 m above the coal seam. Concurrently, the roof subsidence curve shows significant changes, with the maximum roof subsidence reaching 11 m, nearly equivalent to the thickness of the coal seam.

When the working face progresses to 160 m, the separation position expands upward, and the fracture extends laterally. At an advancement of 200 m, distinct zoning phenomena become evident in the caving area, categorized into delamination fracture zone, broken fracture zone, and compaction zone. Both the delamination fracture and compaction zones assume a trapezoidal shape. The separation fracture's height develops to 60 m. Roof subsidence curves indicate subsidence occurring below the separation fracture at the caving area's top, with the maximum subsidence remaining unchanged.

Upon excavating the coal seam to 240 m, the highest point of the separation fracture zone develops to 100 m above the coal seam, and the cantilever rock layer in the goaf undergoes periodic breaks. The delamination fracture zone transitions from trapezoidal to arch-shaped, while the compaction zone adopts a tower-like form. Upon excavating to 280 m, the delamination fracture zone rises to 110 m above the coal seam and expands horizontally. At 320 m excavation, the vertical height of the delamination fracture zone ceases to change, extending only horizontally with the continuous coal seam advancement. Roof subsidence curves from various survey lines reveal that roof subsidence decreases as the vertical height increases.

Visible fractures generally do not occur in rock formations above 110 m from the coal seam in the final stages of mining. Consequently, the height of the caving and fracture zones in the bedrock is determined to be 110 m above the coal seam. The distribution of subsidence curves at different advancement distances is approximately symmetrical, with the center point of the subsidence curve progressively shifting forward as the advancing distance increases.

3.3 Distribution law of expansion coefficient and porosity

The dilatancy coefficient is defined as the ratio of the rock mass volume after fracturing to the volume of the rock mass before fracturing.

where \(K\) is bulking coefficient; \(V{\prime}\) is rock mass volume after rock breaking, m³; \(V\) is rock mass volume before rock breaking, m³.

Constraints exist on the strike and tendency of the working face following the collapse of the broken immediate roof, where bulking is predominantly vertical. Thus, the bulking coefficient after the goaf's immediate roof collapses can approximately equal its vertical bulking. This approximation allows for a simplified analysis of the structural changes post-collapse.

where \(h_{n\sim n + 1}{\prime}\) is the vertical distance between two points before rock breaking, m; \(h_{n\sim n + 1}\) is the vertical distance between two points after rock breaking, m.

The difference in distance between the n and n + 1 measuring points reflects the subsidence difference between the two points. Hence, the initial equation can be modified to account for this spatial variation

where \(\Delta w\) is the sinking difference between two points, m; \(E_{n\sim n + 1}\) is n and n + 1 vertical deformation between two points.

The porosity can be expressed as follows:

Based on the definitions provided, the average bulking factor and porosity variation across each material layer in the working face simulation test are calculated using measuring points, as depicted in Figs. 4 and 5. Given that the expansion coefficient and porosity distribution tend to align, the analysis focuses solely on the changes in porosity.

Distribution law of bulking coefficient of each layer after excavation

Porosity distribution law of each layer after excavation

Figure 5 indicates that the maximum porosity of the layer appears at the overlying strata 120 m away from the working face, and the porosity of the rock layer increases first and then decreases with the advance of the working face. Hence, it can be inferred that the height of the fracture zone in the working face is about 120 m, which is also basically consistent with the 110 m obtained by similar simulation tests. In addition, this figure also illustrates that the maximum porosity of the fractured zone can reach more than 0.3, but in most mining cases, the porosity of the fractured zone is between 0.2 and 0.3, which correlates with the bedding distribution of the rock formation area. Furthermore, it can be observed from the interlayer porosity curve at 20 m of the overlying strata that the porosity at both ends of the cantilever beam area is the highest, about 0.1. The remaining areas experience continuous compaction during the mining process, resulting in porosity far less than 0.1.

4.1 Building the model and parameter calibration

When employing the numerical simulation method to address the issue of roof collapse, the modeler needs to make appropriate assumptions, simplifications, calibrations, and verifications of the model to obtain reliable output. The material constitutive model should be capable of simulating strength reduction after rock collapse to enhance collapse analysis (Le et al. 2017). Due to computing power limitations, the specific model cannot be entirely built to the surface during modeling. Hence, this study simulates the rock pressure of the surface by adjusting the rock density at the top of the model. In the earth model proposed by Dziewonski and Anderson (1981), pressure increases by 3 GPa for every 100 km increase in underground rock strata.

where \(P\) is the pressure at different underground depths, MPa;\(\rho\) is the rock stratum density, kg/m³; \(g\) is the acceleration of Gravity Normal, N/kg; \(h\) is the rock stratum height, m.

Considering the buried depth of the coal seam in the 40204 working face, the rock pressure between the surface and the top of the model is obtained, and the density of the upper rock layer can be calculated by applying it to Eq. (5), facilitating the assignment of model attributes for discrete elements.

Since this experiment observes the failure characteristics of the overburden rock in the advancing direction of the working face, the model is established using PFC2D. The macroscopic mechanical parameters of each rock stratum are presented in Table 1. The model is 500 m long, and the height of the rock layer is 142.99 m, adhering to a 1:1 ratio. In the discrete element model, the minimum radius of particles is 300 mm, and the maximum radius is 600 mm. The constitutive model between particles is the Linear Parallel Bond Model. In the final model, the number of particles is 92,172, as depicted in Fig. 6.

Model diagram of 40204 working face before coal seam mining

In applying constitutive contact to the model, it is essential to calibrate the rock strata's parameters at the specific working face. Previous literature and engineering practices (Yoon 2007) have typically employed a 'trial and error' method for parameter calibration. However, this approach demands considerable computing power and proves time-consuming and inefficient. Therefore, this study adopts the parameter fitting formula developed by Feng et al. (2022) through multiple fitting tests to calibrate the parameters. It can derive the microscopic mechanical parameters required for the model by inserting the macroscopic mechanical parameters of the rock into Eq. (6).

The macroscopic parameters of each rock stratum in working face 40204 were input into Eq. (6) to calculate the microscopic parameters for each stratum in the numerical model.

where \(k_{\text{n}}\) is the micro-normal stiffness of the model, GPa; \(K_{\text{n}}\) is the macro normal stiffness of the model, GPa.

The micro-normal stiffness values were determined as follows: mudstone at 9.75 GPa for ball-ball and 11.8 GPa for ball-facet; sandstone at 13.3 GPa for ball-ball and 20.8 GPa for ball-facet; and Coal 4 at 10.6 GPa for ball-ball and 10.9 GPa for ball-facet.

where \(k_{\text{s}}\) is the micro tangential stiffness of the model, GPa; \(K_{\text{s}}\) is the macro tangential stiffness of the model, GPa.

The micro-shear stiffness values were found to be: for mudstone, 3.1 GPa for ball-ball and 5.2 GPa for ball-facet; for sandstone, 12.9 GPa for ball-ball and 107.5 GPa for ball-facet; and for Coal 4, 3.1 GPa for ball-ball and 3.4 GPa for ball-facet.

Based on these values, the micro-elastic modulus of the constitutive model was calculated using Eq. (8), with the specific modulus distribution illustrated in Fig. 7.

where \(\overline{\varphi }\) is the micro-internal friction angle of the model, º; \(\varphi\) is the macro internal friction angle of the model, º.

Constitutive model diagram with different rock elastic modulus

The micro-internal friction angles for mudstone, sandstone, and Coal 4 were calculated using Eq. (9) to be 41°, 50°, and 28°, respectively. Figure 8 displays the distribution of the internal friction angle for different lithologies.

where \(c_{2}\) is the micro cohesion of the model, MPa;

The constitutive model diagram of applying different internal friction angles of rock strata

Based on Eq. (10), the micro-cohesion values for mudstone, sandstone, and Coal 4 were found to be 0.47, 0.62, and 0.41 MPa, respectively.

where \(c_{\text{b}}\) is the micro-bond strength of the model, MPa; \(c\) is the macro bond strength of the model, MPa.

The micro-bond strengths for mudstone, sandstone, and Coal 4, calculated using Eq. (11), were 5.3, 2.3, and 2.6 MPa, respectively. The cohesion distribution for different lithologies is shown in Fig. 9.

The constitutive model diagram of applying different cohesion of rock strata

4.2 Visual analysis of overburden rock movement

The previous section highlights the simulation experiment's limitations, as it can only simulate overburden failure up to a mining depth of 320 m. This study addresses these limitations by employing numerical simulation methods to model overburden movement and stress change as mining advances. This approach confirms the previous experimental results. Figure 9 depicts the roof displacement at various advancing distances.

Figure 10 demonstrates that as the distance advances from 80 to 440 m, the displacement distribution of roof strata in the caving area exhibits a trapezoidal shape. During the mining process, up to 320 m, the displacement of the roof in most areas of the caving strata ranges from 9 to 12 m, with the maximum displacement reaching 13.7 m at the bottom of the caving strata. The maximum discrepancy in displacement between these findings and those from similar simulations varies from approximately 3.4% to 12%. This indicates that the numerical simulation results largely agree with similar simulation outcomes.

Rock displacement diagram under different advance distances

As the mining distance extends from 320 to 440 m, the normal height of the overlying strata migration area ceases to change, continuing to expand horizontally, and the maximum displacement of the roof gradually stabilizes between 12 and 13 m.

Figure 11 illustrates the stress distribution of rock strata under different advancing distances.

The stress distribution diagram of rock strata under different advancing distances

Figures 11, 12, and 13 reveal that during the mining process, from 80 to 120 m, the stress distribution in the displacement area of the overlying strata maintains a trapezoidal shape. The 'arch-beam' structure appears in the overburden displacement area at 160 m. At an advancing distance of 80 m, the minimum and maximum stresses in the roof displacement area are 7.5 and 10 MPa, respectively. When the advancing distance reaches 120 m, these stresses are 7.5 and 12.5 MPa, respectively.

The distribution of the monitoring system

Working face support layout diagram

It is observed from an analysis of stress distribution during the mining process of the working face that the vertical stress distribution in the rock stratum shows a gradual decrease in stress values and a more dispersed stress distribution with increasing distance from the coal seam. This suggests that the stress distribution follows a normal distribution pattern, gradually reducing the dynamic influence area caused by mining in the vertical direction. The distribution of overburden stress in the goaf forms a symmetrical trapezoidal shape, aligning with patterns observed in analogous simulation tests where overburden fissures develop. These fissures originate at the top and propagate laterally on both sides of the working face. In addition, as the working face progresses in the mining process, the stress reduction and dispersion areas tend to re-concentrate during subsequent continuous mining, leading to the re-compaction of fractures in previously mined areas. The accuracy of this analysis is confirmed through Figs. 2f–h.

5.1 Monitoring system and support distribution

The newly established KJ24 coal mine roof and rock burst monitoring system at Dafosi Mine, part of the Binchang Mining Industry, was officially launched. This system enables real-time online monitoring, transmission, intelligent analysis, and early warning of critical data such as roadway roof separation and anchor cable stress. It enhances the mine's safety monitoring capabilities and assurance coefficient, providing robust and scientific data support for optimizing roadway support parameters. The specific layout is depicted in Fig. 12.

The working face is equipped with 150 hydraulic supports, comprising eight head and tail transition supports and 142 middle supports. It is outfitted with 18 pressure sensor monitoring substations within the middle supports, each overseeing four measurement points on a set of hydraulic supports. The measurement points are distributed approximately every seven units from the #4 frame of the working face to the machine's tail. The sequence of frames includes #4, #12, #20, #28, #36, #44, #52, #60, #68, #76, #84, #92, #100, #108, #116, #124, #132, and #140. The support layout is illustrated in Fig. 13. The collected data are then transmitted to the surface computer and analyzed by technicians.

5.2 Analysis of monitoring results

The exported real-time monitoring data undergoes post-processing, and the segments that conform to the specified advancing distance intervals are compiled into a table, as listed in Table 3.

The stress interval diagram is illustrated in Figs. 14 and 15 using the numerical simulation results and the field monitoring data.

Value range of stress arch under different advancing distance

Value range of stress beam under different advancing distance

Analysis of the stress interval diagram reveals that within the overburden rock's failure stress arch interval, the maximum and minimum stress values decrease as the advancing distance increases. In the overburden rock's stress beam interval, the maximum stress value fluctuates periodically with an increase in the advancing distance, while the minimum value consistently declines. Figure 15 compares the real-time monitoring data from the coal mine site with the stress beam interval range. It indicates that the in-field stress fluctuates within the simulated stress range, validating the simulation study's accuracy.

This research innovatively analyzed the failure characteristics of overlying strata in the working face of specific coal mines, utilizing the coupling of field similarity simulation experiments and discrete element numerical simulations. This analysis provides reference for solving engineering problems. The main conclusions are as follows:

1. (1)

   The similarity simulation experiment results indicated that the fracture zone's height in the 40204 working face is approximately 110 m above the coal seam. The fractures predominantly extend laterally, with separation cracks and compaction areas developing in vertical and horizontal directions. In addition, throughout the entire excavation process, the maximum subsidence of the roof stabilizes at 11 m, corresponding to the coal seam's thickness. The porosity in the area where the separation fracture occurred ranges from 0.2 to 0.3, whereas the porosity of the remaining roof is less than 0.1.

2. (2)

   The caving zone at a mining depth of 200 m presents a trapezoidal shape, composed of an 'arch-beam' structure and can be categorized into three zones: 'delamination fracture zone', 'broken fracture zone', and 'compaction zone'. The numerical simulation displacement results indicated that most roof strata areas displace between 9 and 12 m, with a maximum displacement of 13.7 m. Compared to the similarity simulation's roof subsidence results, the error range is between 3.4 and 12%, which aligns closely with the similarity simulation outcomes.

3. (3)

   Failure stress analysis indicated that the stress distribution in the overburden displacement area remains trapezoidal during mining at depths of 80 to 120 m. The overlying strata's displacement area exhibits an 'arch-beam' structure from a mining depth of 160 m. The maximum and minimum values of the failure stress arch interval of the overburden rock decrease as the advancing distance increases. In the stress beam interval of the overburden rock, the maximum value shows periodic changes with the increase in advancing distance, while the minimum value continues to decrease.