Influence of Incident Orientation on the Dynamic Response of Deep U-Shaped Cavern Subjected to Transient Loading

: In deep rock engineering, caverns are often disturbed by engineering loads from different directions. To investigate the dynamic response of deep U-shaped caverns under different incident orientations, a theoretical solution of the dynamic stress concentration factor along the cavern boundary was derived based on the wave function expansion and conformal mapping method, and the failure characteristics around the cavern were further investigated by PFC2D (Particle Flow Code in two dimensions). As the incident orientation increases from 0° to 90°, the dynamic compressive stress concentration area transforms from both the roof and the floor to the sidewalls, and the peak dynamic stress concentration factor of the roof decreases from 2.98 to −0.20. The failure of the floor converts from dynamic compression shear failure to dynamic tensile failure. Compared to a stress wave incident from the curved boundary, a stress wave incident from the flat boundary causes severer damage. When the stress wave is incident from the sidewall, the cavern with a larger height-to-width ( h / w ) ratio exhibits severer damage. Conversely, the cavern with a smaller h / w ratio tends to fail as the stress wave is incident from the floor. This paper provides a basic understanding of dynamic responses of the deep U-shaped cavern.


Introduction
With the development of underground space and resources, a large number of deep caverns have been established [1][2][3].During the excavation and construction of deep caverns, engineering disasters such as rockburst, spalling, and large deformation occur frequently due to the influence of complex geological environments and stress conditions [4][5][6][7][8][9].These engineering disasters pose a serious threat to the stability of related underground engineering, as well as the safety of personnel and equipment.Studying the disaster mechanism of deep rock engineering is of great significance.
The occurrence of engineering disasters is influenced by various factors, and many scholars have researched these factors extensively.Carter [10] and Zhao et al. [11] conducted laboratory experiments to study the influence of circular opening size on the failure strength around cavities and found that the failure strength of the opening decreased as the opening size increased.Huang et al. [12] compared the failure process of a Dshaped tunnel under different maximum principal stress directions.They found that the failure position of the surrounding rock shifted from the sidewall to the spandrel and finally to the floor as the maximum principal stress direction shifted from vertical to horizontal, and the size of rock fragments around the opening decreased first and then increased.Pan et al. [13] conducted laboratory experiments to investigate the failure modes of tunnels caused by tangential stress concentration.The results revealed that the presence of an open fracture around the tunnel resulted in more extensive failure at the tunnel boundary compared to the tunnel without a nearby fracture.By numerical simulation, Zeng et al. [14] studied the mechanical properties of sandstone specimens containing a hole under uniaxial compression.They found that the shape of the opening significantly affects the stress field and failure mode around the openings.
Aside from the high in situ stress level, the stability of the cavern is also significantly affected by various dynamic disturbances [15][16][17].Kouretzis et al. [18] compared the dynamic responses around the circular cavern subjected to P-waves and S-waves through a theoretical analysis and found that P-waves cause greater tangential stress at the cavern boundary than S-waves.Yi et al. [19] obtained the distribution characteristics of the dynamic stress concentration factor around a lined cavern based on theoretical calculations and pointed out that the incident wave frequency, bonding parameters, and distance of disturbance sources have a significant influence on the dynamic responses around the cavern.Zhao et al. [20] investigated the dynamic stress distribution around a circular cavern by theoretical analysis and numerical simulation.Their findings showed that the existence of the excavation damage zone reduces the dynamic stress concentration around the cavern, and the dynamic stress concentration factor around the intact surrounding rock decreases with the increasing thickness of the damage zone.Saiang [21] analyzed the blast-induced damage zone around the cavern by coupled continuum-discontinuum methods.Their results indicated that the blast-induced damage zone around the cavern at shallow depths is more discrete, while the deep cavern experiences regional failure.
This research has made a great contribution to our understanding of the failure mechanism of deep caverns under dynamic disturbances.In practical engineering, the U-shaped cross-section is a common option, and the dynamic disturbances may derive from various directions of the cavern [22,23].However, there are few studies investigating the dynamic response of uncircular caverns subjected to transient loading, and the failure mechanism around the deep-buried U-shaped cavern under different incident orientations remains unclear.To further investigate the influence of disturbance incident orientation on the transient response of deep U-shaped caverns, the theoretical solution of dynamic stress concentration factors along the cavern boundary under different incident orientations was obtained by the wave function expansion and conformal mapping methods.Furthermore, the failure characteristics and mechanism of deep U-shaped caverns under different incident angles were systematically studied.This work provides insight into failure mechanisms of deep U-shaped caverns subjected to dynamic disturbances with different incident orientations, as well as guidance for the support design of deep rock engineering.

