# Dimensionless Charts for Predicting the Range of Goaf Roof Caving

Authors: Jing Zhang; Rongxing He; Fengyu Ren; Zhihua Ouyang
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

In engineering, the method of charts can provide a convenient query for specific engineering problems. To provide the basis for the potential hazard evaluation and rational governance of the goaf, it is necessary to study the quantitative evaluation for the range of goaf roof caving. Undoubtedly, the charts used to visually query the caving range can simplify the workload of the quantitative evaluation. Therefore, the methods of dimensional analysis, numerical simulation, and linear interpolation are introduced to study the dimensionless charts for predicting the caving range. The dimensionless analysis is used to establish the fuzzy function relationship among the influence factors of the goaf roof caving, and the numerical simulation is used to calculate the dimensionless groups in the fuzzy function. Using the linear interpolation, the dimensionless charts in this work can predict the range of goaf roof caving under more working conditions. The results show that the characteristics of the goaf roof caving corresponding to the dimensionless curves are consistent with the actual situation. With the continuous increase of the goaf spanl, the dimensionless curves of the caving range experience zero growth, rapid growth, and steady growth. The growth degree varies with the fracture spacing S. Especially in the zero growth phase, the duration of the relatively stable state in the overlying strata of the goaf increases with the increase of fracture spacing S. Moreover, based on the case study of Shirengou Iron Mine, the dimensionless charts obtained in this work can predict the range of goaf roof caving under different working conditions, which indicates the findings of this study have certain guiding significance to the treatment of the goaf.


---

## Body

## 1. Introduction

The goaf is formed by manmade excavation or natural geological movement under the ground surface, which seriously threatens the safety production of the mine [1, 2]. The goaf roof caving, one of the main sources of disasters in the goaf, has been one of the research hotspots in the field of underground excavation engineering. Generally, the backfilling materials are used to fill the goaf to prevent the roof caving [3, 4]. Certainly, the measure of inducing roof caving can manage the mine safety hazards by optimizing mining technology, such as the induced caving method [5, 6], which includes the caving control engineering [7] and the certain safety measures [8]. Therefore, the research on the goaf roof caving has important practical significance for the underground excavation engineering.Regardless of preventing or inducing the roof caving to ensure the safety production in the mine, the range of roof caving needs to be researched in advance. Most scholars [9–11] contributed to the study of the mechanism and law of the roof caving, which laid the foundation for predicting the range of goaf roof caving. Regarding the prediction of the range of goaf roof caving, the research methods can be divided into three categories: theoretical mechanics model, physical experiment, and numerical simulation. In the theoretical mechanics model, the caving arch mechanics model [12–14] was the most widely used. Based on that model, the relationship between critical caving span and caving height can be obtained, which can be used to evaluate the state of the goaf roof caving [15]. However, many physical and mechanical parameters affecting the goaf roof caving were neglected in the obtained relationship. Some scholars [16–18] used the physical experiment to conduct the simulation of the goaf roof caving. The equipment for the physical experiment was generally more systematic, such as the physical modeling system for longwall mining [17] and the electro-hydraulic servo multichannel similar material test bench [19]. Therefore, the method of physical experiment was time consuming and costly. In addition, with the upgrading of numerical simulation software, some scholars [20–22] utilized the numerical simulation to study the range of goaf roof caving. The numerical simulation software mainly used the discrete element software, such as UDEC [23], 3DEC [24], and YADE [25]. The results in the numerical simulation were relatively close to the engineering site. However, the numerical simulation was applicable to analyze the goaf roof caving under a special working condition. According to the current research status of the range of goaf roof caving, the method that can not only comprehensively consider the influence factors of the goaf roof caving but also quantitatively evaluate the range of goaf roof caving under different working conditions is lacking.The goaf roof caving is a relatively complex physical phenomenon, which involves many physical and mechanical parameters, such as the size of the excavation area and the properties of the surrounding rock. Therefore, it is difficult to visually predict the range of goaf roof caving. However, the method of charts can provide a convenient query for specific engineering problems. The intuitive information can be obtained in the charts. For example, the support classification chart of theQ system [26, 27] can provide the support method of the excavation area according to the engineering information such as rock mass classification, excavation height, support spacing, and bolt length. Mathews stability graph [28, 29] can evaluate the stope stability through the stability number and hydraulic radius. The stability charts for rock slopes [30, 31] determine the safety factor of the slope according to the slope stability influencing factor. Based on the application of the above chart methods, this paper attempts to utilize the chart method to predict the range of goaf roof caving.The research strategy on the dimensionless charts for the range of goaf roof caving is shown in Figure1. In the following sections, considering that the dimensional analysis [32] is widely utilized to deal with complex physical phenomena containing multiple parameters, the Buckingham pi theorem [33] in the dimensional analysis is first selected to establish the fuzzy function relationship between the range of roof caving and the influencing factors of the goaf roof caving. Then, the discrete element software UDEC is applied to throw light on the fuzzy function relationship by numerical calculation. The large amounts data of the numerical calculation can obtain the dimensionless charts included the range of roof caving. Finally, based on a case study in the Shirengou Iron Mine, the range of goaf roof caving under the different working conditions can be predicted in the dimensionless charts using the method of linear interpolation.Figure 1
Research strategy on the dimensionless charts for the range of goaf roof caving.


## 2. Buckingham Pi Theorem

