Blasting serves as a widely employed technique for rock fragmentation in coal mining operations. However, natural rock formations, including coal deposits, often exhibit various discontinuities such as voids, joints, and faults, which have a substantial influence on the propagation and attenuation of stress waves (Huang et al. 2015; Wang et al. 2022; Raina 2019). Consequently, these discontinuities can significantly alter both the mode of crack propagation and the ultimate blasting effect (He et al. 2021). The comprehension of stress distribution and crack propagation within jointed rock masses becomes imperative when investigating dynamic fracture mechanisms under blasting loads in the realm of coal science.

In recent years, various theoretical studies have been conducted on the propagation of explosion-induced cracks in rock masses. From the perspective of the rock medium’s response, the dominant factor leading to blasting failure is often attributed to the tensile effect resulting from the reflection of blast stress waves (Xie et al. 2017; Yang et al. 2017). Due to the lower tensile strength of rock-like materials compared to compression, circumferential cracks are frequently initiated under tensile stress, particularly in the vicinity of free surface areas (Banadaki and Mohanty 2012; Xu et al. 2021). To examine the mechanism of blasting-induced cracks under the influence of blasting stress waves, various failure criteria have been employed, such as the Mohr-Coulomb criterion and the maximum tensile stress criterion. Ming et al. (2009) and Li and Ma (2010) established the wave propagation equation to study the interaction between shock waves and rock joints, taking into account the influence of confining pressure in deep rock masses. In addition to the aforementioned reflected tension waves, explosive gases have been recognized as playing a significant role in rock blasting (Yang et al. 2016; Yi et al. 2016; Yuan et al. 2019). These gases penetrate into cracks formed by shock waves, expand in the confined spaces, and generate pressure loading, thereby promoting and accelerating crack propagation. Previous studies (McHugh 1983; Daehnke et al. 1997) have indicated that the effect of gas pressure may even be more pronounced than that of stress waves in containment conditions. Furthermore, several dynamic fragmentation models based on fracture mechanics have been proposed to describe the morphology of rock cracks and evaluate the impact of blasting. Torbica and Lapcevic (2014) improved a fragmentation model for rock fractures, enabling the estimation of the density of tensile cracks and the radius of the fracture zone. Dai et al. (2010) and Zhang and Zhao (2013) applied quasi-static fracture mechanics theory to rock testing, providing a means to determine dynamic stress intensity factors. The accumulated knowledge from these studies contributes to a comprehensive understanding of rock blasting and paves the way for improved techniques and safety measures in the field.

However, the prediction and analysis of cracks become more complex when considering the heterogeneity of the rock mass, especially with the inclusion of joints. In response to this challenge, numerous researchers have conducted laboratory experiments on rocks (Zhang et al. 2021; Chi et al. 2019; Wang et al. 2020) or rock-like materials (Yang et al. 2008; Wang et al. 2016) to investigate the process of fragment bursting under various complex conditions. Theoretical studies of rock blasting are generally consistent with the observed crack propagation phenomenon in tests. However, the prediction and analysis of cracks become more complex when considering the heterogeneity of the rock mass, especially with the inclusion of joints. In response to this challenge, numerous researchers have conducted laboratory experiments on rocks (Zhang et al. 2021; Chi et al. 2019; Wang et al. 2020) or rock-like materials (Yang et al. 2008; Wang et al. 2016) to investigate the process of fragment bursting under various complex conditions. Zhang et al. (2021, 2020b) studied the effect of cylindrical granite samples on rock fracture under different stemming conditions. Chi et al. (2019) reported a series of cylinder tests on granite cylinders with a central charge in a closed explosion room to study the pressure and attenuation of shock waves in granite. In Wang’s research (Wang et al. 2020), the damage and degree of damage of small-sized specimens after blasting with an uncoupled charge in the chamber were qualitatively assessed, and the relationship between the uncoupled coefficient and the blasting effect was examined. Furthermore, some experiments (Liu and Qu 2014; Li et al. 2018, 2021; Miranda et al. 2013; Chen et al. 2018) were conducted on jointed rock samples to study the effects of discontinuous structures on blasting failure and the propagation laws of stress waves. Using PAMM as the experimental material, Yang et al. (2008) and Wang et al. (2016) designed experiments on prefabricated jointed rock masses to study the crack propagation behavior at the joint tips and the distribution of the surrounding stress field at the end of defects through the dynamic caustics method. Yang (2008) proposed a field test procedure for jointed rock masses and concluded that joints significantly affect the characteristics of stress wave propagation, such as peak attenuation, frequency spectrum, and spatial variation. Renshu et al. (2016) and Hao et al. (20012001) studied the shielding effect of joints on crack propagation through laboratory experiments and obtained the dynamic stress intensity factor formed at the joints and the propagation trend of the crack propagation speed under the action of blasting stress.

