Progressive failure process of secondary lining of a tunnel under creep effect of surrounding rock

Tunnelling and Underground Space Technology 90 (2019) 76–98

Contents lists available at ScienceDirect

Tunnelling and Underground Space Technology journal homepage: www.elsevier.com/locate/tust

Progressive failure process of secondary lining of a tunnel under creep effect of surrounding rock

T

Guowen Xu , Chuan He, Qinhao Yang, Bo Wang ⁎

Key Laboratory of Transportation Tunnel Engineering, Ministry of Education, Southwest Jiaotong University, Chengdu 610031, Si Chuan, China

ARTICLE INFO

ABSTRACT

Keywords: Creep Secondary lining Progressive failure Discrete element method

During the operation of tunnels, a problem that is encountered often is the cracking of secondary lining caused by the time-dependent deformation of weak surrounding rocks. In this article, the Dujia mountain tunnel, which is situated in a weak and crushing phyllite stratum, was taken as an example to study this issue. First, we conducted a field investigation of the cracking situation of the secondary lining. Then, the similarity model test and particle-discrete-element numerical approach were used to study the failure process of the secondary lining. The results showed that the similarity model test and the results of numerical simulation were consistent with respect to the order in which cracks occurred and the final failure pattern. The progressive failure process can be divided into four stages, i.e., (1) the elastic stage, (2) the initial damage stage, (3) the ultimate limit stage, and (4) the instability stage. For a lateral pressure coefficient (λ) > 1, cracks initially appeared at the left haunch, the right haunch, and the spandrel due to bending damage caused by large-eccentricity compression. For λ = 1, cracks appeared first at the left arch springing because of the compressive and shear damage. For λ < 1, the bending damage first occurred at the spandrel. For the secondary lining with pre-existing cracks, micro-cracks were first generated where pre-existing cracks initially existed, and one or two hinges appeared at the location of the pre-existing cracks when the structure was in its ultimate limit stage. A void behind lining influenced the cracking pattern of the secondary lining in three different ways, i.e., (1) no effect; (2) the order in which cracks emerged was affected, but the final fracture pattern was not; and (3) the order in which cracks emerged and the final fracture pattern were affected.

1. Introduction Among the composite supporting structures, which is used extensively in tunnels that are designed and constructed using the New Austrian Tunnelling Method (NATM), the secondary lining (lining for short) is the main guarantee for their long-term safe operation (Tonon, 2016). Linings are prone to form cracks under the action of various environmental factors, and their cracking patterns can be divided into two categories, i.e., sudden failure and progressive failure. Sudden failure is caused mainly by the sudden increase of a local, concentrated load. For example, a large area of the side wall of Mo tianling tunnel (situated in Chongqing, China) collapsed because of the local increase in water pressure when karst caves behind the lining formed catchment channels after continuous rainfall (Fig. 1). Progressive failure is more common than sudden failure as the result of cracking, aging, and other diseases occurred in linings during their operation. In addition, the extent of the damage to linings increases as time passes. Thus, linings must be reinforced or replaced by grouting, arching, or dismantling

⁎

(Asakura and Kojima, 2003) to maintain their load-bearing capacity. Field observation is an intuitive way to investigate the cracking of linings. Inokuma and Inano (1996) found that the main factors that caused defects in the linings of tunnels on Japanese highways were variations in the ground pressure, degradation of materials, seepage, and frost damage. Sandrone and Labiouse (2011) conducted a survey of the service status of tunnels on the Swiss national highway and summarized the laws and influential factors that contributed to the deterioration of linings. Wang (2010) and Chiu et al. (2017) analyzed the spatial distribution, morphology, and cause of surficial cracks of linings due to the slipping action of adjacent slopes. Xiao et al. (2014) investigated the cracking mechanism of bias tunnels in loose strata. Bian et al. (2016) found that the initial damage to linings and their insufficient thicknesses were the main factors that contribute to the cracking of a bifurcation tunnel at the Huizhou pumped storage power station. The similarity model test can record the cracking and failure process of linings, so researchers are very interested in using it. Singh et al.

Corresponding author. E-mail address: [email protected] (G. Xu).

https://doi.org/10.1016/j.tust.2019.04.024 Received 21 April 2018; Received in revised form 25 February 2019; Accepted 22 April 2019 Available online 29 April 2019 0886-7798/ © 2019 Elsevier Ltd. All rights reserved.

Tunnelling and Underground Space Technology 90 (2019) 76–98

G. Xu, et al.

Nomenclature

cb Ψ N h Eh ηk E Eg ¯c k¯ n kn µc k¯ s

List of symbols M b Eco Ek Eb Ec λ ¯c k¯ s

ks k¯ n µ

bending moment of lining (N·m) width of lining (m) elastic modulus of lining (Pa) elastic modulus of the Kelvin element (Pa) elastic modulus of bolt elastic modulus of sprayed concrete lateral pressure coefficient shear strength of parallel bonded contact (Pa) shear stiffness of parallel bonded contact (Pa·m−1) shear stiffness of particle (Pa·m−1) normal stiffness of smooth joint bond (Pa·m−1) friction coefficient of smooth joint bond

c b

(1997) studied the relationship between the thickness of linings and their bearing capacity under different load patterns. Lee et al. (2002) analyzed the distribution of cracks, cracking loads, and the ultimate bearing capacity of linings based on different loadings. In Japan, the Railway Technical Research Institute (2005) investigated the bearing capacity of the lining of the Shinkansen tunnel with different degrees of cracking. He and Se (2006) systematically studied the influence of construction defects on the morphology of cracks in linings. Lei et al. (2015) studied the mechanical characteristics and failure process of linings under bias loads. Fang et al. (2016) determined the effect of voids on the structural failure of linings in water-rich strata. With the advance of numerical simulation technology, many advanced numerical methods capturing cracking and fracture of materials have been put forward, such as strong discontinuity embedded approach and cracking elements method (cracking is assumed to localize in a discrete surface with initial zero width. Foster et al., 2007; Armero and Garikipati, 1996), extended finite element method (this approach is characterized by the incorporation of additional degrees of freedom through nodal enrichments based on the concept of partition of unity. Borja and Liu, 2010; Hansbo and Hansbo, 2004), phase field method (A phase-field variable is used to describe a smooth transition between the damaged and undamaged phases. Miehe et al., 2010; Nguyen and Wu, 2018), and cracking particles method (The crack is modeled by a set of cracking segments that pass through the nodes. Rabczuk et al., 2009; Caleyron et al., 2012). Thus, as supplements to the similarity model tests, many numerical simulations have been conducted to reveal the cracking patterns of lining. Shi et al. (2001) used the discrete crack approach to determine the cracking modes of linings with vault voids. Huang et al. (2013) studied the distribution and propagation of cracks under the action of bias load, voids, and ground pressure by using the

cohesion of smooth joint bond (Pa) dilation angle of smooth joint bond (°) axial force of lining (N) thickness of lining (m) elastic modulus of the Hoke element (Pa) viscosity coefficient of the Kelvin element (Pa·y) elastic modulus of sprayed concrete-steel arch composite elastic modulus of steel arch normal strength of parallel bonded contact (Pa) normal stiffness of parallel bonded contact (Pa·m−1) normal stiffness of particle (Pa·m−1) friction coefficient between particles shear stiffness of smooth joint bond (Pa·m−1) normal strength of smooth joint bond (Pa) Friction angle of smooth joint bond (°)

extended finite element method. David et al. (2014) analyzed the cracking patterns of linings that were bearing concentrated or distributed loads based on lumped damage mechanics theory. Based on the theory of the stress intensity factor and the finite element method, Gao et al. (2017) calculated the critical cracking thickness of linings. It is apparent that the previous studies focused mainly on the influence of construction defects on the mechanical behavior of linings. However, the time-dependent behavior of the surrounding rock also is an important factor that can lead to the cracking of linings. The lining of the Zhegu mountain tunnel (Meng et al., 2013) on National Road 317 had a large number of cracks after 10 years of operation due to the creep of phyllite. Guan et al. (2008) conducted long-term safety monitoring of the Ureshino tunnel and found that sustained deformation of lining had occurred during its operation; further study (Guan et al., 2009) indicated that the creep of the surrounding rock was the main factor that led to this sustained deformation, and reinforcement methods, including grouting and applying bolts at arch springing, ultimately were used to suppress this deformation. In this article, the Dujia mountain tunnel on the Guanggan highway in China was taken as an example to study the long-term safety of linings under the creep effect of weak rock. First, we conducted a detailed investigation of the formation and distribution of cracks. Then, a large-scale similarity model test was used to study the cracking and failure process of the lining adopting the overloading method. Subsequently, a numerical assessment of the lining was conducted based on the discrete element method, and we used a discrete-element, finite-difference coupling algorithm to determine its long-term safety with the cases of a lining with no faults, a lining with pre-existing cracks, and a lining that had a void behind it. The purpose of this article was to determine and present the progressive failure process of secondary linings under creep loads and provide a basis for the optimization of their structural design and safety assessment during the operational period of tunnels. 2. Engineering overview 2.1. Engineering background The Dujia mountain tunnel (Fig. 2) is located near the Longmen mountain fault, which is the seismogenic fault causing the Wenchuan earthquake in 2008. The tunnel is a two-way, four-lane, separated tunnel (the way from Guanyuan to Yaodu is called T-1, and the reverse way is called T-2) with a length of 1846 m, a width of 13 m, a height of 10.5 m, and a maximum depth of 180 m. Three thrust faults (F1, F2, F3) developed at the site of this tunnel. F1 crossed through the tunnel with an intersection angle of 79°. The tunnel is situated in a weak and crushing phyllite stratum,

