Acta Polytechnica CTU Proceedings doi:10.14311/APP.2018.15.0074 Acta Polytechnica CTU Proceedings 15:74–80, 2018 © Czech Technical University in Prague, 2018 available online at http://ojs.cvut.cz/ojs/index.php/app THE INFLUENCE OF BOUNDARY CONDITIONS ON THE RESPONSE OF UNDERGROUND STRUCTURES SUBJECTED TO EARTHQUAKE Veronika Pavelcová, Tereza Poklopová, 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 prediction of the response of a real underground structure subjected to earthquake. A fully dynamic analysis is carried out in the GEO5 FEM program using the Finite Element Method. Limiting our attention to a two-dimensional analysis we focus on the implementation of special boundary conditions along the vertical edges of the computational model. A simple study is carried out first to show that incorrectly applied boundary conditions may significantly influence the actual design of underground structures loaded by vertically propagating shear waves. This study promotes the combination of so called free-field and static boundary conditions as demonstrated on a simple example. Keywords: Earthquake, shear weave, accelerogram, finite element method, free-field boundary, absorbing boundary, fixed boundary. 1. Introduction Crowded cities, pushing the need for fast and safe transportation, are the main cause of ever faster grow- ing number of underground structures. Such struc- tures, built in densely populated areas, often experi- ence a high degree of complexity. To maintain reliable performance of tunnels during their life time calls for reliable designs. These are typically based on detailed finite element simulations. In this regard, the preliminary design typically adopts a certain sim- plification such as a dimensional reduction from 3D (three-dimensional) to 2D (two-dimensional). While simplified 2D formulations attempting to account for out-of-plane deformations do exit [1–3], the majority of practical engineers rely on a fully 2D plane-strain anal- ysis employing the convergence confinement method (2DCC) either in deterministic [4] or stochastic [5] environment. Owing to its popularity and availability in various geotechnical softwares, e.g. GEO5 Tun- nel [6], the 2DCC analysis is often the choice even if addressing the TBM (tunnel boring machine) based excavation strategies. While in most applications the static loading con- ditions play the principal role in the design of tunnel lining or in assessing a potential failure during tunnel construction and when in service, the dynamic analy- sis in regions prone to the occurrence of earthquake may prove important. For a general and clear intro- duction to this subject the reader is referred to an excellent monograph by I. Towhata [7]. When limit- ing our attention to the determination of lining forces caused by earthquake it is often sufficient to carried out a simplified pseudo-static analysis with the help of either analytical [8, 9] or FEM based methods [10–12]. The principal concern in connection to a fully dy- namic analysis, regardless whether 2D or 3D, is the representation of infinite boundary in the horizon- tal direction. Assuring that horizontally propagating waves are not reflected back to the domain requires implementation of appropriate boundary conditions along vertical boundaries. These are typically in the form of viscous dampers relating the boundary trac- tions to velocities of the propagating waves. Clearly, absorbing these waves at the boundary has the same effect as if transmitting them with no restrictions. The most widely used absorbing/transmitting boundary conditions are those introduced by Lysmer and Kuhlemeyer [13]. If represented by these bound- ary conditions the boundary absorbs the entire wave propagating across this boundary. When modeling the response of structure to earthquake it appears more appropriate to absorb only the portion of en- ergy associated with the perturbation from so called free field conditions, the state preceding any construc- tion. Such perturbation can be attributed to, e.g. the built above-ground or underground structure such as embankment or tunnel. Such boundary conditions, termed henceforth the free-field boundary conditions, were put forward by Zinkiewicz, et al. [14]. To acquire the necessary free field response they build upon the solution of one-dimensional (1D) free field column problem [10, 12, 14]. Point out that free-field bound- ary conditions should become active only if there are other than vertically propagating waves. Thus to properly represent the results of 1D free field column analysis within a 2D environment calls for other types of boundary conditions that assure, in case of verti- 74 http://dx.doi.org/10.14311/APP.2018.15.0074 http://ojs.cvut.cz/ojs/index.php/app vol. 15/2018 Influence of boundary conditions in earthquake analysis cally propagating shear weaves, the stress state the same along every vertical section in the 2D model as derived from 1D analysis. Starting with a brief introduction to a theoretical background in Section 2 we continue in Section 3 with a simple parametric study on the effect of various types of potential boundary conditions applied along the lateral boundaries of the analyzed domain. The most appropriate boundary conditions are then tested on an illustrative example of ground excavation in Sec- tion 4. The essential findings are finally summarized in Section 5. 2. Theoretical background Suppose that a simple homogeneous rectangular do- main is loaded at its boundary by a vertically propa- gating shear weave as schematically shown in Fig. 1(a). The terms uu,ud stand for the upward and downward traveling waves, respectively. We further assume an infinitely stiff base rock so the total motion in the form of outcrop motion can be applied at the bottom boundary. (a) (b) Figure 1. (a) Schematic representation of loading through outcrop motion, (b) Loading and boundary conditions applied to 1D and 2D domains. The displacement field in space and time can then be decomposed as u(x, t) = u0(t) + uR(x, t), (1) where at the bottom boundary it holds u(x, t)|y=0 = u0(t), uR(x, t)|y=0 = 0. (2) The discretized form of the equation of motion then reads MüR(t) + Cu̇R(t) + KuR(t) = −Mü0(t), (3) where M, C, K are the mass, damping and stiffness matrices, respectively. Note that no material damping is assumed so that the damping matrix C collects only the viscous dampers applied along the domain lateral boundaries as shown in Fig. 1(b). Therein, the viscosity cp = √ Eoed/ρ corresponds to the pressure wave velocity whereas cs = √ G/ρ is associated with the shear wave velocity; ρ,Eoed,G stand for the soil density and oedometric and shear moduli, respectively. Further details can be found in [13–15]. The time integration of Eq. (3) is typically performed with the help of Newmark integration scheme. As already mentioned in the introductory part the analysis of 1D free field column problem seen in Fig. 1(b) is needed to provide free-field velocities u̇F FR and shear tractions py = τxy. In the present study, a simple impulse plotted in Fig. 2 was considered in both 1D and 2D analyses. (a) (b) Figure 2. Prescribed impulse: (a) acceleration, (b) displacement. To eliminate any potential errors caused by finite element discretization we assumed an optimal ele- ment size defined such that for the selected time step of 0.01s the shear wave passes through only just one element, i.e. lelem = cs∆t = √ G/ρ∆t =√ 28.6 × 106/1960 · 0.01 = 1.2 m. The potential er- ror arising from discretization if deviating from an optimal element length is clearly visible in Fig. 3. Thus only the optimal element length was adopted when performing all relevant calculations. One such analysis considers a 1D free-field column to provide results to be reproduced by a 2D analysis discussed in the next section. Example of time variation of horizontal displacements of two selected points (on surface and at the center of the domain) is plotted in Fig. 4 for illustration. 75 V. Pavelcová, T. Poklopová, T. Janda, M. Šejnoha Acta Polytechnica CTU Proceedings (a) (b) Figure 3. Influence of finite element discretization. Distribution of horizontal displacement u along verti- cal direction at time: (a) t = 0.4 s, (b) t = 4.0 s. (a) (b) Figure 4. Time variation of horizontal and vertical displacements u,v: (a) on surface, (b) at the center of the domain. As seen, the vertical displacements are clearly zero as only the horizontal displacements in terms of ac- celeration ü0(t) are prescribed to the entire domain. Owing to the fixed boundary uR(0, t) = 0 and no material damping, the wave is free to bounce back to the domain when arriving at the bottom boundary. Also note that the top boundary, free of shear stress, amplifies the propagating shear wave at the bound- ary two times in comparison to the wave within the domain, see e.g. [7]. 3. Parametric study on the influence of lateral boundary conditions To see how well may various boundary conditions rep- resent the free-field response we perform a simple para- metric study. To that end, a rectangular homogeneous domain with the shear modulus G = 28.6 MPa, the Poisson ratio ν = 0.4 and the density ρ = 1960 kg/m3 is considered. As in the previous section we reserve ourselves to fixed boundary at the bottom and loading conditions specified in Fig. 2. In particular, four types of boundary conditions plotted in Fig. 5 are examined in a sequel. (a) (b) (c) (d) Figure 5. Selected boundary conditions along lat- eral edges: (a) fixed vertical displacement v(t) = 0 (BC1), (b) free boundary-no restrictions (BC2), c) free- field boundary conditions (BC3), d) static boundary conditions - py (t) prescribed (BC4). 3.1. Fixed boundary conditions (BC1) With reference to the results provided by free field column analysis we begin by prescribing a kinematic constrain in the vertical direction, Fig. 5(a). Such boundary conditions assure zero vertical displacement along lateral edges while generating, through the reac- tion forces, the expected non-zero shear stress τxy. As evident from Fig. 6, these are the correct boundary conditions producing the same response as the 1D free field column analysis. This is supported by the distribution of the shear stress found at the end of analysis in Fig. 10(a). Unfortunately, such boundary conditions are gen- erally not acceptable as they do not allow for trans- mitting the horizontally propagating wave caused, e.g. by some structural perturbation. Thus application of such boundary conditions could only be justify in cases when the lateral boundaries are sufficiently far 76 vol. 15/2018 Influence of boundary conditions in earthquake analysis from the point of interest, e.g. tunnel, and the re- flected horizontal wave would not interfere with the incoming wave. However, such a requirement may lead to computationally demanding geometrical mod- els. Creating unstructured meshes with coarsening in the vicinity of lateral boundaries as typical of static analysis may not be the solution bearing in mind the mutual connection between the integration time step and the element size. Such a boundary condition is, therefore, generally not recommended. (a) (b) Figure 6. Time variation of horizontal and vertical displacements u,v: (a) on surface, (b) at the center of the domain. 3.2. Free boundary (BC2) Just for illustration we consider next a free boundary with no restrictions. Similar to terrain surface, such a boundary requires the shear stresses along lateral edges be zero. This in turn generates unacceptable vertical displacements within the whole domain and yield a nonuniform distribution of the shear stress as evident in Fig. 10(b). The evolution of ineligible vertical displacements is further seen in Fig. 7(b), which generates a strong deviation of the horizontal displacement in every section of the 2D domain when compared to the 1D free field column analysis. For the sake of clarity, the influence of free boundary on the results pertinent to surface point on the left lateral edge are reported in Fig. 7(a). The free boundary is therefore a totally incorrect boundary condition, even if the lateral boundaries are sufficiently far. (a) (b) Figure 7. Time variation of surface displacements: (a) horizontal displacement u, (b) vertical displacement v at the edge. 3.3. Free-field boundary conditions (BC3) Recall that free-field boundary conditions in Fig. 5(c) should allow for transmitting the structural or ma- terial perturbation only causing, as will be seen in Section 4, e.g. a horizontally propagating shear weave. However, no such perturbation is present in the present example. So both viscous dampers (the horizontal one transmitting the horizontally propagat- ing pressure wave, the vertical one transmitting the horizontally propagating shear wave) should remain inactive during the analysis. However, the formulation of damper does not dis- tinguish between the type and direction of the propa- gating wave. With reference to the previous section it is therefore expected the vertical damper to become active as attempting to damp the difference between the 1D (v = 0) and 2D (v 6= 0) analyses. This is evi- dent from Fig. 8(b) showing the variation of surface vertical displacement being 10 times smaller than in case of free boundary examined in Section 3.2. The results, fully compatible with the 1D free field column analysis, could thus be achieved by setting the shear viscosity cs to infinity, see [15]. Clearly, the infinite viscosity acts in this case as a penalty forcing the velocity of vertical displacement to zero thus in turn generating, similar to Section 3.1, the nonzero shear stresses. Nevertheless, even with damper of finite viscosity the horizontal displacements in Fig 8(a) are repre- 77 V. Pavelcová, T. Poklopová, T. Janda, M. Šejnoha Acta Polytechnica CTU Proceedings sented fairly well. Similarly, the distribution of shear stresses displayed in Fig. 10(c) is considerably im- proved. But keep in mind that the action of vertical damper is in this case purely artificial, as it is activated for purpose, for which it is not actually designed. (a) (b) Figure 8. Time variation of surface displacements: (a) horizontal displacement u, (b) vertical displacement v at the edge. 3.4. Static boundary conditions (BC4) Intuitively, and recalling the discussion in Section 3.1, the most appropriate boundary conditions to repro- duce the 1D free field column analysis are static bound- ary conditions represented by the prescribed shear tractions as shown in Fig. 5(d). It has been demon- strated in [15] that such a boundary conditions yield identical displacements and stresses as produced with the BC1 type of boundary conditions, even if com- bined with the free-field boundary conditions (BC3), which in this case remain inactive. Figure 9. Periodic boundary conditions. It is interesting to point out that similar results would be derived with the application of periodic boundary conditions [14] schematically depicted in Fig. 9. These boundary conditions, widely used in the field of homogenization [16], are essentially the tying constrains which force the homologous displacements on the opposite sides of the analyzed domain the same. [kPa] (a) [kPa] (b) [kPa] (c) Figure 10. Distribution of shear stress τxy at the end of analysis for different lateral boundary conditions: (a) BC1 and BC4, (b) BC2, (c) BC3. 4. Example to illustrate the effect of combined free-field and static boundary conditions To demonstrate the applicability of the proposed boundary conditions we consider a geometrical model plotted in Fig. 11(a) and loaded by the same impulse as used in the previous section, recall Fig. 2. In par- ticular, the effect of perturbation caused by the 90◦ wedge will be examined. In this particular example it is expected that the upward traveling shear wave, when hitting the wedge boundary, will be reflected in the 90◦ angle to become a horizontally propagating shear wave as illustrated in Figs. 11(b,c). To appreciate the action of free-field boundary con- ditions we begin with the static boundary conditions only and follow the time evolution of the vertical and horizontal components of the displacement vector linked to one particular point marked in Fig. 11(a) 78 vol. 15/2018 Influence of boundary conditions in earthquake analysis (coordinates [40,30]). The results appear in Fig. 12. It is clear from Fig. 12(b) that in view of the horizon- tally propagating shear wave the horizontal boundary behaves as a free boundary and reflects the wave back into the domain. Because the bottom boundary is assumed fixed, the wave is trapped in the domain, and in case of no material damping will continue to infinitely move between the two boundaries. The fact that the impulse is applied along the entire bottom boundary promotes interaction of individual waves which in turn is manifested by a continuous rise of the wave amplitude. (a) (b) (c) Figure 11. Testing boundary conditions: (a) geomet- rical model, (b) static boundary conditions assumed only, (c) combination of free-field and static boundary conditions. A considerably different response is observed when properly combining the static boundary condition with the free-field boundary condition, see the red curves in Figs. 12. As evident from Fig. 12(b), the latter boundary condition is able to transmit, damp, the horizontally propagating incoming wave across the lateral edges, so it is not reflected back. The subsequent vibration is again the result of interaction of the neighboring waves. The continuous vibration of the tracked point in the horizontal direction, see Fig. 12(b), is caused by continuously traveling shear wave in the vertical direction caused again by the presence of fixed boundary assumed at the bottom of the model. To bring the results closer to reality one would have to either introduce the material damping or the absorbing boundary condition at the bottom boundary, see e.g. [11, 12, 15]. 5. Conclusion Application of various boundary conditions along ver- tical edges of a rectangular domain was critically exam- (a) (b) Figure 12. Time variation of the displacement of the marked point: (a) horizontal displacement u, (b) ver- tical displacement v. ined in this paper. It was proved first in Section 3 that static boundary conditions are the correct boundary conditions to represent the 1D free field column anal- ysis in 2D simulations. These boundary conditions were then combined in Section 4 with the free-field boundary conditions to allow for properly transmit- ting the horizontally propagating shear weave across the lateral boundaries perturbed in this particular case by an excavated wedge. The ability of free-field boundary conditions to absorb such waves has been confirmed. However, the performance of these absorb- ing conditions in cases of incident angles other than 90◦ is still an open issue and will be examined in our future studies. 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] T. Janda, M. Šejnoha, J. Šejnoha. Application of converegence measurements in 2D analysis of tunnel excavation in comparison with Convergenece Confinement Method. Tunel 19(4):52–57, 2010. [2] T. Janda, M. Šejnoha, J. Šejnoha. Modeling successive excavation within two dimensional finite element mesh. Acta Geodynamica et Geomaterialia 8(1):69–78, 2011. [3] T. Janda, M. Šejnoha, J. Šejnoha. Modeling of soil structure interaction during tunnel excavation: An engineering approach. Advances in Engineering 79 V. Pavelcová, T. Poklopová, T. Janda, M. Šejnoha Acta Polytechnica CTU Proceedings Software 62-63:51–60, 2012. Special Issue dedicated to Professor Zdeněk Bittnar on the occasion of his Seventieth Birthday: Part I. [4] T. Svoboda, D. Mašín. Comparison of displacement field predicted by 2D and 3D finite element modelling of shallow NATM tunnels in clays. Geotechnik 34(2):115–126, 2011. [5] T. Janda, M. Šejnoha, J. Šejnoha. Applying bayesian approach to predict deformations during tunnel construction. International Journal for Numerical and Analytical Methods in Geomechanics 42:1765–1784, 2018. [6] Fine ltd. Geo 5 - Civil Engineering Sowftware. [online]: www.fine.cz. [7] I. Towhata. Geotechnical Earthquake Engineering. Springer-Verlag Berlin Heidelberg, 2008. [8] J. Penzien. Seismically induced racking of tunnel linings. Earthquake engineering and structural dynamics 29:683–691, 2000. [9] Y. Hashash, J. Hook, B. Schmidt, J. Yao. Seismic design and analysis of underground structures. Tunneling and undergrounding Space Technology 16:247–293, 2001. [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. [11] 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. [12] T. Poklopová, V. Pavelcová, T. Janda, M. Šejnoha. Evaluation of real underground structure subjected to earthquake - pseudostatic approach. Acta Polytechnica 2018. Accepted. [13] J. Lysmer, R. Kuhlemeyer. Finite dynamic model for infinite media. Journal of the engineering mechanics division, Proceedings of the American Society of Civil Engineers 95(EM4):859–877, 1969. [14] O. Zinkiewicz, N. Bicanic, F. Shen. Generalized Smith boundary - a transmitting boundary for dynamic computations, vol. 207. INNME Swansea, 1986. [15] 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. [16] M. Šejnoha, J. Zeman. Micromechanics in Practice. WIT Press, Southampton, Boston, 2013. 80 Acta Polytechnica CTU Proceedings 15:74–80, 2018 1 Introduction 2 Theoretical background 3 Parametric study on the influence of lateral boundary conditions 3.1 Fixed boundary conditions (BC1) 3.2 Free boundary (BC2) 3.3 Free-field boundary conditions (BC3) 3.4 Static boundary conditions (BC4) 4 Example to illustrate the effect of combined free-field and static boundary conditions 5 Conclusion Acknowledgements References