Mechanical response of a tunnel subjected to strike-slip faulting processes, based on a multi-scale modeling method

Abstract

The stick-slip action of strike-slip faults poses a significant threat to the safety and stability of underground structures. In this study, the north-east area of the Longmenshan fault, Sichuan, provides the geological background; the rheological characteristics of the crustal lithosphere and the nonlinear interactions between plates are described by Burger’s viscoelastic constitutive model and the friction constitutive model, respectively. A large-scale global numerical model for plate squeezing analysis is established, and the seemingly periodic stick-slip action of faults at different crust depths is simulated. For a second model at a smaller scale, a local finite element model (sub-model), the time history of displacement at a ground level location on the Longmenshan fault plane in a stick-slip action is considered as the displacement loading. The integration of these models, creating a multi-scale modeling method, is used to evaluate the crack propagation and mechanical response of a tunnel subjected to strike-slip faulting. The determinations of the recurrence interval of stick-slip action and the cracking characteristics of the tunnel are in substantial agreement with the previous field investigation and experimental results, validating the multi-scale modeling method. It can be concluded that, regardless of stratum stiffness, initial cracks first occur at the inverted arch of the tunnel in the footwall, on the squeezed side under strike-slip faulting. The smaller the stratum stiffness is, the smaller the included angle between the crack expansion and longitudinal direction of the tunnel, and the more extensive the crack expansion range. For the tunnel in a high stiffness stratum, both shear and bending failures occur on the lining under strike-slip faulting, while for that in the low stiffness stratum, only bending failure occurs on the lining.

Graphical abstract

Keywords

stick-slip action / plate squeezing analysis / multi-scale modeling method / lining cracking / mechanical response

1 Introduction

Sudden brittle fractures or stick-slip actions between plates in the earth’s crust release large amounts of energy and cause earthquakes. After the San Francisco earthquake (1906) of magnitude M_w 7.7, the Wrights tunnel, which crosses the San Andreas fault, was severely damaged and offset by 1.5 m [1]. The Wenchuan earthquake (2008) of magnitude M_w 7.9 generated a surface rupture over 200 km long with a maximum dislocation of 8 to 10 m [2]. The permanent deformation of the stratum caused by fault stick-slip action usually leads to the cracking and destruction of tunnels and can even result in the collapse of entire structures. Therefore, it is important to investigate, in the design and construction stages, the mechanical response of tunnels subjected to strike-slip faulting processes, thus enabling rapid disaster response and efficient traffic control.

Strike-slip faults impose both longitudinal and lateral forces on tunnel structures, severely threatening tunnel structural integrity [3]. To ensure the safe operation of tunnels, numerous related studies, including theoretical analyses, numerical simulations, and laboratory model tests, have been conducted [4–7]. Newmark and Hall [8] proposed an analytical method by adopting a small-displacement model with static soil pressure and static friction force. Based on the Pasternak beam model, Qiao et al. [9] proposed an analytical model for the longitudinal mechanical analysis of tunnels by considering the effects of fault dislocations. Scale model and centrifuge tests have also been conducted using sand and clay samples to determine the propagation and geometric characteristics of fractures [10,11]. Zhou et al. [12] conducted a series of model tests and numerical simulations to investigate the deformation and failure mechanism of an underwater tunnel on the seabed with flexible joints subjected to strike-slip faulting. Qiao et al. [13] analyzed the mechanism of the longitudinal response of a tunnel under normal faulting based on experimental data and numerical simulations. A three-dimensional (3D) numerical discrete-continuum coupling approach was also adopted by Ma et al. [14] to investigate the deformation and failure mechanism of tunnels subjected to normal faulting processes. However, most studies use a constant velocity of displacement or a fixed displacement to simulate the stick-slip action of faults, without adequately considering the spatial and dynamic effects of the faulting in the regional geological environments, resulting in tunnel stress and strain values that were not appropriately simulated.

Large-scale structural analysis of long tunnels is difficult due to limitations of computing power, especially for tunnels embedded in complex regional geological environments. When crossing active strike-slip faults, the stick-slip action induced by plate squeezing is complicated by the strong nonlinear interactions between the hanging wall and footwall of faults, and the process is more likely to cause serious damage to tunnel structures [15]. In addition, if the size difference of adjacent elements is too large in numerical simulations, then a large calculation error is inevitable. The multi-scale modeling method is an advanced analysis approach in finite element methods, so it can be considered to be a potential tool for solving the above key issues, and has been widely used in tunnel engineering. Liu and Wang [16] conducted a three-dimensional numerical analysis of an underground bifurcated tunnel based on multi-scale modeling technology to analyze the branching location, corresponding stress and deformation results in practice. Cao et al. [17] adopted a multi-scale modeling method to simulate water hammer impacts in a 14 km long pressurized water conveyance tunnel, which not only successfully simulated water hammer effects along the full tunnel length, but also obtained the detailed mechanical responses of the tunnel lining. Sun et al. [18] performed weld strength analysis of T-Joint segments of a metro tunnel based on the multi-scale modeling method, establishing a reliable and advanced simulation process. Zhang et al. [19] explored the construction of a 3D geological model near fault zones, based on a multi-scale numerical analysis of in situ stress field of long tunnels. However, few studies have been carried out on the mechanical response analysis of tunnels crossing active faults using a multi-scale modeling method and considering the real regional geological environments and complicated characteristics of stick-slip action.

