2.2.1. Penalty IB method

In the penalty IB method (Kim Reference Kim2003; Kim & Peskin Reference Kim and Peskin2016), also known as the feedback forcing method (Goldstein et al. Reference Goldstein, Handler and Sirovich1993), the coupling force between the fluid and the body can be determined by

(2.5)\begin{equation} {\boldsymbol{F}}\left( {{\boldsymbol{X}},t} \right) = K \int_0^t {\left( {{\boldsymbol{V}}\left( {{\boldsymbol{X}},\tau} \right) - {\boldsymbol{v}}\left( {{\boldsymbol{X}},\tau} \right)} \right){\rm{d}} \tau} + C \left( {{\boldsymbol{V}}\left( {{\boldsymbol{X}},t} \right) - {\boldsymbol{v}}\left( {{\boldsymbol{X}},t} \right)} \right), \end{equation}

where ${\boldsymbol {V}}( {{\boldsymbol {X}}, \tau } )$ and ${\boldsymbol {V}}( {{\boldsymbol {X}},t} )$ denote the velocity of boundary node $\boldsymbol {X}$, ${\boldsymbol {v}}( {{\boldsymbol {X}}, \tau } )$ and ${\boldsymbol {v}}( {{\boldsymbol {X}},t} )$ denote the interpolated fluid velocity at boundary node $\boldsymbol {X}$, i.e. (2.4). Here, $K$ and $C$ are free constants according to the problem being solved. For the inviscid fluid and the rigid body, the free-slip conditions hold at the interface. The boundary coupling force is exerted only in the normal direction of the boundary $\varGamma$, as shown in figure 4(a), and is given as

(2.6)\begin{align} {\boldsymbol{F}}\left( {{\boldsymbol{X}},t} \right) &= \left\{ K \Delta t\left[ {{\boldsymbol{V}}\left( {{\boldsymbol{X}},t} \right) - {\boldsymbol{v}}\left( {{\boldsymbol{X}},t} \right)} \right] \boldsymbol{\cdot} {\boldsymbol{n}}\left( {{\boldsymbol{X}},t} \right) \right.\nonumber\\ &\quad +\left. C \left[ {\boldsymbol{V}}\left( {{\boldsymbol{X}},t} \right) - {\boldsymbol{v}}\left( {{\boldsymbol{X}},t} \right) \right]\boldsymbol{\cdot} {\boldsymbol{n}}\left( {{\boldsymbol{X}},t} \right) \right\}{\boldsymbol{n}}\left( {{\boldsymbol{X}},t} \right) , \end{align}

where $\boldsymbol {n}$ is the unit normal vector of the boundary $\varGamma$, and $\Delta t$ is the time increment. From a physical point of view, it can be regarded as a damped simple harmonic oscillator placed between the boundary and the fluid interface for the consistent motion of the IB (see figure 4b). Hence, $K$ is regarded as the spring stiffness, and $C$ is the damping coefficient. Aquelet, Souli & Olovsson (Reference Aquelet, Souli and Olovsson2006) proposed a coefficient formula to define the two coefficients $K$ and $C$, and the formula for $K$ is expressed as

(2.7)\begin{equation} K=K\left( {{\boldsymbol{X}},t} \right) = {p_f}\frac{{{{\left[ {\mathbb{S}\left( {{\boldsymbol{X}},t} \right)} \right]}^2}}}{{\mathbb{V}\left( {{\boldsymbol{X}},t} \right)}}{K_b}\left( {{\boldsymbol{X}},t} \right) , \end{equation}

where $K_b$ is the bulk modulus of the fluid cell containing the boundary node $\boldsymbol {X}$, and $\mathbb {V}$ is the volume of the fluid cell; $\mathbb {S}$ is the average area of the body surface grids connected to the boundary node $\boldsymbol {X}$; $p_f$ is an empirical coefficient with a range of $0 < {p_f} < 1$. A larger $p_f$ value may cause instabilities, and $p_f$ is generally set to 0.1. Also, $C$ is taken as the critical damping coefficient of the spring system to eliminate numerical oscillations and is calculated by

(2.8)\begin{equation} C=C\left( {{\boldsymbol{X}},t} \right) = 2\sqrt {K\left( {{\boldsymbol{X}},t} \right) \bar{m}\left( {{\boldsymbol{X}},t} \right)} , \end{equation}

where $\bar {m}( {{\boldsymbol {X}},t} )$ is the fluid mass obtained by interpolation at the boundary node, that is,

(2.9)\begin{equation} \bar{m}\left( {{\boldsymbol{X}},t} \right) = \sum\limits_\varOmega {m({\boldsymbol{x}},t)\delta ({\boldsymbol{x}},{\boldsymbol{X}},t)}, \end{equation}

where $m({\boldsymbol {x}},t)$ is concentrated nodal mass corresponding to the component of the lumped mass matrix described in (2.37) in § 2.3, which is the average mass of the connected fluid cells.

Figure 4. Sketch of coupling force calculation in the IB method. (a) The IB (red solid line) and the fictitious fluid interface (green solid line) coincide at $t=t_n$. After one time increment $\Delta t$ without FSI treatment, they (dashed line) deviate from each other with the respective velocities $\boldsymbol {V}$ and $\boldsymbol {v}$. A pair of equal and opposite forces $\boldsymbol {F}$ and $-\boldsymbol {F}$ in the normal direction to the surface are exerted on the fictitious fluid interface and IB, respectively, to counteract their deviation. (b) In the penalty IB method, the coupling force $\boldsymbol {F}$ is determined by a damped oscillator with a zero resting length, which is used to connect the boundary node (red hollow dot) and fictitious fluid point (green solid dot). Here, $K$ and $C$ are artificial empirical coefficients and are regarded as the spring stiffness and damping coefficient from a physical point of view. (c) In the present IB method, the coupling force $\boldsymbol {F}$ is directly calculated based on the velocity difference $\Delta V$ between the boundary node and the fictitious fluid point, where $\bar {m}( {{\boldsymbol {X}},t} )$ is the fluid mass obtained by interpolation at the boundary node.

The resultant force ${{\boldsymbol {F}}_f}$ and moment ${{\boldsymbol {T}}_f}$ exerted by the fluid on the rigid body are given by

(2.10)\begin{equation} \left.\begin{gathered} \displaystyle {{\boldsymbol{F}}_f}\left( t \right) =- \sum\limits_\varGamma {{\boldsymbol{F}} \left( {{\boldsymbol{X}},t} \right)}, \\ \displaystyle {{\boldsymbol{T}}_f}\left( t \right) =- \sum\limits_\varGamma {\left( {{\boldsymbol{X}} - {\boldsymbol{G}}} \right) \times {\boldsymbol{F}}\left( {{\boldsymbol{X}},t} \right)} . \end{gathered}\right\} \end{equation}

With its inherent advantage of momentum conservation, the penalty IB method has been extensively applied to the study of FSIs by many researchers. The details of the numerical scheme can be found in Aquelet et al. (Reference Aquelet, Souli and Olovsson2006), Wang & Guedes Soares (Reference Wang and Guedes Soares2014) and Tian et al. (Reference Tian, Zhang, Liu and Wang2021b).

2.2.2. Improved IB method

