Peridynamics simulation of surrounding rock damage characteristics during tunnel excavation

Tunnelling and Underground Space Technology 97 (2020) 103289

Contents lists available at ScienceDirect

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

Peridynamics simulation of surrounding rock damage characteristics during tunnel excavation

T

⁎

Chenglu Gaoa, Zongqing Zhoua,b,c, , Zhuohui Lia, Liping Lia, Shuai Chenga a

Geotechnical and Structural Engineering Research Center, Shandong University, Jinan 250061, China Key Laboratory of Geotechnical Mechanics and Engineering of Ministry of Water Resources, Yangtze River Scientific Research Institute, Wuhan 430010, China c State Key Laboratory of Mining Disaster Prevention and Control Co-founded by Shandong Province and the Ministry of Science and Technology, Shandong University of Science and Technology, Qingdao 266590, China b

A R T I C LE I N FO

A B S T R A C T

Keywords: Peridynamics Material point dormancy method Excavation damage zone Tunnel excavation

The improved peridynamics is introduced to analyze the stability of surrounding rock in tunnel excavation process. By introducing short-range repulsive force into peridynamic equation of motion, the tensile and compression failure characteristics of rock materials is simulated. In order to simulate the process of tunnel excavation, the material point dormancy method is proposed in this paper. By simulating the distribution characteristics of excavation damage zone (EDZ) during the excavation of a circular tunnel with high stress difference, the relationship between the distribution position of EDZ and the direction of maximum principal stress is revealed. The v-shaped notches developed on the tunnel periphery at the direction of perpendicular to the maximum principal stress. The predicted EDZ distribution characteristics is consistent with the results of previous studies. And the displacement field of surrounding rock after tunnel excavation is also in good agreement with FEM simulation result. The simulation results show that this method not only has good stability, but also has high computational efficiency. By analyzing the deformation, damage and failure characteristics of surrounding rock under different buried depths, different lateral pressure coefficients, different excavation methods and different section shapes, the evolution process of rock instability caused by tunnel excavation and unloading is revealed. Which provides reference for the design and optimization of surrounding rock support scheme in deep buried tunnels.

1. Introduction With the increase in urbanization and development of economic construction all over the world, tunneling has become a preferred construction method for transportation and underground utility systems (Meguid et al., 2008). However, due to the complexity of geological conditions and the limitation of construction technology, geological disasters such as collapse, water inrush and large deformation often occur during tunnel construction, which seriously affect the construction safety and service life of tunnels (Li et al., 2013; Gao et al., 2018; Wang et al., 2019; Song et al., 2019). Tunnel excavation is a typical unloading process, which breaks the balance of the original stress filed and causes stress redistribution. The mechanical, hydraulic and geochemical properties of the surrounding rock will be altered during the process. The degree and extent of the excavation damaged zone (EDZ) is important for the design and construction of deep buried tunnels. EDZ is the zone around an excavation

⁎

where in-situ rock mass properties and conditions have been altered due to stress redistribution (Cai and Kaiser, 2005). EDZ can be mechanically unstable and also form a permeable pathway of groundwater flow which would raise the safety concern of the deep buried tunnels (Li et al., 2019). It is therefore important to characterize the EDZ for the design and construction of these deep buried tunnels. Numerical simulation is an important method to research the damage characteristics of surrounding rock during tunnel excavation. The continuum mechanics method represented by finite element method (FEM) and discontinuous mechanics method represented by discrete element method (DEM) are the two most important numerical simulation methods at present. The finite element method (FEM) (Sterpi and Cividini, 2004; Kasper and Meschke, 2004) or finite difference method (FDM) (Zhao et al., 2015; Senent and Jimenez, 2015) can simulate the deformation and stress distribution of surrounding rock and give the evolution characteristics of plastic zone. But it is difficult to accurately describe the damage state of rock mass due to their assumptions of

Corresponding author at: Geotechnical and Structural Engineering Research Center, Shandong University, Jinan 250061, China. E-mail address: [email protected] (Z. Zhou).

https://doi.org/10.1016/j.tust.2020.103289 Received 8 September 2019; Received in revised form 2 December 2019; Accepted 3 January 2020 0886-7798/ © 2020 Elsevier Ltd. All rights reserved.

Tunnelling and Underground Space Technology 97 (2020) 103289

C. Gao, et al.

dioxide blasting and analyzed the crack propagation mechanism. However, most of current researches are limited to the scale of laboratory test, focusing on the damage mechanism of rock samples. The numerical theory and simulation method of tunnel and other engineering problems based on peridynamics needs to be improved. The rapid development of peridynamics provides a new research idea for solving the stability of surrounding rock in tunnels and other underground engineering. In this study, the material point dormancy method is introduced to the bond-based peridynamics to study the damage characteristics of surrounding rock during tunnel excavation. The simulation results shows that the method has good stability and high efficiency in simulating the deformation and damage of surrounding rock in tunnel excavation process. The deformation and damage mechanism of surrounding rock under different buried depths, different lateral pressure coefficients, different excavation methods and different tunnel section shapes are analyzed. Which provides reference and guidance for the construction safety and design optimization of deep buried tunnels.

continuity and small deformation. Especially for jointed and fractured rock mass, the continuum mechanics method is not applicable because of the singularity of the solution at crack tips. The discrete element method (DEM) (Zhang et al., 2012; Jiang and Yin, 2012; Zhang et al., 2011a, 2011b) treats rock mass as discrete medium completely, and has certain advantages in describing the instability and large deformation of rock mass. But the low calculation efficiency makes it difficult to meet the demand of large-scale engineering calculation (Zhou et al., 2019). Meshless method, such as discontinuous deformation analysis (DDA) (Shi and Goodman, 1989; Do and Wu 2020), the numerical manifold method (NMM) (Shi 1991; He et al., 2017; Ma et al., 2009) and general particle dynamics (GPD) (Zhou et al., 2015; Bi et al., 2016) has been adopted by more and more researchers because of their advantage of not being restricted by grid in simulating material damage. These methods provide important reference for studying the distribution characteristics of EDZ during tunnel excavation. The continuousdiscontinuous method combines the advantages of the two methods in the simulation of rock mass deformation and damage, which has attracted extensive attention from researchers. For example, Zhang et al. (2016) carried out the numerical simulation to investigate the mechanical behavior of caverns during physical tests by a coupled DDAFEM method (NDDA). Morris et al. (2006) simulated the fracture and fragmentation of geologic materials using the combined FEM-DEM method to analyze the dynamic response of elaborate, underground facilities and tunnel systems to shock-wave loading. Wu et al. (2019) used a coupled FEM-DEM method to study the effect of confining pressure on the rock breakage efficiency and fragment size distribution of a TBM cutter. Meanwhile, some researchers have put some improved theories and methods, such as RFPA (Jia and Tang, 2008; Hu and Xu, 2011) and XFEM (Zhou and Chen, 2019; Rivas et al., 2019), to simulate the deformation and damage of rock materials. However, the numerical methods which including these special techniques and improvements still suffer from low efficiency and unsatisfactory simulation accuracy. The damage of rock mass is a process from continuous to discontinuous, the mathematical framework of conventional continuum mechanics is no longer suitable to modeling such problems. This is because that partial derivatives of the spatial coordinates are used to represent the relative displacement and force between the materials, but the partial derivatives are undefined along the discontinuities. Therefore, Silling (2000) proposes a non-local continuum mechanics method, which is known as peridynamics (PD). Compared with the conventional numerical methods, peridynamics replaces the partial differential equations with integral-differential equations to characterize the governing equations, which can solve discontinuities problems such as crack propagation (Silling and Askari, 2005). Because of its advantages in solving discontinuities problems, peridynamics is widely used in predicting the damage and fracture processes of materials, such as membranes and fibers (Silling and Bobaru, 2005), concrete structures (Gerstle et al., 2007), composite laminate (Sun and Huang, 2016), and polycrystalline materials (Meo et al., 2016), etc. Meanwhile, peridynamics shows high stability in solving the dynamic and static failure problems of materials without any external fracture criterion. For example, Ha and Bobaru (2010) studied the dynamic crack propagation and branching in brittle materials and compared the results with experimental results. Diehl et al. (2016) carried out a quantitative study in linear elastic fracture mechanics prediction about the critical traction before crack growth under quasi-static lading based on bond-based peridynamics. The wide and successful application of peridynamics in other fields has attracted the attention of geotechnical researchers. Zhou and Wang (2016) studied the crack propagation and coalescence in pre-cracked rock-like Brazilian disks using the non-ordinary state-based peridynamics. Ai et al. (2019) studied the dynamic mechanical properties and crack propagation law of rock under high strain rate impact loading by Split Hopkinson Pressure Bar (SHPB) tests and numerical simulation of peridynamics. Zhang et al. (2019) established the peridynamics model of rock fracturing under liquid carbon

