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

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 stressseepage-damage coupling program and the damage parameter displacement back analysis program is compiled using C++ language. -en, the program is used to calculate the coupling problem of tunnel elastoplastic stress-seepage-damage. -e results show that the proposed elastoplastic damage constitutive model can well describe the mechanical behavior of rock. -e computational procedure can also simulate practical engineering problems, which can provide specific guidance for site construction.


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. e 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. e 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]. e 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].
ere is a lot of research on the seepage characteristics of tunnels surrounding rocks and underground caverns during construction [3][4][5][6][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. erefore, 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 fluidsolid 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][16][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) ere 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) e 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.
erefore, 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. e global optimization algorithm of difference evolution is introduced to back-analyze the damage parameters. e step iteration, numerical solving methods, and constitutive integration algorithm are studied. e elastoplastic stress-seepage-damage coupling program and intelligent back analysis program are compiled by C++ language. e programs are used to calculate the diversion tunnel of Shuibuya Hydropower Station. e 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.

Rock Permeability.
e 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. e 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.

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.
is interaction is defined as seepage-stress coupling.

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 seepagestress coupling, mainly as follows: the impact of damage on the seepage process, the weakening of water, and the damage induced by osmotic stress. is problem is the seepage-damage coupling problem of rock during the rupture process. e seepage-stress coupling study focuses on the coupling methods for establishing different pore structure systems and describes their applicable conditions. e 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.

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. en, a coupled model describing rock stress-seepage-damage is established. e model can effectively solve the stressseepage-damage coupling problem of the tunnel surrounding rock. e assumptions in this paper are as follows. (1) e skeleton of the saturated body studied is an ideal elastoplastic damage isotropic body and satisfies the assumption of small deformation. (2) e seepage is laminar, following Darcy's law. (3) e seepage fluid is incompressible, and the effects of temperature changes are negligible.

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: e boundary value problem is shown in Figure 1. Equation (1) can be composed of the following three parts: 2 Advances in Civil Engineering where f j is body force, σ ij is the total stress tensor, Ω is the problem-solving domain, n i is the normal border cosine, t j is the surface forces acting on the known boundary surface, Γ 1 is the known force boundary, u m is the known displacement on the border, and Γ 2 is the known displacement boundary. e relationship between stress and strain is different for different constitutive models. e constitutive equation is abbreviated as follows: e geometric equation is as follows: where ε is the strain tensor and u is the displacement vector. e boundary value problem is transformed into a problem of solving the displacement u 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. e principle of effective stress in porous media is as follows: 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.

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: Bring h � (p/c) + z into equation (6); then, where h is the infiltration head, x, y, and z are the space coordinate, t is the time coordinate, k x , k y , and k z is the permeability coefficient of the x, y, and z axes as the main axis direction, S s is the unit storage amount, c 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: e boundary conditions, i.e., the head boundary and the flow boundary, are as follows: where Μ 1 is the known head boundary, f 1 is the known head boundary value, Μ 2 is the known flow boundary, and f 2 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. e equilibrium differential equation under the action of the seepage field embodies the dynamic coupling effect of stress seepage.

Damage Variable and the Damage Evolution Equation.
is 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. e calculation of equivalent plastic strain can be shown as follows: where ε p1 , ε p2 , and ε p3 are three principal plastic strain, respectively. e evolution equation of the corresponding damage variable D is as follows: where the equivalent plastic strain threshold ε p 0 � 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.

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). is 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. e 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. e 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. erefore, 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. e 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]: where n 0 is the initial porosity, K 0 is the initial permeability of the rock, and ε v is the volumetric strain. e permeability coefficient of the rock in the plastic phase is shown as follows [21]: where K M and K D are permeability coefficients of the undamaged rock and fractured rock, respectively, and ε pF v � Dε p v , which is the plastic volumetric strain for the defect phase.

Numerical Solving and Step Iteration
Method for Coupling Model e solution to the elastic-plastic-damage-seepage coupling model of rock is a complex nonlinear problem. e 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.

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: where [K m ] is the total stiffness matrix of the stress field, u { } is the node displacement column vector, Δu is increment displacement, and {f m } is a nodal force vector that includes physical strength, surface strength, and pore pressure of the equivalent load. e overall stiffness matrix K is assembled with the element stiffness matrix. K (e) T is the function of consistent tangent modulus C as follows: where B is the strain matrix, that is, the matrix for finding strain based on displacement and B T 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].
e Drucker-Prager model is widely used in rock materials to describe the characteristics of plastic stressdeformation.
e plastic yield function considering the damaging effect is as follows: (16) 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), J 2 is the second invariant of deviatoric stress, J 2 � 1/2s: s, s(σ) is the deviatoric tensor of stress, s(σ) � σ − p(σ)I, I is the second-order tensor, I � δ ij e i ⊗ e j , 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. e influence of damage on the internal friction angle is very small, so only the effect of damage on cohesion c 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]: where c(ε p ) ′ is the cohesion, c r is the cohesion of rock when it is obviously damaged, and ζ is the material parameter with a value between 0 and 1.