Since the violent water impact is usually inertia dominated, fluid viscosity is currently ignored. Therefore, the free-slip condition (Anderson Reference Anderson1995) holds on the IB for the inviscid fluid and the rigid body. On the IB $\varGamma$, it satisfies

(2.11)\begin{equation} \Delta V \left( {{\boldsymbol{X}},t} \right) = \left[ {{\boldsymbol{V}}\left( {{\boldsymbol{X}},t} \right) - {\boldsymbol{v}}\left( {{\boldsymbol{X}},t} \right)} \right] \boldsymbol{\cdot} {\boldsymbol{n}}\left( {{\boldsymbol{X}},t} \right) = 0 . \end{equation}

Without the FSI treatment, $\Delta V \ne 0$ is self-evident at the next time increment, as shown in figure 4(a). To eliminate the velocity difference $\Delta V$ between the boundary node and the fictitious fluid point (see figure 4c), the novel coupling force $\tilde {\boldsymbol {F}}$ at the boundary node $\boldsymbol {X}$ in one time increment $\Delta t$ is directly defined as

(2.12)\begin{equation} \tilde{\boldsymbol{F}}\left( {{\boldsymbol{X}},t} \right) = {\bar{m}} \left( {{\boldsymbol{X}},t} \right) {\boldsymbol{A}}\left( {{\boldsymbol{X}},t} \right) , \end{equation}

where $\bar {m}( {{\boldsymbol {X}},t} )$ is the fluid mass calculated by (2.9), and ${\boldsymbol {A}}( {{\boldsymbol {X}},t} )$ is defined as

(2.13)\begin{equation} {\boldsymbol{A}}\left( {{\boldsymbol{X}},t} \right) = \frac{{\Delta V ^* \left( {{\boldsymbol{X}},t} \right)}}{{\Delta t}}{\boldsymbol{n}}\left( {{\boldsymbol{X}},t} \right) = \frac{{\left[ {{\boldsymbol{V} ^*}\left( {{\boldsymbol{X}},t} \right) - {\boldsymbol{v} ^*}\left( {{\boldsymbol{X}},t} \right)} \right] \boldsymbol{\cdot} {\boldsymbol{n}}\left( {{\boldsymbol{X}},t} \right)}}{{\Delta t}}{\boldsymbol{n}}\left( {{\boldsymbol{X}},t} \right) , \end{equation}

where the superscript $*$ represents the variables calculated at the next time increment without FSI treatment. Namely, ${\boldsymbol {v} ^*}( {{\boldsymbol {X}},t} )$ denotes the interpolated fluid velocity by solving the momentum equation without taking into account the body-force field. However, since the coupling force is zero, ${\boldsymbol {V} ^*}( {{\boldsymbol {X}},t} )$ is replaced by the boundary node velocity at the current moment, i.e. ${\boldsymbol {V}}( {{\boldsymbol {X}},t} )$. Thus,

(2.14)\begin{equation} \tilde{\boldsymbol{F}}\left( {{\boldsymbol{X}},t} \right) = {\bar{m}} \left( {{\boldsymbol{X}},t} \right) \frac{{\left[ {{\boldsymbol{V}}\left( {{\boldsymbol{X}},t} \right) - {\boldsymbol{v} ^*}\left( {{\boldsymbol{X}},t} \right)} \right] \boldsymbol{\cdot} {\boldsymbol{n}}\left( {{\boldsymbol{X}},t} \right)}}{{\Delta t}}{\boldsymbol{n}}\left( {{\boldsymbol{X}},t} \right) . \end{equation}

In the conventional IB method, the coupling force $\tilde {\boldsymbol {F}}( {{\boldsymbol {X}},t} )$ is distributed to the surrounding fluid grid nodes through the Dirac delta function, as in (2.3), that is,

(2.15)\begin{equation} {\boldsymbol{f}}({\boldsymbol{x}},{\boldsymbol{X}},t) = \tilde{\boldsymbol{F}}\left( {{\boldsymbol{X}},t} \right) \delta ({\boldsymbol{x}},{\boldsymbol{X}},t) = \delta ({\boldsymbol{x}},{\boldsymbol{X}},t) {\bar{m}} \left( {{\boldsymbol{X}},t} \right) {\boldsymbol{A}}\left( {{\boldsymbol{X}},t} \right) . \end{equation}

The Dirac delta function is dependent only on the distance parameters. When the fluid cells near the IB contain free surfaces or multiphase interfaces, the conventional methods may cause discontinuities in the flow velocity due to the uneven mass distribution. This will cause non-physical force oscillations, numerical instability and even interface disturbance.

Substituting (2.9) into (2.12), we obtain

(2.16)\begin{equation} \tilde{\boldsymbol{F}}\left( {{\boldsymbol{X}},t} \right) = {\boldsymbol{A}}\left( {{\boldsymbol{X}},t} \right) \sum\limits_\varOmega {m({\boldsymbol{x}},t)\delta ({\boldsymbol{x}},{\boldsymbol{X}},t)} . \end{equation}

It is noticeable that (2.16) contains the summation term for the mass of surrounding fluid grid nodes. With this feature, we can assign the force to the surrounding fluid nodes as

(2.17)\begin{equation} \tilde{\boldsymbol{f}}({\boldsymbol{x}},{\boldsymbol{X}},t) = m\left( {{\boldsymbol{x}},t} \right){\boldsymbol{A}}\left( {{\boldsymbol{X}},t} \right)\delta \left( {{\boldsymbol{x}},{\boldsymbol{X}},t} \right). \end{equation}

By combining (2.16) and (2.17), we can obtain

(2.18)\begin{align} \tilde{\boldsymbol{F}}({\boldsymbol{X}},t) = {\boldsymbol{A}}\left( {{\boldsymbol{X}},t} \right) \sum\limits_\varOmega { m({\boldsymbol{x}},t)\delta ({\boldsymbol{x}},{\boldsymbol{X}},t)} = \sum\limits_\varOmega {{\boldsymbol{A}}\left( {{\boldsymbol{X}},t} \right) m({\boldsymbol{x}},t)\delta ({\boldsymbol{x}},{\boldsymbol{X}},t)} = \sum\limits_\varOmega {\tilde{\boldsymbol{f}}({\boldsymbol{x}},{\boldsymbol{X}},t)}, \end{align}

which indicates that the novel force distribution method still satisfies linear momentum conservation.

Integrating $\tilde {\boldsymbol {f}}({\boldsymbol {x}},{\boldsymbol {X}},t)$ over the IB $\varGamma$, we obtain the body-force field exerted on the fluid:

(2.19)\begin{equation} \tilde{\boldsymbol{f}}\left( {{\boldsymbol{x}},t} \right) =\sum\limits_\varGamma { \tilde{\boldsymbol{f}} ({\boldsymbol{x}},{\boldsymbol{X}},t)}. \end{equation}

The ratio of the fluid grid size to the boundary grid size in the IB method is generally greater than one to ensure the robustness of the numerical scheme and is set to approximately two in the present numerical scheme. The term ${{{{[ {\mathbb {S}( {{\boldsymbol {X}},t} )} ]}^2}}}/{{\mathbb {V}( {{\boldsymbol {X}},t} )}}$ in (2.7) represents the correction for the difference in grid sizes. Similarly, the dimensionless coefficient $\lambda$ is defined to consider the discrepancy between the fluid grids and boundary grids, that is,

