Modeling hard rock failure induced by structural planes around deep circular tunnels

Accepted Manuscript Modeling hard rock failure induced by structural planes around deep circular tunnels Fan Feng, Xibing Li, Jamal Rostami, Diyuan Li PII: DOI: Reference:

S0013-7944(18)30629-5 https://doi.org/10.1016/j.engfracmech.2018.10.010 EFM 6184

To appear in:

Engineering Fracture Mechanics

Received Date: Revised Date: Accepted Date:

21 June 2018 7 October 2018 8 October 2018

Please cite this article as: Feng, F., Li, X., Rostami, J., Li, D., Modeling hard rock failure induced by structural planes around deep circular tunnels, Engineering Fracture Mechanics (2018), doi: https://doi.org/10.1016/ j.engfracmech.2018.10.010

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Modeling hard rock failure induced by structural planes around deep circular tunnels Fan Feng1,2,3, Xibing Li1,3*, Jamal Rostami2, Diyuan Li1,3 1

School of Resources and Safety Engineering, Central South University, Changsha, Hunan, China

3

Earth Mechanics Institute, Department of Mining Engineering, Colorado School of Mines, Golden, CO,

USA 3

Hunan Key Laboratory of Resources Exploitation and Hazard Control for Deep Metal Mines, Changsha,

Hunan, China Corresponding author: Xibing Li (E-mail: [email protected]) Abstract: Construction of tunnels is often associated with fault or structural features that could affect tunnel stability during the construction phase and service life. Therefore, failure characteristics of hard rock around circular tunnels induced by pre-existing structural features under high in-situ stresses have been the subject of various studies. In the current study, a combined finite element approach, namely ELFEN, has been used for better reflection of entire failure process (including crack initiation, propagation and coalesence) and intrinsic properties of hard rock mass, thus rock heterogeneity, around circular tunnels during excavation unloading process. Parametric analysis which consider the dip angle, location (exposure or not), frictional coefficient of structural planes and lateral pressure coefficient was conducted in detail to reveal the mechanical responses of circular tunnel induced by structural plane under unloading condition. Numerical results indicate that the failure intensity of rock tunnel is a function of both dip angles and frictional coefficients of structural planes. The most critical dip angles of structural planes leading to failure in tunnel rocks largely depend on the frictional coefficient. Also, the results indicate that the released strain energy for the case of exposed structural plane is higher than those not intersecting the tunnel, leading to more

violent rock failure for the former. With the increase of the lateral pressure coefficient, the failure intensity and damage extent around the tunnel is aggravated, especially for the roof and floor of the tunnel. Rock

failure can be categorized as slabbing failure near excavation boundary and shear slip failure, controlled by structural plane. Progressive slabbing failure induced by excavation unloading may activate internal structural planes and the extensive release of energy caused by shear and slip failure may in turn further induce the slabbing failure. Rockburst is more prone to be triggered under such condition. Keywords: Circular tunnel; Structural plane; Rockburst; Crack propagation; Numerical simulation; Energy evolution

1. Introduction

The increasing demand for resources has led to construction of a large number of mining and tunneling projects at depth due to the exhausting of reserves at shallow depths (Li et al., 2018b, Li et al., 2017). Similarly, the need for passing under high mountains and need for access to certain points inside the ground, such as China, Iran and Norway etc., has been the reason for construction of deep tunnels around the globe (Dammyr et al., 2017; Sun et al., 2018; Farrokh and Rostami, 2008; Farrokh and Rostami, 2009). In-situ stress measurements had proved extremely distinct failure phenomenon in deep tunnels or mines such as spalling, zonal disintegration and rockburst, which are threat to the safety of construction workers and can cause substantial damages to the structure and equipment (Li et al. 2018b; Martin and Maybee, 2000; Mortazav and Molladavoodi 2012). This is primarily because excavation of the tunnel leads to stress redistribution of rock mass around the tunnel (Lu et al. 2012; James and Jr 2003). The radial principal stress turn to zero at tunnel wall (lack of lateral support for the wall elements) and high concentration of mainly compressive tangential principal stresses in the periphery (Lu et al. 2012; Engineering Fracture Mechanics). The formation of free surface caused by excavation unloading process also provides the favorable stress conditions for imitation of fracture in rock that is prelude to ground failure. The occurrence of deep rock mass collapse, rockburst and mine earthquake is closely related to the geological structure, rock and rock mass properties, characteristics of induced stress field and excavation path. Fissured rock masses are complex engineering medium that is widely encountered in underground engineering, such as water tunnels and hydro power structures, mining and transportation (Wang et al. 2018). Structural discontinuities and weakened planes play a very important role in the mechanical properties of underground deep rock mass (Feng et al. 2013; Fu et al. 2017). Various studies show that a deep underground chamber is more prone to failure when it approaches a geological discontinuity such as fault, dyke, and contact (Hedley et al., 1992; Snelling et al., 2013; Manouchehrian and Cai 2017; Wasantha et al. 2014; Zhao et al. 2017). Some researchers have focused on the mechanical properties of rock materials with prefabricated fissures to examine the ground behavior in such conditions. Bobet and Einstein (1998) investigated the fracture coalescence in rock-type materials under uniaxial and biaxial compression, with two pre-existing fractures or flaws that were arranged in different geometries. Wong et al (2002) analyzed the crack propagation laws of the prefabricated fractured rock specimen under

uniaxial compression by using RFPA software (Realistic Failure Process Analysis). The analysis of parameters, such as angle and length of the flaws, specimen width and the arrangement of flaw locations, was conducted to examine its influence on the growth and coalescence behavior. Zhao et al. (2015) adopted digital image correlation technology to investigate the deformation and failure characteristics of rock-like materials with a pre-existing single fissure under uniaxial compression, and the crack initiation, propagation and the rock damage were analyzed quantitatively based on the global strain field on the surface of specimen at meso-level. Zhou and Wang (2005) studied the coalescence mechanism of splitting failure of crack-weakened rock masses under compressive loads, and a simplified mechanism of crack propagation, in which the crack grows along the direction of maximum principal compressive stress, was proposed. Yang et al. (2016) investigated the influence of petroleum on the failure pattern of saturated pre-cracked and intact sandstone. They found that the failure for the pre-cracked specimen is strongly influenced by pre-existing flaws, and the fracturing degree is aggravated due to the immersion in petroleum, as stress concentration at a high level will form at the tips of a pre-existing flaw. Other researchers have concentrated on failure characteristics around tunnels or rock samples containing central holes by means of numerical simulation, laboratory tests, theoretical analysis and field monitoring. Li et al. (2017) investigated the effects of pre-fabricated holes with different geometries, including size, shape and inclination angle, on the strength and fracture behavior of prismatic marble specimens. They found that the propagation of tensile cracks surrounding the hole in marble specimens was mainly affected by the nucleation and propagation of strain localization zones. Zhou et al. (2017) studied the fracture coalescence behaviour of rectangular cuboid marble specimens containing a single rectangular cavity or two rectangular cavities. They found that both cavity number and layout have a significant influence on mechanical behaviour of marble specimens. Yang et al. (2018) studied the dynamic stress adjustment and rock damage

during blasting excavation in a deep-buried circular tunnel, and a theoretical model for a two-dimensional (2D) circular excavation was developed to investigate the stress evolution and resultant rock damage arising from millisecond-delay blasting. Li et al. (2014) used PFC2D to investigate the influence of stress path on excavation unloading response, and a mathematical model was applied to characterize the unloading mechanisms of brittle rock under different stress path. Zhong et al. (2018) performed a biaxial compression test on rhyolite to determine the fracture mechanism of the hard and brittle surrounding rock containing natural cracks. Results indicated that the pre-existing natural cracks in the hard and brittle rock played an important role in consuming strain energy by

