Calculation of Elastic Modulus for Fractured Rock Mass Using Dimensional Analysis Coupled with Numerical Simulation

Underground mining activities make the fractures in the natural rock mass develop randomly. The elastic modulus of the fractured rock mass E m is changed with the redistribution and development of the fractures. An equivalent model of fractured rock mass is structured to represent the hydraulic conductivity and the rock mass strain because of the continuum theory. Dimensional analysis is very useful to build relationship of the parameters in complex physical phenomena. Based on the engineering phenomenon of groundwater ﬂowing into the goaf along the fracture in the rock mass, a fuzzy expression among parameters such as the parameter E m , the seepage ﬂow Q , and the exposed area of the goaf A is obtained using dimensionless analysis. To calculate the parameter E m , the fuzzy relationship is then characterized by Darcy’s law and numerical simulation. Under the scripting environment of Python, an automated program to realize the numerical simulation of all scenarios is established, which also provides convenience for drawing the dimensionless ﬂux charts. The results show that the parameter E m can be calculated by the dimensional analysis coupled with numerical simulation. In addition, the parameter E m decreases with the increase of the parameter Q , and the integrity of rock mass is also worse. Finally, a mine example is used to verify the feasibility of the method.


Introduction
e elastic modulus of rock mass is one of the most important mechanical parameters in geotechnical engineering and has important guiding significance for engineering practice [1]. Considering the cutting effect of the structural plane in the rock mass, numerous fractures are randomly formed inside the rock mass. erefore, the study of the elastic modulus of fractured rock mass is more suitable for engineering practice than the complete rock mass. Most scholars carried out many works to obtain more practical elastic modulus, such as empirical estimation method [2][3][4][5][6][7][8], acoustic wave test method [9,10], and in situ test method [11]. However, the empirical estimation method is subject to subjective factors, which leads to the estimation results vary with different engineers. e acoustic wave test method and in situ test method are costly and are obviously interfered with by environmental factors. However, it is worth noting that the above methods indicate that the fractures in the rock mass have a great influence on the elastic modulus, so we attempt to study the fractures of the rock mass to calculate the elastic modulus. e fractures in natural rock mass are chaotic and disorderly [12]. Most scholars [13][14][15][16] research the mechanical behavior of fractured rock mass by establishing an equivalent continuous medium model, that is, to convert discrete fractured rock mass into continuous medium by equivalent hydraulic conductivity methods. e equivalent continuum model simplifies the study of the mechanical behavior of fractured rock mass and provides convenience for studying the elastic modulus of fractured rock. On the other hand, the natural rock mass is generally saturated with groundwater, which produces seepage within the rock mass [17,18]. e seepage characteristics are determined by the fractures. Considering that the seepage flow in the exposed area of the goaf can be measured, we can establish the connection between the elastic modulus and the seepage flow in the goaf by utilizing the fractures of rock mass.
In engineering, the elastic modulus and fractures of the rock mass change due to the influence of mining, and this change involves many physical quantities, not only related to the physical and mechanical parameters of rock mass and the mining parameters but also to the initial state of the fractures. To find the relationship between the elastic modulus and the seepage flow, we attempt to utilize dimensional analysis coupled with numerical simulation, say, finite element method (FEM). Dimensional analysis [19][20][21] is widely used to investigate hypotheses about complex physical phenomena. e Buckingham Pi theorem [22][23][24][25] in dimensional analysis provides a way to form a function of dimensionless variables. As a result, the stress and strain at any circumstances can be easily estimated using FEM following the pattern given in the dimensionless function.
In the following sections, based on the engineering phenomenon of groundwater flow into the goaf, we regard the fractured rock mass as an equivalent continuous medium model and derive the expressions of hydraulic conductivity and rock mass strain under three-dimensional conditions. en, we integrate the physical quantities of this physical phenomenon using the Buckingham Pi theorem and establish the functional relation between the elastic modulus and the seepage flow. Considering that the functional relationship contains the rock mass strain, we use the numerical simulation to solve the strain value in the rock mass and the indirect coupling method to calculate the hydraulic conductivity and the corresponding seepage flow. With the help of the scripting language Python, the FEM input files for all scenarios formed automatically can be executed. Consequently, we can draw the dimensionless flux charts about the elastic modulus conveniently; that is, we can obtain the elastic modulus of fractured rock mass under any working condition. Finally, we take a mine example to verify the accuracy of this method.

