2.1.1 Shock pressure

The shock wave generated by a UNDEX is a high-magnitude, short-duration compressive pressure wave that travels faster than the wave speed of the surrounding medium (Reference ShinShin 2004). This phenomenon occurs instantaneously with the combustion and chemical decomposition of explosives, which produce gaseous byproducts at extremely high density and temperature (Reference ColeCole 1948). The surrounding water is rapidly accelerated, causing a sharp, discontinuous increase in pressure, which then decays exponentially as the gaseous byproducts expand into a bubble. The shock wave radiates outward from the source of the explosion in a spherical wavefront, with its amplitude diminishing due to spherical spreading. As the expansion continues, the shock front velocity decelerates, eventually aligning with the water's acoustic speed (Reference ShinShin 2004).

The characteristics of a UNDEX shock pressure wave, $P(t)$, can be captured by (2.1) (Reference SwisdakSwisdak 1978; Reference ShinShin 2004). The amplitude, ${P_{max}}$, and decay constant, $\theta $, can also be modelled using empirical (2.2) and (2.3). Although empirical relations, these equations mirror the dimensionless equations and solutions that describe strong shock waves in the air (Reference RamamurthiRamamurthi 2021), using explosive energy instead of weight and distance instead of hydrostatic pressure (since these parameters are proportional)

(2.1)\begin{gather}P(t) = {P_{max}}\,{\textrm{e}^{ - t/\theta }}, \end{gather}

(2.2)\begin{gather}{P_{\max }} = {K_P}{\left( {\frac{{{W^{1/3}}}}{R}} \right)^{{\alpha _P}}},\end{gather}

(2.3)\begin{gather}\theta = {K_\theta }{W^{1/3}}{\left( {\frac{{{W^{1/3}}}}{R}} \right)^{{\alpha _\theta }}}.\end{gather}

In these empirical equations, ${K_P}$, ${K_\theta }$, ${\alpha _P}$ and ${\alpha _\theta }$ are empirical constants that are specific to the explosive material, W is the mass of the explosive and R is the distance from the explosion to the point of interest. These equations illustrate the intensity of pressure near the explosion, highlighting that pressures can become exceedingly high as one approaches the explosion source. Similar equations can also be formulated to solve for energy and impulse with the fitting parameters for common explosive materials such as TNT (Reference SwisdakSwisdak 1978). The limitations of using these empirical equations lie in the parameters being tied to a specific explosive material and the range of the data used to determine the fitting parameters.

2.2.1 Bubble dynamics overview

The gas bubbles from UNDEX events are comprised of explosive byproducts and rapidly phase-shifting water vapour engendered by the high temperatures of an exothermic reaction. Historically conceptualized as an oscillating gas globe or sphere in classical literature (Reference ColeCole 1948), this bubble undergoes rapid radial expansion propelled by the exceedingly high pressure contained within. As it expands, the bubble's boundary is propelled outward spherically until it is curtailed by the external fluid pressure, eventually overcoming the internal pressure, causing the boundary velocity to decline and ultimately arrest at the bubble's maximum size.

Subsequently, the bubble enters a collapse phase, signified by the reversal of boundary velocity. This collapse releases a rarefaction wave, or an under-pressure wave, with sufficient force to potentially deform compliant structures. The collapse maintains a near-spherical shape in an unobstructed free-field setting. This infinite fluid domain is absent of any close structures, surfaces or boundaries to influence the flow induced by the bubble's oscillations. The bubble then contracts to a small but finite minimum size (determined by the volume of its condensable and non-condensable gases) before it grows outward again. When the bubble reaches its minimum size, a second shock wave, known as a bubble pulse, is released. Although less intense, this wave has a longer duration than the initial shock wave and results from the abrupt shift in boundary velocity at minimum size. Each successive expansion and collapse cycle diminishes the bubble's maximum radius and the strength of subsequent bubble pulses (Reference ColeCole 1948; Reference ShinShin 2004).

Free-field bubbles also exhibit an upward migration trend driven by buoyant forces from the pressure gradient along the bubble boundary. Larger bubbles experience more pressure change across this boundary, leading to a more pronounced ascension. Figure 1 offers a graphical representation of these general free-field bubble dynamics and the associated pressure changes, as shown by Reference ShinShin (2004).

Figure 1.

Gas bubble growth, migration and bubble pulse (Reference ShinShin 2004).

2.2.2 Free-field bubble empirical modelling

Reference ColeCole (1948) discusses empirical formulas that delineate the peak pressure and duration of UNDEX shock waves and a free-field bubble's maximum radius and period across varying hydrostatic pressures and explosive charge masses. These formulas have found utility in characterizing various explosive materials (Reference SwisdakSwisdak 1978; Reference ShinShin 2004) and accurately predicting their behaviour. Variations of these formulas, adapted to use imperial measurements, are presented in (2.4) and (2.5) (given in the SI-MKS unit system):

(2.4)\begin{gather}{R_{max}} = {K_R}{\left( {\frac{W}{{D + 10}}} \right)^{1/3}}, \end{gather}

(2.5)\begin{gather}T = {K_T}\frac{{{W^{1/3}}}}{{{{(D + 10)}^{5/6}}}}.\end{gather}

In these empirical equations, ${R_{max}}$ signifies the bubble's maximum radius in metres, T is the period of the initial bubble oscillation in seconds, W is the explosive charge mass in kilograms, D represents the depth in metres and ${K_R}$ and ${K_T}$ are empirical coefficients that correlate with the specific type of explosive in use (in SI-MKS units, ${K_R}$ is in ${m^{4/3}}/k{g^{1/3}}$ and ${K_T}$ is in $s({m^{5/6}}/k{g^{1/3}})$).

2.2.3 The energy of underwater shock and bubble