The dimensional analysis [32, 34, 35] is widely utilized to study the general law of the complex physical phenomena, which includes Rayleigh’s method and the Buckingham pi theorem. Rayleigh’s method is used to determine the functional relationships of physical phenomena with fewer variables [32]. The Buckingham pi theorem can be used to research the physical phenomena containing multiple physical quantities [36, 37]. Therefore, the Buckingham pi theorem is selected to deal with the influencing factors of the goaf roof caving. Assuming that a physical phenomenon contains n physical quantities, the functional relationship of the physical phenomenon can be expressed as [33](1)fx1,x2,…,xn=0.If the independent physical quantities in the functional relationship arem (m < n), the independent physical quantities can be regarded as basic physical quantities, and the remaining physical quantities in the functional relationship are regarded as nonbasic physical quantities. Thus, the nonbasic physical quantities can be expressed by a compound quantity of the basic physical quantities, and a new function can be combined in the form of (n-m) dimensionless numbers π. Therefore, equation (1) can be expressed as(2)fπ1,π2,…,πn−m=0.Obviously, a decrease in the number of physical quantities in equation (2) reduces the complexity of the physical phenomenon, which simplifies the theoretical analysis and experimental design alike.The fractures in the natural rock mass are complex. Under the effect of joints and fractures, the goaf roof caving in the macroscopic is the tensile stress exceeds the tensile strength or the shear stress exceeds the shear strength. The microscopic performance is the development of rock mass fractures including new fracture development and the extension of original fracture. Therefore, the goaf roof caving is closely related to the fractures. To research the mechanical behavior of fractured rock mass, Warren and Root [38] introduced three groups of orthogonal fracture systems to represent the complex fracture systems. Then, the equivalent model for fractured rock mass can be established, as shown in Figure 2.Figure 2
The equivalent model for fractured rock mass.The equivalent hydraulic conductivityK of a set of parallel fractures in Figure 2 can be expressed as [39, 40](3)K=gB312μkS,where B is the fracture aperture and S is the fracture spacing. The parameters g and μk are the gravitational acceleration and the kinematic viscosity of water, respectively.According to equation (3), the parameters B and S determine the value of hydraulic conductivity K, which can be used to characterize the fracture development degree of rock mass. It indicates that the parameters B and S are closely related to the goaf roof caving. Taking the x-z plane in Figure 2 as the analysis object, as shown in Figure 3, the original rock stress around the goaf is redistributed by the excavation of the initial goaf. The caving of the goaf roof occurs under the continuously increase of the goaf span l. Using the caving height h to represent the range of goaf roof caving, as shown in Figure 3, the parameter h is related to the goaf span l, the goaf initial height hb, and the goaf depth H. The parameter hb is mainly used to control the compensation space of the goaf roof caving. In addition, the rock density ρ and gravity acceleration g determine the bulk weight of rock mass, which is related to the gravity stress of the goaf roof. The rock mass elastic modulus Em and Poisson's ratio ν are also closely related to the excavation and subsequent caving of the goaf.Figure 3
Schematic diagram of influencing factors of the goaf roof caving.According to the analysis of the influence factors of the goaf roof caving in the fractured rock mass, the physical and mechanical parameters involved in the goaf roof caving include the parametersh, l, hb, H, ρ, g, Em, ν, S, and B. The dimension of each parameter is shown in Table 1. The surrounding rock of the goaf in Table 1 can be any kind of rock or ore. The properties of the surrounding rock vary with the working conditions of the goaf. Based on the Buckingham pi theorem, a fuzzy function relationship among the influence factors of the goaf roof caving can be established, which can be expressed as(4)fh,l,hb,H,ρ,g,Em,ν,S,B=0.Table 1
Unit and dimension of physical quantities involved in the goaf roof caving.
NumberPhysical quantitySymbolUnitDimension1Goaf roof caving heighthmL2Goaf spanlmL3Goaf initial heighthbmL4Goaf depthHmL5Rock densityρkg/m−3ML−36Gravity accelerationgm/s2LT−27Rock mass elastic modulusEmkg/(ms2)ML-1T−28Poisson's ratioν—19Fracture spacingSmL10Fracture apertureBmLThe solutions to equation (5) is derived in Appendix and summarized as follows:(5)hhb=flhb,ρgHνEm,BS,where the dimensionless group h/hb represents the extent of the goaf roof caving. The dimensionless group l/hb represents the ratio of width to height in the initial goaf. The dimensionless group ρgHν/Em represents the basic properties of rock mass at a certain goaf depth. The dimensionless group B/S represents the fracture development degree of rock mass.


## 3. Numerical Simulation

The numerical simulation method for the goaf roof caving can be conducted by the finite element software and the discrete element software. However, the finite element software is not applicable to study the large deformation and displacement of the block system. Therefore, the discrete element software is selected to study the goaf roof caving. The Universal Distinct Element Code (UDEC) is a two-dimensional discrete element program for processing discontinuities. It can be used to simulate the response of discontinuities (such as joints and fractures in the rock mass) under static or dynamic loads. In the discrete element software UDEC, the discontinuous surfaces are treated as boundary surfaces among blocks and the large displacement and rotation of the block along the discontinuity surface are allowed. Using this software, most scholars [41–45] have achieved abundant results in the excavation of fractured rock mass. Therefore, the discrete element software UDEC is selected in this work to study the range of goaf roof caving.


### 3.1. Mechanical Parameters

The Mohr–Coulomb plastic model is selected as the constitutive model in the UDEC. The mechanical parameters of rock mass in the discrete element calculation include the rock mass uniaxial tensile strengthσmt, the rock mass cohesion c1, the angle of internal friction of rock mass φ1, the rock mass volume modulus Km, and the rock mass shear modulus Gm, which can be obtained according to the rock mass strength criterion and continuity assumption.


#### 3.1.1. Rock Mass Uniaxial Tensile Strengthσmt

The rock mass uniaxial tensile strengthσmt can be obtained by the Hoek–Brown criterion. The generalized Hoek–Brown rock mass strength criterion can be expressed as [46](6)σ1=σ3+σcmσ3σc+sai,where σ1 and σ3 are the maximum and minimum principal stresses, respectively. The parameter σc is the rock uniaxial compressive strength. The parameters m, s, and ai are constants related to rock mass materials, which can be expressed as(7)m=miexpGSI−10028−14D,(8)s=expGSI−1009−3D,(9)ai=12+16e−GSI/15−e−20/3,where GSI is the geological strength index and D is a factor, which depends upon the degree of disturbance to which the rock mass has been subjected to blast damage and stress relaxation [46]. The parameters GSI and D can be obtained through field investigation. The parameter mi can be selected according to the Hoek–Brown constant charts of rock mass [47].From equation (9), the value of parameter ai is approximately equal to 0.5 when the parameter GSI ≥ 50. Then, substituting σ1 = 0 into equation (6), the rock mass uniaxial tensile strength σmt can be expressed as(10)σmt=σc2m−m2+4s.


