《1 Self-simulation solution of a strong TNT charge explosion》

1 Self-simulation solution of a strong TNT charge explosion

The point-source explosion theory is one of simplest models for studying of the shock wave of a spherical charge explosion. It is assumed that the mass of detonation product is infinitesimal, the energy released by the charging explosion is limited, and the influence of static air in motion is ignored. Therefore, there is no single parameter that can represent the initial linear scale of the explosion. The distribution of any physical quantity in space is self-similar with time. That is, the gas behind the shock wave close to the center of the explosion satisfies the self-similar condition [1]. Taylor [2], Von Neumann [3], Bach and Lee [4], and Sedov [5] independently obtained the theoretical solutions of the point source strong explosion self-simulation. Compared with the experimental results, the flow field characteristics of the initial stage of an explosion can be well-described [6]. This study, on the basis of the point-source explosion theory, analyzed the self-simulation flow field of the one-dimensional symmetric surface (surface source, line source, and point source) of the plane, cylinder, and sphere.

《1.1 Governing equation》

1.1 Governing equation

The propagation of a blast wave is a highly complicated and uncertain process, but it satisfies the conservation of mass, momentum, and energy. Because of the rapid change of motion considered in explosive mechanics, there is no obvious momentum and energy transfer between the streamlines, so viscosity and heat conduction can be neglected [7].

1.1.1 Governing equation in Euler coordinates

The state quantities of the flow field in Euler coordinates mainly include pressure P, density ρ, particle flow velocity u, internal energy E, and temperature T, all of which are functions of independent variables x, y, z, and t.

Conservation of mass:

where, and u, v, and w are the projection of velocity vector V on the x, y, and z coordinates, respectively.

Conservation of momentum:

where F is the mass force.

Conservation of energy:

Equation of state:

For the one-dimensional symmetric TNT charge explosion, the motion of the medium can be regarded as one-dimensional adiabatic motion, that is, the values in the medium change depending only on time and a single geometric coordinate. If the mass force is removed, namely, F = 0, the governing equations described by Euler coordinates in the Cartesian coordinate system, cylindrical coordinate system, and spherical coordinate system can be uniformly written as follows:

where r is the Euler position coordinate and t is time. N is a constant equal to 0, 1, or 2, corresponding to one-dimensional plane symmetry, one-dimensional cylindrical symmetry, and onedimensional spherical symmetry, respectively.

1.1.2 Governing equation in Lagrange coordinates

In order to improve the accuracy of the local calculation of a shock wave front, the equations can be written in Lagrangian form:

with an added motion relation:

where R is the Lagrangian position coordinate, r is the Euler position coordinate, u is the particle velocity, P is the pressure, E is the internal energy, ρ is the density, ρ₀ is the initial air density of the explosion wave front, t is time, and N is the corresponding exponent of the plane coordinate, cylinder coordinate, and sphere coordinate.

For the convenience of calculation, the mass coordinate equation is introduced:

and. Then, the governing equation in Lagrange coordinates can be transformed into dimensionless form as follows:

where N is a constant equal to 0, 1, or 2, corresponding to dimensionless equations of one-dimensional plane symmetry, one-dimensional cylindrical symmetry, and one-dimensional spherical symmetry, respectively.

《1.2 Equation of state for detonation products of TNT charge in real air》

1.2 Equation of state for detonation products of TNT charge in real air

1.2.1 The dimensionless form state equation in real air

The experimental results show that the ideal gas equation of state is no longer satisfied when the temperature is greater than 2000 K or the shock wave overpressure is greater than or equal to 4 MPa [7]. With increasing temperature and pressure, there will be chemical reactions between various molecules in the air, such as the dissociation of oxygen molecules and the formation of NO and NO2. The air component is different at different temperatures. When the temperature is higher than 8000 K, ionization reactions of various atoms in the air take place, and the air component becomes more complicated. When the temperature reaches 300 000 K, the ionization process of the gas can be completed at any density. The ionization, molecular dissociation, and excited molecular vibration of each atom in the air during the rising temperature require energy consumption. Therefore, at high temperatures, the equation of state of real air is significantly different from that of ideal gas [8,9]. From experiments, Russian scholar Kuznetzov used the statistical physics and theoretical method of shock wave gas dynamics to obtain the thermodynamic function table of air and the calculating table of shock hot lines. The fitting formula of the real air state equation is obtained by fitting the data results by Brode [10]:

where, in formulas (10)–(16), e, θ, π, and η are dimensionless internal energy, dimensionless temperature, dimensionless pressure, and dimensionless density, respectively.

1.2.2 State equation for TNT detonation products

Brode, a scholar who studied the equation of state of the detonation products of TNT charge in detail, found that the equation of state is given in the same form as the real air equation of state [10]:

where, in formula (19), α₁ = 1.76, b₁ = 52.4892, and c₁ = 7.3.

《1.3 Calculation of flow field parameters in the detonation stage of a TNT charge with surface source, line source, and point source》

1.3 Calculation of flow field parameters in the detonation stage of a TNT charge with surface source, line source, and point source

Using the theorems of Π in dimensional analysis, the independent self-simulation solution variables of the strong explosion of a TNT charge are obtained:, charge with a surface source;, charge with a line source; and, charge with a point source. Here, ρ₀ is the initial density, C₀ is the velocity of sound in the atmosphere, E is the energy released by the explosion, r is the spherical coordinates, and t is the time. If the position of the shock wave front is rs, then the value of ξ is ξ₀. Thus, the governing equations of the strong explosion self-simulation solution of the TNT symmetrical charge from the surface source, line source, and point source can be obtained as follows:

Charge of surface source:

Charge of line source:

Charge of point source:

Boundary conditions: f(1) = g(1) = h(1) = 1, at the shock wave front; f(0) = 0 at the blasting center. Here, f is the dimensionless velocity, g is the dimensionless density, and h is the dimensionless pressure.

Formulas (21)–(23) are closed equations, which can be solved. The self-simulation numerical solutions of the strong explosion of a surface source, line source, and point source charge can be obtained by calculation. The results are shown in Figs. 1–3.

《Fig. 1》

Fig. 1. Distribution of parameters of the strong explosion solution for a surface source charge.

《Fig. 2》

Fig. 2. Distribution of parameters of the strong explosion solution for a line source charge.

《Fig. 3》

Fig. 3. Distribution of parameters of the strong explosion solution for a point source charge.

Figs. 1–3 show that the characteristics of the distribution of the strong explosion parameters with changes in the relative positions of the shock wave front are essentially the same, and the variation trend is similar. The particle velocity decreases linearly with decreasing r/r_s. In the process of decreasing the pressure with the change in r/r_s from 1 to 0, the pressure first decreases rapidly. When r/r_s drops to 0.5 until the explosion center, the pressure remains essentially stable at 0.36–0.37. The density quickly decreases and approaches zero with the change in r/r_s.

The experimental results show that the self-simulated solution [7] can well describe the initial stage of an explosion. Therefore, it is accurate and reasonable to use the self-simulation solution of the strong explosion of a non-point source, line source, and point source under the real air condition as the initial conditions for the numerical calculation of the near-region flow field characteristics of a TNT charge explosion.

《2 Flow field characteristics in the near field of a TNT charge explosion shock wave》

2 Flow field characteristics in the near field of a TNT charge explosion shock wave

《2.1 Explosion flow field physical model》

2.1 Explosion flow field physical model

When the detonation wave propagates to the charge surface with the initial calculation time, the front plane of the detonation wave divides the whole flow field into two parts (Fig. 4).

According to the laws of conservation, the following relations must be satisfied for the shock wave front:

Conservation of mass:

Conservation of momentum:

Conservation of energy:

where, in formulas (24)–(26), P₀, u₀, ρ₀, and E₀ are, respectively, the pressure, particle velocity, density, and internal energy of air when it is not disturbed before the shock wave front. P, u, ρ, and E are, respectively, the pressure, particle velocity, density, and internal energy of air when it is disturbed after the shock wave front. D is the shock wave velocity.

