2.1. Discrete element model

When a projectile possesses sufficient impact kinetic energy, its intense interaction with marine floating structures during water entry can induce significant deformation or even fracture of the floating structure. Throughout this process, the originally continuous medium of the structure gradually evolves into a discontinuous system composed of numerous discrete particles. To characterise this transition from continuum to discontinuous failure, the DEM is employed to model the floating structure. In the DEM solution process, particle motion is described through explicit integration of Newton’s equations of motion to capture both translational and rotational dynamics of granular flow. In this framework, particle i realises contact, collision, rotation and other movements in space through contact with another particle j, the solid wall and its gravity. The governing equation for translational motion of particles under external forces in DEM can be expressed as

(2.1) \begin{equation} m_i \frac {{\rm d}{v}_i}{{\rm d}t} =\sum _{i \neq j} {f}_{i,j} + {f}_{i,w} + {f}_{i,f}, \end{equation}

where $m_i$ is the mass of particle $i$ , $v_i$ is the velocity of particle, $f_{i,j}$ is the collision force between particles, $f_{i,w}$ is the collision force between particle $i$ and solid wall, and $f_{i,f}$ is the interaction force between particle $i$ and the fluid.

The rotational motion of the particle can be expressed as

(2.2) \begin{equation} I_i \frac {{\rm d}\boldsymbol{\omega }_i}{{\rm d}t} = \sum _{i \neq j} \big ( {M}_{t,ij} + {M}_{r,ij} \big ) + {M}_{i,w} + {M}_{i,f}, \end{equation}

where $I_i$ is the moment of inertia of the ice particle, ${M}_{t,ij}$ and ${M}_{r,ij}$ are the tangential torque and rolling friction torque generated by the contact between particle i and particle j (with $ i \neq j$ ), ${M}_{i,w}$ is the torque generated by the contact between particle i and solid wall, ${M}_{i,f}$ is the torque of fluid acting on particle i, and $\boldsymbol{\omega }_i$ is the angular velocity of the particle.

In this study, the modelling of the floating structure is based on discrete element particle clusters established through classical mechanical principles. Within the DEM framework, a soft-sphere contact approach is employed to identify interactions between particles or between particles and solid boundaries. The method accounts for deformations between contacting bodies, which are represented by particle overlaps and subsequently used to compute interparticle forces. The interaction forces are calculated using a nonlinear slip-spring-damper system. The Hertz–Mindlin contact model (Renzo & Maio Reference Renzo and Maio2005) is adopted in this work to characterise the interactions between particles as well as between particles and walls, as shown in figure 1. The collision force between two particles or between a particle and a solid wall can be expressed as

(2.3) \begin{equation} F_{c} = F_{n}\boldsymbol{n} + F_{t}\boldsymbol{t}, \end{equation}

where $F_{n}$ and $F_{t}$ are the normal and tangential components of the contact force, respectively, and the normal force $F_{n}$ can be defined as

(2.4) \begin{equation} F_{n} = -K_{n}d_{n} - N_{n}v_{n}, \end{equation}

where $K_{n}$ is the normal spring stiffness, $d_{n}$ is the normal overlapping distance between particles or between particles and solid wall, $N_{n}$ is the normal damping, and $v_{n}$ is the normal relative velocity between particles or between particles and solid wall.

Figure 1. Schematic diagram of particle contact model.

The tangential force can be expressed as

(2.5) \begin{equation} \begin{cases} F_{t} = \left |K_{n}d_{n}\right |\!C_{\!{fs}}d_{t}/\!\left |d_{t}\right |, & \quad \left |K_{t}d_{t}\right | \geqslant \left |K_{n}d_{n}\right |\!C_{\!{fs}}, \\ F_{t} = -\!K_{t}d_{t} - N_{t}v_{t}, & \quad \left |K_{t}d_{t}\right | \leqslant \left |K_{n}d_{n}\right |\!C_{\!{fs}}, \end{cases} \end{equation}

where $K_{t}$ is the tangential spring stiffness, $d_{t}$ is the tangential overlapping distance between particles, $v_{t}$ is the tangential relative velocity between particles or between particles and solid wall, and $C_{\!{fs}}$ is the friction coefficient.

The particle bonding approach enables the extension of discrete particle models to simulate complex geometries and structures. By implementing the parallel bond model to represent mechanical connections between particles, cohesive granular assemblies can be constructed. These assemblies exhibit quasi-continuum elastoplastic mechanical responses when subjected to stresses below the failure threshold (Potyondy & Cundall Reference Potyondy and Cundall2004). Once local stresses exceed the stress limit, the model naturally captures the evolution of crack initiation, bond rupture and fragmented particle behaviour, thereby simulating the complete physical process from intact structure to fragmentation. The parallel bond assumes a bonding disk (or cylinder) with radius $R_b$ centred at the contact point between adjacent particles. The bond radius is typically selected as the smaller radius of the two bonded particles. Forces and moments are transmitted simultaneously through the parallel bond, with normal and tangential forces denoted as $F_n$ and $F_t$ , and normal and tangential moments as $M_n$ and $M_t$ , respectively, as shown in figure 2. The tensile stress $\sigma$ (normal bonding strength) and shear stress $\tau$ (tangential bonding strength) acting on the parallel bond can be expressed as