It is often necessary to model the shock and bubble pressure signatures in terms of their energy output. Measuring the energy output of an explosion event is critical because it helps quantify the loading to structures or systems impacted by the blast. For scientific research, such measurements can lead to a better understanding of explosive dynamics and developing predictive models. The energy terms associated with the total energy in an explosion event typically include the kinetic energy, the thermal energy released, the energy transmitted as shock waves through the medium and the potential energy associated with the gas bubble formation.

Previous work by Reference Arons and YennieArons & Yennie (1948) looked into the energy partition during UNDEX events. These energy components include the irreversible energy flux contributing to the acoustic radiation as the shockwave travels away from the explosion source and the reversible after flow energy stored in the water due to fluid-particle movement during bubble formation and oscillations. The energy transferred to the surrounding medium (water) can be partitioned into various forms, but the work by Reference Arons and YennieArons & Yennie (1948) focuses on the radiated energy and the energy related to the after flow process associated with the bubble dynamics post-explosion, which are the main components of the total energy in the system. The total radiated and reversible energy, E, for a spherical wave, can be found as a function of the pressure time history, P, and fluid properties, as shown in (2.6). Some energy terms are neglected in this equation due to their relatively small contribution to the total energy, including the kinetic energy imparted to the water and the change in internal energy in the mass of water as the shockwave propagates. Furthermore, other energy components not explicitly considered in (2.6) but discussed in the work by Reference Arons and YennieArons & Yennie (1948) include the total irreversible thermal, viscous and turbulent energy losses in water. These energy loss terms become critical for calculating the after flow energy after each subsequent bubble oscillation cycle

(2.6)\begin{equation}E = \frac{{4\pi {R^2}}}{{{\rho _o}{c_o}}}\int {{P^2}\,\textrm{d}t} + \frac{{2\pi R}}{{{\rho _o}}}{I^2}.\end{equation}

In this relationship, the first term represents radiated energy, and the second represents reversible energy. The R term is the radius from the explosive source to the spherical shock front, ${\rho _o}$ the fluid density, ${c_o}$ is the fluid's acoustic wave speed and I is the specific impulse of the pressure wave as defined by (2.7). The limits of integration for (2.6) and (2.7) are relative to when the shock pressure begins and ends

(2.7)\begin{equation}I = \int {P\,\textrm{d}t} .\end{equation}

2.2.4 Near-field bubble dynamics

The interaction of UNDEX bubbles with nearby boundaries in their vicinity, termed ‘near-field’ interactions, carries the potential for severe structural loading. This is attributed to the formation of a re-entrant jet during the bubble's collapse in proximity to structures, which can have a destructive effect akin to a water hammer. This phenomenon occurs as a high-velocity jet of water thrusts into the bubble from the side opposite to the boundary, enabled by the high compliance of the gas within the bubble, which allows the jet to accelerate substantially, as illustrated in figure 2. The result is an intensely localized pressure that impacts a focused area of the structure (Reference Javier, Galuska, Papa, LeBlanc, Matos and ShuklaJavier et al. 2020a). These events are incredibly chaotic, and no reliable analytical models can be used for prediction purposes. Simulations of this phenomenon have been explored computationally, as illustrated in figure 3, and are discussed in more detail in § 5.

Figure 2.

Bubble collapse near a rigid barrier, resulting in jet formation.

Figure 3.

Experimental and simulated maximum and minimum bubble volumes for standoff distances of (a) γ = 1.5, (b) γ = 1.0 and (c) γ = 0.75 (Reference Javier, Galuska, Papa, LeBlanc, Matos and ShuklaJavier et al. 2020a).

Collapses of underwater bubbles and cavities near structures have historically been shown to cause significant deformation and even corrosion. While ‘near field’ lacks a universally accepted definition, it typically refers to a scenario where a bubble's proximity to a boundary – such as a structure or free surface – is close enough to cause consequential fluid–structure interactions or influence the bubble dynamics. In these near-field cases, the structure can inhibit bubble expansion by constraining fluid motion on one side, leading to an asymmetric bubble growth cycle and, potentially, to the formation of a bubble jet. A range of factors influences the trajectory of this jet: the distance between the explosive charge and the structure, known as the standoff distance; the effects of buoyant forces; the structure's compliance, shape and boundary conditions; and the presence and positioning of additional structural elements or free surfaces (Reference Javier, Galuska, Papa, LeBlanc, Matos and ShuklaJavier et al. 2020a; Reference Leger, Matos, Shukla and JavierLeger et al. 2023).

2.2.5 Near-field bubble jetting

The study of gas bubble dynamics traces back to Rayleigh's analytical model of a spherical bubble collapsing in a free-field environment, established in 1917 (Reference RayleighRayleigh 1917). Research on near-field rigid boundary jetting expanded significantly after Reference Konfeld and SuvorovKonfeld & Suvorov's (1944) proposition that bubble jetting was the principal cause of cavitation damage. Benjamin and Ellis's experiments supported this hypothesis of a liquid jet forming during bubble collapse (Reference Benjamin and EllisBenjamin & Ellis 1966), which produced high-speed photographic evidence of bubbles generated by electrical sparks. They also suggested applying the analytical Kelvin impulse concept to evaluate the fluid momentum of an asymmetrical cavitation bubble. An experimental investigation by Reference Naudé and EllisNaudé & Ellis (1960) honed in on the loads caused by bubble jetting and the interaction of the bubble boundary, employing photoelastic materials to visualize the stress conditions induced on the boundary.

Reference Plesset and ChapmanPlesset & Chapman (1970) then advanced a numerical method predicated on a Rayleigh bubble to forecast these loads and utilized the water hammer principle to estimate the likely jet velocities and pressures. Reference Lauterborn and BolleLauterborn & Bolle (1975) validated this computational approach by experimentally determining jet tip velocities through high-speed photography.