After making the upper formula dimensionless, the shock relation corresponding to the single value can be solved using the difference method to solve the partial differential equations (9). In order to keep the calculations continuous at the shock wave contact discontinuity and jump discontinuity, the artificial viscosity method is adopted, which is used to add the artificial viscous term ζ when the pressure term exists. It is written in the dimensionless form ρ = ζ/η. Thus, equations 

《Fig. 4》

Fig. 4. Flow field distribution at the initial calculated time.

(9) (when N = 2) can be written in dimensionless form as follows:

where l is the integral grid of the differential method, n is the integral step length of time, e is the intersection point of the half grid of location l–1/2, and the integral grid of time n + 1.

The artificial viscous term is:

As to the difference form of the energy equation, because of the complexity of the parameters, a new form can be obtained by substituting e on the basis of the dimensionless state equation:

where ξ = lnη.

Thus, the difference form of the energy equation is obtained:

Because the upper expression is an implicit difference form, the quadratic iteration method is used to solve the equation. Establishing the assembly method at the back of the shock wave means first performing shock wave exploration by calculating the maximum value of ξ using the cubic spline interpolation method, and making the location of the maximum value of ξ as the location of the shock wave front λ_s. Thus, the pressure of the shock wave front is obtained. Then, the parameters at the shock wave front are solved with the solver subroutine. The time step of the numerical calculation is determined by the Courant condition and the stability condition of the difference equation in the shock region.

《2.2 Example analysis and result》

2.2 Example analysis and result

2.2.1 Air and TNT charge parameters

The parameters of stationary air are shown in Table 1.

《Table1》

Table1. Parameters of static air.

The parameters of the TNT charge are shown in Table 2.

《Table 2》

Table 2. Parameters of TNT charge.

2.2.2 Results

The following results can be obtained by using the calculation method established in this study:

(1) Overpressure

Fig. 5 shows the curves of the dimensionless peak overpressure ΔPs and the peak dynamic pressure qs with the position of the dimensionless shock wave front. It can be seen from the diagram that the calculated results in this study agree closely with the results of Beker [11].

《Fig. 5》

Fig. 5. Curves of peak overpressure and peak dynamic pressure with shock radius (the solid points are the result of this study, and the hollow points are the result of Beker).

(2) Time-space curves of contact discontinuous surface and main shock wave front

Fig. 6 shows the main shock wave（MS）, secondary shock wave （2S）, and the dimensionless position λ of the contact discontinuity surface change with dimensionless time τ. From Fig. 6, one can see the phenomenon of the secondary shock wave found by Friedman [12], and many subsequent shock waves caused by the reflection of the waves. The influence range of a shock wave is much larger than the diffusion range of the contact discontinuity surface of detonation products from Fig. 6.

《Fig. 6》

Fig. 6. Time-space curve of shock wave and contact discontinuity.

(3) Curve of overpressure peak and dynamic pressure peak varying with the shock wave front position

Figs. 7 and 8 respectively show the curves of the overpressure peak of the main shock wave and the change of the dynamic pressure peak of the main shock wave with the position of the shock wave front (using exponential coordinates) for a shock wave produced by the explosion of a spherical TNT charge in real air with different charges. From Figs. 7 and 8, it can be seen that for the same charge amount, the overpressure peak and dynamic pressure peak of spherical TNT charge decreases with increasing distance. The changes in explosion overpressure and dynamic pressure with different charge amounts conform to the similarity law. With increasing charge amount, the overpressure peak and dynamic pressure peak at the same position increase accordingly.

《Fig. 7》

Fig. 7. Curves of peak explosion overpressure of different spherical TNT charge amount varying with the position of the shock wave front.

《Fig. 8》

Fig. 8. Curves of the dynamic pressure peak of a spherical TNT charge with different charge amount and the position of the shock wave front.

(4) Curve of shock wave arrival time varying with the shock wave front position

Fig. 9 shows the effect of shock wave arrival time on the position of the shock wave front of a spherical TNT charge in real air. It can be seen that with increasing charge amount, the time of the blast wave of a spherical TNT charge to the same position decreases. The variation law of shock wave arrival time of different charge amount is approximately the same, so it can be seen that the similarity law relation is also satisfied.

