Three-dimensional numerical simulation of cavity dynamics of a stone with different spinning velocities

Highlights

• The evolution of cavity morphology is well captured and investigated by three-dimensional numerical simulation.

• The rotation around the y-axis are discovered. And the factors caused the rotation are analyzed.

• The pressure distribution for the stone with spinning velocity is no longer regular and symmetrical, and the factors are analyzed.

• The trajectory of the stone with higher spinning velocity shows a more apparent bend in the o-xy plane.

Abstract

To investigate the complex cavity evolution, hydrodynamic action and motion response of a stone impacting on the water, a calculation model is established with six degrees of freedom. The Large Eddy Simulation (LES) model is used to solve the turbulence based on Computational Fluid Mechanics (CFD). The accuracy of the numerical model is verified by comparing the numerical results with experimental results. The stone is released above the free surface with different initial spinning velocities under a certain attack angle. The cavity dynamics, including the dome, open cavity and splash, are investigated during the impact process. Asymmetric impact cavity and the rotation around the y-axis are discovered. And the factors caused the rotation are analyzed. In particularly, there is a special cavity flow phenomenon found by the simulation that obvious tortuosity appears between first cavity and second cavity for a certain spinning velocity. Moreover, the hydrodynamic action, trajectory and reduce rate of velocity are studied for different spinning velocities. Compared to the non-spinning stone, the pressure distribution for the stone with spinning velocity is no longer regular and symmetrical. The maximum value of the pressure at the initial moment of the impact increases with increasing spinning velocity. The results demonstrate that the deviation in z-axis also increases with increasing spinning velocity. Furthermore, the trajectory of the stone with higher spinning velocity shows a more apparent bend in the o-xy plane. In the primary stage of impact, the spinning velocity of the stone decreases at a faster rate for the stone has a higher initial spinning velocity. Then, a transient increase of the reduce rate appears for the stone with violent tumble.

Introduction