initiating and propagating fractures during biaxial compression. Zhu et al. (2005) studied the progressive fracturing processes around underground excavations under biaxial compression by using RFPA. In their study, lateral coefficients of horizontal stress and tunnel shapes were considered. It was showed that the mechanical parameters of the red sandstone containing two unparallel fissures were weaker than those of intact specimens, but were strongly dependent on fissure angle, and the macroscopic deformation characteristics of red sandstone under unaxial compression were systematically evaluated. It should be noted that while much effort has been dedicated to understand failure properties of rock materials containing pre-existing cracks or central holes, little attention has been paid to the effects of the integrated response of structural plane and excavation unloading effect. Moreover, the frictional effect of pre-existing cracks was seldom considered in previous studies. This could be due to the difficulties in preparation of rock samples in laboratory tests or acquiring the exact structural plane parameters in field investigation. Zhou et al. (2015) analyzed the rockburst mechanisms induced by structural planes in deep tunnels. They stated that weak planes could induce rockburst in tunnels with three possible mechanisms including fault-slip, shear rupture, and buckling. However, their research did not consider the effect of lateral stresses, excavation unloading process, dip angle, and frictional coefficient of structural planes. As such, their conclusions were mainly based on a series of qualitative analysis. In practice, the above-mentioned parameters are of vital importance in controlling the failure characteristics and mechanical behavior of rock mass around deep-buried tunnels (Ashby and Hallam 1986; Zhu et al. 2005; Cao et al. 2016; Wong et al. 2002). Numerical approaches are the main tools used for simulating the fracturing processes of a rock and rock mass, because modeling can greatly reduce the cost of studying rock behavior in various stress states. Manouchehrian and Cai (2017) studied the role of discontinuities around tunnels in rockburst by using the Abaqus. They found when a weak plane was presented in the model, it resulted in more released kinetic energy and higher element velocity, indicating that rock failure could become more violent. Zhang et al. (2013) adopted FLAC3D to investigate the rock mass damage development following two extremely intense rockbursts induced by fault around deep tunnels at Jinping II hydropower station, in southwestern China. However, due to the selected specific numerical code in their study, the entire failure process from continuum to dis-continuum could not be observed. Therefore, it is important to select a suitable numerical tool in order to capture the entire rock failure process including crack initiation, propagation and coalescence

stage. The FDEM approach (also called FEM/DEM combined method), which integrates the FEM (Finite Element Method) and DEM (Discrete Element Method) into a singular tool, has provided a medium for both continuous and discontinuous problems (Vyazmensky et al. 2010; Elmo et al. 2013; Mitelman and Elmo 2014; Hamdi et al. 2014; Li et al. 2018a; Feng et al. 2017). For example, Antolini et al. (2016) adopted FDEM to describe the simulation results of the triggering mechanism and the evolution scenarios of the Torgiovannetto di Assisi rockslide (in central Italy). Cai (2008) used ELFEN software and showed that forming of tunnel surface-parallel fractures and microcracks was attributed to material heterogeneity and the existence of relatively high intermediate principal stress, as well as near zero minimum principal stress. In his numerical models, the Mohr–Coulomb constitutive model with variable crack orientation was adopted to achieve tensile and shear failure in hard rocks. The authors (Feng et al. 2017; Li et al. 2018a) also adopted ELFEN software to investigate the rock failure characteristics by splitting test of circular ring and uniaxial compressive tests, the entire failure process was fully demonstrated by this method, and numerical simulation results showed good agreement with laboratory testing. Note that though FDEM approach has been employed in analyzing rock failure process in the aspects of laboratory scale or field measurements, little attention was paid to the failure mechanisms of rock around deep tunnels affected by near-by faults or weak planes. In the present study, a FDEM approach (ELFEN, Rockfield Software Ltd.) is adopted to simulate the hard rock fracture induced by structural planes around deep-buried circular tunnels, with particular attention being paid to the effect of dip angle, location (exposure to the wall or being sealed and inside the ground), frictional coefficient of structural planes, and lateral pressure coefficient. Characteristics of a typical hard rock (Marble) extracted from Hongling Lead-Zinc mine, located in Inner Mongolia, China was selected for representing material properties in the numerical simulation. Rock heterogeneity is also considered by specifying a certain range of modeling parameters in the current study, and a linear unloading process is adopted to realistically reflect gradual relaxation during excavation of a tunnel. Parametric analysis was carried out to investigate the failure modes, evolution of released strain energy and effective plastic strain around an underground tunnel. The entire failure processes including crack initiation, propagation and coalescence, induced by structural planes, are reproduced to further investigate the rock failure mechanism under different conditions. This paper is intended to provide an insight into the failure mechanism of hard rock, induced by structural planes in the periphery of tunnels to be excavated. 2. Modelling circular tunnel fracturing induced by structural plane

2.1 Brief introduction of fundamental principles of FDEM based on ELFEN The combined finite/discrete approach is a numerical method that combines continuum mechanics principles with dis-continuum algorithms to simulate multiple interacting deformable solids (Elmo et al. 2013; Feng et al. 2017). This approach was proposed by Munjiza (2004), who takes advantage of both techniques and to combine the FEM and DEM into a hybrid finite/discrete modeling methodology. In FEM/DEM, each discrete element is discretized into finite elements, so transition from continuum to dis-continuum is done through a fracture and fragmentation process. If the failure criterion within the intact rock (initially represented as a FEM domain) is met, then a crack is initiated. ELFEN, a software developed by Rockfield Software Ltd. (Rockfield, 2013), contains a variety of constitutive models and is capable of both explicit and implicit analyses of 2D and 3D types. In the present study, ELFEN was adopted to simulate the failure mechanism of hard rock failure induced by structural plane around circular tunnels. An obvious distinction from other numerical codes is its ability to allow the new cracks to cut across the original mesh (intra-element fracturing) rather than along or around the element boundary (inter-element fracturing) merely (Li et al. 2017). Whether the fractures split the elements is determined by the value for“Smallest Element” set under the discrete contact parameters. If this value is the same as or larger than the original mesh size, the fractures will propagate along the element boundaries only, as the minimum size has already been reached. If it is smaller than the original mesh, the fractures may split elements, if conditions dictate, as can be seen in Fig. 1. For a 2D plane strain problem, the Mohr–Coulomb constitutive model with rotating crack was selected because it can simulate rock fracturing due to both tension and compression currently. A detailed description of fundamental principles of FDEM based on ELFEN and the selected constitutive model could refer to Feng et al. (2017), Li et al. (2018a) and Cai (2008). 2.2 Numerical modeling The geometry of the models used in this study is shown in Fig. 2. For 2D simulations, plane strain analysis in ELFEN assumes unit thickness (1 mm in the present study). The rock mass is represented by a single face, with a circular hole representing the tunnel to be excavated. The dimension of the area that is modeled is 20m × 20m, and the tunnel diameter is 4m. A closed structural plane, with a length of 2m, was included in the model located on the top of the circular hole (The structural plane in the sidewall was not considered in order to avoid repetition due to the circular shape of the tunnel). Four parameters were considered in this study: (1) dip angle of structural plane; (2) location of structural plane (exposed or unexposed); (3) frictional coefficient of structural plane; (4) coefficient of lateral stresses. Four dip angles at

30º, 45º, 60º and 90º are modeled. For the position of the structural plane relative to the tunnel, two conditions were simulated, one with the discontinuity exposed with a distance of 0m to the tunnel boundaries and unexposed plane with a distance of 6m to the tunnel boundaries. The frictional coefficients of the structural plane in the models were set to be 0.1, 0.5 and 1. Various in-situ stress states were also modeled with four lateral pressure coefficients k at 0.5, 1, 1.5, and 2, respectively. Four monitoring points in roof, floor and sidewalls of the circular tunnel were selected to monitor the effective plastic strain and tangential stress. Properties of a typical hard rock (Marble) from Hongling Lead-zinc mine was used for numerical simulation. For a realistic representation of the properties of the underground rock mass, the material heterogeneity was considered in the present study. Note that the heterogeneity in the current study is irrelevant to the direction inside the material. It is mainly used for characterizing the discrepancy of physical and mechanical properties in each location or position (Rockfield 2013). Material heterogeneity was applied on an element by element basis, allowing each element to have an independent value assigned for Young’s Modulus, Possion’s ratio, density, tensile strength, and fracture energy. In ELFEN, a random number generator was used to create the values for each property for each element, within specified limits either side of a mean value. The limits were given as a percentage of the mean value. The random number generator also has the ability to “seed” the first number, so that the “random” series of numbers is repeated. The distribution of the Young’s modulus, Possion’s ratio, density, tensile strength and fracture energy are presented in Fig. 3. The scale in the legend indicates the ratio of the each parameter to the average value (Cai 2008). For example, the number 1.135 for the Possion’s ratio at a point inside the model means that the Possion’s ratio at that point is 13.5% larger than the average value of 0.23. The physical properties of the hard rock and the discrete contact parameters used in the modeling are shown in Table 1. The total numbers of unstructured triangle elements established by the model for solving the problem were 4568. Because stress concentration due to unloading process and the slip of the structural plane was expected, very dense mesh was used in this area close to the tunnel boundaries. The mesh dimension for the rest of the regions was set to be larger in order to reduce the computation time. An important feature of ELFEN is that an explicit fracture could be introduced when all the fracture energy was consumed (Feng et al. 2017; Li et al. 2018a). This limit was defined as the amount of energy needed to create continuous crack per unit area. In fact, the fracture energy used for study of hard and brittle materials often ranges from 0.01 to 0.3 N/mm (Klerck 2000; Cai 2013; Hamdi et al. 2014) as noted in