This study selected the north-east area of the Longmenshan fault as the geological background, considering the rheological characteristics of the crustal lithosphere and the nonlinear interactions between plates, to establish a large-scale plate squeezing numerical model for simulating the stick-slip action of strike-slip faults. On this basis, the displacement time history of the fault plane during the stick-slip action was treated as the displacement loading in the local finite element model to evaluate the crack propagation and mechanical response of a tunnel that crosses strike-slip faults. The proposed multi-scale modeling method and the extended finite element theory were used to analyze the cracking process and mechanical response of a tunnel crossing strike-slip faults. The numerical results obtained may provide a reference for relevant studies on the damage and deformation analyses of similar tunnels under strike-slip faulting.

2 Methodology

2.1 Multi-scale modeling

The lengths of tunnels are frequently several kilometres, and some of those tunnels inevitably cross active faults. Their regional geological environments are complicated, especially at the plate boundaries. To simulate the stick-slip action of faults during plates squeezing, a large-scale global model is necessary. Such a model was discretized, with a relatively coarse mesh in the area away from faults and a relatively fine mesh near the faults, and detailed structures of the tunnel were ignored. A mesh sensitivity analysis was conducted to determine a suitable mesh dimension that could balance computational time with the accuracy of results.

A tunnel lining crossing strike-slip faults is associated with various detailed data, including the tunnel size, lining thickness, fault inclination and width, and stratum type. The details of mechanical response cannot be appropriately simulated by using a large-scale global model. A simulation needs to be discretized with a relative fine mesh; a small size of elements should be adopted. As a result, a large-scale tunnel model inevitably leads to an oversized finite element model, that is difficult to apply even on supercomputers.

The scale difference between tunnels and plates is extremely large, causing the computational efficiency of the numerical analysis of the tunnel–plate system to be extremely poor. Therefore, multi-scale modeling, not only considering large-scale, should be used, and there is an inevitable tradeoff between accuracy and computational cost. This work defines multi-scale modeling as the combination of a large-scale global numerical model of plate squeezing and a local sub-model of a tunnel subjected to strike-slip faulting, and the proposed method can decrease the complexity and computational expense of numerical simulations.

2.2 Simulation workflow

Structural analysis of tunnels located in active faults requires a number of specialized tools. In the framework adopted in Refs. [20,21], the stick-slip actions of faults were simulated in the large-scale global numerical model, and relevant key parameters (acceleration, velocity and displacement) of the fault in the stick-slip action were calculated. Subsequently, the displacement time history near the ground on the Longmenshan fault plane obtained during the stick-slip action was taken as the displacement loading in the local finite element model, and the dynamic faulting analysis was solved by a sub-model method based on extended finite element theory. Fig.1 provides a detailed view of the simulation workflow.

3 Large-scale model for plate squeezing analysis

3.1 Configurations and mesh

The numerical simulation of plate squeezing requires the development of a mathematical model of a portion of the earth’s crust with a width of hundreds of kilometres; the depth of the crust is relatively small relative to the width of the plates. Therefore, a three-dimensional finite element model would cost unbearable computing power when considering the highly nonlinear frictional relations between plates. Therefore, the two-dimensional (2D) plane-strain finite element model was used to simulate the stick-slip action of faults caused by plate squeezing, at different crust depths.

The Wenchuan M_w 7.9 earthquake (2008) occurred along the Longmenshan fault, the central part of which provided the geological background for creating the large-scale model (plate squeezing model). Fig.2 presents a simplified plate squeezing model of the Longmenshan fault zone. The strike of the fault was NE 50°, the width and length of the numerical model were taken to be 300 km. To minimize potential boundary effects, the elements’ size near the fault plane was set to be less than 100 m based on a series of sensitivity analyses in the meshing, and then gradually transitioned to larger units with increasing distance, to a maximum size of 5 km. The north-eastern boundary of the Sichuan Basin plate was normally constrained, while a dynamic displacement with a squeezing rate of 5 mm per year [22] was loaded on the south-western boundary of the Xizang Plateau to simulate the plate squeezing. The stress boundary conditions varying with depth were simulated by applying horizontal stress to the large-scale model. According to the in-situ stress test results for locations in North America and South Africa [23], for depths in excess of 1 km, horizontal stress and vertical stress tend to be equalized, as suggested by Heim’s rule. Therefore, the horizontal stress can be given by:

(1){\sigma }_{\mathrm{H}}={\sigma }_{\mathrm{V}}=\sum _{i=1}^n{\gamma }_i\times h_i,

where {\sigma }_{\mathrm{H}} and {\sigma }_{\mathrm{V}} are horizontal and vertical stresses, respectively. {\gamma }_i is the unit weight of the ith stratum, and h_i is the thickness of the ith stratum.

It should be noted that the buried depth of the tunnel is almost negligible relative to the plate width and crustal depth; the displacement time history obtained at the ground (h = 0 km) in the stick-slip action was used as subsequent input displacement loading of the sub-model.

3.2 Burger’s viscoelastic constitutive model of the crustal lithosphere

The viscoelastic characteristics of the lithosphere have important consequences in plate movement analyses. The continental lithosphere generally has a sandwich-like layered rheological structure of the brittle upper crust, flexible lower crust, and strong upper mantle. Hoechner et al. [24] obtained the viscous coefficient of rock strata through the observation of surface morphology and long-term deformation monitoring after earthquakes, reflecting the viscoelasticity of the crustal structure. Burger’s viscoelastic constitutive model (Fig.3), composed of Maxwell and Kelvin models in series, not only reflects the rigid response under an instantaneous force but also the creep characteristics under a long-term force.

The differential equation of the constitutive model governing Burger’s model [25] can be expressed as:

(2)\sigma +p_1\cdot \dot{\sigma }+p_2\cdot \ddot{\sigma }=q_1\cdot \dot{\epsilon }+ q_2\cdot \ddot{\epsilon },

p_1,p_2,q_1 , and q_2 are defined as:

p_1=\frac{{\eta }_{\mathrm{M}}}{E_{\mathrm{M}}}+\frac{{\eta }_{\mathrm{M}}}{E_{\mathrm{K}}}+\frac{{\eta }_{\mathrm{K}}}{E_{\mathrm{K}}},p_2=\frac{{\eta }_{\mathrm{M}}{\eta }_{\mathrm{K}}}{E_{\mathrm{M}}{E}_{\mathrm{K}}},q_1={\eta }_{\mathrm{M}},q_2=\frac{{\eta }_{\mathrm{M}}{\eta }_{\mathrm{K}}}{E_{\mathrm{K}}},

where \sigma and \epsilon are the stress and strain, respectively; \dot{\sigma } and \ddot{\sigma } are the first and second derivatives of stress with respect to time, respectively; \dot{\epsilon } and \ddot{\epsilon } are the first and second derivatives of strain with respect to time, respectively; E_{\mathrm{M}} and {\eta }_{\mathrm{M}} are the elastic modulus and viscosity coefficient in the Maxwell model, respectively; E_{\mathrm{K}} and {\eta }_{\mathrm{K}} represent the elastic modulus and viscosity coefficient in the Kelvin model, respectively.

When the viscoelastic constitutive model is defined in Abaqus [26], Burger’s model parameters are converted into a Prony series. The shear relaxation modulus in Burger’s term after Laplace transformation can be expressed:

(3)G(t)=\frac{G_{\mathrm{M}}}{(\alpha -\beta )}\left[\left(\frac{G_{\mathrm{K}}}{{\eta }_{\mathrm{K}}^{\prime }}-\beta \right)\cdot {\mathrm{e}}^{-\beta t}+\left(\alpha -\frac{G_{\mathrm{K}}}{{\eta }_{\mathrm{K}}^{\prime }}\right)\cdot {\mathrm{e}}^{-\alpha t}\right],

where t is the time, and

\begin{array}{rl}& G_{\mathrm{M}}=\frac{E_{\mathrm{M}}}{2(1+v)},G_{\mathrm{K}}=\frac{E_{\mathrm{K}}}{2(1+v)},{\eta }_{\mathrm{M}}^{\prime }=\frac{{\eta }_{\mathrm{M}}}{2(1+v)},\\ & {\eta }_{\mathrm{K}}^{\prime }=\frac{{\eta }_{\mathrm{K}}}{2(1+v)},\alpha ,\beta =\frac{p_1\mp \sqrt{p_1^2-4p_2}}{2p_2}\end{array}

The shear relaxation modulus is considered as:

(4)G\text{(}t\text{)}=G_{\mathrm{\infty }}+\sum _n^{i=1}{G}_i\cdot {\text{e}}^{-t/{\tau }_i}.

Comparing Eq. (4) to Eq. (3), one can obtain:

(5)\begin{array}{rl}& G_{\mathrm{\infty }}=0,G_1=\frac{G_{\mathrm{M}}}{(\alpha -\beta )}\left(\frac{G_{\mathrm{K}}}{{\eta }_{\mathrm{K}}^{\mathrm{\prime }}}-\beta \right),\\ & G_2=\frac{G_{\mathrm{M}}}{(\alpha -\beta )}\left(\alpha -\frac{G_{\mathrm{K}}}{{\eta }_{\mathrm{K}}^{\mathrm{\prime }}}\right),\\ & {\tau }_1=\frac{1}{\beta },{\tau }_2=\frac{1}{\alpha }.\end{array}

The rest dimensionless variable in Abaqus can be expressed by:

(6)g_1=\frac{1}{(\alpha -\beta )}\left(\frac{G_{\mathrm{K}}}{{\eta }_{\mathrm{K}}^{\prime }}-\beta \right),g_2=\frac{1}{(\alpha -\beta )}\left(\alpha -\frac{G_{\mathrm{K}}}{{\eta }_{\mathrm{K}}^{\prime }}\right),

E_{\mathrm{M}},v, g_1,g_2,t_1,t_2 are parameters of Burger’s viscoelastic constitutive model in Abaqus. The relevant parameters of the lithosphere of the Sichuan Basin and the Xizang Plateau were determined based on related research on the lithosphere of mainland China [27,28], as listed in Tab.1.

3.3 Friction constitutive model of fault planes

stick-slip action is a periodic phenomenon relevant to frictional instability, and its process can be represented by Eq. (7). Under the initial condition, the hanging wall and footwall of the fault are bonded (the bonding state in Eq. (7a)), and the strain energy accumulated on the Longmenshan fault plane increases the squeezing pressure. Once the shear force increases to exceed the static frictional force, the plates slip on the fault plane (the sliding state in Eq. (7b)). With the decrease of frictional resistance, the sliding rate increases sharply, causing the strain energy to be released rapidly. The slip rate decreases in the process of strain energy release, and the friction coefficient increases again, so that the fault plane then tends to return to a bonding state, forming a cycle of sliding-bonding states.