Equivalent Model for Fractured Rock Mass
e equivalent model of fractured rock mass (as shown in Figure 1) approximately treats the rock mass object as a continuous medium and assumes that the normal direction of the fractured rock group corresponds to the direction of coordinate axis [26]. Generally, the normal directions of the surfaces in the above three groups of fractures are arbitrary and not perpendicular to each other because of the anisotropy of rock mass. However, Warren and Root [27] proposed an equivalent model for fractured rock mass that transformed the complex fracture system to a threedimensional state with three groups of mutually orthogonal fracture systems, which greatly simplifies the study of the hydraulic conductivity. In this way, we can analyze the mechanical behavior of fractured rock mass by using the continuum mechanics method and the seepage theory from a macro perspective, regardless of the geometry structure of a particular fracture.
Regardless of the deformation of the material, the hydraulic conductivity of a set of parallel fractures can be defined as [28,29] where g is the gravitational acceleration and μ k is the kinematic viscosity of water. e parameters B and S are the parallel fractures of aperture and spacing, respectively. Equation (1) expresses the hydraulic conductivity in the absence of deformation and for a prescribed initial fracture aperture. Under the condition of orthotropic, the hydraulic conductivity of three groups of basic fractures that the normal direction is consistent with the direction of coordinate vector can be expressed as where K ij is the tensor for hydraulic conductivity and the numbers 1, 2, and 3 refer to the coordinates x, y, and z in Figure 1, respectively. K 1 ij can be expanded as Obviously, under the condition of orthotropic, when i ≠ j and K 1 ij � 0, and when i � j � 1 and K 1 ij � 0, then And equation (2) can be expressed as  Mathematical Problems in Engineering Assuming that the aperture of fractures in the three groups of fractures is expressed by B 1 , B 2 , and B 3 , respectively, and the spacing of fractures in the three groups of fractures is expressed by S 1 , S 2 , and S 3 , respectively. en, According to equation (6), the hydraulic conductivity is positively proportional to the cube of the aperture of fractures, indicating that the aperture of fractures plays a leading role in the change of the hydraulic conductivity. Regardless of the influence of the change of the spacing of fractures, we label the variables of the aperture of fractures as Δu b1 , Δu b2 , and Δu b3 . Assuming that the initial aperture and spacing of the three groups of fractures in the fractured rock mass are the same (namely, , the fractured rock mass equivalent hydraulic conductivity can be expressed as As shown in Figure 2, we select the model in Figure 1 to project onto the x-z plane. Taking a set of fractures in the rock mass as the analysis object, we regard the fractured rock mass as an equivalent continuous rock mass.
From Hooke's law, for a linear rock matrix, where Δε s is the strain in the direction of the applied stress Δσ and E is the elastic modulus of the rock matrix.
erefore, between two fractures in the rock mass at separation S, the displacement Δu s is as follows: Defining the deformation across a single fracture as Δu b , then [26] where k n is the normal fracture stiffness. Total displacement Δu m is obtained from the sum of component displacements as where E m is the elastic modulus of the fractured rock mass. Substituting equations (9) and (10) into equation (11) gives Rearranging equation (11) yields where R m is the elastic modulus reduction ratio [26]. Rearranging equations (10) and (9), Substituting equation (14) into equation (11) gives where Δε is the strain for equivalent continuous rock mass and Δε � (Δσ/E m ). Substituting equation (13) into equation (15) and rearranging the equation yield Generalizing equation (16) to the three-dimensional situation, Substituting equation (17) into equation (7) gives

Mathematical Problems in Engineering
Equation (18) expresses the functional relationship between the hydraulic conductivity and the strain. Meanwhile, the relationship shows that the strain has a great influence on the hydraulic conductivity. When the rock mass is excavated, the groundwater will flow into the excavation area. Darcy's law [30] suggested that the relationship between the seepage flow and the hydraulic conductivity in the anisotropic medium can be expressed as where Q i is the seepage flow. J is the hydraulic gradient.
As to the orthogonally anisotropic medium, expanding the equation (19) can obtain erefore, as to the mined excavations, the total seepage flow of the goaf can be expressed as

Buckingham Pi Theorem
e Buckingham Pi theorem can describe physical processes involving n physical quantities (x 1 , x 2 ,..., x n ), and the physical processes can be expressed as functional relations of the physical quantities involved as follows: Among the n physical quantities involved in the physical process, the n − m nonbasic physical quantities can be represented by the composite quantity of m basic physical quantities and converted into n − m dimensionless parameters, then equation (22) can be expressed as [22,31] e important fact to notice is that the new relation involves fewer variables than the original relation, which simplifies the theoretical analysis and experimental design alike. e development of fractures in rock mass determines the initial hydraulic conductivity and the elastic modulus of the fractured rock mass. Considering that the stress state of the surrounding rock disturbed by the mining activity, a certain displacement can be generated in the normal direction of the fracture surface, which caused the hydraulic conductivity to change. In addition, a large number of goafs of different sizes are formed by the excavation of the mine, and the exposed area of the goaf is directly related to the seepage flow, indicating that there is a fuzzy relationship between seepage flow and elastic modulus of the fractured rock mass. However, the fuzzy relationship involves many parameters (as shown in Figure 3), such as the exposed area of goaf, the depth of goaf, the bulk weight of rock mass, and the aperture and spacing of fractures.
Considering that the fuzzy relationship between the elastic modulus of the fractured rock mass and the seepage flow involves many parameters, and so far there is no clear equation to express the relationship, we utilize the Buckingham Pi theorem to derive the relationship as follows: where Q is the seepage flow, H is the depth of the goaf, A is the exposed area of the goaf, h is the height of the goaf, c is the bulk weight of rock mass, and g is the gravitational acceleration. Table 1 summarizes all physical quantities and corresponding dimensions in the fuzzy relationship. e physical quantities Q, h, and E m are selected as the basic physical quantities, and the remaining physical quantities are regarded as the quantities to be determined. en, where α i , β i , c i , and η i are the parameters to be determined, and i � 1, 2, 3, 4, 5, and 6.
Taking equation (25) as an example, we can obtain From equation (26), we can obtain α i � c i � 0 and η i � −β i ; then, Similarly, π 2 , π 3 , π 4 , π 5 , and π 6 can be expressed as erefore, the dimensionless relationship can be expressed as e values of η 5 , β 1 , β 2 , β 3 , c 4 , and β 6 in equation (29) can change the axis values of the dimensionless flux charts but cannot change the final evaluation results. To make the calculation more convenient, the values of the parameters need to be clear. Considering that the parameter Q is one of the key parameters for calculating the elastic modulus for fractured rock mass, we can take η 5 � 1/2. Referring to the expression of the hydraulic conductivity equation (1) and the dimension of B and S, the expression B/S can be a dimensionless expression. Besides, the magnitude of the parameters h and B is quite different; therefore, we can take β 1 � 1 and β 2 � −1 to obtain the dimensionless expression B/ S. Generally, the value of parameter A is greater than the parameter h 2 ; we take β 3 � −1 to make the dimensionless flux chart clearer. Given that the parameter h can be eliminated by using the dimensionless expression E m /c h to multiply the dimensionless expression h/H, we take c 4 � 1 and β 6 � 1. Finally, we can obtain Considering that the multiplication and division of dimensionless quantities have no effect on the final analysis results, equation (30) can be expressed as Poisson's ratio ] can reflect the lateral deformation of the material. As one of the important factors in the calculation of elastic mechanics, it has a certain effect on the seepage flow and elastic modulus. Considering that its dimension is 1, equation (31) can be expressed as or

Numerical Simulation
Zhengxing Multiple goafs of varying sizes are formed due to the unreasonable mining planning. e maximum exposed area of the goaf is 10000 m 2 , and the minimum exposed area is 2400 m 2 . e maximum height of the goaf is 50 m, and the minimum is 10 m. ere are many developed local fault zone and derived fractures in the mine because of the unreasonable support. Considering that the in situ stress gradually increases with the mining depth, the confining pressure of the goaf is more and more complicated. In addition, the increase of mining depth accelerates the development of fractures in rock mass.
To realize the automatic execution for numerical simulation of all scenarios, we use Python as the scripting environment for modeling and finite element calculation. Gmsh [32][33][34], a 3D modeling software, can conduct finite element mesh division. In addition, the output of this software can be .tmpl file format for Python software to call, so that we can use Python's loop function to change the model parameters in the .tmpl file, so as to achieve automatic model establishment and mesh division. SfePy [35,36] (Simple Finite Element in Python) is an open-source finite element analysis software written in Python, which uses the Galerkin method to solve. At present, the SfePy has mature models in various aspects such as acoustics, fluid mechanics, and elastic mechanics. Use Python to call the SfePy module and load the model created by Gmsh into SfePy for setting boundary condition and finite element calculation. Finally, the dimensionless diagram of elastic modulus can be obtained by using the scientific calculation module of Python to substitute the finite element calculation results into the function relations determined by the Buckingham Pi theorem. e flowchart of numerical simulation is shown in Figure 4.
Generally, numerical simulation simplifies the actual size of irregular goaf for finite element calculation. e modeling and meshing based on Gmsh are shown in Figure 5. e size of the rock mass is 1000 m * 500 m * 1100 m, and the size of the goaf is an unfixed value. We use the parameters of l, w, and h in Figure 3 to control the size of the goaf. e density of the mesh is increased near the goaf.
When the model is loaded into the SfePy, considering that the finite element calculation in this work involves the coupling calculation of elasticity and fluid mechanics, we set two boundary conditions: (1) stress boundary condition: the vertical stress is calculated according to the gravity stress, and the horizontal stress is calculated according to the Poisson effect. e lateral side of the rock mass restricts the horizontal displacement, and the bottom of the model is a fixed constraint. e Earth surface and the goaf are free. (2) Stress boundary condition: the model is a saturated rock mass medium, and a fixed water head value is applied to the boundary around the rock mass. e water head value is the distance from the underground water surface to the goaf. e Earth surface is a free boundary and the pore water pressure is zero. e physical and mechanical parameters of Zhengxing Iron Mine are shown in Table 2. e parameters of c, g, υ, and μ k are fixed values, which are taken from field investigation and laboratory experiment. en, as shown in equation (39) Table 2).     Table 2. e parameter A can be expressed as Once the parameters w and h are determined, the value of A is determined by the span of the goaf l. erefore, the value range of l is 3.0∼160.0 m. e final value of all physical parameters is shown in Table 2.
When the seepage flow into the goaf is calculated, the strain induced by the mining process is first determined by the finite element model. en, the seepage flow Q can be obtained by equation (21). e results of numerical simulation are then substituted into equation (33) to obtain the dimensionless parameters. As the size of the goaf changes, we must process a large amount of data. Consequently, we establish an automated program to handle the onerous task.
is procedure is described in the flowchart of Figure 6, where the parameters B, H, and l are variables, and the parameter B controls the development of fractures in different rock masses, the parameter H controls the gravity stress, and the parameter l controls the exposed area of the goaf. e calculation of the seepage flow is combined with theoretical analysis and numerical simulation. e mathematical model of saturated continuous medium is established by Gmsh, and the finite element numerical calculation results are obtained by SfePy. Using the value of element strain in the numerical simulation results, the hydraulic conductivity can be obtained by equation (18). Combining the hydraulic conductivity and the value of pore water pressure in the numerical simulation results, the seepage flow can be calculated by equation (21). e hydraulic gradient J can be obtained by Darcy's law [30]: where ΔH w is the water head drop of nodal element sections and L is the spacing of nodal element sections. e water head H w can be expressed as where p is the pore water pressure, c w is the bulk density of water, and h w is the position head.

