Pressure hull analysis under shock loading

The hull of high performance submarines must resist underwater shock loading due to exploding torpedoes or depth bombs. An underwater shock involving an initial shock wave andsuccessive bubble pulsating waves is socomplex that a theoretical technique for deriving shock pressure distribution is required for improving simulation efficiency. Complete shock loading is obtained theoretically in this work, and responses of a submarine pressure hull are calculated using ABAQUS USA (Underwater Shock Analysis) codes. In the long run, this deflection and stress data will assist in examining the structural arrangement of the submarine pressure hull.


Introduction
A typical submarine structure system can be deconstructed into four parts: outer shell, pressure hull, sail, and rudder. The pressure hull is typically composed from fore to aft of three conventional pressure hulls: dome, cylinder, and cone. Pressure hulls must resist high water pressure to protect humans and instruments and, thus, require critical strength. During attack, antisubmarine weapons, such as depth bombs and torpedoes, generate powerful explosive forces that can seriously damage a submarine's structure. This work focuses on pressure-hull strength, an important component of a submarine's structural design, under shock loading of the submerged submarine.
When an underwater explosion occurs, the explosive charge in the bomb is rapidly converted into a gas at a very high temperature and pressure. An explosion's initial shock wave is followed by a succession of pressure waves at 10-15% of the initial shock-wave pressure. These waves result from repeated expanding and contracting of the high-pressure gas. The initial shock wave, which has a short duration and a high-pressure variation, can be considered as an extremely brief transient load with respect to its action on a vessel. Conversely, although the shock waves resulting from the bubble pulses are substantially smaller than the initial shock wave, they typically last longer, and result in vibratory submarine motion, since the duration of a bubble-pulse cycle can be close to the fundamental vibration period of a ship. In this work the initial shock wave and bubble-pulsating waves are applied as the total load during underwater shock analysis.
Although finite element method (FEM) packages, such as ABAQUS with an Underwater Shock Analysis (USA) module, can solve for the explosive load -including the initial shock wave and bubble-pulsating waves -massive computing resources are required. This work attempts to obtain theoretically the total incident pressure, and uses this pressure as the input for ABAQUS. Figure 1 presents the methodology flow chart used to solve the underwater shock problem. An initial shock wave (see left-hand column of Fig. 1) can be defined effectively using empirical equations from the literature; however, the bubble-pulsating waves must be extracted from a series of theories (see right-hand column). Then ABAQUS software is utilized for computations. Each theory portrayed in Fig. 1 is described separately in this paper.
This work utilizes the pressure-hull model that simulates the German HDW U209-1400 submarine (Fig. 2). Analysis results are used to evaluate the submarine's ability to resist underwater explosions.

Empirical equation for initial shock wave
The initial underwater shock wave can be approximated using an exponential function that describes the pressure history of a shock wave as observed from fixed location. The shock wave comprises an initial, instantaneous increase in pressure, followed by an exponential decline. This equation is well known as Cole's empirical equation [1]: where P sw (R, t, W ) is the shock-wave pressure at the observation point, R is the distance between observation point and explosion center (m), W is the weight of explosive charge (kg), P max (R, W ) is the maximum pressure (Pa); t is the time elapsed since the arrival of the shock; and θ e (R, W ) is the decay time constant, which is the time of observed pressure decays from P max to P max /e. Peak pressure P max and decay constant θ e are depicted by Eqs (2) and (3), respectively: and where α 1 , β 1 , α 2 and β 2 are constants based on the explosive charge type. For a TNT charge, α 1 = 5.212 × 10 7 , β 1 = 1.18, α 2 = 0.0895 × 10 −3 and β 2 = −0.185.

Bubble-pulse dynamic equation
The bubble-pulse dynamic model can be addressed using two principal approaches: considering the surrounding fluid as non-compressible, or as compressible. In 1917, Rayleigh [2] where a is the bubble radius, P g is the bubble interior pressure, P h = ρgH 0 (static water pressure at the bubble center H 0 ), and c is the sound speed in the fluid. Solving Eq. (4) yields the bubble radius history, where the interior bubble pressure P g comes from the gas state equation (discussed in the following section).