(2.6) \begin{align} \sigma &= -\frac {F_n}{A_b} + \frac {\left |M_{t}\right |}{I_b}R_b, \end{align}

(2.7) \begin{align} \tau &= \frac {\left |F_t\right |}{A_b} + \frac {\left |M_{n}\right |}{J_b}R_b, \end{align}

Figure 2. Parallel bond between particles.

where $A_b$ is the cross-sectional area of the parallel key, $J_b$ is the polar moment of inertia of the parallel key and $I_b$ is the moment of inertia of the parallel key, which can be expressed as

(2.8) \begin{align} A_b &= \pi\! R_b^2,\\[-12pt]\nonumber \end{align}

(2.9) \begin{align} I_b &= \frac {\pi\! R_b^4}{4}, \\[-12pt]\nonumber \end{align}

(2.10) \begin{align} J_b &= \frac {\pi\! R_b^4}{2}. \end{align}

When the tensile stress $\sigma$ or shear stress $\tau$ acting on a parallel bond exceeds its failure strength, the bond will rupture. This bond failure mechanism enables the model to naturally capture crack initiation and damage propagation in the simulated floating ice structure.

In marine environments, floating structures typically exhibit diverse structural compositions and material properties, with their failure behaviour involving complex mixed-mode fracture mechanisms and energy dissipation pathways (Adamson & Dempsey Reference Adamson and Dempsey1998; Schreyer & Munday Reference Schreyer, Sulsky and Munday2006; Paavilainen & Polojärvi Reference Paavilainen, Tuhkuri and Polojärvi2009; Dempsey & Wang Reference Dempsey, Cole and Wang2018; Lilja & Tuhkuri Reference Liu and Ji2021). However, this study focuses on the transient response during the impact water entry of the projectile, with particular emphasis on how structural damage induced by impact affects subsequent hydrodynamic characteristics. Given that the temporal scope of this investigation focuses on the initial impact and fragmentation phases, a fracture criterion based on interparticle bond failure proves sufficient to capture the essential structural damage features while meeting research requirements. This fracture model maintains computational stability while effectively simulating the transition from intact structures to fragmented particulate states, and their influence on cavity evolution and flow field disturbances. Moreover, it possesses the distinct advantage of characterising structural fragmentation and fracture propagation mechanisms, thereby providing an effective numerical approach for investigating the dynamic response of floating structures subjected to projectile impacts.

2.2. Interaction between fluid and particles

The interaction between the fluid and the dispersed solid phase is described using a locally volume-averaged two-phase formulation within the CFD–DEM framework. The total hydrodynamic force ( $\boldsymbol{f}_{\!i,f}$ ) acting on each particle is expressed as the sum of the drag ( $\boldsymbol{f}_{\!d}$ ), pressure-gradient ( $\boldsymbol{f}_{\kern -3pt p}$ ), added-mass ( $\boldsymbol{f}_{\!a}$ ) and buoyancy ( $\boldsymbol{f}_{\!b}$ ) contributions:

(2.11) \begin{align} \boldsymbol{f}_{\!i,f} =\boldsymbol{f}_{\!d}+\boldsymbol{f}_{\kern -3pt p}+\boldsymbol{f}_{\!a}+\boldsymbol{f}_{\!b}, \end{align}

(2.12) \begin{align} \boldsymbol{f}_{\!d} = \frac {1}{2} C_{d} A_{\!p} \rho _{\!f} \big ( \boldsymbol{u}_{\!f} - \boldsymbol{u}_{\!p} \big ) \left | \boldsymbol{u}_{\!f} - \boldsymbol{u}_{\!p} \right |\!, \end{align}

where $A_p$ is the projected particle area, $C_d$ is the drag coefficient, $\boldsymbol{u}_{\!f}$ is the fluid velocity and $\boldsymbol{u}_{\!p}$ is the particle velocity.

The pressure-gradient term accounts for the background flow acceleration:

(2.13) \begin{align} \boldsymbol{f}_{\kern -3pt p} = - V_{\!p} \boldsymbol{\nabla }\!p_{\!f}, \end{align}

(2.14) \begin{align} \boldsymbol{f}_{\!b} = V_{\!p} (\rho _{\!f} - \rho _{p}) \boldsymbol{g}, \end{align}

where $V_{\!p}$ is the particle volume and $p_{\!f}$ is the local fluid pressure.

The added-mass contribution is included to capture the unsteady acceleration of the surrounding fluid:

(2.15) \begin{equation} \boldsymbol{f}_{\!a} = C_{a} V_{\!p} \rho _{\!f} \left ( \frac {{\rm D}\boldsymbol{u}_{\!f}}{{\rm D}t} - \frac {{\rm d}\boldsymbol{u}_{\!p}}{{\rm d}t} \right )\!, \end{equation}

where $C_{a}$ is the virtual mass coefficient.

In computational fluid dynamics, due to the interaction between the gas–liquid interface and particles, the governing equations of the fluid can be written as the volume-averaged continuity and momentum equations:

(2.16) \begin{equation} \frac {\partial (\alpha _{\!f} \rho _{\!f})}{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }(\alpha _{\!f} \rho _{\!f} \boldsymbol{u}_{\!f}) = 0, \end{equation}