2. Excavation damage zone of tunnel surrounding rock The rock damage is caused by the initiation, accumulation and growth of stress-induced cracks or fractures in rocks (Cai and Kaiser, 2005). The stability control of surrounding rock is one of the key problems that tunnel builders must consider. Tunnel excavation is the main factor causing the damage of surrounding rock. Therefore, to study the evolution law and distribution characteristics of EDZ, it is necessary to fully understand the effect of tunnel excavation on surrounding rock unloading. According to the research of Martino and Chandler (2004), the excavation-affected rock surrounding an excavation comprises an excavation damage zone and an excavation disturbed zone. In this paper the EDZ refers to irreversible damage to the rock resulting from the creation of an opening. The rock damage occurs from energy imparted to the rock from the excavation method and redistribution of in-situ stresses around the opening. Stress redistribution around openings causes zones of compressive or tensile stress, which may locally exceed the rock strength. This can lead to micro-or macrocracks in the rock surrounding the excavation and in more severe cases may result in spalling from the excavation surface. In severe cases of high stress concentration from stress redistribution, v-shaped notches may develop, as shown in Fig. 1. The slabbing failure occurs in the compressive region on the tunnel periphery at the direction of perpendicular to the maximum principal stress. In high horizontal in-situ stress environments, the EDZ would be expected to occur in the roof and floor of the tunnel. The EDZ distribution characteristics of surrounding rock is closely related to the magnitude and direction of in-situ stress in deep buried tunnels. Many researchers tried to reveal the distribution law of EDZ and its influence on the stability of tunnel surrounding rock through theoretical and numerical models. Schwartz and Einstein (1980) give an analytical solution based on the relative stiffness between the support and the surrounding ground. Li and Wang (2008) proposed an elastic plane strain solution for stresses and displacements around a lined circular tunnel in an isotropic geologic material due to uniform ground loads and internal pressure. Mortazavi and Molladavoodi (2012) developed a comprehensive damage model to include the material degradation due to compressive and tensile stress fields. Yang et al. (2013) deemed that the excavation process, time dependent behavior of rock mass and the evolution of fractured zones in rock mass are mainly three aspects contribute to the formation of the EDZ around a tunnel and developed an analysis model based on damage evolution mechanics. These developing theoretical models play an important role in understanding the mechanism of EDZ. Combined with damage and fracture mechanics and acoustic emission (AE) and microseismic (MS) measurement techniques, many researchers gave the prediction results of 2

Tunnelling and Underground Space Technology 97 (2020) 103289

C. Gao, et al.

(a) The obseved damage zone (Mortazavi and Molladavoodi, 2012).

(b) The development of the breakout notch Martino and Chandler (2004).

Fig. 1. Notch formed in the experiment tunnel.

acceleration of any material point at x in the reference configuration at time t is expressed as (Dayal and Bhattacharya, 2006; Zhou et al., 2018)

EDZ, as shown in Fig. 2. These EDZ prediction results of tunnel surrounding rock by various methods are consistent. The results suggest that the visible damage and cracking area is mostly concentrated in the place perpendicular to the maximum principal stress, which is resulted in the stress redistribution caused by the excavation unload effect. In the following part, peridynamics is used to study the evolution law and distribution characteristics of EDZ in tunnel excavation process, and the simulation results are expected to be consistent with the above research.

ρu¨ (x , t ) =

∫H

x

