Mechanism of Splitting Failure for High Sidewall Cavern of Hydropower Station Based on Complex Function and Strain Gradient

Abstract

With the rapid development of society, the number of hydropower projects has increased. During the construction of these projects, due to excavation-induced unloading, the high sidewalls of the hydropower station are often subject to splitting failure, which produces many adverse effects on the construction of the cavern. In order to reveal the formation mechanism of splitting failure of hydropower stations, based on the strain gradient theory and elasto-plastic damage theory, we proposed an elasto-plastic damage softening model. Using the ODE45 program in MATLAB, we solved the numerical solution of displacement and stress of circular cavern based on our proposed elastoplastic damage model. Then, we apply the complex function method and use the Schwarz–Christoffel integral formula to obtain the mapping function from the outer domain of the high sidewall cavern to the outer domain of the unit circle. Finally, the elastic-plastic region and displacement distribution of the high sidewall cavern are obtained by mapping the obtained elastic-plastic solution of the circular cavern under the axisymmetric condition. In future research, it is necessary to further study the corresponding relationship between the internal length parameter of the material and its internal microstructure in order to accurately determine the internal length parameter.

1. Introduction and Background

With the rapid development of society, people’s demand for electric energy is increasing, so there are a greater number of hydropower projects. However, with an increase in hydropower projects, it was discovered that some nonlinear deformations and failures occurred in the surrounding rock of caverns during the construction of these projects. These deformations and damages are a serious threat to the construction safety and long-term stability of underground caverns. Among them, the deformation failure mode of the splitting failure of the hydropower projects with a high sidewall is quite different from that of the shallow cavern, showing obvious discontinuity and nonlinear behavior. This splitting failure has been widely recorded and described by many researchers [1,2,3,4], as shown in Figure 1. However, there is no clear understanding of the formation mechanism of splitting failure of high sidewall cavern, and the traditional continuum mechanics theory cannot provide a satisfactory explanation.

Because the failure process of the rock is very complicated, it is difficult to obtain the ideal result if only the classical continuity mechanics method is used to describe it. Therefore, it is necessary to introduce the damage mechanics method in the study of splitting failure of high-side caverns of hydropower stations. Damage mechanics can well simulate the nonlinear behavior of rocks after the peak. In addition, it shows good applicability in simulating the strain-softening characteristics of rocks. Li et al. [5] applied elastic damage mechanics to model the softening behavior of rock masses. Then, they proposed a numerical method based on the maximum tensile stress criterion and strain energy density theory to simulate the zonal fracture phenomenon. Tao et al. [6] applied a damage formulation to simulate softening and modulus reduction to propose a new method to understand the mechanisms and behaviors of high initial stress rocks. Jia et al. [7] studied the damage evolution around deep underground openings under multiaxial stress through 3D numerical tests. The results showed that biaxial stress parallel to the free surface constrained by zero or low confinement stress contributes to spalling or parallel fracture of the slab.

In the compression test of rocks, as the load increases, the internal stress of the rock will rapidly decrease to a lower level after reaching its peak strength, showing obvious post-peak strain softening. At the same time, along with the decline in the bearing capacity of the rock, the deformation and failure of the rocks are also concentrated in one or several narrow-banded areas, which is the phenomenon of localization of the rock strain. In the experimental research on the strain localization of rocks, many scholars [8,9,10,11,12] have given strong evidence that there is an obvious strain gradient effect in the deformation localization zone, and the strain gradient effect plays a controlling role in the localized deformation and failure process of the rock. Fleck et al. [13,14,15,16,17,18,19] carried out a large number of studies related to strain gradient. Based on the concepts of statistical storage and geometric asymmetric dislocation, they introduced plasticity into the strain gradient theory and established a phenomenological strain gradient plasticity theory. Menzel et al. [20] Established geometric linear formulas for high-order gradient plasticity of single crystal and polycrystalline materials based on dislocation and incompatible continuum theory. They found that due to dislocations or incompatible stresses, higher-order strain gradients are introduced into the yield conditions, which is similar to the structure of motion hardening in form, so the back stress follows the law of nonlocal evolution. Therefore, if the strain gradient effect is ignored in the post-peak localized deformation and failure analysis of rocks, it will cause serious deviations between the analysis results and the real situation. Therefore, it is very important to introduce strain gradient theory into the study of splitting failure of hydropower station caverns.

The common tunnel section shape of an underground powerhouse of a hydropower station is not circular, but high sidewall. For non-circular caverns, the commonly used method is to equivalent them to circles with the same area for analysis, but results obtained in this way will produce large errors. In order to reduce the error, some experts and scholars introduce the complex function theory into the study of solving the analytical solutions of stress and displacement of underground non-circular caverns. In order to reduce the error, some scholars [21,22,23,24] have introduced the theory of complex variable function to solve the analytical solution of the stress and displacement of underground non-circular caverns. Deng et al. [25] used the complex variable function solution method to derive the analytical solution of the stress and displacement of the surrounding rock of a highway tunnel of arbitrary shape under the action of supporting pressure and remote site stress, and solved the numerical solution by programming on the MATLAB platform. Lu et al. [26] used the conformal mapping method of the complex function theory to determine the optimal shape of the closed concrete support. Yang et al. [27] used complex variable function theory and a multipolar coordinate system to study the scattering problem of lining circular tunnels in an anisotropic half-space, and gave typical numerical simulation results of surface displacement amplitude under different anisotropic parameters.

Although great progress has been made in studying the microscopic mechanisms of rock specimens [28,29,30,31,32], laboratory or field experiments [33,34,35,36,37,38,39], and numerical simulations [40,41,42,43], so far, the theoretical method of splitting failure of surrounding rock of high side wall caverns under high ground stress is not perfect, and there is no theoretical explanation widely accepted by everyone. It is urgent to establish a more scientific mechanical model to reveal the formation mechanism of splitting failure.

2. Objectives

In this paper, we propose an elasto-plastic damage softening model based on strain gradient theory and elasto-plastic damage theory. Based on this elasto-plastic damage model, the numerical solutions for the displacements and stresses in the circular cavern were solved using the ODE45 program in MATLAB. Then, we applied the complex function method and use the Schwarz–Christoffel integral formula to obtain the mapping function from the outer domain of the high sidewall cavern to the outer domain of the unit circle. Next, the elasto-plastic solution conformal mapping of the circular cavern under axisymmetric conditions is used to obtain the elasto-plastic area and displacement distribution of the cavern with a high sidewall. Taking Pubugou Hydropower Station (Figure 2) as the engineering background, during the construction of the main powerhouse, several splitting cracks appeared on the sidewall, which caused the cracking of the concrete shotcrete on the sidewall. The maximum crack width was 20 mm. The numerical analysis results are basically consistent with the on-site monitoring and model test results. Finally, the variation of tangential stress of high side wall caverns with different height span ratios is obtained by comparison, which provides a theoretical basis for the study of the formation mechanism of splitting failure.

3. Establishment of Constitutive Model

3.1. Virtual Working Principle and Control Equation

In the Toupin–Mindlin strain gradient theory, in addition to the conventional strain tensor ${\epsilon }_{ij}$ and stress tensor ${\sigma }_{ij}$, higher-order strain terms ${\eta }_{ijk}$ and higher-order stress tensors ${\tau }_{ijk}$ that describe the state of the material are also required [11]. The expression of strain and its strain gradient term is:

$${\epsilon }_{ij}=\left(u_{i,j}+u_{j,i}\right)/2, {\eta }_{ijk}=\left(u_{k,ij}+u_{k,ji}\right)/2$$

(1)

Not only strain and higher-order strain terms are symmetric tensors, their corresponding stress and higher-order stress terms are also symmetric tensors. So, suppose that the strain rate and higher-order strain rate are both decomposed into elastic and plastic parts:

$${\dot{\epsilon }}_{kl}={\dot{\epsilon }}_{kl}^e+{\dot{\epsilon }}_{kl}^p, {\dot{\eta }}_{lmn}={\dot{\eta }}_{{}_{lmn}}^e+{\dot{\eta }}_{{}_{lmn}}^p$$

(2)

Similar to the strain gradient plasticity model proposed by previous scholars [17,20,44,45,46,47], it is assumed that the plastic higher-order strain term contributes to the work per unit volume. Since the principle of virtual work is divided into the work of elastic strain, the work of plastic strain, and the work carried out by their higher-order terms, the virtual work of its internal unit volume V can be expressed as:

$$\delta W_i={\displaystyle {\int }_V\left({\sigma }_{ij}\delta {\epsilon }_{ij}^e+q_{ij}\delta {\epsilon }_{ij}^p+{\tau }_{ijk}\delta {\epsilon }_{ij,k}^p\right)}\mathrm{d}V$$

(3)

where $q_{ij}$ is called the micro stress tensor conjugated with plastic strain ${\epsilon }_{ij}^p$ work; ${\tau }_{ijk}$ is the higher-order stress tensor conjugated to the plastic strain gradient ${\epsilon }_{ij,k}^p$ work. Since the influence of $q_{ij}^{\prime }$ and ${\tau }_{ijk}^{\prime }$ on virtual work $\delta W_i$ are omitted in this paper, the above formula can be expressed as:

$$\delta W_i={\displaystyle {\int }_V\left({\sigma }_{ij}\delta {\epsilon }_{ij}+\left(q_{ij}-{\sigma }_{ij}^{\prime }\right)\delta {\epsilon }_{ij}^p+{\tau }_{ijk}\delta {\epsilon }_{ij,k}^p\right)}\mathrm{d}V$$

(4)

Using Gauss’s theorem to deal with the above formula, the following formula can be obtained:

$$\delta W_i={\displaystyle {\int }_S\left({\sigma }_{ij}{n}_j\delta u_i+{\tau }_{ijk}{n}_k\delta {\epsilon }_{ij}^p\right)}\mathrm{d}S-{\displaystyle {\int }_V\left[{\sigma }_{ij,j}\delta u_i+\left({\tau }_{ijk,k}+{\sigma }^{\prime }-q_{ij}\right)\delta {\epsilon }_{ij}^p\right]}\mathrm{d}V$$

(5)

where $n_i$ is the normal vector outward in the direction of the surface S. Since Equation (3) should be applied to any changes in plastic strain, and the second integral on the right side of Equation (5) disappears with any change, two sets of equilibrium equations can be determined:

$$\left\{\begin{array}{l}{\sigma }_{ij,j}=0\\ {\tau }_{ijk,k}+{\sigma }_{ij}^{\prime }-q_{ij}=0\end{array}\right.$$

(6)

According to the principle of virtual work, the corresponding conventional $T_i={\sigma }_{ij}{n}_j$ and high-order boundary conditions $t_{ij}={\tau }_{ijk}{n}_k$ are obtained.

3.2. The Establishment of the Yield Function

For materials that exhibit plasticity, the stress state ${\sigma }_{ij}$ reaches a certain yield condition before it can enter the plastic state. The yield condition can be expressed by independent stress components and hardening parameters:

$$f=f(\sigma ,H)=0$$

(7)

where H is the hardening parameter, which represents the relationship between yield stress and plastic deformation. For the classical elastic-plastic theory, the yield of a point in the material depends on its stress state and hardening parameters; for the theory of plastic strain gradient, the yield of the point in the material is not only related to the stress state and hardening parameters of the point itself, but also the influence of the surrounding points on the hardening parameters. When considering the influence of surrounding points on hardening parameters, point P with coordinate $\left(x,y\right)$ in the overall coordinate system is selected, and the neighborhood with surrounding radius R is taken as the neighborhood affecting point P, and the distance from any point in the neighborhood to the point is r, as shown in Figure 3. In the local coordinate system, the rectangular coordinates and polar coordinates of point P are $\left(r,\theta \right)$ and $\left(\xi ,\eta \right)$, respectively. Therefore, it is necessary to define a weighted average of hardening parameters to represent the influence of each point in the influence area on point P. The expression of the weighted average is:

$$\overline{H\left(x,y\right)}={\displaystyle \underset{\mathrm{\Omega }}{\int }g\left(r\right)h\left(x,y\right)\mathrm{d}\mathrm{\Omega }}$$

(8)

Among them, $g\left(r\right)$ is the material weight function, and the expression is selected as:

$$g\left(r\right)=\left\{\begin{array}{l}C\left(1-r^2/R^2 \right),r<R\\ 0,r\ge R\end{array}\right.$$

(9)

where, R is the radius of influence, $R={\rho }_{0}l$; l is the internal length parameter of the material, reflecting the length of the internal microstructure of the material; ${\rho }_0$ is the influence coefficient, for two-dimensional plane problems ${\rho }_0=={(3/4)}^{1/3}$; C is the weight function parameter, $C=3/\left(\pi R^2\right)$; r is the distance from any point in the neighborhood to that point. The weight function has the characteristics of a single point of decline and the nature of normalization:

$${\displaystyle \underset{\mathrm{\Omega }}{\int }g\left(r\right)\mathrm{d}\mathrm{\Omega }}=1$$

(10)

Use the above weight function to weighted average the hardening parameters at a certain point:

$$\overline{H\left(x,y\right)}={\displaystyle \underset{\mathrm{\Omega }}{\int }g\left(r\right)h(x+\xi ,y+\eta )\mathrm{d}\mathrm{\Omega }}$$

(11)

From this, the h(x,y) value of the internal length parameter l of the material can be obtained, which is expressed as the interaction between the microstructures. In order to obtain h(x + ξ, y + η) worthy of specifically related expressions, perform Taylor expansion of h(x,y) at this point:

$$\begin{array}{l}h(x+\xi ,y+\eta )=h\left(x,y\right)+(\xi \frac{\partial }{\partial x}+\eta \frac{\partial }{\partial y})h\left(x,y\right)+\frac{1}{2}!{(\xi \frac{\partial }{\partial x}+\eta \frac{\partial }{\partial y})}^{2}h\left(x,y\right)\\ +\frac{1}{3}! {(\xi \frac{\partial }{\partial x}+\eta \frac{\partial }{\partial y})}^{3}h\left(x,y\right)+\dots \end{array}$$

(12)

For isotropic materials, the odd-order term is 0 after integration, use $\mathrm{d}\xi \mathrm{d}\eta =r\mathrm{d}r\mathrm{d}\theta$ to transform the polar coordinates and rectangular coordinates in the integration and then take the integration before two items get:

$$\overline{H\left(x,y\right)}={\displaystyle \sum _{n=0,2}{\displaystyle \iint \frac{1}{n!}{\left(\xi \frac{\partial }{\partial x}+\eta \frac{\partial }{\partial y}\right)}^{2}h\left(x,y\right)}}\mathrm{d}\xi \mathrm{d}\eta$$

(13)

Simplify the above formula to obtain:

$$\overline{H\left(x,y\right)}=h\left(x,y\right)+C{\nabla }^{2}h\left(x,y\right)$$

(14)

where C is a constant; ${\nabla }^2$ is in Laplacian form, ${\nabla }^2=\frac{{\partial }^2}{\partial x^2}+\frac{{\partial }^2}{\partial y^2}$. Equation (14) can be understood as the hardening parameter expression considering the effect of strain gradient, and $\overline{H\left(x,y\right)}$ is used to replace $\kappa$ in the yield function of the classical plastic theory. The plastic strain gradient is introduced as an internal variable into the yield function of classical plastic theory. Introduce the hardening parameter considering the gradient term into the yield function:

$$F(\sigma ,H,{\nabla }^{2}H)=0$$

(15)

At the same time, in order to meet the constraints, the yield function should also satisfy the following formula:

$$F(\sigma ,H,{\nabla }^{2}H)\le 0,\lambda \ge 0,\lambda F=0$$

(16)

From this, the yield function of the plastic strain gradient theory is derived as:

$${\left\{\frac{\partial F}{\partial \sigma }\right\}}^T\sigma +\frac{\partial F}{\partial H}H+\frac{\partial F}{\partial {\nabla }^{2}H}{\nabla }^{2}H=0$$

(17)

The hardening modulus is defined as:

$$\kappa (H,{\nabla }^{2}H)=-\frac{1}{\lambda }\frac{\partial F}{\partial H}H$$

(18)

For further calculation, establish a certain relationship between the hardening function and the plastic multiplier:

$$H=\alpha \lambda$$

(19)

where α is a relational constant. Further collation can be obtained:

$${\left\{\frac{\partial F}{\partial \sigma }\right\}}^T\sigma -\kappa \lambda +\alpha \frac{\partial F}{\partial {\nabla }^{2}H}{\nabla }^2\lambda =0$$

(20)

Let $\alpha \frac{\partial F}{\partial {\nabla }^{2}H}=-\kappa l^2$, and further obtain the yield function of the plastic strain gradient theory as:

$${\left\{\frac{\partial F}{\partial \sigma }\right\}}^T\sigma -\kappa \lambda -\kappa l^2{\nabla }^2\lambda =0$$

(21)

In the plane problem, the strengthening parameter H can be expressed by the equivalent plastic strain $E^p$, so the yield model including the strain gradient term can be obtained from the formula:

$$\begin{array}{l}F=F(\sigma )-{\sigma }_y(E^p,{\nabla }^2{E}^p)\\ {\sigma }_y={\sigma }_0+\kappa E^p+\kappa l^2{\nabla }^2{E}^p\end{array}$$

(22)

where $\kappa l^2{\nabla }^2{E}^p$ reflects the influence of plastic strain gradient and material internal length parameters on the yield function; κ is the hardening modulus, $\kappa =\frac{1}{\lambda }\frac{\partial f}{\partial H}H$. Based on the Drucker–Prager strength criterion, the following yield function f and plastic potential function g are constructed:

$$\left\{\begin{array}{l}f=f\left(\mathit{\sigma },H\right)=\alpha I_1\left(\mathit{\sigma }\right)+\sqrt{J_2\left(\mathit{s}\right)}-H\left({\epsilon }_{ij}^p\right)-{\sigma }_y(E^p,{\nabla }^2{E}^p)\\ g=g\left(\mathit{\sigma },H\right)=\beta I_1\left(\mathit{\sigma }\right)+\sqrt{J_2\left(\mathit{s}\right)}-H\left({\epsilon }_{ij}^p\right)-{\sigma }_y(E^p,{\nabla }^2{E}^p)\end{array}\right.$$

(23)

where J₂ is the second invariant of stress deviation; $\alpha$ and $\beta$ are intensity parameters, and their expressions are as follows:

$$\left\{\begin{array}{l}\alpha =\frac{2\sin \phi }{\sqrt{3}\left(3-2\sin \phi \right)}\\ \beta =\frac{2\sin \psi }{\sqrt{3}\left(3-2\sin \psi \right)}\end{array}\right.$$

(24)

where $\phi$ is the friction angle and $\psi$ is expansion angle. When $\phi$ and $\psi$ are equal, it is the associated flow law.

3.3. Constitutive Equation

According to the strain gradient shaping theory, it is necessary to consider introducing the high-order stress related to the strain increment into the dissipation work generated by the high-order strain in a thermodynamically consistent manner [47,48]. The key to constructing the plastic constitutive equation considering the strain gradient is to define the effective stress $\Sigma$ work conjugated with the effective plastic strain rate ${\dot{E}}^p$ of the gradient enhancement to ensure the plastic work rate:

$$\Sigma {\dot{E}}^p=q_{ij}^D{\dot{\epsilon }}_{ij}^p+{\tau }_{ijk}^D{\dot{\epsilon }}_{ij,k}^p$$

(25)

Since $\Sigma {\dot{E}}^p\ge 0$, ${\dot{E}}^p$ can be defined as:

$${\dot{E}}^p=\sqrt{\frac{2}{3}{\dot{\epsilon }}_{ij}^p{\dot{\epsilon }}_{ij}^p+L^2{\dot{\epsilon }}_{ij,k}^p{\dot{\epsilon }}_{ij,k}^p}$$

(26)

where L is the dissipative length parameter. Thus, $\Sigma$ can be defined as:

$$\Sigma =\sqrt{\frac{2}{3}{q}_{ij}^D{q}_{ij}^D+L^{-2}{\tau }_{ijk}^D{\tau }_{ijk}^D}$$

(27)

From this, the dissipative plastic micro-stress tensor and the dissipative plastic high-order stress tensor are obtained, respectively:

$$q_{ij}^D=\frac{2}{3}\frac{\Sigma }{{\dot{E}}^p}{\dot{\epsilon }}_{ij}^p,{\tau }_{ijk}^D=\frac{\Sigma }{{\dot{E}}^p}{L}^2{\dot{\epsilon }}_{ij,k}^p$$

(28)

When there is damage in the material, the Helmholtz free energy function can be defined as the following form with the dissipation potential function according to the first law of thermodynamics:

$$\begin{array}{l}W=\Psi -\Phi =\\ \frac{1}{2}\left({\epsilon }_{ij}-{\epsilon }_{ij}^p\right)C_{ijkl}^e\left({\epsilon }_{kl}-{\epsilon }_{kl}^p\right)+\frac{1}{2}{\eta }_{ijk}^e{\Lambda }_{ijklmn}{\eta }_{lmn}^e+\frac{1}{2}{\epsilon }_{ij}^p{C}_{ijkl}^{pd}{\epsilon }_{kl}^p+\frac{1}{2}G{l}^2{\eta }_{ijk}^p{\eta }_{ijk}^p\\ -\sqrt{\frac{2}{3}{\dot{\epsilon }}_{ij}^p{\dot{\epsilon }}_{ij}^p+L^2{\dot{\epsilon }}_{ij,k}^p{\dot{\epsilon }}_{ij,k}^p}\sqrt{\frac{2}{3}{q}_{ij}^D{q}_{ij}^D+L^{-2}{\tau }_{ijk}^D{\tau }_{ijk}^D}\end{array}$$

(29)

where $C_{ijkl}^e$ is fourth-order elastic tensor; ${\Lambda }_{ijklmn}$ is the sixth-order elastic strain gradient damage tensor, ${\Lambda }_{ijklmn}=l^2{C}_{ijlm}^e{\delta }_{kn}$; $C_{ijkl}^{pd}$ is fourth-order plastic damage tensor, $C_{ijkl}^{pd}=\left(1-D+{\delta }_{ij}D/3\right)C_{ijkl}^p$; G is the shear modulus, $G=\left(1-d\right)E/\left[2\left(1+\nu \right)\right]$; l is the internal length parameter closely related to the micro-cracks and micro-defects in the rocks, and is an inherent parameter of the material; ${\delta }_{kn}$ is a Kronecker symbol; D is the modified damage variable, $D=\delta d$; $\delta$ is the correction parameter of the damage variable; d is the damage variable. Generally speaking, there is no damage when a material is in elastic state, and damage begins to appear when it enters plastic state. So, we define the following expression of the damage evolution equation:

$$d=1-\frac{{\kappa }_i}{\kappa }\left[\left(1-\alpha \right)+\alpha e^{-\beta \left(\kappa -{\kappa }_i\right)}\right]$$

(30)

where ${\kappa }_i$ is the threshold for damage initiation; $\alpha$ and $\beta$ are two material parameters; $C_{ijlm}^e$ is elastic tensor of the isotropic body, and its matrix expression is:

$$C_{ijlm}^e=\left[\begin{array}{cccccc}\frac{E\left(1-\nu \right)}{\left(1+\nu \right)\left(1-2\nu \right)}& \frac{E\nu }{\left(1+\nu \right)\left(1-2\nu \right)}& \frac{E\nu }{\left(1+\nu \right)\left(1-2\nu \right)}& 0& 0& 0\\ \frac{E\nu }{\left(1+\nu \right)\left(1-2\nu \right)}& \frac{E\left(1-\nu \right)}{\left(1+\nu \right)\left(1-2\nu \right)}& \frac{E\nu }{\left(1+\nu \right)\left(1-2\nu \right)}& 0& 0& 0\\ \frac{E\nu }{\left(1+\nu \right)\left(1-2\nu \right)}& \frac{E\nu }{\left(1+\nu \right)\left(1-2\nu \right)}& \frac{E\left(1-\nu \right)}{\left(1+\nu \right)\left(1-2\nu \right)}& 0& 0& 0\\ 0& 0& 0& \frac{E}{2\left(1+\nu \right)}& 0& 0\\ 0& 0& 0& 0& \frac{E}{2\left(1+\nu \right)}& 0\\ 0& 0& 0& 0& 0& \frac{E}{2\left(1+\nu \right)}\end{array}\right]$$

(31)

The plastic tensor $C^p$ can take the following form according to the incremental theory in the elastoplastic theory:

$$C^p=\frac{C^e\frac{\partial g}{\partial \mathsf{\sigma }}{\left(\frac{\partial f}{\partial \mathsf{\sigma }}\right)}^T{C}^e}{{\left(\frac{\partial f}{\partial \mathsf{\sigma }}\right)}^T{C}^e\frac{\partial g}{\partial \mathsf{\sigma }}-\frac{\partial f}{\partial H}{\left(\frac{\partial H}{\partial {\mathsf{\epsilon }}^p}\right)}^T\frac{\partial g}{\partial \mathsf{\sigma }}}$$

(32)

where H is the hardening function. Since rocks have both plastic hardening and plastic softening characteristics (Figure 4), its manifestation is defined as the following exponential form:

$$H=H\left({\epsilon }_{ij}^p\right)=k_1{\epsilon }_{ij}^p\exp \left(k_2{I}_1-k_3{\epsilon }_{ij}^p\right)$$

(33)

where k_i (i = 1, 2, 3) is the material parameter; I₁ is the first stress invariant representing the average stress or hydrostatic stress. The f and g in Equation (32) are the yield function and plastic potential function defined in Section 3.2, respectively.

The elastoplastic damage constitutive equation considering the strain gradient can be derived from the free energy function strain and strain gradient according to the second law of thermodynamics:

$${\sigma }_{ij}=\left\{\begin{array}{l}{C}_{ijkl}^e{\epsilon }_{kl}^e\\ \left(1-D+{\delta }_{ij}D/3\right)C_{ijkl}^p{\epsilon }_{kl}^p-\frac{2}{3}\frac{\Sigma }{E^p}{\epsilon }_{kl}^p\end{array}\right.$$

(34)

$${\tau }_{ijk}=\left\{\begin{array}{l}{l}^2{C}_{ijkl}^e{\delta }_{mn}{\eta }_{lmn}\\ \frac{\left(1-D+{\delta }_{ij}D/3\right)C_{ijkl}^p}{2\left(1+\nu \right)}{l}^2{n}_{ijk}^p-\frac{\Sigma }{E^p}{l}^2{\epsilon }_{kl}^p\end{array}\right.$$

(35)

4. Obtained Mapping Function

4.1. Basic Assumption

Some scholars [21,22,23,24,25,26,27] conducted preliminary research on the elastoplastic solution of non-circular caverns, including elastoplastic boundary and displacement field. Like these scholars, we apply the method of conformal function to solve the displacement and stress of the high side wall cavity. The specific steps are as follows: For the unit circular cavern with a radius of 1 in the physical z plane and the unit circle in the complex plane $\zeta$, the complex variable function method is used to conformally map the outer region of the circular cavern in the z plane to the outer region of the unit circle in the complex plane $\zeta$. In this case, the elastoplastic solution of the circular cavern in the physical plane z is the same as that of the unit circle in the complex plane $\zeta$. They extended it to non-circular caverns, and inversely mapped the plastic zone and displacement distribution on the complex plane $\zeta$ to obtain the plastic zone and displacement distribution on the physical plane.

The result of conformal transformation is to change the singly connected domain to be discussed on the z-plane into the inside or outside of the unit circle centered on the origin on the complex plane $\zeta$, and the boundary line of this area is transformed into the circumference of the unit circle. The position of any point on the $\zeta$ plane is expressed in polar coordinates as:

$$\zeta =\rho e^{i\theta }=\rho \left(\cos \theta +i\sin \theta \right)$$

(36)

The coordinate lines are circumferential and radial. The circumferential line is represented by $\rho =\mathrm{const}$, and the radial line is represented by $\theta =\mathrm{const}$. The two are always orthogonal at each point (Figure 5). According to the nature of conformal transformation, the curves corresponding to the two groups of lines on the z-plane are also orthogonal. The curves on the two groups of z-planes can be regarded as the coordinate lines of the orthogonal curve coordinate system on the z-plane. So we also use $\rho =\mathrm{const}$ and $\theta =\mathrm{const}$ to sign, but the $\rho =\mathrm{const}$ and $\theta =\mathrm{const}$ lines on the z-plane are no longer circular and radial lines. Through the rectangular coordinate system and orthogonal curvilinear coordinate system, the coordinate transformation formulas of displacement and stress components between the rectangular coordinate system in the physical plane and polar coordinate system in the complex plane can be established.

4.2. Coordinate Transformation of Displacement Component

From Figure 5, it can be obtained that:

$$\left\{\begin{array}{l}{u}_x=u_{\rho }\cos \alpha -u_{\theta }\sin \alpha \\ u_y=u_{\rho }\sin \alpha +u_{\theta }\cos \alpha \end{array}\right.$$

(37)

where $\alpha$ is the angle between the $\rho$ axis and the x axis in the z plane; $u_x$ and $u_y$ are the displacements in x and y directions of a point in the cavern in the physical z plane, and $u_{\rho }$ and $u_{\theta }$ are the displacements in $\rho$ and $\theta$ directions of a point in the complex plane $\zeta$. The points in two different planes are connected by conformal mapping function:

$$\begin{array}{cc}\hfill u_x+i{u}_y& =\left(u_p\cos \alpha -u_{\theta }\sin \alpha \right)+i\left(u_p\sin \alpha +u_{\theta }\cos \alpha \right)\hfill \\ & =\left(u_p+i{u}_{\theta }\right)\left(\cos \alpha +i\sin \alpha \right)=e^{i\alpha }\left(u_p+i{u}_{\theta }\right)\hfill \end{array}$$

(38)

It can be seen from the above formula that the corresponding relationship between the displacement components in two coordinate systems can be obtained only when $e^{i\alpha }$ is calculated. Since the $\rho =\mathrm{const}$ and $\theta =\mathrm{const}$ curve families on the z-plane are connected by the $\rho =\mathrm{const}$ and $\theta =\mathrm{const}$ curve families on the $\zeta$ plane through $z=\omega \left(\zeta \right)$, $e^{i\alpha }$ should start from this. Assuming that point B is given an infinitesimal increment dz along the $\rho$ axis on the z-plane, corresponding to the infinitely small increment $\mathrm{d}\zeta$ of point A on the $\zeta$ plane along the radial direction, then dz and $\mathrm{d}\zeta$ are, respectively:

$$\mathrm{d}z=e^{i\alpha }\left|\mathrm{d}z\right|,\mathrm{d}\zeta =e^{i\theta }\left|\mathrm{d}\zeta \right|$$

(39)

Because of $z=\omega \left(\zeta \right)$, we can obtain:

$$\mathrm{d}z={\omega }^{\prime }\left(\zeta \right)\mathrm{d}\zeta ,\left|\mathrm{d}z\right|=\left|{\omega }^{\prime }\left(\zeta \right)\mathrm{d}\zeta \right|=\left|{\omega }^{\prime }\left(\zeta \right)\right|\left|\mathrm{d}\zeta \right|$$

(40)

Then, we can obtain the following results:

$$e^{i\alpha }=\frac{\mathrm{d}z}{\left|\mathrm{d}z\right|}=\frac{{\omega }^{\prime }\left(\zeta \right)\mathrm{d}\zeta }{\left|{\omega }^{\prime }\left(\zeta \right)\right|\left|\mathrm{d}\zeta \right|}=e^{i\theta }\frac{{\omega }^{\prime }\left(\zeta \right)}{\left|{\omega }^{\prime }\left(\zeta \right)\right|}=\frac{\zeta }{\rho }\frac{{\omega }^{\prime }\left(\zeta \right)}{\left|{\omega }^{\prime }\left(\zeta \right)\right|}$$

(41)

Therefore, the corresponding relationship between the rectangular coordinate system in the physical plane and the displacement component in the polar coordinate system in the complex plane is:

$$u_x+i{u}_y=\frac{\zeta }{\rho }\frac{{\omega }^{\prime }\left(\zeta \right)}{\left|{\omega }^{\prime }\left(\zeta \right)\right|}\left(u_p+i{u}_{\theta }\right)$$

(42)

4.3. Coordinate Transformation of Stress Component

According to the axis formula in elastoplastic mechanics, the coordinate transformation formulas of the stress components in the following two different coordinate systems can be obtained:

$$\left\{\begin{array}{l}{\sigma }_{\rho }=\frac{1}{2}\left({\sigma }_x+{\sigma }_y\right)+\frac{1}{2}\left({\sigma }_x-{\sigma }_y\right)\cos 2\alpha +{\tau }_{xy}\sin 2\alpha \\ {\sigma }_{\theta }=\frac{1}{2}\left({\sigma }_x+{\sigma }_y\right)-\frac{1}{2}\left({\sigma }_x-{\sigma }_y\right)\cos 2\alpha -{\tau }_{xy}\sin 2\alpha \\ {\tau }_{\rho \theta }=\frac{1}{2}\left({\sigma }_y-{\sigma }_x\right)\sin 2\alpha +{\tau }_{xy}\cos 2\alpha \end{array}\right.$$

(43)

It is obvious from the above equation that there is ${\sigma }_x+{\sigma }_y={\sigma }_{\rho }+{\sigma }_{\theta }$, and it follows that:

$$\begin{array}{cc}\hfill {\sigma }_{\rho }-{\sigma }_{\theta }+2i{\tau }_{\rho \theta }& =\left({\sigma }_y-{\sigma }_x\right)\cos 2\alpha -{\tau }_{xy}\sin 2\alpha +2i\left[\frac{1}{2}\left({\sigma }_y-{\sigma }_x\right)\sin 2\alpha +{\tau }_{xy}\cos 2\alpha \right]\hfill \\ & =\left({\sigma }_y-{\sigma }_x\right)\left(\cos 2\alpha +i\sin 2\alpha \right)+2i{\tau }_{xy}\left(\cos 2\alpha +i\sin 2\alpha \right)\hfill \\ & =\left({\sigma }_y-{\sigma }_x+2i{\tau }_{xy}\right)e^{2i\alpha }\hfill \end{array}$$

(44)

According to the previous expression of $e^{i\alpha }$ (Equation (25)), $e^{2i\alpha }$ can be determined as:

$$e^{2i\alpha }==\frac{{\zeta }^2}{{\rho }^2}\frac{{\left[{\omega }^{\prime }\left(\zeta \right)\right]}^2}{{\left|{\omega }^{\prime }\left(\zeta \right)\right|}^2}=\frac{{\zeta }^2}{{\rho }^2}\frac{{\left[{\omega }^{\prime }\left(\zeta \right)\right]}^2}{{\omega }^{\prime }\left(\zeta \right)\overline{{\omega }^{\prime }\left(\zeta \right)}}=\frac{{\zeta }^2}{{\rho }^2}\frac{{\omega }^{\prime }\left(\zeta \right)}{\overline{{\omega }^{\prime }\left(\zeta \right)}}$$

(45)

Thus, the following corresponding relations of stress components in the two coordinate systems are obtained:

$${\sigma }_y-{\sigma }_x+2i{\tau }_{xy}=\frac{{\rho }^2}{{\zeta }^2}\frac{\overline{{\omega }^{\prime }\left(\zeta \right)}}{{\omega }^{\prime }\left(\zeta \right)}\left({\sigma }_{\theta }-{\sigma }_{\rho }+2i{\tau }_{\rho \theta }\right)$$

(46)

4.4. Approximate Solution of the Mapping Function

For high sidewall cavern, the mapping function from the outer domain of the cavern to the outer domain of the unit circle often takes the following form:

$$Z=\omega \left(\zeta \right)=c\left(\zeta -\frac{a_1}{\zeta }-\frac{a_{2}i}{2{\zeta }^2}-\cdots -\frac{a_{n-1}}{\left(n-1\right){\zeta }^{n-1}}-\frac{a_{n}i}{n{\zeta }^n}\right)$$

(47)

In the physical plane, the sidewall and bottom plate of the high sidewall cavern have a tendency to curve toward the inside of the cavern. That is, for the straight-line boundary part, when the number of mapping items is small, the boundary of the cavern determined by the mapping function is a curved edge. As the number of terms increases, the accuracy of the mapping becomes higher and higher, and the boundary of the cavity obtained by the mapping function gets closer to a straight line. Usually n = 4, that is, taking 5 terms has sufficient accuracy, so the integral function is taken as:

$$Z=\omega \left(\zeta \right)=c\left(\zeta -\frac{a_1}{\zeta }-\frac{a_{2}i}{2{\zeta }^2}-\frac{a_3}{3{\zeta }^3}-\frac{a_{4}i}{4{\zeta }^4}\right)$$

(48)

where, $a_i\left(i=1,2,3,4\right)$ is a real constant, and there is a corresponding real constant a for hole shapes with different height span ratios a_i, c. Based on the FORTRAN language, a mapping function parameter solving program was written according to the Schwarz–Christoffel formula. The parameters of the mapping function are shown in

Table 1. The parameters of the mapping function.

c	a1	a2	a3	a4
1.512308	0.269926	0.029007	0.366010	−0.158103

.

5. Numerical Results and Discussions

5.1. Basic Assumption

Considering the mechanical calculation model, when the length of the deep-buried horizontal cavern is large, it can be treated as a plane strain problem. Meanwhile, in classical elasto-plastic mechanics, a circular cavern under pressure and regardless of a body force is considered as a cylinder buried in an infinite elastoplastic body subjected to a uniform pressure q. Therefore, we define the following basic assumptions (Figure 6):

• $R_0$ is the radius of the circular cavern, $R_p$ is the radius of the plastic zone.
• The cavern is a deep buried cavern with infinite longitudinal length, and the uniform geostress $p_0$ is applied at the infinite distance of the cross section, regardless of body force.
• The surrounding rock of the cavern is isotropic elastic-plastic material.
• The compressive stress is positive and the tensile stress is negative.

5.2. Derivation of Balance Equation

Based on the basic assumptions, the cavern is an axisymmetric plane strain problem, ignoring body force. Due to the axial symmetry of the mechanical state of the surrounding rock of the deep cavern, it can be known that the displacement of any point in the surrounding rock is only a function of r, so the simplified geometric equation is:

$$\left\{\begin{array}{l}{\epsilon }_{rr}=\frac{\mathrm{d}u}{\mathrm{d}r},{\epsilon }_{\theta \theta }=\frac{u}{r}\\ {\eta }_{rrr}=\frac{{\mathrm{d}}^{2}u}{\mathrm{d}{r}^2},{\eta }_{\theta \theta r}=\frac{1}{r^2}\left(r\frac{\mathrm{d}u}{\mathrm{d}r}-u\right)\\ {\eta }_{r\theta \theta }={\eta }_{\theta r\theta }=\frac{1}{r}\left(\frac{\mathrm{d}u}{\mathrm{d}r}-\frac{u}{2r}\right)\end{array}\right.$$

(49)

According to the basic assumptions ignoring physical strength, the balance equation in the plane strain problem can be simplified to:

$$\frac{\partial {\sigma }_{rr}^{*}}{\partial r}+\frac{1}{r}\left({\sigma }_{rr}^{*}-{\sigma }_{\theta \theta }^{*}\right)=0$$

(50)

In the above equation, ${\sigma }_{rr}^{*}$ and ${\sigma }_{\theta \theta }^{*}$ are the generalized radial and tangential stresses in the surrounding rock after considering the higher order stress terms, respectively, and their expressions are as follows:

$$\left\{\begin{array}{l}{\sigma }_{rr}^{*}={\sigma }_{rr}-\left[\frac{\partial {\tau }_{rrr}}{\partial r}+\frac{1}{r}\left({\tau }_{rrr}-{\tau }_{\theta \theta r}-{\tau }_{r\theta \theta }\right)\right]\\ {\sigma }_{\theta \theta }^{*}={\sigma }_{\theta \theta }-\left[\frac{\partial {\tau }_{rrr}}{\partial r}+\frac{1}{r}\left({\tau }_{\theta r\theta }+{\tau }_{\theta \theta r}+{\tau }_{r\theta \theta }\right)\right]\end{array}\right.$$

(51)

Substituting Equation (50) into Equation (49), the static equilibrium equation with higher-order stress terms can be obtained:

$$\frac{\partial {\sigma }_{rr}}{\partial r}+\frac{{\sigma }_{rr}-{\sigma }_{\theta \theta }}{r}-\frac{{\partial }^2{\tau }_{rrr}}{\partial r^2}+\frac{1}{r}\frac{\partial \left({\tau }_{\theta r\theta }+{\tau }_{\theta \theta r}+{\tau }_{r\theta \theta }-2{\tau }_{rrr}\right)}{\partial r}+\frac{{\tau }_{\theta r\theta }+{\tau }_{\theta \theta r}+{\tau }_{r\theta \theta }}{r}=0$$

(52)

For the problem of thick-walled cylinders, the internal and external boundary conditions are:

$$\left\{\begin{array}{l}{T}_r=-{\sigma }_{rr}+\frac{\partial {\tau }_{rrr}}{\partial r}-\frac{1}{r}\left({\tau }_{rrr}+{\tau }_{\theta \theta r}\right)\left|{}_{r=R_0}\right.=0\\ R_r={\tau }_{rrr}\left|{}_{r=R_0}\right.=0\end{array}\right.$$

(53)

$$\left\{\begin{array}{l}{T}_r={\sigma }_{rr}-\frac{\partial {\tau }_{rrr}}{\partial r}+\frac{1}{r}\left({\tau }_{\theta \theta r}-{\tau }_{rrr}+2{\tau }_{r\theta \theta }\right)\left|{}_{r=R_e}\right.=p_0\\ R_r={\tau }_{rrr}\left|{}_{r=R_e}\right.=0\end{array}\right.$$

(54)

where T_r and R_r are surface traction and high-order surface traction, respectively.

5.3. Derivation of Related Equations for Surrounding Rock of Deep Circular Cavern

The damage is not considered when the surrounding rock is in an elastic state. Combining the previous constitutive equations (Equations (34) and (35)) and the equilibrium equation (Equation (52)) together, the equilibrium equation for the elastic stage of the surrounding rock of the cavern can be obtained:

$${\dot{u}}_{,rrrr}-\frac{11}{5r}{\dot{u}}_{,rrr}-\left(\frac{61}{20r^2}+\frac{\lambda +2G}{5c{l}^2}\right){\dot{u}}_{,rr}+\left(\frac{51}{20r^3}+\frac{\lambda +2G}{5cr{l}^2}\right){\dot{u}}_{,r}-\left(\frac{51}{20r^4}+\frac{\lambda +2G}{5cr{l}^2}\right)\dot{u}=0$$

(55)

In order to solve the above-mentioned fourth-order ordinary differential equation and transform the higher-order ordinary differential equilibrium equation into a system of first-order differential equations, the following auxiliary variables need to be defined first:

$$\left\{\begin{array}{l}{y}_1\left(r\right)=\dot{u}\\ y_2\left(r\right)=\frac{d{y}_1}{dr}\\ y_3\left(r\right)=\frac{d{y}_2}{dr}\\ y_4\left(r\right)=\frac{d{y}_3}{dr}\end{array}\right.$$

(56)

The following expression can be obtained:

$$f(\mathrm{Y},r)=\frac{d}{dr}\left(\begin{array}{c}{y}_1\\ y_2\\ y_3\\ y_4\end{array}\right)=\left(\begin{array}{c}{y}_2\\ y_3\\ y_4\\ k_1{y}_1+k_2{y}_2+k_3{y}_3+k_4{y}_4\end{array}\right)$$

(57)

The corresponding coefficients in the above formula are: $k_1=\frac{51}{20r^4}-\frac{\lambda +2G}{5c{r}^2{l}^2}$, $k_2=\frac{\lambda +2G}{5cr{l}^2}-\frac{51}{20r^3}$, $k_3=\frac{61}{20r^2}+\frac{\lambda +2G}{5c{l}^2}$, $k_4=\frac{11}{5r}$. The boundary conditions at the boundary of the elastic zone ($r=R_e$) can also be obtained:

$$\left\{\begin{array}{l}\left[-5c{l}^2{y}_4-\frac{5c{l}^2}{2r}{y}_3-\left(\lambda +2G+\frac{59c{l}^2}{4r^2}\right)y_2+\left(\frac{25c{l}^2}{2r^3}-\frac{\lambda }{r}\right)y_1\right]\left|{}_{r=R_e}\right.=0\\ c{l}^2\left(5y_3+\frac{4}{r}{y}_2-\frac{13}{4r^2}{y}_1\right)\left|{}_{r=R_e}\right.=0\end{array}\right.$$

(58)

From the previously assumed yield function (Equation (16)), using the method of calculating the stress and displacement of the surrounding rock plastic zone used [49], we can obtain:

$${\sigma }_{rr}-{\sigma }_{\theta \theta }=-\frac{6\alpha {\sigma }_{rr}\left[\left(1-D+{\delta }_{ij}D/3\right)C_{ijkl}^p-\frac{2\Sigma }{3E^p}\right]+2k_1{\sigma }_{rr}\exp \left[k_2{I}_1-k_3{\sigma }_{rr}/\left(\left(1-D+{\delta }_{ij}D/3\right)C_{ijkl}^p\right)-\frac{2\Sigma }{3E^p}\right]}{\left(\sqrt{1-3{\alpha }^2}+3\alpha \right)\left[\left(1-D+{\delta }_{ij}D/3\right)C_{ijkl}^p-\frac{2\Sigma }{3E^p}\right]}$$

(59)

Then, the equation of the plastic damage state of the surrounding rock of the cavern can be obtained:

$$\begin{array}{l}\left[\left(1-D+{\delta }_{ij}D/3\right)C_{ijkl}^p-\frac{2\Sigma }{3E^p}\right]\left(-5c{l}^2\right){\dot{u}}_{,rrrr}+\left[\left(1-D+{\delta }_{ij}D/3\right)C_{ijkl}^p-\frac{2\Sigma }{3E^p}\right]\left(-\frac{11c{l}^2}{r}\right){\dot{u}}_{,rrr}\\ +\left[\left(1-D+{\delta }_{ij}D/3\right)C_{ijkl}^p-\frac{2\Sigma }{3E^p}\right]\left(\lambda +2G-\frac{2G}{r^2}+\frac{19c{l}^2}{2r^2}\right){\dot{u}}_{,rr}+\left[\left(1-D+{\delta }_{ij}D/3\right)C_{ijkl}^p-\frac{2\Sigma }{3E^p}\right]\left(\frac{\lambda +2G}{r}+\frac{21c{l}^2}{4r^3}\right){\dot{u}}_{,r}\\ +\left[\left(1-D+{\delta }_{ij}D/3\right)C_{ijkl}^p-\frac{2\Sigma }{3E^p}\right]\left(\frac{51c{l}^2}{4r^4}-\frac{\lambda }{r^2}\right)\dot{u}+{\left[\left(1-D+{\delta }_{ij}D/3\right)C_{ijkl}^p-\frac{2\Sigma }{3E^p}\right]}^{\prime }\left(-10c{l}^2\right){\dot{u}}_{,rrr}\\ +{\left[\left(1-D+{\delta }_{ij}D/3\right)C_{ijkl}^p-\frac{2\Sigma }{3E^p}\right]}^{\prime }\left(-\frac{15c{l}^2}{r}\right){\dot{u}}_{,rr}+{\left[\left(1-D+{\delta }_{ij}D/3\right)C_{ijkl}^p-\frac{2\Sigma }{3E^p}\right]}^{\prime }\left(\lambda +2G+\frac{8c{l}^2}{r^2}\right){\dot{u}}_{,r}\\ +{\left[\left(1-D+{\delta }_{ij}D/3\right)C_{ijkl}^p-\frac{2\Sigma }{3E^p}\right]}^{\prime }\left(\frac{\lambda }{r}-\frac{51c{l}^2}{4r^3}\right)\dot{u}+{\left[\left(1-D+{\delta }_{ij}D/3\right)C_{ijkl}^p-\frac{2\Sigma }{3E^p}\right]}^{″}\left(-5c{l}^2\right){\dot{u}}_{,rr}\\ +{\left[\left(1-D+{\delta }_{ij}D/3\right)C_{ijkl}^p-\frac{2\Sigma }{3E^p}\right]}^{″}\left(-\frac{4c{l}^2}{r}\right){\dot{u}}_{,r}+{\left[\left(1-D+{\delta }_{ij}D/3\right)C_{ijkl}^p-\frac{2\Sigma }{3E^p}\right]}^{″}\left(\frac{13c{l}^2}{4r^2}\right)\dot{u}=0\end{array}$$

(60)

Substituting the previous static balance equation (Equation (36)) can obtain:

$$\begin{array}{l}\frac{\partial {\sigma }_{rr}}{\partial r}-\frac{6\alpha {\sigma }_{rr}\left[\left(1-D+{\delta }_{ij}D/3\right)C_{ijkl}^p-\frac{2\Sigma }{3E^p}\right]+2k_1{\sigma }_{rr}\exp \left[k_2{I}_1-k_3{\sigma }_{rr}/\left(\left(1-D+{\delta }_{ij}D/3\right)C_{ijkl}^p\right)-\frac{2\Sigma }{3E^p}\right]}{\left(\sqrt{1-3{\alpha }^2}+3\alpha \right)\left[\left(1-D+{\delta }_{ij}D/3\right)C_{ijkl}^p-\frac{2\Sigma }{3E^p}\right]r}\\ -\frac{{\partial }^2{\tau }_{rrr}}{\partial r^2}+\frac{1}{r}\frac{\partial \left({\tau }_{\theta r\theta }+{\tau }_{\theta \theta r}+{\tau }_{r\theta \theta }-2{\tau }_{rrr}\right)}{\partial r}+\frac{{\tau }_{\theta r\theta }+{\tau }_{\theta \theta r}+{\tau }_{r\theta \theta }}{r}=0\end{array}$$

(61)

Define the same auxiliary variable as the elastic segment (Equation (56)), but the expression should be written as follows:

$$f(\mathrm{Y},r)=\frac{d}{dr}\left(\begin{array}{c}{y}_1\\ y_2\\ y_3\\ y_4\end{array}\right)=\left(\begin{array}{l}{y}_2\\ y_3\\ y_4\\ (g_{1}A{y}_4+g_{2}A{y}_3+g_{3}A{y}_2+g_{4}A{y}_1\\ +g_5{A}^{\prime }{y}_4+g_6{A}^{\prime }{y}_3+g_7{A}^{\prime }{y}_2+g_8{A}^{\prime }{y}_1\\ +g_9{A}^{″}{y}_3+g_{10}{A}^{″}{y}_2+g_{11}{A}^{″}{y}_1)/B\end{array}\right)$$

(62)

The corresponding coefficients in the above formula are: $g_1=-\frac{11c{l}^2}{r}$, $g_2=\lambda +2G-\frac{2G}{r^2}+\frac{19c{l}^2}{2r^2}$, $g_3=\frac{\lambda +2G}{r}+\frac{21c{l}^2}{4r^3}$, $g_4=\frac{51c{l}^2}{4r^4}-\frac{\lambda }{r^2}$, $g_5=-10c{l}^2$, $g_6=-\frac{15c{l}^2}{r}$, $g_7=\lambda +2G+\frac{8c{l}^2}{r^2}$, $g_8=\frac{\lambda }{r}-\frac{51c{l}^2}{4r^3}$, $g_9=-5c{l}^2$, $g_{10}=-\frac{4c{l}^2}{r}$, $g_{11}=\frac{13c{l}^2}{4r^2}$. The expressions of A, B and E are as follows:

$$\left\{\begin{array}{l}A=\left(1-D+{\delta }_{ij}D/3\right)C_{ijkl}^p-\frac{2\Sigma }{3E^p}\\ B=\left[\left(5c{l}^2\right)A-\left(\frac{13c{l}^2}{4r^2}\right)E{y}_1+\left(\frac{4c{l}^2}{r}\right)E{y}_2+\left(5c{l}^2\right)E{y}_3\right]\\ E=\frac{l^2{y}_3\left[1-\frac{{\kappa }_i}{\kappa }\left[\left(1-\alpha \right)+\alpha e^{-\beta \left(\kappa -{\kappa }_i\right)}\right]\right]}{{\left[\sqrt{\frac{4}{3}{y}_2^2+\frac{4}{9}{\epsilon }_0^2+l^2\left(y_3^2+3\frac{y_2^2}{r^2}+\frac{3}{2}\frac{y_1^2}{r^4}-4\frac{y_2{y}_1}{r^3}\right)}\right]}^3}\end{array}\right.$$

(63)

The boundary conditions at the inner boundary ($r=R_0$) of the post peak strain softening zone can be obtained:

$$\left\{\begin{array}{l}\left(5c{l}^2\right)A{y}_4+\left(-\frac{5c{l}^2}{2r}\right)A{y}_3+\left(-\frac{11c{l}^2}{r^2}-\lambda -2G\right)A{y}_2+\left(-\frac{c{l}^2}{2r^3}-\lambda \right)A{y}_1\\ {\left.+\left(5c{l}^2\right)A^{\prime }{y}_3+\left(\frac{4c{l}^2}{r}\right)A^{\prime }{y}_2+\left(-\frac{13c{l}^2}{4r^2}\right)A^{\prime }{y}_1\right|}_{r=R_0}=0\\ {\left.\left(5c{l}^2\right)A{y}_3+\left(\frac{4c{l}^2}{r}\right)A{y}_2+\left(-\frac{13c{l}^2}{4r^2}\right)A{y}_1\right|}_{r=R_0}=0\end{array}\right.$$

(64)

5.4. Solution Method for Displacement and Stress

The equilibrium equation (Equation (55)) in the elastic state and the equilibrium equation (Equation (61)) in the plastic damage state are both fourth-order ordinary differential equations (ODE), and the corresponding boundary conditions (Equation (57)) and (Equation (64)) are third-order ordinary differential equations. It is almost impossible to obtain analytical solutions to the above problems. Therefore, the numerical solutions of the above problems are solved by using MATLAB numerical analysis software.

The boundary conditions given in the above problems belong to boundary value problems at the inner and outer boundaries of surrounding rock; when the ODE45 function is used in MATLAB numerical analysis software, it should be first transformed into initial value problems. ODE45 is a special function for solving differential equations in MATLAB numerical analysis software, which adopts the fourth-fifth order Runge–Kutta algorithm. Its fourth-order algorithm provides candidate solutions, and the fifth-order algorithm controls errors. It is an adaptive step size numerical solution of ordinary differential equations and solves Nonstiff ordinary differential equations.

When applying the ODE45 function to solve high-order ordinary differential equations, it is necessary to provide initial solutions for calculation intervals and differential equations, and it is also necessary to convert high-order differential equations into first-order differential equations. The first-order differential equations transformed by higher-order differential equations have been given in the previous section. In order to solve the differential equation, it is still necessary to give the initial solution of the differential equation on the boundary. In the above problem, the surrounding rock in the range of A is in an elastic state. It is required to solve the fourth-order ordinary differential equation in the elastic zone (Equation (55)), and the initial solution on the outer boundary (Equation (58)) must be given, that is, the derivative values of each order of displacement at $r=R_e$ ($y_1\left(R_e\right)=u$, $y_2\left(R_e\right)=u_{,r}$, $y_3\left(R_e\right)=u_{,rr}$ and $y_4\left(R_e\right)=u_{,rrr}$).

From the calculation and analysis of Zhao et al. (2007), it can be seen that the displacement value obtained by applying two different theories (classical theory and strain gradient theory) in the elastic stage is not much different from its first derivative value. Therefore, the displacement value ($u\left|{}_{R_e}\right.$) and the displacement first derivative value ($u_{,r}\left|{}_{R_e}\right.$) at E can be roughly estimated from the analytical solution of the thick-walled cylinder in the classical elastic theory, and then $u\left|{}_{R_e}\right.$ and $u_{,r}\left|{}_{R_e}\right.$ are substituted into the boundary conditions (Equation (58)) to obtain $u_{,rr}\left|{}_{R_e}\right.$ and $u_{,rrr}\left|{}_{R_e}\right.$. In this way, the initial solution on the boundary required to solve the higher-order differential equation is obtained. The analytical solution of the thick-walled cylinder in the classical elastic theory is as follows:

$$\left\{\begin{array}{l}u=\frac{1-\upsilon }{E}\frac{p_a{a}^2-p_b{b}^2}{b^2-a^2}r+\frac{1+\upsilon }{E}\frac{\left(p_a-p_b\right)a^2{b}^2}{b^2-a^2}\frac{1}{r}\\ u^{\prime }=\frac{1-\upsilon }{E}\frac{p_a{a}^2-p_b{b}^2}{b^2-a^2}-\frac{1+\upsilon }{E}\frac{\left(p_a-p_b\right)a^2{b}^2}{b^2-a^2}\frac{1}{r^2}\end{array}\right.$$

(65)

In the above equation, E and $\nu$ are the elastic modulus and Poisson’s ratio of the surrounding rock respectively, a and b are the inner and outer diameters of the cylinder respectively, and $p_a$ and $p_b$ are the pressures on the inner and outer walls of the cylinder, respectively.

Due to the continuity of the deformation of the surrounding rock of the cavern, $r=R_p$ is both the inner boundary of the elastic zone and the outer boundary of the plastic damage zone, that is to say, the surrounding rock conforms to the equilibrium equation of the elastic zone (Equation (39)) and plastic damage zone (Equation (45)). Therefore, the following method is used to determine the size of $R_p$: first, the calculation interval is set as $r\in \left[R_0\hspace{0.33em},\hspace{0.33em}{R}_e\right]$; combined with the initial values of the displacement derivatives $u\left|{}_{R_e}\right.$, $u_{,r}\left|{}_{R_e}\right.$, $u_{,rr}\left|{}_{R_e}\right.$ and $u_{,rrr}\left|{}_{R_e}\right.$ on the outer boundary of the elastic zone are estimated and determined; the ODE45 function in MATLAB is used to solve the differential equations (Equation (57)), and the numerical solution of each point in the range of $R_0\le r\le R_e$ can be obtained. According to the research of Gao et al. [8], find a point that satisfies the Equation (66), and the distance between this point and the center of the roadway is the size of $R_p$.

$$\tilde{\epsilon }=\sqrt{\frac{4}{3}{u^{\prime }}^2+\frac{4}{9}{\epsilon }_0^2+l^2\left({u^{″}}^2+3\frac{{u^{\prime }}^2}{r^2}+\frac{3}{2}\frac{u^2}{r^4}-4\frac{u^{\prime }u}{r^3}\right)}={\epsilon }_f$$

(66)

After the boundary between the elastic zone and plastic damage zone is determined, the displacement derivatives ($u\left|{}_{R_e}\right.$, $u_{,r}\left|{}_{R_e}\right.$, $u_{,rr}\left|{}_{R_e}\right.$ and $u_{,rrr}\left|{}_{R_e}\right.$) on the inner boundary of the elastic zone can be calculated; it is taken as the initial value on the outer boundary of the plastic damage zone, and the calculation interval is set as $r\in \left[R_0\hspace{0.33em},\hspace{0.33em}{R}_p\right]$; the numerical solution ($u$, $u_{,r}$, $u_{,rr}$ and $u_{,rrr}$) of each point in the range of $R_0\le r\le R_p$ can be obtained by solving the differential equations (Equation (62)) with ODE45 function in MATLAB.

Take the derivative values of the displacements ($u\left|{}_{R_0}\right.$, $u_{,r}\left|{}_{R_0}\right.$, $u_{,rr}\left|{}_{R_0}\right.$ and $u_{,rrr}\left|{}_{R_0}\right.$) on the boundary of the plastic damage zone and substitute them into Equations (67) and (68) to check the calculation errors. This calculation error is caused by the estimated initial value ($u\left|{}_{R_e}\right.$, $u_{,r}\left|{}_{R_e}\right.$, $u_{,rr}\left|{}_{R_e}\right.$ and $u_{,rrr}\left|{}_{R_e}\right.$) on the outer boundary of the elastic zone. If the error meets the requirements (${\mathbf{F}}^0\le TOL$), the calculation ends to obtain the numerical solution of the entire calculation area; if the error does not meet the requirements (${\mathbf{F}}^0>TOL$), adjust the initial value on the outer boundary of the elastic zone and repeat the above calculation steps until the error meets the requirements, and obtain the numerical solution of the entire calculation area. The calculation flowchart of the solution process is shown in Figure 7.

$$\left\{\begin{array}{l}{f}_T^0=\left(5c{l}^2\right)A{u}_{,rrr}+\left(-\frac{5c{l}^2}{2r}\right)A{u}_{,rr}+\left(-\frac{11c{l}^2}{r^2}-\lambda -2G\right)A{u}_{,r}\\ +\left(-\frac{c{l}^2}{2r^3}-\lambda \right)Au+\left(5c{l}^2\right)A^{\prime }{u}_{,rr}+\left(\frac{4c{l}^2}{r}\right)A^{\prime }{u}_{,r}{\left.+\left(-\frac{13c{l}^2}{4r^2}\right)Du\right|}_{r=R_0}\\ {\left.f_R^0=\left(5c{l}^2\right)A{u}_{,rr}+\left(\frac{4c{l}^2}{r}\right)A{u}_{,r}+\left(-\frac{13c{l}^2}{4r^2}\right)Au\right|}_{r=R_0}\end{array}\right.$$

(67)

$${\mathbf{F}}^0=\sqrt{{\left(f_T^0\right)}^2+{\left(f_R^0\right)}^2}$$

(68)

5.5. Analysis of Theoretical Calculation Results

In this paper, the radial displacement, radial stress, and tangential stress distribution of the high sidewall cavern are obtained using rock parameters (

Table 2. Rock parameters in theoretical calculation.

Initial Elastic Modulus/GPa	Unit Weight/(KN·m−3)	Compressive Strength/MPa	Tensile Strength/MPa	Peak Strain	Ultimate Strain	Poisson’s Ratio	Lame Constant/GPa	Internal Length Parameter/m
41.50	26.60	128.80	8.00	$9.38\times {10}^{-3}$	$16.68\times {10}^{-3}$	0.27	16.34	0.01

) of the Pubugou Hydropower Station as calculation parameters. The top of the main powerhouse is 360 m from the ground surface and the cross-section size of the cavern cross-section is measures 26.8 m × 70.1 m. Then, the theoretical calculation results are compared and analyzed with the geomechanical model test results [50] to initially verify the applicability of the splitting failure elastic-plastic damage model in the study of splitting damage in a deep cavern. The geomechanical model test program and results are shown in Figure 8 and Figure 9, respectively.

Using the above calculation parameters, the numerical solution of the radial displacement, radial stress, and tangential stress of the rock of the high sidewall cavern is calculated according to the solution process shown in Figure 7. In this section, in order to facilitate the comparison between theoretical calculation results and model test results, all model test results are converted into prototypes based on the principles of similitude. Figure 10 shows the radial displacement distribution diagram of the surrounding rock of the high sidewall cavern. The theoretically calculated values and experimentally measured values of the radial displacement, radial stress, and tangential stress at the corresponding positions of the measuring points around the cavern are listed in

Table 3. Radial displacement around the cavern.

Location	Category	Radial Displacement Value/mm
		0.1 L	0.5 L	L	1.5 L	2 L	3 L
Right sidewall	Model test	63.0	27.0	36.0	9.0	18.0	3.0
	Theoretical calculation	58.30	14.34	30.82	4.82	10.94	0.30
Right hance	Model test	57.0	24.0	30.0	7.5	15.0	3.0
	Theoretical calculation	49.92	12.98	25.48	3.50	9.71	0.30
Vault	Model test	42.0	30.0	24.0	15.0	9.0
	Theoretical calculation	32.75	12.81	10.17	3.24	0.14
Left hance	Model test	57.0	18.0	33.0	7.5	15.0	1.5
	Theoretical calculation	49.92	12.98	25.48	3.50	9.71	0.30
Left sidewall	Model test	64.5	21.0	34.5	6.0	16.5	1.5
	Theoretical calculation	58.30	14.34	30.82	4.82	10.94	0.30

.

It can be seen from Figure 10 that the radial displacement at the high sidewall of the cavern obtained by using the elastic-plastic damage model of splitting failure shows an oscillating attenuation, and the oscillation amplitude gradually decreases with an increase in the distance from the sidewall, which is consistent with the variation measured by the model test. Because the theoretical model ignores the influence of excavation damage, the calculated displacement value is less than the displacement value of the model test. The difference close to the sidewall is as small as no more than 13%, and the value at the trough of the oscillating change is quite different, with the maximum difference being 12.66 mm. Comparing the measured values in Table 3 with the calculated values, it can be seen that they are very close in most positions. It is preliminarily confirmed that the elasto-plastic damage model of splitting failure is applicable to the study of splitting failure of high sidewall cavern.

The stress distribution of rock of the high sidewall cavern was also obtained by using the splitting failure elastoplastic damage model, as shown in Figure 11 and Figure 12.

Table 4. Radial stress around the cavern.

Location	Category	Radial Stress Value/MPa
		0.1 L	0.5 L	L	1.5 L	2 L	3 L
Right sidewall	Model test	−1.32	−6.20	−4.90	−8.61	−7.96	−11.77
	Theoretical calculation	−2.62	−8.33	−4.90	−9.68	−9.13	−12.76
Right hance	Model test	−1.22	−7.24	−4.79	−9.59	−8.89	−12.57
	Theoretical calculation	−3.82	−8.53	−7.07	−10.67	−10.31	−13.45
Vault	Model test	−2.08	−6.94	−8.21	−10.39	−12.63
	Theoretical calculation	−4.50	−8.36	−10.41	−12.41	−14.08
Left hance	Model test	−1.38	−7.57	−4.90	−9.12	−8.21	−11.93
	Theoretical calculation	−3.82	−8.53	−7.07	−10.67	−10.31	−13.45
Left sidewall	Model test	−1.82	−6.60	−4.88	−8.17	−7.98	−11.85
	Theoretical calculation	−2.62	−8.33	−4.90	−9.68	−9.13	−12.76

and

Table 5. Tangential stress around the cavern.

Location	Category	Tangential Stress Value/MPa
		0.1 L	0.5 L	L	1.5 L	2 L	3 L
Right sidewall	Model test	−3.83	−6.94	−20.36	−17.48	−19.38	−16.08
	Theoretical calculation	−4.82	−12.13	−26.29	−22.34	−23.51	−17.28
Right hance	Model test	−2.09	−8.39	−19.91	−16.99	−18.68	−16.39
	Theoretical calculation	−5.33	−12.77	−24.16	−21.16	−22.75	−16.79
Vault	Model test	−23.39	−19.83	−16.4	−15.88	−15.12
	Theoretical calculation	−26.59	−23.49	−20.46	−19.55	−17.62
Left hance	Model test	−2.92	−9.24	−19.35	−17.72	−19.55	−16.18
	Theoretical calculation	−5.33	−12.77	−24.16	−21.16	−22.75	−16.79
Left sidewall	Model test	−2.07	−8.58	−18.71	−17.95	−19.29	−16.97
	Theoretical calculation	−4.82	−12.13	−26.29	−22.34	−23.51	−17.28

present the radial and tangential stresses of theoretical calculations and model tests, respectively. Figure 11 shows the comparison of the radial stress distribution between theoretical calculation and model test. We find that the distribution of radial stresses is similar to that of displacements. The calculated radial stress value is relatively close to the test monitoring result, and the maximum difference does not exceed 2 MPa. As shown in Figure 12, the theoretically calculated tangential stress distribution is also consistent with the model test. When the distance from the sidewall is less than L, the rock is in a plastic state and the tangential stress is released; when the distance from the sidewall is greater than L, the tangential stress is concentrated and the rock is in an elastic state. The tangential stresses at the vault are all greater than the in-suit value (σ_θ = 15.25 MPa). Due to ignoring the excavation disturbance, the calculated value of the tangential stress is generally greater than the test value. When the tangential stress just appears to be concentrated, the gap is large, and the maximum gap between the two is 28.83%.

The numerical solution of the radial displacement and radial stress of the surrounding rock of the high side wall cavern using the splitting failure elastoplastic damage softening model for theoretical calculations shows that the oscillation attenuation occurs alternately with wave crests and wave troughs. This is completely different from the monotonic attenuation obtained by using the classic continuous elastoplastic model and is consistent with the change rule obtained by the geomechanical model test. This indicates that the split failure elastoplastic damage model can explain the split failure of high side wall caverns under high ground stress conditions.

In order to study the formation mechanism of splitting failure of high sidewall caverns of hydropower stations, we compared the distribution of tangential stress around the high sidewall caverns under different height–span ratios, as shown in Figure 13. We find that the stress distribution of the surrounding rock of the high side wall cavern is significantly affected by the height–span ratio. The tangential stress is unevenly distributed along the boundary of the cavern and is symmetrically distributed along the central axis of the cavern. There are obvious stress fluctuations in surrounding rock near the vault, arch toe, and lower corner of the side wall. With an increase in the height–span ratio of the cavern, the tangential stress of the vault and the lower corner becomes larger and larger. The maximum tangential stress at the midpoint of the side wall shows a decreasing trend. The release of the tangential stress increases gradually, and the damage of the surrounding rock increases gradually.

6. Conclusions

Based on the strain gradient theory and continuum damage mechanics, the splitting failure elasto-plastic damage model is established. By using complex variable function mapping, the variations of radial displacement, radial stress, and tangential stress of surrounding rock of high sidewall cavern are obtained. The theoretical calculation results and model test results are compared and analyzed, and the following research conclusions are obtained:

(1)

A yield function with strain gradient is established, which can describe the plastic hardening and softening of quasi-brittle materials such as rock at the same time. Based on this, an elastic-plastic damage softening constitutive model considering strain gradient is established.

(2)

The calculation program is compiled by MATLAB software, and the variations of radial displacement, radial stress, and tangential stress of surrounding rock of circular cavern are calculated. Then, the variations of high side wall cavern under high ground stress are mapped by using the obtained mapping function.

(3)

The theoretical calculation value and the measured value of the model test are in good agreement in both the value and the variations, which confirms the applicability of the splitting failure elasto-plastic damage model.

(4)

Through comparison, the variations of tangential stress in high side wall caverns with different height–span ratios are obtained: as the height-to-span ratio of the cavern increases, the tangential stress release at the sidewall gradually increases, and the damage to the surrounding rock also gradually increases.

7. Recommendations for Future Research

In this paper, an elastoplastic damage softening constitutive model that can simultaneously describe the plastic hardening and softening processes of rock and other quasi-brittle materials considering the strain gradient is established. In addition, we use the model to calculate the displacement and stress values of the split failure of the high side wall cavity. Among them, the internal length parameter of the material is related to the internal microstructure of the rock material and is a key parameter in the strain gradient theory. At present, the method of determining the internal length parameters of materials by measuring the width of the shear band can only be applied to specific models and materials and does not have broad significance. Therefore, it is necessary to further study the corresponding relationship between the internal length parameters and their internal microstructure in order to accurately determine the internal length parameters.