Numerical validation of analytical solutions and their use for equivalent-linear seismic analysis of circular tunnels

Soil Dynamics and Earthquake Engineering 66 (2014) 206–219

Contents lists available at ScienceDirect

Soil Dynamics and Earthquake Engineering journal homepage: www.elsevier.com/locate/soildyn

Numerical validation of analytical solutions and their use for equivalent-linear seismic analysis of circular tunnels S. Kontoe n, V. Avgerinos, D.M. Potts Department of Civil & Environmental Engineering, Imperial College London, South Kensington Campus, London SW7 2AZ, UK

art ic l e i nf o

a b s t r a c t

Article history: Received 12 December 2013 Received in revised form 20 June 2014 Accepted 8 July 2014

The first part of this paper presents an extensive validation of four analytical solutions for the seismic design of circular tunnels. The validation is performed with a quasi-static finite element (FE) model which conforms to the assumptions of the analytical solutions. Analyses are performed for a wide range of flexibility ratios, slippage conditions at soil–lining interface, assuming both drained and undrained behaviour. Based on the numerical predictions the relative merits of the considered analytical solutions are discussed and recommendations are given for their use in design. The second part of this paper explores the use of equivalent linear soil properties in analytical solutions as an approximate way of simulating nonlinearity. The results of equivalent linear site response analyses are used as an input for the analytical solutions. The comparison of the analytical predictions with nonlinear numerical analysis results is very satisfactory. The results of this study suggest that analytical solutions can be used for preliminary design using equivalent linear properties and the corresponding compatible strain as an approximate way of accounting for nonlinear soil response. & 2014 Elsevier Ltd. All rights reserved.

Keywords: Quasi-static finite element analysis Seismic design of tunnels Equivalent linear site response analysis

1. Introduction The seismic design of circular tunnels is currently based on analytical solutions [11,19,4,5,18], deformation based quasi-static finite element (FE) analysis [15,10,22,2]and time-domain FE analysis [15,14,1,6,2,16] which simulate vertical propagation of shear waves. The above mentioned analytical solutions and quasi-static FE analysis consider the lining ground interaction, but they neglect any inertial interaction which can only be examined with a timedomain analysis. The analytical solutions are very attractive tools for preliminary design, as they provide a quick and easy calculation of the seismic design loads in the tunnel lining in terms of axial force and bending moment. Most analytical solutions assume either zero friction (fullslip condition) or full connexion (no-slip condition) between the soil and the lining. Park et al. [18] provides the only analytical solution, to the authors' knowledge, that considers the slippage effect at the soil–lining interface using a spring-type flexibility coefficient. As suggested by many studies [9,22,16] for most cases the interface condition is actually somewhere between full-slip and no-slip and, in that respect, the contribution of the Park et al. [18] solution is significant. In practice though, it is difficult to, a-priori, quantify the

n

Corresponding author. Tel.: þ 44 20 7594 5996; fax: þ 44 20 7594 5934. E-mail addresses: [email protected] (S. Kontoe), [email protected] (V. Avgerinos), [email protected] (D.M. Potts). http://dx.doi.org/10.1016/j.soildyn.2014.07.004 0267-7261/& 2014 Elsevier Ltd. All rights reserved.

degree of slippage and consequently it is difficult to select an appropriate value for the flexibility coefficient. The most widely used and discussed analytical solutions are the generalised method of Hoeg [11] (often referred to as Wang [25]) and the method of Penzien [19]. These two approaches result in identical seismic loads for the case of full-slip assumption along the soil–structure interface, whereas for the no-slip case significant discrepancies are observed [9,15,18]. Hashash et al. [10] and more recently Sedarat et al. [22] compared the results of simple quasistatic finite element analyses with the two analytical solutions in terms of axial force and bending moments for the no-slip case. These studies not only confirmed the differences between the two considered analytical solutions, but they also revealed that the Penzien [19] approach severely underestimates the seismic axial forces. Another important assumption which has been adopted by most analytical solutions is that of dry soil conditions. Naturally, this assumption is not realistic when tunnels are placed in saturated deposits and are subjected to rapid loading. Bobet [4,5] dealt with this limitation and provided a solution for undrained conditions. In the first part of this paper, an extended validation of four analytical solutions [11,19,4,5,18] is presented for various slippage conditions, a wide range of flexibility ratios and for both drained and undrained conditions. The validation is performed using linear elastic quasi-static analyses with the finite element code ICFEP [20]. The second part of this paper explores the use of equivalent linear properties in analytical solutions as an approximate way of

S. Kontoe et al. / Soil Dynamics and Earthquake Engineering 66 (2014) 206–219

simulating soil's nonlinearity, which was previously suggested by Wang [25] and later employed by Hashash et al. [9] and Amorosi and Boldini [1]. A recent study by Argyroudis and Pitilakis [2] shows a poor comparison between analytical solutions and nonlinear quasistatic FE analysis, thus suggesting the exclusive use of fragility curves derived by numerical analysis for preliminary design. This study shows that the results of most analytical solutions when adopting equivalent linear properties compare very favourably with nonlinear FE analysis and provides possible explanations for the discrepancy with the Argyroudis and Pitilakis [2] findings.

2. Problem statement and analytical solutions It is widely accepted [17] that the response of circular tunnels to ground shaking can be described by three main deformation modes: (a) axial deformation along the tunnel, (b) longitudinal bending and (c) ovaling deformation of the tunnel's cross section. The first two deformation modes occur due to seismic waves propagating parallel or at an angle to the tunnel axis, while the latter one is related to seismic shear waves propagating perpendicularly to the tunnel axis. Obviously, in order to simulate rigorously all these three modes of deformation, a three-dimensional model is needed. It is widely

207

accepted though that two-dimensional plane strain models provide a reasonable approximation of the problem, as the most critical mode is the ovaling deformation of the tunnel's cross section [19] which is schematically illustrated in Fig. 1a. In the present study four analytical solutions (generalised method of Hoeg's [11] (often referred to as Wang [25]), Penzien [19], Park et al. [18] and Bobet [5]) are examined and compared with numerical results, concerning the bending moments, axial forces and lining deflection due to ovaling deformation. When employing such analytical solutions it is very important to first appreciate their underlying assumptions:

 The tunnel lining (for the seismic loading) is considered to act as  

an elastic beam under a uniform shear-strain field of amplitude γmax, ignoring any inertial soil–lining interaction. The initial stresses in the ground and the lining are not considered and the excavation is assumed to take place simultaneously with the construction of the lining. The strain is applied under fully drained conditions (note that Bobet [5] gives a separate solution for undrained conditions).

