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

Authors: An Gang-jian; Cao Zheng-zheng; Duan Hong-fei; Li Guo-sheng; Li Ao
Journal: Geofluids
 (2022)

---

## Abstract

In order to reveal the mechanical response of surrounding rock of karst tunnel under stress-damage-seepage coupling effect, a new damage constitutive mechanical model of surrounding rock of karst tunnel under stress-damage-seepage coupling effect 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 them 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 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 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. 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 floor. 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.


---

## Body

## 1. 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 characterized 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–16].The multiple coupling effects exist in the engineering problems of undesirable geological sections, with complex disastrous mechanisms and high construction risks [17–19]. It is obvious that the scholars at home and abroad have achieved many research results in water inrush disaster of karst tunnel [20–23]. However, the mechanical response of surrounding rock of karst tunnel under stress-damage-seepage 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.


## 2. 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(1)σ∗=σ1−D,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(2)fx=mnxnm−1exp−xnm,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(3)D=∫−∞xftdt.The damage variable is yielded in equation (4), by combining equation (2) with equation (3).
(4)D=1−exp−xnm.Hooke’s law is satisfied in the stress-strain relationship of the undamaged part of rock mass in the principal stress state, shown in(5)εi=1Eσi∗−vσj∗+σk∗,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).
(6)σ1=Eε1exp−xnm+vσ2+σ3.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.
(7)σ1=Eε1exp−xnm+2vσ3.According to the principle of effective stress, the relationship is satisfied in(8)σij=σij′+pwδij,where pw 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
(9)σ1=Eε1exp−xnm+2vσ3+1−2vpw.


### 2.2. Parameter Distribution of Representative Cells

Based on the statistical law of mechanical parameters of rock mass, the Weibull distribution probability density function is employed, according to the mean values of physical and mechanical parameters of rock mass, such as elastic modulus and uniaxial (triaxial) compressive strength. A series of data conforming to the Weibull distribution are generated with the MATLAB program, and microelement of rock mass is assigned to corresponding positions, which describes the heterogeneity of rock mass. The Weibull distribution function is obtained in(10)fx=mnxnm−1exp−xnm.Therefore, the cumulative distribution function is obtained in(11)Fx=∫0xftdt=1−e−x/nm.According to the inverse function, ifn and m are the given values, random numbers are generated in
(12)x=n−log1−Fn1/m.The expectation of the Weibull distribution function is obtained in(13)Ex=nΓ1+1m,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.Whenn 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.Figure 1
Weibull probability density function whenn is 2 and m are 1, 2, 3, 5, and 7.Whenm 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.Figure 2
Weibull probability density function whenm is 2 and n are 1, 1.2, 1.5, 2.0, and 3.0.


### 2.3. 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 representative 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, cross-coupling 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(3)
According to the maximum tensile stress theory and the Mohr-Coulomb criterion, the damage of rock mass microunit is judged in each loading step. If the damage of rock mass occurs, the mechanical parameters of rock mass microunit are modified according to the damage theory, which are recalculated under the existing mechanical conditions until no damage occurs in this loading step. The modified and recalculated process is the computational convergence process in the numerical software COMSOL Multiphysics. Therefore, the final result is the true state of rock mass in this loading step(4)
Continue to load step by step, and repeat the previous step until the model totally fails in the last loading step. Since the loading value increases as the loading step increases, the description of (3) is the final damage state of rock mass in this loading step, and the description of (4) is the final damage state of rock mass in the last loading step, corresponding to the maximum loading value


### 2.4. Parameter Heterogeneity of Rock Mass

Based on the damage constitutive mechanical model of rock mass, a numerical simulation model of uniaxial compression test of rock mass is established in COMSOL Multiphysics, which is employed to discuss the influence of elastic modulus, tensile strength, and compressive strength on stress and damage of rock mass in different parameter heterogeneity.(1)
Numerical simulation modelAccording 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 widthw 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.Figure 3
Numerical model of uniaxial compression of rock mass.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 of1mm×1mm 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 strength is 6.04 MPa, the density is 2500 kg/m3, Poisson’s ratio is 0.25, and internal friction angle is 33.7°. The source of the corresponding parameters is the laboratory experiment of rock mass, and the processing and experiment equipment is shown in Figure 4. The selection basis is the national standard of the People’s Republic of China GBT23561.1-2009 Method for determination of physical and mechanical properties of coal and rock.
(2)
Result analysis of numerical simulationFigure 4
Processing and experiment equipment.
(a)
Vertical core drilling machine(b)
Automatic rock cutter machine(c)
Double face grinding machine(d)
Electronic universal testing machineAccording 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.Figure 5
Elastic modulus distribution of microelement of rock mass with different homogeneity.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 Figure6. 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.Figure 6
Mises stress distribution of microelement of rock mass with different homogeneity.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 Figure7. 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. 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.Figure 7
Strain energy density distribution of microelement of rock mass with different homogeneity.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 Figure8. 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.Figure 8
The number of damage points per loading step of rock mass with different homogeneity.


