Numerical Study on the Evolution Characteristics of the Bubble Dynamics and Free Surface Structures in Extremely Shallow Water Explosion

1. Introduction

Underwater explosions in shallow water hold significant strategic and engineering importance in nearshore environments, and their study provides practical value for the assessment of underwater blasting operations and the design of port protection. In shallow water conditions, explosion bubbles often exhibit nonlinear behaviors such as breaking and collapse, with the resulting jets potentially impacting underwater pipelines and seabed structures, while the water spike may threaten vessel navigation safety or coastal engineering projects [1,2,3]. Under extremely shallow water conditions, where the water depth is comparable to the maximum bubble radius scale, the bubble dynamics differ markedly from those underwater explosions in conventional shallow water. In such scenarios, both the water surface and the seabed boundary simultaneously influence the bubble pulsation process, leading to a strongly nonlinear coupling effect among bubble morphology and water spike evolution characteristics. This makes conclusions from traditional shallow water explosion studies not directly applicable, and numerical research under these specific conditions remains insufficient. Therefore, conducting numerical simulations of underwater explosion in extremely shallow water is of great significance for gaining deeper insights into the bubble pulsation mechanisms and water spike generation characteristics.

In terms of experimental research, scholars both domestically and internationally have conducted extensive work on the dynamic behavior of bubbles near free surfaces or solid walls, revealing characteristics such as bubble pulsation, jet formation, and water spike evolution. Early on, Chahine et al. [4,5] experimentally investigated the influence of a free surface on bubble pulsation, finding that the free surface exerts a repulsive effect on the bubble and induces a downward-directed jet of bubbles. Blake et al. [6,7] further explored the role of the free surface on bubble contraction and jet direction, noting that during the collapse phase, the bubble can generate a water spike opposite to the jet direction. Robinson et al. [8] used laser-induced cavitation bubble experiments to study the dynamics of single and multiple bubbles near a free surface, observing that bubbles consistently produce high-speed jets away from the free surface under its influence, rapidly forming spike-shaped water spikes. Zhang et al. [9] employed bubble experiments using a high-voltage spark generator to record interfacial instabilities such as "breaking wrinkles" and "spraying water film" during bubble pulsation near a water surface. Yan et al. [10] conducted bottom-charge underwater explosion experiments, discovering that the maximum bubble radius under such conditions could reach 1.05 times that of a free-field underwater explosion. Xu, Wang et al. [11,12] and Zhang et al. [13] obtained key parameters, such as the bubble radius multiplication factor and equivalent charge coefficient, through experiments involving bottom-sitting and water-soil interface explosions. While the aforementioned experimental studies have largely focused on either a single rigid boundary or a free surface boundary. There is still a lack of systematic experiments on the coupling of bubbles and multiple boundaries under extreme shallow water conditions. Moreover, the experimental methods have limitations in capturing the fine dynamic information within the flow field, making it necessary to employ numerical approaches for more in-depth mechanistic investigation.

Building upon experimental research, various numerical methods have provided effective approaches for revealing the coupling mechanisms between bubbles and the free surface as well as the bottom boundary in extremely shallow water explosions. Dong et al. [2] employed LS-DYNA to analyze the influence of explosion depth on bubble pulsation and load characteristics in shallow water explosions. To address the complex evolution of bubbles near a free surface, Tian et al. [14] revealed the re-closure phenomenon after the breaking of a bubble based on the Eulerian finite element method (EFEM), pointing out that the water film can re-merge due to pressure differences on its two sides. Koukouvinis and Li et al. [15,16] summarized three breaking modes of bubbles under the influence of the free surface by using a compressible two-phase flow solver combined with the volume of fluid (VOF) method. Singh et al. [17] studied the effects of symmetric and asymmetric distributions of multiple bubbles near a free surface on their breaking behavior. Based on an inviscid and incompressible assumption, Liu et al. [18] investigated the influence of distance parameters on the coupling characteristics between bubbles and the free surface in shallow water environments using the finite volume method (FVM). Furthermore, Ming and Tian et al. [19,20] used an improved smoothed particle hydrodynamics (SPH) method to simulate the shock wave, surface disturbance, and bubble expansion process in near-surface underwater explosions, validating the effectiveness of the meshless methods in handling free surfaces with large deformations. Zhang and Li et al. [21,22,23] developed a novel theoretical model for bubble dynamics that integrates multiple influencing factors, including fluid compressibility, viscosity, phase change, bubble interaction, bubble migration, boundary effects, and ambient flow fields. Regarding the seabed boundary that can be treated as a rigid wall, Tian et al. [24] studied the shock wave loads and bubble dynamics of near-bottom underwater explosions, while Li and Zhang et al. [25,26,27,28] systematically analyzed the growth, collapse, and jet behavior of bubbles near solid walls using the boundary element method (BEM).

By comprehensively comparing the aforementioned numerical methods, the EFEM combines the advantage of the Eulerian framework in handling large-deformation interfaces with the numerical accuracy of the finite element method. It can naturally capture strongly nonlinear phenomena in extremely shallow water explosions, such as bubble breaking, water film re-closure, jet formation, and water spike evolution. Therefore, this paper will establish an axisymmetric two-dimensional numerical model based on the EFEM to systematically investigate the influence of key parameters, including water depth, charge weight, and detonation distance, on the bubble pulsation characteristics in extremely shallow water explosions.