#### 3.1.2. Rock Mass Cohesionc1 and Angle of Internal Friction of Rock Mass φ1

The parametersc1 and φ1 can be obtained by the Mohr–Coulomb criterion [48], which states that the relationship between σ1 and σ3 is linear.(11)σ1=1+sinφ11−sinφ1σ3+2c1cosφ11−sinφ1.Equation (11) can be expressed as(12)σ1=σmc+kσ3,where σmc is the rock mass uniaxial compressive strength, σmc = 2c1cosφ1/(1 − sinφ1), and k = (1 + sinφ1)/(1 − sinφ1).The Hoek–Brown criterion and Mohr–Coulomb criterion can predict rock strength well under the lower confining pressure. Hoek and Brown [49] proposed that the predicted results of both criteria are consistent under the condition of 0 <σ3 < 0.25σc. When 0 <σ3 < 0.25σc, the data of the parameters σ1 and σ3 can be calculated using equation (6). Then, using equation (12), the values of σmc and k can be obtained by fitting the data of principal stresses. Under the premise that the value of k is known, the value of φ1 is easy to find.Fromσmc = 2c1cos φ1/(1 − sin φ1),(13)c1=σmc2k.According to equation (13), the rock mass cohesion c1 can be obtained.


#### 3.1.3. Elastic ModulusEm, Volume Modulus Km, and Shear Modulus Gm

The rock mass volume modulusKm and the rock mass shear modulus Gm can be expressed as(14)Km=Em31−2ν,(15)Gm=Em21+ν.The value of Poisson's ratioν is easy to obtain. When the elastic modulus of rock mass Em is obtained, the values of the parameters Km and Gm can be calculated by equations (14) and (15). Therefore, the parameter Em is the key parameter to be determined.The existence of joints and fractures directly affects the integrity of rock mass. Meanwhile, the joints and fractures determine the elastic modulus of rock mass, which the fracture spacingS and fracture aperture B have the most direct effect on the parameter Em. A set of parallel fractures in Figure 3 is selected as the research object (as shown in Figure 4). It is convenient to regard the fractured rock mass as the continuous rock mass to study the mechanical behavior of fractured rock mass. Under the vertical stress Δσ, it can be known from Hooke's law that(16)Δεs=ΔσE,where Δεs is the strain of fractured rock mass and E is the rock elastic modulus.Figure 4
The equivalent continuous model for fractured rock mass.Then, the displacement Δus in the fracture spacing S can be expressed as(17)Δus=SΔεs=SΔσE.Define the deformation of a single fracture is Δub, then [39](18)Δub=Δσkn,where is kn the normal fracture stiﬀness.As shown in Figure4, define the deformation of the continuous rock mass is Δum; then,(19)Δum=Δus+Δub=S+BΔσEm.Substituting equations (17) and (18) into equation (19) gives(20)Em=EknS+BSkn+E.According to equations (6)–(20), the parameters σmt, c1, φ1, Km, and Gm can be calculated. The rock physical and mechanical parameters involved in the calculation can be obtained from the field investigation.


### 3.2. Numerical Calculation Flowchart

As shown in Figure5, the fracture spacing S, the goaf depth H, and the goaf span l are selected as the variables to obtain the dimensionless charts of the goaf roof caving height h. The elastic modulus of rock mass Em can be obtained by equation (20). When the rock physical and mechanical parameters in equations (6)–(20) are known, the parameter h can be obtained by the numerical simulation. Then, the dimensionless groups in equation (5) can be calculated. Subsequently, add the values of dimensionless groups into the dimensionless charts, then one calculation ends. In the next calculation, gradually increase the variables l, H and S for loop calculation. Finally, the whole loop calculation ends until the variable S takes the maximum value.Figure 5
Flowchart of numerical simulation.It is worth noting that the more times of the loop calculation, the higher the accuracy of the prediction of the goaf roof caving height in the dimensionless charts. Based on this point, the flowchart of numerical simulation in Figure5 could increase the complexity of numerical simulation and the complexity of data processing. To overcome this difficult, the flowchart of dimensionless charts is established. As shown in Figure 6, the computer language Python is introduced to call the discrete element software UDEC using the os.system module. And then the numerical calculation results are directly provided to Python. Finally, using the NumPy module and Matplotlib module, the values of dimensionless groups can be calculated and the dimensionless charts can be plotted.Figure 6
Flowchart of dimensionless charts.


## 3.1. Mechanical Parameters

The Mohr–Coulomb plastic model is selected as the constitutive model in the UDEC. The mechanical parameters of rock mass in the discrete element calculation include the rock mass uniaxial tensile strengthσmt, the rock mass cohesion c1, the angle of internal friction of rock mass φ1, the rock mass volume modulus Km, and the rock mass shear modulus Gm, which can be obtained according to the rock mass strength criterion and continuity assumption.


### 3.1.1. Rock Mass Uniaxial Tensile Strengthσmt

The rock mass uniaxial tensile strengthσmt can be obtained by the Hoek–Brown criterion. The generalized Hoek–Brown rock mass strength criterion can be expressed as [46](6)σ1=σ3+σcmσ3σc+sai,where σ1 and σ3 are the maximum and minimum principal stresses, respectively. The parameter σc is the rock uniaxial compressive strength. The parameters m, s, and ai are constants related to rock mass materials, which can be expressed as(7)m=miexpGSI−10028−14D,(8)s=expGSI−1009−3D,(9)ai=12+16e−GSI/15−e−20/3,where GSI is the geological strength index and D is a factor, which depends upon the degree of disturbance to which the rock mass has been subjected to blast damage and stress relaxation [46]. The parameters GSI and D can be obtained through field investigation. The parameter mi can be selected according to the Hoek–Brown constant charts of rock mass [47].From equation (9), the value of parameter ai is approximately equal to 0.5 when the parameter GSI ≥ 50. Then, substituting σ1 = 0 into equation (6), the rock mass uniaxial tensile strength σmt can be expressed as(10)σmt=σc2m−m2+4s.