(2.17) \begin{equation} \frac {\partial (\alpha _{\!f} \rho _{\!f} \boldsymbol{u}_{\!f})}{\partial t} +\boldsymbol{\nabla }\boldsymbol{\cdot }(\alpha _{\!f} \rho _{\!f} \boldsymbol{u}_{\!f} \otimes \boldsymbol{u}_{\!f}) = - \alpha _{\!f} \boldsymbol{\nabla }\!p_{\!f} + \boldsymbol{\nabla }\boldsymbol{\cdot }(\alpha _{\!f} \boldsymbol{\tau }_{\!f}) + \alpha _{\!f} \rho _{\!f} \boldsymbol{g} + \boldsymbol{S}_{\!f}, \end{equation}

where $\rho _{\!f}$ is the density of fluid, $\boldsymbol{u}_{\!f}$ is the velocity of the fluid, $p$ is the pressure, $\boldsymbol{\tau }_{\!f}$ is the viscous shear stress tensor, $\alpha _{\!f}$ is the volume fraction of the fluid in a single fluid unit (fluid porosity), using the assumption of local average (Tsuji, Kawaguchi & Tanaka Reference Tsuji, Kawaguchi and Tanaka1993), and $S_{\!f}$ is the source term after volume correction (Du et al. Reference Du, Bao, Xu and Wei2006).

To account for the multiphase media, the volume of fluid (VOF) method is employed to track interfaces between immiscible fluids. In this approach, different fluid phases share a single set of momentum equations, while their interfaces are captured by solving transport equations for volume fraction variables assigned to each phase. The fundamental constraint requires the sum of volume fractions in each computational cell to equal 1. Among them, the density and viscosity of the fluid are weighted and averaged according to the phase fraction, that is, $\rho _{\!f} = \alpha _1 \rho _1 + \alpha _2 \rho _2, \mu _{\!f} = \alpha _1 \mu _1 + \alpha _2 \mu _2$ , where the subscripts 1 and 2 represent the liquid phase and the gas phase, respectively. In the VOF framework, the transport equation of the liquid phase is given as

(2.18) \begin{equation} \frac {\partial (\alpha _1)}{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }(\alpha _1 \boldsymbol{u}_{\!f}) + \boldsymbol{\nabla }\boldsymbol{\cdot }\left [ \alpha _1 (1 - \alpha _1) \boldsymbol{u}_r \right ] = 0, \end{equation}

where $\boldsymbol{u}_r$ is the compression velocity, which is an artificially applied relative velocity field designed to compress the VOF interface and thereby reduce numerical diffusion.

The first two terms in (2.18) have the same advection form as the volume-fraction form of the continuity equation, indicating that for an incompressible fluid, mass conservation is equivalent to volume conservation; therefore, a volume-fraction variable is introduced. The presence of $\alpha _{\!f}$ in these two terms arises from the local averaging procedure, while $\alpha _1$ serves as an indicator of wet/dry conditions. Within the computational domain, a control volume with $\alpha _1 = 1$ is fully occupied by liquid, whereas a control volume with $\alpha _1 = 0$ corresponds to a dry air region.

For three-phase flow (gas–liquid–particle), the volume fraction of each phase is shown in figure 3. For three-phase flow, the fluid volume fraction can be expressed as

(2.19) \begin{equation} \alpha _{\!f} = 1 - \sum V_{\textit{pi}}/{V_{t}} = (V_{1} + V_{2})/{V_{t}}, \end{equation}

where $V_{\textit{pi}}$ is the volume of the $i$ particle mapped from the particle domain to the corresponding fluid cell, $V_{1}$ and $V_{2}$ are the volumes of the liquid phase and gas phase, respectively, and $V_{t}$ is the total volume of the cell. While $\alpha _{1}$ and $\alpha _{2}$ can be expressed as the volume fractions of their respective volume ratios ( $V_{\!f} = V_{1} + V_{2}$ ), which can be expressed as

(2.20) \begin{equation} \alpha _{1} = V_{1}/{V_{\!f}}, \end{equation}

(2.21) \begin{equation} \alpha _{2} = V_{2}/{V_{\!f}} = 1 - \alpha _{1}. \end{equation}

Figure 3. Volume distribution of the gas–liquid–particle three-phase of the VOF method.

Therefore, the local average phase fraction $\alpha _{\!f}\alpha _{1} = V_{1}/{V_{t}}$ is calculated based on the total volume $V_{t}$ , which represents the transportation volume as described in (2.20).

The fluid–structure interaction computation is initiated by calculating the fluid porosity in each computational cell based on particle positions and fluid grid information. Subsequently, the fluid–particle interaction forces acting on each particle are determined using particle velocities, fluid velocities and the pressure/stress tensor at the current fluid time step. The Newtonian framework, combined with a drag model, computes the particle–fluid interaction forces within each fluid cell. A sub-cycling DEM time-stepping scheme is then implemented, where particle positions and velocities are updated through explicit integration of Newton’s equations using smaller time increments within each CFD time step. The fluid phase’s mass and momentum conservation equations are solved using the computed porosity and volumetric fluid–particle interaction forces in each cell. The PISO (pressure implicit with splitting of operators) algorithm is employed for solving the fluid equations, with the QUICK (quadratic upstream interpolation for convective kinematics) scheme adopted for discretising convective terms. Gradient and divergence terms are evaluated using a second-order finite difference approach. Numerical stability is governed by the Courant number (Co), which determines the maximum allowable time step. To ensure computational stability, the fluid motion within any grid cell must not exceed a characteristic length per time step, necessitating the use of a fraction of the critical time step for stable temporal integration. The coupling frequency between CFD and DEM domains is determined by the ratio of their respective time steps. Figure 4 presents the complete workflow of this coupling procedure.