Gas state equation
The equation describing the gas pressure variation in the bubble following an explosion is called the gas state equation. In 1948, Jones and Miller [4] utilized the assumption of heat-and chemical energy conservation for an entire explosion to develop the Virial Expansion theory that yields the pressure of bursting gas. Later, in 1950, Taylor [5] applied this theory and, assuming that the specific heat ratio γ g is constant during an explosion, obtained the Jones-Miller-Taylor Equation of State (JMT EOS): where K is a parameter varying with dynamite type and charge density, W is the charge weight, and V g is the bubble volume. Table 1 defines the parameters for a TNT charge.

Initial condition for the bubble pulse dynamic equation
To solve the bubble-pulse dynamic equation, Eq. (4), appropriate initial conditions, including initial bubble radius a(0) and explosion speedȧ(0), are needed. In ideal conditions, assume that a charge explodes and transforms into a high-temperature, high-pressure gas with the same mass and volume as that of its original solid state. The gas then pierces the packing shell of the dynamite and explodes into the surrounding water. In this case, a(0) = a c (a c is the charge's radius) andȧ(0) = 0.
However, the bubble-pulse dynamic equation cannot effectively describe an explosion at the initial explosion time due to the complexity of the energy conversion event. Hunter [6] proposed utilizing the condition at the time ranging between 3θ e and 7θ e as the initial condition.
In 1950, J.S. Coles [1] used experimental results to modify Eq. (1) to The applicable time range in Eq. (6) is the time when the initial condition is given, which can be set as 3θ e to 7θ e . (t * is the time when the initial condition is given).
In 1956, Strasberg [7] surveyed acoustic theory for a spherical surface with large movement and solved the relationship between the bubble radius and the shock pressure aṡ and Equation (6) yields maximum pressure P max . Substituting P max into Eqs (7) and (8) produces the initial conditions a(t * ) andȧ(t * ).

Empirical formula for bubble dynamics [8]
For TNT charges, an approximation of the maximum bubble radius (a max ) and the first period of the bubble-pulse wave (T 1 ) can be expressed using the empirical formulae below in Eqs (9) and (10): and W is the weight of the charge in newtons, and D is the depth of an explosion in meters. These two empirical equations can be used for a comparison after extracting the bubble motion histories theoretically.

Pressure variation caused by bubble motion
To allocate pressure variation caused by bubble motion a(t), assumptions are required to simplify the physical phenomenon. The liquid density ρ L and the dynamic viscosity µ L are assumed constant; and the contents of the bubble are assumed homogeneous, that is, to have uniform temperature T B (t) and pressure P B (t) [9].
Outside the bubble, P bu (R, t) and T (R, t) denote pressure and temperature respectively; and u(R, t) denotes the radial velocity outward for the bubble. Conservation of mass requires that where F (t) is the kinematic boundary condition at the bubble surface. In the ideal case of zero mass crossing the bubble surface, it is clear that u(a, t) = da/dt; hence, Equation (12) is typically an effective boundary denotation, even when considering evaporation and condensation of a fluid occurring continually on the interface. Newly produced vapor volume must equal bubble volume variation, that is, 4πa 2 · da/dt; therefore, the mass flow rate of evaporation is where ρ V (T B ) is the saturated vapor density at bubble temperature T B . Expression (13) must equal the mass flow rate of liquid moving inward to the interface and, hence, the inward velocity of liquid can be obtained as (11) and (12) are then further expanded into Eqs (14) and (15): In practical cases, ρ V (T B ) << ρ L and, therefore, the approximate form of Eq. (11) is sufficient and will be employed hereafter.
Assuming seawater is a newtonian liquid, the Navier-Stokes equation for motion in the R direction is Substituting the expression for u from Eq. (11) generates Upon integrating Eq. (17), one obtains: where P ∞ is the pressure far from the explosion; this is commonly the static water pressure at explosion depth. Finally, an equation for water pressure variation caused by bubble pulsating motion can be derived using the condition P bu → P ∞ while R → R ∞ :

The fluid-structure interaction effects
Suppose a condition of a given structure (e.g., a submarine or ship) submerged or semi-submerged in water; the water, as an ideal medium, should transmit linear wave motion [10]. When P i is the incident wave, the diffracted wave due to a rigid body is P d and the reflected wave from the water surface is P p ; the pressure field in the water is represented by Eq. (20): On the wetted surface (S 0 ) of a ship, the boundary condition is where n is the normal unit vector on the structure surface, ρ is the fluid density, u is the surface displacement of the structure andü n is the normal acceleration of the structure surface. At the water surface S p , the characteristic impedances of air and water vary significantly, such that the pressure wave cannot pass through and, therefore, is regarded as zero: At infinity, the following Sommerfeld radiation condition should be satisfied: where r is the distance between the observation point and explosion center, i = √ −1, and k is the wave number. Furthermore, the wave motion should satisfy the following Helmholtz equation: where c is sound speed, and ∇ 2 is the Laplacian operator. In this work, the incident wave P i comprises the initial shock wave and bubble pulsating wave, and can be configured as P i = P sw (R, t, W ) + P bu (R, t). (25)

Coupled fluid-structure interaction equations [11]
To solve the coupled fluid-structure interaction problems, the FEM software ABAQUS combined with the USA module, which is theoretically based on the Doubly Asymptotic Approximations (DAA), is utilized. During the solution process, the USA structural interface (USI) elements should be constructed on the structure surface that is connected with the ship structure and DAA boundary. The incident wave spreads to the DAA boundary, and the USA transforms the wave into pressure acting on the structure surface. Next, ABAQUS performs structure transient analysis and USA solves the coupled fluid-structure interaction equations based on displacement and velocity of nodes.
The coupled fluid-structure interaction equations are and where {q s } =

Model for analysis
The pressure hull model was adopted for the Germany HDW U209-1400. A typical U209 submarine is constructed with a single pressure hull composed of HY80 high-tensile steel. The U209 nominally carries 8 torpedo launch pads and 14 torpedoes; some U209s can launch anti-vessel missiles. Table 2 lists the principal specifications of the U209. Figure 3 presents the finite element model of the pressure hull constructed for this research. The black triads in Fig. 3 are damping boundaries that will be described in Section 13. A conductor tower and two access hatches are included; but the outer shell is omitted, as it does not resist water pressure. Although detailed structural configurations of the U209 cannot be obtained, the pressure hull dimensions including the outside surface, ribs and other stiffeners, were carefully constructed to represent the actuality as closely as possible (Fig. 2). The model has the capability of diving to 400 m based on an analysis considering geometrical nonlinearity and plastic material properties.

Depth bomb
Assume that a depth bomb, the explosion source, has TNT charges weighing 160 kg and a radius of a c = 0.1625 m. The TNT charge density is 1500 kg/m 3 , with parameters K = 77.7E+3, and γ g = 1.27 in Eq. (6). This bomb is exploded at 15.81 m beside the submarine hull while the submarine is at a depth of 250 m. In this scenario, the Keel Shock Factor (KSF) is 0.8, and conforms to the naval vessel design standard. The KSF, which indicates the degree of shock, is defined with Eq. (28): where W is the charge weight, L is the distance between explosion center and ship keel, and θ is the angle between the explosion center, free surface and ship Referring to Fig. 4, θ SUB is employed here for a submarine while θ SHIP is adopted for a surface ship.

Results of shock wave evaluation
The histories of bubble motion and pressure variation for the simulated depth bomb are obtained theoretically (Fig. 5) and by using Eqs (4) and (19), which include the bubble radius, the swelling/contracting velocity, and the acceleration and pressure variation observed at the standoff point -a standoff point is the nearest point to the explosion center on the submarine surface. In Fig. 5, the time origin is the time at which the pressure wave first contacts the submarine surface. Simulation results indicate that the maximum shock pressure reaches 1.6034e+7 Pa and then declines immediately. Very small negative pressure (suction) also exists sometimes. At t = 1.055e-1 sec, the bubble contracts significantly to a minimum radius, followed by high-acceleration swelling that induces another shock wave similar to an explosion.
The empirical values of maximum bubble radius and first period of the bubble pulsating wave, derived by Eqs (9) and (10), are a max = 2.9 m and T 1 = 0.1111 sec, respectively. The theoretical results for maximum bubble radius and first period of the bubble pulsating wave are a max = 2.65 m and T 1 = 0.1055 sec, respectively. The closeness of these two results demonstrates the accuracy of the strategy applied in this work.
Additionally, the LS/DYNA commercial software is used to solve the Jones-Wilkins-Lee Equation of State (JWL EOS) using the same bomb and exploding position settings. Figure 6 presents the history of shock pressure variation compared with theoretical results. Maximum pressure P max (5.0 MPa) calculated by LS/DYNA is higher than that (4.0 MPa) calculated using the theoretical approach, although both peaks occur at almost the same time. The difference between two maximum pressures is thought to be due to the fact that the P max derived by Eq. (19) required the data for bubble radii a,ȧ andä, which were obtained from Eqs (4), (7) and (8). The bubble radius formulae were derived from integration procedures that likely accumulate errors and lead to a smaller bubble radius and, consequently, a smaller pressure than the LS/DYNA's results.

Damping boundary condition for pressure hull
Some clamped boundary conditions for a submarine hull are incorrect for an underwater shock problem. Submarines vibrate and move slightly during a shock wave. Moreover, the fluid surrounding the submarine has added mass and damping effects. In this work, the added water mass is addressed by the ABAQUS USA; however, the damping forces acting on the submarine must be determined.
The total damping coefficient in the heave direction, C HD [12], can be defined by where ρ is the fluid density, B is the ship's width, U is velocity, and K L is the first-order differentiation of the C L − α chart (chart of lift coefficient vs. angle of attack for a submarine cross section). Calculate C HD utilizing the empirically estimated velocity of the submarine subjected to a specific load. This work assumes that damping forces are distributed uniformly on the submarine surface. Total damping forces are simulated using 14 sets of connector elements with damping. These connector elements are placed on three hull sections (front, middle, and back sections) with four groups for each section; the remaining two groups are on the bow and stern. Each group contains three connectors that are pointed toward three space coordinate directions. The damping coefficients (force per relative velocity, N · sec/m) can be reasonably specified in the range 10 2 -10 3 ; here, the value employed is 10 3 . One end of a connector is attached to the hull and the other is fixed in space. This arrangement facilitates finite movement of the pressure hull. The black triads in Fig. 3 present the damping connector locations.
Furthermore, based on preliminary calculated results, when stress is concentrated on the connected point, more numbers of connector can be utilized to substitute one connector while decreasing the damping coefficient for each connector.

Results of hull response
In Fig. 7, the target point indicates the bomb positioned beside the hull, and the arrows indicate shock wave pressure distribution.
According to calculated response results (Figs 8 and 9), the submarine first deforms and shifts in the force direction. At t = 0.005 sec, the instantaneous acceleration is 650 G. As the shock wave decays, acceleration decreases; however, at t = 0.11 sec, acceleration again peaks at −250 G due to the bubble pulse.
Next, the maximum stress (420 MPa) at the standoff point occurs at t = 0.02 sec, at which time the region surrounding the standoff point also attains maximum displacement. Computational results show that the front and back extremities of the hull deform less than the standoff point, which is located at the middle hull. At t = 0.07 sec, the displacement is in the opposite direction: the two extremities have larger displacements than the standoff point, implying that the fundamental-mode vibration period for this hull is 0.1 sec, a short period.
Stress and displacement are described as follows. When t = 0.0 sec, the shock wave has not yet contacted the hull and, thus, stresses and displacements are zero. When t = 0.01 sec (Figs 10 and 11), the pressure hull is affected by the initial shock wave, inducing an elevating stress (590 MPa) near the forward access hatch and yielding the local hull. When t = 0.04 sec (Fig. 12), the stress concentrates at the bottom of conductor tower and the aft access hatch.
When t = 0.12 sec (Fig. 13), with the bubble pulsating wave passing through the hull, stress reaches 575 MPa, causing the plates at the conductor tower bottom to yield. Finally, at t = 0.2 sec (Figs 14 and 15), the stress becomes low and the next pulsating bubble wave comes into play. At this time, total shift and deformation at the midship section is 24 cm; and at the stern, 27 cm.

Conclusion
(1) The FEM packages may be more mature and correct, and take more computing time than the theoretical method when solving for the dynamic bubble phenomenon. In this research, the theoretical method, which was applied to solve the pulsating dynamic bubble equation and the pressure variation, obtained a satisfactory solution compared with a FEM solution, and consumed less computer resources. This scheme is recommended for use during early design stages. (2) Three bomb positions were tested: beside, below, and above the hull. When the bomb explodes below the submarine, it produces the most serious damage to the hull. When the bomb explodes above the hull, the damage is slighter than the other two cases.