previous research works in this field. Therefore, the value of 0.02 N/mm was selected for describing properties of hard and brittle Marble. Penalty parameters were used to evaluate the normal and tangential contact forces. The value of the normal penalty was usually in the range of 0.5E to 2.0E, where E is the Young’s Modulus, and the tangential penalty was about 0.1 of the normal force penalty (Rockfield 2013). The contact field defines the thickness of contact layer which corresponds to the maximum permissible penetration. The value issued in modeling is normally 10-20% of the minimum side length of any element in the mesh. In the present study, the smallest element was set to be 50% of the original mesh, which allows newly-generated cracks (distinguish from structural plane) to propagate either along or split elements. In 2D plane strain problem, Face loading is selected to apply in-situ stresses to the model boundaries, where vertical stress is applied on the top domain while horizontal stresses are applied on both sides of the model. Linear loading over 0.01s is employed to avoid instantaneous loading which results in violent and unwanted fracture propagation in the model. After 0.01s, the face loads in each direction are remained constant until the end of the numerical calculation. The excavation relaxation was specified on the circular tunnel object. A linear relaxation curve was assigned to the circular tunnel representing the excavated material, in this case a linear ramp over 0.005s starting from 0.01s was used. When the tunnel is deactivated, the relaxation was automatically applied to release the stresses according to the specified curve. This avoids the sudden change in stress state which would occur with instantaneous removal of material. The bottom boundary is fixed in the vertical direction. The loading and unloading paths are shown in Fig. 4. Note that

the numerical simulation is terminated till t=0.02s. The in situ stress field of the Hongling Lead-Zinc mine was employed in this study, and the fitting equations of the vertical principal stress, maximum horizontal principal stress and minimum horizontal principal stress are presented, respectively, as follows:  v  0.0279h,    h max  0.0256h  10.202,   h min  0.0233h  1.4743. 

(1)

In the present study, a circular tunnel at a burial depth of 1000m was selected to reflect the deep highly-stressed structure. The horizontal and vertical stress components are σv=27.9MPa, σhmax=35.8MPa, σhmin=24.8MPa according to Eqn. (1). Note that σhmin is not used because of the 2D simulation problem. Parametric analyses were conducted based on specific in-situ stresses except the cases for various lateral pressure coefficients. In the present study, a numerical viscous damping (i.e., proportional to velocity) is introduced in the FDEM simulation. Velocity proportional (viscous) damping can be applied to any type of entity and to some

or all degrees of freedom, as required. It is specified as a percentage of critical damping. ELFEN calculates the critical damping factor automatically, and then multiplies by the factor specified by the user. The factor ranges from 0 (i.e. no damping) to 1 (i.e. critical damping). The damping applied to all the nodes created for entities to which it is assigned in the pre-processor. For example, if all the surfaces in a 2D problem are selected, then every node in the problem will be subjected to point damping (Rockfield, 2013). Actually, a damping value of 0.1 or less should be used for dynamic relaxation simulation of linear problems, so we choose a damping value 0.08 in the current study. During dynamic unloading process, when unloading stress wave arrives at the model boundary, a certain of tensile stress wave might be reflected back from the constrained boundary given the fact that the prescribed boundaries is inappropriate or finite. To tremendously reduce or even eliminate the influence of reflected waves on simulation results, the “non-reflecting boundary” set in ELFEN numerical simulation was employed in the current study. The non-reflecting boundary enables high frequency waves to pass though the boundary rather than being reflected, as in the case of either a prescribed or free boundary (Rockfield, 2013), which guarantees the accuracy and validity of numerical results. 2.3 Model validation A series of models without the structural plane were set up and analyzed for comparison with the model containing structural plane to validate ELFEN software package. Four lateral pressure coefficients of k=0.5, 1, 1.5, and 2 were specified with vertical stress being 30MPa. Fig. 5 presents the final maximum principal stress (also in this case tangential stress) and effective plastic strain contour lines under the condition of each lateral pressure coefficient. After excavation and unloading of the walls, the stresses around the tunnel are disturbed, with the radial principal stress being released and tangential principal stress increasing to reach an equilibrium. By observation of the maximum principal stress contour lines, the stress concentration seems to be mainly scattered in two sides of circular tunnel when k=0.5. When k is raised to 1, the stress distribution around the tunnel tends to be uniform compared to the former case, as anticipated. Further increase in k to 1.5 and 2 will lead obvious stress concentration appearing in tunnel roof and floor. As stated by Li and Weng (2016), when the k is less than 1.0, two sidewalls of the opening are subjected to high compressive stress, which results in high strain energy intensity near the sidewalls. But when the k is larger than 1.0, the strain energy mainly intensifies at the roof and floor. Numerical results show good agreement with the analytical solutions and previous findings (Li and Weng 2016; Cao 2016; Li et al. 2018c). By observation of the newly generated cracks and the effective plastic strain around the circular tunnel,

the evolution law with the increasing k in the present study seems to show discrepancy compared with the case of homogeneity. It can be seen that the cracks distribution is strongly related to effective plastic strain. Only when obvious effective plastic strain occurs, the cracks begin to initiate. However, the effective plastic strain is not symmetrically distributed around the circular tunnel. This is due to the rock heterogeneity that leads to strength anisotropy of each grid point in the model. Therefore, due to the diversity of the strength in each point, the plastic strain (thus cracks) may appear even when the maximum principal stress is less than the average uniaxial compressive strength of the material. This is why the cracks could still be observed around the circular tunnel when k=0.5 and 1. Fig. 6 presents the final tangential stress distribution along the radial direction at four monitoring points. It can be seen that Fig. 6a and b show different trends to Fig. 6c and d, especially near excavation boundary. The reduction of tangential stress near excavation boundaries in both point 5 and 6 illustrate that the rock mass around has entered into plastic behavior. Under such stress state, wall will fracture at point 5 and 6 with visible cracks, while point 7 and 8 are still in stable state (always present a monotonous decreasing trend with the increase of distance). This is due to the material heterogeneity that leads to obvious stress concentration in certain places. Note that when k=1.5 and 2, most of the cracks are observed in the tunnel roof and floor, illustrating the higher stress concentration in these areas caused by excavation unloading effect. Fig. 7 presents the released strain energy-time evolution curves for various lateral pressure coefficients. It can be seen that for k=0.5 and 1, the released strain energy is about 11KJ, which remains a relatively low level. For the lower k, the tangential stress concentration (accumulated strain energy) is lower, so the released strain energy used for fracturing rock at the tunnel should not be very high. Under such condition, fewer cracks are observed around the circular tunnel, and the failure intensity is weaker. However, with the increase in the k, more strain energies are stored in the rock mass. Once unloaded being, the stored strain energy is bound to be released largely due to the stress re-distribution and concentration, which might lead to violent and severe destruction around the tunnel, even the occurrence of rockburst. 3. Results and Discussion 3.1 Influence of dip angle of structural plane In this section, the failure characteristics of an underground tunnel at a depth of 1000m are investigated, with a constant frictional coefficient of 0.1, and the structural plane as assumed to be exposed to the tunnel wall. The dip angle of structural plane varied from 30º to 90º. Fig. 8 and 9 depict the maximum principal stress and effective plastic strain contour lines in these models. It appears that for structural plane with

orientation of 45º and 60º from the horizontal plane the cracks propagate easier. The conditions of dip angle of 30º and 90º are the most difficult for fracture propagation. In the case of θ=30º, only few cracks were generated around the tunnel and at the structural plane tip. These cracks could be caused by the excavation unloading and local stress concentration around the tip. Note that the typical slabbing failure is observed around tunnel boundary, which is similar to the previous findings by Li et al. (2011) and Li et al. (2018b). The fewer generated cracks may illustrate lower influence of structural plane on the tunnel fracturing. With the increase in dip angle (45º), the failure extent gets more violent and the failure range expands, with more cracks generated at the tunnel roof around the structural plane. Under such circumstances, the closed structural plane is activated by the progressive slabbing failure caused by excavation loading. With further concentration of tangential stress, the normal force on the structural plane is gradually reduced, and slide and dislocation is triggered by increasing shear stress. Meanwhile, the activated weak plane may in turn cause more intense slabbing failure, as can be seen in Fig. 8b and 9b. The condition is further aggravated when the dip angle is increased to 60º, with more blocks and thin rock plates ejecting towards the free surface, which is similar to the rockburst phenomenon encountered in deep underground openings (He et al. 2010; Feng et al. 2016). However, when dip angle is increased to 90º, the rock failure around tunnel is reduced. Only local failure is observed in the roof and right side of the tunnel, which should be controlled by the localized rock properties. According to Zhou et al. (2015), excavation unloading will lead slabbing failure near tunnel boundaries. For the condition of dip angle of 90º, shear and slip failure will not be activated by the structural plane in the tunnel wall because the tangential stress is normal to the structural plane, so extensive sliding and dislocation could not be observed in practice. Moreover, the existence of structural plane along the radial direction may reduce the integrity of the thin rock plates. Therefore, the rockburst phenomenon is most unlikely to be observed in this case. Fig. 10 presents the released strain energy-time evolution curves for different dip angles. It can be seen from Fig. 10 that during the elastic loading stage (0-0.01s), no released strain energy is observed for all the cases since the total energy in rock mass is gradually accumulated instead of being released. After 0.01s, the tunnel is excavated with unloading internal forces at tunnel walls. During the unloading stage, the released strain energy begins to increase and remains constant once the surface force is unloaded to 0 completely. A closer look at Fig. 10 shows that with the increase in dip angle, the released strain energy presents various trends of initial increase, followed by a subsequent decrease. The maximum released strain energy (about 500KJ) is registered when the dip angle is equal to 60º, then gradually decreases at 45º and 30º, finally the