Results
Based on the dimensional analysis coupled with numerical simulation, the dimensionless flux charts of the fractured rock mass parameters obtained using the flow chart shown in Figure 6 are shown in Figure 7. e results show that the fuzzy relationship between the parameter Q and the parameter E m can be expressed indirectly. In addition, as shown in Figure 7, the parameter Q increases with the increase in the ratio of B/S at the same height of the goaf h. In other words, the worse the integrity of rock mass, the greater the parameter Q, which is consistent with the engineering practice.

Mathematical Problems in Engineering
Certainly, the relationship among the parameters of Q, A, and E m can be obtained under a certain working condition. Assuming that the height of the goaf h is a constant, we can obtain the relationship between the parameter Q and the parameter A when the value of dimensionless expression (E m /]cH) is given. As shown in Figure 8, when h � 20 m, the parameter Q increases with the exposed area of the goaf A, which is in line with the actual situation. Similarly, the relationship between the parameter Q and the parameter E m is shown in Figure 9. In the same working condition, the parameter Q decreases with the increase of the parameter E m . In fact, the seepage flow Q is measured easily in the field. erefore, we can use the measured parameter Q to estimate the parameter E m under specific working conditions. Certainly, the relationship among the parameters of Q, A, and E m can be obtained under a certain working condition. Assuming that the height of the goaf h is a constant, we can obtain the relationship between the parameter Q and the parameter A when the value of dimensionless expression (E m /]cH) is given. As shown in Figure 8, when h � 20 m, the parameter Q increases with the exposed area of the goaf A, which is in line with the actual situation. Similarly, the relationship between the parameter Q and the parameter E m is shown in Figure 9. In the same working condition, the parameter Q decreases with the increase of the parameter E m . In fact, the seepage flow Q is measured easily in the field. erefore, we can use the measured parameter Q to estimate the parameter E m under specific working conditions. Figure 7 contains the relationship among the physical quantities related to the fractured rock mass parameters of Zhengxing Iron Mine under all possible scenarios. en, we   Table 1 Selecting a goaf at the mining level as the research object, we can measure that the exposed area of the goaf A is 5675 m 2 , the height of the goaf h is 20 m, and the seepage flow of the goaf Q is 16 According to B/S � 1.27 × 10 −4 , we can utilize Figure 7(c) to calculate the elastic modulus of fractured rock mass E m . As shown in Figure 10, the characteristic point G (14.19, 3.46 × 10 −8 ) represents the coordinate point of the exposed area and seepage flow of the goaf. e interval [4.09, 4.14] of the right coordinate axis is divided into 5 equal parts to accurately obtain the corresponding elastic modulus of fractured rock mass E m at G point. Making a corresponding curve through the point G, as shown in the black dashed line in Figure 10, we can obtain that the value of log 10 (E m /]cH)is 4.14. e trend of the curve through the point G is similar to that of the red curve in Figure 10. Considering that the elevation of Earth surface is +50 m, we can obtain the depth of the goaf H is 250 m. Finally, we calculate that the parameter E m is equal to 23.29 GPa.
Most scholars carried out many works to obtain more practical elastic modulus. e empirical estimation method   [3] proposed the expression of the relationship between the ratio E m /E and the rock mass quality index RQD, as follows: Bieniawski [4] suggested that the elastic modulus of rock mass is correlated with the rock mass quality index RMR. When RMR > 50, E m can be expressed as When RMR ≤ 50, Serafim and Pereira [5] suggested that E m can be expressed as Based on the database of rock mass deformation modulus measurements, Hoek and Diederichs [6] proposed that the elastic modulus of rock mass E m can be represented as where GSI is the geological strength index. D is a factor that depends upon the degree of disturbance to which the rock mass has been subjected to blast damage and stress relaxation. e index parameters of rock mass in the empirical estimation method introduced in this work are shown in Table 3. According to the index parameters of the Zhengxing Iron Mine, the rock mass elastic modulus E m corresponding to each empirical formula can be calculated. As shown in Table 3, the value of E m calculated by the empirical formula is close to the calculation result of the method in this work, which indicates that it is feasible to use the dimensional analysis coupled with numerical simulation to calculate the elastic modulus of the fractured rock mass. In addition, when the working conditions are clear, we can directly obtain the value of E m from the dimensionless flux charts according to the seepage flow in the goaf, which facilitates the calculation of the elastic modulus on-site.

