- Email: [email protected]
Control of fault lay-out on seismic design of large underground caverns
Tunnelling and Underground Space Technology 50 (2015) 305–316
Contents lists available at ScienceDirect
Tunnelling and Underground Space Technology journal homepage: www.elsevier.com/locate/tust
Control of fault lay-out on seismic design of large underground caverns Saeid Ardeshiri-Lajimi a, Mahmoud Yazdani a,⇑, Arya Assadi Langroudi b a b
Tarbiat Modares University, Civil and Environmental Engineering Department, Tehran, Iran Newcastle University, School of Civil Engineering and GeoSciences, England, United Kingdom
a r t i c l e
i n f o
Article history: Received 29 August 2014 Received in revised form 12 May 2015 Accepted 7 July 2015 Available online 13 August 2015 Keywords: Fault Underground Plastic zone Serviceability Seismic
a b s t r a c t Although buried structures are generally believed to suffer a lesser degree of damage in the event of earthquake – than that of over-ground structures – significant damage has been widely reported to buried assets after major earthquakes, including the 1995 Kobe and the 2008 Wen-Chuan. Discontinuity is one key feature of rock as the supporting medium around subsurface excavated spaces. Joints, faults and bedding planes influence, by-and-large, the stability of structures made from/into rock. In particular, fault system around underground caverns such as hydropower house has a marked control on assets’ seismic stability. This study builds on the current understanding through vigorous numerical modelling of fault-structure system under seismic excitation. A parametric approach is followed to determine the most critical layout of a single fault crossing a benchmark cavern. Fault system is systematically broken down into several combinations of dips and intersection points with cavern wall. For each case, a nonlinear dynamic analysis is conducted. To simulate the discontinuous medium, the hybrid finite difference – discrete element code CA2 (Continuum Analysis 2 dimensional) is implemented. The work showed that, similar to static conditions, fault influences the seismic stability of underground caverns through a tendency in extending the plastic zones and increasing displacements as well as asymmetric distribution of the latter and the former in rock medium. A 40–50° dip, single-point-intersection-on-crown k0 = 1 fault layout renders the most critical combination from both ultimate and serviceability limit states perspective. Under earthquake loading conditions however, the critical limit states condition took place for single fault intersected the cavern at heel and sidewall. The latter critical condition led to the tensile failure of cavern right sidewall. For faults intersecting the carven crown and having a k0 = 0.5, collapse would be more likely as fault dip increases. Collapse would be less likely with increasing dip for k0 = 0.5 fault crossing the bed and sidewall of caverns. Ó 2015 Elsevier Ltd. All rights reserved.
1. Introduction Common to nearly all reported damages to underground structures in the event of earthquake is the pronounced contribution of nearby faults to structural failure, highlighting the need for a better understanding of the interaction between rock discontinuities and structures they accommodate within (Zhang et al., 2013; Yashiro et al., 2007; Wang et al., 2001). Li (2012) emphasized on the significance of secondary fractures of seismogenic faults in proposing installation of reinforced concrete as secondary lining should tunnel alignment intersect such features. He recommended the lining to mantle the active fault zone and extend beyond, up to a minimum distance of 5 m. The micromechanics of fracture growth
⇑ Corresponding authors. E-mail addresses: [email protected] (S. Ardeshiri-Lajimi), myazdani@ modares.ac.ir (M. Yazdani), [email protected] (A. Assadi Langroudi). http://dx.doi.org/10.1016/j.tust.2015.07.002 0886-7798/Ó 2015 Elsevier Ltd. All rights reserved.
under high geo-energies and implications were addressed within a fractal context in a more recent work of Assadi et al. (2014). Ichimura et al. (2012) contributed to the understanding of explicit links between buried structures and faults through developing a finite element tool, which allows the measurement of seismic response of underground structures. Building on the latter platform, they demonstrated the reliability of fault-structure system analysis as a function of discontinuity geometry and orientation. Perhaps one of the earliest works in appreciation of faults as the origin of subsurface void space collapse is that of Brekke and Selmer-Olsen (1965), which discussed the matter from a geochemical perspective. They examined a number of failure incidents in Norwegian tunnels due to the smectite occurrences in between supporting rock discontinuities. From geomechanical standpoint however, fault-related instabilities include the collapse of large rock mass volumes within the fault zone into neighbouring caverns, particularly underground openings excavated through large thick fault zones. Instability is also likely after slip surface
306
S. Ardeshiri-Lajimi et al. / Tunnelling and Underground Space Technology 50 (2015) 305–316
formation along fault planes and consequent downfall of blocks or wedges – intersected by fault and other minor joints – into adjacent caverns (Nagelhout and Roest, 1997; Brekke and Howard, 1973). Other major contributions to association of instabilities and faults in rock mediums – a pivotal issue that have engaged both practical engineers and researchers – include Heap et al. (2010), Hao and Azzam (2005), Hao and Azzam (2001), Hashash et al. (2001), Tang and Kaiser (1998) and Tang (1997). Another stream of stability analysis of discontinuum fractured mediums concerns determination of joint properties from captured displacements, one latest being the wok of Yazdani et al. (2012). The developed knowledge however continues to lack a concrete theoretical solution to measurement of underground cavern stability. One approach to fill the latter shortfall is rigour benchmark modelling and cross-disciplinary researching, collating and analysing the generated big datasets in the framework of probabilistic analysis. Pursuing such approach shall positively impact the profession in enhancing sustainable and reliable support systems design for future buried infrastructures. Seismic response of underground structures is fundamentally different with that of surficial structures, in part due to their full enclosure in soil/rock medium and their often significant length dimension e.g. road tunnels. For buried structures, hence, to be designed fully resistant to seismic excitation, implementation of numerical modelling is not just common but an inevitable practice. Over the past years, various methods have been developed and extensively applied across the profession, yet, problem oriented numerical approaches are scarce. The relaxation of the latter is deemed here possible through the use of CA2 code (Fakhimi, 1998). The present work aims at evaluating the impacts of faults on induced plastic zones i.e. ultimate limit states as well as displacement re-distribution in subsurface openings i.e. serviceability limit states, in the event of earthquake. The main focus will be on measuring the effects of fault layout on the stability of caverns they cross – a function of fault dip and location. To simulate the benchmark discontinuity scenarios, the finite difference – discrete element code CA2 (for large static and dynamic strains) is used throughout, which particularly benefits in generation of tailored interface elements. The CA2 code is used to model the seismic excitation of a benchmark large underground cavern, seeking the critical position/dip of a single fault crossing the cavern section. That understanding will inform practitioners in taking tailored mitigation measures in lining design with respect to the host environment morphology.
2. Faults and modelling approaches Faults in definition are the commonly encountered large geological discontinuities in hard rock masses. Most severe underground instabilities are found to be closely associated with fault systems present nearby. Hao and Azzam (2001) in an investigation into the influence of fault on stability of large-scale underground caverns under static conditions demonstrated that faults may expand the plastic zones, give rise to displacements and distribute both displacements and plastic points in an asymmetrical fashion across the rock mass and around excavated spaces. The early types of joint elements proposed for simulating rock joints or faults in rock mechanics (Griffiths, 1985; Pande and Sharma, 1979; Goodman et al., 1968) are still widely implemented to model joints and interfaces across a broad range of geomechanical applications. Schweiger et al. (1990) employed the thin-layer joint element model to assay the fault thickness and stiffness effects on underground tunnel stability. Lei et al. (1995) used a contact-friction interface element to investigate the fault influence on the stress redistribution, bending moment and normal force across lining. Steindorfer (1998) utilised the ratio of longitudinal
to vertical displacement to assess the effect of fault zones in the process of tunnel excavation. More recent developments include the works of Bograd et al. (2011) and Mayer and Gaul (2007). Within the present work, joints are modelled as normal and shear springs with normal and tangential stiffness i.e. Kn and Ks, respectively between any pair of planes in contact (see Scheme 1). The model allows the explanation of principle fault characteristics including slipping, opening and compression behaviour. The model offers an appreciation of the fact that faults exhibit zero tension in the normal direction and follow Coulomb’s friction law for shear. The CA2 code employed is a two-dimensional explicit finite difference programme using time-marching procedures at very small steps through which the interface is divided into many differential elements. The incremental relative displacement vector at the contact point is resolved into the normal and shear directions and the total normal (Fn) and shear (Fs) forces are determinable using Eqs. (1) and (2). DtÞ F ðtþ ¼ F nðtÞ K n DU nðtþDt=2Þ L n DtÞ F ðtþ ¼ F sðtÞ K s DU sðtþDt=2Þ L s
ð1Þ
Eq. (1) is expressed in form of matrix as of,
DF n DF s
¼
Kn 0
0 DU n L DU s L Ks
ð2Þ
where the stiffness K n and K s equate stress over displacement, L is the length of interface element, DU n is the displacement increment of the interface element in the normal direction and DU s is the displacement increment of a segment of interface element in the tangential direction. The status of fault is expressed within the model by the following criteria: (a) Tensile yield condition:
F t ¼ rt t where rt is the tensile stress acting on the fault and t is the tensile strength of the fault. For Ft P 0, the interface breaks and the shear and normal forces are set to zero. The tensile strength of faults is generally agreed as of zero. (b) Shear yield condition: Shear stress along faults can be expressed using Mohr–Coulomb failure criterion as of:
F s max ¼ c l þ F n tan u where c is cohesion intercept along the interface, l is effective contact length and u is friction angle of interface. For |Fs| P Fs max, then Fs = Fs max. The rock is modelled as a bonded particle system. The rock grains are assumed as circular cylinders that interact through normal and shear springs. Each discrete element in CA2 is a circular cylinder. The cylinder–cylinder and cylinder–wall (i.e. interface) interactions are modelled by normal and shear springs, normal
Scheme 1. Fault model. The interface represented by side A and side B, which are connected by normal (kn) and shear (ks) stiffness springs. (S = slider, T = tensile strength and LN = length associate with grid-point N).
S. Ardeshiri-Lajimi et al. / Tunnelling and Underground Space Technology 50 (2015) 305–316
and shear bonds, and a friction coefficient. Therefore, the micromechanical constants for cylinder–cylinder or cylinder–wall interactions are the normal and shear spring constants (kn and ks), the normal and shear bonds (nb and sb), and the friction coefficient (l). This simple contact bond model behaves elastically under the applied stresses until the tensile normal force or shear force at a contact point exceeds the normal or shear bond leading to tensile or shear failure at that contact point. A failed contact point loses its bond permanently but still it can carry shear forces through a Coulomb law using the friction coefficient l if it is subjected to a compressive load – also see Fakhimi and Gharahbagh (2009) and Hosseinpour (2008). 3. Modelling and boundary conditions To assess the influence of fault system on seismic stability of underground openings, a benchmark model is adopted – see Fig. 1. The underground opening considered in this model represents many forms of hydroelectric schemes such as underground power houses and surge chambers. As such, adopted model’s cross-section has a crown of arch in shape and sidewalls sizing 20 m wide by 30 m high. To minimize the impact of uninteresting boundary conditions on modelling outputs, the size of model in each direction is adopted eight times the cavern dimension. The fault is extended far enough from the excavation and the fault displacements dissipate at model borders, replicating the remote rock masses. Under static conditions, vertical boundaries are fixed along horizontal direction and zero displacement condition is applied across the base of the model. The effect of overlying (i.e. overhead) rock weight is considered on the model top boundary by a normal stress. Adopting displacement
307
boundary condition instead of stress boundary condition on the vertical boundaries of the model is the premise of obtaining ideal numerical results for the fault-related problem (Hao and Azzam, 2001). A fault with different dips crosses through the rock mass in the vicinity of the cavern i.e. intersects the cavern. The fault strike is presumed parallel or nearly parallel to the alignment of the cavern and hence a plane strain state governs. The latter is well established as the most unfavourable condition in the context of underground caverns stability (Hoek and Brown, 1982). 3.1. Model considerations under static loads Throughout the analysis, the underground cavern is assumed to be at a depth of 300 m below the ground level. The in-situ stress state prior to cavern excavation is sourced from the weight of rock masses at rest. The vertical stress varies linearly with depth while the horizontal stress is calculated as of rh = k0rv, where k0 is the ratio between the horizontal and vertical stresses. Intact rock is considered as an elastic–perfectly plastic material that follows Mohr–Coulomb failure criterion. Slip failure of the fault also follows the Mohr–Coulomb criterion. Input parameters for the analysis purposes are determined in a way to satisfy following requirements: The properties of rock mass and fault need to represent the most unfavourable fault influence, while reflecting real characteristics of the materials. As such, the elastic modulus of the rock mass is taken as of 20 GPa, the internal friction angle is fixed at c. 35°, while small values for tensile strength and cohesion are adopted. These specified geotechnical parameters represent sparse jointed rock masses with low-to-medium strength in nature. Since the emphasis of this analysis is to discuss the effects of key fault parameters, only fault stiffness, tensile strength and cohesion are needed to be specified. The intention is to study a typical weak fault with a weathered, persistent and smooth surface in absence of tensile strength and cohesion and thin thickness. Thereby, shear dilation, strain hardening and softening of fault are neglected in the model. The latter parameters have no control on the general discussion of the problem throughout the analysis. Within the scope of current work, a comparative analysis of the effect of different fault parameters is followed throughout the coming sections, leaving the absolute prediction of a practical problem out of the interest of the present article. The adopted geotechnical parameters for the rock mass and fault render a stable model under static conditions. 3.2. Model considerations under dynamic loads The relation between the elastic modulus in dynamic loading environments, Ed, and the elastic modulus under static loading circumstances, Es, is given by Hayashi and colleagues in 1973, reading:
Ed ¼ ð1:3 1:7ÞEs
ð3Þ
The magnitude of Ed is obtained through field full-scale dynamic loading test, with the Poisson’s ratio md ranging from 0.2 to 0.3. According to Price et al. (1969), the tensile strength increases up to 1.5 times that of the static value in the event of quick loading. The mechanical properties of rock mass and fault used in the present research under static and dynamic conditions are summarized in Table 1. Fig. 3 shows the baseline-corrected response spectra graph i.e. input motion. 3.3. Dynamic analysis
Fig. 1. Sketch of the model for this numerical analysis.
Following the static analysis under the static effective stress conditions, below steps are adopted to undertake the dynamic analysis:
308
S. Ardeshiri-Lajimi et al. / Tunnelling and Underground Space Technology 50 (2015) 305–316
Table 1 Input parameters of rock mass and fault for the numerical modelling. Material properties 3
Density q (kg/m ) Internal friction angle / Cohesion c (MPa) Poisson ratio Normal stiffness kn (GPa/m) Shear stiffness ks (GPa/m) Es (GPa) Ed (GPa) Tensile strength rt (kPa)
Rock mass
Fault
2600 35 2 0.2 –
24 0 – 20 2 – – –
20 25 1.7
Fig. 2. CA2 model boundary conditions for dynamic simulation.
Step 1: In order to use horizontal component of the acceleration record as input event, it is essential to accomplish baseline correction on above-mentioned histories first. As a result, after integrating the acceleration histories twice, the obtained displacement does not reach zero value at the end of the motion. One approach for correcting this error is to add a low frequency wave to earthquake loading so that displacement reaches zero in the end. Step 2: Another issue in dynamic loading which is of considerable importance is the wave propagation in the model. Generally, in performing dynamic analyses the presence of error in the form of wave propagation is likely to happen. The finite difference grid dimensions were selected by taking into account the maximum frequency (f) of the shear wave that the model could logically respond to during earthquake loading. Consequently, this prevents the incorrect propagation of the waves. The frequency is determined by the following equation (Kuchlemeyer and Lysmer, 1973):
f ¼
Cs 10Dl
ð4Þ
where C s is the shear wave velocity of the soil and Dl is the largest grid zone size in the model. Referring to equation, a uniform zone size of 1 m 1 m was selected. Since the lowest shear wave velocity in the model belongs to the soil deposits, the highest admissible frequency for a propagating shear wave is 5 Hz. Therefore, the input earthquake record shall be filtered by a low pass filter to remove frequency components higher than 5 Hz. A frequency of 4 Hz was ultimately selected as the low pass filter for the reduction in shear wave velocity which may occur due to plastic flow during seismic loading. Step 3: Application of dynamic boundary conditions: The boundary conditions used in the dynamic analysis are illustrated in Fig. 2. Quiet boundaries were used on all outside boundaries of model. These boundaries prevent reflection of outgoing seismic waves back into the model. Quiet boundaries were combined with free-field boundaries on the vertical outside boundaries that prevent distortion of vertically propagating plane waves along the boundaries. Dynamic loading was applied at the bottom of model, as propagating vertically upwards (Pakbaz and Yareevand, 2005). The record of accelerogram of 1999 Kocaeli, Turkey earthquake with magnitude of 7.4 in scale of Richter (Fig. 4) has been used. Step 4: Rayleigh damping, which consists of two viscous elements, is conventionally used in the numerical analyses herein. The two elements of Rayleigh damping are both frequency dependent (Lysmer and Kuchlemeyer, 1969). By choosing a mid-frequency at which the combined effects of the two elements cancel out, it is possible to have a damping that is nearly independent of frequency over a fairly wide range of frequencies on either side of the mid-frequency (White et al., 1977). The mid-frequency is usually chosen in the range between the natural frequency of the model and the predominant frequency of the input motion. Rayleigh damping was assigned to each element of the model in the mid-frequency. In analyses that use one of the plasticity constitutive models (e.g. Mohr– Coulomb), a considerable amount of energy dissipation can occur during plastic flow (viscous damping). Thus, for many dynamic analyses that involve large-strain, only a minimal percentage of damping (e.g. 5%) may be required (Itasca, 2002). Therefore, the damping ratio was assumed to be 5% in the analyses. Dynamic analyses were performed and the results were extracted for interpretation and further assessment. 4. Numerical results and discussion The effect of earthquake on underground opening depends on various parameters including peak acceleration, magnitude and duration of earthquake, the relative rigidity between underground opening and ground and in-situ stress state. Furthermore, for
Fig. 3. Acceleration record of kocaeli earthquake.
S. Ardeshiri-Lajimi et al. / Tunnelling and Underground Space Technology 50 (2015) 305–316
309
underground caverns crossed by a single fault, additional parameters including intersection of fault with cross section and mechanical properties of fault have control on the stability. In this section the effects of fault orientation and the relative location of its intersection with cavern perimeter together with the PGA (Peak Ground Accelerations) of earthquake loading are examined on seismic stability of cavern. To quantify the stability analysis, two indices, namely the extended plastic area in rock mass (analogous to ultimate limit state) and the maximum displacement in tunnel wall or roof (analogous to serviceability limit state) have been adopted. Quantified outputs are obtained for a combination of fault dips, PGAs and horizontal to vertical stress ratio (k0). Discussions are lastly built on comparisons between displacements and area of plastic zones before and after application of earthquake loading is undertaken.
Fig. 6. Variation of maximum vertical displacement of roof versus fault dip for different PGA in case 3: fault intersects cavern floor in the middle point and also right cavern wall (k0 = 1; u = 24°).
Fig. 4. Variation of maximum vertical displacement of roof versus fault dip for different PGA in case 1: fault intersects cavern roof (k0 = 1; u = 24°).
Fig. 7. Variation of maximum vertical displacement of roof versus fault dip for different PGA in case 4: fault intersects cavern floor in the corner point (k0 = 1; u = 24°).
4.1. Maximum vertical displacement
Fig. 5. Variation of maximum vertical displacement of roof versus fault dip for different PGA in case 2: fault intersects left cavern wall in the middle point and also the roof (k0 = 1; u = 24°).
Figs. 4–7 illustrate the variation of maximum rock mass vertical displacement against fault dips for varying PGAs and a suite of faults crossing the cavern periphery. A k0 = 1 and fault frictional angle of u = 24° was considered in building the models. The presented values represent displacements induced by earthquake loading at its ultimate excitation and are independent from initial displacement originated from in-situ stresses. The maximum vertical displacement took place at cavern roof, lending evidence to the remarkable impact of fault when intersecting the cavern at a single point on crown. In general, maximum roof settlement occurred for fault dips ranging from 40° to 50°. Other key observations include: (1) Vertical displacement on crown appeared more sensitive to PGA for accelerations greater than 0.15 g. (2) Dip-displacement
310
S. Ardeshiri-Lajimi et al. / Tunnelling and Underground Space Technology 50 (2015) 305–316
Fig. 8. Relationship of the area of rock mass plastic zones varying with fault dips and horizontal stress ratio for case 1, PGA = 0.15 g & friction angle of fault u = 30°.
Fig. 10. Variation of plastic area extent versus fault dip for different PGA in case 1: fault intersects cavern roof (k0 = 1; u = 24°).
Fig. 9. Relationship of the area of rock mass plastic zones varying with fault dips and horizontal stress ratio for case 3, PGA = 0.15 g & friction angle of fault u = 30°.
Fig. 11. Variation of plastic area extent versus fault dip for different PGA in case 2: fault intersects left cavern wall in the middle point and also the roof (k0 = 1; u = 24°).
trend tool a bi-modal distribution when fault intersected the cavern at crown and side wall both. (3) Predicted maximum displacement was 15–20 times greater when fault intersected the crown at a single point rather than any other relative position. Examining the vertical displacement values in Figs. 4–7 infers a 40–50° dip, single-point-intersection-on-crown fault layout renders the most critical combination in terms of serviceability.
the critical dip. Thereby, a 0.15g PGA was adopted throughout this course of analysis. Original (prior to excavation) vertical in-situ stress was assumed to equate the weight of the rock mass below the ground level – which varies linearly with depth. Horizontal in-situ stress was taken as rh = k0 rv, where k0 is the ratio between horizontal stress and vertical stress. k0 was varied at constant vertical stress (i.e. constant buried depth of the benchmark cavern) for every analysis case. Analysis inferred that the critical dip values discussed in Section 4.1 are strongly a function of k0. For Case 1 (where the fault system intersects the crown at single point) plastic zone area gained remarkably high values for k0 = 0.5 and followed an upsurge trend with dip, dissimilar to the unimodal pattern captured for k0 = 1. As such, the conclusion gathered in Section 4.1 may be read as: a 40–50° dip, single-point-intersection-on-crown k0 = 1 fault
4.2. In-situ stress state The control of fault on the plastic zones surrounding the cavern for varying in-situ stress state was investigated through parametric analysis (see Fig. 8 and 9). Overall, PGA was found directly proportional to the area of plastic zone and displacement magnitude, with no potential control on
S. Ardeshiri-Lajimi et al. / Tunnelling and Underground Space Technology 50 (2015) 305–316
311
Fig. 12. Variation of plastic area extent versus fault dip for different PGA in case 3: fault intersects cavern floor in the middle point and also right cavern wall (k0 = 1; u = 24°).
Fig. 14. Relationship of the area of rock mass plastic zones varying with fault dips and fault location, (a) under static loads and (b) under dynamic loads for PGA = 0.15 g, horizontal stress ratio k0 = 1 & friction angle of fault u = 24°. Fig. 13. Variation of plastic area extent versus fault dip for different PGA in case 4: fault intersects cavern floor in the corner point (k0 = 1; u = 24°).
layout renders the most critical combination from a serviceability perspective. For faults intersecting the carven crown and having a k0 = 0.5, collapse would be more likely as fault dip increases. Collapse would be less likely with increasing dip for k0 = 0.5 fault crossing the bed and sidewall of caverns.
reworded as: a 40–50° dip, single-point-intersection-on-crown k0 = 1 fault layout renders the most critical combination from both ultimate and serviceability limit states perspective. Critical dip was measured as of 50° for when fault intersects the cavern at the side wall and in the middle point (Figs. 11 and 12), and 30° for when fault intersects the cavern at the lower floor corner (Fig. 13). 4.4. Effect of fault relative location (position)
4.3. Area of plastic zone For four fault locations, the area of plastic zone – surrounding cavern was plotted against fault dip in Figs. 10–13. Plots were re-produced for four PGAs. Plastic area appeared a function of fault dip, the maximum content of which being strongly a factor of fault location. For faults crossing cavern crown (established as the critical position on Sections 4.1 and 4.2), plastic zone gently expands with fault dip, although the maximum plastic area occurs at a critical 40° dip. To this end, the main line of conclusion in Section 4.2 may be
As depicted in Fig. 1, four possible fault sites – with respect to the cavern – were considered in the current work. To facilitate the discussion and better coherency with the previous arguments, Case 1 is defined for when the fault crosses (is tangent to) the cavern roof. In Case 2, the fault intersects the mid-height of the left sidewall. In Case 3, the fault cuts the cavern crown and through the centre of the bottom surface. In Case 4, the fault emerges out close to the heel of the cavern’s sidewall. For all the cases, the fault dips are assumed within the range of 20–80° for consistent mechanical properties of rock mass and fault.
312
S. Ardeshiri-Lajimi et al. / Tunnelling and Underground Space Technology 50 (2015) 305–316
Fig. 14 plots the plastic area in the rock mass against fault dip for a suite of fault locations. An interesting drawing from the results is the resembling pattern of plastic area – fault dips of Case 1 and Case 3, with Case 2 and Case 4, respectively, where maximum plastic area is attained as the fault dips approach 40° and 60°, respectively. Under seismic excitation and from ultimate limit state perspective (i.e. area of plastic zone), the critical condition took place for single fault intersected the cavern at heel and sidewall. This is dissimilar with the static conditions in which critical conditions were associated with the single fault crossing the cavern’s crown. 4.5. Maximum shear displacement of fault Fig. 15 shows the variation of maximum shear displacement of fault against dip. Independent of fault location, maximum shear displacement took place at 50° dip under all PGAs. 4.6. Comparing static and dynamic results
Fig. 15. Maximum shear displacement of fault.
Maximum values of vertical displacements on cavern roof, horizontal displacements on cavern walls and area of plastic zone under dynamic and static condition are compared in pairs and discussed in Figs. 16–18:
Fig. 16. Comparing static and dynamic values of maximum vertical displacements of cavern roof, for PGA = 0.15 g horizontal stress ratio k0 = 1 & friction angle of fault u = 24°.
S. Ardeshiri-Lajimi et al. / Tunnelling and Underground Space Technology 50 (2015) 305–316
4.6.1. Maximum vertical displacement Fig. 16 shows the variation of maximum rock mass vertical displacement against fault dip under static and dynamic condition for PGA = 0.15 g. A k0 = 1 and u = 24° fault frictional angle was adopted throughout. The annotation ‘static’ on plot represents the displacement values induced by in-situ stresses after excavations, while ‘dynamic’ tag represents displacements induced after earthquake loading at its final moment of excitation, independent from (excluding) the initial displacement produced by in-situ stresses. The annotation ‘sum’ stands for the total of static and dynamic values. Following drawings were made through examining the plots: (1) Vertical displacement of cavern crown is less dependent on dip under static loading conditions than dynamic excitation. (2) Critical dip, in terms of crown vertical displacement i.e. serviceability limit state sits in the 40° to 50° interval, predominantly controlled by strains in dynamic environment. (3) Maximum vertical displacement is likely to take place for fault crossing the top of the cavern. For other fault locations, the maximum value falls at least 10 times shorter than the latter, reinforcing the need for particular mitigation measures to be taken for excavation immediately beneath fault lines. (4) No link between
313
the trend of maximum vertical displacement under static and dynamic conditions could be identified. 4.6.2. Maximum horizontal displacement Figs. 17 and 18 illustrate the variation of maximum horizontal displacements of cavern left and right sidewall against fault dip under static and dynamic (PGA = 0.15 g) conditions. A k0 = 1 and u = 24° fault frictional angle was adopted throughout. Plots suggest no correlation between static and dynamic values. Due to the tensile failure of cavern right sidewall on Case 4, horizontal displacement-dip plot could not possibly be generated. Key observations include: (1) For fault location intersecting the left side-wall, horizontal displacement on the left-wall under static loading appeared to be strongly a function of dip. In contrary, left wall horizontal displacement followed a fairly plateau trend against dip for faults with no intersection with the left side-wall. Similar pattern follows for the right side-wall. (2) No link could be established between the fault location (i.e. Case) and the occurrence of maximum horizontal displacement. In Fig. 17, maximum horizontal displacement occurred in Case 1 when fault crossed the cavern crown, while in Fig. 18, maximum horizontal
Fig. 17. Comparing static and dynamic values of maximum horizontal displacements of cavern left wall, for PGA = 0.15 g horizontal stress ratio k0 = 1 & friction angle of fault u = 24°.
314
S. Ardeshiri-Lajimi et al. / Tunnelling and Underground Space Technology 50 (2015) 305–316
displacement occurred in Case 3 for fault crossing the mid-height of the right side-wall.
4.7. Area of plastic zone Fig. 19 illustrates the variation of plastic area within the rock mass around the excavated space against fault dips for static and dynamic conditions and for four fault positions (intersection of fault with cavern periphery). Plots suggest the existence of a link between the fault-induced area of rock mass plastic zone under earthquake loading and static loading. Whilst the plots represent a PGA of 0.15 g for the dynamic analysis, same trends were re-produced for PGAs of 0.1 g, 0.20 g and 0.25 g. This was confirmed further in Ardeshiri (2008), where input PGA values were extended to six. Eq. (5), illustrates the established relationship between area of plastic zones before and after earthquake loading. The equation was successfully validated for two other earthquake waves i.e. Loma Prieta and Duze, with different predominate period values (Ardeshiri, 2008).
ðAPa Þd ¼ ð1 þ 0:25e10:92amax Þ ðAPa Þs
ð5Þ
where ðApa Þd is the area of plastic zones after earthquake loading, amax is maximum earthquake acceleration and ðApa Þs is area of plastic zones under in-situ stress condition.
5. Conclusion
Fig. 18. Comparing static and dynamic values of maximum horizontal displacements of cavern right wall, for PGA = 0.15 g horizontal stress ratio k0 = 1 & friction angle of fault u = 24°.
In this paper, CA2 code was used to revisit, in a systematic manner, the effects of fault lay-out on the dynamic and static response of a benchmark cavern excavated into a rock mass. Numerical analysis indicates that the influence of a fault to the underground cavern stability is different from its inherent properties. Similar to static conditions, fault influences the seismic stability of large underground caverns through a tendency in extending the plastic zones and increasing displacements as well as asymmetric distribution of the latter and the former in rock medium. This however is strongly a function of fault dip and location (i.e. lay-out), together with the situ stress state. Fault dip in relation to the underground cavern is proved to have key control on cavern’s response to earthquake forces. The present work confirmed a nearby 40° and 50° fault dip to be most critical. For a suite of k0 and the typical configuration of cavern adopted in the present work, an earlier expectation was met in the attained decreasing area of rock mass plastic zone with an increase in the horizontal stress ratio. The significance of fault in seismic stability of underground cavern was further highlighted in the high dependency of serviceability e.g. vertical displacement of cavern crown on dip under dynamic excitation rather than static stresses. When fault system crosses the cavern section, sharp edges could form on cavern section. Faced to earthquake waves, these sharp edges can generate a state of tensile failure. Given the local nature of such damage on cavern section however, sufficient flexibility could provide the required resistance of cavern against seismic loads. An equation for determining the area of the earthquake-induced plastic zones was proposed, which informs engineers in improving the design of support system in seismic regions. The work can be extended by interested workers in bringing into parametric analysis other key parameter, in addition to fault parameters, including PGV (Peak Ground Velocity), frequency content of earthquake wave, depth of excavated cavern.
S. Ardeshiri-Lajimi et al. / Tunnelling and Underground Space Technology 50 (2015) 305–316
315
Fig. 19. Comparing static and dynamic area of plastic zones values of cavern, for PGA = 0.15 g horizontal stress ratio k0 = 1 & friction angle of fault u = 24°.
References Ardeshiri, S., 2008. Influence of Discontinuities on Seismic Behaviour of Large Underground Caverns. Master of Science Tarbiat Modares Univfersity – TMU. Assadi, A., Jefferson, I., O’Hara-dhand, K., Smalley, I., 2014. Micromechanics of quartz sand breakage in a fractal context. Geomorphology 211, 1–10. Bograd, S., Reuss, P., Schmidt, A., Gaul, L., Mayer, M., 2011. Modeling the dynamics of mechanical joints. Mech. Syst. Signal Process. 25, 2801–2826. Brekke, T.L., Howard, T.R., 1973. Functional Classification of Gouge Materials from Seams and Faults in Relation to Stability Problems in Underground Openings. US Bureau of Mines (AROA) – Con. No. H0220022. Brekke, T.L., Selmer-Olsen, R., 1965. Stability problems in underground constructions caused by montmorillonite-carrying joints and faults. Eng. Geol. 1, 3–19. Fakhimi, A., 1998. Theory and User Manual of CA2 Computer Programe – Report no. 262. Building and Housing Research Centre, Tehran, Iran. Fakhimi, A., Gharahbagh, E.A,. 2009. Discrete element modeling of the influence of void size and distribution on the mechanical behavior of rock. In: Diederichs, M., Grasselli, G., (Eds.), Proceedings of the 3rd CANUS Rock Mechanics Symposium. Toronto, Canada. Goodman, R.E., Taylor, R.L., Brekke, T.L., 1968. A model for the mechanics of jointed rock. J. Soil Mech. Foundation Division, ASCE 94, 637–659. Griffiths, D.V., 1985. Numerical modelling of interfaces using conventional finite elements. In: Kawamoto, T., Ichikawa, Y., (Eds.), Proceedings of the 5th International Conference on Numerical Methods in Geomechanics, Nagoya, Japan. pp. 837–844. Hao, L.H., Azzam, R.,2001. Analysis of fault effects on the stability of large underground caverns. In: ISRM International Symposium – 2nd Asian Rock Mechanics Symposium. Beijing, China.
Hao, Y.H., Azzam, R., 2005. The plastic zones and displacements around underground openings in rock masses containing a fault. Tunnell. Underground Space Technol. 20, 49–61. Hashash, Y.M.A., Hook, J.J., Schmidt, B., I-Chiang Yao, J., 2001. Seismic design and analysis of underground structures. Tunnell. Underground Space Technol. 16, 247–293. Heap, M.J., Faulkner, D.R., Meredith, P.G., Vinciguerra, S., 2010. Elastic moduli evolution and accompanying stress changes with increasing crack damage: implications for stress changes around fault zones and volcanoes during deformation. Geophys. J. Int. 183, 225–236. Hoek, E., Brown, E.T., 1982. Underground Excavations in Rock. Institute of Mining and Metallurgy, London. Hosseinpour, H., 2008. Numerical and experimental evaluation of the effect of an oversize particle on direct shear test results. New Mexico Bureau Geol. Min. Resour. Ichimura, T., Hori, M., Quinay, P.E., Wijerathne, M.L.L., Suzuki, T., Noguchi, S., 2012. Comprehensive numerical analysis of fault-structure systems – computation of the large-scale seismic structural response to a given earthquake scenario. Earthquake Eng. Struct. Dyn. 41, 795–811. Itasca, F., 2002. Fast Lagrangian Analysis of Continua – Version 4.0 user guide. Thrasher Square East, 708: Itasca Consulting Group Inc. Kuchlemeyer, R.L., Lysmer, J., 1973. Finite element method accuracy for wave propagation problems. J. Soil Mech. Found. Div., ASCE 99, 421–427. Lei, X.Y., Swoboda, G., Zenz, G., 1995. Application of contact – friction interface element to tunnel excavation in faulted rock. J. Comput. Geotech. 17, 349–370. Li, T., 2012. Damage to mountain tunnels related to the Wenchuan earthquake and some suggestions for aseismic tunnel construction. Bull. Eng. Geol. Environ. 71, 297–308. Lysmer, J., Kuchlemeyer, R.L., 1969. Finite dynamic model for infinite media. J. Eng. Mech. Div., ASCE 95, 859–877.
316
S. Ardeshiri-Lajimi et al. / Tunnelling and Underground Space Technology 50 (2015) 305–316
Mayer, M.H., Gaul, L., 2007. Segment-to-segment contact elements for modelling joint interfaces in finite element analysis. Mech. Syst. Signal Process. 21, 724– 734. Nagelhout, A.C.G., Roest, J.P.A., 1997. Investigating fault slip in a model of an underground gas storage facility. Int. J. Rock Mech. Min. Sci. 34, 212. Pakbaz, M.C., Yareevand, A., 2005. 2-D analysis of circular tunnel against earthquake loading. Tunnell. Underground Space Technol. 20, 411–417. Pande, G.N., Sharma, K.G., 1979. On joint interface elements and associated problems of numerical ill-conditioning. Int. J. Numer. Analyt. Methods Geomech. 3. Price, D.G., Malkin, A.B., Knill, K.J., 1969. Foundations of multi-storey blocks on coal measures with special reference to old workings. Q. J. Eng. Geol. 1, 271–322. Schweiger, H.F., Haas, W., Handle, E., 1990. A thin layer element for modelling joints and faults. In: Rossmanith, H.P., (Ed.), Proceedings of International Conference on Jointed and Faulted Rock, 1990 Vienna, Austria. pp. 559–564. Steindorfer, A., 1998. Short term prediction of rock mass behaviour in tunelling by advanced analysis of displacement monitoring data. In: Schuvert. Riedmuller, Semprich (Eds.), Schriftenreihe der Gruppe Geotechnik, Heft 1. Tang, C., 1997. Numerical simulation of progressive rock failure and associated seismicity. Int. J. Rock Mech. Min. Sci. 34, 249–261.
Tang, C.A., Kaiser, P.K., 1998. Numerical simulation of cumulative damage and seismic energy release during brittle rock failure—part I: fundamentals. Int. J. Rock Mech. Min. Sci. 35, 113–121. Wang, W.L., Wang, T.T., Su, J.J., Lin, C.H., Seng, C.R., Huang, T.H., 2001. Assessment of damage in mountain tunnels due to the Taiwan Chi–Chi Earthquake. Tunnell. Underground Space Technol. 16, 133–150. White, W., Lee, I.K., Valliappan, S., 1977. Unified boundary for finite dynamic models. J. Eng. Mech. Div., ASCE 103, 949–964. Yashiro, K., Kojima, Y., Shimizu, M., 2007. Earthquake damage to tunnels in Japan and case studies of railway tunnels in the 2004 Niigataken-Chuetsu earthquake. Quarterly Report Railway Tech. Res. Inst. 48, 136–141. Yazdani, M., Sharifzadeh, M., Kamrani, K., Ghorbani, M., 2012. Displacement-based numerical back analysis for estimation of rock mass parameters in Siah Bisheh powerhouse cavern using continuum and discontinuum approach. Tunnell. Underground Space Technol. 28, 41–48. Zhang, J.W., Mei, Z.R., Quan, X.J., 2013. Failure characteristics and influencing factors of highway tunnels damage due to the earthquake. Disaster Adv. 6, 142–150.
Contents lists available at ScienceDirect
Tunnelling and Underground Space Technology journal homepage: www.elsevier.com/locate/tust
Control of fault lay-out on seismic design of large underground caverns Saeid Ardeshiri-Lajimi a, Mahmoud Yazdani a,⇑, Arya Assadi Langroudi b a b
Tarbiat Modares University, Civil and Environmental Engineering Department, Tehran, Iran Newcastle University, School of Civil Engineering and GeoSciences, England, United Kingdom
a r t i c l e
i n f o
Article history: Received 29 August 2014 Received in revised form 12 May 2015 Accepted 7 July 2015 Available online 13 August 2015 Keywords: Fault Underground Plastic zone Serviceability Seismic
a b s t r a c t Although buried structures are generally believed to suffer a lesser degree of damage in the event of earthquake – than that of over-ground structures – significant damage has been widely reported to buried assets after major earthquakes, including the 1995 Kobe and the 2008 Wen-Chuan. Discontinuity is one key feature of rock as the supporting medium around subsurface excavated spaces. Joints, faults and bedding planes influence, by-and-large, the stability of structures made from/into rock. In particular, fault system around underground caverns such as hydropower house has a marked control on assets’ seismic stability. This study builds on the current understanding through vigorous numerical modelling of fault-structure system under seismic excitation. A parametric approach is followed to determine the most critical layout of a single fault crossing a benchmark cavern. Fault system is systematically broken down into several combinations of dips and intersection points with cavern wall. For each case, a nonlinear dynamic analysis is conducted. To simulate the discontinuous medium, the hybrid finite difference – discrete element code CA2 (Continuum Analysis 2 dimensional) is implemented. The work showed that, similar to static conditions, fault influences the seismic stability of underground caverns through a tendency in extending the plastic zones and increasing displacements as well as asymmetric distribution of the latter and the former in rock medium. A 40–50° dip, single-point-intersection-on-crown k0 = 1 fault layout renders the most critical combination from both ultimate and serviceability limit states perspective. Under earthquake loading conditions however, the critical limit states condition took place for single fault intersected the cavern at heel and sidewall. The latter critical condition led to the tensile failure of cavern right sidewall. For faults intersecting the carven crown and having a k0 = 0.5, collapse would be more likely as fault dip increases. Collapse would be less likely with increasing dip for k0 = 0.5 fault crossing the bed and sidewall of caverns. Ó 2015 Elsevier Ltd. All rights reserved.
1. Introduction Common to nearly all reported damages to underground structures in the event of earthquake is the pronounced contribution of nearby faults to structural failure, highlighting the need for a better understanding of the interaction between rock discontinuities and structures they accommodate within (Zhang et al., 2013; Yashiro et al., 2007; Wang et al., 2001). Li (2012) emphasized on the significance of secondary fractures of seismogenic faults in proposing installation of reinforced concrete as secondary lining should tunnel alignment intersect such features. He recommended the lining to mantle the active fault zone and extend beyond, up to a minimum distance of 5 m. The micromechanics of fracture growth
⇑ Corresponding authors. E-mail addresses: [email protected] (S. Ardeshiri-Lajimi), myazdani@ modares.ac.ir (M. Yazdani), [email protected] (A. Assadi Langroudi). http://dx.doi.org/10.1016/j.tust.2015.07.002 0886-7798/Ó 2015 Elsevier Ltd. All rights reserved.
under high geo-energies and implications were addressed within a fractal context in a more recent work of Assadi et al. (2014). Ichimura et al. (2012) contributed to the understanding of explicit links between buried structures and faults through developing a finite element tool, which allows the measurement of seismic response of underground structures. Building on the latter platform, they demonstrated the reliability of fault-structure system analysis as a function of discontinuity geometry and orientation. Perhaps one of the earliest works in appreciation of faults as the origin of subsurface void space collapse is that of Brekke and Selmer-Olsen (1965), which discussed the matter from a geochemical perspective. They examined a number of failure incidents in Norwegian tunnels due to the smectite occurrences in between supporting rock discontinuities. From geomechanical standpoint however, fault-related instabilities include the collapse of large rock mass volumes within the fault zone into neighbouring caverns, particularly underground openings excavated through large thick fault zones. Instability is also likely after slip surface
306
S. Ardeshiri-Lajimi et al. / Tunnelling and Underground Space Technology 50 (2015) 305–316
formation along fault planes and consequent downfall of blocks or wedges – intersected by fault and other minor joints – into adjacent caverns (Nagelhout and Roest, 1997; Brekke and Howard, 1973). Other major contributions to association of instabilities and faults in rock mediums – a pivotal issue that have engaged both practical engineers and researchers – include Heap et al. (2010), Hao and Azzam (2005), Hao and Azzam (2001), Hashash et al. (2001), Tang and Kaiser (1998) and Tang (1997). Another stream of stability analysis of discontinuum fractured mediums concerns determination of joint properties from captured displacements, one latest being the wok of Yazdani et al. (2012). The developed knowledge however continues to lack a concrete theoretical solution to measurement of underground cavern stability. One approach to fill the latter shortfall is rigour benchmark modelling and cross-disciplinary researching, collating and analysing the generated big datasets in the framework of probabilistic analysis. Pursuing such approach shall positively impact the profession in enhancing sustainable and reliable support systems design for future buried infrastructures. Seismic response of underground structures is fundamentally different with that of surficial structures, in part due to their full enclosure in soil/rock medium and their often significant length dimension e.g. road tunnels. For buried structures, hence, to be designed fully resistant to seismic excitation, implementation of numerical modelling is not just common but an inevitable practice. Over the past years, various methods have been developed and extensively applied across the profession, yet, problem oriented numerical approaches are scarce. The relaxation of the latter is deemed here possible through the use of CA2 code (Fakhimi, 1998). The present work aims at evaluating the impacts of faults on induced plastic zones i.e. ultimate limit states as well as displacement re-distribution in subsurface openings i.e. serviceability limit states, in the event of earthquake. The main focus will be on measuring the effects of fault layout on the stability of caverns they cross – a function of fault dip and location. To simulate the benchmark discontinuity scenarios, the finite difference – discrete element code CA2 (for large static and dynamic strains) is used throughout, which particularly benefits in generation of tailored interface elements. The CA2 code is used to model the seismic excitation of a benchmark large underground cavern, seeking the critical position/dip of a single fault crossing the cavern section. That understanding will inform practitioners in taking tailored mitigation measures in lining design with respect to the host environment morphology.
2. Faults and modelling approaches Faults in definition are the commonly encountered large geological discontinuities in hard rock masses. Most severe underground instabilities are found to be closely associated with fault systems present nearby. Hao and Azzam (2001) in an investigation into the influence of fault on stability of large-scale underground caverns under static conditions demonstrated that faults may expand the plastic zones, give rise to displacements and distribute both displacements and plastic points in an asymmetrical fashion across the rock mass and around excavated spaces. The early types of joint elements proposed for simulating rock joints or faults in rock mechanics (Griffiths, 1985; Pande and Sharma, 1979; Goodman et al., 1968) are still widely implemented to model joints and interfaces across a broad range of geomechanical applications. Schweiger et al. (1990) employed the thin-layer joint element model to assay the fault thickness and stiffness effects on underground tunnel stability. Lei et al. (1995) used a contact-friction interface element to investigate the fault influence on the stress redistribution, bending moment and normal force across lining. Steindorfer (1998) utilised the ratio of longitudinal
to vertical displacement to assess the effect of fault zones in the process of tunnel excavation. More recent developments include the works of Bograd et al. (2011) and Mayer and Gaul (2007). Within the present work, joints are modelled as normal and shear springs with normal and tangential stiffness i.e. Kn and Ks, respectively between any pair of planes in contact (see Scheme 1). The model allows the explanation of principle fault characteristics including slipping, opening and compression behaviour. The model offers an appreciation of the fact that faults exhibit zero tension in the normal direction and follow Coulomb’s friction law for shear. The CA2 code employed is a two-dimensional explicit finite difference programme using time-marching procedures at very small steps through which the interface is divided into many differential elements. The incremental relative displacement vector at the contact point is resolved into the normal and shear directions and the total normal (Fn) and shear (Fs) forces are determinable using Eqs. (1) and (2). DtÞ F ðtþ ¼ F nðtÞ K n DU nðtþDt=2Þ L n DtÞ F ðtþ ¼ F sðtÞ K s DU sðtþDt=2Þ L s
ð1Þ
Eq. (1) is expressed in form of matrix as of,
DF n DF s
¼
Kn 0
0 DU n L DU s L Ks
ð2Þ
where the stiffness K n and K s equate stress over displacement, L is the length of interface element, DU n is the displacement increment of the interface element in the normal direction and DU s is the displacement increment of a segment of interface element in the tangential direction. The status of fault is expressed within the model by the following criteria: (a) Tensile yield condition:
F t ¼ rt t where rt is the tensile stress acting on the fault and t is the tensile strength of the fault. For Ft P 0, the interface breaks and the shear and normal forces are set to zero. The tensile strength of faults is generally agreed as of zero. (b) Shear yield condition: Shear stress along faults can be expressed using Mohr–Coulomb failure criterion as of:
F s max ¼ c l þ F n tan u where c is cohesion intercept along the interface, l is effective contact length and u is friction angle of interface. For |Fs| P Fs max, then Fs = Fs max. The rock is modelled as a bonded particle system. The rock grains are assumed as circular cylinders that interact through normal and shear springs. Each discrete element in CA2 is a circular cylinder. The cylinder–cylinder and cylinder–wall (i.e. interface) interactions are modelled by normal and shear springs, normal
Scheme 1. Fault model. The interface represented by side A and side B, which are connected by normal (kn) and shear (ks) stiffness springs. (S = slider, T = tensile strength and LN = length associate with grid-point N).
S. Ardeshiri-Lajimi et al. / Tunnelling and Underground Space Technology 50 (2015) 305–316
and shear bonds, and a friction coefficient. Therefore, the micromechanical constants for cylinder–cylinder or cylinder–wall interactions are the normal and shear spring constants (kn and ks), the normal and shear bonds (nb and sb), and the friction coefficient (l). This simple contact bond model behaves elastically under the applied stresses until the tensile normal force or shear force at a contact point exceeds the normal or shear bond leading to tensile or shear failure at that contact point. A failed contact point loses its bond permanently but still it can carry shear forces through a Coulomb law using the friction coefficient l if it is subjected to a compressive load – also see Fakhimi and Gharahbagh (2009) and Hosseinpour (2008). 3. Modelling and boundary conditions To assess the influence of fault system on seismic stability of underground openings, a benchmark model is adopted – see Fig. 1. The underground opening considered in this model represents many forms of hydroelectric schemes such as underground power houses and surge chambers. As such, adopted model’s cross-section has a crown of arch in shape and sidewalls sizing 20 m wide by 30 m high. To minimize the impact of uninteresting boundary conditions on modelling outputs, the size of model in each direction is adopted eight times the cavern dimension. The fault is extended far enough from the excavation and the fault displacements dissipate at model borders, replicating the remote rock masses. Under static conditions, vertical boundaries are fixed along horizontal direction and zero displacement condition is applied across the base of the model. The effect of overlying (i.e. overhead) rock weight is considered on the model top boundary by a normal stress. Adopting displacement
307
boundary condition instead of stress boundary condition on the vertical boundaries of the model is the premise of obtaining ideal numerical results for the fault-related problem (Hao and Azzam, 2001). A fault with different dips crosses through the rock mass in the vicinity of the cavern i.e. intersects the cavern. The fault strike is presumed parallel or nearly parallel to the alignment of the cavern and hence a plane strain state governs. The latter is well established as the most unfavourable condition in the context of underground caverns stability (Hoek and Brown, 1982). 3.1. Model considerations under static loads Throughout the analysis, the underground cavern is assumed to be at a depth of 300 m below the ground level. The in-situ stress state prior to cavern excavation is sourced from the weight of rock masses at rest. The vertical stress varies linearly with depth while the horizontal stress is calculated as of rh = k0rv, where k0 is the ratio between the horizontal and vertical stresses. Intact rock is considered as an elastic–perfectly plastic material that follows Mohr–Coulomb failure criterion. Slip failure of the fault also follows the Mohr–Coulomb criterion. Input parameters for the analysis purposes are determined in a way to satisfy following requirements: The properties of rock mass and fault need to represent the most unfavourable fault influence, while reflecting real characteristics of the materials. As such, the elastic modulus of the rock mass is taken as of 20 GPa, the internal friction angle is fixed at c. 35°, while small values for tensile strength and cohesion are adopted. These specified geotechnical parameters represent sparse jointed rock masses with low-to-medium strength in nature. Since the emphasis of this analysis is to discuss the effects of key fault parameters, only fault stiffness, tensile strength and cohesion are needed to be specified. The intention is to study a typical weak fault with a weathered, persistent and smooth surface in absence of tensile strength and cohesion and thin thickness. Thereby, shear dilation, strain hardening and softening of fault are neglected in the model. The latter parameters have no control on the general discussion of the problem throughout the analysis. Within the scope of current work, a comparative analysis of the effect of different fault parameters is followed throughout the coming sections, leaving the absolute prediction of a practical problem out of the interest of the present article. The adopted geotechnical parameters for the rock mass and fault render a stable model under static conditions. 3.2. Model considerations under dynamic loads The relation between the elastic modulus in dynamic loading environments, Ed, and the elastic modulus under static loading circumstances, Es, is given by Hayashi and colleagues in 1973, reading:
Ed ¼ ð1:3 1:7ÞEs
ð3Þ
The magnitude of Ed is obtained through field full-scale dynamic loading test, with the Poisson’s ratio md ranging from 0.2 to 0.3. According to Price et al. (1969), the tensile strength increases up to 1.5 times that of the static value in the event of quick loading. The mechanical properties of rock mass and fault used in the present research under static and dynamic conditions are summarized in Table 1. Fig. 3 shows the baseline-corrected response spectra graph i.e. input motion. 3.3. Dynamic analysis
Fig. 1. Sketch of the model for this numerical analysis.
Following the static analysis under the static effective stress conditions, below steps are adopted to undertake the dynamic analysis:
308
S. Ardeshiri-Lajimi et al. / Tunnelling and Underground Space Technology 50 (2015) 305–316
Table 1 Input parameters of rock mass and fault for the numerical modelling. Material properties 3
Density q (kg/m ) Internal friction angle / Cohesion c (MPa) Poisson ratio Normal stiffness kn (GPa/m) Shear stiffness ks (GPa/m) Es (GPa) Ed (GPa) Tensile strength rt (kPa)
Rock mass
Fault
2600 35 2 0.2 –
24 0 – 20 2 – – –
20 25 1.7
Fig. 2. CA2 model boundary conditions for dynamic simulation.
Step 1: In order to use horizontal component of the acceleration record as input event, it is essential to accomplish baseline correction on above-mentioned histories first. As a result, after integrating the acceleration histories twice, the obtained displacement does not reach zero value at the end of the motion. One approach for correcting this error is to add a low frequency wave to earthquake loading so that displacement reaches zero in the end. Step 2: Another issue in dynamic loading which is of considerable importance is the wave propagation in the model. Generally, in performing dynamic analyses the presence of error in the form of wave propagation is likely to happen. The finite difference grid dimensions were selected by taking into account the maximum frequency (f) of the shear wave that the model could logically respond to during earthquake loading. Consequently, this prevents the incorrect propagation of the waves. The frequency is determined by the following equation (Kuchlemeyer and Lysmer, 1973):
f ¼
Cs 10Dl
ð4Þ
where C s is the shear wave velocity of the soil and Dl is the largest grid zone size in the model. Referring to equation, a uniform zone size of 1 m 1 m was selected. Since the lowest shear wave velocity in the model belongs to the soil deposits, the highest admissible frequency for a propagating shear wave is 5 Hz. Therefore, the input earthquake record shall be filtered by a low pass filter to remove frequency components higher than 5 Hz. A frequency of 4 Hz was ultimately selected as the low pass filter for the reduction in shear wave velocity which may occur due to plastic flow during seismic loading. Step 3: Application of dynamic boundary conditions: The boundary conditions used in the dynamic analysis are illustrated in Fig. 2. Quiet boundaries were used on all outside boundaries of model. These boundaries prevent reflection of outgoing seismic waves back into the model. Quiet boundaries were combined with free-field boundaries on the vertical outside boundaries that prevent distortion of vertically propagating plane waves along the boundaries. Dynamic loading was applied at the bottom of model, as propagating vertically upwards (Pakbaz and Yareevand, 2005). The record of accelerogram of 1999 Kocaeli, Turkey earthquake with magnitude of 7.4 in scale of Richter (Fig. 4) has been used. Step 4: Rayleigh damping, which consists of two viscous elements, is conventionally used in the numerical analyses herein. The two elements of Rayleigh damping are both frequency dependent (Lysmer and Kuchlemeyer, 1969). By choosing a mid-frequency at which the combined effects of the two elements cancel out, it is possible to have a damping that is nearly independent of frequency over a fairly wide range of frequencies on either side of the mid-frequency (White et al., 1977). The mid-frequency is usually chosen in the range between the natural frequency of the model and the predominant frequency of the input motion. Rayleigh damping was assigned to each element of the model in the mid-frequency. In analyses that use one of the plasticity constitutive models (e.g. Mohr– Coulomb), a considerable amount of energy dissipation can occur during plastic flow (viscous damping). Thus, for many dynamic analyses that involve large-strain, only a minimal percentage of damping (e.g. 5%) may be required (Itasca, 2002). Therefore, the damping ratio was assumed to be 5% in the analyses. Dynamic analyses were performed and the results were extracted for interpretation and further assessment. 4. Numerical results and discussion The effect of earthquake on underground opening depends on various parameters including peak acceleration, magnitude and duration of earthquake, the relative rigidity between underground opening and ground and in-situ stress state. Furthermore, for
Fig. 3. Acceleration record of kocaeli earthquake.
S. Ardeshiri-Lajimi et al. / Tunnelling and Underground Space Technology 50 (2015) 305–316
309
underground caverns crossed by a single fault, additional parameters including intersection of fault with cross section and mechanical properties of fault have control on the stability. In this section the effects of fault orientation and the relative location of its intersection with cavern perimeter together with the PGA (Peak Ground Accelerations) of earthquake loading are examined on seismic stability of cavern. To quantify the stability analysis, two indices, namely the extended plastic area in rock mass (analogous to ultimate limit state) and the maximum displacement in tunnel wall or roof (analogous to serviceability limit state) have been adopted. Quantified outputs are obtained for a combination of fault dips, PGAs and horizontal to vertical stress ratio (k0). Discussions are lastly built on comparisons between displacements and area of plastic zones before and after application of earthquake loading is undertaken.
Fig. 6. Variation of maximum vertical displacement of roof versus fault dip for different PGA in case 3: fault intersects cavern floor in the middle point and also right cavern wall (k0 = 1; u = 24°).
Fig. 4. Variation of maximum vertical displacement of roof versus fault dip for different PGA in case 1: fault intersects cavern roof (k0 = 1; u = 24°).
Fig. 7. Variation of maximum vertical displacement of roof versus fault dip for different PGA in case 4: fault intersects cavern floor in the corner point (k0 = 1; u = 24°).
4.1. Maximum vertical displacement
Fig. 5. Variation of maximum vertical displacement of roof versus fault dip for different PGA in case 2: fault intersects left cavern wall in the middle point and also the roof (k0 = 1; u = 24°).
Figs. 4–7 illustrate the variation of maximum rock mass vertical displacement against fault dips for varying PGAs and a suite of faults crossing the cavern periphery. A k0 = 1 and fault frictional angle of u = 24° was considered in building the models. The presented values represent displacements induced by earthquake loading at its ultimate excitation and are independent from initial displacement originated from in-situ stresses. The maximum vertical displacement took place at cavern roof, lending evidence to the remarkable impact of fault when intersecting the cavern at a single point on crown. In general, maximum roof settlement occurred for fault dips ranging from 40° to 50°. Other key observations include: (1) Vertical displacement on crown appeared more sensitive to PGA for accelerations greater than 0.15 g. (2) Dip-displacement
310
S. Ardeshiri-Lajimi et al. / Tunnelling and Underground Space Technology 50 (2015) 305–316
Fig. 8. Relationship of the area of rock mass plastic zones varying with fault dips and horizontal stress ratio for case 1, PGA = 0.15 g & friction angle of fault u = 30°.
Fig. 10. Variation of plastic area extent versus fault dip for different PGA in case 1: fault intersects cavern roof (k0 = 1; u = 24°).
Fig. 9. Relationship of the area of rock mass plastic zones varying with fault dips and horizontal stress ratio for case 3, PGA = 0.15 g & friction angle of fault u = 30°.
Fig. 11. Variation of plastic area extent versus fault dip for different PGA in case 2: fault intersects left cavern wall in the middle point and also the roof (k0 = 1; u = 24°).
trend tool a bi-modal distribution when fault intersected the cavern at crown and side wall both. (3) Predicted maximum displacement was 15–20 times greater when fault intersected the crown at a single point rather than any other relative position. Examining the vertical displacement values in Figs. 4–7 infers a 40–50° dip, single-point-intersection-on-crown fault layout renders the most critical combination in terms of serviceability.
the critical dip. Thereby, a 0.15g PGA was adopted throughout this course of analysis. Original (prior to excavation) vertical in-situ stress was assumed to equate the weight of the rock mass below the ground level – which varies linearly with depth. Horizontal in-situ stress was taken as rh = k0 rv, where k0 is the ratio between horizontal stress and vertical stress. k0 was varied at constant vertical stress (i.e. constant buried depth of the benchmark cavern) for every analysis case. Analysis inferred that the critical dip values discussed in Section 4.1 are strongly a function of k0. For Case 1 (where the fault system intersects the crown at single point) plastic zone area gained remarkably high values for k0 = 0.5 and followed an upsurge trend with dip, dissimilar to the unimodal pattern captured for k0 = 1. As such, the conclusion gathered in Section 4.1 may be read as: a 40–50° dip, single-point-intersection-on-crown k0 = 1 fault
4.2. In-situ stress state The control of fault on the plastic zones surrounding the cavern for varying in-situ stress state was investigated through parametric analysis (see Fig. 8 and 9). Overall, PGA was found directly proportional to the area of plastic zone and displacement magnitude, with no potential control on
S. Ardeshiri-Lajimi et al. / Tunnelling and Underground Space Technology 50 (2015) 305–316
311
Fig. 12. Variation of plastic area extent versus fault dip for different PGA in case 3: fault intersects cavern floor in the middle point and also right cavern wall (k0 = 1; u = 24°).
Fig. 14. Relationship of the area of rock mass plastic zones varying with fault dips and fault location, (a) under static loads and (b) under dynamic loads for PGA = 0.15 g, horizontal stress ratio k0 = 1 & friction angle of fault u = 24°. Fig. 13. Variation of plastic area extent versus fault dip for different PGA in case 4: fault intersects cavern floor in the corner point (k0 = 1; u = 24°).
layout renders the most critical combination from a serviceability perspective. For faults intersecting the carven crown and having a k0 = 0.5, collapse would be more likely as fault dip increases. Collapse would be less likely with increasing dip for k0 = 0.5 fault crossing the bed and sidewall of caverns.
reworded as: a 40–50° dip, single-point-intersection-on-crown k0 = 1 fault layout renders the most critical combination from both ultimate and serviceability limit states perspective. Critical dip was measured as of 50° for when fault intersects the cavern at the side wall and in the middle point (Figs. 11 and 12), and 30° for when fault intersects the cavern at the lower floor corner (Fig. 13). 4.4. Effect of fault relative location (position)
4.3. Area of plastic zone For four fault locations, the area of plastic zone – surrounding cavern was plotted against fault dip in Figs. 10–13. Plots were re-produced for four PGAs. Plastic area appeared a function of fault dip, the maximum content of which being strongly a factor of fault location. For faults crossing cavern crown (established as the critical position on Sections 4.1 and 4.2), plastic zone gently expands with fault dip, although the maximum plastic area occurs at a critical 40° dip. To this end, the main line of conclusion in Section 4.2 may be
As depicted in Fig. 1, four possible fault sites – with respect to the cavern – were considered in the current work. To facilitate the discussion and better coherency with the previous arguments, Case 1 is defined for when the fault crosses (is tangent to) the cavern roof. In Case 2, the fault intersects the mid-height of the left sidewall. In Case 3, the fault cuts the cavern crown and through the centre of the bottom surface. In Case 4, the fault emerges out close to the heel of the cavern’s sidewall. For all the cases, the fault dips are assumed within the range of 20–80° for consistent mechanical properties of rock mass and fault.
312
S. Ardeshiri-Lajimi et al. / Tunnelling and Underground Space Technology 50 (2015) 305–316
Fig. 14 plots the plastic area in the rock mass against fault dip for a suite of fault locations. An interesting drawing from the results is the resembling pattern of plastic area – fault dips of Case 1 and Case 3, with Case 2 and Case 4, respectively, where maximum plastic area is attained as the fault dips approach 40° and 60°, respectively. Under seismic excitation and from ultimate limit state perspective (i.e. area of plastic zone), the critical condition took place for single fault intersected the cavern at heel and sidewall. This is dissimilar with the static conditions in which critical conditions were associated with the single fault crossing the cavern’s crown. 4.5. Maximum shear displacement of fault Fig. 15 shows the variation of maximum shear displacement of fault against dip. Independent of fault location, maximum shear displacement took place at 50° dip under all PGAs. 4.6. Comparing static and dynamic results
Fig. 15. Maximum shear displacement of fault.
Maximum values of vertical displacements on cavern roof, horizontal displacements on cavern walls and area of plastic zone under dynamic and static condition are compared in pairs and discussed in Figs. 16–18:
Fig. 16. Comparing static and dynamic values of maximum vertical displacements of cavern roof, for PGA = 0.15 g horizontal stress ratio k0 = 1 & friction angle of fault u = 24°.
S. Ardeshiri-Lajimi et al. / Tunnelling and Underground Space Technology 50 (2015) 305–316
4.6.1. Maximum vertical displacement Fig. 16 shows the variation of maximum rock mass vertical displacement against fault dip under static and dynamic condition for PGA = 0.15 g. A k0 = 1 and u = 24° fault frictional angle was adopted throughout. The annotation ‘static’ on plot represents the displacement values induced by in-situ stresses after excavations, while ‘dynamic’ tag represents displacements induced after earthquake loading at its final moment of excitation, independent from (excluding) the initial displacement produced by in-situ stresses. The annotation ‘sum’ stands for the total of static and dynamic values. Following drawings were made through examining the plots: (1) Vertical displacement of cavern crown is less dependent on dip under static loading conditions than dynamic excitation. (2) Critical dip, in terms of crown vertical displacement i.e. serviceability limit state sits in the 40° to 50° interval, predominantly controlled by strains in dynamic environment. (3) Maximum vertical displacement is likely to take place for fault crossing the top of the cavern. For other fault locations, the maximum value falls at least 10 times shorter than the latter, reinforcing the need for particular mitigation measures to be taken for excavation immediately beneath fault lines. (4) No link between
313
the trend of maximum vertical displacement under static and dynamic conditions could be identified. 4.6.2. Maximum horizontal displacement Figs. 17 and 18 illustrate the variation of maximum horizontal displacements of cavern left and right sidewall against fault dip under static and dynamic (PGA = 0.15 g) conditions. A k0 = 1 and u = 24° fault frictional angle was adopted throughout. Plots suggest no correlation between static and dynamic values. Due to the tensile failure of cavern right sidewall on Case 4, horizontal displacement-dip plot could not possibly be generated. Key observations include: (1) For fault location intersecting the left side-wall, horizontal displacement on the left-wall under static loading appeared to be strongly a function of dip. In contrary, left wall horizontal displacement followed a fairly plateau trend against dip for faults with no intersection with the left side-wall. Similar pattern follows for the right side-wall. (2) No link could be established between the fault location (i.e. Case) and the occurrence of maximum horizontal displacement. In Fig. 17, maximum horizontal displacement occurred in Case 1 when fault crossed the cavern crown, while in Fig. 18, maximum horizontal
Fig. 17. Comparing static and dynamic values of maximum horizontal displacements of cavern left wall, for PGA = 0.15 g horizontal stress ratio k0 = 1 & friction angle of fault u = 24°.
314
S. Ardeshiri-Lajimi et al. / Tunnelling and Underground Space Technology 50 (2015) 305–316
displacement occurred in Case 3 for fault crossing the mid-height of the right side-wall.
4.7. Area of plastic zone Fig. 19 illustrates the variation of plastic area within the rock mass around the excavated space against fault dips for static and dynamic conditions and for four fault positions (intersection of fault with cavern periphery). Plots suggest the existence of a link between the fault-induced area of rock mass plastic zone under earthquake loading and static loading. Whilst the plots represent a PGA of 0.15 g for the dynamic analysis, same trends were re-produced for PGAs of 0.1 g, 0.20 g and 0.25 g. This was confirmed further in Ardeshiri (2008), where input PGA values were extended to six. Eq. (5), illustrates the established relationship between area of plastic zones before and after earthquake loading. The equation was successfully validated for two other earthquake waves i.e. Loma Prieta and Duze, with different predominate period values (Ardeshiri, 2008).
ðAPa Þd ¼ ð1 þ 0:25e10:92amax Þ ðAPa Þs
ð5Þ
where ðApa Þd is the area of plastic zones after earthquake loading, amax is maximum earthquake acceleration and ðApa Þs is area of plastic zones under in-situ stress condition.
5. Conclusion
Fig. 18. Comparing static and dynamic values of maximum horizontal displacements of cavern right wall, for PGA = 0.15 g horizontal stress ratio k0 = 1 & friction angle of fault u = 24°.
In this paper, CA2 code was used to revisit, in a systematic manner, the effects of fault lay-out on the dynamic and static response of a benchmark cavern excavated into a rock mass. Numerical analysis indicates that the influence of a fault to the underground cavern stability is different from its inherent properties. Similar to static conditions, fault influences the seismic stability of large underground caverns through a tendency in extending the plastic zones and increasing displacements as well as asymmetric distribution of the latter and the former in rock medium. This however is strongly a function of fault dip and location (i.e. lay-out), together with the situ stress state. Fault dip in relation to the underground cavern is proved to have key control on cavern’s response to earthquake forces. The present work confirmed a nearby 40° and 50° fault dip to be most critical. For a suite of k0 and the typical configuration of cavern adopted in the present work, an earlier expectation was met in the attained decreasing area of rock mass plastic zone with an increase in the horizontal stress ratio. The significance of fault in seismic stability of underground cavern was further highlighted in the high dependency of serviceability e.g. vertical displacement of cavern crown on dip under dynamic excitation rather than static stresses. When fault system crosses the cavern section, sharp edges could form on cavern section. Faced to earthquake waves, these sharp edges can generate a state of tensile failure. Given the local nature of such damage on cavern section however, sufficient flexibility could provide the required resistance of cavern against seismic loads. An equation for determining the area of the earthquake-induced plastic zones was proposed, which informs engineers in improving the design of support system in seismic regions. The work can be extended by interested workers in bringing into parametric analysis other key parameter, in addition to fault parameters, including PGV (Peak Ground Velocity), frequency content of earthquake wave, depth of excavated cavern.
S. Ardeshiri-Lajimi et al. / Tunnelling and Underground Space Technology 50 (2015) 305–316
315
Fig. 19. Comparing static and dynamic area of plastic zones values of cavern, for PGA = 0.15 g horizontal stress ratio k0 = 1 & friction angle of fault u = 24°.
References Ardeshiri, S., 2008. Influence of Discontinuities on Seismic Behaviour of Large Underground Caverns. Master of Science Tarbiat Modares Univfersity – TMU. Assadi, A., Jefferson, I., O’Hara-dhand, K., Smalley, I., 2014. Micromechanics of quartz sand breakage in a fractal context. Geomorphology 211, 1–10. Bograd, S., Reuss, P., Schmidt, A., Gaul, L., Mayer, M., 2011. Modeling the dynamics of mechanical joints. Mech. Syst. Signal Process. 25, 2801–2826. Brekke, T.L., Howard, T.R., 1973. Functional Classification of Gouge Materials from Seams and Faults in Relation to Stability Problems in Underground Openings. US Bureau of Mines (AROA) – Con. No. H0220022. Brekke, T.L., Selmer-Olsen, R., 1965. Stability problems in underground constructions caused by montmorillonite-carrying joints and faults. Eng. Geol. 1, 3–19. Fakhimi, A., 1998. Theory and User Manual of CA2 Computer Programe – Report no. 262. Building and Housing Research Centre, Tehran, Iran. Fakhimi, A., Gharahbagh, E.A,. 2009. Discrete element modeling of the influence of void size and distribution on the mechanical behavior of rock. In: Diederichs, M., Grasselli, G., (Eds.), Proceedings of the 3rd CANUS Rock Mechanics Symposium. Toronto, Canada. Goodman, R.E., Taylor, R.L., Brekke, T.L., 1968. A model for the mechanics of jointed rock. J. Soil Mech. Foundation Division, ASCE 94, 637–659. Griffiths, D.V., 1985. Numerical modelling of interfaces using conventional finite elements. In: Kawamoto, T., Ichikawa, Y., (Eds.), Proceedings of the 5th International Conference on Numerical Methods in Geomechanics, Nagoya, Japan. pp. 837–844. Hao, L.H., Azzam, R.,2001. Analysis of fault effects on the stability of large underground caverns. In: ISRM International Symposium – 2nd Asian Rock Mechanics Symposium. Beijing, China.
Hao, Y.H., Azzam, R., 2005. The plastic zones and displacements around underground openings in rock masses containing a fault. Tunnell. Underground Space Technol. 20, 49–61. Hashash, Y.M.A., Hook, J.J., Schmidt, B., I-Chiang Yao, J., 2001. Seismic design and analysis of underground structures. Tunnell. Underground Space Technol. 16, 247–293. Heap, M.J., Faulkner, D.R., Meredith, P.G., Vinciguerra, S., 2010. Elastic moduli evolution and accompanying stress changes with increasing crack damage: implications for stress changes around fault zones and volcanoes during deformation. Geophys. J. Int. 183, 225–236. Hoek, E., Brown, E.T., 1982. Underground Excavations in Rock. Institute of Mining and Metallurgy, London. Hosseinpour, H., 2008. Numerical and experimental evaluation of the effect of an oversize particle on direct shear test results. New Mexico Bureau Geol. Min. Resour. Ichimura, T., Hori, M., Quinay, P.E., Wijerathne, M.L.L., Suzuki, T., Noguchi, S., 2012. Comprehensive numerical analysis of fault-structure systems – computation of the large-scale seismic structural response to a given earthquake scenario. Earthquake Eng. Struct. Dyn. 41, 795–811. Itasca, F., 2002. Fast Lagrangian Analysis of Continua – Version 4.0 user guide. Thrasher Square East, 708: Itasca Consulting Group Inc. Kuchlemeyer, R.L., Lysmer, J., 1973. Finite element method accuracy for wave propagation problems. J. Soil Mech. Found. Div., ASCE 99, 421–427. Lei, X.Y., Swoboda, G., Zenz, G., 1995. Application of contact – friction interface element to tunnel excavation in faulted rock. J. Comput. Geotech. 17, 349–370. Li, T., 2012. Damage to mountain tunnels related to the Wenchuan earthquake and some suggestions for aseismic tunnel construction. Bull. Eng. Geol. Environ. 71, 297–308. Lysmer, J., Kuchlemeyer, R.L., 1969. Finite dynamic model for infinite media. J. Eng. Mech. Div., ASCE 95, 859–877.
316
S. Ardeshiri-Lajimi et al. / Tunnelling and Underground Space Technology 50 (2015) 305–316
Mayer, M.H., Gaul, L., 2007. Segment-to-segment contact elements for modelling joint interfaces in finite element analysis. Mech. Syst. Signal Process. 21, 724– 734. Nagelhout, A.C.G., Roest, J.P.A., 1997. Investigating fault slip in a model of an underground gas storage facility. Int. J. Rock Mech. Min. Sci. 34, 212. Pakbaz, M.C., Yareevand, A., 2005. 2-D analysis of circular tunnel against earthquake loading. Tunnell. Underground Space Technol. 20, 411–417. Pande, G.N., Sharma, K.G., 1979. On joint interface elements and associated problems of numerical ill-conditioning. Int. J. Numer. Analyt. Methods Geomech. 3. Price, D.G., Malkin, A.B., Knill, K.J., 1969. Foundations of multi-storey blocks on coal measures with special reference to old workings. Q. J. Eng. Geol. 1, 271–322. Schweiger, H.F., Haas, W., Handle, E., 1990. A thin layer element for modelling joints and faults. In: Rossmanith, H.P., (Ed.), Proceedings of International Conference on Jointed and Faulted Rock, 1990 Vienna, Austria. pp. 559–564. Steindorfer, A., 1998. Short term prediction of rock mass behaviour in tunelling by advanced analysis of displacement monitoring data. In: Schuvert. Riedmuller, Semprich (Eds.), Schriftenreihe der Gruppe Geotechnik, Heft 1. Tang, C., 1997. Numerical simulation of progressive rock failure and associated seismicity. Int. J. Rock Mech. Min. Sci. 34, 249–261.
Tang, C.A., Kaiser, P.K., 1998. Numerical simulation of cumulative damage and seismic energy release during brittle rock failure—part I: fundamentals. Int. J. Rock Mech. Min. Sci. 35, 113–121. Wang, W.L., Wang, T.T., Su, J.J., Lin, C.H., Seng, C.R., Huang, T.H., 2001. Assessment of damage in mountain tunnels due to the Taiwan Chi–Chi Earthquake. Tunnell. Underground Space Technol. 16, 133–150. White, W., Lee, I.K., Valliappan, S., 1977. Unified boundary for finite dynamic models. J. Eng. Mech. Div., ASCE 103, 949–964. Yashiro, K., Kojima, Y., Shimizu, M., 2007. Earthquake damage to tunnels in Japan and case studies of railway tunnels in the 2004 Niigataken-Chuetsu earthquake. Quarterly Report Railway Tech. Res. Inst. 48, 136–141. Yazdani, M., Sharifzadeh, M., Kamrani, K., Ghorbani, M., 2012. Displacement-based numerical back analysis for estimation of rock mass parameters in Siah Bisheh powerhouse cavern using continuum and discontinuum approach. Tunnell. Underground Space Technol. 28, 41–48. Zhang, J.W., Mei, Z.R., Quan, X.J., 2013. Failure characteristics and influencing factors of highway tunnels damage due to the earthquake. Disaster Adv. 6, 142–150.