90º. Simulation results discussed in this section suggest that the dip angle of structural plane do play an important role in affecting the failure behavior of circular tunnel. 3.2 Influence of frictional coefficient of structural plane In this section, three different frictional coefficients for structural plane of 0.1, 0.5 and 1 were considered in the analysis, respectively. Fig. 11 to 14 show the numerical results of entire failure process of tunnels with various frictional coefficients at dip angles of 30º, 45º, 60º and 90º for the offsetting with exposed structural plane. Fig. 15 presents the released strain energy-time evolution curves with various frictional coefficients at different dip angles. Crack initiation, propagation and coalescence around the circular tunnels are depicted in detail. By observation of these plots, it is found that the crack is initiated before the internal surface load is unloaded to 0, which means the rock mass has been damaged during the unloading process. The crack initiation location is always found in the roof tunnel approaching the excavation boundary, irrespective of the frictional coefficients and the dip angles. This was expected since the maximum tangential stress is concentrated in the tunnel roof and floor based on the applied in-situ stresses. The main reason for lack of crack initiation at the tunnel floor should be the location of structural plane. With further unloading, the cracks begin to propagate more consistently. However, a deep insight into Fig. 11 to 15 shows that the failure intensity and damage range may show distinct variation trend in practice. In the case of θ=30º dip angle, the failure intensity present a monotonous increasing trend with the increase of frictional coefficient. When μ (for convenience, we use μ to illustrate the frictional coefficient) is 0.1, only slabbing failure with visible cracks are observed around the tunnel. With the increase in μ, the final failure is gradually aggravated and the released strain energy is increased, typically for the case of μ=1. It can be seen from Fig. 15a that the corresponding maximum released strain energy value reached 130KJ, which is about 2-3 times of the former two conditions with lower friction coefficient. In fact, the higher the released strain energy, the easier the cracks propagate and coalesce, and the more serious the failure intensity and severity. Additionally, Fig. 11c clearly present the slabbing failure coupled with shear slip failure, with more rock thin plates and anomalous rock bocks collapsed from the roof tunnel around the structural plane. The influence of the structural plane is maximized with μ of 1 when dip angle is 30º. However, for the case where θ=45º, an obvious distinction can be made for the variation of failure intensity compared to that of θ=30º. Numerical results indicate that the most favorable frictional coefficient for crack propagation and nucleation around the structural plane is 0.5. The failure intensity around the circular tunnel presents a trend of initial increase, followed by a subsequent drop. Fig. 15b shows the

corresponding released strain energy values for the case of μ=0.1, 0.5 and 1 are 132.3KJ, 324.5KJ and 104.9KJ. By observation of Fig. 12, it is found that more cracks are generated when μ=0.5, and these cracks are mainly distributed in the roof and floor of the tunnel. When looking at the failure characteristics influenced by frictional coefficients at dip angle of 60º, the variation of failure intensity and extent of a monotonous decreasing trend with the increasing frictional coefficients can be deducted. Fig. 15c shows that the corresponding released strain energy values for the case of μ=0.1, 0.5 and 1 are 510.4KJ, 85.4KJ and 68.3KJ respectively. This means the structural plane with μ=0.1 is the easiest to initiate crack propagation and coalescence (intense failure), while the structural plane with μ=1 is the most difficult for crack propagation (slight failure). Fig. 13 shows that only slabbing failure accompanied by rock blocks and rock plates are observed in the roof when μ=0.5 and 1, while the failure intensity and range is most serious under the condition of μ=0.1. Fig. 14 presents the failure process of circular tunnel with various μ at dip angle of 90º. It can be seen that the failure intensity and the damage extent present no obvious correlation with frictional coefficient of structural plane. Moreover, the tunnel stability is well maintained, with less cracks generated. Fig. 15d shows that the released strain energy values for the case of μ=0.1, 0.5 and 1 are 35.8KJ, 41.1KJ and 40.2KJ respectively. This means the frictional coefficient has little influence on the tunnel failure under such condition. This could be due to lack of shear stresses along the structural plane at this configuration, combined with the addition of normal stresses at the joint surface. Based on the above-mentioned analysis, it may be inferred that the failure characteristics of hard rock around circular tunnels is the integrated response of dip angle and frictional coefficient. Fig. 16 presents the maximum released strain energy distribution influenced by frictional coefficient at different dip angles in the current study. It can be observed that the favorable dip angle leading to rock failure (crack propagation and nucleation) around the tunnel heavily depend on the frictional coefficient. Actually, these two factors are not isolated in affecting the dynamic responses of the underground tunnel. According the description of Wong et al. (2002), the crack initiated from the flaw inclining at 25º is the most difficult to propagate, and the flaws with orientation near 60º to the axial loading direction is easier to occur. In their study, the dip angle denotes the angle between the pre-existing flaws and the loading direction (thus maximum principal stress), and a frictional coefficient of about 0.2 was selected for numerical simulation. It can be inferred that our findings are similar to their results, where for a relatively lower frictional coefficient (0.1-0.2), the favorable dip angle of structural plane leading to rock failure also correspond to a higher value (about 60º). According to

the analysis of tunnel stability using analytical solutions by Ashby and Hallam (1986), the most favorable flaw angle for crack nucleation is equal to  =1/ 2 tan 1 (1/  ) . Moreover, Wong and Chau (1998) also investigated crack coalescence laws in a rock-like material containing two cracks. In their laboratory tests, the frictional coefficients are in the range of 0.6–0.9, and their experiments showed that the most favorable angles θ of the structural plane for crack nucleation were in the range of 25–30º from axial load. Similarly, their conclusions were also analogous to our results. Note that previous research works mainly concentrated on the influence of either dip angle or frictional coefficient on the rock failure, which did not consider the integrated responses and mutual interactions between the two factors. Nevertheless, previous research also validates our findings and conclusions. 3.3 Influence of location of structural plane In underground tunnel engineering, some structural planes intersect tunnel boundaries (exposed) while others are concealed in the rock (unexposed) but close to the proximity of tunnel walls. To investigate the influence of structural plane location on the rock failure characteristics, the unexposed structural plane around the circular tunnel was simulated to make comparison with the results of the modeling of exposed structural plans that was discussed in section 3.1. The frictional coefficient of the structural plane was set to be 0.1. Fig. 17 and 18 present the maximum principal stress and effective plastic strain contour lines of model containing unexposed structural plane at different dip angles. In this condition, the failure intensity and extent of damage in the rock around the circular tunnel are basically the same as that of exposed condition, which present a trend of initial increase, then subsequent decrease with the increase of dip angle. The newly generated cracks (effective plastic strain) are mainly distributed in the crown area. This seems to be due to the consistent fracturing at the tunnel roof that causes the stress re-distribution around the tunnel. Fig. 19 presents the released strain energy-time evolution curves for different dip angles of unexposed structural plane and shows the sequence of the released strain energy values to be Ed60º> Ed45º> Ed30º> Ed90º.

Combined with released strain energy values, one can conclude that only local slabbing failure can be found for the case of θ=30º and 90º. The obvious shear slip failure induced by structural plane is observed for the case of θ=45º and 60º. For the latter one (60º), the rock blocks ejection and extensive roof collapse phenomena could be seen in Fig. 17c and 18c, illustrating the severity of failure due to the shear slip effect. Fig. 20 shows the maximum released strain energy value of the model containing exposed and exposed structural plane at different dip angles. It can be inferred that the released strain energy under unexposed