\{\begin{array}{l}\tau \le \mu \times {\sigma }_{\text{n}},\mu ={\mu }_{\text{s}},\mathrm{b}\mathrm{o}\mathrm{n}\mathrm{d}\mathrm{i}\mathrm{n}\mathrm{g}\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{t}\mathrm{e},(7\mathrm{a})\\ \tau >\mu \times {\sigma }_{\text{n}},\mu ={\mu }_{\text{d}},\mathrm{s}\mathrm{l}\mathrm{i}\mathrm{d}\mathrm{i}\mathrm{n}\mathrm{g}\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{t}\mathrm{e},(7\mathrm{b})\end{array}

where τ is the tangential stress of the fault plane, μ is the friction coefficient, μ_d is the dynamic friction coefficient, μ_s is the static friction coefficient, and {\sigma }_{\mathrm{n}} is the normal contact stress.

The friction constitutive model plays an important role in the instability analysis of fault interfaces [29], and the static and dynamic friction coefficients are defined based on the slip-weakening, slip rate-weakening, and rate-state dependent laws. Among these laws, the slip-weakening law is limited when applied to the entire cycle simulation of fault slip because it does not consider the mechanism of relocking and strength recovery after fault instability and hence cannot explain the recurrence of stick-slip actions. The rate-state dependent law is derived in a laboratory environment with low slip rates (10⁻⁷ to 1 mm/s), which are very different from the magnitudes (measured in m/s) for actual large earthquakes [30]. Moreover, the rate-state dependent law [31,32] requires additional parameters that are difficult to determine. Only a static-kinetic exponential decay definition needs to be determined to construct the decay function of the friction coefficient in the slip rate-weakening model.

(7)\mu ={\mu }_{\mathrm{d}}+\left({\mu }_{\mathrm{s}}-{\mu }_{\mathrm{d}}\right){\mathrm{e}}^{-dc|v|},

where v is the sliding rate, and dc is the decay coefficient. Based on previous studies [33,34], μ_s, μ_d, and dc were set to 0.1, 0.6, and 0.2, respectively.

3.4 strike-slip faulting analysis subjected to plates squeezing

Fig.4 illustrates the time histories of shear stress on the Longmenshan fault plane at depths of 1, 5, 10, and 15 km below the surface. It can be seen that the shear stress increases gradually under the action of plates squeezing, and it drops drastically when reaching a certain magnitude, which corresponds to a stick-slip action on the fault plane. The cycle of the increase and decrease of shear stress indicates the seemingly periodic stick-slip action. At different crustal depths, the shear stress threshold reached by each stick-slip action is different. The shear stress threshold ranges from 35.9 to 53.4 MPa for h = 1 km, 102.4 to 233.4 MPa for h = 5 km, 204.9 to 499.4 MPa for h = 10 km, and 268.9 to 821.8 MPa for h = 15 km. Furthermore, the stress accumulation and release process on the fault plane in stick-slip action differs slightly after the first stick-slip action; that is, the recurrence interval is relatively stable. Tab.2 lists the key parameters of the fault subjected to plate squeezing at different depths. It can be seen that the larger the crustal depth, the later the first stick-slip action, and the longer the recurrence interval. At the crustal depth of 15 km, the stick-slip action interval is 3234 years, basically consistent with the recurrence interval of 2233 years to 4167 years for the Longmenshan Tectonic Zone provided by Zhang et al. [35] based on field investigations and seismic inversions.

3.5 Time histories of displacement on the Longmenshan fault plane at the ground

Since the buried depth of a tunnel is almost negligible relative to the plate width and crustal depth, as mentioned in Subsection 3.1, the displacement time history of the fault plane obtained at the ground (h = 0 km) in a stick-slip action is herein approximated as the displacement loading applied to the sub-model of a tunnel subjected to stick-slip action. Fig.5 shows the time histories of acceleration and velocity on the fault plane at the ground. Similarly to the several stick-slip action illustrated in Fig.4, the thresholds of acceleration and velocity for h = 0 km in Fig.5(a) are not uniform, range from −1.15 to −2.28 m/s², and 0.10 to 0.13 m/s (Fig.5(b)), respectively. Unlike the variation of shear stress, the acceleration and velocity at each instant between adjacent stick-slip action are very small (nearly equal to 0), indicating that the acceleration and velocity are induced by stress release and have little relation with stress accumulation. Fig.6(a) shows the time histories of accumulated displacement in the process of plates squeezing. It can be seen that the displacement is accumulated in each stick-slip action, and the average accumulated displacement in a stick-slip action is about 50 cm. The time history of incremental displacement in a certain stick-slip action is presented in Fig.6(b), the incremental displacement is accumulated to reach a maximum value of 49.19 cm at an increasing rate from fast to slow. A small rebound occurs when the incremental displacement reaches the maximum value; this may be induced by the pulse of acceleration or velocity in the stick-slip action.

4 Sub-model of a tunnel subjected to strike-slip faulting

4.1 Numerical model of a tunnel subjected to strike-slip faulting

Based on the observation results of surface rupture and focal mechanism of the main earthquake, combined with the geophysical data of Longmenshan and the Sichuan Basin [2,36], the earthquake-related strike-slip faults in the Longmenshan area may be characterized by a dip angle of 50°–60°. On this basis, a sub-model of the tunnel subjected to the strike-slip faulting with a fault dip angle of 60° was established in combination of the extended finite element method.

Fig.7 shows the established three-dimensional finite element model. The overall dimensions of the numerical model were 200 m × 100 m × 100 m, and the inner width and height of the tunnel were 12.5 and 9.6 m, respectively, with a lining thickness of 0.5 m, as shown in Fig.7. The distance from the top of the tunnel to the ground surface was 40 m. To minimize the boundary effects, a series of sensitivity analyses were conducted to ensure that the boundary effects on the mechanical response of the tunnel were negligible by adopting the geometries in Fig.7. In this study, the strata of hanging wall (moving block) and footwall (fixed block) were simulated by eight-node reduced-integration brick elements (C3D8R), while the tunnel was modeled by eight-node brick elements (C3D8). Since the effect of fault width was neglected in the large-scale model, a contact surface without considering fault width was used for the fault rupture simulation in the sub-model, with reference to previous experimental and numerical analyses [37,38]. To improve computational efficiency and ensure computational accuracy, much finer meshes were applied to the elements of the stratum close to the fault plane, where the mesh size in the longitudinal direction was set to 2.5 m. Similarly, the mesh size of the tunnel was set to 1.0 m in the longitudinal direction.

The time history of incremental displacement (Fig.6(b)) obtained in the large-scale model was considered to be the displacement loading at the nodes on the external surfaces on the hanging wall (moving block) in the sub-model, for consideration of the dynamic effects of the fault movement on the mechanic response of the tunnel, as illustrated in Fig.7 and Fig.8. The maximum of the incremental displacement was about 50 cm (49.19 cm in Fig.6(b)). Wang [39] conducted a series of laboratory tests on the static response of a tunnel a tunnel subjected to strike-slip faulting processes and obtained the typical structural deformation and failure modes. The reliability of the numerical model can be verified to some extent by comparing with the experimental results in Ref. [39].

In the local finite element model, the lateral boundaries of both the hanging wall and footwall were fixed in the x direction, while the bottom boundaries were fixed in the y direction. The time history of incremental displacement was loaded on the squeezed side of the hanging wall, while the boundary on the opposite side of the footwall was fixed in the z direction. The other boundaries of the footwall and hanging wall were free, as shown in Fig.8.

A hard contact using the Lagrange multiplier method was assumed in the normal direction, whereas a friction contact using the penalty method was assumed in the surface-parallel direction to track the interface behavior on the Longmenshan fault plane. According to BSI [40], the friction coefficient was herein assumed to be 0.4. The concrete lining and the stratum were assumed to be in full contact, in line with the parameters in Dadashi et al. [41]. The elastic-perfectly plastic Mohr–Coulomb constitutive model was adopted for simulating the three types of stratums, while the linear elastic constitutive model was used to simulate the C35 class concrete tunnel lining. Their basic parameters were determined according to the guidelines for “Design of Highway Tunnel” [42], as listed in Tab.3. The damaged evolution could be formulated in terms of fracture energy, with reference to published studies [43,44]; the fracture energy value of the concrete was herein assumed to be 80 N/m. The tensile strength of the concrete was 1.57 MPa.

There are several damage generation criteria in ABAQUS to predict the propagation of cracks. The maximum principal stress criterion was selected herein to analyze the cracking characteristics of the tunnel lining. It can be formulated by:

(8)f=\left\{\frac{⟨{\sigma }_{max}⟩}{{\sigma }_{max}^0}\right\},

where {\sigma }_{max} and {\sigma }_{\mathrm{mmx}}^0 are the maximum and the maximum principal stress, respectively, ⟨⟩ is the Macaulay bracket, signifying that a pure compressive stress state cannot generate damage. f = 1 indicates the initiation of damage.

4.2 Comparison of cracking propagation with the experimental results

A series of experimental model tests with a geometric similarity ratio of 1:100 was conducted by Wang [39] to investigate the cracking process and mechanical response of the stratum and lining under strike-slip faulting. The density, elastic modulus, and Poisson ratio of the actual tunnel lining in Ref. [39] were 2400 kg·m⁻³, 30 GPa, and 0.2, respectively, and the matching figures for the C35 class concrete in this study were ρ = 2500 kg·m⁻³, E = 31.5 GPa, v = 0.2. The elastic modulus and Poisson ratio of the original stratum in Ref. [39] were 11 GPa and 0.28, matching the material parameters of stratum I in this study, E = 10 GPa, v = 0.3. Therefore, the rationality of this sub-model method can be validated by comparing the numerical results with the experimental results. Please refer to Wang [39] for more details about the model test.

Fig.9 shows that an initial crack first occurred at the central section on the left wall, corresponding to the 1# crack in this numerical simulation. As the stratum was subjected to further squeezing, another crack occurred at the inverted arch and extended to the central section of the right wall, corresponding to the 2# crack in this numerical simulation. The 1# and 2# cracks extended constantly under strike-slip faulting to form a closed circumferential crack when they were connected. This phenomenon can be observed in both the model test of Wang [39] and this numerical simulation. It is undeniable that the numerical results of this study varied somewhat from those of the model test of Wang [39] due to the differences in the shape and size of the tunnel, displacement loading, and fault characteristics. However, on the whole, the numerical results were consistent with the experimental results of Wang [39], verifying the rationality of this numerical simulation.

4.3 Lining cracking propagation process for different strata

Fig.10 shows the crack propagation process of the tunnel lining in different strata. To clarify the expression, Fig.10 also illustrates a straight line of the oblique-cut tunnel cross-section as the projection of the fault plane of the tunnel. The gray bold arrows in Fig.10 indicate the direction of hanging wall dislocation. The central axis of the tunnel cross-section is considered as the boundary, the squeezed side of the hanging wall tunnel represents the positive side of the Z-axis, and it is termed the P side of the model, whereas the negative side is termed the N side. Regardless of the stratum type, the following conclusions can be drawn. 1) An initial crack (2# crack) occurred at the inverted arch in the footwall on the squeezed side (P side) of the fault plane, and then spread to the negative side (N side) at a certain angle to form a new inclined crack. 2) The initial crack also extended upwards to the central section and spandrel after reaching the arch foot on the P side, but this phenomenon was not obvious on the N side. 3) When viewed from the top, the propagation direction of the cracks formed at the inverted arch and vault was at an acute angle with the presupposed displacement direction of the hanging wall; that is, the inclination directions of 2# crack and newly formed crack were consistent.

Differences also existed in the lining crack propagation process in different strata as can be seen by comparing Fig.10(a)–Fig.10(c). The following conclusions can be drawn. 1) For the tunnel in stratum I, in addition to the initial crack (2# crack) at the inverted arch on the P side, an inclined crack (1# crack) also occurred at the tunnel vault on the N side and expanded to connect 2# crack, nearly forming a ring crack. 2) For the tunnel in stratum II, 1# crack did not occur, but another crack (3# crack) occurred at the central section in the footwall on the P side and expanded to connect 2# crack to form a “V”-shape crack. 3) For the tunnel in stratum III, 1# and 3# crack did not occur; only 2# crack occurred at the inverted arch in the footwall on the P side and expanded along both ends of the crack. Through preliminary analysis, we could conclude that these differences may be attributed to the effect of the stratum stiffness on the local shear, squeezing, and bending action of the strike-slip fault on the tunnel lining, thus affecting the formation and propagation of lining cracks.

Fig.11 presents the time history of crack opening (2# crack) at the inverted arch for different strata. The maximum crack openings in strata I, II, and III were 10.3, 8.3, and 1.8 cm, respectively, decreasing with decreasing stratum stiffness. The crack opening herein was calculated by extracting the time histories of relative displacement at the cracks. The distribution range of longitudinal cracks along the tunnel became wider with decreasing stratum stiffness, while the acute angle between the propagation direction of the cracks and the longitudinal direction of the tunnel decreases and the extension length became longer with decreasing stratum stiffness.

5 Mechanical response of the tunnel subjected to strike-slip faulting

5.1 Lining stress analysis of initial cracking for different strata

It can be observed from Fig.12 that the shear stress of the 2# crack at the inverted arch in the footwall on the P side was far less than the maximum shear stress. By analyzing the position and inclined direction of the 2# crack, we could conclude that this crack was induced by the bending moment of the tunnel rotating around the vertical direction (Y-axis) under the strike-slip faulting of the hanging wall. Since the thickness direction of the inverted arch of a small radius and large width was approximately parallel to the vertical direction (Y-axis), the central section and spandrel of the lining had large radius, causing the longitudinal bending moment to mainly act on the inverted arch. Therefore, regardless of the stratum type, an initial crack (2# crack) always occurred at the inverted arch in the footwall on the squeezed side (P side) of the fault plane under strike-slip faulting.

Since initial cracks were determined based on the maximum principal stress principle, it can be seen from Fig.12 that the position of initial lining cracks corresponds to the position of a maximum of the maximum principal stress. The maximum principal stress appeared at the tunnel vault on the N side for stratum I. However, the maximum principal stress appeared at the inverted arch in the footwall on the squeezed side for strata II and III. The maximum principal stress at the tunnel vault on the N side for stratum I was 1.64 MPa, which exceeded the tensile strength of concrete of 1.57 MPa, leading to the occurrence of 1# crack. For strata II and III, the maximum principal stresses at the tunnel vault on the N side were 1.56 and 1.11 MPa, which were smaller than the tensile strength of concrete of 1.57 MPa, resulting in no cracks at this position. For the same position at the tunnel vault, the maximum principal stresses for strata I, II, and II were 1.64, 1.56, and 1.11 MPa, respectively, indicating that the stress concentration decreased with decreasing stratum stiffness.

The above analyses of shear stress and maximum principal stress reasonably explain the mechanism of initial cracking described in Subsection 4.3. It is further concluded that for the tunnel in the high-stiffness stratum, both shear and bending failures occur on the lining under strike-slip faulting, while for that in the low stiffness stratum, only bending failure occurs on the lining.

5.2 Lining deformation analysis for different strata

Fig.13 shows the deformation state of the tunnel and stratum under strike-slip faulting. It can be observed that both sides (P side and N side) of the tunnel had obvious squeezing deformation, and large dislocation of the stratum could also be observed in the hanging wall and footwall. Fig.14 shows the deformation at the spandrel node along the longitudinal direction of the tunnel under strike-slip faulting. In combination with Fig.13(a), it can be seen that the spandrel of the tunnel on the N side is subjected to severe lateral squeezing due to the fault dislocation between the hanging wall and footwall. This led to a significant lateral tension at the spandrel of the tunnel and a stress concentration of maximum principal stress at this position. For the tunnel in stratum I, the deformation at the spandrel along the longitudinal direction had an obvious sharp change near the fault plane, indicating that the tunnel near the fault plane was subjected to a large transverse shear force. However, for the tunnel in stratum II and stratum III, the deformation curves of the spandrel along the longitudinal direction tunnel were gentle, in agreement with the conclusion in Subsection 5.1 that the stress concentration at the spandrel on the N side decreased with decreasing stratum stiffness. Therefore, it can be concluded that the large dislocation of the stratum between the hanging wall and footwall during strike-slip faulting resulted in the lateral tensile stress at the spandrel of the tunnel on the N side, due to the local bending moment. For the tunnel in the high stiffness stratum (stratum I), the tunnel lining near the fault plane cracked first due to the significant constraint effect of the stratum.

Fig.15 shows the deformation at the foot node along the longitudinal direction of the tunnel under strike-slip faulting. It can also be seen that the direction of longitudinal bending of the tunnel rotating around the vertical direction (Y axis) is transitional on the fault plane. According to the analysis of the lining cracking process provided in Subsection 4.3, the deformation indices of the tunnel and stratum were obtained, as listed in Tab.4. It can be seen that the longitudinal bending deformation ranges of the tunnel for strata I, II, and III were 35, 54, and 88 m, respectively, and the maximum squeezing deformations of the stratum for strata I, II, and III were 9.60, 21.46, and 24.12 m, respectively. Also, the angles between the crack and longitudinal direction of the tunnel for strata I, II, and III were 60°, 57°, and 37°, respectively. Once again, this phenomenon () illustrates that, with the decrease in the stratum stiffness, the crack expansion range and the maximum squeezing deformation of the stratum become larger, while the included angle between the crack and the longitudinal direction of the tunnel decreases.

6 Conclusions

In this study, the conditions in the north-east zone of the Longmenshan (Sichuan) fault were considered as the geological background. A large-scale model was established to investigate the stick-slip action of the fault subjected to plate squeezing and to obtain the dynamic parameters of the fault in the stick-slip action. The dynamic displacement obtained in the large-scale model was considered as the displacement loading boundary condition in the sub-model to consider the effects of regional geological environments and the dynamic effects of the fault movement on the mechanic response of the tunnel subjected to stick-slip action of strike-slip faulting. By means of this proposed multi-scale modeling method, the following conclusions can be drawn.

1) The large-scale model incorporating Burger’s viscoelastic constitutive model and friction constitutive model can successfully simulate the stick-slip action of the fault induced by plate squeezing. The stick-slip action are seemingly periodic; the greater the crustal depth, the later the first strike-slip occurrence and the longer the recurrence interval.

2) Regardless of the stratum stiffness, an initial crack first occurs at the inverted arch of the tunnel in the footwall on the squeezed side under strike-slip faulting. A new crack also occurs at the tunnel vault on the non-squeezed side in stratum I, but not in strata II and III.