2. Numerical Model and Method

2.1. Numerical Model

Bubble pulsation is influenced by the attraction of rigid walls, leading to the generation of high-speed jets directed toward the wall. Simultaneously, the presence of a free surface induces the formation of jets away from the free surface and water spikes moving in the opposite direction. These phenomena interact and couple with each other, collectively determining the characteristics of bubble pulsation. This study focuses on extremely shallow water environments, paying particular attention to the complex dynamic characteristics exhibited by bubbles in such explosions under the influence of different boundaries (the free surface and the rigid seabed). Therefore, it is of vital importance to reveal the influence patterns of key parameters such as water depth, charge weight, and the relative position of the explosive on the bubble pulsation and free surface evolution mechanisms in extremely shallow water explosions.

To systematically study the above issues, this paper defines the dimensionless parameters related to the seabed water depth H, the explosives depth d and the charge weight W as follows:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}\lambda =H/R_{\mathrm{m}}\\ \gamma =d/H\\ \delta =\sqrt{{\displaystyle \frac{\rho g{R}_{\mathrm{m}}}{P_0}}}\end{array},\right.\end{array}$$

(1)

in which $\lambda$ is the water depth parameter, defined as the ratio of the actual water depth H to the maximum bubble radius $R_{\mathrm{m}}$. The condition $\lambda \to \infty$ corresponds to the deep-water scenario, whereas $\lambda \le 1.5$ represents the extremely shallow water scenario of primary interest in this study. $\gamma$ is the position parameter, defined as the ratio of the distance from the charge to the free surface d to the water depth H, a lager d indicates a weaker influence of the free surface on bubble dynamics. $\delta$ is the buoyancy parameter, characterizing the relative importance of gravity during the bubble pulsation process and serving as a key factor governing the uplift height of the water spike at the free surface. Herein, the ambient pressure is expressed as $P_{\infty }=P_0+\rho gd$, where $P_0$ is taken as 101325 Pa. The maximum bubble radius $R_{\mathrm{m}}$ can be calculated from the charge weight W using the Cole’s empirical formula [29]:

$$R_{\mathrm{m}}=3.38{\left({\displaystyle \frac{W}{d+10.3}}\right)}^{1/3}.$$

(2)

Considering the axisymmetric nature of the extremely shallow water explosion problem, a two-dimensional axisymmetric computational model is established as shown in Figure 1. The model employs a cylindrical coordinate system $(r,z)$, with the the origin O set at the center of the computational domain’s bottom boundary. Here, r and z denote the radial and vertical direction, respectively. The computational domain has a total radial extent of $20R_{\mathrm{m}}$, exploiting symmetry, the model’s radial dimension is taken as $10R_{\mathrm{m}}$, with a vertical height of $10R_{\mathrm{m}}$. The bottom and axisymmetric boundaries of the computational domain are configured as a rigid wall, while the top and side boundaries are all set as non-reflecting boundary conditions to simulate infinite far-field effects. Meanwhile, a hybrid mesh strategy combining uniform and stretched meshes is adopted in the computational domain to further reduce the computational scale. Specifically, more finely uniform meshes are applied only within the region of [$0,3R_{\mathrm{m}}$] × [$0,3R_{\mathrm{m}}$], while the mesh size is gradually increased (stretched) in the remaining areas.

2.2. Eulerian Finite Element Method

The EFEM based on operator splitting technique provides an effective numerical framework for solving the bubble dynamics in underwater explosion problems. To accurately capture the highly nonlinear evolution of the gas–liquid interface during the explosion process(such as bubble breaking, fusion and water film re-closure behaviors), the VOF method is introduced in this work. This approach reconstructs and tracks the moving interface by defining and updating the volume fraction of the target fluid in each computational cell, and has significant advantages in dealing with such large interface deformation problems.

During underwater explosion, both the Reynolds number $Re$ and the Weber number $We$ are sufficiently large, indicating that inertial forces dominate, and the hydrodynamic pressure difference is much greater than the additional pressure induced by surface tension. Moreover, the bubble pulsation period is relatively short, allowing insufficient time for full turbulence development. For the reasons mentioned above, the fluid medium can reasonably be simplified as an inviscid, compressible ideal medium with negligible surface tension effects [14,30,31,32]. Within this framework, the governing equation system consists of the volume fraction transport equation together with the conservation equations of mass, momentum, and energy, as follows:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{\displaystyle \frac{\partial {\alpha }_i}{\partial t}}+\nabla ·\left({\alpha }_i\mathbf{u}\right)={\alpha }_i\nabla ·\mathbf{u}\hfill \\ {\displaystyle \frac{\partial \left({\alpha }_i{\rho }_i\right)}{\partial t}}+\nabla ·\left({\alpha }_i{\rho }_i\mathbf{u}\right)=0\hfill \\ {\displaystyle \frac{\partial \left(\overline{\rho }\mathbf{u}\right)}{\partial t}}+\nabla ·(\overline{\rho }\mathbf{u}\otimes \mathbf{u})+\nabla p=\overline{\rho }\mathbf{g}\hfill \\ {\displaystyle \frac{\partial \left({\alpha }_i{\rho }_i{e}_i\right)}{\partial t}}+\nabla ·\left({\alpha }_i{\rho }_i{e}_i\mathbf{u}\right)+p{\alpha }_i\nabla ·\mathbf{u}=0\hfill \end{array},\right.\end{array}$$

