# Coupled Hydrologic-Mechanical-Damage Analysis and Its Application to Diversion Tunnels of Hydropower Station

Authors: Annan Jiang; Xiurong Yang; Mengfei Xu; Tengfei Jiang
Journal: Advances in Civil Engineering
 (2021)

---

## Abstract

Since the traditional model cannot sufficiently reflect the multifield coupling problem, this paper established an elastoplastic stress-seepage-damage analysis model considering the seepage field, stress field, and damage field. Simultaneously, the elastoplastic damage model involves many parameters and is difficult to determine. An inverse analysis program is compiled based on the differential evolution algorithm, and the surrounding rock damage parameters are inverted. Finally, the elastoplastic stress-seepage-damage coupling program and the damage parameter displacement back analysis program is compiled using C++ language. Then, the program is used to calculate the coupling problem of tunnel elastoplastic stress-seepage-damage. The results show that the proposed elastoplastic damage constitutive model can well describe the mechanical behavior of rock. The computational procedure can also simulate practical engineering problems, which can provide specific guidance for site construction.


---

## Body

## 1. Introduction

During tunnel construction on the steep slope of the dam abutment of the large-scale water conservancy and hydropower project in Southwest China, the shallow and deep rock mass stress redistribution caused by excavation usually leads to damage to the surrounding rock. The degree of damage is related to various factors, such as excavation methods, physical and mechanical properties of the rock mass, initial geostress field, and natural fracture distribution. During the tunnel excavation process, the rock mass permeability changes significantly with the initiation, expansion, and penetration of rock mass fissures. The rock mass is low-permeability before failure, and the seepage-stress coupling effect is not obvious. However, with the initiation and expansion of the fissures, the evolution process and interaction of the stress field, seepage field, and damage field inside the rock mass become very significant [1]. The deformation and failure of rock mass under the stress-seepage coupling are not only a frontier and hot issue in the development of basic science but also a key scientific problem to be solved in applied research [2].There is a lot of research on the seepage characteristics of tunnels surrounding rocks and underground caverns during construction [3–7]. Due to the stress-seepage coupling problem in the hydropower project and the flood-rich tunnel during construction, the disaster will lead to greater losses. Therefore, many scholars have conducted research on tunnel seepage problems in water-rich region projects. Liu et al. [8, 9] proposed a coupled seepage-erosion water inrush model based on classical theories of solute transport and fluid dynamics in porous media to investigate the characteristics of seepage-erosion properties. Wu et al. [10, 11] studied the characteristics of water flow and optimization of escape routes after water inrush in a post-open karst tunnel. Zhang et al. [12] studied the effect of Longsheng Reservoir on the seepage effect of the adjacent Kuzi Village Tunnel (located in Ulanchap, Inner Mongolia, North China).Because of the uncertainty of parameters of the fluid-solid coupling calculation, the back analysis method is studied by researchers. Based on the Levenberg–Marquit method of complex variable differentiation method, Liu et al. [13] established a multiparameter inversion method for the seepage-stress coupling problem with displacement information as a known quantity. Based on the full coupling analysis method for solving the coupling problem between stable seepage field and elastic displacement field, Wang Yuan [14] proposed a parameter inversion method for seepage and stressed full static coupling of a fractured rock mass.With the deepening of underground engineering activities, the research on the water-force coupling of the rock mass is more and more influenced by the damage evolution and pore fluid flow in rock mass [15–17]. Rich achievements have been made in the existing research about fluid-solid coupling calculation. However, there are also some problems to be solved. (1) There are many studies on the coupling of rock elastic brittle damage and seepage, and the model of the coupling between rock plastic damage and seepage is rare. (2) Existing models generally use a constant permeability coefficient, and the permeability coefficient changes less with the damage. (3) The inversion of coupling parameters for the coupling of the damage model and seepage involves complex optimization problems. However, the general inverse analysis optimization algorithm is easily limited to the local optimum and does not easily converge to the optimal global solution.Therefore, in this paper, the damage model method with a full coupling of the seepage field and the stress field is used based on the DP criterion and variable permeability coefficient method. The global optimization algorithm of difference evolution is introduced to back-analyze the damage parameters. The step iteration, numerical solving methods, and constitutive integration algorithm are studied. The elastoplastic stress-seepage-damage coupling program and intelligent back analysis program are compiled by C++ language. The programs are used to calculate the diversion tunnel of Shuibuya Hydropower Station. The distribution of stress field, seepage field, and damage field of the tunnel surrounding rock and the distribution curve of tunnel surrounding rock permeability coefficient are analyzed.


## 2. The Basic Theory of Elastoplastic Stress-Seepage-Damage Coupling of Rock

### 2.1. Rock Permeability

The permeability of a rock refers to the ease with which a gas, liquid, or ion passes through a rock. Due to the inherent porosity of the rock, there is a phenomenon that the liquid or gas migrates from the high pressure to the low pressure. The permeability of the rock mainly depends on the pore structure of the rock and the performance of the aggregate. Rock materials contain pores and cracks of various sizes, so porosity is one of the main factors affecting permeability.


### 2.2. Coupling Properties

#### 2.2.1. Seepage-Stress Coupling

Under the action of water pressure, the seepage of water acts on the rock with effective stress, which affects the crack of the rock stress. At the same time, the change of the rock stress field often leads to the fissure closure or expansion, affecting the penetration performance of fissures. Consequently, the seepage field is redistributed as the permeability of the fracture changes. This interaction is defined as seepage-stress coupling.


#### 2.2.2. Seepage-Damage Coupling

With the deepening of the understanding of the seepage coupling problem, people gradually realize that the damage and crack propagation have a significant effect on the seepage-stress coupling, mainly as follows: the impact of damage on the seepage process, the weakening of water, and the damage induced by osmotic stress. This problem is the seepage-damage coupling problem of rock during the rupture process. The seepage-stress coupling study focuses on the coupling methods for establishing different pore structure systems and describes their applicable conditions. The seepage-damage coupling is a combination of the above model and commonly used numerical calculation software, introducing medium fracture and damage judgment criteria, embedding in the characterization equation of the seepage-damage coupling of the medium to destroy the expansion zone, and investigating the seepage-damage coupling behavior in the engineering.


## 2.1. Rock Permeability

The permeability of a rock refers to the ease with which a gas, liquid, or ion passes through a rock. Due to the inherent porosity of the rock, there is a phenomenon that the liquid or gas migrates from the high pressure to the low pressure. The permeability of the rock mainly depends on the pore structure of the rock and the performance of the aggregate. Rock materials contain pores and cracks of various sizes, so porosity is one of the main factors affecting permeability.


## 2.2. Coupling Properties

### 2.2.1. Seepage-Stress Coupling

Under the action of water pressure, the seepage of water acts on the rock with effective stress, which affects the crack of the rock stress. At the same time, the change of the rock stress field often leads to the fissure closure or expansion, affecting the penetration performance of fissures. Consequently, the seepage field is redistributed as the permeability of the fracture changes. This interaction is defined as seepage-stress coupling.


### 2.2.2. Seepage-Damage Coupling

With the deepening of the understanding of the seepage coupling problem, people gradually realize that the damage and crack propagation have a significant effect on the seepage-stress coupling, mainly as follows: the impact of damage on the seepage process, the weakening of water, and the damage induced by osmotic stress. This problem is the seepage-damage coupling problem of rock during the rupture process. The seepage-stress coupling study focuses on the coupling methods for establishing different pore structure systems and describes their applicable conditions. The seepage-damage coupling is a combination of the above model and commonly used numerical calculation software, introducing medium fracture and damage judgment criteria, embedding in the characterization equation of the seepage-damage coupling of the medium to destroy the expansion zone, and investigating the seepage-damage coupling behavior in the engineering.


## 2.2.1. Seepage-Stress Coupling

Under the action of water pressure, the seepage of water acts on the rock with effective stress, which affects the crack of the rock stress. At the same time, the change of the rock stress field often leads to the fissure closure or expansion, affecting the penetration performance of fissures. Consequently, the seepage field is redistributed as the permeability of the fracture changes. This interaction is defined as seepage-stress coupling.


## 2.2.2. Seepage-Damage Coupling

With the deepening of the understanding of the seepage coupling problem, people gradually realize that the damage and crack propagation have a significant effect on the seepage-stress coupling, mainly as follows: the impact of damage on the seepage process, the weakening of water, and the damage induced by osmotic stress. This problem is the seepage-damage coupling problem of rock during the rupture process. The seepage-stress coupling study focuses on the coupling methods for establishing different pore structure systems and describes their applicable conditions. The seepage-damage coupling is a combination of the above model and commonly used numerical calculation software, introducing medium fracture and damage judgment criteria, embedding in the characterization equation of the seepage-damage coupling of the medium to destroy the expansion zone, and investigating the seepage-damage coupling behavior in the engineering.


## 3. Stress-Seepage-Damage Coupling Model for Rock

