Microsoft Word - numero_30_art_63 L. Bing et alii, Frattura ed Integrità Strutturale, 30 (2014) 526-536; DOI: 10.3221/IGF-ESIS.30.63 526 Time-history simulation of civil architecture earthquake disaster relief- based on the three-dimensional dynamic finite element method Liu Bing College of Mechanical and Electronic Engineering, China University of Petroleum, Qingdao 266580, China College of Mechanical and Electronic Engineering, Shandong University of Science and Technology, Qingdao 266510, China. liubing5195@163.com Qi Yaoguang, Du Jiyun College of Mechanical and Electronic Engineering, China University of Petroleum, Qingdao 266580, China ABSTRACT. Earthquake action is the main external factor which influences long-term safe operation of civil construction, especially of the high-rise building. Applying time-history method to simulate earthquake response process of civil construction foundation surrounding rock is an effective method for the anti-knock study of civil buildings. Therefore, this paper develops a civil building earthquake disaster three-dimensional dynamic finite element numerical simulation system. The system adopts the explicit central difference method. Strengthening characteristics of materials under high strain rate and damage characteristics of surrounding rock under the action of cyclic loading are considered. Then, dynamic constitutive model of rock mass suitable for civil building aseismic analysis is put forward. At the same time, through the earthquake disaster of time-history simulation of Shenzhen Children’s Palace, reliability and practicability of system program is verified in the analysis of practical engineering problems. KEYWORDS. Earthquake disaster; Three-dimensional dynamic finite element method (FEM); Civil architecture; Time-history analysis method. INTRODUCTION nderground structure seismic research involves many subjects such as seismology, dynamics, rock mechanics, the theory of spray anchor bracmg and so on and it has always been a hot research topic at home and abroad. This paper from the viewpoint of engineering application, focuses on civil building earthquake disaster dynamic process simulation and system application development, carried on the comprehensive study of the civil construction dynamic time-history seismic analysis method and developed the three-dimensional dynamic finite element numerical simulation system of civil building earthquake disaster process. The explicit central difference method is adopted to realize the incremental solution to basic differential equation of kinematics. Strengthening characteristics of materials under high strain rate and damage characteristics of surrounding rock under the action of cyclic loading are considered and dynamic constitutive model of rock mass suitable for civil building aseismic analysis is put forward. Finally, a three-dimensional dynamic simulation of the earthquake disaster process was carried out of Shenzhen Children’s Palace. The results show that calculation speed can satisfy the requirement of engineering analysis program. The law is correct and calculation result is reasonable. Therefore, it can be applied to earthquake disaster process simulation in actual underground engineering. U L. Bing et alii, Frattura ed Integrità Strutturale, 30 (2014) 526-536; DOI: 10.3221/IGF-ESIS.30.63 527 THE EXPLICIT FINITE ELEMENT SOLUTION OF THE SEISMIC WAVE FIELD henzhen Children’s Palace, which is under design is located at the foot of Lotus Hill. Its seismic response is mainly affected by the seismic wave propagation, which is a question of near field wave. Solution of seismic wave field in the area of engineering is the key to the earthquake disaster process simulation of the Children’s Palace. For three- dimensional simulation for complex engineering, finite element method is generally used to realize the integral solution of wave equation. The solution ways of finite elements can be divided into two kinds of implicit and explicit. The Explicit Central Difference Method Based on system kinematics basic theorem, wave equation can be expressed as differential equation of motion of Incremental Lagrange format [11] as follows:   int ext dampMa f f f (1) In the equation, a is the acceleration of mesh nodes; int , extf f are internal forces and external forces of mesh nodes respectively; dampf is the damping force; M is the mass matrix. The lumped mass matrix form is adopted in the program and is defined as follows:       T M ρ dV (2) In the equation,   is the matrix ofi , i selects 1 within the area of i and 0 outside the area. In the program, for hexahedral 8 node unit, the volume of node region of i is 1 8 of its neighboring units. For ease of matrix inversion, local damping force is adopted, defined dampf as follows:   int damp extf α f f signa (3) In the equation, signa is sign of grid node's speed and direction.  is the local damping coefficient,   D . D is the critical damping ratio. For geotechnical materials, the scope of D is generally 2%~5%. For the explicit solution of Eq. (1) there are many methods. And central difference method with second order accuracy is adopted in this paper, as shown in Fig. 1: 0t 1nt 2/1nt nt nt 2/1nt 1nt 2/1 nt2/1 nt Figure 1: Timeline of central difference method Assume that displacement, velocity and acceleration of 1 20,  nt , t , ,t are known, then, displacement, velocity and acceleration of 1nt can be calculated by central difference approximation, as follows: S L. Bing et alii, Frattura ed Integrità Strutturale, 30 (2014) 526-536; DOI: 10.3221/IGF-ESIS.30.63 528     1 1 intn ext dampa M f f f (4)   1 2 1 2  n / n- / na a Δta (5)     1 2 1 21 n nn na a Δt a (6)           1 2 1 2 1 2 12n / n- / n / n nt t t / , t t t (7) In the equation, 1n na , a are respectively displacement of 1n nt ,t ;  1 2 1  n n / na , a , a are respectively velocity of  1 2 1n n / nt ,t ,t ; na , 1na is respectively accelerated speed of 1n nt , t . The explicit central difference method calculation process can be described as follows: (1) According to the initial condition 0a , 0a , we can calculate that    0 1 0 1 2 0,  damp n- /a M f f a a . (2) For each time step: 1 2n /a is calculated by Eq. (5); displacement 1na of 1nt is calculated by Eq. (6); accelerated speed 1na is calculated by Eq. (4). Finite element explicit calculation is a conditional stability algorithm which is often constrained by the solution stability and calculation time consuming. This is because to ensure the stability of calculation, we generally need to select the smaller calculation time step t to increase system cycle count and extend the calculation time. Therefore, on the premise of ensuring stability calculation, we need to choose the larger t . Considering the physical meaning of the wave finite element, t can be regarded as the minimum duration in a computational time step, that wave propagation distance is not more than any one unit size. It can be calculated as follows:     min 0 80 0 98e e e l t α . α . C (8) In the equation, eC is the unit wave velocity; el is the ratio of the largest area between element volume and unit 6 surface. The system time step of Eq. (8) only considers the seismic wave transmission stability in calculation model rather than artificial boundary conditions and computational stability requirements of supporting measures such as anchor bars to the system time step. Therefore, in the practical engineering calculation, according to different situations, we need to make adjustments of t . Node Internal Force Solution In the explicit solution of type (1), the key is the solution to node internal force fint. In incremental variational type calculation, we first calculate unit stress increment and update the unit stress. Then the updated nodes stress is integral to the unit nodes to obtain internal forces. Deformation description In the Delta Lagrange explicit calculation, deformation is described by strain rate ije :                       1 1 2 2 i i j i j ij ij ij j j i i i v v v v v e Ω D x x x x x (9) In the equation, ijΩ is the spin rate tensor; ijD is the deformation rate tensor; iv is the velocity component; jx is the weight coordinates, ,i j =1, 2, 3. L. Bing et alii, Frattura ed Integrità Strutturale, 30 (2014) 526-536; DOI: 10.3221/IGF-ESIS.30.63 529 Stress update The numerical algorithm of integral rate constitutive equation is called as stress update algorithm. Stress  ijσ t dt of t dt can be obtained by ij as follows:        ij ij ijσ t dt σ t dt (10a) The Cauchy stress rate can be expressed as:  ij      ij ik jk jk ikΩ Ω (10b) In the equation,  ij is the Jaumann rate,  ij  ijkl klC D , ijklC is the elastic constitutive tensor. Stress integral  ijσ t dt obtained through Eq. (10a) is performed to the unit nodes, and then node’s internal force vector intf is obtained. THE DYNAMIC ELASTOPLASTIC DAMAGE CONSTITUTIVE OF THE SURROUNDING ROCK The Dynamic Elastoplastic Damage Constitutive Model of Civil Construction Foundation Surrounding Rock he dynamic elastoplastic damage constitutive model of civil construction foundation surrounding rock is the basis of simulation of the Children’s Palace earthquake disaster process. On the basis of theoretical study, we need to combine with a large number of indoor and outdoor tests, and establish a scientific and practical dynamic damage constitutive model of rock. Civil construction foundation surrounding rock shows the dynamic strengthening characteristics and fatigue damage properties of the surrounding rock. The research results [5] showed that: Under dynamic loading, the surrounding rock strain rate increases. And elastic modulus of surrounding rock materials and damage resistance are improved. While under cyclic loading conditions, there appears fatigue damage of surrounding rock and elastic modulus of surrounding rock materials and destruction resistant properties are reduced. The two factors are in opposition to each other and coexist at the same time. It should be considered in the constitutive model as follows: (1) The study on rock dynamic strengthening characteristics is mainly divided into high strain rate    2 110 s and middle strain rate      4 1 1 110 10s s study. Among which, high strain rate is mainly to solve the problem of shock blasting. Because the loading ways are simple, there are relatively more tests. While Children’s Palace earthquake effect strain rate is generally at the range of     6 1 1 110 10s s , belonging to medium strain rate and there are less tests. From the existing research results [5, 14 ~ 17] we can get the following conclusions:  Rock dynamic strength and modulus of elasticity increase with the increase of strain rate.  Rock dynamic strength increases with strain rate. This is because that cohesion increases with strain rate while internal friction angle is not affected by strain rate.  The dynamic characteristics of medium strain rate of the rock are affected by the way of loading, loading rate, rock type and test equipment. All kinds of test results are with large discreteness, so consistent rules are difficult to conclude based on experiments. (2) There are many rock damage constitutive researches. Damage description is usually with tensor description form of first, second and fourth order. In view of the aeoplotropism of seismic wave propagation, this paper adopts second order damage tensor to describe the damage of surrounding rock in the earthquake. Considering the above two kinds of characteristics of surrounding rock in the earthquake, the dynamic elastic modulus of surrounding rock of dynamic constitutive model can by expressed as:  0 E P E (11) In the equation, E is the static elastic modulus of surrounding rock; T L. Bing et alii, Frattura ed Integrità Strutturale, 30 (2014) 526-536; DOI: 10.3221/IGF-ESIS.30.63 530  P is the strain rate function greater than 1. The Mohr-Coulomb yield criterion is adopted and assumes that the compressive stress is negative and tensile stress is positive. Line AB is the shear yield criterion  0sf , line BC is the tensile yield criterion  0tf , that is:               * * 1 3 * 3 2 ( )s t t f N TQ c N f (12a) Among which, N       1 sin 1 sin (12b) In the equation, ,c are respectively cohesion and internal friction angle; * 1 σ , * 3 σ are respectively first and third effective principal stress;  t is ultimate tensile strength,  max / tan t c ; T is the damage effect coefficient. Under the second order tensor form,    2 2 21 2 31T D D D , 1 2 3, ,D D D are damage coefficient of principal stress direction, solved by the damage evolution equation. Within every calculation step,  , T Q is constant. According to the non-correlation of surrounding rock, corresponding potential function to yield criterion is as follows:       * $ 1 3 * 3 s ψ t g σ σ N g σ (13a) Among which,       1 sin 1 sin N (13b) In the equation,  is the surrounding rock dilatancy angle. In the calculation, element stress applies the effective force form considering damage:                   1 1 2 2 3 3 1 0 0 0 1 0 0 0 1 * σ D h σ σ D h σ D h (14) Under the tension cases, assume  1.0h ; in press conditions,  0.2h . The damage process of surrounding rock is along with the micro fissure development and micro porosity increase and decrease of strength of materials. In the larger sense, damage of surrounding rock is caused by tension. According to some practical observation results, most of the surrounding rock tensile damage is caused by the actual strain exceeding the limit tensile strain of surrounding rock. For the three-dimensional situation, it can be considered that element stress enters plasticity and first principal tensile strain exceeds ultimate tensile strain. Then damage will occur, that is:                     1 1 0 0 0 s s ij or f D h and f (15) In the equation,   ijh is the three-dimensional damage evolution equation, L. Bing et alii, Frattura ed Integrità Strutturale, 30 (2014) 526-536; DOI: 10.3221/IGF-ESIS.30.63 531 1 is the first principal tensile strain of the surrounding rock,   is the ultimate tensile strain of the surrounding rock. For the situation of the ultimate tensile strain of untested surrounding rock in practical engineering, we adopt specific value of rock tensile strength sR and elasticity modulus E , that is:      sR kE (16) In the equation, k is the safety coefficient. In the damaged area, surrounding rock stress and physical parameters decrease with the increase of slant plastic strain. The greater the cumulative plastic deformation is, the higher the damage degree of surrounding rock will be. In the adjacent area in front of the ultimate strength of surrounding rock material, micro crack damage increase is very rapid. This kind of change rate of damage can usually be described by exponential function. Therefore, three-dimensional damage evolution equation can be expressed as follows:        1 exp 1, 2, 3p pi i iD R i (17) In the equation,  pi is the slant plastic strain in the principal strain direction, R is damage constant of the material. For the purposes of showing the state of surrounding rock damage in the post-processing, second order damage tensor is converted to a first-order damage scalar when output the results file, expressed as   2 2 21 2 3D D D D . Stress Correction and Damage Calculation Unit stress calculation is the core of the explicit finite element method. In the elastic and plastic damage calculation of the tn+1 moment, the gauss point elastic stress is calculated using unit material initial parameters and whether unit is in a state of damage is judged (when    0D tn , unit damage coefficient of the moment is calculated). If the unit is in a state of damage, each gaussian point stress is revised as the effective stress. On this basis, considering damage effect on yield surface and carry out yield judgment. In case of yield, plastic stress correction is required. Also, we need to judge whether principal tensile strain exceeds the ultimate tensile strain. If there is exceeding, it's considered unit damage and damage coefficient is calculated. This stress correction method is applicable to the situation that yield function  nf is the linear function. Based on the principle that plastic stress should flow on the yield surface, it is required that plastic stress of tn + 1 moment as follows:     0n nf σ σ (18) For the linear function, Eq. (18) can be converted to:       0n n*f σ f σ (19a)          0n*f f f (19b) In the equation, * (·)f is calculated by yield function minuses to the constant term  0nf . Based on this transformation, linear yield function plastic stress correction formula was deduced by service manual [20] as follows:             N I i i i n g S (20) In the equation,  is the plastic flow factor; iS is the elastic constitutive; L. Bing et alii, Frattura ed Integrità Strutturale, 30 (2014) 526-536; DOI: 10.3221/IGF-ESIS.30.63 532  Ii ,  N i are respectively elastic assumed stress and the revised plastic stress. Within each computing time step t, in the constitutive model,  Q is the constant and yield function is the linear function. Eq. (20) can be applied to carry out plastic correction of unit effective stress considering the damage. For the plastic shear failure, there is:                     * 1 1 1 2 * 2 2 * 3 3 2 1 2 1 N I N I N I s ψ s ψ s ψ σ σ λ α α N σ σ λ α N σ σ λ α N α (21a)                * * 1 3 1 2 1 2 ,s I I s ψ f σ σ λ N N N (21b) Among which,     1 2 4 2 , 3 3 K G K G (21c) In the equation, K , G are respectively bulk modulus and shear modulus;  *In is the equivalent elastic assumed stress after damage correction. For the ultimate tensile damage, there is:                  * * 1 1 3 * * 2 2 3 3 2 1 2 1 Nt I I Nt I I Nt t t t α σ σ σ σ α α σ σ σ σ α σ σ (22) In the damage calculation, damage coefficient iD should be calculated according to the damage evolution equation. Computation formula is as follows:                      ( 1) ( ) * 2 ( 1) ( 1) 01 exp ( ) p p s i tn i tn I i p p i tn i tn F D R (23) Among which,          0 1 ( 1) 2( 1) 3( 1) 3p p p ptn tn tn (24) ENGINEERING PROJECTS Establishment of the Model he Children’s Palace subject has a maximum altitude of 37.25 m and a maximum span of 24 m. Its foundation is in the dolomite limestone of Lower Cambrian system. The rock is of hard property and fissure is of medium development. There are mainly III category and partial IV category of surrounding rock. Physical and mechanical parameters are shown in Tab. 1. This example temporarily is regardless of the anchor rod, anchor rope, lining and other supporting measures. Since this project did not carry out live stress test, initial geostatic stress field calculation in this paper considers self-weight stress field affected by river valley and excavation slope. And local damping is adopted and set as 0.157 1. This project did not carry out the test of rock mass deformation modulus and cohesion-strain rate test. Considering engineering safety, the P, Q value is set as 1. T L. Bing et alii, Frattura ed Integrità Strutturale, 30 (2014) 526-536; DOI: 10.3221/IGF-ESIS.30.63 533    3/ g cm E GPa  c MPa         t MPa 2.7 4.0 0.3 0.8 45 43 0.3 Table 1: The Children's Palace foundation rock physics parameters Considering that the Children’s Palace is close to the earth surface, seismic wave field inversion model in the engineering area and the Children’s Palace cavity analysis model can use the same model. Finite element model is established including administrative building, draft tube gate chamber, the tailrace tunnel, bus tunnel and diversion tunnel. And it is a total of 272 304 hexahedron eight nodes unit containing 301112 nodes. The axis of the main building is the y direction X axis and y axis are perpendicular to the downstream, z-axis coincides with geodetic coordinates. The calculation scope along three direction x, y, z axis are respectively 176.0, 108.0, 141.0 m. The Children’s Palace finite element calculation model is shown in Fig. 2. To reflect the three dimensional input features, seismic wave in this paper applied three-dimensional earthquake acceleration time history recorded by Wolong TV station in WenChuan earthquake and make an interception of 20 s for calculation. This project filter frequency range is 0.1~10.0 Hz. Corresponding maximum grid feature sizes is not greater than 7.0 m. After filtering, baseline correction were carried out on original recorded seismic waves, measured acceleration time history curve is as shown in Fig. 2. Seismic wave incident reference point is the origin of coordinates, as shown in Fig. 2 (a), incident direction vector is (1, -1, 1). Note: look for a new perspective of (a) and change (b) into an architectural renderings of The Children’s Palace (shown in Fig. 3), keeping the draft tube gate chamber, bus tunnel, diversion tunnel, tailrace tunnel. And make administration main workshop into the hall. (a) (b) Figure 2: The finite element calculating model for the children's palace: (a) FEM model; (b) Excavation model Cavity Earthquake Disaster Simulation On the basis of rock damage and stress disturbance around Children’s palace foundation, datum plane wave and power attenuation coefficient β would be of inversion. The process of earthquake disaster in Children's Palace group will be calculated and simulated. There are about 27 myriad of hexahedral eight nodes units in this model�dynamic calculation lasts about 25 s (there is seismic wave input in former 20 s). FLAC3D general numerical analysis software is adopted about 10.8h, then Fig. 4 is obtained. The analysis of calculation results is as following: 1) As can be seen from the Fig. 4, the displacement time history law is basically in line with inversion input seismic waves. Because of the airport surface effect, peak ground displacement of each monitoring point has experienced a sharp increase compared with input seismic waves. Among which, the main building roof arch monitoring point A increased by 80%, L. Bing et alii, Frattura ed Integrità Strutturale, 30 (2014) 526-536; DOI: 10.3221/IGF-ESIS.30.63 534 upstream wall monitoring point B increased by 55%, downstream sidewall monitoring point D increased by 45%, slope monitoring point F increased by 58%. The amplification effect of free face is unfavorable for surface surrounding rock stability. When the 25s of calculation is finished, along with layman of seismic wave and damping dissipation, kinetic energy of the whole of model system is basically to 0 and reaches a steady state. Figure 3: 3D model for the Children's Palace. 0.20 0.00 0.04 0.08 0.12 0.16 0 5 10 15 20 25 t/s Crown of measuring point A Upstream sidewall point B Downstream sidewall point D The slope point F Figure 4: Displacement time history curve of dynamic calculation process monitoring. Under the impact of wave time difference, space contour, plastic damage and damping dissipation, displacement time history curve of the monitoring point shows difference in quantity and there appears the relative motion. This is also an important reason for instability of the children's palace under the action of earthquake. From Fig. 4, we can see that largest relative displacement between main workshop upstream sidewall and downstream sidewall is 6 cm; which is 50% of input waveform amplitude load value. Considering install the backfill of concrete to the factory units can slow down the situation to a certain extent. 2) As can be seen from the Fig. 5, under the action of seismic wave load, third principal stress of surrounding rock is in the cycle of severe shock. This kind of cyclic loading results in a sharp increase in damage of surrounding rock and can cause fatigue failure of surrounding rock. The three curves show that after the calculation unit stress state and the initial stress state are quite similar. It suggests that after the earthquake surrounding rock stress state did not produce large disturbance. Due to unit stress state in the earthquake changing with time, stress state cannot reflect the stability of surrounding rock objectively. The author thinks that damage degree and stress state should be adopted together to describe degradation degree and the stability of surrounding rock in the earthquake. To sum up, in the earthquake disaster of Shenzhen Children’s palace, cave surrounding rock produced large absolute and relative displacement (upstream and downstream sidewall 6cm). In the process of seismic wave propagation, surrounding rock stress disturbance is strong. Some areas were in the alternative changing of elastic and plastic stress state and caused fatigue damage of surrounding rock. With the development of time, cumulative plastic strain increases, which causes irreversible damage to the surrounding rock material and reduce of the bearing capacity of surrounding rock and surrounding rock stability is influenced. At the 13″ moment, part of rock whole interval was filled with plastic damage area. After the calculation of 25″ moment, cave surrounding rock plastic damage zone was widened and damage degree of surrounding rock was significantly improved. Under earthquake loading effect, Children’s palace foundation rock mass may lose stability in the earthquake. L. Bing et alii, Frattura ed Integrità Strutturale, 30 (2014) 526-536; DOI: 10.3221/IGF-ESIS.30.63 535 Figure 5: Time-history curves of 3rd principal stress at typical area of surrounding rock. Figure 6: Time-history curves of 3rd plastic damage elements of global model. We could see from Fig. 6 that, with the spread of seismic wave, part of the surrounding rock has been in the alternate change state between elastic stress and plastic stress for a long time. At the moment of 13 s, the number of whole model’s plasticity damage has reached peak, there are 56 545 units as a whole. The biggest moment of plastic damage zone is consistent with seismic displacement peak time. At the time of 13 s, most of the sidewall rocks in Children’s Palace have entered the plastic damage area (see Fig. 6). The depth of local cracking area is as high as 7m, plastic damage zone depth is as much as 20m, Children’s Palace has the possibility of instability. CONCLUSION his paper is to expound the finite element simulation theory of earthquake disaster in Children’s Palace and develop three-dimensional dynamic finite element numerical simulation system in the earthquake disaster process of Children’s Palace through starting from engineering analysis. The dynamic constitutive model of rock mass adaptive to seismic analysis of Children’s Palace is put forward through adopting central difference method and considering the enhancement of material parameter, as well as fatigue damage of material under circulating load effect. The earthquake disaster simulation is performed on Children’s Palace in Shenzhen through applying software. It has T L. Bing et alii, Frattura ed Integrità Strutturale, 30 (2014) 526-536; DOI: 10.3221/IGF-ESIS.30.63 536 indicated that computation speed of procedure in this paper could satisfy the requirement of engineering analysis, and this procedure has performed relatively real and reasonable simulation on earthquake disaster of Children’s Palace. Due to that article’s length is limited, there inevitably exists vulnerable spot of initial defects and structure on other perspective. The architecture is regarded as continuum medium field through applying damage mechanics. The material containing various microfracture and microdefect is indistinctly deemed as continuous medium with damage field. The occurring and development of damage is considered as the evolution process of damage. The appropriate damage variable is introduced to describe physical and mechanical properties of this continuum. Heterogeneity and mesoscopic structure characteristics are grasped. The corresponding numerical model and theoretical analysis is constructed at micro and macro level respectively. Analysis about the initiation extension of building micro cracks, formation of macroscopic fracture, adhesive failure, damage fracture criterion as well as the influence of mesoscopic composition on macro elastic modulus and so on is performed. The failure mechanism of structure is further demonstrated, and it has certain promotive role in perfecting material and structural design. ACKNOWLEDGEMENTS he issue of national sci-tech support plan -‘The application demonstration of industrializing oil equipment's research and development platform’, No (2012BAH26F04). REFERENCES [1] Liu, J. B., Du, Y. X., Wang, Z. Y., 3D viscous-spring artificial boundary in time domain, Earthquake Engineering and Engineering Vibration, 5(1) (2006) 93–102. [2] Du, X. l., Zhao, M., Wang, J. T., Artificial stress boundary condition of the near field wave simulation, Acta Mechanica Sinica, 38(1) (2006) 49–56. [3] Li, X. J, Lu, T., Explicit finite element analysis of earthquake response for underground caverns of hydropower stations, Journal of Hydroelectric Engineering, 28(5) (2009) 41–46. [4] Zhang, X. Z., Xie, L. L., Problems in numerical solution of complex open system by using explicit finite element method, Earthquake Engineering and Engineering Vibration, 25(2) (2005) 10–15. [5] Qian, Q. H., Qi, C. Z., Dynamic strength and dynamic fracture criteria of rock and rock mass, Journal of Tongji University (Natural Science), 36(12) (2008) 1599–1605. [6] Qi, C. Z., Qian, Q. H., Physical mechanism of dependence of material strength on strain rate for rock-like material, Chinese Journal of Rock Mechanics and Engineering, 22(2) (2003) 177–181. [7] Zhao, J., Li, H. B., Estimating the dynamic strength of brittle rock using Mohr-Coulomb and Hoek-Brown criteria, Journal of Rock Mechanics and Engineering, 22(2) (2003) 171–176. [8] Li, H. B., Zhu, L., Lu, T., Dynamic response analysis of large underground excavations in jointed rock, Chinese Journal of Rock Mechanics and Engineering, 27(9) (2008) 1757–1766. [9] Chen, J. Y., Hu, Z. Q., Lin, G., 3D seismic response study on large scale underground group caverns, Chinese Journal of Geotechnical Engineering, (2001). T