(3)

where ${\alpha }_i$ is the volume fraction of the i-th fluid phase, ${\rho }_i$ and $e_i$ represent the density and specific internal energy of this phase, respectively ($i=0,1,2$ correspond to the water, air, and TNT explosion product phase), $\mathbf{u}=(u_r,u_z)$ is the velocity vector, p is the fluid pressure, $\mathbf{g}$ is the gravitational acceleration vector, and $\overline{\rho }$ represents the average density of the mixed cell, defined as $\overline{\rho }=\sum {\alpha }_i{\rho }_i$.

For the axisymmetric model established in this paper, the three-dimensional governing equations are simplified to a two-dimensional axisymmetric form, which can significantly improve the solution efficiency while maintaining the calculation accuracy. In the cylindrical coordinate system $(r,z)$, the expression for the velocity divergence is:

$$\nabla ·\mathbf{u}={\displaystyle \frac{\partial u_r}{\partial r}}+{\displaystyle \frac{\partial u_z}{\partial z}}+{\displaystyle \frac{u_r}{r}}.$$

(4)

To close the system of governing equations, corresponding equations of state (EOS) need to be introduced to describe each fluid phase. In this work, the Tammann EOS [33,34] is employed for the water and air media, and its expression is:

$$p=\rho e(\gamma -1)-\gamma P_w,$$

(5)

where $\gamma$ and $P_w$ represent the specific heat ratio and reference pressure of the fluid, respectively, and their relevant values can be found in Table 1.

The explosion product of TNT is described by the Jones-Wilkens-Lee (JWL) EOS [35,36], which can be expressed as:

$$p=A\left(1-{\displaystyle \frac{\omega \beta }{R_1}}\right){\mathrm{e}}^{-R_1/\beta }+B\left(1-{\displaystyle \frac{\omega \beta }{R_2}}\right){\mathrm{e}}^{-R_2/\beta }+\omega \rho e_{in},$$

(6)

in which $\beta =\rho /{\rho }_0$ is the ratio of the explosion product density to the explosive’s initial density, and A, B, $R_1$, $R_2$, $\omega$, $e_{in}$ are constants related to the type of explosive. The specific constant values of TNT are shown in Table 2.

The core idea of the operator splitting technique in the EFEM is to separate the transport terms in the numerical calculation, and then split the solution process within each time step into Lagrangian phase and Eulerian phase to ensure numerical stability and computational accuracy under conditions involving large density ratios and strong physical discontinuities. For the sake of convenience in description, Eq. 3 can be written uniformly as [37]:

$${\displaystyle \frac{\partial \mathbf{E}}{\partial t}}+\nabla ·\left(\mathbf{Eu}\right)=\mathbf{S},$$

(7)

with $\mathbf{E}$ representing the conserved variables such as mass, momentum, energy, and volume fraction, and $\mathbf{S}$ denoting the source term.

Based on the operator splitting technique, Eq. 7 can be decoupled into

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}{\displaystyle \frac{\partial \mathbf{E}}{\partial t}}+\mathbf{E}\nabla ·\mathbf{u}=\mathbf{S}\hfill & \mathrm{for}\mathrm{Lagrangian}\mathrm{phase}\hfill \\ {\displaystyle \frac{\partial \mathbf{E}}{\partial t}}+\mathbf{u}·\nabla \mathbf{E}=\mathbf{0}\hfill & \mathrm{for}\mathrm{Eulerian}\mathrm{phase}\hfill \end{array}.\right.\end{array}$$

(8)

The Lagrangian phase is employed to obtain the changes in physical quantities induced by pressure work, external forces, etc. After remapping the results onto the mesh, they are substituted into the Eulerian phase for solution, thereby completing the calculation of convective effects. This splitting scheme effectively maintains numerical stability and conservation throughout the computation.

2.3. Nondimensionalization

To facilitate subsequent calculations and the discussion and analysis of results, the maximum bubble radius $R_{\mathrm{m}}$, water density ${\rho }_w$, and ambient pressure $P_{\infty }$ are selected as the basic reference scales for length, density, and pressure respectively, that is, $L_{ref}=R_{\mathrm{m}}$, ${\rho }_{ref}={\rho }_w$ and $P_{ref}=P_{\infty }$. Then the scales of other relevant quantities can be naturally derived, such as $U_{ref}=\sqrt{P_{\infty }/{\rho }_w}$ for velocity, $T_{ref}=R_{\mathrm{m}}\sqrt{{\rho }_w/P_{\infty }}$ for time, $A_{ref}=P_{\infty }/\left({\rho }_w{R}_{\mathrm{m}}\right)$ for accelerated velocity, and so on. Subsequently, the subscript * will be uniformly used to represent the dimensionless variable, that is, ${(•)}^{*}=(•)/{(•)}_{ref}$.

2.4. Convergence Test and Model Validation

