Survey on Experimental and Numerical Approaches to Model Underwater Explosions

Abstract

The ability of predicting material failure is essential for adequate structural dimensioning in every mechanical design. For ships, and particularly for military vessels, the challenge of optimizing the toughness-to-weight ratio at the highest possible value is essential to provide agile structures that can safely withstand external forces. Exploring the case of underwater explosions, the present paper summarizes some of the fundamental mathematical relations for foreseeing the behavior of naval panels to such solicitation. A broad state-of-the-art survey links the mechanical stress-strain response of materials and the influence of local reinforcements in flexural and lateral-torsional buckling to the hydrodynamic relations that govern the propagation of pressure waves prevenient from blasts. Numerical simulation approaches used in computational modeling of underwater explosions are reviewed, focusing on Eulerian and Lagrangian fluid descriptions, Johnson-Cook and Gurson constitutive materials for naval panels, and the solving methods FEM (Finite Element Method), FVM (Finite Volume Method), BEM (Boundary Element Method), and SPH (Smooth Particle Hydrodynamics). The confrontation of experimental tests for evaluating different hull materials and constructions with formulae and virtual reproduction practices allow a wide perception of the subject from different yet interrelated points of view.

1. Introduction

Acknowledging how to properly soften the effects of impact-related damage is an imperative design guideline in shipbuilding. Specially for military applications, underwater explosions (UNDEX) prevenient from subsea blast loads can infer irreparable structural impairment to vessels, where the dynamic response of ships depends on the influence of several parameters, such explosive power, blast distance, hull panel composition, and reinforcements. The mechanics involved in ship collisions and impact can thus be looked at from two inter-dependent perspectives [1], where the external mechanics describe the hydrodynamics surrounding the vessel and the magnitude with which they affect it, and the internal mechanics that regard how the materials respond to these forces and dissipate energy through strain.

The necessary comprehension of the materials behavior goes beyond traditional steels, whereas modern designs utilize advanced metallic alloys [2], composite structures (either purely polymeric, as in the case of small and medium sized boats [3], or in a sandwich layout composed by metal sheets and foam cores [4]) and wood [5]. Also, the usage of innovative materials, such as polymers reinforced by natural fibers, have been addressed in the literature, on one hand presenting inferior mechanical properties to carbon or glass fiber, and on the other having enhanced sustainability [6]. In addition, the diverse existent hull architecture requires the designer not only to master different materials, but also to understand the response of the vessel itself as a composition of the several interlinked structural elements built with those materials, such as panels, stiffeners, and girders. The way in which these structures are realized and arranged ultimately determines the resistance of the ship and define if reinforcements like stiffeners and girders apply, and what is the optimal geometry for plate cutouts (e.g., for hatches) to avoid undesired elastoplastic deformations to buckling [7,8,9].

External mechanics, in turn, involves a series of differential hydrodynamic relations that describe the propagation of intermittent pressure waves caused by explosions defining the velocity, temperature, period, and force with which hulls are hit. The consequent fluid-structure interactions with relevant severity do not only involve liquid matter; gas bubbles generated by the blasts can also cause noticeable damage given the high pressure and temperature they may present [10,11].

The relevance of UNDEX for military purposes is clear, as highlighted in studies dated from 1990 [12], which explicitly express the cold war as the impetus for such research initiative, stating as border condition for the case study a soviet missile impacting a double-hulled structure. The aforementioned work defined through a series of experiments the strongest parameters that influence UNDEX effects as being: enthalpy of detonation, number of moles of gas, molecular weight of gas, and solid phase density.

The strategic usage of the arctic ocean as an operational area for submarines has also had its importance reflected by studies such as Barash’s [13], which analyze the effect of UNDEX beneath artic ice, not only to improve the safety of those vessels to enemy blasts, but also to allow them to have quick access to the surface by breaking ice packs using explosives. Potentially harmful debris originated from the blasts were addressed by Bryant [14], concluding that little or no metallic debris are deposited in the medium during early bubble oscillations, whereas they are actually transported by the bubbles and released in a certain area as the oscillations terminate.

Early analyses on the physics of the pressure waves were also performed. Snay [15] confirmed that the pulse shape of waves varies during their propagation in the water, where the steepness of the pressure at the head of the wave rises until the formation of an impact front. However, the author states the validity of approximating UNDEX to acoustic waves (which do not change their shape) for the solution of interaction problems, modelling them as small-amplitude constant-sound-velocity waves.

The scientific way of approaching UNDEX problems itself is an issue with distinct propositions in precursor military-grade reports. O’Daniel et al. [16] evaluated the mechanics of bubble jets, suggesting that such analysis could be conducted by performing three experiments of bubble jet strike on vertical targets in different scales: a small-scale one to study the phenomenology of the bubble jet; a mid-scale one to measure the loads applied by these jets on targets; and a large-scale one for the structural response of targets. Miller [17] categorized three fundamental steps: comparison of structural response to conventional and nuclear pressure profiles, followed by finite element analysis, and then by writing an appropriate computational numerical code to describe the experiments. Naturally, since then, computer-aided simulations have been immensely improved, as shown hereby.

