Failure Mechanism and Numerical Simulation of Splitting Failure for Deep High Sidewall Cavern Under High Stress

With the increase of the depth of the underground engineering, the phenomenon of splitting failure of the deep rock will appear, which is very different from the shallow cavern. In order to reveal the formation mechanism of splitting failure, mechanical model test and numerical simulations of splitting failure were carried out respectively. Using the Pubugou Hydropower Station as the engineering background, a three-dimensional (3D) geomechanical model test was conducted relying on a high stress three-dimensional load test system. The splitting failure phenomenon of high sidewall cavern was observed, and the oscillation variations of displacement and stress were measured. Based on strain gradient theory and continuum damage mechanics, an elastic–plastic damage softening model for splitting damage was established. The relationship between rock damage and energy dissipation was analyzed. Based on the strain energy density theory, the splitting failure criterion based on the strain gradient is established. A numerical analysis method for splitting damage was proposed, and a regional disintegration calculation program was developed based on a commercial finite element code. The numerical simulation results are in basic agreement with the 3D geomechanical model test.


Introduction
As a renewable energy source, hydropower has received increasing attention in China. Many large hydropower projects are under construction or have been completed. A series of complex discontinuous deformations were detected during the construction of these projects. This phenomenon poses a serious threat to the construction safety and long-term stability of high sidewall cavern. Among them, splitting failure is a special failure mode which is often encountered.
The splitting failure of cavern with high sidewall is an engineering failure phenomenon caused by the release of strain energy stored in rock mass due to excavation-induced unloading, which results in the initiation, expansion and penetration of cracks in surrounding rock. The splitting failure of hard rock has been extensively recorded and described by many researchers. Ortlepp et al. (1994Ortlepp et al. ( , 1997 considered that splitting failure is a form of failure caused by excavation unloading under high geostress conditions. Splitting cracks is generally parallel to the direction of maximum tangential stress. As the splitting failure develops, a v-shaped notch will eventually form. Martin and Maybee (2000) observed a field statistical analysis of 178 pillar failures in Canadian hard-rock mines. The study indicates that when the width-toheight ratio of hard rock pillars is less than 2.5, the dominant failure mode is progressive slabbing and spalling which eventually leads to an hour-glass shape. Hibino (2001) and Yoshida et al. (2004) conducted onsite observations on 16 high sidewall caverns of the Japanese power system. Monitoring results show that the opening deformation of splitting cracks in the cavern accounted for a considerable part of the total deformation of the surrounding rock.
The damage mechanics method has good applicability to simulate the stiffness degradation and strain softening characteristics of rock materials, and can simulate the nonlinear behavior of rock materials after the peak more comprehensively to be widely used in investigating the failure responses of rock caverns. Li and Feng (2013) 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 zonal disintegration. Tao et al. (2013) 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. (2012) studied the damage evolution around deep underground openings under multiaxial stress through 3D numerical tests. The results show that biaxial stress parallel to the free surface constrained by zero or low confinement stress contributes to spalling or parallel fracture of the slab.
According to the classical continuum theory, the stress at a point of the material is only a function of the strain and deformation history at that point, independent of the strain in the region of the adjacent zone. However, the classical continuum theory is not strictly satisfied due to the imperfect nature of the rock material. The strain gradient effect does not control the deformation and strength of geotechnical materials in their elastic phase, but when the material reaches its peak strength, the strain gradient becomes large to become the cause of strain localization initiation and development. In recent years, the strain gradient theory that introduces internal length parameters of materials have gradually been developed as a nonlocal theory. Strain gradient theory can be used to describe the local and post-peak strain softening of rock (Steinmann, 1994;Zbib and Aifantis, 1988a, b). The local strain in the rock is usually accompanied by strain softening. It is precisely because the influence of strain gradient is ignored that the previous numerical analysis methods are difficult to explain the problems of splitting failure (Zhang et al. 2017). Based on the strain gradient theory and plastic deformation theory, Gao et al. (2019) established the elastic-plastic damage softening model and elaborated the formation mechanism and development law of zonal disintegration.
With its advantages of visualization, intuition and realism, the geomechanical model test is of great significance to the study of damage mechanism of underground caverns (Chen et al. 2013a, b, c;Zhang et al. 2013Zhang et al. , 2017Gao et al. 2018). The method can also simulate excavation and unloading under real three-dimensional high ground stress conditions, reproduce the damage pattern of the cavern, and capture the variations in deformation and stress (Zhang et al., 2019a, b).
Although great progress has been made in studying the microscopic mechanisms of rock specimens (Tang et al. 2001;Liu et al. 2002;Fang and Harrison 2002;Sahouryeh et al. 2002;Dyskin et al. 2003), laboratory or field experiments (Bauch and Lempp 2004;Wong et al. 2004;Li et al. 2011;Jiang et al. 2017;Liu et al. 2017;Gong et al. 2012Gong et al. , 2018 and numerical simulations (Zhu et al. 2010;Wang et al. 2016;Sharif et al. 2019;Wowk et al. 2019), so far there is no widely accepted theoretical explanation for splitting failure, and further improvement of theoretical methods is still needed. In this paper, the 3D geomechanical model test of splitting failure is firstly introduced, and the variations of displacement and stress around the high sidewall cavern are given. Based on the strain gradient theory and continuous damage theory, an elasticplastic damage model for splitting failure is established. According to the strain energy density theory, the energy damage criterion for splitting failure is proposed. In order to consider the strain gradient effect in the numerical simulation, a high-order hexahedral cell is constructed using Hermite interpolation function, and the shape function and stiffness matrix of the cell are derived. The numerical analysis method of regional rupture was proposed, and the corresponding calculation program was developed based on the ABAQUS platform. The results of numerical simulation and model test are in general agreement. The reliability of the damage softening model and the numerical analysis method of regional rupture are proved. Finally, on the basis of theoretical analysis and numerical simulation, the formation mechanism of deep high sidewall cavern splitting failure is described.
2 Three Dimensional Geomechanical Model Test

The Scheme of the Geo-mechanical Model
We take the main powerhouse project of the Pubugou Hydropower Station as the engineering background. The underground powerhouse cavern is located in the granite rock mass on the left bank downstream of the dam axis. During the construction of the main powerhouse, several splitting cracks appeared on the sidewall, causing cracking of the concrete shotcrete on the sidewall. The overlying rock thickness is 360 m, and cross-section size of the main powerhouse is 26.8 m 9 70.1 m. According to the similarity principle (Zhang et al. 2019c) and considering the scale of the main powerhouse and the size of the test system (0.7 m 9 0.7 m 9 0.7 m), the geometrical similarity of the model test is selected as C L = C P /C M = 300.
The extent and in-situ stress simulated by the model test are shown in Fig. 1.
We chose IBSCM (iron-barites-silica cementation material) as the analogous material (Zhang et al. 2019a, b, c). The physical and mechanical parameters of the granite and the analogous material and its mixture proportion are provided in Tables 1 and 2, respectively. Figure 2a shows the three typical monitoring sections to accurately reflect the deformation and stress changes of surrounding rock during the testing, in which section I and III are stress monitoring sections (section I is radial stress monitoring section and section III is tangential stress monitoring section), and section II is deformation monitoring section. As shown in Fig. 2b, five measuring lines are arranged for each monitoring section, and each measuring line is arranged with six measuring points. Due to the size of the model, in the vertical measurement Line C is arranged five measuring points.

Construction, Excavation and Overloading of the Geo-mechanical Model
To ensure the integrity and homogeneity of the surrounding rock, the geomechanical model is paved, compacted and air-dried layer by layer (Zhang et al. 2017). The body of the model is divided into 10 layers of 70 mm each. When the design elevation is reached, the grooves are cut and the measuring sensors are buried. To avoid delamination at the adjacent interfaces of the model, the surface of the previous layer of material must be chiseled rough and moistened with alcohol before the next layer is constructed. Boundary loads were applied to the model using a 3D loading system (Fig. 3). The model body was loaded in increments of 0.1 times the original in-suit stress until the design value was reached. After the design stress was reached, the load was held constant for 24 h to produce an initial stress field in the model that was consistent with the original stress field. Due to the limitations of the cavern dimensions, the nine excavation levels were combined into three layers in the actual project, as shown in Fig. 4. Each layer was excavated in 17 steps, each with an excavation length of 40 mm (equivalent to 12 m in the prototype). When the excavation of one footage was completed, the excavation was paused and the data was observed. When the data monitored by the sensors is stable, the next step of excavation is started until the entire excavation is completed. Once the cavern excavation is complete, the load is held constant until the monitored data is stable and then recorded.

Model Test Results
As shown in Fig. 5, the sidewall of cavern is seriously damaged. There are significant longitudinal splitting cracks on both sides of the model cavern, but no obvious damage to the vault. This is a typical damage of splitting failure in high sidewall cavern. The length of the cracks near the high sidewalls is slightly less than the height of the sidewalls. The farther away from the sidewall, the shorter the cracks are. The distance from the outermost crack to the sidewall is about 20.8 m, which is consistent with the measured data at the project site.
In this section, all the model test results have been converted to the prototype according to the principles of similitude. After excavation, the deformation of the cavern perimeter all developed toward the cavern (Fig. 6). The displacement of the nearest measurement Cementing agent as a percentage of aggregate weight (%) Iron powder: barite powder: Silica sand 6 5.5 1:1:0.5 Fig. 2 Layout of monitoring sections a and monitoring sensors b (Unit: mm). Section I: the radial stress monitoring section; Section II: the radial deformation monitoring section; Section III: the tangential stress monitoring section point is relatively large, indicating that the area close to the cavern wall is a fractured area of the surrounding rock in the traditional sense, which is consistent with the damage of the surrounding rock observed in the test. The displacement of the measurement point farthest from the cavern sidewall is much smaller than  that of other measurement points on the line, indicating that the excavation disturbance in the high sidewall cavern at 3L from the cavern has almost no effect. The displacements at the sidewall show oscillatory attenuation, which is completely different from that around the shallow cavern. The displacements at the vaults all show a monotonic decay trend, which is completely different from that at the sidewall. After excavation, due to the stress release around the high sidewall cavern chamber, the deformation was all toward the cavern, and the radial stress around the cavern decreased (Fig. 7). With the increase of the distance from the sidewall of cavern, the damage of surrounding rock decreases. The change in radial stress is consistent with the radial displacement. The tangential stress variation shown in Fig. 8 is relatively different from the displacement and radial stress variation. The tangential stresses here are close to zero due to severe damage near the high sidewall resulting in stress relief. Then the tangential stress increases gradually with distance until it reaches a maximum at L. Since the maximum tangential stress shifts from the periphery of the cavern wall to the junction of the elastic-plastic zone when surrounding rock enters the plastic state. It means that the surrounding rock in the range of L is in a plastic state. The surrounding rock in the range of L to 3L is in the elastic state, and the fluctuation of oscillatory attenuation is small, and the variation between two adjacent points is about 15%.

Virtual Working Principle and Control Equation
In Toupin-Mindlin strain gradient theory, in order to describe the state of a material, in addition to Euler strain tensor e ij and Cauchy stress tensor r ij , strain gradient term g ijk and the higher order stress tensor s ijk are also needed (Zhang et al. 2017). The expressions of Euler strain and strain gradient are as follows: In the above equation, both the Eulerian strain e ij and the strain gradient g ijk are symmetric tensors, as are the corresponding Cauchy stress r ij and higher order stress s ijk . It is further assumed that the strain rate and strain gradient rate can be decomposed into elastic part and plastic part respectively: Similar to the strain gradient plasticity model proposed by previous scholars (Mentzel and Steinmann 2000;Fleck and Hutchinson 2001;Gurtin 2000Gurtin , 2002Gurtin , 2003Gudmundson 2004), it is assumed that the plastic strain gradient contributes to the work per unit volume. Therefore, the principle of virtual work should include the work contribution of elastic strain and the work done by plastic strain and its gradient. The virtual work dW i of internal unit volume V can be expressed as: where q ij is called the micro-stress tensor, which is conjugate to the work of plastic strain e p ij ; s ijk denotes the higher-order stress tensor, which is conjugate to the work of the plastic strain gradient e p ij;k . It should be Fig. 7 Variation of the radial stress around the cavern noted that only the deviatoric parts of q ij and s ijk contributes to the virtual work dW i . Therefore, q ij and s ijk in this paper are the omission of q 0 ij and s 0 ijk respectively. The Eq.
(3) can be expressed as follows: Then the following formula can be obtained by applying Gauss' theorem: where n i is the normal vector outward in the direction of the surface S. Equation (3) should be valid for any change in displacement and plastic strain. The first integral on the right side of Eq. (5) can be identified as external virtual work. The second integral on the right side of Eq. (5) should be vanished with arbitrary variations, so two sets of equilibrium equation can be determined: Then, the corresponding conventional T i ¼ r ij n j and higher order t ij ¼ s ijk n k boundary conditions can be obtained according to the principle of virtual work.

Constitutive Equation
In a thermodynamically consistent manner, the dissipative work generated by higher-order strains is introduced by employing higher-order stress related to strain increments (Gudmundson 2004;Gurtin and Anand 2005). The key step in constructing such a constitutive equation is to define the effective stress R work conjugated with the gradient-enhanced effective plastic strain rate _ E p to ensure the plastic work rate: Since the effective stress work R _ E p is always positive, the effective plastic strain rate can be defined as: where L is the dissipative length scale. Thus, the effective stress can be defined as: Then the dissipative plastic micro-stress tensor and dissipative plastic high-order stress tensor can be easily obtained: According to the first law of thermodynamics, when there is damage present in the material, the Helmholtz free energy function is defined with the dissipative potential function as: where C ed ijkl is elastic damage tensor and is a fourth-order tensor, C ed ijkl = 1 À d ð ÞC e ijkl ; K ijklmn is elastic damage tensor considering strain gradients and is a six-order tensor, K ijklmn = 1 À d ð Þl 2 C e ijlm d kn ; C pd ijkl is plastic damage tensor and is a fourth-order tensor, ; l is the internal length parameter of the material, which is closely related to the microcracks and microdefects inside the material and is an inherent property of the material; d kn is a Kronecker symbol. In the above parameters, D is the modified damage variable, D ¼ Dd, D is the damage variable correction parameter or initial damage number; d is the damage variable. In this paper, the damage evolution law is selected to describe the mechanical properties of rock under high in-suit stress. It is considered that there is no damage in the material of elastic state, and damage begins to appear in the plastic state when the rock yield. The expression of damage evolution equation is as follows: where j i indicates the threshold for damage initiation; a and b are material parameters. The C e ijlm described above is the elastic tensor of an isotropic body, and its matrix expression is: According to the incremental theory in the classical elastoplastic theory, the plastic tensor C p can take the following form: where H is the hardening function, which is a function of plastic strain. In order to reflect the characteristics of plastic hardening and plastic softening of rock materials (Fig. 9), its expression form takes the following exponential form: where k i (i = 1,2,3) is the material parameter; I 1 is the first stress invariant, which can represent the average stress or hydrostatic stress. The f and g in Eq. (14) are the yield function and the plastic potential function, respectively. Considering the influence of the intermediate principal stress, the yield function f and plastic potential function g constructed based on the Drucker-Prager strength criterion are as follows: where J 2 is the second invariant of stress deviation; a and b are intensity parameters, and their expressions are as follows: where u and w represent the friction angle and expansion angle of the rock, respectively. When u ¼ w, that is, a ¼ b, it is the associated flow law. According to the second law of thermodynamics, the stress tensor and higher-order stress tensor can be derived from the free energy function strain and strain gradient respectively: 4 Energy Damage Failure Criterion

Rock Destruction and Energy Dissipation
The deformation damage of rocks is essentially a damage evolution process of energy dissipation, so it will be easier to reflect the characteristics of the mechanical response of rocks by studying the relationship between the energy change law and rock damage during the damage process. The damage of rocks is a process of progressive deterioration of materials due to the gradual increase of damage, and microscopic cracks are generated and expanded until penetration. Therefore, the study of rock rupture from the perspective of microscopic mechanics, the analysis of the damage law of tiny units of rocks, and the study of unit damage through the unit energy dissipation and energy discrimination criteria can show the whole process of rock rupture in a more systematic and complete way. Rocks exhibit significant post-peak strain softening under high stress conditions, but a criterion for determining whether damage has occurred in rocks needs to be given. Based on the damage characteristics of the material and the different loading conditions, many different damage criteria have been proposed, Fig. 9 Stress-strain curve of uniaxial compression (Gao et al. 2019) which can generally be divided into four main categories: stress or strain-based damage criteria, energy-based damage criteria, damage-based damage criteria and empirical damage criteria. In essence, the deformation and damage of rocks is a process of continuous energy dissipation, therefore, if the energy change in the process of rock damage can be studied in depth, the strength criterion based on the energy change will be able to describe the process and mechanism of rock damage more accurately.

Element Failure Criterion
In order to obtain the distribution of failure zone in the surrounding rock, it is necessary to give a criterion to determine the failure of element. In the numerical simulation, two criteria are used to determine the failure of elements: maximum tensile strain criterion and energy damage failure criterion. The maximum tensile strain of the unit is calculated according to Eq. (20). If the maximum tensile strain e max is greater than or equal to the ultimate tensile strain e tu , the element is subjected to tensile failure. In order to maintain the integrity and continuity of the entire calculation during the numerical simulation, a small residual elastic modulus E C (E C ¼ 0:05E, E is elastic modulus) is applied to the rock mass element where the tensile failure occurs.
where e x , e y and e z are strains in X, Y and Z directions respectively. If the maximum tensile strain of the element is lower than the ultimate tensile strain (e max \e tu ), the unit will not undergo tensile failure. In this case, the energy damage failure criterion based on strain gradient is adopted as the criterion of element failure. Considering the influence of strain gradient term, the formula for calculating the strain energy density of the element is proposed as follows: The calculation formulas of elastic strain energy density dW=dV ð Þ s , critical strain energy density dW=dV ð Þ c and ultimate strain energy density dW=dV ð Þ u are as follows: The energy damage failure criterion is determined by the change in strain energy density dW=dV ð Þ . When dW=dV ð Þ \ dW=dV ð Þ s ¼ R e s 0 r ij de ij (e s is the yield strain in the process of loading), the material is in the linear elastic phase and no damage occurs inside the element. When dW=dV ð Þ! dW=dV ð Þ s ¼ R e s 0 r ij de ij , the material enters into plastic damage, indicated by the damage variable d, and micro-cracks appear and begin to expand. When dW=dV ð Þ! dW=dV ð Þ c ¼ R e o 0 r ij de ij (e o is the residual strain in the process of loading), the microcracks in the material develop into macroscopic cracks, and the rock mass enters the residual stage, and there is still bearing capacity. When dW=dV ð Þ! dW=dV ð Þ u ¼ R e u 0 r ij de ij (e u is the ultimate strain in the process of loading), the unit is destroyed b Fig. 10 The flow chart of zonal disintegration calculation program and can no longer bear the load. As with the tensile failure, residual elastic modulus Es is imparted to the failure element to maintain the integrity and continuity of the overall calculation.

Development of the Simulation Program
We develop a program based on the finite element method to simulate splitting failure in the surrounding rock of cavern with high sidewall under high stress. The program is developed based on the elastic damage softening model and the energy damage failure criterion which are established from the strain gradient theory and the continuous damage mechanics. In this section, we provide an introduction to the simulation program and compare the simulation results with those obtained from the geo-mechanical model test.
In numerical simulations based on elastic damage softening models, the effect of strain gradient must be taken into account. However, in general finite element analysis, the cell shape function has only first-order continuity and the second-order derivative of the displacement within the cell has a value of 0 after interpolation. This cell analysis does not reflect the effects of strain gradients. Therefore, it is necessary to construct elements with higher order continuity, i.e. elements with C 1 -order continuity, and to require that the nodal displacements and their first-order derivatives maintain continuity. High-order hexahedral cells were constructed using Hermite interpolation functions, and the shape functions and stiffness matrices of the cells were derived (Zhang 2015). This provides a technical basis for considering strain gradient effects in numerical simulations. The simulation program was developed using the finite element software ABAQUS as a platform.
In the elastic-plastic damage softening model, the total strain at any point is decomposes into the Eulerian strain e ij and the higher order strain g ijk , which correspond to the Cauchy stress r ij and the higher order stress s ijk , respectively. In general, the Cauchy stress r ij and Eulerian strain e ij each have 6 independent components, while the higher order stress s ijk and higher order strain g ijk each have 18 independent components, after taking symmetry into account in the 3D condition.
The relationship between Cauchy stress and Eulerian strain is given by: where r ¼ ð r xx r yy r zz r xy r xz r yz Þ T ð27Þ e ¼ ð e xx e yy e zz c xy c xz c yz Þ T where C is a 6 9 6 matrix divided into an elastic matrix C e (Eq. (13)) and a plastic matrix The relationship between higher order stress and higher order strain in the elastic state is given by:  123 g ¼ g xxx g yyx g zzx g xyx g yzx g xzx g xxy g yyy g zzy À g xyy g yzy g xzy g xxz g yyz g zzz g xyz g yzz g xzz Á T where K is an 18 9 18 matrix, K ¼ l 2 D 1 D 2 D 3 D 4 ! . Wherein D 1 , D 2 , D 3 , D 4 as follows: In the above matrix, l is the internal length parameter of the material, n i i ¼ 1 ; 5 ð Þis the elastic constant associated with the strain gradient. According to previous scholars' research (Mindlin 1965;Eshel and Rosenfeld 1975;Yang et al. 2002;Lam et al. 2003), their values are The relationship between higher order stress and higher order strain in the plastic phase can be obtained from Eq. (20) as: Figure 10 shows the simulation program flow chart.

Setup of Numerical Model and Selection of Parameters
The scope of the numerical model is consistent with the scope of the physical model discussed in Sect. 2. As shown in Fig. 11, the numerical model has a scale of 210 m along X, Y and Z. The 3D computational grid consists of 60,000 elements and 67,769 nodes. The external loads applied on the numerical model are shown in Fig. 11. According to the laboratory experiments (Zhang et al. 2021), the initial elastic modulus is E = 41:50 GPa, the Poisson's ratio is t = 0:27, the density rock is q = 26:60 KN m À3 , the compressive strength is r c = 128:80 MPa, the tensile strength is r t = 8:00 MPa, the Lame constant is G ¼ 16:34 GPa, the internal length parameter of the material is l ¼ 0:01 m. The boundary conditions and cavern construction process in the numerical model are consistent with the physical model.

Simulation Results
The above external loads were applied to the numerical model shown in Fig. 11. After the calculations are balanced, the displacements and strains at each node in the model are set to 0, thus establishing the initial insuit stress field of the model. The cavern was excavated to obtain the displacement, strain and stress fields of the cavern envelope. The calculated results were transformed into a column coordinate system to obtain the distribution of radial displacement, radial strain and radial stress of surrounding rock (Fig. 12). As shown in Fig. 12, spaced rupture zones occur in the surrounding rock. The rupture zones are distributed on both sides of the cavern and parallel to the cavern sidewalls. Three rupture zones appear around the cavern chamber in terms of displacement, strain and stress. Near the rupture zones, deformation is concentrated in the middle of the cavern sidewalls and stresses are parallel to the height of the cavern sidewalls. This is consistent with Fig. 5, which shows a vertical fracture parallel to the sidewall in the vicinity of the cavern. However, the damage to surrounding rock from the model test is greater than that from the numerical simulation. This may be due to the fact that the intrinsic model used in the numerical calculations does not take into account the damage effect of excavation on the surrounding rock.
In order to facilitate the comparative analysis between the numerical simulation results and the model test results, the radial displacement data of the numerical model at the corresponding locations of the test model monitoring points were read out, as shown in Table 3. A comparison of the model test and radial displacement numerical simulation results is shown in Fig. 13 for Line A as an example. From the numerical model, it can be seen that the radial displacement of the surrounding rock shows an oscillatory pattern, n 1 þ n 4 þ n 5 0 n 2 n 2 þ n 5 n 1 þ n 2 0 0 0 0 0 n 4 0 0 0 n 5 0 0 0 2 n 2 0 n 3 þ n 4 n 3 2n 2 þ n 3 0 0 0 0 2 n 2 þ n 5 0 n 3 n 3 þ n 4 2n 2 þ n 3 0 0 0 0 2 n 1 þ 2n 2 0 2n 2 þ n 3 2n 2 þ n 3 2n 1 þ 4n 2 þ n 3 þ n 4 þ 2n which is characterised by alternating peaks and troughs, with the amplitude of the oscillations decreasing with distance from the wall. These characteristics are consistent with the radial displacement variation pattern measured by the model test. The reliability of the elastic damage softening model and numerical analysis method based on the ABAQUS platform is verified. The results of the numerical simulation show that the excavation unloading leads to the rupture of the surrounding rock. The stress waves generated by the rock rupture propagate outwards from the cavern walls and form a new stress field in the cavern envelope. The rock, after superimposed with the original rock stress field, may rupture again, forming new rupture and non-ruptured zones.

Conclusions
Splitting failure in deeply buried high sidewall cavern is a specific and regular strain localization phenomenon. The study of splitting failure should start with the strain localization and strain softening of the rock behind the peak. The effect of strain gradients should also be considered.
Based on strain gradient theory and continuum damage mechanics, an elastic-plastic damage softening model for splitting failure was developed. The relationship between rock damage and energy dissipation was analyzed. Based on the strain energy density theory, a strain gradient-based splitting failure criterion was established.
To account for strain gradient effects, an eight-node high-order hexahedral cell was constructed. The shape functions and stiffness matrices of the higher order cells are derived. A numerical analysis method for splitting failure is proposed. A splitting failure calculation program based on ABAQUS is also developed. The numerical simulation results based on ABAQUS and 3D geomechanical model tests are in general agreement. The reliability of the splitting failure elastic-plastic damage softening model and the numerical analysis method of splitting failure were effectively verified.  Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.