In recent years, the rapid growth of the coal–electricity–chemical industry has resulted in significant economic progress. However, this growth has also led to the accumulation of large quantities of coal-based solid waste, such as gasification slag, gangue, fly ash, and slag, which pose challenges for major coal production areas. The disposal of these waste materials has caused land degradation and environmental deterioration, hindering the balanced development of coal resource extraction and environmental protection in mining areas (Zhang et al., 2021; Li et al., 2022). This situation contradicts the principles of ecological civilization, which emphasizes the harmonious coexistence between humans and nature, as well as sustainable development. Filling mining technology offers significant technical advantages in preventing and controlling rock instability, reducing solid waste emissions, and improving coal recovery rates. As a result, filling mining plays a crucial role in achieving green mining practices.

At present, many deep mines are actively exploring new paths that can not only protect the surface environment and buildings but also achieve safe and efficient mining. The “mining filling” integrated technology is an effective method to solve the aforementioned problems (Tu et al., 2021). The coordinated layout of underground coal mining and filling space in deep mining and filling integrated mines can fully leverage the technical advantages of the two types of mining processes and simultaneously meet the engineering needs of rock migration control, rock burst prevention and control, roadway surrounding rock control, water resource protection, and coal seam pressure relief and permeability enhancement (Zhang et al., 2019). At present, many domestic experts and scholars have conducted extensive research on the optimization layout of underground coal mining and filling space. For example, Zhang Qiang et al. studied the solid intelligent filling mining method in coal mines, providing the concept, characteristics, and difficulties of the solid intelligent filling mining method. Combining mining geological conditions, static size parameters of filling equipment mechanisms, and dynamic trajectory parameters during operation, a typical interference state judgment and coordination equation for the mechanism during the filling process was constructed, forming a mechanism for regulating the interference state of the filling equipment mechanism (Zhang et al., 2022); Wu Yabin et al. used numerical simulation methods to analyze the impact of overlying rock movement and surface deformation under three different mining and filling sequences in mines. The results showed that as the mining depth increased, the difference in surface deformation caused by the mining and filling sequence decreased (Wu et al., 2020). Based on theoretical analysis and practical experience, Sun Kai et al. derived formulas for calculating the key parameters of mining and filling process coordination, such as the number of separated tunnels between two alleys, the number of simultaneously filled tunnels in one group, and the number of circulating mining and filling tunnel groups. From this, they proposed a design method for coordinating and optimizing the mining and filling process of longwall tunnels by tunnel cementation filling (Sun et al., 2019); Ma Liqiang et al. proposed the “parallel mining and filling” and “continuous mining and filling” filling and water conservation coal mining methods (Ma et al., 2017). During the mining process, all the mining branches within the entire mining block are first divided into multiple mining stages. Once all the mining branches within one stage are mined and filled, the next stage of mining and filling of the mining branches is carried out until all the mining branches are mined and filled, ultimately achieving complete replacement of the filling body and coal; Lu Bin et al. proposed the short-wall gangue cementation filling technology and studied the material ratio, support tunnel parameters, evolution of surrounding rock stress field, and surface deformation laws (Lu et al., 2017); Ma Liqiang and Li Yongliang proposed the wall-type continuous mining and water filling coal mining method. First, the entire mining block was designed with a longwall system for coal transportation, material transportation, and ventilation systems. Afterward, the mining roadway (wide roadway) within the block was divided into several mining stages, during which skip mining was carried out on the mining roadway. A filling slurry underground mixing and ground mixing conveying system suitable for different coal seam dip angles was proposed, as well as corresponding working face compaction filling methods, forming a complete set of continuous mining and filling cementitious filling coal mining technology (Ma Li et al., 2019; Li et al., 2022); Wang Feifei and Ji Yu established an AHP-Fuzzy evaluation model for the safety status of mining areas and a safety analysis model of “mining area → overlying rock and surface.” They used the elastic foundation beam theory to establish a differential equation for the deflection of the roof, clarifying the relationship between the elastic foundation coefficient of the filling body and the direct elastic modulus of the top and the subsidence of the roof (Jiang et al., 2018; Ji et al., 2023). Tang Weijun, Chang Qingliang, and Sun Xiaoguang used theoretical analysis, numerical simulation, and on-site measurement methods to study the “activation” mechanism of the overlying rock, the law of roof movement, and the factors affecting the stability of the overlying rock when filling and recovering coal pillars. They also predicted the range of damage and water inrush of the bottom plate during paste filling mining, revealing the mechanical mechanism of controlling the bottom plate damage during paste filling mining. The numerical simulation software FLAC was used. We studied the surface subsidence and distribution of elastic–plastic zones on the roof and floor after mining the lower layer of coal with different filling body strengths and analyzed the rock pressure and roof stability of the working face during filling mining (Sun et al., 2007; Chang et al., 2016; Tang et al., 2017). Zhou et al. analyzed the effects of raw material mix ratio, pore size, porosity, and surface structure on the performance of the filling body by a single-factor test. The influence of the mine water environment was simulated, the immersion and dry–wet cycle tests were carried out on the field filling materials, and corrosion was carried out with different durations. Then, the influence of mine water corrosion on the mechanical properties of filling materials was studied by the uniaxial compression test. In order to improve the shear and tensile properties of the filling material, three different lengths of glass fiber were added to the filling material. It is concluded that the average peak strength of the cement-based backfill material increases with the increase in the length of the glass fiber. (Zhou et al., 2021a; Zhou et al., 2021b; Zhou et al., 2022). A green thixotropic cement filling material containing molybdenum tailings was prepared by Gao et al. A thixotropic coefficient method based on the correlation between shear stress and time was proposed to evaluate the thixotropy of the paste. The rheological properties, hydrophilicity, segregation rate, and compressive strength of the oxyphilic paste filling material were studied by the orthogonal test. (Gao et al., 2023). Chen et al. designed experiments based on the response surface method to explore the effects of gasification slag, cement, and water on the performance of paste filling materials (Chen et al., 2022). Based on the response surface method test design and data analysis, Sun Guoqing et al. simulated the roof subsidence law of the paste filling face under different filling body stiffnesses, top coal stiffnesses, support stiffnesses, and filling body strengths by using the finite difference software UDEC and analyzed the influence of single-factor and multi-factor interactions on the roof stability of the paste filling face and its reasons (Sun et al., 2016). Li et al. explored the impact of backfill mining on rock burst prevention and control and used a comprehensive research method of theoretical analysis, numerical simulation, and on-site monitoring to study the influence of controlling backfill rate on the stress of backfill working faces and the evolution of overlying rock spatial structure (Li et al., 2023). In order to reasonably and accurately determine the strength required for filling materials to fill the goaf, Han et al. proposed on-site measurement based on sensors and conducted on-site measurements in Jinchuan Mine. Formulas based on the thick plate theory were derived to verify the measurement results (Han et al., 2022). Ma et al. studied the mechanical properties of coal gangue gypsum filling materials through indoor experiments and analyzed the critical conditions for bending and fracture of the coal seam roof. The discrete element numerical simulation software was used to study the process of roof fracturing, and a reasonable filling ratio of 95% was determined to ensure the stability of the roof. At the same time, the movement of overlying strata and the law of rock pressure reduction were studied through numerical simulation and on-site measurements (Ma et al., 2022).