(2.20)\begin{equation} \lambda \left( {{\boldsymbol{X}},t} \right) = \frac{{\mathbb{S}\left( {{\boldsymbol{X}},t} \right)}}{{{{\left[ {\mathbb{V}\left( {{\boldsymbol{X}},t} \right)} \right]}^{2/3}}}} . \end{equation}

Adding the dimensionless coefficient $\lambda$ to (2.19), the novel fluid body-force field can be obtained as

(2.21)\begin{align} \tilde{\boldsymbol{f}}\left( {{\boldsymbol{x}},t} \right) & =\sum\limits_\varGamma {\lambda \left( {{\boldsymbol{X}},t} \right) \tilde{\boldsymbol{f}}({\boldsymbol{x}},{\boldsymbol{X}},t)}\nonumber\\ & = m\left( {{\boldsymbol{x}},t} \right) \sum\limits_\varGamma { \frac{{ \mathbb{S}\left( {{\boldsymbol{X}},t} \right)}}{{{{\left[ {\mathbb{V}\left( {{\boldsymbol{X}},t} \right)} \right]}^{2/3}}}} \frac{{\left[ {{\boldsymbol{V}}\left( {{\boldsymbol{X}},t} \right) - {\boldsymbol{v} ^*}\left( {{\boldsymbol{X}},t} \right)} \right] \boldsymbol{\cdot} {\boldsymbol{n}}\left( {{\boldsymbol{X}},t} \right)}}{{\Delta t}} {\boldsymbol{n}}\left( {{\boldsymbol{X}},t} \right) \delta \left( {{\boldsymbol{x}},{\boldsymbol{X}},t} \right) }. \end{align}

Therefore, the novel coupling force is calculated as

(2.22)\begin{align} \tilde{\boldsymbol{F}}({\boldsymbol{X}},t) &= \lambda \left( {{\boldsymbol{X}},t} \right) {\boldsymbol{A}}\left( {{\boldsymbol{X}},t} \right) \sum\limits_\varOmega {m({\boldsymbol{x}},t)\delta ({\boldsymbol{x}},{\boldsymbol{X}},t)}\nonumber\\ & = \frac{{ \mathbb{S}\left( {{\boldsymbol{X}},t} \right)}}{{{{\left[ {\mathbb{V}\left( {{\boldsymbol{X}},t} \right)} \right]}^{2/3}}}} \frac{{\left[ {{\boldsymbol{V}}\left( {{\boldsymbol{X}},t} \right) - {\boldsymbol{v} ^*}\left( {{\boldsymbol{X}},t} \right)} \right] \boldsymbol{\cdot} {\boldsymbol{n}}\left( {{\boldsymbol{X}},t} \right)}}{{\Delta t}}{\boldsymbol{n}}\left( {{\boldsymbol{X}},t} \right) \sum\limits_\varOmega {m({\boldsymbol{x}},t)\delta ({\boldsymbol{x}},{\boldsymbol{X}},t)}. \end{align}

According to (2.10), the resultant forces and moments on the rigid body can also be obtained as

(2.23)\begin{equation} \left.\begin{gathered} \displaystyle {{\boldsymbol{F}}_f}\left( t \right) =- \sum\limits_\varGamma {\tilde{\boldsymbol{F}}\left( {{\boldsymbol{X}},t} \right)}, \\ \displaystyle {{\boldsymbol{T}}_f}\left( t \right) =- \sum\limits_\varGamma {\left( {{\boldsymbol{X}} - {\boldsymbol{G}}} \right) \times \tilde{\boldsymbol{F}}\left( {{\boldsymbol{X}},t} \right)} .\end{gathered}\right\} \end{equation}

2.2.3. Comparison between penalty IB method and present IB method

Although the calculations of body force in the above two methods are completely different, as one is a spring system and the other is directly calculated based on velocity difference $\Delta V$, they lead to expressions with very similar structures. In the penalty IB method, the body force $\boldsymbol {f}$ is deduced by (2.3) and (2.6), and the formula can be reorganized as

(2.24)\begin{equation} {\boldsymbol{f}}\left( {{\boldsymbol{x}},t} \right) = \sum\limits_\varGamma {\left[ {K\left( {{\boldsymbol{X}},t} \right)\Delta t + C\left( {{\boldsymbol{X}},t} \right)} \right]\Delta V\left( {{\boldsymbol{X}},t} \right) {\boldsymbol{n}}\left( {{\boldsymbol{X}},t} \right)\delta \left( {{\boldsymbol{x}},{\boldsymbol{X}},t} \right)} . \end{equation}

In the improved IB method, the body force $\tilde {\boldsymbol {f}}$ is rewritten in a similar form as

(2.25)\begin{equation} \tilde{\boldsymbol{f}}\left( {{\boldsymbol{x}},t} \right) = m\left( {{\boldsymbol{x}},t} \right) \sum\limits_\varGamma {\frac{{\lambda \left( {{\boldsymbol{X}},t} \right)}}{{\Delta t}}\Delta V ^* \left( {{\boldsymbol{X}},t} \right) {\boldsymbol{n}}\left( {{\boldsymbol{X}},t} \right)\delta \left( {{\boldsymbol{x}},{\boldsymbol{X}},t} \right)} . \end{equation}

Observing the two expressions for the body force, it can be found that both of them are proportional to velocity difference $\Delta V$, but the coefficients are different. This also conforms to the core of the IB method, which is that the fluid velocity near the IB is reconstructed with the body force added into the momentum equation of the fluid for consistent motion of the FSI interface.

In the classical IB method, the coupling force is distributed to the surrounding fluid grid nodes through the Dirac delta function as (2.3). When the fluid is a single homogeneous material, it works well. However, for multiphase flows with IB, the mass difference of fluid nodes at the interface, i.e. the fluid cell containing the moving interface, as shown in figure 3, may cause excessive nodal acceleration. The discontinuity in the velocity of the flow field inevitably leads to non-physical force oscillation. For example, for the high-speed water entry in the present work, due to the high density ratio of water and air at the interface, the concentrated nodal mass may differ by hundreds of times. The general approaches of the smooth distribution function and the implicit scheme usually cannot deal with this issue of uneven mass distribution. Tian et al. (Reference Tian, Zhang, Liu and Wang2021b) corrected (2.3) by adding the term of mass weight, that is,

(2.26)\begin{equation} {\boldsymbol{f}}\left( {{\boldsymbol{x}},t} \right) = \sum\limits_\varGamma {{\boldsymbol{F}}\left( {{\boldsymbol{X}},t} \right)\frac{{m\left( {{\boldsymbol{x}},t} \right)\delta \left( {{\boldsymbol{x}},{\boldsymbol{X}},t} \right)}}{{\displaystyle\sum\nolimits_k {m\left( {{{\boldsymbol{x}}_k},t} \right)\delta \left( {{{\boldsymbol{x}}_k},{\boldsymbol{X}},t} \right)} }}}, \end{equation}

where the subscript $k$ indicates the nodes of the fluid grid overlapping the boundary node. The accelerations induced by body force for Tian's method and the improved IB method are compared in (2.27) and (2.28):

