Effect of Seismic Frequency Spectra on Surrounding Rock Damage Evolution of Large Underground Caverns

To explore the distribution and evolution of the surrounding rock of underground caverns of hydropower stations, the dynamic compression tests for granite sample were carried out and the dynamic constitutive relation suitable to computer programming was established in terms of strain-based static elastic-plastic damage constitutive relation. By inputting five earthquakes of various frequency spectra, the fully nonlinear dynamic analysis was performed for some large underground caverns in the gorge and mountain region in Sichuan Province in China with the explicit difference method. +e simulation results show that, as the input seismic waves become strong gradually, the damage zone generally extends deeply from the middle of side walls.+e damage zone reaches the maximum value in areas immediately after the seismic acceleration peak arrives and keeps constant thereafter. According to the initial analysis, if the dominant frequency of input wave approaches the natural frequency of underground caverns, the maximum area of damage zone will increase. +ose findings may provide practical data for earthquake-resistant design and construction.


Introduction
In the past, it was believed that the damage caused by earthquakes was not serious for underground facilities compared to surface structures, and minimal attention was focused on the seismic study of underground facilities until the 1990s when earthquakes occurred in Kobe (Japan, 1995), Chi-Chi (Taiwan, 1999), and Kocaeli (Turkey, 1999), destroying selected underground structures [1].e collapse of the Daikai subway station was the first collapse of an urban underground structure due to earthquake forces rather than ground instability.During the Ms 8.0 Wenchuan earthquake which occurred on May 12, 2008, eleven tunnels of the Dujiangyan to Wenchuan highway suffered different degrees of damage [2], and the safe operation of the two hydropower stations in Yingxiuwan and Gengda was threatened [3].In addition, the surrounding rock of underground caverns of Yingxiuwan Hydropower Station in Southern China deformed seriously during the 2008 Wenchuan earthquake, the lining cracked, and the bottom concrete plate was displaced, which paralyzed the hoisting machine service [4].Hence, it is essential to investigate the seismic dynamic response of the engineering rock mass of underground caverns.ree methods, namely, theoretical study, experimental study, and numerical modelling, are commonly used for studying this subject.e theoretical study does not easily consider complicated geological cases [5][6][7].
e experimental study is sometimes limited by the existing test techniques and scale effects [8,9].e limitations of these two methods motivated the rapid development of the numerical modeling which provides an economical, convenient, and efficient approach.
In recent years, many researchers have investigated the damage mechanism and influence factors of tunnels or caverns under seismic loading [10][11][12][13][14]. Zhang applied the modified DDA method to study the seismic response of the underground houses of the Dagangshan Hydropower Station in Western China, and valuable results were obtained [15].Cui et al. [16][17][18][19] analyzed the seismic performance of a rock joint with a modified continuously yielding model by DEM software and also studied the control effect of a large geological discontinuity on the seismic response and stability of underground rock caverns.Fu et al. [20,21] used the discontinuous deformation analysis method to analyze the seismic wave propagation across a rock mass of a large rock cavern complex.
Furthermore, the reasonable seismic input is vital to the dynamic response of underground rock caverns in the mountainous region.Most earthquakes applied at the numerical or analytic models were obtained directly or indirectly from the seismic records measured on or near the ground.However, some seismic data measured at various depths show that the peak acceleration in deep position is generally less than that on the ground.
e Japan code stipulates that the acceleration at 20 m depth is 1/2-1/3 of that on the ground.Hu and Xie studied the seismic propagation in Hosokura mining region in Japan and found that the average horizontal peak acceleration at about 400 m depth is approximately 0.65 times that on the ground [22].In addition, the mountain, especially the mountain peak, can amplify the acceleration.Fortunately, most actual mountainous seismic records were measured at mountain foot, which is little influenced by the amplification effect.In mountainous region, the so called "ground" is assumed to be at the river surface based on the fact that most seismic records were measured at the mountain foot.us, for many underground caverns of hydropower stations in mountainous region in Southern China, the numerical model bottom is always hundreds of meters deep from the river surface, and the depth attenuation should be considered in order to simulate the dynamic response under seismic loading in a relatively reasonable way.
In addition, there are abundant microfractures, pores, or cracks within surrounding rock of underground caverns, so the growth of the damage zone under seismic loading should be also considered.Based on the static damage accumulation model provided by Zhu et al. [23], a dynamic damage accumulation model is proposed which is assumed to be proportional to the static one.is paper established a 3D numerical simplified model to analyze the damage zone evolution under seismic loadings of various Fourier spectra.