The impact of objects on water has been investigated over a hundred years (Bodily et al., 2014; Kiara et al., 2017; R. Zhao and O. Faltinsen, 1993; Truscott and Techet, 2009; Wu et al., 2004; Xia et al., 2019). Most of these works have mainly focused on hydrodynamics force (Glasheen and Mcmahon, 1996; Ma et al., 2016), cavity evolution (Mansoor et al., 2014; Wei and Hu, 2015), and trajectory Xia et al., 2019 [11] of the objects passing through the water surface. In recent decades, skipping stone, where the object rebounds from the water surface, has attracted the attention of scholars from all over the world. Besides being a recreational and competitive sport (Most skips of a skimming stone, Guinness World Records 2013, UNITED STATES EMPORIUM PA, n.d.), skipping stone has also been widely used in military engineering, such as the surface landing of seaplanes (Khabakhpasheva and Korobkin, 2013) and the skipping flight of amphibious Unmanned Aerial Vehicle (UAV) (Paul, 2001). The most famous military application during World War II is the ‘Dambuster’ (http://www.dambusters.org.uk/), which succeeded in destroying the target by using the principle of skipping stone to avoid obstacles. Later, Johnson (1998; 1975) gave a brief history of the idea and technology of the projectile bouncing before it reaches its target. The phenomenon of skipping stone also can be found in the nature. Glasheen and McMahon (1996) used a high-speed camera to record and analyze the whole running process of basilisk lizard on the water surface. Therefore, it is necessary to have a deeper investigation of skipping stone.

In the available public literature, skipping stone has been established based on experiments (Bocquet and Clanet, 2006; Clanet et al., 2004; Rosellini et al., 2005). Bocquet and Clanet (2006) found that the skipping stone must rotate at a minimum speed to maintain stability, which is called gyroscopic effect. Clanet et al. (2004) found the optimum water entry angle called ‘magic angle’ was 20°. Rosellini et al. (2005) then extended to multiple bounces under further experimental work, uncovering four different cavity evolution modes due to different motion responses of the stone.

A computational model was proposed under the condition of high spin, and the physical principle of skipping and the source of dissipation during rebounding were analyzed (Babbs, 2019; Bocquet, 2003). Nagahiro and Hayakawa (2005) established the three-dimensional numerical calculation model of the skipping based on SPH method without analysis of the motion process, and deduced the motion differential equation of the stone. Truscott et al. (2014) described the physics of the rebounding behavior of stones and spheres, and emphasized how the elasticity and flexibility of objects promote the skipping. Smith and Liu (2017) examined the conditions under which a body is able to either emerge from or sink deeper into water, with particular respect to the incident angle, rotation, scaled gravity and buoyancy of the body.

The published research are mainly obtained through experiments (Clanet et al., 2004; Bocquet and Clanet, 2006; Hewitt et al., 2011; Truscott et al., 2012, 2014) and theoretical models (Bocquet, 2003; Nagahiro and Hayakawa, 2005; Yabe et al., 2005; Hicks and Smith, 2011; Khabakhpasheva and Korobkin, 2013; Liu, 2017). In experimental study, some data are difficult to capture, such as the cavity morphology at different viewing angles and pressure distribution on the disk. Simulating this phenomenon numerically is a necessary way to improve these results. Yan and Monaghan (2017) simplified the motion in two dimensions and used SPH method to investigate the change of velocity and position. As a result, the three-dimensional cavity dynamics and the effect of spinning velocity were ignored in their calculation. Therefore, a numerical calculation model in three dimensions is established based on Volume of Fluid (VOF) multiphase flow model. And the effects of spinning velocity on the cavity dynamics, hydrodynamic action and motion response of a stone impacting on the water surface, or more exactly the impact of a circular plate are studied in the present paper.

This paper is organized as follows: Section 2 shows the numerical approach and computational details for the simulated model. Section 3 represents the validation of numerical sensitivity and validation test. The effects of spinning velocity on the cavity evolution, hydrodynamic action and motion response are reported in Section 4.1 to Section 4.3, respectively. Finally, Section 5 draws the main conclusions.

Compared with the Reynolds-averaged Navier-Stokes (RANS), LES was found to give a better overall quantitative agreement with the experimental data (CHENG et al., 2003; Xie and Castro, 2006) for the turbulent flow at a relatively high Reynolds number. For the Reynolds number defined by Re = U₀R/υ is of the order of 10⁵ for the velocity in the present paper, where υ is the viscosity of water, LES is set to simulate the motion of the water.

The Sub-grid Scale (SGS) model proposed by Smagorinsky (1963) is used in the present study, which takes the sub-grid motions into account. For the VOF model for air/water two-phase flows, the multiphase fluid components are assumed to share the same velocity and pressure. The main governing equations consist of the mass and momentum conservation equations,

$$\frac{\partial \rho }{\partial t}+\frac{\partial \left(\rho u_j\right)}{\partial x_j}=0$$

$$\frac{\partial \left(\rho u_i\right)}{\partial t}+\frac{\partial \left(\rho u_i{u}_j\right)}{\partial x_j}=-\frac{\partial p}{\partial x_i}+\frac{\partial }{\partial x_j}\left(\mu \frac{\partial u_i}{\partial x_j}\right)$$

where, u_i, i ∈ [1, 2, 3] denotes the velocity component in three directions (x, y and z) in Cartesian coordinate system, respectively. p is the mixture pressure. The laminar viscosity, μ, and the mixture density, ρ, are defined as

$$\mu ={\alpha }_a{\mu }_a+(1-{\alpha }_a){\mu }_w$$

$$\rho ={\alpha }_a{\rho }_a+(1-{\alpha }_a){\rho }_w$$

where, subscripts w and a refer to the water and the air, respectively. α_a is the volume fraction of air in each cell, and (1-α_a) corresponds to the water volume fraction in each cell.

Applying a Favre-filtering operation to Eqs. (1) and (2) gives the LES equations:

$$\frac{\partial \rho }{\partial t}+\frac{\partial \left(\rho \overline{u_j}\right)}{\partial x_j}=0$$

$$\frac{\partial \left(\rho \overline{u_i}\right)}{\partial t}+\frac{\partial \left(\rho \overline{u_i}\overline{u_j}\right)}{\partial x_j}=-\frac{\partial \overline{p}}{\partial x_i}+\frac{\partial }{\partial x_j}\left(\mu \frac{\partial \overline{u_i}}{\partial x_j}\right)-\frac{\partial {\tau }_{ij}}{\partial x_j}$$

where, $\overline{u_i}$,i ∈ [1, 2, 3] denotes the filtered velocity component in three directions (x, y and z) in Cartesian coordinate system, respectively. $\overline{p}$ is the filtered pressure, and τ_ij is the SGS stresses and given by the following equation,

$${\tau }_{ij}=\rho \left(\overline{u_i{u}_j}-{\overline{u}}_i{\overline{u}}_j\right)$$

By Boussinesq approximation, introducing the sub-grid eddy viscosity ${\upsilon }_t$ yields the following equations for the SGS stress,

$${\tau }_{ij}-\frac{1}{3}{\delta }_{ij}{\tau }_{kk}=-2{\upsilon }_t{\overline{S}}_{ij}$$

$${\overline{S}}_{ij}=\frac{1}{2}\left(\frac{\partial {\overline{u}}_i}{\partial x_j}+\frac{\partial {\overline{u}}_j}{\partial x_i}\right)$$

where δ_ij is the Kronecker delta and ${\overline{S}}_{ij}$ denotes the strain rate tensor in the resolved field. ${\upsilon }_t$ is a function of ${\overline{S}}_{ij}$ and the sub-grid length l, and can be calculated by the following equation,

$${\upsilon }_t=l^2\left|{\overline{S}}_{ij}\right|$$

$$l=C_s\overline{\Delta }$$

$$\left|{\overline{S}}_{ij}\right|=\sqrt{2{\overline{S}}_{ij}{\overline{S}}_{ij}}$$

where $\overline{\Delta }$ denotes the grid filter width, and $C_s=0.2$is the Smagorinsky constant.

Schematic of the skipping problem considered in the current study is displayed in Fig. 1. One end of the stone with mass M is just above the air-water interface. Referring to the experiment of Rosellini et al. (2005), the stone with a radius of R = 25mm and a thickness of h = 2.75mm (where h ≪ R) has an attack angle of α. The stone has a translation velocity U₀ and a spinning velocity Ω Ω n, where n is the unit vector normal to the stone surface. The direction of motion is define by the impact angle β. The simulations are conducted with the stone (s) to water (w) density ratio: ρ_s / ρ_w ≈ 2.7, and the mass M of the stone is 14.6g. The investigation of the motion is based on the two following coordinate systems:

Inertial coordinate system: The coordinate system o-xyz is established at the first contact point o between the stone and water surface, where the x-axis is perpendicular to the water surface and o-yz is the fixed horizontal plane.

Body-fitted coordinate system: The coordinate system o₁-x₁y₁z₁ is established at the centroid o₁ of the stone, where the x₁-axis coincides with the central axis of the stone. And the coordinate translates with the stone but does not rotate with the stone.

To solve the problem of skipping stones, a numerical model is established by the STAR-CCM+ software. PISO (Pressure Implicit with Splitting of Operators) algorithm is used to solve the coupling between pressure and velocity in the computational process, and the second-order discrete scheme is used for transient simulations. The overset mesh is created between fluid region and the rigid stone by using the overset mesh boundary, which achieves the filed data exchange between the fluid and solid regions shown in Fig. 2(a). The dynamics of the fluid is solved by the governing equations. (1) ~ (12) listed in the present paper based on the VOF model and the LES model. Moreover, the Dynamic Fluid Body Interaction (DFBI) module is used to simulate the motion of the stone in response to pressure and shear forces the fluid exerts, and to additional forces we define (such as the gravity). When the DFBI is added to the stone, the corresponding six DOF solver is created. Then, fluid forces, moments, and gravitational forces on the stone is computed based on the six DOF solver. And the forces and moments acting on the stone are used to compute the translational motion of the center of mass of the body and the angular motion of the orientation of the body by solving the governing equations of rigid body. Furthermore, a free-slip wall condition is applied to the four side walls and bottom of the domain, and gradient of volume fraction is considered to be zero. A stagnation inlet condition is applied to the upper side of the domain, implying the velocity gradient is set to be zero. The boundary condition for the stone surface is taken to be no-slip boundary where pressure gradient along the direction perpendicular to the wall is zero.

The size of computational domain and grid structures is illustrated in Fig. 2(b). The length and width of domain are 39R and 12R, and the depth of water is 10R. The numerical model is solved based on an unstructured body-fitted grid. To satisfy with the requirements of LES model for solving and capturing steep gradients accurately, the trimmer mesher model for high-quality grids, the surface remesher model, and the prism layer mesher model are employed. A refined mesh, is applied close to the surface of the stone in both axial and radial directions, the region of movement and water surface. Moreover, to ensure the non-dimensional distance y⁺ < 1 for the LES model(Georgiadis et al., 2010), the boundary layer is added to the stone surface marked with red rectangle as shown in Fig. 2(b) by the prism layer mesher model.