(2.27)$$\begin{gather} {\boldsymbol{a}}\left( {{\boldsymbol{x}},t} \right) =\frac{{\boldsymbol{f}}\left( {{\boldsymbol{x}},t} \right)}{m\left( {{\boldsymbol{x}},t} \right)} = \sum\limits_\varGamma {{\boldsymbol{F}}\left( {{\boldsymbol{X}},t} \right)\frac{{ \delta \left( {{\boldsymbol{x}},{\boldsymbol{X}},t} \right)}}{{\displaystyle\sum\nolimits_k {m\left( {{{\boldsymbol{x}}_k},t} \right)\delta \left( {{{\boldsymbol{x}}_k},{\boldsymbol{X}},t} \right)} }}}, \end{gather}$$

(2.28)$$\begin{gather}\tilde{\boldsymbol{a}}\left( {{\boldsymbol{x}},t} \right) =\frac{\tilde{\boldsymbol{f}}\left( {{\boldsymbol{x}},t} \right)}{m\left( {{\boldsymbol{x}},t} \right)} = \sum\limits_\varGamma {\frac{{\lambda \left( {{\boldsymbol{X}},t} \right)}}{{\Delta t}}\Delta V^* \left( {{\boldsymbol{X}},t} \right) {\boldsymbol{n}}\left( {{\boldsymbol{X}},t} \right)\delta \left( {{\boldsymbol{x}},{\boldsymbol{X}},t} \right)} . \end{gather}$$

We find that $\tilde {\boldsymbol {a}}$ in (2.28) is dependent only on the velocity difference $\Delta V$ and not the concentrated nodal masses of the fluid grid. In addition, Tian's method is based on the penalty IB method and indeed causes non-physical force oscillation. The applicability of Tian's method in the problem of high-speed water impact needs to be examined further (see § 3 for the detailed discussion). Additionally, it is also dependent on the empirical coefficients, which is crucial for the stability and robustness of the numerical scheme.

The present IB method provides two major advantages over the previous penalty IB methods: it can tackle the FSI problems of multiphase flow with high density ratios, and also removes the uncertainties related to the artificial coefficients in the penalty IB methods. Therefore, the improved IB method has good applicability and robustness. The comparative analysis of the above methods for high-speed water entry is detailed in § 3.

2.2.4. Internal treatment of the body

In the present numerical scheme, there are no special treatments for the fluid inside the IB, leaving the inner fluid free to develop. Previous research indicated that the fluid inside the closed IB has an influence on the dynamics of the rigid body (Suzuki & Inamuro Reference Suzuki and Inamuro2011). In all cases of water entry carried out in the present work, the inner fluid is set as air, of which the density is much smaller than that of the solid. Therefore, in this case, the influence is negligible.

2.3.1. Governing equations of fluid

According to previous research on the hydrodynamics of high-speed water entry, the characteristics of high velocity and transient and strong impact at the initial stage are observed, so the fluid viscosity, surface tension and heat conduction can be ignored. Therefore, the Cauchy stress tensor ${\boldsymbol {\sigma }}$ in the Navier–Stokes equations (Anderson Reference Anderson1995) can be simplified to its isotropic part $- p{\boldsymbol {\delta }}$, that is, the nonlinear Euler equations of inviscid and initially irrotational compressible single-phase flow:

(2.29)\begin{equation} \frac{{\partial {\boldsymbol{E}}}}{{\partial t}} + \boldsymbol{\nabla} \boldsymbol{\cdot} \left( {{\boldsymbol{Ev}}} \right) = {\boldsymbol{S}}, \end{equation}

where $\boldsymbol {E}$ is the conservation variable of the fluid, $\boldsymbol {v}$ is the fluid velocity and on the right side of the equation, $\boldsymbol {S}$ is the source term. In the coordinate system of $( {xyz} )$, as shown in figure 2, the conservation variable of the fluid and the source term are

(2.30)\begin{equation} {\boldsymbol{E}} = \left[ {\begin{array}{@{}c@{}} \rho \\ {\rho {v_i}}\\ {\rho {e_{{{in}}}}} \end{array}} \right], \quad{\boldsymbol{S}} = \left[ {\begin{array}{@{}c@{}} 0\\ { - (\boldsymbol{\nabla} p)_i + \rho {g_i} + {f_i}}\\ { - p\boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol{v}}} \end{array}} \right] ,\quad i = 1, 2, 3, \end{equation}

where the subscript $i$ represents three components, corresponding to the three axis directions of the coordinate system $( {xyz} )$. Also, $\rho$ is the fluid density, ${e_{{{in}}}}$ is the internal energy per unit mass of the fluid, $p$ is the pressure, $v_i$ is the component of fluid velocity $\boldsymbol {v}$ and $g_i$ is the component of gravitational acceleration $\boldsymbol {g}$. The gravity direction is specified along the axial direction of $z$, and $\boldsymbol {g} = ( 0, 0, -9.8)\ \textrm {{m}}\ \textrm {{s}}^{-2}$. In addition, $f_i$, the component of ${\boldsymbol {f}}$, is the body force from the IB.

The volume of fluid method (Hirt & Nichols Reference Hirt and Nichols1981) is used to deal with multiphase flow in the process of water entry, and the advection equation is

(2.31)\begin{equation} \frac{{\partial \alpha_j}}{{\partial t}} + {\boldsymbol{v}} \boldsymbol{\cdot} \boldsymbol{\nabla} \alpha_j = 0 , \end{equation}

where $\alpha _j$ is the volume fraction of the fluid phase and satisfies $0 \leqslant {\alpha _j} \leqslant 1$ and $\sum \nolimits _j {{\alpha _j}} = 1$ in one fluid cell. Here, $j$ is the fluid phase number, i.e. $j=1,2$, representing the water phase and the air phase, respectively. In a mixed fluid cell, the fluid is assumed to be a homogeneous mixture of fluid phases, and the fluid density and pressure are computed as a function of $\bar {\rho } = \sum {{\alpha _j}{\rho _j}}$ and $\bar {p} = \sum {{\alpha _j}{p_j}}$.

The advection equation (2.31) is reorganized and added to (2.29) to compose the Euler equations for multiphase compressible flow, in which the conservation variable of the fluid $\boldsymbol {E}$ and source term $\boldsymbol {S}$ are

(2.32)\begin{align} {\boldsymbol{E}} = \left[ {\begin{array}{@{}c@{}} {{\alpha _j}}\\ {{\alpha _j}{\rho _j}}\\ {\bar{\rho} {v_i}}\\ {{\alpha _j}{\rho _j}{e_{{{in}}_j}}} \end{array}} \right],\quad {\boldsymbol{S}} = \left[ \begin{array}{@{}c@{}} {\alpha _j} \boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol{v}}\\ 0\\ - (\boldsymbol{\nabla} \bar{p})_i + \bar{\rho}{g_i} + {f_i}\\ - {\alpha _j}\bar{p} \boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol{v}} \end{array} \right],\quad i = 1, 2, 3,\ j = 1,2 . \end{align}

In (2.29), the number of unknown variables in the equations is greater than the number of equations. Therefore, to make the Euler equations closed, the equation of state (EOS) is used. It is assumed that air is an ideal gas, and the $\gamma$-law EOS (Fedkiw et al. Reference Fedkiw, Aslam, Merriman and Osher1999) used for an ideal gas is given as

(2.33)\begin{equation} p_g = \rho_g {e_{{{in}}_g}} \left( {\gamma_g - 1} \right) , \end{equation}