Based on the experimental results of the evolution of seepage law during rock failure, this section introduces the elastoplastic damage constitutive relation, damage variable, permeability coefficient evolution equation, and effective stress based on the elastoplastic theory, damage mechanics, and classical seepage mechanics. Then, a coupled model describing rock stress-seepage-damage is established. The model can effectively solve the stress-seepage-damage coupling problem of the tunnel surrounding rock.The assumptions in this paper are as follows. (1) The skeleton of the saturated body studied is an ideal elastoplastic damage isotropic body and satisfies the assumption of small deformation. (2) The seepage is laminar, following Darcy’s law. (3) The seepage fluid is incompressible, and the effects of temperature changes are negligible.


### 3.1. Rock Mass Mechanics Field Equation

A unit body is taken out at any point in the object, and each stress component on any one of the unit bodies should satisfy the static balance condition. Establish a balanced differential equation in a three-dimensional Cartesian coordinate system as follows:(1)∂σij∂xij+fj=0,i,j=1,2,3.The boundary value problem is shown in Figure1.Figure 1
Boundary value problem.Equation (1) can be composed of the following three parts:(2)Lu=σij,i+fj=0onΩ,σijni=t¯jonΓ1,um=um¯onΓ2,where fj is body force, σij is the total stress tensor, Ω is the problem-solving domain, ni is the normal border cosine, t¯j is the surface forces acting on the known boundary surface, Γ1 is the known force boundary, um¯ is the known displacement on the border, and Γ2 is the known displacement boundary.The relationship between stress and strain is different for different constitutive models. The constitutive equation is abbreviated as follows:(3)σij=fεijordσij=dfεij.The geometric equation is as follows:(4)εij=12ui,j+uj,i=ui,j,where ε is the strain tensor and u is the displacement vector.The boundary value problem is transformed into a problem of solving the displacementu while satisfying the boundary condition constraints. When the displacement u is solved, the strain and stress states can be solved by the deformation equation and the constitutive equation.The principle of effective stress in porous media is as follows:(5)σij′=σij−δijαp,0≤α≤1,where σij′ is the effective stress tensor (pressure is positive and pull is negative), p is the pore water pressure, α is the equivalent pore pressure coefficient, and δij is the Kronecker symbol.


### 3.2. Seepage Equation

Assuming that water is incompressible, according to Darcy’s law, the continuous equation of seepage under passive unsteady conditions [18] is as follows:(6)∂∂xkx∂h∂x+∂∂yky∂h∂y+∂∂zkz∂h∂z=Ss∂h∂t.Bringh=p/γ+z into equation (6); then,(7)∂∂xkx∂p∂x+∂∂yky∂p∂y+∂∂zkz∂p∂z=γSs∂h∂t,where h is the infiltration head, x,y,andz are the space coordinate, t is the time coordinate, kx,ky,andkz is the permeability coefficient of the x, y, and z axes as the main axis direction, Ss is the unit storage amount, γ is the water weight, and p is the pore water pressure.In the conditions with the free surface under nonpressure unsteady seepage, the initial condition is as follows:(8)h|t=0=h0x,y,z,0.The boundary conditions, i.e., the head boundary and the flow boundary, are as follows:(9)h|Μ1=f1x,y,z,t,kn∂h∂n|Μ2=f2x,y,z,t,where Μ1 is the known head boundary, f1 is the known head boundary value, Μ2 is the known flow boundary, and f2 is the known flow boundary value.According to the effective stress principle (equation (5)) and the equilibrium condition (equation (1)), the equilibrium differential equation based on the principle of effective stress can be obtained. The equilibrium differential equation under the action of the seepage field embodies the dynamic coupling effect of stress seepage.


### 3.3. Damage Variable and the Damage Evolution Equation

This study considers the equivalent plastic strainε¯p as the evolution process of the rock damage variable [15]. Experimental research shows that, with the increase of ε¯p, the damage aggravates gradually. Damage variable D is a nonlinear function ε¯p. It can be expressed as the exponential function of equivalent plastic strain.The calculation of equivalent plastic strain can be shown as follows:(10)ε¯p=23εp1−εp22+εp2−εp32+εp3−εp12,where εp1, εp2, and εp3 are three principal plastic strain, respectively.The evolution equation of the corresponding damage variableD is as follows:(11)D=1−exp−κε¯p−ε¯0p,where the equivalent plastic strain threshold ε¯0p=0, namely, the equivalent plastic strain produced damage evolution and κ is the normal number from the test. Equation (11) shows that, with the increase of cumulative plastic strain, the damage evolution eventually stabilized.Combining the simultaneous stress field (equation (2)), the percolation differential equation (equation (6)), the damage evolution equation (equation (11)), and the corresponding initial and boundary conditions (equations (8) and (9)), the coupling model of the stress field, damage field, and seepage field can be established.


## 3.1. Rock Mass Mechanics Field Equation

A unit body is taken out at any point in the object, and each stress component on any one of the unit bodies should satisfy the static balance condition. Establish a balanced differential equation in a three-dimensional Cartesian coordinate system as follows:(1)∂σij∂xij+fj=0,i,j=1,2,3.The boundary value problem is shown in Figure1.Figure 1
Boundary value problem.Equation (1) can be composed of the following three parts:(2)Lu=σij,i+fj=0onΩ,σijni=t¯jonΓ1,um=um¯onΓ2,where fj is body force, σij is the total stress tensor, Ω is the problem-solving domain, ni is the normal border cosine, t¯j is the surface forces acting on the known boundary surface, Γ1 is the known force boundary, um¯ is the known displacement on the border, and Γ2 is the known displacement boundary.The relationship between stress and strain is different for different constitutive models. The constitutive equation is abbreviated as follows:(3)σij=fεijordσij=dfεij.The geometric equation is as follows:(4)εij=12ui,j+uj,i=ui,j,where ε is the strain tensor and u is the displacement vector.The boundary value problem is transformed into a problem of solving the displacementu while satisfying the boundary condition constraints. When the displacement u is solved, the strain and stress states can be solved by the deformation equation and the constitutive equation.The principle of effective stress in porous media is as follows:(5)σij′=σij−δijαp,0≤α≤1,where σij′ is the effective stress tensor (pressure is positive and pull is negative), p is the pore water pressure, α is the equivalent pore pressure coefficient, and δij is the Kronecker symbol.


## 3.2. Seepage Equation

Assuming that water is incompressible, according to Darcy’s law, the continuous equation of seepage under passive unsteady conditions [18] is as follows:(6)∂∂xkx∂h∂x+∂∂yky∂h∂y+∂∂zkz∂h∂z=Ss∂h∂t.Bringh=p/γ+z into equation (6); then,(7)∂∂xkx∂p∂x+∂∂yky∂p∂y+∂∂zkz∂p∂z=γSs∂h∂t,where h is the infiltration head, x,y,andz are the space coordinate, t is the time coordinate, kx,ky,andkz is the permeability coefficient of the x, y, and z axes as the main axis direction, Ss is the unit storage amount, γ is the water weight, and p is the pore water pressure.In the conditions with the free surface under nonpressure unsteady seepage, the initial condition is as follows:(8)h|t=0=h0x,y,z,0.The boundary conditions, i.e., the head boundary and the flow boundary, are as follows:(9)h|Μ1=f1x,y,z,t,kn∂h∂n|Μ2=f2x,y,z,t,where Μ1 is the known head boundary, f1 is the known head boundary value, Μ2 is the known flow boundary, and f2 is the known flow boundary value.According to the effective stress principle (equation (5)) and the equilibrium condition (equation (1)), the equilibrium differential equation based on the principle of effective stress can be obtained. The equilibrium differential equation under the action of the seepage field embodies the dynamic coupling effect of stress seepage.


## 3.3. Damage Variable and the Damage Evolution Equation

This study considers the equivalent plastic strainε¯p as the evolution process of the rock damage variable [15]. Experimental research shows that, with the increase of ε¯p, the damage aggravates gradually. Damage variable D is a nonlinear function ε¯p. It can be expressed as the exponential function of equivalent plastic strain.The calculation of equivalent plastic strain can be shown as follows:(10)ε¯p=23εp1−εp22+εp2−εp32+εp3−εp12,where εp1, εp2, and εp3 are three principal plastic strain, respectively.The evolution equation of the corresponding damage variableD is as follows:(11)D=1−exp−κε¯p−ε¯0p,where the equivalent plastic strain threshold ε¯0p=0, namely, the equivalent plastic strain produced damage evolution and κ is the normal number from the test. Equation (11) shows that, with the increase of cumulative plastic strain, the damage evolution eventually stabilized.Combining the simultaneous stress field (equation (2)), the percolation differential equation (equation (6)), the damage evolution equation (equation (11)), and the corresponding initial and boundary conditions (equations (8) and (9)), the coupling model of the stress field, damage field, and seepage field can be established.


## 4. Permeability Characteristics in the Process of Rock Damage