It should be mentioned that the dynamic interaction can be of importance when (a) the dimensions of the tunnel's cross section

Fig. 1. (a) Ovaling deformation due to vertically propagating shear waves (adapted from [17]), (b) corresponding seismic shear loading and (c) equivalent static loading (adapted from [18]).

208

S. Kontoe et al. / Soil Dynamics and Earthquake Engineering 66 (2014) 206–219

are comparable to the wavelengths of the predominant frequency of the earthquake loading, (b) the tunnel is relatively shallow and (c) for stiff structures in soft soil. The loading in the tunnel under static conditions depends on the stiffness of the lining, the overburden pressure and the coefficient of earth pressure at rest, K0. For the seismic case, the free-field shear stress replaces the in-situ overburden pressure and the at-rest coefficient of earth pressure is assigned a value of (  1). Hence the loading induced by shear waves is assumed as being equivalent to a compressive and tensile principal stress acting at 451 with respect to the horizontal axis (as shown in Fig. 1b and c). Clearly due to this assumption, regarding the seismically induced stress state in the field, all these analytical solutions predict maximum ovaling deformation at the major and minor axes at θ¼ 7 451 from the horizontal axis. The main parameters that govern the predicted loads by the analytical solutions are the compressibility, C, and the flexibility, F, ratios and the shear-strain γmax. The shear strain is calculated either as the ratio of the particle velocity to the shear wave propagation velocity or, preferably, from a site response analysis as the maximum shear strain at the level of the tunnel axis. The compressibility, C, and the flexibility, F, ratios are means to quantify the relative stiffness of the medium compared to that of the lining. The compressibility ratio is a measure of the extensional stiffness of the medium relative to that of the lining and the flexibility ratio is a measure of the flexural stiffness C¼

Em ð1  ν2l Þr El tð1 þ νm Þð1  2νm Þ

ð1Þ

F¼

Em ð1 ν2l Þr 3 6El Ið1 þ νm Þ

ð2Þ

where Em and El are Young's moduli of the medium and the lining respectively, νm and νl are Poisson's ratios of the medium and the lining respectively, r and t are the radius and the thickness of the lining, respectively and I is the moment of inertia per unit width of the lining. Values of the flexibility ratio F greater than 1 suggest that the tunnel lining has lower stiffness than the surrounding medium. Hence, values of F-1 imply that the tunnel undergoes identical deformations to that of an unlined tunnel without resisting the ovaling deformation. The analytical solutions examine separately the cases of full and no-slip condition along the interface of ground and lining. Slip at the interface is generally possible depending on the construction method, the relative stiffness at the soil–lining interface and the severity of the seismic loading, while as many studies suggest [9,22,16] in most cases the reality is somewhere in between full-slip and no-slip. For completeness the four analytical solutions are presented concisely in the Appendix.

3. Numerical validation of analytical solutions 3.1. Numerical model and parameters considered The validation of the analytical solutions was based on a simple quasi-static model which was created with the finite element code ICFEP [20]. Quasi-static finite element analysis can either be forced based, in which seismically induced inertia forces are approximated as a constant body force, or deformation based, in which the mesh is subjected to shear deformation, as schematically illustrated in Fig. 2. In this study a deformation based quasistatic methodology is adopted, as this approach is consistent with the assumptions of the analytical solutions and it does not suffer from the inherent limitations of the forced based approach (e.g. sensitivity on the dimensions of the model) [13].

Fig. 2. Schematic representation of the mesh configuration in a quasi-static deformation based analysis consistent with simple-shear conditions. Table 1 Soil parameters and corresponding flexibility ratios. Soil's Young's modulus Εm (kPa)

Flexibility ratio F

2018.2 7427.1 16,145.8 35,843.7 350,455.7 1,120,000.0 2,240,000.0 3,000,000.0

0.12 0.46 1.00 2.22 21.7 69.4 138.7 185.8

The assumed model represents a soil layer of thickness of H¼ 50 m overlying bedrock. The vertical boundaries of the mesh are 50 m away from the centre of the tunnel. The soil was considered to behave elastically with Poisson's ratio νm ¼0.25, while its Young's modulus, Em , was varied parametrically in order to capture a wide range of flexibility ratios F (Table 1). The modelled tunnel has a circular cross section of radius r¼ 3 m and its centre is located 15 m below the ground surface. The lining's width is t¼0.3 m and was modelled with elastic beam elements of Young's modulus El ¼24.8 GPa and Poisson's ratio, νl ¼0.2. Prior to all quasi-static analyses, a static analysis was undertaken to establish the initial stresses acting on the lining. During the static analysis, displacements were restricted in both directions along the bottom mesh boundary and horizontal displacements were restricted along the side boundaries. The tunnel construction was modelled using the convergence-confinement method which is described in detail by Potts and Zdravković [21]. The tunnel excavation was performed in 10 increments and the linings were constructed at 50% of stress relaxation (i.e. increment 5) prior to the completion of excavation. In order to be consistent with the assumptions of the analytical solutions (which ignore the initial stresses in the ground and the lining), in all cases the “seismic increment” was obtained by subtracting the loads computed at the end of the excavation from those at the end of the quasi-static analysis. For the modelling of the slippage, interface elements of zero thickness were used [8] with a very high normal stiffness (KN ¼ 107 kN/m3) and varying shear stiffness depending on the degree of slippage. In particular for the full-slip case a nominal value of shear stiffness (Ks ¼1 kN/m3) was used. In order to be consistent with the assumptions of the analytical solutions, the vertical and horizontal displacements were restricted along the bottom boundary, while the vertical displacement was also restricted along the side and top boundaries of the model. Furthermore, a uniform displacement u was applied along the top boundary, while a triangular displacement distribution was applied along the vertical boundaries of the mesh (see Fig. 2). The magnitude of the uniform applied displacement u was u ¼0.126 m corresponding to a shear strain: γ max ¼

u 0:126 m ¼ ¼ 0:00252 H 50 m

ð3Þ

S. Kontoe et al. / Soil Dynamics and Earthquake Engineering 66 (2014) 206–219

209

Fig. 3. Comparison of calculated thrust (a) and bending moments (b) by analytical solutions and numerical analysis for the full-slip condition.