Pioneering studies of gas bubbles in proximity to other types of boundaries, including free surfaces and compliant boundaries, were conducted by Reference GibsonGibson (1968). His experiments analysed the differences in jet behaviour and bubble–structure interactions at various explosive standoff distances, demonstrating that greater standoff distances diminish the potential for jetting and interaction with the structure.

The implications of standoff distances on bubble jetting were further examined in a small-scale UNDEX study by Reference Javier, Galuska, Papa, LeBlanc, Matos and ShuklaJavier et al. (2020a). Their research delved into the pressures experienced by specimens due to bubble jetting against the surface of a rigid boundary. Both experimental approaches and numerical simulations were employed to capture the phenomenon. High-speed photography was utilized to record the bubble jetting, and the results were compared with numerical solutions for three dimensionless standoffs. The findings are shown in figure 3, which illustrates that reduced standoff distances (represented by the dimensionless parameter $\gamma = \textrm{Standoff}/{R_{max}}$) lead to more pronounced jetting towards the rigid boundary.

2.2.6 Damage of shock waves and explosions

Underwater explosion shock waves are recognized for their potential to damage structures, which is relevant in far-field and near-field scenarios. Reference KennardKennard (1944) conducted a seminal analytical investigation into the shock loading of air-backed plates and diaphragms of infinite size. His predictions encompassed damage assessment, surface cavitation and the impact of baffles on structural distortion. The studies by Reference TaylorTaylor (1963) and Reference FriedlanderFriedlander (1941) delved into the reflected pressures and structural reactions of plates with varying degrees of rigidity and compliance when subjected to shock waves.

It has been established that vaporous surface cavities can develop on certain structures due to underwater shock, particularly in situations characterized by low hydrostatic pressure and high specimen velocities (Reference Pavlov and GalievPavlov & Galiev 1977). These surface cavities engage with the structure, exacerbating structural deformation. Galiev conducted parallel experimental (Reference GalievGaliev 1979) and theoretical (Reference GalievGaliev 1987) inquiries, positing that the impulse imparted by a shock wave is more critical than its magnitude in determining the resulting deformation. His numerical analyses forecasted the initiation of cavitation stemming from shock wave reflections. In a subsequent study, Galiev examined the phenomena of liquid jets striking submerged targets (Reference GalievGaliev 1995), predicting the genesis of a cavitation zone akin to that identified by Friedlander (Reference FriedlanderFriedlander 1941). He noted this zone's dependency on the wavelength of the incident pulse, which could amplify peak deflection by a factor of two to three. Following this line of inquiry, experimental work was expanded to encompass a range of elastic-plastic plate materials and configurations, such as circular and rectangular plates (Reference GalievGaliev 1997).

2.2.7 Damage from near-field gas bubbles

Ramajeyathilagam et al. undertook a series of UNDEX shock loading experiments and computational damage assessments for various structures, including thin rectangular plates (Reference Ramajeyathilagam, Vendhan and RaoRamajeyathilagam et al. 2000; Reference Ramajeyathilagam and VendhamRamajeyathilagam & Vendham 2004), cylindrical plates (Reference Ramajeyathilagam, Vendham and RaoRamajeyathilagam et al. 2001) and panels with differing curvatures (Reference Ramajeyathilagam and VendhamRamajeyathilagam & Vendham 2003). These tests were likely influenced by bubble effect loading due to the substantial explosive charges and the relatively close detonation distances. The researchers analysed postmortem specimens to gauge plastic deformation and rupture. They compared their findings with a theoretical rupture strain model integrating effective plastic strain (Reference Ramajeyathilagam and VendhamRamajeyathilagam & Vendham 2004) to enhance simulation rupture prediction. The failure modes identified in these panels, depicted in figure 4, predominantly featured edge tearing and central petalling. The research emphasized the criticality of incorporating strain rate-dependent material models in UNDEX experimental simulations. Suresh replicated these experiments using a coupled fluid–structure interaction (FSI) model in LS-DYNA, addressing both non-rupture (Reference Suresh and RamajeyathilagamSuresh & Ramajeyathilagam 2020) and rupture scenarios (Reference Suresh and RamajeyathilagamSuresh & Ramajeyathilagam 2021). Reference Riley, Paulgaard, Lee and SmithRiley et al. (2010) experimented with air-backed circular and square thin plates subjected to near-contact explosions. These scenarios were modelled using LS-DYNA to compare the central ruptures with actual test results, as shown qualitatively in figure 5, without accounting for bubble effects. Reference Hung, Hsu and Hwang-FuuHung et al. (2005) employed experimental pressure data within a finite element analysis framework to predict the response behaviours of specimens.

Figure 4.

Failure modes of mild steel plates: (a) large deformation (mode I), (b) tensile tearing (mode II), (c) central rupture (mode III) and (d) combined shear failure and tensile tearing (modes II–III) (Reference Ramajeyathilagam and VendhamRamajeyathilagam & Vendham 2004).

Figure 5.

Final petal formation in numerical and experimental contact charges for (a) 350WT steel and (b) A653 steel (Reference Riley, Paulgaard, Lee and SmithRiley et al. 2010).

The collapse of UNDEX bubbles generates additional shock waves, a phenomenon studied by Reference Arons, Slifko and CarterArons et al. (1948), Reference AronsArons (1948). These secondary shock emissions, while less intense in peak magnitude than the initial wave, persist longer. The researchers proposed using point impulse integrations of pressure–time plots to evaluate the potential damage from each shock wave. Their findings confirmed that UNDEX gas bubbles oscillate, dissipating their energy gradually. Chapman's study (Reference Chapman1985) provided insights into predicting the peak and duration of these pulses following multiple collapses, also considering the effects of buoyancy, hydrostatic pressure and bubble migration on this dynamics. These studies illustrate that the impulse magnitude of the bubble collapse is significant and should be considered during near-field UNDEX problems.

2.2.8 Challenges and limitations