Fig. 1. Sudden destruction of secondary lining. 77

Tunnelling and Underground Space Technology 90 (2019) 76–98

G. Xu, et al.

Guang yuan Fig. 2. Location of the Dujia mountain tunnel.

excavation method with reserved core soil was used, as shown in Fig. 6. The layout of the components of supporting system and their parameters are shown Fig. 6 and Table 1.

which easily can be softened by water (Fig. 3), and it contains welldeveloped weak planes. In addition, the tunnel is affected by periodic aftershocks of the Wenchuan earthquake. Thus, the combination of soft rock and aftershocks multiplies the problem of large deformation (Fig. 5) during its construction. The deformation of the tunnel had the following characteristics: (1) Rock collapse, which was concentrated in the areas from the haunches to the vaults and occurred frequently after the unstable, large deformation. (2) The rate of deformation of the tunnel increases as the water content of the surrounding rock increases, especially in the fault zone. (3) Secondary collapses could occur in some areas. For example, after a collapse occurred in sections between chainage K 15 + 160 and K 15 + 172 (Fig. 4) and although it was dealt with properly, another collapse occurred in this area due to aftershocks. To prevent the undesirable deformation of the tunnel, the tri-bench

Fig. 3. Hydrostatic disintegration test of phyllite: (a) sample; (b) sample immersed in water for about 1 min. 78

Tunnelling and Underground Space Technology 90 (2019) 76–98

G. Xu, et al.

and weathered secicite phyllite

Mileage K15+024 250 Chaniage Buried depth/m

200 150

Yao du

weathered secicite phyllite

Chaniage Mileage K16+910

fault fracture zone

100 50 0 -50

K15+450

K15+650

K15+900

K16+250

K16+550

Chaniage

Fig. 4. Longitudinal profile of the Dujia mountain tunnel.

were mainly distributed from chainage K 15 + 270 to K 15 + 450 and from chainage K 15 + 490 to K 15 + 770. (2) For sections in which there were more cracks, the widths of the cracks generally were greater than 0.5 mm, and some of them were larger than 2 mm. For other sections, the widths generally were less than 0.5 mm. (3) Six cracks were selected randomly to determine their depths, and Fig. 9 and Table 2 show the results of the drilling core. It was concluded that larger widths of the cracks were indicative of greater depths. The widths and depths of the cracks were the smallest at chainage K 15 + 508 and were larger at chainage K 15 + 693 and K 15 + 735. 2.3. Factors that cause cracks The geological radar test (Fig. 10) was conducted to detect the concrete reinforcement and the thicknesses of linings, as well as the distances between the steel arches and the defective areas behind linings. In addition, the bounce test was used to determine the strength of linings. The results showed that the defective areas (including voids and looseness) were sparse. The distance between the steel arches, the reinforcement of linings, and the strength of linings met the design requirements, indicating that construction defects were not the main reason causing cracking of linings. During the construction of the tunnel, concrete strain gauges were installed in 14 typical sections to analyze the inner forces of the lining. During the operation of the tunnel, data were transmitted to a central control system through automatic acquisition equipment and an optical fiber communication network, and these data are used in the evaluation of safety of linings. Fig. 11 shows the arrangement of the monitoring sections and the locations of the sensors. The monitoring data of each section had the same characteristics. Taking chainage K 16 + 230 as an example, the internal forces of the lining increased with the passage of time (Fig. 12), which indicated that

Fig. 5. Large deformation: (a) collapse; (b) shear failure of steel arches; (c) damage of preliminary supports; (d) cracking and leaking of secondary lining.

2.2. Cracking feature of linings Four years after the tunnel operation, many cracks appeared in the lining. Fig. 7 shows the distribution of cracks, and Fig. 8 shows the typical morphology of the cracks. The figures show that: (1) The cracks were mainly circumferential and longitudinal cracks that occurred on the haunches, vaults, and spandrels, and they were concentrated at certain sections. Specifically, cracks in T-1 were distributed in sections from chainage ZK 15 + 100 to ZK 15 + 350 and from chainage ZK 15 + 640 to ZK 15 + 800, and cracks in T-2

1

1

2

2 2m

3

3 2m

4

4 2m

(a)

(b)

Fig. 6. Supporting structures and construction method: (a) sketch map; (b) tunnel face. 79

Tunnelling and Underground Space Technology 90 (2019) 76–98

G. Xu, et al.

Table 1 Parameters of supporting structures. Advanced support

Double-layered small pipe with length of 5 m and diameter of 4.2 cm

Preliminary support

Secondary linings

Shotcrete/cm

Cable

Steel fabric

Steel arch

Concrete/cm

reinforcement

24

Fully grouted cable with diameter, length, axial and circumferential intervals of 2.2 cm, 4 m, 1 m and 0.6 m, respectively.

Diameter: 8 mm Interval: 2 cm

I 22b with interval of 40 cm

60

HPB 335 with diameter of 2 cm and interval of 40 cm

the crushed phyllite has a strong creep effect. Thus, it can be inferred that the time-dependent deformation of the surrounding rock was the main reason causing cracking of linings of this tunnel.

(1) Elastic stage (section AB): For a loading level is less than 9.6 MPa, the displacement remained stable while the internal force increased as the loading level increased. (2) Initial damage stage (section BC): The loading level at point C was 12.0 MPa. In this stage, both the deformation and force increased rapidly. Cracks appeared on the inner vault and left spandrel, and the lining spalled near the left arch springing. (3) Ultimate limit stage (section CD): cracks appeared at the right spandrel and left arch springing, and longitudinal cracks coalesced at the vault. In this stage, the growth rate of the deformation and the force accelerated. The maximum failure load was 14.2 MPa. (4) Instability stage (section after point D): the left arch springing was crushed, and the crack at the left and right spandrel coalesced through the lining. The structure had no load-bearing capability at this stage. The ultimate load was 15.2 MPa.

3. Similarity model test 3.1. Test device The test device (Fig. 13) was placed horizontally and consisted of a soil box, a foundation, horizontal jacks, vertical jacks, a cover plate, a reaction frame, and a control system. The soil box are 3.6 m long, 3.6 m wide, and 0.3 m high. The distances from center of the tunnel to the loading plate were equal to 2.5 times the equivalent diameter of the tunnel. The stress applied on each horizontal loading plate was controlled separately to allow anisotropic loading by horizontal jacks on each side. Four vertical jacks on the cover plate were used to maintain the plain strain state of test soil. In the test, we increased the stress on the loading plate gradually to approximately simulate the continuous increase of external loads on the lining due to the creep effect of surrounding rock, and the lateral pressure coefficient was chosen as 0.5.

4. Creep constitutive model and parameter inversion 4.1. Creep constitutive model (1) Creep constitutive model of surrounding rock

3.2. Similar materials and test equipment

For weak rock, when the external load is lower than a certain threshold, a recoverable viscoelastic deformation occurs, and, the growth rate of deformation decreases until it tends to a fixed value as time passes. This creep process is called the primary and secondary creep stage. When the external load exceeded the threshold, first, an irreversible deformation occurred with a constant growth rate; then, the rate accelerated during what is called the tertiary creep stage. Sterpi and Gioda (2009) found that the combination of the general Kelvin creep model and the Mohr-Coulomb mode can effectively describe this typical creep behavior of rock (Fig. 17). Thus, in this article, we used this model to describe the time-dependent behavior of the surrounding rock. The constitutive equation of the model for σ < σs is:

The similarity ratios of geometrical size and gravity were 25 and 1. According to the theories of scale model similitude (Gohl, 1991):

C =C =C =1

CE = C = Cc = 25

(1)

