#### Abstract

It is of great significance for the protective design of submarine to study the influences of coverings on the damage characteristics of single and double cylindrical shells subjected to underwater contact explosions. The SPH models of single and double cylindrical shells coated with foam silicone rubber are established to analyze shockwave propagation, damage characteristics, and elastoplastic responses, which provides reasonable parameters of covering position and thickness. The results can be concluded as follows: the superposition of multiple waves may cause the inhomogeneity and discontinuity; for the single cylindrical shell with inner or outer coverings, the damage mode is mainly tensile and shear failure is caused by detonation waves and detonation products; compared with out-covering approach, the in-covering approach has better antishock performance; the best protective effect comes out when the thickness of covering is close to that of the shell; as for the double cylindrical shell without interlayer water, the destruction of inner shell mainly results from the puncture of high-speed fragments from the outer shell, so for the outer shell, out-covering is a better choice; however, since the interlayer water is very effective in protecting the inner shell, in-covering will be better for the inner shell.

#### 1. Introduction

With the development of precision-guided weapons, the probability of naval structures being attacked/damaged by underwater contact explosion increases gradually. Therefore, special materials or structures are always adopted to resist the severe shock. As a submerged vehicle, submarine is always coated with multiple acoustic coverings [1] in order to improve the concealment, and thus it is of great significance to study the influences of the coverings and their antishock performance. Numerous researches on the covering performance were carried out and can be classified into mainly two categories: those on the acoustic performance such as the vibration and acoustic radiation of finite submerged cylindrical shell [2] and the acoustic radiation of laminated composite shells [3], and those on the antishock performance such as the dynamic response of coated double-hull subjected to underwater explosion [4], antishock characteristics of the sandwich structure subjected to underwater explosion [5, 6], the theoretical and experimental researches on antishock characteristics of a new cavity structure [7–9], the analyses of structural dynamic responses and failure modes based on finite element program [10], and the antishock performance of plate-frame covering [11, 12]. Owing to the high-cost and poor-portability of experimental research and the limitation of theoretical research, the numerical research is more popular. Due to the presence of large-deformation, moving interface, multiphase mixing, and other challenges, traditional methods based on mesh come across many difficulties in solving such problems [13–15]. However, as the oldest mesh-free method, characterized by Lagrangian particle property [16], SPH method has great advantages in dealing with these problems [17, 18].

Foam silicone rubber is chosen as the covering material since it has no significant influence on acoustic performance. SPH method is applied to study the influences of the covering on the damage characteristics of single and double cylindrical shells. Firstly, the feasibility of SPH method to solve strong discontinuity problems will be discussed in terms of dynamic continuity and motion continuity; and then the coated single and double cylindrical shells will be modeled with SPH method, the results of which will be compared with those of AUTODYN to verify the effectiveness of the present SPH model; on the basis of these, the parameters of covering position and covering thickness will be considered to discuss the influences on the structural damage of single and double cylindrical shells; more reasonable parameters of the coverings will be presented in terms of antishock performances to give a reference for the structural design of submarines.

#### 2. Theoretical Background

##### 2.1. SPH Approximation Equations

In SPH method, the approximation of function and its derivative can be discretized as [19]

Thus, the conservation of mass, momentum and energy in SPH can be expressed as (without regard to body force) [19] where , , , , , , and denote density, mass, velocity, energy, time, coordinates, and stress, respectively; and indicate the direction along axis; is the smoothed function of a pair of particles and , the cubic spline function is applied in present paper; represents artificial viscosity [19].

##### 2.2. Constitutive Equation

As for fluid, the stress is composed of isotropic pressure and viscous shear stress , namely, [19] where can be obtained from equation of state (EOS), Mie-Gruneisen EOS is used for water [20]; is the viscous shear stress; is the strain rate. As for Mie-Gruneisen EOS of water [20], is initial density of water; is sound velocity; is density ratio before and after the perturbation; is volume correction coefficient; , , and are fitting coefficients; is initial energy. The parameters for EOS of water are given in Table 1.

As for solid, the stress is composed of isotropic pressure and deviatoric stress tensor , namely, [19]

The solutions of and are explained as follows.

*(1) The Solution of *. The deviatoric stress can be drawn from the stress rate by integration, and yet can be gained from Jaumann stress rate [19], which can be noted as
where is shear modulus, is the strain rate, and is the rotation tensor rate.

