Influence of blasting load directions on tunnel stability in fractured rock mass

Xiohn Li, Zheming Zhu,*, Meng Wng,**, Yun Shu, Shui Deng, Dingjun Xio

a Key Laboratory of Deep Underground Science and Engineering (Ministry of Education), College of Architecture and Environment, Sichuan University, Chengdu,610065, China

b Shock and Vibration of Engineering Materials and Structures Key Laboratory of Sichuan Province, Southwest University of Science and Technology, Mianyang,621010, China

Keywords:Blasting Tunnel stability Crack propagation Failure mode Numerical simulation

ABSTRACT Tunnels in fractured rock masses are typically damaged by dynamic disturbances from various directions.To investigate the influence of blasting load directions on the stability of a tunnel with a precrack nearby, blasting tests were conducted on the physical models of an external crack around a tunnel(ECT)in this study.Failure modes of the tunnels were analysed based on stress wave theory.The Riedel-Hiermaier-Thoma (RHT)material model was employed to perform the numerical simulations on ECT models.Stress distribution around the tunnels and final failure patterns of the tunnels were characterised.The results show that, under blasting loads, the pre-crack propagates and then new cracks initiates on the incident side of the tunnel.These cracks extend towards each other and eventually coalesce.Blasting load directions significantly influence the ultimate failure mode of the tunnel in the fractured rock masses.The new cracks on the shadow side of the tunnel appear at different positions when the blasting stress waves come from various directions.The results are meaningful to the analysis of tunnel stability and optimisation of the tunnel support scheme.

1.Introduction

Tunnels, roadways, and cavities are usually excavated in complex geological environments.The fractured rock mass encountered generally contains many weak planes, such as natural joint and crack, excavation-induced crack, and pass-through fault (Lei et al., 2017; Di et al., 2020; Han et al., 2020; Xu et al., 2020; Zhao et al., 2020), as shown in Fig.1.These weak planes have a great influence on the stability of underground structures.In addition,tunnels may undergo dynamic disturbances induced by impact,earthquake,and blasting(Zhang and Zhao,2013;Shen et al.,2014;Mobaraki and Vaghefi, 2015; Huang et al., 2017).Under these disturbances,especially under blasting loads with high amplitude and strain rate, cracks near the tunnel will initiate, propagate, and eventually lead to tunnel failure, which can seriously threaten the safety of personnel,delay the project progress,and cause economic loss (Zhou et al., 2012; Liu et al., 2017; Weng et al., 2017a).When blasting loads come from various directions,tunnel failure may be characterised by various failure patterns.Therefore,when there is a crack near a tunnel,it is important to identify the failure mode of a tunnel subjected to blasting loads from various directions,which is also crucial for safe construction and support design.

In recent years, researchers have studied the factors that cause rock failure using various methods, including theoretical analyses(Timoshenko and Goodier, 1951; Fraldi and Guarracino, 2010),laboratory experiments (Gong et al., 2019; Luo et al., 2019), and numerical simulations (Zhu et al., 2005; Li et al., 2018a).These researches have revealed that the failure mode of rock is closely related to its mechanical properties and the existence of cavities.However, in most of these studies, it is assumed that the rock is intact,which has not considered the influence of natural or artificial cracks in the actual rock mass.Suorineni et al.(1999) pointed out that,when applying the open stope mining method,faults increase the risk of stope instability and damage.The angle between the fault and the stope is an important factor that affects the stope stability.Bruneau et al.(2003) reported the shaft damage induced by a fault.Jia and Tang(2008)used rock failure process analysis tounderstand the failure mechanism of a tunnel in a fractured rock mass and revealed that the joint angle can significantly influence the failure mode of the tunnel.Huang et al.(2013)conducted model tests and numerical simulations to investigate the effects of weak interlayers on tunnel stability.Their results showed that a weak interlayer can cause asymmetry in the stress distribution and aggravate tunnel failure.

Fig.1.Schematic diagram of a tunnel excavated in a fractured rock mass.

The above studies mainly focused on the analyses of the stability of tunnels or cavities under static or quasi-static loads.In addition,dynamic disturbance is another important factor that can induce rock failure (Martino and Chandler, 2004; Zhu et al., 2010; Wang et al., 2020).Huang and Wang (1999) and Yan et al.(2015)demonstrated that the rockburst around a tunnel is the result of coupled static and dynamic forces.Su et al.(2017) studied the mechanism of rockbursts induced by dynamic disturbances.Weng et al.(2017b) conducted the split Hopkinson pressure bar (SHPB)experiments to study the failure modes of rock specimens containing mini-tunnels and demonstrated that their failure modes are dominated by tensile cracks when the pre-stresses are low.Liu et al.(2014)studied the failure of rock specimens with holes subjected to coupled static and dynamic loads.Zhu et al.(2014) investigated zonal disintegration around a tunnel using AUTODYN under dynamic loads.Li and Weng (2016) used LS-DYNA software to investigate the energy evolutions of the rock mass at different locations around the tunnel under different lateral pressure coefficients,and analysed the effects of the stress wave amplitude and direction on the failure pattern of the tunnel.Li and Li(2018)and Li et al.(2018b)determined that,when a plane wave is incident on a circular tunnel,cracks appear on the incident and shadow sides of the tunnel and extend parallel to the direction of stress wave propagation.These studies(e.g.Liu et al.,2014; Zhu et al.,2014;Li and Weng,2016;Weng et al.,2017b;Li and Li,2018;Li et al.,2018b)showed that,when no cracks exist,newly developed cracks appear around the tunnel or cavity owing to the stress concentration after the stress wave arrives.Nevertheless, when cracks exist near the tunnel or cavity, the stress concentration around the tunnel or cavity can affect the crack propagation path and ultimately influence its failure mode.The cracks initiate and propagate when subjected to dynamic loads.Some researchers used ingenious models that contain cracks and tunnels to study the failure modes of tunnels under dynamic loads(Zhou et al.,2018;Tang et al.,2019).Zhou et al.(2018)found that the failure patterns of cracked tunnels subjected to dynamic loads are quite different from those under static loads.However, the influence of blasting loads on tunnel failure has not been explored in these studies.Blasting loads are not the same as impact loads in terms of loading rate and amplitude(Liu et al., 2020).Regardless of the drop weight test or the SHPB test,the generated wave is only a plane or a one-dimensional(1D)P-wave.The blasting stress wave is a spherical or cylindrical wave,which results in a large difference in the mechanical behaviours of rock between the rock under impact and that subjected to blasting loads.Therefore,studying the failure mode of a tunnel with a crack nearby under blasting loads is of great importance.Cracks initiate and expand under blasting loads and are very likely to propagate to the tunnel face, eventually leading to tunnel failure.In addition,crack propagation behaviours on the incident and shadow sides of the tunnel may exhibit different modes.

