1 Introduction

With the development of precision-guided weapons,the probability of warship attacked by underwater contact explosion increases gradually. The physical process of underwater contact explosion is quite complex^{[1, 2, 3, 4]},involving steam-water mixture of high temperature and high pressure,water jet of high velocity,structural responses of nonlinearity,etc. The structures may be regarded as “fluid” in such high temperature and high pressure environment^{[3, 4]}. Therefore,a slight change in the energy radiation may directly affect the structural damage mode. The charge,as the energy source of underwater contact explosion,takes great effects on the energy radiation. Therefore,the studies on the charge parameters of underwater contact explosion such as charge shapes and detonation modes are of great significance. As more and more scholars pay attention to underwater contact explosion,many numerical algorithms^{[4, 5, 6, 7, 8, 9, 10]} have been proposed. However,some of them are still confined to empirical formulas and related standards when dealing with contact explosion loads,where the applied loads on structures are only drawn from a simple relation of the charge weight and the specified distance. While the process of charge detonation and the interactions with its surrounding media are seldom considered. Based on the axisymmetric smoothed particle hydrodynamics (SPH) method^{[11, 12]},an underwater explosion model is established in this paper to study the influence of charge parameters on near-field loads.

The SPH method featured by Lagrangian particles^{[5, 13]} has great advantages for solving nonlinear problems of large-deformation,splash,multiphase mixing,moving interface,etc. It can simulate the real process of charge detonation,and is more accurate and convictive when dealing with underwater contact explosion^[13]. However,the simulations of underwater contact explosion with the SPH method are mostly restricted to two dimensional models,and the researches based on three dimensional models are rarely published due to the expensive program. In this paper,some improvements are made based on Omang’s axisymmetric SPH method^[11] with the given symmetry of charge. Due to the large gradient of density in underwater contact explosion,the density approximation suitable for multiphase problems is introduced. For the presence of pressure fluctuation and poor accuracy near the symmetric axis^{[12, 14]},a new approach of dynamic boundary particles is proposed to deal with the symmetric axis. Then,the influences of charge parameters on near-field loads are investigated,and a simple model of underwater contact explosion is established to verify the above analyses. 2 Theoretical background

The present axisymmetric SPH method is established based on Omang’s model^[11]. In this section,the cores of the axisymmetric SPH will be detailed in terms of the kernel function,the governing equations,etc. The constitutive relation in the axisymmetric coordinates will be also introduced. Then,a new approach for dealing with the symmetric axis is proposed. 2.1 Axisymmetric kernel function

In cylindrical coordinates,the integral of function f(r) is defined as follows^{[11, 12]}:

where W_(q) is the B-spline kernel function^[13]. Define the function as the axisymmetric kernel function. Then,we have When 0 ≤ q≤ 2,W_(q) ≠0. Therefore,the definition domain of φ is reduced to be where Meanwhile,the function W_(q) is an even function. Therefore,the axisymmetric kernel function becomes^{[11, 12]} Therefore,the axisymmetric function is Similarly,the approximations of derivatives ▽f(r,z) and ▽· f(r,z can be obtained as follows: The partial derivatives of W_c versus r,z,r′,and z′ are obtained as follows^{[11, 12]}: It is obvious that the partial derivatives of the axisymmetric kernel function W_c are asymmetric versus r and r′ but symmetric versus z and z′^[11],i.e., 2.2 Axisymmetric governing equation

The process of underwater explosion is transient accompanied with strong shock. Thus,the fluid can be regarded as inviscid,and the whole process is adiabatic. Due to the presence of gas,liquid,and solid phases,the continuity equation suitable for multiphase flow is applied here. The discretized equations are presented as follows:

where Eqs. (12)-(15) denote the continuity equation,the momentum equation,and the energy equation,respectively. The symbols ρ,υ,e,σ,p,S,and ε indicate the density,the velocity,the internal energy,the stress tensor,the isotropic pressure,the deviatoric stress,and the strain tensor,respectively. α and β represent the direction of coordinates. Ⅱ_ij is the artificial viscosity^[13]. The last term in the energy equation only appears in solid material. 2.3 Axisymmetric constitutive relation

The resultant stress tensor σ^αβ of solid material is composed of two parts^[13]: the isotropic pressure p and the deviatoric stress S,which is expressed as follows:

where the deviatoric stress S can be obtained from the integration of the deviatoric stress rate. The expressions of is where G,^αβ,and denote the shear modulus,the strain rate tensor,and the rotation rate tensor,respectively. The expressions of ^αβ and are as follows: Therefore,the discretized forms in the axisymmetric SPH method are The Johnson-Cook yield model^[15] is introduced to evaluate the dynamic yield stress. The von Mises yield criterion^[16] is applied to judge whether the solid material yields or not. If the second-order invariant of the deviatoric stress J₂ exceeds σ²_y/3,the modified deviatoric stress of the particles can be treated as follows: 2.4 Symmetric axis treatment

To prevent the particles penetrating the symmetric axis non-physically,Zhang et al.^[12] proposed the mirror particle method,which showed the effectiveness of preventing particles from penetrating the axis. However,it may cause serious pressure fluctuation near the symmetric axis. In theory,the axisymmetric SPH method should show perfect results near the symmetric axis even without any treatment. However,since the symmetric axis is always positioned at the boundary,with the accumulation of approximation errors,the penetration always comes into being,especially in the problems of strong shock. Therefore,it is not beyond the imagination that a number of boundary particles are distributed near the symmetric axis to prevent the penetration. A new treatment for dynamic boundary particles is proposed in the present paper. The dynamic boundary particle does not participate in the density approximation of real SPH particles,but provides artificial resistance to the real SPH particles in order to prevent them from penetrating the symmetric axis. The results prove that the instability problem of pressure near the symmetric axis can be overcome easily.

The properties of dynamic boundary particles are

where r'_i,z'_i,υ'_ri,υ'_zi,υ'_i,m'_i,p'_i,e'_i,and h'_i denote the parameters of dynamic boundary particles,while r_i,z_i,υ_ri,υ_zi,υ_i,m_i,p_i,e_i,and h_i denote the parameters of real SPH particles near the symmetric axis. When r_ii' 6 h_i is satisfied between the real SPH particle i and the mirror particle i',it is considered that the SPH particle tends to penetrate the symmetric axis. At the moment,an artificial resistance is applied to the real SPH particle by the dynamic boundary particle as follows: where β is an adjustable parameter,and c is the sound velocity.

The underwater explosion model of a spherical charge is established to test the effects of the present treatments for the symmetric axis. The mirror particle and the dynamic boundary particle are positioned near the symmetric axis in the model,respectively. The numerical results are shown in Fig. 1. It is obvious that the pressure around the detonation gas in the proposed method seems to be more uniform,and the fluctuation mentioned in Refs.^{[12, 14]} are overcome well. Moreover,the shape of the explosion bubble is closer to a sphere. In order to show the differences more clearly,the relative peak pressure is defined as p = p/p_∑-1,where p denotes the peak pressure measured in a given point and p_∑ is the average value of the peak pressure of the given radius.

The result of p varying with the angle θ in a distance of R = 0.5m apart from the charge center is shown in Fig. 2. The lower end of the symmetric axis is regarded as θ = 0◦. It is not hard to draw that the relative pressure of the present method fluctuates around zero,and the biggest deviation is within 2%; while the relative pressure of the mirror particle method is serious,the biggest deviation is up to 16% near the symmetric axis and −10% at θ = 90◦,and the explosion bubble shape is not like a sphere. In conclusion,the treatment of the symmetric axis in the axisymmetric SPH model is a key problem,which not only affects the numerical stability but also the numerical accuracy. Last but not least,the axisymmetric SPH model established in this paper has greatly improved the accuracy with the dynamic boundary particle at the symmetric axis and has met the high-accuracy requirement of the present study,which is also verified in the following numerical verification.

3 Results and discussion

Two kinds of underwater explosion models,i.e.,spherical charge and cylindrical charge,are established to study the influences of charge parameters on near-field loads,including the peak pressure and impulse. The spherical and cylindrical TNT charges of the same mass are detonated in spherical water. The spherical charge has a radius of r = 0.085 m,and is discretized into 459 particles; while the cylindrical charge has a length-diameter ratio l: 2r = 1:1,1:1.5,1:2 on the premise of the same mass,and the length l and the radius r are variable. The particle spacing of the above models is △x = 0.005 m.

The underwater explosion loads are evaluated in a certain distance R and a certain range of θ,as shown in Fig. 3,θ ∈ [0◦,180◦]. The measure points are positioned at the interval of 10◦.

The spherical charge model is selected as the reference condition. p₀,I₀ and p,I indicate the peak pressure and impulse in the reference condition and the investigation condition,respectively. Thus,the relative pressure p = p/p₀-1 and the relative impulse I = I/I₀-1 are defined. The JWL state equation^[17] is used for TNT detonation products,and the Mie-Grüneisen state equation^[18] is used for water. The propagation process of shock wave is shown in Fig. 4.

3.1 Validation of axisymmetric SPH model

In order to validate the numerical accuracy of the present axisymmetric SPH method,the underwater explosion model of free field is established to compare with experimental results. The experiment is carried out at China Academy of Engineering Physics. A spherical charge of 8 kg,as shown in Fig. 5,is detonated in the explosion tank. The same model is established with 252 869 SPH particles including 104 TNT particles,the diameter of the spherical water field is d = 10 m. The pressure of shockwave varying with time at a distance of 4 m apart from the charge center is shown in Fig. 6,which indicates that the result of the axisymmetric SPH agrees well with that of the experiment (noted as “EXP”). There is a pity that the pressure of shock wave takes a longer time to climb and dip. The reason may come from the viscosity and energy dissipation of real fluid.

3.2 Influence of charge shapes

It is of great significance to study the influence of charge shapes on near-field loads. On the premise of a same charge mass,the cylindrical charge models of different length-diameter ratios

are established. The peak pressure and impulse varying with the distance and the angle are investigated; the point-source detonation is adopted in all cases to avoid the effect of detonation modes. The pressure nephogram of the cylindrical charge with the length-diameter ratio l : 2r = 2:1 is shown in Fig. 7. It is obvious that the circumferential pressure is nonuniform at the same radius,and at the detonation end,the pressure is weakened but enhanced at the other side.

The dimensionless radius R/r is defined here to evaluate the distance between the measure point and the charge center,where r indicates the radius of the spherical charge,and R is the distance of a measure point apart from the charge center. The relative peak pressure p and the relative impulse I versus R/r are shown in Figs. 8 and 9,respectively,the arrow length represents the amplitude and the direction indicates a positive or negative value. The positive value means pointing outward while the negative value indicates pointing inward.

In Figs. 8(a)-8(c),as for R/r = 6.47,when the length-diameter ratio is l : 2r = 1 : 1,the enhanced or weakened effects of the peak pressure are not obvious; as the ratio increases,the influence of the charge shape increases; when l : 2r = 2:1,the peak pressure is enhanced at the interval 0° < θ < 125° with the maximum value up to 20%,especially in the radial direction,while is weakened seriously at the interval 125° < θ < 180° and the amplitude in the normal direction of the detonation end is reduced up to 25% maximally. However,with the increase in the distance R/r,the influence reduces,as shown in Figs. 8(c),8(f),and 8(i). In Fig. 9,when the length-diameter ratio increases from l: 2r= 1:1 to l: 2r = 2 : 1,the impulse I is uniform.

To sum up,the peak pressure of shock wave is greatly affected by the charge shapes in the near field,while the impulse is less affected. In a certain range,the bigger the ratio l : 2r is,the more obvious the directivity of the peak pressure will be. Moreover,in underwater contact explosion,the structural damage is greatly affected by the peak pressure. Therefore,the effect of the charge shape in the underwater explosion of the near field should not be overlooked. 3.3 Influence of detonation modes

The energy radiation of the charge is significantly affected by detonation modes. In general,the cylindrical charge is detonated at a point of one end. However,in some engineering practice,a reasonable detonation mode is seldom considered.

The cylindrical charge with the length-diameter ratio l : 2r = 2 : 1 is taken as an example to study the influence of detonation modes. The peak pressure of different distances is shown in Figs. 10(a)-10(c). Compared with the point-source detonation in Figs. 8(c),8(f),and 8(i),the law of the plane-source detonation around the charge is similar,but the amplitude in either the weakened or the enhanced direction is raised by about 5% around the charge. As the dimensionless radius R/r increases,the differences decrease gradually.

The impulse of different distances around the cylindrical charge of l : 2r = 2 : 1 is shown in Figs. 10(d)-10(f). Compared with the point-source detonation in Figs. 9(c),9(f),and 9(i),the impulse of the plane-source detonation increases up to about 5% when 90° < θ < 180°,while the amplitude is similar at 0° < θ < 90°. Therefore,the plane-source detonation can increase the energy of the shock wave radiating into water,i.e.,compared with the point-source detonation,the energy of the plane-source detonation is easier to radiate into water.

The shock wave pressure varying with the time of the cylindrical charge with l : 2r = 2 : 1 is shown in Fig. 11. The measure points are positioned at R/r = 6.47. In the left figure,the point-source detonation is adopted,while the plane-source detonation is adopted in the right one. It can be drawn from the curves why the impulse of shock wave is uniform. Though the peak pressure at θ = 0° is much bigger than that of θ = 180°,the pressure at θ = 0° decreases fast. Therefore,the impulse around the charge is similar.

In conclusions,the near-field load of the plane-source explosion has a greater directivity. Compared with the point-source explosion,the peak pressure of the plane-source detonation can increase up to about 5% in the near field and at almost all angles. Besides,the impulse can increases about 5% only in the angles near the detonation plane,while is similar to the point-source detonation at other angles. Therefore,the detonation modes of the near field underwater explosion deserve attentions. 3.4 Numerical verification with underwater contact explosion model

In order to validate the above conclusions,an axisymmetric model of underwater contact explosion is established. The influences of the charge parameters on underwater contact explosion are verified through the deflection of the steel plate,the velocity of the fragment,etc.

The damage of the steel plate in contact explosion generally goes through three stages: dishing,discing,and petalling^{[19, 20, 21]}. The simulation of underwater contact explosion with the axisymmetric SPH method in the present research is only limited to the dishing and discing stages. The numerical model,a circular steel plate with a radius of 0.8m and a thickness of 0.024m is subjected to the shock load of the cylindrical TNT charged 1.26 kg with the lengthdiameter ratios l : 2r = 2 : 1 and l : 2r = 1 : 1 in a hemispherical water field,respectively. The Mie-Grüneisen equation in solid mechanics^[22] and the Jackson-Cook constitutive relation^[15] are applied for steel,the detailed parameters are shown in Tables 1 and 2. In order to validate the influences of the charge parameters (charge shapes and detonation modes),four typical cases are arranged here: the point-source detonation at the upper end with l : 2r = 2 : 1,the point-source detonation at the lower end with l : 2r = 2 : 1,the plane detonation at the lower end l : 2r = 2 : 1,and the plane detonation at the lower end with l : 2r = 1 : 1. If the equivalent plastic strain of a steel particle reaches 0.3,the particle turns to be failure,and it would never be considered in the numerical model.

The shock wave pressure profiles in water and Mises stress in the steel plate at t = 0.15ms are shown in Fig. 12. It is obvious that the propagation of the shock wave and the response of the steel plate comply with the physical law well. The round crevasse and the deflection of the steel plate at t = 0.24ms are shown in Fig. 13. Considering the influence of charge shapes,the explosion load around the charge has a great directivity. Therefore,the responses of the circular steel plate in weakened and enhanced zones are different,which are illustrated in Figs. 12,13(a),and 13(b). It can be seen that the steel plate is positioned in the weakened zone of peak pressure for the case of point-source detonation at the upper end while is in the enhanced zone for the case of point-source detonation at the lower end. The energy radiated to the water in Fig. 13(a) is larger than that in Fig. 13(b). Therefore,the deflection in Fig. 13(a) is smaller than that in Fig. 13(b),which proves the above conclusions in Section 3.2. Considering the influence of detonation modes,as shown in Figs. 12,13(b),and 13(c),the failure of the plate with the plane-source detonation is more serious,the damage radius is larger but the deflection along the radial direction is slightly smaller,which indicates the stronger directivity of the plane-source detonation. Considering Figs. 13(c) and 13(d),because of the directivity of energy radiation,the damage radius and the deflection of the charge shaped with l : 2r = 1 : 1 are slightly larger,which may be well explained from the peak pressure and the impulse from Figs. 8,9(a),and 9(c). Besides,the energy absorbed by the plate in Fig. 13(d) is larger than that in Fig. 13(c) due to the smaller length-diameter ratio of Fig. 13(d). To sum up,with the same charge mass,the radius of discing in underwater contact explosion is related to charge shapes and detonation modes.

In order to show the influences of charge parameters further,several curves are presented. The deflection curves varying with time at the point of R = 0.116m apart from the plate’s center are shown in Fig. 14. The slope of the curves shows the deformation velocity. As for l : 2r = 2 : 1,it is obvious that the slop,which indicates the extents of structural deformation,is the smallest for the case of the point-source detonation at the upper end and the largest for the case of the plane-source detonation at the lower end. As for l : 2r = 1 : 1,the deflection is larger than that for l : 2r = 2 : 1 (see Figs. 13(c) and 13(d)). Though the impulses in Figs. 13(c) and 13(d) (shown in Figs. 8 and 9) in the direction of the steel plate seem to be similar,the contact area of TNT and the plate in Fig. 13(d) is larger and the time of the shock wave impact on the plate is longer,which will cause larger deflection for the case of plane-source detonation at the lower end with l : 2r = 1 : 1. At the same time,the larger detonation plane for the case of plane-source detonation at the lower end with l : 2r = 1 : 1 also promotes the result. The displacement curves of the fragment from the center of the steel plate are shown in Fig. 15. Similarly,the plane-source detonation drives the fragment faster,and the charges of different length-diameter ratios also show different results. In conclusion,the detonation modes and charge shapes have great influences on the structural damage.

4 Conclusions

The axisymmetric SPH method is introduced to deal with near-field underwater explosion,in which the density approximation equation suitable for multiphase flow and the constitutive relation of solid material based on axisymmetric coordinates are adopted. Through the improvement of the symmetric axis,the influences of the charge parameters on near-field shock loads in underwater explosion are investigated. Finally,an underwater contact explosion model is established to verify the above analyses,and the conclusions are as follows:

(i) The dynamic boundary particle can improve the numerical stability and accuracy near symmetric axis.

(ii) The influence of charge shapes on near-field loads in underwater contact explosion should not be overlooked. Compared with spherical charge,the peak pressure of cylindrical charge with different length-diameter ratios shows a strong directivity in the near field. On the premise of a same charge,in a certain range,the larger the length-diameter ratio is,the greater the peak pressure along the axial direction is weakened and the peak pressure along the radial direction is enhanced,which may be up to 25% and 20%,respectively. The impulse around the charge is uniform,which is less affected by charge shapes.

(iii) The influence of detonation modes on both peak pressure and impulse of the near field is also considerable. Compared with the point-source detonation,the peak pressure of the plane-source detonation around the charge may be raised by 5%. Besides,the impulse near the detonation end may be raised about 5%,but similar to the point-source detonation at other angles.

(iv) The differences of structural damage parameters,such as the radius of crevasse,the deflection,and the velocity of fragments,have proved that the influences of charge parameters should not be ignored in underwater contact explosion. The improved axisymmetric SPH method can be applied to other dynamic simulations.