3) For the tunnel in the high-stiffness stratum (strata I and II), both shear and bending failures occur on the lining under strike-slip faulting, while for the tunnel in the low-stiffness stratum (stratum III), only bending failure occurs on the lining.

4) With the decrease in the stratum stiffness, the crack expansion range and maximum squeezing deformation of the stratum increase, while the angle between the crack and the longitudinal direction of the tunnel decreases.

In this work, the mechanical response of a tunnel subjected to stick-slip actions of strike-slip fault was preliminarily studied to provide a reference of structural analysis for similar tunnels subjected to strike-slip faulting. The effects of construction joint, steel, fault width, and other factors on the dynamic performance of the tunnel were not taken into account. Furthermore, the width and expansion range of cracks may vary somewhat from an actual situation. The above issues will be considered emphatically in future work.

References

[1] Kontogianni V A, Stiros S C. Earthquakes and seismic faulting: Effects on tunnels. Turkish Journal of Earth Sciences, 2003, 12(1): 153–156

[2] Liu-Zeng J, Zhang Z, Wen L, Tapponnier P, Sun J, Xing X, Hu G, Xu Q, Zeng L, Ding L, Ji C, Hudnut K W, van der Woerd J. Co-seismic ruptures of the 12 May 2008, Ms 8. 0 Wenchuan earthquake, Sichuan: East–west crustal shortening on oblique, parallel thrusts along the eastern edge of Tibet. Earth and Planetary Science Letters, 2009, 286(3–4): 355–370