## 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(1)σ∗=σ1−D,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(2)fx=mnxnm−1exp−xnm,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(3)D=∫−∞xftdt.The damage variable is yielded in equation (4), by combining equation (2) with equation (3).
(4)D=1−exp−xnm.Hooke’s law is satisfied in the stress-strain relationship of the undamaged part of rock mass in the principal stress state, shown in(5)εi=1Eσi∗−vσj∗+σk∗,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).
(6)σ1=Eε1exp−xnm+vσ2+σ3.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.
(7)σ1=Eε1exp−xnm+2vσ3.According to the principle of effective stress, the relationship is satisfied in(8)σij=σij′+pwδij,where pw 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
(9)σ1=Eε1exp−xnm+2vσ3+1−2vpw.


## 2.2. Parameter Distribution of Representative Cells

Based on the statistical law of mechanical parameters of rock mass, the Weibull distribution probability density function is employed, according to the mean values of physical and mechanical parameters of rock mass, such as elastic modulus and uniaxial (triaxial) compressive strength. A series of data conforming to the Weibull distribution are generated with the MATLAB program, and microelement of rock mass is assigned to corresponding positions, which describes the heterogeneity of rock mass. The Weibull distribution function is obtained in(10)fx=mnxnm−1exp−xnm.Therefore, the cumulative distribution function is obtained in(11)Fx=∫0xftdt=1−e−x/nm.According to the inverse function, ifn and m are the given values, random numbers are generated in
(12)x=n−log1−Fn1/m.The expectation of the Weibull distribution function is obtained in(13)Ex=nΓ1+1m,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.Whenn 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.Figure 1
Weibull probability density function whenn is 2 and m are 1, 2, 3, 5, and 7.Whenm 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.Figure 2
Weibull probability density function whenm is 2 and n are 1, 1.2, 1.5, 2.0, and 3.0.


## 2.3. 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 representative 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, cross-coupling 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(3)
According to the maximum tensile stress theory and the Mohr-Coulomb criterion, the damage of rock mass microunit is judged in each loading step. If the damage of rock mass occurs, the mechanical parameters of rock mass microunit are modified according to the damage theory, which are recalculated under the existing mechanical conditions until no damage occurs in this loading step. The modified and recalculated process is the computational convergence process in the numerical software COMSOL Multiphysics. Therefore, the final result is the true state of rock mass in this loading step(4)
Continue to load step by step, and repeat the previous step until the model totally fails in the last loading step. Since the loading value increases as the loading step increases, the description of (3) is the final damage state of rock mass in this loading step, and the description of (4) is the final damage state of rock mass in the last loading step, corresponding to the maximum loading value


## 2.4. Parameter Heterogeneity of Rock Mass

Based on the damage constitutive mechanical model of rock mass, a numerical simulation model of uniaxial compression test of rock mass is established in COMSOL Multiphysics, which is employed to discuss the influence of elastic modulus, tensile strength, and compressive strength on stress and damage of rock mass in different parameter heterogeneity.(1)
Numerical simulation modelAccording 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 widthw 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.Figure 3
Numerical model of uniaxial compression of rock mass.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 of1mm×1mm 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 strength is 6.04 MPa, the density is 2500 kg/m3, Poisson’s ratio is 0.25, and internal friction angle is 33.7°. The source of the corresponding parameters is the laboratory experiment of rock mass, and the processing and experiment equipment is shown in Figure 4. The selection basis is the national standard of the People’s Republic of China GBT23561.1-2009 Method for determination of physical and mechanical properties of coal and rock.
(2)
Result analysis of numerical simulationFigure 4
Processing and experiment equipment.
(a)
Vertical core drilling machine(b)
Automatic rock cutter machine(c)
Double face grinding machine(d)
Electronic universal testing machineAccording 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.Figure 5
Elastic modulus distribution of microelement of rock mass with different homogeneity.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 Figure6. 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.Figure 6
Mises stress distribution of microelement of rock mass with different homogeneity.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 Figure7. 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. 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.Figure 7
Strain energy density distribution of microelement of rock mass with different homogeneity.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 Figure8. 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.Figure 8
The number of damage points per loading step of rock mass with different homogeneity.


