Mechanical Response of Surrounding Rock of Karst Tunnel under Stress-Damage-Seepage Coupling Effect

In order to reveal the mechanical response of surrounding rock of karst tunnel under stress-damage-seepage coupling e ﬀ ect, a new damage constitutive mechanical model of surrounding rock of karst tunnel under stress-damage-seepage coupling e ﬀ ect is established, which is calculated by the COMSOL Multiphysics in this paper. When the mechanical parameters are assigned to microelements of rock mass by the Weibull distribution function, the larger the m value is, the more homogeneous the rock mass is. The variation trend of strain energy density with m is similar to equivalent stress, which increases ﬁ rstly and then decreases. The number of damage points increases with the increase in loading step and decreases rapidly after reaching the peak value and then remains a small number in the later loading stage. With the increase in γ , the stress range expands to the rock mass above the vault and below the ﬂ oor; the stress value increases signi ﬁ cantly, and the surrounding rock of karst tunnel is closer to strength limit, leading to the damage of rock mass. With the increase in γ , the area of the damage area in the upper part of the vault becomes larger, and most of the rock mass below the bottom plate is damaged; the damage area is semicircular, which indicates that both places are damaged by shearing action, resulting in the developed ﬁ ssures. Besides, there are the distribution characteristics of “ high value on both sides with a peak value and low value in the middle position ” in the permeability distribution, and high permeability is located at the arch foot, and the low permeability is located at the ﬂ oor. The larger the value of γ , the larger the permeability. The research achievements provide an important theoretical basis for prediction and treatment for dynamical disaster of karst tunnel.


Introduction
The wide distribution of karst landscapes poses a serious threat for the safety construction of transportation engineering in China [1,2]. Underground engineering in karst areas encounters more complex geological conditions, with highly developed karst fissures. The underground construction inevitably produces disturbance to the surrounding rock, inducing karst fissures in the karst water-rich area [3,4]. It is likely to cause engineering accidents such as water and mud inrush [5]. As one of the key geological dynamical disasters in tunnel construction, karst water inrush is char-acterized by the strong suddenness, fast catastrophic evolution, and large impact range [6,7].
It is important to analyze the formation and evolutionary mechanism of karst water inrush disaster and grasp the critical disastrous conditions and accurately predict and control karst water inrush to ensure the safe construction of the tunnel project. Tian et al. [8] employed the flow attenuation analysis method to study the relationship between the characteristics of the porosity pipeline underground river system and the characteristics of flow attenuation. Huang et al. [9] conducted an in-depth discussion on the water inrush disaster mechanism of tunnel filling karst pipeline and found that the groundwater has a great influence on the safety of karst pipeline. Zhang et al. [10] believed that the reason for the water inrush and mud inflow of the tunnel is that there is a weak "water barrier" between the tunnel bottom and the karst pipeline. Shen [11] analyzed the effects of different pressures on pipeline karst phenomena. In the process of construction of the Yiwan railway, karst water inrush disaster occurred frequently and seriously affected the construction of the tunnel project. Besides, many scholars summarized the effective and successful engineering experiences for handling with water inrush disaster in tunnel [12][13][14][15][16].
The multiple coupling effects exist in the engineering problems of undesirable geological sections, with complex disastrous mechanisms and high construction risks [17][18][19]. It is obvious that the scholars at home and abroad have achieved many research results in water inrush disaster of karst tunnel [20][21][22][23]. However, the mechanical response of surrounding rock of karst tunnel under stress-damageseepage coupling effect is not considered and researched systematically. In this paper, the parameter heterogeneity of rock mass is introduced in the damage constitutive mechanical model of rock mass, and the COMSOL Multiphysics is employed to verify the reliability and accuracy of the mechanical model. Besides, the damage evolution and seepage law of rock mass of karst tunnel are researched by establishing the numerical simulation model of surrounding rock of karst tunnel. The research results provide a theoretical support and judgment basis for the design, construction, and operation of karst tunnel.