The abovementioned scholars have conducted a detailed analysis on the migration law, fracture characteristics, and deformation mechanism of the goaf roof in the filling working face. However, there is relatively little research on the spatiotemporal continuity of mining and filling replacement in the “mining filling coordination” operation of the working face. This article takes the 110904 working face of Renjiazhuang Mine as the background. Combined with theoretical analysis, the research methods of numerical simulation and on-site experiments have deeply analyzed the reasonable filling step distance that meets the mining and filling needs of the working face, as well as the main control factors that affect the filling effect. Good experimental results have been achieved on site, and the research results are helpful for promoting the efficient coordination and cooperation of underground coal mining and filling systems, guiding the planning and layout of underground mining and filling space, and dynamic adjustment. Realizing the integrated, safe, efficient, and green mining of “mining selection filling” in deep mines is of positive significance and can provide theoretical and technical support for the formulation of mining filling coordination plans for the same type of working face.

Renjiazhuang Coal Mine employs a mixed development mode, utilizing inclined shafts and vertical shafts. The initial design of the mine divides it into five mining areas: 11, 12, 13, 21, and 23. The mining process starts with the central 11 and 21 mining areas, followed by the connection of the northern 12 mining area and the southern 13 and 23 mining areas. Finally, the mining areas in the west wing of the anticline are extracted. The mine commenced operations in 2008, with the current focus on the 11 mining area, while the 12 mining area remains the continuous mining area. As depicted in Figure 1, mining operations of the 110904 working face in the 9th coal seam of the 11 mining area are being conducted at Renjiazhuang Coal Mine. The 110904 working face is situated adjacent to the already completed 110902 working face on the left and the 110906 working face, which is scheduled for mining by the end of 2023, on the right.