In recent years, various numerical methods, including the Finite Element Method (FEM) (Liu et al. 2019; Chao et al. 2020), Finite Difference Method (FDM) (Zhu et al. 2007, 2008; Wang et al. 2017), Discrete Element Method (DEM) (Hajibagherpour et al. 2020), Smoothed Particle Dynamics (SPH) (Gharehdash et al. 2020), hybrid methods (Yuan et al. 2018; Wang et al. 2021), and non-local damage models (Liu et al. 2022; Wang et al. 2023), have been widely employed to study the dynamic failure response of rock masses. Researchers such as Zhang et al. (2020a) and Siamaki et al. (2018) have utilized the particle-based DEM (PFC) methodology to simulate rock blasting and analyze the effect of blasting vibration on the degradation of joint shear strength. While discrete methods are suitable for modeling cracks, they require additional parameter calibration and are challenging to extend to elastoplastic behaviors. Continuous methods such as FEM and FDM make up for the deficiencies of discrete element methods in the complex constitutive model, allowing them to reproduce the blast processing of complex structures and the nonlinear characteristics of materials (Dai et al. 2022). Furthermore, hybrid methods combine the advantages of both continuous and discrete methods, making them a valuable tool for studying rock blasts (Wang and Konietzky 2009). In this work, we focus on the FEM approach implemented in the commercial software LS-DYNA, which has become one of the promising numerical tools for studying rock blast damage (Hongtao et al. 2015; Ainalis et al. 2017; Zhang et al. 2016; Yi et al. 2021). For this numerical method, a reasonable damage evolution law and constitutive model are crucial in dealing with discontinuous structures such as nucleation, propagation, coalescence, and fragmentation. At present, a series of constitutive models (e.g., RHT, HJC) have been implemented in the finite element framework to conduct numerical studies on various explosion-induced failures, such as the blasting problem under in-situ stress (Jayasinghe et al. 2019), porous blasting processes (Xie et al. 2017), explosion-induced crack propagation (Yi et al. 2018), and damage evolution and spatial distribution (Wang et al. 2019).

Despite the noteworthy contributions in the field, studying the blasting response of jointed rock mass is still a challenging subject, especially when it comes to the influence of joints on stress wave propagation and the propagation of rock crack damage, which is rarely considered. Hence, this study aimed to address this knowledge gap by conducting a series of numerical experiments on crack propagation in rock masses with varying combinations of double joints. The objective was to investigate the mechanical behavior of crack propagation and the distribution of stress waves in jointed rock masses under blasting loads. By analyzing changes in the stress field and failure patterns caused by different joint arrangements, we reveal the general law of explosive crack propagation and failure models. The results of this study provide a reference for understanding the mechanical and failure characteristics of fissured rock mass defects under blasting loads.

This paper is organized as follows: In Sect. 2, we briefly state the different constitutive equations of the numerical methodology and algorithm of fluid-solid coupling. In Sect. 3, we establish three series of RHT granite models and simulate the failure processes of jointed rocks with different joint arrangements under blasting. Furthermore, we quantitatively explore the effects of joint length, inclination, and distance from the blasting hole on the mechanical and failure characteristics of jointed granite specimens in Sect. 4. Finally, we draw conclusions and summarize our findings in Sect. 5.

2.1 Constitutive models for blasting simulation

For numerical studies of blasting problems involving nonlinearities and large deformations, choosing a reasonable material model is crucial to the numerical results. LS-DYNA software provides many constitutive models for simulating rocks or rock-like materials. In this numerical study, the RHT model and the JWL equation of state will be employed to model the mechanical response of rock under blast loading. These models will be discussed below.

2.1.1 RHT constitutive model for rock material

The RHT constitutive model was proposed by Riedel, Hiermaier, and Thomas in 1998, based on previous research on the mechanical properties of concrete (Sparks and Menzies 1973; Rostasy et al. 1984; John et al. 1992), which embeds pressure-related elastic limit surface equations, failure surface equations, and residual strength surface equations (Shu et al. 2022).The model incorporates various parameters and equations to describe the fracture process, taking into account factors such as stress state, strain rate, and material properties. As illustrated in Fig. 1, the material is assumed to be elastic when the pressure P is lower than the pressure at which the material pores begin to crush \(P_{\text{crush}}\). Once P exceeds \(P_{\text{crush}}\), the pore collapses, resulting in a decrease in the effective volume modulus of the material. If the pressure \(P_{\text{crush}}\) is between P and the pressure at which the material’s pores are compacted \(P_{\text{comp}}\), the material is assumed to be plastic. Moreover, the material will only yield if the load exceeds its yield limit, given a specific stress state and loading rate. It has been widely applied in studies related to impact and explosion dynamics, allowing researchers to better understand and predict the fracture behavior of brittle materials under dynamic loading conditions.

In the RHT model, the equations of state for \(p-\alpha\) compaction of rock masses are:

where \(P_{R}\) represents the pressure of the equation of state (EOS) used in the RHT model, with material constants \(B_{0}\) and \(B_{1}\), initial porosity of the rock \(\alpha _{0}\), initial density of rock \(\rho _{0}\), internal energy per unit mass e, and volumetric strain \(\mu\). To determine the polynomial coefficients, the following expressions can be used: \(A_{1}\), \(A_{2}\), and \(A_{3}\):

where, \(c_{0}\) is the sound speed at ambient pressure and temperature; \(T_{1}\) is the material constant, and k is the empirical constant of the material.