《Fig. 9》

Fig. 9. Curves of shock wave arrival time of different spherical TNT charge amount varying with the position of the shock wave front.

(5) Curve of shock wave overpressure impulse varying with the shock wave front position

Fig. 10 shows the curve of the overpressure impulse of the explosion shock wave of a spherical TNT charge with different charge amount and the position of the shock wave front. From Fig. 10, one can see that for the same charge amount, the impulse of a spherical TNT charge decreases gradually with increasing position of the shock wave front. When a certain position is reached, there is a gentle or even slightly increased phase in this area. Then, attenuation begins because of the increase in barotropic impulse caused by the secondary shock catching up with the main shock wave. On the other hand, with increasing charge, the impulse at the same position increases gradually.

《Fig. 10》

Fig. 10. Curve of overpressure impulse of an explosion shock wave of different spherical TNT charge amount with the position of the shock wave front.

(6) Curve of positive pressure time of a shock wave with the shock wave front position

Fig. 11 shows the curve of the time of positive pressure acting on the shock wave of a spherical TNT charge with different charge amount and the position of the shock wave front. One can see from Fig. 11 that for the same charge amount, the variation in the shock wave with the distance of the shock wave front is similar to that in Fig. 10, and there is a turning point in the time of positive pressure action, which decreases with increasing distance. However, at a certain point it starts to increase with distance. This is mainly due to the existence of secondary shock waves. As the shock wave continues to propagate outward, the peak overpressure of the main shock wave decreases,and the longer the main shock wave propagates, the smaller the positive pressure action time will be. Because the secondary shock wave continues to catch up, when it arrives at a certain position, the positive pressure of the main shock wave is not reduced to negative pressure, and the secondary shock wave catches up, thus increasing the duration of the barotropic action. For a TNT charge of 1–50 kg, the positions of the turning point are not different, approximately 4 m; however, for a large-volume charge ball of 100–1000 kg, the distance from the turning point to the center of the charge increases with increasing charge amount. In addition, for a ball with 1–50 kg charge, the distance from the turning point to the center of detonation increases with increasing amount of charge, and at the fixed position point, the time of positive pressure action increases with increasing charge amount before and after the turning point; however, for a charge ball with 100–1000 kg charge, the time of positive pressure action also increases with increasing charge amount before the turning point. However, after the turning point, the positive pressure time decreases with increasing charge (Fig. 11).

《Fig. 11》

Fig. 11. Curves of the response time of shock wave positive pressure for different spherical TNT charge amount varying with the position of the shock wave front.

《3 Conclusions》

3 Conclusions

Because there is a great difference between propagation characteristics in in the near field of a TNT charge explosion shock wave in real air and that in ideal air, the real air state equation is used to describe the characteristics of the near region of the explosion flow field. From the theory established in this study, combined with the artificial viscosity technique and shock wave assembly technology, the parameters of the explosion flow field in real air of a TNT charge are calculated. The results agree well with the previous research results, and the following conclusions are obtained:

(1) After the explosion of a TNT charge, a strong shock wave with an initial overpressure of 60 MPa is produced, and the intensity decreases gradually with the continuous propagation of the shock wave; there is a contact discontinuity surface that propagates continuously after the main shock wave.

(2) The most typical physical characteristic of the near region of the explosion is the existence of a detonation product region. In the region with a radius less than three times the radius of the charge, the pressure intensity is far greater than the pressure value of the main shock wave.

(3) When the detonation wave reaches the interface between the charge and the air, it forms an inward shock wave as well as an outward-propagating air shock wave, and the inward propagating shock wave is strengthened after the collision at the center of the detonation. Then, it continues to propagate outward and forms a new shock wave that propagates inward and outward after reaching the contact surface and forms the second and third shock waves in the negative-pressure region of the fixed-point overpressure waveform.

Therefore, the method proposed in this paper can well simulate and describe the propagation of a TNT charge in in the near field of a TNT charge explosion shock wave. The research results of this paper can play an important role in the study of explosion mechanics and protection engineering.