where Cν, Cε, Cφ, CE, Cσ, and Cc represent the similarity ratios of passion ratio, strain, inner friction angle, elastic modulus, stress, and cohesion, respectively. The surrounding rock was simulated by a mix of barite powders, fly ash, fine sands, and oils in certain proportions (Table 3). The bolt was simulated by aluminum wires with diameter of 2.5 mm, and its parameters were obtained based on the similarity ratio of tensile stiffness. The steel arch was modeled by aluminum wires with diameter of 4 mm, and its mechanical parameters were calculated according to the similarity ratio of bending stiffness. The lining and the shotcrete were simulated with a mixture of plaster, water, and diatomite in certain proportions to meet the similarity ratios of CE and Cσ. Reinforcement was modeled using iron material with a specific diameter according to the similarity ratio of bending stiffness. A lining that was 60 cm thick was selected to investigate the failure process. Five earth pressure cells were laid on the outer part of vault, spandrel, haunch, arch springing, and base slab. Strain gauges were installed on the inner and outer surface of the lining. The differential displacement gauges were set inside the tunnel, as shown in Fig. 14. The internal forces of lining were obtained based on strains of the inner and outer surfaces of lining described in Appendix I.

=

Ek

1

e

Ek k

t

+

(2)

Eh

where Eh is the elastic modulus of the Hoke element, Ek is the elastic modulus of the Kelvin element, and ηk is the viscosity coefficient of the Kelvin element. For σ ≥ σs, the equation is:

=

Ek

1

e

Ek k

t

+

Eh

+

p

(3)

where εp is the plastic strain. (2) Creep constitutive of anchorage zone Fig. 18 shows the axial force of the bolts in a typical section of the Dujia mountain tunnel. It was determined that the bolts bear small tensile forces. This was especially the case for the bolt at the haunch, where the maximum value was only 4.5 kN, which is far less than the designed value. The reason is that the relaxation area of the surrounding rock was larger than the length of the bolts, causing the shear forces at the contact surfaces between the bolts and the rock mass to be

3.3. Test results Figs. 15 and 16 show the test results, and the figures show that the failure process of the lining can be divided into the following four stages: 80

Tunnelling and Underground Space Technology 90 (2019) 76–98

G. Xu, et al. ZK15+ 050

070

090

110

130

150

170

190 200

220

240

260

280

300

320

340

360

380

400

420

440

460

480

500

520

540

560

580

600

Left haunch Vault Right haunch ZK15+ 200

Left haunch Vault Right haunch ZK15+ 400

Left haunch Vault Right haunch width >2 mm

widith:0.5~2mm

width <0.5mm

(a) ZK15+ 600

620

640

660

680

700

720

740

760

780

820

840

860

880

900

920

940

960

980 ZK16+ 000

800

Left haunch Vault Right haunch ZK15+ 800

Left haunch

Fig. 8. Typical morphology of cracks: (a) circumferential crack; (b) crossed crack; (c) longitudinal crack; (d) oblique crack.

Vault Right haunch ZK16+ 000

020

040

060

100

080

120

160

140

180

200

viscoelastic deformation, the general Kelvin element and the bolt element (called the B-K model, H.J. Zhao et al., 2012; T.B. Zhao et al., 2012) were used to describe this zone’s viscoelastic behavior. When the anchorage zone undergoes plastic deformation, the improved MohrCoulomb model was used to describe its mechanical behavior as follows:

Left haunch Vault Right haunch width <0.5mm

widith:0.5~2mm

width >2 mm

(b) K15+ 030

050

070

110

090

130

150

170

190

210

230

Left haunch

fs =

Vault 270

290

310

330

350

370

390

410

430

450

470

490

510

530

550

570

590

610

630

Left haunch

c =c+

Vault Right haunch

width <0.5mm

width:0.5mm-2mm

650

670

690

710

730

750

770

790

850

870

890

910

930

950

970

990 ZK16+ 010 030

810

830

Left haunch Vault Right haunch K15+ 830

=

Left haunch Vault 050

070

090

110

130

150

170

190

210

+

1

1

Vault width <0.5mm

width:0.5mm-2mm

exp

where = , = modulus of the bolt. For σ ≥ σs,

230

Left haunch Right haunch

1

Eh + Eb Eh + Ek

Right haunch K16+ 030

d 0 s (1/2 + /180) sin(45 + /2) 4 3 sa sc

(5)

where d0 is the diameter of the bolts; σs is the tensile strength of the bolts; and sa, sc are the axial and circumferential intervals of the bolts, respectively; c, φ are the cohesion and internal friction angle of surrounding rock, respectively. The constitutive equation of the B-K model: For σ < σs,

width >2mm

(c) K15+ 630

(4)

1 + sin

250

Left haunch Vault Right haunch K15+ 430

+ 2c N

where N = 1 - sin , c′ is the cohesion of the anchorage zone, and ϕ′ is the internal friction angle of the anchorage zone. The internal friction angle of the anchorage zone was set equal to that of the surrounding rock (Yang and Zhang, 2003). The cohesion of the anchorage zone can be calculated as follows (Torres, 2009):

Right haunch K15+ 230

3N

1

=

width >2mm

(d)

1

+

1

1

t

Eh Eb Eh + Eb

exp

t

(6)

+ Eb ,

+

= Eh + Eb , and Eb is the elastic

p

(7)

(3) Implementation of the creep constitutive model of the anchorage zone into PFC 2D

Fig. 7. Distribution of cracks on secondary lining: (a) ZK 15 + 050 ∼ ZK 15 + 600; (b) ZK 15 + 600 ∼ ZK 16 + 200; (c) K 15 + 030 ∼ YK 15 + 630; (d) K 15 + 630 ∼ YK 16 + 230.

In the particle discrete element software PFC 2D, the macroscopic creep behavior of a material is determined by setting the creep contact between particles, and the plasticity of the material is determined by adjusting the plastic parameters of bonding between particles. Fig. 20a and b show the general Kelvin contact and the B-K contact, respectively. The implementation of the general Kelvin contact in PFC 2D can be referred to Wang’s research (Wang et al., 2009, Appendix II). In this article, the finite

small. Therefore, the bolts were not built explicitly in the numerical approach. Instead, their effects were considered by strengthening the mechanical parameters of the surrounding rock in the anchorage zone, as shown in Fig. 19. To be specific, when the anchorage zone undergoes 81

Tunnelling and Underground Space Technology 90 (2019) 76–98

G. Xu, et al.

Fig. 9. Depths of cracks.

consistent with the laboratory test results in terms of the failure morphology, stress-strain curve, and peak strength (Zhang et al., 2010, Appendix II). The BPM model was used to represent the bonding behavior between concrete particles. In BPM model, the micro-parameters include the size distribution of the particles, the radius multiplier ( ¯ ), the normal and shear strength of bonds ( ¯c and ¯c ), the normal and shear stiffness of particles (k n and k s ), the normal and shear stiffness of bonds (k¯ n and k¯ s ), and the friction coefficient between particles (μc). The BPM model was used to simulate the sprayed concrete-steel arch composite (Table 5), and its equivalent elastic modulus is (Torres and Diederichs, 2009):

Table 2 Geometric parameters of cracks. Drilling position

Width/mm

Depth/cm

Photo

K K K K K K

0.2 0.5 0.4 1.4 0.6 1.1

2.4 7.0 5.0 10.5 5.5 13.5

Fig. Fig. Fig. Fig. Fig. Fig.

15 + 508 15 + 596 15 + 631 15 + 693 15 + 715 15 + 735

9a 9b 9c 9d 9e 9f

difference scheme of the B-K contact was developed (see Appendix III) and implemented into PFC 2D. An example was presented to validate this proposed contact. In this example (Fig. 20c), two balls were put together with the overlapping dimension of 0.01 m. The position and rotation of both balls were constrained. The micro-parameters of the B-K contact are shown in Table 4. The relaxation of contact force between balls is shown in Fig. 20c. The figure shows that the contact force decreased exponentially as the time increased, and the theoretical results and calculated results were consistent with each other.

E = Ec +

Ag Eg Sc

(8)

where E, Ec and Eg are the elastic modulus of sprayed concrete-steel arch composite, sprayed concrete and steel arch, respectively; Ag and Sc are the cross-sectional area of steel arch and sprayed concrete, respectively. Specimen with width of 5 cm and height of 10 cm was used to calibrate the mechanical parameters of BPM (Fig. 21a). The typical fracture pattern is shown in Fig. 21b. Comparison between the stressstrain curves of secondary lining (C30 concrete) and sprayed concretesteel arch composite (C15 concrete) obtained from lab test and numerical simulation are displayed in Fig. 21c-d. It shows that the simulation results are well consistent with the laboratory test results in terms of the stress-strain curve and peak strength. Micro-parameters of BPM for secondary lining and sprayed concrete-steel arch composite are shown in Table 5. Fig. 22 shows the numerical model. The size of the model was 80 × 80 m with a discrete area of 30 × 30 m. The thicknesses of the sprayed concrete and the lining were 24 and 60 cm, respectively. The thickness of the anchored area was 4 m. In the discrete area, the average radius of the rock particles was 9 cm, and the total number of rock particles was 45,632. The average radius of the concrete particles was 10 mm, and the total number of concrete particles was 16,340.