Figure 4. Flow chart of fluid–particle interaction.

2.3. Dual-grid method

During grid division, local flow characteristics and resolution requirements must be considered. When projectiles interact with floating structures composed of discrete particles, the relationship between particle size and local fluid grid resolution becomes critical for ensuring both solution accuracy and numerical stability in coupled simulations. To properly capture particle–solid and particle–fluid interaction behaviours, it is generally recommended that particle dimensions should not exceed the scale of local fluid grid cells (Jiang et al. Reference Jiang, Guo, Yu, Hua, Lin, Wassgren and Curtis2021; Cui et al. Reference Cui, Harris, Howarth, Zealey, Brown and Shepherd2022; Mathews et al. Reference Mathews, Khan, Sharma, Kumar and Kumar2022; Zhu et al. Reference Zhu, He, Zhao, Vowinckel and Meiburg2022). However, when dealing with relatively large particles, coarser grids may be employed in surrounding regions to reduce computational costs, which may reduce the resolution of local flow details. In the present study, a balanced approach is adopted to maintain computational efficiency and the analysability of fragmented structures. The particle size is deliberately maintained above a minimum threshold to guarantee that post-impact floating structures develop representative particulate configurations with observable dynamic responses. Simultaneously, the water entry process of the projectile involves complex gas–liquid interface deformation and cavity evolution phenomena, necessitating refined grids for accurate interface dynamics resolution. This inherent conflict between particle-scale requirements and fluid grid refinement creates fundamental challenges in simultaneously satisfying the resolution needs of both physical processes within a single grid system.

To address the scale disparity between particle dimensions and fluid grid resolution, this study introduces a dual-grid methodology. The approach employs a nested coarse–fine grid structure that facilitates bidirectional information transfer between particles and fluids, significantly enhancing simulation efficiency while maintaining physical consistency. The schematic diagram of the dual-grid algorithm is shown in figure 5. Based on classical volume-averaging theory, the method first processes particle distributions on the coarse grid to construct a spatially continuous particle-phase field. This field estimates local porosity, drag source terms and particle-induced momentum transport effects. This process can be regarded as a reverse interpolation operation that projects discrete Lagrangian particle data onto an Eulerian coarse-grid framework for coupling term computation. Following the derivation of averaged particle-phase variables on the coarse scale, the coupling information is subsequently mapped onto the fine-scale fluid resolution grid to capture local flow features such as free-surface deformation and cavity boundaries. The fluid governing equations are solved within this support domain, with initial conditions provided by the coarse-scale averaged field. The projection of Lagrangian variables ensures source-term consistency at the averaging scale, while the support-domain grid can be adaptively refined based on local flow requirements, typically independent of particle distribution characteristics. Upon obtaining the fine-scale fluid solutions, these results are back-projected to update the coarse-grid averaged field, completing an information feedback loop that enables efficient multiscale, multiphysics interaction and dynamic synchronisation between particles and fluids.

Figure 5. A cluster group of continuous grid elements coupled by the dual-grid method. Larger grid elements are used to calculate the interaction between the particle beam and the liquid phase. Refined grid elements are then employed to solve the fluid equations. Upon the interaction calculation, the simulation distributes the effects of volume fraction, momentum, energy source terms and other transport quantities across the component grid elements.

Based on the aforementioned dual-grid coupling framework, this study employs a three-dimensional bilinear interpolation scheme for information transfer between coarse grids ( $\varOmega _{c}$ ) and fine grids ( $\varOmega _{k}$ ), ensuring accurate mapping of physical quantities across different grid scales. To maintain interpolation accuracy and numerical consistency, the fluid domain is discretised using structured hexahedral grids. For each fine-grid node, the interpolation references data from the eight vertex nodes of its encompassing coarse-grid cell, as illustrated in figure 6. Specifically, a given fine-grid node is located within a coarse cell defined by eight vertices with spatial coordinates $(x_{i}, y_{i}, z_{i})$ and $i=1$ –8. These vertices correspond to physical variable values $\varphi _{\textit{coarse}}(x_{i},y_{i},z_{i})$ on the coarse grid. The interpolated value at the fine-grid node is computed through a weighted average, in which the contribution of each coarse grid vertex is determined according to the distance weight $w_{i}$ , which effectively reflects the spatial influence of coarse grid vertices on fine-grid nodes. However, for fine-grid nodes positioned on the interface between two or more adjacent coarse cells, a multi-cell weighted averaging method is adopted to guarantee the continuity of the field across coarse-cell boundaries. The fine-scale quantity is computed as

(2.22) \begin{equation} \varphi _{\textit{fine}}(x,y,z) = \frac {\sum _{i=1}^{N_c} w_i \varphi _{\textit{coarse}}(x_i,y_i,z_i)}{\sum _{i=1}^{N_c} w_i}, \end{equation}

(2.23) \begin{equation} w_i = \frac {1}{| \boldsymbol{r}_{\!f} - \boldsymbol{r}_c |^n}, \end{equation}