condition is always lower than that of exposed condition, which means the failure intensity and damage extent is also weaker than the exposed condition. The distinction becomes rather obvious for θ=45º and 60º. This could be explained by presence of intact rock between the excavation boundaries and the structural plane. It is the intact rock that dampens the activation effect of the structural plane. Therefore, the shear and slip may be difficult to start. In addition, under such circumstance, the location of structural plane is away from the excavation boundary, meaning the degree of tangential stress concentration around the unexposed structural plane is lower than the exposed one, and the shear stress acting on the structural plane is also reduced. This also leads to conclusion that rockburst is less likely to occur if there are the unexposed structural planes around the hard rock tunnel under the same stress conditions. According to the description of Zhou et al. (2015), the intensity of shear rupture burst evoked by unexposed inclined structural plane may be higher than that of slip burst induced by exposed plane under the same condition. In their statement, the failure process is accompanied with both shear failure from intact rock mass between the structural plane and the excavation boundary and the shear and slip failure from structural plane, which leads to a higher energy release for the unexposed condition. Actually, although the failure process involves both shearing of intact rock and shearing-off the asperities or filling materials of structural plane, the dynamic response is greatly weakened because of reduced tangential stress concentration around the structural plane and reduced disturbance from the slabbing failure blocked by intact rock. For the exposed structural planes, the mutual interaction between the slabbing failure induced by excavation unloading and the shear and slip failure activated by structural plane is further accelerated and enhanced, so the failure intensity and damage extent is more severe and violent. 3.4 Influence of lateral pressure coefficient on rock failure Underground openings are seldom located under hydrostatic pressure, so it is necessary to investigate the influence of the ratio of horizontal and vertical in situ stresses on rock failure, at presence of structural plane. In this section, the conditions under various k with different dip angles of unexposed structural plane are modeled. The in-situ stress states are identical to the conditions discussed in section 2.3. Fig. 21 to 24 plot the failure process of circular tunnel with various k values at different dip angles of structural plane. The numerical results show that the failure intensity and damage extent are strongly dependent on k. When k is 0.5, fewer cracks are observed around the excavation boundaries. The local slabbing failure can only be found at right side of and top of the tunnel (the left side presents no failure might be attributed to the rock heterogeneity). When k is raised to 1, the failure intensity shows no obvious distinction compared to the case

of k=0.5 since the degree of stress concentration around the tunnel is still kept at a lower level. It is also found that the visible cracks at the right side of the tunnel are reduced while more cracks are generated at the roof tunnel close to structural plane tip. When k is 1.5, the failure intensity gets worse, and most of the cracks are distributed in tunnel roof and floor, especially in the tunnel crown. With the increase in k above 1, the tangential stress will concentrate around the roof and floor. The slabbing failure near excavation boundaries is further promoted and extended inside the walls. Under such circumstance, the structural plane may be activated by slabbing failure. The shear stress acting on the structural plane will lead to the extensive shear and slip failure, and it may in turn promote the slabbing failure, as can be seen in Fig. 21, 22c, 23c and 24c. Both failure modes (shear and slabbing) may interact and promote each other in a mutual manner. The condition is further exacerbated when k is increased to 2. It is inferred that rockburst is more likely to take place under the condition of higher lateral pressure coefficient. Our results are similar to the previous findings on this topic (Li et al. 2018c; Li and Weng 2016). Fig. 25 plots the released strain energy-time evolution curves at different dip angles and lateral pressure coefficients of unexposed structural plane. Fig. 26 shows the maximum released strain energy distribution at different dip angles and lateral pressure coefficients of unexposed structural plane. It can be observed that, on the one hand, the released strain energy values always remain a lower level (in the range of 10-20KJ) when k is equal or less than 1. However, there is an obvious increase in released strain energy when k is greater than 1. The released strain energy all presents a monotonous increasing trend for a constant dip angle with the increase of k. On the other hand, , the released strain energy presents a trend of initial increase, followed by decreasing values with the increase of dip angle for a constant k. In order to better illustrate the influence of k and structural plane on the failure extent at different locations around the circular tunnel, the effective plastic strain-time evolution curves for various k at different angles of structural plane are plotted in Fig. 27. The maximum effective plastic strain distribution of each monitoring point is also plotted in Fig. 28. For monitoring point 5 (at crown), when k=0.5, the effective plastic strain is always kept at a lower level for various dip angles. This coincides with the stress distribution around a circular tunnel at crown. With the increase of lateral pressure coefficient (k=1, 1.5 and 2), the effective plastic strain begins to increase gradually, and in the case of θ=60º it shows the maximum

value. The sequence of effective plastic strain for a constant λ is always: Sp60º> Sp45º Sp30º Sp90º. For monitoring point 6 (rightside of the tunnel), the effective plastic strain remains at relatively higher level (greater than 0.006) when k=0.5, but it decreases to less than 0.004 when k is increased to 1,

indicating that the plastic strain is restrained due to the changes in in-situ stress conditions. With further

increase in k, the plastic strain values present a monotonous increasing trend. This should be due to progressive fracturing in roof and floor that causes the stress re-distribution on tunnel walls. The higher effective plastic strain values for a constant k always correspond to θ=45º and 60º. For monitoring point 7 (tunnel invert), the values are always 0 when k=0.5. The effective plastic strain for the case of θ=30º, 45º and 90º when k=1 is also equal to 0, except for θ=60º. The different effective plastic strain values between point 5 and 7 under the same stress state should be due to the location of the structural plane and the rock heterogeneity. The effective plastic strain is significantly increased only when k is greater than 1. For monitoring point 8 (left side of the tunnel), no effective plastic strain is observed when k is equal or less than 1. The relatively higher rock strength, compared to rock at the right side of the tunnel (due to the rock heterogeneity), makes the left side of the tunnel relatively stable. When k is increased to 1.5, the effective plastic strain begins to increase, but still remains relatively low at various dip angles. However, an obvious increase can be noticed when k is raised to 2, especially for the case of θ=45º and 60º, indicating

that this area would be fractured under such condition. By observation of Fig. 22d and 23d, it is found that more cracks are generated in the left side of the tunnel, and the failure intensity is even more severe than the right side. Since the structural plane inclines to the left side, the visible cracks in the top left periphery rock are more prone to be generated once λ is raised to higher level. Under such circumstance, the tunnel boundaries near monitoring points 8 are much more disturbed by progressive failure and finally lead to obvious failure, though the left rock strength is higher than the right rock strength.

The above-mentioned analysis and observations suggest that the extent of failure at different locations around the circular tunnel should be assessed based on the analysis of the combined effects of

lateral pressure coefficient, configuration of structural plane, and rock heterogeneity. 4 Conclusions A combined finite/discrete element approach using ELFEN software package was adopted to study the hard rock failure around deep circular tunnel at the presence of structural planes of weakness. A linear relaxation path was used for realistic simulation of the excavation-unloading process, and the rock heterogeneity was also considered in the present study. The dip angle, location and frictional coefficient of the structural plane as well as the lateral pressure coefficient were changed in the modeling to investigate their influence on the failure characteristics of the circular tunnel. Numerical results show that rock

heterogeneity may cause asymmetric failure in the periphery rock of the circular tunnel. The failure modes can be classified into slabbing failure induced by excavation unloading and shear and slip failure induced by the activation of structural plane. However, these two types of failure modes seem to interact in some cases. While the gradual slabbing failure may contribute to the stress concentration and dynamic response around the structural plane, the activated structural plane begins to shear and slip and may in turn further causes slabbing failure near excavation boundaries. Under such condition, the extensive and violent coupled reaction may lead to the serious failure and dynamic response, even the rockburst. Based on the parametric analysis, the following conclusions have been obtained: (1) The favorable dip angle leading to rock failure (crack propagation and nucleation) around the tunnel depend on the frictional coefficient. These two factors are not isolated in controlling the dynamic responses of the underground tunnel. When frictional coefficient is 0.1, the rock failure around the tunnel is most likely when a structural plane at dip angle of 60º is present, while the dip angle of 45º and 30º are more active at frictional coefficient of 0.5 and 1, respectively. The dip angle of 90º has little influence on rock failure, where local slabbing failure composed of fewer visible cracks is only observed near tunnel boundaries. (2) The failure intensity and damage extent with an unexposed structural plane around the periphery of circular tunnel is not so severe compared to the case where the tunnel wall is intersected by an exposed plane of discontinuity. Hence, under the same condition, possibility of encountering rock burst is lower with an unexposed structural plane. (3) The failure intensity and extent of damage (released strain energy) presents a monotonic increase with the increase of lateral pressure coefficient for a constant dip angle. On the other hand, failure intensity follows a trend of initial increase, then decrease with the increase of dip angle for a constant lateral pressure coefficient. Rockburst is more likely under the condition of higher lateral pressure coefficient. Analysis of effective plastic strain in each monitoring point suggests that failure extent at different locations around the circular tunnel is not only influenced by the lateral pressure coefficient, but also depend on the configuration of structural plane and rock heterogeneity. The FDEM modeling is deemed as an effective and convenient way to study the failure processes of rocks subject to different stress state, stress pathways and occurrence of structural planes in underground engineering involving tunnels at high depth. Acknowledgments The authors would like to acknowledge the financial supports from the Key projects of National Natural