The insights gained from most studies that analyse small, spark or laser-induced gas bubbles under low hydrostatic pressure conditions do not fully translate to understanding the damage potential from UNDEX events. The shock wave generated by an explosive event poses a risk of structural damage and can inhibit bubble growth. Underwater explosion bubbles are typically larger, facing more intense internal pressures and temperatures and experiencing greater pressure gradients, which lead to stronger buoyancy forces. Whereas small cavitation bubbles may only result in pitting on near-field materials, these effects might not be consistent when scaled up to larger structures or when other forms of damage are considered. Moreover, the potential energy contained in larger UNDEX bubbles can render even thick metal specimens like aluminium or steel relatively compliant in the face of the extreme velocities and pressures generated by bubble jetting.

3.1.1 Polymer-coated plates

Polyurea for coating applications is favoured for its rapid curing, chemical resistance, fire, corrosion resistance and high stiffness under loading rates (Reference Roland, Twigg, Vu and MottRoland et al. 2007). In air blast scenarios, it has successfully reduced deformation in steel and composite plates (Reference Yi, Boyce, Lee and BalizerYi et al. 2006; Reference Tekalur, Shukla and ShivakumarTekalur et al. 2008; Reference Mohotti, Fernando, Weerasinghe and RemennikovMohotti et al. 2021). LeBlanc and others have explored damage mitigation in composite structures, showing that coating thickness and placement significantly influence underwater blast resistance. For instance, curved plates with thicker back-face coatings displayed reduced deformations (Reference Leblanc, Gardner and ShuklaLeblanc et al. 2013; Reference Leblanc and ShuklaLeblanc & Shukla 2015). Li et al. found that front-face coatings on isotropic aluminium plates were more effective than back-face coatings in far-field explosions (Reference Li, Chen, Zhao, Cao, Jiang, Xiao and FangLi et al. 2019), while LeBlanc et al. observed that air-backed flat composite plates coated on the back face with polyurea resisted rupture under similar conditions, as illustrated in figure 6 (Reference LeBlanc, Shillings, Gauch, Livolsi and ShuklaLeBlanc et al. 2016). In contrast, Dai's study on metal plates suggested that front-face polyurea coatings offer better mitigation, highlighting the role of bond strength between the polyurea and the substrate (Reference Dai, Wu, An and LiaoDai et al. 2018). Liu et al. observed that steel plates with front-face coatings had lower final deformation than those with back-face coatings (Reference Liu, An, Wu, Liao, Zhou, Xue and GuoLiu et al. 2022b), illustrated in figure 7, and further research showed that dual-faced polyurea applications could be more effective than single-sided ones, illustrated in figure 8 (Reference Liu, An, Niu, Zhang, Feng, Li and WuLiu et al. 2022a).

Figure 6.

Centre point out-of-plane displacements of composite plates, adjusted by areal weight (AWR), with and without polyurea coating (Reference LeBlanc, Shillings, Gauch, Livolsi and ShuklaLeBlanc et al. 2016).

Figure 7.

Postmortem images for three polyurea coating methods on steel plates (D = 8 mm): (a) the bare steel plate; (b) front-face polyurea coating (4 mm); (c) back-face polyurea coating (4 mm); (d) deformation profiles of the steel plates (Reference Liu, An, Wu, Liao, Zhou, Xue and GuoLiu et al. 2022b).

Figure 8.

Maximum deformation as a function of polyurea coat thickness for plates with various coating methods (Reference Liu, An, Niu, Zhang, Feng, Li and WuLiu et al. 2022a).

3.1.2 Polymer-coated cylinders

Similar investigations have been conducted on hollow cylinders, with studies examining the effects of rubber and polyurea coatings under UNDEX conditions. Kwon's numerical analysis on aluminium and steel cylinders with rubber coatings revealed that, while these coatings can trap shockwave energy, increasing the rubber's shear modulus and thickness could improve energy mitigation (Reference Kwon, Bergersen and ShinKwon et al. 1994). Gauch's experiments on carbon fibre/epoxy hollow cylinders coated with polyurea showed significant damage reduction, particularly with thicker coatings, as illustrated in figure 9 (Reference Gauch, LeBlanc and ShuklaGauch et al. 2018). Nayak developed a machine learning algorithm to analyse the behaviour of polyurea-coated composite cylinders, concluding that such coatings’ protective effects largely depend on their thickness (Reference Nayak, Lyngdoh, Shukla and DasNayak et al. 2022).

Figure 9.

Exterior view of cylinder damage – 2.54 cm charge standoff, (a) uncoated, (b) thin coating (2.34 mm thickness), (c) thick coating (3.04 mm thickness) (Reference Gauch, LeBlanc and ShuklaGauch et al. 2018).

These studies collectively indicate that polymer coatings, especially polyurea, are a promising solution for mitigating damage from UNDEX. The effectiveness of these coatings, however, varies based on factors like coating thickness, placement and the specific material properties of a front-face application. Contrary to some previous findings, the data also indicate that a dual-faced polyurea application has a superior mitigative impact than coating solely on the front or the back face. These results illustrated that polymer coating systems may be an effective solution to mitigate damage from UNDEX. However, the optimization of the coating locations depends on the structural performance and the shock loading condition.

3.1.3 Sandwich composites

As previously mentioned, sandwich composites are increasingly recognized for their effectiveness in mitigating underwater shock, particularly during UNDEX events. These composites, comprising high-strength face sheets and a lightweight core (such as polymeric foams or honeycomb structures), offer a unique combination of strength and energy absorption. The core material is critical in dissipating the energy from shock waves, thereby reducing peak stresses and potential damage (Reference Wanchoo, Matos, Rousseau and ShuklaWanchoo et al. 2021). Much like the differences with coating systems, no universal solution would be ideal for dissipating every type of UNDEX loading. The loading case's specifics must be considered when designing a sandwich structure. There is also much ongoing work on novel core designs, such as auxetic and lattice structures for shock mitigation due to technological advancements in manufacturing (Reference Wanchoo, Matos, Rousseau and ShuklaWanchoo et al. 2021).