Figure 6. Interpolation processing of grid nodes in the dual-grid method.

where $N_c$ denotes the number of adjacent coarse grids, $ \boldsymbol{r}_{\!f}$ and $ \boldsymbol{r}_c$ denote the position vectors of the fine-node and coarse-grid centres, respectively, and $ n$ is a weighting exponent. The value of $ n$ is 2 (weighted by inverse distance square), which is based on its balance between locality and smoothness, and is consistent with the commonly used spatial interpolation methods.

Through the aforementioned weighted averaging procedure, physical variables at fine-grid nodes can be derived from the values at the eight vertex nodes of their encompassing coarse-grid cell. The interpolation process involves two key steps: first, determining the position of the target fine-grid node within the coarse grid (i.e. identifying its host coarse cell); second, extracting the spatial coordinates $(x_{i}, y_{i}, z_{i})$ and corresponding variable values $\varphi _{\textit{coarse}}(x_{i}, y_{i}, z_{i})$ from these eight vertices to compute the interpolation weights $w_{i}$ , ultimately yielding the fine-grid nodal value $\varphi _{\textit{fine}}(x,y,z)$ .

To achieve bidirectional information transfer between fine and coarse grids, the fluid variables obtained from fine-grid computations must be projected back onto the coarse grid to update the averaged field at the coarse scale. This coarsening procedure is implemented through a control-volume weighted averaging approach. The average physical quantity $\overline {\varphi }_{\textit{coarse}}$ on each coarse grid ( $\varOmega _{c}$ ) is calculated. This is determined by collecting the contributions of all fine grids ( $\varOmega _{\!f}$ ),

(2.24) \begin{equation} \overline {\varphi }_{\textit{coarse}} = \frac {\sum _{\varOmega _{c}}\! \varphi _{\textit{fine}}\! V_{\textit{fine}}}{\sum _{\varOmega _{c}}\! V_{\textit{fine}}}, \end{equation}

where $\overline {\varphi }_{\textit{coarse}}$ is the average physical quantity on each coarse grid and $V_{\textit{fine}}$ is the volume of the fine grid. For a given coarse grid node, its value is determined by the central value $\overline {\varphi }_{\textit{coarse}}$ of all adjacent coarse grids sharing the node:

(2.25) \begin{equation} \varphi _{\textit{coarse}}(x_{i}, y_{i}, z_{i}) = \frac {\sum _{N_{c}}\! w_{i} \overline {\varphi }_{\textit{coarse}}}{\sum _{N_{c}}\! w_{i}}. \end{equation}

Note that (2.25) reconstructs coarse-grid nodal values from cell-centred averages for use in the next interpolation step (2.23). While the inverse-distance weighting does not guarantee exact reversibility between nodal and cell-centred values, global conservation is maintained via the volume-weighted averaging in (2.24) after each fine-grid solution, because the sum of fine-grid fluxes over a coarse cell equals the coarse-cell flux. To minimise artificial diffusion of sharp gradients, the fine grid is concentrated in regions of high flow activity (impact zone, cavity interface) and a gradient-limited interpolation is employed when mapping from coarse to fine grids. The method has been validated on benchmark problems with steep gradients and shows negligible loss of local accuracy for the present impact water entry simulations.

The dual-grid system, through the aforementioned interpolation and weighted-averaging procedures, establishes bidirectional variable transfer between coarse- and fine-scale grids, ensuring effective coupling between particulate and fluid phases while demonstrating remarkable precision preservation in high-resolution computational domains. This methodology constructs averaged fields on coarse grids to reduce computational expense, while simultaneously resolving local fluid dynamics through fine-scale meshing, thereby maintaining global numerical stability without compromising local accuracy. By synergistically combining volume-averaging theory with multipoint interpolation techniques, the approach provides a robust and efficient solution strategy for multiscale multiphase flow problems, particularly well suited for numerical simulations involving complex boundary deformation, particle fragmentation and strong fluid–structure interactions.

3.1. Physical model

Under the CFD–DEM dual-grid algorithm framework, a physical model for the impact water entry of a projectile is established, as shown in figure 7. The dimensions of the projectile are scaled based on the currently employed marine detectors (Li et al. Reference Li, Sun, Yao, Wang and Li2020; Jiang et al. Reference Jiang, Guo, Yu, Hua, Lin, Wassgren and Curtis2021; Hu, Wei & Wang Reference Hu, Wei and Wang2023) and modelled as a rigid body, featuring a conical-nosed cylindrical geometry with a diameter of $D_{0} = 0.016\,\text{m}$ , a body length of $L_{\!p} = 3.375D_{0}$ , a nose cone angle of $\theta = 120^\circ$ (the conical head is processed by arc and the radius of arc is 0.0008 m) and a cross-sectional area of $S_{p} = \pi (D_{0})^{2}/4$ . The projectile is constructed from aluminium with a density of $\rho _{p} = 2700\,\text{kg}\,\text{m}^{-3}$ , an elastic modulus of $E_{\!p} = 27\,\text{GPa}$ and a strength of $150\,\text{MPa}$ . The computational domain measured $L_{\!f} \times W_{\!f} \times H_{\!f} = 20D_{0} \times 20D_{0} \times 42D_{0}$ , with a water depth of $H_{w} = 32D_{0}$ . Pressure outlet boundary conditions are imposed on the domain boundaries. The density and viscosity of water are $\rho _{1} = 1000\,\text{kg}\,\text{m}^{-3}$ and $\mu _{1} = 0.00179\,\text{Pa}\boldsymbol{\cdot }\text{s}$ , while the density and viscosity of air are $\rho _{2} = 1.294\,\text{kg}\,\text{m}^{-3}$ and $\mu _{2} = 0.000017\,\text{Pa}\boldsymbol{\cdot }\text{s}$ , respectively.