The 9th coal seam of Renjiazhuang Mine is 4 m thick, mainly comprising dark coal. It is semi-dark, powdery, and occasionally lumpy and has a uniform structure and weak asphalt luster. The thickness of the coal seam is generally stable. There is a phenomenon of local thinning and thickening. The coal seam contains gangue, and the coal seam structure is more complicated. The uniaxial compressive strength of the coal seam is 1.83 MPa. After more than 0.05 m of the 9th coal seam is removed from the gangue, the average coal thickness is 3.59 m, and the overall thickness of the coal seam is stable. The 110904 working face is located in the east gate of Renjiazhuang Coal Mine Industrial Square, which is approximately 580.9 m eastward along the 120.5° direction, approximately 264.5 m east–west width, and approximately 1,452 m southward along the 207.5° direction. The ground elevation is + 1,324 m to + 1,328 m. The strike length of the working face is 1,435 m, the dip length is 265 m, and the average dip angle of the coal seam is 14°. The average thickness of the coal seam is 4.2 m, the bulk density is 1.46 t/m³, and the mining method is longwall mining. With reference to the 208 boreholes in the 110901 working face of the adjacent working face of the 110904 working face with the 11th mining area and the 9th coal seam, combined with the on-site technical data of the 110904 working face in Renjiazhuang, the comprehensive histogram of the 110904 working face rock stratum, rock stratum characteristics, and roof and floor conditions are shown in Figure 2.

Based on the theory of beam structure in the mechanics of materials, combined with the comprehensive histogram of the 110904 working face and the characteristic table of rock strata, the theoretical analysis and calculation of the limit caving step distance are carried out. The stress state of the main roof of the working face is assumed to be the fixed beam breaking model, and the stress analysis diagram of the model is shown in Figure 3.

The normal stress of any point in the clamped beam is as follows: σ = M y J z . (1)

In Formula 1, M represents the section bending moment in kN∙m, y represents the distance to the central axis in meters, and J_z represents the central axis section distance. The section distance of the beam is denoted as follows: J z = h 3 12 . (2)

The stress at each position in the structure is as follows: σ = 12 M y h 3 , (3) τ x y = 3 2 Q z h 2 − 4 y 2 h 3 . (4)

In Formula 3 and Formula 4, σ is the normal stress, MPa; h is the thickness, m; τ _xy is the shear stress, MPa; and Q _z is the shear force, kN.

According to the mechanical characteristics of the fixed beam model, the maximum bending moment and tensile stress are as follows: M Max = − q L 2 12 , (5) σ Max = − q L 2 2 h 2 . (6)

According to the model structure, when the external force reaches the tensile strength of the rock mass, the rock layer begins to break and collapse. Based on this, the following conclusions are drawn: R T = − q L 2 2 h 2 , (7) l c = h 2 R T q . (8)