Whether the stress should be updated or not is determined by Mises yield criterion [21], in which the Mises stress is compared with yield strength to judge the stress state of particle . If , particle is in elastic state, and there is no need to modify the stress component. On the contrary, if , particle is in plastic state, it is necessary to modify the stress component as [21] where means the yield strength, Johnson-cook model is applied to calculate the yield strength of the steel [22]. Meanwhile, the yield strength of the foam silicone rubber can be drawn from the dynamic constitutive equation [23].

In Johnson-Cook constitutive model of steel [22], is the static yield strength; and are strain hardening coefficients; is strain rate strengthening coefficient; is thermal softening coefficient; is equivalent plastic strain. , where is the second invariant of deviator stress tensor; is equivalent plastic strain rate; is reference strain rate; is dimensionless temperature which corresponds to ; and , where is room temperature and is melting temperature. , where is specific internal energy, is initial specific internal energy, and is specific heat; is shear modulus of steel. Some of the parameters are listed in Table 2. In the dynamic constitutive equation of the foam silicone rubber [23], is equivalent plastic strain rate; is reference strain rate; is equivalent plastic strain; , , , , , , and are fitting coefficients. The detailed parameters are shown in Table 3.

*(2) The Solution of *. The pressure of detonation products, steel and foam silicone rubber, can be obtained from Jones-Wilkins-Lee (JWL) EOS, Mie-Gruneisen EOS for solid, and Murnaghan EOS, respectively [24–27]. For JWL EOS for explosive gas [24], is initial density; is detonation velocity; , , , , and are experimental fitting coefficients; is detonation energy per unit mass. Then in Mie-Gruneisen EOS for steel [25], is initial density; is Gruneisen parameter; is linear participation coefficient of impact velocity and particles velocity; is slope; is initial energy of steel; is Poisson’s ratio. Finally, as for Murnaghan EOS for foam silicone rubber [26, 27], is initial density; is linear participation coefficient of impact velocity and particles velocity; is slope of Hugoniot curve [26]. And some parameters are listed in Tables 4, 5, and 6.

##### 2.3. Continuity Condition of Discontinuity Surface

According to the requirement of continuum on wavefront, the particle displacement is continuous, but its derivatives, such as velocity and strain, may be discontinuous for shock problems. In this section, the essential conditions of stress and strain on discontinuity surface, namely, the conditions of dynamic continuity and motion continuity [28], will be deduced with SPH method.

Given a finite region , as shown in Figure 1 [28], and denote the volume and the surface. The region is divided into two parts by a propagating shockwave , the volume and surface of the two parts are , and , , respectively; the wavefront velocity in the direction of outward normal is assumed to be , the corresponding particle velocity and normal velocity of the two parts are , and , , respectively.

The function is assumed to be continuous and differentiable in and , yet it may be discontinuous on . The particle position is a function of time, and the density is a scalar function of and , namely, .

###### 2.3.1. Dynamic Continuity Condition

With the propagation of shock wave, not only the function but also including and would change. Therefore, considering the derivative relations [28] and the particle approximation of function and its derivative, the above equation in can be discretized as Similarly, the equation in can be obtained as follows:

*(1) Replace ** with Function ** in (8) and (9), namely, *. Assume that and tend to the wavefront and approach zero, and thus, both the boundaries and tend to . Since the outward normal and are opposites when , should be in (8). Then, adding (8) to (9),

Due to the mass conservation in any volume and the arbitrariness of the wave front , (10) can be written as

*(2) Replace ** in (8) and (9) as **, *. Similarly, (8) is added to (9) as follows:

As a result of the opposite directions between and when , by applying Green formula, the momentum conservation equation regardless of the external force [19] can be derived as [28]

Substituting the above equations into (12), we have the following:

Owing to the arbitrariness of wavefront , we have

Substituting (11) into the above equation yields

Applying the particle approximation, the following equation can be drawn: which is the dynamic continuity condition [28].

###### 2.3.2. Motion Continuity Condition

The discontinuity surface turns to be after from . Therefore, attend , the discontinuous values on the discontinuity surface can be given as [28]

Then during , the rate of discontinuous value with respect to time is [28] as follows:

If , the equation above turns to be [28] as follows:

In general, the propagation velocity is much larger than the particle velocity and , and thus the particle velocity in (7) can be ignored [28] as follows:

Discretizing the two equations above, subtract (23) from (22) as follows:

Replace with the displacement , due to the continuity of the particle displacement, namely, , so we get

Therefore, restoring the discretized equation (25), in Cartesian coordinate, it turns to be as follows:

According to the Helmholtz velocity decomposing theorem [29], in the case of small-deformation, the deformation can be divided into linear deformation and angular deformation. Therefore, the motion continuity condition can be derived as [28]

In addition, we use a slight penalty force of Lennard-Jones model to solve interface problem, and the molecular force is so slight that it just prevents particles’ penetration. When particles on both sides of an interface tend to penetrate, in the case where , there would be a molecular force acting on two approaching particles, where is cutoff radius and it is generally close to the initial spacing of particles and is the distance between particles and . The direction of the is along the center line of the particles, preventing the particle from penetrating;, , and are set parameters. Consequently, we use the ratio of to to deal with the interface; this can guarantee good numerical stability without errors caused by the smoothing length.

#### 3. Contact Explosion of Single Cylindrical Shell

##### 3.1. Numerical Model and Parameters

The model of a coated single-hull submarine subjected to the attack of torpedo can be simplified as Figure 2, and (a) indicates no-covering, (b) denotes in-covering, and (c) represents out-covering. The relevant parameters and cases are listed in Tables 7 and 8. The origin of coordinates is positioned at the center of ring-shaped water. The column TNT is detonated from the center. The TNT dosage is about 12.8 kg, a popular amount in most model tests. The entire model is discretized with nonuniform particle spacing, 0.004 m for the cylindrical shell and 0.008 m for TNT and water, and some verification for the initial particle spacing is explained in Section 3.2. The numbers of different kinds of particles are listed as Table 8. The test point A is placed at (0, 2.1, 0).

**(a)**

**(b)**

**(c)**

##### 3.2. Numerical Verification

In order to verify the effectiveness of the present SPH method, the no-covering model of underwater contact explosion is selected, the SPH result will be compared with that of AUTODYN.

The propagations of shockwave obtained from different methods are shown in Figure 3 and the pressure of test point A is shown in Figure 4. It is obvious that the propagation laws of shockwave agree well with each other; moreover, the curves obtained from different methods also show good agreements, which all prove the effectiveness and accuracy of the present SPH method to solve the problems of underwater contact explosion.

**(a)**

**(b)**

The pressure of test point A in the cases of different initial particle spacing is shown in Figure 5. We note that the ratio of initial particle spacing between water and shell is and discuss the different cases (resp., = 1 : 1, 2 : 1, 3 : 1, and 4 : 1). The following figures are pressure-time curves of test point A in different cases and the verified model. As Figure 5 shows, the pressure-time curves in cases where = 1 : 1 and = 2 : 1 are basically consistent, but the case where = 1 : 1 requires a significantly larger amount of computation than the case where = 2 : 1. However, there are significant pressure noises at the interface in cases where = 3 : 1 and = 4 : 1 and serious numerical oscillation occurs when = 4 : 1. Therefore, = 2 : 1 is chosen.

##### 3.3. Results and Discussion

###### 3.3.1. Shockwave Propagation

After TNT detonation, the shockwaves propagates and overlaps in multilayer media, which will cause the inhomogeneity in time and space. According to the principle of impedance matching, the results above are in good accordance with the physics law, and the numerical results are fairing and smooth, which verifies that the dynamic and motion continuity conditions are satisfied in the SPH method.

For the single cylindrical shell, the pressure of point A is shown in Figure 6. It is obvious that the peak pressure of the out-covering is smaller than those of the other two cases. The reason may be that a rarefaction wave is reflected when shock wave transmits from TNT of high-impedance to rubber of low-impedance, which will result in low pressure.

The pressure of the in-covering with different covering thickness is shown in Figure 7. Define the dimensionless thickness as the ratio of covering thickness and the thickness of steel. When , the so-called “double peak” or “multimodality” of pressure is observed. With the increase of , the peak pressure goes up. When , it is up to the maximum, and does not increase with the increase of covering thickness.

###### 3.3.2. Damage Characteristics of Single-Layer Cylindrical Shell

*Damage Characteristics of Different Covering Positions*. The failure modes of shell subjected to underwater explosion can be divided into three modes [30]: Mode 1: plastic large-deformation; Mode 2: tensile tearing in outer fibers at the support; Mode 3: transverse shear failure at the support.

