Acta Polytechnica CTU Proceedings doi:10.14311/APP.2019.23.0038 Acta Polytechnica CTU Proceedings 23:38–43, 2019 © Czech Technical University in Prague, 2019 available online at https://ojs.cvut.cz/ojs/index.php/app EVALUATION OF UNDERGROUND STRUCTURES SUBJECTED TO SEISMIC LOADS Jan Pruška Czech Technical University in Prague, Faculty of Civil Engineering, Department of Geotechnics, Thákurova 7, 166 29 Prague, Czech Republic correspondence: Pruska@fsv.cvut.cz Abstract. The paper is focused on the evaluation of the effect of earthquakes on underground structures. Free-field analysis is one solution of this task common mainly in engineering tunnelling practice, but it has some rather simplified aspects (e.g. equivalent shear strain is constant). Pseudo- static finite element calculation combines free-field analysis and the advantages of a FEM model. Dynamic effects are introduced in the form of displacements prescribed along the vertical boundaries of the FEM model in a usually static manner. This approach also implies constant material parameters for the geological profile in the horizontal direction, an arbitrary geometry of excavation, soil structure interaction and description of share waves as a time-dependent 1D analysis of the so called free-field column. Moreover, there is shown an example comparing pseudo-static FEM analysis with an analytical method. Finally, the advantages of the pseudo-static FEM method are presented. Keywords: Finite element method, earthquake, 1D free field analysis, dynamics, pseudo static analyses. 1. Introduction The response of underground structures subject to seismic load can be solved in several ways, which can be roughly classified into four categories: laboratory testing, analytical solutions, pseudo static finite ele- ment calculations and full dynamic analysis. (1.) Centrifuge modelling techniques allow well for control testing conditions regarding investigations of soil structure interaction in simulated seismic events. However, there are several difficulties for a dynamic centrifuge model such as accurate modelling of the boundary conditions or the scaling law for dynamic velocity. (2.) Most simplified methods were derived from closed form solutions but involving neglect of some fac- tors. In general, the methodology of the analytical method is that the seismic loads for an underground structure are characterized in terms of strains im- posed on the structure by the surrounding ground or their interaction. Although this is the most con- servative method having only a limited application, it does provide a first-order estimate of the liming stresses so to potentially serve as an accuracy check on more complex calculations. (3.) A pseudo-static finite element calculation pro- vides a solution for any shape of underground structure in a generally non-homogeneous layered rock/soil mass with a potentially nonlinear material response. The basic idea of this approach is that the maximum shear deformation in a free-field con- dition of the layer represents the maximum dynamic earthquake stress in this layer and can be an input to the numerical model as a boundary condition. The shear strain equivalent to the actual dynamic loading can be calculated from the particle velocity estimated either in a simplified manner based on limited information or from the solution of a one dimensional simulation of free-field response to an actual earthquake. (4.) Full dynamic analyses are based on a direct ap- plication of recorded accelerations to the actual computational model of the tunnel where the soil and tunnel responses are mechanically coupled and analysed via numerical modelling, such as finite element methods. 2. Analytical method The methodology of the analytical method for a seis- mic loading design is based on the principle that the static condition has to include the additional loading imposed by seismicity. In the case of underground structures, the seismic loads are described in terms of strains imposed on the structure by the surrounding ground mass or their interaction. 2.1. Seismic coefficient method This method was originally developed for dynamic analysis of over ground buildings. Inertial forces due to earthquakes F1 of the structure and F2 of the upper soil mass are described by the following equations: F1 = a g · Q = Kc · Q, (1) F2 = ηc · Kh · m1 · g. (2) This method is not suitable for tunnels with large diameters because this method does not reflect the real dynamic property of tunnel structure. 38 https://doi.org/10.14311/APP.2019.23.0038 https://ojs.cvut.cz/ojs/index.php/app vol. 23/2019 Underground structures subjected to seismic loads 2.2. Free – field deformation approach Tunnel structures are assumed to move together with the surrounding soil/rock mass during earthquakes and therefore undergo the following deformations: • ovaling, • axial, • curvature. The ovaling deformation of a circular tunnel lining is primarily caused by seismic waves that propagate perpendicular to the longitudinal axis of the tunnel (Figure 1). The vertically propagating shear waves produce usually the most critical ovaling distortion of the lining. The free field deformation approach assumes that the deformation of the underground structure should conform to the deformation of the soil in the free field under seismic shaking. In these simplified methods described in [1–4] the soil – struc- ture interaction and the effects of the waves are ne- glected. Moreover the maximum shear strain γmax (Equation 6) induced by shear waves propogating in vertical planes is constant along the depth z. This may lead to an overestimate of tunnel lining deformations especially in loose rock and soils. Figure 1. Ovaling deformation of a Circular Cross Section. In the case of a circular tunnel in the homogeneous medium the ovaling deformation is determined by the following steps: (1.) according to the Seismic Report we get the Pick Ground Acceleration (3) and the site-specific Pick- Ground-Acceleration (4): agR = 40% · Ss, (3) amax,s = S · agR (4) (2.) and then we determine the pick acceleration at the depth of the tunnel: az,max = C · S · agR, (5) (3.) az,max max is used to determine maximum shear deformation in free-field condition as shown in equa- tions (6) and (7): γmax = Vz Cs , (6) Vz = k · az,max (7) (4.) and finally we obtained the ovaling deformation of the circular tunnel: ∆dff = γmean d 2 . (8) 2.3. Soil-structure interaction approach The kinematic interaction between the tunnel and the surrounding soil/rock mass is used to improve the accuracy of the widely used free-field methods. Gener- ally, analytical solutions for estimating soil-structure interaction for tunnels are based on the following as- sumptions: • the tunnel is circular; • the soil is an infinite, elastic, homogeneous and isotropic medium; • the tunnel lining is generally an elastic tube under plane strain conditions or it is based on the theory of an elastic beam on an elastic foundation; • the interface between the soil and the tunnel lining has a full-slip or a no-slip condition; • the soil-structure interaction effects operate in a quasi-static manner, ignoring any inertial interac- tion effect. There exist many approaches in the technical liter- ature describing seismic design and analysis of under- ground structures using a soil structure interaction approach. For example: Billota et al. [5] discussed four different closed-form elastic solutions to evaluate the maximum shear strain γmax (all approaches are based on the equilibrium of the soil column from the surface to the given depth z), while expressions for displacements and internal lining forces in the circular tunnels excavated in a rock mass with a near fault are presented in [6], and closed-form expressions to calculate the circular tunnel liner forces due to the compressive seismic P-wave propagation are described by Kouretzis et al. [7]. In addition an analytical solution for estimating internal forces in the lining of a circular tunnel in the homogeneous medium is published in [8, 9]. The ovaling deformation is determined by the assumed shear strain ∆d (Figure 1). We get for parameters ∆dff, ∆dlin the follow-on formulas: ∆dlin = Rn∆dff, ∆dff = γmean d 2 , Rn = 4(1 − νm) αn + 1 ,αn = 12El(5 − 6) d3Gm(1 − ν2l ) . (9) It is evident that the material and geometric parame- ters Rn and αn take into account the rigidity of the subsoil-lining system. Expressing the curvature of an elliptical excavation using ∆dlin and parameters describing soil structure interaction then the moment in the lining (Figure 2) can be written using the polar 39 Jan Pruška Acta Polytechnica CTU Proceedings coordinate ϕ (positive clockwise from horizontal) as: M(ϕ) = 6El∆dlin d2(1 − ν2l ) cos (2(ϕ + φ 4 )). (10) With reference to Figure 2 (c), the normal force on the curved element can be expressed as: N(ϕ) = − 1 R dM2 dϕ + PnR (11) and after rearranging using equation (10): N(ϕ) = 48ElI∆dlin d3(1 − ν2l ) cos (2(ϕ + φ 4 )) + Rpn (12) Figure 2. Internal forces and contact stress on the differential element. 3. pseudo static FEM Implementation of free-field analysis results into the static FEM analyses can be made in several ways. Generally, the shear strains in the 2D solution without lining have to be in conformity with 1D free-field analyses and the shear deformation in the free-field deformation is applied by prescribed displacements at the boundaries of the numerical model. In the first method, we consider a homogeneous soil/rock mass with an effective shear modulus G. We determine the maximum shear deformation in free-field condition γmax (6). The value of γmax cor- responds to the maximum horizontal displacement applied in the boundaries of the numerical model (Figure 3), calculated by equation: ∆xmax = γmax · hm. (13) The second method considers the layered soil/rock mass (Figure 4) where some sort of homogenization is needed. Adopting the Voigt assumption of constant deformation at each point of the solid we obtained effective shear modulus GhomV : GhomV = h1 h G1 + h2 h G2 (14) and average shear deformation of the numerical model γmean: γmean = h1 h γ1 + h2 h γ2. (15) Figure 3. Homogenized mass and constant shear strain. Figure 4. Layered mass and constant shear strain . For the layered mass the better method is based on the setup constant corresponding movements by parts – each layer has a different shear deformation – Figure 5. The prescribed shear deformation values should be in the ratio of the stiffness of the materials in layers to ensure the continuity of the stress at the layer interface: γ1 γ2 ≈ G2 G1 (16) Figure 5. Layered mass with different values γ at layers 1 and 2. 4. Full dynamic calculations using FEM The numerical solution of dynamic load effects using the finite element method is based on the principle of virtual displacements, which express the external volume load by the effects of inertial forces and leads to equations of differential motion in the form: Mü(t) + Cu̇(t) + Ku(t) = F (t). (17) 40 vol. 23/2019 Underground structures subjected to seismic loads The assembly of the stiffness matrix C is usually a very complicated task thus the Rayleigh formulation is used i.e. as a linear combination of the mass tensor M and the stiffness tensor K: C = αM + βK. (18) Coefficients α and β have been calculated according to the double frequency method that avoids an overesti- mate of damping throughout the considered frequency range [10]. Assuming that the first natural frequency is the least damping, the equation (17) can be decom- posed to the natural vibrations. Using the damping ratio ξ we get a relationship for individual frequencies (for real systems with a large number of degrees of freedom): ξi = 1 2 ( α ωi + βωi). (19) In order for there to be a damping ratio ξ minimal for ω =ωi it must fulfill the next equation: dξ dω = 1 2 (− α ωi + β). (20) After substitution ω =ωi into equation (20) we get α = ω21β (21) and then from equation (20) β = ξ1 ω1 . (22) In common practice, equation (17) is solved by direct integration using the Newmark integration method. The discrete form of equation (17) already includes the influence of boundary conditions and the effect of loading in the form of a prescribed accelera- tion of displacements at the lower limit. 5. Case exemplification The FEM calculation results will be compared with the results of the analytical calculation published in [11]. The geometry of the model and the finite element mesh is shown in Figure 6 and Table 1 contains the parameters of the subsoil and the tunnel lining. The seismic loading has PGA 2 ms−2 and the wave period is 0.01 s. All FEM calculations were made by Geo5 Tunnel software [12]. For the sake of clarity the rep- resentative extreme values obtained for ϕ = φ4 have been set – see Table 2. For illustration the resulting moment in the lining is also shown in Figure 7. The two approaches in terms of the absolute per- centage of the difference in the last column of Table 2 can be mutually assessed. It is evident that in the case of a homogeneous mass and a circular excavation both methods give comparable results. The exception is the value of normal force in the case without a considering of the interaction between the lining and the soil (effect of the contact stress). If we take into account the value of contact stress in the equation (18), the resulting values ∆dff a ∆dlin of normal force for both methods are almost identical. Type Unit Value Unit Loading γmean 0.0105 [-] Soil Gm = Ghom 8.65 [MPa] νm 0.3 [Pa] ρm 2000 [kg/m3] Lining γmean 0.0105 [-] νl 0.3 [-] Gl 12.5 [GPa] d 6 [m] h 0.3 [m] I 0.00225 [m4] Table 1. Input parameters of the numerical model Figure 6. Numerical model. 6. Conclusions The paper is focused on the evaluation of the seismic effects on the tunnel lining by free-field analysis and pseudo static calculation using FEM. Generally, both calculation methods give comparable results in the case of a homogeneous soil environment and ignoring soil structure interaction. However, from a practical point of view, the FEM combined with the 1D free- field analysis yields clear advantages: • It allows a possibility to derive the load for any time of the earthquake (measured or synthetically generated accelerograms). • It analyses any shape of excavation (lining). • It takes into account the potentially non-linear sub- soil response. • It take into consideration the layered rock/soil mass depending on the "real" material parameters. Another indisputable advantage is taking into ac- count the soil structure interaction whose neglect in the case of an analytical method may lead to a sig- nificant overestimation of the normal force. Such a non-realistic value of the normal force leads to a lower degree of reinforcement in the effects of the bending moment for the considered cross section. 41 Jan Pruška Acta Polytechnica CTU Proceedings Variable FEM Analytical Difference [%] ∆ dff [mm] 31.3 31.5 0.6 ∆dlin [mm] 31.0 36.1 14.1 M [kNm/m] –367.6 –422.8 13.1 V [kN/m] 244.5 281.0 13.0 N [kN/m] –380.8 –562.0 32.2 N∗ [kN/m] –380.8 –390.1∗ 2.4∗ * The influence of cont. nor. stresses was taken into the calc. Table 2. Output parameters. Figure 7. Moment in the lining. List of symbols a earthquake acceleration an ratio between the stiffness of lining and massive az,max pick acceleration at the depth of the tunnel d diameter of the excavation g gravitational acceleration h thickness of the lining k ratio of PGV to PGA m1 weight of soil mass u actual displacement vector C ground motion at depth z C damping matrices Cs apparent propagation velocity of S-wave D tunnel radius El Young modulus of the lining F vector of nodal loading G shear modulus of soil/rock Gm elasticity modulus of the massive I moment of inertia of the lining K stiffness matrices Kv seismic coefficient in vertical direction Kh seismic coefficient in horizontal direction M mass matrices Q weight of structure PGA peak ground acceleration R reduction factor of depth influence R radius of excavation Rn lining – massive racking ratio under normal loading S coefficient based on the types of elastic response spectra suggested in Eurocode 8 Ss short period spectral acceleration Vs peak ground velocity α parameter that determines the influence of mass in the damping of systems αM damping proportional to displacements β parameter that determines the influence of the stiffness in the damping of the system βK damping proportional strain rate γmean average shear strain of the model ηc integrate response coefficient related to the tunnel depth, soil property and project significance νm Poisson number of the massive νl Poisson number of the lining σv Poisson number of the lining ω angular frequency ω1 angular frequency of the first natural mode ωi angular frequency of i-th mode Acknowledgements The article was produced with financial support from the Competence Centres Program of the Technology Agency of the Czech Republic (TA CR), project no. TE01020168 Centre for Effective and Sustainable Transport Infrastruc- ture (CESTI). References [1] J.-N. Wang. Seismic design of tunnel – A simple state of the art design approach. Parsons, Brinckerhoff, Quade & Douglas, Inc. 1993. Monograph 7. [2] M. S. Power, D. Rosidi, J. Kaneshiro. Strawman: screening, evaluation, and retrofit design of tunnels. National Center for Earthquake Engineering Research, Buffalo, New York. Report Draft. Vol. 3. [3] J. Penzien, C. Wu. Stresses in linings of bored tunnels. International Journal of Earth-quake Engineering and Structural Dynamics 3(27):283–300, 1998. [4] Y. M. A. Hashash, J. J. Hook, B. Schmidt, J. I.-C. Yao. Seismic design and analysis of under-ground structures. Tunnelling and Underground Space Technology (16):247–293, 2001. [5] E. Bilotta, G. Lanzano, G. Russo, et al. Pseudo-static and dynamic analyses of tunnels in transversal and longitudinal directions. In Proceedings of 4th International Conference on Earthquake Geotechnical Engineering. ICEGE, Thessaloniki, 2007. 42 vol. 23/2019 Underground structures subjected to seismic loads [6] M. Corigliano, L. Scandella, C. Lai, R. Paolucci. Seismic analysis of deep tunnels in near fault conditions: a case study in southern italy. Bulletin of Earthquake Engineering 9(4):975–995, 2011. [7] G. Kouretzis, K. Andrianopoulos, S. Sloan, J. Carter. Analysis of circular tunnels due to seismic p-wave propagation, with emphasis on unreinforced concrete liners. Computers and Geotechnics 55:187–194, 2014. [8] AASHTO. Technical manual for design and construction of road tunnels – civil elements. First Edition. American Association of State Highway Transportation Officials, 2010. [9] C. Hunk, J. Monses, N. Munfah, J. Wisniewski. Technical manual for design and construction of road tunnels – civil elements. First Edition. Federal Highway Administration, 2009. Report FHWA-NHI-10-034. [10] G. Lanzo, A. Pagliaroli, B. D’Elia. Influenza della modellazione di rayleigh dello smorzamento viscoso nelle analisi di risposta sismica locale. In Proceedings of 11th National Conference L’ingegneria Sismica in Italia. 2004. [11] D. Kučera. Posouzení geotechnické konstrukce na seismické zatížení. Diplomová práce. ČVUT v Praze, Fakulta stavební, 2017. In Czech. [12] GEO 5 FEM. Online, available from: https: //www.fine.cz/geotechnicky-software/mkp-tunel/. 43 https://www.fine.cz/geotechnicky-software/mkp-tunel/ https://www.fine.cz/geotechnicky-software/mkp-tunel/ Acta Polytechnica CTU Proceedings 23:33–38, 2019 1 Introduction 2 Analytical method 2.1 Seismic coefficient method 2.2 Free – field deformation approach 2.3 Soil-structure interaction approach 3 pseudo static FEM 4 Full dynamic calculations using FEM 5 Case exemplification 6 Conclusions List of symbols Acknowledgements References