4.2. Inversion of parameters (1) Numerical inversion model In the numerical model, the surrounding rock was simulated using a discrete-continuous coupling method (Indraratna et al., 2015). That is, discrete particles were used to simulate the rock mass within the excavation influence area, and continuous elements were used to simulate rock mass outside the influence area to improve the computational efficiency. The lining was characterized by discrete particles. Concrete is composed of aggregate, cement mortar, and various interfaces. It is assumed that the particle size follows a Gaussian distribution without taking into account the actual distribution of each component. The simulation results of this simplified approach were

Fig. 10. Field test: (a) geological radar detection; (b) bouncing test. 82

Tunnelling and Underground Space Technology 90 (2019) 76–98

G. Xu, et al.

Fig. 11. Structural health monitoring: (a) Layout of monitoring sections; (b) sensor placement of each section; (c) strain gauge.

Fig. 12. Internal forces of secondary lining of typical monitoring section (K 16 + 230): (a) bending moments; (b) axial forces.

(2) Inversion method and results The micro-parameters of the rock were determined as follows: (1) The plastic parameters of B-K contact were calibrated according to the field tests results (Xu et al., 2016). (2) For the viscoelastic parameters of the B-K contact, an intelligent algorithm based on the PFC + PSO (particle swarm optimization) method was used to perform the inversion based on the inner force of the lining shown in Fig. 12. The detailed description can be found in Yang et al. (2017). The short-term and creep micro-parameters of rock are shown in Tables 6 and 7. The mechanical parameters of finite differential meshes are obtained through the micro-parameters of surrounding rock (Table 8). Good agreement between the numerical simulation and the measured results was observed in terms of the inner force of the lining (Fig. 23), which indicated that the inversed parameters can be used for subsequent analysis.

Horizontal jacks Vertical jacks

Horizontal jacks

Fig. 13. Test device. Table 3 Mechanical parameters of surrounding rock. Type

Elastic modulus E0/GPa

Density γ/kN·m−3

Cohesion c/MPa

Internal friction angle φ/°

Prototype Model Similarity ratio

0.5–1 0.025 25

19.2–19.5 19.4 1

3–4.2 0.15 25

9–11 10.5 1

5. Long-term safety analysis of the lining The lining can be regarded as a three-fold, statically-indeterminate structure. The degree of indeterminacy becomes one when a plastic hinge appears, and then it becomes zero (indeterminate structure) when another plastic hinge appears. Furthermore, more plastic hinges will be generated, and the structure will become unstable if additional damage occurs. In this stage, the structure is in its ultimate bearing state 83

Tunnelling and Underground Space Technology 90 (2019) 76–98

G. Xu, et al.

Fig. 14. Test preparation: (a) direct shear test on similar material of surrounding rock; (b) uniaxial compression test on similar material of lining; (c) placing earth pressure cells; (d) Steel arches and bolts; (e) shotcrete; (f) arrangement of strain gauges; (g) arrangement of displacement gauges.

and may undergo large deformation at any time (Zhang, 2012). The subsequent analysis of the failure process of the lining was based on this theory, as shown in Fig. 24. Three geo-stress cases were studied in this calculation, i.e., case I, λ (the lateral pressure coefficient, λ = σy/ C

B

A

D

1800 A

vault inverted arch spandrel haunch arch springing

20 10

-10

D

900 600 300

0

2

4

6

8

σ y /MPa

10

12

14

0

16

0

2

4

6

8

σ y /MPa

(a) B

C

800 A

D

12

Moment/kN.m

8000

4000

B

vault inverted arch spandrel haunch arch springing

600

vault inverted arch spandrel haunch arch springing

12000

10

14

16

(b)

16000 A

Axial force/kN

C

vault inverted arch spandrel haunch arch springing

1200

0

-20

B

1500 Contact pressure/kPa

Displacement/mm

30

σx) = 1, σx = σy = 8 MPa; case II, λ = 2, σx = 4 MPa, σy = 8 MPa; case III, λ = 0.5, σx = 8 MPa, σy = 4 MPa.

400 200

C

D

12

14

0 -200 -400 -600

0

0

2

4

6

8 10 σ y/MPa

12

14

-800

16

(c)

0

2

4

6

8

σ y/MPa

10

(d)

Fig. 15. Experimental results: (a) displacement; (b) contact pressure; (c) axial force; (d) bending moment. 84

16

Tunnelling and Underground Space Technology 90 (2019) 76–98

G. Xu, et al.

Fig. 16. Failure pattern of lining: (a) sketch; (b) photos.

5.1. Intact lining

which ran through the lining, appeared at the left and right spandrels and at the left haunch. When the lining eventually lost its load-bearing capacity (point d), numerous macro-cracks appeared. For λ = 1, micro-cracks also were distributed sporadically after tunnel operation. At point b, micro-cracks gathered at the left arch springing and the right spandrel. As time passed, three macrocracks appeared, one at the right spandrel and one each at the left and right arch springing (point c). The widths of the three macrocracks increased as time passed, while no new macro-cracks appeared (point d). For λ = 0.5, micro-cracks were distributed sporadically before point a. At point b, micro-cracks appeared on left and right spandrel. At point c, three macro-cracks appeared, one at the vault, and one each at the left and right spandrels. At point d, the damage degree of the three macro-cracks increased, and new macro-cracks appeared on the left arch springing.

Fig. 25 shows the failure process of different geo-stress cases for lining that has no initial damage and is in good contact with the surrounding rock. It can be concluded from Fig. 25 that: (1) For different geo-stress cases, the deformations of linings were significantly different. When λ = 2, the haunch convergence was 5.24 cm and the vault displacement was 0.97 cm at failure; when λ = 1, these values were 3.84 and 3.02 cm, respectively; when λ = 0.5, these values were 2.52 and 5.19 cm, respectively. (2) The failure processes also were significantly different for different geo-stress cases. For λ = 2, micro-cracks were distributed sporadically in the lining after tunnel operation. As time increased, micro-cracks began to increase constantly, and they accumulated at left and right haunches (point b). At point c, three macro-cracks,

(t)

<

s

ek

(t)

M-C

t

(t) p

>

(t)

Eh

σ

s

ηk Ek

σf

ek

(t)

Generalized Kelvin model

eh

t

(a)

(b) Fig. 17. Creep feature of weak rock: (a). creep curve; b constitutive model.

85

M-C

σ

Tunnelling and Underground Space Technology 90 (2019) 76–98

G. Xu, et al.

2.4 2.0

5.6

1.3 3.2

25.5

1.1 3.7

14.3

15.2 7.9

0.4 0.8

5.2

3.0

2.5

2.4 1.3

4.5

(a)

0.5

(b)

Fig. 18. Distribution of axial forces of bolts: (a) installing rock bolts; b Axial forces (units: kN).

Eb

σ

It can be concluded that the failure process of the lining could be divided into four stages, i.e., the elastic stage (before point a); the initial damage stage (point a to point b); the ultimate bearing stage (point b to point c); and the instability stage (point c to point d). The result was similar to the experimental findings.

M-C

ηk

Eh

σ

σf

(3) For λ > 1, micro-cracks initially appeared at the left haunch, the right haunch, and the inside part of the spandrel as the result of the bending damage caused by large-eccentricity compression. Thus, it is suggested that the circumferential reinforcement of the upper sections of the haunch should be strengthened to improve the bending capacity of secondary lining in this case. For λ = 1, cracks initially appeared at the left arch springing as a result of the

Ek

M-C

B-K

Fig. 19. Creep model of the anchorage rock.

m1

m1

Kmn

Kmn

Ckn

Kkn

m1 m2

Kbn

Kks fs

Ckn

Kks

Kkn

m1

m2

m2

fs K ms

Kms

Cks

m2

Cks

(a)

Kbs

(b) 2.0

Contact force/MPa

1.8

Numerical simulation analytical solution

1.6

1.4

1.2

1.0

0

10

Time/d

20

30

(c)

Fig. 20. Model validation: (a) general Kelvin creep contact; (b) B-K creep contact; (c) computational and theoretical results of normal contact force of B-K creep contact. 86

Tunnelling and Underground Space Technology 90 (2019) 76–98

G. Xu, et al.

compression-shear damage. Thus, it is suggested that the hooping at the arch springing should be strengthened to improve the shear bearing capacity of the secondary lining. For λ < 1, the bending damage initially occurred at the spandrel, so it is recommended that the circumferential reinforcement be enhanced at and above the spandrel. (4) The failure process and the cracking position of the lining obtained by numerical simulation were similar to the results of the similarity model test. Preliminarily, it was proved that the overloading method can approximately simulate the mechanical behavior of the lining when exposed to the creep effect of the surrounding rock macroscopically. However, the load acting on the lining increased at an equal rate in overloading method, which hardly can reflect the three creep stages of rock.