Heterogeneity Theory of Parameter
Distribution of Rock Mass 2.1. Damage Constitutive Model of Rock Mass. Rock mass is a complex material containing joints and fissures, and the stress-strain relationship shows the nonlinearity characteristics under the action of external force. Therefore, it is of great significance to establish a damage constitutive model of rock mass to analyze the mechanical property of rock mass during the damage process. According to the Lemaitre equivalent strain hypothesis, the basic relational expression of damage constitutive model of rock mass is expressed in where σ * is the effective stress of the undamaged part of rock mass, σ is the nominal stress, and D is the damage variable. In order to describe the heterogeneity characteristics of rock mass, it is important to assign mechanical parameters of rock mass, such as elastic modulus, strength, and permeability. Assuming that the Weibull distribution rule is obeyed, the probability density function is obtained in where m represents the heterogeneity degree of rock mass, x is a random distribution variable of rock mass parameters, and n is the average value of rock mass parameters.
Assuming that the probability density function is obeyed in the failure probability of rock mass microelement, the damage variable is obtained in The damage variable is yielded in equation (4), by combining equation (2) with equation (3).
Hooke's law is satisfied in the stress-strain relationship of the undamaged part of rock mass in the principal stress state, shown in where v is Poisson's ratio, ε is the strain, and E is the elastic modulus.
Based on the Weibull distribution, the damage constitutive model of rock mass is obtained in equation (6), by combining equations (1) and (4) with equation (5).
The second principal stress is equal to the third principal stress in the conventional triaxial compression test of rock mass, and equation (7) is obtained.
According to the principle of effective stress, the relationship is satisfied in where p w is the pore water pressure, α is the Biot coefficient, σ ij ′ is the effective stress of the undamaged part of rock mass considering the pore water pressure, σ ij is the total stress, and δ ij is the Kronecker symbol. The stress measured by the test is the effective stress of undamaged part of rock mass considering the pore water pressure, which is substituted into equation (8). Based on the Weibull distribution, the damage constitutive model of rock mass in seepage is obtained in Therefore, the cumulative distribution function is obtained in According to the inverse function, if n and m are the given values, random numbers are generated in The expectation of the Weibull distribution function is obtained in where x is the independent variable, representing the generated physical and mechanical parameter; n is the scale parameter; m is the shape parameter, representing the rock heterogeneity; and Γ is the gamma function. When n is 2 and m are 1, 2, 3, 5, and 7, respectively, the Weibull distribution function is shown in Figure 1. When m is 1, the Weibull distribution function is an exponential function. The larger the shape parameter m is, the more random numbers are concentrated in the middle of the curve.
When m is 2 and n are 1, 1.2, 1.5, 2, and 3, respectively, the Weibull distribution function is shown in Figure 2. When the shape parameter m is constant, the shapes of the five curves are similar, but becoming gentle with the increase in the scale parameter n, and the abscissa corresponding to the highest point of the curve increases.