Science Foundation (grant number 41630642), and State Key Research Development Program of China (grant number 2016YFC0600706). References Antolini, F., Barla, M., Gigli, G., Giorgetti, A., Intrieri, E., Casagli, N., 2016. Combined finite–discrete numerical modeling of runout of the Torgiovannetto di Assisi rockslide in central Italy. International Journal of Geomechanics 16(6): 1-16. Ashby, M.F., Hallam, S.D., 1986. The failure of brittle solids containing small cracks under compressive stress states. Acta Metallurgica 34(3):497–510. Bobet, A., Einstein, H.H., 1998. Fracture coalescence in rock-type materials under uniaxial and biaxial compression. International Journal of Rock Mechanics and Mining Sciences 35(7):863–888. Cai,.M., 2008. Influence of intermediate principal stress on rock fracturing and strength near excavation boundaries—Insight from numerical modeling. International Journal of Rock Mechanics and Mining Sciences 45: 763-772. Cai, M., 2013. Fracture initiation and propagation in a Brazilian disc with a plane interface: a numerical study. Rock Mechanics and Rock Engineering 46(2): 289–302. Cao, W.Z., Li,X.B., Tao, M., Zhou, Z.L., 2016. Vibrations induced by high initial stress release during underground excavations. Tunnelling and Underground Space Technology 53: 78-95. Dammyr, Ø., Nilsen, B., Gollergger, J., 2017. Feasibility of tunnel boring through weakness zones in deep Norwegian subsea tunnels. Tunnelling and Underground Space Technology 69: 133-146. Elmo, D., Stead, D., Eberhardt, E., Vyazmensky, A., 2013. Applications of Finite/Discrete Element Modeling to Rock Engineering Problems. International Journal of Geomechanics 13(5): 565-580. Farrokh, E., Rostami, J., 2008. Correlation of tunnel convergence with TBM operational parameters and chip size in the Ghomroud tunnel, Iran. Tunnelling and Underground Space Technology 23(6): 700-710. Farrokh, E., Rostami, J., 2009. Effect of adverse geological condition on TBM operation in Ghomroud tunnel conveyance project. Tunnelling and Underground Space Technology 24(4): 436-446. Feng, F., Li, X.B., Li, D.Y., Wang, S.F., 2016. Mechanism and control strategy of buckling rockbursts of orthotropic slab. Chinese Journal of Geotechnical Engineering 39(7): 1302-1311.(in Chinese) Feng, F., Li, X.B., Li, D.Y., 2017. Modeling of failure characteristics of rectangular hard rock influenced by sample height-to-width ratios: A finite/discrete element approach. Comptes Rendus Mecanique 345: 317-328.

Feng, X.T., Chen, B.R., Zhang, C.Q., 2013. Mechanism, warning, and dynamic control of rockburst development processes. Science Press, Beijing (in Chinese) Fu, J.W., Zhang X. Z., Zhu, W.S., Chen, K., Guan, J.F., 2017. Simulating progressive failure in brittle jointed rock masses using a modified elastic-brittle model and the application. Engineering Fracture Mechanics 178: 212-230. Hamdi, P., Stead, D., Elmo, D., 2014. Damage characterization during laboratory strength testing:

A

3D-finite-discrete element approach. Computers and Geotechnics 60: 33-46. He, M.C., Miao, J.L., Feng, J.L., 2010. Rockburst process of limestone and its acoustic emission characteristics under true-triaxial unloading conditions. International Journal of Rock Mechanics and Mining Sciences 47: 286–298. Hedley, D., 1992. Rockburst handbook for Ontario hard rock mines. Ottawa, Canada. James M.A., Jr, J.C.N. 2003. The effect of crack tunneling on crack growth: experiments and CTOA analyses. Engineering Fracture Mechanics 70: 457-468. Klerck, P.A., 2000. The finite element modelling of discrete fracture in quasi-brittle materials. University of Wales Swansea, Wales. Li, D.Y., Li, C.C., Li, X.B., 2011. Influence of sample height-to-width ratios on failure mode for rectangular prism samples of hard rock loaded in uniaxial compression. Rock Mechanics and Rock Engineering 44(3): 253–267. Li, D.Y., Zhu, Q.Q., Zhou, Z.L., Li, X.B., Ranjith, P.J., 2017a. Fracture analysis of marble specimens with a hole under uniaxial compression by digital image correlation. Engineering Fracture Mechanics 183: 109-124. Li., X.B., Cao, W.Z., Zhou, Z.L., Zou, Y., 2014. Influence of stress path on excavation unloading response. Tunnelling and Underground Space Technology 42: 237-246. Li, X.B., Feng, F., Li, D.Y., 2018a. Numerical simulation of rock failure under static and dynamic loading by splitting test of circular ring. Engineering Fracture Mechanics 188: 184-201. Li, X.B., Feng, F., Li, D.Y., Du, K., Ranjith, P.G., Rostami, J., 2018b. Failure characteristics of granite influenced by sample height‑ to‑ width ratios and intermediate principal stress under true‑ triaxial unloading

conditions.

Rock

Mechanics

and

Rock

Engineering

https://doi.org/10.1007/s00603-018-1414-4 Li, X.B., Li, C.J., Cao, W.Z., Tao, M., 2018c. Dynamic stress concentration and energy evolution of

deep-buried tunnels under blasting loads. International Journal of Rock Mechanics and Mining Sciences 104, 131-146. Li, X.B., Weng, L., 2016. Numerical investigation on fracturing behaviors of deep-buried opening under dynamic disturbance. Tunnelling and Underground Space Technology 54: 61–72. Lu, W.B., Yang, J.H., Yan, P., Chen, M., Zhou, C.B., Yi, L., Li, J., 2012. Dynamic response of rock mass induced by the transient release of in-situ stress. International Journal of Rock Mechanics and Mining Sciences 53(9):129–141. Manouchehrian, A., Cai, M., 2017. Analysis of rockburst in tunnels subjected to static and dynamic loads. Journal of Rock Mechanics and Geotechnical Engineering 9: 1031-1040. Martin, C.D., Maybee, W.G., 2000. The strength of hard-rock pillars. International Journal of Rock Mechanics and Mining Sciences 37: 1239–1246. Mitelman, A., Elmo,D., 2014. Modelling of blast-induced damage in tunnels using a hybrid finite-discrete numerical approach. Journal of Rock Mechanics and Geotechnical Engineering 6: 565-573. Mortazavi, A., Molladavoodi, H., 2012. A numerical investigation of brittle rock damage model in deep underground openings. Engineering Fracture Mechanics 90: 101-120. Munjiza A., 2004. The combined finite-discrete element method. John Wiley & Sons, UK. Rockfield, 2013. ELFEN Explicit/Implicit Manual, Version 4.7.1. Rockfield Software Limited, West Glamorgan, UK. Snelling, P.E., Godin, L., McKinnon, S.D., 2013 The role of geologic structure and stress in triggering remote seismicity in Creighton mine, Sudbury, Canada. International Journal of Rock Mechanics and Mining Sciences 58:166-179. Sun, X.M., Chen, F., Miao, C,Y., Song, P., Li, G., Zhao, C.W., Xia, X., 2018 Physical modeling of deformation failure mechanism of surrounding rocks for the deep-buried tunnel in soft rock strata during the excavation. Tunnelling and Underground Space Technology 74: 247-261. Vyazmensky, A., Elmo, D., Stead, D., 2009. Role of Rock Mass Fabric and Faulting in the Development of Block Caving Induced Surface Subsidence. Rock Mechanics and Rock Engineering 43: 533-556. Wang, Y,L., Tang, J.X., Dai, Z.Y., Yi, Ting., 2018. Experimental study on mechanical properties and failure modes of low-strength rock samples containing different fissures under uniaxial compression. Engineering Fracture Mechanics 197: 1-20. Wasantha, P.L.P., Ranjith, P.G., Shao, S.S., 2014. Energy monitoring and analysis during deformation of