Depth Attenuation in Mountainous Region
e Dagangshan Hydropower Station in mountainous Western China was taken as the case study, and the mountain peak is about 600 m higher than valley bottom.One horizontal seismic record measured in the Panzhihua region in 2008 Wenchuan earthquake was selected as the input earthquake wave for the numerical model, as shown in Figure 1. is seismic record is about 300 km away from the epicentre.Its Fourier spectrum is shown in Figure 2, and the frequencies more than 5 Hz are filtered.
e design peak acceleration for the Dagangshan dam site under the exceedance probability of 3% in 50 years is 4.403 m/s 2 .
To approximately obtain the acceleration attenuation with depth, a simple model with no consideration of caverns was established.e mountain top was averagely flattened to simply simulate the mountain condition.e seismic wave was vertically input at the bottom of the model, with free-field boundary.
e peak acceleration of Panzhihua earthquake is assumed to be 1 m/s 2 for the purpose of simplification, and the damp ratio is 0.05.e FLAC3d was used to simulate the acceleration propagation, and the results are shown in Figure 3. e results show that the horizontal peak acceleration at the bottom is about 72 percent of that at the 2 Advances in Materials Science and Engineering valley, which means the dynamic response will be overestimated if the depth attenuation is neglected.

Damage Accumulation Model
e damage constitutive model is the foundation for evaluating the safety of rock masses that surround underground caverns during an earthquake.Zhu et al. provided a static elastic-plastic damage constitutive model to describe the deformation and failure of brittle rocks such as granite [23].
e model is expressed as where E s is the initial static elastic modulus, I is the unit matrix, and ω s is the static damage third-order tensor.If the damage variables in three principal directions are independent of each other, then the damage evolution equation is where ε i is the principal strain, i � 1, 2, and 3, ε s0 is the initially cracked static strain, and ε sc is the static failure strain.
Assume that the rock strain consists of the elastic strain and crack strain, namely, where ε 1 is the axial strain, ε 3 is the lateral strain, ε e 1 is elastic axial strain, ε e 3 is the elastic lateral strain, ε f 1 is the axial strain due to crack propagation, and ε f 3 is the lateral strain due to crack propagation.ε e 1 and ε e 3 can be obtained by Hooker's law, and the crack strain can be determined by subtracting the crack strain from the elastic strain.
e ω s1 − ε 1 and ω s3 − ε 3 curves under the static test of 10 −5 •s −1 strain ratio for granite drilled from the Dagangshan hydropower site are shown in Figure 4, which indicates that the crack stain mainly comes from the lateral direction.us, it is assumed that the only elastic strain occurs along the ε 1 direction, neglecting the axial brittle damage.In addition, the test results show that the lateral failure strain is not sensitive to the surrounding pressures and loading ratios.
us, a simple and practical failure criterion was established as follows: Assume the dynamic elastic-plastic damage constitutive model is similar to the static one, and thus is written as where E d is the initial dynamic elastic modulus.
According to the assumption above, the dynamic constitutive curve is similar to the static one except that their peak stresses and initial elastic modulus are different, as shown in Figure 5. us, it is inferred that there is some relationship between their damage variables.
Figure 5 shows the geometric relation as us, (BC/AC) � (E s /E d ), and then multiplied by For the purpose of simplification, it is assumed here that the lateral failure strain of bristle rock is not affected by the strain ratio, and the strain amplification coefficients are identical for ε do ≤ ε ≤ ε dc , and hence -0.006 -0.004 -0.002 0 0.002 0.004 0.006 0.008 0.01 Advances in Materials Science and Engineering 3 e relationship among rock peak stress, initial elastic modulus, and strain ratio is established by the dynamic test results of Dagangshan granite.
e threshold stain for dynamic damage is assumed to not be affected by the strain ratio, and then it can be determined by the following equation: e formula above can be used to numerically simulate the damage evolution of Dagangshan underground caverns by using the embedded FISH language of the FLAC 3D software.