The equivalent plastic strain of coated shell is shown in Figure 8. As for (a), it is observed that the damage process experiences plug failure, dent, and rolling. As for (b), the shell firstly develops a pit with a greater plastic deformation and becomes ruptured at about ms. The size of the fragment is close to the contact area of TNT. As for (c), rupture occurs at the center of the shell first. With the increase of deflection, the tensile failure arises. Consequently, it is concluded that the main failure mode of no-covering shell is plastic large-deformation and tensile failure, while the mode will be tension failure and shear failure if the shell has coverings. The size of fragment is related to the contact area. In addition, the differences between in-covering and out-covering are initial rupture positions and initial rupture time.

**(a)**

**(b)**

**(c)**

*Damage Characteristics of Different Covering Thicknesses*. The equivalent plastic strain of different covering thicknesses is shown in Figure 9. For (a), (b), and (c) where = 0, 0.5 and 0.75, respectively, the main failure mode is plastic large-deformation and tensile failure, and the crevasse tends to roll; as for (d) where , the tension failure and shear failure should be the main failure modes, and the size of fragment is close to the contact area of TNT; as for (e) and (f) where = 1.25 and 1.5, which are similar to , the ruptures come out at different positions. In a word, when the covering thickness is close to that of steel, the protective effect is better.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

*Crevasse Radius*. The crevasse radius of different coated shells is shown in Figure 10. The dimensionless radius is defined as the ratio of crevasse radius to the TNT radius. As for Case 1, the crevasse radius increases through time, on the contrast, there is a little rebound process in Case 2 and Case 3. The turning point of Case 3 further proves that the main failure mode of out-covering is shear failure. In brief, the covering may not always take positive effects on the shell, better protective effect is found in the case of in-covering.

The crevasse radius of different covering thicknesses is shown in Figure 11. When , the crevasse radius increases with time first, and decreases with the increase of covering thickness. When , the crevasse radius undergoes a little rebound after the peak. When = 1.25, 1.5, a turning point turns out which also proves that shear failure is one of the main failure modes; in addition, the thicker the covering is, the larger will be. To sum up, the better protective effect appears when the covering thickness is equal to shell thickness.

#### 4. Contact Explosion of Double Cylindrical Shell

##### 4.1. Numerical Model and Parameters

Based on the researches of single cylindrical shell, the influences of covering position and interlayer water on the damage characteristics of the double cylindrical shell are studied. The basic model is shown in Figure 12; the detailed parameters and setting cases are listed in Tables 9 and 10. The column TNT is detonated from the center. The TNT dosage is about 12.8 kg, a popular amount in most model tests. The origin of coordinates is also positioned at the center of the model. The two test points are A (0, 2.1, 0) and B (0, 0.73, 0). In case (a), there is no interlayer water while case (b) has full interlayer water. The circular region represents the main research object. The whole model is discretized with nonuniform particle spacing, 0.003 m for the cylindrical shell and 0.006 m for TNT and water. And the numbers of different kinds of particles are listed in Table 10.

**(a)**

**(b)**

##### 4.2. Results and Discussion

###### 4.2.1. Propagation of Shockwaves

The propagation of shockwaves in the case of in-covering on inner surface is shown in Figure 13. Shock wave propagates from water to inner shell at ms; at the same time, a reflected shock wave is generated on the lower surface of inner shell. When arriving at outer shell, a reflected shock wave is generated on the upper surface of outer shell, which is shown in the rectangular region. At ms, not only reflects a shock wave on the upper surface of outer shell (shown in the rectangular region) but also generates a shock wave (shown in the elliptical region) transmitting to the water; the pressure in the circular region is larger because of the superposition of , , , and . In the process, it is obvious that the shock wave velocity in steel is larger than that in water or rubber, yet the velocity in water is close to that in rubber because of similar impedances. The impedance mismatch causes the complexity of shock wave propagation in different media.

**(a)**

**(b)**

**(c)**

The pressure at test point A is shown in Figure 14. As for 1, due to the air-backed outer shell, a rarefaction wave is generated, which will decrease the pressure at test point A. However, as for 2, the rarefaction wave takes little effects on the test point. So the peak pressure of 2 is larger than that of 1. In addition, as for 2, the second peak pressure appears at ms because there is a reflected shock wave transmitting from the inner shell to water.

