# A Method for Determining the Thicknesses of the Reserved Protective Layers in Complex Goafs Using a Joint-Modeling Method of GTS and Flac3D

Authors: Xin Zhao; Dianshu Liu; Shenglin Li; Meng Wang; Shuaikang Tian; Xiaoming Niu; Lin Mo
Journal: Advances in Civil Engineering
 (2021)

---

## Abstract

In this study, a C-ALS underground cavity scanner was used to detect the shapes of mining goafs. In addition, GTS software was adopted to establish a three-dimensional geological model based on the status of the stopes, geological data, and mechanical parameters of each rock mass and to analyze the roof areas of the goafs. In regard to the morphology of the study area, based on a thin plate theory and the obtained field sampling data, a formula was established for the thicknesses of the reserved protective layers in the goafs. In addition, a formula for the thicknesses of the protective layers in the curved gobs was obtained. The thickness formula of the protective layers was then successfully verified. The detection results showed that the roof shapes of the goafs in the Yuanjiacun Iron Mine were mainly arc-shaped, and the spans of the goafs were generally less than 20 m. The stability of the arc-shaped roofs was found to be greater than that of the plate-shaped roofs. Therefore, by reducing the thicknesses of the protective layers in mining goafs, the ore recovery rates can be increased on the basis of safe production conditions. The formula of the thickness of the security layers obtained through the thin plate theory was revised based on the statistical results of the roof shapes of the goafs and then combined using GTS and FLAC3D. The modeling method successfully verified the stability of the mined-out areas. It was found that the verification results were good, and the revised formula was able to improve the recovery rate of the ore under the conditions of meeting safe production standards. Also, it was found that the revised formula could be used in the present situation. At the same time, it was also determined that the complexity of the rock masses obstructed the full identification of the joints and fissures in the present orebodies. Therefore, it is necessary to incorporate C-ALS underground cavity scanners to regularly observe the shapes of the goafs in order to ensure that stability and safety standards are maintained.


---

## Body

## 1. Introduction

Cave mining practices have the advantages of simple production technology, low production costs, and production safety and are thus widely used in the field of metal ore mining [1]. However, as mining depths increase, the surface subsidence areas continue to expand. This has resulted in a variety of geological disasters [2]. It has been found that, with the depletion of easily accessible resources, an increasing number of metal ore mines have transferred to deep mining methods. However, in order to protect ecological environments and ensure the safety of mining personnel, during the stable rock formation stages, horizontal isolation pillars of sufficient size are required. When designing the sizes of these pillars, it is believed that the larger the thickness is, the safer the area will be and the more beneficial it will be to the ecological environment. However, if the sizes of the isolated pillars become too large, resources will be lost unnecessarily. According to the statistical data, the recovery rates of metal ore mines in China are only approximately 40% [3].Determining the optimal size of the isolated pillars is vital for safe mining practices in metal ore mines. At the present time, the methods used for determining the thicknesses of horizontally separated coal pillars mainly include load transfer line intersection methods [4, 5]; fractured arch theory methods [6]; the KB Rupenieit theory estimation method [7]; the fixed beam theory method [8, 9]; and voussoir beam theory [10]. Sofianos and Kapenis [11] and Diederichs et al. [12] further discussed the critical conditions of roof sliding and instability in goaf. These theories have been used in previous production and research projects and have played major roles in the field’s development. With the increased complexity of production conditions with deeper depths, the boundary conditions, spatial effects, and breaking rules of horizontally isolated mine pillars have also become increasingly complicated. At the present time, the applicability of two-dimensional mathematical models remains limited. It has been found that, for isolated mine pillars with complex shapes, researchers are more often considering using plate and shell theories for their analysis processes [13–15]. Although the applications of the above-mentioned theoretical calculation methods have played certain roles in practice, some shortcomings still remain. One of the main factors affecting the calculations of the thicknesses of the safety layers is the selection of the rock layer mechanical parameters. The other main factors that affect the calculations are the morphology of the rock masses. After caving has been completed, the morphology of the roof areas will be irregular. However, the selections of the rock layer parameters adopt certain reference values. In regard to the shapes of the isolated pillars, the rectangular plate shapes are generally defaulted [16]. As a result, large errors from the actual site may occur, which makes it difficult to generalize the obtained conclusions, and large amounts of mineral resources may potentially be wasted.In view of the aforementioned issues, this study was based on the field data of the Yuanjiacun Iron Mine. The rock layer parameters were measured by comprehensively using both field and laboratory instruments, and a model for predicting the safe thicknesses of isolated mining pillars was established in combination with an elastic plate theory. A cavity scanner was used to first detect the mining goaf and then analyze the roof shape of the goaf. A numerical simulation analysis was performed based on the on-site detection results. Then, based on the obtained results, this study’s prediction model for the safe thicknesses of isolated mining pillars was revised.


## 2. Analysis of the Occurrence Conditions of Mining Goafs

### 2.1. Engineering of the Geological Conditions

The Taiyuan Iron and Steel Yuanjiacun Iron Mine [17] is located in Liangjiazhuang Township, Lanxian County, Lvliang City, Shanxi Province. A highway exists from the mining area to the county seat, which is located 20 km away. It is also situated near the proposed Taixing Railway site, as illustrated in Figure 1. The iron ore deposits of the Yuanjiacun Iron Mine belong to the category of Anshan style sedimentary metamorphic magnetic red mixed iron ore, with a geological grade of TFe 32.05%. The associated beneficial and harmful components are both known to be low and belong to the acidic ore deposits. The mine has the following characteristics: a high degree of exploration; large scale; concentrated reserves; simple structural and hydrogeological conditions; small stripping ratio; and suitability for large-scale open-pit mining procedures. However, due to the many different types of ore in the mining area, as well as the high hardness and complex compositions, it is recommended that the sizes of the protective layers of the goafs should be reasonably reserved. These issues have not yet been resolved.Figure 1
Location of the Yuanjiacun Iron Mine [18].


### 2.2. Morphology Analysis of the Mining Goaf Roof Areas

During actual production activities, goaf roof areas are mainly arc-shaped, with small numbers of semicircular and near-planar shapes observed. At the present time, the thickness analysis methods used to examine the security layers of goaf roof areas assume that the roofs are the same as those examined using plane analysis methods, which is different from the actual situations. Therefore, certain information gaps have occurred. The currently used calculation formulae are not always accurate. Moreover, it has been established that curve-shaped roofs are more stable than the plate-shaped roofs. When using traditional calculation methods, it is assumed that curved goafs will be the same as the plate goafs when the sizes of the goafs are similar. In reality, the safe thicknesses are smaller than those of the plate-shaped goafs [9]. Therefore, if traditional methods are used to calculate the reserved thicknesses, major waste of mineral resources may potentially occur.In order to improve the practicability of the calculation formulae of the safe thicknesses of the goaf roofs, it is planned to further modify the roof shapes of the goaf in mines in order to achieve more reliable results. In this study, the C-ALS underground cavity scanner was used to detect the goafs in the study area for a period ranging from 2012 to 2016. During the examination process, detailed statistical analysis was performed. A total of 145 goafs were scanned, of which approximately 60 were roadway and chamber-type goafs, and approximately 85 were larger goafs. Among the various types of goafs, it was observed that the larger goafs were the most unstable. These large goafs had the highest potential of collapse and were also the types which required urgent research inquiries. The survey results showed that, in the larger goafs, the roofs were mainly arc-shaped. The statistical results showed that the goafs’ spanning ranges were mainly less than 20 m.Figure2 [19] shows a representative three-dimensional shape of a goaf. Figure 2(a) shows the No. 27 goaf, in which the top of the goaf is obviously an irregular semicircular shape. Figure 2(b) illustrates the No. 42 goaf, in which the top of the goaf can be seen to display an irregular arc-shape. Figure 2(c) shows the No. 52 goaf, where the top of the goaf displays an irregular semicircular formation. As can be seen in Figure 2(d), the top of the No. 54 goaf was an inclined plate shape. Finally, Figure 2(e) shows that the roof area of the No. 59 goaf was both inclined and arc-shaped.Figure 2
Morphology of the goaf [19].
(a)(b)(c)(d)(e)


### 2.3. Equipment and Method for Obtaining Rock Samples

