Multiphysics Coupling Model of Rock Mass considering Damage and Disturbance and Its Application

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. *e results show that the improved model reflects the extent and range of mining-induced failure more accurately and fits well with the actual detection. *ese 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.


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.e 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.en, 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.e 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.e 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.
erefore, 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.rough 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.

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).e 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, F is the volume force tensor, and u is the displacement tensor.
It is generally believed that the internal seepage of a rock mass is incompressible and obeys Darcy's law: where v is the seepage speed, Q m is the source item, p is the seepage pressure, μ is the dynamic viscosity coefficient of fluid, ρ is the fluid density, g is the gravitational acceleration, and K is the permeability tensor.e 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, C is the elastic stiffness tensor of rock mass, α is the effective stress coefficient of Biot, and ε v 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.erefore, it is necessary to improve the conventional multifield coupling model.

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.
e 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.e rock mass becomes an ideal elastoplastic body, which is inconsistent with real-world behavior.
erefore, 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.e 2 Advances in Civil Engineering expression of damage variable, D, is defined by the following stages: where ε e is the elastic strain, ε e0 is the mean value of elastic strain, ε p is the equivalent plastic strain, and s 1 and s 2 are correction coefficients.e 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 E 0 and c 0 are initial modulus of elasticity and cohesion and E(D) and c(D) are modulus of elasticity and cohesion considering damage.

Improved Multifield Coupling Model.
e damage variable, D, 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 C(D) K(D), we can reflect the effect of stress and seepage damage.In addition, αp in the rock constitutive equation reflects the effective stress effect of seepage, and the bulk strain source, zε v /zt, reflects the rock mass effect of deformation on seepage continuity.e 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 k 0 is the initial permeability, ω is the correlation coefficient (0.35 MPa −1 ), and σ ′ i 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.erefore, 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 I 1 is the first invariant of the stress tensor, J 2 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.e plastic strain increment is calculated by the associated flow rule: where c is the plastic multiplier, I � δ ij e i ⊗ e j , and s � σ − (I 1 /3)I.
e loading-unloading criterion is described as follows: e piecewise linear hardening function is used to approximate the nonlinear hardening function, and the expression is described as follows: where κ c 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.
erefore, 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.

Blasting Equivalent and Expression of Strength
Heterogeneity.
e 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.e method of uncoupled charge is usually used in deep hole blasting in a mine, and the equivalent load P s of blasting is calculated as follows [17]: where ρ 0 is the density of explosives, v b is the explosive velocity, c is the isentropic exponent of explosives, d c is the charge diameter, d b is the hole diameter, l e is the total length of charge, l b is the hole length, n is the pressure increase coe cient, and M and N are the attenuation coe cients of the blasting load.e Weibull probability density function [18] is used to characterize the heterogeneity of rock mass strength parameters: where ξ m are the rock mass strength parameters (such as elastic modulus, cohesive force, etc.), ξ m is the mean value of the strength parameter, and m is the correlation coe cient of inhomogeneity.
As shown in Figure 1, when the value of the nonhomogeneous coe cient, m, is higher, the value of the intensity parameter becomes more concentrated and tends to the mean value.

Failure Criterion.
For rock mass engineering, the concept of break approaching degree is introduced to quantitatively indicate the failure degree of rock mass units.
e 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.rough a series of geotechnical engineering applications, Zhang et al. [20] obtained the criterion of failure approaching degree: when the calculated value of F s is greater than 2, the rock mass element will be destroyed.

Numerical Calculation Programming
e improved multi eld coupling model is constructed by a set of di erential equations that include multiple dependent variables, while considering the heterogeneity of rock mass strength parameters and seepage anisotropy.e model is highly nonlinear, and the coupling relation is very complex.It is di cult to obtain an analytical solution by means of series variation or integration.erefore, the improved multi eld coupling model is solved using a numerical method.COMSOL Multiphysics is a nite element numerical analysis software speci cally designed to solve multi eld coupling problems.It also provides powerful programming functions.Based on the COMSOL numerical platform, a program for improving the multi eld 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 coe cient and stress, the de nition of the damage variable, the e ective stress, and the source of volume strain.
e validity and rationality of the improved multi eld coupling model and its calculation program will be veri ed by practical application.