Table 4 Micro-parameters of B-K contact. Kelvin element

Elastic element

Bolt

Ek/GPa

ηk/GPa·y

Eh/GPa

Eb/GPa

1e8

1e8

1e8

1e8

Table 5 Micro-parameters of BPM for secondary lining and sprayed concrete-steel arch composite. Rmin/mm

Rmax/Rmin

Ec/GPa

kn/ks

E¯c /GPa

8 k¯ n /k¯ s 2.5

1.5 ¯c /MPa

26(22) ¯c /MPa

2.5 μc

26(22) ρ/kg·m−3

18 ± 3.6(12 ± 3)

18 ± 3.6(12 ± 3)

0.5

Lateral pressure coefficient has great influence on the results. Thus, it is important to determine its value in real engineering practices. Insitu stress measurement (such as, stress relief method, stress recovery method and hydraulic fracturing method. Huang et al., 2014; Jiang et al., 2016; Moayed et al., 2012) is widely adopted to obtain accurate

2330

Notes: The value of parameters in the parentheses are those of sprayed concrete-steel arch composite. For other parameter, the values of secondary lining are the same as sprayed concrete-steel arch composite.

(a)

(b)

35

20

30 Simulation Lab test

Simulation Lab test

15

20

UCS/MPa

UCS/MPa

25

15 10

10

5

5 0

0

1

2

3

4

5

6

axial strain/10-3

0

0

1

2

3

4

5

axial strain/10-3

(c)

(d)

Fig. 21. Comparison between lab test and numerical simulation results: (a) specimen; (b) typical fracture pattern; (c) stress-strain curve of C30 concrete; (d) stressstrain curve of sprayed concrete-steel arch composite.

87

Tunnelling and Underground Space Technology 90 (2019) 76–98

G. Xu, et al.

Anchorage zone

arch

Sprayed concrete-steel arch composite

spandrel

hauch

Secondary lining

arch springing invert

Fig. 22. Numerical model.

5.2. The influence of pre-existing cracks

Table 6 Short-term micro-parameters of surrounding rock. Rmin/cm

Rmax/Rmin

Ec/GPa

kn/ks

E¯c /GPa

7 k¯ n /k¯ s 2.5

1.66 ¯c /MPa

0.5 ¯c /MPa

2.5 μc

0.5 ρ/kg·m−3

2 ± 0.4(2.3 ± 0.5)

3 ± 0.6(3.4 ± 0.7)

0.5

Pre-existing cracks are prone to form in the lining during construction. Pre-existing cracks were modeled using the smooth joint model (SJM) (Vallejos et al., 2016). In the SJM, the micro-parameters include the normal and shear stiffness of contact (k¯ n and k¯ s ), normal strength and cohesion of contact ( c and cb ), friction coefficient, and the friction angle between contact (μ and φb). The micro-parameters of preexisting cracks are shown in Table 9, and the numerical model is shown in Fig. 26. Taking λ = 2 as an example, we studied the failure processes of a lining with one pre-existing crack at the right spandrel; a lining with two pre-existing cracks, one at the right spandrel and one at the arch springing; and a lining with three pre-existing cracks, i.e., one each at the right spandrel, the haunch, and the arch springing. In this calculation, the damage degree (the ratio of the thickness of cracked section to that of lining) was 0.5. Fig. 27 shows the failure process of lining with pre-existing cracks. It can be concluded form Fig. 27 that:

2640

Notes: The value of parameters in the parentheses are those of anchorage zone. For other parameters, the values of non-anchorage rock and anchorage zone are the same. Table 7 Creep micro-parameters of surrounding rock. Kelvin element

Elastic element

Bolts

Ek/GPa

ηk/GPa·y

Eh/GPa

Eb/GPa

1

1.3

2.1

0.88

(1) The existence of pre-existing cracks had little influence on the deformation feature of the lining. For the three cases, the haunch convergences and vault displacements were all about 4.3 and 2 cm, respectively. (2) The pre-existing cracks had a great influence on the failure process of the lining. In the elastic stage, micro-cracks formed around the pre-existing cracks, which is different from a lining with no preexisting cracks. In the initial damage stage, micro-cracks penetrated through the cross-section of the lining along the pre-existing cracks. In the ultimate bearing stage, macro-cracks formed at the vault and the left and right spandrels for the lining with one pre-existing crack, or distributed at the left spandrel and around the two preexisting cracks for the lining with two or three pre-existing cracks. In the instability stage, the damage degree of the three macrocracks increased, and new macro-cracks appeared at the arch springing.

Table 8 Mechanical parameters of finite differential meshes. Kelvin element

Elastic element

Bulk modulus

Ek/GPa

ηk/GPa·y

Em/GPa

K/GPa

1

1.3

2.1

0.38

stress field information. However, only several discrete points can be conducted to measure in-situ stress for specific project, and these data can only reflect the local in-situ stress field near the measuring points. Thus, inversion method becomes the main way to determine the overall distribution of in-situ stress in engineering area. That is, a three-dimensional numerical model considering topographic, stratigraphic, geomorphological and lithological characteristics of engineering area is established firstly. Then, the geo-stress field in this area is calculated by using various inversion method (such as, multiple linear regression, genetic algorithm, particle swarm optimization and surrogate model accelerated random search algorithm. H.J. Zhao et al., 2012; T.B. Zhao et al., 2012; Zhang et al., 2016; Djurhuusa and Aadnøyab, 2003) based on the in-situ stress measure results of individual points.

The phenomena regarding influence of pre-existing cracks on failure pattern of lining is related to the mechanical properties of pre-existing cracks and the mechanical behavior of lining where pre-existing cracks are located, that is,

88

Tunnelling and Underground Space Technology 90 (2019) 76–98

G. Xu, et al.

Fig. 23. Inversion results: (a) bending moments; (b) axial forces.

5.3. The influence of voids Voids are likely to form behind linings, resulting in poor contact between the linings and the surrounding rock. Thus, we analyzed the failure process of a lining with a void behind the vault, the right spandrel, and the right haunch. Fig. 28a shows the geometric features of the void. It is a semicircle with angle α of 18° (taking void at vault as an example, α is the angel between 9° on both sides from the symmetrical axis of the tunnel section) and diameter of 0.5 m. Fig. 28b shows the numerical model of void. Taking lining with a void behind the vault as an example, Fig. 29a shows the evolutionary curves of vault settlement, haunch convergence, and micro-cracks versus time. Fig. 29b-d show the evolution of cracks when λ = 2, 1, and 0.5, respectively. The following conclusions were based on the figure:

rigid arch

two hinged arch

three hinged arch Fig. 24. Failure process of secondary lining.

(1) For cracked section with damage degree of 0.5, it has certain ability to restrain the creep deformation of surrounding rock for its relatively large stiffness. While compared to the complete section, the effective area of section bearing external force is small, which makes pre-existing cracks propagate under external loads easily. (2) It can be seen form Figs. 12 and 15 that, pre-existing crack at spandrel or arch springing is in the compressive zone of lining, while pre-existing crack at haunch is in the tensile zone of lining. Thus, for pre-existing crack at haunch, it propagated mainly due to the tensile breakage of contacts at the tip of crack. When it penetrated the section, unloading occurred in particles around the crack. Therefore, there were few micro-cracks in particles near the crack when the lining was damaged. While for pre-existing crack at spandrel or arch springing, it propagated mainly due to the compressive-shear breakage of contacts at the tip of crack. When it penetrated the section, the creep-load born by the cracked section was transferred from cracked section to the nearby complete section. Thus, there were many micro-cracks in particles near the crack when the lining was damaged.

(1) Because of the absence of a constraint near the void, the value of haunch convergence and the direction of vault deformation of lining were quite different from those with no void. For λ = 2 and λ = 1, the haunch convergences were 8.5 and 4.8 cm, respectively, whereas these values were 5.24 and 3.84 cm when there was no void. For λ = 0.5, the haunch convergences were roughly equal for linings with and without a void. The vault deformation developed towards the hollow direction for all cases. (2) The existence of void had a great influence on failure process of the lining. In the elastic stage and the initial damage stage, the evolutions of cracks were similar for the linings with and without a void. In the ultimate bearing stage, one of the three penetrating cracks occurred near the void for each geo-stress cases, which was different from the lining with no void. Specifically, for λ = 2, the locations of the penetrating cracks were the vault and the left and right haunches. For λ = 1, the locations were the vault and the left and right arch springing. For λ = 0.5, the locations were the vault, the left arch springing, and the right spandrel. Thus, it was concluded that a hinge would form below the void when the lining lost its load-bearing capacity, which was the most notable influence of

89

300

4 c

200

2 100

b

40 30 20 10

Displacement/cm

d

6 cracking rate/count·y-1

haunch convergence cumulative cracks cracking rate

vault settlement haunch convergence

cumulative cracks cracking rate