In the porous continuum model and the equivalent continuum model of rock mass seepage, the rock mass is regarded as a uniform medium composed of skeleton particles and pores (fractures). This structural feature makes the microscopic geometry of the rock mass change after the load or disturbance of the rock medium and the skeleton particles will be rearranged, resulting in changes in the porosity and permeability of the rock mass medium. The rock mass permeability coefficient should be considered as a variable in a fully coupled analysis, which is usually a function of porosity, strain, or stress. Guang-Ting [19] systematically summarized and introduced the three research methods of rock seepage stress coupling characteristics and reviewed the rationality and application of various results.The permeability coefficient-strain (or stress) equation is an indispensable governing equation for numerical analysis of seepage stress coupling. According to the different stress states of the rock mass, the permeability coefficient is defined as the function of stress and damage, and the evolution law of rock mass permeability coefficient in the elastic and plastic stages is realized. Under actual conditions, once the rock mass material yields and breaks, the rock mass permeability coefficient will increase sharply with the expansion and penetration of the original crack and the generation of a large number of new cracks. Therefore, some scholars have introduced the concept of area reduction rate or sudden jump coefficient to simulate the permeability change law of rock damage yield stage, but this does not directly affect the impact of the damage process on rock mass permeability change. Compared with current research, this study more realistically reflects the change of the permeability coefficient in the elastoplastic stress-seepage-damage coupling analysis of rock mass.The relational expression between the coefficient of permeability of rock mass and volume strain can be obtained based on the elastic stage Kozeny–Carman equation [20]:(12)K=K01+εv/n031+εv,where n0 is the initial porosity, K0 is the initial permeability of the rock, and εv is the volumetric strain.The permeability coefficient of the rock in the plastic phase is shown as follows [21]:(13)K=1−DKM+DKD1+εvpF3,where KM and KD are permeability coefficients of the undamaged rock and fractured rock, respectively, and εvpF=Dεvp, which is the plastic volumetric strain for the defect phase.


## 5. Numerical Solving and Step Iteration Method for Coupling Model

The solution to the elastic-plastic-damage-seepage coupling model of rock is a complex nonlinear problem. The difficulty of solving is mainly reflected in the relationship between elastoplasticity, damage, seepage calculation, and stress-damage-seepage interaction of rock mass. It is quite difficult to solve a problem involving so many nonlinear factors by only one iteration. After the above factors are summarized, this study iteratively solves them in a certain order, and many nonlinear problems can be solved in order to achieve the ultimate realization of the coupled model.Based on solid and seepage finite element theory, elastoplastic constitutive integral theory, and step-by-step iterative coupling method, the elastoplastic stress-seepage-damage program is compiled in this study.


### 5.1. Elastoplastic Damage Finite Element and Constitutive Integration Algorithm

Discrete equations (2)–(5) to obtain a finite element equation with displacement as an unknown or its incremental form can be obtained as follows:(14)Kmu=fmorKmΔuΔu=Δfm,where Km is the total stiffness matrix of the stress field, u is the node displacement column vector, Δu is increment displacement, and {fm} is a nodal force vector that includes physical strength, surface strength, and pore pressure of the equivalent load.The overall stiffness matrixK is assembled with the element stiffness matrix. KTe is the function of consistent tangent modulus C as follows:(15)KTe=∫ΩeBTCBdv,where B is the strain matrix, that is, the matrix for finding strain based on displacement and BT is the matrix transpose of B.For rock, plasticity mainly refers to frictional sliding between internal cracks or joint surfaces and damage refers to the occurrence and expansion of internal cracks. Plastic damage coupling has two meanings: (1) interact with each other through their potential functions (and loading functions); (2) interact with each other through their consistency conditions, i.e., the evolution of the two internal variables of plasticity and damage interact [21].The Drucker–Prager model is widely used in rock materials to describe the characteristics of plastic stress-deformation. The plastic yield function considering the damaging effect is as follows:(16)fpσ,ε¯p,D=J2sσ+ηpσ−1−Dξcε¯p,where pσ is the mean normal stress, pσ=1/3trσ, trσ is the trace of stress tensor (i.e., the algebraic sum of all principal stresses), J2 is the second invariant of deviatoric stress, J2=1/2s:s, sσ is the deviatoric tensor of stress, sσ=σ−pσI, I is the second-order tensor, I=δijei⊗ej, i,j=1,2,3, δij is the Kronecker symbol, cε¯p is the cohesion under the damage, D is the damage variable, and η and ξ are the material parameters.The influence of damage on the internal friction angle is very small, so only the effect of damage on cohesionc is considered. With the accumulation of damage, the plastic strain increases, and the cohesion decreases gradually. It can be described as a power function [22]:(17)cε¯p=cε¯p′−cε¯p′−crDζ,where cε¯p′ is the cohesion, cr is the cohesion of rock when it is obviously damaged, and ζ is the material parameter with a value between 0 and 1.


#### 5.1.1. Constitutive Integration Algorithm

The displacement increment can be obtained according to equation (14), and the strain increment can be calculated from displacement. Each iteration step calculates the stress increment from the given strain increment. The solution of the stress is related to the selection of the constitutive model.The study uses the constitutive model of Drucker–Prager (equation (16)). In the nonlinear finite element algorithm, each iteration step calculates the stress increment from the given strain increment. However, the updated values of stress and internal variables during the solution process are prone to not satisfy the yield function, resulting in inaccurate calculation results. This paper adopts the return mapping algorithm proposed by Simo [23], which has two processes of elastic prediction and plastic correction, and the algorithm is shown in Figure 2. Elastic prediction is to calculate the stress state σn+1trial according to the overall strain. If the stress state is outside the yield surface, then plasticity correction is driven by incremental plasticity factor, and the elastic trial stress returns to the yield surface so that plastic consistency is re-established in the updated state; then, σn+1 is obtained.Figure 2
The general form of the return mapping algorithm. (a) Strain softening model. (b) Strain hardening model.
(a)(b)


#### 5.1.2. Consistent Tangent Modulus

When it returns to the smooth conical surface, the consistent tangent modulus is derived as follows:(18)Cdep=2GDnaId+2GDnbT⊗T−2GDnKDnc+KDndI⊗I,where Cdep is the consistent tangent modulus, GDn and KDn are damage shear modulus and damage bulk modulus, respectively, Id is volume tensor, T is the second-order unit tensor parallel to the elastic predicted strain, a, b, c, and d are coefficients, εdn+1et is the strain during elastic prediction at the time of tn+1, Dn is the damage variable at the time of tn, and η and ζ are the material parameters:(19)a=1−Δγ2εdn+1et,b=Δγ2εdn+1et−GDnGDn+KDnηη¯+1−Dnζ2H,c=ηT⊗I+η¯I⊗T,d=1−KDnηη¯A.When it returns to a sharp point, the consistent tangent modulus is derived as follows:(20)Cdep=KDn1−KDnKDn+αβ1−DnHI⊗I,where H is the hardening modulus. The study adopts the associated flow rule, α and β are the parameters related to the internal friction angle, and they are chosen according to the required approximation to the Mohr–Coulomb criterion.


### 5.2. Seepage Finite Element Method

Disperse the seepage field into combinations of a finite number of units. Disperse equation (6) to solve the seepage field in the hydraulic head function h of finite element basic equations, and the forms are(21)Ksh=fs,where [Ks] is seepage matrix, h is hydraulic head column vector, and fs is a free term column vector. Thus, the original solution of the partial differential equation is used instead of solving algebraic equations.Based on the variational principle, the finite element solution program for seepage is compiled. The preprocedure process divides the mesh by means of ANSYS and then converts it into a seepage finite element program to solve the input file. The postprocessing is displayed by Tecplot, and the solution result of the seepage finite element program is converted into a format that can be displayed by the Tecplot software.


### 5.3. Elastic-Plastic Stress-Seepage-Damage Coupling Iteration

The coupling mechanism between the groundwater seepage field and stress field in the rock mass is a relatively complex dynamic process. The interaction between the stress field and the seepage field is linked by changes in the permeability of the rock mass. When the rock mass is disturbed and the permeability is changed, the two mechanisms corresponding to the stress field and the seepage field are repeated to achieve a dynamic stable state.There is much research related to the numerical methods of rock stress-seepage coupling which has given different coupling methods. There are two main methods, i.e., step-by-step iteration and one-time coupling. In this paper, the step-by-step iterative method is used to achieve elastoplastic stress-damage-seepage coupling.Under the initial stress state of the rock mass, the stress field, deformation field, and damage field of the elastoplastic damaged rock mass can be obtained by the incremental iteration calculation in Section5.1. The volumetric strain of the element is obtained from the calculation of the deformation field. Then, under the updated stress state, the permeation coefficient matrix of each unit is calculated according to the obtained volume strain (equations (12) and (13)), and the updated permeation coefficient matrix is subjected to the calculation of the seepage finite element in Section 5.2 to obtain the seepage field.Finally, the pore water pressure of the joint is calculated according to the calculation of the hydraulic head in Section5.2, and it is brought back to the calculation of the mechanical field of Section 5.1 by the principle of effective stress. The reciprocating progress is continued until the stress field, the damage field, and the seepage field are obtained twice before and after satisfying the convergence criterion. The elastoplastic stress-seepage-damage coupling procedure is shown in Figure 3.Figure 3
Elastic-plastic stress-seepage-damage coupling process.


## 5.1. Elastoplastic Damage Finite Element and Constitutive Integration Algorithm