Solution
Steps of Numerical Simulation. COMSOL Multiphysics is finite element numerical software for solving partial differential equations and is widely employed to analyze the interaction of multiphysics fields. Besides, COMSOL Multiphysics is developed by MATLAB, and both programming languages are the same. Therefore, the problems of rock damage under the coupled action of stress and seepage fields are studied by COMSOL Multiphysics and MATLAB software. The solution steps to build and analyze the numerical simulation results are set.
(1) The geometric model is established in combination with the actual situation and divided into representa-tive units, namely, rock microunits, and then, the physical and mechanical parameters of the corresponding microunits, such as elastic modulus, tensile strength, and compressive strength, are generated from the Weibull distribution (2) The physical and mechanical parameters are assigned to the geometric model through the interpolation function. Besides, the material parameters, crosscoupling terms, initial conditions, and boundary conditions of the multiphysical fields are set, and the mesh is divided reasonably; thus, the incremental iteration method is employed to solve the problem   (1) Numerical simulation model According to the requirements of the uniaxial compression test of rock mass, the standard cylinder rock sample is established, which is simplified as a plane strain problem in COMSOL Multiphysics. The width w of the numerical model is 50 mm, and the height H is 100 mm. Therefore, the numerical model of uniaxial compression of rock mass is shown in Figure 3.
The left and right boundaries are set as free boundaries in the solid mechanic field, and the bottom boundary is set as roller support, and the displacement is applied on the top boundary. In consideration of the nonuniformity of rock mass, according to the REV theory, the numerical model of rock mass is divided into 5000 rock microelements with a size of 1 mm × 1 mm and a total of 5151 nodes. Besides, the parameters of microelement nodes are generated by the Monte-Carlo random simulation method and Weibull distribution function embedded in MATLAB. The mechanical parameters of rock mass are given by the average mechanical parameters of rock mass. Specifically, the elastic modulus is 34.5 GPa, the compressive strength is 100 MPa, the tensile   (2) Result analysis of numerical simulation According to the Weibull distribution function, m is a parameter characterizing heterogeneity characteristic of rock mass. In this section, the numerical simulation model of uniaxial compression of rock mass is set when the value of m is 2, 3, and 5, respectively. In order to research the effect of the Weibull distribution function on the heterogeneity characteristic of rock mass, the elastic modulus is taken as an example. The average value of elastic modulus of rock mass is 34.5 GPa. The generated random number is assigned to each microelement of rock mass, and the elastic modulus distribution of microelement of rock mass with different homogeneity is obtained in Figure 5. It is obvious that the maximum value of elastic modulus decreases and the minimum value increases with the increase in m; therefore, the distribution range of elastic modulus is shrinking. Besides, the color of the cloud image is closer to the middle value of legend in Figure 5, which indicates that the elastic modulus of microelements of rock mass tends to be the average value. When the value of m becomes large, the legend range shrinks and the rock mass tends to be homogeneous.
The three principal stresses are considered in the Mises stress, which is suitable for a reasonable yield criterion of rock mass. The Mises stress distribution of microelement of rock mass with different homogeneity is shown in Figure 6. It is obvious that the maximum value of Mises stress increases with the increase in homogeneity value, when homogeneity value is 2, 3, and 4. Besides, the Mises stress value in homogeneity value m = 5 is less than that in m = 4 and larger than that in m = 3, and the Mises stress value in homogeneity value m = 1 is less than that in m = 3 and larger than that in m = 2.
The strain energy density is the energy stored by elastic deformation of rock mass per unit volume under external force. The strain energy density distribution of microelement of rock mass with different homogeneity is shown in Figure 7. It is obvious that the strain energy density of rock mass increases firstly and then decreases with the increase in homogeneity value, similar to that of Mises stress.

Geofluids
Within a certain range of homogeneity value, the larger the homogeneity value of rock mass is, the higher the bearing capacity is; when the homogeneity value reaches a certain value, the bearing capacity of rock mass decreases. Based on the damage constitutive model of rock mass, the damage evolution and mechanical characteristics of rock mass are analyzed by counting the number of damage points in the numerical simulation model of rock mass under different homogeneity, which is shown in Figure 8. The x-axis is the loading step, while the y -axis is the number of damage points in rock mass. When homogeneity value m is larger, the number of peak damage points is larger, and the loading step reaching the peak value lags behind. It is obvious that the more homogeneous the rock mass is, the less likely it is destroyed. In case of macroscopic failure, more rock elements reach the strength limit. Before loading to step 4 and step 5, the displacement at the top of the model increases with the increase in loading step, and the damage inside the rock mass increases and gradually decreases after reaching the peak value of the damage point. In the last several loading steps, the number of damage points changes slightly relative to the total number of nodes and remains stable.

Damage Evolution and Seepage Law of Rock Mass of Karst Tunnel
According to the damage constitutive model of rock mass microelement and geological condition of karst tunnel, the damage evolution and seepage law of rock mass of karst tunnel under the action of in situ stresses and confined aquifer are researched, by establishing the numerical simulation model of surrounding rock of karst tunnel.

Numerical Simulation Model of Surrounding Rock of
Karst Tunnel. The numerical simulation model of surrounding rock of karst tunnel is established in COMSOL Multiphysics, shown in Figure 9. The karst tunnel is located in the middle of the 50 m × 50 m strata, with straight wall lining, and the section structure is not lined. The width of the bottom plate is 5 m, the height to the vault is 6 m, and the radius of the semicircular arch is 2.5 m. The vertical displacement is limited in the bottom boundary of the numerical model, the lateral initial in situ stress σ x is applied to the left and right boundaries, and the additional load F t represents the self-weight of rock mass. Therefore, taking the proportionality coefficient γ = σ x /F t , four numerical simulation schemes are set, and the proportionality coefficient is 0.1, 0.5, 1, and 2, respectively. Besides, the upper boundary is the bottom of the confined aquifer, and the pore water pressure P w is 4.6 MPa; the top line segment is the permeable boundary, and the right line segment, left line segment, and bottom line seg-ment are the impervious boundary. The air-side surface of karst tunnel is connected to the atmosphere, and the air pressure is 0.1 MPa. In order to ensure the calculation accuracy, the model is divided into free triangles, with a total of 4626 units. The line segment m-n for monitoring settlement is set from the aquifer to the vault, and the line segment k-j at 0.5 m below the bottom plate is set in the numerical simulation model.

