Effect of Different Tunnel Distribution on Dynamic Behavior and Damage Characteristics of Non-Adjacent Tunnel Triggered by Blasting Disturbance

: When a blasting is executed near two tunnels, the blasting wave will trigger a dynamic response and damage to the tunnels. Depending on the tunnel distribution, the path of the blasting wave to the remote non-adjacent tunnels will change. The aim of this study is to analyze the effect of the tunnel distribution on the dynamic response characteristics of a remote non-adjacent tunnel. Numerical models of two tunnels were established by PFC2D and three different tunnel distributions were considered. The two tunnels were divided into the adjacent tunnel and the non-adjacent tunnel according to their relative distance to the blasting source. The dynamic stress evolution, damage characteristics and the evolution of strain energy of the non-adjacent tunnel were initially analyzed. The results show that the stress wave amplitude of the non-adjacent tunnel is closely related to the tunnel distribution, but only near the sidewalls of the non-adjacent tunnel is the stress wave waveform sensitive to the tunnel distribution. The larger the tunnel dip, the more severe the damage to the non-adjacent tunnel. In addition, as the tunnel dip increases, the maximum strain energy densities (SEDs) in the roof, ﬂoor and sidewalls of the non-adjacent tunnel exhibit different trends. The inﬂuence of the wavelength of the blasting wave is further discussed. It is shown that the dynamic stress ampliﬁcation factor and damage degree around the non-adjacent tunnel is usually positively correlated with the wavelength of the blasting wave. Moreover, the release of strain energy around the non-adjacent tunnel has a positive correlation with the wavelength. The SED variations in different areas around the non-adjacent tunnel also exhibit different trends with the increase of tunnel dip.