This value of shear strain is identical to one adopted by Hashash et al. [9] in their validation study. 3.2. Drained analysis In this section, the results of a series of numerical analyses for both full-slip and no-slip conditions are compared with the four analytical solutions. Fig. 3 presents the comparison for the full-slip assumption in terms of lining thrust and bending moment. The resulting thrust and bending moment for all analytical methods coincide, while the analytical results are also in excellent agreement with the numerical results for all examined flexibility ratios. It should be noted that Sedarat et al. [22] performed a similar comparison for three flexibility ratios and two lining thickness (t ¼0.36 m and t¼1.26 m). They concluded that the results of extended Hoeg and Penzien for the full-slip condition in terms of thrust loads are comparable with the numerical results only for the thick lining and for low values of the flexibility ratio (up to 3.2). In the present study although the lining thickness (t¼ 0.3 m) is similar to the thin lining case of Sedarat et al. [22] (which showed poor agreement with the analytical solutions) the computed thrust values herein are in excellent agreement with the analytical results for all examined flexibility ratios. It should be noted that Sedarat et al. [22] modelled the full-slip case with a frictionless contact, adopting an elastoplastic law for the soil– lining interface. This prohibits the development of tensile normal contact stresses and leads to the observed discrepancies, for the full-slip case, with the numerical results of the present study and the analytical solutions. While the modelling of the interface condition in the present study is consistent with the assumptions of analytical solutions, the use of an elastoplastic law for the interface condition is more realistic providing that it is accompanied by an elastoplastic model for the surrounding soil. Among the existing analytical solutions, Park et al. [18] is the only one that considers the slippage effect at the soil–lining interface employing a spring-type flexibility coefficient (D). According to Park et al. [18], while the flexibility coefficient can theoretically vary from D ¼0 (corresponding to the no-slip case) to D-1 (representing the full-slip case), in practice it is limited to a narrower range. Fig. 4 compares the calculated thrust and bending moments with the Park et al. [18] solution for D ¼10  6–10  4

which corresponds to Ks values of 106–104 kN/m3 for the shear stiffness of the interface elements. The numerical results compare very well with the analytical solution both in terms of thrust and bending moments for all the considered flexibility coefficients. The flexibility coefficient D has a more pronounced effect on the resulting thrust values rather than the bending moments with the effect becoming progressively larger as the flexibility ratio (F) increases. Figs. 5 and 6 present the comparison for the no-slip case again in terms of lining thrust and bending moment. In Fig. 5 the results for the lower examined flexibility ratios are shown (up to F ¼2.22), while Fig. 6 presents the results for the higher flexibility ratios. Concerning the thrusts, Hoeg, Park et al. and Bobet solutions match perfectly the numerical solution for the lower flexibility ratios (Fig. 5a), while the comparison is also very good for the larger ones (Fig. 6a). However Penzien's solution gives comparable results with the other approaches only for very low flexibility ratios, while as the flexibility ratio increases the Penzien method results diverge from the numerical solution considerably, reaching differences of 28.3% and 46.3% for flexibility ratios F ¼1.0 and F ¼2.22 respectively. The underestimation of the thrust values by this approach becomes even worse for greater flexibility ratios, as the predicted thrust values are lower by at least an order of magnitude compared to the other approaches and the numerical solution. This is in agreement with the findings of Hashash et al. [10] and Sedarat et al. [22]. With regards to the bending moments, the comparison between Hoeg's or Penzien's approaches and the numerical solution is very good only for very low flexibility ratios F o1.0, while the results of Park et al.'s and Bobet's solutions are matching the numerical ones for even larger flexibility ratios (up to F ¼2.22 in Fig. 5b). For flexibility ratios F 421.7 (Fig. 6b) the numerical solution predicts a very small change in bending moments, while the analytical solutions result in an increase of its value with increasing flexibility ratio (with a reduced rate though). The percentage differences of the bending moments between analytical and numerical solutions are summarised in Table 2. The superiority of the Park et al. [18] and Bobet [5] solutions is clear as the maximum observed deviation, for the flexibility ratios examined herein, is almost half of the other two approaches (15% compared to 27% deviation of the other two solutions).

210

S. Kontoe et al. / Soil Dynamics and Earthquake Engineering 66 (2014) 206–219

Fig. 4. Comparison of calculated thrust (a) and bending moments (b) by Park et al. [18] and numerical analysis for various slippage levels.

Fig. 5. Comparison of calculated thrust (a) and bending Moments (b) by analytical solutions and numerical analysis for the no-slip condition and for F¼ 0.125, 0.46, 1.0 and 2.22.

In order to compare the thrusts and moments (for the no-slip assumption) of the analytical procedures with the numerical solution for a greater range of flexibility ratios, all the above calculations and analyses were repeated for an increased tunnel radius of r ¼5 m. All the remaining parameters concerning both the soil and the lining were kept the same. The larger radius led to an increase of all the examined flexibility ratios. The results are summarised in Fig. 7, while a logarithmic axis is used because the range of the investigated flexibility ratios varied from less than 1 to almost 650. The results, both for thrusts (Fig. 7a) and for bending moments (Fig. 7b), confirm the observations of Figs. 5 and 6 and show that the discrepancy between analytical and numerical solutions in terms of bending moment increases with increasing flexibility ratio.

Furthermore Fig. 8 compares the results of the analytical and numerical solutions in terms of normalised lining deflection at θ¼ 7451 from the horizontal axis, which is defined as follows: Δdl Δdl ¼7 Δdfreefield rγ max

ð4Þ

where Δdl is the lining deflection, Δdf reef ield is the corresponding diametric strain for a circular cross section for non-perforated ground, r is the tunnel radius and γmax is the maximum free-field shear-strain at the level of the tunnel axis. In both cases (i.e. fullslip and no-slip) the discrepancy between the numerical predictions and the analytical solutions is significant for flexibility ratios less or equal to one, while the comparison improves for higher values of flexibility ratio. For the full-slip case Hoeg's, Penzien's

S. Kontoe et al. / Soil Dynamics and Earthquake Engineering 66 (2014) 206–219

211

Fig. 6. Comparison of calculated thrust (a) and bending moments (b) by analytical solutions and numerical analysis for the no-slip condition and for F¼ 21.7, 69.4, 138.7 and 185.8.

Table 2 Percentage differences (%) between calculated maximum bending moments by analytical solutions and numerical analysis for the no-slip assumption. F

Extended Hoeg

Penzien [19]

Park et al [18]

Bobet [5]

0.12 0.46 1.00 2.22 21.7 69.4 138.7 185.7

11.93 13.63 15.64 18.04 22.89 25.12 26.85 27.63

0.19 3.86 7.97 12.88 22.06 24.86 26.72 27.53

 1.28  0.62 0.34 1.51 5.18 9.34 13.35 15.33

 1.28  0.62 0.34 1.51 5.18 9.34 13.35 15.33