In the present study, massive rock samples were obtained, and two tests were conducted in the Rock Mechanics Laboratory of the Changsha Mine Research Institute Co., Ltd. For drilling, cutting, smoothing, and air drying. Then, such basic data as the geometric dimensions and quality of the rock specimens were measured. The main test equipment used was as follows:(1)
A JW/4 drilling prototype(2)
TY-450 rock sample cutting machine(3)
SMD-150 double-end grinder(4)
MTS-815 fully digital hydraulic servo rigidity testing machine, along with some auxiliary supporting equipment(5)
30 t universal material pressure testing machine, along with some auxiliary supporting equipment(6)
TP-1000A electronic balance with an accuracy of 0.01 grams(7)
A point load meterThe rock blocks taken at the study site were first formed into cylindrical cores using a drilling prototype and then made into test specimens of a certain specification using a cutting machine and a grinder. The uniaxial compressive strength test samples had an aspect ratio of approximately 2 : 1. That is to say, the size was Φ 50 × 100 mm. The Brazilian test method has previously been recommended by the International Society of Rock Mechanics for testing rock tensile strength. In China, this method has been listed in the national standard “Standard for test methods of engineering rock mass (GB/T50266-99)” and the “Industry standard specification for rock test of water conservancy and hydropower engineering (SL264-2001).” In China, cylindrical samples with a height diameter ratio of 0.5–1.0 are often used for testing rock tensile strength. Therefore, the Brazilian Method tensile strength test sample was processed into a cylindrical sample with a diameter of approximately 50 mm and a height of approximately 35 mm. The height-to-diameter ratio was approximately 0.7 : 1. During the processing, it was found that the samples were too hard in some places to complete the sampling procedures. In those cases, irregular test blocks were selected at the site for point load testing. This study’s processed samples are shown in Figures3 and 4, and the number of samples is detailed in Table 1.Figure 3
Samples processed.
(a)(b)Table 1
Number of samples.
ProjectNo. 1 oxide oreNo. 1 primary oreNo. 2 oxidized oreNo. 10 primary oreNo. 10 magnetiteNo. 10 oxidized oreUniaxial compression test122210—45Brazil split71818435Point load test—153926——Figure 4
test block for partial point load test.
(a)(b)During this study’s experimental research processes, uniaxial compressive strength (UCS), elastic modulus, and Poisson’s ratio tests were completed using a 250-ton fully digital hydraulic servo rigid rock mechanics test system (MTS-815) provided by Rock Mechanics Laboratory of Changsha Mine Research Institute Co., Ltd. The loading test system (MTS-815) is shown in Figure5.Figure 5
MTS-815 test system.The results of the two rock mechanical parameter tests were synthesized, and various rock mechanical parameters were successfully summarized. The final results are listed in Table2.Table 2
Summary table of rock mechanics test results.
Lithologyγ (kN m−3)σ (MPa)σt (MPa)E (GPa)vc (MPa)φ (°)No. 1 oxide ore33.82120.3411.017113.860.24918.2156.33No. 1 primary ore35.0372.736.5232.200.2210.8956.66No. 2 oxidized ore32.54107.8210.37113.840.3116.7255.54No. 10 primary ore36.374182.979.13686.350.21120.4464.81No. 10 magnetite34.43134.377.6073.800.1815.9863.24No. 10 oxidized ore—160.346.21——15.7867.73γ is the bulk density, σ is the compressive strength, σt is the tensile strength, E is the modulus of elasticity, v is Poisson’s ratio, c is the cohesion, and φ is the internal friction angle.


### 2.4. Conversion of the Mechanical Parameters of the Rock Masses

The mechanical parameters obtained in laboratory can be used to estimate the cohesive force. Aydan et al. [20], Hoek and Brown [21], and others have proposed the estimation method. The Simpson Empirical Method is commonly used in China. According to the Simpson Empirical Method, it is believed that there is a relationship between the cohesion Cm of rock masses and the cohesion CR of the complete rock blocks, and the fracture density can be written as follows:(1)Cm=0.114e−0.48i−2+0.02CR.The friction angle in a rock body should be between the friction angle of the complete rock block and the discontinuity in-plane friction angle. Previously, experienced rock mechanic experts had summarized large amounts of experimental data and concluded that the friction coefficient in a complete rock ranges between 1.10 and 1.20. [22] The rock cohesion will be reduced by 10 to 20 times that of the cohesion of the rock mass [22]. In the present study [21, 23], by using the above-mentioned method, the rock mechanical parameters detailed in Table 2 were reduced to the rock mass mechanical parameters, as shown in Table 3.Table 3
Mechanical parameters of rock mass [24–27].
Lithologyc (MPa)φ (°)γ (kN/m3)E (MPa)μC (MPa)φ (°)σt (MPa)RMRSimpsonExperience reductionRMRHoek–BrownExperience reductionNo. 1 oxide ore0.3–0.40.400.9135°–45°26.87°46.9°33.8283940.2490.5537.90.50No. 1 primary ore0.2–0.30.240.5425°–35°24.07°47.2°35.0353680.220.3437.10.30No. 2 oxidized ore0.3–0.40.360.8435°–45°26.87°46.3°32.5483940.310.5237.70.48No. 10 primary ore0.3–0.40.531.0235°–45°29.76°54.0°36.37111930.210.6341.30.57No. 10 magnetite0.3–0.40.410.8035°–45°29.76°52.7°34.43111930.180.5240.80.48No. 10 oxidized ore0.3–0.40.410.7935°–45°29.76°56.4°36.37111930.210.5242.10.48


## 2.1. Engineering of the Geological Conditions

The Taiyuan Iron and Steel Yuanjiacun Iron Mine [17] is located in Liangjiazhuang Township, Lanxian County, Lvliang City, Shanxi Province. A highway exists from the mining area to the county seat, which is located 20 km away. It is also situated near the proposed Taixing Railway site, as illustrated in Figure 1. The iron ore deposits of the Yuanjiacun Iron Mine belong to the category of Anshan style sedimentary metamorphic magnetic red mixed iron ore, with a geological grade of TFe 32.05%. The associated beneficial and harmful components are both known to be low and belong to the acidic ore deposits. The mine has the following characteristics: a high degree of exploration; large scale; concentrated reserves; simple structural and hydrogeological conditions; small stripping ratio; and suitability for large-scale open-pit mining procedures. However, due to the many different types of ore in the mining area, as well as the high hardness and complex compositions, it is recommended that the sizes of the protective layers of the goafs should be reasonably reserved. These issues have not yet been resolved.Figure 1
Location of the Yuanjiacun Iron Mine [18].


## 2.2. Morphology Analysis of the Mining Goaf Roof Areas

During actual production activities, goaf roof areas are mainly arc-shaped, with small numbers of semicircular and near-planar shapes observed. At the present time, the thickness analysis methods used to examine the security layers of goaf roof areas assume that the roofs are the same as those examined using plane analysis methods, which is different from the actual situations. Therefore, certain information gaps have occurred. The currently used calculation formulae are not always accurate. Moreover, it has been established that curve-shaped roofs are more stable than the plate-shaped roofs. When using traditional calculation methods, it is assumed that curved goafs will be the same as the plate goafs when the sizes of the goafs are similar. In reality, the safe thicknesses are smaller than those of the plate-shaped goafs [9]. Therefore, if traditional methods are used to calculate the reserved thicknesses, major waste of mineral resources may potentially occur.In order to improve the practicability of the calculation formulae of the safe thicknesses of the goaf roofs, it is planned to further modify the roof shapes of the goaf in mines in order to achieve more reliable results. In this study, the C-ALS underground cavity scanner was used to detect the goafs in the study area for a period ranging from 2012 to 2016. During the examination process, detailed statistical analysis was performed. A total of 145 goafs were scanned, of which approximately 60 were roadway and chamber-type goafs, and approximately 85 were larger goafs. Among the various types of goafs, it was observed that the larger goafs were the most unstable. These large goafs had the highest potential of collapse and were also the types which required urgent research inquiries. The survey results showed that, in the larger goafs, the roofs were mainly arc-shaped. The statistical results showed that the goafs’ spanning ranges were mainly less than 20 m.Figure2 [19] shows a representative three-dimensional shape of a goaf. Figure 2(a) shows the No. 27 goaf, in which the top of the goaf is obviously an irregular semicircular shape. Figure 2(b) illustrates the No. 42 goaf, in which the top of the goaf can be seen to display an irregular arc-shape. Figure 2(c) shows the No. 52 goaf, where the top of the goaf displays an irregular semicircular formation. As can be seen in Figure 2(d), the top of the No. 54 goaf was an inclined plate shape. Finally, Figure 2(e) shows that the roof area of the No. 59 goaf was both inclined and arc-shaped.Figure 2
Morphology of the goaf [19].
(a)(b)(c)(d)(e)


## 2.3. Equipment and Method for Obtaining Rock Samples