Conclusions
e aperture and spacing of fractures, which are closely related to the hydraulic conductivity, have the most direct influence on the elastic modulus of the fractured rock mass E m . e proposed equivalent model is useful to analyze the mechanical behavior of fractured rock mass and conceptually represents the relationship between the hydraulic conductivity and the strain. Considering that the joint stiffness parameters are elusive, the hydraulic conductivity can be obtained conveniently by introducing a modulus reduction ratio R m . e dimensionless flux charts developed in this work integrate most physical quantities involved in the physical process of groundwater flow into goaf. With the help of dimensional analysis coupled with numerical simulation, the elastic modulus of the fractured rock mass E m can be calculated. When the working conditions are specific, the parameter E m can be obtained by querying the relationship diagram of the seepage flow Q and the parameter E m . Taking the Zhengxing Iron Mine as an example, the elastic modulus for the fractured rock mass calculated by the method in this paper is close to the calculation result of the existing empirical formula, which shows that it is feasible to use the method to calculate the parameter E m . In addition, the Python introduced in this work, which is a useful tool to handle the data generated under different working conditions, provides an automated program to realize the automatic execution for numerical simulation of all scenarios.
In conclusion, the method of the dimensional analysis coupled with the numerical simulation throws new light on the reasonable calculation of the parameter E m . Further field measurements and observations contribute to the straightforward practical use of the method.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that they have no conflicts of interest.