### 3.1.2. Rock Mass Cohesionc1 and Angle of Internal Friction of Rock Mass φ1

The parametersc1 and φ1 can be obtained by the Mohr–Coulomb criterion [48], which states that the relationship between σ1 and σ3 is linear.(11)σ1=1+sinφ11−sinφ1σ3+2c1cosφ11−sinφ1.Equation (11) can be expressed as(12)σ1=σmc+kσ3,where σmc is the rock mass uniaxial compressive strength, σmc = 2c1cosφ1/(1 − sinφ1), and k = (1 + sinφ1)/(1 − sinφ1).The Hoek–Brown criterion and Mohr–Coulomb criterion can predict rock strength well under the lower confining pressure. Hoek and Brown [49] proposed that the predicted results of both criteria are consistent under the condition of 0 <σ3 < 0.25σc. When 0 <σ3 < 0.25σc, the data of the parameters σ1 and σ3 can be calculated using equation (6). Then, using equation (12), the values of σmc and k can be obtained by fitting the data of principal stresses. Under the premise that the value of k is known, the value of φ1 is easy to find.Fromσmc = 2c1cos φ1/(1 − sin φ1),(13)c1=σmc2k.According to equation (13), the rock mass cohesion c1 can be obtained.


### 3.1.3. Elastic ModulusEm, Volume Modulus Km, and Shear Modulus Gm

The rock mass volume modulusKm and the rock mass shear modulus Gm can be expressed as(14)Km=Em31−2ν,(15)Gm=Em21+ν.The value of Poisson's ratioν is easy to obtain. When the elastic modulus of rock mass Em is obtained, the values of the parameters Km and Gm can be calculated by equations (14) and (15). Therefore, the parameter Em is the key parameter to be determined.The existence of joints and fractures directly affects the integrity of rock mass. Meanwhile, the joints and fractures determine the elastic modulus of rock mass, which the fracture spacingS and fracture aperture B have the most direct effect on the parameter Em. A set of parallel fractures in Figure 3 is selected as the research object (as shown in Figure 4). It is convenient to regard the fractured rock mass as the continuous rock mass to study the mechanical behavior of fractured rock mass. Under the vertical stress Δσ, it can be known from Hooke's law that(16)Δεs=ΔσE,where Δεs is the strain of fractured rock mass and E is the rock elastic modulus.Figure 4
The equivalent continuous model for fractured rock mass.Then, the displacement Δus in the fracture spacing S can be expressed as(17)Δus=SΔεs=SΔσE.Define the deformation of a single fracture is Δub, then [39](18)Δub=Δσkn,where is kn the normal fracture stiﬀness.As shown in Figure4, define the deformation of the continuous rock mass is Δum; then,(19)Δum=Δus+Δub=S+BΔσEm.Substituting equations (17) and (18) into equation (19) gives(20)Em=EknS+BSkn+E.According to equations (6)–(20), the parameters σmt, c1, φ1, Km, and Gm can be calculated. The rock physical and mechanical parameters involved in the calculation can be obtained from the field investigation.


## 3.1.1. Rock Mass Uniaxial Tensile Strengthσmt

The rock mass uniaxial tensile strengthσmt can be obtained by the Hoek–Brown criterion. The generalized Hoek–Brown rock mass strength criterion can be expressed as [46](6)σ1=σ3+σcmσ3σc+sai,where σ1 and σ3 are the maximum and minimum principal stresses, respectively. The parameter σc is the rock uniaxial compressive strength. The parameters m, s, and ai are constants related to rock mass materials, which can be expressed as(7)m=miexpGSI−10028−14D,(8)s=expGSI−1009−3D,(9)ai=12+16e−GSI/15−e−20/3,where GSI is the geological strength index and D is a factor, which depends upon the degree of disturbance to which the rock mass has been subjected to blast damage and stress relaxation [46]. The parameters GSI and D can be obtained through field investigation. The parameter mi can be selected according to the Hoek–Brown constant charts of rock mass [47].From equation (9), the value of parameter ai is approximately equal to 0.5 when the parameter GSI ≥ 50. Then, substituting σ1 = 0 into equation (6), the rock mass uniaxial tensile strength σmt can be expressed as(10)σmt=σc2m−m2+4s.


## 3.1.2. Rock Mass Cohesionc1 and Angle of Internal Friction of Rock Mass φ1

The parametersc1 and φ1 can be obtained by the Mohr–Coulomb criterion [48], which states that the relationship between σ1 and σ3 is linear.(11)σ1=1+sinφ11−sinφ1σ3+2c1cosφ11−sinφ1.Equation (11) can be expressed as(12)σ1=σmc+kσ3,where σmc is the rock mass uniaxial compressive strength, σmc = 2c1cosφ1/(1 − sinφ1), and k = (1 + sinφ1)/(1 − sinφ1).The Hoek–Brown criterion and Mohr–Coulomb criterion can predict rock strength well under the lower confining pressure. Hoek and Brown [49] proposed that the predicted results of both criteria are consistent under the condition of 0 <σ3 < 0.25σc. When 0 <σ3 < 0.25σc, the data of the parameters σ1 and σ3 can be calculated using equation (6). Then, using equation (12), the values of σmc and k can be obtained by fitting the data of principal stresses. Under the premise that the value of k is known, the value of φ1 is easy to find.Fromσmc = 2c1cos φ1/(1 − sin φ1),(13)c1=σmc2k.According to equation (13), the rock mass cohesion c1 can be obtained.


## 3.1.3. Elastic ModulusEm, Volume Modulus Km, and Shear Modulus Gm