Discrete equations (2)–(5) to obtain a finite element equation with displacement as an unknown or its incremental form can be obtained as follows:(14)Kmu=fmorKmΔuΔu=Δfm,where Km is the total stiffness matrix of the stress field, u is the node displacement column vector, Δu is increment displacement, and {fm} is a nodal force vector that includes physical strength, surface strength, and pore pressure of the equivalent load.The overall stiffness matrixK is assembled with the element stiffness matrix. KTe is the function of consistent tangent modulus C as follows:(15)KTe=∫ΩeBTCBdv,where B is the strain matrix, that is, the matrix for finding strain based on displacement and BT is the matrix transpose of B.For rock, plasticity mainly refers to frictional sliding between internal cracks or joint surfaces and damage refers to the occurrence and expansion of internal cracks. Plastic damage coupling has two meanings: (1) interact with each other through their potential functions (and loading functions); (2) interact with each other through their consistency conditions, i.e., the evolution of the two internal variables of plasticity and damage interact [21].The Drucker–Prager model is widely used in rock materials to describe the characteristics of plastic stress-deformation. The plastic yield function considering the damaging effect is as follows:(16)fpσ,ε¯p,D=J2sσ+ηpσ−1−Dξcε¯p,where pσ is the mean normal stress, pσ=1/3trσ, trσ is the trace of stress tensor (i.e., the algebraic sum of all principal stresses), J2 is the second invariant of deviatoric stress, J2=1/2s:s, sσ is the deviatoric tensor of stress, sσ=σ−pσI, I is the second-order tensor, I=δijei⊗ej, i,j=1,2,3, δij is the Kronecker symbol, cε¯p is the cohesion under the damage, D is the damage variable, and η and ξ are the material parameters.The influence of damage on the internal friction angle is very small, so only the effect of damage on cohesionc is considered. With the accumulation of damage, the plastic strain increases, and the cohesion decreases gradually. It can be described as a power function [22]:(17)cε¯p=cε¯p′−cε¯p′−crDζ,where cε¯p′ is the cohesion, cr is the cohesion of rock when it is obviously damaged, and ζ is the material parameter with a value between 0 and 1.


### 5.1.1. Constitutive Integration Algorithm

The displacement increment can be obtained according to equation (14), and the strain increment can be calculated from displacement. Each iteration step calculates the stress increment from the given strain increment. The solution of the stress is related to the selection of the constitutive model.The study uses the constitutive model of Drucker–Prager (equation (16)). In the nonlinear finite element algorithm, each iteration step calculates the stress increment from the given strain increment. However, the updated values of stress and internal variables during the solution process are prone to not satisfy the yield function, resulting in inaccurate calculation results. This paper adopts the return mapping algorithm proposed by Simo [23], which has two processes of elastic prediction and plastic correction, and the algorithm is shown in Figure 2. Elastic prediction is to calculate the stress state σn+1trial according to the overall strain. If the stress state is outside the yield surface, then plasticity correction is driven by incremental plasticity factor, and the elastic trial stress returns to the yield surface so that plastic consistency is re-established in the updated state; then, σn+1 is obtained.Figure 2
The general form of the return mapping algorithm. (a) Strain softening model. (b) Strain hardening model.
(a)(b)


### 5.1.2. Consistent Tangent Modulus

When it returns to the smooth conical surface, the consistent tangent modulus is derived as follows:(18)Cdep=2GDnaId+2GDnbT⊗T−2GDnKDnc+KDndI⊗I,where Cdep is the consistent tangent modulus, GDn and KDn are damage shear modulus and damage bulk modulus, respectively, Id is volume tensor, T is the second-order unit tensor parallel to the elastic predicted strain, a, b, c, and d are coefficients, εdn+1et is the strain during elastic prediction at the time of tn+1, Dn is the damage variable at the time of tn, and η and ζ are the material parameters:(19)a=1−Δγ2εdn+1et,b=Δγ2εdn+1et−GDnGDn+KDnηη¯+1−Dnζ2H,c=ηT⊗I+η¯I⊗T,d=1−KDnηη¯A.When it returns to a sharp point, the consistent tangent modulus is derived as follows:(20)Cdep=KDn1−KDnKDn+αβ1−DnHI⊗I,where H is the hardening modulus. The study adopts the associated flow rule, α and β are the parameters related to the internal friction angle, and they are chosen according to the required approximation to the Mohr–Coulomb criterion.


## 5.1.1. Constitutive Integration Algorithm

The displacement increment can be obtained according to equation (14), and the strain increment can be calculated from displacement. Each iteration step calculates the stress increment from the given strain increment. The solution of the stress is related to the selection of the constitutive model.The study uses the constitutive model of Drucker–Prager (equation (16)). In the nonlinear finite element algorithm, each iteration step calculates the stress increment from the given strain increment. However, the updated values of stress and internal variables during the solution process are prone to not satisfy the yield function, resulting in inaccurate calculation results. This paper adopts the return mapping algorithm proposed by Simo [23], which has two processes of elastic prediction and plastic correction, and the algorithm is shown in Figure 2. Elastic prediction is to calculate the stress state σn+1trial according to the overall strain. If the stress state is outside the yield surface, then plasticity correction is driven by incremental plasticity factor, and the elastic trial stress returns to the yield surface so that plastic consistency is re-established in the updated state; then, σn+1 is obtained.Figure 2
The general form of the return mapping algorithm. (a) Strain softening model. (b) Strain hardening model.
(a)(b)


## 5.1.2. Consistent Tangent Modulus

When it returns to the smooth conical surface, the consistent tangent modulus is derived as follows:(18)Cdep=2GDnaId+2GDnbT⊗T−2GDnKDnc+KDndI⊗I,where Cdep is the consistent tangent modulus, GDn and KDn are damage shear modulus and damage bulk modulus, respectively, Id is volume tensor, T is the second-order unit tensor parallel to the elastic predicted strain, a, b, c, and d are coefficients, εdn+1et is the strain during elastic prediction at the time of tn+1, Dn is the damage variable at the time of tn, and η and ζ are the material parameters:(19)a=1−Δγ2εdn+1et,b=Δγ2εdn+1et−GDnGDn+KDnηη¯+1−Dnζ2H,c=ηT⊗I+η¯I⊗T,d=1−KDnηη¯A.When it returns to a sharp point, the consistent tangent modulus is derived as follows:(20)Cdep=KDn1−KDnKDn+αβ1−DnHI⊗I,where H is the hardening modulus. The study adopts the associated flow rule, α and β are the parameters related to the internal friction angle, and they are chosen according to the required approximation to the Mohr–Coulomb criterion.


## 5.2. Seepage Finite Element Method

Disperse the seepage field into combinations of a finite number of units. Disperse equation (6) to solve the seepage field in the hydraulic head function h of finite element basic equations, and the forms are(21)Ksh=fs,where [Ks] is seepage matrix, h is hydraulic head column vector, and fs is a free term column vector. Thus, the original solution of the partial differential equation is used instead of solving algebraic equations.Based on the variational principle, the finite element solution program for seepage is compiled. The preprocedure process divides the mesh by means of ANSYS and then converts it into a seepage finite element program to solve the input file. The postprocessing is displayed by Tecplot, and the solution result of the seepage finite element program is converted into a format that can be displayed by the Tecplot software.


## 5.3. Elastic-Plastic Stress-Seepage-Damage Coupling Iteration

The coupling mechanism between the groundwater seepage field and stress field in the rock mass is a relatively complex dynamic process. The interaction between the stress field and the seepage field is linked by changes in the permeability of the rock mass. When the rock mass is disturbed and the permeability is changed, the two mechanisms corresponding to the stress field and the seepage field are repeated to achieve a dynamic stable state.There is much research related to the numerical methods of rock stress-seepage coupling which has given different coupling methods. There are two main methods, i.e., step-by-step iteration and one-time coupling. In this paper, the step-by-step iterative method is used to achieve elastoplastic stress-damage-seepage coupling.Under the initial stress state of the rock mass, the stress field, deformation field, and damage field of the elastoplastic damaged rock mass can be obtained by the incremental iteration calculation in Section5.1. The volumetric strain of the element is obtained from the calculation of the deformation field. Then, under the updated stress state, the permeation coefficient matrix of each unit is calculated according to the obtained volume strain (equations (12) and (13)), and the updated permeation coefficient matrix is subjected to the calculation of the seepage finite element in Section 5.2 to obtain the seepage field.Finally, the pore water pressure of the joint is calculated according to the calculation of the hydraulic head in Section5.2, and it is brought back to the calculation of the mechanical field of Section 5.1 by the principle of effective stress. The reciprocating progress is continued until the stress field, the damage field, and the seepage field are obtained twice before and after satisfying the convergence criterion. The elastoplastic stress-seepage-damage coupling procedure is shown in Figure 3.Figure 3
Elastic-plastic stress-seepage-damage coupling process.


## 6. Damage Parameters’ Inversion Based on Differential Evolution Algorithm