RHT model (Borrvall and Riedel 2011): a Schematic description of the \(p-\alpha\) EOS; b Stress limit surfaces and loading scenario

The RHT model also includes a description of the hydrostatic pressure using the polynomial Hugoniot curve and the \(p-\alpha\) compaction relationship, as shown in Fig. 1b. To describe the model’s yield surfaces, compressive strength, regularized yield functions, and William-Warnke functions are utilized (Riedel et al. 1999):

where, \(\sigma _{y}^{*}\) denotes the standard yield function; \(f_{c}\) stands for the uniaxial compressive strength; \(R_{3}\) represents the William-Warnke function, and \(\theta _{l}\) denotes the pulse angle. Furthermore, the equation includes the dynamic strain rate increase factor, \(F_{r}\), the normalization pressure, \(P_{0}^{*}\), defined as \(P_{0}^{*} = P_{0}/f_{c}\), hydrostatic pressure, \(P_{0}\), strain rate \(\varepsilon _{p}\), and the effective plastic strain rate \({\dot{\varepsilon }}_{p}\).

Thus, Failure surfaces indicate (Johnson and Holmquist 1994):

where, \(\sigma _{f}^{*} = \sigma _{f}/f_{c}\) is the normalized strength; \(A_1\) and \(N_1\) are the failure surface parameters.

The residual strength surface describes the strength of the broken concrete, which is related to the confining pressure state of the material. So, the residual strength surfaces can be described as follows:

where, B is the residual failure surface constant; M is the residual failure surface index, and \(P^{*}\) is the normalized hydrostatic pressure.

As the stress state of a material approaches its ultimate strength on the destruction surface, any further inelastic deformation or plastic strain results in a gradual accumulation of damage. In the event of destruction, the amount of plastic strain can be used to determine the extent of the damage. It should be noted that the plastic strain in case of failure is represented by:

where \(\varepsilon _{p}^{f}\) represent the plastic strain at failure, \(\varepsilon _{p}^{m}\) is the minimum damage residual strain, \(P_{t}^{*}\) is the failure end pressure, \(D_{r}\) is a damage variable for the RHT model, \(D_{1}\) and \(D_{2}\) injury constants.

Finally, the damage variable \(D_{r}\) is defined as the accumulation of plastic strain as:

2.1.2 Constitutive parameters of coupling material air

As a wildly used coupling material, the air model is described as a linear elastic model in this research. It uses the material with the keyword \(*\)MAT\(\_\)NULL, the equation of state is represented by the keyword EOS\(\_\)LINEAR\(\_\)POLYNOMIAL (Hallquist et al. 2006), and the pressure is determined by the following form:

where \(\mu =\rho /\rho _{0}-1\), \(\rho\) is the density, \(\rho _{0}\) is the starting density, \(E_{0}\) is the energy density of the beginning, \(C_{0}\), \(C_{1}\), \(C_{2}\), \(C_{3}\), \(C_{4}\), \(C_{5}\), and \(C_{6}\) are parameters in the equation as listed in Table 1.

2.1.3 Dynamic loading produced by explosives

During blasting, explosives undergo rapid chemical reactions that can be described by various equations of state. One commonly used model is the Jones-Wilkins-Lee (JWL) EOS, which accurately predicts pressure changes resulting from blasting (Hallquist et al. 2006; Alia and Souli 2006). In this study, we employed the JWL\(\_\)EOS to describe the blasting process, and the blast pressure was calculated using the following equation:

where, V is the relative volume; E is the initial internal energy per unit volume of explosives, A, B, \(R_{1}\), \(R_{2}\), and \(\omega\) are all parameters of the JWL equation of state as listed in Table 2.

2.1.4 Material parameters of joint

The keyword \(*\)MAT\(\_\)PLASTIC\(\_\)KINEMATIC describes isotropic hardening and follow-up hardening plastic models, commonly used as the constitutive structure of rock in rock mass blasting simulations (Hallquist et al. 2006), and the radius of the yield surface under this model is the initial yield strength plus hardening reads:

\(\sigma _{y}\) is yield strength, \(\sigma _{0}\) is the initial yield strength, \(E_{p}\) is a plastic hardening modulus, \(\varepsilon _{eff}^{p}\) is an effective plastic strain.

The modulus of plastic hardening is given by the following equation

E is the modulus of elasticity and the \(E_{t}\) is the tangent modulus.

According to the different \(\beta\) values, it can be used to describe various hardening models: when \(\beta = 0\), it is follow-up hardening, the yield surface size is unchanged, and it moves in the direction of plastic strain; When \(\beta = 1\), it is isotropic hardening, the position of the yield surface is unchanged, and the size changes with strain, at \(0 \le \beta \le 1\), it is mixed hardening. Moreover, the Cowper–Symonds model is used in the \(*\)MAT\(\_\)PLASTIC\(\_\)KINEMATIC model to consider the effect of strain rate as follows:

where, p and c are the parameters related to the strain rate, which is determined by the strain rate characteristics of the material; \({\dot{\varepsilon }}\) is the load strain rate.

Then the \(*\)MAT_PLASTIC_KINEMATIC model can be expressed as follows:

The specific parameters for joint material are listed in Table 3.