and Park et al.'s solutions give identical results, while the Bobet [5] solution gives lower values for the entire range of flexibility ratios. For the no-slip case and for F 41 the solutions of Park et al. [18] and Bobet [5] coincide and are in better agreement with the numerical predictions than the remaining two solutions. It should be noted that Hashash et al. [10] observed a good agreement between Hoeg's and Penzien's approaches and the numerical solution in terms of normalised lining deflection examining the no-slip case (percentage differences up to 15.5%), but they considered a limited range of flexibility ratios (F ¼9.63, 16.2 and 18.58) within which the present study also shows a satisfactory comparison (12% for Hoeg and 11% for Penzien for F ¼21.71). 3.3. Undrained analysis Bobet [4] extended the analytical solutions for dry ground to fully saturated ground under undrained conditions, for both full slip and no-slip conditions and expressed the bending moment and the thrust as a function of the quasi-static free-field vertical stress ðσ v Þ and horizontal stress ðσ h Þ. Later, Bobet [5] presented the equations as a function of the far-field shear stress (as shown in the Appendix) and further extended the methodology for the case of rectangular tunnels. In order to validate the expressions proposed by Bobet [5] for undrained conditions, the set of analyses for a tunnel of radius r ¼3.0m was repeated with the water table at the ground surface and undrained conditions. Fig. 9 presents the comparison for the

full-slip case, showing a good agreement between the analytical and the numerical solution both in terms of thrust and bending moments with the corresponding percentage differences reaching almost 6% for the thrust and 10% for the bending moment for the high flexibility ratios (Table 3). The deviation between the numerical and analytical solution is more significant for the no-slip case as depicted in Fig. 10 and Table 3. The numerical analysis predicts larger values of thrust and bending moment than Bobet [5] for low flexibility ratios, while this trend is reversed for F421.7. While the comparison in terms of bending moments is within acceptable level the differences in terms of thrust values are significant. It should be noted that the expressions of Bobet [5] for dry ground and saturated ground under undrained conditions give identical results if a Poisson ratio, νm ¼ 0:5 is adopted.

4. Numerical evaluation of analytical solutions for nonlinear response The use of the analytical solutions in seismic design usually involves the adoption of the maximum shear modulus and the maximum free-field shear-strain as input parameters. This can lead to erroneous estimation of the lining loads, as it is well established that soil behaves nonlinearly during seismic shaking and can experience significant stiffness degradation. To improve this shortfall of the analytical solutions, it has been suggested [25] to perform first an equivalent-linear site response analysis of the stratigraphy and to then use as input parameters for the analytical solutions the equivalent linear (converged value) of shear modulus and the corresponding “effective” shear-strain (which is a fraction of γ max ) at the level of the tunnel axis. In order to verify the validity of the suggested approach, the quasi-static finite element analyses were repeated adopting a nonlinear model using for the input deformation the effective shear-strain of each equivalent analysis at the depth of the tunnel axis. 4.1. Equivalent linear analysis and nonlinear quasi-static analysis To estimate the equivalent linear shear stiffness and the compatible free-field shear-strain (γ ef f ) at the level of the tunnel

212

S. Kontoe et al. / Soil Dynamics and Earthquake Engineering 66 (2014) 206–219

Fig. 7. Comparison of calculated thrust (a) and bending moments (b) by analytical solutions and numerical analysis for the no-slip condition and for a wide range of flexibility ratios (for tunnel with radius r¼ 5.0 m).

Fig. 8. Comparison of normalised lining deflection calculated by analytical solutions and numerical analysis for a wide range of flexibility ratios. (a) Full-slip, and (b) No-slip.

axis, site response analyses of a uniform one dimensional soil column of 50 m depth were performed with the software EERA [3]. The column was subjected to the north–south component of the Griffith observatory motion recorded during the 1994 Northridge earthquake. This record was selected because of its very rich frequency content which is expected to engage a wide range of modes of the problem. A parametric analysis was undertaken varying the intensity of the ground motion and the soil type, aiming to cover a wide range of shear strains. In particular, the input motion was scaled to maximum peak ground accelerations of 0.2 g and 0.4 g and two soil types, a clean sand with plasticity index PI¼0 and a medium plasticity clay with PI¼50, were examined. The adopted stiffness degradation and damping curves were based on the expressions of Darendeli [7] and were also subsequently used to calibrate the nonlinear model of the FE quasi-static analysis. Table 4 summarises the examined cases

along with the results of the site response analysis in terms of equivalent shear stiffness and compatible free-field shear-strain (γ ef f ) at the level of the tunnel axis. It should be noted that the ratio of γ ef f to γ max was taken equal to 0.6 for all site response analyses, which according to the empirical expression of Idriss and Sun [12] implies an earthquake of moment magnitude M ¼7.0. Subsequently, for all cases shown in Table 4, nonlinear elastic quasi-static finite element analyses were carried out employing the boundary conditions and the lining parameters described previously for the linear case and assuming no-slip between the soil and the lining. These nonlinear elastic analyses were carried out with a cyclic nonlinear model for the soil which adopts a hyperbolic function for the backbone curve as described in [23, 24]. The model follows the Masing rules for unloading and reloading, but this feature of the model was not used herein as only monotonic loading was imposed for the purposes of this

S. Kontoe et al. / Soil Dynamics and Earthquake Engineering 66 (2014) 206–219

213

Fig. 9. Comparison of calculated thrust (a) and bending moments (b) by Bobet [5] and numerical analysis for undrained conditions and full-slip assumption.

Table 3 Percentage differences (%) between calculated maximum thrust and bending moments by Bobet [5] and numerical analysis for undrained conditions. F

0.12 0.46 1.00 2.22 21.7 69.4 138.7 185.7

Full-Slip

No-Slip

Tmax

Mmax

Tmax

Mmax

 3.27  0.65 1.30 3.05 5.59 5.88 5.70 5.59

 0.15 2.00 3.79 5.55 8.84 9.93 9.70 9.68

 360.74  149.66  63.89  30.57  6.75 0.85 5.56 7.67

 14.78  9.24  5.45  2.38 1.77 7.43 13.49 16.35