The damage parameters inverse problem turned into the optimization problem of constraint [24]:(22)minfx1,x2,…,xn=1m∑i=1mYi0−Yi2.The constraint condition:(23)xil≤xi≤xiu,i=1,2,…,n,where Yi0 is the observed displacement, Yi is the FEM calculated value by the seepage-stress-damage coupling model, m is the number of observations, xi is the damage parameter, and xil and xiu are the upper bound and lower bound on xi.The differential evolution (DE) algorithm is a global optimization algorithm proposed by Rainer Storn and Kenneth Price [24]. It is a new algorithm after genetic algorithm, ant colony algorithm, and particle swarm algorithm. It has great advantages in search success rate and computational efficiency, no requirement for the initial value, less controlled variables, fast convergence, good adaptability, optimization for multivariable complex problems, and no need to encode and decode operating.DE algorithm includes generating initial population, mutation and crossover operation, and selection operation. The specific principle is as follows.


### 6.1. Generate Initial Population

DE constitutesR-dimension solution vectors with the optimization problem directly. Each solution vector is basically individual evolution. Let the number of solution vectors in the Gth population solution be NP. Solution vectors are xi, G = (x1, x2, …xn), i = 1, 2, …, NP, where G means the group evolution of every generation and i means the individual location in the population. The initial population is produced according to the method of random that uniform distribution is in the solution space:(24)xi,1=xjL+randxjU−xjL,i=1,2,…,NP,where rand∈[0, 1] is a random number, xi,1 is the solution vector of first-generation, j = 1, 2, …, R, R is the number of the dimensions of solution vector, and xjU and xjL are the upper and lower bounds of the jth component, respectively.


### 6.2. Mutation Operation

The mutation operation uses a different strategy, which uses the difference vector between individuals in the population to perturb the individual to achieve individual variation. The size of the difference vector can be automatically adjusted according to the individual distribution within the population, and the adaptation is good. For each target vectorxi, G in the Gth generation, each vector individual contains R components, and the solution of the variance vector νi, G + 1 is(25)vi,G+1=xr1,G+F⋅xr2,G−xr3,G,where r1, r2, r3∈[1, 2, …, NP] are not the same random integers each other and not equal to i and F∈[0, 1] is the mutagenic factor that one of the main control parameters in the algorithm has to adjust the step magnitude of the vector difference.


### 6.3. Interlace Operation

Interlace operation is in order to increase the diversity of the population. The new test vectorui,G+1 can be calculated by hybridizing the target vector xi, G, and the variation vector νi, G + 1 according to the following rules:(26)ui,G+1=u1i,G+1,u2i,G+1,…,uDi,G+1,uji,G+1=vi,G+1,ifrandbj≤CRj=rnbri,xi,G,ifrandbj>CR&j=rnbri,where randbj ∈ [0, 1] is random decimals corresponding with the jth component, j = 1, 2, …, R, CR ∈ [0, 1] is crossed factors that is another one of the main control parameters, and rnbri is the coefficient corresponding with the ith vector, random integers among 1, 2, …, R.


### 6.4. Selecting Operation

The greedy search method is used to select whether to select the test vectorui,G+1 as the target vector of the (G + 1)th generation.Compare test vectorui,G+1 with target vector xi, G, and choose vector ui,G+1 if ui,G+1 corresponding the smaller target value. Otherwise, retain xi,G if xi,G is corresponding the smaller one. New population xi,G+1 can be generated after selection as follows:(27)xi,G+1=ui,G+1,iffui,G+1≤fxi,G,xi,G,iffui,G+1>fxi,G.The DE algorithm adopts real-coded and regards equation (22) as a fitness function in this paper. The back analysis calculation process is shown in Figure 4.Figure 4
Calculation process of back analysis.In this study, the self-developed elastoplastic stress-seepage-damage finite element solution program is named RMAST in the intelligent displacement back analysis program, which can calculate the displacementYi of equation (22) and complete the inversion of damage parameters.


## 6.1. Generate Initial Population

DE constitutesR-dimension solution vectors with the optimization problem directly. Each solution vector is basically individual evolution. Let the number of solution vectors in the Gth population solution be NP. Solution vectors are xi, G = (x1, x2, …xn), i = 1, 2, …, NP, where G means the group evolution of every generation and i means the individual location in the population. The initial population is produced according to the method of random that uniform distribution is in the solution space:(24)xi,1=xjL+randxjU−xjL,i=1,2,…,NP,where rand∈[0, 1] is a random number, xi,1 is the solution vector of first-generation, j = 1, 2, …, R, R is the number of the dimensions of solution vector, and xjU and xjL are the upper and lower bounds of the jth component, respectively.


## 6.2. Mutation Operation

The mutation operation uses a different strategy, which uses the difference vector between individuals in the population to perturb the individual to achieve individual variation. The size of the difference vector can be automatically adjusted according to the individual distribution within the population, and the adaptation is good. For each target vectorxi, G in the Gth generation, each vector individual contains R components, and the solution of the variance vector νi, G + 1 is(25)vi,G+1=xr1,G+F⋅xr2,G−xr3,G,where r1, r2, r3∈[1, 2, …, NP] are not the same random integers each other and not equal to i and F∈[0, 1] is the mutagenic factor that one of the main control parameters in the algorithm has to adjust the step magnitude of the vector difference.


## 6.3. Interlace Operation

Interlace operation is in order to increase the diversity of the population. The new test vectorui,G+1 can be calculated by hybridizing the target vector xi, G, and the variation vector νi, G + 1 according to the following rules:(26)ui,G+1=u1i,G+1,u2i,G+1,…,uDi,G+1,uji,G+1=vi,G+1,ifrandbj≤CRj=rnbri,xi,G,ifrandbj>CR&j=rnbri,where randbj ∈ [0, 1] is random decimals corresponding with the jth component, j = 1, 2, …, R, CR ∈ [0, 1] is crossed factors that is another one of the main control parameters, and rnbri is the coefficient corresponding with the ith vector, random integers among 1, 2, …, R.


## 6.4. Selecting Operation

The greedy search method is used to select whether to select the test vectorui,G+1 as the target vector of the (G + 1)th generation.Compare test vectorui,G+1 with target vector xi, G, and choose vector ui,G+1 if ui,G+1 corresponding the smaller target value. Otherwise, retain xi,G if xi,G is corresponding the smaller one. New population xi,G+1 can be generated after selection as follows:(27)xi,G+1=ui,G+1,iffui,G+1≤fxi,G,xi,G,iffui,G+1>fxi,G.The DE algorithm adopts real-coded and regards equation (22) as a fitness function in this paper. The back analysis calculation process is shown in Figure 4.Figure 4
Calculation process of back analysis.In this study, the self-developed elastoplastic stress-seepage-damage finite element solution program is named RMAST in the intelligent displacement back analysis program, which can calculate the displacementYi of equation (22) and complete the inversion of damage parameters.


## 7. Engineering Application

### 7.1. General State of Shuibuya Hydropower Station

The Shuibuya Water Conservancy Project is located in the middle section of Qingjiang River in Badong County, Hubei Province, and is one of the important power stations for the development of the Qingjiang River Basin, as shown in Figure5. The geological engineering conditions of the surrounding rock of the underground cavern are very complicated, and the tailwater tunnel section is particularly serious, passing through a variety of soft and hard phase rock formations. Therefore, it is very important to analyze the stability of the surrounding rock in the tailwater tunnel during the construction process.Figure 5
Underground cavern tunnel structure diagram. (a) Shuibuya hydropower station. (b) Underground powerhouse.The right bank diversion underground power station of Shuibuya Water Conservancy Project is located on the right bank of the dam site NE30°. There are 4 power stations installed with a unit capacity of 400 MW and a total installed capacity of 1600 MW. Power station buildings include diversion canals, water inlets, diversion tunnels, main powerhouses, installation sites, busbars, tailwater tunnels, tailwater platforms, tailraces, 500 kV substations, traffic tunnels, ventilation, and pipeline tunnels, and off-site drainage holes.


### 7.2. The Damage Parameters’ Back Analysis