Constitutive Integration Algorithm.
e 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. e solution of the stress is related to the selection of the constitutive model. e 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. is 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 σ trial n+1 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.

Consistent Tangent Modulus.
When it returns to the smooth conical surface, the consistent tangent modulus is derived as follows: where C ep d is the consistent tangent modulus, G(D n ) and K(D n ) are damage shear modulus and damage bulk modulus, respectively, I d is volume tensor, T is the second-order unit tensor parallel to the elastic predicted strain, a, b, c, and d are coefficients, ε et dn+1 is the strain during elastic prediction at the time of t n+1 , D n is the damage variable at the time of t n , and η and ζ are the material parameters: When it returns to a sharp point, the consistent tangent modulus is derived as follows: (20) where H is the hardening modulus. e 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.

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 where [K s ] is seepage matrix, h { } is hydraulic head column vector, and f s is a free term column vector.
us, 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. e preprocedure process divides the mesh by means of ANSYS and then converts it into a seepage finite element program to solve the input file. e 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.

Elastic-Plastic Stress-Seepage-Damage Coupling Iteration.
e coupling mechanism between the groundwater seepage field and stress field in the rock mass is a relatively complex dynamic process. e 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.
ere is much research related to the numerical methods of rock stress-seepage coupling which has given different coupling methods. ere are two main methods, i.e., stepby-step iteration and one-time coupling. In this paper, the Advances in Civil Engineering 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 Section 5.1. e volumetric strain of the element is obtained from the calculation of the deformation field.
en, 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 Section 5.2, and it is brought back to the calculation of the mechanical field of Section 5.1 by the principle of effective stress. e 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. e elastoplastic stress-seepage-damage coupling procedure is shown in Figure 3.

Damage Parameters' Inversion Based on
Differential Evolution Algorithm e damage parameters inverse problem turned into the optimization problem of constraint [24]: e constraint condition: where Y 0 i is the observed displacement, Y i is the FEM calculated value by the seepage-stress-damage coupling model, m is the number of observations, x i is the damage parameter, and x l i and x u i are the upper bound and lower bound on x i . e 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. e specific principle is as follows. where G means the group evolution of every generation and i means the individual location in the population. e initial population is produced according to the method of random that uniform distribution is in the solution space: where r and ∈[0, 1] is a random number, x i,1 is the solution vector of first-generation, j � 1, 2, . . ., R, R is the number of the dimensions of solution vector, and x j U and x j L are the upper and lower bounds of the jth component, respectively.

Mutation Operation.
e mutation operation uses a different strategy, which uses the difference vector between individuals in the population to perturb the individual to achieve individual variation. e 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 vector x i, G in the Gth generation, each vector individual contains R components, and the solution of the variance vector ] i, G + 1 is   Advances in Civil Engineering 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.