2.2 Algorithm of fluid–solid coupling

The finite element method (FEM) is a powerful and robust numerical technique for simulating the progressive failure of rock materials and revealing the fracture and fracture mechanisms of rocks under dynamic loading (Zhang et al. 2019; Liu et al. 2021). In FEM, there are three basic algorithms for three-dimensional elements: Lagrangian, Eulerian, and Arbitrary Lagrange-Eulerian Algorithm (ALE), which are controlled by THE ELFORM keyword in \(*\)SECTION\(\_\)SOLID. The Lagrange algorithm attaches elements to the material, enabling deformation to be generated with the flow of the material. However, numerical difficulties can arise when the structure undergoes large deformation, particularly due to high-stress gradients at the failure interface, which often causes the simulation to break down before analysis is complete because of convergence issues associated with topological changes. In contrast, the Euler algorithm is suitable for computational analysis of fluid mechanics and explosion of explosives in the rock mass, but it results in a longer and more complicated calculation process for the model. Therefore, to overcome the limitations of the Lagrangian and Eulerian formulations, the ALE algorithm, which couples both formulations, is used for fluid–structure interaction analysis in this research.

The basic principles of the ALE algorithm involve performing one or several Lagrange time-step calculations, during which the element mesh deforms with the flow of material, followed by an ALE time-step calculation consisting of (1) maintaining the boundary conditions after the object’s deformation and remeshing the interior, and (2) transporting the element variables and node velocity vectors from the deformed mesh to the repartitioned mesh.

Algorithm of the fluid-solid coupling with explicit time integration in ALE

As shown in Fig. 2, the rock is solid, and explosives and air are fluids. Therefore, The explosive and air parts are controlled by ALE formulation based on the common nodes, while the rock belongs to the Lagrange algorithm. If only explosives and rocks are set up, the rocks will not be damaged, and the model is also a failure. The coupling between explosives, air, and rocks is the formation of a domain of air between rocks and explosives as a medium. The energy waves generated by the explosion of explosives are transmitted to the rock through the air medium, thus acting on the rock and causing destructive damage to rock. Currently, coupling algorithms have been widely used in finite element analysis, especially in the explosion of structures. The interaction effect between structural elements and fluids using the LAGRANGE algorithm is achieved by defining the keyword \(*\)CONSTRAINED\(\_\)LAGRANGE\(\_\)IN\(\_\)SOLID, and the numerical model of the target structure and fluid and the element grid can overlap. The methods of coupling the target structure to the fluid are: acceleration, velocity, penalty function, and other constraints to accomplish the interaction between the fluid and the solid.

In this section, we explore the blasting response of a jointed rock mass using the constitutive model and parameters introduced in Sect. 2. First, we present a benchmark of granite specimens with a single borehole for numerical validation. Then, we conduct a series of numerical experiments with two joints to investigate the effect of joint arrangements. All numerical experiments are carried out under a two-dimensional plane strain assumption.

3.1 Calibration of the RHT model

For the first case, we conducted a numerical simulation of single-hole blasting on a standard square specimen made of granite. The specimen has a length of 300 mm and is discretized into 22,500 quadratic quadrilateral elements with uniform size (\(h_{e} = 2\) mm). Additionally, we created a circular air domain consisting of approximately 17,663 elements and an explosion source consisting of 4 elements, as shown in Fig. 3. To prevent the propagation of reflected tensile cracks, we applied non-reflecting boundary conditions around the model.

Schematic diagram for borehole blasting experiment: a geometric and boundary conditions, b ultimate failure pattern predicted by numerical test, and c the failure pattern observed in experiment of the Barre granite (Banadaki and Mohanty 2012)

The geometric setup for the borehole blasting experiment, including the rock, explosives, and air, is depicted in Fig. 3a. Microscopic parameters of the simulated granite’s RHT model are determined in Table 4 based on previous laboratory experiments using the split Hopkinson pressure bar (SHPB) test of granite (Li 2016).

The predicted crack growth paths under explosive loading conditions are shown in Fig. 3b. As depicted in the figure, a crushed zone emerges around the borehole, and four radial cracks propagate towards the outer boundary of the rock specimen. Moreover, the ultimate failure patterns of the Barre granite in the blasting test obtained from the previous experimental observations (Banadaki and Mohanty 2012) are compared to the numerical results, as demonstrated in Fig. 3c. It can be seen that the presented numerical results are in good agreement with the previous experimental findings, which suggests that our numerical model and its parameters are capable of reproducing the mechanical and failure behaviors of granite.

3.2 Setup for the blasting experiment

In this numerical study, parallelepiped granite specimens with dimensions of \(400 \, {\text{mm}} \times 500 \, {\text{mm}} \times 2 \, {\text{mm}}\) are used. The filling medium is also granite. Two prefabricated joints are created in the intact granite specimen by replacing the material model. The front joint (Joint A) is located at the center of the specimen and has a width of 4 mm. As shown in Fig. 4a, the length of Joint A is denoted by L, the distance from the center of Joint A is denoted by D, and \(\alpha\) represents the joint inclination angle. Another horizontal joint (Joint B) is located 30 mm away from Joint A to study the dynamic propagation of rock defects under diffraction waves. Its length and width are fixed at 80 mm and 4 mm, respectively. An explosive (2.7 g/cm\(^3\) and 4 mm in diameter) is placed on the left side of Joint A to simulate the formation of blast loads. To eliminate the stress wave reflection from the outer boundaries, a non-reflection boundary is adopted in the outer boundaries of the model. In the numerical simulation, 50,000 8-node elements for a total of 400,000 nodes with uniform size (\(h = 2\) mm) are used to discretize the computational domain, as shown in Fig. 4b. Four monitoring points (A1, A2, B1, B2) are set at the wing end element with two prefabricated joints.