The section of the dam toe plate is selected as the calculated section and calculated according to the plane strain model.The calculation range is centered on the diversion tunnel and is 321.5 m wide and 186.6 m high. TheX-axis is forwardly directed to the side of the mountain, and the Y-axis is oriented vertically upward. The tunnel diameter equals 8.5 m, and tunnel center locates on the level line of y and equals 50 m height. There are 1#, 2#, 3#, and 4# from left to right, respectively. The bottom edge of the calculation domain is fixed, and the normal constraints on both sides. The simplified calculation profile is shown in Figure 6. The model is divided into 459 nodes and 817 units. The initial hydraulic head height is set at y = 80 m, and the left and right sides exert a hydraulic head pressure that varies in a gradient along the gravity direction. The two sides of the model and the perimeter of the tunnel are permeable boundaries. The bottom of the model is the impervious boundary.Figure 6
Model of the tunnel of the finite element method.Next, determine the elastoplastic mechanical parameters of the surrounding rock. Due to a large number of damage constitutive parameters, considering that the damage range of surrounding rock in tunnel excavation is limited, this paper only uses the DP damage model for the argillaceous limestone formation and uses the traditional ideal elastoplastic DP model for the limestone formation. Considering that the surrounding rock is only in the argillaceous limestone formation, only the stratum damage parameters are inverted. According to the results of preliminary surveys and laboratory test, the limestone bulk densityγ = 25 kN/m3, elasticity modulus E = 1.9 GPa, Poisson ratio μ = 0.23, cohesion c = 0.25 MPa, internal frictional angle Φ = 25°, initial porosity of surrounding rock e = 0.038, and the initial permeability kx = ky = 9.71 × 10−3 m/d. The argillaceous limestone bulk density γ = 27 kN/m3, elasticity modulus E = 2.7 GPa, Poisson ratio μ = 0.25, cohesion c = 0.46 MPa, internal frictional angle Φ = 22°, initial porosity of surrounding rock e = 0.033, and the initial permeability kx = ky = 6.48 × 10−3 m/d.Then, substitute the top arch displacement and horizontal convergence displacement of the four tail water holes in the field monitoring and numerical calculation into the objective function (equation (22)), and the damage parameters are used as the optimization variables. Next, carry damage parameters inversion according to the calculated displacement inversion calculation process shown in Figure 4. The iterative curve of CR = 0.9 and F = 0.5∼0.8 is shown in Figure 7.Figure 7
Iterative curves of different variation factors and crossover factors. (a) CR = 0.9 and F = 0.5∼0.8. (b) F = 0.9 and CR = 0.5∼0.8.
(a)(b)Two important parameters in the DE algorithm: mutation factor F and crossover factor CR are studied. The DE/rand/1/bin difference strategy is selected, the mutation factorF = 0.9, the crossover factor CR = 0.5∼0.8, and the crossover factor CR = 0.9. The iterative curve when the variation factor F = 0.5∼0.8 is shown in Figure 7. It can be seen from the above figure that, under the premise of selecting the fixed difference strategy, the variation of the crossover factor CR and the variation factor F have a certain influence on the accuracy and convergence speed of the inversion results, but they can converge to the optimal solution faster.The range of surrounding rock damage parameters is as follows. The cohesion with the damagecr ranges from 0 to 460 MPa, which is less than the range of cohesion c. The parameter ζ ranges from 0 to 1, and the parameter κ is obtained from experiments and is generally 1 to 5000. The displacement inversion parameters calculate results are shown in Table 1.Table 1
Results of back analyzed parameters.
Damage parameterscr/kPaζκInversion results60.1210.684473.167In order to fully explain the calculation function and applicability of the program, the calculation is carried out in two cases: (1) the elastoplastic damage mechanics field calculation is performed separately without considering the seepage effect; (2) the stress-seepage-damage constitutive model established above is used to carry out the calculation of fluid-solid coupling. At the same time, for the actual seawater backflow phenomenon, three water level heights are considered, which are 50 m, 80 m, and 100 m, respectively.


### 7.3. Analysis of Calculation Results considering and Not considering Seepage

The specific calculations and calculation results are as follows:(1)
The elastoplastic damage mechanics field calculation of the surrounding rock and the lining after tunnel excavation is carried out without considering the seepage.The stress cloud diagram after tunnel excavation of the holes without lining and lining is shown in Figures8 and 9, respectively. It can be seen that the x and y direction stresses around the tunnel after the lining is significantly increased compared with the stress excavation of the pores, and stress concentration occurs.Figure 8
X- and Y-direction stress calculation without lining. (a) X-direction stress. (b) Y-direction stress.
(a)(b)Figure 9
X- and Y-direction stress calculation with lining. (a) X-direction stress. (b) Y-direction stress.
(a)(b)Figure10 is a diagram of the plastic zone around the hole after tunnel excavation. It can be found that the range of the plastic zone of the tunnel is significantly reduced due to the influence of the tunnel support.Figure 10
The plastic zone with different support patterns. (a) Without lining after excavation. (b) With lining after excavation.
(a)(b)Figure11 describes the excavation of the pores and the damage cloud pattern. The damaged area of the rock mass caused by excavation is mainly distributed on the left and right sides of the cavern. The damage of the rock mass in the left and right edges of the cavern is particularly serious, and the damage variable value D reaches 0.9, which causes the rock bearing capacity of the corresponding area to decrease significantly.Figure 11
The damage nephogram with different support patterns. (a) Without lining after excavation. (b) With lining after excavation.
(a)(b)This study selects 4 monitoring points above the four tunnel vaults to compare the displacement settlement values of the excavation and lining excavation. The settlement curve of the detection points is shown in Figure12. It can be seen that the surface settlement after the lining is smaller than the surface settlement after the excavation of the pores, and the settlement value of each tunnel is also different. The larger the thickness of the soil above the tunnel, the larger the settlement value of the tunnel.(2)
When calculating the stress-damage-seepage flow of the tunnel underlining, the key factor affecting the coupling between the seepage field and the stress field is the permeability coefficient of the surrounding rock. When the surrounding rock conditions are poor, there are a large number of joint fissures or the pores of the rock and soil which are large, the permeability coefficient of the surrounding rock tends to be relatively large, and the coupling between the seepage field and the stress field will be stronger. Among them, groundwater contributes a lot to the deformation of the overlying stratum of the tunnel. Therefore, if the coupling effect between the seepage field and the stress field is not considered, there would be a large error in the calculation result.Figure 12
The measure points’ settlement curve with different support patterns.Figure13 is a comparison of the amount of surface settlement considered with or without stress-seepage-damage. It can be seen that the surface settlement value considering the stress-seepage-damage coupling effect is larger than the surface without considering the coupling effect, indicating that the seepage flow cannot be ignored for the deformation and failure of the rock mass.Figure 13
The measure points coordinates’ settlement curve for considering coupling and without coupling.


### 7.4. Analysis of Calculation Results of Different Hydraulic Heads

Tunnel excavation will destroy the aquifer structure of the surrounding rock and expose part of the groundwater channel, causing a sharp change in the hydrodynamic condition and the balanced state of the surrounding rock mechanics. Groundwater or other water bodies with which it is hydraulically connected and stored energy are turned from a relatively static state to a flowing state under neutral action. Groundwater flows into the tunnel through the seepage channel and enters the tunnel. The groundwater level decreases, and the pore water pressure decreases correspondingly, forming a reduced area adjacent to the tunnel area. The outside of the tunnel is a slowly changing area, as shown in Figure14(a).Figure 14
The flow vector figures with different support patterns. (a) Without lining after excavation. (b) With lining after excavation.
(a)(b)After tunnel lining construction, because the lining concrete material has strong impermeability, the permeability coefficient is much smaller than the stratum, which can block the drainage of groundwater into the tunnel, and the amount of water in the tunnel decreases, as shown in Figure14(b). The permeability coefficient of the lining is much smaller than that of the surrounding rock, and the water blocking effect is more obvious. Therefore, the amount of water inflow after tunnel lining construction is much smaller.After the tunnel excavation, the pore water pressure of the surrounding rock is continuously consumed, and the groundwater penetrates into the cave, causing the change of the seepage field. Finally, the distribution shape of the seepage field similar to the seepage funnel centering on the tunnel excavation area is formed. Figure15 shows the water pressure distribution of the pores around the hole at different hydraulic head heights. The influence of the lining on the pore water pressure of the surrounding rock is not obvious.Figure 15
The pore water pressure nephogram under different hydraulic heads. (a) Hydraulic head (H)s = 50 m. (b) Hydraulic head (H)s = 80 m. (c) Hydraulic head (H)s = 100 m.
(a)(b)(c)The plastic zone and the damaged zone under different hydraulic heads are shown in Figures16 and 17, respectively. As the hydraulic head Hs increases, the plastic zone and the damage zone of the tunnel gradually increase. In the rock mass stress-seepage-damage coupling model, the cause of rock mass damage evolution is the dynamic change of osmotic water pressure and rock mass stress in the rock mass. This study believes that the rock mass stress is coupled by the dynamic hydraulic gradient, which disturbs the stress distribution of the rock mass and leads to the damage evolution of the rock mass.Figure 16
The plastic zone under different hydraulic heads. (a) Hydraulic head (H)s = 50 m. (b) Hydraulic head (H)s = 80 m. (c) Hydraulic head (H)s = 100 m.
(a)(b)(c)Figure 17
The damage zone under different hydraulic heads. (a) Hydraulic head (H)s = 50 m. (b) Hydraulic head (H)s = 80 m. (c) Hydraulic head (H)s = 100 m.
(a)(b)(c)Based on the elastoplastic stress-seepage-damage constitutive model established above, considering the influence of seepage, the stress field and seepage field of the tunnel can be calculated. The tunnels are numbered 1#, 2#, 3#, and 4# from left to right. In order to consider the change of pore water pressure caused by different thicknesses of the soil, the pore water pressure of the monitoring points above the four tunnels is extracted, as shown in Figure18. It is apparent that due to the location of each tunnel and the thickness of the overlying soil, the pore water pressure of the tunnel is different. Moreover, after the lining is applied, there is a large pore water pressure in the tunnel, so the larger the thickness of the overlying soil, the larger the pore water pressure.Figure 18
Different thickness of the overlying soil pore water pressure.Figure19 shows the distribution curve of the surrounding rock permeability coefficient along with the horizontal and vertical directions of the tunnel surrounding rock after tunnel excavation. Then, Figure 19(a) is a schematic diagram of the tunnel monitoring point. Based on the measuring points in Figure 19(a), the relationship between the permeability coefficient of the surrounding rock and the surrounding rock distance can be analyzed, as shown in Figures 19(b) and 19(c). It can be seen that the stress redistribution caused by tunnel excavation has a nonnegligible influence on the permeability coefficient of the surrounding rock of the tunnel. From the permeability coefficient of line AB and line IJ in Figure 19(b), it can be seen that the farther the tunnel is, the smaller the horizontal permeability coefficient is. From the survey lines CD, EF, and GH, it can be seen that the permeability coefficient of the surrounding rock of the tunnel increases first and then decreases. Moreover, the permeability coefficient of the vertical position of the tunnel gradually decreases with increasing distance from the circumference of the tunnel (Figure 19(c)).Figure 19
Permeability coefficient of horizontal and vertical direction in surrounding rock. (a) Sketch map. (b) Horizontal direction. (c) Vertical direction.
(a)(b)(c)Therefore, the closer the surrounding rock is, the larger the permeability coefficient is. As the depth of burial increases, the permeability coefficient around the tunnel also increases. The permeability coefficient of surrounding rock is obviously affected by the damage of rock mass, which indicates that the mechanical-hydraulic-damage (MHD) coupling model used in this model can not only reflect the evolution of damage but also reflect the relationship of the damage-permeability coefficient. So, it is more reasonable than the mechanical-hydraulic (MH) model of the constant permeability coefficient.In summary, the excavation of the rock mass causes serious damage to the left and right edge areas of the cavern. The bearing capacity of the rock mass decreases, the permeability coefficient increases rapidly, and the seepage flow increases. Therefore, in the case of high water pressure, these areas are most prone to water inrush, so the corresponding reinforcement measures should be adopted. In the actual project, the variation of the stress field of the rock mass and the distribution of the damage zone can be used to delineate the change zone of permeability, which provides a basis for the antiseepage design. On the contrary, when the permeability of the rock mass changes significantly, the structure must be destroyed, and the damage extent and failure mode of the rock mass can be determined.