Numerical Model and Boundary Conditions
e three caverns of Dagangshan Hydropower Station all have a gate arch shape, with axis direction NE55.A threedimensional model is established for three adjacent underground caverns by means of FLAC 3D software [24], as shown in Figure 6.
e model dimensions are 600 m × 420 m × 200 mm, and the caverns are about 390 to 520 m deep.e top of the model is applied by the corresponding overlying pressure in accordance with the model top depth and rock density.No reinforcement support is considered in this model.e zone around the underground caverns is mainly composed of relatively fresh and intact granite, with several faults.However, if the faults and fractures are considered in the model, it is very hard to deal with the dynamic simulation due to the much time cost.us, the whole rock mass is considered intact rock mass, and the rock property parameters used for numerical simulation are listed in Table 1.
e reflection of seismic waves on the boundaries will heavily affect the dynamic simulation of underground caverns, so it is important to apply proper artificial boundaries around the model.us, the free-field boundary is used, and the seismic waves are applied at the model bottom [25].Considering the depth attenuation, the peak input acceleration is 4.043 × 0.72 � 2.911 m/s 2 under the exceedance probability of 3% in 50 years.

Input Seismic Wave
e underground cavern is often affected by low-frequency compositions of the earthquake [10], so the five earthquake records which contain rich low-frequency components are adopted as the input waves.In order to compare the dynamic response of underground caverns under seismic waves of various frequency spectra, Panzhihua earthquake record, Kobe earthquake, and three artificial earthquakes are selected for numerical simulation.
e Kobe earthquake happened in Japan in 1995.
e other three waves are artificially made by triangle series superposition method.
e Panzhihua earthquake, Kobe earthquake, and three artificial earthquakes differ greatly in acceleration-time history and Fourier spectrum.
A seismic wave is usually a series of discrete acceleration points, and it should be filtered and rectified before used for numerical simulations [26].Higher frequency will need smaller mesh dimension, larger grid number, and longer calculation time under the same computation accuracy condition.In this study, the frequencies larger than 8 Hz are removed from the seismic wave, which greatly reduces the number of grids and improves the computation efficiency.

Results and Discussion
e above five earthquakes were input into the numerical model for simulation separately.For each earthquake, the damage zone evolution of the surrounding rock of underground caverns was carried out under three typical moments which were selected based on the great change in acceleration.

Panzhihua Earthquake.
e distributions of the damage zone of the surrounding rock under Panzhihua earthquake at the 1.522 s moment, 3.862 s moment, and 10.31 s moment   16 shows the extension of damage zones of powerhouse, transformer chamber, and tailrace surge chamber with time.At the 1.522 s moment, the damage zone of powerhouse is 778.1 m 2 in area, that of transformer chamber 235.6 m 2 , and that of tailrace surge chamber 770.9 m 2 , totaling about 1785 m 2 .At the 3.862 s moment, the damage zone of powerhouse is 959.9 m 2 in area, that of transformer chamber 355.9 m 2 , and that of tailrace surge chamber 913.9 m 2 , totaling about 2229 m 2 .At the 10.31 s moment, the damage zone of powerhouse is 1010 m 2 in area, that of transformer chamber 380.1 m 2 , and that of tailrace surge chamber 1032.9 m 2 , totaling about 2422 m 2 .After the 10.31 s moment, the damage zones of the surrounding rock around three caverns remain constant.e evolution of damage zones around three caverns with time is shown in Figure 18.
Figures [15][16][17][18] indicate that the damage zone near the middle of side walls of powerhouse and tailrace surge chamber slightly extends inward as the input acceleration increases.When the acceleration approaches the peak, the damage zone   Advances in Materials Science and Engineering extends greatly.e damage zone is generally located near the tops and sides of the tailrace surge chamber and powerhouse as well as the bottom and sides of the transformer chamber.
e damage zone area reaches the maximum when the peak acceleration arrives, and then remains constant.