The rock mass volume modulusKm and the rock mass shear modulus Gm can be expressed as(14)Km=Em31−2ν,(15)Gm=Em21+ν.The value of Poisson's ratioν is easy to obtain. When the elastic modulus of rock mass Em is obtained, the values of the parameters Km and Gm can be calculated by equations (14) and (15). Therefore, the parameter Em is the key parameter to be determined.The existence of joints and fractures directly affects the integrity of rock mass. Meanwhile, the joints and fractures determine the elastic modulus of rock mass, which the fracture spacingS and fracture aperture B have the most direct effect on the parameter Em. A set of parallel fractures in Figure 3 is selected as the research object (as shown in Figure 4). It is convenient to regard the fractured rock mass as the continuous rock mass to study the mechanical behavior of fractured rock mass. Under the vertical stress Δσ, it can be known from Hooke's law that(16)Δεs=ΔσE,where Δεs is the strain of fractured rock mass and E is the rock elastic modulus.Figure 4
The equivalent continuous model for fractured rock mass.Then, the displacement Δus in the fracture spacing S can be expressed as(17)Δus=SΔεs=SΔσE.Define the deformation of a single fracture is Δub, then [39](18)Δub=Δσkn,where is kn the normal fracture stiﬀness.As shown in Figure4, define the deformation of the continuous rock mass is Δum; then,(19)Δum=Δus+Δub=S+BΔσEm.Substituting equations (17) and (18) into equation (19) gives(20)Em=EknS+BSkn+E.According to equations (6)–(20), the parameters σmt, c1, φ1, Km, and Gm can be calculated. The rock physical and mechanical parameters involved in the calculation can be obtained from the field investigation.


## 3.2. Numerical Calculation Flowchart

As shown in Figure5, the fracture spacing S, the goaf depth H, and the goaf span l are selected as the variables to obtain the dimensionless charts of the goaf roof caving height h. The elastic modulus of rock mass Em can be obtained by equation (20). When the rock physical and mechanical parameters in equations (6)–(20) are known, the parameter h can be obtained by the numerical simulation. Then, the dimensionless groups in equation (5) can be calculated. Subsequently, add the values of dimensionless groups into the dimensionless charts, then one calculation ends. In the next calculation, gradually increase the variables l, H and S for loop calculation. Finally, the whole loop calculation ends until the variable S takes the maximum value.Figure 5
Flowchart of numerical simulation.It is worth noting that the more times of the loop calculation, the higher the accuracy of the prediction of the goaf roof caving height in the dimensionless charts. Based on this point, the flowchart of numerical simulation in Figure5 could increase the complexity of numerical simulation and the complexity of data processing. To overcome this difficult, the flowchart of dimensionless charts is established. As shown in Figure 6, the computer language Python is introduced to call the discrete element software UDEC using the os.system module. And then the numerical calculation results are directly provided to Python. Finally, using the NumPy module and Matplotlib module, the values of dimensionless groups can be calculated and the dimensionless charts can be plotted.Figure 6
Flowchart of dimensionless charts.


## 4. Case Study

### 4.1. Background

The ore deposit in the Shirengou Iron Mine is an Anshan-type magnetite deposit with a surface elevation of +145 m. The ore body is produced in layers with an inclination of 65°–75°. The surrounding rock is hornblende gneiss. Above the −60 m level, the shallow-hole mining method was adopted for the stable ore rock. The shallow-hole mining method had left many goafs of varying sizes, forming an intricate goaf group in space. The open stope subsequent filling method is used to mine the ore body below the −60 m level, and the section height is 15 m. At present, the Shirengou Iron Mine is mined to −210 m level. As shown in Figure7, the roof of the goaf left by the illegal mining induced below the level of −210 m penetrated to the −210 m level and is intersected with the goaf of the operation area. Thus, a large irregular goaf (referred to as M2 main goaf in Figure 7) is formed.Figure 7
Right view of the 3D laser scanning of the M2 main goaf.It is measured that the height of the M2 main goaf is about 120 m, and the maximum width is about 102 m. Considering that the M2 main goaf is in an unstable state, it is necessary to carry out the potential hazard evaluation and rational governance to ensure the safety mining of the surrounding ore body. At this time, a quantitative evaluation for the range of the M2 main goaf caving is of important guiding significance to the safety measures of the M2 main goaf.


### 4.2. Calculation

To predict the range of goaf roof caving under different working conditions, this paper takes the Shirengou Iron Mine as the engineering background. As for the layered ore body, the goaf can be regarded as a 2-D model [50, 51]. The numerical model for the goaf roof caving is established (as shown in Figure 8). The field investigation shows that the fracture spacing S is mainly concentrated in 1∼4 m, the joint occurrence is mostly 234°∠0°, the goaf depth H is mainly concentrated in the range of 200∼450 m, and the goaf span l is 10∼102 m. Therefore, the range of the variables S, H, and l in Figure 5 can be set to 1∼4 m, 100∼500 m, and 0∼100 m, respectively. The excavation step distance of the goaf is set to 2 m, which can obtain massive data of the goaf roof caving.Figure 8
The schematic diagram of numerical simulation.The Mohr–Coulomb criterion is used as the elastic-plastic criterion in the discrete element calculation. The left and right boundaries of the numerical model are subjected to deformation constraints to limit the horizontal displacement. The bottom boundary of the numerical model is subjected to fixed constraints to limit the horizontal and vertical displacement, and the top is the free boundary. The numerical model just considers the gravity stress and neglects the tectonic stress and mining stress. The gravitational accelerationg is 9.8 m/s2. Physical and mechanical parameters of rock mass and joint are obtained by field investigation and relevant equations, which are shown in Table 2. The range of goaf roof caving corresponding to the excavation step distance can be obtained according to the flowchart of numerical simulation (as shown in Figure 5).Table 2
Physical and mechanical parameters of the Shirengou Iron Mine.
ParameterValueMethodGeological strength index GSI76Field investigationDisturbance factorD0.7Rock material constantmi24Gravitational accelerationg (m/s2)9.8Initial goaf heighthb (m)15Hydraulic conductivityK (m/s)1.68 × 10−8Kinematic viscosity of water (20°C)μk (m2/s)1.01 × 10−6Rock densityρ (kg/m3)3.58Rock uniaxial compressive strengthσc (MPa)99.44Rock elastic modulusE (GPa)80.3Poisson's ratioν0.21Joint normal fracture stiﬀnesskn (GPa)9.67Field investigation and simulation experienceJoint tangential fracture stiffnessks (GPa)6.53Joint cohesionc2 (MPa)1.67Joint internal friction angleφ2 (°)22.35Joint tensile strengthσt2 (MPa)0.23Fracture apertureB (m)2.73 × 10−4Equation (3)Rock mass material constantm6.42Equation (7)Rock mass material constants0.031Equation (8)Rock mass material constanta0.5Equation (9)Rock mass tensile strengthσmt (MPa)0.48Equation (10)Rock mass cohesionc1 (MPa)4.97Using equation (11) to fit the principal stress dataRock mass internal friction angleφ1 (°)40.98Fracture spacingS (m)1.0–4.0Field investigationGoaf depthH (m)100–500Goaf spanl (m)0–100.0Rock mass elastic modulusEm (GPa)Equation (20): the value of Em varies with the parameter SRock mass volume modulusKm (GPa)Equation (14): the value of Km varies with the parameter EmRock mass shear modulusGm (GPa)Equation (15): the value of Gm varies with the parameter Em


