Acta Polytechnica CTU Proceedings doi:10.14311/APP.2018.15.0088 Acta Polytechnica CTU Proceedings 15:88–93, 2018 © Czech Technical University in Prague, 2018 available online at http://ojs.cvut.cz/ojs/index.php/app MODELING OF UNDERGROUND STRUCTURES SUBJECTED TO EARTHQUAKE BY COMBINING 1D FREE-FIELD AND 2D PSEUDO-STATIC ANALYSES Tereza Poklopová, Veronika Pavelcová, Tomáš Janda, Michal Šejnoha∗ Czech Technical University in Prague, Faculty of Civil Engineering, Thákurova 7, 166 29 Prague 6, Czech Republic ∗ corresponding author: sejnom@fsv.cvut.cz Abstract. The paper deals with the evaluation of tunnel construction subjected to earthquake using a pseudostatic analysis combined with the finite element method. The initial stage of calculation is concerned with the design of computational model and subsequent analysis of the actual excavation process. The case study of the response of a selected construction to earthquake follows next. To that end, the so called 1D free filed dynamic analysis is performed first to generate appropriate loading conditions in terms of a layered-wise constant shear strain. Therein, two particular boundary conditions, termed the fixed and absorbing boundary, are examined. The corresponding loading conditions are finally introduced in a 2D plane strain analysis to yield the internal forces developed in the tunnel lining. The results clearly show inadequacy of the fixed boundary and promote the use of absorbing boundary conditions for the present soil profile. Keywords: Earthquake, shear weave, accelerogram, finite element method, free-field, pseudostatic, absorbing boundary, fixed boundary, lining forces. 1. Introduction Design of underground structures loaded by earth- quake is typically based on analytical solutions uti- lizing pseudostatic methods. This may yield less ac- curate or rather conservative results. In addition, adopting simplified methods generally calls for a num- ber of restrictions such as the assumption of homo- geneous elastic medium and a circular or rectangular frame structure representing the tunnel. If these as- sumptions apply, one may easily follow approaches described, e.g. in [1, 2]. The effect of earthquake is then simply represented by a static load due to shear strain derived from a velocity of the propagating shear wave and the associated maximum velocity of vibrat- ing soil grains. This causes, e.g. ovaling or racking of a tunnel cross-section, which in turn serves as a basis for the stress analysis of the tunnel lining [2–4]. In this step, the soil-structure interaction is, however, neglected. Overcoming the above restrictions offers an appli- cation of the finite element method [5, 6]. From the reliability of structural design point of view it seems most appropriate the use of a fully dynamic analy- sis [7, 8]. Nevertheless, such an approach typically suffers from the lack of input parameters, notwith- standing the computational complexity. A suitable alternative, at least in the stage of a preliminary de- sign, is then seen in exploitation of the previously mentioned pseudostatic method built upon the use of FEM and one-dimensional (1D) fully dynamic free field analysis generating the load, see e.g. [3, 4]. This computational strategy will also be be examined in this paper with several recommendations for its prac- tical application. We open the discussion is Section 2 by briefly intro- ducing the assumed computational model and general steps preceding the pseudostatic analysis. The basic principles of free field analysis (FFA) will be presented next in Section 3. The actual pseudostatic analysis with emphases on the application of various types of boundary conditions is outlined in Section 4. The essential findings are finally summarized in Section 5. 2. Calculation steps To introduce the subject we first present two com- putational models of a circular tunnel built in a soft soil about 10m above the base rock. The material parameters of the two layers are listed in Table 1. While the response of clayey soil is assumed to be well represented by the Mohr-Coulomb (MC) material model, the limestone layer is considered as elastic. In [4] the material parameters of the MC model, not provided by the laboratory measurements including the increasing magnitude of the Young modulus with depth, were tuned to render the resulting distribution of the lining internal forces at the end of excavation similar to those provided by the Cam clay model. 2.1. Excavation step - selection of geometrical model No doubt that the response of soil subjected to earth- quake depends on the current state of stress generated 88 http://dx.doi.org/10.14311/APP.2018.15.0088 http://ojs.cvut.cz/ojs/index.php/app vol. 15/2018 Modeling of underground structures subjected to earthquake Soil layer E [MPa] ν [-] c [kPa] ϕ [◦] γ [kN/m3] Clayey soil 15-40 0.4 10 25 19.6 Limestone rock 630 0.3 - - 21.7 Table 1. Material parameters of subsoil layers. Element size Number of Settlement [mm] Lining forces Plastic strain [m] elements Terrain Crest Invert Nmax [kN/m] Mmax [kNm/m] [%] 0.25-2 4840 16.4 53.5 -26.8 952.8 129.8 0.91 0.5-4 2322 16.6 53.5 -26.3 958.4 135.5 0.85 0.5-8 1700 16.9 54.2 -26.6 948.6 138.8 0.86 1 - 8 1172 16.7 52.8 -24.7 935.2 165.6 0.59 Table 2. Evolution of settlement and lining internal forces as a function of mesh density. by the preceding actions. Here we limit our attention to the effect of excavation only. While the actual tunnel was built using the tunnel boring machine, the modeling of the excavation process exploited the clas- sical convergence confinement method (CC). To bring it closer to the, perhaps more suitable, volume loss method, the excavation sequence was proposed such as to generate final terrain settlement less than 2 cm. To that end, two geometrical models in Fig. 2 were created to address the influence of the rock layer in static analysis on the one hand, and, on the other hand, to allow for considering two particular boundary conditions used in the dynamic analysis discussed next in Section 3. (a) (b) Figure 1. Geometrical model: a) Model-1 - one layer, b) Model-2 - two layers. The excavation analysis was carried out in three consecutive steps to obtain the initial distribution of internal forces, see [4] for details. The evolution of settlement together with the evolution of inter- nal forces as a function of the mesh density were adopted in search for the optimal mesh. The results presented in Table 2, see also Fig. 2 for illustration, promoted the mesh with 6-node triangular elements of the element size ranging from 0.5 m to 8 m for the pseudostatic analysis. Point out that both mod- els in Fig. 1 provided essentially the same estimates of the initial distribution of internal forces thus giv- ing seemingly the preference to less computationally demanding Model-1. 2-0.254-0.58-0.5 8-1 0 20 40 60 80 100 120 140 160 180 0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 5500 M m ax (k N m /m ) n elements Figure 2. Evolution of lining moment as a function of mesh density. 2.2. Pseudostatic analysis - loading conditions As mentioned in the introductory part, the pseudo- static analysis considers the kinematic boundary con- ditions to generate the layer-wise uniform shear strain corresponding to free field conditions, the state prior to any construction. In case of a single wave propa- gating through a homogeneous material the free-field shear strain is provided by γ = vs cs , (1) where cs = √ G/ρ is the velocity of the shear wave and vs is the velocity of the soil particle; ρ,G stand for the soil density and elastic shear modulus, respectively. The particle velocity can be estimated based on the location and magnitude of earthquake, type of soil and the depth of tunnel [2]. However, the above assumptions are generally not met particularly if assuming a record of a real earth- quake such as the one in Fig. 3. It will be shown in the next section that in such a case it is more appro- priate to derive the shear strain directly from a fully dynamic analysis of a 1D free field column. 3. Free field dynamic analysis If a geotechnical construction can be analyzed in a two-dimensional (2D) environment under plane-strain conditions, the free field dynamic analysis reduces to the solution of a 1D column problem, see Fig. 4(a), 89 T. Poklopová, V. Pavelcová, T. Janda, M. Šejnoha Acta Polytechnica CTU Proceedings 0 5 10 15 20 25 30 35 40 45 50 a (m /s 2 ) t (s) 0.5 -0.5 -1.5 -2.5 1.5 2.5 (a) 0 5 10 15 20 25 30 35 40 45 50 v ( m /s ) t (s) 0.3 0.1 -0.1 -0.3 -0.5 0.5 (b) Figure 3. Adopted earthquake records: a) accelera- tion, b) velocity. having the same soil profile as the actual 2D com- putational model and being subjected to the same loading and boundary conditions. As it might not always be fulfilled by the 2D model, the free field conditions assume parallel and homogeneous layers in the horizontal direction. FIXED B.C. 𝑢(𝑦, 𝑡) = 𝑢(0, 𝑡) + 𝑢𝑅(𝑦, 𝑡) 𝑢𝐼 𝑢𝑂 𝑢𝑅(0, 𝑡) = 0 𝑢𝐼0 = 𝑢𝐼 + 𝑢𝑂 𝑦 ABSORBING B.C. 𝑢(𝑦, 𝑡) = 𝑢𝐼0(𝑡) + 𝑢𝑅(𝑦, 𝑡) 𝑢𝐼0 𝑢𝑂 𝑢𝑂(0, 𝑡) = 𝑢𝑅(0, 𝑡) 𝜌𝑐𝑠(𝑢�̇�(0, 𝑡) − 𝑢𝐼0̇ (𝑡)) (a) (b) rock outcrop motion on terrain: 𝑢(𝑦, 𝑡) = 2𝑢𝑢 𝑢𝑢 𝑢𝑑 𝑢𝑢 = 𝑢𝑑 𝑢𝐼0 = 2𝑢𝑢 fixed b.c. absorbing b.c. 𝑢𝑑 𝑢𝐼0 = 𝑢𝑢 (c) Figure 4. a) Free field column with fixed boundary, b) Free field column with absorbing boundary, c) Loading and boundary conditions. Turning now our attention to boundary conditions, two scenarios may occur. First, the soft soil layer containing a construction might be separated by an infinitely stiff base rock. The computational model is then truncated at the interface of the two layers and the resulting boundary is termed the fixed boundary. Such a model appears in Fig. 4(a) and corresponds to the 2D model in Fig. 1(a). However, we often encounter a situation where the model must be trun- cated within a layer. Such a boundary, depicted in Fig. 4(b), requires the outgoing wave uO to be trans- mitted across this boundary. Because the layer, con- taining the boundary, is assumed to be infinitely long, this requirement is essentially the same as if absorbing this wave at the boundary so it cannot be reflected back to the analyzed domain. This boundary will therefore be termed the absorbing boundary and will be linked to the 2D model in Fig. 1(b). Introduction of these two types of boundary conditions into the FEM analysis is outlined in the next section. As for loading, we consider a real acceleration in Fig. 3 in the form of a rock outcrop motion, see Fig. 4, where uu, ud correspond to the upward and downward traveling waves, respectively. It turns out, see also the next section, that in case of fixed boundary conditions the total motion1 will be adopted when employing the fixed boundary, whereas only one half of the total motion2, corresponding to the incoming wave uI , will be prescribed in case of absorbing boundary. 3.1. Theoretical background We begin by writing the total displacement field u(y, t) = uu(y, t) + ud(y, t) = uI(y, t) + uO(y, t) = uI0(t) + uR(y, t), (2) where uI0(t) = u(0, t) is the displacement field pre- scribed at the bottom boundary at y = 0, and uR(y, t) = {uR(y, t), vR(y, t)} is the searched rela- tive displacement vector. Going back to Fig. 1(a) it becomes clear that in case of fixed boundary it holds uR(0, t) = 0 at the bottom boundary. (3) The fixed boundary does not therefore require any special treatment. So we focus now on absorbing boundary and briefly describe the derivation of the absorbing boundary condition presented in Fig. 4(b). For more details the interested reader is referred to [9]. Examining Eq. (2) suggests u ′ O(y, t) = u ′ R(y, t)− u ′ I(y, t), (4) where (f ′) stands for the space derivative with respect to y. As shown in [9, 10] we may further write u ′ O(y, t) = 1 cs u̇O(y, t) (5) u ′ I(y, t) = − 1 cs u̇I(y, t), (6) 1The actual magnitude of the record being equal to some of the incoming and outgoing waves as the outgoing wave is reflected back to the domain. 2This assumption is valid for the case study in Fig. 4(b). 90 vol. 15/2018 Modeling of underground structures subjected to earthquake where (ḟ) represents the time derivative. Substituting Eqs. (5) and (6) into Eq. (4) gives for y = 0, see also Fig. 4(b), u ′ O(0, t) = u ′ R(0, t) + 1 cs u̇I0(t) = 1 cs u̇R(0, t). (7) Considering the shear wave only, Eq. (7) combined with the Hook law provides the static boundary con- dition in the form τ(0, t) = Gu ′ R(0, t) = ρcs (u̇R(0, t)− u̇I0(t)) (8) to be satisfied at the bottom boundary. Note that this condition corresponds to a viscous damper with the viscosity η = ρcs introduced at y = 0, recall Fig. 4(b). In the absence of material damping the principal of virtual work reads δüT RMüR + δuT RKuR + δu̇T RCu̇R|y=0 = −δüT RMüI0 + δu̇T RCu̇I0|y=0. (9) where M,K represent the mass and stiffness matrices and the damping matrix C, active at the bottom boundary only, is provided by C = [ ρcs 0 0 ρcp ] , (10) where cp = √ Eoed/ρ is the pressure wave velocity with Eoed being the oedometric modulus. The discretized system of governing equations of motion resulting from Eq. (9) is typically solved with the help of Newmark method [5]. 3.2. Loading conditions - shear strain The time variation of the shear strain u ′(y, t) was found by employing Eq. (9). To that end, the accel- eration üI0 and velocity u̇I0 records in Fig. 3 were adopted. Note that in case of fixed boundary the damping term in Eq. (9) drops out. It is also worth mentioning that Eq. (8) was de- rived assuming the elastic behavior. In this regard, the elastic moduli G,Eoed corresponding to an un- loading/reloading modulus should be employed. The deformation modulus for the clayey soil in Table 1 was therefore replaced by a dynamic modulus and set equal to 80 MPa for the entire layer. In the light of subsequent analysis the maximum values of the shear strain u′(y, t) developed along the vertical direction of the analyzed domain during the duration of earthquake were collected. A graphical representation for the two types of boundary condi- tions appears in Fig. 5. The corresponding particle velocity, to be potentially introduced in Eq. (1), is plotted for the sake of completeness. We immediately observe that the maximum shear strain is in case of fixed boundary considerably bigger in comparison to absorbing boundary. This result, attributed to the interference of the incoming wave and the one being reflected back on the fixed boundary (infinitely stiff base rock) will be analyzed in more details in the next section. The resulting kinematic boundary conditions for the two models are plotted in Fig. 6. Model-1 FIXED B.C. -49 -39 -29 -19 -9 z (m ) γmax (-) γmean = 0.023 0.0 0.01 0.02 0.03 0.04 0.05 -49 -39 -29 -19 -9 z (m ) vs,max (m/s) vs,mean = 2.802 m/s 0.0 1.0 2.0 3.0 4.0 5.0 tunnel axis (a) Model-2 -59 -49 -39 -29 -19 -9 z (m ) γmax (-) γmean= 3.541∙10-3 γmean= 0.630∙10-3 0.0 0.002 0.004 0.006 0.008 -59 -49 -39 -29 -19 -9 z (m ) vs,max (m/s) vs,mean= 0.335 m/s vs,mean= 0.549 m/s 0.00 0.20 0.40 0.60 0.80 1.00 tunnel axis ABSORBING B.C. (b) Figure 5. Distribution of maximum shear velocity and shear strain: a) fixed boundary, b) absorbing boundary. 4. Pseudostatic analysis This section reports on the application of simplified pseudostatic analysis to two geometrical models in Fig. 1 being linked two types boundary conditions described in the previous section. The pseudostatic analysis represents a new calculation stage in which the analyzed domain is loaded by the prescribed dis- placements compatible with the average shear strain provided by the free field analysis, see Fig. 6. γ = 0.023 1 127.0 mm γ = 0,0035 γ = 0.00063 179.8 mm 6.3 mm γ = 0.0035 (a) (b) Figure 6. Kinematic loading conditions: a) fixed boundary, b) absorbing boundary. The value of the maximum displacement associated with the fixed boundary immediately suggests inade- quacy of this type of boundary condition in the present analysis. Herein, the Young modulus of the base rock equal to 630 MPa, recall Table 1, is far from infin- ity particularly if compared to the value of 80 MPa adopted for the clayey soil. This can be supported by calculating the value of particle velocity following 91 T. Poklopová, V. Pavelcová, T. Janda, M. Šejnoha Acta Polytechnica CTU Proceedings the procedure outlined in [2]. In our particular case, see [4] for details, the particle velocity vs would be approximately equal to 0.99 m/s. This is about three times bigger than the average velocity reported in Fig. 5(a) for the fixed boundary and comparable with the one found for the absorbing boundary in Fig. 5(b). The suitability of using the absorbing boundary con- ditions will further be supported by inspecting the values of lining internal forces in the next section. But before proceeding, a word of caution is in place as to the application of material models in pseudo- static analysis. This becomes clear by examining the results plotted in Fig. 7. These were derived by con- sidering the computational Model-1 with only one layer of soil free of any tunnel (free field conditions) and loaded by the prescribed constant shear strain γ = 0.023, recall Fig. 5(a). To be consistent with the free field conditions the same value of shear strain should be reproduced at every point within the do- main3. This is certainly not true in the first case, see Figs. 7(a,b) assuming the soil with a variable stiffness and its response driven by the MC material model. As expected such a loading may generated an unaccept- able plastic strains at the vicinity of terrain surface where the compressive stresses are close to zero. Thus the nonlinear constitutive model, if at all, should only be applied in the vicinity of the tunnel. As further seen in 7(c) the use of elastic model while still keeping the variable stiffness assumed in excavation analysis does not provide an acceptable result. As evident from Fig. 7(d), that is only found if considering a truly homogeneous (one modulus for the entire layer) elastic soil, thus the material setting used previously in FFA. (a) (b) (c) (d) Figure 7. a) Total shear strain - variable stiffness and MC model, b) Equivalent plastic strain - variable stiffness and MC model, c) Total shear strain - variable stiffness and elastic model, d) Total shear strain - constant stiffness with the layer and elastic model. Such conditions then yield the expected distribution of shear strains for the two computational models 3This is not exactly true as the static boundary condition requires the shear stress and thus also the shear strain be zero at the terrain boundary [9]. as evident from Fig. 8. Clearly, the perturbation caused by the presence of tunnel is confined to its close vicinity. γmax = 4.04% γmax = 0.64% (a) (b) Figure 8. Distribution of shear strain: a) Model- 1, loading from Fig. 6(a), b) Model-2, loading from Fig. 6(b). 4.1. Lining forces due to earthquake To further appreciate the differences in the response provided by the two types of boundary conditions we examine the resulting values of the lining internal forces at the end of the pseudostatic analysis. For illustration we plot in Fig. 10(a) the distribution of mo- ments found at the end of excavation step, thus prior to the introduction of surcharge due to earthquake. The impact of earthquake, in the light of pseudostatic analysis, is clearly visible in Figs. 10(b,c). -1500 -1000 -500 0 500 1000 1500 0 π/4 π/2 3π/4 π 5π/4 3π/2 7π/4 2π M ( k N m /m ) ϕ (rad) End of excavation After pseudostatic analysis 1380.8 1393.8 1384.2 1393.6 65.0 90.4 64.8 90.3 90.3 (a) -300 -200 -100 0 100 200 300 0 π/4 π/2 3π/4 π 5π/4 3π/2 7π/4 2π M ( k N m /m ) ϕ (rad) End of excavation After pseudostatic analysis 229.5 231.3 231.9 230.2 90.3 65.0 90.4 90.3 64.8 (b) Figure 9. Distribution of moments prior to and after running the pseudostatic analysis: a) fixed boundary, b) absorbing boundary. Even more clear illustration of the influence of the selected type of boundary is seen in Figs. 9 showing variation of moments as a function of the circumfer- ential angle. The use of fixed boundary conditions is clearly in- correct as the maximum moment in Fig. 9(a) exceeds 10 times the value of the moment at the end of exca- vation step. The normal force, not shown, becomes even positive at some locations with the largest value over 1600 kN/m. Given the depth of the tunnel, such 92 vol. 15/2018 Modeling of underground structures subjected to earthquake -90,4 6 4 ,8 -90,3 90,3 -6 5 ,0 90,4 -85,8 -1384,2 13 93 ,6 -90,5 90,5 1380,8 -1 39 3, 8 85,8 -90,4 -231,9 23 0, 2 -88,8 88,8 229,5 -2 31 ,3 90,4 (a) (b) (c) Figure 10. Distribution of moments: a) at the end of excavation step, b) at the end of pseudostatic analysis - fixed boundary, c) at the end of pseudostatic analysis - absorbing boundary. a result is highly improbable. Much reasonable re- sults are achieved with the use of absorbing boundary as displayed in Fig. 9(b) completely eliminating the positive normal forces, see [4] for further details. The resulting distributions of moment increments caused by earthquake and pertinent to two types of boundary conditions are finally compared in Fig. 11 further highlighting the inadequacy of the fixed bound- ary conditions, at least for the present geometrical and material settings. -1500 -1000 -500 0 500 1000 1500 0 π/4 π/2 3π/4 π 5π/4 3π/2 7π/4 2π Δ M ( k N m /m ) ϕ (rad) Model-1 Model-2 216.1 221.4 215.8 216.6 1385.5 1388.5 1383.3 1384.0 Figure 11. Distribution of moments comparing the impact of fixed and absorbing boundary conditions. 5. Conclusions The present paper was concerned with the application of pseudostatic analysis to address the effect of earth- quake on the evolution of internal forces within the lining of a circular tunnel. Unlike in most practical approaches the loading conditions for the final 2D analysis were found from a fully dynamic analysis of a 1D free field column problem subjected to a real earthquake. In this regard, two types of boundary conditions prescribed at the bottom boundary were examined. While fixed boundary is typically adopted in cases when a soft soil layer is supported by an in- finitely stiff base rock, the absorbing boundary must be introduced if the depth of computation model ap- pears within a layer, regardless whether soft or stiff. Such a boundary condition then allows for completely absorbing the incoming wave so it is not reflected back into to the analyzed domain as in case of fixed boundary. Unless the based rock is truly stiff, the fixed boundary should be avoided, as it leads, unlike the absorbing boundary, to unacceptable results. This was shown herein for one particular example. Acknowledgements The support provided by the SGS project No. SGS18/037/OHK1/1T/11 and the TAČR project no. TE01020168 is gratefully acknowledged. References [1] J. Penzien. Seismically induced racking of tunnel linings. Earthquake engineering and structural dynamics 29:683–691, 2000. [2] Y. Hashash, J. Hook, B. Schmidt, J. Yao. Seismic design and analysis of underground structures. Tunneling and undergrounding Space Technology 16:247–293, 2001. [3] D. Kučera. Using finite element method in practical design of geotechnical structure under seismic loading. Master’s thesis, Czech Technical University in Prague, Faculty of Civil Engineering, 2017. In Czech. [4] T. Poklopová. Evaluation of real underground structure subjected to earthquake - pseudostatic analysis. Bachalor’s thesis, Czech Technical University in Prague, Faculty of Civil Engineering, 2018. In Czech. [5] O. Řeřicha, Z. Bittnar. Finite element method in dynamic analysis of structures. Nakladatelství technické literatury, 1981. In Czech. [6] Z. Bittnar, J. Šejnoha. Numerical methods in structural engineering. ASCE Press, 1996. [7] 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. [8] V. Pavelcová, T. Poklopová, T. Janda, M. Šejnoha. The influence of boundary conditions on the response of underground structures subjected to earthquake. Acta Polytechnica 2018. Accepted. [9] O. Zinkiewicz, N. Bicanic, F. Shen. Generalized Smith boundary - a transmitting boundary for dynamic computations, vol. 207. INNME Swansea, 1986. [10] D. Kučera. Posouzení geotechnické konstrukce na seismické zatížení metodou konečných prvku. Master’s thesis, Czech Technical University in Prague, Faculty of Civil Engineering, 2017. 93 Acta Polytechnica CTU Proceedings 15:88–93, 2018 1 Introduction 2 Calculation steps 2.1 Excavation step - selection of geometrical model 2.2 Pseudostatic analysis - loading conditions 3 Free field dynamic analysis 3.1 Theoretical background 3.2 Loading conditions - shear strain 4 Pseudostatic analysis 4.1 Lining forces due to earthquake 5 Conclusions Acknowledgements References