In the present study, massive rock samples were obtained, and two tests were conducted in the Rock Mechanics Laboratory of the Changsha Mine Research Institute Co., Ltd. For drilling, cutting, smoothing, and air drying. Then, such basic data as the geometric dimensions and quality of the rock specimens were measured. The main test equipment used was as follows:(1)
A JW/4 drilling prototype(2)
TY-450 rock sample cutting machine(3)
SMD-150 double-end grinder(4)
MTS-815 fully digital hydraulic servo rigidity testing machine, along with some auxiliary supporting equipment(5)
30 t universal material pressure testing machine, along with some auxiliary supporting equipment(6)
TP-1000A electronic balance with an accuracy of 0.01 grams(7)
A point load meterThe rock blocks taken at the study site were first formed into cylindrical cores using a drilling prototype and then made into test specimens of a certain specification using a cutting machine and a grinder. The uniaxial compressive strength test samples had an aspect ratio of approximately 2 : 1. That is to say, the size was Φ 50 × 100 mm. The Brazilian test method has previously been recommended by the International Society of Rock Mechanics for testing rock tensile strength. In China, this method has been listed in the national standard “Standard for test methods of engineering rock mass (GB/T50266-99)” and the “Industry standard specification for rock test of water conservancy and hydropower engineering (SL264-2001).” In China, cylindrical samples with a height diameter ratio of 0.5–1.0 are often used for testing rock tensile strength. Therefore, the Brazilian Method tensile strength test sample was processed into a cylindrical sample with a diameter of approximately 50 mm and a height of approximately 35 mm. The height-to-diameter ratio was approximately 0.7 : 1. During the processing, it was found that the samples were too hard in some places to complete the sampling procedures. In those cases, irregular test blocks were selected at the site for point load testing. This study’s processed samples are shown in Figures3 and 4, and the number of samples is detailed in Table 1.Figure 3
Samples processed.
(a)(b)Table 1
Number of samples.
ProjectNo. 1 oxide oreNo. 1 primary oreNo. 2 oxidized oreNo. 10 primary oreNo. 10 magnetiteNo. 10 oxidized oreUniaxial compression test122210—45Brazil split71818435Point load test—153926——Figure 4
test block for partial point load test.
(a)(b)During this study’s experimental research processes, uniaxial compressive strength (UCS), elastic modulus, and Poisson’s ratio tests were completed using a 250-ton fully digital hydraulic servo rigid rock mechanics test system (MTS-815) provided by Rock Mechanics Laboratory of Changsha Mine Research Institute Co., Ltd. The loading test system (MTS-815) is shown in Figure5.Figure 5
MTS-815 test system.The results of the two rock mechanical parameter tests were synthesized, and various rock mechanical parameters were successfully summarized. The final results are listed in Table2.Table 2
Summary table of rock mechanics test results.
Lithologyγ (kN m−3)σ (MPa)σt (MPa)E (GPa)vc (MPa)φ (°)No. 1 oxide ore33.82120.3411.017113.860.24918.2156.33No. 1 primary ore35.0372.736.5232.200.2210.8956.66No. 2 oxidized ore32.54107.8210.37113.840.3116.7255.54No. 10 primary ore36.374182.979.13686.350.21120.4464.81No. 10 magnetite34.43134.377.6073.800.1815.9863.24No. 10 oxidized ore—160.346.21——15.7867.73γ is the bulk density, σ is the compressive strength, σt is the tensile strength, E is the modulus of elasticity, v is Poisson’s ratio, c is the cohesion, and φ is the internal friction angle.


## 2.4. Conversion of the Mechanical Parameters of the Rock Masses

The mechanical parameters obtained in laboratory can be used to estimate the cohesive force. Aydan et al. [20], Hoek and Brown [21], and others have proposed the estimation method. The Simpson Empirical Method is commonly used in China. According to the Simpson Empirical Method, it is believed that there is a relationship between the cohesion Cm of rock masses and the cohesion CR of the complete rock blocks, and the fracture density can be written as follows:(1)Cm=0.114e−0.48i−2+0.02CR.The friction angle in a rock body should be between the friction angle of the complete rock block and the discontinuity in-plane friction angle. Previously, experienced rock mechanic experts had summarized large amounts of experimental data and concluded that the friction coefficient in a complete rock ranges between 1.10 and 1.20. [22] The rock cohesion will be reduced by 10 to 20 times that of the cohesion of the rock mass [22]. In the present study [21, 23], by using the above-mentioned method, the rock mechanical parameters detailed in Table 2 were reduced to the rock mass mechanical parameters, as shown in Table 3.Table 3
Mechanical parameters of rock mass [24–27].
Lithologyc (MPa)φ (°)γ (kN/m3)E (MPa)μC (MPa)φ (°)σt (MPa)RMRSimpsonExperience reductionRMRHoek–BrownExperience reductionNo. 1 oxide ore0.3–0.40.400.9135°–45°26.87°46.9°33.8283940.2490.5537.90.50No. 1 primary ore0.2–0.30.240.5425°–35°24.07°47.2°35.0353680.220.3437.10.30No. 2 oxidized ore0.3–0.40.360.8435°–45°26.87°46.3°32.5483940.310.5237.70.48No. 10 primary ore0.3–0.40.531.0235°–45°29.76°54.0°36.37111930.210.6341.30.57No. 10 magnetite0.3–0.40.410.8035°–45°29.76°52.7°34.43111930.180.5240.80.48No. 10 oxidized ore0.3–0.40.410.7935°–45°29.76°56.4°36.37111930.210.5242.10.48


## 3. Mechanical Analysis of the Recommended Safety Thicknesses of Goaf Protective Layers

### 3.1. Instability Assumptions and a Mechanical Equation for Goaf Protective Layers

Mining activities are carried out in three-dimensional spaces. However, after orebodies have been mined, open-field stripe or space pillar mining techniques are often adopted. The suspended roof rock layers are supported by the surrounding pillars, and a fixed three-dimensional board is formed above the active mining area. The stability of stope roofs can be assessed by the strength calculation results of the slab structures [28]. In the present study, the width of the stope roof was denoted by Lx; the length was Ly; the thickness was indicated as h; and the load concentration was q, as shown in Figure 6.Figure 6
Stress analysis of the fixed roof.Before a top plate breaks, the four sides of the visible plate are fixed supports. Therefore, in accordance with the theory of thin plates with small deformations in elastic mechanics, a rectangular thin plate with wideLx and length Ly will be fixedly supported by the uniformly distributed load q. Since the displacement will be 0, it can be assumed that the deflection function of the midplane of the plate will be as follows [14]:(2)wx,y=∑m∑nwmn41−−1mcos2mπxLx1−−1ncos2nπyLy,m=n=1,3,5,....The Rayleigh-Ritz Method is a mathematical variational problem of functional extremes [29, 30]. It is known to be an approximate method for solving mechanical problems using the principle of minimum potential energy or the principle of minimum residual energy. Therefore, for a rectangular thin plate fixed on four sides, the strain energy formula will be as follows:(3)U=D2∬∇2ω2dxdy.The total potential energy will be(4)I=U−∬qωdxdy=D2∬∂2ω∂x2+∂2ω∂y22dxdy−∬qωdxdy.Then, by substituting equation (2) into equations (3) and (4), the following can be obtained:(5)∂I∂Cn=0.Subsequently, the undetermined coefficientwmn can be obtained by equation (5).Then, in order to simplify the calculation, it can be assumed thatm = n = 1, and then(6)w11=qLx4Ly4Dπ43Lx4+Ly4+2Lx2Ly2.Therefore, by substituting equation (6) into equation (2), the following deflection function can be obtained:(7)w=qLx4Ly4Dπ43Lx4+Ly4+2Lx4Ly4cos2πxLxcos2πyLy,where D=Eh3/121−v2 is the bending stiffness of the plate and E and v represent the elastically deformed modulus and Poisson’s ratio of the plate material, respectively.According to the theory of elasticity, the relationship between the stress and the deflection in the plate can be deduced as follows:(8)σx=Ez1−v2∂2w∂x2+v∂2w∂y2,σy=Ez1−v2∂2w∂y2+v∂2w∂x2,τxy=Ez1−v∂2w∂x∂y.Then, by substituting equation (7) into equation (8), the stress expression in the plate can be obtained:(9)σx=ALy2cos2πxLxcos2πyLy+vLx2cos2πyLycos2πxLx,σy=ALx2cos2πyLycos2πxLx+vLy2cos2πxLxcos2πyLy,τxy=1−v2Asin2πxLxsin2πyLy,where A=24Lx2Ly2qz/π23Lx4+Ly4+2Lx2Ly2h3.The stress value at any point in the plate can be obtained using equation (9). The calculation of the maximum principal stress in the plate will be as follows:(10)σ1=12qLx2Ly2Ly2+vLx2π23Lx4+Ly4+2Lx2Ly2h2,σ2=12qLx2Ly2Lx2+vLy2π23Lx4+Ly4+2Lx2Ly2h2,σ3=0.


### 3.2. Instability Criteria for the Goaf Protective Layers