Theoretical Investigation for the Dynamic Responses of U-Shaped Caverns under Different Orientations
As shown in Figure 1, it is assumed that a deep-buried U-shaped cavern located in the homogenous elastic medium and subjected to a transient loading, where h and w are the height and width of the cavern, and a = w/2 is the radius of the semi-circle arch.The incidence angle is θ0, and the vertical and horizontal in situ stresses are σv and σh, respectively.According to the superposition principle of elasticity theory [24][25][26], the dynamic responses of a deep-buried cavern under transient loading can be divided into the static part induced by in situ stress (Figure 1b) and the dynamic part induced by transient loading (Figure 1c).For one-dimensional problems, non-periodic disturbance can be regarded as the boundary condition, and the wave equations can be solved by the traveling wave method, variable separating method, or integral transform method [27].However, there are a lot of mathematical difficulties in directly solving the wave equations of two-dimensional problems with complex boundary conditions.Comparatively, the wave function expansion method is more suitable for solving two-dimensional wave equations.In this paper, the conformal mapping method is employed to obtain the stress field around the U-shaped cavern.To simplify the calculation, the cavern is symmetrical about the x-axis.
As presented in Figure 1d, the exterior of the cavern section in the z-plane is transformed to the exterior of a unit circle in the ζ-plane by using the following mapping function: The static stress field of a U-shaped cavern under in situ stress can be obtained by complex variable theory [28,29].This section aims to obtain the theoretical solution of dynamic responses of the U-shaped cavern subjected to transient stress waves (Figure 1c).

The Steady-State Responses of U-Shaped Caverns under Harmonic Stress Waves
The transient response of a cavern under dynamic disturbance can be calculated by superposing the steady responses of harmonics of all frequencies.Thus, the dynamic response of a U-shaped cavern subjected to harmonic P-wave was investigated first.
The incident harmonic P-wave in the ζ-plane can be described as where φ0 is the amplitude of incident wave, α is the P-wave number, () W  is the conjugate of the mapping function ( ) W  , θ1 is the angle between incident direction of the Pwave and x-axis, , ω is the circular frequency of the incident stress wave, and t is the time.
When the incident P-wave propagates to the boundary of the cavern, the reflected waves, P-wave , the reflected waves satisfy the wave equation [30,31]: where β is the S-wave number, and ( ) are the first derivative of the mapping function and its conjugate, respectively.To satisfy the Sommerfeld radiation condition, the reflected waves can be written as (1) where Hn (1) is the first kind of Hankel function, and an and bn are unknown coefficients.The total wave field in the surrounding rock can be expressed as In the ζ-plane, the radial stress r  , tangential stress   , and shear stress r  around the U-shaped cavern can be written as where λ and μ are the Lame constants of the surrounding rock.Since there is no external force on the boundary of the cavern, the corresponding boundary condition is Substituting Equations ( 9) and ( 11) into Equation ( 12) and making use of the orthogonality of is e  ( ) , the undetermined coefficients an and bn can be determined by solving a set of equations [32]: The elements in these equations can be found in Appendix A. This paper focuses on the stress distribution on a cavern boundary where only tangential stress exists.Thus, the dynamic stress concentration factor is defined as where is the stress amplitude of the incident P-wave.In this section, the geometric dimensions of the U-shaped cavern are h = w = 2.5 m.Based on the conformal mapping method, the corresponding mapping function which has been calculated by our previous work can be expressed by the following equation [33]: ( ) 1.348 0.141 0.008 0.114 0.108 0.038 0.010 Taking the red sandstone in Linyi, Shandong, China, as an example, the density, elastic modulus, and Poisson's ratio of the surrounding rock are 2440 kg/m 3 , 18.60 GPa, and 0.21, respectively, and the Lame constants are λ = 5.57GPa and μ = 7.69 GPa [33].With the different incident orientations (θ0 = 0°, 30°, 45°, 60°, 90°, and 270°) and wave numbers (αa = 0.1, 1.0, and 5.0), the dynamic stress concentration factors along the U-shaped cavern boundary under harmonic P-waves are illustrated in Figure 2. When the stress wave is horizontally incident from the left side of the cavern (θ0 = 0°), the stress concentrations occur at the roof and floor of the cavern.With increasing incident angle (θ0 = 30°, 45°, and 60°), the stress concentration zones transform from the roof and floor to the left spandrel and right bottom corner.When the stress wave is vertically incident (θ0 = 90° and 270°), the stress concentrations can be observed around both sides of the cavern.As αa = 0.1, the stress distribution of θ0 = 90° is identical to that of θ0 = 270°.With increasing wave number, there are significant differences in the dynamic stress concentration factors between θ0 = 90° and 270°, which indicates that the geometry of the incident boundary would affect the dynamic stress distribution around the cavern.
Comparing the dynamic responses with different wave numbers, the dynamic stress concentration factors of αa = 0.1 are larger than those of αa = 1.0 and 5.0.When the wave number αa approaches zero, the frequency of the incident wave is quite small, and the surrounding rock of the cavern is approximately under a biaxial static loading state, as shown in Figure 3a, where  is the Poisson's ratio of the surrounding rock.Figure 3b demonstrates the static stress concentration factors along the cavern boundary based on complex variable theory.With an increase in θ0 from 0° to 90°, the stress concentration factor at the roof decreases from 3.1 to −0.1, and the factor at the midpoint of floor decreases from 1.41 to −0.54, while the factor at the midpoint of the left sidewalls increases from −0.4 to 2.0.As shown in Figures 2 and 3, the stress distribution of αa = 0.1 agrees with the stress distribution under static biaxial compression, which validates the accuracy of the theoretical method of this section.