The size of mesh within the computational domain directly affects the accuracy of multiphase interface capturing. To evaluate the influence of it on the computational results, a mesh convergence analysis was performed. The minimum mesh size $L_{\mathrm{e}}$ of the computational domain was set to 0.03$R_{\mathrm{m}}$, 0.02$R_{\mathrm{m}}$ and 0.01$R_{\mathrm{m}}$, respectively, with the water depth parameter $\lambda =0.224$, position parameter $\gamma =0.191$, and buoyancy parameter $\delta =0.0164$. Figure 2 illustrates the time history of the water spike height under three types of $L_{\mathrm{e}}$. As shown in the figure, when the mesh size decreases from 0.03$R_{\mathrm{m}}$ to 0.02$R_{\mathrm{m}}$, the computational results change significantly. However, further decreases to 0.01$R_{\mathrm{m}}$ shows that the water spike height curve closely matches the results of $L_{\mathrm{e}}=0.02R_{\mathrm{m}}$, indicating that the computational results have essentially converged. Considering both computational accuracy and efficiency, the mesh size of $L_{\mathrm{e}}=0.02R_{\mathrm{m}}$ is adopted for subsequent studies.

After determining the grid resolution, it is necessary to further verify the accuracy of the numerical method in simulating bubble dynamics and free surface evolution. The spark-generated bubble experiment near the free surface conducted by Zhang et al. [9] was selected for comparative validation. In this experiment, the maximum bubble radius is $R_{\mathrm{m}}=27.75$mm, the bubble period is $T=4.95$ms, and the distance parameter from the initial bubble center to the free surface is ${\gamma }_f=d/R_{\mathrm{m}}=0.85$. It is determined that the numerical simulation parameters corresponding to this paper remain consistent with the settings used in the grid convergence analysis.

Figure 3 compares the experimental snapshots and numerical simulation results of bubble morphology and free surface evolution at several typical time instants. During the initial stage (2.18ms and 3.18ms), the bubble generates a downward jet under the repulsive effect of the free surface. Subsequently, the jet penetrates the lower surface of the bubble at 4.64ms, and a obvious water spike has already formed on the free surface. As the bubble continues to pulsate downward at 8.68ms, some distinct "wrinkles" appears around the water spike, and the early wrinkles will splash to form the water skirt. At 16.36ms, the water spike rises continuously at high speed, causing the water spike to break, and the surrounding "wrinkles" also begin to splash. The above results are in good agreement with the experimental observations, which indicates that the model established in this work successfully reproduces the entire bubble dynamics process observed in the experiment and can effectively simulate the complex nonlinear coupling effect between the bubble and the free surface. Figure 4 presents a comparison of the water spike height variation over time. Throughout the pulsation period, the numerical simulation curve shows good agreement with the experimental data, especially during the main stages of jet development, when the two almost completely coincided. This quantitative comparison further confirms that the numerical method employed in this study provides reliable accuracy in predicting the water spike dynamics induced by bubbles near the free surface.

Based on the comparative analysis and validation presented above, the EFEM with VOF interface capturing adopted in this work can accurately and stably simulate the strongly coupled dynamic processes between an underwater explosion bubble and the free surface in an extremely shallow water environment, laying a reliable methodological foundation for subsequent parametric research.

3. Results and Discussion

3.1. Underwater Explosion Bubble Dynamics in Shallow Water

To deeply investigate the dynamic evolution of the underwater explosion bubbles and their coupling characteristics with the free surface in extremely shallow water environments $\lambda \le 1.5$, we selected the water depth parameter $\lambda$, the position parameter $\gamma$, and the buoyancy parameter $\delta$ as control variables. A total of 13 numerical cases were designed, with details provided in Table 3. Among them, Case 0 is set as the baseline condition to obtain the basic evolution process of bubble dynamics in shallow water explosions. Subsequently, three series of cases (Case 1–4, Case 5–8, and Case 9–12) were constructed by respectively fixing two of the parameters while varying the third based on Case 0. This approach aims to quantitatively reveal the influence of each individual parameter on the bubble morphology and free surface evolution in shallow water explosions, thereby elucidating the nonlinear coupling mechanism between the bubble and the free surface.

Figure 5 shows the contours of the dimensionless vertical velocity $v^{*}$ and pressure $p^{*}$ during the bubble evolution under the standard case (Case 0), fully revealing the coupled dynamic process between the bubble and the free surface in an extremely shallow water explosion.