[3] Huang Z, Fu H, Zhang J, Chen W, Shi Y. Structural damage evaluation method for metro shield tunnel. Journal of Performance of Constructed Facilities, 2019, 33(1): 04018097

[4] Ranjbarnia M, Zaheri M, Dias D. Three-dimensional finite difference analysis of shallow sprayed concrete tunnels crossing a reverse fault or a normal fault: A parametric study. Frontiers of Structural and Civil Engineering, 2020, 14(4): 998–1011

[5] Sabagh M, Ghalandarzadeh A. Centrifugal modeling of continuous shallow tunnels at active normal faults intersection. Transportation Geotechnics, 2020, 22: 100325

[6] Wang Q, Geng P, Chen J B, He C. Dynamic discrimination method of seismic damage in tunnel portal based on improved wavelet packet transform coupled with Hilbert–Huang transform. Mechanical Systems and Signal Processing, 2023, 188: 110023

[7] Wang Q, Geng P, Chen C J, Chen J B, He C. Determination of seismic response of reinforced tunnel portal slope using shaking table tests. Tunnelling and Underground Space Technology, 2023, 136: 105072

[8] NewmarkN MHall W J. Pipeline design to resist large fault displacement. In: Proceedings of US National Conference on Earthquake Engineering. Ann Arbor, MI: Department of Energy Office of Scientific and Technical Information, 1975, 416–425