Interlace Operation.
Interlace operation is in order to increase the diversity of the population. e new test vector u i,G+1 can be calculated by hybridizing the target vector x i, G, and the variation vector ] i, G + 1 according to the following rules: G+1 , if (randb(j) ≤ CR‖j � rnbr(i)), where randb(j) ∈ [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 rnbr(i) is the coefficient corresponding with the ith vector, random integers among 1, 2, . . ., R.

Selecting
Operation. e greedy search method is used to select whether to select the test vector u i,G+1 as the target vector of the (G + 1)th generation.
Compare test vector u i,G+1 with target vector x i, G , and choose vector u i,G+1 if u i,G+1 corresponding the smaller target value. Otherwise, retain x i,G if x i,G is corresponding the smaller one. New population x i,G+1 can be generated after selection as follows: e DE algorithm adopts real-coded and regards equation (22) as a fitness function in this paper. e back analysis calculation process is shown in Figure 4. In this study, the self-developed elastoplastic stressseepage-damage finite element solution program is named RMAST in the intelligent displacement back analysis program, which can calculate the displacement Y i of equation (22) and complete the inversion of damage parameters.

General State of Shuibuya Hydropower Station.
e 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 Figure 5. e 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. erefore, it is very important to analyze the stability of the surrounding rock in the tailwater tunnel during the construction process. e right bank diversion underground power station of Shuibuya Water Conservancy Project is located on the right bank of the dam site NE30°.
ere 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.

e Damage Parameters' Back Analysis.
e section of the dam toe plate is selected as the calculated section and calculated according to the plane strain model. Advances in Civil Engineering e calculation range is centered on the diversion tunnel and is 321.5 m wide and 186.6 m high. e X-axis is forwardly directed to the side of the mountain, and the Y-axis is oriented vertically upward. e tunnel diameter equals 8.5 m, and tunnel center locates on the level line of y and equals 50 m height. ere are 1#, 2#, 3#, and 4# from left to right, respectively. e bottom edge of the calculation domain is fixed, and the normal constraints on both sides. e simplified calculation profile is shown in Figure 6. e model is divided into 459 nodes and 817 units. e 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. e two sides of the model and the perimeter of the tunnel are permeable boundaries. e bottom of the model is the impervious boundary.
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 c � 25 kN/m 3 , 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 k x � k y � 9.71 × 10 − 3 m/d. e argillaceous limestone bulk density c � 27 kN/m 3 , 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 k x � k y � 6.48 × 10 − 3 m/d. en, substitute the top arch displacement and horizontal convergence displacement of the four tail water holes in the field monitoring and numerical calculation into the  Figure 7.
Two important parameters in the DE algorithm: mutation factor F and crossover factor CR are studied. e DE/ rand/1/bin difference strategy is selected, the mutation factor F � 0.9, the crossover factor CR � 0.5∼0.8, and the crossover factor CR � 0.9. e 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. e range of surrounding rock damage parameters is as follows. e cohesion with the damage c r ranges from 0 to 460 MPa, which is less than the range of cohesion c. e parameter ζ ranges from 0 to 1, and the parameter κ is obtained from experiments and is generally 1 to 5000. e displacement inversion parameters calculate results are shown in Table 1.
In 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.

Analysis of Calculation Results considering and Not considering Seepage.
e specific calculations and calculation results are as follows: (1) e elastoplastic damage mechanics field calculation of the surrounding rock and the lining after tunnel excavation is carried out without considering the seepage.
e stress cloud diagram after tunnel excavation of the holes without lining and lining is shown in Figures 8 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 10 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 11 describes the excavation of the pores and the damage cloud pattern. e damaged area of the rock mass caused by excavation is mainly distributed on the left and right sides of the cavern. e 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.
is study selects 4 monitoring points above the four tunnel vaults to compare the displacement settlement values of the excavation and lining excavation. e settlement curve of the detection points is shown in Figure 12. 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. e 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 Advances in Civil Engineering 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. erefore, 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 13 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.

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. e groundwater level decreases, and the pore water pressure decreases correspondingly, forming a reduced area adjacent to the tunnel area. e outside of the tunnel is a slowly changing area, as shown in Figure 14(a).
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 Figure 14(b). e permeability coefficient of the lining is much smaller than that of the surrounding rock, and the water blocking effect is more obvious. erefore, 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. Figure 15 shows the water pressure distribution of the pores around the hole at different hydraulic head heights. e influence of the lining on the pore water pressure of the surrounding rock is not obvious. e plastic zone and the damaged zone under different hydraulic heads are shown in Figures 16 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.
is 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.
Based on the elastoplastic stress-seepage-damage constitutive model established above, considering the influence   Advances in Civil Engineering of seepage, the stress field and seepage field of the tunnel can be calculated. e 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 Figure 18. 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 19 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. en, 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)    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)). erefore, 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. e 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. e bearing capacity of the rock mass decreases, the permeability coefficient increases rapidly, and the seepage flow increases. erefore, 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.

Conclusions
rough the research in this paper, the following conclusions are obtained: (1) is 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-seepagedamage model of rock is established. e numerical solution of the coupled model is given, and the corresponding program is compiled by C++ language. e established elastoplastic stress-seepagedamage coupling model of rock is applied to the tunnel simulation. e 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. ere is a certain gap between the results of the coupled model of elastoplastic damage model without considering seepage and considering seepage. erefore, it is necessary to consider the deformation, damage, and failure characteristics of the seepage in the tunnel stability analysis.
(2) e 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. e 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. e 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. e return mapping algorithm has the characteristics of accuracy and stability. e 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. e inversion results show that the program has good accuracy and stability. (4) e algorithm in this paper is based on the continuous mechanical medium model and the linear strength criterion, DP criterion. e 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  16 Advances in Civil Engineering or discrete element numerical algorithm based on nonlinear strength criterion. is is the direction that needs further research in the future.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.