In terms of the dynamic responses of rocks subjected to blasting loads,some researchers ascertained crack initiation toughness and identified the crack propagation path and failure mode of rocks(He et al., 2018; Li et al., 2018a; Liu et al., 2020).Using the physical models containing simulated tunnels, Guo et al.(2016) conducted relevant studies on crack propagation behaviour under blasting loads.However, the polymethyl methacrylate material used is relatively isotropic and homogeneous, which may not mimic the mechanical behaviour of an anisotropic and heterogeneous material such as rock under blasting loads.In addition,blasting tests are expensive and difficult to conduct.Most blasting tests have focused on the basic mechanical properties of rock.Although this is also vital to rock engineering, more attention should be paid to the failure mode of tunnels in practical engineering.In particular,when cracks exist near a tunnel, the failure mode of the tunnel under blasting loads is essential to support design.To date,there is a lack of related research, and thus the aim of the present study is to fill this gap.

In this study,the physical models of an external crack around a tunnel(ECT)were built,and blasting tests were conducted on ECT models.Crack propagation gauges (CPGs) were utilised to record the crack propagation data,including the crack initiation time and growth path.The experimental results were analysed based on the theory of stress wave propagation and reflection.The Riedel-Hiermaier-Thoma (RHT) material model parameters of green sandstone were then determined.With the RHT model adopted,LSDYNA software was used to conduct the numerical study, and the reliability of the numerical simulation was verified.The displacement trend line method was employed to identify the failure pattern around a tunnel subjected to blasting loads from various directions.

2.Experimental study

2.1.ECT model and measuring system

To understand the crack propagation behaviour and failure mode of a tunnel under blasting loads that come from various directions when cracks exist near the tunnel, a simplified ECT model was designed to conduct blasting tests,as shown in Fig.2a.Because the stress wave reflected from the model boundary can significantly influence the experimental results,the model should be sufficiently large to reduce this influence.The model size is 600 mm × 600 mm × 15 mm (length × width × thickness).A blasthole with 8 mm in diameter is located at 200 mm and 300 mm away from the left and the lower boundaries of the model, respectively.The length of the pre-crack was 70 mm, and the distance between the left crack tip and the centre of the blasthole is 20 mm.An inverted U-shaped tunnel is excavated in the model with a circular arch radius of 25 mm,sidewall height of 25 mm,and width of 50 mm.The centre of the arch is 365 mm and 300 mm away from the left and lower boundaries of the model,respectively.As shown in Fig.2a,the distances from the blastholecentre to the left, upper, and lower boundaries of the model are 200 mm, 300 mm, and 300 mm, respectively.The distance between the blasthole centre and tunnel centre is 165 mm.The Pwave velocity in the sandstone adopted in this study is 2563 m/s.Therefore,the time for the blasting stress wave to propagate from the blasthole to the left boundary and the time for the reflected stress wave to propagate from the left boundary to the pre-crack tip are longer than that for the blasting stress wave to arrive at the tunnel.This means that the crack initiation and propagation around the area between the pre-crack and tunnel are not affected by the reflected stress wave.

Fig.2.Sketches of(a)the ECT model and(b)the blast loading system.The unit for the dimensions in (a) is mm.

To facilitate the investigation of the influence of blasting load directions on the failure mode of the tunnel,instead of altering the position of the blasthole,the angle between the tunnel symmetrical axis and the horizontal direction (α) is changed.The investigated angles are 0?, 30?, 45?, 60?, 90?, 120?, 135?, 150?, and 180?.Two models are prepared for each case to ensure that at least one test is successful.

The model is made of green sandstone, which has been widely used to investigate the dynamic failure of rocks(e.g.Li et al.,2018a;Zhou et al.,2018;Liu et al.,2020).The mechanical parameters of the green sandstone used in the blasting tests are listed in Table 1.

Table 1Mechanical parameters of green sandstone.

Blasting tests were conducted in the Blast Laboratory at the Southwest University of Science and Technology,as shown in Fig.3.In this test, an electric detonator with a diameter of 6.8 mm was placed at the centre of the blasthole,and the coupling medium was air.A pedestal was placed under the model to prevent the detonator from moving up and down.In addition,a circular rubber ring was used to restrict the detonator from moving left and right.The measuring system included a strain amplifier, a constant voltage supply,a data collection system,and a workstation.A strain gauge(SG)was pasted on the blasthole wall to trigger the data collection system, which means that the data collection system starts to collect data when the strain monitored by the strain gauge exceeds the threshold value.