Damage Evolution of Surrounding Rock of Karst Tunnel.
The Mises stress distribution of rock mass of karst tunnel under four working conditions is shown in Figure 10. When γ is less than 1, the equivalent stress of the side walls on both sides of karst tunnel is relatively large, especially in the arch shoulder and arch foot. The maximum values shown in the legend are located at the arch foot. With the increase in γ, the maximum stress increases, and the larger range of the side wall stress gradually decreases. When γ is greater than 1, the Mises stress of the tunnel vault and floor is larger; as γ increases, the stress range expands to the rock mass above the vault and below the floor; the stress value increases significantly, and the surrounding rock of karst tunnel is closer to strength limit, leading to the damage of rock mass.
The damage evolution of rock mass of karst tunnel is shown in Figure 11. When γ = 0:1, shear damage occurs in the spandrel, arch foot, and bottom center of karst tunnel; when γ = 0:5, the damage area decreases, indicating that the increase in horizontal initial stress alleviates the damage development to a certain extent. When γ = 1, the area with    The permeability is an important mechanical parameter in the seepage field. The permeability distribution of rock mass of karst tunnel is shown in Figure 12. When γ is small, namely, σ x < F t , k is larger at the tunnel spandrels and arches (green and red), and k is smaller at the vault and floor (dark blue). As γ increases, the k of the side walls on both sides of karst tunnel decreases and the range expands. When γ increases to 1, the k of the vault and the straight wall is reversed from the previous case (the colors are swapped), but the k remains at a large value at the foot of the vault. When γ continues to increase, namely, σ x is significantly larger than F t , the k of the vault and floor increases, extending to a wider rock layer.
The permeability from the aquifer to the surrounding rock of the vault of karst tunnel is reflected by the permeability in the m-n line. The permeability distribution at different depths in the m-n line is shown in Figure 13. It is obvious that the permeability increases with depth. When σ x < F t at the vault of karst tunnel, the permeability reaches a peak and then decreases; when σ x > F t , the permeability continues to increase.
In order to research the permeability of rock mass under the tunnel floor, a horizontal k-j line is set, and the permeability at different positions is obtained in Figure 14. The distribution characteristics of "high value on both sides with a peak value and low value in the middle position" are obtained, and high permeability is located at the arch foot and the low permeability is located at the floor. Besides, the larger the value of γ is, the larger the permeability is.

Conclusions
(1) When the mechanical parameters are assigned to microelements of rock mass by the Weibull distribution function, the larger the m value is, the more homogeneous the rock mass is. The variation trend of strain energy density with m is similar to equivalent stress, which increases firstly and then decreases. The more homogeneous the rock is, the more the damage points are, and the peak value lags behind. The number of damage points increases with the increase in loading step and decreases rapidly after reaching the peak value and then remains a small number in the later loading stage (2) With the increase in γ, the stress range expands to the rock mass above the vault and below the floor; the stress value increases significantly, and the surrounding rock of karst tunnel is closer to strength limit, leading to the damage of rock mass. With the increase in γ, the area of the damage area in the upper part of the vault becomes larger, and most of the rock mass below the bottom plate is damaged; the damage area is semicircular, which indicates that both places are damaged by shearing action, resulting in the developed fissures (3) The permeability at different depths in the m-n line increases with depth. When σ x < F t at the vault of karst tunnel, the permeability reaches a peak and then decreases; when σ x > F t , the permeability continues to increase. Besides, there is the distribution characteristics of "high value on both sides with a peak value and low value in the middle position" in the permeability distribution, and high permeability is located at the arch foot and the low permeability is located at the floor. The larger the value of γ is, the larger the permeability is

Data Availability
The data used to support the findings of this study are included within the article.

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