study. The model was calibrated to the stiffness degradation curves that were used for the site response starting with the elastic shear stiffness (Gmax) values shown in Table 4. For each analysis the imposed deformation was calculated based on Eq. (3) using the compatible free-field shear-strain (γ ef f ) at the level of the tunnel axis shown in Table 4. In order to be consistent with the assumptions of the analytical solutions, the simulation of the excavation and the tunnel construction was undertaken assuming linear elastic behaviour (using Gmax) and the nonlinear model was activated only for the quasi-static part of the analysis. In all cases the “seismic increment” was obtained by subtracting the loads computed at the end of the excavation from those at the end of the quasi-static analysis. Fig. 11 shows the thrust and bending moment variation with flexibility ratio resulting from the nonlinear FE analyses for all cases. It should be noted that the flexibility ratios have been calculated using the equivalent linear values of shear modulus shown in Table 4. The thrust values seem to reach, with some fluctuations, a plateau for high flexibility ratios, while for the linear case they were monotonically increasing with increasing flexibility ratio (see Fig. 7a). On the other hand, the bending moments reduce with increasing flexibility ratio, while in the linear case they were reaching a plateau for the high flexibility ratio (see Fig. 7b). It should be noted though that, while all linear elastic analyses were performed for the same value of shear-strain, in the nonlinear case the imposed shear-strain value varies with

flexibility ratio. The higher the flexibility ratio, the stiffer the soil and hence the lower the value of effective shear-strain that results from the site response analysis (see Table 4). This explains the differences observed in the trends of linear and nonlinear analyses, while it indicates that using the elastic properties can result in overestimation of the lining loads for high flexibility ratios. 4.2. Comparison with analytical solutions Figs. 12 and 13 compare the nonlinear analyses results with the analytical solutions for PGA values of 0.2 g and 0.4 g respectively. As previously explained Geq and the compatible free-field shearstrain (γ ef f ) at the level of the tunnel axis were used as input parameters for the analytical solutions. Similar to the linear case (see Figs. 5 and 6) Park et al. [18] and Bobet [5] solutions compare very well with the nonlinear FE analyses, Hoeg's solution shows larger differences in terms of bending moments, while Penzien [19] underestimates the thrust values. The corresponding percentage differences between analytical solutions and numerical analyses are shown in Tables 5 and 6 for PGA values equal to 0.2 g and 0.4 g respectively. Overall the comparison between the analytical solutions and numerical analysis for the thrust values deteriorates with increasing nonlinearity (i.e. higher PGA level and lower flexibility ratios), but it is always within acceptable levels (excluding Penzien for thrust predictions). Furthermore Fig. 14 shows the comparison between numerical analysis and analytical solutions in terms of normalised lining deflection. Again the agreement between the two methodologies (i.e. analytical and numerical) is good with the Hoeg and Penzien solutions showing slightly better comparison with the numerically obtained deflections. These observations confirm that analytical solutions with equivalent linear shear stiffness and the corresponding compatible shearstrain at the tunnel axis level as input parameters can be successfully used for the preliminary seismic design of circular tunnels. The proposed approach can lead to more accurate estimation of the lining loads as soil nonlinearity is approximated in a consistent manner. Recently, Argyroudis and Pitilakis [2] performed a similar comparison between analytical solutions and FE quasi-static analysis observing a poor agreement. Therefore they concluded that numerically derived fragility curves should be used for preliminary

214

S. Kontoe et al. / Soil Dynamics and Earthquake Engineering 66 (2014) 206–219

Fig. 10. Comparison of calculated thrust (a) and bending moments (b) by Bobet [5] and numerical analysis for undrained conditions and no-slip assumption.

Table 4 Parameters considered in site response analyses and corresponding results in terms of equivalent shear modulus and effective shear strain at the level of the tunnel axis. PGA (g)

Gmax (kPa)

Geq (kPa)

F

Geq/Gmax

γ ef f (%)

0

0.2

50

0.2

0

0.4

50

0.4

140,182.3 448,000.0 896,000.0 1,200,000.0 14,337.48 140,182.3 448,000.0 896,000.0 1,200,000.0 140,182.3 448,000.0 896,000.0 1,200,000.0 140,182.3 448,000 896,000 1,200,000

72,811.27 273,580.5 657,803.7 1,008,380.0 8364.264 95,935.35 336,672.5 812,796.8 1,075,343.0 38,412.16 239,162.0 559,027.3 860,149.2 79,797.1 289,267.6 748,747.8 1,017,085.0

11.27 42.36 101.85 156.14 1.30 14.85 52.13 125.85 166.50 5.95 37.03 86.56 133.18 12.36 44.79 115.94 157.48

0.52 0.61 0.73 0.84 0.58 0.68 0.75 0.91 0.90 0.27 0.53 0.62 0.72 0.57 0.65 0.84 0.85

0.0465 0.0310 0.0167 0.0083 0.0850 0.0527 0.0367 0.0102 0.0117 0.1465 0.0437 0.0291 0.0183 0.0905 0.0637 0.0207 0.0182

PI

Fig. 11. Thrust and bending moment variation with flexibility ratio calculated with nonlinear quasi-static FE analysis.

S. Kontoe et al. / Soil Dynamics and Earthquake Engineering 66 (2014) 206–219

215

Fig. 12. Comparison of calculated thrust (a) and bending moments (b) by analytical solutions and nonlinear numerical analysis for PGA ¼0.2 g (no-slip condition).

Fig. 13. Comparison of calculated thrust (a) and bending moments (b) by analytical solutions and nonlinear numerical analysis for PGA ¼0.4 g (no-slip condition).

Table 5 Percentage differences (%) between calculated maximum thrust and bending moments by analytical solutions and nonlinear numerical analysis for PGA ¼ 0.2 g. PI

F

Extended Hoeg Penzien

Park et al.

Bobet

Tmax

Tmax

Table 6 Percentage differences (%) between calculated maximum thrust and bending moments by analytical solutions and nonlinear numerical analysis for PGA ¼ 0.4 g. PI

Tmax 0

11.27 15.47 42.36 13.11 101.85 7.39 156.14 4.04 50 1.30 5.06 14.85 9.43 52.13 7.66 125.85 2.87 166.50 3.40

Mmax

Tmax

Mmax

Mmax

17.01 22.31 27.27 28.51 14.81 21.13 25.86 28.30 28.74

 352.15  1336.53  3115.64  4548.19  40.40  507.30  1727.19  3890.11  4808.00

15.42 15.70  2.01 15.70  2.01 21.87 13.16 5.04 13.16 5.04 27.09 7.39 12.88 7.39 12.88 28.39 4.03 15.71 4.03 15.71 7.82 6.24  1.19 6.24  1.19 19.94 9.62 2.98 9.62 2.98 25.52 7.70 9.69 7.70 9.69 28.16 2.86 14.74 2.86 14.74 28.63 3.38 16.22 3.38 16.22

F

Extended Hoeg Penzien

Park et al.

Bobet Tmax

Mmax

0

5.95 37.03 86.56 133.18 50 12.36 44.79 115.94 157.48

Tmax

Mmax

Tmax

Mmax

Tmax