In Formula 7 and Formula 8, h is the layer thickness, m; R _T is the tensile strength, MPa; and q is the uniform load.

It can be seen from Eq. 8 that the initial fracture step of the main roof is related to the tensile strength of the rock stratum, the thickness of the rock stratum, and the load above it. The determination of the load above the main roof is an important part of the initial breaking step. It is generally believed that the load above the overlying rock of the working face is evenly distributed in the original state (Li et al., 2022; Wang et al., 2022; Song et al., 2023). In addition, combined with the “key layer theory” proposed by academician Qian Minggao, it is believed that two discriminant conditions of the key layer are obtained from the characteristics of the key layer.

The first condition is that if the control range of the first stratum reaches the second stratum, the second stratum is the key stratum of the next stratum, which must satisfy the principle of Formula 8. q n + 1 1 < q n 1 . (9)

In Formula 9, (q _n+1)₁ is the load of the No. 1 layer of rock, kN, when the No. n + 1 layer acts on the 1 layer; (q _n)₁ is the load of the No. 1 layer of rock, kN, when the No. n layer acts on the 1 layer.

The second condition, according to Formula 9, to determine the hard rock layer is that the lower fracture distance must be less than the upper fracture distance. l j < l j + 1 . (10)

In Formula 10, l _j is the breaking distance of layer j.

According to the theory of key strata, it is assumed that n and m are two strata (n > m). If L_n > L_m > L_l, the n layer is the main key layer, and the m and l layers are the sub-key layers. If the j layer does not satisfy Formula 10, the discrimination is continued after recalculation.

According to the current classical surrounding rock control theory, the load on the rock stratum is calculated. According to the principle of composite beam, the load of overlying strata can be calculated by the following formulas: q n = γ n h n , (11) q m n = E n h n 3 γ n h n + γ n + 1 h n + 1 + . . . + γ m h m E n h n 3 + E n + 1 h n + 1 3 + . . . + E m h m 3 . (12)

In Formula 12, q _n is the load of the n_th layer. q _mn is the load on the n_th layer when the m_th layer acts on the n_th layer, where m > n. γ is the volume force, h is the thickness, and E is the elastic modulus.

Through the analysis of the characteristic parameters of the rock strata in Figure 2, the key strata of each rock stratum are distinguished. According to the detection data of the coal mine technology department and the fully mechanized mining team, the basic roof of the target working face in this paper is mudstone. The following calculation can be obtained.

The load of the first layer itself is q₁ = 31.68 kPa, and the effect of the second layer of rock on the first layer is (q₂)₁ = 1.88 kPa, (q₂)₁ < q₁, indicating that the second layer of rock has no effect on the first layer of rock; then, the second layer of mudstone is the first key layer. The load of the second layer (mudstone) is q₂ = 205.02 kPa, the effect of the third and fourth layers on the second layer is (q₃)₂ = 243.4 kPa and (q₄)₂ = 238.0 kPa, (q₄)₂ < (q₃)₂, indicating that the fourth layer has no effect on the second layer, and the fourth layer is the second key layer. Based on the aforestated calculation, the influence of the third layer on the load of the second rock layer should be considered. The load on the main roof is 0.25 MPa. Combining it with Formula 8, the initial ultimate fracture step of the main roof of the 110904 face is calculated to be 26.32 m.

The periodic weighting step of the basic roof is often determined by the cantilever fracture of the basic roof, according to the mechanics of materials, σ = M y J . Here, the maximum bending moment is M Max = q L 2 2 (L is the ultimate span of the cantilever beam), Y is h 2 (h is the thickness of the rock layer), σ is the ultimate tensile strength R _T, L = h R T 3 q , and it is compared with the ultimate span of the main roof, L b = h 2 R T q . The periodic weighting step is equivalent to 1/2.45 of L_b, and after calculation, the periodic caving step of the main roof of the 110904 working face is 10.74 m. According to the operation rules, the cutting depth of each knife in the 110904 working face is 0.865 m. Four knives of coal are cut in one shift, 12 knives of coal can be cut in the periodic caving step, and three shifts of work can be taken. The mining and filling coordination scheme of three-mining and one-filling can be adopted, that is, filling once every three shifts of coal mining is completed (Cui et al., 2021; Wang et al., 2021; Jin et al., 2022; Xu et al., 2023; Zhang et al., 2023).