According to the H Tresca Yield Criterion, when shear yielding occurs at a dangerous point on a roof, the principal stress at that point will satisfy the following formula:(11)σ1−σ3=2τmax.Therefore, by substituting equation (10) into the above formula, the following can be obtained:(12)τmax=6Lx2Ly2Ly2−Lx21−vqπ23Lx4+Ly4+2Lx2Ly2h2,where τmax represents the maximum shear stress in the roof rock layer, MPa; H is the thickness of the overburden acting on the roof, m; h denotes the thickness of the roof rock layer, m; q indicates the overburden rock load; and Lx and Ly represent the width and length of the study area, respectively, where Lx = min (Lx, Ly).Therefore, by using equation (7), whether or not the roof can withstand the shear force can be determined. The obtained value can then be compared with the tensile strength and shear strength of the rock masses of the roof. Subsequently, if one of the items reaches or exceeds the strength of the rock masses, it can be assumed that the roof will break. The fracture criterion of a stope roof can be written as follows:(13)maxσx,σy>σT,(14)τmax>τ.In the formula, [σ] and [τ] indicate the tensile strength and shear strength of the rock masses, MPa.


### 3.3. Calculations and Analyses of the Recommended Safe Thicknesses of the Protective Layers in Mining Goafs

In accordance with the known occurrences of mined-out areas in the Yuanjiacun Iron Mine, this study’s analysis process considered that the maximum load which the security layers must bear included the gravity of the overburden and the open-pit mining processes. This study included the following equipment: (1) WK20 electric shovel (850 t); (2) Bellas 7530 cars (371 ∗ 2 t); and (3) D8 N bulldozer (38 t). The calculations were obtained using the following formula:(15)q=γh+GS,=2.44MPa,where γ indicates the overburden layer bulk density, 34 kN/m3; H denotes the thickness of overlying strata, 70 m; G indicates the ground equipment weight, 1.63 × 106 N; and S is the minimum goaf area, 25 m2.Then, by substituting equation (15) into equations (13) and (14), the maximum principal stress of the roof of the goaf when security layers with different side lengths and different thicknesses were left could be analyzed.In this study, in accordance with the listed mechanical parameters of various rock masses under different span conditions, the thicknesses of the various rock mass security layers in the study area were determined as shown in Table4.Table 4
Calculation results of the theoretical analysis method.
Rock characterSpan of goaf (m)101520253035404550No. 1 oxide ore4.56.99.412.115.018.221.524.828.5No. 1 primary ore5.68.511.715.018.722.626.831.536.0No. 2 oxidized ore5.27.910.713.716.820.223.827.231.0No. 10 primary ore4.06.18.310.813.416.319.322.525.8No. 10 magnetite4.06.28.511.113.916.920.123.527.2No. 10 oxidized ore4.67.09.612.415.318.521.825.429.0Then, in accordance with the obtained data shown in Table3, combined with equations (13) and (14), the relationships between the safe thicknesses of the goaf roofs and the sizes of the goafs under the different lithology conditions were successfully determined in the Yuanjiacun Iron Mine (Figure 7).Figure 7
Relationship between the safe thickness of the goaf roof and the span of the goaf.The data shown in Table4 were calculated using equations (13) and (14). However, the calculations obtained using equations (13) and (14) were found to be complex and inconvenient to use on-site. At the same time, as can be seen in Figure 7, that the thicknesses of the roofs and the spans were in a function of one variable, and the fitting relationship was good. These results were found to have good practical value. The fitting equation of roof thicknesses and spans was as follows:(16)No.1oxide ore:h=0.6b−2.34,No.1primary ore:h=0.76b−3.24,No.2oxidized ore:h=0.64b−1.99,No.10primary ore:h=0.54b−2.33,No.10magnetite:h=0.57b−2.77,No.10oxidized ore:h=0.61b−2.37.


## 3.1. Instability Assumptions and a Mechanical Equation for Goaf Protective Layers

Mining activities are carried out in three-dimensional spaces. However, after orebodies have been mined, open-field stripe or space pillar mining techniques are often adopted. The suspended roof rock layers are supported by the surrounding pillars, and a fixed three-dimensional board is formed above the active mining area. The stability of stope roofs can be assessed by the strength calculation results of the slab structures [28]. In the present study, the width of the stope roof was denoted by Lx; the length was Ly; the thickness was indicated as h; and the load concentration was q, as shown in Figure 6.Figure 6
Stress analysis of the fixed roof.Before a top plate breaks, the four sides of the visible plate are fixed supports. Therefore, in accordance with the theory of thin plates with small deformations in elastic mechanics, a rectangular thin plate with wideLx and length Ly will be fixedly supported by the uniformly distributed load q. Since the displacement will be 0, it can be assumed that the deflection function of the midplane of the plate will be as follows [14]:(2)wx,y=∑m∑nwmn41−−1mcos2mπxLx1−−1ncos2nπyLy,m=n=1,3,5,....The Rayleigh-Ritz Method is a mathematical variational problem of functional extremes [29, 30]. It is known to be an approximate method for solving mechanical problems using the principle of minimum potential energy or the principle of minimum residual energy. Therefore, for a rectangular thin plate fixed on four sides, the strain energy formula will be as follows:(3)U=D2∬∇2ω2dxdy.The total potential energy will be(4)I=U−∬qωdxdy=D2∬∂2ω∂x2+∂2ω∂y22dxdy−∬qωdxdy.Then, by substituting equation (2) into equations (3) and (4), the following can be obtained:(5)∂I∂Cn=0.Subsequently, the undetermined coefficientwmn can be obtained by equation (5).Then, in order to simplify the calculation, it can be assumed thatm = n = 1, and then(6)w11=qLx4Ly4Dπ43Lx4+Ly4+2Lx2Ly2.Therefore, by substituting equation (6) into equation (2), the following deflection function can be obtained:(7)w=qLx4Ly4Dπ43Lx4+Ly4+2Lx4Ly4cos2πxLxcos2πyLy,where D=Eh3/121−v2 is the bending stiffness of the plate and E and v represent the elastically deformed modulus and Poisson’s ratio of the plate material, respectively.According to the theory of elasticity, the relationship between the stress and the deflection in the plate can be deduced as follows:(8)σx=Ez1−v2∂2w∂x2+v∂2w∂y2,σy=Ez1−v2∂2w∂y2+v∂2w∂x2,τxy=Ez1−v∂2w∂x∂y.Then, by substituting equation (7) into equation (8), the stress expression in the plate can be obtained:(9)σx=ALy2cos2πxLxcos2πyLy+vLx2cos2πyLycos2πxLx,σy=ALx2cos2πyLycos2πxLx+vLy2cos2πxLxcos2πyLy,τxy=1−v2Asin2πxLxsin2πyLy,where A=24Lx2Ly2qz/π23Lx4+Ly4+2Lx2Ly2h3.The stress value at any point in the plate can be obtained using equation (9). The calculation of the maximum principal stress in the plate will be as follows:(10)σ1=12qLx2Ly2Ly2+vLx2π23Lx4+Ly4+2Lx2Ly2h2,σ2=12qLx2Ly2Lx2+vLy2π23Lx4+Ly4+2Lx2Ly2h2,σ3=0.


## 3.2. Instability Criteria for the Goaf Protective Layers

According to the H Tresca Yield Criterion, when shear yielding occurs at a dangerous point on a roof, the principal stress at that point will satisfy the following formula:(11)σ1−σ3=2τmax.Therefore, by substituting equation (10) into the above formula, the following can be obtained:(12)τmax=6Lx2Ly2Ly2−Lx21−vqπ23Lx4+Ly4+2Lx2Ly2h2,where τmax represents the maximum shear stress in the roof rock layer, MPa; H is the thickness of the overburden acting on the roof, m; h denotes the thickness of the roof rock layer, m; q indicates the overburden rock load; and Lx and Ly represent the width and length of the study area, respectively, where Lx = min (Lx, Ly).Therefore, by using equation (7), whether or not the roof can withstand the shear force can be determined. The obtained value can then be compared with the tensile strength and shear strength of the rock masses of the roof. Subsequently, if one of the items reaches or exceeds the strength of the rock masses, it can be assumed that the roof will break. The fracture criterion of a stope roof can be written as follows:(13)maxσx,σy>σT,(14)τmax>τ.In the formula, [σ] and [τ] indicate the tensile strength and shear strength of the rock masses, MPa.


## 3.3. Calculations and Analyses of the Recommended Safe Thicknesses of the Protective Layers in Mining Goafs