The data during crack propagation were obtained using CPGs.The CPG used in the blasting tests was BKX3-17.8CY-21-W.As shown in Fig.4,the effective monitoring range is 44 mm×20 mm(length×width),and the resistance is about 3 Ω.A total of 21 wires were used.The distance between two adjacent wires is 2.2 mm.These wires have different widths, and thus their resistances are not identical.The first wire has the lowest resistance,while the 21st wire has the highest resistance.In the experiment,the first wire of the CPG should be aligned with the pre-crack tip and firmly stuck along the crack propagation path.After the detonator was fired,the pre-crack initiated and propagated under the blasting load,and the CPG wires broke.As a result, the total resistance of the CPG changed,and the voltage signals of the CPG monitored during crack propagation varied.Generally,when CPG wires break one by one,a steplike voltage signal is produced.According to these steplike voltage signals,the crack initiation time and the time at which each wire was broken can be determined.The crack propagation velocity can then be calculated based on the distance and breaking time between two wires.

2.2.Blasting test results

Fig.5 shows the blasting test results.It can be seen from the figure that failure modes of ECT models are different when blasting loads come from various directions.The overall failure zones can be divided into three main regions for analysis: the vicinity of the blasthole, the incident side of the tunnel, and the shadow side of the tunnel.

After the detonator was fired, the strength of the shock wave induced by the blasting was extremely high.When its strength exceeds the dynamic compressive strength of the rock, a crushed zone is formed near the blasthole (Kutter and Fairhurst, 1971),resulting in the formation of a cavity with a certain radius,as shown in Fig.5.Because the formation of this zone consumes considerable energy, the shock wave decays quickly.Outside the crushed zone,the shock wave attenuates into a compressive stress wave, and its strength is no longer greater than the dynamic compressive strength of the rock.Therefore, the crushed zone is usually small.Under the compressive stress wave, rock particles are pushed outwards in the radial direction,thereby generating tensile stress in the tangential direction.Because the dynamic tensile strength of the rock is much lower than the dynamic compressive strength,radial cracks are generated.Fig.5 shows that 1-3 major long cracks appeared and propagated to the model boundaries.

2.2.1.Failure modes on the incident side of the tunnel

As can be seen from the red dashed boxes in Fig.5, the crack propagation path is not a straight line, but a curved shape.However,the blasthole is located at the extended line of the pre-crack;consequently, the blasting stress wave propagates parallel to the pre-crack.If there is no tunnel,it can be regarded as a mode-I crack,and the crack propagation path would be almost a straight line(Liu et al., 2020).The crack propagation path is deflected, whichillustrates that the existence of the tunnel has a certain influence on the pre-crack propagation around the tunnel, and this effect is related to the blasting load direction.For these cases, pre-cracks were all initiated, but they propagated to various positions of tunnels when the angle α was different.When α = 0?, i.e.the blasting stress wave propagates perpendicular to the floor of the tunnel,the pre-crack eventually propagates to around the midpoint of the floor.As α increases to 30?, 45?, and 60?, the pre-cracks deflect towards the corner of the tunnel.For the case of α = 45?,the pre-crack finally propagates to the corner of the tunnel, and when α is 30?and 60?,the pre-crack extends to a position that has a certain distance from the corner.For α=90?,i.e.the blasting stress wave propagates vertically to the sidewall of the tunnel, the precrack finally propagates to the spandrel.When α = 120?, 135?,and 150?,the pre-cracks ultimately propagate to the collision point between the extended line of the pre-crack and the tunnel.This also indicates that,for these three cases,the blasting load direction has a limited effect on the failure mode of the ECT model.When α=180?,i.e.the blasting stress wave propagates to the roof of the tunnel, the pre-crack eventually propagates to the roof.The above scenarios are mainly due to the superposition of the incident and reflected stress waves; complicated stress concentration occurs at different positions of the tunnel, resulting in different pre-crack propagations and tunnel failure modes, which will be further analysed in Section 3.

Fig.3.Blasting test and the measuring system.

The voltage signals monitored by the CPG for various values of α are shown in Fig.6.The steplike voltage signal indicates that the CPG wires are broken.The derivative of the voltage signal versus time was calculated,and the time corresponding to its extremum is regarded as the breaking time of the CPG wire, which is also the time when the crack propagated to the position where the CPG wire is broken.Because the distance between adjacent wires of the CPG is fixed(2.2 mm),the corresponding crack propagation velocity can be obtained based on the breaking time of the wires.According to our previous research (Li et al., 2018a; Zhou et al., 2018; Liu et al.,2020), for the propagation of mode-I cracks, CPG wires will be cut from the first wire to the 21st wire to record 21 steps of voltage signals.In the present blasting tests,Fig.5 shows that,for each case,all the CPG wires were cut.If the pre-crack extended steadily from the crack tip to the tunnel face,there were 21 voltage signal steps.However, as α increased from 0?to 180?, the measuring system recorded 10, 8, 7, 11, 9, 14, 15, 14, and 12 steps of voltage signals instead of 21 steps, which indicates that the 21st wire of the CPG broke before the pre-crack expanded to it,i.e.the last broken wire of the CPG is not the 21st wire.Therefore, it can be inferred that,under blasting loads, the pre-crack initiates and propagates, and a new crack will initiate and propagate from the tunnel.Finally,they coalesce with each other where the deflection occurs in the crack propagation path.Because this process is extremely fast and a highspeed camera was not available, the process will be further analysed later in conjunction with the numerical simulations in Section 4.