Built upon classic mathematical models, numerical simulation software is a well-established time-saving tool for yielding precise UNDEX results, leaving aside the need for resource costly experimental tests. In marine engineering, accurate virtual reproductions can be done for UNDEX provoked pressure waves [18] and bubbles [19] over vessels, as well as water blast on sandwich panels [4] and cavitation in propellers [20]. For that, a variety of constitutive material models and numerical solving techniques may be used. Given the importance of understanding how external and internal mechanics work and relate, the present work briefly summarizes the fundamental numerical expressions of each one, giving also a glance in how they are usually represented in numerical simulations.

Portraying the three independent but correlated parts of UNDEX described by the navy as shown in Brenner et al. [21], i.e., the responses of water, of the vessel structure, and of the materials to the shock waves, Figure 1 depicts through a flow chart scheme the aspects approached in this work in an orderly manner. In other words, the current state-of-the-art explains how the pressure impulse propagates through water and with which magnitude it hits the ship (taking into account pressure waves and bubbles), how the ship structure responds to this impulse, and in which way the material components of the hull absorb this energy, deform, and eventually fail.

The characterization of such phenomena is hence made through analytical models, and directions on how to reproduce all the three steps by numerical simulation are also given. Some fundamental notions are provided on Eulerian, Lagrangian, and Aleatory Lagrangian-Eulerian (ALE) fluid descriptions on Finite Element and Volume Methods (FEM and FVM, respectively), Boundary Element Method (BEM), and Smoothed Particle Hydrodynamics (SPH) numerical modelling techniques, and on Johnson-Cook’s and Gurson’s constitutive materials.

Acknowledging that the study of UNDEX phenomena is an intricate and wide field of study with deeper aspects to be considered, some of its basic yet most important parameters are hereby summarized and briefly explained, intending to allow the comprehension of fundamental notions of the subject based on a thorough literature survey.

2. Analytical Models

2.1. Hydrodynamics and Fluid-Structure Impact

Amongst the challenges in ship design, there are several forces to be taken into account that may act either separately or combined over the structure originated from gravity, wind pressure, friction, freight, and hydrostatic lifting force, which can infer massive shear and torsional efforts given the large size of a ship and due to its weight distribution. Moreover, fluid-structure interactions (FSI), such as pressure waves, gas bubbles, cavitation, and underwater explosions also have the potential to cause significant damage, given that water transmits explosive energy much more efficiently than air [16].

Differently from solid-particle impacts, where a specific damage is locally caused by a particularly shaped projectile [22], in UNDEX, all the extension of the hull reachable by the same pressure wave is approximately equally impaired. These spherical waves originate a mass flow behind them (i.e., the afterflow), which leads to the formation of bubbles that can be regular UNDEX bubbles, steam bubbles in the case of nuclear explosions, or cavitation bubbles in the case of explosions that do not provoke gases [15].

Bubbles generated by UNDEX, for instance, which constitute hot air masses that can hit a vessel in more than 10 consecutive waves [23], can reach impact pressures of up to 15% of the wave pressure [11], causing significant damage and even a potential and dangerous match in resonance frequencies with the ship. Furthermore, Cui et al. [24] experimentally analyzed the mechanics of these bubbles, finding out that the pressure peaks they induce on the impacted body is highly dependent on their shape before collapse, given that the less spherical and more asymmetrical it is, the smaller is the pressure due to the influence of splashing jet water. The shape of bubbles, as shown by Gong et al. [19], is directly related to the position of the explosive charge in relation to the hull.

The response to underwater explosions is a topic of interest for military purposes, given that even if the hull is not directly hit by a torpedo or collided to other structures, the sole propagation of energy from a far blast center through water is capable of causing significant damage to the ship. This energy release can include gas bubbles at nearly 3000 °C and pressures up to 5 GPa, besides solid particles made from lead or alumina, for instance, which altogether may impact the hull’s surface at nearly the velocity of sound [10].

It is possible to quantify the pressure in the proximity of the wave front P(t) if the charge is detonated at less than 100 m in depth [25] by the expression of Equation (1) [10]:

$$P\left(t\right)=P_m\{\begin{array}{c}{e}^{-t/\theta }ift<\theta \\ \frac{0.368\theta }{t}if\theta <t<t_1\end{array}$$

(1)

where the time constant $\theta$ can be defined by Equation (2) [10]:

$$\theta =R_0\{\begin{array}{c}0.45\overline{r}{10}^{-3}if\overline{r}\le 30\hfill \\ \frac{3.5}{c}\sqrt{log\overline{r}-0.9}if\overline{r}>30\hfill \end{array}$$

(2)

where $\overline{r}$ is the ratio between the distance of the explosion center and the quotient of the measured location R by the initial radius of the spherical explosive R₀. The peak pressure P_m (MPa) and the time decay constant $\theta$ (ms) can also be defined as a function of the charge load W and the standoff distance S (Equations (3) and (4), respectively) [26]:

$$P_m=\alpha {\left(\frac{W^{1/3}}{S}\right)}^{\beta }$$

(3)

$$\theta =0.058W^{1/3}{\left(\frac{W^{1/3}}{S}\right)}^{-0.22}$$

(4)

where $\alpha$ and $\beta$ are constants, which highlights that it is important to notice the discrepancies found in the literature for them. While [26,27] adopt $\beta$ as 1.13, [18] considers 1.81. By its turn, $\alpha$ is adopted as 148.93 in [26], 52.16 in [27], and 29.9 in [18]. In addition, [26] considers Equation 3 valid when the measuring point is located until 10 times the explosion radius, while [18] states that the adequate tolerance distance is up to 6 times the charge radius.

When the underwater wave passes through a certain point in space, it is then submitted to a transient pressure P(t) and displaced in a velocity v(t) in the direction of the flow. Considering a spherical flow, it is possible to estimate and correlate both into the same time-dependent Equation (5) [26], where the first term stands for the velocity of a plane wave, and the second, called “after flow” term, is important for large time intervals and for when the measuring point is close to the explosions, i.e., the standoff is low.

$$v\left(t\right)=\frac{P\left(t\right)}{\rho c}+\frac{1}{\rho S}{{\displaystyle \int }}_0^{t}p\left(t\right)dt$$

(5)

Although the strongest shock wave comes along with the first UNDEX impact, it is important to acknowledge that the explosion damage is characterized by repetitive impacts originated in the explosion center, which intensity decays exponentially when hitting a marine structure [28], and the second impact can even be the most damaging one due to the bubble pulse, especially to air-backed hull panels [29]. This concept is illustrated in Figure 2.

2.2. Structural Response

Structural elements are meant to be strategically assembled in the construction of ships in places where they can aggregate normal, shear, or torsional resistance to the vessel, preferably in profiles and shapes that favor their mechanical properties and account for the lowest weight attainable. The main structures are longitudinal and transversal girders and stiffeners, providing enhanced resistance to underwater or contact explosions compared to unstiffened vessels [30]. These can grant noticeable economy without compromising strength and durability if stiffeners are adequately positioned at regular gaps, while dividing the hull shell into lattices in which the intersection point between transversal and longitudinal girders are the most resistant zones [31].

Hung et al. [30] have demonstrated experimentally the importance of stiffeners by performing underwater explosion tests on 3 cylindrical specimens, one being unstiffened, one externally, and one internally stiffened with welded rings. The latter showed much smaller strain, pressure, and fluid acceleration inside the coupon. Ming et al. [18] have complemented this approach by numerically and experimentally studying the reaction of plates reinforced with “T” profile stiffeners placed in the opposite side of an underwater explosion, demonstrating the damage steps of the reinforced plate. The study states that after both the plate and the stiffeners are initially damaged due to pressure waves, the panel in the lattice is torn, leading to the collapse of the longitudinal and transversal stiffeners joints, weakening the whole structure significantly.

Gordo et al. [32] have developed a numerical method to estimate the ultimate longitudinal strength of the hull girder, yielding accurate confrontations against similar approaches and the real ship failure. The compression strength of stiffened plate columns at a determined strain is defined by two independent expressions that describe flexural buckling (${\varphi }_F$) and lateral-torsional buckling (or tripping) (${\varphi }_T$). The first is described by Equation (6):

$${\varphi }_F\left({\epsilon }_i\right)={\varphi }_{JO}\left({\epsilon }_i\right)\frac{A_s+w_{e}bl}{A_s+bl}$$

(6)

where ${\varphi }_{JO}$ is based on the Johnson-Ostenfeld formulation that considers inelastic effects during buckling [33], and the second term of the equation stands for the actual degradation of the plate due to compression loading. $A_s$ is the sectional area, b is the breadth of plating between longitudinal stiffeners, l is the length of the stiffeners, and $w_e$ is the effective width of the plate given by Equation (7):

$$w_e=max\left(-1,min\left[1,{\epsilon }_i\right]\right)\left(\frac{2}{\beta }-\frac{1}{{\beta }^2}\right)$$

(7)

where $\beta$ is the slenderness of the plate equal to $\beta =\left(b/l\right).{\epsilon }_v$, in which the second term represents the instantaneous strain. In turn, the tripping strength can be calculated by Equation (8):

$${\varphi }_T\left({\epsilon }_i\right)={\varphi }_{Tmin}\frac{{\epsilon }_t}{{\epsilon }_i}\frac{A_s+w_{e}bl}{A_s+bl}$$

(8)

The first term stands for the maximum elastic tripping stress, and the second to the strain of maximum load. It is important to highlight that the resistance of the girder can be affected by residual stress and external agents such as corrosion [32], and that although both influence the behavior of the material, they should not be taken into account simultaneously, given that while residual stress acts mainly in the beginning of the operational life of the ship, corrosion is negligible in the beginning and increases with time. To avoid lateral-torsional bending in ships, the influence of the residual stress over the ultimate bending moment is quantified in Equation (9):

$$M_{ur}=\left(1-0.3{\overline{\sigma }}_r\right)M_u$$

(9)

where ${\overline{\sigma }}_r$ is the residual stress of the panels in compression and $M_{ur}$ and $M_u$ are, respectively, the ultimate bending moment with and without residual stress.

Innovative mechanical reinforcement techniques that go beyond classic metallic structures have also been studied, such as rubber coated ships [34], which show the ability of damping the high-frequency response from pressure waves, helping to keep the integrity of the vessel to crashes and close-range underwater explosions that may cause elastic and even plastic deformations to the side structures. However, this coating is not so effective regarding low-frequency whipping motions caused by gas bubbles.

2.3. Material Response

The first important material feature that comes to mind when one needs a vessel design to withstand impact and wear solicitations is high mechanical resistance thresholds [35,36]. However, in addition to being tough, building materials must also be as lightweight as possible to consume less fuel, leave room for more load bearing capacity and to provide the ability to quickly overcome inertia when maneuvering, especially in the case of military vessels, to make them difficult targets.

The study of metallic phases and alloys figure as an important resource to understand and improve the resistance of military naval structures to underwater explosions [37], especially in terms of the high inherent density attributed to metals. Latourte et al. [2] compared experimentally and numerically high-strength martensitic and austenitic alloys, investigating their deformation and fracture characteristics, where the first was designed to present higher fracture toughness and maximum strength, while the latter has a better uniform ductility, avoiding premature necking and consequent localized failure. The interesting outcome is that the martensitic steel was found out to sustain underwater fluid-structure impacts better, where the force magnitude needed to infer failure is 25% higher than for austenitic steel (although the austenitic performance can be improved by tempering) [2].

Besides traditional monolithic metals, the coupling of different materials in sandwich layouts have presented interesting features relevant for ship design, like the flexibility in the core constitution, which can lead to a stiffer overall material in the through-the-thickness direction, or can soften and dampen impact and decrease the force transmission to the supporting structure [4]. Furthermore, sandwich panels present an enhanced performance against underwater fluid shocks when compared to single-constituent elements of equal weight [38,39].

Fan et al. [40] investigated the response of sandwich panels with aluminum sheet skin and a honeycomb core to underwater blasts, comparing the results with monolithic plates of equivalent mass, proving a better performance of the sandwich both in terms of deformation resistance and secondary pressure wave intensity. The first advantage can be furthermore improved by increasing the equivalent thickness of the composite. As for composites with polymeric matrix, Gong et al. [19] studied the transient behavior of glass-epoxy composites to UNDEX, showing that although a composite hull might amortize the effects of bubble impacts better than steel, it becomes more susceptible to global mechanical effects, not only local ones as would be the case for metals.

In an experimental underwater explosion test, it is also possible to identify the energy absorption capacity of a certain material by analyzing its deflection, once it has been proved that during plastic regime deformation, the ability of the material to absorb energy is proportional to the square of its deflection [23]. Another important aspect is the specimen geometry to be adopted in a bulge test, given that parabolic and spherical shapes absorb practically the same energy, while conical and hyperbolic absorb only half of that considering the same bulge depth [27]. Regarding experimental tests, it is possible to predict the maximum von Mises stress for both thin circular (${\sigma }_c$) [41] and rectangular plates (${\sigma }_r$) [42], as shown respectively in Equations (10) and (11), considering that it is an explosion with low intensity (for specimens farther from the charge more than 10 times its explosion radius) and that the deformation is realized within the elastic regime.

$${\sigma }_c=\sqrt{\frac{6E{\rho }_P{P}_m^2{\chi }^{2/\left(1-x\right)}}{{\rho }^2{c}^2\left(1-\upsilon \right)}}$$

(10)

$${\sigma }_r=0.867\sqrt{\frac{14E{\rho }_P{P}_m^2{\chi }^{2/\left(1-x\right)}}{{\rho }^2{c}^2}}$$

(11)

where E is the elastic modulus and ${\rho }_P$ is the density of the plate material, $\rho$ is the density of water, c is the velocity of sound in the water medium, P_m is the peak pressure, $\chi$ is the angle of incidence to the plate, and $\upsilon$ is the Poisson’s ratio. In order to indicate the level of damage caused by a shock, the shock factor index (SF) is used as displayed in Equation (12) [42]. To verify whether the shock energy can cause yielding to a thin air backed plate, the yield shock factor (SF_y) would have to reach the limiting value expressed by Equation (13) [42].

$$SF=0.445\frac{\sqrt{W}}{S}$$

(12)

$$S{F}_y=Y\frac{1}{\sqrt{\eta }}{\sigma }_y\sqrt{t},$$

(13)

where W is the TNT equivalent of charge quantity (kg), S is the standoff from the explosion center, t is the thickness of the plate, ${\sigma }_y$ is the yield stress of the plate material, Y is the yield factor equivalent to 2.212 × 10⁻⁹ for circular plates (Y_c), and 1.997 × 10⁻⁹ for rectangular ones (Y_r), and $\eta$ stands for the coupling factor, which is a function of the incidence angle as shown in Equation (14).

$$\eta =4{\chi }^{\left(1+\chi \right)/\left(1-\chi \right)}$$

(14)

It is also possible to determine the yielding shock factor (SF_y) as a function of the necessary charge quantity W (Equation (15)) or the standoff S (Equation (16)) and a time constant $\theta$ (ms).

$$W={\left(\frac{SF}{0.445}\right)}^{0.5946}{\left(\frac{\theta }{96.6\times {10}^{-6}}\right)}^{2.7026}$$

(15)

$$S={\left(\frac{0.445}{SF}\right)}^{0.7027}{\left(\frac{\theta }{96.6\times {10}^{-6}}\right)}^{1.3513}$$

(16)

As for the outcome material damage, it mainly depends on the material, boundary conditions, and loading rate considered. Ming et al. [18] have shown through experiments the particular case of clamped flat metal plates, in which damage can be visually categorized in three subsequent steps. The formation of a localized protuberance in the point of maximum stress immediately before the material fracture takes place (“bulging”), leads to detachment of a small metal fragment giving place to a hole (“discing”), whose size is determined by the plate thickness, charge weight, and yield stress [26,43] that continues to grow until the material dissipates this energy through the propagation of radial and equally spaced cracks, deforming the metal in petals (“petaling”), as detailed in Figure 3. A particular aspect noticed by this study is that while the hole size is sensitive to the peak pressure, the deflection near the whole is strongly affected by the impulse of the explosion.

The hole radius R found in the discing step was calculated by Rajendran et al. [43], as described in Equation (17):

$$R=\sqrt{2\beta W{E}_i/\pi t{\sigma }_y{\epsilon }_f}$$

(17)

where $\beta$ is the ratio of effective work, E_i is the inner energy per unit mass of the explosive, and ${\epsilon }_f$ is the fracture strain. Before that, while the material is in the bulging stage, its round-shaped deformation can be calculated by Equation (18) [43], also represented by Figure 4, where R_S and R_C are, respectively, the radii of the sphere and the circle; t is the thickness of the plate; and D_S is the depth of the instantaneous indentation over the plate:

$$R_S=\frac{R_C{}^2}{2D_S}+\frac{D_S}{2}$$

(18)

Once the mechanical response of metals to UNDEX is highly dependent on the strain rate of the solicitation, it is vital to illustrate that this effect must be taken into account in calculations for non-idealistic results. Rearranging the stress-strain relations found in [44], this behavior of metallic plates can be summarized in Figure 5:

It is possible to see that metallic plates tend to absorb less energy and behave in a more brittle manner when the strain rate effects are acknowledged in stress-strain calculations. This is made by considering a strain rate factor n as a constant associated to stress. The factor n defines the property demeaning level of metallic plates, being directly proportional to the impulse load and inversely proportional to the thickness and yielding stress of the plate. This relation is valid for both circular and rectangular plates [27]. As demonstrated by Jones [44], the influence of n can be calculated by rearranging stress-strain as a relation between the deflection at the center of the plate ($\delta$), the thickness (t), and a modified damage parameter (${\varphi }_M$) that embraces both loading and plate dimensions, through Equations (19) and (20) for circular ($n_c$) and rectangular plates ($n_r$), respectively.

$${\left(\frac{\delta }{t}\right)}_c=0.817\frac{{\varphi }_M}{\sqrt{n_c}}$$

(19)

$${\left(\frac{\delta }{t}\right)}_r=0.95\left[{\left(1+0.6637\frac{{\varphi }_M{}^2}{n_r}\right)}^{1/2}-1\right]$$

(20)

In those relations, the strain rate factor can be calculated as shown by Equations (21) and (22), where D and q are material parameters, I is the impulse, t is the thickness of the plate, ${\sigma }_y$ is the yielding stress, ${\rho }_p$ is the material density, and R is the radius of the circular hole bored.

$$n_c=1+{\left(\frac{I^2}{3{\rho }_p^2{t}^{2}DR}{\left(\frac{{\rho }_p}{3{\sigma }_y}\right)}^{1/2}\right)}^{1/q}$$

(21)

$$n_r=1+0.0357{\left(\frac{I^2}{t^2\sqrt{{\sigma }_y}}\right)}^{1/q}$$

(22)

As seen above, small-scale UNDEX experimental tests are routinely carried out by many authors to allow the perception of how materials would behave in operational conditions. However, especially in the case of UNDEX, where high-magnitude forces directly affect materials and structures, involving a complex array of fracture mechanisms, the ability to scale those experiments to the larger and thicker panels that are actually applied in ship building and are subjected to higher stresses is relevant.

However, the intrinsic difficulty to reproduce these experiments in a larger-scale that is closer to reality is clear: high cost, safety concerns, logistics, access to high amounts of explosives, and confidentiality issues make these real-sized tests practically exclusive to national Defense Departments. Therefore, as stated by Cui et al. [24], small-charge experiments remain the best way to study UNDEX phenomena.

Jacob et al. [45] analyzed the scaling effect in regular test scale experiments. Regarding the relation between deflection and impulse, it is reported that the exposed plate area does not influence the deformation, but thickness does. Also, the charge diameter is relevant to determine the plate response. As for the effects of charge height and diameter, the impulse increases as either height or diameter increase independently of the thickness of the plate. This study provides some graphic extrapolation of data to predict the outcome of not-so-practicable tests, but it advises the usage of numerical simulations as a reliable tool to reproduce even operational-scale conditions, meaning that there is no point in performing extensive, dangerous, and costly large-scale experiments.

3. Numerical Simulations

Diffused in most engineering fields, numerical simulations represent a very useful tool for faithfully predicting the behavior of materials to various mechanical solicitations, constituting a valid and time-saving approach in designing structures and materials to avoid rework and ensure safety at an optimal cost. Thus, some of the fundamental simulation aspects regarding fluid description, structural, and materials responses to UNDEX will be analyzed.

3.1. Fluid Description Algorithms

In numerical simulations, the fluid dynamics can be described by two main approaches, the Eulerian and the Lagrangian. Both can be used to reproduce free fluid surface and deforming wall boundaries, but present some differences. In the first, the properties of the flow are given for defined spatial coordinates as a function of time, by observing how the properties of a certain point in space change due to the flow that passes through it. On the other hand, the Lagrangian method depicts the flow as a large number of individual particles whose motion is described, following them and tracking their property variations over time. These different methods are reproduced within numerical simulations by either setting a fixed grid in space (Eulerian) or nodes that move according to the velocity field in a meshless setup (Lagrangian).

Because of the intrinsic high computational cost to simulate fluid interfaces proper of fixed-grid Eulerian-based models, and the high distortions for violent fluid-structure interactions of the Lagrange method [46], both typical characteristics of UNDEX, the combination of the two into the Arbitrary Lagrangian Eulerian approach (ALE) is quite common in this field of study. It considers a moving grid to provide enhanced interfacial precision and to optimize the usage of computational resources for the calculations [47]. The high CPU cost to simulate fixed mesh problems has been addressed by several researchers. Wang et al. [48] were able to cut down this effect by defining optimal ratios between the radius of the charge and the side length of the mesh elements through a dimensionless variable, whose recommended value was found to range between 3 and 6 for most cases.

Even with its limitations for complex simulations, the description of compressible flow by Eulerian algorithms is still popular in literature. Liu et al. [49] developed a continuous UNDEX simulation model that embraces both shock wave generation and bubble motion stages consecutively to overcome the usual two-step routine. Through this method, the high-pressure bubbles become smaller and weaker than the real ones with errors up to 10%, most likely related to intense interaction of the reflection wave with the expanding bubble in the proximity of the wall. Hu et al. [50] managed to obtain satisfactory values for compressible multi-fluid flows prevenient from UNDEX (gas and liquid water) through a sharp-interface (i.e., considering discontinuous material properties across the interface). Accordingly, Ma et al. [51] successfully presented an extended version of a known set of one-dimensional equations for compressible multiphase fluids to two and three dimensions.

Purely Lagrangian techniques have also been effectively used to numerically reproduce UNDEX, such as in Ming et al. [18], where a faithful fracture pattern on metallic panels was achieved, and Zhang et al. [52], where the penetration of metal jet on steel plates was accurate, thus, proving the interfacial efficiency of this method. Aiming to supply the faults and couple the benefits of both methods, ALE has also been applied by authors like Jafarian et al. [53] to allow the simulation of UNDEX compressible flow and cavitation (formation, development, and collapse of bubbles), using a single fluid by Petrov et al. [54] to take into account cavitation and rarefaction waves that propagate through the liquid, and by Wardlaw [55] to validate simulations at one, two, and three dimensions and with a variable number of moving boundaries.

3.2. Numerical Modelling

The Eulerian grid is applicable to solving methods like Finite Element Method (FEM) [2,28,30,31,36] and Finite Volume Method (FVM) [48,51,55,56], since both calculate the values of the fluid properties at discrete places on a meshed geometry. In FEM, these discretized spatial units are defined as elements and have constant properties, being mostly used to model solid ship structures and materials, although fluid domains can also be formulated [57]. Prusty et al. [31] and Gupta et al. [28] are good examples of the precision with which FEM models can reproduce the response of metallic stiffened panels to UNDEX. Also, FEM constitutes the best technique through which the hull girder strength can be assessed by a progressive collapse analysis and its consequent non-linearity [56], as Chung et al. [58] demonstrated by evaluating the structural response of a catamaran, even taking into account the strain rate effects of the material. Not confined to monolithic structures, Tillbrook et al. [37] demonstrated that beams made of sandwich material compositions can have their dynamic response to UNDEX adequately described by FEM, which is a recognized standard approach to model ship structures by naval classification agencies and the navy [59].

In turn, FVM supports unstructured meshes of discretized units named cells, through which matter is allowed to flow following conservative physical laws. Unlike in elements where boundary conditions may be applied in surfaces or nodes, in FVM these are directly applied in the control volume within the cell, thus not influencing neighbor matter units. Wang et al. [48] make use of FVM to define both water and explosive charge materials, allowing them to flow over the mesh cells, and Ma et al. [51] applied the FVM to model one, two, and three-dimensional problems with multiple fluids efficiently, proving the validity of this method. Although FVM is not a popular choice to model either structures or fluids, military-grade validations of its usage have been carried by studies such as Wardlaw’s [55].

Other techniques with interesting resources, such as Boundary Element Method (BEM) [19,60,61,62] (also used for propeller hydrodynamic calculations [63]) and Smooth Particle Hydrodynamics (SPH) [18,52] have proven to be efficient, having their own intrinsic advantages.

BEM, also known as BIM (Boundary Integral Method), discretizes only the edges of the control volumes formed by a grid. Due to this grid simplicity, BEM allows the calculation of more complex integral equations with a reduced computational power, potentially resulting in highly accurate outputs. It is particularly popular to model problems involving bubble impact, as it provides efficient algorithms to represent its transient shape and pressure from the moment of the explosion until the collapse. Wang et al. [60] proposed numerical modifications to include not only the compressibility over the bubble to simulate its damped oscillation, but also the effect of oscillation cycles subsequent to the first one before it breaks down into smaller bubbles. Zhang et al. [61] introduced improvements to the regular 3D BEM model, enhancing the accuracy and stability at large deformations of the toroidal bubble phase. Li et al. [62] studied the pressure field caused by a collapsing bubble by the auxiliary function method, finding out that its dynamic total pressure can be decomposed and quantified in two parts: one regarding the pressure gradient between the gas within the bubble and surrounding liquid, and another related to its motion. Coupled approaches have also been explored for hydrodynamic problems, such as RANS-BEM (Reynolds-averaged Navier-Stokes) [20] and BEM-FEM. For example, the latter, as studied by Gong et al. [19], allowed the simultaneous simulation of the physics of bubble growth, contraction, and collapse, and the consequent behavior of the hull. The shape transformation of the bubble from its generation until impact, as generically represented in these studies, is displayed in Figure 6.

However, the discretization of a domain in elements might lead to difficulties to deal with free-surfaces, deformable boundaries, moving interfaces, and large fluid deformations, and assuring a satisfactory mesh quality may often be a tricky and time-consuming process. Within this context, Lagrangian based methods, such as SPH, present an interesting alternative to overcome the barriers imposed by a grid system [64]. SPH is a robust method for intricate hydrodynamic problems [46] that represents a set of particles with material properties that interact with each other in the flow. This method is attributed to numerical simulations, in which the analysis deals with large deformations, given that it has a meshless nature. Zhang et al. [52] applied SPH to simulate a charge detonation, the metal-jet formation, and the penetration on a steel plate, agreeing with experimental results. Although impacts on the air-water interface constitute a limitation of the SPH-FEM coupling because of its instabilities in modelling particular water physics [65], it is a valid and precise approach for UNDEX. Ming et al. [18] demonstrated through the comparison of experiments using SPH and SPH-FEM methods that both are accurate to reproduce UNDEX. De Vuyst et al. [66] have actually proved through coupled SPH-FEM simulations the interesting concept that considering the case of having either one or multiple underwater explosions for a same charge quantity, the latter is less harmful to the steel than the first, as it generally results in smaller final deformation and might avoid cracks that could be formed by a single blast.

Given the aforementioned considerable amount of input on numerical modelling techniques, a similar approach to Hirdaris et al. [67] is reflected in

Table 1. Taxonomy of hydrodynamic modelling methods as a function of their discretization dependency.

Level	Method	Key Features	References
1	FEM	Mesh-dependent Conventionally used to simulate the behavior of materials and ship structures Can be used to model nonlinear impact Impractical for violent flows If large deformations are considered, the mesh must be too fine and computations become slow	Latourte et al. [2] Gupta et al. [28] Hung et al. [30] Prusty et al. [31] Tillbrook et al. [37] Holtmann et al. [59] Wang et al. [68]
2	FVM	Mesh dependent, but allows unstructured meshes Cells obey conservative laws Boundary conditions are applied noninvasively within the volume cell	Wang et al. [48] Ma et al. [51] Wardlaw [55] Rigo et al. [56]
3	BEM	It is highly accurate due to the simplicity of the grid that allows the adoption of complex integral equations Velocity potential in fluid domain represented by a distribution of sources over the mean wet body surface Suitable for bubble formation, growth and collapse Precise for problems involving stress concentrations May give unreliable values for added mass and damping at irregular shock frequencies	O’Daniel et al. [16] Gong et al. [19] Wang [60] Zhang et al. [61] Li et al. [62] DNV [69]
4	SPH	Meshless The Lagrangian nature tracks the mass of material particles and difficulties boundary conditions setup Advection can be calculated Adequate for free-surface, interfacial and violent flow Exact and simultaneous conservation of mass, momentum, angular momentum, energy and entropy Computationally expensive due to higher number of neighbor elements Easy to parallelize. So, even if it is costlier than other methods, it has the potential to be faster	Ming et al. [18] Zhang et al. [46] Zhang et al. [52] Liu et al. [64] Hughes et al. [65] De Vuyst et al. [66] Jones et al. [70]

to summarize the key features of each method addressed in this subsection, aiming to provide a clearer visualization of their particularities. The taxonomy for classification adopted is based on the level of mesh dependency they present to yield accurate outputs, as some methods always include discretization, others partially, and others none.

It is vital to state that a taxonomic descriptive table was preferred to point out the peculiarities of each method, because it is not possible to indicate in a simple manner which method is more appropriate in each occasion; after all, the method choice is subjective for the user, and takes into account their affinity with specific software, the CPU power available, the level of precision requested, and the boundary conditions of the problem to be analyzed. Hence, the same method can present diverse levels of complexity depending on the situation. Besides, through the scientific review hereby carried, it was possible to notice that all methods are capable of providing valid results. Despite this, the survey showed a literature inclination to FEM to reproduce vessel structures, and to BEM and SPH to model the fluids, be they gas (prevenient from bubbles), liquid, or both.

Overcoming the academic sphere, naval classification agencies and official military reports also support these methods, confirming their legitimacy to simulate UNDEX. Holtmann from Det Norske Veritas [59] advises FEM to evaluate the shock resistance, accelerations, strains, and stresses of hull girders for providing a fast and efficient approach. Wang et al. from American Bureau of Shipping make use of nonlinear FEM to design ship structures for ice loads [68]. Warlaw [55] from the U.S. Navy considers the division of the UNDEX simulation domain into FVM cells applied to problems with different dimensions and border conditions. Det Norske Veritas [69] sees BEM as the most common method to solve potential flow problems, whereas O’Daniel et al. [16] from the U.S. Army faces this method as the most suitable to simulate the behavior of bubbles. Also, Jones et al. [70] from the Department of Defense of Australia recognizes SPH as appropriate for modeling several fluid problems of particular military interest.

3.3. Constitutive Models

The techniques used to model impact solicitations differ depending on the type of material considered [71,72,73], and particular sets of equations must be adopted to form a software-embedded constitutive model, which generally requires calibration of variables. Often, given the wide application of metals in shipbuilding, the models used for UNDEX make use of relations developed for metallic materials, which are generally highly dependent on strain rate effects [74].

Accordingly, ship design relies on such resource to properly characterize the vessel itself and all its structural components to efforts intrinsic of navigation, as well as to external input forces such as UNDEX. A variety of studies deal with different strategies to model the aforementioned problem, exhibiting the versatility with which it can be done by using diverse approaches and modeling techniques.

Characterizing a common ground among most authors, the application of the Johnson-Cook constitutive model [75] for the ship hull or single panels is widely utilized, coupled with equations of state in numerical codes, allowing its usage for UNDEX loading. Furthermore, it accounts for equivalent plastic strain, strain rate, and the influence of temperature, while requiring only five parameters that can be found by means of experimental tests, as described in Equation (23):

$$\sigma =\left(A+B{\epsilon }_{eff}^n\right)\left(1+Cln{\epsilon }_{ref}\right)\left(1-T_h^m\right)$$

(23)

where ${\epsilon }_{eff}$ is the effective plastic strain, ${\epsilon }_{ref}$ is the same strain but at a reference rate ${\epsilon }_0=1s^{-1}$ and $T_h$ is the homologous temperature calculated by the quotient of the difference between the material temperature and the room temperature by the difference between the melting temperature and the room temperature. As for the five parameters, namely A, B, C, m, and n, the first two and the latter are obtained from tension tests, while C and m are defined through split Hopkinson pressure bar tests. This model was used by Kong et al. [76], which numerically reproduced a blast load detonated inside a multi-layer protective and stiffened naval structure, demonstrating in this case that stiffened plates are severely damaged by fragments which penetrate, causing diverse spots of crack initiation culminating into crack propagation. On relatable UNDEX studies [18,28], stiffened plates were also modeled, making use of the Johnson-Cook criterion.

Another relevant constitutive model used in this field, although secondary compared to Johnson-Cook’s for being less accurate to represent UNDEX, is Gurson’s [77]. It can be calibrated by the realization of mechanical tests [2] and describes stress flow and rupture of materials depending on void growth. A description of the damage growth rate ($\dot{f}$) as a function of the void volume fraction ($f$) and the plastic strain rate ($D^P$) in a pure shear stress case (i.e., in the kk direction) follows (Equation (24)):

$$\dot{f}=\left(1-f\right)D_{kk}^P$$

(24)

In this case, the limitation assumed by this model is that no damage is predicted when strain is under zero mean stress, depending exclusively on the formation of voids. Aiming to address this issue and include the effect of shear-induced fracture also in the absence of voids, a modification of this model was introduced by Nahshon et al. [78], which adds the numerical constant $k_{\omega }$ to define the rate of damage evolution in a shear-predominant stress state. The evolution of the model is thus presented by the incorporation of a second term on Equation (24) resulting in Equation (25):

$$\dot{f}=\left(1-f\right)D_{kk}^P+k_{\omega }f\omega \left(\sigma \right)\frac{s_{ij}{D}_{ij}^P}{{\sigma }_e}$$

(25)

where $\omega \left(\sigma \right)$ describes the relation between a third invariant of stress with effective stress [78], $s_{ij}$ is the stress deviator $s_{ij}={\sigma }_{ij}-1/3{\sigma }_{kk}{\delta }_{ij}$, and ${\sigma }_e$ stands for the effective stress. This enhanced approach was validated by Xue et al. [79], where the calibration of the damage parameters inherent of this constitutive model were able to reproduce tension and shear induced failures on steel coupons. Latourte et al. [2] used it to numerically describe the behavior of high performance steel alloys to UNDEX.

As for sandwich panels, studies such as the ones from Tilbrook et al. [37,80], considering anisotropic foam cores and elastically ideal skin sheets, provide particular constitutive models allying the mechanics of core compression with the bending response of sandwich beams. These two impulsive load responses can be described in four ways, depending on whether they are analyzed separately (if the accelerations of the sheets differ) or together, and if the core densification is total or partial.

4. Discussion

Addressing the mechanics embraced by underwater explosions, the present survey underlines essential mathematical relations from the hydrodynamics of subsea explosions, to the propagation of high pressure waves and bubbles, to the behavior shown by materials and structures commonly used in shipbuilding. The validity of these methods has been checked through the high correlation of purely analytical formulations with several experimental works, although some discrepancies among authors have also been identified, such as the constants used for determining peak pressure as a function of the charge load.

Likewise, the variety of computational numerical resources applied in this area have permitted the constitution of advanced research studies that are able to mirror UNDEX, conditions considering complex wave force transmission mechanisms and progressively accurate material responses in parallel, thus ratifying the forefront role played by virtual simulations in related engineering problems.

Moreover, the confrontation of the late efforts in improving and perfecting the precise reproducibility of constitutive models to real situations, with the reported advanced material characterization and hydrodynamic studies, show that this traditional but ongoing research field has kept its evolving pace and still presents a remarkable development potential.

Among the possible future directions for this area of study, it is hereby encouraged that a comparative study considering different numerical simulation techniques (such as FVM, BEM, and SPH) is carried, adopting the same boundary conditions, aiming to narrow the current subjectivity in selecting the modelling method to a more grounded approach, crossing desired output accuracy with computational power available. Regarding the analytical side, it would be considerably relevant if new experimental studies could confront the response of materials for hull panels to UNDEX to those in the sometimes-conflicting studies available in the literature. In this case, special attention should be given to fiber reinforced composite materials for their growing importance in the nautical industry over the past years, considering their behavior to underwater explosions has not yet been properly explored.