bedded-sandstone: Use of acoustic emission. Ultrasonics 54(1): 217-226. Wong, R.H.C., Tang, C.A., Chau, K.T., Lin, P., 2002. Splitting failure in brittle rocks containing pre-existing flaws under uniaxial compression. Engineering Fracture Mechanics 69: 1853-1871. Yang, H.Q., Liu, J.F., Wong, L.N.Y., 2017. Influence of petroleum on the failure pattern of saturated pre-cracked and intact sandstone. Bulletin of Engineering Geology and the Environment https://doi.org/10.1007/s10064-017-1075-7 Yang, J.H., Jiang, Q.H., Zhang, Q.B., Zhao, J., 2018. Dynamic stress adjustment and rock damage during blasting excavation in a deep-buried circular tunnel. Tunnelling and Underground Space Technology 71: 591-604. Zhang, C.Q., Feng, X.T., Zhou, H., Qiu, S.L., Wu, W.P., 2013. Rockmass damage development following two extremely intense rockbursts in deep tunnels at Jinping II hydropower station, southwestern China. Bulletin of Engineering Geology and the Environment 72: 237-247. Zhao, C., Tian, J.S., Song, T.H., Zhao, C.F., Bao, C., 2015. Crack propagation and damage of rock under uniaxial compression based on global strain field analysis. Chinese Journal of Rock Mechanics and Engineering 34(4): 763-769. (in Chinese) Zhao, J.J., Zhang, Y., 2017. Studies on rock failure of layered rock in underground mining-face and control techniques. Geomechanics and Geophysics for Geo-Energy and Geo-Resources 4(3): 405-414. Zhong, Z.B., Deng, R.G., Lv, F., Fu, X.M., Yu, J., 2018. Fracture mechanism of naturally cracked rock around an inverted U-shaped opening in a biaxial compression test. International Journal of Rock Mechanics and Mining Sciences 103: 242-253. Zhou, H., Meng, F.Z., Zhang, C.Q., Hu, D.W., Yang, F.J., Lu, J.J., 2015. Analysis of rockburst mechanisms induced by structural planes in deep tunnels. Bulletin of Engineering Geology and the Environment 74: 1435-1451. Zhou, X.P., Wang, J.H., 2005. Study on the coalescence mechanism of splitting failure of crack-weakened rock subjected to compressive loads. Mechanics Research Communications 32:161-171. Zhou, Z.L., Tan, L.H., Cao, W.Z., Zhou, Z.Y., Cai, X., 2017. Fracture evolution and failure behaviour of marble specimenscontaining rectangular cavities under uniaxial loading. Engineering Fracture Mechanics 184: 183-201. Zhu, W.C., Liu, J., Tang, C.A., Zhao, X.D., Brady, B.H., 2005. Simulation of progressive fracturing processes around underground excavations under biaxial compression. Tunnelling and Underground

Space Technology 20: 231-247.

Failure direction

Intra-element fracturing

Inter-element fracturing

Fig.1 Schematic figure of illustration of crack insertion procedure in ELFEN (Li et al. 2018a)

P0

P0

5

Linear unloading kP0

kP0

20m

kP0

8

6 7

Circular tunnel 20m

Bottom Fixed in y direction Structural plane

Excavation relaxation x

o

Dense grid

y Pt ≠0

Sparse grid

Pt＝0

Fig.2 Model geometry and meshes of the circular tunnel containing structural plane

kP0

(b) Possion’s ratio

(a) Young’s modulus

(d) Tensile strength

(c) Density

(f) Fracture energy

Fig.3 Distribution of various rock parameters in the model (Before excavation unloading)

P

Load factor

Face loading (MPa)

1

0

0.01

0.015

0.02 t (s)

0

(a) Model boundaries

0.01

0.015

(b) Internal tunnel surface

Fig. 4 Loading and unloading paths

0.02 t (s)

(a) k=0.5

(b) k=1

(c) k=1.5

(d) k=2

Fig.5 Maximum principal stress and effective plastic strain contour for various k (without structural plane)

80 70

60

Tangential stress (Mpa)

Tangential stress (Mpa)

70

50

40

30

60 50 40 30

20

20

0

2

4

6

8

0

2

Distance to point 5 along radial direction (m)

(a)

6

8

(b) 70

Ta n g e n t isat lr e s sM(p a)

60

Tangential stress (Mpa)

4

Distance to point 6 along radial direction (m)

50

40

30

60

50

40

30

20

20

0

2

4

6

Distance to point 7 along radial direction (m)

(c)

8

0

2

4

6

8

Distance to point 8 along radial direction (m)

(d)

Fig. 6 The tangential stress distribution along radial direction

350

Released strain energy (KJ)

300 250

k=0.5 k=1 k=1.5 k=2

200 150 100 50 0 0.0050

0.0075

0.0100

0.0125

0.0150

0.0175

0.0200

0.0225

t (s)

Fig. 7 Released strain energy-time evolution curves for various k (without structural plane)

(a) θ=30°

(b) θ=45°

(c) θ=60°

(d) θ=90°

Fig. 8 The maximum principal stress contour line of model containing exposed structural plane at different dip angles

(a) θ=30°

(c) θ=60°

(b) θ=45°

(d) θ=90°

Fig.9 The plastic strain contour line of model containing exposed structural plane at different dip angles

Released strain energy (KJ)

500

400

θ=30° θ=45° θ=60° θ=90°

300

200

100

0 0.005

0.010

0.015

0.020

t (s)

Fig.10 Released strain energy-time evolution curves for different dip angles of exposed structural plane

t=0.0137s

t=0.0149s

t=0.0143s

t=0.02s

(a) μ=0.1

t=0.0138s

t=0.0154s

t=0.0148s

t=0.02s

(b) μ=0.5 Shear and slip failure induced by unloading

Slabbing failure

t=0.0137s

t=0.0151s

t=0.0162s

Coupled extensive failure (slabbing and shear)

t=0.02s

(c) μ=1 Fig.11 Failure process of circular tunnel with various μ at dip angle of 30º

t=0.0131s

t=0.0149s

t=0.0161s

t=0.02s

t=0.0162s

t=0.02s

(a) μ=0.1

t=0.0138s

t=0.0150s

(b) μ=0.5

t=0.0135s

t=0.0149s

t=0.0162s

t=0.02s

(c) μ=1 Fig.12 Failure process of circular tunnel with various μ at dip angle of 45º

t=0.0139s

t=0.0168s

t=0.0152s

t=0.02s

(a) μ=0.1

t=0.0132s

t=0.0163s

t=0.0144s

t=0.02s

(b) μ=0.5

t=0.0135s

t=0.0145s

t=0.0155s

t=0.02s

(c) μ=1 Fig.13 Failure process of circular tunnel with various μ at dip angle of 60º

t=0.0140s

t=0.0149s

t=0.0158s

t=0.02s

t=0.0155s

t=0.02s

(a) μ=0.1

t=0.0136s

t=0.0148s

(b) μ=0.5

t=0.0137s

t=0.0146s

t=0.0158s

(c) μ=1 Fig.14 Failure process of circular tunnel with various μ at dip angle of 90º

t=0.02s

350

μ=0.1 μ=0.5 μ=1

160

Released strain energy(KJ)

Released strain energy(KJ)

200

120

80

40

300 250 200 150 100 50

0 0.005

0.010

0.015

0 0.005

0.020

0.010

0.015

0.020

t (s) (b) θ=45°

t (s) (a) θ=30° 600

50

μ=0.1 μ=0.5 μ=1

500

Released strain energy(KJ)

Released strain energy(KJ)

μ=0.1 μ=0.5 μ=1

400

300

200

100

0 0.005

0.010

0.015

t (s) (c) θ=60°

0.020

40

μ=0.1 μ=0.5 μ=1

30

20

10

0 0.005

0.010

0.015

0.020

t (s) (d) θ=90°

Fig.15 Released strain energy-time evolution curves for various μ at different dip angles

600 Released strain energy (KJ)

θ=30° 500

θ=45° θ=60°

400

θ=90° 300 200

100 0 0

0.25

0.5 Frictional coefficient μ

0.75

1

Fig.16 Maximum released strain energy distribution influenced by μ at different dip angles

(a) θ=30°

(a) θ=30°

(b) θ=45°

(c) θ=60°

(d) θ=90°

Fig.17 The maximum principal stress contour line of models containing unexposed structural plane at different dip angles

(a) θ=30°

(b) θ=45°

(c) θ=60°

(d) θ=90°

Fig.18 The effective plastic strain contour line of model containing unexposed structural plane at different dip angles

Released strain energy (KJ)

350 300 250

θ=30° θ=45° θ=60° θ=90°