Introduction
In underground engineering excavation, a large number of tunnels are excavated in a limited area to reduce the workload and cost [1][2][3]. These tunnels are important structures to ensure the safety of underground engineering. However, with the increase of excavation depth and the number of engineering structures, the stress concentration around the engineering structures becomes more obvious and a large amount of strain energy is stored around them. These initial strain energies can induce a high risk of rock burst, collapse and fracture [4][5][6]. In addition, the dynamic disturbance in underground engineering is another major cause of engineering disasters [7][8][9][10]. In underground tunnels which suffer different distribution, the propagation of stress waves is also complex and diverse, which makes it very difficult to predict and prevent disasters. Therefore, it is necessary to understand the dynamic behavior and damage characteristics of underground tunnels with different distributions.
In recent decades, the static mechanical behavior of underground tunnels has attracted extensive attention. Some scholars have calculated the stress distribution around the tunnel by the elastic mechanics and complex function methods [11,12]. They generally believe that the lateral pressure coefficient and the tunnel cross-section shapes are important factors affecting stress distribution. For example, Kirsch first obtained the stress distribution function (Kirsch's solution) around the circular tunnel under different lateral pressure conditions based on elastic mechanic theory [12]. Subsequently, the stress distribution formula of the elliptical tunnel was obtained by using the conformal mapping method of complex function and the classical Kirsch's solution [12]. Recently, some scholars have also solved the stress distribution of specific shape tunnels (such as a rectangular and semi-circular tunnel) using numerical regression analysis, complex function theory and numerical simulation [13,14]. For example, Exadaktylos and Stavropoulou [13] adopted the complex function method to calculate the stress distribution of multiple shapes of underground tunnels with rounded corners and further verified the theoretical solution by the FLAC3D numerical simulation. Zhao et al. [14] also obtained the initial function for solving the stress distribution around the rectangular tunnel by defining a coefficient related to the height-width ratio. On the other hand, many scholars have studied the failure mechanism of underground tunnels [15][16][17][18]. For example, Zhu et al. [15] analyzed the influence of lateral pressure coefficient on the failure of a U-shaped tunnel based on RFPA numerical simulation. The results showed that, when the lateral pressure coefficient is 1, the roof and floor would suffer shear failure, while when the lateral pressure coefficient is 4, the sidewalls would suffer tensile failure. Gong et al. [16] also studied the rock burst mechanism of hard rock tunnels with a circular cross-section. They found that the rock burst process of deep hard rock tunnels has a typical time effect which can be divided into four stages and the spalling can further be developed into the rock burst. Si et al. [17] conducted a series of triaxial compression tests and investigated the spalling mechanism of sidewalls of D-shaped tunnels. The results show that the spalling will be inhibited and the depth of the V-notch will be reduced under higher lateral pressure.
Recently, some scholars have paid attention to the dynamic mechanical behavior and failure characteristics of underground tunnels [19][20][21][22][23]. In general, the dynamic disturbance or unloading disturbance can prompt the rapid deformation and energy conversion of surrounding rock in the form of stress waves. Li et al. [19] have indicated that the release of strain energy is closely related to the failure modes of underground tunnels. The result shows that, when roof spalling is induced, the release of strain energy will last for a long time, but when a violent strain rock burst occurs, massive strain energy will be released instantaneously. Si et al. [20] conducted the triaxial unloading compression test and investigated the strength-weakening effect of unloading and unloading rate on finegrained granite. They found a lower unloading rate more conducive to the improvement of the bearing strength and the storage of elastic energy. In addition, some scholars have also analyzed the influence of disturbance location, in-situ stress, disturbance amplitude and disturbance duration on the dynamic stability of the underground tunnel [24][25][26][27]. For example, Qiu et al. [25] carried out a series of physical model tests on deep tunnels and they found that, even when the disturbance distance is the same but the disturbance dip is different, the dynamic stability of deep tunnels varies greatly. Li et al. [26] studied the effect of stress wave wavelength on the failure model of the deep tunnel and found that the short stress wave tends to cause spalling and the length stress wave tends to cause rock burst. Zhu et al. [27] analyzed the failure mechanism of the deep tunnel under different disturbance amplitudes. The results show that, the larger the disturbance amplitude, the more serious the dynamic failure. In addition, Kulynych et al. [28] studied the action process of gaseous products of explosive on rock fracturing behavior and found that the effect of rock mass disturbance would gradually decrease with increased relaxation of defect formation. Slashchov et al. [29] analyzed the relationship between the emanation activity of radon decay products in mining tunnels and the geological dislocations and believed that the hidden tectonic disturbances and high-stress concentrations can be determined by monitoring the gas-dynamic processes in mining tunnels. Arnau et al. [30] compared the dynamic behavior of a double-decker circular tunnel under train-induced loads with that of a simple tunnel and the results showed that the soil response would be significantly different across the frequency range studied. Behshad et al. [31] analyzed the influence of Dynamic Vibration Absorbers (DVA) on the vibration reduction of a double-deck circular railway tunnel and found that DVA can effectively reduce the total energy flow acting on the tunnel when trains pass by.
In underground engineering, numerous tunnels are usually excavated in a limited area, as shown in Figure 1. These tunnels are locally distributed horizontally, vertically and obliquely. During the excavation or internal blasting of these tunnels, existing tunnels will inevitably be affected by the adjacent working area, including the redistribution of static stress and the dynamic response. Usually, the impact of blasting disturbance in the adjacent working area is particularly significant. Numerous previous studies have confirmed the impact of blasting disturbance or dynamic disturbance on adjacent tunnels, which usually causes local stress surges, triggering rock bursts or surrounding rock spalling [19,[24][25][26]. The object of these studies is usually a single tunnel or the tunnel near the disturbance source, but little attention has been paid to multiple tunnels, especially the remote nonadjacent tunnels. Limited studies have analyzed the dynamic behavior of multiple tunnels under dynamic or unloading disturbances [32,33]. For example, Feldgun et al. [33] analyzed the dynamic behavior of a rectangular existing tunnel induced by the internal blasting inside another horizontal parallel tunnel. The results show that the left sidewall of the existing tunnel will bend under the blasting disturbance. Li et al. [32] also analyzed the effect of unloading waves caused by neighboring tunnel excavation on the existing tunnel. The results show that the unloading wave can cause a strong dynamic response of the existing tunnel and the increase of the unloading rate can amplify the dynamic effect. However, due to the multiple transmission and reflection of stress waves among multiple tunnels, the dynamic stability and stress concentration of the non-adjacent tunnel will not only depend on the amplitude and wavelength of the initial incident wave, but also depend on the interaction with the adjacent tunnel. Therefore, the tunnel distribution is also an important factor, which can affect the dynamic response and fracturing behavior around the tunnels triggered by blasting disturbance [34,35]. For this reason, the main object of this study is to analyze the influence of tunnel distribution on the dynamic response and stability of the non-adjacent tunnel under blasting disturbance. Therefore, the numerical models with the two tunnels were established by the particle flow code (PFC2D) method and different tunnel distributions were also considered. In addition, the stress evolution, energy evolution and failure characteristics around the two tunnels caused by blasting disturbance were analyzed. The effect of stress wave wavelength on the dynamic behavior of a non-adjacent tunnel is further discussed.

Description of PFC2D
The discrete element method PFC2D is widely used in rock engineering because it has the advantage of synchronous microcracks display and does not need to consider the convergence issue in calculations. PFC2D provides several typical modeling methods to simulate different mechanical properties of materials, such as the contact bond (CB) model, parallel bond (PB) model, smooth joint (SJ) model, flated joint (FJ) model, etc. Generally, the CB model is used to simulate the soil material or the bulk material, because it can only transfer the force and not moments through the contact between particles. The PB model is suitable for rock-like materials because it can effectively transmit forces and moments [37,38]. The SJ model is usually used to simulate structural planes such as joints, cracks and bedding-in materials. FJ mode is an increasingly popular new modeling method for rock-like materials because it provides a more realistic ratio between the tensile strength and compression strength of materials. However, because the FJ model will generate too many micro-cracks, it also has certain disadvantages in observing the failure mode of the rock mass. Therefore, the PB model was chosen to simulate the mechanical behavior of rock mass in this study. As shown in Figure 2a, a series of non-uniform-sized rigid particles can be seen as the basic constituent element of the PB model, the linear contact is arranged between the particles and a parallel bond can be conceived as a finitesize concrete cemented between particles. The mechanical behavior between particles can be assumed to be produced by a series of linear springs, viscous dashpots and two specific boned elements (with the normal and shear strengths c  and c  ) [38], as shown in Figure  2b.

Description of PFC2D
The discrete element method PFC2D is widely used in rock engineering because it has the advantage of synchronous microcracks display and does not need to consider the convergence issue in calculations. PFC2D provides several typical modeling methods to simulate different mechanical properties of materials, such as the contact bond (CB) model, parallel bond (PB) model, smooth joint (SJ) model, flated joint (FJ) model, etc. Generally, the CB model is used to simulate the soil material or the bulk material, because it can only transfer the force and not moments through the contact between particles. The PB model is suitable for rock-like materials because it can effectively transmit forces and moments [37,38]. The SJ model is usually used to simulate structural planes such as joints, cracks and bedding-in materials. FJ mode is an increasingly popular new modeling method for rock-like materials because it provides a more realistic ratio between the tensile strength and compression strength of materials. However, because the FJ model will generate too many micro-cracks, it also has certain disadvantages in observing the failure mode of the rock mass. Therefore, the PB model was chosen to simulate the mechanical behavior of rock mass in this study. As shown in Figure 2a, a series of non-uniform-sized rigid particles can be seen as the basic constituent element of the PB model, the linear contact is arranged between the particles and a parallel bond can be conceived as a finite-size concrete cemented between particles. The mechanical behavior between particles can be assumed to be produced by a series of linear springs, viscous dashpots and two specific boned elements (with the normal and shear strengths σ c and τ c ) [38], as shown in Figure 2b. For the linear contact, the moment is zero and the contact force can be iterated by: where the subscripts n and s denote normal and shear direction, respectively. The subscripts l and d denote linear spring and dashpot contributions, respectively. It is worth noting that the contact activation state is determined by the comparison between the contact gap (gc) and the reference gap (gr). When gc > gr, the linear contact is inactivated and the calculations of force-displacement is ignored. When gc  gr, the linear contact is activated and the force components   For the linear contact, the moment is zero and the contact force can be iterated by: where the subscripts n and s denote normal and shear direction, respectively. The subscripts l and d denote linear spring and dashpot contributions, respectively. It is worth noting that the contact activation state is determined by the comparison between the contact gap (g c ) and the reference gap (g r ). When g c > g r , the linear contact is inactivated and the calculations of force-displacement is ignored. When g c ≤ g r , the linear contact is activated and the force components F l n , F l s , F d n , F d s of linear contact can be upgraded by the following equations: where (F l s ) 0 denotes the shear force of linear contact at the last step. ∆U s represents the relative shear displacement increment. µ denotes the friction coefficient. m c denotes the effective mass of linear contact. v n and v s denote the relative normal velocity and shear velocity, respectively. k n and k s denote the normal and shear stiffness of linear contact, respectively.
For the parallel bond, there is a linear relationship between force and displacement. In addition, the bending moment M b is also calculated. The force components F n , F s and the bending moment M b can be iterated by: where A is the area of the parallel bond cross-section. I is the moment of inertia of the parallel bond. ∆θ b is the relative rotation increment. k n and k s denote the normal and shear stiffness of parallel bond, respectively. The maximum normal and shear stresses acting on the parallel bond can be calculated: where R is the radius of the bond cross-section. The failure state of the parallel bond can be determined by a comparison between these stresses and the tensile and shear strengths where the tensile strength is preset manually and the shear strength can be updated by the Mohr-Coulomb criterion: where c is the cohesion and φ is the friction angle.

Modelling Procedure and Calibration of PB Model Parameters
Calibration of the numerical model should be performed to obtain appropriate microparameters of the PB model for matching macroscopic behaviors between the numerical model and experiment. In this work, the granite specimens from Linglong gold mine were tested by the uniaxial compression method and the corresponding numerical specimens were established to calibrate these results. A series of trial and error procedures are executed by adjusting the microscopic parameters of the numerical model so that the macroscopic mechanical characteristics of the numerical model, including UCS, elastic modulus and Poisson, have a good match with that of actual granite. The empirical uniaxial compression strength (UCS) is about 158.45 MPa, the elastic modulus is 32.3 GPa and the Poisson's ratio is 0.258. The corresponding numerical results are listed in Table 1. As shown in Table 1, the error between the numerical results and the real granite is less than 5%, indicating that the calibrated numerical model can well represent the real granite. The stress-strain curves and failure modes of the numerical model and granite are also compared in Figure 3 [24]. These results also show that the numerical model is in good agreement with the testing result The obtained calibrated micro-parameters of the PB model are listed in Table 2.
The in-situ stress of Linglong Gold Mine was considered in this work. The tunnels are assumed to be excavated along the direction of minimum horizontal principal stress. Thus, the linear functions of the maximum horizontal principal stress and vertical principal stress are presented as follows [39] σ hmax = 0.4612 + 0.0588h (13) σ v = −0.4683 + 0.0316h (14) where σ hmax and σ v are the maximum horizontal principal stress and vertical principal stress, respectively. h is the depth. The numerical model with dimensions of 48 m × 24 m was established by the calibrated PB model, as shown in Figure 4. The radius of particles is in the range of 0.06 m-0.096 m. The viscous boundary condition is set to reduce the reflected wave at the boundary [40]. The depth of 1200 m was investigated. Correspondingly, the horizontal and vertical stresses applied to the model boundaries are 71.02 MPa and 37.45 MPa, respectively. Two tunnels excavated before these boundary stresses were loaded. The tunnel modelling process in this study involves the main factor: tunnel distribution. As shown in Figure 4, two circular tunnels with the same diameter (4 m) are considered. The distance between the two tunnel centers is set to 8 m. The tunnel dip β is defined as the angle between the two tunnel centers and the horizontal direction. As shown in Figure 4a, the tunnel distribution can be obtained by changing the tunnel dip. In this study, three different tunnel dips were considered: 0 • , 45, • and 90 • .
Two tunnels excavated before these boundary stresses were loaded. The tunnel modelling process in this study involves the main factor: tunnel distribution. As shown in Figure 4, two circular tunnels with the same diameter (4 m) are considered. The distance between the two tunnel centers is set to 8 m. The tunnel dip  is defined as the angle between the two tunnel centers and the horizontal direction. As shown in Figure 4a, the tunnel distribution can be obtained by changing the tunnel dip. In this study, three different tunnel dips were considered: 0, 45, and 90.  The real blasting wave is extremely complex and it is almost impossible to completely reconstruct it through numerical simulation. For this reason, some scholars proposed the simplified triangular wave method to replace the real blasting wave [24,41,42]. Their results show that this simplified triangular wave can effectively reflect the effects of blasting. To this end, the triangular wave method was adopted in this study. Before the blasting, a borehole with a radius of 0.3 m was first excavated, then a triangle stress wave was applied to the borehole periphery, as shown in Figure 5. The tunnel close to the borehole can be regarded as the adjacent tunnel and the tunnel far away from the borehole is the nonadjacent tunnel. The distance between the borehole location and the adjacent tunnel is set as 4 m. The peak stress of the triangle stress wave is 3 GPa, the rising time tr is 250 s and the total time tm is 1250 s. The time ratio k = tm/tr, which is defined as the ratio of rising time to total time, is set to 5. Four measuring circles with a radius of 0.3 m are arranged around the non-adjacent tunnel. The measuring circles near the non-adjacent tunnel are named B1, B2, B3 and B4, as shown in Figure 4b. The real blasting wave is extremely complex and it is almost impossible to completely reconstruct it through numerical simulation. For this reason, some scholars proposed the simplified triangular wave method to replace the real blasting wave [24,41,42]. Their results show that this simplified triangular wave can effectively reflect the effects of blasting. To this end, the triangular wave method was adopted in this study. Before the blasting, a borehole with a radius of 0.3 m was first excavated, then a triangle stress wave was applied to the borehole periphery, as shown in Figure 5. The tunnel close to the borehole can be regarded as the adjacent tunnel and the tunnel far away from the borehole is the non-adjacent tunnel. The distance between the borehole location and the adjacent tunnel is set as 4 m. The peak stress of the triangle stress wave is 3 GPa, the rising time t r is 250 µs and the total time t m is 1250 µs. The time ratio k = t m /t r , which is defined as the ratio of rising time to total time, is set to 5. Four measuring circles with a radius of 0.3 m are arranged around the non-adjacent tunnel. The measuring circles near the non-adjacent tunnel are named B1, B2, B3 and B4, as shown in Figure 4b.

Dynamic Stress Characteristics
When a blast is executed, the nearby surrounding rock will directly break into pieces and the far-field surrounding rock will accumulate or release deformation under the attenuated stress wave, causing unexpected damage. Generally, the stress evolution of surrounding rock induced by the blasting disturbance is mainly manifested in two aspects: dynamic stress and static stress. In PFC2D, the stress recorded by the measuring circle actually includes the static part and dynamic part. Figure 6 presents the typical stress-

Dynamic Stress Characteristics
When a blast is executed, the nearby surrounding rock will directly break into pieces and the far-field surrounding rock will accumulate or release deformation under the attenuated stress wave, causing unexpected damage. Generally, the stress evolution of surrounding rock induced by the blasting disturbance is mainly manifested in two aspects: dynamic stress and static stress. In PFC2D, the stress recorded by the measuring circle actually includes the static part and dynamic part. Figure 6 presents the typical stresstime curve recorded by measuring circles. It can be seen from Figure 6 that the static stress before blasting (t < 0 µs) is usually a stable value and the static stress after blasting will be stabilized again. For expedient analysis, the dynamic part is separated from the superimposed curve in the subsequent sections. Furthermore, it should be noted that the tensile stress is positive and the compressive stress is negative.

Dynamic Stress Characteristics
When a blast is executed, the nearby surrounding rock will directly break into pieces and the far-field surrounding rock will accumulate or release deformation under the attenuated stress wave, causing unexpected damage. Generally, the stress evolution of surrounding rock induced by the blasting disturbance is mainly manifested in two aspects: dynamic stress and static stress. In PFC2D, the stress recorded by the measuring circle actually includes the static part and dynamic part. Figure 6 presents the typical stresstime curve recorded by measuring circles. It can be seen from Figure 6 that the static stress before blasting (t < 0 s) is usually a stable value and the static stress after blasting will be stabilized again. For expedient analysis, the dynamic part is separated from the superimposed curve in the subsequent sections. Furthermore, it should be noted that the tensile stress is positive and the compressive stress is negative. In the periphery of the tunnel, the radial stress is usually very small, so only the tangential stress is discussed in this study. In this section, the dynamic tangential stress was examined. Figure 7 presents the tangential stress waves around non-adjacent tunnels under different tunnel dips. As shown in Figure 7, these stress curves generally undergo one or multiple peaks, of which the compression peak is less than zero and the tensile peak is In the periphery of the tunnel, the radial stress is usually very small, so only the tangential stress is discussed in this study. In this section, the dynamic tangential stress was examined. Figure 7 presents the tangential stress waves around non-adjacent tunnels under different tunnel dips. As shown in Figure 7, these stress curves generally undergo one or multiple peaks, of which the compression peak is less than zero and the tensile peak is greater than zero. It can be found these peaks exhibit different variation trends at different locations, as follows: (1) For zone B1, when β = 0 • and 45 • , the maximum tensile peaks are generally greater than the maximum compression peaks, so attention should be paid to the tensile failure of surrounding rock in the vicinity of this zone. When β = 90 • , the tensile peak is not obvious and the maximum compression peak is 180.4 MPa, which is far greater than the maximum tensile stress. The result shows that, when β = 90 • , the compression failure tends to occur near zone B1. (2) For zone B2, the maximum tensile peak is generally greater than the maximum compression peak, which indicates that the tensile failure tends to be near this zone. (3) For zone B3, the stress amplitudes of these curves are commonly smaller than those of other zones. When β = 0 • and 45 • , the maximum peak stress is tensile and when β = 90 • , the maximum stress is compressive. The result is similar to that of zone B1. (4) For zone B4, the result is similar to that of zone B2, in which the maximum stress is generally tensile. In addition, it should be noted that, when β = 90 • , the tensile peak is also not obvious. (5) Generally, the stress wave waveform will not change significantly in zones B2 and B4, but the stress amplitude will. In zones B1 and B3, both the waveform and amplitude of the stress wave will change. In addition, it can be observed that the first peaks of zones B1, B3 and B4 decrease first and then increase and the first peaks of zone B2 increase first and then decrease.
(4) For zone B4, the result is similar to that of zone B2, in which the maximum stress is generally tensile. In addition, it should be noted that, when  = 90, the tensile peak is also not obvious. (5) Generally, the stress wave waveform will not change significantly in zones B2 and B4, but the stress amplitude will. In zones B1 and B3, both the waveform and amplitude of the stress wave will change. In addition, it can be observed that the first peaks of zones B1, B3 and B4 decrease first and then increase and the first peaks of zone B2 increase first and then decrease.   Figure 8 presents the microcrack distribution under different tunnel dips. As shown in Figure 8, when β = 0 • , the microcracks were mainly distributed in the right sidewall of the adjacent tunnel and almost no microcracks were formed in the vicinity of the non-adjacent tunnel. When β = 45 • , there are some microcracks distributed at the roof, right sidewall and floor of the non-adjacent tunnel and the vicinity of the adjacent tunnel, while almost no microcracks are generated between the adjacent and non-adjacent tunnels. Besides, it can be observed that the failure on the right sidewall of the non-adjacent tunnel is spalling. When β = 90 • , there are a large number of microcracks generated at the roof, right sidewall and floor of the non-adjacent tunnel. Obviously, the larger the tunnel dip is, the more severe the damage to the non-adjacent tunnel caused by blasting disturbance. In addition, it can also be found that, when β = 90 • , an obvious penetrating failure zone merges between adjacent and the non-adjacent tunnels. Based on these results, it can be inferred that, when the tunnel dip exceeds a critical value, the disaster hazard induced by blasting will be more severe, because the cascading failure tends to occur between multiple tunnels. Figure 9 shows the development process of this microcrack. As shown in Figure 9, the microcrack development process can be generally divided into four stages: no crack generation stage (t < 0.1 ms), rapid microcrack growth stage (t = 0.1-2.5 ms), slow microcrack growth stage (t = 2.5-5.6 ms) and stable microcrack stage (t < 5.6 ms). In the no crack generation stage (t < 0.1 ms), there are almost no new microcracks formations in the models. In the rapid microcrack growth stage (t = 0.1-2.5 ms), the microcracks will increase rapidly and the damage caused by this stage is also the most serious. In this stage, the microcrack curves of different tunnel dip almost coincide, because the damage of surrounding rock is mainly concentrated around the borehole and the left sidewall of the adjacent tunnel. In the slow microcrack growth stage (t > 2.5 ms), the damage begins to occur around the non-adjacent tunnel as the stress wave propagation. Therefore, the three microcrack curves gradually separate. In this stage, the growth rate of the microcrack gradually slows down. In the stable microcrack stage (t > 5.6 ms), the total number of microcracks in the surrounding rock gradually tends to be stable and the damage of the surrounding rock is basically completed. Figure 8 presents the microcrack distribution under different tunnel dips. As shown in Figure 8, when  = 0, the microcracks were mainly distributed in the right sidewall of the adjacent tunnel and almost no microcracks were formed in the vicinity of the nonadjacent tunnel. When  = 45, there are some microcracks distributed at the roof, right sidewall and floor of the non-adjacent tunnel and the vicinity of the adjacent tunnel, while almost no microcracks are generated between the adjacent and non-adjacent tunnels. Besides, it can be observed that the failure on the right sidewall of the non-adjacent tunnel is spalling. When  = 90, there are a large number of microcracks generated at the roof, right sidewall and floor of the non-adjacent tunnel. Obviously, the larger the tunnel dip is, the more severe the damage to the non-adjacent tunnel caused by blasting disturbance.

Damage Characteristics
In addition, it can also be found that, when  = 90, an obvious penetrating failure zone merges between adjacent and the non-adjacent tunnels. Based on these results, it can be inferred that, when the tunnel dip exceeds a critical value, the disaster hazard induced by blasting will be more severe, because the cascading failure tends to occur between multiple tunnels.  Figure 9 shows the development process of this microcrack. As shown in Figure 9, the microcrack development process can be generally divided into four stages: no crack generation stage (t < 0.1 ms), rapid microcrack growth stage (t = 0.1-2.5 ms), slow microcrack growth stage (t = 2.5-5.6 ms) and stable microcrack stage (t < 5.6 ms). In the no crack generation stage (t < 0.1 ms), there are almost no new microcracks formations in the models. In the rapid microcrack growth stage (t = 0.1-2.5 ms), the microcracks will increase rapidly and the damage caused by this stage is also the most serious. In this stage, the microcrack curves of different tunnel dip almost coincide, because the damage of surrounding rock is mainly concentrated around the borehole and the left sidewall of the adjacent tunnel. In the slow microcrack growth stage (t > 2.5 ms), the damage begins to occur around the non-adjacent tunnel as the stress wave propagation. Therefore, the three microcrack curves gradually separate. In this stage, the growth rate of the microcrack gradually slows down. In the stable microcrack stage (t > 5.6 ms), the total number of microcracks in the surrounding rock gradually tends to be stable and the damage of the surrounding rock is basically completed.

Evolution Characteristics of Strain Energy
The blasting disturbance usually causes the accumulation, release and dissipation of the strain energy. Some studies believe that the evolution of strain energy is a major cause of disasters and an important feature reflecting the stability of surrounding rock. [4,[43][44][45]. For example, Li et al. [4] assessed the rock burst characteristics around a tunnel based on an energy index, namely strain energy density (SED). The results show that as a large amount of strain energy is released, the strain rock burst will occur on the floor and corner. Luo and Gong [43] also assessed the established invariable feature of the ultimate internal elastic index based on the law of energy release and dissipation during rock failure. To accurately evaluate the evolution characteristics of strain energy, the strain energy density (SED) is applied in this work because it is not affected by the volume of the measuring area. According to PFC2D, the strain energy stored in linear contact and parallel bond can be obtained:

Evolution Characteristics of Strain Energy
The blasting disturbance usually causes the accumulation, release and dissipation of the strain energy. Some studies believe that the evolution of strain energy is a major cause of disasters and an important feature reflecting the stability of surrounding rock. [4,[43][44][45]. For example, Li et al. [4] assessed the rock burst characteristics around a tunnel based on an energy index, namely strain energy density (SED). The results show that as a large amount of strain energy is released, the strain rock burst will occur on the floor and corner. Luo and Gong [43] also assessed the established invariable feature of the ultimate internal elastic index based on the law of energy release and dissipation during rock failure. To accurately evaluate the evolution characteristics of strain energy, the strain energy density (SED) is applied in this work because it is not affected by the volume of the measuring area. According to PFC2D, the strain energy stored in linear contact and parallel bond can be obtained: Ik n ) (16) where E s and E s denote the strain energies of linear contact and parallel bond, respectively. F n and F s denote the normal and shear force of linear contact, respectively. F n and F s denote the normal and shear force of parallel bond, respectively. Therefore, the SED can be given by: where m is the total number of linear contacts in the measuring area. A is the area of the measuring area. Figure 10 shows the SED-time curve of typical areas under different tunnel dips. As shown in Figure 10, before the blasting wave arrives (t < 0 µs), some strain energy is accumulated in the surrounding rock and can be regarded as the initial SED. During the blasting, the SED will increase sharply, then release rapidly and finally reach an approximately stable value (final SED). In addition, it can be found that, for zones B1 and B3, the maximum SEDs tend to decrease first and then increase with the tunnel dip. However, for zone B2, the maximum SEDs increase first and then decrease with the tunnel dip, which is contrary to the cases of zones B1 and B3. For zone B4, the maximum SEDs are positively correlated with tunnel dip. Furthermore, it can be found that, after the blasting wave path, the strain energy does not recover to the initial value, indicating that the blasting disturbance causes some irreversible deformation and even some severe damage.
Mathematics 2022, 10, x FOR PEER REVIEW 13 of 20 mulated in the surrounding rock and can be regarded as the initial SED. During the blasting, the SED will increase sharply, then release rapidly and finally reach an approximately stable value (final SED). In addition, it can be found that, for zones B1 and B3, the maximum SEDs tend to decrease first and then increase with the tunnel dip. However, for zone B2, the maximum SEDs increase first and then decrease with the tunnel dip, which is contrary to the cases of zones B1 and B3. For zone B4, the maximum SEDs are positively correlated with tunnel dip. Furthermore, it can be found that, after the blasting wave path, the strain energy does not recover to the initial value, indicating that the blasting disturbance causes some irreversible deformation and even some severe damage.

Discussion
In the past decades, some scholars have realized that the stress wave wavelength is an important factor affecting and controlling the failure characteristics of surrounding rock [30,46,47]. For this reason, we will focus on the influence of stress wave wavelength on the dynamic response and damage characteristics of the non-adjacent tunnel. Gener-

Discussion
In the past decades, some scholars have realized that the stress wave wavelength is an important factor affecting and controlling the failure characteristics of surrounding rock [30,46,47]. For this reason, we will focus on the influence of stress wave wavelength on the dynamic response and damage characteristics of the non-adjacent tunnel. Generally, the different stress wave wavelengths can be obtained by adjusting the duration of the stress wave. Therefore, in this section, the blasting wave amplitude and rise time are kept the same as those in Section 2.2 and three time ratios k = t m /t r (2.5, 5 and 10) are set.
Usually, the damage behavior of the surrounding rock is controlled by its total stress. However, the total stress of the surrounding rock is not only related to the dynamic stress but also closely related to the static stress. In order to accurately evaluate the dynamic effect caused by blasting disturbance, some scholars generally define a dynamic stress amplification factor (DSAF) [40,48]. For this reason, a dynamic stress amplification factor is applied to this work: where ϕ is the dynamic stress amplification factor. σ(t) is the total stress and σ s0 is the static stress before blasting. Figure 11 presents the dynamic stress amplification factor around a non-adjacent tunnel. As shown in Figure 11, for zone B1, the DSAFs increase slightly at first and then increase sharply with the tunnel dip β. For zone B2, the DSAFs tend to increase first and then decrease with the tunnel dip. Especially, when the time ratio k = 10, the DSFA decreases slightly at first. For zones B3 and B4, the DSAFs tend to decrease initially and then increase with the tunnel dip, but for the case of k = 2.5, it can be found that the DSAFs monotonically increase with the tunnel dip. On the other hand, the DSAF generally increases with the time ratio k, which means that the longer the stress wave wavelength, the greater the total maximum stress in the surrounding rock. Figure 12 presents the microcrack distribution in surrounding rock under different wavelengths. When β = 0 • , there are almost no microcracks around the non-adjacent tunnel, but when the wavelength increases to a certain extent (e.g., k = 10), the microcracks will gradually extend from the adjacent tunnel to the non-adjacent tunnel, such as the microcrack C1. This phenomenon shows that, when the wavelength exceeds a certain critical value, the instability of adjacent tunnels may lead to instability in non-adjacent tunnels. When β = 45 • and β = 90 • , it is clear that, with the increase of time ratio k, the damage around the non-adjacent tunnel becomes more and more serious. Besides, the interaction between the two tunnels seems to be more obvious under the higher wavelength. For example, when β = 45 • , no penetrating failure zone forms between the two tunnels under the smaller time ratios (e.g., k = 2.5), but does do so under larger time ratios (e.g., k = 10). Certainly, when β = 90 • , the penetrating failure zones were also formed between the two tunnels, but it is obvious that the penetrating failure zones increase with the wavelength. These results show that the long stress wave is more likely to cause damage to the non-adjacent tunnel than the short stress wave. Therefore, the designing and protecting strategy of underground engineering can be optimized. For example, for dynamic disturbance with a short wavelength, the roof and floor of a non-adjacent tunnel should be protected. However, for dynamic disturbance with a long wavelength, the areas between the adjacent tunnel and non-adjacent tunnel should also be monitored to determine whether there is a trend of penetrating failure.
around a non-adjacent tunnel. As shown in Figure 11, for zone B1, the DSAFs increase slightly at first and then increase sharply with the tunnel dip . For zone B2, the DSAFs tend to increase first and then decrease with the tunnel dip. Especially, when the time ratio k = 10, the DSFA decreases slightly at first. For zones B3 and B4, the DSAFs tend to decrease initially and then increase with the tunnel dip, but for the case of k = 2.5, it can be found that the DSAFs monotonically increase with the tunnel dip. On the other hand, the DSAF generally increases with the time ratio k, which means that the longer the stress wave wavelength, the greater the total maximum stress in the surrounding rock.  nel, but when the wavelength increases to a certain extent (e.g., k = 10), the microcracks will gradually extend from the adjacent tunnel to the non-adjacent tunnel, such as the microcrack C1. This phenomenon shows that, when the wavelength exceeds a certain crit- ical value, the instability of adjacent tunnels may lead to instability in non-adjacent tunnels. When  = 45 and  = 90, it is clear that, with the increase of time ratio k, the damage around the non-adjacent tunnel becomes more and more serious. Besides, the interaction between the two tunnels seems to be more obvious under the higher wavelength. For example, when  = 45, no penetrating failure zone forms between the two tunnels under the smaller time ratios (e.g., k = 2.5), but does do so under larger time ratios (e.g., k = 10).
Certainly, when  = 90, the penetrating failure zones were also formed between the two tunnels, but it is obvious that the penetrating failure zones increase with the wavelength. These results show that the long stress wave is more likely to cause damage to the nonadjacent tunnel than the short stress wave. Therefore, the designing and protecting strategy of underground engineering can be optimized. For example, for dynamic disturbance with a short wavelength, the roof and floor of a non-adjacent tunnel should be protected. However, for dynamic disturbance with a long wavelength, the areas between the adjacent tunnel and non-adjacent tunnel should also be monitored to determine whether there is a trend of penetrating failure. The rapid release of strain energy is an important index of rock burst intensity. Generally, more release of strain energy will cause severe rock burst. To further evaluate the strain energy evolution law before and after blasting, the SED variation is defined, which is the difference between initial SED and final SED. The value of SED variation greater than zero represents the accumulation of strain energy and the value less than zero represents the release of strain energy. Figure 13 presents the variation of strain energy density before and after blasting. For zone B1, the SED variations tend to initially increase and then decrease with the tunnel dip . Especially, when k = 10, the SED variation directly decreases with the tunnel dip . This is the reason that the accumulated deformation of zone B1 tends to increase with the tunnel dip , but when the tunnel dip approaches a specific value, some damage will occur around the zone B1, which leads to the partial deformation recovery of the zone B1 and the decrease in SED variation. For example, when  = 45 and k = 2.5 or 5, there is almost no damage in zone B1, so the corresponding SED variation is positive; but when  = 45 and k = 10, there is obvious damage in zone B1 (as shown in Figure 12), so the corresponding SED variation is negative. The rapid release of strain energy is an important index of rock burst intensity. Generally, more release of strain energy will cause severe rock burst. To further evaluate the strain energy evolution law before and after blasting, the SED variation is defined, which is the difference between initial SED and final SED. The value of SED variation greater than zero represents the accumulation of strain energy and the value less than zero represents the release of strain energy. Figure 13 presents the variation of strain energy density before and after blasting. For zone B1, the SED variations tend to initially increase and then decrease with the tunnel dip β. Especially, when k = 10, the SED variation directly decreases with the tunnel dip β. This is the reason that the accumulated deformation of zone B1 tends to increase with the tunnel dip β, but when the tunnel dip approaches a specific value, some damage will occur around the zone B1, which leads to the partial deformation recovery of the zone B1 and the decrease in SED variation. For example, when β = 45 • and k = 2.5 or 5, there is almost no damage in zone B1, so the corresponding SED variation is positive; but when β = 45 • and k = 10, there is obvious damage in zone B1 (as shown in Figure 12), so the corresponding SED variation is negative.
For zone B2, as the tunnel dip β increases, the damage to surrounding rock will be more and more serious. Therefore, the SED variation tends to decrease with the tunnel dip β (such as in the case of k = 5 and 10). Especially, when k = 2.5, the SED variation exhibits a trend of decreasing first and then increasing. This is because, when the tunnel dip is 90 • , the total energy obtained by zone B2 is greater than those of the other dips. Another reason is that when k = 2.5, the stress wavelength is too short and there is not enough damage in zone B2. The combination of these two causes leads to the insufficient release of strain energy. Therefore, the SED variation in the case of β = 90 • and k = 2.5 is the largest.
For zones B3 and B4, the SED variations decrease with the tunnel dip β. In addition, it is worth noting that, for zones B2 and B4, the SED variation decreases with the time ratio k. The results suggest that the strain energy release of the roof and floor will increase with the wavelength. For zones B1 and B3, there is no simple positive or negative correlation between the SED variation and wavelength.
In summary, the stress evolution, energy evolution and damage characteristics of the surrounding rock are closely related to the tunnel dip and stress wave wavelength. In practical engineering, the effect of the nearby blasting disturbance can be predicted to a certain extent based on the existing wavelength information of the blasting wave and tunnel distribution information. For example, in two successive blasting activities with different charge lengths, the information of the second blasting, such as the damage and energy evolution characteristics of the roof and floor of the non-adjacent tunnel, can be effectively predicted based on the first blasting. Therefore, it is very necessary to evaluate the dynamic response and damage of non-adjacent tunnels with different tunnel distributions. It should be noted that the real rock stratum may contain a large number of random discontinuities (not considered in this study). According to previous studies [49,50], these discontinuities often affect the mechanical behavior of surrounding rock widely. Therefore, in our subsequent research, the properties of these discontinuities, including size, density, distribution and cohesiveness, will be further considered to reveal the dynamic behavior of the tunnels.
is the difference between initial SED and final SED. The value of SED variation greater than zero represents the accumulation of strain energy and the value less than zero represents the release of strain energy. Figure 13 presents the variation of strain energy density before and after blasting. For zone B1, the SED variations tend to initially increase and then decrease with the tunnel dip . Especially, when k = 10, the SED variation directly decreases with the tunnel dip . This is the reason that the accumulated deformation of zone B1 tends to increase with the tunnel dip , but when the tunnel dip approaches a specific value, some damage will occur around the zone B1, which leads to the partial deformation recovery of the zone B1 and the decrease in SED variation. For example, when  = 45 and k = 2.5 or 5, there is almost no damage in zone B1, so the corresponding SED variation is positive; but when  = 45 and k = 10, there is obvious damage in zone B1 (as shown in Figure 12), so the corresponding SED variation is negative. For zone B2, as the tunnel dip  increases, the damage to surrounding rock will be more and more serious. Therefore, the SED variation tends to decrease with the tunnel dip  (such as in the case of k = 5 and 10). Especially, when k = 2.5, the SED variation exhibits a trend of decreasing first and then increasing. This is because, when the tunnel dip is 90°, the total energy obtained by zone B2 is greater than those of the other dips. Another reason is that when k = 2.5, the stress wavelength is too short and there is not enough damage in zone B2. The combination of these two causes leads to the insufficient release of strain energy. Therefore, the SED variation in the case of  = 90 and k = 2.5 is the largest.
For zones B3 and B4, the SED variations decrease with the tunnel dip . In addition, it is worth noting that, for zones B2 and B4, the SED variation decreases with the time ratio k. The results suggest that the strain energy release of the roof and floor will increase with the wavelength. For zones B1 and B3, there is no simple positive or negative correlation between the SED variation and wavelength. In summary, the stress evolution, energy evolution and damage characteristics of the surrounding rock are closely related to the tunnel dip and stress wave wavelength. In practical engineering, the effect of the nearby blasting disturbance can be predicted to a certain extent based on the existing wavelength information of the blasting wave and tunnel distribution information. For example, in two successive blasting activities with dif-

Conclusions
This work extends the existing research on the mechanical problem of multiple tunnels under static load to the dynamic problem, because blasting disturbance often triggers a different mechanical response from that under static load. In addition, it is also different from previous studies in that the focus of this work has shifted from a single tunnel or the tunnel near the disturbance source to a remote non-adjacent tunnel. To study the dynamic behavior of the non-adjacent tunnel, the numerical models of two tunnels with different distributions were established by the particle flow code (PFC2D).
The dynamic stress evolution around the non-adjacent tunnel is initially examined. It can be found that the tunnel distribution can affect the dynamic stress response around the non-adjacent tunnel. In the roof and floor of the non-adjacent tunnel, the stress wave waveform is commonly unchanged, while the stress amplitude will change obviously. In the sidewalls of the non-adjacent tunnel, both the waveform and amplitude of the stress wave will change obviously. In addition, the damage of the surrounding rock is closely related to the tunnel distribution. Generally, the damage around the non-adjacent tunnel increases with the tunnel dips. When β = 0 • , there is no damage around the non-adjacent tunnel, but when β = 90 • , there is obvious damage around the non-adjacent tunnel and the penetrating failure forms between the two tunnels. On the other hand, the evolution of strain energy was examined. In general, the strain energy density (SED) will undergo a stage of rapid accumulation and release and the maximum strain energy density in different areas around the non-adjacent tunnel will show different trends with the increase in tunnel dip.
The dynamic response and damage characteristic of the non-adjacent tunnel caused by the blasting wave with different wavelengths was further examined. It can be found that the tunnel dip has an obvious influence on the dynamic stress amplification factor (DSAF) and the effect of tunnel dip changes with wavelength. Generally, the DSAF around the nonadjacent tunnel increases with the wavelength. By observing the distribution characteristics of microcracks around non-adjacent tunnel under different wavelength, it can be found that the long wavelength is more likely to induce rock mass damage than the short wavelength. Subsequently, the SED variation before and after blasting was further analyzed. In the right sidewall of the non-adjacent tunnel, the SED variations tend to increase first and then decrease with the tunnel dip. In the roof of the non-adjacent tunnel, the SED variations tend to decrease with the tunnel dip. In the left sidewall and floor, the SED variations tend to decrease with the tunnel dip. In addition, the longer the wavelength, the more conducive to the strain energy release of the roof and floor. Overall, this study provides some insight into the dynamic behavior of local twin tunnels. Some prediction and analysis, including stability analysis of surrounding rock, vibration assessment of surrounding rock and optimization and control of blasting method, can be preliminarily given.