4

0

10

20 30 Time/year

40

0

50

0

50

d 400

40

300 c

2 b

a

0

500

200 100

Cumulative cracks

50

400

vault settlement

Cumulative cracks

Displacement/cm

6

30 20 10

cracking rate/count·y-1

Tunnelling and Underground Space Technology 90 (2019) 76–98

G. Xu, et al.

a

0

0

10

20 Time/year

0 30

0

(a)

(b)

(c)

(d) Fig. 25. Failure processes of secondary lining under different geo-stress cases: (a) evolutional curves of cracks; (b) λ = 2; (c) λ = 1; (d) λ = 0.5. (The figures from left to right in Fig. 25b–d corresponding to points a, b, c, d in Fig. 25 a and this is the same for Fig. 27, Fig. 29). Table 9 Micor-parameters of pre-existing crack. Parameter

Value

Normal stiffness, k¯n (GPa·m−1) Shear stiffness, k¯ s (GPa·m−1)

6

Friction coefficient, μ Dilation angle, ψ (°) Tensile strength, σc (MPa) Cohesion, cb (MPa) Friction angle, φb (°)

14

0.4 0 0 0 0

the vault void on lining failure.

Fig. 26. Secondary lining with pre-existing cracks.

Fig. 30 shows the failure pattern of a lining with a void behind the right spandrel and the right haunch. Figs. 25, 29, and 30 show that there are three types of impacts that a void can have on the failure pattern of a lining:

(1) No effect. For example, for a lining with λ = 0.5 and a void behind the right haunch, the order of the appearance and the positions of the macro-cracks were the same as those of a lining without a void when λ = 0.5.

90

Tunnelling and Underground Space Technology 90 (2019) 76–98

0

10

4 200 c

2 100

b

40 30 20 10

cumulative cracks cracking rate

4

10

20 Time/year

0

30

0

0

0

10

d

300

200

2 b

a

a

0

100

300

400

vault settlement haunch convergence

c

100

50 40 30 20 10

cracking rate/count·y-1

b

20

d

6

Displacement/cm

c

2

30

cumulative cracks cracking rate

50

cracking rate/count·y-1

200

40

400

vault settlement haunch convergence

Cumulative cracks

4

300

6

Displacement/cm

Displacement/cm

d

50

cracking rate/count·y-1

400

vault settlement

haunch convergence cumulative cracks cracking rate

Cumulative cracks

6

Cumulative cracks

G. Xu, et al.

a

20 Time/year

0

30

0

0

0

10

20 Time/year

0

30

0

(a)

(b)

(c)

(d) Fig. 27. The influence of pre-existing cracks on the failure pattern of lining: (a) evolutional curves of cracks; (b) lining with one pre-existing crack; (c) lining with two pre-existing cracks; (d b) lining with three pre-existing cracks.

α

Ο

α α

(a)

(b)

Fig. 28. A void behind lining: (a) Schematic diagram of void; (b) numerical model.

91

Tunnelling and Underground Space Technology 90 (2019) 76–98

c

2 b

100

0

20 10

4 200 c

0

0

30 20

cumulative cracks cracking rate

10

20 Time/year

0

30

0

-4

0

10

20 30 Time/year

0

40

(a)

50 40 30

200 c

0

100

b

10 0

300

2

a

a

a

-2

100

b

40

400

-2

0

10

20 Time/year

30

0

20 10

cracking rate/count·y-1

200

300

d

vault settlement haunch convergence

d

Displacement/cm

4

30

4

50

cracking rate/count·y-1

6

40

400

vault settlement haunch convergence cumulative cracks cracking rate

Cumulative cracks

300

Displacement/cm

d

8

50

cracking rate/count·y-1

cumulative cracks cracking rate

8 Displacement/cm

400

vault settlement haunch convergence

Cumulative cracks

10

Cumulative cracks

G. Xu, et al.

0

(b)

(c)

(d) Fig. 29. The influence of a void behind vault on the failure pattern of lining: (a) evolutional curves of cracks; (b) λ = 2; (c) λ = 1; (d) λ = 0.5.

(2) A void can affect order in which macro-cracks appear but not affect the final failure pattern of the lining. Taking λ = 1 as an example, for a lining without void, the order of the appearance of macrocracks was right spandrel, followed by left and right arch springing. For a lining with a void behind the right spandrel, micro-cracks appeared earlier at the right arch springing than at the left arch springing. (3) A void can affect both the order of appearance of macro-cracks and the failure pattern of the lining. Taking λ = 1 as an example, for lining with a void behind the right haunch, macro-cracks appeared near the void first and then were distributed at the left arch springing, the right haunch, and the base slab. However, for a lining with no void, the macro-cracks appeared at the right spandrel first and then were distributed at the right spandrel, the left and right

arch springing. 6. Conclusion The similarity model test and a discrete-finite difference coupling method were used to study the progressive failure process of secondary lining. The conclusions are as follows: (1) The field investigation indicated that the time-dependent behavior of broken phyllite was the most important factor that caused the cracking of secondary lining of the Dujia mountain tunnel. (2) The similarity model test and numerical simulation results showed that the failure process of secondary lining can be divided into four stages, i.e., the elastic stage, initial damage stage, bearing capacity

92

Tunnelling and Underground Space Technology 90 (2019) 76–98

G. Xu, et al.

③ ②

③

②

①

(3) For λ > 1, initially, cracks appeared at the left haunch, right haunch, and inside part of the spandrel as the result of bending damage caused by large-eccentricity compression. Thus, it is suggested that the circumferential reinforcement of the upper sections of the haunch should be strengthened to improve the bending capacity of the secondary lining in this case. For λ = 1, cracks first appeared at the left arch springing, which was the result of compression-shear damage. Thus, it is suggested that the hooping at the arch springing should be strengthened to improve the shear bearing capacity of the secondary lining. For λ < 1, the bending damage first occurred at the spandrel, so it is recommended that the circumferential reinforcement be enhanced at and above the spandrel. (4) Pre-existing cracks had significant influence on the ultimate failure pattern of secondary lining. In each case, macro-cracks appeared near the pre-existing cracks first, and one or two hinges appeared near the pre-existing cracks when the lining lost its load-bearing capacity. (5) There are three types of impacts that a void can have on the failure pattern of lining, i.e., (1) no effect; (2) an effect on the order in which cracks appear while not affecting the final fracture pattern of lining; and (3) effects on the order in which cracks appear and the final fracture pattern of lining.

①

④

① ② ③

④

③

②

①

③

②

②

③

①

①

④

Acknowledgements

④

This research was supported by the National key research and development program of China (Grant No. 2016YFC0802210). Fig. 30. Influence of a void behind right spandrel and haunch on the failure pattern of lining ( =2, 1, 0.5 from top to bottom and the number ①-④ in the figure are the occurring order of cracks.)

Conflict of Interest The authors declare that they have no conflicts of interest.

stage, and instability stage. Appendix I

The inner force of lining in model test is calculate using the following method (PRC, 2010.). The secondary lining is in the state of small eccentric compression, and the state can be divided into three categories according to the value of ei:

N

e

e'

N

e

ei

α1σc

as'

As'

As

b

As'

As

As'

x

σs'As'

σsAs

σs'As'

σsAs

b

b h0 - as' h

α1σc

α1σc

σs'As'

As

as

e'

ei

ei

σs As

N

e

e'

as

h0- as' h

(a)

x

as'

as

(b) Fig. I.1. The different state of small eccentric compression: (a) state I; (b) state II; (c) state III. 93

h0- as' h

(c)

as'

Tunnelling and Underground Space Technology 90 (2019) 76–98

G. Xu, et al.

(1) State I: Reinforcement As is beyond the equivalent compressive section and it is in tension (Fig. I.1a). (2) State II: Reinforcement As is in the equivalent compressive section and it is in compression (Fig. I.1b). (3) State III: The whole section is in compression (Fig. I.1c). It can be obtained according to the equilibrium of force and moment:

N=

1 c bx

M=

1 c bx

+

s As

x + 2

h0

(I-1)

s As s As

(h 0

as )

(I-2)

where α1 is the coefficient of equivalent rectangular stress diagram, α1 = 1, b is the width of lining, x is the thickness of lining under compression, h0 is the effective thickness of lining, As, AS′ are the area of reinforcement, as′ is the distance from the center of reinforcement to the edge of section, σc, σs′ and σs are the stress of concrete, reinforcement AS′ and As. The following assumptions were made: a. The strain of reinforcement is equal to the strain of concrete at the edge of section near the reinforcement. b. x is described as:

x=

c

(

c

+ | t |)

h0

(I-3)

where εc is the strain at the edge of compressed section, εt is the strain at the edge of tensile or compressed section. Thus,

N=

M=

1 Ec0 c bh 0

1 Ec0 c bh 0

c

(

c

+ | t |)

(

c

+ | t |)