In addition, Freund (1990) pointed out that the crack propagation velocity cannot exceed the Rayleigh wave velocity, which is lower than the shear wave velocity(Viktorov,1967).Therefore,the crack extension lags behind the stress wave propagation.Related experimental studies (Li et al., 2018a; Liu et al.,2020) have shown that the crack propagation velocities of sandstone under blasting loads are approximately several hundred meters per second.Zhang and Zhao (2014) also pointed out that the crack propagation velocities of sandstone are in the range of 300-650 m/s.The research by Li et al.(2018a) indicated that the average crack propagation velocity for sandstone subjected to a blasting load is 529.3 m/s.The P-wave velocity of the green sandstone used in the experiment was 2563 m/s.It is assumed here that the lag time from when the blasting stress wave reaches the pre-crack tip to its initiation is equal to the lag time from the wave arriving at the tunnel face to the occurrence of the new crack initiated from the tunnel.It takes about 17.2 μs for the P-wave to propagate a distance of 44 mm (the effective length of the CPG),and the pre-crack propagation velocity is conservatively taken as 300 m/s.When the pre-crack initiates and propagates to a distance of 5.16 mm,it can be calculated that a new crack that initiates from the tunnel has not occurred.Therefore, the breaking time of the first three wires of the CPG is accurate,which means that these three wires were cut in sequence as a result of the pre-crack propagation.It is not possible to determine whether the following wires are cut as a result of the propagation of the pre-crack or the new crack initiated from the tunnel.The initiation time of the pre-crack and average crack propagation velocity are listed in Table 2.

Table 2Data for pre-crack propagation.

The above results illustrate that the blasting load directions have only a minor effect on the initiation time of the pre-crack.The precrack is the earliest to initiate for α = 90?, while it is the latest to initiate for α=150?.The difference in the pre-crack initiation time is within 18.6 μs, which may be caused by measurement errors during the blasting tests.The average propagation velocity of thepre-crack was also maintained in a relatively narrow range, but it exhibits no obvious trend.

Fig.4.Sketch of the CPG.

Fig.5.Failure modes of ECT models:(a)α=0?,(b)α=30?,(c)α=45?,(d)α=60?,(e)α = 90?, (f) α = 120?, (g) α = 135?, (h) α = 150?, and (i) α = 180?.

2.2.2.Failure modes on the shadow side of the tunnel

The failure mode on the shadow side of the tunnel exhibits a distinct regularity.When α is in the range of 0?-60?,under blasting loads, new cracks initiate from the tunnel middle, i.e.the collision point between the extended line of the pre-crack and the shadow side of the tunnel,and gradually propagate to the boundaries.As α increased to 90?, a crack appeared at the corner of the tunnel instead.Similar phenomena occur for α = 120?, 135?, 150?, and 180?, where new cracks initiate from the corners of tunnels and propagate to the boundaries, indicating that, under these conditions, a high stress concentration will occur at the corner of the tunnel where is most likely to damage.

3.Blasting stress wave propagation and reflection

The above experimental results can be analysed from the perspective of stress wave propagation and reflection.When the detonator is fired,cylindrical waves are generated in the rock, and when cylindrical waves propagate for a certain distance, they can be regarded as plane waves(Li and Li,2018;Li et al.,2018b).P-and S-waves are induced by blasting;however,to simplify the analysis,only the plane P-wave is considered.

As shown in Fig.7,when the plane P-wave(Pi)reaches the crack tip, it induces the diffracted P-wave (Pd), diffracted S-wave (Sd),Rayleigh wave (Rp), and head wave (De Hoop,1958).Pdand Sdare two cylindrical waves diverging outwards from the crack tip, and the stress induced by them attenuates sharply with the propagation distance.The head wave is a plane shear wave formed by the convergence of Sd,and it has high strength(Yue et al.,2019).Rpis a surface wave that propagates in the area near the crack surface.This can cause the particles near the crack surface to undergo elliptical motion, producing an obvious displacement perpendicular to the crack surface, and thus has a significant influence on crack initiation and propagation.

So far,no clear conclusion can be reached on the effects of these four waves on crack initiation and propagation,and it is difficult to identify the influence of each wave on the behaviour of crack initiation and propagation.Erdogan and Sih (1963) proposed the maximum circumferential tensile stress theory, which converts stresses near a crack tip into radial and tangential stresses, as shown in Fig.8.Thereafter, the crack initiation and expansion behaviour can be clarified by analysing the radial and tangential stresses, which makes it convenient to analyse such problems.When the stress wave is incident parallel to a crack,a symmetrical stress field is generated near the upper and lower surfaces of the crack.Once the circumferential stress exceeds its threshold, the crack starts to initiate and propagate, which explains why, during the early stage of the pre-crack propagation, its path is almost a straight line.

As the stress wave continues to propagate and arrive at the tunnel face, it is reflected.After the reflected and incident stress waves are superimposed, a complex stress field is formed around the incident side of the tunnel.This stress field is analysed as follows.

The detonation pressure P(t) induced by blasting acts on the blasthole wall, and the stress distribution in the rock can be calculated theoretically.Here, we assume that the rock medium conforms to the basic assumptions in elastic mechanics,i.e.the rock is continuous,completely elastic,homogeneous,and isotropic.The governing equation of the stress wave in the medium in polar coordinates can be expressed as follows (Miklowitz,1984):

where a is the radius of the blasthole, r is the distance from the blasthole, φ(r,t) is the displacement potential function, andis the radial stress.P(t) can be obtained by(Blair, 2007):

where PVNis the blasthole pressure after blasting, n and γ are the integers and pressure decay parameters, respectively.

The stress in the medium can then be expressed as

Fig.6.Curves of voltage signals: (a) α = 0?, (b) α = 30?, (c) α = 45?, (d) α = 60?, (e) α = 90?, (f) α = 120?, (g) α = 135?, (h) α = 150?, and (i) α = 180?.

where λ and μ are the Lame constants,represents the tangential stress,andis the normal stress in the out-of-plane direction.

The Laplace transform can be used to solve these equations,and the general solution is as follows:

where m is the Laplace transform parameter; H(t - t′)is the Heaviside function; and K0and K1are the second-type zero-order and first-order modified Bessel functions, respectively.