In accordance with the known occurrences of mined-out areas in the Yuanjiacun Iron Mine, this study’s analysis process considered that the maximum load which the security layers must bear included the gravity of the overburden and the open-pit mining processes. This study included the following equipment: (1) WK20 electric shovel (850 t); (2) Bellas 7530 cars (371 ∗ 2 t); and (3) D8 N bulldozer (38 t). The calculations were obtained using the following formula:(15)q=γh+GS,=2.44MPa,where γ indicates the overburden layer bulk density, 34 kN/m3; H denotes the thickness of overlying strata, 70 m; G indicates the ground equipment weight, 1.63 × 106 N; and S is the minimum goaf area, 25 m2.Then, by substituting equation (15) into equations (13) and (14), the maximum principal stress of the roof of the goaf when security layers with different side lengths and different thicknesses were left could be analyzed.In this study, in accordance with the listed mechanical parameters of various rock masses under different span conditions, the thicknesses of the various rock mass security layers in the study area were determined as shown in Table4.Table 4
Calculation results of the theoretical analysis method.
Rock characterSpan of goaf (m)101520253035404550No. 1 oxide ore4.56.99.412.115.018.221.524.828.5No. 1 primary ore5.68.511.715.018.722.626.831.536.0No. 2 oxidized ore5.27.910.713.716.820.223.827.231.0No. 10 primary ore4.06.18.310.813.416.319.322.525.8No. 10 magnetite4.06.28.511.113.916.920.123.527.2No. 10 oxidized ore4.67.09.612.415.318.521.825.429.0Then, in accordance with the obtained data shown in Table3, combined with equations (13) and (14), the relationships between the safe thicknesses of the goaf roofs and the sizes of the goafs under the different lithology conditions were successfully determined in the Yuanjiacun Iron Mine (Figure 7).Figure 7
Relationship between the safe thickness of the goaf roof and the span of the goaf.The data shown in Table4 were calculated using equations (13) and (14). However, the calculations obtained using equations (13) and (14) were found to be complex and inconvenient to use on-site. At the same time, as can be seen in Figure 7, that the thicknesses of the roofs and the spans were in a function of one variable, and the fitting relationship was good. These results were found to have good practical value. The fitting equation of roof thicknesses and spans was as follows:(16)No.1oxide ore:h=0.6b−2.34,No.1primary ore:h=0.76b−3.24,No.2oxidized ore:h=0.64b−1.99,No.10primary ore:h=0.54b−2.33,No.10magnetite:h=0.57b−2.77,No.10oxidized ore:h=0.61b−2.37.


## 4. Amendments to the Security Layers of the Goaf Roofs

### 4.1. Corrections for Recommended Safe Thicknesses of the Goaf Protection Layers

In accordance with the previous research results, numerical simulations were performed on the rock mass mechanical parameters of the No. 1 primary orebody in the study area. The goaf roofs were arc-shaped, flat, and semicircular at the different spans of 10 m, 20 m, and 30 m, respectively. The circle and safety factors of the same roof thicknesses were calculated and analyzed. In this study, when the goaf span was 10 m and the roof thickness was 10 m (greater than the thickness of the required security layer), the cloud diagram of the displacement maps of the different goaf forms was as shown in Figure8. Also, when the span of the area was 20 m and the thickness of the roof was 12.5 m (approximately equal to the thickness of the required security layer), the cloud diagram of the displacement maps of the roof of the different goaf forms resembled that shown in Figure 9. When the span of the goaf was 30 m and the thickness of the roof was 15 m (less than the required security layer thickness), the cloud diagram of the displacement maps of the different goaf forms was as shown in Figure 10. Among the cloud diagrams, Figures 8(a), 9(a), and 10(a) were curved roofs; Figures 8(b), 9(b), and 10(b) were plate-shaped roofs; and Figures 8(c), 9(c), and 10(c) were semicircular top plate roofs.Figure 8
Cloud displacement map of the goaf with a span of 10 m and a roof thickness of 10 m.
(a)(b)(c)Figure 9
Cloud displacement map of the goaf with a span of 20 m and a roof thickness of 12.5 m.
(a)(b)(c)Figure 10
Cloud displacement map of the goaf with a span of 30 m and a roof thickness of 15 m.
(a)(b)(c)This study completed a comparison of the displacement maps shown in Figures8–10. The results revealed that when the thicknesses of the protective layers in the goafs were greater than the minimum thickness of the protective layer (Figure 8), the top plate displacements of the three forms were within the safe range. The displacements of the slab-shaped gob were found to be the largest. This was followed by the semicircular shaped roofs, and the displacements of the arc-shaped slabs were the smallest. When the remaining thicknesses of the gob protection layer in the goaf were approximately equal to the minimum thickness of the required protection layer (Figure 9), the displacements of the slab-shaped gob were the largest. The semicircular type roofs were second, and the displacements of the curved roofs were the smallest. In those cases, the roofs of the plate-shaped goafs were found to be more prone to instability. It was also observed that when the thicknesses of the protective layers in the goafs were less than the minimum thicknesses of the protective layer (Figure 10), the displacements at the top of the roofs and the semicircular roofs were large, and danger of instability existed in the goafs. The displacements of the arc roofs were maintained within a small range. In summary, when compared with the flat roof types, it was found that the curved-shaped roofs had displayed a better ability to maintain the stability of the goafs. That is to say, when the roofs of the goafs were curved, the thicknesses of the required protective layers were smaller than those required for the goafs with plate-shaped roofs.The minimum safety factors of the recommended thicknesses of the goaf protection layers were obtained by the reduced strength method in Flac2D, the safety factor is the ratio of the thickness of the roof to the depth of the plastic deformation zone, and the statistics are shown in Table 5. The following conclusions were made in accordance with the data shown in the table: (1) Under the same security layer thickness conditions, it was found that, for the curved and semicircle roof types, the formations of the pressure arches in the roof shapes were beneficial to the stability of the goaf roofs. (2) When the heights of the arc roofs were greater than 1/4 of the span of the goaf, the different radians were found to have only small impacts on the safety factors of the goafs, for example, less than 2%. (3) When the span of the goaf was 10 m, the safety factors were increased by 6% when the roof of the goaf was greater than the thickness of the security layer. When the span of the goaf was 20 m and the roof of the goaf was close to the thickness of the security layer, the safety factor was increased by 10%. In addition, when the span of the goaf was 30 m and the roof of the goaf was smaller than the thickness of the security layer, it was determined that the safety factor was increased by 19%.Table 5
Analysis results of safety factors of the roof in different goaf shapes.
Gob spanGoaf patternPlateOvalSemicircle102.272.412.39201.691.871.92301.661.981.97Therefore, it was ascertained from the statistical results that the roof of the goafs of Yuanjiacun Iron Mine in Taigang was mainly arc-shaped. The current method of retaining layer protection usually involves calculations based on slab shapes. However, when compared with the actual production situation in the study area, that type of calculation procedure may result in unnecessary losses of mineral resources. Therefore, by revising the equation of the relationship between the safe thickness of the goaf in the roof and the span of the goaf, the recovery rates of the minerals can potentially be improved.


### 4.2. Corrections in the Calculation Method for the Protective Layers of Gobs with Curved Roofs

In the present study, the formula for the thickness of the reserved security layer in gobs with curved roofs was revised according to the statistical results of the goaf roof shapes and the analysis results of the safety factors of the goafs. The rewritten formula was as follows:(17)No.1oxide ore:h=0.57b−0.38,No.1primary ore:h=0.68b−0.51,No.2oxidized ore:h=0.59b−0.35,No.10primary ore:h=0.53b−0.44,No.10magnetite:h=0.56b−0.59,No.10oxidized ore:h=0.58b−0.34,where h represents the thickness of the protective layer of the goaf, m, and b denotes the span of the goaf, m.


## 4.1. Corrections for Recommended Safe Thicknesses of the Goaf Protection Layers

