Abstract
A brief literature review of numerical studies on excavation damage zone (EDZ) is conducted to compare the main numerical methods on EDZ studies. A hybrid finite-discrete element method is then proposed to model the EDZ induced by blasts. During the excavation by blasts, the rock mass around the borehole is subjected to dynamic loads, i.e., strong shock waves crushing the adjacent rocks and high-pressure gas expanding cracks. Therefore, the hybrid finite-discrete element method takes into account the transition of the rock from continuum to discontinuum through fracture and fragmentation, the detonation-induced gas expansion and flow through the fractured rock, and the dependence of the rock fracture dynamic behaviour on the loading rates. After that, the hybrid finite-discrete element method is calibrated by modelling the rock failure process in the uniaxial compression strength (UCS) test and Brazilian tensile strength (BTS) test. Finally, the hybrid finite-discrete element method is used to model the excavation process in a deep tunnel. The hybrid finite-discrete element method successfully modelled the stress propagation and the fracture initiation and propagation induced by blasts. The main components of the EDZ are obtained and show good agreements with those well documented in the literature. The influences of the initial gas pressure, in situ stress, and spacing between boreholes are discussed. It is concluded that the hybrid finite-discrete element method is a valuable numerical tool for studying the EDZ induced by blasts in deep tunnels.
1. Introduction
Blast is widely employed to develop an underground space, such as tunnels, mines, and shafts. Due to the blasting impact, the original state of rock mass around a tunnel is inevitably changed in the form of creation of new fractures, closure and opening of preexisting fractures, and redistribution of stress [1]. The disturbance zone is often referred interchangeably to as excavation damaged zone (EDZ), damaged rock zone (DRZ) [2, 3], or blast-induced damaged zone (BIDZ) [4]. As EDZ affects the overall response of the near-field rock mass and plays a negative role in the stability of a tunnel, more considerations are needed for this zone during the design stage. Additionally, it is of paramount significance to better understand the mechanical behaviour of EDZ and its influences on the performance of excavations since this will benefit the underground geostructures, e.g., tunnelling and drifting projects, on which the EDZ is highly concerned.
Generally speaking, the blast-induced EDZ is formed during the blasting interaction of an individual borehole. Figure 1 shows the fracture process by the blast of a single borehole in an infinite rock mass (Figures 1(a)–1(d)) and a rock mass with one free surface (Figure 1(e)). It can be seen from Figure 1 that the fracture pattern consists of three zones, i.e., crushed zone, cracked zone, and fragment formation zone. If free surfaces are present, failure at the free surface due to the reflected tensile wave will occur.
(a)
(b)
(c)
(d)
(e)
While the explosive is detonated, a great heat and high-density gas coupled with an extremely high-pressure pulse, i.e., explosive wave, is generated. Then, dilatational waves at sonic velocities in the rock are produced as the high-pressure pulse transmits in the rock [6]. A crushed zone which lies in the vicinity around a borehole is then formed by the high pressure, which exceeds severalfold the compressive strength of the rock, exerting on immediately around the borehole (Figure 1(a)). A crushed zone is characterized by intensive crushed and shattered fragments. The compressive failures or shear failures dominate the crushed zone but there are also rare tensile failures. While the wave attenuates, a cracked zone immediately follows the crushed zone in which the compressive failures decrease while the tensile failures increase. In the cracked zone, the fracture types range from severe crushing through plastic deformation to partial fracturing and the rapidly increasing particle size as the radial distance increases [5]. The fragment formation zone is mainly formed by the tensile failures as the compressive stress rarely meets the critical stress of the rock [7]. Nevertheless, the tensile stress is still large enough to cause radial fractures due to the much smaller tensile strength of rock compared with the compressive strength. Meanwhile, apart from the explosive wave, a high gas pressure, generated simultaneously while the chemical charge is detonated, also promotes and accelerates the fracture propagation, especially at the fragment formation zone. In the case of a free surface presented as shown in Figure 1(e), after forming the crushed zone, as the compressive stresses continue to propagate outwards and reach the free surface, the compressive waves are reflected as tensile waves. Then, tensile failures occur at the vicinities of the free surface. Joints and preexisting cracks produced during the blasting also play the role of the free surface. During excavating by blasts, the fracture process of the rock around every single borehole follows the steps depicted in Figure 1. However, the final results of the blast-induced EDZ will depend on coalescent and interaction of the fractures induced by the blasts of all boreholes.
In order to better understand the EDZ, many researchers have focused their study on the characterizations and classification of the EDZ [3, 4, 8]. On the basis of knowledge provided by Sato et al. [1], Saiang and Nordlund [4], Kuzyk and Martino [8], and Martino and Chandler [9], a conceptual model of EDZ in a tunnel due to excavation by blasts is illustrated in Figure 2. The EDZ, lying on the immediate vicinity of the tunnel wall (the blue area in Figure 2), consists of an inner damage zone and an outer damage zone. The mechanical and hydraulic properties are dramatically changed in the inner damaged zone while these properties are changed gradually in the outer damaged zone [3]. Outside the damaged zone is a disturbed zone in which there is no damage but the stresses are redistributed. Outside the disturbed zone is an undisturbed zone where the rock is in its in situ state. Additionally, Kuzyk and Martino [8] defined a failed zone which is part of the inner damage zone but the rocks in the failed zone have physically detached from the opening wall.
Nowadays, instead of studying EDZ experimentally, many researchers turn to study the blast-induced damaged zone by various numerical methods since the current computer technologies are able to complete large-scale numerical calculation in a short time. Based on the assumption of material continuity, numerical methods are classified into three types, i.e., continuous method, discontinuous method, and hybrid method. The continuous method, e.g., finite element method and finite difference method, mainly focuses on modelling the stress distribution of EDZ after excavation, blast-induced damage before rock fragmentation, and rock mass displacement. Sato et al. [1] used the finite element method (FEM) to analyse the rock mass displacement from a drift wall and compare it with the displacement measurement. They found that EDZ formed by mechanical excavation was narrower than that by blasting. Sheng et al. [10] used the 2D/3D elastoplastic finite element to analyse the stress redistribution and predicted the plastic yield zones in the three-gorge dam permanent ship lock slopes during excavation. It is concluded that the predicated plastic zones generally coincide with the EDZ obtained from field testing and monitoring. Liu et al. [11] developed the rock and tool interaction code (R-T2D) on the basis of the rock failure process analysis (RFPA) model to analyse the crack propagation by blast in the TASQ Tunnel of Äspö Hard Rock Laboratory. In their research, the main phenomena of rock by blasts, e.g., crushed zone and long radial cracks, are modelled. Wang et al. [12] used RFPA to model the effects of heterogeneity and anisotropy on EDZ. Four models with different excavation shapes are modelled, and it is indicated that the anisotropy and heterogeneity of a rock mass play a significant role in its engineering stability. Saiang and Norlund [4] employed the continuous method FLAC to evaluate the effects of the strength and stiffness of the blast-induced damaged zone and other relevant parameters on induced boundary stresses and ground deformation in the near-field rock mass in a shallow tunnel. They concluded that the presence of the blast-induced damage zone does affect the behaviour of the boundary stresses and ground deformation.
As mentioned above, the continuum method can capture the characters, e.g., stress distribution and rock mass displacements, of EDZ before damage occurs in an excavation process or the influence of the surroundings after excavation. However, the continuous method is limited to only sufficiently large volumes compared to the scale of the problems and individual rock discontinuity is not taken into account [13]. In the case of underground excavation, rock masses comprise joints, faults, and in situ stress-induced preexisting cracks. Moreover, the excavated rock masses are broken into fragments while cracks are produced in newly formed tunnels. Therefore, the discontinuity of the rock masses should not be ignored for the blasting simulation. Since the discontinuity of the rock mass is not taken into account, the continuous method has limited us to better understanding the EDZ formation process in terms of fracture initiation and propagation, especially for the blast-induced EDZ. To overcome the limitation, the discontinuous method is then used to model EDZ. Jiang et al. [14] developed an expanded distinct element method (EDEM) based on the distinct element method code UDEC to carry out excavation simulations of a deep underground cavern. The excavation process is modelled, and the positions of crack initiation are investigated. The functions of rock bolts on controlling the deformation of rock masses and the generation of new cracks are studied by the proposed discontinuous method. The KTH-SKI research team [15] used Particle Flow Code (PFC) based on the discrete element method (DEM) to assess the time-dependent effects on the EDZ around a heat-releasing nuclear waste emplacement drift in fractured rock. The research indicated that the time-dependent effects played a very limited impact on mechanical behaviour for this particular simulation problems.
The discontinuous method provides unique advantages for modelling the behaviour of the rock mass which is generally considered losing continuity during its fracture. However, the large computational demand tends to limit its applicability to small-scale problems. Under this scenario, hybrid methods, which take the advantages of both continuous and discontinuous method, are adopted to study the rock mass stability influenced by EDZ. Cai et al. [16] used the coupled numerical method, FLAC/PFC, to capture acoustic emission (AE) and microseismic (MS) events for locating and evaluating the underground excavation induced rock mass degradation or damage. In FLAC/PFC, the rock surrounding an AE sensor is modelled using PFC while the remaining rock mass is simulated by FLAC. Thus, the calculation time can be significantly reduced. Saiang [17] used a coupled continous-discontinuous method FLAC and PFC2D to analyse the behaviour of the blast-induced damage zone. It is found that the failures mapped in the hybrid method, i.e., FLAC/PFC2D, for shallow excavations were consistent with the yield elements in Phase 2. It is worth noting that, in the FLAC/PFC, the continuous and discontinuous materials are confined in specified zones, respectively. The transition from continuum to discontinuum is not fully achieved, which limited the method to model the large-scale fracture and fragmentation induced by blasts.
To naturally model the crack initiation and propagation and resultant fragmentation, the continuum to discontinuum should be freely transited. More recently, a combined finite-discrete method, originally proposed by Munjiza [18], was used to study the fracture and fragmentation of rock mass induced by impact loads [13], the fracturing behaviour of unsupported circular excavations in laminated rock masses under various in situ stress conditions [19], and damage process and failure mechanisms around underground openings in clay shales [20]. In addition, the authors also used a hybrid finite-discrete element method to study the blast-induced excavation damaged zone in the top heading of deep tunnels [21]. Recently, a grain-based finite-discrete element method is proposed to investigate the crack initiation and propagation mechanism in brittle rocks [22]. The hybrid finite-discrete element method can naturally model the transition from continuum to discontinuum through failure, fracture, and fragmentation of the continua [13, 18], which makes it a promising method for EDZ study.
According to the brief literature review, rare numerical methods are able to model the entire formation process of the blast-induced EDZ, i.e., stress propagation, crack initiation and propagation, and separation, between the newly formed tunnel wall and excavated rock masses. Therefore, in this study, a hybrid finite-discrete element method is employed to model the EDZ formation process in a deep tunnel to gain insight into the rock crack initiation and propagation during excavation by blasts.
2. Hybrid Finite-Discrete Element Method
A hybrid finite-discrete element model to be used for modelling EDZ formation process by blasts can include a single discrete element body or a number of interactive discrete bodies, and each of which is of general shape and size and can be discretized into finite elements. The interactions of discrete bodies are model by the discrete element method while the deformability, fracture, and fragmentation of discrete bodies are modelled by the finite element method. The essential components of the hybrid finite-discrete element method consist of contact detection and interaction among those individual discrete bodies, deformability, and transition from continuum to discontinuum through fracture and fragmentation of individual discrete bodies, temporal integration scheme, and computation fluid dynamics. Among them, the transition from continuum to discontinuum through fracture and fragmentation is the key component, which is used to model the deformability, fracture, and fragmentation of rock by blasts. The transition from continuum to discontinuum in the hybrid finite-discrete element method occurs through three fracture modes, i.e., mode-I (i.e., tensile failure), mode-II (i.e., shear failure), and mixed-mode I-II [18, 23, 24]. However, it only occurs in the mode-I fracture mode in the original open-source combined finite-discrete element library Y2D. Thus, mode-II and mixed-mode are introduced in detail in Section 2.1. Moreover, the computation fluid dynamics is another key component, which is used to model the gas pressure as surface loads on the rock around the borehole and its flow into subsequently generated fractures. Thus, Section 2.2 briefly introduces computation fluid dynamics implemented into hybrid finite-discrete element method. Furthermore, for the rock fracture and fragmentation process during rock blasting, the rock is under impact load. Since the behaviour of the rock is sensitive to loading rate, the empirical relationship between physical-mechanical properties of rock and the loading rate is taken into account in the hybrid finite-discrete element method and it is briefly introduced in Section 2.3.
2.1. Transition from Continuum to Discontinuum
The hybrid finite-discrete element numerical model for modelling EDZ formation process during blasting in a deep tunnel consists of a number of interactive discrete bodies with general shape and size. The interaction of the discrete elements is governed by contact law: a stiffness in the normal direction and a stiffness and friction angle in the tangential direction. The discrete elements are then discretized into finite elements to analyse the transition from continuum to discontinuum of discrete bodies, i.e., deformability, fracture, and fragmentation.
The finite elements are bonded together in the discrete bodies before fracturing and the fractures are assumed to coincide with the element surfaces during fracturing. Thus, the transition from continuum to discontinuum occurs through the separation of the bonded finite elements. A bonding stress is involved during separating the elements, which is regarded as the function of the size of separation . At any point on the surface of a crack, the separation of can be divided into two components in the following equation:where and are the unit vectors in the normal and tangential directions, respectively, of the surface at such a point and and are the magnitudes of the components of in the normal and tangential directions, respectively.
Figure 3(a) illustrates the relationship of bonding stress and opening displacement of element surface for the mode-I fracture (i.e., tensile fracture) while Figure 3(b) shows the relationship of bonding stress and sliding displacement of element surface for mode-II fracture (i.e., shear fracture).
(a)
(b)
In terms of mode-I fracture, i.e., tensile fracture, as shown in Figure 3(a), if the opening displacement of element surfaces, i.e., the separation in the normal direction , reaches the critical displacement, i.e., , prescribed by the tensile strength of the element, , the tensile damage is assumed to occur. As the opening displacement of element surfaces increases, i.e., , the bonding stress gradually decreases. When the opening displacement of element surfaces is beyond the ultimate opening displacement , tensile failure has completed.
The relationship of bonding stress and opening displacement of the element surface for tensile failure is indicated by the following equation:where is a damage variable between 0 and 1, is damage function described in the mechanical damage model [25], is the tensile strength, and all other parameters have the same meanings as those introduced before.
As for the mode-II fracture, i.e., shear fracture, as shown in Figure 3(b), if the sliding displacement of the adjacent element surface reaches the critical value, i.e., , defined by the shear strength of the element, , the damage is assumed to occur. As the sliding displacement of adjacent element surface increases, i.e., , the bonding stress gradually decreases. As the sliding displacement of element surfaces exceeds a residual opening displacement , the shear failure has completed and the bonding stress becomes a purely frictional resistance. Equation (3) is the relationship of sliding displacement of adjacent element surfaces and the bonding stress for shear failure:where is a damage variable between 0 and 1, is damage function described in the mechanical damage model [25], is the joint residual friction angle, is the shear strength, and all other parameters have the same meanings as those introduced before.
The ultimate opening displacement for tensile failure is determined by the mode-I release rate which is equal to the area under the curve of bonding stress and opening displacement depicted by equation (4) while the residual sliding displacement depends on the mode-II fracture energy release rate which is described by the following equation:
In the case of mixed-mode I-II fracture, if a criterion (i.e., equation (6)) involving both the opening and sliding parameter is satisfied, the fracture occurs:
The mode-I and mode-II are pure tensile failure and shear failure, respectively. The mixed-mode I-II describes the rock failure under combination of the tensile stress and shear stress. There is rarely pure compressive failure. The compressive stress at any point of the surface of a crack or element boundary can be divided into stress in normal direction and tangential direction. Thus, the compressive stress in tangential direction will result in the sliding of the adjacent elements along the element boundaries, which eventually results in shear failure. The compressive stress in the normal direction can lead to the element penetration, i.e., adjacent element will penetrate to each other in the normal direction along the element boundaries. The compressive stress in the normal direction can be obtained by the penalty function method, which is used to describe the contact force of two bodies penetrate each other.
After fracture and fragmentation, the central difference explicit time integration scheme is applied in the hybrid finite-discrete method to integrate the equations of motion of either the initially discrete element or the discrete elements formed by the fracture and fragmentation algorithm.
The central difference explicit time integration scheme assumes thatwhere , , and are the displacement, velocity, and acceleration, respectively; is the sum of the body forces, contact forces, and external loads together with any damping forces, and is the mass associated; is the current time while is the time increment.
As an explicit scheme, there is no need for stiffness matrices to be assembled or stored in the step-by-step solution. The central difference explicit integration scheme obtains the solutions for forces, stresses, and displacements at time using the equilibrium condition at the last time steps, i.e., at time . It should be noted that the time step adopted for the central difference explicit integration scheme must be smaller than a critical value for a stable scheme. The approximate critical value for the time step can be obtained by the following equation:where is the smallest elemental length, and are Lame constants, and is the density of geomaterials.
2.2. Detonation-Induced Gas Expansion and Flow through the Fracturing Rock Mass
For modelling the EDZ formation process during blasting, it is required to model the detonation-induced gas expansion and flow through the fracturing rock mass process. As the explosive is detonated, the borehole is full of gaseous detonation productions at very high pressure and temperature. Then, the pressure is exerted immediately on the wall of the borehole and generating radial compressive stress [26]. Since the pressure exerted on the borehole wall is much higher than the rock strength, the rock becomes yielded and extensively broken or crushed. The borehole pressure generated by the detonated explosive describes its expansion work during rock breakage process, and this pressure indicates the transfer of the explosive energy into the rock [26]. This is first section of rock blasting process, i.e., detonation-induced gas produces exerting high pressure on the borehole wall and generating radial compressive stress. In the first section, the borehole pressure is the most important information in successfully modelling the explosive performance and predicting the blasting results. Instead of modelling the chemical reaction in the borehole, many researchers used the pressure-time histories generated from empirical equations to model the detonation-induced pressure, i.e., borehole pressure.
As this research is not aimed at study, the detonation theories but the rock fracture and fragmentation by blasts, the entire chemical reaction process of explosive in a borehole is not simulated. However, a gas pressure-time history curve generated from commercial finite element code AUTODYN modelling the chemical reaction in the blast through the Jones-Wilkens-Lee (JWL) equation of state (EOS) in equation (11) is used to model the interaction between detonation production and surrounding rocks:where is the instantaneous pressure at any time, , , , , and are the material constants, and are the densities of the explosive and the detonation products, respectively, and is the specific energy.
The second section of the rock blasting process is the detonation gas with high-pressure flowing through the rock cracks. During the second section, the detonation exerts pressures on the cracks. Thus, the fractures propagate to form the crushed zone and cracks out of the crushed zone. To model this process, a relatively simple model for gas flow, the voids created by the fracturing and fragmenting solid, is adopted [18]. For this model, the gas is present only within a predefined zone, i.e., a circle around the borehole [18, 23]. The cracks in the gas zone are approximated by constant area ducts while the detonated gas is approximated by steady-state flow of ideal gas. Thus, the spatial and temporal distribution of gas pressure exerting on cracks can be achieved by calculating the gas pressure in the constant area ducts.
2.3. Effect of the Loading Rate on the Dynamic Behaviour of Rock
Studies have shown that the strengths of the rock increase with the increasing loading rate (e.g., [27, 28]). Thus, the effect of loading rate significantly influences the fracture behaviour of rock. In the excavation by blasts, the rocks surrounding boreholes are subjected to strong shock waves and high gas pressure. The loading rate is much higher than that for static loads. Therefore, the effect of the loading rate on the dynamic deformation and fracture behaviour of rock should be taken into account due to the dynamic loads, i.e., high gas pressures ranging from 0.5 GP to 50 GP [29] and strong shock waves with a velocity around 5500 m/s.
In order to take the loading rate into account, a semilog formula in equation (12) proposed by Zhao [28] is implemented to the hybrid finite-discrete element method to describe the relationship between the strength and the loading rate:where is the dynamic uniaxial compressive strength (MPa), is the dynamic loading rate (MPa/s), is the quasistatic loading rate (approximately ), is the uniaxial compressive strength at the quasistatic loading rate (MPa), and is a material parameter, which is 11.9 for the Bukit Timah granite [28].
As mentioned in Section 2.1, the fracture energy release rate is implemented to regulate the mode-I and mode-II damage processes. In the hybrid finite-discrete method, the fracture energy release rate is assumed to increase with the loading rate according to equation (12). Thus, the implementation of the fracture energy release rate in the hybrid method can model the tensile and shear failure under dynamic loading rates.
3. Calibration of the Hybrid Finite-Discrete Element Method
The calibration of the hybrid finite-discrete element method in modelling of the rock fracture and fragmentation processes has been done in detail in the authors’ previous research [23]. In our previous research, a single-borehole blast in a rock mass with a free surface, two-borehole simultaneous blasts, and two-borehole consecutive blasts have been modelled and the obtained results have been compared with those documented in the literature. The hybrid finite-discrete element method has well modelled the stress wave propagation, reflection and attenuation, the detonation gas flow through the fracturing rock, the formation of the crushed, cracked and elastic zones, the long radial crack initiation and propagation, the tensile failure caused by the stress wave reflection at the free surface, and the effect of various properties of the explosive and rock on rock fracture process in rock blast [23]. Thus, instead of modelling rock blasting, to calibrate the hybrid finite-discrete element method in this research, the rock fracture processes in uniaxial compression strength (UCS) test and Brazilian tensile (BTS) test are modelled using the proposed method, and the obtained results are compared with those conducted in the rock mechanics laboratory and documented in the literature.
A series of rock fracture tests under static loads is taken to obtain the rock parameters and rock fracture patterns. Figure 4 illustrates the rock sample for the Brazilian tensile strength (BTS) tests and uniaxial compressive strength (UCS) test while Figure 5 demonstrates the rock fracture patterns of the tests which will be used to compare with the model results. Table 1 shows the rock properties obtained from the static test which is then used as input parameters for the hybrid finite-discrete element method for modelling the UCS and BTS tests.
(a)
(b)
(c)
(d)
(a)
(b)
3.1. Modelling Rock Fracture during the Uniaxial Compressive Strength Test
Figure 6 depicts the geometrical and numerical models constructed for the UCS test following the standards of the ISRM suggested the method for rock characterization, testing, and monitoring: 2007–2014 [30]. For the specimen, as shown in Figure 6, the diameter of the top and bottom of the sample is 50 mm, and the height is 100 mm. It can be seen from Figure 6(a) that the test is simplified as plane stress problems and only the vertical section is considered.
(a)
(b)
The rock specimen is placed between two plates. During the UCS test, a constant displacement increment of 0.1 m/s is applied on the top-loading plate for modelling the rock failure process while both top- and bottom-loading plates are fixed in the horizontal direction. As shown in Figure 6(b), the numerical model is discretized by triangular elements and it can be seen that a relatively big mesh size is adopted. The reason why a relative big mesh size is applied is that the calculation time can be significantly reduced as the mesh size is significantly related to the calculation time. The time step for the UCS test is .
The material properties of the rock specimen follow those in Table 1. The mode-I and mode-II fracture energy release rates are essential parameters for hybrid models, which are not directly obtained by tests in this research. However, it can be calculated by the following equations:where and are the mode-I and mode-II fracture energy release rates and and are the mode-I and mode-II fracture toughness.
For the mode-I fracture toughness , a relationship between tensile strength and mode-I fracture toughness for various rock materials as shown in equation (15) is adopted to estimate it [31]:
In addition, Ingraffea, Whittaker et al., and Al-Shayea et al. indicated that model-II fracture toughness is usually larger than mode-I fracture toughness [32–34]. According to Brazilian disc tests for limestone rock, the ratio of is 1.72 [35], while it is 1.64 obtained from numerical modelling of the notched Brazilian disc test for granite materials [36]. Accordingly, the mode-II fracture toughness is assumed to be 1.7 times larger than the mode-I fracture toughness. Thus, the mode-I and mode-II fracture energy release rates can be estimated on the basis of the tensile strength and equations (13)–(15). The penalty terms of the rock specimen and loading plates are assumed to be equal to the elastic modulus of the rock specimen and half of that of the steel, respectively.
Figure 7(a) illustrates the minor principal stress evolution process during the UCS test, and the value of the minor principal stress is indicated by the colour legend on the top side of Figure 7(a). Figure 7(b) shows the crack initiation and propagation with the increasing displacement of the top plate, and the loading times are also indicated. In the UCS test, the red colour represents shear failure while the blue represents tensile failure and model boundaries. Figure 8 shows the stress-displacement curve of the UCS test, and the alphabets on the curve correspond to those in Figure 7.
(a)
(b)
As the top plate contacts the rock specimen, the load is applied on the top-loading plate. Then, the load is transferred to the upper part of the rock specimen through the contact between the top-loading plate and the rock specimen (Figure 7(a)-A). The arrow in Figure 7(a)-A indicates the stress propagation direction. The stress induced by the applied load from the top-loading plate propagates from the upper part of the rock specimen to the lower part of the rock specimen (Figure 7(a)-B). As the stress researches the contact of the bottom of the rock specimen and the bottom-loading plate, a load from the bottom-loading plate is applied on the lower part of the rock specimen due to the interaction of the rock specimen and the bottom-loading plate. The arrow in Figure 7(a)-C indicates the reflection of the stress from the bottom plate. The stress is then reflected from the contact of the rock specimen and the bottom-loading plate. After the stress propagates between the top and bottom for a few times, relatively uniform stress distribution in the specimen is formed (Figure 7(a)-E). Before the stress meets the relative stress equilibrium, there is no fracture (Figure 7(b)-A–D) and the stress-displacement curve is approximately linear (Figure 8-A–E). At point E in Figure 8, the stress almost meets its peak stress and a shear crack occurred at the top-right corner of the rock specimen (Figure 7(b)-E). With the top-loading plate continues to move, a shear fracture at the middle of the specimen is produced (Figure 7(b)-F) and the stress-loading displacement curve approximately meets its peak (Figure 8-F). With the increase of the displacement of the top-loading plate, some shear fractures with an angle of approximate 45° occur and some vertical tensile failures due to the expansion are produced (Figure 7(b)-G–K). Meanwhile, the top load drops dramatically from their peaks to the bottom as shown in Figure 8-F–L, which indicates that the specimen loses its ability to carry loads.
Figure 9 compares the finally fracture pattern of the modelled and the experimental results. For the modelled results, the blue colour represents the rock mass while the white colour represents the fractures. In general, the modelled result agrees with the experimental result in terms of the fracture pattern. In both the two results, some fragments are produced on the upper part of the specimen due to the high-stress concentration at the top part of the specimen, and some fractures are produced with a certain angle. According to the stress-displacement curve, hybrid finite-discrete element method can reproduce the failure process of the UCS test: a nearly linear loading increase region (Figure 8-A–F), a big jump from the peak value to bottom (Figure 8-F–H), and a postfailure region (Figure 8-H–L). In addition, according to the stress-displacement curve (Figure 8), the compressive strength, i.e., the peak stress, is 107 MPa. The ratio between the modelled result (107 MPa) and the input value (103 MPa) is 1.04, which indicates that the hybrid finite-discrete element method can well model the rock behaviour of the rock. However, according to the peak stress (107 MPa) and the corresponding strain (0.25 mm/100 mm), the Young`s modulus is approximately 42.8 GPa, which is much lower than the input value. It might be caused by the relatively big mesh size adopted in the UCS test model. The strain of the specimen is mainly caused by the opening or sliding of the adjacent finite element and the penetration of the elements to each other. The big mesh size indicates the big finite elements adopted. As compared with the geometrical model, the finite elements are much larger. Thus, the deformation of the model might be influenced. According to the modelled result, the strain of the model is relatively larger. Thus, the calculated elastic modulus is smaller than the input value.
(a)
(b)
3.2. Modelling Rock Fracture during the Brazilian Tensile Strength Test
Figure 10 depicts the geometrical and numerical models constructed for the Brazilian tensile strength test following the standards of the ISRM suggested the method for rock characterization, testing, and monitoring: 2007–2014 [30]. As shown in Figure 10, the test is simplified as a plane stress problems and only the vertical section is considered. The time step for BTS test modelling is . During the BTS test, a constant displacement increment of 0.1 m/s is applied on both the top- and bottom-loading plates for modelling the rock failure process while both top- and bottom-loading plates are fixed in the horizontal direction.
(a)
(b)
Figure 11 visually illustrates the modelled failure processes in terms of the distribution of the minor principal stress and the initiated cracks during the BTS test while the corresponding force and displacement curves are shown in Figure 12 in which the alphabets correspond to those in Figure 11. The magnitude of the stress in Figure 11(a) can be referred to as the legend on the top side of Figure 11(a). In the BTS test, tensile failure is marked using the red colour while the compressive cracks and the boundaries are marked using the blue colour as shown in Figure 11(b).
As the top- and bottom-loading plates contact the rock specimen, high-stress concentrations are immediately induced at the loading areas (Figure 11(a)-A). As the loading plates further move toward each other, the stress concentrations at the loading areas become more intensive (Figure 11(a)-B, C). During this period, there are no fractures produced (Figure 11(b)-A–D) and the force-loading displacement curve is approximately linear (Figure 12-A–C). As the plates continue to move, the loading force researches its peak (Figure 12-D) and a tensile crack along the vertical diameter is produced (Figure 11(b)-E, F). After that, the forces drop quickly as the sample loses its ability to carry load (Figure 12-E-F). With the displacements of the two loading plates increasing, the crack propagates towards the top- and bottom-loading points (Figure 12-G–I). In addition, some shear fractures also occur at the contact of the loading palates and rock specimen during the high-stress concentration at the loading areas (Figure 11(b)-J-K). Finally, the rock specimen is split into two halves along the vertical diameter (Figure 11(b)-K) and the force drop to a very low level (Figure 12-F–K).
According to the loading-displacement curve (Figure 12), the tensile strength of the rock specimen can be calculated as follows:
The ratio of the modelled result (12.78 MPa) to the input value (10.86 MPa) is 1.18. The modelled result is a little bit larger than the experimental result.
Figure 13 compares the hybrid finite-discrete element modelled result and experiential result for the BTS test. For the modelled result (Figure 13(a)), the blue colour represents rock material while the white colour represents fractures. For the laboratory test, there is a crack along the loading diameter due to the tensile stress and some fragments at the loading area due to the high compressive stress. The modelled fracture pattern (Figure 13(a)) includes both the tensile failure along the vertical diameter and some shear failure induced by the high compressive strength at the contact area of the disc and the loading plates. It should be noted that the main crack is relatively straight along the loading diameter, which is mainly caused by the tensile stress distributed along the loading diameter as illustrated in Figure 11(a). However, the mesh orientation is also the main reason. The rock only can crack along the joint elements in the hybrid finite-discrete element method. Since a relatively coarse mesh is adopted to save the calculation time, the joint elements (or the edges of the finite elements) can be connected along the vertical direction. Thus, the crack can propagate along the loading straight diameter.
(a)
(b)
The hybrid finite-discrete element method has well reproduced the rock failure process of the BTS, test and the modelled result agrees well with the laboratory test. In addition, the hybrid finite element method can capture the rock mechanical behaviour during the BTS test since the obtained force-loading displacement curve indicates a typical brittle material fail process.
4. Modelling the Excavation Process and Resultant EDZ Formation Process by Blasts
In order to reduce the extension of damage due to blasting, perimeter blasting or smooth blasting is widely employed for excavation, e.g., Saiang and Nordlund [4]. Therefore, in this section, the smooth blasting technique is adopted for the excavation of a deep tunnel.
In this section, the numerical model for modelling the EDZ formation process is simplified from a deep tunnel called the TASQ tunnel, which was excavated by drill and blast and specially designed for a rock mechanics experiments, the Aspo Pillar Stability Experiment (APSE) in Sweden [37].
As shown in Figure 14, the tunnel is closed to the hoist shaft and ventilation shafts and located at the −450 m level. There are three different phases for the excavation work as shown in Figure 14(b). The first phase is an ordinary tunnel excavated by full-face blasting [11]. The third phases include a top heading and a bench. The second phase is the transition phase from the first to third phases. A vertical section of the TASQ tunnel is separated into top heading and bench, and the dimensions are illustrated in Figure 15. To mitigate the influence of the boundary, a 22 m × 22 m square is taken as the boundary. Figure 16 illustrates the drilling pattern for the tunnel excavation. It is separated into two parts, i.e., top heading and bench. The spacing between boreholes is 0.4 m while the burden for the perimeter holes is 0.5 m which is somewhat tighter pattern than normal. This research aimed at illustrating the potential use of the hybrid finite-discrete element method in modelling excavation of tunnels. Thus, to save the calculated time, only the last step of the blasting is considered as illustrated in Figure 16. The red boxes in Figure 16 indicate the last steps of the excavation for the top heading and the bench. Thus, the hybrid finite-discrete element method is used to model the last step of the excavation to obtain the EDZ formation process according to the drilling pattern.
4.1. Modelling Excavation Process and Resultant EDZ Formation Process by Blasts in the Top Heading of the Tunnel
Due to the high denotation velocity of the explosive along the length of the borehole, it is supposed that the explosive is detonated at the same time along the length of the borehole, and the rock fracture mechanisms along the length of the borehole are similar in each section along the borehole. Thus, for the excavation by blasts as illustrated in Figure 16, it is can be simplified as a two-dimensional plane strain problem. As a matter of fact, the two-dimensional model has been widely employed to study the blasting mechanism under both experimental condition and numerical environment [5, 7].
Figure 17(a) demonstrates the geometrical model for the top heading of the tunnel. As can be seen, there are 46 boreholes with diameters of 0.1 m set around the top heading of the tunnel. The blasting boundary is 0.5 m while the distance between the boreholes is 0.4 m. Three points, i.e., points A, B, and C, which are located at 0.2 m, 1.5 m, and 3.5 m far from the borehole in the horizontal direction, respectively, are selected to analyse the evolution of the minor principal stress, at inner damage zone, outer damage zone, and disturbed zone, respectively. As shown in Figure 17(b), triangular elements are used to discretize the model. The dense elements are used in the interested areas, i.e., the potential EDZ area around the tunnel, to model the crack initiation and propagation, while the coarse element is employed far from the tunnel to model the stress propagation. The model is comprised of 3008 finite elements and 1564 nodes. The number of elements between two adjacent boreholes is approximately 60. The time step size is set to be to make sure the central difference explicit time integration scheme used in the hybrid finite-discrete element method is stable. It should be noted that before the last step of the excavation, the in situ stress around the tunnel wall and the blasting of boreholes beside the lasting step around the tunnel wall have already produced some fractures around the desirable tunnel wall. As for this research, the model is simplified as the rock mass has been excavated and only the rock mass around the tunnel wall is going to be excavated by the last step of the blasting, i.e., smooth blasting. Only the in situ stress-induced fractures on the rock mass of the last step excavation are considered. The rock mass in Sweden is generally classed as good quality or above average in terms of engineering classification [17]. In addition, the rock mass is dominated by plutonic Äspö diorite and Ävrö granite. Thus, the general granite rock mass properties as shown in Table 2 are adopted for the tunnel models. The in situ stresses around the model boundaries in both horizontal and vertical directions follow those obtained by Staub et al. [37] and can be found in Table 3. In the following modelling, the bottom boundary is fixed in both horizontal and vertical directions as illustrated in Figure 17(b). The in situ stresses in both the vertical and horizontal directions are applied on the other three boundaries (Figure 17(b)) of the models first to obtain deep tunnel stress condition and to consider the fractures induced by the in situ stresses in both the horizontal and vertical conditions. After the in situ stresses achieve equilibrium, the initial gas pressures are then applied on the borehole walls to simulate the stress propagation and crack initiation and propagation. In this research, in order to save the occupational time, only the last step of the exaction by blasts is considered.
(a)
(b)
4.1.1. Stress Propagation and Fracture Initiation and Propagation by Blast in the Top Heading
Figure 18 visually illustrates the temporal distributions of the in situ stress and the blast-induced minor principal stress modelled by hybrid finite-discrete element method during the smooth blasting in the top heading of the tunnel. The magnitude of the stress in Figure 18 is represented by the colour and can be referred to legend shown on the top of Figure 18. In order to achieve the stress equilibrium environment around the tunnel, the in situ stresses are applied on the three boundaries first, i.e., major principal stress lies in the horizontal direction and minor principal stress lies in the vertical direction. Figure 19 visually demonstrates the fracture initiation and propagation and rock fragmentation process during the smooth blasting.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
As shown in Figures 18(a) and 18(b), the stresses initiate from the three boundaries, i.e., top, left, and right boundary, and propagate towards the tunnel. At 0.04 ms, the stresses reach the tunnel. Then, the stresses reflect from the free surface, i.e., the tunnel wall. After the interaction of stresses from the three boundaries and those reflected from the free surface for a while, the stress distribution (Figure 18(c)) keep almost static for a while, which indicates that the stress equilibrium achieves. A few cracks are observed (Figure 19(b)) around the tunnel wall due to the in situ stress in the process of achieving equilibrium.
After that, the explosives in the 46 boreholes are detonated simultaneously and high-pressure gas in each borehole is created and interacted with the borehole wall. Then, the intense stress wave is produced immediately and propagates radially out of the borehole as illustrated in Figure 18(d). Meanwhile, some fractures are produced as shown in Figure 19(c) due to the interaction of the high-pressure gas and the borehole wall. However, no obvious crush zones around boreholes are created due to a relatively lower initial gas pressure adopted for the smooth blasting. Instead, cracks zones are produced immediately around borehole walls and expand radially by the high gas pressure. The cracks from adjacent boreholes reach each other at 0.16 ms (Figure 19(d)). Consequently, a cracked line formed which connects the boreholes (Figure 19(d)). As the gas enters the crack line, the gas pressure then mainly expands the cracked line paralleling to the free surface (Figures 18(d) and 19(d)) for the excavation of the tunnel. As the stress wave continually propagates radially and the high-pressure gas expands the cracks (Figures 18(e) and 18(f)), long cracks radially from the borehole propagate further. Meanwhile, long cracks can be observed (Figures 19(e) and 19(f)). While the stress wave reaches the tunnel wall, the compressive stress can be reflected and turn to be tensile stress. Since the tensile strength of the rock is much smaller than the compressive strength, the rock is much easier to break by tensile stress compared with compressive strength. Thus, it can be seen that due to the reflected tensile stress from the tunnel wall, the rock mass near the tunnel wall breaks into fragments (Figures 19(f) and 19(g)). Finally, the rock mass is cut along the borehole line to form a new tunnel wall and the excavated rock mass breaks into fragments. As the gas releases to air from the cracks and the strong shock wave reduces to a threshold of the rock strengths, no more cracks can be produced (Figures 18(i) and 19(i)). However, the rock fragments will drop to the floor to form a fragment pile. Since the designed new tunnel wall along the borehole has been produced, the drop process of the rock fragments is not simulated to save the calculation time.
Figure 20 depicts the newly formed tunnel wall and the EDZ induced by smooth blasting. The newly formed tunnel wall is denoted by the red dash line while between the red dash line and black dash line is the blast-induced EDZ. It can be seen that there are three zones, i.e., inner damage zone, outer damage zone, and stress redistribution zone. Rock mass in the inner damage zone is characterized by visually blast-induced fine and dense cracks while the outer damage zone is featured by long cracks. The modelled result agrees with that in the literature reviewed in the Introduction Section and Figure 2. Additionally, Figure 20 demonstrates that the fine cracks are mainly concentrated on the immediate area of the newly formed tunnel wall which agrees with Golshani et al. [39] who indicated that blast-induced damage is mostly concentrated in the vicinity of the side wall of the circular opening. Martino and Chandler [9] presented the same conclusion from the measurements of ultrasonic wave velocities in the damage zone surrounding the excavation at Underground Research Laboratories (URLs) of Atomic Energy of Canada, and they concluded that the depth of cracks of EDZ is at least 1 m into the rock mass. The modelled result by the hybrid finite-discrete element method (Figure 20) illustrates that the inner damage zone is approximately 0.3 m and the out damage zone is roughly 1.5 m. The modelled size of EDZ generally shows agreements with the literature [9]. However, it is should be noted that the sizes of EDZ highly depend on the rock properties, initial gas pressure, and in situ stress.
4.1.2. Evolution of Minor Principal Stresses at Monitoring Points
Figure 21 records the evolution of the minor principal stresses at specific monitoring points A, B, and C located in the inner damage zone, the outer damage zone, and the disturbed zone, respectively, induced by the smooth blasting in the top heading of the tunnel.
In this modelling, the in situ stresses are firstly applied on boundaries before the explosives in the boreholes are detonated. Thus, the closer for the minoring points to the boundaries, the earlier the stresses initiate. At around 0.5 ms, the stresses at points C, B, and A reach their first peaks, successively, and the stresses at these points keep constant roughly until the blast wave achieved, which means that the in situ stress in the rock mass achieves equilibrium. Since the boreholes are detonated at 1.45 ms, the stress at point A increases rapidly at around 1.47 ms to its peak and drops sharply. Finally, stress at point A drops to around zero. The stress at point B also follows the same trend as stress at point A. Stress at point B rises quickly to its peak and drops gradually to the bottom. After 2.25 ms, the stress at point B becomes a constant value. While the stress at point B is about to decrease, the stress at point C begins to increase to its peak and then gradually drops to a constant value too. The final constant stresses at points C, B, and A indicate that stresses there achieve another equilibrium after the smooth blasting. It is also illustrated from Figure 21 that the amplitudes of the stress curves at the three points are dependent on the distances from the monitoring points to the borehole. Consequently, the closer to the borehole, the higher the stress curve amplitude is.
In short, due to the in situ stresses are applied on the boundaries first and the boreholes are ignited at 1.45 ms, the monitored stresses are significantly related to the distances to the boundaries before the in situ stress achieves equilibrium. Then, the stresses at the monitoring points keep almost constant before the explosives in the boreholes are detonated. After 1.45 ms, the stresses are significantly influenced by the smooth blasting and the stresses fluctuation times at the monitoring points due to the rock blasting depend on the distance between the monitoring points and the nearest borehole. Furthermore, the constant values of stresses at monitoring points from 0.5 ms to 1.45 ms demonstrate that the stress in the model achieves equilibrium. Therefore, according to the variables at the monitoring points, the modelled results of excavation by blasts are reasonable.
4.2. Modelling of the Excavation Process and Resultant EDZ Formation Process by Blasts in the Bench of the Deep Tunnel
In this section, the EDZ formation process in the bench of the tunnel is modelled. Figure 22 demonstrates the geometrical model and the numerical model. Figure 23 visually illustrates the minor principal stress evolution process while Figure 24 depicts the fracture initiation and propagation process during the smooth blasting at the bench of the tunnel.
(a)
(b)
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
As shown in Figure 22(a), the boreholes with diameters of 0.1 m are set around the bench. The distance among boreholes is 0.4 m, and burdens for boreholes are 0.5 m. As illustrated in Figure 22(b), triangular elements are employed to discretize the model, and dense elements are employed in the interested area. In this model, 3566 elements and 1877 nodes are employed. The bottom is fixed in both x and y direction while other boundaries are set in situ stress, e.g., in the horizontal direction and in the vertical direction. The properties of the rock follows that in Table 2. The in situ stresses are applied to the boundaries to achieve stress equilibrium before the explosives in boreholes are detonated. To make sure the central difference explicit time integration scheme stable, the time step size of is adopted. The stress propagation, crack initiation, and propagation in the bench of the tunnel during blasting are similar to that in Section 4.1. As the in situ stress is applied on the left, right, and top boundaries, the in situ stress propagates to the tunnel and meets the tunnel wall at 1.40 ms as shown in Figure 23(b). After the stress redistribution, the in situ stress equilibrium is achieved as indicated in Figure 23(c). Then, the explosives in the boreholes are ignited at 2.05 ms. The contained chemical charge in the borehole reacts rapidly, which results in high-pressure gas. While the high-pressure gas expands and interacts with the rock materials surrounding the borehole, intense stress waves from each borehole are emitted into the rock through the borehole walls and propagate radially out of the boreholes. Due to the limited distance between two boreholes, i.e., 0.4 m, the stress waves from the adjacent boreholes overlap and interact with each other (Figure 23(d)). Meanwhile, the stresses reach the free surface, i.e., the original tunnel wall of the bench. During this period, in the immediate vicinity of the borehole, cracks are produced due to stress wave exceeding the rock strength as illustrated in Figure 24(c). However, in order to obtain a desirable newly formed contour in smooth blasting, a relatively lower pressure gas of 0.5 GPa is adopted to avoid extensive damage. Consequently, a crushed zone is not formed. Instead, long cracks are produced directly, which interact with cracks from adjacent boreholes. While the stress wave continues to propagate (Figures 23(e) and 23(f)), boreholes are connected by the produced cracks and the big chamber is created. Due to the high-pressure gas expanding, cracks mainly propagate along the curve path approximately parallel to the contour of the bench (Figure 24(d)). The compressive stresses produced from the boreholes are reflected from the free surface (Figures 23(f)–23(h)). As a result, tensile stress at the vicinity of the free surface is produced. As the stress wave attenuates with time going by, high-pressure gas plays an important role in producing cracks and expands exiting cracks, which results in more cracks in the excavated rock mass as can be seen from Figures 24(f)–24(h). Finally, the excavated rock mass is lifted up and pushed away by the high-pressure gas as illustrated in Figure 24(i). It should be noted that long cracks at the two bottom corners are produced and propagate upwards probably due to the stress concentration at those two corners.
Figure 25 depicts the final contour around the tunnel. The immediate vicinity of the contour, where cracks are shown, is the EDZ. It is comprised of some attached rocks on the newly formed tunnel wall, intensive damage, and long cracks disturbed along the contour of the bench. Again, the modelling results show a good agreement with the well-documented literature reviewed in Introduction Section and Figure 2.
5. Discussion
In terms of excavation underground, there are lots of factors influencing the crack initiation and propagation, eventually influencing the shape and size of the EDZ. Among all the factors, the in situ stress, initial gas pressure, and borehole spacing are of great importance. Thus, these factors are discussed in this section.
5.1. Effect of the In Situ Stress
In this section, the last step of excavation for the top heading and bench of a tunnel are studied in two boundary conditions with and without in situ stresses. An initial gas pressure of 500 MPa is adopted for the two scenarios, i.e., with and without in situ stresses, to model the crack initiation and propagation. Figure 26 compares the modelled results at 1.5 ms after the explosives in the boreholes are detonated. For all the modelled results, the clear blast-induced contours, i.e., newly formed tunnel walls, are obtained. It is found that, under the two scenarios, the shapes of newly formed tunnels are similar and the rock fragment sizes and size distribution in excavated rock mass are of little difference.
(a)
(b)
(c)
(d)
Besides, there are two obvious differences between the modelled results in the two cases. Firstly, the most obvious difference is that the cracks in the scenario of without in situ stress are much longer than that in the case of with in situ stress probably due to the in situ stress increasing the resistance of the rock to fracture. Secondly, the cracks mainly propagate radially in the case of without in situ stress while under the scenario of with in situ stress, the cracks approximately propagate horizontally. It seems that the crack propagation in the vertical direction is suppressed due to the in situ stress.
Therefore, it is concluded that the in situ stress plays an important role in the excavation by blasts. For the same initial gas pressure, in situ stress suppresses the propagation of cracks. Thus, due to the suppression of the in situ stress, the resultant EDZ is relatively smaller compared with the case without in situ stress.
5.2. Effect of the Initial Gas Pressure for Top Heading
In the modelling of the crack propagation by blasts in both the top heading and bench of the tunnel, three initial gas pressures, i.e., 100 MPa, 500 MPa, and 1 GPa, are applied in the boreholes to model the fracture and fragmentation processes during excavation by blasts. Figure 27 depicts the modelling results of both top heading and bench at 3.00 ms with three initial gas pressures, respectively.
(a)
(b)
It is found that the excavation highly depends on the initial gas pressure. In the case of 100 MPa (Figures 27(a)-i and 27(b)-i), although a clear contour is formed, the burden for both top heading and bench are not broken into fragments and the initial gas pressure is not able to push the burdens away. As a result, the tunnel fails to be excavated by smooth blasting due to low initial gas pressure. In the case of 500 MPa (Figures 27(a)-ii and 27(b)-ii), the interaction of cracks from adjacent boreholes forms a big chamber and the high-pressure gas in it expands to separate the burden and the newly formed tunnel wall. Moreover, due to relatively high-pressure gas adopted, the burdens are broken into fragments. In the case of 1 GPa, relatively big damaged zones are caused. Dense and long cracks are observed. More shear failures are produced, which is denoted by red colour.
Therefore, it is concluded that it is essential to choose an appropriate initial gas pressure not only to form a new tunnel wall as designed but also to avoid producing big EDZ.
5.3. Effect of the Spacing between Boreholes
For smooth blasting, the spacing between boreholes plays an important role. As indicated by Khoshrou and Mohanty [40], the spacing between boreholes should not exceed 0.8 times the burden for effective smoothing blasting. To demonstrate the effect of spacing between boreholes, two spacings, i.e., 0.4 m and 0.8 m, are adopted for smoothing blasting, and the modelled results are compared. In Figures 28(a) and 28(b), on the left half of the tunnel, the blue represents the rock mass while the white indicates the cracks. On the right half of the tunnel, the lines represent cracks while the white represents rock mass. On the left half of the two models, only the major cracks with relatively large separations are shown. Figure 28 demonstrates the results at 0.3 ms when the blasting processes have almost finished and the fragments are going to fall down. In the case of smoothing blasting with 0.4 m spacing, it produces a newly tunnel wall which passes through the boreholes and approximately parallel to the original tunnel wall as shown in the left half of Figure 28(a). However, in the case of smooth blasting with 0.8 m spacing, as shown in Figure 28(b), the produced tunnel wall is not as smooth as that with closer spacing, i.e., 0.4 m. In addition, the burden is broken into bigger fragments (Figure 28(b)) compared with that for low spacing smoothing blasting (Figure 28(a))
(a)
(b)
It is easier for boreholes to join together during crack propagation for lower spacing smooth blasting. Therefore, the quick joining of boreholes in a row permits release of excess gaseous energy through the cracks joining the boreholes and prevents extension of cracks in the periphery wall and the cumulative force generated from each borehole help to the move the burden rock [41]. That is why lower spacing smooth blasting produces a regular tunnel wall.
It also can be seen from the right half of the tunnel in Figure 28 that lower spacing smoothing blasting (Figure 28(a)) produces more small fragments. There are more boreholes for lower spacing blasting and the boreholes and the cracks initiating from them provided free surfaces for adjacent borehole blasting. Therefore, more explosive and free surface help to produce fine fragments.
The EDZ produced by lower spacing blasting is smaller than longer spacing blasting because for the lower spacing blasting, more gaseous is used to move the burden rock instead of producing long cracks. Therefore, the spacing should be considered in the design stage to achieve a desirable tunnel wall.
6. Conclusion
In this study, the hybrid finite-discrete element method is implemented to model the rock fracture process in a deep tunnel and the resultant excavation damage zone (EDZ) induced by smooth blasting. The previous study on numerical modelling EDZ is firstly reviewed to compare the advantages and disadvantages of numerical methods available in terms of simulating stress distribution and redistribution, crack initiation, and propagation during EDZ formation process around a tunnel. Then, the main component of hybrid finite-discrete element method, i.e., transition from continuum to discontinuum, which makes hybrid method superior to the continuum-based method and discontinuum-based method, is introduced. Additionally, the detonation-induced gas expansion and flow through fracturing rock and the dependence of the rock dynamic behaviour on the loading rate are taken into account. After that, the hybrid finite-discrete element method is implemented to model the rock fracture process in the uniaxial compression strength (UCS) test and the Brazilian tensile strength (BTS) test. The modelled results are compared with the experimental results to calibrate the hybrid finite-discrete element method. Afterwards, the hybrid finite-discrete element method is employed to study the rock fracture process and the EDZ formation process induced by blasts in a deep tunnel. In the end, a variety of factors influencing the propagation of the rock fracture, the final shape of the newly formed tunnel wall, and the shapes and sizes of EDZ are discussed. According to the research mentioned above, it is concluded as follows:(i)The hybrid finite-discrete element method has well modelled the transition of rock from continuum to discontinuum through the blast-induced fracture and fragmentation. The hybrid method has even well modelled the whole EDZ formation process by smoothing blasting, which includes jointing adjacent boreholes by radial cracks from boreholes, producing fragments in the burden rock mass, and moving the burden rock mass by high-pressure gas. The comparison of the modelled outer damage zone and the inner damage zone of EDZ with those well documented in the literature indicates that the hybrid finite-discrete element method can reproduce the EDZ formation process and capture the main characteristics of EDZ.(ii)The ability of hybrid finite-discrete element method modelling the excavation process by blasts is because of the key component, i.e., transition from continuum to discontinuum through rock fracture and fragmentation, which makes the hybrid finite-discrete element method superior to the traditional continuum-based finite element method and discontinuum-based discrete element method.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This research was partly supported by the Scientific Research Fund Project of the Yunnan Provincial Department of Education (grant no. 2020J0051), funding for new arrival lecturers of the Kunming University of Science and Technology (grant no. KKSY201867017), and a two-year visiting PhD scholarship provided by the China Scholarship Council (CSC), which are greatly appreciated.