Abstract
The generation of stress waves induced by explosions underground is governed by material nonlinear responses of materials surrounding explosions and affected by source region mediums and local structures. A nonlinear finite element (NFE) method can simulate the generation efficiently. However, the calculation using the NFE to observational distances, where motions are elastic, is computationally challenging. In order to tackle this problem, we present a subsection numerical simulating method for forward modelling the generation and propagation of stress waves with a hybrid method coupling the NFE and a linear finite element (LFE). The subsection idea is developed based on previous works; calculating steps of the subsection method as well as techniques of passing motions from a source region to an elastic region are discussed. 3D numerical simulations of stress wave propagation in rock generated by decoupled explosion underground with two methods for comparison are carried out. The accuracy of the subsection method is demonstrated with simulated results. The demand of PC memory and the calculating time are investigated. The subsection method provides another approach for modeling and understanding the generation and propagation of explosion-induced stress waves, though, currently, studies are preliminary.
1. Introduction
The forward modeling of blasting stress waves has become an important approach to investigate the character of explosions underground as well as effects of blast and shock on infrastructures. The generation and propagation of blasting stress waves are quite complex issues. During the generation of stress waves, explosions cause irreversible nonlinear behavior in the surrounding material (rock). After the detonation of an explosive charge, a strong shock wave generates and propagates outward from the cavity. The strong waves cause inelastic deformations in the surrounding rock. Rock first flows plastically; then, it is stressed beyond its brittle failure limit and becomes granulated; then, it is stressed to the point where radial cracks grow but failure is not reached; then, it deforms nonlinearly where preexisting cracks slide but do not grow; and finally it deforms elastically [1]. Some nonlinear responses (e.g., material strength and damage) have strong impacts on the observed far-field seismic motions [2].
On the other hand, source region mediums and local structures have effects on the generation of blasting stress waves. For example, Bungum et al. [3] found that there are significant conversions of compressional energy to shear energy due to influences of the geologic heterogeneity and the local structures in the source region, after investigating waveform data of rock bursts and mining blasts and 3D finite-difference calculations on mining blasts.
Analytical methods of wave propagation, such as ray theory [4] and its modifications (e.g., the Kirchhoff-Helmholtz methods [5], the generalized ray method [6], and the Born approximation method [7]), mode matching methods, and the reflectivity method [8], can get theoretical results of the wave propagation under the point-source approximation of an explosion; however, these methods are only applicable to simple problems.
With the further investigation of the wave propagation and the increasing computational resources, many different numerical methods have been developed, such as, for example, the finite-difference methods (Geller and Takeuchi [8]; Kristek and Moczo [9]), the Fourier pseudospectral methods (Carcione [10]; Fornberg [11]), boundary integral equation methods and boundary element methods (Bouchon and Sánchez-Sesma [12]), and spectral element methods (Komatitsch et al. [13]; Chaljub et al. [14]). These methods have individual advantages during modeling particular wave phenomena with sufficient accuracy, but they have difficulties in modeling the generation of blasting stress because of irreversible nonlinear behavior in the rock of the explosion source.
The NFE methods can simulate the generation and attenuation of the shock in the rock surrounding the cavity efficiently taking irreversible nonlinear behaviors of the rock into account. For example, Rougier et al. [15] modeled the North Korean nuclear explosion on 25 May 2009. Recently, Xu et al. [16] preformed simulations on buried explosions and showed that the elastic source spectra simulated had been consistent with that predicted by the empirical source model of Mueller and Murphy [17] and Stevens and Day [18].
However, the NFE methods require extremely fine temporal and spatial discretization to model the nonlinear behavior under shock waves. Therefore, the simulating scale is only limited to several hundred meters from the blasting center due to the limitation of PC memory and acceptable time, which makes it computationally challenging to calculate the response on the scales required to compare calculations with observed ground-motion data recorded at tens of kilometers away from a source, which is very important to investigate explosion source characters and evaluate the effects of blast and shock on infrastructures.
Hybrid methods which couple two or more modelling methods provide feasible approaches to tackle this problem. For example, Moczo et al. [20] used a coupling scheme with a discrete wave number method for the source and path to a finite element method in order to incorporate the viscoelastic effects of the valley filled with a low-velocity sediment. Ma et al. [21] combined the finite element and staggered-grid finite-difference methods to simulate elastic P-SV wave propagation with complicated boundary conditions. De Martin et al. [22] simulated the wave propagation sequentially from the finite-difference method to finite element method domains by transferring the wave field information along the boundary. Ducellier and Aochi [23] combined finite element and 4th-order finite-difference techniques to model SH and P-SV seismic wave propagation in a 2D elastic medium with irregular surface topography. Wen [24] developed a hybrid method for SH wave field simulation including 2D localized heterogeneous structures by combining generalized ray theory and finite-difference method and studied shear velocity structures in the lowermost mantle beneath the central Pacific and South Atlantic Oceans. Wang et al. [25, 26] modelled the seismic wave propagation for P-SV wave field in a whole Moon model with a hybrid method based on pseudospectral method and finite-difference method and presented a parallel hybrid algorithm for two-dimensional global SH wave field simulation. Stevens and O’Brien [27] modelled the generation of blasting stress waves in a source region with a three-dimensional nonlinear finite element to get motions in a monitoring surface and then used the representation theorem to propagate motions to an elastic region. These hybrid methods involve either using analytical solutions in simplified layered structures or combining two linear numerical methods with various computational advantages; consequently, the detailed nonlinear responses of rock generated by shock waves are taken into account hardly. Recently, Rodgers et al. [28] and Xu et al. [16] developed a hybrid modeling approach with a hydrodynamic-to-elastic coupling in three dimensions. Near-source motions are computed with Eulerian hydrodynamics code (GEODYN) with adaptive mesh refinement and seismic waves are modeled with an elastic finite-difference code (WPP). The hybrid scheme provides a useful pathway to investigate the full-scale wave propagation from the source to the receivers. However, complete understanding of the wave-propagation spectrum from the near field to far field is a challenging issue.
In this paper, we presented a subsection method which combines the nonlinear finite element (NFE) method which applies to the rock in a source region and the linear finite elements (LFE) method which applies to an elastic region to simulate the stress wave. In order to verify this method, we carry out 3D numerical simulations of stress wave generated by the subsection method and a one-way method for comparison. The accuracy of the subsection method and the demand of PC memory and the calculating time are investigated.
2. Subsection Method
2.1. Subsection Idea
Yan and Zeng [29] presented the subsection idea that the whole region was divided into two regions, an inelastic region and an elastic region, and different methods were applied in different regions. However, the dimension of inelastic region determined by empirical formulas or blasting experiments [30] was approximate, which caused some problems. For example, the rock response in the elastic region was elastic in simulations but in fact not all the rock response in the elastic region was elastic. Therefore, in present study, we divide the whole region into the source region and the elastic region based on the rock dynamic responses and characters of the NFE and the LFE methods. The source region includes the entire nonlinear response region of rock and also includes part of the elastic response region. The NFE method is applied in the resource region, while the LFE method is applied in the elastic region.
Figure 1 shows the partition of regions. The region in the circle denotes the source region. The region between the rectangle and the circle means the elastic region, the star in the center denotes explosive source, and the circle denotes the monitoring surface between the two regions. The radius of the circle is larger than the elastic radius. We pass the motions from source region to elastic region via displacements of nodes on the monitoring surface. In simulations, the source region can also be set as rectangle which contains the entire nonlinear region to pass data conveniently.
2.2. Calculating Steps of Subsection Method
We carry out a subsection simulation with the following steps.(1)The elastic radius is calculated by empirical formulas or blasting experiments.(2)3D finite element models of source region are set up.(3)Node displacements on the monitoring surface are calculated with the NFE method and saved.(4)If necessary, the monitoring surface data are resampled, as permitted by the required resolution.(5)3D finite element models of elastic region are set up.(6)Node displacements calculated in Step are passed to the corresponding nodes of the elastic region.(7)The propagation of stress waves in the elastic region is modeled with the LFE method.
3D source region finite element models in Step can take detailed structures of the source region, such as the topography, the geometry of cavity, and the geometry and placement of charge, into account.
In Step , we use LS_DYNA 971 to model shock-wave motions. Some thermodynamically consistent nonlinear constitutive relations in LS_DYNA 971, such as the MAT_111 [31], and the user defined material model [32] can simulate the cavity formation, yielding, porous compaction, tensile failure, bulking, and the damage of rock during the propagation and attenuation of shock waves.
Passing motions from the source region to the elastic region with displacements in Step is consistent with the transformation strain theorem presented by Eshelby [33] and Aki and Richards [34]. Figure 2 illustrates the generation process of the stress waves according to the transformation strain theorem. The red region in Figure 2 denotes an explosion source, the green region denotes the surrounding material (rock), and the blue circle denotes the interface of the explosion source and the surrounding material. Firstly, the explosion source is incised from the whole region at time , and then a traction is put on the outside surface of the explosion source to keep it having the original shape at time . Secondly, another traction is also put on the outside surface of the explosion source to make it have the same transformation strain at time which is not incised from the whole region. Thirdly, we put the explosion source undergoing the transformation strain back to the surrounding region at time and remove the second traction. Finally, the stress waves in the elastic region from time to time are generated by the transformation strain. The displacements in Step from time to are the integral of the transformation strain at the normal to the monitoring surface.
3. Validation of Subsection Method
In order to investigate the fidelity and reliability of the subsection method, the generation and the propagation of the stress waves of explosion underground are modeled with both the subsection method and the one-way method. In present study, “the one-way method” refers to modelling the generation and propagation of the stress waves beginning with the explosion in the whole region mentioned above continuously only with the NFE method.
3.1. Finite Element Model
A decoupled explosion in rock is simulated with two methods. The whole region is rectangular with the dimension of . One ton of explosive charge is set in a ellipse charge cave. The radius of the charge is 0.527 m. The center of the explosive charge sphere is offset 1 m from the center of the ellipse.
Figure 3 shows 1/2 of the 3D element model of the subsection method. All of materials are modeled by solid elements. The whole simulating region is divided into the source region (Figure 3(a)) and the elastic region (Figure 3(c)). Figure 3(b) is the zoom of the explosive source element model in the middle of Figure 3(a). The dimension of the source region is and that of the elastic region is . 530800 elements and 547600 elements are used, respectively. The elastic radius of rock is taken as 20 m based on empirical formulas conservatively; so all of the inelastic responses are included in the source region. In the subsection model, we first simulate the generation of shock waves in the source region and then the propagation of stress waves in the elastic region. In Figure 3(a), the size of mesh transits from fine to coarse in radius directions. The fine element size of rock is the order of millimeter magnitude to capture the front of shock waves. In simulations, the NFE method is applied. The source region consisted of the rock which is modelled by Lagrange mess and the air and high explosive which is modelled by the Euler mesh. The interaction between the Lagrange mess and Euler mesh can be used to simulate the nonlinear response. In Lagrange mess, the coordinates move with the material, and, in the Euler mesh, the grid is fixed and materials are allowed to flow through. The Lagrange mesh imposes a geometric constraint on the Euler mesh whereas the Euler mesh provides a pressure boundary to the Lagrange mesh. The boundary conditions of the Euler mesh and the Lagrange mesh are set as an outflow boundary. In the elastic region in Figure 3(c), the rock of the elastic region is only modelled by the Lagrange mesh and the LFE method applied.
(a) The element model of the source region
(b) The element model of the explosive source
(c) The element model of the elastic region
Figure 4 shows 1/2 of 3D element models of the one-way method. The source region and the elastic region are connected together. In simulations, 1078400 elements are used in the whole region. In the one-way method we begin with the detonation of the explosive charge and then simulate the propagation of stress waves continuously in the whole region.
3.2. Material Model and Parameters
In the numerical model, the rock in the source region is modeled by the plastic kinematic material model [31]. On the other hand, the rock in the elastic region is modeled by an isotropic elastic material model. The plastic kinematic material model is expressed as follows:where is the initial yield stress; is the yield stress; is the hardening parameter; is the strain rate; and are Cowper-Symonds parameters; is the effective plastic strain; is the hardening modulus; and is the tangent modulus.
Air is modeled by an ideal gas. The pressure is expressed aswhere is a constant, is air density, and is the specific energy.
High explosives, TNT, are modeled using the Jones-Wilkins-Lee EOS. The pressure can be written aswhere is the hydrostatic pressure; is the specific volume; is the specific internal energy; and , and are material constants.
Material models and corresponding parameters in both methods are listed in Table 1. For the convenience of comparison following, we also list regions, element number, methods, memory demand, and calculating time in Table 1.
3.3. Comparisons between Two Methods
Simulated pressure radiation patterns in 1/2 of the elastic region are shown in Figure 5. Figures 5(a), 5(c), and 5(e) are pressure patterns simulated with the one-way method at s, 0.14 s, and 0.18 s after detonation, respectively, while Figures 5(b), 5(d), and 5(f) are pressure patterns simulated with the subsection method at s, 0.14 s, and 0.18 s, respectively.
(a) Radiation patterns with one-way method at s
(b) Radiation patterns with subsection method at s
(c) Radiation patterns with one-way method at s
(d) Radiation patterns with subsection method at s
(e) Radiation patterns with one-way method at s
(f) Radiation patterns with subsection method at s
Comparing the gradient of pressure radiation patterns in Figure 5(a) with that in Figure 5(b), the gradient in Figure 5(c) with that in Figure 5(d), and the gradient in Figure 5(e) with that in Figure 5(f), we find that simulated pressure radiation patterns with two methods have consistence at three time points, respectively. Elastic waves in the elastic region propagate continuously and each front wave at any time point is a wave source of the next time point; consequently, the consistence in Figure 5 indicates that the pressure radiation patterns simulated with two different methods have consistency in the elastic region. However, there are small differences between two sets of figures. The impossible reason is that the time step in two methods is different. The time step in simulations is determined by the minimum characteristic length of elements. For the demand of the shock-wave simulation, the magnitude of the minimum characteristic length of the whole region in the one-way method is a millimeter, while that of the elastic region in the subsection method is a meter considering the demand of the elastic-wave simulation. Therefore, the magnitude of the time step corresponding is a microsecond and a millisecond, respectively. The differences of three orders of magnitude will cause the differences in accumulating calculating errors in pressures, which need further investigation in our following study.
Figure 6 shows the simulated velocity of particle motion curve of node 537321 which is 800 m away from the center of the ellipse with two methods. Figures 6(a), 6(b), and 6(c) are those in the , , and directions, respectively. The solid line in each figure denotes the simulated velocity-time curve with the one-way method and the dot dash line denotes the simulated velocity-time curve with the subsection method. The rising and falling time, the maximum and the minimum, and the width of the impulses of solid line in each figure have good consistency with those of the corresponding dot dash lines. The slight difference between the solid line and the dot line in each figure is resulted in the difference of the time step in two methods mentioned above possible. The good consistency of simulated wave forms shows that we can get the same simulated wave forms with both methods.
(a) The velocities of node 537321
(b) The velocities of node 537321
(c) The velocities of node 537321
In the one-way method, we simulate the generation and propagation of stress wave continuously in whole region including both the source region and the elastic region using the NFE method, whereas in the subsection method we first model the generation of stress waves in the source region using the NFE method to get the displacement time of nodes in the monitoring surface, then pass the displacement to the corresponding node in the monitoring surface of elastic region, and then model the propagation of stress waves of the elastic region using the LFE method. The consistency of both pressure radiation patterns and velocity of particle motion in the elastic region with different methods indicates that the subsection method has the same reliability as the one-way method during the simulating propagation of stress waves in the elastic region generated by explosion.
In simulations, we must calculate stress, strain, strain rate, and so forth of 1,078,400 elements at each time step in the one-way method, while in the subsection method we only first calculate those of 530800 elements in the source region and then calculate those of 547600 elements in the elastic region. Though we implement two calculations in the subsection method, the memory demand is much lower than that in one-way method (8 GB in the one-way method and 4 GB in the subsection method, listed in Table 1), almost one second. The subsection method can reduce the demand for the memory which is limited in a PC.
The calculating time with both methods in WinXP64 operating system and 4 Inter(R) i7-2600 (3.4 GHz) processors PC is also listed in Table 1. The total time cost in the one-way method is 14 hr. 07 min. 46 sec., whereas that in the subsection method is 9 hr. 56 min. 50 sec. The comparison of the total time is 70.4%; that is, almost three-tenths of time is saved.
The reason why the subsection method reduces the calculating time and the memory demand is that we need to calculate total elements of the whole region at each time with a small time step in the one-way method, whereas we only calculate elements of the source region with the same time step firstly and then elements of the elastic region with a large time step in the subsection method.
4. Conclusion
This paper presents a subsection forward method to simulate the generation and propagation of stress waves induced by explosions underground. Firstly, we make an improvement on the subsection idea. We divide the whole region into a source region and an elastic region. The NFE method is applied in the source region, while the LFE method is applied in the elastic region. Motions are passed from the source region to the elastic region based on the continuity of displacements in the monitoring surface between the two regions. Secondly, we put forward 6 calculating steps to perform the subsection modelling and make the explanation on the propagation of motions between the two regions with the transformation strain theorem. Thirdly, we carry out 3D numerical simulations of stress wave by the subsection method and the one-way method for comparison. Simulated results show that pressure patterns and node velocities by the subsection method have good consistency with those by the one-way method. Simulated results also indicate that the PC memory demand and calculating time by the subsection method are decreased obviously compared with those by the one-way method because total elements in the whole region are calculated at each time with a small time step in the one-way method, whereas we only calculate elements of the source region with the same time step firstly and then elements of the elastic region with a large time step in the subsection method. These studies are preliminary; however, results are promising. Further efforts will be made to investigate accumulating calculating errors with different time steps.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This research is supported by the National Natural Science Foundation of China (no. 41374005 and no. 11402297). This paper is published with the approval of the National University of Defense Technology of China.