Durbin (1974) conducted a numerical inversion of Eq.(4).Durbin’s method and Dubner’s method used by Yi et al.(2018) were compared by Tao et al.(2020), who determined that the results obtained by these methods are consistent.

The result of Durbin’s method can be expressed as follows:

where α is any real number (0 ≤α ≤Re(m)) and T is the time interval (0 ≤t ≤T/2).

By using Eqs.(1)-(5), the stress distribution inside the rock under a blasting load was obtained.The radius of the blasthole in our experiment was 4 mm,and the dynamic stress evolution in the green sandstone is shown in Fig.9.

Fig.6.(continued).

Fig.7.Wavefronts when a plane P-wave reaches a crack tip.

Fig.8.Stress distribution around a crack tip.

It can be seen from Fig.9 that, after blasting, the radial and tangential stresses always change from compressive stress to tensile stress, and they attenuate rapidly with time.In the radial direction, the radial compressive stress quickly reaches the extremum and then changes to tensile stress.In addition, in the vicinity of the blasthole, the maximum value of the radial compressive stress was greater than the extremum of the radial tensile stress.With increasing distance from the blasthole, the extremum of the radial compressive stress gradually approaches the maximum value of the radial tensile stress.In the tangential direction, after the tangential stress transforms from compressive to tensile, the extremum of the tangential tensile stress becomes higher than the maximum tangential compressive stress near the blasthole.As the distance from the blasthole increases, the extremum of the tangential compressive and tensile stresses gradually become equal.Therefore, the dominant stress is the compressive stress in the radial direction.Generally, for blasting problems, more attention is paid to the tensile stress in the tangential direction.Fig.9 shows that the peak value of the tangential compressive stress is much lower than that of the radial compressive stress.Therefore,the tensile stress can be regarded as the dominant stress in the tangential direction.This is consistent with the results reported by Yang et al.(2018).

The stress field around the incident side of the tunnel when the blasting stress wave passed through the pre-crack is shown in Fig.10.Here,the cases of α=0?and 120?are selected for analysis;other cases can also be analysed based on the illustration.

The incident stress wave can be decomposed into radial compressive stress and tangential tensile stress, as shown in Fig.10a.For α = 0?, this is the case in which the stress wave is incident perpendicular to the free surface, i.e.the floor of the tunnel.The incident stress wave was reflected off the floor of the tunnel.The reflected stress wave was in the direction opposite to the incident stress wave.The radial compressive stress was reflected as the radial tensile stress,and the tangential tensile stress was reflected as the tangential compressive stress.For cases in which the stress wave is obliquely incident on the tunnel, the propagation direction of the reflected stress wave is related to α.As shown in Fig.10c and d,when α=120?,the left side of the tunnel vault and the left sidewall are the free surfaces.The components of the stress wave in the radial and tangential directions will follow the same law as that for normal incidence of the stress wave,i.e.the radial compressive stress is reflected as the radial tensile stress,and the tangential tensile stress is reflected as the tangential compressive stress.These reflected and incident stress waves will be superimposed, and eventually, a complicated stress concentration will develop around the tunnel.When the superimposed tangential tensile stress exceeds the dynamic tensile strength of the rock, cracks occur around the tunnel.This is also consistent with the studies by Li and Li(2018)and Li et al.(2018b),who found that cracks appear around the tunnel owing to the stress concentration under the plane P-wave.It should be noted that the reflection and diffraction of blasting stress waves at the free surface and the finalsuperimposed stress field are extremely complicated.A certain degree of simplification was performed in the above discussion to qualitatively analyse the problem.

Fig.9.(a) Radial and (b) tangential stresses at different distances from the blasthole.

Fig.10.Stress evolution around the incident side of the tunnel: (a) Propagation of normal incident stress wave, (b) Reflection of normal incident stress wave, (c) Propagation of oblique incident stress wave, and (d) Reflection of oblique incident stress wave.

As explained in Section 2.2.2,crack propagation lags behind the stress wave propagation.Therefore, the pre-crack starts to propagate,and before it reaches the tunnel face,a new crack also initiates from the tunnel as a result of the superposition of stress waves.They eventually coalesce with each other.

In summary, the fracture process of rock around the tunnel starts with the initiation and propagation of a pre-crack outside the tunnel under blasting loads.Then, when the blasting stress wave arrives at the tunnel,a new crack initiates and propagates from the incident side of the tunnel.Finally, the two cracks coalesce with each other.This explains the deflection of the crack propagation path, as shown in Fig.5.

To verify the above analysis, two strain gauges SG1 and SG2 were pasted on the back of the model for the case α = 150?, as shown in Fig.11a.Fig.11b shows the strain histories recorded by the strain gauges.It can be seen that the fracture times of SG1 and SG2 are t1=144.05 μs and t2= 108.05 μs, respectively.SG2 fractured first, which means that, before the pre-crack propagated to the tunnel face, a new crack initiated from the tunnel and propagated towards the incident direction of the blasting stress wave and finally coalesced with the pre-crack.This proves that the analysis above is accurate and reliable.

4.Numerical study

To further analyse the crack propagation behaviour and failure modes of the tunnel under blasting loads, the finite element software LS-DYNA was used for numerical simulation.The model shown in Fig.2, which contains rock, explosive, and air materials,was established.The average mesh size was 1 mm, and the mesh was appropriately refined near the blasthole.To prevent the grid near the blasthole from being deformed significantly after blasting and causing the numerical calculation to terminate incorrectly,the Fluid-Structure Interaction method in LS-DYNA was applied in this simulation.The explosive and air parts are controlled by Arbitrary Lagrangian-Eulerian (ALE) formulation.The rock was set as the Lagrangian part.This method can prevent the rock from deforming significantly during the calculation process to ensure the robustness of the calculation.