Engineering
Situation.An underground room mining area is located between a −400 m level and −360 m level section.e stope ore reserves are approximately 22,200 t, the lower ore average angle is 55 °, and the ore grade (Pb + Zn) is 14%.e 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.e 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.e lower part is a nonpillar mining at-type chamber.e ore is transported by remote LHD, and the cavity is subsequently back lled.

Simulation Parameters and Conditions.
Several parameters (mainly from eld 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 4 Advances in Civil Engineering blasting numerical simulation.Based on gravity stress and tectonic stress, the initial geostress eld is obtained by using the multiple regression method.e expression of the heterogeneity of rock mass strength parameters is achieved by de ning 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.

Simulation Results, Analysis, and Veri cation
(1) Stress distribution characteristics: the stope is perpendicular to the direction of the ore body (northsouth 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.e mechanical response and failure characteristics of stope will directly a ect the stoping of the adjacent stope.erefore, the typical pro le 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 eld distribution obtained by numerical simulation is also heterogeneous.e 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.
e stress concentration area is mainly located at the top and bottom of the stope, and the maximum principal stress is 30.2MPa. e 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.
is causes tensile failure to occur.
(2) Seepage distribution characteristics: a typical seepage characteristic pro le is extracted from the strike direction, as shown in Figure 7.In the gure, the arrow represents the direction of uid ow, and the digital label corresponds to the seepage contour.ere is a redistribution of seepage pressure around the stope because of excavation unloading.As the

Advances in Civil Engineering
goaf is approached, the water pressure decreases rapidly and tends to zero.Under the action of water head pressure di erence, the ow is obviously owing to the goaf, and the isoline of seepage pressure bends to the goaf, and nally presents the distribution similar to the precipitation funnel type.
(3) Damage e ect analysis: according to the analysis results of stress and seepage, the characteristic points in the direction of mining in uence 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 di erent.After excavation, the elastic modulus of the measuring points decreases by 8.7%, with an average decrease of 6.4%.e results of the numerical analysis and eld 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 re ects the rationality of damage on the strength parameters of rock mass deterioration.As a result, a damage variable is de ned.(4) Failure analysis and test veri cation: the failure is the result of coupling of stress, seepage, damage, disturbance, and so on.e failure characteristics of mining were analyzed and compared with the 6 Advances in Civil Engineering measured results so that the validity of the improved multi eld coupling model can be veri ed.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 reect the damage degree of the unit.erefore, a rock mass failure scenario can be quantitatively judged based on the damage approaching degree function de ned in the previous paper.Using CMS (cavity monitoring system) to carry out eld measurements in the mined out area, a three-dimensional measurement model was obtained (Figure 10).At the same time, the conventional multi eld coupling and excavation simulation under the single stress eld were carried out.e section with the most severe failure degree is extracted and compared with the numerical simulation results, as shown in Figure 11.
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.
e most serious damage distance from the design boundary reaches to 1.98 m. e numerical simulation of the control and improvement of the multi eld model by considering the heterogeneity e ect of blasting damage, strength parameters, and the failure form of the simulation also showed asymmetry.
e damage range is signi cantly greater than the conventional single and multi eld coupling model of the stress eld and has a high degree of agreement with the measured results.
e failure range of the conventional multi eld coupling and single eld model is obviously smaller than the measured range, and the damage presents symmetry.
e comparison results verify the e ectiveness and rationality

Conclusion
(1) Aiming at the shortcomings of the conventional multi eld coupling model, the phase de nition 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 multi eld coupling model that takes damage and disturbances into account was established, and a corresponding numerical calculation program was compiled.(2) e improved multi eld coupling model was used to analyze the response characteristics of excavation stress, seepage, and the failure of a stope.e obtained results were compared with other numerical simulation conditions and eld measurement results.e results show that the improved multi eld coupling model can re ect the extent of mining damage more accurately than the conventional coupling and single eld model and has a high degree of agreement with eld measured result.is veri es the validity and rationality of the improved model and its calculation program.

Figure 11 :
Figure 11: Comparison of numerical failure characteristics with measured ones.