## 3. 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.


### 3.1. 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 Figure9. The karst tunnel is located in the middle of the 50m×50m 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 Ft represents the self-weight of rock mass. Therefore, taking the proportionality coefficient γ=σx/Ft, four numerical simulation schemes are set, and the proportionality coefficient is 0.1, 0.5, 1, and 2, respectively.Figure 9
The numerical simulation model of surrounding rock of karst tunnel.Besides, the upper boundary is the bottom of the confined aquifer, and the pore water pressurePw is 4.6 MPa; the top line segment is the permeable boundary, and the right line segment, left line segment, and bottom line segment 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.


### 3.2. 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 Figure10. 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.Figure 10
The Mises stress distribution of rock mass of karst tunnel.
(a)
γ=0.1(b)
γ=0.5(c)
γ=1(d)
γ=2The damage evolution of rock mass of karst tunnel is shown in Figure11. 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 large Mises stress corresponds to the damaged vault and arch foot. Therefore, with the increase in γ, 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. The fissures in the rock mass are developed and the lithology is weak, and the groundwater in aquifer accelerates the rapid development of damaged area, thereby causing the water inrush disaster.Figure 11
The damage evolution of rock mass of karst tunnel.


### 3.3. Seepage Law of Surrounding Rock of Karst Tunnel

The permeability is an important mechanical parameter in the seepage field. The permeability distribution of rock mass of karst tunnel is shown in Figure12. When γ is small, namely, σx<Ft, 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 Ft, the k of the vault and floor increases, extending to a wider rock layer.Figure 12
The permeability distribution of rock mass of karst tunnel.The permeability from the aquifer to the surrounding rock of the vault of karst tunnel is reflected by the permeability in them-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<Ft at the vault of karst tunnel, the permeability reaches a peak and then decreases; when σx>Ft, the permeability continues to increase.Figure 13
The permeability distribution at different depths in them-n line.In order to research the permeability of rock mass under the tunnel floor, a horizontalk-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.Figure 14
The permeability distribution at different locations in thek-j line.


## 3.1. 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 Figure9. The karst tunnel is located in the middle of the 50m×50m 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 Ft represents the self-weight of rock mass. Therefore, taking the proportionality coefficient γ=σx/Ft, four numerical simulation schemes are set, and the proportionality coefficient is 0.1, 0.5, 1, and 2, respectively.Figure 9
The numerical simulation model of surrounding rock of karst tunnel.Besides, the upper boundary is the bottom of the confined aquifer, and the pore water pressurePw is 4.6 MPa; the top line segment is the permeable boundary, and the right line segment, left line segment, and bottom line segment 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.


## 3.2. 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 Figure10. 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.Figure 10
The Mises stress distribution of rock mass of karst tunnel.
(a)
γ=0.1(b)
γ=0.5(c)
γ=1(d)
γ=2The damage evolution of rock mass of karst tunnel is shown in Figure11. 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 large Mises stress corresponds to the damaged vault and arch foot. Therefore, with the increase in γ, 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. The fissures in the rock mass are developed and the lithology is weak, and the groundwater in aquifer accelerates the rapid development of damaged area, thereby causing the water inrush disaster.Figure 11
The damage evolution of rock mass of karst tunnel.


## 3.3. Seepage Law of Surrounding Rock of Karst Tunnel

The permeability is an important mechanical parameter in the seepage field. The permeability distribution of rock mass of karst tunnel is shown in Figure12. When γ is small, namely, σx<Ft, 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 Ft, the k of the vault and floor increases, extending to a wider rock layer.Figure 12
The permeability distribution of rock mass of karst tunnel.The permeability from the aquifer to the surrounding rock of the vault of karst tunnel is reflected by the permeability in them-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<Ft at the vault of karst tunnel, the permeability reaches a peak and then decreases; when σx>Ft, the permeability continues to increase.Figure 13
The permeability distribution at different depths in them-n line.In order to research the permeability of rock mass under the tunnel floor, a horizontalk-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.Figure 14
The permeability distribution at different locations in thek-j line.


## 4. Conclusions

(1)
When the mechanical parameters are assigned to microelements of rock mass by the Weibull distribution function, the larger them 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 them-n line increases with depth. When σx<Ft at the vault of karst tunnel, the permeability reaches a peak and then decreases; when σx>Ft, 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



---

*Source: 6879808-2022-08-13.xml*