4.1.Material models

4.1.1.RHT model for rock

For numerical simulation, it is crucial to determine whether the material model selected is reliable for solving nonlinear,large-deformation problems such as blasting (Chao et al., 2020).LS-DYNA provides many material models that can be used to simulate rock or rocklike materials, such as the Johnson-Holmquist model (Ma and An, 2008), the continuous surface cap model (Tao et al., 2012), the Cowper-Symonds model (Zhao et al., 2017), and the Karagozian-Case model (Liu et al., 2018b,2019).The Cowper-Symonds model is efficient, but may be too simple to characterise the dynamic response of the rock subjected to blasting loads.Determining the parameters of the continuous surface cap model is complicated (Li and Weng,2016).The Johnson-Holmquist model can describe the damage of the material in the compressed state, but it is not good atdescribing the tensile damage of the material (Liu et al., 2018a).In contrast, the Karagozian-Case model can simulate the tensile damage of the material perfectly, but its modelling of the compressive damage of the material is insufficient (Brannon and Leelavanichkul, 2009).Recently, the RHT material model has been widely used to solve rock blasting problems (Borrvall and Riedel, 2011; Xie et al., 2017).It considers compressive and tensile damages, strain rate effect, confining pressure effect, strain hardening, and damage softening effects.Therefore, the RHT model was utilised in this study.

Fig.11.Method for characterising the crack propagation behaviour for α = 150?: (a) Position of strain gauges, and (b) Histories of strains versus time.

Related studies(Yang et al.,2017;Zhao et al.,2017)have shown that rocklike materials exhibit different strain rate effects under compression and tension.Therefore,for the strain rate effect of the RHT model,the strain rate strength factors under compression and tension (Fcr(˙εp) and Ftr(˙εp), respectively) should be calculated separately:

where ˙εpis the strain rate;andare the reference strain rates under compression and tension, respectively (here taken as= 3 × 10-5s-1and= 3 × 10-6s-1); andandare the current strain rates under compression and tension,respectively.γcand γtcan be determined from the continuity requirements;however, they are relatively difficult to obtain.Xie et al.(2017)fitted the data in Zhang and Zhao (2013) and pointed out that,when=30 s-1,γc=0.512;and when=10-1s-1,γt=2.4.βcand βtare the material constants that can be calculated by

where fcis the uniaxial compressive strength (in MPa).βcand βtwere calculated as 0.043 and 0.045, respectively.The strain rate dependence can then be expressed as

where Fr(˙ε) is the strain strength factor, p is the pressure, and ftis the uniaxial tensile strength (in MPa).

The P-α compaction equation of state (EOS) is given by

where P(ρ,e) is the EOS pressure; B0and B1are the material constants; α0is the initial porosity of the material; ρ0is the initial density of the material;e is the internal energy per unit mass;εvisthe volumetric strain; and A1, A2, and A3are the polynomial coefficients.These coefficients can be calculated according to the equations given by Xie et al.(2017):

where Vsis the velocity of sound at ambient pressure and temperature, and M is the material constant.A1, A2, and A3can be calculated as 13.37 GPa,19.68 GPa,and 8.5446 GPa, respectively.

The failure surface parameters A and N can be obtained by

Table 3RHT parameters of green sandstone.

4.1.2.Parameters for the explosive

In our experiments, the No.8 instantaneous electric detonator produced by the Yahua Group was used, and the explosive washexogen.In the simulation, *MAT_HIGH_EXPLOSIVE_BURN was used to model the explosive, and the Jones-Wilkens-Lee (JWL)EOS was adopted to characterise the pressure change during blasting.The JWL EOS is expressed as (Kury et al.,1965):

where P is the detonation pressure; V is the relative volume; E is the internal energy per unit volume;and A,B,R1,R2,and ω are the material constants.The parameters utilised in the simulation,which were obtained from the literature(Banadaki,2010),are listed in Table 4.

4.2.Numerical results

Rock elements can be blanked out when their damage level exceeds a certain threshold to facilitate the analysis of the pre-crack propagation behaviour and failure mode of the tunnel.For this threshold, various values have been obtained by different researchers(e.g.Yang et al.,2017;Yi et al.,2018).The threshold needs to be determined according to the specific experimental results.In this study, the elements whose damage levels exceed 0.4 were blanked out, and the final results are shown in Fig.12.Enlarged views around the tunnels are indicated by red dashed frames.The black dashed lines represent the crack propagation paths.The numerical results are consistent with the experimental results shown in Fig.5.

Under blasting loads, for cases of different values of α, all precracks were first initiated and propagated, and the crack propagation paths are almost straight lines during the early stage, which is consistent with the results shown in Fig.5 and the discussion in Section 3.In other words,during the early stage,the pre-crack initiation and propagation are mainly affected by the incident blasting stress wave.

On the incident side of tunnels,cracks deflect to the nearest corners of tunnels when α=30?,45?,and 60?;whereas,in other cases,crack propagation paths remain at the middle of the tunnels.For α=0?and 180?,the blasting stress wave is incident vertically on the floor and roof of the tunnel,respectively.Owing to the symmetry of the free surface, the stress field after the reflected wave and the incident wave are superimposed should also be symmetrical.Therefore,the crack propagation paths on the incident side of the tunnel under these two conditions are approximately straight lines.However, the experimental results show that the cracks deflect, which may be caused by the heterogeneity of the rock.In addition, for all cases of various α,new cracks were initiated from the shadow side of the tunnel and propagated to the boundaries of the model.When α ≤60?, the new cracks always initiated from the middle of the tunnels,whereas in other cases,they always appeared at the corners of the tunnels.