In accordance with the previous research results, numerical simulations were performed on the rock mass mechanical parameters of the No. 1 primary orebody in the study area. The goaf roofs were arc-shaped, flat, and semicircular at the different spans of 10 m, 20 m, and 30 m, respectively. The circle and safety factors of the same roof thicknesses were calculated and analyzed. In this study, when the goaf span was 10 m and the roof thickness was 10 m (greater than the thickness of the required security layer), the cloud diagram of the displacement maps of the different goaf forms was as shown in Figure8. Also, when the span of the area was 20 m and the thickness of the roof was 12.5 m (approximately equal to the thickness of the required security layer), the cloud diagram of the displacement maps of the roof of the different goaf forms resembled that shown in Figure 9. When the span of the goaf was 30 m and the thickness of the roof was 15 m (less than the required security layer thickness), the cloud diagram of the displacement maps of the different goaf forms was as shown in Figure 10. Among the cloud diagrams, Figures 8(a), 9(a), and 10(a) were curved roofs; Figures 8(b), 9(b), and 10(b) were plate-shaped roofs; and Figures 8(c), 9(c), and 10(c) were semicircular top plate roofs.Figure 8
Cloud displacement map of the goaf with a span of 10 m and a roof thickness of 10 m.
(a)(b)(c)Figure 9
Cloud displacement map of the goaf with a span of 20 m and a roof thickness of 12.5 m.
(a)(b)(c)Figure 10
Cloud displacement map of the goaf with a span of 30 m and a roof thickness of 15 m.
(a)(b)(c)This study completed a comparison of the displacement maps shown in Figures8–10. The results revealed that when the thicknesses of the protective layers in the goafs were greater than the minimum thickness of the protective layer (Figure 8), the top plate displacements of the three forms were within the safe range. The displacements of the slab-shaped gob were found to be the largest. This was followed by the semicircular shaped roofs, and the displacements of the arc-shaped slabs were the smallest. When the remaining thicknesses of the gob protection layer in the goaf were approximately equal to the minimum thickness of the required protection layer (Figure 9), the displacements of the slab-shaped gob were the largest. The semicircular type roofs were second, and the displacements of the curved roofs were the smallest. In those cases, the roofs of the plate-shaped goafs were found to be more prone to instability. It was also observed that when the thicknesses of the protective layers in the goafs were less than the minimum thicknesses of the protective layer (Figure 10), the displacements at the top of the roofs and the semicircular roofs were large, and danger of instability existed in the goafs. The displacements of the arc roofs were maintained within a small range. In summary, when compared with the flat roof types, it was found that the curved-shaped roofs had displayed a better ability to maintain the stability of the goafs. That is to say, when the roofs of the goafs were curved, the thicknesses of the required protective layers were smaller than those required for the goafs with plate-shaped roofs.The minimum safety factors of the recommended thicknesses of the goaf protection layers were obtained by the reduced strength method in Flac2D, the safety factor is the ratio of the thickness of the roof to the depth of the plastic deformation zone, and the statistics are shown in Table 5. The following conclusions were made in accordance with the data shown in the table: (1) Under the same security layer thickness conditions, it was found that, for the curved and semicircle roof types, the formations of the pressure arches in the roof shapes were beneficial to the stability of the goaf roofs. (2) When the heights of the arc roofs were greater than 1/4 of the span of the goaf, the different radians were found to have only small impacts on the safety factors of the goafs, for example, less than 2%. (3) When the span of the goaf was 10 m, the safety factors were increased by 6% when the roof of the goaf was greater than the thickness of the security layer. When the span of the goaf was 20 m and the roof of the goaf was close to the thickness of the security layer, the safety factor was increased by 10%. In addition, when the span of the goaf was 30 m and the roof of the goaf was smaller than the thickness of the security layer, it was determined that the safety factor was increased by 19%.Table 5
Analysis results of safety factors of the roof in different goaf shapes.
Gob spanGoaf patternPlateOvalSemicircle102.272.412.39201.691.871.92301.661.981.97Therefore, it was ascertained from the statistical results that the roof of the goafs of Yuanjiacun Iron Mine in Taigang was mainly arc-shaped. The current method of retaining layer protection usually involves calculations based on slab shapes. However, when compared with the actual production situation in the study area, that type of calculation procedure may result in unnecessary losses of mineral resources. Therefore, by revising the equation of the relationship between the safe thickness of the goaf in the roof and the span of the goaf, the recovery rates of the minerals can potentially be improved.


## 4.2. Corrections in the Calculation Method for the Protective Layers of Gobs with Curved Roofs

In the present study, the formula for the thickness of the reserved security layer in gobs with curved roofs was revised according to the statistical results of the goaf roof shapes and the analysis results of the safety factors of the goafs. The rewritten formula was as follows:(17)No.1oxide ore:h=0.57b−0.38,No.1primary ore:h=0.68b−0.51,No.2oxidized ore:h=0.59b−0.35,No.10primary ore:h=0.53b−0.44,No.10magnetite:h=0.56b−0.59,No.10oxidized ore:h=0.58b−0.34,where h represents the thickness of the protective layer of the goaf, m, and b denotes the span of the goaf, m.


## 5. Validation of Calculation Methods for the Goaf Protection Layers

### 5.1. GTS and FLAC3D Joint-Modeling Method

FLAC3D is widely used in the study of mine disasters, including roadway support, roof stability, surface subsidence, and other issues [31–34]. In the present study, in order to verify the reliability of formula (16), a combined modeling method of GTS and FLAC3D was adopted in conjunction with the No. 27 goaf. The specific modeling process is shown in Figure 11. First, the goaf shape was detected using a C-ALS underground cavity scanner. Then, GTS software was used to establish a three-dimensional geological model, according to the status of the stope, geological data, and mechanical parameters of each rock mass. After the modeling had been completed, the node information and unit information in the grid model was read and saved as a .txt file. Then, the fish language was adopted to read the saved node and unit information of the model into the computer’s memory. At the same time, it was also saved as a FLAC3D file containing the node and unit data. The boundary conditions of the model were displacement boundary conditions. The displacements of the left, right, and lower sides of the model were set to 0, and the upper part was the open displacement boundary conditions. After defining parameters such as mesh shape and size, the software was adapted to model division. This method is able to efficiently process complex and irregular models. Finally, the grid data were read through a FLAC3D grid import interface, and a FLAC3D mesh model was generated.Figure 11
GTS and FLAC3D joint-modeling steps.


### 5.2. Parameter Selection

The material properties required by FLAC3D were divided into two types: elastic deformation properties and strength properties. The determination of the material properties was also a key part of the modeling process, as well as a difficult point in the modeling. These were not the same as the material properties of the rock masses caused by certainty.In the field of the underground engineering of surrounding rock during mining processes, the testing of rock material properties is generally performed under laboratory conditions. The differences between the laboratory results and the field data are generally obvious. Therefore, the data available for the modeling processes of numerical simulations are very limited. It is recommended that the material property parameters used for modeling should thereby be determined through the derivation of the existing material properties and based on the actual situations of a site and some previous experience. In addition, the following two are commonly used in the calculations and derivation of bulk modulus and shear modulus. The formula is as follows:(18)K=E31−2v,G=E21+v,where v and E are Poisson’s ratio and the elastic modulus of the rock mass.According to the existing geological data of the examined mine, combined with the results of the previous mechanical tests of the rock masses, stability calculations of the mined-out areas were carried out in this study. In the numerical simulations, the characteristic parameters of the material were based on the experimental data. The model mainly considered four types of lithology: orebodies; No. 2 mining site; schist; and diabase.The selection of physical and mechanical parameters of each rock mass is shown in Table3. The mechanical parameters used for model calculation are shown in Table 6.Table 6
Physical and mechanical parameters of rock mass.
LithologyTensile strength (MPa)Cohesion (MPa)Internal friction angle (°)Density (kN m−3)Young’s modulus (GPa)Poisson’s ratioNo. 10 primary ore0.480.5242.132.6511193.00.17Chlorite schist0.520.6731.928.3610000.00.20Diabase0.821.0535.029.4612000.00.21


### 5.3. Stability Analysis of the No. 27 Goaf