In this numerical research, the stress wave propagation process is accomplished by setting a coupling of rock, air, explosives, and joints. Since the discontinuous structural planes in natural rock masses usually appear in the form of closed fractures, filler-based joints are used that are first dug out and then filled with medium.

Schematic diagram for borehole blasting experiment. a Geometric and boundary conditions and b the finite element discretization

In order to study the fracture characteristics of joint rock mass under blasting load, a total of 13 study cases were studied, and the following three groups of parameters were considered, as listed in Table 5.

1. (1)

   When the inclination angle \(\alpha\) of Joint A is fixed at \(90^{\circ }\) and blast hole distance D is 50 mm, increase the length L of Joint A from 0 to 80 mm with an interval of 20;

2. (2)

   Joint length \(L = 60\) mm and blast hole distance \(D = 50\) mm maintained constant, change the inclination angle of Joint A \(\alpha\) as \(0^{\circ }\), \(30^{\circ }\), \(45^{\circ }\), \(60^{\circ }\), and \(90^{\circ }\);

3. (3)

   Fix the joint length L as 80 mm and inclination angle \(\alpha\) as \(90^{\circ }\), take blast hole distance D in 25 mm, 50 mm, 70 mm, 90 mm, and 100 mm.

3.3 Results of numerical experiment

3.3.1 Fracture patterns and stress distribution

For the case of \(L = 80\) mm, \(\alpha = 90^{\circ }\), and \(D = 50\) mm, the evolution of crack growth paths, maximum shear stress, and maximum principal stress fields during blasting is presented in Fig. 5. The blasting fracture process is characterized by four stages, namely the compressive, crushing, propagating, and failure stages. As shown in Fig. 5a, the shock wave generated by the explosion initially causes extensive plastic deformation in the radial direction due to the high strength of the rock, resulting in compressive failure. This deformation leads to the formation of a compressive crushed zone near the blast hole in stage I (t = 15.998 μs). As the blasting loads increase (t = 15.998–23.796 μs), the stress wave propagates to Joint B, leading to a reflected tensile wave. The low tensile strength of the rock induces plastic tensile deformation, forming a tensile failure zone between the blast hole and Joint B. Stress waves that are not affected by the joints propagate along the specimen, forming the main crack, which propagates radially. When the stress wave reaches the tip of Joint B, a difference in deformation on both sides leads to a deformation gradient at the joint ends, causing a tensile shear crack. Three radial and two wing cracks propagate along the axial direction in stage III (t = 67.996 μs) as the stress wave attenuates. Finally, the stress waves and damaged cracks reach the surface, and the cracks branch after propagating a certain distance due to the coalescence of the stress waves, as shown in Fig. 5. The transmitted wave leads to the production of two bifurcated cracks at the right end of Joint B.

The spatiotemporal distribution of a crack growth paths; b maximum shear stress fields (unit: MPa); and c maximum principal stress fields (unit: MPa) at different stages

The data from the monitoring points in the granite specimen during the failure process are shown in Fig. 6. Four representative stages of the failure process are plotted in the history curves. The results show that points A1 and A2 have similar trends due to their symmetry. When the stress wave reaches Joint A, the monitoring points A1 and A2 both produce rightward horizontal displacements and reach peak values at 25.6 μs. As the stress wave propagates, the horizontal displacement at these points changes from positive to negative, indicating the formation of a reflected tensile wave. However, the displacements of points B1 and B2 show a monotonically increasing trend. Figure 6b shows the horizontal velocity of the monitoring points. The velocity of points A1 and A2 reaches a peak value of 2.04 (μm/μs) at 13.9 μs and then drops to state II due to the rebound effect. The stress wave then converges again, and the velocity of the two measuring points increases significantly, reaching the peak value again. The same trend can also be observed for monitoring points B1 and B2, but with a time hysteresis.

Evolution of a horizontal displacement, b horizontal velocity, and c horizontal acceleration at the monitoring points in numerical simulations

The effective stress and maximum principal stress at the monitoring points recorded in the simulation are plotted in Fig. 7a, b, respectively. The stress attenuation is clearly visible, and the history curves of the two stresses have the same characteristic. The stress peaks at about 2.0 MPa at point A, while at point G, the peak stress drops to below 0.2 MPa. This indicates the shielding effect of joints on the propagation of stress waves.The recorded strain energy and kinetic energy history are shown in Fig. 7c. The elastic strain energy rises rapidly at the initial load due to the high pressure of the blasting load. With further loading, the kinetic energy decreases gradually, while the strain energy tends to a constant value of 0.6 kJ.