The case of α=150?is taken as an example to analyse the failure process around the incident side of the tunnel,as shown in Fig.13.Fig.13a shows that the pre-crack initiates at 69.978 μs.At 112 μs,a crack initiated from the tunnel appears, but at this time, the precrack has not yet propagated to the tunnel face.Then, the two cracks propagate towards each other and finally coalesce.This is consistent with the discussion in Section 2.2.2.The initiation time of the pre-crack obtained from the experiment was 73.8 μs,and the fracture time of SG2 was 108.05 μs.The errors between the numerical and experimental results were 5.18% and 3.66%, respectively, which were both within a reasonable range.Therefore, the numerical results fit well with the experimental results.The numerical results of the crack initiation time for different α values are listed in Table 5.The results indicate that the pre-crack initiationtimes for different values of α are essentially the same.It can be inferred that the difference in the initiation time of the pre-crack measured in the experiment may be caused by the inhomogeneity of the rock and the experimental error.The pre-crack initiation was mainly affected by the incident blasting stress wave,and the reflected wave did not influence the crack initiation behaviour.As α increases from 0?to 45?, the time for a crack to initiate from the incident side of the tunnel gradually decreases.When α=45?,the time is the earliest;then,as α increases,the time gradually increases and finally remains stable at about 112 μs.The time for the crack initiated from the shadow side of the tunnelreaches its peak at 309 μs when α=45?,and the difference in the time is not too large for other cases.

Table 4The parameters of the explosive(Banadaki, 2010).

Table 5Numerical results for crack initiation time.

Fig.12.Numerical results for various α: (a) α = 0?, (b) α = 30?, (c) α = 45?, (d) α = 60?, (e) α = 90?, (f) α = 120?, (g) α = 135?, (h) α = 150?, and (i) α = 180?.

Fig.13.Numerical results for α = 150?: (a) Initiation of the pre-crack, (b) Crack initiated from the tunnel, (c) Crack propagation, and (d) Crack coalescence.

Fig.14.Stress distribution around the tunnel in polar coordinates.

Crack initiation and propagation are typically characterised by the maximum circumferential tensile stress theory.In addition,because the tensile strength of the rock is much lower than its compressive strength, it is generally believed that tunnel failure is mainly caused by the tangential tensile stress.Therefore, more attention was paid to the tangential stress distribution around the tunnel.The tangential stress distribution around the tunnel can be viewed by establishing a local coordinate system at the centre of the tunnel,as shown in Fig.14.Fig.15 shows the tangential tensile stress evolution for α = 150?.As shown in Fig.15a, the maximum tangential tensile stress occurs at the pre-crack tip at 69.978 μs,and thus the pre-crack starts to initiate.Then, at 112 μs, the large tangential tensile stress is concentrated around the incident side of the tunnel, which results in the initiation of a new crack from the tunnel.There is always a high tangential tensile stress between the pre-crack and the tunnel,which causes the pre-crack and the crack initiated from the tunnel to propagate towards each other and finally coalesce.In addition, Fig.15 shows that a large tangential tensile stress concentration develops around the corner of the tunnel.Therefore, it can be inferred that the corner is likely to damage on the shadow side of the tunnel.

4.3.Failure modes

The analysis in Section 3 indicates that,after the blasting stress wave reaches the pre-crack, it can be regarded as a pure mode-I crack at the initial stage during its propagation.However, when the blasting stress wave propagates to the tunnel face and is reflected,a complex stress field is formed.The crack propagation was deflected, indicating that the crack pattern also changed.The RHT material model cannot distinguish between tensile damage and compressive damage, which makes it impossible to determine whether the crack pattern is a tensile or shear failure.Therefore,the displacement trend line method proposed by Zhang and Wong(2013) was utilised to confirm the failure mode.According to the displacement trend line,the crack pattern can be divided into three categories: tensile failure, mixed tensile-shear failure, and shear failure, as shown in Fig.16.For the tensile failure of rock, the displacement trend lines are separated from each other,indicating a relative tensile displacement (see Fig.16a).The shear failure of rock is mainly characterised by the displacement trend lines converging towards each other,resulting in a strong shearing effect(see Fig.16c).In addition, the displacement trend lines may have the characteristics of relative tension and shear, which can be considered as mixed tensile-shear failure (see Fig.16b).

The failure modes of rock around the tunnel in the numerical simulation are shown by red dashed frames in Fig.12.It is clear that,at the initial stage of pre-crack initiation and propagation,its failure mode is manifested as a tensile failure, which means that, at thisstage, crack initiation and propagation are mainly affected by the incident stress wave, while the reflected wave has limited effects.As the crack continues to propagate forwards, the effect of the reflected wave gradually increases.It can be concluded that, when α=0?,45?,90?,120?,150?,and 180?,the crack propagation on the incident side of the tunnel exhibits a tensile failure mode.For the case of α = 30?, the incident side of the tunnel presents shear failure.When α=60?,the same zone exhibits mixed tensile-shear failure.Similar phenomena were also observed for α = 135?.Nevertheless,all new cracks that initiated from the shadow side of the tunnels were characterised as tensile failure.

Fig.15.Evolution of the tangential tensile stress for α = 150?: (a) t = 69.978 μs, (b) t = 112 μs, (c) t = 115.99 μs, and (d) t = 118.97 μs.

4.4.Stress analysis

To clarify the danger level at different positions of the tunnel for cases of various α values, such as the roof, spandrel, sidewall,corner, and floor, stresses at monitoring points A-G, as shown in Fig.14, were derived.By establishing a polar coordinate system at the centre of the tunnel, the tangential stress σθθat any position around the tunnel can be obtained by

where σθθis the tangential stress;σx,σy,and τxyare the x-stress,ystress and shear stress,respectively; and θ is the polar angle.

