Dimensionless Charts for Predicting the Range of Goaf Roof Caving

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 span l, 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.


Introduction
e 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]. e 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]. erefore, 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][10][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][13][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][17][18] used the physical experiment to conduct the simulation of the goaf roof caving.
e 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]. erefore, the method of physical experiment was time consuming and costly. In addition, with the upgrading of numerical simulation software, some scholars [20][21][22] utilized the numerical simulation to study the range of goaf roof caving. e numerical simulation software mainly used the discrete element software, such as UDEC [23], 3DEC [24], and YADE [25]. e 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. e 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. erefore, 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. e intuitive information can be obtained in the charts. For example, the support classification chart of the Q 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. e 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. e research strategy on the dimensionless charts for the range of goaf roof caving is shown in Figure 1. 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. en, the discrete element software UDEC is applied to throw light on the fuzzy function relationship by numerical calculation. e 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.

Buckingham Pi Theorem
e 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]. e Buckingham pi theorem can be used to research the physical phenomena containing multiple physical quantities [36,37]. erefore, 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] f x 1 , x 2 , . . . , x n � 0. (1) If the independent physical quantities in the functional relationship are m (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. us, 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 π. erefore, equation (1) can be expressed as 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. e 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. e microscopic performance is the development of rock mass fractures including new fracture development and the extension of original fracture. erefore, 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.
en, the equivalent model for fractured rock mass can be established, as shown in Figure 2. 2 Mathematical Problems in Engineering e equivalent hydraulic conductivity K of a set of parallel fractures in Figure 2 can be expressed as [39,40] where B is the fracture aperture and S is the fracture spacing. e 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. e 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 h b , and the goaf depth H. e parameter h b 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. e rock mass elastic modulus E m and Poisson's ratio ] are also closely related to the excavation and subsequent caving of the goaf.
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 parameters h, l, h b , H, ρ, g, E m , ], S, and B. e dimension of each parameter is shown in Table 1. e surrounding rock of the goaf in Table 1 can be any kind of rock or ore. e 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 e solutions to equation (5) is derived in Appendix and summarized as follows: where the dimensionless group h/hb represents the extent of the goaf roof caving. e dimensionless group l/hb represents the ratio of width to height in the initial goaf. e dimensionless group ρgH]/Em represents the basic properties of rock mass at a certain goaf depth. e dimensionless group B/S represents the fracture development degree of rock mass.

Numerical Simulation
e 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. erefore, the discrete element software is selected to study the goaf roof caving. e 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][42][43][44][45] have achieved abundant results in the excavation of fractured rock mass. erefore, the discrete element software UDEC is selected in this work to study the range of goaf roof caving.  e Mohr-Coulomb plastic model is selected as the constitutive model in the UDEC. e mechanical parameters of rock mass in the discrete element calculation include the rock mass uniaxial tensile strength σ mt , the rock mass cohesion c 1 , the angle of internal friction of rock mass φ 1 , the rock mass volume modulus K m , and the rock mass shear modulus G m , which can be obtained according to the rock mass strength criterion and continuity assumption.

Rock Mass Uniaxial Tensile Strength σ mt .
e rock mass uniaxial tensile strength σ mt can be obtained by the Hoek-Brown criterion. e generalized Hoek-Brown rock mass strength criterion can be expressed as [46] where σ 1 and σ 3 are the maximum and minimum principal stresses, respectively. e parameter σ c is the rock uniaxial compressive strength. e parameters m, s, and a i are constants related to rock mass materials, which can be expressed as 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]. e parameters GSI and D can be obtained through field investigation. e parameter m i can be selected according to the Hoek-Brown constant charts of rock mass [47]. From equation (9), the value of parameter a i is approximately equal to 0.5 when the parameter GSI ≥ 50. en, substituting σ 1 � 0 into equation (6), the rock mass uniaxial tensile strength σ mt can be expressed as

Rock Mass Cohesion c 1 and Angle of Internal Friction of
Rock Mass φ 1 . e parameters c 1 and φ 1 can be obtained by the Mohr-Coulomb criterion [48], which states that the relationship between σ 1 and σ 3 is linear.
en, 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 According to equation (13), the rock mass cohesion c1 can be obtained.

Elastic Modulus E m , Volume Modulus K m , and Shear Modulus G m .
e rock mass volume modulus K m and the rock mass shear modulus G m can be expressed as e value of Poisson's ratio ] is easy to obtain. When the elastic modulus of rock mass E m is obtained, the values of the parameters K m and G m can be calculated by equations (14) and (15). erefore, the parameter E m is the key parameter to be determined. e 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 spacing S and fracture aperture B have the most direct effect on the parameter E m . 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 where Δε s is the strain of fractured rock mass and E is the rock elastic modulus. en, the displacement Δu s in the fracture spacing S can be expressed as Define the deformation of a single fracture is Δu b , then [39] where is k n the normal fracture stiffness. As shown in Figure 4, define the deformation of the continuous rock mass is Δu m ; then, Substituting equations (17) and (18) into equation (19) gives According to equations (6)- (20), the parameters σ mt , c 1 , φ 1 , K m , and G m can be calculated. e rock physical and mechanical parameters involved in the calculation can be obtained from the field investigation. Figure 5, 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. e elastic modulus of rock mass E m 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. en, 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.

Numerical Calculation Flowchart. As shown in
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 Figure 5 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.

Background.
e ore deposit in the Shirengou Iron Mine is an Anshan-type magnetite deposit with a surface elevation of +145 m. e ore body is produced in layers with an inclination of 65°-75°. e surrounding rock is hornblende gneiss. Above the − 60 m level, the shallow-hole mining method was adopted for the stable ore rock. e shallowhole mining method had left many goafs of varying sizes, forming an intricate goaf group in space. e 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 Figure 7, 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. us, a large irregular goaf (referred to as M2 main goaf in Figure 7) is formed.
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.

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]. e numerical model for the goaf roof caving is established (as shown in Figure 8). e 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. erefore, 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. e excavation step distance of the goaf is set to 2 m, which can obtain massive data of the goaf roof caving. Mathematical Problems in Engineering 5 e Mohr-Coulomb criterion is used as the elasticplastic criterion in the discrete element calculation. e left and right boundaries of the numerical model are subjected to deformation constraints to limit the horizontal displacement. e 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. e numerical model just considers the gravity stress and neglects the tectonic stress and mining stress. e gravitational acceleration g is 9.8 m/s 2 . Physical and mechanical parameters of rock mass and joint are obtained by field investigation and relevant equations, which are shown in Table 2. e 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).

Results.
When the goaf span is excavated to 100 m, the caving ranges of the goaf roof under different fracture spacing are shown in Figures 9-12. e 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.
As shown in Figure 13, 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 Figures 13(a) example, the caving height of the goaf roof h increases with the increase of the goaf span l. e results show that there is a rapid increase process between the points O j and R j (j � 1, 2, 3, 4). e 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. e speed of caving can appear short acceleration when the caving of overlying strata is occurred.
Under the premise that the goaf span l 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.

Application.
e 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 depth H is 355 m, the goaf span l is 72 m, and the fracture spacing S is 1.4 m. e 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]). e 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.
As shown in Figure 7, 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    (11)    unstable state. Occasionally, the M2 main goaf has the sound of caving rubble according to the field investigation. is 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. erefore, 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.

Summary and Conclusions
e 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. e 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. e 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 span l, the dimensionless curves of the caving range in the dimensionless charts experience zero growth, rapid growth and steady growth. e 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. e 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.