The evolution process of a effective stress, b maximum principal stress, and c rock energy in granite under the blasting conditions

The peak stress response at various locations under different front joints is as follows: at the A1 end, it occurs within the range of 5–10 μs; at the A2 end, it manifests at approximately 7–8 μs; at the B1 end, it peaks between 19 and 24 μs; and at the B2 end, it reaches its maximum around 38–39 μs. Similarly, for different front joint angles, the peak stress at the A1 end is observed within 4–6 μs, while at the A2 end it ranges from 18 to 20 μs. At the B1 end, it occurs between 20 and 24 μs, and at the B2 end, it is observed at approximately 37–40 μs. In terms of blasting hole spacing, the peak stress at the A1 end of the joint spans from 6 to 23 μs, at the A2 end it ranges from 7 to 24 μs, at the B1 end it falls between 45 and 50 μs, and at the B2 end it is recorded within the range of 57–60 μs.

3.3.2 Blasting with different joint length

To study the influence of the length L of Joint A on mechanical and failure characteristics, we compared the numerical results of five jointed granite specimens with the same configuration except for L, as shown in Fig. 8. It can be observed from Fig. 8 that four main cracks in the shape of a cross are initiated and grow. When the main crack propagates in the opposite direction of Joint A, the crack is hardly affected by the joint, and a branching phenomenon is inclined to occur at the adjacent boundary. A similar phenomenon can be observed for those main cracks that propagate perpendicular to the joint direction, except that the crack propagation direction is deflected as the length of Joint A increases. Crack initiation is observed at different time intervals within the system: the crack in the crushed zone occurs at 3.89 μs, followed by cracks at 6 μs, the joint zone at 12 μs, and finally at 18 μs.

A significant difference was found when the main crack impinges on the prefabricated joints. When the length of Joint A is set to 0 mm (i.e., no front joint exists), the main crack propagates directly along Joint B and bifurcates at the end of Joint B. When the length of the joint is 20 mm, the stress wave generates a reflected tensile wave and a transmitted wave when it reaches Joint A. The reflected tensile wave induces many microcracks and forms a triangular ring crack zone (Friedrich and Schlarb 2011) on the joint surface, while the transmitted wave allows the main crack to further propagate to Joint B in the original direction. As the joint length increases to 40 mm, this phenomenon becomes more pronounced: the ring crack zone becomes larger and connects with the blast hole. Furthermore, as shown in Fig. 8c, two horizontal cracks propagate from the end of Joint A. When the joint length is increased to 60 or 80 mm, the transmission of the main crack becomes less evident, although the bifurcation of the crack can still be observed at the end of Joint B. The triangular ring crack zone formed by the microcracks gradually increases into a rectangle during this stage. Furthermore, the inclination of the wing cracks becomes larger with increasing joint length due to the deflection of the stress wave propagation direction.

Ultimate failure patterns after blasting of rock specimens with different joint length L: a \(L = 0\) mm, b \(L = 20\) mm, c \(L = 40\) mm, d \(L = 60\) mm, and e \(L = 80\) mm

3.3.3 Blasting with different joint inclination

Crack propagation and joint interaction are simulated for five different inclination angles (\(0^{\circ }\), \(30^{\circ }\), \(45^{\circ }\), \(60^{\circ }\), and \(90^{\circ }\)) of Joint A, as shown in Fig. 9. When Joint A and Joint B are in contact (\(\alpha = 0^{\circ }\)), the main crack propagates along the joints and bifurcates at the end, similar to the result in the case of \(L = 0\) mm. For an inclination angle of \(30^{\circ }\), more radial cracks are formed due to the uniform energy distribution, as shown in Fig. 9. The upper wing end of Joint A is located in the crushed zone, and no evident crack propagation is observed. However, the main crack is shifted towards Joint A under the influence of this joint end, while the stress wave releases from the lower wing end, creating a more extended wing crack. The main crack generated through the transmitted wave is deflected in the joint’s vertical direction due to the effect of Joint A, but a small crack is still observed from Joint A to Joint B. Increasing the inclination angle \(\alpha\) leads to a transition of the main crack from deflection to penetration, as observed from Fig. 9. Results from the case with an inclination of \(30^{\circ }\) are almost identical to those observed from the case of \(\alpha = 60^{\circ }\) due to rotational symmetry. However, no additional main cracks were found in the case of \(\alpha = 45^{\circ }\). As the number of main cracks increases with the loading rate, it can be inferred that when the inclination angle is at a certain angle, the stress wave reflected by Joint A will enhance the intensity of shock waves around the blast hole.

Ultimate failure patterns after blasting of rock specimens with different inclination angle \(\alpha\) of Joint A: a \(\alpha = 0^{\circ }\), b \(\alpha = 30^{\circ }\), c \(\alpha = 45^{\circ }\), d \(\alpha = 60^{\circ }\), and e \(\alpha = 90^{\circ }\)

3.3.4 Blasting with different distance from joint