3.1.4 Polymer-modified structures

Polymer-modified structures offer enhanced mitigation against blast loads due to their improved toughness and energy absorption properties. The modifications in the matrix resin, including liquid rubbers and pre-formed rubber particles, enhance the overall response of the composite structure to dynamic loadings such as UNDEX loadings. These changes improve the structural damping and energy dissipation under below-failure shock loadings and exhibit limited crack growth under high-intensity shock conditions. Studies have shown that these viscoelastic materials become more dissipative with increasing loading frequency, although they may become hard and non-dissipative at very high frequencies (Reference Wanchoo, Chaudhary, Matos and ShuklaWanchoo et al. 2023). While these modifications significantly increase the impact resistance and shock absorption capabilities, a notable limitation is the base resin's reduction in modulus and strength due to the inclusion of these modifiers (Reference Wanchoo, Chaudhary, Matos and ShuklaWanchoo et al. 2023).

3.2.1 Impedance mismatching

Impedance mismatching is a technique that leverages the difference in acoustic impedance between different materials to enhance stress wave reflection. By strategically choosing materials for the protective layer that significantly differ in impedance from the underlying structure, the stress waves resulting from an explosion may be effectively reflected away. This mismatch minimizes the transmission of shock wave energy into the protected structure, thereby reducing the potential damage. It is particularly effective when the shock wave's energy can be redirected rather than absorbed. Hence, this mitigation approach has the opposite effect of polymer coating.

3.2.2 Protective cladding

Protective or sacrificial cladding involves the use of materials that can absorb and dissipate the energy of the blast wave, thus protecting the underlying structure. This technique reduces the stress transferred to the structure and increases the loading duration, spreading the impact. The cladding materials often undergo plastic deformation or destruction, sacrificing themselves to protect the primary structure (Reference Wanchoo, Pandey, Leger, LeBlanc and ShukaWanchoo et al. 2024). Due to the cladding destruction, this mitigation approach is exceptionally effective but may be limited to a single use, unlike polymer coatings, which may be designed for multi-use. Materials used for cladding are chosen for their energy absorption capabilities and could include advanced composites, foams or specialized alloys.

3.2.3 Geometrical and structural features

Geometrical arrangements in design play a crucial role in deflecting or redirecting blast waves away from critical areas of a structure. By strategically designing the shape and orientation of structures, the impact of a blast wave can be mitigated. This might involve curved surfaces (Reference Leger, Matos, Shukla and JavierLeger et al. 2023), angled panels or specific structural contours that channel the energy of the blast wave in a less harmful direction. This approach is advantageous in large-scale structures where direct absorption of blast energy is impractical. Other new areas that include this type of energy mitigation include metamaterials and metastructures with substructural-scale features incorporated for energy mitigation purposes (Reference He and FanHe & Fan 2021).

3.2.4 Blast-wave disrupters

Blast-wave disrupters are structures or materials placed between the source of the blast and the protected structure. These disrupters act as barriers or baffles, breaking up and diffusing the energy of the blast wave. They can be specially designed structures or materials that fragment or deform under blast pressure, disrupting the blast wave's continuity and intensity. The effectiveness of these disrupters depends on their material properties, design and placement relative to the protected structure and the anticipated blast source (Reference Wanchoo, Matos, Rousseau and ShuklaWanchoo et al. 2021).

3.3.1 Mitigating of bubble jetting

Gibson's study (Reference Gibson1968) delved into the nuances of jetting behaviours associated with different boundary conditions, making significant observations. Particularly noteworthy were the jetting patterns that deviated from free surfaces when subjected to minor buoyant forces – typically observed in smaller bubbles with brief collapse durations. Moreover, in scenarios involving flexible boundaries, small standoff distances led to the absence of jetting, while increased standoff distances resulted in jetting that veered away from the boundary. This behaviour highlights the potential of flexible boundaries to significantly reduce the risk of structural damage caused by the collapse of near-field bubbles.

Gibson suggested the use of flexible coatings on rigid surfaces to create a boundary that would be more compliant. An increase in compliance can neutralize or even reverse the direction of jetting. Reference Gibson and BlakeGibson & Blake (1982) expanded upon this idea by suggesting amendments to the Kelvin impulse formulation, aiming to predict the jetting direction when dealing with elastic boundaries conceptualized as a spring–mass damper system. Their research compared the jetting behaviours of plain flexible rubber sheets against rubber-coated rigid boundaries. The findings revealed that while the flexible and coated surfaces exhibited repulsive jetting, the stiffer rubber options displayed annular jetting. This latter type of jetting led to the bubble splitting and subsequent formation of linear jets shooting in opposite directions post-collapse.

The culmination of their findings and methodologies can be found in their 1987 publication (Reference Blake and GibsonBlake & Gibson 1987). Visual aids in their work offer a clearer understanding of these concepts: figure 10 outlines a parameter space that correlates the impact of the dimensionless standoff, γ, with the buoyancy parameter, δ, on the direction of bubble jetting against rigid boundaries. Furthermore, figure 11 illustrates the varied jetting behaviours when interacting with an array of compliant surfaces, demonstrating both annular jetting, shown in (b,c), and repulsive jetting, depicted in (a–d).

Figure 10.

The $\gamma$–$\delta $ parameter space for buoyant vapour bubbles (Reference Gibson and BlakeGibson & Blake 1982).

Figure 11.

Examples of the interaction of pulsating bubbles with various deformable surfaces (Reference Gibson and BlakeGibson & Blake 1982).