In LS-DYNA, a negative stress value represents a compressive stress, whereas a positive value represents a tensile stress.The extreme values of the tangential tensile stress at different positions for identical α values are shown in Fig.17.The figure illustrates that the extreme values of the tangential tensile stress around the tunnel appear at different locations for a certain α,which means that, when the blasting load comes from a certain direction, danger levels at different locations of the tunnel differ.When the blasting stress wave passes through the pre-crack andis incident perpendicular to the floor of the tunnel (α = 0?),position A at the shadow side of the tunnel is most likely to be damaged.For α = 30?, failure is prone to occurring at position E2on the incident side of the tunnel.When α = 45?, position B1is most likely to be damaged.In addition, high tangential tensile stresses appear at A, C1, D1, and E2; therefore, these positions are also prone to fracturing and require more attention.For α = 60?,the dangerous places are E2, B1, C1, D1, and E1.When α = 90?, the failure is most likely to occur at C2,and B1,C1,D1,and E1are likely to be damaged.As α increases to 120?,cracks can easily appear at C2and E1.When α = 135?, E1on the shadow side of the tunnel is most likely to be damaged.There is also a high tangential tensile stress at B2on the incident side of the tunnel,which is also likely to be damaged.For the case of α = 150?, the final failure mode is similar to that at α =135?, i.e.positions E1and B2are most likely to be damaged.When α=180?,failure easily occurs at position A on the incident side of the tunnel.

Fig.16.Three failure modes based on displacement trend line method (Zhang and Wong, 2013): (a) Tensile failure, (b) Mixed tensile-shear failure, and (c) Shear failure.

Fig.17.Extreme tangential tensile stresses at different positions for identical α values: (a)α = 0?, (b)α = 30?, (c)α = 45?, (d)α = 60?, (e)α = 90?, (f)α = 120?, (g)α = 135?, (h)α = 150?, and (i) α = 180?.

Fig.18 illustrates that danger levels at the same location of the tunnel are different when blasting loads come from various directions.For the tunnel roof A, there is a great tangential tensile stress when α=0?and 180?,and damage is most likely to occur.At position B1, damage is likely to occur when α = 45?, 60?, and 90?.For position C1, the extreme value of the tangential tensile stress gradually increases as α increases from 0?to 90?,and after reaching its maximum value when α = 90?, it then gradually declines with increasing α.Position D1is most likely to be damaged when α=60?.Position E1is most damage-prone when α ranges from 60?to 150?,since there is a high tangential tensile stress.At position F1, when α = 135?, damage is most likely to occur.Position G is damageprone for cases of α = 0?, 30?,135?,150?, and 180?.Failure easily occurs at position F2when α = 30?.Position E2is more damageprone when α = 30?, 45?, and 60?.Positions D2and C2are most likely to be fractured when α = 90?.Failures may easily appear at position B2when α = 120?and 135?.

Based on the above analyses, the position that is most likely to be damaged when the blasting load comes from a certain direction can be determined.The tunnel can then be supported in a costeffective manner.

5.Conclusions

Tunnels in fractured rock masses are usually affected by dynamic disturbances from various directions, which will have a significant influence on the stability of tunnels.To understand the failure modes of tunnels under blasting loads when there are cracks near the tunnel, ECT models were built in this study.A series of blasting tests was conducted by applying the ECT models and the CPG monitoring system.The final failure modes of the ECT modelsare discussed based on stress wave theory.By adopting the RHT material model, a numerical simulation was performed using LSDYNA to verify the experimental results and to clarify the stress distribution around the tunnel.The displacement trend line method was used to identify failure patterns.From this study, the following conclusions can be drawn:

Fig.17.(continued).

(1) Under blasting loads, the pre-crack first initiates and propagates,and a new crack initiates from the incident side of the tunnel.The two cracks propagate towards each other and finally coalesce.At the early stage, the pre-crack propagates as a mode-I crack, and the blasting load direction generally has no effect on the initiation time of the pre-crack.

Fig.18.Extreme tangential tensile stresses for various α values at an identical position(monitoring point):(a)A,(b)B1,(c)C1,(d)D1,(e)E1,(f)F1,(g)G,(h)F2,(i)E2,(j)D2,(k)C2,and(l) B2.

(2) The incident direction of the blasting load significantly influences the tunnel stability.For the incident side of the tunnel,when α = 30?, 45?, and 60?, pre-crack propagation is always deflected towards the nearest corners of the tunnels, and propagates to the middle of the tunnels for other angles.New cracksinitiateontheshadowsidesofthetunnels.Whenα<90?,the new crack always initiates at the middle of the tunnels,whereas for α ≥90?,it appears at the corner of the tunnels.

Fig.18.(continued).

(3) The RHT material model used in this study can well characterise the failure process of green sandstone subjected to blasting loads,and the numerical results are consistent with the experimental results.

(4) On the incident side of the tunnel, crack propagation is characterised as tensile failure for cases α=0?,45?,90?,120?,150?,and 180?.When α=30?,the failure pattern of the crack is shear failure, whereas for α = 60?and 135?, it exhibitsmixed tensile-shear failure.On the shadow side of the tunnel,new cracks can be characterised as tensile failures for various values of α.

(5) As the blasting stress wave comes from various directions,the danger levels at various positions of the tunnel are diverse, depending on the tangential tensile stress around the tunnel.The greater the tangential tensile stress is, the more likely the damage to occur.

Declaration of competing interest

The authors wish to confirm that there are no known conflicts of interest associated with this publication and there has been no significant financial support for this work that could have influenced its outcome.

Acknowledgments

This research was funded by the National Natural Science Foundation of China(Grant No.U19A2098), the open fund of MOE Key Laboratory of Deep Underground Science and Engineering(Grant No.DESEYU202101), and the Sichuan Science and Technology Program (Grant No.2021YJ0511).