Acta Polytechnica CTU Proceedings doi:10.14311/APP.2020.26.0064 Acta Polytechnica CTU Proceedings 26:64–70, 2020 © Czech Technical University in Prague, 2020 available online at https://ojs.cvut.cz/ojs/index.php/app INFLUENCE OF MATERIAL DAMPING ON THE RESPONSE OF UNDERGROUND STRUCTURES SUBJECTED TO EARTHQUAKE Veronika Pavelcová, Tereza Poklopová, Michal Šejnoha∗, Tomáš Janda Czech Technical University in Prague, Faculty of Civil Engineering, Department of Mechanics, Thákurova 7, 166 29 Prague 6, Czech Republic ∗ corresponding author: sejnom@fsv.cvut.cz Abstract. The paper describes a finite element simulation of the response of a real underground structure subjected to earthquake using GEO5 FEM program. It concentrates on the influence of material damping with respect to a specific type of boundary condition prescribed at the bottom of the analyzed domain. It is seen that considering material damping is inevitable particularly in case of so called fixed boundary conditions to arrive at meaningful results. This is demonstrated on an artificial earthquake generated according to a design spectrum defined in Eurocode 8. A viscous damping ratio combined with the results of eigenvalue analysis is used to derive parameters of Rayleigh damping for three specific scenarios promoting the approach based on the lowest natural frequency as sufficiently accurate for the present task. Keywords: Earthquake, tunnel, material damping, accelerogram, finite element method, absorbing boundary, fixed boundary. 1. Introduction Numerical modeling of geotechnical structures sub- jected to earthquake typically distinguishes between damping associated with the application of radiation boundary condition [1] and intrinsic material damp- ing associated with the energy dissipation within the material caused by various mechanisms including inter- particle sliding and friction and pore fluid viscosity [2]. The former type of damping has been addressed in detailed in [3, 4]. The present contribution is therefore an extension of our previous work with emphasis on the effect of material damping. Although valid for a single degree of freedom the complex structural analysis usually grounds on the knowledge of a single value damping ratio. Such a parameter is typically estimated experimentally from a hysteric stress-strain loop as the energy loss in a cycle [5]. While some new measurement technique arise [6] to eliminate uncertainty in the determination of geometric damping thus increasing reliability in the damping ratio measurement, the determination of damping ratio for all mode shapes is still far from trivial. This is why most practical applications rely on a single value damping ratio even in case of multi degree of freedom tasks combined with the concept of Rayleigh proportional damping [7]. This computa- tional strategy is adopted also in the present study focused on the effect of material damping in com- bination with various types of boundary conditions introduced at the bottom of the computational model. The remainder of the paper is organized as follows. Starting with a theoretical background in Section 2 we continue in Section 3 with the definition of material damping. The theoretical part of the paper is finalized in Section 4 addressing the construction of artificial accelerograms. Details of the computational models including the finite element mesh refinement are pro- vided in Section 5. The overview of computational results starts in Section 6 with modal analysis to pro- vide data for the derivation of parameters of Rayleigh damping. The results of free field (FF) analysis are presented in Section 7 followed by the summary of the essential outcomes delivered by two-dimensional (2D) analyses reflecting the presence of a tunnel. 2. Governing equations A detailed derivation of the set of governing equations driving the earthquake analysis in geotechnical engi- neering can be found, e.g., in [1, 3, 8]. Therein, the formulation was limited to damping associated with absorbing boundary conditions. When taking into account the material damping the discretized form of the equation motion with quiet boundary assumed at the bottom of the model be- comes, see also Fig. 1, MüR + Cu̇R +KuR + CBBu̇R|y=0 + CLBu̇R|x=0,L = −MüI0 − Cu̇I0 + CBBu̇I0|y=0 + CLBu̇FFR |x=0,L (1) −Rτ |x=0 +Rτ |x=L, assuming the total displacement u(x, t) be split into a relative displacement uR(x, t) and the displacement uI0(t) applied to the whole domain such that u(x, t) = uR(x, t) + uI0(t). (2) In Eq. (1), M,K are the mass and stiffness matrices of the domain, respectively. The damping matrix C 64 https://doi.org/10.14311/APP.2020.26.0064 https://ojs.cvut.cz/ojs/index.php/app vol. 26/2020 Influence of material damping in earthquake analysis Figure 1. Boundary and loading conditions for both FF and 2D analysis. corresponds to material damping only, while matri- ces CBB and CLB reflect the effect of artificial viscous dampers applied along the absorbing boundaries (the superscripts BB and LB identify the bottom and lat- eral boundaries; henceforth the absorbing BB bound- ary will be called the quiet boundary and denoted as QBB whereas the notation FBB will be adopted to identify the so called fixed boundary [3, 8]). Finally, the vector Rτ stores the nodal forces associated with the prescribed boundary shear stress. Recall that in case of QBB the vector uI0(t) stands for the incoming wave only. Should the top soil be supported by an infinitely stiff bedrock we may in- troduce at the bottom boundary the FBB conditions suggesting a total reflection of the incoming wave back into the domain. The applied load than corresponds to the total motion u0(t) so that u(x, t) = uR(x, t) + u0(t), (3) and Eq. (1) simplifies as MüR + Cu̇R +KuR + CLBu̇R|x=0,L = −Mü0 − Cu̇0 + CLBu̇FFR |x=0,L (4) −Rτ |x=0 +Rτ |x=L, For simplicity, the present study assumes the incoming wave be equal to one half of the total motion. 3. Rayleigh damping Material damping may have a strong impact on the response of an underground structure when subjected to earthquake. However, a physically supported for- mulation of damping matrix C is not an easy task mainly because of the lack of experimental evidence. The proportional Rayleigh damping is then a typical choice in most practical applications. In such a case, the damping matrix is expressed as a linear combina- tion of the mass and stiffness matrices in the form, see e.g. [9, 10], C = αM + βK. (5) Accepting Eq. (5) allows us to write the damped natural frequency ωdi associated with the i-th orthog- onal mode shape φi of undamped system in terms of the viscous damping ratio ξi and the i-th natural frequency of undamped system ωi 2ωdi = 2ξiωi = φi TCφi = α+ βω2 i , (6) because φi TMφi = 1, φiTKφi = ω2 i . (7) Thus if two viscous damping ratios ξi are known for the pair of natural frequencies ωi, the values of param- eters α, β readily follow from Eq. (6). However, most often only a single value of parameter ξ is provided. Then, the hypothetical assumption of the first natu- ral frequency being least damped is usually accepted. Therefore, dξ dω = 1 2 ( − α ω2 + β ) = 0. (8) Substituting ω = ω1 into Eqs. (6) and (8) then gives α = ξ1ω1, β = ξ1 ω1 . (9) 4. Synthetic accelerograms A finite element analysis of geotechnical structures subjected to earthquake-induced loading requires spec- ifying the time history of free field accelerations. Eu- rocode 8 [11] allows for using either recorded or ar- tificial accelerograms. The artificial accelerograms have to be generated to match the elastic response spectrum defined within Eurocode 8. An example of such a spectrum used in the present study is plotted in Fig. 2. The algorithm for generating artificial accelero- grams utilized in this paper is mostly inspired by the algorithm proposed in [12] and consists of the following steps: 65 V. Pavelcová, T. Poklopová, M. Šejnoha, T. Janda Acta Polytechnica CTU Proceedings 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 Time period T [s] 0.5 1.0 1.5 2.0 2.5 3.0 3.5 A cc el er at io n a [m / s2 ] ω1ω3 ωRS Design response spectrum Computed response spectrum Figure 2. Design response spectrum, computed re- sponse spectrum, location of frequencies to calculate parameters of Rayleigh damping. (1.) A Fourier spectrum with constant spectral am- plitudes and random phases is generated. (2.) Time history of acceleration is obtained from the spectrum by means of Fourier transform. (3.) Elastic responses of a single degree-of-freedom systems with natural periods corresponding to fre- quencies used in the Fourier spectrum are computed for the time history. (4.) Factors of the computed elastic response accel- eration and the design response acceleration are evaluated. (5.) The spectral amplitudes of the original Fourier spectrum are divided by the factors while the phases are kept. (6.) Steps 2-5 are repeated for the updated Fourier spectrum until the computed elastic response spec- trum matches the design elastic response spectrum within 10% relative error, see Fig. 2. Although the generated signal matches the target elastic response spectrum, it is stationary and lacks the characteristic stages of a typical earthquake signal. To take the rise, the strong motion and the decay periods into account, the amplitude of the stationary signal is modulated by an envelope function [13] E(t) = atb exp(−ct), (10) with parameters b = −ε ln η 1 + ε(ln ε− 1) , (11) c = b εTw , (12) a = ( exp(1) εTw )b , (13) where Tw denotes the duration of earthquake. Peak of the envelope function occurs at fraction ε of the earthquake duration and the original amplitude is reduced to η fraction at time TW . An example of 0.0 2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0 Time t [s] −1.5 −1.0 −0.5 0.0 0.5 1.0 A cc el er at io n a [m / s2 ] Stationary Modulated Figure 3. An example of stationary (original) and modulated (with rise, strong motion and decline stages) accelerograms. stationary and modulated synthetic accelerogram is shown in Figure 3. 5. Model description Two computational models were created to demon- strate the effect of material damping in relation to the selected boundary conditions applied at the bottom of the domain. The first model in Fig. 4(a) consists of a soft layer supported by an infinite layer of a relatively stiff bedrock. Truncating the model within the rock layer thus requires application of absorbing bound- ary conditions ensuring that the incoming wave is fully damped and not reflected back into the domain. The second model in Fig. 4(b) assumes in turn the bedrock be infinitely stiff thus replacing this layer with the fixed boundary. Both layers are assumed linearly elastic with material properties listed in Table 1. Soil Density ρ 1960 [kg/m3] Poisson ratio ν 0.40 [-] Elastic modulus E 80.00 [MPa] Rock Density ρ 2170 [kg/m3] Poisson ratio ν 0.30 [-] Elastic modulus E 630.00 [MPa] Table 1. Material parameters. As in [3] the tunnel being 6 m in diameter is located in the soil layer at a depth of 37.35 m. The tunnel lining is modeled by linearly elastic beam elements. 6. Rayleigh damping parameters As demonstrated in Section 3, defining parameters α, β in Eq. (5) requires the knowledge of a pair of natural frequencies, Eq. (6), or at least of the first natural frequency ω1 should Eq. (9) be used. 66 vol. 26/2020 Influence of material damping in earthquake analysis Point B Point A2Point A1 0,00 -6,00 -12,00 -18,00 -24,00 -30,00 -36,00 -42,00 -48,00 -54,00 -59,00 0, 00 -1 2, 00 -6 ,0 0 -1 8, 00 -2 4, 00 -3 0, 00 -4 0, 00 36 ,0 0 24 ,0 0 12 ,0 0 30 ,0 0 40 ,0 0 6, 00 (a) Point B Point A2Point A1 0,00 -6,00 -12,00 -18,00 -24,00 -30,00 -36,00 -42,00 -49,00 0, 00 -6 ,0 0 -1 2, 00 -1 8, 00 -2 4, 00 -3 0, 00 -4 0, 00 -6 ,0 0 -1 2, 00 -1 8, 00 -2 4, 00 -3 0, 00 -4 0, 00 36 ,0 0 24 ,0 0 12 ,0 0 30 ,0 0 40 ,0 0 6, 00 (b) Figure 4. Computational model - a) quiet bound- ary conditions (QBB), b) fixed boundary conditions (FBB). Eurocode 8 suggests proportional damping with a single value of ξ = 0.05 (proportional damping of 5%). Thus the way of calculating parameters α, β deter- mines, how much will individual natural frequencies be damped. The following three options are examined in this study • Type 1: The first natural frequency is assumed to be the least damped frequency. • Type 2: The least damped frequency domain is be- tween the first and the third natural frequency. This option results from the recommendation provided in the PLAXIS user manual [14]. • Type 3: The least damped frequency domain is between the first natural frequency and the most dominant frequency ωRS = 12.5664 rad/s of the design response spectrum displayed in Fig. 2. The required natural frequencies were found from a standard eigenvalue analysis. The first three natu- ral frequencies corresponding to the shear dominated mode shapes were extracted from the list of calculated frequencies based on the value of a significant modal participation factor [15] for the x direction, recall the coordinate system displayed in Fig. 1. Three analyses in particular were carried out. The first one assumed the free field column (1DFF), the second one adopted the 2D model in Fig. 4(a) (2DT) and the third one considered the same model but for the state prior to excavation (2DI). In all calculations, the bottom boundary was fixed. In 2D models the lateral bound- aries were additionally constrained in the y-direction to recover the mode shapes compatible with the free field column analysis. The results are stored in Ta- ble 2. The values in parentheses represent the actual order when considering all modes. ω1 ω2 ω3 1DFF 3.7547 (1) 11.5307 (3) 20.1191 (4) 2DT 3.7481 (1) 11.4440 (5) 19.5180 (10) 2DI 3.7505 (1) 11.3844 (5) 19.4204 (10) Table 2. Calculated natural frequencies. As can be seen, all analyses provided almost iden- tical results. Thus to determine parameters α, β in Table 3 the results from 2DI analysis were employed. Damping ω ξ α β Type 1 ω1 0.05 0.1875 0.0133 Type 1b ω1 0.10 0.3750 0.0266 Type 2 ω1 + ω3 0.05 0.3143 0.0043 Type 3 ω1 + ωRS 0.05 0.2888 0.0061 Table 3. Rayleigh damping parameters. Figure 5 shows how the natural frequencies are damped for the three types of damping with ξ = 0.05. Figure 5. Viscous damping ratio as a function of the natural period (frequency). 7. 1D free field analysis The influence of material damping was first addressed in the light of one-dimensional free-field column anal- ysis. This is in fact also the first step in a two- dimensional analysis needed in generating the time variation of free-field velocities and boundary shear stresses. It also allows us to test the influence of finite element mesh, which is in this case quite non-uniform as seen in Fig. 4. The 1DFF mesh corresponds to the refinement along the edges of the 2D domains and certainly is not optimal neither from the material, nor from the integration time step point view, which was 67 V. Pavelcová, T. Poklopová, M. Šejnoha, T. Janda Acta Polytechnica CTU Proceedings set to 0.01 s in accordance with the adopted accel- eration record. Only the soil layer was assumed to be damped. All the results presented in Figs. 6 - 9 show variation of a relative displacement along the x-direction uRx observed on the terrain (Point A1). 0 5 10 15 20 Time [s] -0,12 -0,08 -0,04 0,00 0,04 0,08 0,12 u R x [ m ] No damp. Damp.T 1 ACC Fixed BC Point A1 Figure 6. Time variation of horizontal displacement uRx monitored at Point A1. Blue line is the scaled variation of the prescribed acceleration. Figure 6 compares variation of damped and un- damped displacement uRx considering the fixed bound- ary conditions, Fig. 4(b). The scaled acceleration is added just to illustrate the time interval of strong mo- tion. The results suggest that appropriate damping is absolutely essential in case of fixed boundary, because in undamped model the vertically propagating shear waves reflect not only from the terrain free boundary but also from the fixed bottom boundary. Therefore, in the undamped model the motion does not cease even after the input signal diminishes. The influence of the selected type of damping, recall Section 6 and Table 3, is examined in Figs. 7 and 8. While some differences in the model response can be observed in Fig. 8 in the case of QBB, the results plot- ted in Fig. 7 considering FBB suggest almost negligible effect of the type of damping, at least in the present example. It is also seen in Fig. 8 that unlike the fixed boundary, the influence of damping in the case of QBB is far less pronounced. This is attributed to the fact that the downward traveling wave is intentionally damped at the bottom boundary having the same effect as allowing the downward traveling waves to leave the examined model at the truncated boundary, so the motion eventually stops even without material damping. To directly compare the results derived for FBB and QBB we first plot in Fig. 9 the time variation of the horizontal displacement uRx monitored at various points to see that the displacement at the interface of the two materials practically coincides with the motion monitored at the bottom boundary. This indicates that the vertically propagating wave travels very fast through the rock layer so the difference between the two models dimensions in Fig. 4 can safely be neglected. This comparison of the response of damped model with FBB and undamped model with QBB is plotted in Fig. 10. 0 5 10 15 20 Time [s] -0,10 -0,05 0,00 0,05 0,10 0,15 u R x [ m ] No damping Damp. type 1 Damp. type 2 Damp. type 3 Point A1 Fixed BC Figure 7. Influence of the type damping - FBB. 0 2 4 6 8 10 Time [s] -0,04 -0,02 0,00 0,02 0,04 0,06 u R x [ m ] No damping Damping type 1 Damping type 2 Damping type 3 Point A1 Quiet BC Figure 8. Influence of the type damping - QBB. It can be observed that even with 5% damping, the model with fixed bottom boundary yields non- negligible displacements even beyond 20 s of monitor- ing where the prescribed motion comes to zero. The figure also shows the importance of correctly prescrib- ing only the incoming wave for the model with QBB. Adopting the total motion would yield the displace- ment at the terrain exceeding the motion found in the model with FBB, which is unrealistic. 8. 2D analysis of tunnel Taking into account the results discussed in the previ- ous section the 2D analyses were limited to the QBB model in Fig. 4(a) assuming no damping, while only the damping type 1 was considered when examin- ing the FBB model in Fig. 4(b). For both models, the horizontal displacements plotted in Fig. 11 were monitored at points A2 and B. Similar conclusions to what has been observed in Section 7 can be drawn here suggesting again the need for material damping with fixed boundary. Because choosing the particular value of damping ratio has mostly no rigorous support, we compared the results pertinent to 5% of damping to those derived for 10% of damping. Again the 5% of damping still overes- timates the displacements at the end of earthquake, while the results found for 10% of damping resemble the response of QBB model at least at the terrain boundary. In the vicinity of tunnel lining, the motion is, however, overdamped. 68 vol. 26/2020 Influence of material damping in earthquake analysis 0 5 10 15 20 Time [s] -0,04 -0,02 0,00 0,02 0,04 0,06 u R x [ m ] Terrain Soil layer mid point Interface Bottom No damping Point A1 Figure 9. Time variation of horizontal displacement uRx monitored at various points - QBB. 0 5 10 15 20 Time [s] -0,10 -0,05 0,00 0,05 0,10 0,15 u R x [ m ] Fixed BC (Damping type 1) Quiet BC (No damp. total motion) Quiet BC (No damp. incom. wave) Point A1 Figure 10. Comparison of damped model with FBB and undamped model with QBB. A smooth variation of damped displacements also suggests that particularly high frequencies of the mo- tion are damped, which agrees well with Fig. 5. An- other observation follows from a close inspection of Figs. 10 and 11(a) showing only a minor difference between the response provided by 1DFF and 2DT analyses, mainly because of a relatively stiff lining. This thus supports the results obtained already from modal analysis in Section 6. 9. Conclusions The present paper examined the influence of mate- rial damping in earthquake analysis particularly with reference to the boundary conditions prescribed at the bottom boundary of the model. The following conclusions can be proposed: • Considering material damping in case of QBB has only a minor effect on the predicted response. • On the other hand, the FBB model should not be considered without material damping as rather er- roneous results can be obtained. • The way of determining the Rayleigh damping pa- rameters has, at least in the present study, a rela- tively minor effect. Acknowledgements The support provided by the SGS project No. SGS19/032/OHK1/1T/11 and the TAČR project No. 0 5 10 15 20 Time [s] -0,10 -0,05 0,00 0,05 0,10 u R x [ m ] Quiet BC - ξ=0% Fixed BC - ξ=5% Fixed BC - ξ=10% Point A2 (a) 0 5 10 15 20 Time [s] -0,03 -0,02 -0,01 0,00 0,01 0,02 u R x [ m ] Quiet BC - ξ=0% Fixed BC - ξ=5% Fixed BC - ξ=10% Point B (b) Figure 11. Time variation of relative displacement at selected points during the period of strong motion: a) point A, b) point B. TE01020168 is gratefully acknowledged. References [1] O. Zinkiewicz, N. Bicanic, F. Shen. Generalized Smith boundary - a transmitting boundary for dynamic computations, vol. 207. INNME Swansea, 1986. doi:10.1007/978-3-7091-2828-2_3. [2] A. Ashmawy, R. Salgado, S. Guha, V. Drnevich. Soil damping and its use in dynamic analyses. In Third International Conferences on Recent Advances in Geotechnical Earthquake Engineering and Soil Dynamics. Missouri University of Science and Technology, 1995. Paper No. 1.13. [3] V. Pavelcová, T. Poklopová, T. Janda, M. Šejnoha. The influence of boundary conditions on the response of underground structures subjected to earthquake. Acta Polytechnica, CTU Proceedings 15:74–80, 2018. doi:10.14311/APP.2018.15.0074. [4] J. Pruška, M. Šejnoha, T. Janda, et al. Analysis of underground structures subjected to earthquake using finite element method. In 14th International Conference entitled "Underground Construction Prague 2019". ITA-AITES Czech Tunnelling Association, 2019. June 3-5, Prague. [5] I. Towhata. Geotechnical Earthquake Engineering. Springer-Verlag Berlin Heidelberg, 2008. [6] S. Badsar, M. Schevenels, W. Haegeman, et al. Determination of the material damping ratio in the soil from SASW tests using the half-power bandwidth method. Geophysical Journal International 182:1493– 1508, 2010. doi:10.1111/j.1365-246X.2010.04690.x. 69 https://doi.org/10.1007/978-3-7091-2828-2_3 https://doi.org/10.14311/APP.2018.15.0074 https://doi.org/10.1111/j.1365-246X.2010.04690.x V. Pavelcová, T. Poklopová, M. Šejnoha, T. Janda Acta Polytechnica CTU Proceedings [7] A. Zidan. Modeling of soil damping for seismic ground response by nonlinear finite element analysis. In Third International Conference on Advances in Civil and Structural Engineering - CSE. Institute of Research Engineers and Doctors, USA, 2015. [8] V. Pavelcová. Evaluation of real underground structure subjected to earthquake - fully dynamic analysis. Bachalor’s thesis, Czech Technical University in Prague, Faculty of Civil Engineering, 2018. In Czech. [9] O. Řeřicha, Z. Bittnar. Metoda konečných prvků v dynamice konstrukcí. Nakladatelství technické literatury, 1981. [10] Z. Bittnar, J. Šejnoha. Numerical methods in structural engineering. ASCE Press, 1996. [11] Comité Européen de Normalisation. Eurocode 8: Design of structures for earthquake resistance, 2004. [12] A. Kumar. Software for generation of spectrum compatible time history. In Proceedings of 13th World Conference on Earthquake Engineering. Canadian Association for Earthquake Engineering, 2004. [13] D. Boore. Stochastic simulation of high frequency ground motion based on seismological models of radiated spectra. Bulletin of Seismological Society of America 73(6):1865–1894, 1983. [14] PLAXIS manual, http://www.plaxis.nl/. [15] T. Irvine. Effective modal mass and modal participation factors Revision H, 2013. 70 Acta Polytechnica CTU Proceedings 26:64–70, 2020 1 Introduction 2 Governing equations 3 Rayleigh damping 4 Synthetic accelerograms 5 Model description 6 Rayleigh damping parameters 7 1D free field analysis 8 2D analysis of tunnel 9 Conclusions Acknowledgements References