where the subscript $g$ indicates the gas phase; $\gamma _g$ is the ratio of the specific heats. For air, the reference densities are $\rho _{g_0} =1.29\ \textrm {{kg}}\ \textrm {{m}}^{-3}$ and $\gamma _g=1.4$. The EOS commonly used for water includes the Tait EOS, stiffened gas (SG) EOS and Mie–Grüneisen EOS. The SG EOS (Saurel et al. Reference Saurel, Le Métayer, Massoni and Gavrilyuk2007) is used in this work, and it can be expressed as

(2.34)\begin{equation} p_l = \rho_l {e_{{{in}}_l}} \left( {\gamma_l - 1} \right) - \gamma_l {P_w} , \end{equation}

where the subscript $l$ indicates the liquid phase; $\gamma _l$ and $P_w$ are parameters obtained from the shock Hugoniot experiment. Typically, for water, the reference density $\rho _{l_0} =1000.0\ \textrm {{kg}}\ \textrm {{m}}^{-3}$, $\gamma _l=7.15$ and $P_w=3.309 \times 10^8$ Pa.

Thus, (2.29) and (2.32)–(2.34) constitute the governing equation system and can be solved for compressible multiphase flow.

2.3.2. Multiphase Eulerian finite element method

In CFD techniques, there are many methods to solve the above governing equations of fluid motion, such as SPH, FVM and finite difference method (FDM). In this work, the multiphase EFEM (Benson Reference Benson1992a; Benson & Okazawa Reference Benson and Okazawa2004; Liu et al. Reference Liu, Zhang, Tian and Wang2019) based on the operator splitting technique (OST) is used to solve the above system of (2.29). The core idea is to separate the convection term from the Euler equations by using operator splitting to avoid the numerical instability caused by handling the convection term using the Galerkin method. With the OST treatment, (2.29) is divided into two equations, that is, (2.35a) and (2.35b), which are solved sequentially in one time increment:

In the first step, (2.35a) can be solved after some mathematical transformations. Substituting the continuity equation into the momentum equation and the energy equation and after some manipulations, (2.35a) can be transformed into

(2.36) \begin{equation} \left. {\begin{array}{l@{}} \displaystyle \dfrac{\partial {\rho _j}}{{\partial t}} =- {\rho _j} \boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol{v}} \\ \displaystyle \bar{\rho}\dfrac{{\partial {v _i} }}{{\partial t}} =- (\boldsymbol{\nabla} \bar p)_i + \bar{\rho} {g_i} + {f_i} \\ \displaystyle {\rho _j} \dfrac{{\partial {e_{{{in}}_j}} }}{{\partial t}} =- \bar{p} \boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol{v}} \end{array}},\quad i = 1, 2, 3,\ j = 1,2.\right\} \end{equation}

It can be found that the only difference between (2.36) and the standard Lagrangian formulation is the type of the time derivative. Therefore, it can be solved by the classical explicit finite element method. In finite element formulations, the velocity component $v_i$ is known at the nodes, while the variables of $\alpha _j$, $\rho _j$ and $e_{{{in}}_j}$ are at the centre of the elements. The accelerations of the nodes are calculated by solving the momentum equation though the Galerkin method (Wu & Gu Reference Wu and Gu2012), and its discrete form is

(2.37)\begin{align} &\sum\limits_k { \left( \int_{{\varOmega ^k}} {{{\bar{\rho} }^k}\varPhi _M^k\varPhi _N^k\,{\rm{d}}{\varOmega ^k}} \right) } \dot{v}_i \nonumber\\ &\quad =\sum\limits_k \left[ \int_{{\varOmega ^k}} {{{\bar{\rho} }^k}\varPhi _M^k{g_i}\,{\rm{d}}{\varOmega ^k}} + \int_{{\varOmega ^k}} {{{\bar{p}}^k}{{( {\boldsymbol{\nabla} \varPhi _M^k})}_i}\,{\rm{d}}{\varOmega ^k}} - \int_{{\varGamma ^k}} {\bar{p}\varPhi _M^k{{\boldsymbol{n}}_i}\,{\rm{d}}{\varGamma ^k}} \right] + {f_i} , \end{align}

where the superscript $k$ represents the fluid element number, $M$ and $N$ are both the total number of nodes of the fluid element $\varOmega$, $\varPhi _M^k$ and $\varPhi _N^k$ are the shape functions of element $k$, $\varOmega ^k$ denotes the element $k$, $\varGamma ^k$ indicates the element boundary of $k$ and $\boldsymbol {n}_i$ represents the component of the outer unit normal of $\varGamma ^k$ in the direction $i$. Then, the velocity field is computed and updated, so the density and energy of each phase are updated by

(2.38) \begin{equation} \left. {\begin{array}{c@{}} \displaystyle {\dot{\rho}_j^k} =- \dfrac{{\rho _j^k}}{{{V^k}}}\int_{{\varGamma ^k}} {\varPhi _N^k {{\boldsymbol{v}}} \boldsymbol{\cdot} {\boldsymbol{n}}}\, {\rm{d}}{\varGamma ^k}\\ \displaystyle {\dot{e_{{in}}} _j^k} =- \dfrac{{{{\bar{p}}^k}}}{{m_j^k}}\int_{{\varGamma ^k}} {\varPhi _N^k{{\boldsymbol{v}}} \boldsymbol{\cdot} {\boldsymbol{n}}} \,{\rm{d}}{\varGamma ^k} \end{array}},\quad j=1,2,\right\} \end{equation}

where $\boldsymbol {n}$ is a unit vector perpendicular to the external surface of $\varGamma ^k$, $V$ is the element volume and $m_j$ is the mass of the $j$ fluid phase.

In the second step, the form of (2.35b) is consistent with the advection equation. Therefore, the transport flux of the conservation variable $\boldsymbol {E}$ between adjacent elements is calculated to solve (2.35b). Actually, the physical time is not advanced during this step. In each time increment, the transport calculation can be split into a sequence of one-dimensional transports (Benson & Okazawa Reference Benson and Okazawa2004). The structured grids are decomposed into one-dimensional transport calculations in three directions. Before solving the advection equation, the interface of multiphase flow is constructed by the piecewise-linear interface construction algorithm (Rider & Kothe Reference Rider and Kothe1998). Then, the volume of material transported between adjacent elements is calculated, and the element-centred variables (mass and internal energy) and the node-centred variable (momentum) are updated by the monotone upwind scheme conservation laws scheme (Van Leer Reference Van Leer1977) and half-index shift algorithm (Benson Reference Benson1992b), respectively. After advection, the pressure is updated by the EOS, and the time increment is also updated. Then, the numerical implementation enters the next loop.

Regarding the simplicity and stability of the technique, the numerical scheme of EFEM has been extensively discussed by many researchers (Benson Reference Benson1992a, Reference Benson1997; Tian et al. Reference Tian, Zhang, Liu and Tao2021a), so the details will not be discussed here.

2.4.1. Motion equations of rigid body