c

+ Er c As

(h 0

h0

Er t As c

2(

c

+ | t |)

(I-4)

) + Er c As (h 0

as )

(I-5)

where Eco is the elastic modulus of lining, Er is the elastic modulus of reinforcement. Appendix II The literatures, Zhang et al. (2010) and Wang et al. (2009), are important to support some conclusions in this manuscript. Thus, the main contents of these two papers are provided by translating them. For the paper of Zhang et al. (2010): In the numerical model, the BPM was used to represent the bonding behavior between concrete particles. The particle size follows a Gaussian distribution without taking into account the actual distribution of each component. The comparison between simulation results of this simplified approach and laboratory test results in terms of peak strength, stress-strain curve and the failure morphology are displayed in Table II.1 and Fig. II.1. Table II.1 shows that, for cubic or cuboid specimen, the peak strength of simulation result agrees well with that of the lab test result. The experimental stress-strain curve of cubic or cuboid specimen is in good agreement with the simulated curve, and the curves can be divided into three typical stages as follows: (1) Elastic stage (OA): the specimen deformed elastically, i.e., the relationship between strain and stress was linear. (2) Plastic stage (AB): the stress underwent non-linear growth with the increase of strain, and the slope of curve was smaller than that in elastic stage. (3) Strain softening stage (BC): After the specimen reached its peak strength, the stress tended to decrease to its residual strength, while the strain continued to increase. The failure morphology of cubic specimen obtained through numerical simulation or lab test shows that crack first appeared in the interior of the specimen, and then it propagated to the corner of the specimen near the loading surface. The failure morphology of cuboid specimen obtained through numerical simulation or lab test shows that the inclination angle of the main crack was about 60 degrees when the specimen broke. For the paper of Wang et al. (2009). The implementation of the general Kelvin contact in PFC 2D is conducted and verified by Wang et al. (2009). The derived expression of general Kelvin contact is:

f t+1 =

2Ck Km + K m2 t t + 1 u 2Ck + 2Km t

ut +

2K k t 2Ck ukt + ft 2Ck + Kk t 2Ck Km + K m2 t

(II-1)

Table II.1 Peak strength. Category of specimen

Dimension of specimen/(cm × cm)

Peak strength/MPa Lab test

Simulation

Cubic concrete Cuboid concrete

15 × 15 15 × 30

24.1 26.67

25.8 27.6

94

Tunnelling and Underground Space Technology 90 (2019) 76–98

G. Xu, et al.

30 25 20 UCS/MPa

Lab test Simulation

B

15

A

10

C

5 0

O 0

1

2

3

4

5

6

axial strain/10-3

(a) 30 25

Lab test Simulation

UCS/MPa

20 15 10 5 0

0

1

2

3

4

5

6

axial strain/10-3