Figure 7. Physical model of impact water entry. The entire domain is divided into three regions: the air domain, the water domain and the particle domain. Air and water form an interface (liquid surface), with the flat ice composed of discrete particles arranged in layers following a hexagonal close-packed structure. The size and physical properties of each particle are uniform.

In marine environments, particularly in polar and high-latitude regions, the free surface is often populated with various floating structures, including floating ice, marine debris and natural drift objects. These floating structures, commonly encountered in polar exploration, marine deployment and underwater operations, constitute a non-negligible pre-existing medium during the water entry process of the projectile. Addressing this practical concern, this study incorporates floating ice as a representative floating structure into the modelling and analysis of impact water entry. In terms of modelling, recent studies have employed DEM to simulate the mechanical response of ice, typically by constructing particle systems with bonding models to reconstruct intact ice structures for numerical analysis of their physical parameters and failure behaviour (Long, Ji & Wang Reference Long, Ji and Wang2019; Pradana & Qian Reference Pradana and Qian2020; Zhang et al. Reference Zhang, Tao, Wang, Ye and Guo2021). Although these investigations have primarily focused on solid mechanics without addressing strong fluid–structure interactions involving fluid, they have established critical modelling foundations for this study, particularly in particle parameter selection and contact model formulation. Within the optimised CFD–DEM coupling framework, floating ice is represented by discrete particle assemblies, where mechanical continuity is maintained through parallel bond models to effectively characterise ice deformation and fracture under impact loading. To balance numerical stability with physical representativeness, the ice dimensions and material parameters are derived from actual polar ice properties, with appropriate geometric and physical scaling based on similarity principles (Yang et al. Reference Yang, Sun, Zhang, Yang and Lubbad2024a ,Reference Yang, Sun, Zhang, Yang and Lubbad b ).

The floating ice is modelled with fundamental dimensions of $L_{l} \times L_{w} \times L_{t} = 12.5D_{0} \times 12.5D_{0} \times 1D_{0}$ , a density of $\rho _{\textit{ice}} = 900\,\text{kg}\,\text{m}^{-3}$ , a bending strength of $\sigma _{f,\textit{ice}} = 0.81\,\text{MPa}$ , a compressive strength of $\sigma _{c,\textit{ice}} = 2.56\,\text{MPa}$ and an elastic modulus of $E_{\textit{ice}} = 0.7\,\text{GPa}$ . From a discrete element perspective, the ice is discretised into elemental particles whose material properties govern the bulk behaviour of the floating ice. The physical parameters of the ice have an adaptive relationship with the particle parameters. The particle-scale parameters set in the numerical simulation are based on the particle-bond mapping framework established by Long et al. (Reference Long, Ji and Wang2019) with a particle density of $\rho _{\textit{ip}} = 900\ \text{kg}\,\text{m}^{-3}$ , normal bonding strength of $\sigma = 1.2\ \mathrm{MPa}$ , tangential bonding strength of $\tau = 0.8\ \mathrm{MPa}$ , friction coefficient of $C_{\!{fs}} = 0.2$ and elastic modulus of $E_{\textit{ip}} = 0.5\ \mathrm{GPa}$ . A detailed validation of these particle parameter selections will be presented in subsequent sections.

Initially, the ice remains stationary on the liquid level. The velocity of the projectile is denoted as $u_{\!p}$ , with an initial water entry velocity of $u_{\!p0} = 20\,\text{m s}^{-1}$ , and the $\textit{Fr} = u_{\!p0}/\sqrt {\textit{gD}_{0}} = 50$ . Characteristic parameters are dimensionless in the numerical simulations, including the hydrodynamic force $F_{n}^{*} = F_{h}/0.5\rho _{w}c_{\textit{ice}}^{2}S_{p}$ , collision force $F_{c}^{*} = F_{c}/0.5\rho _{w}c_{\textit{ice}}^{2}S_{p}$ , water entry speed $u_{\!p}^{*} = u_{\!p}/c_{\textit{ice}}$ , movement time $t^{*} = t/(D_{0}/c_{\textit{ice}})$ and displacement $L^{*}_{\textit{pi}} = L_{\textit{pi}} / D_{0}\ (i = x, y, z)$ . Since the impact of the projectile on the floating ice involves stress wave propagation, the characteristic velocity ( $c_{\textit{ice}}$ ) used in the dimensionless treatment is defined as the stress wave propagation velocity in ice.

3.2. Numerical method verification

The fluid domain is discretised using structured hexahedral grids, as shown in figure 8, to balance computational efficiency with the data mapping regularity between coarse and fine grids in the dual-grid framework. This grid topology offers superior alignment consistency and controllable resolution, thereby enhancing both coupling efficiency and convergence stability in the numerical simulations.

Figure 8. Grid division of the fluid domain.