In the initial stage ($t^{*}=0.05-0.18$), the explosion bubble expands rapidly due to its internal pressure being significantly higher than the ambient pressure, forming a local high-pressure zone in the surrounding water. The shock wave undergoes multiple reflections between the free surface and the bottom boundary. The velocity field is dominated by inertial flow, and the free surface above the bubble has already begun to lift and deform. As the pressure difference inside and outside the bubble gradually decreases, its expansion rate tens to slow down. By $t^{*}=0.71$, the free surface is significantly raised, and the lower part of the bubble tends to flatten near the bottom boundary. Simultaneously, under the combined effect of the repulsion of the free surface and the attraction of the bottom boundary, the fluid above the bubble is violently accelerated to form the prototype of a high-speed jet, and a pressure concentration appears at the head of the jet, indicating a strong load impact. During the stage of $t^{*}=1.42-1.87$, water accumulates at the free surface to form a water spike. Meanwhile, the high-speed jet penetrates the lower surface of the bubble, imposing a significant jet impact load on the bottom boundary and generating radial jets along the bottom boundary, causing the bubble to develop into a circular pulsation in all directions. As the bubble further contracts and collapses to a smaller volume, a shear layer develops between the velocity-disturbed water around the bubble and the fluid in the root region of the water spike due to their velocity difference. This triggers interfacial instability, causing initial wrinkled structures to begin forming in the root region of the water spike. In the later evolution stage ($t^{*}=3.17-3.73$), the water spike continues to rise and its shape gradually becomes thinner. Under the action of multiple bubble pulsations, the wrinkles around water spike increase. The earlier-formed wrinkles splash, creating the water skirt that gives the water spike a crown-shaped morphology. Subsequently, the splashing phenomenon of the free surface wrinkles becomes increasingly intense, and new wrinkles also formed between the water spike and the water skirt. This case clearly illustrates the evolution of a series of nonlinear phenomena induced by the explosion bubble during its multiple pulsations in shallow water, including the water jet, water spike, free surface wrinkles, and their splashing.

Figure 6 shows the multi-period variation of bubble equivalent radius $R^{*}$ with time for this case. Combined with the flow field analysis results in Figure 5 , four characteristic time instants($t_1$, $t_2$, $t_3$, and $t_4$)can be further extracted to respectively characterize the key dynamic phenomena during the bubble pulsation process. Here, $t_1$ and $t_4$ correspond to the moments after $R^{*}$ reaches local maximum in the first and third bubble pulsation periods, respectively, while $t_2$ and $t_3$ represent the transitional stages between the three periods. It can be seen that at $t_1$, a high-speed jet begins to form, and the free surface simultaneously rises significantly. . By $t^{*}=t_2$, a significant water spike has developed on the free surface, the bubble undergoes contraction and collapse, the jet penetrates the bubble and impacts the bottom boundary, and the initial wrinkles begin to appear at the root of the water spike. At $t_3$, the bubble contracts and collapses again, and the number of wrinkles around the water spike increases. The earlier wrinkles splash to form the water skirt, resulting in crown-shaped water spike. When $t^{*}=t_4$, the wrinkling and splashing phenomena on the free surface become even more pronounced. These four characteristic time instants capture the main flow field structures during the evolution of shallow water explosion bubbles, providing a clear reference for the stages of dynamics in subsequent parametric studies.

3.2. Bubble Dynamics with Different Water Depth Parameters

In extremely shallow water explosions, the water depth parameter $\lambda$ is the core dimensionless parameter characterizing the intensity of environmental constraints in a finite water depth. When $\lambda$ is small, the bubble contacts both the free surface and the bottom boundary during the initial expansion phase, resulting in strongly nonlinear bubble pulsation. In this subsection, the numerical cases are set up by varying the magnitude of the water depth parameter while keeping the position and buoyancy parameters unchanged (Cases 1–4), with specific parameters listed in Table 3.

Figure 7 shows a comparison of the time histories of the bubble equivalent radius under different water depth parameters. When $\lambda =0.5$ (Case 1), the maximum bubble radius is twice the water depth. The bubble contacts the free surface during its initial expansion phase, causing it to rupture. The internal gas and energy of the bubble are thus released.

Consequently, the first bubble pulsation cycle in Case 1 is incomplete and is considered an extreme case of interface breaking, which cannot be directly compared and discussed with the other cases. Furthermore, due to the energy dissipation caused by the bubble rupture, its subsequent pulsation periods are significantly shortened. For the remaining cases, in the initial stage of the first pulsation period, the influence of boundary conditions is minimal, and the bubble morphology is essentially the same across different cases. When the bubble volume reaches its peak in different cases, a smaller water depth parameter results in a stronger constraining effec from the free surface and bottom boundary on the bubble expansion, which leads to a reduction in the peak value of the bubble equivalent radius. Additionally, the bubble collapse velocity gradually increases, thereby shortening the bubble pulsation period. Meanwhile, the minimum equivalent bubble radius in the first period decreases as the water depth parameter increases. It also can be observed that the subsequent bubble pulsation period still exhibits a positive correlation with the water depth parameter. As energy dissipates, the maximum bubble equivalent radius in each period gradually decreases. Due to stronger nonlinear coupling between the bubble, the free surface, and the bottom boundary in the later stages, the correlation between the minimum bubble equivalent radius and the water depth parameter becomes insignificant.

Figure 8 presents contours of the dimensionless vertical velocity $v^{*}$ and pressure $p^{*}$ at four characteristic time instants for different water depth parameters (Cases 1-4). Comparing the contours at $t_1$, for $\lambda =0.5$, the bubble contacts the bottom boundary during its expansion phase, resulting in a nearly flattened bottom.