Kobe Earthquake.
e distributions of damage zone of the surrounding rock under Kobe earthquake at the 4.389 s moment, 4.782 s moment, and 7.115 s moment are shown in Figures 19-21, respectively.e evolution of damage zones around three caverns with time is shown in Figure 22.It is   increases.e damage zone is generally located near the tops and sides of the tailrace surge chamber and powerhouse as well as the bottom and sides of the transformer chamber.At the 4.782 s moment, a continuous damage zone occurs in the region between the transformer chamber and tailrace surge chamber, which does not happen under the action of other earthquakes.

ree Artificial
Earthquakes.e evolution of damage zones around three caverns with time under the first artificial earthquake is shown in Figure 23.e damage zone near the middle of side walls of powerhouse and the tailrace surge chamber strongly extends inward as the input acceleration increases.However, the extension area under the first artificial earthquake is less than that under the Kobe earthquake.e damage zone is generally located near the tops and sides of the tailrace surge chamber and powerhouse as well as the bottom and sides of the transformer chamber.e evolution of damage zones around three caverns with time under the second artificial earthquake is shown in Figure 24.e total damage zone areas of surrounding rock around three caverns at 1.979 s, 6.35 s, and 7.2 s moments are 1750 m 2 , 1799 m 2 , and 1805 m 2 , respectively.e damage zone extends little under the second artificial earthquake.
e evolution of damage zones around three caverns with time under the third artificial earthquake is shown in Figure 25.e damage zone gradually extends inward as the input acceleration increases.
e damage zone area reaches the maximum immediately after the peak acceleration arrives and then keeps constant.Generally speaking, the damage zone extends a little under the third artificial earthquake.e damage zone is generally located near the tops and sides of the tailrace surge chamber and powerhouse as well as the bottom and sides of the transformer chamber.
6.4.Discussion. Figure 26 shows the maximum damage zone areas under various earthquakes and no earthquake, which indicates that the maximum damage zone areas vary greatly with Fourier spectrum in the case of the same peak acceleration.Among the five earthquakes, Kobe earthquake caused the greatest damage zone with area of 4378 m 2 , while the second artificial earthquake resulted in the least damage zone with area of 1805 m 2 .
e natural frequency of underground caverns is solved to be approximately 1.9 Hz.Advances in Materials Science and Engineering 2.3 Hz and 3.2 Hz, Kobe earthquake 1.7 Hz, first artificial earthquake 1.4 Hz and 2.9 Hz, second artificial earthquake 2.6 Hz and 4.3 Hz, and third artificial earthquake 3.4 Hz. erefore, the maximum damage zone area will increase when the predominant frequency of underground caverns approaches the natural frequency of underground caverns and vice versa.

Conclusions
In this paper, taking the Dagangshan Hydropower Station in China as a study case, the dynamic response of underground caverns was numerically simulated under various Fourier spectra without consideration of reinforced supports such as anchor.e following results are drawn: (1) Even if the acceleration attenuation effect with depth was considered, a continuous damage zone occurs in the region between the transformer chamber and tailrace surge chamber under the Kobe earthquake, which may result in severe failure of caverns.(2) e damage zone near the middle of side walls of the powerhouse and tailrace surge chamber gradually extends inward as the input acceleration increases.e damage zone area reaches the maximum when the peak acceleration arrives and then keeps constant.
(3) e maximum damage zone area varies greatly with Fourier spectrum in the case of the same peak acceleration.e maximum damage zone area will increase when the predominant frequency of underground caverns approaches the natural frequency of underground caverns and vice versa.

Figures 7 -
show that the Panzhihua earthquake has the predominant frequencies of 2.3 Hz and 3.2 Hz, Kobe earthquake 1.7 Hz, first artificial earthquake 1.4 Hz and 2.9 Hz, second artificial earthquake 2.6 Hz and 4.3 Hz, and third artificial earthquake 3.4 Hz.e five

Figure 17 :
Figure 17: Damage zone from 10.31 s to end.

Figure 18 :
Figure 18: Damage zone area with time.

Figure 21 :
Figure 21: Damage zone from 7.115 s to end.
show that the Panzhihua earthquake has the predominant frequency of

Figure 24 :
Figure 24: Damage zone area with time under the second artificial earthquake.

Figure 25 :)Figure 26 :
Figure 25: Damage zone area with time under the second artificial earthquake. Figure 2: Fourier spectrum of Panzhihua earthquake.

Table 1 :
Mechanical parameters for numerical simulation.