22.69 16.99 12.13 7.78 14.48 12.41 4.53 3.78

10.88 18.63 23.12 27.78 18.46 23.57 27.98 28.54

 154.93  1120.89  2570.96  3858.79  392.64  1419.93  3576.62  4590.82

7.99 18.11 22.90 27.65 17.01 23.16 27.83 28.43

23.07  8.93 23.07  8.93 17.04 0.37 17.04 0.37 12.14 7.45 12.14 7.45 7.77 14.32 7.77 14.32 14.70  0.27 14.70  0.27 12.46 6.67 12.46 6.67 4.53 14.11 4.53 14.11 3.77 15.79 3.77 15.79

Mmax

Mmax

216

S. Kontoe et al. / Soil Dynamics and Earthquake Engineering 66 (2014) 206–219

Fig. 14. Comparison of normalised lining deflection calculated by analytical solutions and numerical analysis for a wide range of flexibility ratios (no-slip condition). (a) PGA ¼ 0.2g, (b) PGA ¼0.4g.

design instead of analytical solutions. They attributed the observed discrepancies mainly to the following reasons:

 The use of an elastoplastic constitutive model for their FE anal 

ysis which can lead to redistribution of stresses around the tunnel lining and plastic deformations. The simulation of the excavation which essentially modifies the initial stress conditions around the tunnel prior to the quasistatic analysis. The application of a deformation variation, with depth while an average value of strain at the tunnel axis is considered for the analytical solutions.

In order to compare numerical analysis with the analytical solutions, within the context of a validation exercise, is important to model as closely as possible the assumptions of the analytical solutions. Once the validation exercise is completed, a further comparison can be carried out with more sophisticated numerical modelling to assess the validity of the assumptions of the analytical solutions. Indeed, the use of an elastoplastic model is not compatible with the assumptions of the analytical solutions, while the nonlinear modelling of the excavation process results to a stress-field around the tunnel which differs from the assumption of the analytical solutions of a uniform strain field around the tunnel. This was further investigated herein by repeating the FE analysis simulating nonlinearly the excavation process and the tunnel construction (in all previous nonlinear quasi-static analyses the excavation/construction was simulated linearly). Fig. 15 shows the thrust and bending moment variations around the tunnel lining for PI¼ 0, PGA ¼0.2 g and Gmax ¼448.0 MPa predicted by numerical analysis both for the case of linear and nonlinear excavation/construction and the corresponding no-slip Bobet [5] prediction. It should be clarified that the loads generated during excavation/construction have been subtracted to facilitate the comparison with the analytical solution. Clearly the nonlinear simulation for construction/excavation leads to a non-symmetric distribution of the loads around the tunnel lining, while the predicted response diverges considerably from the analytical solution in [5]. For both numerical analyses, the maximum response is predicted at θ ¼ 45o to the horizontal axis.

A further important reason for the discrepancy between the numerical results and the analytical solutions in [2] is the adopted boundary conditions in their quasi-static model which are not in agreement with the racking boundary conditions which are assumed by the analytical solutions. In particular, vertical displacement was not restricted along the top boundary of their model (as shown in Fig. 7a of [2]). Furthermore the variation of strain with depth which was obtained by the site response analysis was imposed to the quasi-static FE model instead of a uniform strain corresponding to the level of the tunnel axis. A direct consequence of their adopted boundary conditions is that the numerical model did not predict the maximum lining response at θ ¼ 7 45o to the horizontal axis (shown in Fig. 7b of [2]).

5. Conclusions The first part of this paper examined parametrically the reliability of four analytical solutions for the evaluation of seismically induced loading and deflection in circular lined tunnels under plane strain quasi-static conditions. This was done for a wide range of flexibility ratios and a variation of slippage conditions. Based on this parametric investigation the following conclusions can be drawn:

 Assuming full-slip conditions between the soil and the lining,





the results of the analytical methods, in terms of thrust loads and bending moments, are in perfect agreement with the numerical predictions. The numerical results also compared very well with the Park et al. [18] solution when the slippage effect between the lining and the soil was considered. In agreement with previous studies [10,15,22], the Penzien [19] solution for the no-slip condition underestimates the thrust values and thus its use should be avoided. The thrust loads predicted by the extended Hoeg [11], the Park et al. [18] and Bobet [5] approaches are in perfect agreement with the numerical solution. The bending moments resulting from the Park et al. [18] and Bobet [5] methods for the no-slip case are also in very good agreement with the numerical results. The bending moments

S. Kontoe et al. / Soil Dynamics and Earthquake Engineering 66 (2014) 206–219

217

Fig. 15. Variation of seismic thrust (a) and seismic bending moment around the tunnel lining (PI ¼ 0, PGA ¼ 0.2 g, Gmax ¼448.0 MPa).





predicted by Hoeg's method start diverging from the numerical solution even for very low flexibility ratios (i.e. for Fo1), while for the Park et al. [18] and Bobet [5] methods this divergence is minor and more pronounced only for higher flexibility ratios (F4 21.7, reaching 15.3% for F¼185.7). Generally, the deviations between analytical and numerical solutions in terms of bending moments are increasing with increasing flexibility ratio showing a trend of stabilisation for larger flexibility ratios. In terms of normalised lining deflection, the discrepancy between the numerical predictions and the analytical solutions is significant for flexibility ratios less or equal to one, while the comparison improves for higher values of flexibility ratio. For the full-slip case Hoeg's, Penzien's and Park et al.'s solutions give identical results, while Bobet's solution gives lower values for the entire range of flexibility ratios. For the no-slip case and for large flexibility ratios the solutions in [5,18] are in better agreement with the numerical predictions than the remaining two solutions. The solution in [5] for undrained conditions was found to be in perfect agreement with the numerical predictions for the fullslip case. For the no-slip case the numerical analysis predicts larger values of thrust and bending moment than in [5] for low flexibility ratios, while this trend is reversed for F4 21.7.

The second part of this paper investigated the use of equivalent linear properties in analytical solutions as an approximate way of simulating nonlinearity which was previously suggested by Wang [25] and later employed by Hashash et al. [9] and Amorosi and Boldini [1]. Equivalent linear site response analyses were undertaken varying the intensity of the ground motion and the soil type covering a wide range of shear strains. The resulting equivalent stiffness (Geq) and effective shear strain at the level of the tunnel axis were used as input parameters for the analytical results. The analytical solutions were validated for the no-slip case against FE quasi-static nonlinear analyses subjected to the same value of shear-strain, but starting the loading with the elastic properties (Gmax).