Figure 10 illustrates the ultimate failure patterns resulting from blasting at five different blast hole distances D from Joint A, specifically 25 mm, 50 mm, 70 mm, 90 mm, and 100 mm. As shown in Fig. 10, the propagation direction of the vertical main crack is deflected by the reflected tensile wave generated by Joint A. This effect gradually diminishes as the blast distance increases. During this process, the inclination angle of the wing crack continuously decreases due to the reduction in blast hole and Joint A tip angle. Furthermore, a significant difference exists in the results of the ring cracks zone between the blasting hole and Joint A at different blast hole distances D. At a distance of 25 mm, the explosion shock wave interacts with the reflected tensile wave generated by Joint A, resulting in violent microcracks. As this distance increases, the intensity of damage gradually decreases. When the distance reaches 70 mm, the interaction between these two waves begins to decay, demonstrating the spatial independence of the crushed zone and the ring cracks zone. Moreover, when the blasting distance is set to 25 mm, the crushed zone is observed at 12 μs, followed by the appearance of the main crack at 22 μs. As the blasting distance increases, the development time of these two cracks is observed to be delayed. Specifically, at a blasting distance of 100 mm, a crushed zone is observed at 30 μs, with the main crack appearing at 34 μs.

Ultimate failure patterns after blasting of rock specimens with different blast hole distance D: a \(D = 25\) mm, b \(D = 50\) mm, c \(D = 70\) mm, d \(D = 90\) mm, and e \(D = 100\) mm

4.1 Ultimate failure mode

Previous research (Hino 1956; Friedrich and Schlarb 2011) has identified four failure modes near the blast hole that contribute to the rock failure morphology in numerical simulations: the Crushed zone, Ring cracks zone, Main crack (Long radial crack), and short Radial crack. The present study investigated the effects of joints on the blast-induced failure patterns of jointed granite specimens, and three typical ultimate failure patterns were identified, as shown in Fig. 11.

When a joint is located adjacent to the blast hole, the reflected tensile waves generated by the joint can significantly increase the degree of fragmentation near the blast. This effect is manifested by the deflection of the main crack towards the direction of the reflected wave, and the formation of a complete ring cracks zone between the crushed zone and the joint, as shown in Fig. 11a. It is worth noting that the dynamic disturbance of the stress wave at the main crack tip during the propagation process usually bifurcates the main crack, but this effect is suppressed by the reflected tensile wave. When the joint is far from the blast hole, the main crack will no longer deflect since the reflected wave is dissipated in a long propagation path. In this case, the ring cracks zone is formed only on the joint surface. However, there is still a weakened reflected wave that interacts with the crushed zone and forms a short radial crack, as shown in Fig. 11b. When the joint is arranged at a specific angle, additional main cracks can be observed under the disturbance of the stress wave reflected from the joint. Furthermore, due to the effective blast hole distance gradually increasing, the joint only forms a ring cracks zone on the side close to the blast hole.

Three different ultimate failure patterns near the blast hole under the influence of a adjacent joint, b distant joint, and c inclined joint

On the other hand, the failure modes of the main crack impinging and the response of the rear joint are summarized in Fig. 12. It can be observed that when the main crack impinges on the joint surface, a portion of the stress wave can transmit through the surface, resulting in the penetration of the main crack. As the inclination of the joint increases, the main crack will deflect according to the wave interface refraction. However, the propagation of the main crack behind the joint is weaker than in other directions, whether in deflection or penetration. The existence of joints blocks a part of the blasting energy in the form of reflected waves, thus creating a shield effect (Hao et al. 2001; Hutchinson 1987). In addition, some slight secondary cracks initiated from the joints instead of the main crack. Furthermore, the original shear force caused by the difference in stress wave paths on both sides of the joint end will produce a wing crack at the joint ends. According to the experimental results, it was found that the propagation direction of the wing crack is consistent with the direction of the joint end and the blast hole.

Failure modes and responses of prefabricated joints to main crack impact with a vertical joint and b inclined joint

4.2 Effect of joint arrangements

Figure 13 displays the peak displacement, velocity, and acceleration of the monitoring points for granite specimens with varying joint lengths L. The peak displacement of both joint ends decreases as joint length increases, but this effect is limited as most of the displacement is caused by global separation of the specimen after main crack propagation rather than local deformation under shock waves. In contrast, the velocity and acceleration show more significant changes. Comparing point B for cases with \(L=0\) and \(L=20\), a shield effect can be observed, where the peak acceleration at the joint end can be reduced from 9.8 to 7.5 μm/μs² by adding a vertical joint. This suggests that vertical joints can effectively reduce the impact of shock waves at joint ends compared to horizontal joints. When joint length increases to 40 mm, the impact disturbance of Joint B can be reduced to one-tenth of the direct impact case.

Effect of the joint length L on a peak displacement, b peak vibration velocity, and c peak acceleration of the monitoring points under the blasting conditions

Figure 14 shows the recorded kinetic energy and total potential energy history of five granite specimens for joint length study. The energy history curves in each model are consistent in the stage before the stress wave is propagated to the joints (0–10 μs). The kinetic energy decreases after reaching a peak value of 0.914 kJ, a process involving dissipation for crack formation and absorption by the boundary. As shown in Fig. 14a, the larger the L value, the smaller the decrease in the kinetic energy of the rock. The total energy of the rock reaches saturation after a certain period of accumulation and tends to be stable, as seen in Fig. 14b. The rock’s total energy does not change much for different joint lengths.