## 7.1. General State of Shuibuya Hydropower Station

The Shuibuya Water Conservancy Project is located in the middle section of Qingjiang River in Badong County, Hubei Province, and is one of the important power stations for the development of the Qingjiang River Basin, as shown in Figure5. The geological engineering conditions of the surrounding rock of the underground cavern are very complicated, and the tailwater tunnel section is particularly serious, passing through a variety of soft and hard phase rock formations. Therefore, it is very important to analyze the stability of the surrounding rock in the tailwater tunnel during the construction process.Figure 5
Underground cavern tunnel structure diagram. (a) Shuibuya hydropower station. (b) Underground powerhouse.The right bank diversion underground power station of Shuibuya Water Conservancy Project is located on the right bank of the dam site NE30°. There are 4 power stations installed with a unit capacity of 400 MW and a total installed capacity of 1600 MW. Power station buildings include diversion canals, water inlets, diversion tunnels, main powerhouses, installation sites, busbars, tailwater tunnels, tailwater platforms, tailraces, 500 kV substations, traffic tunnels, ventilation, and pipeline tunnels, and off-site drainage holes.


## 7.2. The Damage Parameters’ Back Analysis

The section of the dam toe plate is selected as the calculated section and calculated according to the plane strain model.The calculation range is centered on the diversion tunnel and is 321.5 m wide and 186.6 m high. TheX-axis is forwardly directed to the side of the mountain, and the Y-axis is oriented vertically upward. The tunnel diameter equals 8.5 m, and tunnel center locates on the level line of y and equals 50 m height. There are 1#, 2#, 3#, and 4# from left to right, respectively. The bottom edge of the calculation domain is fixed, and the normal constraints on both sides. The simplified calculation profile is shown in Figure 6. The model is divided into 459 nodes and 817 units. The initial hydraulic head height is set at y = 80 m, and the left and right sides exert a hydraulic head pressure that varies in a gradient along the gravity direction. The two sides of the model and the perimeter of the tunnel are permeable boundaries. The bottom of the model is the impervious boundary.Figure 6
Model of the tunnel of the finite element method.Next, determine the elastoplastic mechanical parameters of the surrounding rock. Due to a large number of damage constitutive parameters, considering that the damage range of surrounding rock in tunnel excavation is limited, this paper only uses the DP damage model for the argillaceous limestone formation and uses the traditional ideal elastoplastic DP model for the limestone formation. Considering that the surrounding rock is only in the argillaceous limestone formation, only the stratum damage parameters are inverted. According to the results of preliminary surveys and laboratory test, the limestone bulk densityγ = 25 kN/m3, elasticity modulus E = 1.9 GPa, Poisson ratio μ = 0.23, cohesion c = 0.25 MPa, internal frictional angle Φ = 25°, initial porosity of surrounding rock e = 0.038, and the initial permeability kx = ky = 9.71 × 10−3 m/d. The argillaceous limestone bulk density γ = 27 kN/m3, elasticity modulus E = 2.7 GPa, Poisson ratio μ = 0.25, cohesion c = 0.46 MPa, internal frictional angle Φ = 22°, initial porosity of surrounding rock e = 0.033, and the initial permeability kx = ky = 6.48 × 10−3 m/d.Then, substitute the top arch displacement and horizontal convergence displacement of the four tail water holes in the field monitoring and numerical calculation into the objective function (equation (22)), and the damage parameters are used as the optimization variables. Next, carry damage parameters inversion according to the calculated displacement inversion calculation process shown in Figure 4. The iterative curve of CR = 0.9 and F = 0.5∼0.8 is shown in Figure 7.Figure 7
Iterative curves of different variation factors and crossover factors. (a) CR = 0.9 and F = 0.5∼0.8. (b) F = 0.9 and CR = 0.5∼0.8.
(a)(b)Two important parameters in the DE algorithm: mutation factor F and crossover factor CR are studied. The DE/rand/1/bin difference strategy is selected, the mutation factorF = 0.9, the crossover factor CR = 0.5∼0.8, and the crossover factor CR = 0.9. The iterative curve when the variation factor F = 0.5∼0.8 is shown in Figure 7. It can be seen from the above figure that, under the premise of selecting the fixed difference strategy, the variation of the crossover factor CR and the variation factor F have a certain influence on the accuracy and convergence speed of the inversion results, but they can converge to the optimal solution faster.The range of surrounding rock damage parameters is as follows. The cohesion with the damagecr ranges from 0 to 460 MPa, which is less than the range of cohesion c. The parameter ζ ranges from 0 to 1, and the parameter κ is obtained from experiments and is generally 1 to 5000. The displacement inversion parameters calculate results are shown in Table 1.Table 1
Results of back analyzed parameters.
Damage parameterscr/kPaζκInversion results60.1210.684473.167In order to fully explain the calculation function and applicability of the program, the calculation is carried out in two cases: (1) the elastoplastic damage mechanics field calculation is performed separately without considering the seepage effect; (2) the stress-seepage-damage constitutive model established above is used to carry out the calculation of fluid-solid coupling. At the same time, for the actual seawater backflow phenomenon, three water level heights are considered, which are 50 m, 80 m, and 100 m, respectively.


## 7.3. Analysis of Calculation Results considering and Not considering Seepage

The specific calculations and calculation results are as follows:(1)
The elastoplastic damage mechanics field calculation of the surrounding rock and the lining after tunnel excavation is carried out without considering the seepage.The stress cloud diagram after tunnel excavation of the holes without lining and lining is shown in Figures8 and 9, respectively. It can be seen that the x and y direction stresses around the tunnel after the lining is significantly increased compared with the stress excavation of the pores, and stress concentration occurs.Figure 8
X- and Y-direction stress calculation without lining. (a) X-direction stress. (b) Y-direction stress.
(a)(b)Figure 9
X- and Y-direction stress calculation with lining. (a) X-direction stress. (b) Y-direction stress.
(a)(b)Figure10 is a diagram of the plastic zone around the hole after tunnel excavation. It can be found that the range of the plastic zone of the tunnel is significantly reduced due to the influence of the tunnel support.Figure 10
The plastic zone with different support patterns. (a) Without lining after excavation. (b) With lining after excavation.
(a)(b)Figure11 describes the excavation of the pores and the damage cloud pattern. The damaged area of the rock mass caused by excavation is mainly distributed on the left and right sides of the cavern. The damage of the rock mass in the left and right edges of the cavern is particularly serious, and the damage variable value D reaches 0.9, which causes the rock bearing capacity of the corresponding area to decrease significantly.Figure 11
The damage nephogram with different support patterns. (a) Without lining after excavation. (b) With lining after excavation.
(a)(b)This study selects 4 monitoring points above the four tunnel vaults to compare the displacement settlement values of the excavation and lining excavation. The settlement curve of the detection points is shown in Figure12. It can be seen that the surface settlement after the lining is smaller than the surface settlement after the excavation of the pores, and the settlement value of each tunnel is also different. The larger the thickness of the soil above the tunnel, the larger the settlement value of the tunnel.(2)
When calculating the stress-damage-seepage flow of the tunnel underlining, the key factor affecting the coupling between the seepage field and the stress field is the permeability coefficient of the surrounding rock. When the surrounding rock conditions are poor, there are a large number of joint fissures or the pores of the rock and soil which are large, the permeability coefficient of the surrounding rock tends to be relatively large, and the coupling between the seepage field and the stress field will be stronger. Among them, groundwater contributes a lot to the deformation of the overlying stratum of the tunnel. Therefore, if the coupling effect between the seepage field and the stress field is not considered, there would be a large error in the calculation result.Figure 12
The measure points’ settlement curve with different support patterns.Figure13 is a comparison of the amount of surface settlement considered with or without stress-seepage-damage. It can be seen that the surface settlement value considering the stress-seepage-damage coupling effect is larger than the surface without considering the coupling effect, indicating that the seepage flow cannot be ignored for the deformation and failure of the rock mass.Figure 13
The measure points coordinates’ settlement curve for considering coupling and without coupling.