200 150 100 50 0 0.0050

0.0075

0.0100

0.0125

0.0150

0.0175

0.0200

0.0225

t (s) Fig.19 Released strain energy-time evolution curves for different dip angles of unexposed structural plane

Released strain energy (KJ)

600 expsosed structural plane unexposed structural plane

500 400

300 200 100 0

20

40

60 Dip angle θ (ﾟ)

80

100

Fig.20 Maximum released strain energy containing exposed and exposed structural plane at different dip angles

t=0.0135s

t=0.0136s

t=0.0124s

t=0.0148s

t=0.0144s

t=0.0163s

t=0.02s

t=0.02s

t=0.02s

(a) k=0.5

(b) k=1

(c) k=1.5

t=0.0126s

t=0.0158s

t=0.02s

(d) k=2

Fig.21 Failure process of circular tunnel with various k at dip angle of 30º containing unexposed structural plane

t=0.0139s

t=0.0136s

t=0.0131s

t=0.0120s

t=0.0155s

t=0.0147s

t=0.0142s

t=0.0149s

t=0.02s

t=0.02s

t=0.02s

(a) k=0.5

(b) k=1

(c) k=1.5

t=0.02s

(d) k=2

Fig.22 Failure process of circular tunnel with various k at dip angle of 45º containing unexposed structural plane

t=0.0139s

t=0.0138s

t=0.0132s

t=0.0124s

t=0.0144s

t=0.0148s

t=0.0152s

t=0.0154s

t=0.02s

t=0.02s

t=0.02s

(a) k=0.5

(b) k=1

(c) k=1.5

t=0.02s

(d) k=2

Fig.23 Failure process of circular tunnel with various k at dip angle of 60º containing unexposed structural plane

t=0.0140s

t=0.0137s

t=0.0127s

t=0.0125s

t=0.0148s

t=0.0144s

t=0.0148s

t=0.0151s

t=0.02s

t=0.02s

t=0.02s

(a) k=0.5

(b) k=1

(c) k=1.5

t=0.02s

(d) k=2

Fig.24 Failure process of circular tunnel with various k at dip angle of 90º containing unexposed structural plane

900

k=0.5 k=1 k=1.5 k=2

400

800

Released strain energy(KJ)

Released strain energy(KJ)

500

300

200

100

700

k=0.5 k=1 k=1.5 k=2

600 500 400 300 200 100

0 0.005

0.010

0.015

0 0.005

0.020

t (s) (a) θ=30°

800

0.015

t (s) (b) θ=45°

0.020

350

k=0.5 k=1 k=1.5 k=2

Released strain energy(KJ)

Released strain energy(KJ)

1000

0.010

600

400

200

300 250

k=0.5 k=1 k=1.5 k=2

200 150 100 50

0 0.005

0.010

t (s) 0.015 (c) θ=60°

0.020

0 0.005

0.010

t (s) 0.015 (d) θ=90°

Fig.25 Released strain energy-time evolution curves for various k at different dip angles containing unexposed structural plane

0.020

Released strain energy (KJ)

1000

θ=30° θ=45° θ=60° θ=90°

800 600 400 200 0 0

0.5

1 1.5 Lateral coefficient k

2

Fig.26 Maximum released strain energy distribution for various λ at different dip angles containing unexposed structural plane

0.018

0.015

t (s)

0.008 0.007 0.006 0.005 0.004 0.003 0.002 0.001 0.000 0.005

0.010

0.015

0.008 0.007 0.006 0.005 0.004 0.003 0.002 0.001 0.000 0.005

Effective plastic strain(ε)

Effective plastic strain(ε)

0.009

0.010

0.015

t (s)

0.008 0.007 0.006 0.005 0.004 0.003 0.002 0.001 0.000 0.005

0.010

0.015

t (s)

0.020

t (s)

0.009

0.006

0.003

0.015

t (s)

0.025 0.020 0.015 0.010 0.005

0.015

0.015

0.012

0.009

0.006

0.003

0.000 0.005

0.010

0.015

t (s)

0.025 0.020 0.015 0.010 0.005

0.010

0.020

0.015

t (s)

0.03

0.02

0.01

0.015

t (s)

0.015

0.010

0.005

(d) θ=90º

0.015

0.015 0.010 0.005 0.000 0.005

0.010

0.015

0.020

0.015

0.020

t (s)

0.07 0.06 0.05 0.04 0.03 0.02 0.01

0.010

0.10

0.08

0.06

0.04

0.02

0.010

0.015

0.020

t (s)

0.020

0.010

0.020

0.00 0.005

0.020

0.025

0.000 0.005

0.025

t (s)

0.04

0.010

0.030

0.00 0.005

0.020

0.05

(c) θ=60º

0.018

t (s)

0.020

0.030

0.00 0.005

0.020

0.015

0.035

(b) θ=45º

0.030

t (s)

0.010

0.040

0.000 0.005

0.020

0.035

0.010

0.005

(a) θ=30º

0.012

0.010

0.010

0.000 0.005

0.020

0.015

0.000 0.005

0.020

Effective plastic strain(ε)

Effective plastic strain(ε)

t (s)

0.015

0.018

0.000 0.005

0.020

0.010

Effective plastic strain(ε)

0.000 0.005

0.020

0.015

Effective plastic strain(ε)

0.010

0.003

0.020

Effective plastic strain(ε)

0.001

0.006

0.025

Effective plastic strain(ε)

0.002

Effective plastic strain(ε)

0.003

0.009

Effective plastic strain(ε)

0.004

0.012

Effective plastic strain(ε)

0.005

B B B B

0.035

0.030

0.015

Effective plastic strain(ε)

0.006

0.000 0.005

Effective plastic strain(ε)

Effective plastic strain(ε)

0.007

Effective plastic strain(ε)

Effective plastic strain(ε)

0.008

0.020

0.035 0.030 0.025 0.020 0.015 0.010 0.005 0.000 0.005

t (s)

0.010

0.015

0.020

t (s)

Fig.27 Effective plastic strain-time evolution curves for various k at different dip angles containing unexposed structural plane. Note: The columns from left to right indicate the case of λ=0.5, 1, 1.5 and 2, respectively.

0.014

λ=0.5

Effective plastic strain(ε)

Effective plastic strain(ε)

0.1

λ=1

0.08

λ=1.5 λ=2

0.06

0.04 0.02

λ=0.5 λ=1

0.012

λ=1.5

0.01

λ=2

0.008

0.006 0.004 0.002 0

0 20

30

40

50

60

70

80

90

20

100

30

40

50

Dip angle θ (ﾟ) (a) point 5

λ=1

0.02

λ=1.5 λ=2

0.015 0.01 0.005 0 30

40

50

60

70

80

0.016

λ=0.5

70

Dip angle θ (ﾟ) (c) point 7

80

90

100

Effective plastic strain(ε)

Effective plastic strain(ε)

0.025

20

60

90

100

Dip angle θ (ﾟ) (b) point 6 λ=0.5 λ=1 λ=1.5

0.012

λ=2 0.008

0.004

0 20

30

40

50

60

70

80

90

100

Dip angle θ (ﾟ) (d) point 8

Fig.28 Maximum effective plastic strain for each monitoring point as a function of λ at different dip angles containing unexposed structural plane

Table 1. Material properties adopted in the present numerical study Name Young’s modulus (E/GPa) Possion’s ratio (μ)

Marble 35.81 (Variance: upper 25%, lower 20%) 0.23 (Variance: upper 15%, lower 5%)

Shear modulus (G/GPa) Density (Ns2/mm4) Tensile strength (σt/MPa)

14.55 2.61e-09 (Variance: upper 20%, lower 20%) 3.1 (Variance: upper 15%, lower 25%)

Plastic strain

0.0

0.001

0.003

0.005

Cohesion (c/MPa)

20

10

2

1

Frictional angle (φ/ﾟ)

32

35

40

45

Fracture energy (Gf/N/mm)

0.02 (Variance: upper 15%, lower 25%)

Normal penalty (Pn/N/mm2)

36000

Tangential penalty (Pt/N/mm2)

3600

Friction of newly generated cracks (γ)

0.5

Frictional coefficient of structural plane (μ)

0.1, 0.5, 1

Mesh element size (m)

0.5 and 0.3

Smallest element size (m)

0.15

Contact type

Node-Edge

Contact field

0.1

Highlights



Simulation of tunnel failure induced by structural plane is used by FDEM approach.



Rock heterogeneity and excavation unloading process was fully considered.



Dip angle, location, frictional coefficient and lateral coefficient were studied.



The crack propagation and failure intensity is related to selected parameters.



Slabbing failure and shear slip failure may interact and promote mutually.