The top of the bubble pushes and squeezes the water layer between it and the free surface. After the bubble breaking, this part of the liquid is accelerated upward by the explosion products, while the remaining gas forms a cavity. Over time, the water film surrounding the cavity closes, forming an upward-moving water spike, while the action of gravity contributes to the formation of a downward jet. When $\lambda =1.0$, the maximum bubble radius is approximately equal to the water depth. Influenced by the free surface, the top of the bubble is vertically stretched, causing the initial water spike with sharp characteristics. As the water depth parameter increases to 1.25 and 2.0, the influence of the free surface and bottom boundary on the bubble diminishes, causing it to assume a more spherical shape. Correspondingly, the tip of the initial water spike becomes increasingly rounded. By time $t_2$, after the jet penetrates the bubble, the water spike for $\lambda =0.5$ is primarily formed by the closure of the cavity’s water film. Overall, it appears as a narrow, tall spike shape. Furthermore, due to the initial bubble breaking, its height is lower than that of $\lambda =1.0$. From the subsequent cases, it can be observed that a larger water depth parameter leads to a lower water spike height, a broader overall morphology, and a smaller vertical velocity at its tip. From times $t_3$ and $t_4$, it is visually evident that smaller water depth parameters (0.50 and 1.00) cause the kinetic energy generated by bubble pulsation to be highly concentrated in the vertical direction. After the water spike forms, it continues to push the remaining water body upward. The velocity difference between this flow and the initially formed water spike generates early-stage wrinkles. Later, after these wrinkles splash to form the water skirt, other subsequently generated wrinkles in the surrounding area are not prominent. For $\lambda =$ 1.25 and 2.00, the water layer between the bubble and the boundaries is thicker, resulting in a flatter water spike morphology. The influence range of the bubble pulsation is larger, which perturbs the free surface more extensively. Consequently, a significant number of distinct wrinkles appear after the water skirt has formed. When $\lambda =2.0$, new wrinkles appear between the water spike and the water skirt at $t_4$.

Figure 9 shows the relationship curves between the characteristics (vertical velocity $v_{jet}^{*}$ and radial width $w_{jet}^{*}$) of the initial jet and the water depth parameter before the jet penetrates the bubble.

It can be observed that the jet width gradually increases as $\lambda$ increases. When the water depth parameter is small, the water layer between the bubble and the free surface is thinner. The energy generated by the underwater explosion readily leads to the rapid formation of a water spike at the free surface, while the attractive effect of the bottom boundary on the bubble is stronger, resulting in a higher jet velocity and a narrower morphology. The jet velocity exhibits a counterintuitive increase when $\lambda =2.00$. This is because the repulsive effect of the free surface on the bubble plays a significant role under this condition, and it combines with the attractive effect of the bottom boundary, thereby increasing the jet velocity pointing downward.

Figure 10 illustrates the variation of water spike height over time for different water depth parameters.

It can be visually observed that a larger water depth parameter corresponds to a lower upward velocity and a lower final height of the water spike, while the onset of its fallback also occurs earlier. When $\lambda =0.50$, the free surface undergoes breakup at the initial moment. The water spike only begins to form after the cavity’s water film recloses, which causes a delayed onset in the variation of its height.

3.3. Bubble Dynamics with Different Position Parameters

In extremely shallow water explosions, the position parameter $\gamma$ is the key parameter governing the coupling between the explosion bubble, the free surface, and the bottom boundary. This subsection further investigates the influence patterns of different position parameters on bubble dynamics by employing the control variable method, with the case setups detailed in Table 3.

From the evolution of the bubble equivalent radius shown in Figure 11,when $\gamma =0.1$, the explosive charge is located relatively close to the free surface, and the top of the explosion bubble attains a higher velocity, rapidly pushing and squeezing the thin water film between the bubble and the free surface, ultimately leading to bubble breaking.

Consequently, subsequent bubble pulsation periods in this case are not pronounced. Comparing the remaining four cases, as the position parameter increases, the maximum bubble equivalent radius within each period gradually increases, and the bubble pulsation period also progressively lengthens. According to potential flow theory, a larger position parameter reduces the constraining effect of the free surface on the bubble. This favors bubble expansion, allowing it to achieve a larger equivalent radius, while simultaneously slowing down the bubble contraction process, thereby extending the pulsation period.

As shown in Figure 12, at time $t_1$, when $\gamma =0.1$, the bubble breaking causes part of the liquid to be ejected upward with the explosion products. The closure of the water film leads to the accumulation of more liquid along the axis, forming a water spike that moves upward at high speed. When $\gamma =0.3$, the free surface does not rupture due to the bubble’s expansion. Instead, it moves upward at high speed following the bubble’s upper surface, forming a sharp water spike. For water depth parameters of 0.7 and 0.9, the bubble pulsation is more strongly influenced by the bottom boundary, and the bottom of the bubble becomes flatter. A larger position parameter results in a greater relative contact area between the bubble and the bottom boundary, causing the bubble to tend toward a hemispherical shape, while the water spike morphology gradually transitions from a sharp spike to a flatter form. At time $t_2$, it can be observed that, except for $\gamma =0.1$, a larger position parameter corresponds to a lower water spike height, a broader overall morphology, and a lower velocity at the tip of the water spike. The water spike for $\gamma =0.1$ is primarily formed after the water film closure. The bubble’s energy is partially dissipated, and its interaction with the spike is brief, resulting in a lower water spike height at this moment compared to $\gamma =0.3$.. The water spike for $\gamma =0.3$ is generated by the local high pressure between the bubble and the free surface, acquiring a high initial velocity and forming a narrow, tall sharp water spike. Furthermore, a high-pressure region can still be observed at the root of the water spike, promoting its continued upward growth, expansion, and the generation of wrinkles at the root of the water spike. When the position parameters are 0.70 and 0.90, the bubble primarily couples with the bottom boundary. The water spike morphology becomes flatter, and its height gradually decreases as the position parameter increases. At times $t_3$ and $t_4$, no wrinkles are observed around the water spike for $\gamma =0.1$, and part of the liquid in the upward water spike begins to fall back, causing the overall shape to resemble a fountain-shaped water spike. When $\gamma =0.3$, the wrinkles at the root of the water spike undergo splashing and eventually develop into a water skirt. As the position parameter increases to 0.70 and 0.90, the water layer between the bubble and the free surface becomes thicker, resulting in a gentler water spike morphology. The pulsation of bubbles, when reflected by the bottom of the water, disturbs the free surface and generates numerous wrinkles. For $\gamma =0.9$, the wrinkles between the water spike and the water skirt splash to form a new water skirt, the height of which already exceeds that of both at $t_4$.