The Transient Response of U-Shaped Caverns under Different Orientations
Based on the steady-state response of a harmonic wave, the transient response of the U-shaped cavern can be obtained by Fourier transform: where p(ω) is the Fourier transform of the incident stress wave P(t), and ( )  is the admittance function of the harmonic wave.
In this paper, a blasting load which is simplified as a triangular wave is chosen as the incident stress wave: where Pd is the amplitude of the incident stress wave, tr is the rising time, and ts is the total time.
The corresponding dynamic response of the U-shaped can be expressed as [34] ( ) ( ) where R(ω) is the real part of the admittance function.
With the different incident orientations (θ0 = 0°, 45°, and 90°, ts/tr = 5), the dynamic stress concentration factor evolutions at the midpoint of the roof, floor, and left sidewall are illustrated in Figure 4.When θ0 = 0° and tr = 4 ms, the peak dynamic stress concentration factor of the cavern roof is 2.98.The shorter rising time of the incident wave, the peak of the roof becomes smaller.With increasing incident angle, the peak dynamic stress concentration factor of the roof decreases.When θ0 = 90° and tr = 4 ms, the negative peak dynamic stress concentration factor of the cavern roof is −0.20.This indicates that the cavern roof converts from dynamic compressive stress concentration to dynamic tensile stress concentration with the incident angle increasing from 0° to 90°.The tangential stress of the cavern floor exhibits the same evolution process.However, the peak values of the dynamic stress concentration factor of the floor are smaller than that of the roof.In contrast, the dynamic stress concentration factor of the left sidewall increases with increasing incident angle.When θ0 = 0°, the tensile stress concentration occurs on the left sidewall.It transforms to compressive stress concentration as θ0 increases to 45° and 90°.The incident orientation affects the distribution characteristics of dynamic stress concentration areas around the cavern.Comparing the peak value of dynamic stress concentration factor among the roof, floor, and left sidewall under different incident angles, the maximum value occurs at the cavern roof when θ0 = 0°, tr = 4 ms, while the minimum value occurs at the cavern floor when θ0 = 90°, tr = 4 ms.It indicates that the compressive stress concentration at the curved surface (roof) is the most obvious, while the tensile stress concentration at the flat surface (floor) is the most evident.