The pressure at test point B in the case of full interlayer water is shown in Figure 15. It turns out that the most notable feature is multimodality and the peak pressure decays with time. It is not beyond imagination that the superposition of incident wave, reflected wave, and transmitted wave appears at inner shell, outer shell, and the water between them.

###### 4.2.2. Damage Characteristics of Single Cylindrical Shell

*Crevasse Mode*. Failure modes 1–3 are plastic large-deformation, tensile failure, and shear failure, respectively. Moreover, the shear failure can also be classified to three categories 3a–3c. The failure modes of double cylindrical shell coated on different surfaces are shown in Figures 16 and 17.

**(a)**

**(b)**

**(c)**

**(d)**

As for the outer shell, the damage process experiences plug failure, dent, and rolling. As for the inner shell, in the case of no interlayer water and no-covering, the failure of inner shell is caused by high-speed fragments from outer shell, which is called “3a,” as shown in Figure 16(a); in the case of full interlayer water and no-covering, the crack appears at an angle of 45° from the normal direction, which is called “3c,” as shown in Figure 16(b); in the case of full interlayer water and in-covering of inner shell, only plastic large-deformation comes out, as shown in Figure 16(c); in the case of full interlayer water and out-covering of inner shell, the crack appears in the center of inner shell in normal direction, which is called “3b,” as shown in Figure 16(d). To sum up, when the interlayer water is full, the in-covering of the inner shell will better protect the inner shell and the interlayer water can reduce the damage caused by fragments.

*Crevasse Radius*. The crevasse variation of the outer shell without interlayer water is shown in Figure 18. In case 1, without covering, the crevasse increases as time goes on. However, in other cases, there is a little rebound when up to a peak and tends to be steady finally. Obviously, the out-covering provides better protection to the outer shell.

The crevasse variation of the outer shell with full interlayer water is shown in Figure 19. In cases 6, 7, and 8, the crevasses increase through time, and the coincidence indicates that the protective effect on outer shell is not obvious when inner shell is coated. Nonetheless, in cases 9 and 10, the crevasse is smaller, which shows that the covering on the outer shell is beneficial for the outer shell. In addition, the water has certain influences on the outer shell because the water would participate in energy absorption.

The crevasse variation of the inner shell without interlayer water is shown in Figure 20. The turning points account for the fracture of the inner shell. In case 9 with in-covering on the inner shell, the crevasse is larger than case 8 without covering, which indicates that the covering does not always take positive effects.

The crevasse variation of the inner shell with full interlayer water is shown in Figure 21. Except cases 9 and 10, the curves appear to overlap, which indicates that the protection is not obvious. It seems that the effect of in-covering on the inner shell is better; however, considering the operation in practice, out-covering on inner shell would be a practical choice. In addition, comparing the failure modes and the sizes of crevasse between cases with/without interlayer water, it is found that the inner shell can be protected by the water from the damage caused by fragments.

#### 5. Conclusions

An SPH method with mesh-free and Lagrange properties is applied to solve extremely dynamic problems of cylindrical shell subjected to underwater contact explosion in this paper. The influences of coverings on damage characteristics of cylindrical shells are investigated. Through the parametric analyses of covering position and covering thickness, the following conclusions can be drawn.(1)The propagation laws of shockwaves in multilayer media correspond with the theory of impedance matching, and the continuity conditions of discontinuity surface are verified by the smooth SPH results; the results of SPH agree well with those of AUTODYN, which verifies the feasibility and effectiveness of the present SPH method.(2)For single cylindrical shell, the failure mode will be large-deformation and tensile failure when the covering is thinner than the shell, and tension failure and shear failure otherwise. The antishock performance does not necessarily increase with covering thickness; better protective effects appear when the covering thickness is close to that of the shell.(3)As for double cylindrical shell, the obvious protective effect is found when there is out-covering on the outer shell or in-covering on the inner shell. However, considering the operability, the out-covering on inner shell is also a good choice when the interlayer is full of water.(4)The damage of single shell is mainly caused by detonation wave and detonation products; as for double shell without interlayer water, the inner shell is mainly disrupted by high-speed fragments. Moreover, a certain protective effect will be provided by the interlayer water.

#### Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

#### Acknowledgments

This work is supported by the Excellent Young Scientists Fund (51222904), the National Defense Basic Scientific Research (B2420133001), and the National Security Major Basic Research Program of China (613157).