[9] Qiao Y F, Tang J, Zhang X D. Longitudinal mechanical response of tunnels under reverse faulting and its analytical solution. In: Proceedings of ISRM International Symposium-Asian Rock Mechanics Symposium. Beijing: ISRM, 2021, 861(2): 022015

[10] Kiani M, Akhlaghi T, Ghalandarzadeh A. Experimental modeling of segmental shallow tunnels in alluvial affected by normal faults. Tunnelling and Underground Space Technology, 2016, 51: 108–119

[11] Sabagh M, Ghalandarzadeh A. Centrifuge experiments for shallow tunnels at active reverse fault intersection. Frontiers of Structural and Civil Engineering, 2020, 14(3): 731–745

[12] Zhou G, Sheng Q, Cui Z, Wang T, Ma Y. Investigating the deformation and failure mechanism of a submarine tunnel with flexible joints subjected to strike-slip faults. Journal of Marine Science and Engineering, 2021, 9(12): 1412

[13] Qiao Y, Tang J, Liu G, He M. Longitudinal mechanical response of tunnels under active normal faulting. Underground Space, 2022, 7(4): 662–679

[14] Ma Y, Sheng Q, Zhang G, Cui Z. A 3D discrete-continuum coupling approach for investigating the deformation and failure mechanism of tunnels across an active fault: A case study of xianglushan tunnel. Applied Sciences, 2019, 9(11): 2318

[15] Zeng G, Geng P, Guo X, Li P, Wang Q, Ding T. An anti-fault study of basalt fiber reinforced concrete in tunnels crossing a strike-slip fault. Soil Dynamics and Earthquake Engineering, 2021, 148: 106687

[16] Liu X, Wang Y. Three dimensional numerical analysis of underground bifurcated tunnel. Geotechnical and Geological Engineering, 2010, 28(4): 447–455