In this study, the structural deformation of the body during water entry is ignored, and the object is assumed to be a rigid body since the hydrodynamic loads during high-speed water entry are primarily investigated. The local coordinate system fixed on the gravity centre $\boldsymbol {G}$ of the body is established to represent the position and orientation of the rigid body, as shown in figure 2. It is specified that the translational velocity and angular velocity at the gravity centre $\boldsymbol {G}$ are $\boldsymbol {V}_G$ and $\boldsymbol {\omega }$, respectively. The 6-DOF motion equations of a rigid body (Fossen Reference Fossen1994), including translation and rotation equations, are as follows:

(2.39)$$\begin{gather} M( {{\dot{\boldsymbol{V}}}_G^b + {\boldsymbol{\omega }^b} \times {\boldsymbol{V}}_G^b}) = {\boldsymbol{F}}_{0}^b , \end{gather}$$

(2.40)$$\begin{gather}{{\boldsymbol{\mathsf{J}}}_G}{\dot{\boldsymbol{\omega}}^b} + {\boldsymbol{\omega }^b} \times ( {{\boldsymbol{\mathsf{J}}}_G}{\boldsymbol{\omega }^b}) = {\boldsymbol{T}}_{0}^b , \end{gather}$$

where the superscript $b$ represents the variables in the LCS, $M$ is the mass of the rigid body and $\boldsymbol{\mathsf{J}}_G$ is the inertia matrix with respect to the gravity centre $\boldsymbol {G}$ of the rigid body in the LCS,

(2.41)\begin{equation} {{\boldsymbol{\mathsf{J}}}_G} = \left[ {\begin{array}{@{}ccc@{}} {{J_{xx}}} & {{J_{xy}}} & {{J_{xz}}}\\ {{J_{yx}}} & {{J_{yy}}} & {{J_{yz}}}\\ {{J_{zx}}} & {{J_{zy}}} & {{J_{zz}}} \end{array}} \right] , \end{equation}

constructed from the moments $( {{J_{xx}},{J_{yy}},{J_{zz}}} )$ and products $( {{J_{xy}},{J_{xz}},{J_{yz}}} )$ of inertia of the body. In the LCS, the three axes $( {{J_{xy}},{J_{xz}},{J_{yz}}} )$ are principal axes of revolution of the rigid body, so ${J_{xy}} = {J_{xz}} = {J_{yz}} = 0$. Also, ${\boldsymbol {F}}_{0}^b$ and ${\boldsymbol {T}}_{0}^b$ are the resultant forces and moments acting on the rigid body in the LCS, respectively. They are calculated in the GCS and then transformed into the LCS, which are expressed as

(2.42)\begin{equation} \left.\begin{gathered} {\boldsymbol{F}}_{0}^b = {\boldsymbol{\mathsf{R}}}\left( {{\boldsymbol{F}}_{f}} + M{\boldsymbol{g}}\right) ,\\ {\boldsymbol{T}}_{0}^b = {\boldsymbol{\mathsf{R}}}{{\boldsymbol{T}}_{f}}, \end{gathered}\right\} \end{equation}

where $\boldsymbol{\mathsf{R}}$ is the rotation matrix from the GCS to the LCS.

2.4.2. Unit quaternions for 6-DOF system

Before solving the motion equations (2.39) and (2.40), it is necessary to determine the relationship between the LCS and the GCS. Although Euler angles are most commonly used to represent the rigid body motion due to the convenience in implementation and visualization, a phenomenon called gimbal lock sometimes occurs. Unit quaternions (Chelnokov Reference Chelnokov1983; Featherstone Reference Featherstone2008) are used instead of Euler angles to represent spatial orientations and rotations of rigid bodies in three-dimensional space, which have several advantages over Euler angles, including low computation cost, high accuracy, avoidance of trigonometric calculation and the absence of singularity.

Unit quaternions are a normalized set of quaternions and are generally represented in the form

(2.43)\begin{equation} {\boldsymbol{Q}} = {\left( {{q_0},{q_1},{q_2},{q_3}} \right)^{\rm{T}}} = {q_0} + {q_1}{\boldsymbol{i}} + {q_2}{\boldsymbol{j}} + {q_3}{\boldsymbol{k}} , \end{equation}

where ${\boldsymbol {i}},{\boldsymbol {j}},{\boldsymbol {k}}$ are the basic quaternions with $\boldsymbol {i}^2=\boldsymbol {j}^2=\boldsymbol {k}^2=\boldsymbol {i}\boldsymbol{\cdot}\boldsymbol {j}\boldsymbol{\cdot}\boldsymbol {k}=-1$, and ${q_0},{q_1},{q_2},{q_3}$ are real numbers. The four numbers of unit quaternions are called Euler parameters and satisfy

(2.44)\begin{equation} {q_0}^2+{q_1}^2+{q_2}^2+{q_3}^2=1 . \end{equation}

When the LCS coincides with the GCS, ${\boldsymbol {Q}}$ is equal to ${( {1,0,0,0} )^\textrm {{T}}}$. For the initial condition of high-speed water entry in figure 2, the initial water entry angle of the slender body is $\theta _0$. Correspondingly, the initial unit quaternions are set as

(2.45)\begin{equation} {\boldsymbol{Q}} ={\left( \cos \frac{\theta_0}{2}, \sin \frac{\theta_0}{2}, \sin \frac{\theta_0}{2}, \sin \frac{\theta_0}{2} \right)^{\rm{T}}} . \end{equation}

According to the rigid body dynamics of mechanisms (Featherstone Reference Featherstone2008), the rotation matrix $\boldsymbol{\mathsf{R}}$ from the GCS to the LCS can be derived by unit quaternions ${\boldsymbol {Q}}$, that is,

(2.46)\begin{align} {\boldsymbol{\mathsf{R}}} = \left[ {\begin{array}{@{}ccc@{}} {{q_0}^2 + {q_1}^2 - {q_2}^2 - {q_3}^2} & {2 {{q_1}{q_2} + 2 {q_0}{q_3}} } & {2 {{q_1}{q_3} - 2 {q_0}{q_2}} }\\ {2 {{q_1}{q_2} - 2 {q_0}{q_3}}} & {{q_0}^2 - {q_1}^2 + {q_2}^2 - {q_3}^2} & {2 {{q_2}{q_3} + 2 {q_0}{q_1}} }\\ {2 {{q_1}{q_3} + 2 {q_0}{q_2}} } & {2 {{q_2}{q_3} - 2 {q_0}{q_1}} } & {{q_0}^2 - {q_1}^2 - {q_2}^2 + {q_3}^2} \end{array}} \right] . \end{align}

After solving the motion equations (2.39) and (2.40), unit quaternions are updated from $t_n$ to $t_{n+1}$, which is expressed as

(2.47)\begin{equation} {\boldsymbol{Q}}_{t_{n+1}} = {\boldsymbol{Q}}_{t_{n}} + \frac{1}{2} \Delta t \left[ {\begin{array}{@{}cccc@{}} 0 & { - {\omega _x}} & { - {\omega _y}} & { - {\omega _z}}\\ {{\omega _x}} & 0 & {{\omega _z}} & { - {\omega _y}}\\ {{\omega _y}} & { - {\omega _z}} & 0 & {{\omega _x}}\\ {{\omega _z}} & {{\omega _y}} & { - {\omega _x}} & 0 \end{array}} \right]{\boldsymbol{Q}}_{t_{n}} , \end{equation}

where $( {{\omega _x},{\omega _y},{\omega _z}} )$ is the component of $\boldsymbol {\omega }^b$.