If the mining and filling step distance is greater than the periodic caving step distance, the filling space will have difficulty in meeting the requirements of the established filling ratio, thus affecting the space–time continuity of the mining and filling cycle. It is difficult to describe the law of roof migration after filling by theoretical calculation, and it cannot be simply considered that the periodic caving step obtained by theoretical calculation is the maximum filling step. Therefore, on the basis of theoretical calculation, it is very important to further study the law of roof migration in filling the working face and optimize the mining–filling synergistic process parameters through numerical simulation.

Based on the previous calculation, the reasonable mining–filling synergy step is analyzed using 3DEC numerical simulation (Cao et al., 2023). The 3DEC software includes a range of constitutive models, with the Mohr–Coulomb failure model being the most suitable for analyzing the linear failure surface of shear failure in rock mechanics. The transformation formula is as follows: f s = σ 1 − σ 3 1 + sin ⁡ φ 1 − sin ⁡ φ + 2 c 1 + sin ⁡ φ 1 − sin ⁡ φ . (13)

In Formula 13, σ ₁ is the maximum principal stress, MPa; σ ₃ is the minimum principal stress, MPa; φ is the angle of internal friction, °; and c is the force of cohesion, MPa.

If f _s < 0, shear yielding occurs.

In order to prevent the Mohr–Coulomb model from losing its meaning, this is simplified by assuming that the yield surface is extended to the σ ₃ = σ ^t region. The minimum principal stress is not greater than the tensile strength. If σ ₃ > σ ^t , tensile yielding begins to occur. σ^t cannot exceed σ₃, and the maximum value σt max is determined by the following formula: σ Max t = c tan ⁡ φ . (14)

In the Mohr–Coulomb model, the post-peak characteristics of the material must be considered, that is, the response of the material after the initial damage. This chapter uses the Mohr–Coulomb model to simulate the dilatancy effect. In addition, the elastic material in this simulation is assumed to be isotropic. Its K and G can best describe the basic characteristics of the material, so regarding these material properties, the transformation formula is as follows: K = E 3 1 − 2 ν , (15) G = E 2 1 + ν , (16) where G is the shear modulus, GPa; E is the elastic modulus, GPa; and ν is the Poisson ratio.

The 3DEC built-in joint constitutive model mainly includes two types, as detailed in Table 1. The contact Coulomb slip model is the most frequently used engineering model.

In summary, this paper selects the Mohr–Coulomb failure criterion as the constitutive of the numerical model, and the joint constitutive selects the surface contact Coulomb slip model.

Based on the field engineering background of the 110904 working face in Renjiazhuang Mine, this paper aims to determine the limit caving step distance and periodic caving step distance of the goaf roof in the working face through theoretical analysis and calculation. Additionally, the influence of the filling ratio and 1d strength of the filling body on the filling effect is analyzed using numerical simulation. Furthermore, a hanging pipe filling test was conducted in the 110904 working face, leading to the following conclusions:

1) The key strata in the goaf roof are the second and fourth layers. The main roof has a limit caving step distance of 26.32 m and a periodic caving step distance of 10.74 m. The mining–filling coordination scheme is initially defined as “three-mining and one-filling,” which satisfies the engineering technical requirements. The increase in filling ratio and 1d strength of the filling body can result in roof subsidence on the filling side, roof subsidence on the unfilled side, lateral deformation of the filling body, and contact stress between the roof and filling body, all showing a downward trend. The filling ratio is positively correlated with the 1d strength of the filling body and the overall filling effect. The impact of changing the filling ratio on the filling effect is greater than that of changing the 1d strength of the filling body. The best filling effect is achieved when the filling ratio is 80% and the early compressive strength of the filling body is 2.5 MPa.