[17] Cao Y, Wang P, Jin X, Wang J, Yang Y. Tunnel structure analysis using the multi-scale modeling method. Tunnelling and Underground Space Technology, 2012, 28: 124–134

[18] Sun L F, Chen Z L, Li J C, Lv X L, Shang Q. Weld strength analysis of T-joint segments of the metro crossing passage by the shield method based sub-model. Materials Science Forum, 2019, 971: 27–35

[19] Zhang X, Fu P, Zhou C, Yin J, Lu W. Construction of 3D geological model and multi-scale numerical analysis of in-situ stress field of long tunnel in fault development area. ISRM International Symposium-Asian Rock Mechanics Symposium, 2021, 861(6): 062053

[20] Krishnan V R, Hui C Y, Long R. Finite strain crack tip fields in soft incompressible elastic solids. Langmuir, 2008, 24(24): 14245–14253

[21] Lucht T. Finite element analysis of three dimensional crack growth by the use of a boundary element sub model. Engineering Fracture Mechanics, 2009, 76(14): 2148–2162

[22] Meyer B, Tapponnier P, Bourjot L, Metivier F, Gaudemer Y, Peltzer G, Guo S M, Chen Z T. Crustal thickening in Gansu–Qinghai, lithospheric mantle subduction, and oblique, strike-slip controlled growth of the Tibet plateau. Geophysical Journal International, 1998, 135(1): 1–47

[23] BrownE THoek E. Underground Excavations in Rock. Boca Raton, FL: CRC Press, 1980

[24] HoechnerASobolev S VEinarssonIWangR. Investigation on afterslip and steady state and transient rheology based on postseismic deformation and geoid change caused by the Sumatra 2004 earthquake. Geochemistry, Geophysics, Geosystems, 2011, 12(7): n/a

[25] FindleyW NDavis F A. Creep and Relaxation of Nonlinear Viscoelastic Materials. Chelmsford: Courier Corporation, 2013

[26] DassaultS. ABAQUS 6.14 analysis User’s Manual. Providence, RI: Dassault Systems, USA, 2014

[27] Clark M K, Royden L H. Topographic ooze: Building the eastern margin of Tibet by lower crustal flow. Geology, 2000, 28(8): 703–706

[28] Hu J, Xu X, Yang H, Wen L, Li G. S receiver function analysis of the crustal and lithospheric structures beneath eastern Tibet. Earth and Planetary Science Letters, 2011, 306(1–2): 77–85

[29] Scholz C H. Earthquakes and friction laws. Nature, 1998, 391(6662): 37–42

[30] Heaton T H. Evidence for and implications of self-healing pulses of slip in earthquake rupture. Physics of the Earth and Planetary Interiors, 1990, 64(1): 1–20

[31] Carlson J M, Langer J S. Mechanical model of an earthquake fault. Physical Review A: General Physics, 1989, 40(11): 6470–6484

[32] Aagaard B T, Hall J F, Heaton T H. Characterization of near-source ground motions with earthquake simulations. Earthquake Spectra, 2001, 17(2): 177–207

[33] Yuan J, Zhu S B. FEM simulation of the dynamic processes of fault spontaneous rupture. Chinese Journal of Geophysics, 2014, 57(1): 138–156

[34] Ma L, Tao W, Zhang Y, Zeng M, Zheng Q. The numerical simulation study of the earthquake cycles and the dynamic evolutionary processes on the Longmen Shan Fault. Chinese Journal of Geophysics, 2018, 61(5): 1824–1839

[35] ZhangP ZXu X WWenX ZRanY K. Slip rates and recurrence intervals of the Longmen Shan active fault zone, and tectonic implications for the mechanism of the May 12 Wenchuan earthquake, 2008, Sichuan, China. Chinese Journal of Geophysics, 2008, 51(4): 1066–1073 (in Chinese)

[36] Xu X, Wen X, Yu G, Chen G, Klinger Y, Hubbard J, Shaw J. Coseismic reverse-and oblique-slip surface faulting generated by the 2008 Mw 7. 9 Wenchuan earthquake. China Geology, 2009, 37(6): 515–518

[37] Liu X, Li X, Sang Y, Lin L. Experimental study on normal fault rupture propagation in loose strata and its impact on mountain tunnels. Tunnelling and Underground Space Technology, 2015, 49: 417–425

[38] Li H, Li X, Yang Y, Liu Y, Ma M. Structural stress characteristics and joint deformation of shield tunnels crossing active faults. Applied Sciences, 2022, 12(7): 3229

[39] WangT Q. Research on the forced influence of cross-active fault tunnel under strike-slip dislocation. Thesis for the Master’s Degree. Shenyang: Shenyang University of Technology, 2021 (in Chinese)

[40] BritishStandards Institution. Code of Practice for Temporay Works Procedures and The Permssible Stress Design of Falsework: Incorporating Corrigendum No. 1. London: BSI, 2008

[41] Dadashi E, Noorzad A, Shahriar K, Goshtasbi K. Hydro-mechanical interaction analysis of reinforced concrete lining in pressure tunnels. Tunnelling and Underground Space Technology, 2017, 69: 125–132

[42] JTG/TD70-2010. Guidelines For Design of Highway Tunnel. Beijing: Ministry of Transport of the People’s Republic of China, 2010

[43] Giner E, Sukumar N, Tarancón J E, Fuenmayor F J. An Abaqus implementation of the extended finite element method. Engineering Fracture Mechanics, 2009, 76(3): 347–368

[44] Min B, Zhang C, Zhang X, Wang H, Li P, Zhang D. Cracking performance of asymmetric double-arch tunnels due to the voids behind linings. Thin-Walled Structures, 2020, 154: 106856

RIGHTS & PERMISSIONS

Higher Education Press