However, the drawback of unit quaternion representation is the lack of intuitive visualization. To address this issue, Euler angles can be derived based on unit quaternions, and yaw angle $\psi$, pitch angle $\theta$ and roll angle $\phi$ are, respectively, expressed as

(2.48)\begin{equation} \left.\begin{gathered} \displaystyle \psi = \arctan \left[ {\frac{{2\left( {{q_1}{q_2} + {q_0}{q_3}} \right)}}{{{q_0}^2 + {q_1}^2 - {q_2}^2 - {q_3}^2}}} \right],\\ \displaystyle \theta =- \arcsin \left( {2\left( {{q_1}{q_3} - {q_0}{q_2}} \right)} \right),\\ \displaystyle \phi = \arctan \left[ {\frac{{2\left( {{q_2}{q_3} + {q_0}{q_1}} \right)}}{{{q_0}^2 - {q_1}^2 - {q_2}^2 + {q_3}^2}}} \right]. \end{gathered}\right\} \end{equation}

3.2.1. Early stage of water entry

The numerical results of the pressure field and cavity shape for several typical instants at the early stage of water entry are given in figure 12. Figure 12 clearly shows that the body gradually crosses the water surface and enters the water, and a high-pressure area around the nose is formed. With the movement of the body, the pressure decreases gradually to approximately 0.95 MPa at $t = 2\ \textrm {{ms}}$. When $t = 8\ \textrm {{ms}}$, the free surface breaks and splashes up, and the splash velocity is approximately $65\ \textrm {m}\ \textrm {s}^{-1}$, which is much higher than the initial velocity of water entry. Meanwhile, the air is entrained, and a cavity is gradually developed. At $t = 20\ \textrm {{ms}}$, a strongly lateral impact occurs on the body, and the initial cavity is formed. During the whole process, the free surface evolves smoothly, and the pressure field is smooth and stable. The cylinder is subjected to severe water impact during the process of water entry. The axial accelerations and normal accelerations at the gravity centre obtained by different grid resolutions are plotted in figure 13 for the convergence test of the present method. The numerical results converge well with decreasing grid sizes.

Figure 12. Snapshots of the cavity shape and the pressure distribution at four time instants ($t = 2\ \textrm {ms}$, 8 ms, 14 ms and 20 ms) during oblique water entry. The solid black line denotes the water–air interface. The inset at the bottom left corner is the cavity shape at the corresponding time, coloured by the velocity value.

Figure 13. Convergence test with four different grid resolutions. (a) Acceleration along the $x_b$-axis, (b) acceleration along the $z_b$-axis.

To test the performance of the improved IB method proposed here, the penalty IB method and Tian's method are also applied to simulate the above oblique water-entry problem with the same grid resolution $R/\Delta x = 30$ and the same initial conditions. Figure 14(a–c) provides the cavity shapes obtained through the three methods at $t = 20\ \textrm {{ms}}$. These cavity shapes seem to be identical in general. As a further comparison of the mentioned IB methods, the measuring point for impact pressure is arranged at the centre of the cylinder nose. The results for impact pressure are plotted in figure 14(d–f), where the same sampling rate for different methods is presented. It is noticeable that the time evolution of pressure by the penalty IB method shows violent oscillations, especially around the pressure peak. The main reason is that the mass differences of the fluid grid nodes around the IB make the velocity field discontinuous and thus result in non-physical oscillations of pressure. For the result from Tian's method, the pressure curve also shows some oscillations, but they are less than those of the penalty IB method. This indicates that it does partly suppress the non-physical oscillations by adding the term of mass weight, but the oscillations are still obvious around the pressure peak. However, the non-physical oscillations are significantly suppressed in the result obtained by the improved IB method. This can be seen more intuitively from figure 14(g–i), which gives the pressure fields with contour lines at $t = 4\ \textrm {{ms}}$ corresponding to the peak time of impact pressure in figure 14(d–f). There are some notable pressure fluctuations for the results from the penalty IB method and Tian's method, especially in the high-pressure region near the cylinder nose. In contrast, the pressure field is quite smooth for the improved IB method.

Figure 14. Comparison of numerical results by three methods: penalty IB method (a,d,g), Tian's method (b,e,h) and improved IB method (c,f,i). (a–c) Cavity shape at $t = 20 \ \textrm {{ms}}$, (d–f) time evolution of pressure at centre of the flat nose, (g–i) pressure field with contour line at $t = 4\ \textrm {{ms}}$, where white line denotes pressure contour line and black line denotes water–air interface.

3.2.2. Later stage of water entry

Since the period of the projectile crossing the free surface is too short, the later stage of water entry is discussed further to test this method. The accuracy of the improved IB method is further examined by comparison with the experimental results. The sensors are arranged on the head of the projectile in figure 9(c). According to rigid body dynamics (Fossen Reference Fossen1994), the acceleration at this location needs to be calculated by

(3.1)\begin{equation} {\boldsymbol{A}}^H = \dot{\boldsymbol{V}}_G^b + \dot{\boldsymbol{\omega }}^b \times {\boldsymbol{X}}_S^b , \end{equation}

where ${\boldsymbol {X}}_S^b$ indicates coordinates of the location of the accelerometers in the LCS, ${\boldsymbol {X}}_S^b = (0.3, 0.0, 0.0)^\textrm {{T}}$ m. In figure 15, the results obtained by the present method are compared with the experimental results in terms of axial acceleration, normal acceleration, pitch angular velocity and cavity shapes. The comparison of axial acceleration is quite good during the whole water-entry process, even at the steep peak. After the peak, the axial acceleration decays exponentially and gradually tends towards flatness. In terms of normal acceleration and pitch angular velocity, they are in good agreement with the experimental results, although small differences occur, which may be primarily caused by the minor uncertain disturbances in the initial conditions of the present experiment. However, these minor disturbances are not taken into account in the numerical simulation. The pitching motion depends upon the water-entry initial conditions, and minor disturbances in the initial condition will induce significant changes in the pitching motion of the projectile. This can also be confirmed by the experimental results that the colour band of the pitch angular velocity in figure 15(c) is relatively wider than those in figures 15(a) and 15(b) (the reader can refer to Appendix A and supplementary material for more details). The phenomenon of tail slamming caused by the differences in pitching motion can be more intuitively seen in 15(d), which compares the cavity evolutions from experimental and numerical results. The deformation of the free surface, shape of the cavity and movement of the body are very close. However, there are some discrepancies in the later stage after tail slamming. These might be caused by the violent impact, which can induce unstable angular velocities of the body and thus affect the cavity evolution. Conversely, the cavity evolution will affect the body movement. The effect of tail slamming can also be confirmed by the inflection point of normal acceleration at approximately $t = 150\ \textrm {{ms}}$ in figure 15(b). Although the peak induced by tail slamming is smaller than that induced by water impact in the early stage, the duration is significantly longer than the former. In addition, tail slamming can also significantly affect the motion trajectory of the projectile. Therefore, it is necessary to pay attention to the effects of tail slamming. The impact loads of tail slamming for different nose types are discussed further in § 4. Overall, the present numerical method for body acceleration, pitch angular velocity and cavity shape shows good agreement with the experimental results throughout the whole water-entry process. In addition, the present numerical method can also be applied to the oblique water-entry problem of a high-speed supercavitating vehicle, see Appendix D.