f (u (x ', t ) − u (x , t ), x ' − x ) dVx ' + b (x , t ),

(1)

where ρ is the mass density of x, u¨ is the secondary derivation of u, i.e. acceleration, Hx is a neighborhood of x, u is the displacement vector field, b is a prescribed body force density field, and f is a pairwise force vector per unit volume squared that x ' exerts on x. Within a limited distance, the material points interact with each other by the pairwise force vector f on the bonds. Thus, if the distance between the two material points exceeds a fixed distance, there is no interact between them. The limited distance δ , which is called the horizon (Fig. 3), as follows

3. Bond-based peridynamics framework 3.1. Peridynamics theory

Hx = {x ' ∈ R: ‖x ' − x‖ ≤ δ }.

There are two formulations of peridynamics: bond-based peridynamics and state-based peridynamics. And the state-based peridynamics can also be divided into ordinary and non-ordinary state-based peridynamics (Diana and Casolo, 2019; Silling et al., 2007; Warren et al., 2009). The bond-based peridynamics was firstly proposed by Silling (2000), which supposes that the discrete material points interact with each other by a pairwise force vector. The direct physical interaction between the material points is called a bond, or, in the special case of an elastic interaction to be defined a spring (Silling and Askari, 2005). In the bond-based peridynamics, the bonds are independent of each other, which results in the limitation of a fixed Poisson’s ratio (Le et al., 2014). In order to overcome the limitation, the state-based peridynamics was proposed by Silling et al. (2007). However, an instability problem, known as a “zero-energy model” is related to the state-based peridynamics (Yaghoobi and Chorzepa, 2017). The bondbased peridynamics has been certified that it can be used to simulate the fracture process of rock-like materials. For example, Ha et al. (2015) investigated the fracturing responses of a single flaw embedded in rocklike materials under compression. Lee et al. (2017) used peridynamic simulations to determine the extent of coalescing damage and identify the underlying causes. Zhou and Shou (2017) introduced tangent bonds into the bond-based peridynamic theory to investigate the initiation, propagation, and coalescence of cracks in brittle rock materials subjected to compressive loads. In this pater, the EDZ of surrounding rock in tunnel excavation process is studied by using the improved bondbased peridynamics. Because the bond-based peridynamics is the most simple and efficient method to simulate the deformation and damage of rock-like materials, and it meets the requirements of rock material calculation when the influence of Poisson’s ratio is ignored. The peridynamic equation of motion at a point in a homogeneous body can be expressed in the form of Newton’s Second Law and the

(2)

The pairwise force function f can be expressed as a function of the relative position vector ξ and the relative displacement vector η at time t (Kilic et al., 2009; Wang et al., 2019)

f (u (x ', t ) − u (x , t ), x ' − x ) = f (ξ , η)

(3)

where ξ is the relative position vector between the two material points in the reference configuration

ξ = x ' − x,

(4)

and η is the relative displacement vector between the two material points at time t

η = u (x ', t ) − u (x , t ),

(5)

and note that η + ξ represent the current relative position vector between the material points. The pairwise force function assures conservation of linear momentum

f (−η , −ξ ) = −f (η , ξ )

∀ η, ξ ,

(6)

and conservation of angular momentum

(η + ξ ) × f (η + ξ ) = 0

∀ η, ξ .

(7)

The equation means that the force vector between these material points satisfies Newton’s Third Law and is parallel to their current relative position vector and has the magnitude

f (η , ξ ) = f (η , ξ )

η+ξ , ‖η + ξ ‖

where f is a scalar function. 3

(8)

Tunnelling and Underground Space Technology 97 (2020) 103289

C. Gao, et al.

(a) MS model (Cai and Kaiser, 2005)

(b) Cohesion weakening-frictional strengthening (CWFS) model (Hajiabdolmajid and Kaiser, 2003)

(c) Finite element method model (Diederichs, 2007)

(d) Discrete element method model (Li et al., 2014)

(e) Damage model based on fracture mechanics (f) Damage evolution model (Yang et al., 2013) (Homand-Etienne et al., 1998) Fig. 2. Comparison of the EDZ prediction results of tunnels (Cai and Kaiser, 2005; Diederichs, 2007; Hajiabdolmajid and Kaiser, 2003; Homand-Etienne et al., 1998; Li et al., 2014; Yang et al., 2013).

A peridynamics material is referred to as being microelastic if the net work done on any material point x ' due to interaction with another fixed point x as x ' moves along any closed path is zero (Wang et al., 2018). For a microelastic situation, the pairwise force of the peridynamics material is defined as

f (η , ξ ) =

∂w (η , ξ ) ∂η

∀ η, ξ

s=

||η + ξ || − ||ξ || . ||ξ ||

(10)

However, when the deformation stretch exceeds the critical stretch s0 , the bond between the two material points breaks irreversibly and permanently. A scalar-valued function μ is thus introduced to describe the pairwise force of a prototype microelastic brittle (PMB) material (Silling and Askari, 2005)

(9)

f (η , ξ ) = μ (t , ξ ) cs,

where w is the micropotential which is the energy in a single bond and has dimensions of energy per unit volume squared. The bond stretch s is defined by Silling and Askari (2005)

(11)

where c is the micro-modulus denoting the elastic stiffness of a pairwise bond, and μ is defined as follows 4

Tunnelling and Underground Space Technology 97 (2020) 103289

C. Gao, et al.

Fig. 5. The relationship between force and stretch of bond in PMB material model (Silling and Askari, 2005).

Fig. 3. Horizon of point x (Sun and Huang, 2016).

1 if s (t , ξ ) < s0 μ (t , ξ ) = ⎧ ⎨ otherwise ⎩0

(12)

Therefore, local damage φ at a material point is introduced to illustrate the level of bond failure and is written as (Wang et al., 2018)

φ (x , t ) = 1 −

∫Hx μ (x , t , ξ ) dVξ . ∫Hx dVξ

(13)

Fig. 6. Schematic diagram of effects of short-range repulsive force.

Note that the local damage ranges from 0 to 1. When φ = 1, all the interactions initially associated with the point have been eliminated, while φ = 0, means that all interactions are intact. The value of local damage represents the possibility of internal crack formation within a body (Madenci and Oterkus, 2014). For example, initially a material point interacts with all materials in its horizon, as shown in Fig. 4(a); thus, the local damage has a value of zero. However, the creation of a crack terminates half of the interactions within its horizon, resulting in a local damage value of one-half, as shown in Fig. 4(b). With the propagation of the crack, all of the interactions are terminated within its horizon, and the local damage has a value of one in this case, as shown in Fig. 4(c).

2E Aδ2 9E ⎨ πhδ3 ⎪ 12E ⎩ πδ 4

⎧ ⎪ c=

for 2D for 3D

(14) '

in which A is the area of material point x , E is the elastic modulus and h is the thickness of the 2D disk. For the bond-based formulation, the critical stretch s 0 can be obtained from the critical energy release rate of the material, G0 as (Silling and Askari, 2005; Oterkus and Madenci, 2017)

s0 = 3.2. Peridynamics constitutive model

⎧ ⎪ ⎨ ⎪ ⎩

πG0 3κδ

for 2D

5G0 9κδ

for 3D

(15)

where κ is the bulk modulus of the material.

The simplest way to introduce failure into a constitutive model is by allowing bonds to break when they are stretched beyond a predefined limit. As shown in Fig. 5, after bond failure, there is no interaction force sustainable in the bond. And once a bond fails, it is failed forever. This makes the model has history dependent. According to the research of Oterkus et al. (2017), the material parameter c is specified as

(a) All interactions are intact with no damage

for 1D

3.3. Short-range repulsive forces In the preceding section, the material points interact only through bond forces, and the bond failure is judged by using the critical stretch. Therefore, the above method has a good effect in solving the tensile failure of materials, but is not good in solving the compression failure problems. This is because, the material points will be overlapped in the compression process, which is not consistent with the compression

(b) Half of the terminated interactions create a crack Fig. 4. Local damage. 5

(c) All interactions are terminated

Tunnelling and Underground Space Technology 97 (2020) 103289

C. Gao, et al.

Fig. 7. Diagrammatic of material point dormancy method.

ds = min {0.9‖x −x '‖, 1.35(rs + rs' )},

(17)

where rs is defined as the radius of material points. Note that shortrange forces are only repulsive, never attractive. Therefore, the peridynamic equation of motion describing the compression process of materials is established as

ρu¨ (x , t ) =

∫H

x

(f (u (x ', t ) − u (x , t ), x ' − x ) + fr (η , ξ )) dVx ' + b (x , t ). (18)

By improving the peridynamic equation of motion, the overlapping of material points in the compression process is limited, and the material points change from overlapped to incompressible, as shown in Fig. 6. 3.4. Material point dormancy method In order to simulate the excavation effect of underground space such as tunnel, the material point dormancy method is proposed. As shown in Fig. 7, the material points in the designated excavation area will be dormant to eliminate the interactions with the active material points during the tunnel excavation simulation. And at the same time, the interaction between the dormant material points in the excavation area will be eliminated too. However, the interaction between active material points in the normal region are unaffected. This procedure can be represented by introducing a scalar valued function as

Fig. 8. Geometry size and boundary conditions of circular tunnel model.

1 if x or x ' is in active ψ (x , x ', t ) = ⎧ , ' ⎨ ⎩ 0 if x or x is in dormant

(24)

Substituting from Eq. (24) into Eq. (11) results in an improved peridynamics pairwise force function which can be used to describe the excavation effect of materials, as

f (η , ξ ) = ψ (x , x ', t ) μ (t , ξ ) cs.

(25) '

That is, if the material points x or x is in the excavation area, the pairwise force between the two material points is set to zero. On the macro scale, all of the material points in the excavation area are cut off from the other areas to realize the material excavation simulation. The material points in the excavation area become dormant, but not actually nulled form the model matrix. This method has the advantages of invariant data structure, avoiding matrix recombination and ensuring the computational efficiency. In addition, the dormant material points have the possibility of recovery when simulating the process of material filling and generate. In the following, the above method will be used to simulate the distribution law of surrounding rock EDZ under the action of tunnel excavation and unloading.

Fig. 9. Displacement variation of key monitoring points of surrounding rock.

failure of rock materials. To prevent material points from overlapping, short-range repulsive forces fr are introduced into Eq. (1) as follows (Parks et al., 2008)

fr (η , ξ ) =

η+ξ cs min 0, (||η + ξ || − ds ) , ||η + ξ || δ

{

}

(16)

where 6

Tunnelling and Underground Space Technology 97 (2020) 103289

C. Gao, et al.

(a)

(b)

(c)

(d)

(f) Damage index

(e)

Fig. 10. Distribution characteristics of EDZ under different size of material points spacing.

numbers, Vj is the volume of the node j. In Eq. (19), the summation runs over all points within the neighborhood {||xj − xi || ≤ δ } , and the acceleration at the left side can be expressed as explicit central difference form

3.5. Numerical solution method The peridynamic equation of motion is an integral-differential equation, which is not usually amenable for analytical solutions. Therefore, its solution is constructed by using numerical techniques for spatial and time integrations (Foster et al., 2010; Liu et al., 2017). In a reference configuration, a body is uniformly discretized into nodes with a certain volume (Silling and Askari, 2005). Different from the finite element, all nodes form a grid allowing deformation without restriction. After discretization, the integral in Eq. (18) can be replaced by a finite sum, that is

u¨ in =

uin + 1 − 2uin + uin − 1 , Δt 2

where Δt is a constant time step that satisfies the stability condition given by Silling and Askari (2005)

Δt <

2ρ N

∑ j = 1 || x

c

j − xi ||

ρu¨ in =

∑ (f (ujn − uin, xj − xi) + fr (ηn , ξ n)) Vj + bin j

(20)

Vp

, (21)

where N is the total number of nodes inside the horizon of the node xi. After determining the acceleration of a material point at the nth time step from Eq. (21), the velocity and displacement at the next time step

(19)

where n is the time step number, the subscripts denote the node 7

Tunnelling and Underground Space Technology 97 (2020) 103289

C. Gao, et al.

Fig. 11. Distribution law of EDZ under different horizon size.

artificial damping to the system, the solution is guided to the steadystate solution as fast as possible. Therefore, the adaptive dynamic relaxation method (Kilic and Madenci, 2010) is adopted to solve the quasi-static problems in the paper.

can be obtained by employing explicit forward and backward difference techniques in two steps, respectively. The first step determines the velocity at the (n + 1)th time step using the known acceleration and the known velocity at the nth time step as

u̇ in + 1 = u¨ in Δt + u̇ in.

(22)

4. Simulation and prediction of EDZ during tunnel excavation

th

The second step determines the displacement at the (n + 1) time step using the velocity at the (n + 1)th time step from Eq. (22) and the known displacement at the nth time step as

uin + 1 = u̇ in + 1 Δt + uin.

4.1. Numerical model Assuming that a circular tunnel is excavated in a homogeneous rock stratum, where there are no discontinuities such as joints, cracks and bedding plane. The surrounding rock deforms freely without any support after tunnel excavation. Based on the above assumptions, a square numerical model with a side length of 40 m is established. The tunnel is located in the central of the square model with a radius of 5.5 m, as shown in Fig. 8. The numerical model is subjected to symmetrical loads,

(23)

Although the equation of motion of the peridynamics theory is in dynamic form, it can still be applicable to solve quasi-static or static problems by using a dynamic relaxation technique, which is based on the fact that the static solution is the steady-state part of the transient response of the solution (Kilic and Madenci, 2010). By introducing an 8

Tunnelling and Underground Space Technology 97 (2020) 103289

C. Gao, et al.

(a) Horizontal displacement (PD)

(b) Horizontal displacement (FEM)

(c) Vertical displacement (PD)

(d) Vertical displacement (FEM)

Fig. 12. Distribution law of displacement field.

of 2600 kg/m3. Based on the above mentioned prototype microelastic brittle (PBM) material model and adaptive dynamic relaxation solve method, the numerical study of EDZ during tunnel excavation is carried out. The simulation process is divided into two stages. Firstly, the initial in-situ stress balance is calculated. That is, the deformation and stress state of strata under the boundary conditions are calculated before tunnel excavation. Then the tunnel excavation simulation is carried out by means of material point dormancy method under the condition of stress balance. In the current simulation, 1000 iterations are calculated in the stage of stress balance, and another 1000 iterations are calculated in the stage of stress re-balance after tunnel excavation. Four key monitoring points are set at tunnel vault, bottom and left and right side walls. The vertical displacement of tunnel vault (signed as A) and bottom (signed as B) and the horizontal displacement of tunnel left side wall (signed as C) and right side wall (signed as D) are monitored respectively. Fig. 9 shows the variation of displacement of the key monitoring points with the calculation steps in the simulation process. It can be seen that the in-situ stress reaches the balance condition at the end of 1000th calculation step, and the displacement of surrounding rock is suddenly changed due to tunnel excavation, and reaches a new balance state again with iterative calculation. The deformation of surrounding rock represented by the key monitoring points shows a tendency of convergence into the tunnel inside, which is consistent with the actual situation. The influence of material points spacing and horizon size on surrounding rock EDZ will be investigated below.

Fig. 13. Geometry size and boundary conditions of horseshoe tunnel model.

with 6.9 MPa in vertical direction and 34.5 MPa in horizontal direction. A tunnel excavation model with a lateral pressure coefficient of 5.0 and a high stress difference is established. The stresses are applied separately on the virtual boundaries, which is composed of three layers of material points. The thickness of the numerical model is one layer of material points, so the model can be regarded as a plane strain problem. The physical and mechanical parameters of the numerical model are selected as elastic modulus of 18 GPa, Poisson’s ratio of 1/3 and density

4.2. The influence of material points spacing on EDZ In order to investigate the influence of material points spacing on the distribution characteristics of EDZ during tunnel excavation. Five 9

Tunnelling and Underground Space Technology 97 (2020) 103289

C. Gao, et al.

Fig. 14. Damage and deformation of surrounding rock under different buried depths (λ = 1.0).

EDZ in the five numerical models are similar. The EDZ is concentrated in the compressive region on the tunnel periphery at the direction of perpendicular to the maximum principal stress. This is consistent with the phenomenon observed in underground research laboratory and previous studies. But with increasing the number of material points in the calculation model, the distribution of EDZ becomes more concentrated and the contour becomes smoother. The damage index decreases with increasing the number of material points, and finally tends to be stable. This shows that the accuracy of the peridynamics model can be improved with decreasing the material points spacing, but it cannot be increased indefinitely, which is in agreement with δ-convergence of peridynamics theory. Considering the calculation cost and accuracy comprehensively, the material point spacing can be selected as Δ = L/400 , which meets the simulation requirements.

numerical models with material points spacing of Δ = L/200 , Δ = L/300 , Δ = L/400 , Δ = L/500 and Δ = L/600 are respectively simulated using the improved peridynamics and proposed material points dormancy method. The horizon size is specified as δ = 3Δ . In order to characterize the damage degree of surrounding rock in tunnel excavation process, the damage index ψ is defined as the ratio of local damage to the initial integrity of the surrounding rock, stated as follows:

ψ=

∫ φ (x , t )dVx ∫ dVx

(26)

The distribution characteristics of EDZ under different material points spacing is shown in Fig. 10(a)–(e). And the variation of damage index of surrounding rock under different material points spacing is shown in Fig. 10(f). It is obvious that the distribution characteristics of 10

Tunnelling and Underground Space Technology 97 (2020) 103289

C. Gao, et al.

Fig. 14. (continued)

Fig. 16. Damage indexes of surrounding rock under different buried depths (λ = 1.0).

Fig. 15. Displacement variation curve of key points under different buried depths (λ = 1.0).

observed in actual engineering. But when the horizon size is too small or too large, some calculation errors will occur in the simulation results. This is because the small horizon size cannot be able to capture crack propagation behavior. And the discrepancy between analytical and numerical solutions becomes larger when the horizon size increases due to the excessive wave dispersion (Silling and Askari, 2005; Madenci and Oterkus, 2014). The damage indexes between different horizon sizes are not much difference. According to the distribution characteristics of EDZ and the variation curve of damage index, the horizon size can be selected as δ = 3Δ or 4Δ , which meets the requirements of current numerical simulation.

4.3. The influence of horizon size on EDZ In order to investigate the character of m -convergence of the improved peridynamics, five numerical simulations with horizon size of δ = 2Δ , δ = 3Δ , δ = 4Δ , δ = 5Δ and δ = 6Δ are respectively carried out. The horizon size is specified as δ = 0. 3 m . The distribution characteristics of EDZ under different horizon sizes is shown in Fig. 11(a)–(e). And the variation of damage index of surrounding rock under different horizon size is shown in Fig. 11(f). The distribution positon and shape of EDZ under different horizon sizes are similar and also in good agreement with the phenomenon 11

Tunnelling and Underground Space Technology 97 (2020) 103289

C. Gao, et al.

Fig. 17. Damage and deformation of surrounding rock under different lateral pressure coefficient (h = 500 m).

4.4. Distribution characteristics of displacement field

rock deformation is considered. By comparing the distribution characteristics of horizontal and vertical displacement of two numerical simulation methods, it can be seen that PD and FEM are consistent in predicting the deformation trend of surrounding rock during tunnel excavation. Under the condition of high horizontal stress, the left and right side walls of the tunnel tend to converge into the inside of tunnel. The tunnel vault and bottom are under compression stress and tend to expand outwards. However, FEM cannot predict the damage of surrounding rock. Therefore, the simulation results of the two methods are different in the EDZ of tunnel vault and bottom. It is unique advantage of PD to predict the generation, expansion and distribution of EDZ during tunnel excavation. In this section, the influence of material points spacing and horizon size on EDZ distribution characteristics in tunnel excavation process has been investigated. As well as the simulation results between PD and

The deformation of surrounding rock during tunnel excavation can also be simulated by finite element method (FEM). As one of the most common and mature numerical simulation methods, FEM can predict the deformation characteristics of surrounding rock after tunnel excavation. Suppose a tunnel excavation model with the same size and parameters as the above PD simulations. The FEM simulation is also divided into two stages. After the in-situ stress balance, the tunnel excavation simulation is carried out. The large general software COMSOL Multiphysics is used for the FEM simulation. The distribution characteristics of the displacement field of PD and FEM simulation results is shown in Fig. 12. Where the displacement value is the calculated result after tunnel excavation minus the calculated result of in-situ stress balance. That is, only the influence of tunnel excavation on surrounding 12

Tunnelling and Underground Space Technology 97 (2020) 103289

C. Gao, et al.

Fig. 17. (continued)

Fig. 18. Displacement variation law of key points under different lateral pressure coefficients (h = 500 m).

Fig. 19. Damage index of surrounding rock under different lateral pressure coefficients (h = 500 m).

FEM has been contrastive analyzed. Which proved that the improved peridynamics and material point dormancy method are reliable and applicable in simulating tunnel excavation. In the next section, the evolution law of damage and deformation of surrounding rock during tunnel excavation under different conditions were discussed.

length of 40 m is established. The horseshoe-shaped tunnel is located in the central of the square model with a width of 11 m and a high of 8.5 m, as shown if Fig. 13. The numerical model is uniformly dispersed to 400 × 400 material points, and the horizon size is selected as δ = 3Δ. In order to simulate the in-situ stress state, the stress boundary conditions corresponding to the buried depth and lateral pressure coefficient of the tunnel are applied to the four sides of the model. The vertical in-situ stress is simplified to the dead weight of the overlying rock mass: σ1 = ρgh, where ρ is the density of rock mass, g is the gravitational acceleration and h is the depth of overlying rock mass. The horizontal in-situ stress is simplified to the product of vertical in-situ

5. Damage and deformation of surrounding rock under different conditions Assuming that a more common horseshoe tunnel is excavated in a homogeneous rock stratum. A square numerical model with a side 13

Tunnelling and Underground Space Technology 97 (2020) 103289

C. Gao, et al.

(a) Full cross-section method

(b) Two-bench method

(c) Three-bench method

Fig. 20. Schematic diagram of different tunnel excavation methods.

Fig. 21. Damage and deformation of surrounding rock under different tunnel excavation methods (h = 700 m, λ = 1.0).

14

Tunnelling and Underground Space Technology 97 (2020) 103289

C. Gao, et al.

Fig. 21. (continued)

lateral pressure coefficient is specified as 1.0. Fig. 14 shows the numerical simulation results of surrounding rock damage and deformation under different buried depths. It can be seen that the EDZ of surrounding rock is distributed differently around the tunnel under different buried depths. When the buried depth is small, the EDZ is mainly concentrated in the side walls. With the increase of buried depth, the EDZ is expanded to the vault and bottom of the tunnel and the damage degree also increased. By comparing the deformation of surrounding rock under different buried depths, it can be seen that the failure zone of surrounding rock keeps expanding with the increase of buried depth, and the micro-cracks in the failure zone around the tunnel gradually form an annular macro-crack. Four key monitoring points in tunnel vault, bottom and side walls are used to monitor the variation of displacement of surrounding rock with the increase of buried depth. The variation of vertical displacement of key points in tunnel vault (signed as A) and bottom (signed as B) and the horizontal displacement of key points in tunnel left side wall (signed as C) and right side wall (signed as D) is shown in Fig. 15. It is

stress and lateral pressure coefficient: σ2 = λσ1, where λ is the lateral pressure coefficient. The physical and mechanical parameters of the numerical model are selected as elastic modulus of 6.4 GPa, Poisson's ratio of 1/3 and density of 2600 kg/m3. The simulation process is divided into two steps. Firstly, the initial in-situ stress balance calculation is carried out, and then the tunnel excavation is simulated. A series of numerical simulation were carried out using the improved peridynamics. The damage and deformation characteristics of surrounding rock during tunnel excavation under different conditions, such as buried depth, lateral pressure coefficient, excavation method and tunnel section shape, were studied.

5.1. The influence of tunnel buried depth Firstly, the damage and deformation characteristics of surrounding rock during tunnel excavation under different buried depths are studied. It is assumed that there are five numerical models with buried depth of 300 m, 400 m, 500 m, 600 m and 700 m, respectively. And the 15

Tunnelling and Underground Space Technology 97 (2020) 103289

C. Gao, et al.

can be see that when lateral pressure coefficient changes from less than 1.0 to greater than 1.0, the damage degree of surrounding rock in tunnel excavation process increases substantially. This is mainly caused by the damage of surrounding rock in tunnel vault and bottom. Therefore, it is necessary to strengthen the support of tunnel vault and bottom when tunneling in the stratum with high lateral pressure coefficient. 5.3. The influence of tunnel excavation method Next, the variation of surrounding rock stability during tunnel excavation with different excavation methods are studied using the above numerical method. Assuming that a tunnel is excavated in a homogeneous rock stratum by three excavation methods: full cross-section method, two-bench method and three-bench method, as shown in Fig. 20. The buried depth of the tunnel is 700 m and the lateral pressure coefficient is 1.0. The physical and mechanical parameters of the model are same as mentioned above. The simulation results of damage and deformation of surrounding rock for each excavation step are shown in Fig. 21. As can be seen in Fig. 21, there is little difference in the damage and deformation characteristics of surrounding rock after the completion of tunnel excavation between different excavation methods. However, with the increase of excavation steps, the influence of tunnel excavation disturbance on surrounding rock decreases. The damage degree of surrounding rock decreases with the number of excavation steps, as shown in Fig. 22. Meanwhile, it can be seen that the damage degree of surrounding rock is low when the tunnel excavation section is small. But, with the increase of tunnel excavation section, the damage of surrounding rock increases obviously. This is because the effect of support was not considered in the numerical simulations. Therefore, timely support is necessary in the process of tunnel excavation. The simulation examples show that the proposed material point dormancy method can accurately describe the progressive influence of tunnel excavation on surrounding rock damage.

Fig. 22. Damage indexes of surrounding rock under different tunnel excavation methods (h = 700 m, λ = 1.0).

observed that, the displacement of the key monitoring points of the surrounding rock has an increasing trend with the increase of buried depth. The displacement growth rate of the key monitoring points increases gradually. Fig. 16 shows the variation of damage degree of surrounding rock in tunnel excavation process with the increase of buried depth. Obviously, the damage degree of surrounding rock increases with increasing the buried depth.

5.2. The influence of lateral pressure coefficient Secondly, the damage and deformation characteristics of surrounding rock under different lateral pressure coefficients in tunnel excavation process are studied. Five numerical models with lateral pressure coefficients of 0.50, 0.75, 1.00, 1.25 and 1.50 were selected. And the buried depth of tunnel in the numerical simulations is specified as 500 m. The simulation results of damage and deformation of surrounding rock under different lateral pressure coefficients are shown in Fig. 17. Obviously, when the lateral pressure coefficient is less than 1.0, the EDZ is mainly distributed in tunnel side walls. The EDZ in tunnel vault and bottom are developed and expanded gradually with the increase of the lateral pressure coefficient. When the lateral pressure coefficient is greater than 1.0, the EDZ is mainly concentrated in tunnel vault and bottom. This is because the EDZ is significantly affected by the magnitude and direction of in-situ stress, and is mainly distributed in the direction perpendicular to the maximum principal stress. The lateral pressure coefficient not only represents the ratio of horizontal and vertical in-situ stresses, but also controls the direction of the maximum principal stress. Therefore, when the lateral pressure coefficient changes from 0.5 to 1.5, the direction of the maximum principal stress changes from vertical to horizontal and the location of EDZ is also changed. By comparing the deformation characteristics of surrounding rock and the development trend of cracks under different lateral pressure coefficients, it can be seen that the micro-cracks in tunnel vault and bottom gradually increase and expand. And the large cracks are gradually formed in the failure zones with the increase of the lateral pressure coefficient. Meanwhile, it can be seen from Fig. 18 that with the increase of the lateral pressure coefficient, the vertical displacement of tunnel vault and bottom decreases firstly and then increases. And the horizontal displacement of side walls increases firstly, then decreases and then increases. In addition, the change law of surrounding rock damage index under different lateral pressure coefficients is shown in Fig. 19. It

5.4. The influence of tunnel section shape Finally, the stability of tunnel surrounding rock with different section shapes is studied with above numerical methods and models. In the current simulations, three kinds of tunnel section shapes are adopted: circular, horseshoe and rectangular. The buried depth of the tunnels is 700 m, and the lateral pressure coefficient is 1.0. The damage and deformation of surrounding rock under different tunnel section shapes is shown in Fig. 23. And the damage indexes of the three simulation models is shown in Fig. 24. As shown in Figs. 23 and 24, it can be seen that the tunnel excavation with different section shapes has a great influence on the stability of surrounding rock, and the distribution of EDZ is also very different. Due to the regular shape of the circular tunnel, the EDZ is bisymmetric and evenly distributed around the tunnel. Since the rectangular tunnel has four right angles, the EDZ is symmetrically distributed along the four sides of surrounding rock, and the damage degree is the most serious. The damage degree of the horseshoe tunnel is in the middle of the above two conditions, because its section profile is relatively smooth. The distribution characteristics of micro-cracks in tunnel surrounding rock is the same as the above EDZ. Obviously, the stability of surrounding rock in circular tunnel is the best, and the rectangular tunnel is the worst. Therefore, the reasonable design of tunnel section shape is important for the stability of surrounding rock in deep buried tunnels, and the stability of surrounding rock can be improved by optimizing the section shape of tunnels. 6. Discussions In this paper, the material point dormancy method, based on the 16

Tunnelling and Underground Space Technology 97 (2020) 103289

C. Gao, et al.

Fig. 23. Damage and deformation of surrounding rock under different tunnel section shapes (h = 700 m, λ = 1.0).

(2007) proposed state-based peridynamics and Zhou and Shou (2017) introduced tangent bonds into the bond-based peridynamics. But these computational techniques add complexity to the peridynamics. This paper is mainly devoted to solving the application problem of peridynamics in simulating the damage and deformation of surrounding rock in underground engineering. Therefore, the most simple and efficient bond-based peridynamics is chosen to carry out the research, and the results prove that it meets the requirements of rock material calculation when the influence of Poisson’s ratio is ignored. In addition, the numerical models in this paper adopt the isotropy and uniformity hypothesis, without considering the influence of complex geological structures such as joint surface and bedding plane. The gravity effect of the model is neglected, because the tunnel is buried at great depth and the damage evolution is mainly affected by stress redistribution and tunnel excavation.

improved bond-based peridynamics, is proposed to simulate the damage and deformation characteristics of surrounding rock during tunnel excavation. By comparing with previous experimental observation and numerical simulations, the EDZ distribution characteristics obtained by peridynamics in this paper is consistent with the previous research results. The predicted EDZ is mainly distributed on the tunnel periphery at the direction of perpendicular to the maximum principal stress. However, some assumptions and simplifications are used in the numerical models. Firstly, the rock materials simulated in this paper have fixed Poisson’s ratio due to the inherent limitation of bond-based peridynamics. The Poisson’s ratio is fixed at 1/3 in 2D model, and 1/4 in 3D model. Which limits the application for bond-based peridynamics in various complex rock materials. In order to overcome the restriction, researchers have made some useful efforts. For example, Silling et al. 17

Tunnelling and Underground Space Technology 97 (2020) 103289

C. Gao, et al.

(4) Different excavation methods and section shapes also have certain influence on the distribution of EDZ and deformation characteristics of surrounding rock. Multi-step excavation decreases the disturbance of tunneling to surrounding rock. Due to the high smoothness of the circular, the damage degree of surrounding rock of the circular tunnel is minimal. In order to optimize the design scheme and improve the stability of surrounding rock, these factors should be taken into account in deep strata excavation. CRediT authorship contribution statement Chenglu Gao: Writing - original draft. Zongqing Zhou: Conceptualization, Methodology, Writing - review & editing. Zhuohui Li: Formal analysis, Investigation. Liping Li: Resources, Supervision, Project administration, Funding acquisition. Shuai Cheng: Validation, Data curation. Declaration of Competing Interest Fig. 24. Damage indexes of surrounding rock under different tunnel section shapes (h = 700 m, λ = 1.0).

We declared that we have no financial and personal relationships with other people or organizations that can inappropriately influence our work. There is no professional or other personal interest of any nature or kind in any product, service and/or company that could be construed as influencing the position presented in, or the review of, the manuscript entitled Peridynamics Simulation of Surrounding Rock Damage Characteristics during Tunnel Excavation.

Although the present model is not an exact solution, it is helpful for the assessment of excavation disturbed or damage zone. It provides important reference for tunnel and other underground engineering construction safety. In the future work, the study will focus on the stability of surrounding rock and evolution of geological hazards during construction of tunnel and other underground engineering. Moreover, a more accurate numerical model considering multiple geological conditions based on the state-based peridynamics will be established to solve more complex engineering problems.

Acknowledgements Much of the work presented in this paper is supported by the Changjiang River Scientific Research Institute (CRSRI) Open Research Program (CKWV2018468/KY), the National Natural Science Foundation of China (51709159, 51911530214, 51679131), Shandong Provincial Key R&D Program of China (2019GSF111030), and the State Key Laboratory of Mining Disaster Prevention and Control Co-founded by Shandong Province and the Ministry of Science and Technology (MDPC201802).

7. Conclusions In this study, the numerical simulations of tunnel excavation process based on improved peridynamics are carried out. The characteristics of damage and deformation are analyzed to describe the stability of surrounding rock during tunnel excavation. According to the results of numerical simulations, following conclusions are obtained:

Appendix A. Supplementary material

(1) Peridynamics has been proved to be able to simulate the damage evolution and progressive failure process of rock materials. The simulation results shown that the compression failure behavior of rock materials can be accurately described by improving the peridynamic equation of motion through introducing a short-range repulsive force. (2) The EDZ is mainly caused by stress redistribution and tunnel excavation. The v-shaped notches developed on the tunnel periphery at the direction of perpendicular to the maximum principal stress. The phenomenon can be well simulated by introducing the material point dormancy method into peridynamics models. The predicted EDZ distribution law is consistent with the results of previous studies. And the displacement field of surrounding rock after tunnel excavation is also in good agreement with FEM simulation result, which proved the reliability of the proposed method. (3) The distribution characteristics of EDZ is closely related to the magnitude and direction of in-situ stress field. The damage degree of surrounding rock increases with increasing the buried depth of tunnel. And EDZ is mainly distributed in tunnel side walls when the vertical in-situ stress is the maximum principal stress. On the contrary, EDZ is mainly distributed in the vault and bottom of tunnel when the horizontal stress is the maximum principal stress. Therefore, the lateral pressure coefficient affects the main distribution position of EDZ. It is necessary to strengthen the support of tunnel vault and bottom when the lateral pressure coefficient is large than 1.0.

Supplementary data to this article can be found online at https:// doi.org/10.1016/j.tust.2020.103289. References Ai, D., Zhao, Y.C., Wang, Q.F., Li, C.W., 2019. Experimental and numerical investigation of crack propagation and dynamic properties of rock in SHPB indirect tension test. Int. J. Impact. Eng. 126, 135–146. Bi, J., Zhou, X.P., Qian, Q.H., 2016. The 3D numerical simulation for the propagation process of multiple pre-existing flaws in rock-like materials subjected to biaxial compressive loads. 2016, 1611–1627. Cai, M., Kaiser, P.K., 2005. Assessment of excavation damaged zone using a micromechanics model. Tunn. Undergr. Sp. Tech. 20, 301–310. Dayal, K., Bhattacharya, K., 2006. Kinetics of phase transformations in peridynamic formulation of continuum mechanics. J. Mech. Phys. Solids 54, 1811–1842. Diana, V., Casolo, S., 2019. A bond-based micropolar peridynamics model with shear deformability: Elasticity, failure properties and initial yield domains. Int. J. Solids. Struct. 160, 201–231. Diederichs, M.S., 2007. Mechanistic interpretation and practical application of damage and spalling prediction criteria for deep tunnelling. Can. Geotech. J. 44, 1082–1116. Diehl, P., Franzelin, F., Pfluger, D., Ganzenmuller, G.C., 2016. Bond-based peridynamics: a quantitative study of mode I crack opening. Int. J. Fracture 201, 157–170. Do, T.H., Wu, J.H., 2020... Simulation of the inclined jointed rock mass behaviors in a mountain tunnel excavation using DDA. Comput. Geotech. 117, 103249. Foster, J.T., Silling, S.A., Chen, W.W., 2010. Viscoplasticity using peridynamics. Int. J. Numer. Meth. Eng. 81, 1242–1258. Gao, C.L., Li, S.C., Wang, J., Li, L.P., Lin, P., 2018. The rick assessment of tunnels based on grey correlation and entropy weight method. Geotech. Geol. Eng. 36, 1621–1631. Gerstle, W., Sau, N., Silling, S., 2007. Peridynamic modeling of concrete structures. Nucl. Eng. Des. 237, 1250–1258.

18

Tunnelling and Underground Space Technology 97 (2020) 103289

C. Gao, et al.

deformation analysis for forward modelling. Int. J. Numer. Anal. Met. 13, 359–380. Silling, S.A., 2000. Reformulation of elasticity theory for discontinuities and long-range forces. J. Mech. Phys. Solids. 48, 175–209. Silling, S.A., Askari, E., 2005. A meshfree method based on the peridynamic model of solid mechanics. Comput. Struct. 83, 1526–1535. Silling, S.A., Bobaru, F., 2005. Peridynamic modeling of memebranes and fibers. Int. J. Nonlin. Mech. 40, 395–409. Silling, S.A., Epton, M., Wechner, O., Xu, J., Askari, E., 2007. Peridynamic states and constitutive modeling. J. Elasticity 88, 151–184. Song, S.G., Li, S.C., Li, L.P., Shi, S.S., Zhou, Z.Q., Liu, Z.H., Shang, C.S., Sun, H.Z., 2019. Model test study on vibration blasting of large cross-section tunnel with small clearance in horizontal stratified surrounding rock. Tunn. Undergr. Sp. Tech. 92, 103013. Sterpi, D., Cividini, A., 2004. A physical and numerical investigation on the stability of shallow tunnels in strain softening media. Rock. Mech. Rock. Eng. 37 (4), 277–298. Sun, C.Y., Huang, Z.X., 2016. Peridynamic simulation to impacting damage in composite laminate. Compos. Struct. 138, 335–341. Wang, J., Li, S.C., Li, L.P., Lin, P., Xu, Z.H., Gao, C.L., 2019a. Attribute recognition model for risk assessment of water inrush. Bull. Eng. Geol. Environ. 78, 1057–1071. Wang, Q., Wang, Y., Zan, Y.F., Lu, W., Bai, X.L., Guo, J., 2018. Peridynamics simulation of the fragmentation of ice cover by blast loads of an underwater explosion. J. Mar. Sci. Tech. 23 (1), 52–66. Wang, Y.T., Zhou, X.P., Kou, M.M., 2019b. Three-dimensional numerical study on the failure characteristics of intermittent fissures under compressive-shear loads. Acta. Geotech. 14, 1161–1193. Warren, T.L., Silling, S.A., Askari, A., Wechner, O., Epton, M.A., Xu, J.F., 2009. A nonordinary state-based peridynamic method to model solid material deformation and fracture. Int. J. Solids. Struct. 46, 1186–1195. Wu, Z.J., Zhang, P.L., Fang, L.F., Liu, Q.S., 2019. Numerical study of the effect of confining pressure on the rock breakage efficiency and fragment size distribution of a TBM cutter using a coupled FEM-DEM method. Tunn. Undergr. Sp. Tech. 88, 260–275. Yaghoobi, A., Chorzepa, M.G., 2017. Higher-order approximation to suppress the zeroenergy mode in non-ordinary state-based peridynamics. Comput. Struct. 188, 63–79. Yang, H.Q., Huang, D., Yang, X.M., Zhou, X.P., 2013. Analysis model for the excavation damage zone in surrounding rock mass of circular tunnel. Tun. Undergr. Sp. Tech. 35, 78–88. Zhang, C.S., Chu, W.J., Liu, N., Zhu, W.S., Hou, J., 2011a. Laboratory tests and numerical simulation of brittle marble and squeezing schist at Jinping II hydropower station, China. J. Rock. Mech. Geotech. Eng. 3 (1), 30–38. Zhang, Q.B., He, L., Zhu, W.S., 2016. Displacement measurement techniques and numerical verification in 3D geomechanical mdoel tests of an underground cavern group. Tunn. Undergr. Sp. Tech. 56, 54–64. Zhang, Y.N., Deng, J.R., Deng, H.W., Ke, B., 2019. Peridynamics simulation of rock fracturing under liquid carbon dioxide blasting. Int. J. Damage. Mech. 28 (7), 1038–1052. Zhang, Z.X., Hu, X.Y., Scott, K.D., 2011b. A discrete numerical approach for modeling face stability in slurry shield tunneling in soft soils. Comput. Geotech. 38, 94–104. Zhang, Z.X., Xu, Y., Kulatilake, P.H.S.W., Huang, X., 2012. Physical model test and numerical analysis on the behavior of stratified rock masses during underground excavation. Int. J. Rock. Mech. Min. 49, 134–147. Zhao, K., Bonini, M., Debernardi, D., Janutolo, M., Barla, G., Chen, G.X., 2015. Computational modelling of the mechanized excavation of deep tunnels in weak rock. Comput. Geotech. 66, 158–171. Zhou, X.P., Bi, J., Qian, Q.H., 2015. Numerical simulation of crack growth and coalescence in rock-like materials containing multiple pre-existing flaws. Rock. Mech. Rock. Eng. 48, 1097–1114. Zhou, X.P., Chen, J.W., 2019. Extended finite element simulation of step-path brittle failure in rock slopes with non-persistent en-echelon joints. Eng. Geol. 250, 65–88. Zhou, Z.Q., Ranjith, P.G., Yang, W.M., Shi, S.S., Wei, C.C., Li, Z.H., 2019. A new set of scaling relationships for DEM-CFD simulations of fluid-solid coupling problems in saturated and cohesiveless granular soils. Comput. Part. Mech. 6 (4), 657–669. Zhou, X.P., Shou, Y.D., 2017. Numerical simulation of failure of rock-like material subjected to compressive loads using improved peridynamic method. Int. J. Geomehc. 17 (3), 04016086. Zhou, X.P., Shou, Y.D., Berto, F., 2018. Analysis of the plastic zone near the crack tips under the uniaxial tension using ordinary state-based peridynamics. Fatigue Fract. Eng. M. 42 (5), 1159–1170. Zhou, X.P., Wang, Y.T., 2016. Numerical simulation of crack propagation and coalescence in pre-cracked rock-like Brazilian disks using the non-ordinary state-based peridynamics. Int. J. Rock. Mech. Min. 89, 235–249.

Hajiabdolmajid, V., Kaiser, P., 2003. Brittleness of rock and stability assessment in hard rock tunneling. Tun. Undergr. Sp. Tech. 18, 35–48. Ha, Y.D., Bobaru, F., 2010. Studies of dynamic propagation and crack branching with peridynamics. Int. J. Fracture 162, 229–244. Ha, Y.D., Lee, J., Hong, J.W., 2015. Fracturing patterns of rock-like materials in compression captured with peridynamics. Eng. Fract. Mech. 144, 176–193. He, J., Liu, Q.S., Wu, Z.J., 2017. Creep crack analysis of viscoelastic material by numerical manifold method. Eng. Anal. Bound. Elem. 80, 72–86. Hu, J., Xu, N.W., 2011. Numerical analysis of failure mechanism of tunnel under different confining pressure. Procedia Eng. 26, 107–112. Homand-Etienne, F., Hoxha, D., Shao, J.F., 1998. A continuum damage constitutive law for brittle rocks. Comput. Geotech. 22 (2), 135–151. Jiang, M.J., Yin, Z.Y., 2012. Analysis of stress redistribution in soil and earth pressure on tunnel lining using the discrete element method. Tunn. Undergr. Sp. Tech. 32, 251–259. Jia, P., Tang, C.A., 2008. Numerical study on failure mechanism of tunnel in jointed rock mass. Tunn. Undergr. Sp. Tech. 23, 500–507. Kasper, T., Meschke, G., 2004. A 3D finite element simulation model for TBM tunneling in soft ground. Int. J. Numer. Anal. Met. 28, 1441–1460. Kilic, B., Agwai, A., Madenci, E., 2009. Peridynamic theory for progressive damage prediction in center-cracked composite laminates. Compos. Struct. 90, 141–151. Kilic, B., Madenci, E., 2010. An adaptive dynamic relaxation method for quasi-static simulations using the peridynamic theory. Theor. Appl. Fract. Mec. 53, 194–201. Le, Q.V., Chan, W.K., Schwartz, J., 2014. A two-dimensional ordinary, state-based peridynamic model for linearly elastic solids. Int. J. Numer. Meth. Eng. 98, 651–657. Lee, J., Ha, Y.D., Hong, J.W., 2017. Crack coalescence morphology in rock-like material under compression. Int. J. Fract. 203, 211–236. Li, S.C., Gao, C.L., Zhou, Z.Q., Li, L.P., Wang, M.X., Yuan, Y.C., Wang, J., 2019. Analysis on the precursor information of water inrush in karst tunnels: a true triaxle model test study. Rock. Mech. Rock. Eng. 52, 373–384. Li, S.C., Wang, M.B., 2008. Elastic analysis of stress–displacement field for a lined circular tunnel at great depth due to ground loads and internal pressure. Tun. Undergr. Sp. Tech. 23, 609–617. Li, S.C., Zhou, Z.Q., Li, L.P., Xu, Z.H., Zhang, Q.Q., Shi, S.S., 2013. Risk assessment of water inrush in karst tunnels based on attribute synthetic evaluation system. Tunn. Undergr. Sp. Tech. 38, 50–58. Liu, M.H., Wang, Q., Lu, W., 2017. Peridynamic simulation of brittle-ice crushed by a vertical structure. Int. J. Nav. Arch. Ocean. 9, 209–218. Li, X.B., Cao, W.Z., Zhou, Z.L., Zou, Y., 2014. Influence of stress path on excavation unloading response. Tun. Undergr. Sp. Tech. 42, 237–246. Madenci, E., Oterkus, E., 2014. Peridynamic Theory and its Applications. Springer Science, New York. Ma, G.W., An, X.M., Zhang, H.H., Li, L.X., 2009. Modeling complex crack problems using the numerical manifold method. Int. J. Fract. 156, 21–35. Martino, J.B., Chandler, N.A., 2004. Excavation-induced damage studies at the underground research laboratory. Int. J. Rock. Mech. Min. 41, 1413–1426. Meguid, M.A., Saada, O., Nunes, M.A., Mattar, J., 2008. Physical modeling of tunnels in soft ground: a review. Tunn. Undergr. Sp. Tech. 23, 185–198. Meo, D.D., Zhu, N., Oterkus, E., 2016. Peridynamic modeling of granular fracture in polycrystalline materials. J. Eng. Mater. Tech. 138, 041008-1–16. Morris, J.P., Rubin, M.B., Block, G.I., Bonner, M.P., 2006. Simulations of fracture and fragmentation of geologic materials using combined FEM/DEM analysis. Int. J. Impact. Eng. 33, 463–473. Mortazavi, A., Molladavoodi, H., 2012. A numerical investigation of brittle rock damage model in deep underground openings. Eng. Fract. Mech. 90, 101–120. Oterkus, S., Madenci, E., 2017. Peridynamic modeling of fuel pellet cracking. Eng. Fract. Mech. 176, 23–27. Oterkus, S., Madenci, E., Oterkus, E., 2017. Fully coupled poroelastic peridynamic formulation for fluid-filled fractures. Eng. Geol. 225, 19–28. Parks, M.L., Lehoucq, R.B., Plimpton, S.J., Silling, S.A., 2008. Implementing peridynamics within a molecular dynamics code. Comput. Phys. Commun. 179, 777–783. Rivas, E., Parchei-Esfahani, M., Gracie, R., 2019. A two-dimensional extended finite element method model of discrete fracture networks. Int. J. Numer. Meth. Eng. 117 (13), 1263–1282. Schwartz, C.W., Einstein, H.H., 1980. Simplified analysis for groud-structure interaction in tunneling. In: Proceedings of the 21st US Symposium on Rock Mechanics. University of Missouri-Rolla, pp. 787–796. Senent, S., Jimenez, R., 2015. A tunnel face failure mechanism for layered ground, considering the possibility of partial collapse. Tunn. Undergr. Sp. Tech. 47, 182–192. Shi, G.H., 1991. Manifold method of material analysis, Trans. 9th Army Conf. on Appl. Math. and Comp. Rep No. 92-1. Shi, G.H., Goodman, R.E., 1989. Generalization of two-dimensional discontinuous

19