The No. 27 goaf was located in the No. 10 orebody, and there were four north-south rock interlayers near the goaf. The lithology was diabase and schist, and the inclination of the rock layer was approximately 70°. Considering the difficulty of perforation construction, it was finally achieved during the blasting treatments of the empty No. 27 goaf roof at the 1,680 m level. The actual maximum mining range was 24.7 m. According to equation (16), the required safe thickness of the roof of the goaf of the No. 10 orebody could be calculated as follows: h = 0.59 ∗ 24.7–0.49 = 14 m. This was measured after mining at the ground elevation of 1,695 m, with the thickness of the roof of the No. 27 goaf ranging between 26.4 and 33.6 m (average of approximately 30 m). At the ground elevation of 1,680 m, the safety thickness of the roof of the No. 27 goaf was between 11.4 and 18.6 m, with an average thickness of 15 m. The empty areas and rock mass model are shown in Figure 12. The basic size parameters of the No. 27 goaf are detailed in Table 7.Figure 12
Model of goaf and rock mass model.Table 7
Basic parameters of No. 27 goaf.
Maximum span (m)Maximum height (m)Highest elevation of the goaf (m)Lowest elevation of goaf (m)Shadow area (m2)Volume (m3)24.715.71668.61647.4386.72647.0In the present study, a 3D solid model was generated using SURPAC software according to the 3D laser scanning point cloud images, and GTS was used to divide the grid to form the model shown in Figure13. Then, numerical simulation software was used to analyze the stability of the roof of the No. 27 goaf.Figure 13
The FLAC3D model of No. 27 goaf.It was determined that the top of the goaf roof was prone to damages according to the previous research on the stability of the goaf. Therefore, based on the results of the three-dimensional laser scanning, this study compared and calculated the stability of the steps at the 1,680 m and 1,695 m points. The displacement cloud maps (Figure14), stress cloud maps (Figure 15), and plastic area cloud maps (Figure 16) of the area were obtained. Among the figures, Figure 14(a) details a displacement cloud map of a typical north-south profile at the ground level of 1,695 m, and Figure 14(b) details a displacement cloud map of a typical east-west profile at the ground level of 1,695 m. Figure 15(a) details a stress cloud diagram of a typical north-south profile at a ground elevation of 1,680 m, and Figure 15(b) displays a stress cloud diagram of a typical north-south profile at a ground elevation of 1,695 m. Figure 16(a) details a cloud map of a plastic failure zone at a typical elevation of 1,680 m below the ground, while a cloud map of a plastic failure zone in a typical north-south profile at 1,695 m below ground level is shown in Figure 16(b).Figure 14
Displacement cloud diagram of No. 27 goaf. (a) Displacement cloud map of typical north-south profile at the ground level of 1,695 m. (b) Displacement cloud map of typical east-west profile at the ground level of 1,695 m.
(a)(b)Figure 15
Stress cloud diagram of No. 27 goaf. (a) Principal stress cloud diagram of a typical north-south cross section at a ground elevation of 1,680 m. (b) Principal stress cloud diagram of a typical north-south cross section at a ground elevation of 1,695 m.
(a)(b)Figure 16
Cloud diagram of plastic failure in goaf 27. (a) Cloud map of a typical north-south profile plastic failure area at a ground elevation of 1,680 m. (b) Cloud map of a typical north-south profile plastic failure area at a ground elevation of 1,695 m.
(a)(b)After analyzing the results of the numerical simulation calculations (Figure14), it was found that, after the formation of the goaf, the maximum subsidence area of the overburden layer was mainly concentrated at the top of the goaf, within the level of 1695 m. The overall maximum settlement was less than 1 mm, and the displacement was very small. These findings indicated that the roof area of the goaf was generally stable. As can be seen in Figure 14, the position of the floor of the goaf displayed a significant displacement, mainly as a bottom drum. The analysis results at the 1,695 m level (Figure 14) revealed that the maximum kick drum volume was less than 2.5 mm in both states, which was considered to be a stable state.It can be seen from the analysis of the stress cloud diagram (Figure15) that, after the formation of the goaf, compressive stress concentrations occurred at both the top and bottom of the goaf, and the maximum compressive stress was 2,500 kPa and 20,028 kPa, respectively. The location of the tensile stress concentration appeared to be at the top of the goaf. The maximum tensile stress of the roof at the level of 1,695 m in the goaf was 550 kPa. The maximum tensile stress of the roof at the level of 1,680 m in the goaf was 498 kPa, as shown in Figure 15(a). When comparing the mechanical parameters of the rock masses (Table 3), it was determined that the maximum tensile stress was in the goaf. The 1,695 m level was close to the tensile strength of the rock masses (Figure 15(b)), and at the 1,680 m level of the goaf, it was much smaller than the tensile strength of the rock masses. The maximum compressive stress was observed to be mainly concentrated on the sides of the goaf, when compared with that of the rock masses detailed in Table 3. It was indicated from the body mechanics parameters that the maximum value was far smaller than the compressive strength of the goaf and would not cause compression failure of the rock masses. In addition, it was found that there was localized tensile stress on the roof of the goaf close to the tensile strength of the rock masses. However, the range was smaller, which met the requirements of production safety.In the analysis results of the cloud images of plastic failure in the No. 27 goaf (Figure16), it was observed that when the local surface elevation was 1,680 m (Figure 16(a)), local plastic failure had occurred at the end of the goaf roof. The failure mode was determined to be tensile-shear failure, and the plastic failure range was approximately 4 m. The integrity of the bottom of the goaf was observed to be good. When the local surface elevation was 1,695 m (Figure 16(b)), there were no obvious plastic failure zones on the roof of the goaf. The local plastic failure occurred at the bottom of the goaf, and the shallow part of the bottom of the goaf experienced tensile failure. The deeper parts were subjected to shear failures, and the range of plastic failure was approximately 2 m. It was also observed that there was a small range of plastic failure at the bottom of the No. 27 goaf. However, it had not affected the overall stability of the goaf. The top section of the goaf had an elevation of 1,680 m. There was a small range of plastic damage at the end of the goaf where the thickness of the roof had reached as high as 30 m. However, the impact was not ideal. Overall, the level of stability of the No. 27 goaf was determined to be high.This study’s three-dimensional numerical simulation analysis results showed that the No. 27 goaf was stable at the ground elevations of 1,695 m and 1,680 m. At the level of 1,680 m, the minimum roof thickness of the No. 27 goaf was 11.4 m. Therefore, with a safety factor of 2 and the upper equipment load, it was basically in a critical state. However, if there were no other unknown unfavorable factors, it was believed that the safety reserve could basically maintain the stability of the roof. In addition, this study had performed the numerical simulations without considering the safety factors and the combined effects of multiple unfavorable factors. It was also possible to break the existing equilibrium state of the goaf. At the same time, the complexity of the rock masses prevented this study from fully identifying the joints and fissures existing in the orebodies, which were also important factors for the instability of the goaf. In actual production activities, it is necessary to incorporate C-ALS underground cavity scanners to regularly observe the shapes of the goafs in order to ensure that safety standards are maintained.


## 5.1. GTS and FLAC3D Joint-Modeling Method

FLAC3D is widely used in the study of mine disasters, including roadway support, roof stability, surface subsidence, and other issues [31–34]. In the present study, in order to verify the reliability of formula (16), a combined modeling method of GTS and FLAC3D was adopted in conjunction with the No. 27 goaf. The specific modeling process is shown in Figure 11. First, the goaf shape was detected using a C-ALS underground cavity scanner. Then, GTS software was used to establish a three-dimensional geological model, according to the status of the stope, geological data, and mechanical parameters of each rock mass. After the modeling had been completed, the node information and unit information in the grid model was read and saved as a .txt file. Then, the fish language was adopted to read the saved node and unit information of the model into the computer’s memory. At the same time, it was also saved as a FLAC3D file containing the node and unit data. The boundary conditions of the model were displacement boundary conditions. The displacements of the left, right, and lower sides of the model were set to 0, and the upper part was the open displacement boundary conditions. After defining parameters such as mesh shape and size, the software was adapted to model division. This method is able to efficiently process complex and irregular models. Finally, the grid data were read through a FLAC3D grid import interface, and a FLAC3D mesh model was generated.Figure 11
GTS and FLAC3D joint-modeling steps.


## 5.2. Parameter Selection

The material properties required by FLAC3D were divided into two types: elastic deformation properties and strength properties. The determination of the material properties was also a key part of the modeling process, as well as a difficult point in the modeling. These were not the same as the material properties of the rock masses caused by certainty.In the field of the underground engineering of surrounding rock during mining processes, the testing of rock material properties is generally performed under laboratory conditions. The differences between the laboratory results and the field data are generally obvious. Therefore, the data available for the modeling processes of numerical simulations are very limited. It is recommended that the material property parameters used for modeling should thereby be determined through the derivation of the existing material properties and based on the actual situations of a site and some previous experience. In addition, the following two are commonly used in the calculations and derivation of bulk modulus and shear modulus. The formula is as follows:(18)K=E31−2v,G=E21+v,where v and E are Poisson’s ratio and the elastic modulus of the rock mass.According to the existing geological data of the examined mine, combined with the results of the previous mechanical tests of the rock masses, stability calculations of the mined-out areas were carried out in this study. In the numerical simulations, the characteristic parameters of the material were based on the experimental data. The model mainly considered four types of lithology: orebodies; No. 2 mining site; schist; and diabase.The selection of physical and mechanical parameters of each rock mass is shown in Table3. The mechanical parameters used for model calculation are shown in Table 6.Table 6
Physical and mechanical parameters of rock mass.
LithologyTensile strength (MPa)Cohesion (MPa)Internal friction angle (°)Density (kN m−3)Young’s modulus (GPa)Poisson’s ratioNo. 10 primary ore0.480.5242.132.6511193.00.17Chlorite schist0.520.6731.928.3610000.00.20Diabase0.821.0535.029.4612000.00.21


## 5.3. Stability Analysis of the No. 27 Goaf