To maintain numerical stability during time advancement, the Courant number ( $Co = U\!\Delta t/\Delta x$ , where $U\!\Delta t$ represents the fluid flow distance and $\Delta x$ is the grid length) must be maintained below 1. The grid resolution is characterised by $D_{0}/\Delta x$ , where $D_{0}$ serves as the characteristic length scale.

To evaluate the rationality of the fluid domain meshing strategy and ensure numerical convergence, a grid-independence study is conducted using a canonical free water-entry case. This case is selected because it employs exactly the same fluid solver settings, boundary conditions, free-surface treatment and grid generation strategy as the ice-breaking water entry simulations, while avoiding DEM coupling. The absence of discrete elements substantially reduces the computational cost, allowing a more systematic assessment of mesh refinement effects under identical fluid conditions and ensuring that the selected mesh resolution is both physically accurate and computationally efficient. During validation, the background grid size is systematically adjusted in the vicinity of the projectile and along its water entry trajectory, with primary focus on analysing the convergence behaviour of the body’s velocity response under different grid resolutions. Figure 9 illustrates the velocity variations during the impact water entry for various grid resolutions. The numerical results demonstrate that as grid refinement progresses, the velocity components in all directions exhibit convergent behaviour, indicating satisfactory grid independence. Balancing numerical accuracy with computational cost, a grid resolution of $D_{0}/\Delta x=16$ is selected for this study, corresponding to a minimum grid size of 0.001 m near the free surface. This configuration adequately satisfies the resolution requirements for both free surface tracking and cavity interface capture.

Figure 9. Velocity of the projectile during water entry under different grid resolutions. The free water-entry condition maintains consistency with the impact water entry, $\textit{Fr}=50$ . Velocity components along the $x$ , $y$ and  $z$ axes of the Cartesian coordinate system are obtained to characterise the directional velocity distributions.

To validate the numerical accuracy and physical reliability of the CFD–DEM coupling framework with the dual-grid method, a laboratory experiment of impact water entry is designed, as shown in figure 10. The experiment is performed in a $1.5\,\text{m} \times 0.8\,\text{m} \times 0.9\,\text{m}$ glass water tank maintained at low ambient temperature to effectively retard ice phase transition processes.

Figure 10. Experimental device for impact water entry: (a) schematic diagram of the experiment; (b) experimental equipment.

The projectile is launched using a pneumatic ejection system, where the gas filling volume is controlled through a preset air pressure and instantaneously released via an electromagnetic valve. The gas content is adjusted to achieve an initial impact velocity of $u_{\!p0} = 20\,\text{m}\,\rm s^{-1}$ , with the projectile’s geometric dimensions and material parameters maintained consistent with those specified in § 3. The water entry process is captured using a FASTCAM SA-5 high-speed imaging system operating at 1500 fps to ensure sufficient temporal resolution for the dynamic process. A 50 W directional LED light source, positioned behind the water tank and coupled with a flexible diffuser screen, provides uniform backlight illumination to minimise measurement errors caused by interfacial refraction. The spatial relationships among the camera, light source and water surface are calibrated to maintain stable focus in the imaging area and ensure clear visualisation of the water entry cavity. Synchronisation between the launching device and imaging system is achieved through a central control computer.

The ice specimens used in the experiment are prepared by the artificial freezing method to ensure their geometric consistency and physical controllability. The ice specimens used for impact water entry experiments have standardised dimensions of $L_l \times L_w \times L_t = 12.5D_0 \times 12.5D_0 \times 1D_0$ , with a density of $\rho _{\textit{ice}} = 900\,\text{kg}\,\text{m}^{-3}$ , a bending strength of $\sigma _{f,\textit{ice}} = 0.81\,\text{MPa}$ , a compression strength of $\sigma _{c,\textit{ice}} = 2.56\,\text{MPa}$ and an elastic modulus of $E_{\textit{ice}} = 0.7\,\text{GPa}$ . Before impact water entry, each ice specimen is positioned to float freely on the water surface with its geometric centre precisely aligned with the rotational axis of the projectile, ensuring ideal axisymmetric normal impact conditions. The initial water entry velocity is set to $u_p = 20\,\text{m}\,\text{s}^{-1}$ .

Numerical simulations maintained identical geometric parameters and initial motion conditions for the projectile as specified in § 3, with particular attention given to the construction of the floating ice model at the free surface. The ice’s dimensions and material properties correspond to the experimental conditions, where the thickness $L_t = 1D_0$ , the number of particle layers is 8 for ice thickness (the total number of particles is 63 744) and the corresponding single particle size is $D_{\textit{ip}} = 0.0024\,\text{m}$ , based on the particle size calculation method in the supplementary material available at https://doi.org/10.1017/jfm.2026.11203. In the optimal grid strategy adopted in the fluid domain, the minimum grid size of the background area is $0.001\,\text{m}$ , which is smaller than the size of particles, so the dual-grid method is adopted to resolve the matching problem between the fluid grid and the particle size.