Figure 13 presents the results of characteristic variables related to the jet under different position parameters. It can be observed that, except for the exceptional case of bubble breaking of $\gamma =0.1$, the jet velocity gradually decreases and the jet width progressively increases as the position parameter rises. This indicates that the presence of the free surface is a key condition for generating a high-speed, narrow jet, while avoiding the exceptional case of bubble breaking.

Figure 14 illustrates the variation of water spike height under different position parameters. For the case with $\gamma =0.1$, bubble breaking leads to significant energy dissipation, consequently, the water spike formed after water film closure cannot attain a great height. For the remaining cases, it is visually evident that a larger position parameter corresponds to lower water spike velocity and height, as well as a faster rate of fallback. When the position parameters are 0.70 and 0.90, the bubble primarily interacts with the bottom boundary, resulting in a lower initial velocity and a limited final height of the water spike. The water spike height exhibits some fluctuation with the bubble pulsation.

3.4. Bubble Dynamics with Different Buoyancy Parameters

During the coupling process between the explosion bubble and the free surface as well as the bottom boundary, in addition to the water depth and position parameters affecting bubble pulsation, the buoyancy parameter $\delta$ is also a factor that cannot be ignored. This subsection continues to select cases with different buoyancy parameters for numerical investigation. The specific parameters are listed in Table 3. Among them, Case 9 with $\delta =0$ considers the effect of a zero-gravity environment on bubble pulsation, while the remaining computational conditions are the same as those in Case 10 with $\delta =0.2$.

Figure 15 shows a comparison of the time histories of the bubble equivalent radius under different buoyancy parameters.

During the initial expansion phase of the first pulsation period, the bubble radius varies little with changes in the buoyancy parameter. As the bubble expands to its maximum volume, it can be seen that, except for Case 9, the maximum equivalent radius of the bubble increases as the buoyancy parameter increases. This is because the zero-gravity condition in Case 9 imposes a weaker constraining effect on bubble expansion during this period. In the bubble collapse stage, the collapse velocity gradually accelerates with an increase in the buoyancy parameter, the pulsation period shortens accordingly, and the minimum equivalent radius of the bubble decreases as the buoyancy parameter increases. Entering the second pulsation period, the maximum bubble equivalent radius continues to show a similarly strong correlation with the buoyancy parameter. However, as the nonlinear coupling between the bubble and the free surface as well as the bottom boundary intensifies in the later stages, causing the bubble to transition into toroidal pulsation, the correlation between the pulse width and the minimum equivalent radius in this period with the buoyancy parameter becomes insignificant. The same holds true for the subsequent periods of evolution.

Figure 16 presents contours of the dimensionless vertical velocity $v^{*}$ and pressure $p^{*}$ at four characteristic time instants of $\delta =$ 0.0, 0.2, 0.6, and 0.8. Comparing the results at time $t_1$, it can be observed that under the condition of $\delta =0$, the bubble expansion is influenced by the repulsive effect of the free surface and the attractive effect of the bottom boundary. Shortly after the bubble volume reaches its maximum, its upper surface becomes depressed, forming an initial downward jet. Concurrently, the surrounding water flow begins to converge towards the center, and the initial form of a water spike takes shape on the free surface. As the buoyancy parameter increases to 0.2, 0.6, and 0.8, the buoyancy effect becomes significantly stronger, inducing a gradual upward migration of the bubble. This reduces the contact area between the bubble bottom and the bottom boundary. Furthermore, a larger buoyancy parameter elevates the position of the initial jet and the water spike, and correspondingly delays the formation timing of the jet. At time $t_2$, the water spike on the free surface becomes quite pronounced. It can be seen that a larger buoyancy parameter results in a lower water spike height and a lower velocity at its tip. For $\delta =$ 0.0, 0.2, and 0.6, the water jet impacts the bottom boundary, creating a high-pressure zone and generating radial jets along the boundary, this phenomenon becomes less distinct as the buoyancy parameter increases. However, when $\delta =0.8$, the strong buoyancy effect causes the bubble to neck and eventually split between the free surface and the bottom boundary. Subsequently, upward and downward jets are formed vertically, acting on the respective boundaries, thereby subjecting the free surface to the high-pressure load of jet impact as well. Overall, a larger buoyancy parameter causes the bubble pulsation closer to the middle region of the free surface and the bottom boundary. Times $t_3$ and $t_4$ illustrate the evolution patterns of the bubble pulsation and water spike in the later stages. Influenced by the different pressure and velocity fields generated after the first bubble pulsation, distinct wrinkles form at the root of the water spike by time $t_3$. Early wrinkles near the water spike splash and develop into the water skirt, and the distance of this skirt from the axis of the water spike increases with a larger buoyancy parameter. For buoyancy parameters of 0.6 and 0.8, the bubble’s absolute size is increased, and the buoyancy effect is enhanced. While pulsating, the bubble also migrates upward towards the water surface, carrying more energy. This promotes the faster transformation of early wrinkles into the water skirt and encourages the generation of new wrinkles between the water skirt and the main water spike. When $\delta =0.8$, the new wrinkles between the water skirt and the water spike exhibit a more distinct vertically upward velocity. Simultaneously, a larger buoyancy parameter causes height of the depression between the water spike and the water skirt to decrease, bringing it closer to the bottom bubble. This makes the region more susceptible to the influence of the high-pressure zones generated by bubble pulsation, facilitating the transformation of the new wrinkles into the new water skirt. By time $t_4$, the new water skirt in Case 12 has already impinged upon the main water spike at high speed.