The No. 27 goaf was located in the No. 10 orebody, and there were four north-south rock interlayers near the goaf. The lithology was diabase and schist, and the inclination of the rock layer was approximately 70°. Considering the difficulty of perforation construction, it was finally achieved during the blasting treatments of the empty No. 27 goaf roof at the 1,680 m level. The actual maximum mining range was 24.7 m. According to equation (16), the required safe thickness of the roof of the goaf of the No. 10 orebody could be calculated as follows: h = 0.59 ∗ 24.7–0.49 = 14 m. This was measured after mining at the ground elevation of 1,695 m, with the thickness of the roof of the No. 27 goaf ranging between 26.4 and 33.6 m (average of approximately 30 m). At the ground elevation of 1,680 m, the safety thickness of the roof of the No. 27 goaf was between 11.4 and 18.6 m, with an average thickness of 15 m. The empty areas and rock mass model are shown in Figure 12. The basic size parameters of the No. 27 goaf are detailed in Table 7.Figure 12
Model of goaf and rock mass model.Table 7
Basic parameters of No. 27 goaf.
Maximum span (m)Maximum height (m)Highest elevation of the goaf (m)Lowest elevation of goaf (m)Shadow area (m2)Volume (m3)24.715.71668.61647.4386.72647.0In the present study, a 3D solid model was generated using SURPAC software according to the 3D laser scanning point cloud images, and GTS was used to divide the grid to form the model shown in Figure13. Then, numerical simulation software was used to analyze the stability of the roof of the No. 27 goaf.Figure 13
The FLAC3D model of No. 27 goaf.It was determined that the top of the goaf roof was prone to damages according to the previous research on the stability of the goaf. Therefore, based on the results of the three-dimensional laser scanning, this study compared and calculated the stability of the steps at the 1,680 m and 1,695 m points. The displacement cloud maps (Figure14), stress cloud maps (Figure 15), and plastic area cloud maps (Figure 16) of the area were obtained. Among the figures, Figure 14(a) details a displacement cloud map of a typical north-south profile at the ground level of 1,695 m, and Figure 14(b) details a displacement cloud map of a typical east-west profile at the ground level of 1,695 m. Figure 15(a) details a stress cloud diagram of a typical north-south profile at a ground elevation of 1,680 m, and Figure 15(b) displays a stress cloud diagram of a typical north-south profile at a ground elevation of 1,695 m. Figure 16(a) details a cloud map of a plastic failure zone at a typical elevation of 1,680 m below the ground, while a cloud map of a plastic failure zone in a typical north-south profile at 1,695 m below ground level is shown in Figure 16(b).Figure 14
Displacement cloud diagram of No. 27 goaf. (a) Displacement cloud map of typical north-south profile at the ground level of 1,695 m. (b) Displacement cloud map of typical east-west profile at the ground level of 1,695 m.
(a)(b)Figure 15
Stress cloud diagram of No. 27 goaf. (a) Principal stress cloud diagram of a typical north-south cross section at a ground elevation of 1,680 m. (b) Principal stress cloud diagram of a typical north-south cross section at a ground elevation of 1,695 m.
(a)(b)Figure 16
Cloud diagram of plastic failure in goaf 27. (a) Cloud map of a typical north-south profile plastic failure area at a ground elevation of 1,680 m. (b) Cloud map of a typical north-south profile plastic failure area at a ground elevation of 1,695 m.
(a)(b)After analyzing the results of the numerical simulation calculations (Figure14), it was found that, after the formation of the goaf, the maximum subsidence area of the overburden layer was mainly concentrated at the top of the goaf, within the level of 1695 m. The overall maximum settlement was less than 1 mm, and the displacement was very small. These findings indicated that the roof area of the goaf was generally stable. As can be seen in Figure 14, the position of the floor of the goaf displayed a significant displacement, mainly as a bottom drum. The analysis results at the 1,695 m level (Figure 14) revealed that the maximum kick drum volume was less than 2.5 mm in both states, which was considered to be a stable state.It can be seen from the analysis of the stress cloud diagram (Figure15) that, after the formation of the goaf, compressive stress concentrations occurred at both the top and bottom of the goaf, and the maximum compressive stress was 2,500 kPa and 20,028 kPa, respectively. The location of the tensile stress concentration appeared to be at the top of the goaf. The maximum tensile stress of the roof at the level of 1,695 m in the goaf was 550 kPa. The maximum tensile stress of the roof at the level of 1,680 m in the goaf was 498 kPa, as shown in Figure 15(a). When comparing the mechanical parameters of the rock masses (Table 3), it was determined that the maximum tensile stress was in the goaf. The 1,695 m level was close to the tensile strength of the rock masses (Figure 15(b)), and at the 1,680 m level of the goaf, it was much smaller than the tensile strength of the rock masses. The maximum compressive stress was observed to be mainly concentrated on the sides of the goaf, when compared with that of the rock masses detailed in Table 3. It was indicated from the body mechanics parameters that the maximum value was far smaller than the compressive strength of the goaf and would not cause compression failure of the rock masses. In addition, it was found that there was localized tensile stress on the roof of the goaf close to the tensile strength of the rock masses. However, the range was smaller, which met the requirements of production safety.In the analysis results of the cloud images of plastic failure in the No. 27 goaf (Figure16), it was observed that when the local surface elevation was 1,680 m (Figure 16(a)), local plastic failure had occurred at the end of the goaf roof. The failure mode was determined to be tensile-shear failure, and the plastic failure range was approximately 4 m. The integrity of the bottom of the goaf was observed to be good. When the local surface elevation was 1,695 m (Figure 16(b)), there were no obvious plastic failure zones on the roof of the goaf. The local plastic failure occurred at the bottom of the goaf, and the shallow part of the bottom of the goaf experienced tensile failure. The deeper parts were subjected to shear failures, and the range of plastic failure was approximately 2 m. It was also observed that there was a small range of plastic failure at the bottom of the No. 27 goaf. However, it had not affected the overall stability of the goaf. The top section of the goaf had an elevation of 1,680 m. There was a small range of plastic damage at the end of the goaf where the thickness of the roof had reached as high as 30 m. However, the impact was not ideal. Overall, the level of stability of the No. 27 goaf was determined to be high.This study’s three-dimensional numerical simulation analysis results showed that the No. 27 goaf was stable at the ground elevations of 1,695 m and 1,680 m. At the level of 1,680 m, the minimum roof thickness of the No. 27 goaf was 11.4 m. Therefore, with a safety factor of 2 and the upper equipment load, it was basically in a critical state. However, if there were no other unknown unfavorable factors, it was believed that the safety reserve could basically maintain the stability of the roof. In addition, this study had performed the numerical simulations without considering the safety factors and the combined effects of multiple unfavorable factors. It was also possible to break the existing equilibrium state of the goaf. At the same time, the complexity of the rock masses prevented this study from fully identifying the joints and fissures existing in the orebodies, which were also important factors for the instability of the goaf. In actual production activities, it is necessary to incorporate C-ALS underground cavity scanners to regularly observe the shapes of the goafs in order to ensure that safety standards are maintained.


## 6. Conclusions

In this research study, the shapes of the mined-out areas were detected using a C-ALS underground cavity scanner. Then, a three-dimensional geological model was established based on the status of the stopes, geological data, and mechanical parameters of the various rock masses using GTS software. The shapes of the goaf roofs in the mined-out areas were obtained from the analysis results. A thin-layer theory was used to establish a formula for the thicknesses of the reserved protective layers in goafs. In addition, a formula for the thicknesses of the protective layers of the curved roofs was be obtained using a modification process. Finally, the formula for the thicknesses of the protective layers for curved roofs was verified through field practice. The following conclusions are obtained in this study:(1)
A detailed statistical analysis of the detected goafs in the study area for the period ranging from 2012 to 2016 was performed using a C-ALS underground cavity scanner. A total of 145 goafs were scanned, of which about 60 were tunnel-type and chamber-type goafs. There were approximately 85 large-scale goafs identified. The detection results showed that the roofs were mainly arc-shaped, and the spans of the goafs were generally less than 20 m. The stability levels of the arc-shaped roofs were greater than those of the plate-shaped roofs. It was found that the thickness of the roofs improved the ore recovery rates on the basis of ensuring safety during the production processes.(2)
A more concise GTS and FLAC3D joint-modeling method was obtained through the acquired research results. The goaf morphology was obtained using a C-ALS underground cavity scanner, and a three-dimensional geological model was established using GTS software. The saved node and element information of the model was read into the computer memory. At the same time, it is saved as a FLAC3D file containing node and element data. Finally, the grid data were imported through a FLAC3D grid import interface for the purpose of reading the grid data and generating a FLAC3D mesh model.(3)
In the present study, based on the statistical results of the roof shapes of the mined-out areas and the analysis results of the safety factors of the mined-out areas, the formula for obtaining the thicknesses of the reserved security layers in arc-shaped gobs was revised. Then, on-site verifications of the No. 10 orebody were carried out based on the revised formula. The stability of the No. 27 goaf was analyzed using a combined modeling method of GTS and FLAC3D. The results indicated that the No. 27 goaf had good stability. The protective layers which were reserved for the gob in the arc roof thickness formula were as follows:(19)No.1oxide ore:h=0.57b−0.38,No.1primary ore:h=0.68b−0.51,No.2oxidized ore:h=0.59b−0.35,No.10primary ore:h=0.53b−0.44,No.10magnetite:h=0.56b−0.59,No.10oxidized ore:h=0.58b−0.34.(4)
The complexity of the rock masses made it impossible to fully identify the joints and fissures in the orebodies, which were also important factors for the instability of the goaf. Therefore, in actual production activities, it is necessary to incorporate C-ALS underground cavity scanners in order to regularly observe the goaf topography in active mining areas in order to ensure that safe production standards are met.



---

*Source: 8862339-2021-01-31.xml*