Figure 15. Comparisons between experimental and numerical results. (a) Acceleration along the $x_b$-axis, (b) acceleration along the $z_b$-axis, (c) angular velocity around the $y_b$-axis and (d) cavity evolution at four moments. Experimental data are shown by shaded regions from 5 repeated experiments (see Appendix A and supplementary material available at https://doi.org/10.1017/jfm.2023.120 for more details).

4.1.1. Early stage of water entry

For the early stage, the projectile strikes and crosses the free surface. The lower side of the projectile touches the free surface first during oblique water entry. Therefore, in addition to the axial drag force, the normal drag force is the focus of hydrodynamic loads. The time evolutions of the axial and normal drag coefficients (defined by (2.1) and (2.2)) of the projectile are plotted in figure 20. For the axial drag force, the peak value of the case of flat nose is approximately three times that of the hemispherical nose, while the pulse width of the former is much smaller, at a third to half of the latter case. For the normal drag force, the case of a flat nose is very small, but the case of a hemispherical nose is very large. This is determined by the pressure distribution on the nose during the water-entry process. Figures 21(a) and 21(b) display the pressure field of the flat nose and hemispherical nose when crossing the free surface, respectively. For the flat nose projectile, when penetrating the free surface, a high-pressure area is formed at the wetted surface of the nose, and the centre of the high-pressure area gradually moves from the edge of the nose to the centre. During this process, the pressure consistently acts in the direction of the normal of the flat nose. Therefore, the axial drag force increases sharply until the centre of the high pressure moves to the centre of the nose, and the peak forms when $t$ is approximately 0.55 ms. For the hemispherical nose projectile, a similar high-pressure area is also formed and gradually moves from the edge to the pole of the hemisphere. However, because the normal direction of the hemispherical nose changes continuously, a normal drag force is also generated along with the axial force. At the moment the nose contacts the water surface, the peak of the normal drag force occurs. In the following, the axial drag force continues to increase, while the normal force gradually decays from the peak. When the centre of the high-pressure area moves to the pole of the hemisphere, the peak of the axial drag force also occurs.

Figure 20. Drag coefficients for different projectiles at early stage of water entry. (a) Drag coefficient on the $x_b$-axis. (b) Drag coefficient on the $z_b$-axis.

Figure 21. Water impact process for three different projectiles during he early stage of water entry: (a) $t = 0.0$ ms, 0.25 ms, 0.4 ms, 0.55 ms for flat nose, (b) $t = 0.0$ ms, 0.1 ms, 0.25 ms, 0.55 ms for hemispherical nose, (c) $t = 0.0$ ms, 0.05 ms, 0.25 ms, 0.55 ms for truncated hemispherical nose. Straight arrows denote main hydrodynamic force exerted by fluid on objects, and curved arrows represent change in pitch moment about the $y_b$-axis.

During oblique water entry, the hydrodynamic pressure is normal to the wetted surface of the projectile nose, and the resultant force may produce a pitch moment about the gravity centre $\boldsymbol {G}$, which directly leads to the hydrodynamic pitch phenomenon. It is noticeable in figure 21 that pitch moments in the opposite direction are produced for different noses. The time evolutions of the pitch moment (defined as the $y_b$-axis component of $\boldsymbol {T}_{0}^b$ and denoted as $T_{yb}$) for different projectiles are given in figure 22. We can clearly see that the directions of the pitch moment of the two cases are opposite. For the flat nose case, the moment is positive and reaches the peak at $t = 0.25\ \textrm {{ms}}$. This is because the centre of the high-pressure area is consistently located at the lower side of the central axis, which passes the gravity centre. Thus, the moment formed by the pressure in the normal direction of the flat nose with respect to the gravity centre is always clockwise, that is, it is positive in this process. In contrast, the high-pressure area acts in the normal direction of the hemisphere, and an anticlockwise moment with respect to the gravity centre is formed, until, at $t = 0.1\ \textrm {{ms}}$, the peak of the moment occurs. The moment peak of the hemispherical nose case is 2.5 times that of the flat nose case. However, the pulse widths of the pitch moment in the two cases are basically the same.

Figure 22. Pitch moment about the $y_b$-axis vs time for three different projectiles at early stage of water entry. Time labels (hollow circles in different colours) correspond to time sequences of snapshots in figure 21, respectively.

4.1.2. Later stage of water entry

As discussed above, the pitch moment in different directions is generated and significantly affects the trajectory. The tail slamming phenomenon of the two cases is shown in figures 23(a) and 23(b). The times of tail slamming for the flat nose case and the hemispherical case are 48 ms and 15 ms, respectively. This is strongly related to the moment formed at the early stage. Because a greater pitch moment is formed for the hemispherical nose case, the projectile begins to slam at the tail more quickly. In addition, the tail of the projectile with a hemispherical nose impacts the lower wall of the cavity, but on the upper wall for the projectile with a flat nose. The axial and normal drag coefficients for the two cases during the whole process of water entry are plotted in figure 24. It can be seen that there is almost no sudden change during tail slamming for the axial drag force. However, tail slamming causes sudden changes in the normal drag force. Especially for the hemispherical nose case, the peak of the normal drag force induced by tail slamming is almost half of the peak when crossing the water surface. Moreover, the pulse width of the tail slamming force is much larger. In comparison, the tail slamming force of the flat nose case is gentle, and a small sudden change occurs. This indicates that tail slamming is also very dangerous for projectiles with hemispherical noses regarding the normal force.

Figure 23. Tail slamming phenomenon at later stage of water entry: (a) flat nose, (b) hemispherical nose, (c) truncated hemispherical nose. Corresponding moments are 41 ms, 12 ms and 54 ms, respectively.

Figure 24. Drag coefficients for different projectiles during water entry process. (a) Drag coefficient on the $x_b$-axis. (b) Drag coefficient on the $z_b$-axis.

As the projectile is a slender body, the severe normal force, which is perpendicular to the revolving axis of the slender body, will inevitably lead to a bending tendency. Therefore, the section bending moment is an important factor in the hydrodynamic loads of water entry, in addition to the axial and normal drag forces. The time evolutions of the bending moment $M_{yb}$ (defined as the $y_b$-axis component of $\boldsymbol {M}^b$) at the cross-section with respect to the gravity centre $\boldsymbol {G}$ are given in figure 25(a). The definition of the bending moment can be found in Appendix C. For the flat nose case, the bending moment reaches the maximum peak at the moment of impacting the water surface and then decreases gradually. When tail slamming occurs, a slight change takes place. However, for the hemispherical nose case, the amplitude of the bending moment is greater than that of the flat nose case. After the initial impact on the water surface, the bending moment begins to increase and reaches the first peak when the centre of the high-pressure area moves to the pole of the hemisphere. After crossing the free surface, the bending moment gradually increases as the projectile rotates and reaches the maximum peak before tail slamming. When tail slamming occurs, the bending moment is suddenly unloaded and then decays when the projectile rebounds back by the cavity wall. This intensive moment can cause the projectile to move violently and threaten the stability of the thin structure of the vehicle, even causing large deformation or break off.

Figure 25. (a) Time evolutions of bending moment at cross-section with gravity centre $\boldsymbol {G}$ for projectiles with different noses. (b) Time evolutions of pitch angular velocity for projectiles with different noses.