(b) Fig. II.1. Lab test and simulation result: (a) Cubic concrete specimen; (b) Cuboid concrete specimen (For each subfigure, from left to right, stress-strain curve, simulated fracture pattern, and experimental fracture pattern.

where u is the displacement of contact, f is the contact force, t is the length of time step, t is the value of current time step, t + 1 is the value of next time step, CK is the viscosity of Kelvin element, Kk is elastic modulus of Kelvin element, Km is elastic modulus of Elastic element. An example was presented to validate this proposed contact. In this example (Fig. II.2), two balls were put together with the overlapping dimension of 0.05 m. The position and rotation of both balls were constrained. The micro-parameters of the contact are: CK = 1e7 Pa·d, Kk = Km = 1e7 Pa. The relaxation of contact force between balls is shown in Fig. II.2. The figure shows that the contact force decreased exponentially as the time increased, and the theoretical results and calculated results were consistent with each other.

Fig. II.2. Normal contact force calculated from PFC result and analytical result.

95

Tunnelling and Underground Space Technology 90 (2019) 76–98

G. Xu, et al.

Appendix III (1) Constitutive Equation For B-K model (see Fig. III.1):

Eb

σ

Bolt element

ηk

Eh

σ

Ek Generalized Kelvin element

Fig. III.1. B-K model.

=

=

b

=

b

(III-1)

r

+

(III-2)

r

where ε, εb, εr are total strain, strain of bolt element and strain of generalized kelvin element; σ, σb, σr are total stress, stress of bolt element and stress of generalized kelvin element. The generalized kelvin element meets the following formulation: h k r r

= Eh h = Ek k + = h= k = h+ k

k k

(III-3)

where σk, εk, Ek, ηk are stress, strain, elastic modulus and viscosity coefficient of kelvin element, respectively. σh, εh, Eh are stress, strain, elastic modulus of Hoke element, respectively.. The constitutive equation is:

+ p1

r

= p2

r

r

+ p3

(III-4)

r

where p1 = E + E ; p2 = h k For bolt element, k

b

= Eb

Ek Eh ; Eh + Ek

p3 =

k Eh

Eh + Ek

. (III-5)

b

where Eb is elastic modulus of bolt. The constitutive equation of B-K model is: (III-6)

+ p1 = (p2 + Eb ) + (p3 + p1 Eb) For

0,

=

= 0 . Thus, the creep equation is:

+ p1 = (p2 + Eb )

(III-7)

0

The initial relation is: t|t = 0

= (Eb + Eh )

(III-8)

0

Thus, the creep relation of B-K model is:

=

0 [(

where

)e

1 p1 t

= Eb + Eh ,

+ ]

(III-9)

= p2 + Eb

(2) Finite difference scheme The B-K contact model is developed and applied into PFC 2D. In order to program the creep model, its two-dimensional finite difference scheme is derived. For Kelvin element,

uk =

Ek uk ± F1 (III-10)

k

where uk is the displacement of Kelvin element, F1 is the contact force, ± denotes normal and tangential direction. The central finite difference of Eq. (III-10) is:

96

Tunnelling and Underground Space Technology 90 (2019) 76–98

G. Xu, et al.

ukt+ 1

ukt t

Ek (ukt+ 1 + ukt ) F t + 1 + F1t ± 1 2 2

1

=

k

(III-11)

where t is the length of time step, t is the value of current time step, and t + 1 is the value of next time step. Thus, uk can be obtained:

1 t t+1 (F1 + F1t ) Bukt ± A 2 k

ukt+ 1 =

(III-12)

where

A=1+

Ek t 2 k

B=1

Ek t 2 k

(III-13)

For Hoke element,

F1 = ±Eh uh

(III-14)

F1 = ±Eh uh

(III-15)

The central finite difference of the Hoke element is:

±

F1t + 1

F1t

= ±Eh

t

uht+ 1

uht

(III-16)

t

Thus, uh can be obtained:

uht+ 1 = ±

F1t + 1 F1t + uht Eh

(III-17)

Total displacement of generalized Kelvin model is: (III-18)

u = uk + uh The central finite difference of Eqs. (III-18) is:

ut = ukt+ 1

ut + 1

ukt + uht+ 1

(III-19)

uht

Thus, the force of generalized kelvin contact is:

F1t + 1 = ±

2 2

k Eh k

+ Ek Eh t

+ (Eh + Ek ) t

[ut + 1

ut +

2Ek tukt 2 k + Ek t

2

+ (Ek

k

2

k Eh

Eh ) t

+ Ek Eh t

F1t ]

(III-20)

For the bolt element:

F2 = ± Eb ub

(III-21)

F2 = ± Eb ub

(III-22)

For B-K model:

F2 = F

(III-23)

F1 ub = u

Thus:

F1t + 1 F1t

= ±F t+1

= ±F t

(III-24)

Eb ut + 1

((III-25)

Eb ut

The force of B-K contact is:

F t+1 = ±

2

k

1 + (Eh + Ek ) t

[2

k (Eh

- [2

+ Eb) + (Ek Eh + Eb Eh + Eb Ek ) t ] ut + 1

k (Eh + Eb ) + (Ek Eh

Eb Eh + Ek Eb) t ] ut

+ 2Ek Eh tukt + 2

(Ek

k

Ek ) tF t

(III-26)

Borja, R.I., Liu, F., 2010. Finite deformation formulation for embedded frictional crack with the extended finite element method. Internat. J. Numer. Methods Engrg. 82, 773–804. Caleyron, F., Combescure, A., Faucher, V., Potapov, S., 2012. Dynamic simulation of damage-fracture transition in smoothed particles hydrodynamics shells. Int. J. Numer. Methods Eng. 90, 707–738. Chiu, Y.C., Lee, C.H., Wang, T.T., 2017. Lining crack evolution of an operational tunnel influenced by slope instability. Tunn. Undergr. Space Technol. 65, 167–178. Djurhuusa, J., Aadnøyab, B.S., 2003. In situ stress state from inversion of fracturing data from oil wells and borehole image logs. J. Petrol. Sci. Eng. 38, 121–130. Fang, Y., Guo, J.N., Grasmick, J., et al., 2016. The effect of external water pressure on the

References Armero, F., Garikipati, K., 1996. An analysis of strong discontinuities in multiplicative finite strain plasticity and their relation with the numerical simulation of strain localization in solids. Int. J. Solids Struct. 33, 2863–2885. Asakura, T., Kojima, T., 2003. Tunnel maintenance in Japan. Tunn. Undergr. Space Technol. 18, 161–169. Bian, K., Liu, J., Xiao, M., et al., 2016. Cause investigation and verification of lining cracking of bifurcation tunnel at Huizhou Pumped Storage Power Station. Tunn. Undergr. Space Technol. 54, 123–134.

97

Tunnelling and Underground Space Technology 90 (2019) 76–98

G. Xu, et al. liner behavior of large cross-section tunnels. Tunn. Undergr. Space Technol. 65, 167–178. Foster, C.D., Borja, R.I., Regueiro, R.A., 2007. Embedded strong discontinuity finite elements for fractured geomaterials with variable friction. Int. J. Numer. Meth. Engng 72, 549–581. David, L.N., Amorim, D.F., Sergio, P.B., et al., 2014. Simplified modeling of cracking in concrete: Application in tunnel linings. Eng. Struct. 70, 23–35. Huang, H.W., Liu, D.J., Xue, Y.D., et al., 2013. Numerical analysis of cracking of tunnel linings based on extended finite element. Chinese J. Geotech. Eng. 35, 266–275 (in Chinese). Gao, S.M., Sun, Y., Wang, W., et al., 2017. Propagation of cracks in the secondary lining of tunnels subjected to asymmetrical pressure and a safety evaluation. Indian Geotech J. 47, 1–12. Gohl, W., 1991. Response of Pile Foundations to Simulated Earthquake Loading: Experimental and Analytical Results Ph.D. Dissertation. University of British Columbia. Guan, Z.C., Jiang, Y.J., Tanabashi, Y., et al., 2008. A new rheological model and its application in mountain tunnelling. Tunn. Undergr. Space Technol. 23, 292–299. Guan, Z.C., Jiang, Y.J., Tanabashi, Y., et al., 2009. Rheological parameter estimation for the prediction of long-term deformations in conventional tunneling. Tunn. Undergr. Space Technol. 24, 250–259. Hansbo, A., Hansbo, P.A., 2004. Finite element method for the simulation of strong and weak discontinuities in solid mechanics. Comput. Methods Appl. Mech. Engrg. 193, 3523–3540. He, C., Se, J., 2006. Maintenance and Reinforcement of Highway Tunnel. China Communications Press, Peking (in Chinese). Huang, M.Q., Wu, A.X., Wang, Y.M., Han, B., 2014. Geostress measurements near fault areas using borehole stress-relief method. Trans. Nonferrous Met. Soc. China 24, 3660–3665. Inokuma, A., Inano, S., 1996. Road tunnels in Japan: deterioration and countermeasures. Tunn. Undergr. Space Technol. 11, 305–309. Indraratna, B., Ngo, N.T., Rujikiatkamjorn, C., et al., 2015. Coupled discrete element–finite difference method for analysing the load-deformation behaviour of a single stone column in soft soil. Comput Geotech. 63, 267–278. Japan railway technical research institute, 2005. Countermeasure of Design Method to Deal With Tunnel Deformation. Japan Railway Technical Research Institute, Tokyo (in Japanese). Jiang, J.D., Liu, Q.S., Xu, J., 2016. Analytical investigation for stress measurement with the rheological stress recovery method in deep soft rock. Int. J. Min. Sci. Tech. 26, 1003–1009. Lee, D.H., Keun, H., Kim, Y.G., 2002. A study on mechanical behavior and cracking characteristic of tunnel lining by model experiment. Geosystem Eng. 5, 104–112. Lei, M.F., Peng, L.M., Shi, C.H., et al., 2015. Model test to investigate the failure mechanisms and lining stress characteristics of shallow buried tunnels under unsymmetrical loading. Tunn. Undergr. Space Technol. 46, 64–75. Meng, L.B., Li, T.B., Jiang, Y., et al., 2013. Characteristics and mechanisms of large deformation in the Zhegu mountain tunnel on the Sichuan-Tibet highway. Tunn Undergr Sp Tech. 37, 157–164. Miehe, C., Hofacker, M., Welschinger, F., 2010. A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits. Comput. Methods Appl. Mech. Engrg. 199, 2765–2778. Ministry of housing and urban-rural development of PRC, 2010. Code for Design of

Concrete Structure. China Communication Publisher Ltd., Beijing. Moayed, R.Z., Izadi, E., Fazlavi, M., 2012. In-situ stress measurements by hydraulic fracturing method at Gotvand Dam site. Iran. Turkish J. Eng. Env. Sci. 36, 179–194. Nguyen, V.P., Wu, J.Y., 2018. Modeling dynamic fracture of solids with a phase-field regularized cohesive zone model. Comput. Methods Appl. Mech. Engrg. 340, 1000–1022. Rabczuk, T., Song, J.H., Belytschko, T., 2009. Simulations of instability in dynamic fracture by the cracking particles method. Eng. Fract. Mech. 76, 730–741. Singh, B., Goel, R.K., Jethwa, J.L., et al., 1997. Support pressure assessment in arched underground openings through poor rock masses. Eng. Geol. 1997 (48), 59–81. Sterpi, D., Gioda, G., 2009. Visco-plastic behaviour around advancing tunnels in squeezing rock. Roc. Mech. Roc. Eng. 42, 319–339. Sandrone, F., Labiouse, V., 2011. Identification and analysis of Swiss National Road tunnels pathologies. Tunn. Undergr. Space Technol. 26, 374–390. Shi, Z.H., Ohtsu, M., Suzuki, M., et al., 2001. Numerical analysis of multiple cracks in concrete using the discrete element approach. Struct. Eng. 127, 1085–1091. Torres, C.C., 2009. Analytical and numerical study of the mechanics of rock bolt reinforcement around tunnels in rock masses. Rock Mech Roc Eng. 42, 175–228. Torres, C.C., Diederichs, M., 2009. Mechanical analysis of circular liners with particular reference to composite supports. For example, liners consisting of shotcrete and steel sets. Tunn. Undergr. Sp. Tech. 24, 506–532. Tonon, F., 2016. Sequential excavation, NATM and ADECO: What they have in common and how they differ. Tunn. Undergr. Space Technol. 25, 245–265. Vallejos, J.A., Suzuki, K., Brzovic, A., et al., 2016. Application of Synthetic Rock Mass modeling to veined core-size samples. Int J Rock Mech Min Sci. 81, 47–61. Wang, T.T., 2010. Characterizing crack patterns on tunnel linings associated with shear deformation induced by instability of neighboring slopes. Eng. Geol. 115, 80–95. Wang, T., LU, Q., Li, Y., et al., 2009. Development of contact model in particle discrete element method. Chinese J. Rock Mech. Eng. 28, 4040–4045 (in Chinese). Xiao, J.Z., Dai, F.C., Wei, Y.Q., et al., 2014. Cracking mechanism of secondary lining for a shallow and asymmetrically-loaded tunnel in loose deposits. Tunn. Undergr. Space Technol. 43, 232–240. Xu, G.W., He, C., Wang, Y., et al., 2016. Study on the safety performance of cracked secondary lining under action of rheological load. China Civ. Eng. J. 49, 114–123 (in Chinese). Yang, F.J., Zhang, C.Q., Zhou, H., et al., 2017. The long-term safety of a deeply buried soft rock tunnel lining under inside-to-outside seepage conditions. Tunn. Undergr. Sp. Tech. 67, 132–146. Yang, S.S., Zhang, B.S., 2003. The influence of bolt action force to the mechanical property of rocks. Rock Soil Mech. 24, 279–282 (in Chinese). Zhang, S.L., 2012. Study on health diagnosis and technical condition assessment for tunnel lining structure. Dissertation. Beijing Jiaotong University, China (in Chinese). Zhao, H.J., Ma, F.S., Xu, J.M., Guo, J., 2012. In situ stress field inversion and its application in mining-induced rock mass movement. Int. J. Rock Mech. Min. Sci. 53, 120–128. Zhao, T.B., Tan, Y.L., Liu, S.S., et al., 2012. Analysis of rheological properties and control mechanism of anchored rock. Rock Soil Mech. 33, 1730–1734 (in Chinese). Zhang, Z.J., Liu, J., Hu, W., et al., 2010. Two-dimensional simulation of concrete material fracturing by discrete element method. J. Hydroelectric Eng. 29, 22–27 (in Chinese). Zhang, S.R., Hu, A.K., Wang, C., 2016. Three-dimensional inversion analysis of an in situ stress field based on a two-stage optimization algorithm. J. Zhejiang Univ. Sci. A 17 (10), 782–802.

98