Numerical Model
The particle flow code (PFC), which is based on the discrete element method, has been widely used to model failure behavior in rock engineering [35,36].Compared with other numerical software, the complex fracture criterion is not required, and the crack evolution of rock materials can be observed directly in PFC [37,38].The dynamic stress concentration factors along the cavern boundary under different incident angles were calculated in the above section by theoretical analysis, without considering the influence of in situ stress.In addition, the theoretical solution is restricted to the elastic stage of the surrounding rock.To further investigate the dynamic responses and the failure characteristics around the deep U-shaped cavern under transient loading, a two-dimensional numerical model was established by PFC2D.As presented in Figure 5, a 15 × 15 m model was constructed in this section, and a U-shaped cavern with the geometric dimensions h = w = 2.5 m was set up at the center of the model.The model has approximately 30000 particles with radii ranging from 0.03 to 0.06 m.The flat-joint model was used in this paper to better simulate the mechanical and cracking behaviors of natural rock [39,40].According to the mechanical parameters of red sandstone, the corresponding microscopic parameters of the particle and the flat-joint model are listed in Table 1, which have been calibrated in our previous work [33].The particle boundary was used in the model, and the viscous boundary was set up to prevent the reflected wave emerging from the boundary [41], as shown by the green particles in Figure 5.The buried depth of the cavern is 1000 m, and the corresponding vertical stress is 26.1 MPa according to the stress field in China [42].Firstly, the confining stresses are applied on the model boundary to simulate the in situ stress field.The loading rate of the in situ stress is 5.5 × 10 −4 MPa/step, and the local damping coefficient is 0.7 for static stress loading.Maintaining the in situ stress, a triangular plane P-wave is applied on the left side of the model subsequently.To observe the dynamic response of the surrounding rock, the local damping coefficient of the model is modified to 0.0.