Reference Shima, Gibson and BlakeShima et al. (1988) extended the exploration of the bubble dynamics by examining composite specimens, specifically by adding foam and a viscoelastic layer to brass plates. Their investigations indicated that the phenomena of attractive and repulsive jetting are influenced by the standoff distance and the foam rubber layer's thickness. In parallel, Reference Duncan, Milligan and ZhangDuncan et al. (1996) developed a numerical model that simulates bubble dynamics and the energy transfer to elastic coatings. His model showed good concordance with the experimental results of Reference Gibson and BlakeGibson & Blake (1982) as well as with Reference Shima, Gibson and BlakeShima et al. (1988).

Reference Brujan, Nahen, Schmidt and VogelBrujan et al. (2001) experimentally delved into the effects of elastic moduli on jetting behaviours, uncovering that attractive jetting occurred at very close standoffs. They also identified a substantial range where bubble splitting was prevalent and, beyond that, a range indicative of repulsive jetting. Significantly, they observed that the velocity of jets in bubble-splitting scenarios was notably higher than that in cases of attractive jetting, which implies a potential for greater damage compared with rigid boundaries and highlights possible adverse outcomes from using flexible coatings.

5.1.1 Analytical background

Early explorations by Reference ColeCole (1948), Reference EzraEzra (1973), Reference KeilKeil (1961) and Reference KennardKennard (1944) shed light on the interactions between shock waves and plane plates, complemented by Chertock's equations (Reference ChertockChertock 1953) on how water alters vibration modes and frequencies in submerged solids. Taylor's work (Reference TaylorTaylor 1963), a milestone in FSI modelling, introduced a method to calculate pressure pulse reflection and plate deflection. Building on this, Reference TemperleyTemperley (1950) investigated cavitation during pressure pulse interactions, and Reference Xue and HutchinsonXue & Hutchinson (2004) provided a comparative analysis of FSI in different plate backing scenarios. The work of Reference McMeeking, Spuskanyuk, He, Deshpande, Fleck and EvansMcMeeking et al. (2008) expanded Taylor's model, focusing on varied sandwich core topologies and cavitation FSI, followed by several studies (Reference Ramajeyathilagam and VendhamRamajeyathilagam & Vendham 2004; Reference Deshpande and FleckDeshpande & Fleck 2005; Reference Espinosa, Lee and MoldovanEspinosa et al. 2006; Reference Rajendran, Paik and KimRajendran et al. 2006; Reference Teich and GebbekenTeich & Gebbeken 2013) delving into more complex loadings on air-backed plates. Reference Liu and YoungLiu & Young (2008) refined Taylor's model for water-backed plates within the acoustic range. Despite these advancements, there remains a notable gap in understanding dynamic plate bending response and the effects of near-field fluid-particle velocities post-pressure wave impact, a gap highlighted by Reference Arons and YennieArons & Yennie (1948) and Reference KeilKeil (1961). Recent models, including those by Reference Wang, Liang and LiuWang et al. (2013), Reference Junger and FeitJunger & Feit (1986 and Reference KimKim (2010) have yet to fully address these factors, particularly in the context of impulse loads, underscoring the need for further research in this area.

5.1.2 Shock-to-structure interaction models

The study of FSIs of large plates and nearby dynamic pressure pulses has been widely studied using Taylor's model (Reference TaylorTaylor 1963). However, Taylor's model cannot accommodate various pressure pulse shapes since it was specifically designed for a planar exponential pressure pulse. During near-field UNDEX problems, non-planar wave shapes need to be considered. Taylor's model was also designed for air-backed structures; hence, it did not account for transmitted pressures. It is critical to include transmitted pressures to broaden the applicability of an FSI model in the underwater environment. Recently, Taylor's model was validated in cases without transmitted pressures through air-backed shock loading scenarios (Reference SneddonSneddon 1946). Then, Taylor's model was further modified by including transmitted pressures, representing water-backed structure interactions (Reference SneddonSneddon 1946). A further change to Taylor's model came with integrating the plate-bending dynamics into the solution framework. This was achieved by incorporating Sneddon's (Reference Kishore, Senol, Parrikar and ShuklaKishore et al. 2020) solution for plate response to pressure loads, allowing the model to simulate structural responses beyond a single degree of freedom mass. Kishore et al.'s work (Reference SneddonSneddon 1946) also introduced near-field pressure pulse effects through the after flow fluid velocity formulation originally proposed by Reference Arons and YennieArons & Yennie (1948). This addition has enabled the model to provide a more comprehensive understanding of the interaction between structures and dynamic pressures, leading to improved predictions of pressures and displacements compared with earlier models, although still limited to infinitely large elastic plates (Reference SneddonSneddon 1946).

These progressive enhancements in modelling FSI are outlined in table 1, showcasing the continual evolution in understanding and predicting the complex dynamics of FSI. The key limitation of these models is that they do not account for boundary conditions, cavitation effects or any phase changes in the fluid, which are often present during UNDEX. Hence, these models apply to large structures for modelling the shock–structure interactions for events where the medium does not change throughout the interaction, such as during far-field explosions or similar events.

Table 1.

Summary of analytical shock interaction models.

${p_i}$, incident pressure; ${p_r}$, reflected pressure; ${p_t}$, transmitted pressure; w, plate deflection; m, plate's areal mass density; ${\rho _w}$, water's volumetric density; c, water's sound speed; $\mu $, plate's in-air frequency; b, plate stiffness parameter.

5.1.3 Analytical solutions

The solution step involves first finding the reflected pressure and then, subsequently, all other quantities (Reference Kishore, Senol, Parrikar and ShuklaKishore et al. 2020). Upon applying the equation of continuity across the incident side and the transmitted pressure side of the plate, the governing solution of the model by Reference Kishore, Senol, Parrikar and ShuklaKishore et al. 2020 is obtained as follows:

(5.1)\begin{gather}{p_i}(t) - {p_r}(t) + \frac{c}{R}\int_0^t {({p_i}(\tau ) - {p_r}(\tau ))\,\textrm{d}\tau } = \lambda \int_0^t {{{\dot{p}}_r}(\tau )\varOmega \textrm{(}t - \tau )\,\textrm{d}\tau } , \end{gather}

(5.2)\begin{gather}\varOmega (y) = 1 - \frac{2}{\pi }\,\textrm{ta}{\textrm{n}^{ - 1}}\left( {\frac{k}{y}} \right),\end{gather}

(5.3)\begin{gather}\lambda = \frac{{\pi {s^2}{\rho _w}c}}{{8\rho hb}}.\end{gather}

Where, as described in table 1, ${\rho _w}$ is the density of water, c is the speed of sound in water, ${p_i}(t)$ is the incident pressure, ${p_r}(t)$ is the reflected pressure (at the plate's centre), R is the standoff distance, $\rho $ is the density of the plate, $2h$ is the total plate thickness, s is the radial decay constant of the Gaussian pressure distribution on the plate and k is a combination of parameters that equates to ${s^2}/4b$. Furthermore, the plate stiffness parameter, b, can be calculated form the flexural rigidity of the plate, D, using the following (5.4) and (5.5). Lastly, $\lambda $ represents a FSI parameter that combines the impedance ratio of water to the thin plate and is proportional to the circular area where the bulk of the pressure acts:

(5.4)\begin{gather}{b^2} = \frac{D}{{2\rho h}},\end{gather}

(5.5)\begin{gather}D = \frac{{2E{h^3}}}{{3(1 - {\nu ^2})}}.\end{gather}

5.1.4 Future analytical modelling

Future FSI modelling efforts should focus on a more nuanced understanding of wavefronts and structural interactions. This includes addressing the complexity of the wavefront shapes and their effects on structural response, the impact of structural boundaries in this response and the inclusion of more complex material models. Further work is needed to refine models to accurately capture after flow velocities and their impact on plate deflections and transmitted pressures. This will enhance the fidelity of models in predicting the intricate interactions in underwater environments.

5.2.1 Computational background

With advancing computational resources, including both hardware and sophisticated software codes, simulations of UNDEX and the resulting load on submerged structures have become a potent and cost-effective means for predictive modelling. Various theoretical formulations have been incorporated into both commercial and government-developed codes since the early 1960s (Reference NohNoh 1964). This period marks when computational resources first became capable of handling such complex simulations, although the underlying theories of these codes date back even further.

These codes can be broadly categorized into two hydrocode types: Eulerian and Lagrangian. The fundamental distinction lies in the treatment of materials: in Eulerian hydrocodes, materials flow through a static grid of cells, whereas in Lagrangian hydrocodes, the mesh deforms in response to the load. For both types, energy conservation within the system is a key focus. Additionally, some codes feature a coupling interface that enables fully integrated FSI simulations. This discussion is intended to offer a concise overview of several widely used hydrocode formulations and explore how these methods are applied in simulating UNDEX and their effects on submerged structures.

5.2.2 Eulerian fluid solver overview

Eulerian-based computational hydrocodes are commonly utilized for simulating multi-fluid, inviscid flows, with applications spanning air/underwater shocks, laminar/turbulent flows, weather forecasting and complex fluid mixing. These Eulerian hydrocodes adeptly handle the challenges posed by significant distortions, twisting and volume loss in Lagrangian elements due to extensive deformations and severe loading. A fixed grid with regularly shaped cells is employed in Eulerian formulations, through which the various fluids navigate. This approach involves discretizing the fluid domain into a volume of regularly shaped cells that form the computational grid. Such a set-up enables the fluid to move through the grid while the grid cells remain static, allowing fluid characteristics like velocity, density and pressure to be tracked within the cells.

The fluid's motion through these cells is dictated by fluid dynamics’ conservation laws and equations, including the Euler and Navier–Stokes equations. These laws facilitate the simulation of fluid movement resulting from pressure differences and external forces, such as gravity, across discrete time steps. Mass, momentum and energy transfer calculations are performed at the junctions between cells to ensure conservation. The flow dynamics within the discretized fluid domain is determined by boundary conditions, which can be rigid boundaries (like walls and symmetry conditions) or flow boundaries (such as open faces for inflow or outflow). Eulerian solvers are recognized for their efficient scalability during parallel processing, adaptability for various transient and high-speed flow scenarios and different time scales. However, achieving high resolution in complex models often necessitates very fine meshes, leading to large model sizes and considerable demands on computational resources. Stability issues may also arise if the time-stepping methodology and size are not cautiously selected.

The origins of Eulerian-based fluid solvers date back to the 1960s, with Noh's pioneering Eulerian–Lagrangian formulation, which provided one of the initial time-dependent approaches for numerical simulations of a Lagrangian body moving through an Eulerian fluid domain. Subsequent computational developments have been made at various institutions, including the Department of Energy's creation of the KRAKEN code in the 1970s (Reference DeBarDeBar 1974), refined later for explosives and armour strength simulations (Reference Mandell, Adams, Holian, Addessio, Baumgardner and MossoHolian et al. 1989). A thorough review of hydrocodes up to the 1990s can be found in the work (Reference BensonBenson 1990). More recently, these codes have expanded to include applications in space and astrophysics modelling, as detailed by Trac in (Reference Trac and PenTrac & Pen 2003).

5.2.3 Coupled fluid–structure interaction methods

In many structural applications, understanding how a structure within a fluid domain interacts with the fluid forces it encounters is crucial for predicting the structure's response. Computational methods known as FSI solvers have been developed to meet these simulation needs. Fluid–structure interaction simulations are typically employed to resolve issues where fluid flow impacts a structure, exerting forces upon it, and in turn, the structure applies forces back into the fluid. An example includes the detonation of an explosive charge near a submerged structure, where the blast pressure loads the structure, and conversely, the structure influences the shockwave.

These numerical methods utilize coupled Eulerian–Lagrangian (CEL) and arbitrary Lagrangian–Eulerian (ALE) approaches. Fluid–structure interaction coupling treats the structure and fluid as separate entities linked by a coupling interface. The fluid is modelled using the Eulerian method in the simulations, while the structure is modelled using the Lagrangian approach. They operate independently and exchange crucial variables after each integration step. The structure's surface within the fluid is delineated by a stair-step boundary technique, and interface elements that are either singly or doubly wetted define the contact with the fluid. This method is elaborated with benchmark examples in Reference Wardlaw, Luton, Renzi and KiddyWardlaw et al. (2004). The ALE method models the structure with Lagrangian mesh coordinates that interact with the fluid, while the fluid mesh is Eulerian. Periodic re-meshing reconciles these domains into a single reference frame. Detailed descriptions of ALE formulations can be found in Reference Le Tallec and MouroLe Tallec & Mouro (2000), Reference Souli, Ouahsine and LewinSouli et al. (2000), Reference Legay, Chessa and BelytschkoLegay et al. (2006) and Reference NobleNoble et al. (2017).

A common challenge with CEL and ALE methods is the appearance of non-physical spurious oscillations near the material interfaces within the fluid. These artifacts arise from issues like material density smearing and abrupt changes in the Equation of state at the interface (e.g., water-to-air transition). To tackle these issues, a novel approach known as the ghost fluid method was introduced in the late 1990s by Reference Fedkiw, Aslam, Merriman and OsherFedkiw et al. (1999). This level-set method for multiphase flows achieves sharp delineation at surface discontinuities by using a layer of ghost cells to maintain continuous pressure and velocity profiles at the discontinuity. The ghost method has been refined to enhance its precision and adaptability. Reference Wang, Tang and LiuWang et al. (2008) developed an adaptive ghost method, blending the ghost fluid method with adaptive meshing to harness the strengths of both. Specifically for simulating underwater implosions, Farhat has further refined the ghost method into the ghost fluid method for the poor (Reference Farhat, Rallu and ShankaranFarhat et al. 2008). This method leverages arbitrary equations of state with advanced multi-step time integration techniques, proving effective in handling large density and pressure discontinuities at interfaces.

5.2.4 Underwater explosion simulations

The methods previously discussed model full fluid domains numerically and computationally. Depending on the mesh size and resolution, these methods can lead to sizable models that require substantial computational resources to achieve accurate results. In contrast, analytical modelling methods for predicting structural responses to dynamic shock loading provide more time-efficient solutions for UNDEX scenarios, although they come with their own limitations. Geers introduced the doubly asymptotic approximation (DAA) for transient analysis of structures within a fluid (Reference GeersGeers 1978), using differential equations to approximate the behaviour of wetted surfaces under dynamic pressure conditions. This method obviates the need for discretizing the fluid domain, substantially reducing computational costs. The original implementation was later used to study the steady-state vibration of submerged structures (Reference Geers and FelippaGeers & Felippa 1983). The DAA method underpins the underwater shock analysis code (Reference DeRuntz, Geers and FelippaDeRuntz et al. 1978) for modelling spherical shock impacts on structures. For modelling surface ship responses to UNDEX, Reference ShinShin (2004) used a finite element method for the ship and near-field fluid coupled with the DAA approach to minimize the fluid domain. Lastly, Reference RileyRiley (2010) summarized nine unique analytical techniques for predicting the behaviour of underwater gas bubbles from explosions, correlating models for bubble radius over time.

5.2.5 Explosion and bubble dynamics simulations

5.2.6 Limitations

A key limitation of computational modelling is the computational expanse of high-reliability models. This expanse is extreme for underwater explosion simulations when the underwater bubble effects, such as near-field explosion events, must be considered. This is because of the difference in time scales: the shock interaction event happens on the microsecond scale, the first bubble interaction will happen within tens of milliseconds, and the entire bubble–structure interaction may be in the hundreds of milliseconds scale. The event may be on the second or higher scale for larger explosives. It is widely understood that total simulation time is highly influenced, if not dominated, by the smallest time events. Hence, multiple events in different time scales will not be computationally time-efficient. Recent computation modelling of such events for simple rigid structures and small charge sizes (Reference Javier, Galuska, Papa, LeBlanc, Matos and ShuklaJavier et al. 2020a) took tens of thousands of CPU hours to simulate up to the first bubble interaction. A high-reliability model for an entire near-field explosion event that includes bubble oscillation behaviour would be a massive undertaking, and it is not a feasible modelling solution.

Other limitations of explosion computational modelling are challenges with the fluid and solid solvers. For fluids, most commercial solvers do not handle fluid phase change (such as cavitation effects) well, if at all. This is a huge limitation for near-field explosions, but there is much ongoing work trying to address such limitations with fluid solvers (Reference Tian, Zhang, Liu and TaoTian et al. 2021; Reference Yu, Liu, Wang, Wang, Zhou and MaoYu et al. 2022, Reference Yu, Wang, Sheng and Hao2023a; Reference Yu, Song and ChoiYu, Song & Choi 2023b). Modelling the dynamic and failure behaviour of complex materials such as composites is also challenging for solids. Recent work has addressed some of these challenges by including the interplay behaviour of composites in the modelling (Reference Ulbricht, Boldini, Zhang and PorfiriUlbricht et al. 2024). Material models, such as LS-DYNA's MAT213, have been developed to capture the dynamic viscoelastic-temperature behaviour of anisotropic materials. Nonetheless, to develop these more accurate modelling approaches and material models, a huge undertaking for testing and characterizing the unique behaviour specific to each composite system is often required to validate modelling parameters. This review does not discuss all these intricacies and novelties with computational simulations. These topics warrant additional considerations and their dedicated review when developing computation models.