Computational domain and boundary conditions
A numerical simulation of the water-entry process was conducted, replicating the conditions of the physical experiment. A two-dimensional numerical model was constructed in the vertical plane, mirroring the layout and dimensions of the physical model. The computational region was defined by the beam, tank sidewalls, tank bottom, and upper boundary. Wall boundary conditions were applied to the beam, tank sidewalls, and tank bottom, while a pressure outlet condition was imposed on the upper boundary.
The computational domain has a length of x/b = 5 and a height of y/b = 3.75. The water depth is y/b = 1.75, and the air domain above it occupies a height of y/b = 2. A structured grid is employed for discretization to ensure both computational accuracy and efficiency. The total number of grid cells in this setup is approximately 300,000, providing a balance between computational precision and practicality. In the x-direction, the grid is refined near critical regions with a grid size of Δx_fine/b = 7.5e-3, while in the majority of the domain, a coarser grid size of Δx_coarse/b = 1.25e-2 is used to maintain computational efficiency. In the y-direction, the grid refinement strategy varies depending on the fluid domain. For the air domain in the y-direction, Δy_fine/b = 1.25e-3 near high-gradient areas, transitioning to Δy_coarse/b = 1.25e-2. In the water domain, Δy_fine/b = 1.25e-3 is maintained near critical features, gradually changing to Δy_coarse/b = 5e-3. Additionally, to verify grid independence, multiple simulations with varying grid densities were conducted. It was found that when the grid size was less than a certain threshold of Δy_fine/b = 1.25e-3 in the most critical regions, the variations in simulation results were within a preset tolerance of 1%. The configuration of the computational domain and simulation grid is illustrated in Fig.8a1 and a2.
In the simulation, the dynamic mesh method was employed to replicate the relative motion between the beam and the tank, while the layering method was utilized for mesh updating. To mimic the flange’s movement, the lower section of the tank was designated as the moving zone, descending at a velocity equivalent to the physical beam’s speed. The tank sidewalls were designated deformation zones, where mesh splitting and deformation processes occurred.
The fluid comprises both water and air, with water being treated as incompressible and air being considered compressible. The fluid dynamics were described using the compressible two-phase form of the Reynolds-averaged Navier–Stokes equations, with turbulence effects accounted for via the renormalization group k–ε turbulence model. The RANS equations, in their differential form, consist of the mass conservation equation Eq.(1), the momentum conservation equations Eq.(2 ~ 3), and the energy conservation equation Eq.(4).
$$ \frac{{\partial \rho }}{{\partial t}} + \nabla \cdot (\rho \overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\rightharpoonup}$}}{V} ) = 0 $$
(1)
$$ \frac{{\partial (\rho u)}}{{\partial t}} + \nabla \cdot (\rho u\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\rightharpoonup}$}}{V} ) = - \frac{{\partial p}}{{\partial x}} + \frac{{\partial \tau _{{xx}} }}{{\partial x}} + \frac{{\partial \tau _{{yx}} }}{{\partial y}} + \rho f_{x} $$
(2)
$$ \frac{{\partial (\rho v)}}{{\partial t}} + \nabla \cdot (\rho v\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\rightharpoonup}$}}{V} ) = - \frac{{\partial p}}{{\partial y}} + \frac{{\partial \tau _{{xy}} }}{{\partial x}} + \frac{{\partial \tau _{{yy}} }}{{\partial y}} + \rho f_{y} $$
(3)
$$ \begin{gathered} \frac{\partial }{{\partial t}}\left[ {\rho \left( {e + \frac{{V^{2} }}{2}} \right)} \right] + \nabla \cdot \left[ {\rho \left( {e + \frac{{V^{2} }}{2}} \right)\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\rightharpoonup}$}}{V} } \right] = \rho \mathop q\limits^{} + \frac{\partial }{{\partial x}}\left( {k\frac{{\partial T}}{{\partial x}}} \right) \hfill \\ \quad + \frac{\partial }{{\partial y}}\left( {k\frac{{\partial T}}{{\partial y}}} \right) - \frac{{\partial \left( {up} \right)}}{{\partial x}} - \frac{{\partial \left( {vp} \right)}}{{\partial y}} \hfill \\ \quad + \frac{{\partial \left( {u\tau _{{xx}} } \right)}}{{\partial x}} + \frac{{\partial \left( {u\tau _{{yx}} } \right)}}{{\partial y}} + \frac{{\partial \left( {v\tau _{{xy}} } \right)}}{{\partial x}} + \frac{{\partial \left( {v\tau _{{yy}} } \right)}}{{\partial y}} \hfill \\ \quad + \rho \overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\rightharpoonup}$}}{f} \cdot \overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\rightharpoonup}$}}{V} \hfill \\ \end{gathered} $$
(4)
where\(\nabla =\frac{\partial }{{\partial x}}\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\rightharpoonup}$}} {i} +\frac{\partial }{{\partial y}}\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\rightharpoonup}$}} {j} \),\( \overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\rightharpoonup}$}}{V} = u\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\rightharpoonup}$}} {i} +v\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\rightharpoonup}$}} {j}\),\( \overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\rightharpoonup}$}}{f} = f_{x} \overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\rightharpoonup}$}} {i} +f_{y}\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\rightharpoonup}$}} {j}\).
The volume of the fluid model was employed to accurately track the interface between water and air, taking into account a density of 998.2kg/m3 for water and using the ideal gas state equation Eq.(5) to determine the density of air.
$$p=\rho RT$$
(5)
For the numerical solution, an explicit time discretization scheme was applied, while spatial discretization employed the second-order upwind scheme. Additionally, the velocity and pressure fields were coupled using the SIMPLE algorithm24.
The time step utilized for our simulations has been carefully chosen with a minimum value of 1e-6 to ensure accuracy and stability of the numerical solutions. Furthermore, to maintain the computational integrity and avoid potential issues associated with unphysical phenomena, we have controlled the courant number strictly within a limit of 2. This approach ensures that the propagation of information through the domain remains physically meaningful and aligned with the characteristics of the fluid flow being simulated.
Computational domain and simulation mesh. (a1) Computational domain; (a2) Simulation mesh.
Validation and verification
Air velocity verification before entry
Throughout the water entry process, we monitored the x-direction flow velocity at positions x1/b = 0.25 and x2/b = 0.5 along the underside of the flange’s leading edge. Prior to the flange’s lower surface making contact with the water, a phenomenon occurs where as the gap between the flange and the water surface diminishes, the air beneath the flange’s underside is expelled horizontally. Assuming a uniform distribution of the air discharge velocity beneath the flange along the y-axis direction and maintaining a stable water surface, we can derive the formula for calculating the air discharge velocity at x1/b and x2/b through a straightforward analysis (Eq.6).
$$v(i,t)/\sqrt {gb} =\frac{{{v_{\text{0}}}{x_{\text{i}}}\sqrt {gb} }}{{h(t)}}$$
(6)
In the equation, v(i, t) represents the flow velocity at point i at time t, and xi denotes the horizontal distance from the point to the center of the flange. v0 represents the water entry velocity, and h(t) represents the distance between the liquid surface and the bottom of the flange at time t.
Figure9 presents a comparison between the numerical simulation results for flow velocities at x1/b and x2/b and the theoretical solution. It is evident from the figures that for tv0/b < 0.575, the numerical solution aligns with the results of the numerical simulation, affirming the accuracy of the numerical model for the wind field. However, for tv0/b > 0.575, the theoretical solution shows a gradual increase, while the numerical solution gradually diminishes and converges toward 0. The discrepancy between the theoretical solution and the numerical model can be attributed to the theoretical assumption that the water surface remains at that level. However, as the gap between the flange and the water surface diminishes, the wind speed at the water surface gradually increases, causing the water surface to become progressively unstable, resulting in the formation of ripples that impede the flow of gases, as depicted in Fig.10. At tv0/b = 0.564, the wind speed at x1/b is 0.429, and at x2/b, it reaches 0.969, with the liquid surface remaining relatively calm. At tv0/b = 0.576 the wind speed at x1/b surges to 1.53, and at x2/b, it increases to 3.1. This leads to the gradual elevation of the water surface along the flange’s central axis to both sides. This phenomenon occurs due to the gradual increase in wind speed from the center toward the periphery, causing a corresponding reduction in pressure at the flange’s lower surface from the center to the edges, resulting in a deformation of the liquid surface. At tv0/b = 0.583, the wind speed at x1/b is 0.661, while at x2/b, it is 0.954. At this point, the flange descends below the average water surface level, and the lower surface gas directly interacts with the water surface from the sides, inducing ripples on the liquid surface. At tv0/b = 0.598, the air beneath the flange’s lower surface becomes trapped, forming an air pocket.
Air discharge velocity at x1/b and x2/b.
Impact stage gas trapping process at the bottom of the flange. (a) tv0/b = 0.564; (b) tv0/b = 0.576 ; (c) tv0/b = 0.583; (d) tv0/b = 0.598.
Water entry process pressure verification
The measured bottom pressure during the water entry process and the corresponding pressure at the numerical simulation points are shown in Fig.11. The graph shows that the numerical simulation pressure trend closely matches that of the physical model. However, there are certain differences in the peak impact force and pressure during the deceleration process. The primary reason for this disparity is that the physical model exhibits a shorter duration but a greater peak during the impact process. Additionally, due to limitations in the stiffness of the supporting structure and hydraulic cylinders in the physical model, energy absorption occurs during impact, resulting from structural deformation and other factors. Therefore, the impact load in the mathematical model is greater than the impact load in the physical model.
Comparison of physical (wp) and numerical (sp) pressures during water entry.
To validate the accuracy of the peak impact force calculations in our model, the same simulation method was applied to numerically simulate the physical measurement of impact pressure during the free-fall water entry process of a flat plate, as reported in reference25. In the referenced literature, a physical model consisting of a flat plate with a length of b = 0.2m and a weight ranging from 15.3 to 15.8kg was released in free fall from a height of h/b = 2 above the water surface into a rectangular water tank with a depth of h/b = 3.
A two-dimensional numerical model was constructed in the vertical plane to reflect the layout and dimensions of the physical model. The calculation area is defined by the flat plate, tank sidewalls, tank bottom and tank upper boundary. Wall boundary conditions are applied to the plate, tank sidewalls and tank bottom, while pressure outlet conditions are applied to the upper boundary. To discretize the computational domain, a structured grid is used in the flat plate area with a grid size of △x/b = 2.5e-3. During the simulation process, the overset grid was used for dynamic grid calculations. The fluids involved in the calculations include water and air, where water is considered incompressible and air is considered compressible. The fluid dynamics are described by the compressible two-phase form of the Reynolds-averaged Navier–Stokes equations, and the turbulence effects are explained by the renormalization group k–ε turbulence model. Considering that the density of water is 998.2kg/m3, a fluid volume model was used to accurately track the water‒air interface, and the ideal gas equation of state was used to determine the density of air.
Pressure measurement points were arranged at the bottom of the model’s longitudinal cross section. The calculated impact pressure comparison is illustrated in Fig.12, where sp represents the numerical simulation results and wp represents the physical measurement results from the referenced literature. The peak impact pressure error is within 10%. Therefore, it can be concluded that utilizing the CFD method for simulating the load during an object’s water entry process is feasible.
Comparison of physical (wp) and numerical (sp) pressures during water entry.
Simulation conditions
The hydrodynamic force experienced by the beam during its water-entry process depended on two key parameters: the flange width (b) and the entry velocity (v). Initially, the water tank had a still water depth of h/b = 1.75, corresponding to a water level of y/b = -0.75, while the bottom of the flange was situated at y/b = 0. The flange traveled d/b = 1.25 at a constant velocity, with a deceleration rate of a/g = 0.01. In total, 25 simulations were conducted, encompassing five different flange widths (0.2, 0.4, 0.6, 0.8, and 1.0m) and five distinct exit velocities ranging from 0.017 to 0.2m/s.
Data analysis
Total hydrodynamic load
The variation in the total hydrodynamic load throughout the water-entry process is illustrated in Fig.13 for gb/v2 = 3394.464. In accordance with the flow patterns observed during the flange water-entry process, the hydrodynamic loads were categorized into three primary stages. During the first stage, the hydrodynamic load experienced a sudden increase followed by a rapid decrease within an extremely short time, forming a pulse peak. This load was denoted as the impact force, which was generated when the flange collided head-on with the water surface during the impact phase. The peak impact was Ft/(0.5ρv2bL) = 1124.5 in the second stage, and the hydrodynamic load gradually increased, resulting in the formation of a second peak, with a peak load Ft/(0.5ρv2bL) = 605.5. The third stage corresponded to the deceleration phase. As the flange decelerated to a halt, the water flow retained its original velocity. The upper water flow exerted a downward thrust on the flange, while the lower water flow adhered to the flange, leading to a sudden reduction in the overall hydrodynamic load, with an amplitude decrease of approximately Ft/(0.5ρv2bL) = 129.7.
Total hydrodynamic load evolution during the water-entry process.
Virtual buoyancy
In this study, virtual buoyancy was defined as the buoyancy exerted on the beam, assuming that the water in the tank was quiet at a given time during the entire water-entry process. Figure14 shows the virtual buoyancy variation corresponding to a test with gb/v2 = 3394.464. The virtual buoyancy force consisted of three straight lines with different slopes, which corresponded to the flange entering the water, the web entering the water, and remaining stationary in the water after entering the water.
Virtual buoyancy evolution during the water entry process.
Additional hydrodynamic loads
The additional hydrodynamic load is defined as the disparity between the total hydrodynamic load and the virtual buoyancy. Figure15 provides an illustration of the additional hydrodynamic load throughout the water-entry process. During the impact stage, the additional hydrodynamic load achieved its initial maximum value of 1124.5, designated the impact load. In the subsequent cavity stage, the additional hydrodynamic load peaked again at 423.9, termed the cavity load. Finally, in the braking stage, the additional hydrodynamic load decreased to a minimum of -86.5, which is designated the braking load. Both the impact load and cavity load held significant importance in the context of facility protection.
Additional hydrodynamic load evolution during the water-entry process.
Prediction of impact loads
In the study of impact forces, Vanden26 noted that the viscosity and gravity of the fluid may have important effects near the contact point. Korobkin & Pukhnachov27 used the perturbation method to also illustrate that the influence of viscosity should be properly considered at the contact point. This study combines previous research and suggests that the impact force Fi/L per unit length of the flange entering the water mainly depends on the following variables:
$${F_{\text{i}}}/L=f(v,b,\rho ,g,\mu )$$
(7)
where ρ is the water density, g is the acceleration due to gravity, and µ is the dynamic viscosity of the water. Dimensional analysis leads to the following mathematical expression:
$$\frac{{{F_{\text{i}}}}}{{\frac{1}{2}\rho {v^2}bL}}=f(\frac{{gb}}{{{v^2}}},\frac{\mu }{{\rho vb}})$$
(8)
Multiple regression for the normalized impact load, Fi/(1/2ρv2bL), yields:
$$\frac{{{F_{\text{i}}}}}{{\frac{1}{2}\rho {v^2}bL}}=33.22{(\frac{{gb}}{{{v^2}}})^{0.38}}{(\frac{\mu }{{\rho vb}})^{ - 0.05}}$$
(9)
With a correlation coefficient of r2 = 0.94. A comparison of the Fi/(1/2ρv2bL) values predicted using Eq.(4) and the calculated values is illustrated in Fig.16. The results were in close agreement with the line of perfect agreement, and the maximum relative deviation from the calculated values was approximately 10%.
Calculated versus predicted normalized impact loads.
Prediction of the cavity loads
Hydraulic considerations have demonstrated that the cavity load per unit length, Fc/L, depends mainly on the following variables26:
$${F_{\text{c}}}/L=f(v,b,\rho ,g)$$
(10)
Dimensional analysis leads to the following mathematical expression:
$$\frac{{{F_{\text{c}}}}}{{\frac{1}{2}\rho {v^2}bL}}=f(\frac{{gb}}{{{v^2}}})$$
(11)
Multiple regression for the normalized cavity load, Fc/(1/2ρv2bL), yields:
$$\frac{{{F_{\text{c}}}}}{{\frac{1}{2}\rho {v^2}bL}}=2.07{(\frac{{gb}}{{{v^2}}})^{0.66}}$$
(12)
with a correlation coefficient of r2 = 0.99. A comparison of the Fc/(1/2ρv2bL) values predicted using Eq.(7) and the calculated values is illustrated in Fig.17. The results were in close agreement with the line of perfect agreement, and the maximum relative deviation from the calculated values was approximately 10%.
Calculated versus predicted normalized cavity loads.