Failure Characteristics of U-Shaped Caverns under Different Incident Orientations
The theoretical solution points out that the incident orientation of the stress wave affects the dynamic stress distribution of the surrounding rock, leading to different failure modes of the U-shaped cavern.Figure 6 presents the failure modes of the deep-buried Ushaped cavern under different incident angles (θ0 = 0°, 30°, 45°, 60°, 90°, and 270°), where the blue segments are shear microcracks and the red segments are tensile microcracks.In this section, the vertical stress and horizontal stress are σv = σh = 26.1 MPa, and the amplitude, the rising time, and the ratio of total time to rising time of the incident stress wave are Pd = 60 MPa, tr = 2 ms, and ts/tr = 5, respectively.Figure 6 presents the failure modes around the U-shaped cavern under different incident orientations.When the stress wave is horizontally incident (θ0 = 0°), the theoretical solution suggests that the compressive stress concentrations occur at the roof and floor of the U-shaped cavern.Therefore, dynamic compressive-shear failure occurs around the roof and floor.As shown in Figure 6a, the microcrack density around the floor is much higher than that of the roof, resulting in the emergence of a larger damage zone around the floor.In addition, spallation occurs around the incident sidewall due to the reflection of the stress wave.Compared with the damage zone around the floor, the damage range around the sidewall is much smaller.
When the stress waves are obliquely incident, the dynamic compressive stress concentration area transfers to the vicinity of the left spandrel and right bottom corner.The failure mainly occurs on the cavern floor.When θ0 = 30°, the number of microcracks of the surrounding rock decreases compared to θ0 = 0°.There are no macro failures on the left sidewall and roof, and the damaged area around the cavern floor decreases also.With the increasing incident angle (θ0 = 45° and 60°), the number of microcracks and the damage depth near the cavern floor increases.When θ0 = 60°, the dynamic compressive stress concentration results in fragment exfoliation at the left spandrel of the cavern.
When the stress waves are vertically incident, the geometric differences between the roof and floor of the cavern contribute to the different failure modes with θ0 = 90° and 270°.Figure 6e presents the failure mode of the cavern as the stress wave is incident from the cavern floor.Due to the reflection of the stress wave at the cavern floor, a larger number of tensile microcracks is generated, and the spallation of a wedge-shaped slab can be observed near the floor.In addition, the dynamic compressive-shear failure occurs around the sidewall.However, when the stress wave is incident from the cavern roof, the number of microcracks of the model decreases, and only a small fragment separates from the incident side, as shown in Figure 6f.Moreover, the dynamic compressive shear failures both occur on the sidewall when θ0 = 90° and 270°, but the damage extent of θ0 = 90° is larger than that of θ0 = 270°.
Comparing the damage extent of the surrounding rock with different incident angles, it was found that the order of the total number of microcracks from most to least is as follows: 90° > 0° > 60° > 45° > 270° > 30°.Moreover, Zaid et al. [43] compared the failure patterns of different cavern shapes under blast loading and found that the damage extent of rectangular caverns is more severe than that of circular caverns.This indicates that stress waves incident from the flat boundary cause more severe damage around the cavern than that from the curved boundary.Accordingly, the flat boundary should be focused on in practical engineering, and the lining structure with the inverted arch is suggested to reduce the damage degree of the surrounding rock.
As the incident angle of stress waves increases from 0° to 90°, the depth of the failure zone around the cavern floor decreases first and then increases.To further investigate the failure mechanism of the floor area under different incident angles, the evolution of tangential stress and kinetic energy at the midpoint of the floor are illustrated in Figure 7.With increasing incident angle, the peak value of tangential stress at the midpoint of the floor decreases, which means that the dynamic compressive stress concentration around the floor reduces.This is consistent with the theoretical solutions in the above section.Therefore, the damage depth around the floor decreases as the incident angle increases from 0° to 30°.When the stress wave is obliquely incident, the tangential compression and radial tension both occur around the cavern floor.Figure 7b presents the evolution of kinetic energy of the midpoint of the floor under different incident angles.With increasing incident angle, the peak value of kinetic energy increases first and then converges to a constant value.In contrast, the kinetic energy contributed by the tangential compression decreases with increasing incident angle.It indicates that the intensity of the reflected tensile stress along the radial direction increases with the increasing incident angle.This phenomenon is consistent with the theoretical results of Xu et al. [44].They also pointed out that the amplitude of reflected longitudinal waves increases with the increasing angle between the direction of stress wave incidence and the incident surface.Consequently, as the incident angle increases from 30° to 90°, the thickness of the separated rock slab increases also.

Discussion
In practical engineering, the height-to-width (h/w) ratio of a U-shaped cavern varies with the different spatial requirements and the complex geological conditions.To investigate the influence of h/w ratio on the stress distribution and the failure characteristics around a U-shaped cavern, three h/w ratios (h/w = 1.0, 1.5 and 2.0) were specified in this section.The width of the caverns remains constant (w = 2.5 m), while the heights of the caverns are h = 2.5, 3.75, and 5.0 m, respectively.The vertical and horizontal in situ stress are σv = σh = 26.1 MPa, and the amplitude, the rising time, and total time of the incident stress wave are Pd = 45 MPa, tr =2 ms, and ts =10 ms, respectively.Figure 8 presents the failure modes of the deep-buried U-shaped cavern with different h/w ratios under dynamic disturbance.When the stress wave is horizontally incident (θ0 = 0°), the number of tensile microcracks and the thickness of the separated rock slab around the incident sidewall increase with the increasing h/w ratio.Due to the dynamic compressive stress concentration, particle ejection occurs around the cavern floor as h/w = 1.5 and 2.0.When the stress wave is incident from the cavern floor (θ0 = 90°), the number of tensile microcracks around the floor (incident side) decreases with increasing h/w ratio, and spalling occurs at the floor when h/w = 1.0.Moreover, the tensile and shear microcracks generate around the sidewall.With the larger h/w ratio, the microcracks around the sidewall extend to deeper areas of the surrounding rock.Figure 9 presents the evolution of the tangential stress and kinetic energy of the sidewall and floor.When θ0 = 0°, the peak value of the kinetic energy of the midpoint of the left sidewall increases with the increasing h/w ratio.It indicates that the intensity of the reflected tensile stress wave around the incident sidewall increases with the increasing h/w ratio.Thus, the separated slab in the sidewall also increases and is accompanied by the release of more residual kinetic energy.Before the stress wave propagates to the boundary of the cavern, the static tangential stress at the midpoint of the floor increases with the increasing h/w ratio.Therefore, a U-shaped cavern with a larger h/w ratio exhibits more considerable peak tangential stress at the floor.Since the size of the incident side (floor) of the cavern remains constant, the peak kinetic energy at the incident side is generally consistent as the stress wave is incident from the floor.However, spalling only occurs in the floor with h/w = 1.0.When h/w = 1.5 and 2.0, the cavern floor maintains stability after the stress wave passes through the cavern.It indicates that the ultimate strength of the floor increases with an increasing h/w ratio.Moreover, the peak value of the tangential stress at the midpoint of the sidewall decreases with the increasing h/w ratio, as shown in Figure 9d.When h/w = 1.0, the peak tangential stress of the sidewall exceeds the dynamic compressive strength of the surrounding rock, resulting in failure of the sidewall.
Current research mainly concentrates on the failure characteristics of caverns with a single h/w ratio [12,45].In fact, there are significant differences in failure patterns between caverns with different h/w ratios.Comparing the damage characteristics around the Ushaped cavern under different incident angles, the cavern with the larger h/w ratio exhibits severer damage when the stress wave is incident from the sidewall.On the contrary, when the stress wave is incident from the cavern floor, the damage extent of a cavern with a smaller h/w ratio is more notable.To reduce the damage of the surrounding rock, the disturbances incident from the high sidewall need to be minimized as far as possible.

