MODELING OF UNDERGROUND STRUCTURES SUBJECTED TO EARTHQUAKE BY COMBINING 1 D FREE-FIELD AND 2 D PSEUDO-STATIC ANALYSES

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.


Introduction
Design of underground structures loaded by earthquake is typically based on analytical solutions utilizing pseudostatic methods.This may yield less accurate or rather conservative results.In addition, adopting simplified methods generally calls for a number of restrictions such as the assumption of homogeneous elastic medium and a circular or rectangular frame structure representing the tunnel.If these assumptions 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 vibrating 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][3][4].In this step, the soil-structure interaction is, however, neglected.
Overcoming the above restrictions offers an application 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 analysis [7,8].Nevertheless, such an approach typically suffers from the lack of input parameters, notwithstanding the computational complexity.A suitable alternative, at least in the stage of a preliminary design, 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 practical application.
We open the discussion is Section 2 by briefly introducing 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.

Calculation steps
To introduce the subject we first present two computational 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.

Excavation step -selection of geometrical model
No doubt that the response of soil subjected to earthquake depends on the current state of stress generated Soil layer 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 classical 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. 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 internal 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 models in Fig. 1 provided essentially the same estimates of the initial distribution of internal forces thus giving seemingly the preference to less computationally demanding Model-1.

Pseudostatic analysis -loading conditions
As mentioned in the introductory part, the pseudostatic analysis considers the kinematic boundary conditions 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 propagating through a homogeneous material the free-field shear strain is provided by where c s = G/ρ is the velocity of the shear wave and v s 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 earthquake such as the one in Fig. 3.It will be shown in the next section that in such a case it is more appropriate to derive the shear strain directly from a fully dynamic analysis of a 1D free field column.

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), having the same soil profile as the actual 2D computational 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.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 truncated within a layer.Such a boundary, depicted in Fig. 4(b), requires the outgoing wave u O to be transmitted across this boundary.Because the layer, containing 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 u u , u d 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 u I , will be prescribed in case of absorbing boundary.

Theoretical background
We begin by writing the total displacement field u(y, t) = u u (y, t) where u I0 (t) = u(0, t) is the displacement field prescribed at the bottom boundary at y = 0, and u R (y, t) = {u R (y, t), v R (y, t)} is the searched relative displacement vector.Going back to Fig. 1(a) it becomes clear that in case of fixed boundary it holds u R (0, t) = 0 at the bottom boundary. ( 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 where (f ) stands for the space derivative with respect to y.As shown in [9,10] we may further write where ( ḟ ) represents the time derivative.Substituting Eqs. ( 5) and ( 6) into Eq.( 4) gives for y = 0, see also Fig. 4(b), Considering the shear wave only, Eq. ( 7) combined with the Hook law provides the static boundary condition in the form τ (0, t) = Gu R (0, t) = ρc s ( uR (0, t) − uI0 (t)) (8) to be satisfied at the bottom boundary.Note that this condition corresponds to a viscous damper with the viscosity η = ρc s introduced at y = 0, recall Fig. 4(b).
In the absence of material damping the principal of virtual work reads where M, K represent the mass and stiffness matrices and the damping matrix C, active at the bottom boundary only, is provided by where c p = E oed /ρ is the pressure wave velocity with E oed 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].

Loading conditions -shear strain
The time variation of the shear strain u (y, t) was found by employing Eq. ( 9).To that end, the acceleration üI0 and velocity uI0 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 derived assuming the elastic behavior.In this regard, the elastic moduli G, E oed corresponding to an unloading/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 conditions 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.

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 displacements compatible with the average shear strain provided by the free field analysis, see Fig. 6.The value of the maximum displacement associated with the fixed boundary immediately suggests inadequacy 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 infinity 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 the procedure outlined in [2].In our particular case, see [4] for details, the particle velocity v s 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 conditions 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 pseudostatic analysis.This becomes clear by examining the results plotted in Fig. 7.These were derived by considering 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 domain3 .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 unacceptable 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.Such conditions then yield the expected distribution of shear strains for the two computational models as evident from Fig. 8. Clearly, the perturbation caused by the presence of tunnel is confined to its close vicinity.

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 moments 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).
-  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 circumferential angle.
The use of fixed boundary conditions is clearly incorrect as the maximum moment in Fig. 9(a) exceeds 10 times the value of the moment at the end of excavation 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 a result is highly improbable.Much reasonable results 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 boundary conditions, at least for the present geometrical and material settings.

Conclusions
The present paper was concerned with the application of pseudostatic analysis to address the effect of earthquake 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 infinitely stiff base rock, the absorbing boundary must be introduced if the depth of computation model appears 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.

Figure 2 .
Figure 2. Evolution of lining moment as a function of mesh density.

Figure 4 .
Figure 4. a) Free field column with fixed boundary, b) Free field column with absorbing boundary, c) Loading and boundary conditions.

Figure 5 .
Figure 5. Distribution of maximum shear velocity and shear strain: a) fixed boundary, b) absorbing boundary.

Figure 7 .
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 strainconstant stiffness with the layer and elastic model.

Figure 9 .
Figure 9. Distribution of moments prior to and after running the pseudostatic analysis: a) fixed boundary, b) absorbing boundary.

Figure 10 .
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.

Figure 11 .
Figure 11.Distribution of moments comparing the impact of fixed and absorbing boundary conditions.

Table 1 .
Material parameters of subsoil layers.

Table 2 .
Evolution of settlement and lining internal forces as a function of mesh density.