As can be seen from Figure 17, before the jet penetrates the bubble, the jet velocity gradually increases while the jet width gradually decreases as the buoyancy parameter increases. Since the water depth parameter and the position parameter remain constant, the relative positions of the bubble to the free surface and the bottom boundary are unchanged. With the fluid Reynolds number held constant, the increase in the buoyancy parameter causes the inertial effects to become dominant. This results in a downward jet that is more concentrated near the axis, attaining higher velocities and a narrower jet width.

Figure 18 shows the variation of water spike height over time for different buoyancy parameters. At the initial stage, the changes in water spike height are largely similar with different buoyancy parameters. As time goes by, a larger buoyancy parameter corresponds to a lower water spike height and an earlier onset of the falling-back phenomenon. For $\delta =$ 0.0, 0.2, and 0.4, the water spike continues to rise at a certain speed, with the ascent velocity increasing as the buoyancy parameter decreases. When $\delta =0.6$, the water spike height reaches its maximum value during the second bubble pulsation period and then begins to descend. For $\delta =0.8$, the water spike height starts to decrease within the first bubble pulsation period and falls back to the water surface during the third period. As the buoyancy parameter increases, the absolute size of the bubble also grows,and the energy has a greater impact on the water spike at the free surface, leads to a faster jet velocity. After penetrating the bubble, the jet exerts a downward dragging force on the water spike, prompting it to fall back earlier.

4. Conclusions

This paper establishes an axisymmetric two-dimensional numerical model based on EFEM suitable for shallow water explosions. It systematically investigates the influence of the water depth parameter $\lambda$, position parameter $\gamma$, and buoyancy parameter $\delta$ on the dynamic evolution process of bubble dynamics and free surface structures such as water spike and wrinkles in extremely shallow water explosions, thereby revealing the nonlinear coupling mechanisms among them. The following conclusions can be drawn:

(1) Due to the nonlinear coupling between the bubble, the free surface, and the bottom boundary in the subsequent stages, the significant influence of different parameters on bubble dynamics is primarily manifested in the first few pulsation periods. Overall, with the increase of the water depth parameter $\lambda$ and the position parameter $\gamma$, the maximum equivalent bubble radius increases, and the corresponding pulsation period is prolonged. An increase in the buoyancy parameter $\delta$ also enlarges the maximum equivalent bubble radius but accelerates the bubble collapse process, thereby shortening the pulsation period. It should be noted that when $\delta =0$, its constraining effect on the bubble is reduced, resulting in a maximum equivalent bubble radius slightly higher than that in the scenario with $\delta =0.2$.

(2) Under smaller water depth and position parameters, the underwater explosion bubble is closer to the free surface, which can easily lead to the rupture of both the bubble and the free surface, resulting in energy dissipation. As the water depth parameter $\lambda$ increases, the jet width also increases, while the jet velocity first decreases and then increases. This is due to the progressively enhanced repulsive effect of the free surface and its superposition with the attractive effect of the bottom boundary. For the position parameter $\gamma$, excluding cases of bubble breaking, an increase in $\gamma$ reduces the jet velocity and increases the jet width. A larger buoyancy parameter results in a higher jet velocity and a narrower jet width.

(3) The expansion of a bubble in shallow water explosion generates a water spike at the free surface. Simultaneously, the pulsation process of the bubble induces the formation of wrinkles on the free surface, with early wrinkles splashing to form the water skirt. Excluding cases of bubble breaking, as the three parameters increase, the morphology of the water spike transitions from tall to short and from narrow to wide, while the rate of its fallback accelerates. Furthermore, when the three parameters are relatively large, the distance between the water skirt and the axis of the water spike also increases. As the bubble pulsation continues, more wrinkles form around the water skirt. Among them, new wrinkles emerging between the water spike and the early water skirt also splash, forming a new water skirt. When $\gamma =0.9$, the height of the newly formed water skirt exceeds that of the original water spike and the early water skirt. When $\delta =0.8$, the newly formed water skirt impacts the main water spike at high speed.