A systematic comparative analysis is conducted between numerical simulations and experimental results of impact water entry. Figure 11 presents the dynamic responses during the impact process, including cavity evolution, velocity and displacement. The simulation results obtained using the dual-grid optimised CFD–DEM algorithm demonstrate agreement with experimental data in terms of trajectory, velocity drop characteristics and overall water entry behaviour. The numerical simulations exhibit certain discrepancies with experimental results in terms of the quantity and morphology of floating ice fragments generated during impact-induced breakage. These differences primarily stem from the inherent randomness of ice fracture processes and the simplified treatment of particle geometry in discrete element modelling. Nevertheless, the numerical approach effectively captures the overall motion characteristics and spatial distribution patterns of the fragments, demonstrating that the proposed coupling method possesses strong physical representational capabilities for modelling dynamic fracture processes. Comparative analysis with conventional grid simulations reveals that the absence of dual-grid refinement leads to symmetric cavity structures after impact, failing to reproduce the cavity fragmentation and collapse observed experimentally due to insufficient particle–fluid coupling. Both experimental observations and dual-grid simulations demonstrate that cavity development is significantly influenced by ice fragments, resulting in irregular boundaries and discontinuous cavity structures. This highlights the critical role of fractured particles in the coupled body–ice–fluid system. Conventional CFD–DEM approaches with unmatched grid scales fail to adequately resolve two-way fluid–particle interactions, which significantly undermines the accuracy of cavity evolution and the flow field.

Figure 11. Comparison between numerical simulations and experimental results ( $\textit{Fr}=50$ , size of ice $L_{l} \times L_{v} \times L_{\!f} = 12.5D_{0} \times 12.5D_{0} \times 1D_{0}$ ). (a) Experimental results of the impact water entry process. (b) Numerical simulation results of the impact water entry process with the dual-grid method. (c) Numerical simulation results of the impact water entry process without the dual-grid method. (d) Comparison of dimensionless vertical velocity components during the impact water entry process. (e) Comparison of dimensionless vertical displacement components during the impact water entry process. Panels (a)–(c) show front views ( $x$ – $y$ plane observations). In panels (b) and (c), upper parts show global impact processes while lower parts display cavity evolution by hiding particle phases.

To further evaluate the differences in ice structure damage evolution and broken ice particle dynamics under different algorithms, this study conducted a comparative analysis of the failure patterns and particle motion characteristics of floating ice during impact water entry, as shown in figure 12. In the experiment, a flexible capture net is installed at the tank bottom to capture fragmented ice pieces, which are subsequently retrieved and reconstructed immediately after each test to examine final fracture morphology. Figure 12(a) presents the observed ice failure modes. When employing the dual-grid method, despite minor discrepancies in fracture locations compared with experimental results, considering inevitable experimental disturbances (such as temperature gradients, water surface fluctuation, etc.), the simulations successfully reproduced the primary crack propagation paths and their topological features, with fracture counts matching experimental observations. In contrast, simulations without dual-grid refinement exhibited more branched cracks and premature structural failure, along with significantly enhanced particle dispersion phenomena at both the ice–fluid interface and the lower surface of the penetration zone. These observations are corroborated by figures 11(c) and 13, demonstrating that conventional grid resolution fails to resolve local fluid–particle coupling in critical failure regions, leading to amplified contact forces and consequent unphysical particle detachment and structural damage propagation. In addition, the expansion of the damage range in the central area is also related to the excessive grid refinement. The fundamental reason is the numerical error caused by the mismatch between fluid grid and particle size. In the traditional coupling method, the fine grid makes the fluid solver too sensitive to the micro-scale velocity/pressure gradient at the particle boundary. This sensitivity overestimates the momentum transfer rate of fluid to particles and thus calculate the excessive fluid–particle interaction force. These non-physical strong forces drive the excessive particle splashing and abnormal damage propagation observed in the simulation. The introduced dual-grid coupling framework effectively enhances the physical resolution of particle–fluid interactions, validating both the necessity and efficacy of this approach for simulating water impact problems involving marine floating structures.

Figure 12. Crack evolution of ice during impact water entry ( $\textit{Fr}=50$ , size of ice $L_{l} \times L_{v} \times L_{\!f} = 12.5D_{0} \times 12.5D_{0} \times 1D_{0}$ ). (a) Crack evolution of ice with the dual-grid method. (b) Crack evolution of ice without the dual-grid method. The crack evolution process of the ice is observed from the top view ( $x$ – $z$ plane). In panel (a), the crack position is indicated by a red line at 0.015 s, when the crack in the ice is fully developed. At the centre of the circumferential crack, the ice has completely fractured and is submerged in water as small-scale crushed particles. The capture net is unable to capture these small particles, leaving only a large-scale circumferential crack. In panel (b), particles surrounding the ice are small-scale crushed ice particles resulting from fluid forces. From the $x$ – $z$ plane perspective, this part includes particles that are splashing obliquely upward into the air and those falling obliquely downward into the water.

Figure 13. Development of crushed ice particles under the condition of fully developed ice cracks (t = 0.0015 s). (a) Development of ice cracks with the dual-grid method. (b) Development of ice cracks without the dual-grid method. (c) Accumulated position of crushed ice in the plane mapping with the dual-grid method. (d) Accumulated position of crushed ice in the plane mapping without the dual-grid method. Panels (a)–(d) are bottom views, which are the perspective of the $x$ – $z$ plane from bottom to top. Panels (c)–(d) show the positions of the crushed ice particles in the water mapped onto the $x$ – $z$ plane. The overlapping of particles represents different depth positions in the water; the more compact the particle accumulation, the greater the density of crushed ice particles at that depth.