Abstract
Aiming at the deficiency of the conventional multiphysics coupling model, the deterioration of strength parameters was considered by defining elastoplastic damage variables, and the heterogeneity of strength parameters was expressed by the Weibull distribution function. In addition, the relation between effective stress and the anisotropic permeability matrix was established, and the blast was transformed into a load boundary condition. On this basis, an improved multiphysics coupling model that considered damage and disturbance was constructed, while a corresponding finite element calculation program was developed. Taking an excavation stope as the object, the characteristics of the mining-induced stress, seepage, and failure were analyzed by an improved multiphysics coupling model and compared with actual detection data. The results show that the improved model reflects the extent and range of mining-induced failure more accurately and fits well with the actual detection. These results are compared to the conventional multiphysics coupling model and a single physics model. It is indicated that the improved multiphysics coupling model and corresponding calculation program are effective and rational.
1. Introduction
Mineral resources are the material basis for developing a national economy and safeguarding national security. With the rapid industrialization and urbanization process in China, the demand for mineral resources is significant, resulting in the depletion of shallow resources, and many mines have now advanced to a deep mining stage. Deep mining is accompanied by complex mining conditions. Rock mass engineering is in a coupled system composed of a stress field, seepage field, and temperature field, and rock deformation and its failure mechanism is very complex. The conventional single field (usually stress field) analysis has major limitations, and the analysis of rock mass response and failure under multifield coupling has not only become an important subject of rock mechanics research [1], but also has important practical significance.
For most rock mass engineering applications, the variation in temperature gradient is small, and the study of multifield coupling is mainly focused on the coupling of the stress field and seepage field. Baghbanan and Jing [2] simulated the coupling process of seepage and stress under a different fracture distribution and pressure coefficient using the discrete element software known as UDEC and obtained the rule of stress variation on permeability and seepage path. Figueiredo et al. [3] established the function relationship between rock mass porosity and isotropic volume strain. Then, the quantitative relationship between rock mass deformation and permeability was determined, and the numerical analysis of a fluid-solid coupling process in a deep rock mass project was realized. Based on pore elasticity theory, Zhang and Wang [4] deduced the relationship between permeability coefficient and stress variation and determined the range of mining failure in a mine. A fully coupled mathematical model of a seepage field and stress field was proposed by Wang et al. [5], and the two developments of a corresponding numerical calculation program and the application analysis of the fluid solid response characteristics of a hydropower station were realized. According to the actual construction situation of a tunnel [6], the analysis model of tunnel excavation under the action of fluid structure interaction was established by using elastoplastic theory, which provided technical reference for the design and construction of a tunnel.
At present, multifield coupling research is usually based on the Biot consolidation theory, which assumes that the deformation of rock is linearly elastic and the seepage obeys Darcy’s law. The conventional multifield coupling analysis has good practicability, which can realize reliable analysis and solve problems in engineering, but there are also several shortcomings. On the one hand, conventional multifield coupling research does not consider the effect of damage, which means that elastic modulus, cohesion, and coupling strength parameters are constant in the calculation process and that the research will not reflect the nonlinear elastic compression or expansion of microdefects caused by deformation, an important feature of strain hardening or softening stage. The rock mass will behave as an ideal elastoplastic body, which is inconsistent with a real-world scenario. Several scholars have gradually focused on the influence of damage and built the corresponding multifield coupling model [7, 8], but the assumption is that damage occurs only in the plastic stage, which ignores the existence of elastic damage. On the other hand, the engineering disturbance (especially blasting) time is very short compared to the geological action, but the effect and influence of blasting on rock engineering is very strong. In the area of multifield coupling research, the disturbance effect is seldom considered. In addition, the conventional multifield coupling study also lacks the consideration of the important characteristics such as rock mass strength heterogeneity and seepage anisotropy.
Therefore, in this study, the definition of staged damage variables and the equivalent calculation of blasting were taken into account. Considering the heterogeneity of rock mass strength parameters and seepage anisotropy, an improved multifield coupling model that considers damage and disturbances was established, and the corresponding numerical calculation program was compiled. Through numerical analysis and comparison of the excavation response and failure characteristics of a stope, the effectiveness of the improved multifield coupling model and its calculation program was verified.
2. Conventional Multifield Coupling Model
In order to simplify the expression of the formula, this study adopts the abstract notation of a tensor (black body representation). The definition of the rock stress tensor sigma, strain tensor for, based on the assumption of small deformation, equilibrium equations and geometric equations of rock mass are as follows:where is the Laplace operator, is the volume force tensor, and is the displacement tensor.
It is generally believed that the internal seepage of a rock mass is incompressible and obeys Darcy’s law:where is the seepage speed, is the source item, is the seepage pressure, is the dynamic viscosity coefficient of fluid, is the fluid density, is the gravitational acceleration, and is the permeability tensor.
The conventional multifield coupling method is mainly based on the Biot consolidation equation [9], that is, the coupling of stress and seepage through effective stress principle and volume strain source term:where is the effective stress tensor, is the elastic stiffness tensor of rock mass, is the effective stress coefficient of Biot, and is the volumetric strain of rock mass.
Compared to the single stress field analysis, the multifield coupling model considers the interaction between seepage and stress, which can reflect a real-world scenario more effectively. However, the conventional multifield coupling model does not consider the influence of damage, disturbance, rock mass heterogeneity, and other important factors, it is still not sufficient to solve the deformation and failure analysis of a rock mass under complex conditions. Therefore, it is necessary to improve the conventional multifield coupling model.
3. Improved Multifield Coupling Model
3.1. Damage Variable Definition
Damage is the existence and evolution of microdefects, and the deterioration of material stiffness and strength is caused by macroscopic damage. The conventional multifield coupling model does not consider the influence of damage, which will lead to the omission of nonlinear elastic rock mass deformation and the softening or hardening of plastic strain. The rock mass becomes an ideal elastoplastic body, which is inconsistent with real-world behavior. Therefore, the corresponding damage variables are defined for the elastic and plastic stages, respectively. Referring to the related research [10, 11], the damage evolution of the elastic stage is related to elastic strain, and the damage evolution of plastic stage is related to the equivalent plastic strain. The expression of damage variable, , is defined by the following stages:where is the elastic strain, is the mean value of elastic strain, is the equivalent plastic strain, and and are correction coefficients.
The experimental results of Yu et al. [12] show that the damage has obvious influence on the elastic modulus and cohesion, but has little influence on the internal friction angle. According to the Lemaitre strain equivalence principle, the relationship between the strength parameters of rock mass and the damage degradation is obtained:where and are initial modulus of elasticity and cohesion and and are modulus of elasticity and cohesion considering damage.
3.2. Improved Multifield Coupling Model
The damage variable, , is introduced into the conventional multifield coupling model, and the improved multifield coupling model at the elastic stage is obtained as follows:
In the improved multifield coupling model, by considering the damage elastic stiffness tensor and permeability tensor , we can reflect the effect of stress and seepage damage. In addition, in the rock constitutive equation reflects the effective stress effect of seepage, and the bulk strain source, , reflects the rock mass effect of deformation on seepage continuity. The coupling action of stress, seepage, and damage has taken into account in the improved model. At the same time, seepage usually shows anisotropic characteristics [13], the relationship between permeability and effective stress is as follows:where is the initial permeability, is the correlation coefficient (0.35 MPa−1), and is the effective principal stress.
After entering the plastic stage, the improved multifield coupling model lies in the difference of the constitutive equations of the rock mass, that is, the stress and strain are no longer subjected to one-to-one correspondence and are related to the loading-unloading path. Therefore, it is necessary to use the incremental theory to express the relationship between stress and strain. In this study, the Drucker–Prager yield model (DP model) [14, 15] which considers three principal stresses are used, and the damage variable is introduced into the DP model, while the yield function and potential function are in turn:where is the first invariant of the stress tensor, is the second invariant of the tensor of partial stress, and are internal friction angle and expansion angle of rock mass, and is the hardening function.
The plastic strain increment is calculated by the associated flow rule:where is the plastic multiplier, .
The loading-unloading criterion is described as follows:
The piecewise linear hardening function is used to approximate the nonlinear hardening function, and the expression is described as follows:where is plastic hardening or softening internal variables of a rock mass [16].
Compared to the conventional yield model, the modified yield model takes into account the influence of damage variables. In the calculation process, the strength parameters will deteriorate with the evolution of damage. Meanwhile, the piecewise linear function is used to approximate the nonlinear hardening function. Therefore, the improved model can overcome the deficiency of rock mass as an ideal elastoplastic body in conventional multifield coupling, thus reflecting the strain hardening or softening characteristics of the rock mass yielding stage.
3.3. Blasting Equivalent and Expression of Strength Heterogeneity
The effect and influence of blasting on rock engineering is very strong, but it is seldom considered in a multifield coupling study. Blasting can be converted into an action load by an equivalent calculation based on the charge configuration in a given engineering application, and it is applied to the boundary of a free surface in a numerical simulation to reflect the effect of blasting on the multifield coupling process. The method of uncoupled charge is usually used in deep hole blasting in a mine, and the equivalent load of blasting is calculated as follows [17]:where is the density of explosives, is the explosive velocity, is the isentropic exponent of explosives, is the charge diameter, is the hole diameter, is the total length of charge, is the hole length, is the pressure increase coefficient, and and are the attenuation coefficients of the blasting load.
The Weibull probability density function [18] is used to characterize the heterogeneity of rock mass strength parameters:where are the rock mass strength parameters (such as elastic modulus, cohesive force, etc.), is the mean value of the strength parameter, and is the correlation coefficient of inhomogeneity.
As shown in Figure 1, when the value of the nonhomogeneous coefficient, , is higher, the value of the intensity parameter becomes more concentrated and tends to the mean value.
3.4. Failure Criterion
For rock mass engineering, the concept of break approaching degree is introduced to quantitatively indicate the failure degree of rock mass units. The physical meaning is that the distance between the most unfavorable stress path to the yield surface and the distance between the most stable reference point and the most unfavorable stress path along the most unfavorable stress path to the yield surface in the same Rhodes angle is calculated as [19]where is the stress Rhodes angle.
Through a series of geotechnical engineering applications, Zhang et al. [20] obtained the criterion of failure approaching degree: when the calculated value of is greater than 2, the rock mass element will be destroyed.
4. Numerical Calculation Programming
The improved multifield coupling model is constructed by a set of differential equations that include multiple dependent variables, while considering the heterogeneity of rock mass strength parameters and seepage anisotropy. The model is highly nonlinear, and the coupling relation is very complex. It is difficult to obtain an analytical solution by means of series variation or integration. Therefore, the improved multifield coupling model is solved using a numerical method. COMSOL Multiphysics is a finite element numerical analysis software specifically designed to solve multifield coupling problems. It also provides powerful programming functions. Based on the COMSOL numerical platform, a program for improving the multifield coupling model is developed by using the Matlab M programming language. In Figure 2, the programming codes of some cores are listed, which mainly relate to the function relation between permeability coefficient and stress, the definition of the damage variable, the effective stress, and the source of volume strain.
The validity and rationality of the improved multifield coupling model and its calculation program will be verified by practical application.
5. Engineering Application
5.1. Engineering Situation
An underground room mining area is located between a −400 m level and −360 m level section. The stope ore reserves are approximately 22,200 t, the lower ore average angle is 55°, and the ore grade (Pb + Zn) is 14%. The regional strata are D3ta, and the main lithology is deep gray thick-bedded limestone. A total of six blasting schemes are planned for the stope. The mining method is column back lateral caving method without a bottom hole (shown in Figure 3), which is composed of the upper hole down a drilling formation, and the stope is disposed of hollow holes as the initial compensation space, layer by layer formed by ring groove blasting area, the remaining ore body is mined by lateral cave blasting. The lower part is a nonpillar mining flat-type chamber. The ore is transported by remote LHD, and the cavity is subsequently backfilled.
5.2. Simulation Parameters and Conditions
Several parameters (mainly from field research group and laboratory test [21]) are listed in Table 1, and the dynamic strength parameters of rock mass are obtained by [22] and used for blasting numerical simulation. Based on gravity stress and tectonic stress, the initial geostress field is obtained by using the multiple regression method. The expression of the heterogeneity of rock mass strength parameters is achieved by defining the Weibull probability density function (Figure 4). According to the actual charge situation and the equivalent calculation method (Formula (14)), the split blasting of the stope is equivalent to the load time history curve (Figure 5), which is applied to the excavation free surface as a boundary condition.
5.3. Simulation Results, Analysis, and Verification
(1)Stress distribution characteristics: the stope is perpendicular to the direction of the ore body (north-south strike), the east and west ends are the ore rock interfaces, and the two northern and southern ore bodies are adjacent to the ore body. The mechanical response and failure characteristics of stope will directly affect the stoping of the adjacent stope. Therefore, the typical profile of the maximum and minimum principal stress in the direction of direction is intercepted, and the stress distribution characteristics are analyzed (the tensile stress is positive, and the compressive stress is negative).
As shown in Figure 6, because of the homogeneity of the rock mass strength parameters, the stress field distribution obtained by numerical simulation is also heterogeneous. The excavation causes the redistribution of stress around the stope, which is mainly caused by the stress drop and stress concentration which is caused by unloading in several areas. The stress concentration area is mainly located at the top and bottom of the stope, and the maximum principal stress is 30.2 MPa. The unloading area is mainly located on the stope side, and the maximum principal stress is 5.81 MPa. Because of unloading and blasting, a certain range of tensile stress zones is formed, and the maximum tensile stress value is 3.03 MPa (Figure 6), which exceeds the tensile strength (2.35 MPa) of the surrounding rock. This causes tensile failure to occur.(2)Seepage distribution characteristics: a typical seepage characteristic profile is extracted from the strike direction, as shown in Figure 7. In the figure, the arrow represents the direction of fluid flow, and the digital label corresponds to the seepage contour. There is a redistribution of seepage pressure around the stope because of excavation unloading. As the goaf is approached, the water pressure decreases rapidly and tends to zero. Under the action of water head pressure difference, the flow is obviously flowing to the goaf, and the isoline of seepage pressure bends to the goaf, and finally presents the distribution similar to the precipitation funnel type.(3)Damage effect analysis: according to the analysis results of stress and seepage, the characteristic points in the direction of mining influence in the strike direction (Figure 8) were selected, and the change characteristics of their elastic modulus during excavation were monitored. As shown in Figure 9, the deformation strength of the surrounding rock shows a downward trend due to the damage caused by excavation. Because of the heterogeneity of the strength parameters, the initial values of the elastic modulus of the monitoring points are different. After excavation, the elastic modulus of the measuring points decreases by 8.7%, with an average decrease of 6.4%. The results of the numerical analysis and field test of Zhu et al. [23] and Ji [24] (excavation area of elastic modulus damage decreases by an average of 4%–14%) are consistent, which reflects the rationality of damage on the strength parameters of rock mass deterioration. As a result, a damage variable is defined.(4)Failure analysis and test verification: the failure is the result of coupling of stress, seepage, damage, disturbance, and so on. The failure characteristics of mining were analyzed and compared with the measured results so that the validity of the improved multifield coupling model can be verified. In numerical analysis, the distribution and coalescence of the plastic zone are often used for judging the failure of rock mass. In essence, the plastic characterization is that the element enters the state of yielding and irreversible deformation, which is not enough to reflect the damage degree of the unit. Therefore, a rock mass failure scenario can be quantitatively judged based on the damage approaching degree function defined in the previous paper. Using CMS (cavity monitoring system) to carry out field measurements in the mined out area, a three-dimensional measurement model was obtained (Figure 10). At the same time, the conventional multifield coupling and excavation simulation under the single stress field were carried out. The section with the most severe failure degree is extracted and compared with the numerical simulation results, as shown in Figure 11.
(a)
(b)
As shown in Figure 11, the red contour is the measured failure boundary, and the numerical simulation results only show the area where the break approaching degree is greater than two (damage zone). Mining failure is mainly distributed on both sides of the side and presents an asymmetric shape. The most serious damage distance from the design boundary reaches to 1.98 m. The numerical simulation of the control and improvement of the multifield model by considering the heterogeneity effect of blasting damage, strength parameters, and the failure form of the simulation also showed asymmetry. The damage range is significantly greater than the conventional single and multifield coupling model of the stress field and has a high degree of agreement with the measured results. The failure range of the conventional multifield coupling and single field model is obviously smaller than the measured range, and the damage presents symmetry. The comparison results verify the effectiveness and rationality of the improved multifield coupling model and the calculation program.
6. Conclusion
(1)Aiming at the shortcomings of the conventional multifield coupling model, the phase definition of elastoplastic damage variables and the equivalent of blasting load were taken into account, while the heterogeneity and anisotropy of the permeability of rock strength parameters were also considered. An improved multifield coupling model that takes damage and disturbances into account was established, and a corresponding numerical calculation program was compiled.(2)The improved multifield coupling model was used to analyze the response characteristics of excavation stress, seepage, and the failure of a stope. The obtained results were compared with other numerical simulation conditions and field measurement results. The results show that the improved multifield coupling model can reflect the extent of mining damage more accurately than the conventional coupling and single field model and has a high degree of agreement with field measured result. This verifies the validity and rationality of the improved model and its calculation program.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was supported by the National Key R&D Program of China during the Thirteenth Five-Year Plan period (2017YFC0602901); National Natural Science Foundation of China (51274250); and the Fundamental Research Funds for the Central Universities of Central South University (2017zzts204).