### 4.3. Results

When the goaf span is excavated to 100 m, the caving ranges of the goaf roof under different fracture spacing are shown in Figures9–12. The results show that the caving height of goaf roof increases with the increase of the goaf depth H under the conditions of the same fracture spacing S and goaf span l. When the goaf depth H and goaf span l are same, the caving height of goaf roof decreases with the increase of the fracture spacing S. Using the goaf roof caving height corresponding to the excavation step distance to calculate the dimensionless group in equation (5), the dimensionless chart for the range of goaf roof caving (as shown in Figure 13) can be obtained.Figure 9
Caving range of the goaf roof (whenl = 100 m, S = 1 m). (a) H = 100 m. (b) H = 200 m. (c) H = 300 m. (d) H = 400 m. (e) H = 500 m.
(a)(b)(c)(d)(e)Figure 10
Caving range of the goaf roof (whenl = 100 m, S = 2 m. (a) H = 100 m. (b) H = 200 m. (c) H = 300 m. (d) H = 400 m. (e) H = 500 m.
(a)(b)(c)(d)(e)Figure 11
Caving range of the goaf roof (whenl = 100 m, S = 3 m. (a) H = 100 m. (b) H = 200 m. (c) H = 300 m. (d) H = 400 m. (e) H = 500 m.
(a)(b)(c)(d)(e)Figure 12
Caving range of the goaf roof (whenl = 100 m, S = 4 m. (a) H = 100 m. (b) H = 200 m. (c) H = 300 m. (d) H = 400 m. (e) H = 500 m.
(a)(b)(c)(d)(e)Figure 13
Dimensionless charts for the caving range of goaf roof. (a)S = 1 m. (b) S = 2 m. (c) S = 3 m. (d) S = 4 m.
(a)(b)(c)(d)As shown in Figure13, the overlying strata of the goaf can remain relatively stable when the goaf span is small. However, the duration of the relatively stable state is closely related to the fracture spacing S. As shown in Figures 13(a)–13(d), the points O1∼O4 are the starting points of the goaf roof caving when the parameter H is 500 m, respectively. Obviously, the degree of difficulty of the goaf roof caving in the initial stage increases with the increase of fracture spacing S, which indicates the duration of the relatively stable state increases with the increase of fracture spacing S.With the continuous increase of the goaf span, the exposed area of the goaf expands gradually. It results in the occurrence of the overlying strata instability under the action of the gravity stress and internal concentrated stress. As shown in Figures13(a)–13(d), taking the dimensionless curve corresponding to the parameter H is 500 m as an example, the caving height of the goaf roof h increases with the increase of the goaf span l. The results show that there is a rapid increase process between the points Oj and Rj (j = 1, 2, 3, 4). The reason for this phenomenon is related to the accumulation of internal strain energy and the continuous development and expansion of fractures during the period when the overlying strata maintain a relatively stable state. The speed of caving can appear short acceleration when the caving of overlying strata is occurred.Under the premise that the goaf spanl is no longer expanded, a relatively stable natural balance arch can be formed when the roof caving reaches a certain height, and the roof caving could not occur for a long time. However, when the span of the goaf is increasing, the stress balance state of the roof rock is destroyed continuously, and the caving could continue to occur. As shown in Figures 13(a)–13(d), the overlying strata still caving with increase of the goaf span when the rapid growth period of the roof caving is over. At this time, the growth rate of the caving height is gradually smaller, but the scale of the caving is relatively large.


### 4.4. Application

The dimensionless charts obtained by the excavation step distance 2 m contain the goaf roof caving under lots of working conditions. Using the linear interpolation method, the range of goaf roof caving under other working conditions can be predicted. Taking the M2 main goaf of −210 m level in the Shirengou Iron Mine as an example, the goaf depthH is 355 m, the goaf span l is 72 m, and the fracture spacing S is 1.4 m. The dimensionless charts corresponding to S = 1 m and S = 2 m in Figure 13 are used to predict the range of goaf roof caving under this working condition.Defineπ1 = l/hb, π2 = h/hb, and π3 = log10(Em/ρgHν). The feature point of the goaf roof caving range predicted by the dimensionless chart is denoted as P(π1, π2, π3). When the goaf depth H is 355 m and the goaf span l is 72 m, as shown in Figure 14, the coordinates of feature points predicted by the corresponding dimensionless charts corresponding to S = 1 m and S = 2 m are P1 (4.8, 3.62, 3.52) and P2 (4.8, 3.08, 3.77), respectively. Using the dimensionless group h/hb, the ranges of goaf roof caving predicted by the two dimensionless charts are 54.3 m and 46.2 m, respectively. It can be known from the linear interpolation that the predicted range of the goaf roof caving is 51.06 m when the parameter S = 1.4 m.Figure 14
Estimation for the caving range of goaf roof. (a)S = 1 m. (b) S = 2 m.
(a)(b)As shown in Figure7, the caving height of the M2 main goaf roof in the −210 m level is 46.66 m, which is less than the range of roof caving predicted by the dimensionless charts in this work. In fact, the M2 main goaf roof is always in an unstable state. Occasionally, the M2 main goaf has the sound of caving rubble according to the field investigation. This phenomenon indicates that the M2 main goaf is in a state of caving activity, and the range of roof caving predicted by the dimensionless charts in this work is reasonable. Therefore, the dimensionless charts obtained by the dimensional analysis coupled with numerical simulation can preliminarily quantitatively evaluate the caving range of goaf roof under different working conditions, which has certain guiding significance to the treatment of the goaf.


## 4.1. Background

The ore deposit in the Shirengou Iron Mine is an Anshan-type magnetite deposit with a surface elevation of +145 m. The ore body is produced in layers with an inclination of 65°–75°. The surrounding rock is hornblende gneiss. Above the −60 m level, the shallow-hole mining method was adopted for the stable ore rock. The shallow-hole mining method had left many goafs of varying sizes, forming an intricate goaf group in space. The open stope subsequent filling method is used to mine the ore body below the −60 m level, and the section height is 15 m. At present, the Shirengou Iron Mine is mined to −210 m level. As shown in Figure7, the roof of the goaf left by the illegal mining induced below the level of −210 m penetrated to the −210 m level and is intersected with the goaf of the operation area. Thus, a large irregular goaf (referred to as M2 main goaf in Figure 7) is formed.Figure 7
Right view of the 3D laser scanning of the M2 main goaf.It is measured that the height of the M2 main goaf is about 120 m, and the maximum width is about 102 m. Considering that the M2 main goaf is in an unstable state, it is necessary to carry out the potential hazard evaluation and rational governance to ensure the safety mining of the surrounding ore body. At this time, a quantitative evaluation for the range of the M2 main goaf caving is of important guiding significance to the safety measures of the M2 main goaf.


## 4.2. Calculation

To predict the range of goaf roof caving under different working conditions, this paper takes the Shirengou Iron Mine as the engineering background. As for the layered ore body, the goaf can be regarded as a 2-D model [50, 51]. The numerical model for the goaf roof caving is established (as shown in Figure 8). The field investigation shows that the fracture spacing S is mainly concentrated in 1∼4 m, the joint occurrence is mostly 234°∠0°, the goaf depth H is mainly concentrated in the range of 200∼450 m, and the goaf span l is 10∼102 m. Therefore, the range of the variables S, H, and l in Figure 5 can be set to 1∼4 m, 100∼500 m, and 0∼100 m, respectively. The excavation step distance of the goaf is set to 2 m, which can obtain massive data of the goaf roof caving.Figure 8
The schematic diagram of numerical simulation.The Mohr–Coulomb criterion is used as the elastic-plastic criterion in the discrete element calculation. The left and right boundaries of the numerical model are subjected to deformation constraints to limit the horizontal displacement. The bottom boundary of the numerical model is subjected to fixed constraints to limit the horizontal and vertical displacement, and the top is the free boundary. The numerical model just considers the gravity stress and neglects the tectonic stress and mining stress. The gravitational accelerationg is 9.8 m/s2. Physical and mechanical parameters of rock mass and joint are obtained by field investigation and relevant equations, which are shown in Table 2. The range of goaf roof caving corresponding to the excavation step distance can be obtained according to the flowchart of numerical simulation (as shown in Figure 5).Table 2
Physical and mechanical parameters of the Shirengou Iron Mine.
ParameterValueMethodGeological strength index GSI76Field investigationDisturbance factorD0.7Rock material constantmi24Gravitational accelerationg (m/s2)9.8Initial goaf heighthb (m)15Hydraulic conductivityK (m/s)1.68 × 10−8Kinematic viscosity of water (20°C)μk (m2/s)1.01 × 10−6Rock densityρ (kg/m3)3.58Rock uniaxial compressive strengthσc (MPa)99.44Rock elastic modulusE (GPa)80.3Poisson's ratioν0.21Joint normal fracture stiﬀnesskn (GPa)9.67Field investigation and simulation experienceJoint tangential fracture stiffnessks (GPa)6.53Joint cohesionc2 (MPa)1.67Joint internal friction angleφ2 (°)22.35Joint tensile strengthσt2 (MPa)0.23Fracture apertureB (m)2.73 × 10−4Equation (3)Rock mass material constantm6.42Equation (7)Rock mass material constants0.031Equation (8)Rock mass material constanta0.5Equation (9)Rock mass tensile strengthσmt (MPa)0.48Equation (10)Rock mass cohesionc1 (MPa)4.97Using equation (11) to fit the principal stress dataRock mass internal friction angleφ1 (°)40.98Fracture spacingS (m)1.0–4.0Field investigationGoaf depthH (m)100–500Goaf spanl (m)0–100.0Rock mass elastic modulusEm (GPa)Equation (20): the value of Em varies with the parameter SRock mass volume modulusKm (GPa)Equation (14): the value of Km varies with the parameter EmRock mass shear modulusGm (GPa)Equation (15): the value of Gm varies with the parameter Em


## 4.3. Results

When the goaf span is excavated to 100 m, the caving ranges of the goaf roof under different fracture spacing are shown in Figures9–12. The results show that the caving height of goaf roof increases with the increase of the goaf depth H under the conditions of the same fracture spacing S and goaf span l. When the goaf depth H and goaf span l are same, the caving height of goaf roof decreases with the increase of the fracture spacing S. Using the goaf roof caving height corresponding to the excavation step distance to calculate the dimensionless group in equation (5), the dimensionless chart for the range of goaf roof caving (as shown in Figure 13) can be obtained.Figure 9
Caving range of the goaf roof (whenl = 100 m, S = 1 m). (a) H = 100 m. (b) H = 200 m. (c) H = 300 m. (d) H = 400 m. (e) H = 500 m.
(a)(b)(c)(d)(e)Figure 10
Caving range of the goaf roof (whenl = 100 m, S = 2 m. (a) H = 100 m. (b) H = 200 m. (c) H = 300 m. (d) H = 400 m. (e) H = 500 m.
(a)(b)(c)(d)(e)Figure 11
Caving range of the goaf roof (whenl = 100 m, S = 3 m. (a) H = 100 m. (b) H = 200 m. (c) H = 300 m. (d) H = 400 m. (e) H = 500 m.
(a)(b)(c)(d)(e)Figure 12
Caving range of the goaf roof (whenl = 100 m, S = 4 m. (a) H = 100 m. (b) H = 200 m. (c) H = 300 m. (d) H = 400 m. (e) H = 500 m.
(a)(b)(c)(d)(e)Figure 13
Dimensionless charts for the caving range of goaf roof. (a)S = 1 m. (b) S = 2 m. (c) S = 3 m. (d) S = 4 m.
(a)(b)(c)(d)As shown in Figure13, the overlying strata of the goaf can remain relatively stable when the goaf span is small. However, the duration of the relatively stable state is closely related to the fracture spacing S. As shown in Figures 13(a)–13(d), the points O1∼O4 are the starting points of the goaf roof caving when the parameter H is 500 m, respectively. Obviously, the degree of difficulty of the goaf roof caving in the initial stage increases with the increase of fracture spacing S, which indicates the duration of the relatively stable state increases with the increase of fracture spacing S.With the continuous increase of the goaf span, the exposed area of the goaf expands gradually. It results in the occurrence of the overlying strata instability under the action of the gravity stress and internal concentrated stress. As shown in Figures13(a)–13(d), taking the dimensionless curve corresponding to the parameter H is 500 m as an example, the caving height of the goaf roof h increases with the increase of the goaf span l. The results show that there is a rapid increase process between the points Oj and Rj (j = 1, 2, 3, 4). The reason for this phenomenon is related to the accumulation of internal strain energy and the continuous development and expansion of fractures during the period when the overlying strata maintain a relatively stable state. The speed of caving can appear short acceleration when the caving of overlying strata is occurred.Under the premise that the goaf spanl is no longer expanded, a relatively stable natural balance arch can be formed when the roof caving reaches a certain height, and the roof caving could not occur for a long time. However, when the span of the goaf is increasing, the stress balance state of the roof rock is destroyed continuously, and the caving could continue to occur. As shown in Figures 13(a)–13(d), the overlying strata still caving with increase of the goaf span when the rapid growth period of the roof caving is over. At this time, the growth rate of the caving height is gradually smaller, but the scale of the caving is relatively large.


## 4.4. Application

The dimensionless charts obtained by the excavation step distance 2 m contain the goaf roof caving under lots of working conditions. Using the linear interpolation method, the range of goaf roof caving under other working conditions can be predicted. Taking the M2 main goaf of −210 m level in the Shirengou Iron Mine as an example, the goaf depthH is 355 m, the goaf span l is 72 m, and the fracture spacing S is 1.4 m. The dimensionless charts corresponding to S = 1 m and S = 2 m in Figure 13 are used to predict the range of goaf roof caving under this working condition.Defineπ1 = l/hb, π2 = h/hb, and π3 = log10(Em/ρgHν). The feature point of the goaf roof caving range predicted by the dimensionless chart is denoted as P(π1, π2, π3). When the goaf depth H is 355 m and the goaf span l is 72 m, as shown in Figure 14, the coordinates of feature points predicted by the corresponding dimensionless charts corresponding to S = 1 m and S = 2 m are P1 (4.8, 3.62, 3.52) and P2 (4.8, 3.08, 3.77), respectively. Using the dimensionless group h/hb, the ranges of goaf roof caving predicted by the two dimensionless charts are 54.3 m and 46.2 m, respectively. It can be known from the linear interpolation that the predicted range of the goaf roof caving is 51.06 m when the parameter S = 1.4 m.Figure 14
Estimation for the caving range of goaf roof. (a)S = 1 m. (b) S = 2 m.
(a)(b)As shown in Figure7, the caving height of the M2 main goaf roof in the −210 m level is 46.66 m, which is less than the range of roof caving predicted by the dimensionless charts in this work. In fact, the M2 main goaf roof is always in an unstable state. Occasionally, the M2 main goaf has the sound of caving rubble according to the field investigation. This phenomenon indicates that the M2 main goaf is in a state of caving activity, and the range of roof caving predicted by the dimensionless charts in this work is reasonable. Therefore, the dimensionless charts obtained by the dimensional analysis coupled with numerical simulation can preliminarily quantitatively evaluate the caving range of goaf roof under different working conditions, which has certain guiding significance to the treatment of the goaf.


## 5. Summary and Conclusions

The physical phenomenon of the goaf roof caving involves many physical quantities. Only considering the relationship between a single or a few physical quantities is not conducive to the study of the goaf roof caving. The method of the dimensional analysis coupled with numerical simulation in this work can be a useful technique to arrange the physical quantities and simplify the study of the goaf roof caving.The dimensionless charts in this work are convenient for the direct acquisition of the range of goaf roof caving. Using the linear interpolation method, the range of goaf roof caving under different working conditions can be predicted, which can be used as a reference for the potential hazard evaluation and rational governance of the goaf. In addition, with the continuous increase of the goaf spanl, the dimensionless curves of the caving range in the dimensionless charts experience zero growth, rapid growth and steady growth. The growth degree varies with the fracture spacing S.In conclusion, the dimensionless charts in this work can comprehensively consider the influence factors of the goaf roof caving and can quantitatively evaluate the range of goaf roof caving under different working conditions, which have a certain referential value for similar engineering problems. However, the limitation of this study is that it is temporarily applicable to the excavation of layered rock masses or ore bodies. The application of the excavation of other strata with different properties will be studied in the future work. More field applications contribute to the modification and supplementation of the dimensionless charts.



---

*Source: 7382550-2021-08-02.xml*