## 7.4. Analysis of Calculation Results of Different Hydraulic Heads

Tunnel excavation will destroy the aquifer structure of the surrounding rock and expose part of the groundwater channel, causing a sharp change in the hydrodynamic condition and the balanced state of the surrounding rock mechanics. Groundwater or other water bodies with which it is hydraulically connected and stored energy are turned from a relatively static state to a flowing state under neutral action. Groundwater flows into the tunnel through the seepage channel and enters the tunnel. The groundwater level decreases, and the pore water pressure decreases correspondingly, forming a reduced area adjacent to the tunnel area. The outside of the tunnel is a slowly changing area, as shown in Figure14(a).Figure 14
The flow vector figures with different support patterns. (a) Without lining after excavation. (b) With lining after excavation.
(a)(b)After tunnel lining construction, because the lining concrete material has strong impermeability, the permeability coefficient is much smaller than the stratum, which can block the drainage of groundwater into the tunnel, and the amount of water in the tunnel decreases, as shown in Figure14(b). The permeability coefficient of the lining is much smaller than that of the surrounding rock, and the water blocking effect is more obvious. Therefore, the amount of water inflow after tunnel lining construction is much smaller.After the tunnel excavation, the pore water pressure of the surrounding rock is continuously consumed, and the groundwater penetrates into the cave, causing the change of the seepage field. Finally, the distribution shape of the seepage field similar to the seepage funnel centering on the tunnel excavation area is formed. Figure15 shows the water pressure distribution of the pores around the hole at different hydraulic head heights. The influence of the lining on the pore water pressure of the surrounding rock is not obvious.Figure 15
The pore water pressure nephogram under different hydraulic heads. (a) Hydraulic head (H)s = 50 m. (b) Hydraulic head (H)s = 80 m. (c) Hydraulic head (H)s = 100 m.
(a)(b)(c)The plastic zone and the damaged zone under different hydraulic heads are shown in Figures16 and 17, respectively. As the hydraulic head Hs increases, the plastic zone and the damage zone of the tunnel gradually increase. In the rock mass stress-seepage-damage coupling model, the cause of rock mass damage evolution is the dynamic change of osmotic water pressure and rock mass stress in the rock mass. This study believes that the rock mass stress is coupled by the dynamic hydraulic gradient, which disturbs the stress distribution of the rock mass and leads to the damage evolution of the rock mass.Figure 16
The plastic zone under different hydraulic heads. (a) Hydraulic head (H)s = 50 m. (b) Hydraulic head (H)s = 80 m. (c) Hydraulic head (H)s = 100 m.
(a)(b)(c)Figure 17
The damage zone under different hydraulic heads. (a) Hydraulic head (H)s = 50 m. (b) Hydraulic head (H)s = 80 m. (c) Hydraulic head (H)s = 100 m.
(a)(b)(c)Based on the elastoplastic stress-seepage-damage constitutive model established above, considering the influence of seepage, the stress field and seepage field of the tunnel can be calculated. The tunnels are numbered 1#, 2#, 3#, and 4# from left to right. In order to consider the change of pore water pressure caused by different thicknesses of the soil, the pore water pressure of the monitoring points above the four tunnels is extracted, as shown in Figure18. It is apparent that due to the location of each tunnel and the thickness of the overlying soil, the pore water pressure of the tunnel is different. Moreover, after the lining is applied, there is a large pore water pressure in the tunnel, so the larger the thickness of the overlying soil, the larger the pore water pressure.Figure 18
Different thickness of the overlying soil pore water pressure.Figure19 shows the distribution curve of the surrounding rock permeability coefficient along with the horizontal and vertical directions of the tunnel surrounding rock after tunnel excavation. Then, Figure 19(a) is a schematic diagram of the tunnel monitoring point. Based on the measuring points in Figure 19(a), the relationship between the permeability coefficient of the surrounding rock and the surrounding rock distance can be analyzed, as shown in Figures 19(b) and 19(c). It can be seen that the stress redistribution caused by tunnel excavation has a nonnegligible influence on the permeability coefficient of the surrounding rock of the tunnel. From the permeability coefficient of line AB and line IJ in Figure 19(b), it can be seen that the farther the tunnel is, the smaller the horizontal permeability coefficient is. From the survey lines CD, EF, and GH, it can be seen that the permeability coefficient of the surrounding rock of the tunnel increases first and then decreases. Moreover, the permeability coefficient of the vertical position of the tunnel gradually decreases with increasing distance from the circumference of the tunnel (Figure 19(c)).Figure 19
Permeability coefficient of horizontal and vertical direction in surrounding rock. (a) Sketch map. (b) Horizontal direction. (c) Vertical direction.
(a)(b)(c)Therefore, the closer the surrounding rock is, the larger the permeability coefficient is. As the depth of burial increases, the permeability coefficient around the tunnel also increases. The permeability coefficient of surrounding rock is obviously affected by the damage of rock mass, which indicates that the mechanical-hydraulic-damage (MHD) coupling model used in this model can not only reflect the evolution of damage but also reflect the relationship of the damage-permeability coefficient. So, it is more reasonable than the mechanical-hydraulic (MH) model of the constant permeability coefficient.In summary, the excavation of the rock mass causes serious damage to the left and right edge areas of the cavern. The bearing capacity of the rock mass decreases, the permeability coefficient increases rapidly, and the seepage flow increases. Therefore, in the case of high water pressure, these areas are most prone to water inrush, so the corresponding reinforcement measures should be adopted. In the actual project, the variation of the stress field of the rock mass and the distribution of the damage zone can be used to delineate the change zone of permeability, which provides a basis for the antiseepage design. On the contrary, when the permeability of the rock mass changes significantly, the structure must be destroyed, and the damage extent and failure mode of the rock mass can be determined.


## 8. Conclusions

Through the research in this paper, the following conclusions are obtained:(1)
This paper establishes an elastoplastic damage constitutive model based on the Drucker–Prager yield criterion and uses equivalent plastic strain to characterize the evolution of rock damage variables. According to the dynamic evolution formula of permeability coefficient, when the rock is in the elastoplastic state, the elastoplastic stress-seepage-damage model of rock is established. The numerical solution of the coupled model is given, and the corresponding program is compiled by C++ language. The established elastoplastic stress-seepage-damage coupling model of rock is applied to the tunnel simulation. The results show that the existence of a seepage field exacerbates the damage and failure of the surrounding rock, and the damage evolution of the stress field in turn affects the change of permeability coefficient. There is a certain gap between the results of the coupled model of elastoplastic damage model without considering seepage and considering seepage. Therefore, it is necessary to consider the deformation, damage, and failure characteristics of the seepage in the tunnel stability analysis.(2)
The engineering application shows that the coupled model can truly reflect the complex macroscopic failure of rock materials through the interaction of stress, seepage, and damage. The programmed program can simulate the coupling characteristics of groundwater seepage field, stress field, and damage field and provide an analysis method for engineering construction with serious impact on groundwater. The algorithm of this paper has advantages over the elastic-plastic damage of the simple stress field or the ideal elastic-plastic fluid-solid coupling algorithm.(3)
Based on the principle of differential evolution algorithm, an intelligent back analysis method is established for the problems involving many parameters and is difficult to measure in the coupled model. The return mapping algorithm has the characteristics of accuracy and stability. The Newton–Raphson method is used in the iteration to obtain the convergence rate of the approximate square. Besides, the differential evolution algorithm has no requirement for the initial value, the search success rate is high, the convergence speed is fast, and no encoding and decoding operations are needed, which is convenient for practical application. Combining the advantages of these two algorithms, the damage parameters in the coupled model are inverted. The inversion results show that the program has good accuracy and stability.(4)
The algorithm in this paper is based on the continuous mechanical medium model and the linear strength criterion, DP criterion. The algorithm of this paper is not suitable for the fractured jointed rock mass or the fractured jointed rock mass that conforms to the nonlinear strength criterion. For jointed rock mass engineering, it is necessary to explore the Hoek–Brown constitutive damage model or discrete element numerical algorithm based on nonlinear strength criterion. This is the direction that needs further research in the future.



---

*Source: 8341528-2021-02-24.xml*