The conclusions drawn from the comparison between analytical solutions and numerical analysis is similar to the linear case. Park et al. [18] and Bobet [5] predictions, in terms of thrust and bending moments, compare well with the numerical analysis results. Overall the comparison for the thrust values deteriorates with increasing nonlinearity (i.e. higher PGA level and lower flexibility ratios), but it is always within acceptable levels. With regards to normalised tunnel deflection, the agreement between the analytical solutions and the numerical approach is also satisfactory with the Hoeg and Penzien solutions showing slightly better comparison with the numerically obtained deflections. The findings of this study confirm that analytical solutions with equivalent linear shear stiffness and the corresponding compatible shear-strain at the tunnel axis level as input parameters can be successfully used for the preliminary seismic design of circular tunnels. The proposed approach can lead to more accurate estimation of the lining loads as soil nonlinearity is approximated in a consistent manner. Analytical solutions are based, of course, on numerous assumptions and simplifications and therefore detailed design should involve a more rigorous approach. For the detailed design, quasi-static deformation based FE nonlinear analysis can be employed as this approach is capable of dealing with some of the limitations of the analytical solutions (e.g. realistic modelling of the excavation and tunnel construction). Ideally though, timedomain finite element analysis should be used employing an appropriate constitutive model which can capture fundamental facets of soil behaviour under seismic loading.

Appendix This Appendix summarises the equations employed by the four analytical methods for the calculation of axial force (thrust), bending moment and normalised lining deflection for the full and no-slip assumption at the interface. The equations are presented using a consistent notation for all methods, where Em and El are Young's moduli of the medium and the lining respectively, νm and νl are Poisson's ratios of the medium and the lining

218

S. Kontoe et al. / Soil Dynamics and Earthquake Engineering 66 (2014) 206–219

respectively, r and t are the radius and the thickness of the lining, respectively, θ is the angular location of the tunnel lining, I is the moment of inertia per unit width of the lining, γmax is the maximum free-field shear-strain at the level of the tunnels axis and C and F are the compressibility and flexibility ratios defined by Eqs. (1) and (2) respectively.

 υm Þ 3El Ið3  4υm Þ n where R ¼ 7 4ð1 ðα þ 1Þ and α ¼ r 3 G ð1  ν2 Þ m

Δdl ¼R rγ max

l

ð14Þ

Park et al. method [18] Generalised Hoeg's method (often referred to Wang [25])

(a) Full-slip  1 Em π r 2 γ max cos 2 θ þ M ¼  K1 6 ð1 þ νm Þ 4

ð5Þ

 1 Em π rγ max cos 2 θ þ T ¼  K1 6 ð1 þ νm Þ 4

ð6Þ

 νm Þ where K 1 ¼ ð2F12ð1 þ 5  6νm Þ

Δdl

1 ¼ 7 K 1 Fγ max 3 Δdf reef ield

ð7Þ

   ð1  νm ÞEm γ max r 2 4Em π cos 2 θ þ ð1  2νm ÞC þ2 þ D 4 rð1 þ νm Þ ð1 þν ÞΔ'' m

where ð8Þ

Δ'' ¼ CFð1  2νm Þ þFð3  2νm Þ ð2F þ 5  6νm ÞEm þCð2:5  8νm þ 6ν2m Þ þ 6  8νm þ 2D rð1 þ νm Þ   Δdl 4ð1  νm Þ 8Gs FD ¼ ð1  2νm ÞFC þ 2F þ r Δdf reef ield Δ''

where K2 ¼ 1þ

M¼ 

ð16Þ

(b) No-slip Em π rγ T ¼ K 2 cos 2ðθ þ Þ 4 2ð1 þ νm Þ max

Park et al. [18] introduced the shear flexibility coefficient, D in order to improve the modelling of the slip between the lining and the soil. The equations that they propose for the calculation of the seismic induced thrusts and bending moments are given in relation to this shear coefficient, with D ¼0 and D-1 representing the no-slip and full-slip cases respectively:    ð1  νm ÞEm γ max r 4Em π T¼ cos 2 θ þ 2F þ ð1  2νm ÞC þ 4D '' 4 rð1 þ ν Þ m ð1 þ νm ÞΔ ð15Þ

ð17Þ

Fð1  2νm Þð1  CÞ  21 ð1  2νm Þ2 þ 2
F½ð3  2νm Þ þ ð1  2νm ÞC þ C 52  8νm þ 6νm 2 þ 6  8νm Bobet method [5]

It should be noted that no solutions are proposed for the calculation of the bending moment and the normalised lining deflection for the no-slip assumption, but the method suggests using the corresponding equations for the full-slip assumption.

(a) Full-slip – dry soil T¼

Penzien [19] method

12ð1  νm Þ 3ð5  6νm Þ þ ð1  νm ÞF '

Gm γ max r sin 2θ

ð18Þ

M¼T r with F ' ¼ Em r 3 ð1 ν2m Þ=El Ið1  ν2l Þ Δdl 2U r ¼ Δdf reef ield rγ max

(a) Full-slip  3E IRn r γ max π cos 2 θ þ T ¼  l3 2 4 2 r ð1  νl Þ M¼ 

 3El IRn r γ max π cos 2 θ þ 4 2 r 2 ð1  ν2l Þ

ð9Þ

ð10Þ

where υm Þ αn ¼ 3 El I ð5  6 υm Þ=2 r 3 Gm ð1  ν2l Þ, Rn ¼ 7 4ðαð1n  þ 1Þ , and Gm is the shear modulus of the medium. Δdl ¼ Rn r γ max

ð11Þ

(b) No-slip T¼

 3El I R Δdf reef ield π cos 2 θ þ 2 3 4 r ð1  νl Þ

 3El I RΔdf reef ield π cos 2 θ þ M¼  2 2 4 2 r ð1  νl Þ

ð19Þ

ð20Þ

where U r is the radial displacement of the lining, which ignoring the thickness of the lining is defined as 1 þ νm ½1 þ 4ð1  νm ÞC 1 þ C 2 Gm γ max r sin 2θ Em ' 3  ð1  νm ÞF ' with C 1 ¼  , C 2 ¼  3ð1  2νm ÞF þ ð1  νm Þ 3ð5  6νm Þ þ ð1  νm ÞF ' 3ð5  6νm Þ þ ð1  νm ÞF ' Ur ¼

ð21Þ

(b) No-slip – dry soil T ¼ ð1 C 2 ÞGm γ max r sin 2θ

ð22Þ

1 M ¼  ð1 þC 1 þ C 2 ÞGm γ max r 2 sin 2θ 2

ð23Þ