Conclusions
In this paper, based on the theory of elasto-dynamics and conformal mapping, the dynamic stress concentration factor around the U-shaped cavern under transient stress waves with different incident orientations was calculated.Furthermore, a numerical model was established by using PFC2D code to investigate the corresponding failure characteristics.Conclusions could be drawn as follows: (1) With increasing incident orientation of the stress wave from 0° to 90°, the dynamic compressive stress concentration area transforms from both the roof and floor to the sidewall, while the dynamic tensile stress concentration area transforms from the sidewall to both the roof and floor.When the wave number approaches zero, the surrounding rock of the cavern is approximately under a biaxial static loading state.The dynamic stress concentration factor of αa = 0.1 is more significant than those of αa = 1.0 and 5.0.(2) With increasing incident orientation of the stress wave from 0° to 90°, the failure of the floor changes from compression shear failure to tensile failure.Compared to a stress wave incident from the curved boundary, a stress wave incident from the flat boundary causes severer damage around the cavern.The most serious damage occurs around the cavern when the dynamic disturbance is incident from the cavern floor.(3) When the incident orientation of the stress wave is 0° (horizontal direction), the peak kinetic energy of the incident sidewall and the peak dynamic compressive stress of the floor increase with increasing h/w ratio of the cavern.A cavern with a larger h/w ratio is more prone to failure, showing more residual kinetic energy.When the stress wave is incident from the floor, the dynamic responses of the cavern present a reverse trend with increasing h/w ratio.(4) The theoretical solutions based on wave function expansion match the modelling results well.In a real-world project, the theoretical method presented in this paper is possibly beneficial to predict the potential risk area of the surrounding rock.
Author Contributions: Methodology, software, and writing, L.L.; validation and funding acquisition, X.L.; writing-review and editing, Z.L.; software, S.P.All authors have read and agreed to the published version of the manuscript.
Funding: This research was funded by the National Natural science Foundation of China, grant number 51927808.
Data Availability Statement: Data will be made available on request.

Conflicts of Interest:
The authors declare no conflicts of interest.

Figure 1 .
Figure 1.Theoretical model of a deep-buried U-shaped cavern under the P-wave.

Figure 3 .
Figure 3. (a) Model of a U-shaped cavern under static biaxial compression; (b) static stress concentration factors along the cavern boundary under different θ0.

Figure 4 .
Figure 4. Theoretical results of dynamic stress concentration factors of a U-shaped cavern under transient loading.

Figure 5 .
Figure 5.The numerical model of a deep U-shaped cavern under transient loading.

Figure 7 .
Figure 7. Dynamic response of the midpoint of the floor with different incident orientations: (a) tangential stress; (b) kinetic energy.