The evolution process of a the rock kinetic energy and b rock total potential energy in granite specimens under blasting conditions with different joint length L

Figure 15 presents the effects of joint inclination \(\alpha\) on the peak displacement, vibration velocity, and acceleration of monitoring points for granite under blasting conditions. It can be observed in Fig. 15a that point A1 shows singular displacement because it is located in the crushed zone. As the joint inclination \(\alpha\) increases, the peak displacement of the wing end of Joint A decreases steeply, while the peak displacement of the end of Joint B remains relatively stable. When \(\alpha\) is set to \(0^{\circ }\), the velocity and acceleration of each point are the largest, indicating that a joint along the direction of stress wave propagation is disadvantageous. As the angle of Joint A increases, the peak velocity and acceleration both gradually decrease. Moreover, an increase in the inclination angle causes point A1 to gradually move away from the blast hole, while the opposite occurs for point A2. Therefore, the attenuation of point A1 is more significant than that of point A2 during angle changes, as shown in Fig. 15c. It can be concluded that a joint perpendicular to the propagation direction of the stress wave can provide the best shield effect for the rear joint.

Effect of the joint inclination \(\alpha\) on a peak displacement, b peak vibration velocity, and c peak acceleration of the monitoring points under the blasting conditions

To further demonstrate the shield effect of joints, Fig. 16 plots the peak stress of four monitoring points with different joint inclination \(\alpha\). It can be seen that the peak stress at the measuring point of Joint A end ranges from 2.08 to 3.81 MPa. In contrast, the peak stress observed on Joint B is only 0.48 MPa, approximately one-fifth of that on Joint A. In addition, unlike the findings in Fig. 15, the peak stress of different monitoring points is not affected by joint inclination \(\alpha\).

The peak stress of monitoring points with respect to joint inclination \(\alpha\)

4.3 Effect of distance from joint

To gain a better understanding of the effects of blast hole distance D on the mechanical behavior of jointed granite specimens under blasting conditions, Fig. 17 shows the peak displacement, vibration velocity, and acceleration of the monitoring points. It can be observed that the displacement of the jointed end increases with increasing blast hole distance. A snap-back caused by the reflected tensile wave can be observed at the end of joint A when the blast hole distance is set to 25 mm, resulting in a reverse displacement of 40 μm. For the case where D is set as 50 mm, the displacements of four points are similar, indicating that the two joints and the nearby rock mass are displaced as one fragment.

Effect of the blast hole distance D on a peak displacement, b peak vibration velocity, and c peak acceleration of the monitoring points under the blasting conditions

The monitoring results of peak stress in the cases with different blast hole distances are depicted in Fig. 18. It can be observed that as the blast hole distance D increases from 25 mm, the peak stress of points A1 and A2 decreases, while the values of points B1 and B2 remain stable. The results obtained in the study indicate that the effect of blasting load on Joint A decreases with increasing blasting distance, while Joint B is not affected by this process in any significant way.

Effect of the blast hole distance D on the peak stress of granite specimens

In this study, the mechanical behavior of crack propagation and distribution of the stress wave in jointed rock mass under blasting load were numerically explored with the RHT model. The RHT model parameters were calibrated by comparing the result with a single-hole blasting test, based on previous experiments on granite. Subsequently, we performed simulations on three groups of double-jointed numerical specimens with varying joint arrangements, including joint length (L), joint inclination angle (\(\alpha\)), and blast hole distance (D), while subjecting them to blasting loads. By analyzing the numerical calculation results, the following conclusions can be drawn:

1. (1)

   The numerical calculations provide valuable insights into the various ultimate failure patterns observed in granite models subjected to blast loading. In addition to the main crack and crushed zone observed in single-hole blasting tests, the numerical simulations reveal the formation of ring crack zones and short radial cracks between the crushed zone and the joint. These formations occur due to the influence of the reflected tensile wave generated by the joint, which disrupts the propagation direction and number of main cracks. Furthermore, the propagation of the shock wave through the front joint gives rise to wing cracks, main cracks, and secondary cracks behind the joints in the numerical results. Upon encountering the front joint, the main crack deflects and penetrates depending on the joint’s inclination. While a significant portion of the energy is absorbed or reflected by the front joint, a fraction of the stress wave still reaches the subsequent joint and induces secondary cracks. Moreover, shear stress at the joint’s end initiates wing cracks that propagate along the blast hole and the direction of the joint tip.

2. (2)

   The numerical test reveals significant vibration isolation effects of vertical joints on stress waves. The presence of longer joints exhibits a noticeable shield effect on rear joints. The inclination angle of Joint A influences the extent of the failure zone: larger inclination angles result in increased reflection of stress waves, leading to more pronounced reflected tension waves and the promotion of microcracks in the ring cracks zone. Moreover, by adjusting the inclination of Joint A, the number of main cracks can be controlled. Furthermore, the effect of the blasting load on Joint A diminishes as the blasting distance increases, while Joint B remains unaffected during this process. Simultaneously, the ring cracks zone induced by the reflected tensile wave diminishes as the stress wave propagation distance increases.