ð1  νm Þ C ' þ ð1  νm Þ  ½ð1  νm ÞC ' þ 43=F ' where C 1 ¼  2 ð1  ν2 ÞC ' þ ð1  νm Þð3  2νm Þ þ ½ð1  νm Þð5  6νm ÞC ' þ 4ð3  4νm Þ3=F ' 2

ð12Þ

m

C2 ¼

1 ð1  νm ÞC ' 2  C 1 ½ð1  νm ÞC ' þ 4νm  3 ð1  ν ÞC ' þ 2 m

ð13Þ

with C ' ¼ Em rð1 ν2m Þ=El Að1  ν2l Þ ð1  ν2l Þ

and

F ' ¼ Em r 3 ð1  ν2m Þ=El I

S. Kontoe et al. / Soil Dynamics and Earthquake Engineering 66 (2014) 206–219

and 1 þ νm Ur ¼ ½1  2ð1  νm ÞC 1  C 2 Gm γ max r sin 2θ Em

ð24Þ

(c) Full-slip – undrained conditions T¼

6 6 þ ð1  νm ÞF '

Gm γ max r sin 2θ

M¼T r Ur ¼

ð25Þ ð26Þ

1 þ νm ½1 þ 2C 1 þ C 2 Gm γ max r sin 2θ Em

ð27Þ

with ' ' C 1 ¼  3  ð1  νm ÞF and C 2 ¼  ð1  νm ÞF ' 6 þ ð1  νm ÞF 6 þ ð1  νm ÞF ' (d) No-slip – undrained conditions T ¼  ð1  2C 2 ÞGm γ max r sin 2θ

ð28Þ

1 M ¼  ð1 þ 2C 1 þ 2C 2 ÞGm γ max r 2 sin 2θ 2

ð29Þ

with C1 ¼ 

ð1  νm Þ2 C ' þ ð1  νm Þ  ½ð1 νm ÞC ' þ 43=F ' ½2 þ ð1  ν ÞC ' ½ð1 ν Þ þ6=F '  m

C2 ¼

m

1 ð1  νm Þ2 C '  12=F ' 2 ½2 þ ð1  ν ÞC ' ½ð1 ν Þ þ 6=F '  m

m

and Ur ¼

1 þ νm ½1  2C 1  2C 2 Gm γ max r sin 2θ Em

ð30Þ

References [1] Amorosi A, Boldini D. Numerical modeling of the transverse dynamic behavior of circular tunnels in clayey soils. Soil Dyn Earthq Eng 2009;29:1059–72. [2] Argyroudis SA, Pitilakis KD. Seismic fragility curves of shallow tunnels in alluvial deposits. Soil Dyn Earthq Eng 2012;35:1–12. [3] Bardet JP, Ichii K, Lin CH. EERA, a computer program for equivalent linear earthquake site response analysis of layered soils deposits. Los Angeles: University of Southern California; 2000.

219

[4] Bobet A. Effect of pore water pressure on tunnel support during static and seismic loading. Tunn Undergr Space Technol 2003;18:377–93. [5] Bobet A. Drained and undrained response of deep tunnels subjected to farfield shear loading. Tunn Undergr Space Technol 2010;25:21–31. [6] Corigliano M, Scandella L, Lai CG, Paolucci R. Seismic analysis of deep tunnels in near field fault conditions: a case study in Southern Italy. Bull Earthq Eng 2011;9(4):975–95. [7] Darendeli MB. Development of a new family of normalised modulus reduction and material damping curves. ([Ph.D. thesis]). Austin: University of Texas; 2001. [8] Day RA, Potts DM. Zero thickness interface elements-numerical stability and application. Int J Numer Anal Methods Geomech 1994;18:689–708. [9] Hashash YMA, Hook JJ, Schmidt B, Yao JI-C. Seismic design and analysis of underground structures. Tunn Undergr Space Technol 2001;16:247–93. [10] Hashash YMA, Park D, Yao JI-C. Ovaling deformations of circular tunnels under seismic loading, and update on seismic design and analysis of underground structures. Tunn Undergr Space Technol 2005;20:435–41. [11] Hoeg K. Stresses against underground structural cylinders. J Soil Mech Found Div ASCE 1968;94(4):833–58. [12] Idriss IM, Sun JI. User's manual for SHAKE91. Davis: Center for Geotechnical Modeling, Department of Civil Engineering, University of California; 1992. [13] Kontoe S, Pelecanos L, Potts DM. An important pitfall of pseudo-static finite element analysis. Comput Geotech 2013;48:41–50. [14] Kontoe S, Zdravković L, Potts DM, Menkiti CO. On the relative merits of simple and advanced constitutive models in dynamic analysis of tunnels. Geotechnique0016-8505 2011;61:815–29. [15] Kontoe S, Zdravković L, Potts DM, Menkiti CO. Case study on seismic tunnel response. Can Geotech J 2008;45(12):1743–64. [16] Kouretzis GP, Sloan SW, Carter JP. Effect of interface friction on tunnel liner internal forces due to seismic S- and P- wave propagation. Soil Dyn Earthq Eng 2013;46:41–51. [17] Owen GN, Scholl RE. Earthquake engineering of large underground structures. (Report no. FHWA/RD-80/195). Washington, D.C: Federal Highway Administration and National Science Foundation; 1981. [18] Park K-H, Tantayopin K, Tontavanich B, Owatsiriwong A. Analytical solution for seismic-induced ovaling of circular tunnel lining under no-slip interface conditions: a revisit. Tunn Undergr Space Technol 2009;24:231–5. [19] Penzien J. Seismically induced racking of tunnel linings. Int J Earthq Eng Struct Dyn 2000;29:683–91. [20] Potts DM, Zdravković LT. Finite element analysis in geotechnical engineering: theory. London: Thomas Telford; 1999. [21] Potts DM, Zdravković LT. Finite element analysis in geotechnical engineering: application. London: Thomas Telford; 2001. [22] Sedarat H, Kozak A, Hashash YMA, Shamsabadi A, Krimotat A. Contact interface in seismic analysis of circular tunnels. Tunn Undergr Space Technol 2009;24:482–90. [23] Taborda D. Development of constitutive models for application in soil dynamics. ([Ph.D. thesis]). London: Imperial College; 2011. [24] Taborda D, Kontoe S, Zdravković L & Potts DM 2009. Application of cyclic nonlinear elastic models to site response analysis. In: Proceedings of the 1st international symposium on computational geomechanics, Juan-les-Pins, France; 29 April–1 May 2009. p. 956–66. [25] Wang JN. Seismic design of tunnels: a state-of-the-art approach. (Monograph). New York: Parsons, Brinckerhoff, Quade and Douglas Inc.; 1993.