The entry of projectiles into water has been of interest to many scientists and engineers, being crucial to a wide range of engineering applications. The water entry problem is a nonlinear and unsteady phenomenon involving complicated multi-phase flow problems and fluid–solid interaction. Many scientists have been studying water entry problems in various conditions through experimental methods and numerical methods. In this paper, three-dimensional numerical simulations of the water entry problem are carried out. The multiphase flow weakly compressible smoothed particle hydrodynamics model is adopted and three-phase interaction is analyzed using pairwise force smoothed particle hydrodynamics. Dynamic boundary condition and rigid body coupling are introduced for interaction between fluid and solid. Spheres with different wetting characteristics entering water at small Reynolds numbers are investigated. Our results show good agreement with the theoretical models from previous studies into the splashing behavior of spheres. The physics of the different splashing behaviors is discussed in detail.

The entry of a solid projectile into water is a ubiquitous phenomenon that is encountered in everyday life (Ding et al., 2015; Truscott et al., 2014). However, the evolution and formation of the resulting splash, cavity, and jet show complex variations depending on the impact conditions, the material properties, and the body geometries (Aristoff and Bush, 2009; Duez et al., 2007; Gilbarg and Anderson, 1948; Kim and Park, 2019; Marston et al., 2016; May and Woodhull, 1950; Techet and Truscott, 2011; Truscott and Techet, 2009; Truscott et al., 2014; Wei and Hu, 2014; and 2015). Such variations influence the projectile motion into the water, and this is a motivation for investigating the physics behind these phenomena, which have potential applications in a broad range of areas (Kim and Park, 2019). At large scales, well-known applications include ship slamming, sea-landing planes, the Wallis bomb, skipping stones, ballistics, and rowing, while applications at small scales include the biomimetics of bio-locomotion and film coatings (Truscott et al., 2014).

The splash and the cavity formed from the entry of a projectile into water were first investigated in 1897 by Worthington and Cole (1897), who used single-spark photography to capture the air cavity formed behind a solid sphere vertically entering water. Since the advent of modern high-speed camera techniques, which allow the instantaneous capture of high-resolution images, many experiments have been conducted investigating these phenomena (Truscott et al., 2014). It was found that a horizontal jet traveling about 30 times faster than the entry speed is ejected right after the projectile impact (Thoroddsen et al., 2004). At the same time, a splash crown is also generated above the free surface after the impact. Aristoff and Bush (2009) quantitatively modeled the growth and formation of this crown and successfully predicted experimental results. They also investigated the evolution and formation of the cavity behind the object below the free surface. They categorized the air-entraining cavities into four groups according to the generation of the pinching-off processes at low Bond numbers (Bo = ΔρgL2/σ, with ρ, g, L, and σ representing the density of liquid, the gravitational acceleration, the diameter of the object, and the surface tension of liquid, respectively). They also proposed a regime map for the water entry of a spherical projectile with a hydrophobic surface based on theoretical and experimental results. Subsequent studies have also shown that the shape of the cavity and its growth exhibit complex variations depending on the liquid properties (density, viscosity, and surface tension) (Aristoff and Bush, 2009), the wettability of the object’s surface (Duez et al., 2007), the solid–liquid density ratio (Aristoff et al., 2010), the object geometry (Wei and Hu, 2015; Kim and Park, 2019), and the initial impact conditions (velocity, angle, and transverse spin) (Truscott and Techet, 2009; Techet and Truscott, 2011). The majority of previous experimental studies have been performed at high capillary numbers (Ca = μU/σ, with μ, U, and σ representing the viscosity of fluid, the impact velocity of the projectile, and the surface tension of liquid, respectively), which tend to be accompanied by cavity formation behind the entering projectile. Among these, some studies have focused on the critical Ca numbers that lead to cavity formation. Duez et al. (2007) experimentally demonstrated that the cavity-formation tendency of a spherical projectile is affected by the contact angle (θc). This study found that the critical Ca number of hydrophilic spheres has a weak dependency on the contact angle, while that of hydrophobic spheres is strongly proportional to θc3. This was strongly supported by consistency between the theoretical model and the experimental results of this study. Subsequent studies have found that a cavity can be created slightly below the critical Ca limit within a narrow transition region.

The recent advancement of computing technology and development in computational fluid dynamics (CFD) techniques enables us to simulate various water entry problems, which involves unsteady complex multiphase flow and fluid–solid interaction (FSI). Many researchers have introduced grid-based numerical methods, including the Boundary Element Method (BEM) (Xu et al., 2008), Finite Element Method (FEM) (Do-Quang and Amberg, 2009), and Finite Volume Method (FVM) (Krastev et al., 2018), to study the water entry problems. Do-Quang and Amberg (2009) conducted numerical simulations on the water entry of a projectile using a two-dimensional (2D) axi-symmetric finite element method and introduced a phase-field model to include the effects of capillary forces. They obtained results that matched well with the theory presented by Duez et al. (2007). In their study, a solid spherical projectile was modeled as a liquid with a very high viscosity, which enabled control of the surface conditions of the object. Kintea et al. (2016) investigated the water entry of a rotating sphere and non-rotating sphere with different impact angles using the finite volume method. The Volume Of Fluid (VOF) is utilized to model the free surface flow, and the spherical object is assumed to be a perfectly rigid sphere. Krastev et al. (2018) coupled the VOF model with a Laplacian mesh-motion solver to capture the free surface and numerically studied the asymmetric water impact of a two-dimensional wedge.

Smoothed particle hydrodynamics (SPH) is a meshless Lagrangian numerical method that was first developed in astrophysics by Gingold and Monaghan (1977) and is gaining popularity in many science and engineering fields (Zhang et al., 2017; Jo et al., 2019; and Ye et al., 2019). Due to its Lagrangian nature, SPH is very effective for handling phenomena with large deformations, such as free surfaces and multiple phases. Therefore, SPH has been used as a powerful tool for studying water entry problems by a number of recent researchers. Oger et al. (2006) first adapted the weakly compressible smoothed particle hydrodynamics (WCSPH) to simulate the water entry problem. A two-dimensional wedge-entry problem is studied using a variable smoothing length technique. In a subsequent study, they coupled FEM with WCSPH to deal with the fluid–solid interaction and captured the free surface by using a single-phase solver (Fourey et al., 2010). Skillen et al. (2013) investigated the water entry of a 2D cylinder and wedge using improved single-phase Incompressible SPH (ISPH). Lind et al. (2015) performed numerical simulations of the impact of a 2D rigid plate by coupling ISPH for the water phase and WCSPH for the air phase to study the effect of two phases. Gong et al. (2016) conducted numerical simulations of water entry of a 2D wedge using two-phase WCSPH and obtained consistent results with their experimental results. Recently, Sun et al. (2018) carried out 2D and three-dimensional (3D) numerical simulations using single-phase WCSPH. The Adaptive Particle Refinement (APR) technique and δ+-SPH scheme are adopted to obtain accurate results. Most of the previous studies of the water entry problem are 2D cases and the 3D studies with quantitative analysis are insufficient. In particular, 3D numerical analyses that deal with the complex two-phase flow characteristics and the effect of the surface wettability are very rare.

In the present study, 3D multi-phase simulations of the projectile water entry are performed for a wide range of surface wettability conditions using smoothed particle hydrodynamics. It demonstrates that the proposed numerical method reproduces the shapes and regimes of cavity formation in an adequate way. This paper is organized as follows: In Sec. II, the governing equations and numerical formulations are explained with boundary conditions. The continuity and momentum equations are reformulated based on a normalized-density to handle high-density ratios in multi-phase flows. Then, a surface tension model based on a physics-based pair-wise force is adopted to treat various surface wettability conditions. In Sec. III, the presented numerical models are validated with theoretical predictions. In addition, the results of the simulations are presented and validated by comparisons with the experimental observations and the theoretical predictions. The characteristics of the flow and the cavity formation are investigated in detail in a quantitative and qualitative manner. Finally, the effect of anisotropy in surface wettability on a sphere trajectory is simulated. The motion of the spherical projectile affected by the anisotropic flow is verified and analyzed qualitatively. Section IV summarizes and concludes this paper with final discussions.

In SPH, the fluid system is represented by Lagrangian nodes (or particles). The physical properties and the quantities of the nodes are approximated by mathematical interpolations using the neighboring nodes (or particles). The general SPH approximation is mathematically expressed by (Gingold and Monaghan, 1977)

fri=jfjrjWijrirj,hVj=jfrjWij|rirj|,hmjρj,
(1)

where r denotes the point vector, fri is a function at the position of center particle i, and subscript j represents the neighboring particles of the center particle i. Vj, mj, and ρj are the volume, mass, and density of the jth particle, respectively. Wij|rirj|,h denotes a kernel function that determines the weighting between particles i and j. The quantity h refers to the smoothing length that is chosen to be 2.5Δx in our study, where Δx is the initial particle spacing.

In the similar way, the gradient, divergence, and Laplacian of a function can be approximated as (Liu and Liu, 2010).

fri=jfrjiWijρ|rirj|,hmjρj,
(2)
fri=jfrjfriiWij|rirj|,hmjρj,
(3)
2fri=j2frifrjrij2rijiWij|rirj|,hmjρj.
(4)

The basic governing equations for an isothermal, incompressible flow in a Lagrangian frame are the mass conservation in Eq. (5) and the momentum conservation in Eq. (6),

dρdt=ρu,
(5)
dudt=Pρ+ν2u+g+fs,
(6)

where t is time and P, ν, g, and fs are the pressure, kinematic viscosity, gravitational acceleration vector, and surface tension force.

Equation (5) is a commonly used SPH continuity equation form, but the density can also be obtained directly from the SPH summation of the mass of all particles in the support domain (Monaghan, 1992). From the general SPH approximation in Eq. (1), the density estimation can be written as

ρi=jmjWij.
(7)

The density estimation using Eq. (7) is the most commonly used formulation in the WCSPH and able to handle the hydrodynamics problem with single-phase flow. However, the non-physical pressure force can be generated at the interface between multi-phase flows with a large-density ratio due to the over-smoothing of the density. This non-physical pressure force induces serious error/instability near the interface (Hu and Adams, 2006). Therefore, the reformulated density estimation based on the normalized-density is adopted to handle multiphase flow without numerical instability (Park et al., 2020),

ρi=ρ0,ijmjρ0,jWij,
(8)

where ρ0 is the reference density of the fluid. The momentum conservation in Eq. (6) is discretized with the SPH approximation as Morris et al. (1997),

DuDti=jmjPi+PjρiρjiWij+j4mjρiρjμiμjμi+μjrijiWijrij2+η2uiuj+g+fs.
(9)

Each term represents the pressure force, viscous force, gravitational force and surface tension force. The surface tension model is explained in detail in Subsection II C.

In the WCSPH, the fluid is assumed to be weakly compressible, which allows for slight compressibility. The pressure is explicitly calculated with the equation of state (EOS) based on the compressibility of density. The most common EOS used in WCSPH is the Tait equation of state written as (Batchelor and Batchelor, 1967)

P=c02ρ0γρρ0γ1,
(10)

where c0 is the speed of sound and γ denotes the polytropic constant that determines the sensitivity of the pressure to compressibility and normally have the value of 7 for water.

The surface tension is one of the most important forces in multi-phase flow, particularly at the small scale. Physically, the surface tension is generated by the irregular distribution of molecules near the interface, and the imbalance of the molecular force induces discontinuous pressure jump at the interface with non-zero curvature. The most commonly used surface tension models in SPH are the continuum surface force (CSF) model and the pairwise force (PF) model. The CSF model was proposed by Brackbill et al. (1992) for the grid-based methods. Morris (2000) first introduced the CSF model in SPH, and several algorithms are proposed to calculate the normal vector and curvature at the interface (Adami et al., 2010),

fs=σρiκini,
(11)

where σ is the surface tension and κi and ni denote the curvature and surface-normal vector.

The PF model was proposed by Tartakovsky and Meakin (2005) based on the microscopic origin of surface tension. In the natural world, a non-uniform distribution of molecules generates a force imbalance near the phase interface, and thus, the stress tensor becomes anisotropic. Therefore, a certain amount of energy must be added when the surface area is increased by a unit length, which is the mechanical definition of surface tension (Rowlinson and Widom, 2002),

σαβ=+[TτzTn(z)]dz,
(12)

where σαβ is the surface tension between phases α and β, Tτ and Tn are the tangential and normal stresses, respectively, and the z direction is normal to the interface. Using the pairwise force terms from the microscopic origin of surface tension, the PF model induces anisotropy in the force equilibrium of the particles near the interphase region. Several forms of PF models are suggested, and in the present study, we used the model below (Tartakovsky and Meakin, 2005),

fij=sijcos1.5πkhrjirjirji,
(13)

where sij is a parameter representing the surface characteristics. This force generates attraction between the particles at a distance and repulsion in proximity. The surface tension and the contact angle are controlled by sij. Tartakovsky and Panchenko (2016) mathematically derived this parameter for macroscopic wetting characteristics (surface tension and contact angle) using the Hardy formula and the mechanical definition of surface tension in Eqs. (14) and (15) (Hardy, 1982),

σαβ=λnα2sαα+nβ2sββ2nαnβsαβ,
(14)
cosθ0=sββ̄sαᾱ2sβs̄+2sαs̄sαᾱ+sββ̄2sαβ̄.
(15)

Here, n is the number density, λ is a constant that can be calculated mathematically, and the subscripts α, β, and s denote the fluid phase α, β and solid phase, respectively. The validity of these formulations are checked by several numerical validations (Bandara et al., 2013; Nair and Pöschel, 2018; Tartakovsky and Panchenko, 2016; and Arai et al., 2020).

Many researchers have proposed different boundary treatment techniques in the SPH method. Some researchers exerted the forces on the fluid particles adjacent to the boundary particles to prevent non-physical penetration (Monaghan and Kajtar, 2009). The mirror boundary condition or dynamic boundary condition (DBC) is also commonly used for solid boundary treatment (Bierbrauer et al., 2009; Crespo et al., 2007). In this paper, the DBC is applied to model the solid boundary and the rigid object. In this model, the boundaries are composed of fixed particles that follow the same density estimation and EOS as the fluid particles. Without special numerical treatment, the penetration of fluid particles is prevented and the computational accuracy is improved by compensating support domain of fluid particles. The momentum and energy conservations were checked from a single-particle impinging test and dam break simulation by Crespo et al. (2007).

A no-slip boundary condition, which is commonly used in fluid dynamics simulations, can cause stress-singularity problems near the three-phase interface (fluid-fluid-solid) when solving contact-line dynamics at the interface (Ren and E, 2007; Bonn et al., 2009). However, it could be forced to some extent in other areas while not forcing the no-slip boundary condition in the vicinity of the three-phase line. In this perspective, the boundary conditions are set as follows, with reference to previous studies. The velocity of the solid particles is forced to be zero to apply the no-slip boundary condition only at fluid–solid interfaces, instead of forcing the no-slip boundary condition in the entire region. Therefore, the slip motion of three-phase line is allowed and the no-slip boundary condition is implicitly satisfied in the other areas. Under this boundary condition, the dynamic motion of three-phase flows are validated in several papers for single-phase and multi-phase flows (Kordilla et al., 2013; Shigorina et al., 2017; Tartakovsky and Meakin, 2005; and Tartakovsky and Panchenko, 2016).

The rigid objects with isotropic wettability are submerged at a constant velocity, and forces applied to solid body from fluid are not considered. Due to the short simulation domain and the small size of the objects, it is reasonable to assume that the rigid objects have a constant velocity. However, in the case of a rigid object with anisotropic wettability, asymmetric flow is induced and it has a great influence on the behavior of the object. In this case, a two-way SPH-rigid body coupling is introduced in order to reflect the momentum applied to the sphere from the flow. The sphere is assumed to be a perfect rigid object without deformation. Figure 1 represents a conceptual sketch of the two-way SPH-rigid body coupling method. The fluid particles interact with all neighboring particles and the rigid particles interact only with neighboring fluid particles. The total momenta applied to rigid particles are reflected in the linear and angular motion of the rigid object. Based on the Newton’s second law of motion, the linear and angular motion equation of the rigid object are written as

v=v0+fimidt,
(16)
w=w0+I1τdt,
(17)

where v, v0, w, and w0 are the velocity, initial velocity, angular velocity, and initial angular velocity, fi denotes a total force acting on the ith solid particle from neighboring fluid particles, and I and τ are the momentum of inertia tensor and torque vector, respectively.

FIG. 1.

The 2D schematic of the two way SPH-rigid body coupling.

FIG. 1.

The 2D schematic of the two way SPH-rigid body coupling.

Close modal

In this section, the verifications of wetting characteristics modeling are carried out. Due to the surface tension, the pressure jump is generated between the curved fluid–fluid interface, which is called the Laplace pressure, and it reads

PinPout=σκ=σ1R1+1R2,
(18)

where R1 and R2 are the principal radii of curvature. The pressure distributions of droplets in equilibrium are calculated in order to confirm whether the expected pressure jump is well reproduced using the numerical models. Figure 2 shows the particle configuration of the 1 mm sized droplets at the equilibrium state using the different surface tension models with the pressure distribution across the x-axis direction. The density ratio of the inner phase and outer phase is 1.0, and the colored area in the graph means ±15% of the theoretical pressure jump. As shown in Fig. 2, the Laplace pressure jump is well reproduced using both models. A narrow gap at the interface has shown in the PF model, which was reported by Tartakovsky and Panchenko (2016). The large pressure jumps are shown in the PF model due to its microscopic origin, but it was found that the model follows the macroscopic behavior when the curvature radius is larger than the kernel support (Arai et al., 2020). In the recent study by Howard and Tartakovsky (2020), it is found that the PF model shows better convergence properties than the CSF for the droplet collision problem. In addition, the PF method does not require any additional numerical models for normal vector calculation and adhesion force, which are very strong points compared to the CSF model (Arai et al., 2020). Therefore, the PF model is applied in this study.

FIG. 2.

The SPH 3D snapshots and 2D cross-sectional view of the particle distribution and pressure profile of droplets (R0 = 1.0 mm) at the equilibrium state using different surface tension models.

FIG. 2.

The SPH 3D snapshots and 2D cross-sectional view of the particle distribution and pressure profile of droplets (R0 = 1.0 mm) at the equilibrium state using different surface tension models.

Close modal

The adaptability of the PF model is demonstrated through the prediction of the Laplace pressure jump at the interface of a spherical droplet with different spatial resolutions, smoothing resolutions, sizes, and the density ratios. The pressure profile of a droplet at the equilibrium state with different numerical conditions is shown in Fig. 3. In Fig. 3(a), pressure distributions of droplets with the radius of 1 mm are calculated with different density ratios. It has been shown that the model works well at the multi-phase problem with high-density ratios. The theoretical pressure jumps are indicated with solid lines and ±15% of the theoretical pressure differences are indicated with colored regions. In Fig. 3(b), pressure distributions of the same sized droplets (R0 = 1.0 mm) with different spatial resolution are shown. The pressure differences agree well with the theoretical Laplace pressure for low resolution but get unstable with larger resolutions. The same calculations are performed with a larger smoothing resolution (h/Δx = 2.5) in Fig. 3(c). Stable pressure distributions that well agree with the theoretical Laplace pressure are shown regardless of the spatial resolution. In Fig. 3(d), pressure distributions of the droplets with different radii with the same spatial resolution (Δx = 0.1 mm) and smoothing resolution (h/Δx = 2.5) are shown. The theoretical pressure jumps are well reproduced in both cases. From the results, the smoothing resolution is set to the large one (h/Δx = 2.5) for the water entry simulations.

FIG. 3.

The pressure profile of droplets at the equilibrium state with different numerical conditions: (a) density ratios, (b) spatial resolution (h/Δx = 1.5), and (c) spatial resolution (h/Δx = 2.5) (d) radius.

FIG. 3.

The pressure profile of droplets at the equilibrium state with different numerical conditions: (a) density ratios, (b) spatial resolution (h/Δx = 1.5), and (c) spatial resolution (h/Δx = 2.5) (d) radius.

Close modal

The contact angle and wetting characteristics can be controlled by varying the relative pairwise force parameters according to the pairwise force smoothed particle hydrodynamics (PFSPH) theory. A droplet is placed in a solid substrate, and the contact angle is measured at equilibrium. The relative parameters are initially set so that the static contact angles are adjusted between 30° and 150° at intervals of 30°. In Fig. 4, the particle images of droplets with different contact angles at equilibrium are displayed. In the top part of the figure, the expected contact angles and the measured contact angles at x and y directions are compared. The results show that the wetting characteristics can be well reproduced.

FIG. 4.

The SPH snapshots of droplets with different static contact angles.

FIG. 4.

The SPH snapshots of droplets with different static contact angles.

Close modal

Figure 5 shows the input model of the water entry simulation. The simulation domain is divided into three regions: liquid, gas, and projectile. White, yellow, and blue particles denote gas, solid (projectile), and liquid particles, respectively. In the simulation, the solid projectile is introduced to the water vertically with a predefined speed. The Reynolds number (Re = ρvD/μ) is specified in the range 30–400. The density and viscosity of the liquid are set to be 1000 kg/m3 and 0.005 Pa s, respectively. The gas density and the viscosity are set to be 1 kg/m3 and 10−5 Pa s, respectively. The surface tension of the liquid phase is 0.073 N/m and the contact angle between the liquid and the sphere is adjusted from 0° (complete-wetting) to 180° (non-wetting). The properties of the solid object are set to be the same as the liquid phase, and the particles are initially uniformly distributed. In order to determine the appropriate particle spacing, the water entry of a hydrophilic sphere (θ = 0°) is simulated with three different spatial resolutions (D/Δx = 15, 20, 25, and 30). The full wetting spherical object with the constant speed of 1.5 m/s is completely submerged without the cavity. The relative position of the three-phase line is tracked until the sphere is completely submerged in water and plotted in Fig. 6. At the smallest resolution (D/Δx = 15), some staircase effect has been observed. The results of the rest three cases are in good agreement, and therefore, the particle spacing is set to be D/Δx = 20 in the subsequent simulations.

FIG. 5.

Initial arrangement of particles.

FIG. 5.

Initial arrangement of particles.

Close modal
FIG. 6.

Time evolution of the relative position of the triple line.

FIG. 6.

Time evolution of the relative position of the triple line.

Close modal

Figure 7 shows snapshots of the simulation for hydrophilic and hydrophobic spheres. Images in the upper row show a hydrophilic sphere (θ = 0°) submerging at a constant speed of 1.5 m/s, and images in the lower row show a hydrophobic sphere (θ = 120°) submerging at a constant speed of 1.5 m/s. In the case of the hydrophilic sphere, as the projectile is submerged in the water, a large stream surrounding the sphere is formed, which causes the flow to converge along its surface. After the projectile is completely submerged, the surface above the sphere rises from the converged flow and a thick jet is formed at the top center. As the body sinks, the accumulated flow over the sphere and the dragged flow near the rear pole induce a stagnation point. The flow field and the behavior of the free surface are similar with the FEM results by Do-Quang and Amberg (2009). Due to the different range of Reynolds number, a thinner jet is generated in their study. In addition, a tiny bubble is captured at the rear end of the projectile because they performed very high-resolution FEM simulation with a 2D axisymmetric mesh. However, the overall progression of the free surface and the time evolution of flow field were very similar to our study. The generation of the stagnation point and jet is also well explained in Kim and Park (2019). Conversely, the simulation with a hydrophobic sphere showed very different results. Small and complicated vortices are formed near the three-phase contact line as the body submerges, as shown in Fig. 8. These vortices destruct the progress of the three-phase interface and induce the formation of a cavity. The generated cavity expands downward as the sphere enters the water and is narrowed by hydraulic pressure. It can be seen in Fig. 7 that the air velocity increases as the width of the air cavity decreases. Accordingly, the pressure is reduced by conservation of energy, and further reduction of the cavity width is promoted until pinch-off (Rayleigh–Plateau instability). The similar sequential process of the flow field is shown in Do-Quang and Amberg (2009). In their study, the generation of small vortices near the contact line and its role on stability of the triple line are well explained (see videos 1–3 of the supplementary material).

FIG. 7.

Snapshots of simulations of the entry of a projectile into water.

FIG. 7.

Snapshots of simulations of the entry of a projectile into water.

Close modal
FIG. 8.

Vorticity and streamline visualizations of (a) a submerging hydrophilic sphere (θ = 0°) at 1.5 m/s and (b) a submerging hydrophobic sphere (θ = 120°) at 1.5 m/s (t* = 1.5).

FIG. 8.

Vorticity and streamline visualizations of (a) a submerging hydrophilic sphere (θ = 0°) at 1.5 m/s and (b) a submerging hydrophobic sphere (θ = 120°) at 1.5 m/s (t* = 1.5).

Close modal

In Fig. 9, the dimensionless position of the three-phase line over dimensionless time is plotted for different velocities and contact angles. The position of the triple line is non-dimensionalized with an angle relative to the ball. Filled symbols represent spheres with a constant impact speed of 1.0 m/s, and empty symbols represent spheres with an constant impact speed of 1.5 m/s. The result confirms that the progress of the triple line is slowed as the hydrophobicity of the ball is increased under conditions of the same capillary number. Regardless of the capillary number, the non-dimensional positions of the triple line at 0°, 30°, and 60° change consistently with non-dimensional time. However, inconsistent result is shown in the case of 90° and deviation from the trend can be seen, especially for at a velocity of 1.5 m/s. This case is near the critical capillary number, where the three-phase line starts to be instable. Accordingly, it starts to gradually deviate from the general tendency, and further increase in the capillary number leads to greater instability of the triple line. The instable triple line allows the invasion of gas, which results in cavity formation and a complete departure from the general tendency. The motions of the triple line of projectiles with the wide range of wettability are well shown and explained by Do-Quang and Amberg (2009).

FIG. 9.

Position of the triple lines of falling spheres with different wettabilities over non-dimensional time.

FIG. 9.

Position of the triple lines of falling spheres with different wettabilities over non-dimensional time.

Close modal

Figure 10 shows comparisons of the critical Ca number between the simulation results and the theoretical regime map proposed by Duez et al. (2007). The solid line is the theoretical prediction of the critical Ca number that separates the area without and with cavity formation. Empty symbols denote submerging objects with no cavity, and filled symbols denote submerging objects with a cavity. As shown in this figure, the critical Ca numbers of the simulations agree well with the theoretical line and the previous studies (Duez et al., 2007; Li et al., 2019). For the hydrophilic sphere (θc < 90°), the critical Ca number shows a weak dependence on the contact angle, and for the hydrophobic sphere (θc > 90°), the critical Ca number dramatically decreases with increasing contact angle. Figures 10(a) and 10(d) show the simulation results, and Fig. 10(e) shows the experimental observation of three different cavity shapes in different regimes. The triangle, asterisk, and diamond symbols denote the different regimes of cavity type (quasi-static, shallow seal, and deep seal). Figure 10(a) shows a hydrophilic sphere submerging without a cavity, and Fig. 10(b) shows the impact of a hydrophobic sphere with a small Ca number. A shallow cavity is formed, the size of which is of the order of the capillary length. As the sphere is immersed in the water, the free surface is deformed by the surface tension and dragged downward, resulting in a pinch-off at a depth of the order of the capillary length. This cavity with small air entrainment is called a “quasi-static cavity,” which is also shown in Fig. 10(e). In Figs. 10(c) and 10(d), the hydrophilic and hydrophobic spheres both submerge in the water with large air cavities. The cavities expand downward because of the large inertia of the sphere, and little Rayleigh–Taylor instability is observed near the free surface. The cavity shown in Fig. 10(d) is classified as belonging to the regime between the shallow-seal region and the deep-seal region, which are shown in Fig. 10(e). The exact regime cannot be clearly classified because of the short simulation domain. The shapes of the cavities formed by the hydrophobic projectiles in Figs. 10(b) and 10(d) are consistent with the previous study performed by Aristoff and Bush (2009), which categorized the cavity types generated by hydrophobic spherical objects according to the conditions of the Froude and Bond numbers.

FIG. 10.

Critical Ca numbers for spherical projectiles of different wettabilities and different types of cavity shapes with experimental observation: (a) hydrophilic sphere (θ = 0°) at 1.5 m/s,1 (b) hydrophobic sphere (θ = 150°) at 0.15 m/s, (c) hydrophilic sphere (θ = 30°) at 2.0 m/s, (d) hydrophobic sphere (θ = 120°) at 1.5 m/s, and (e) experimental observation. Reproduced with permission from Aristoff, J. M. and Bush, J. W. M., J. Fluid Mech. 619, 45–78 (2009). Copyright 2009 Cambridge University Press.

FIG. 10.

Critical Ca numbers for spherical projectiles of different wettabilities and different types of cavity shapes with experimental observation: (a) hydrophilic sphere (θ = 0°) at 1.5 m/s,1 (b) hydrophobic sphere (θ = 150°) at 0.15 m/s, (c) hydrophilic sphere (θ = 30°) at 2.0 m/s, (d) hydrophobic sphere (θ = 120°) at 1.5 m/s, and (e) experimental observation. Reproduced with permission from Aristoff, J. M. and Bush, J. W. M., J. Fluid Mech. 619, 45–78 (2009). Copyright 2009 Cambridge University Press.

Close modal

A simulation of the water entry of a sphere with a hydrophobic surface (θ = 150°) on its right hemisphere and a hydrophilic surface (θ = 0°) on its left hemisphere is performed. As mentioned in Sec. II D, a two-way SPH-rigid body coupling is introduced to reflect the momentum applied to the spherical object from the fluid. The density of rigid particles is five times of fluid particles and it is submerged at an initial speed of 1.5 m/s. It is confirmed that asymmetric cavities are created, as shown in Fig. 11. Asymmetric progression of the triple line induces asymmetrical flow and vortices, and the path of the ball is bent toward the hydrophilic side. The motion of the submerging ball and the shape of the cavity are similar to the experimental results (Truscott and Techet, 2009) as in Figs. 11 and 12. Figure 12(a) presents a series of snapshots of experimental study by Truscott and Techet (2009). The sphere projectile with a half-and-half surface treatment is introduced with 1.72 m/s. The left hemisphere is coated with the hydrophobic surface coating and the right hemisphere is hydrophilic. In Fig. 12(b), the numerical simulation in our study is visualized from top-view. The similar flow characteristics between both cases are shown, but no more quantitative analysis is discussed in this paper (see video 4 of the supplementary material).

FIG. 11.

Vorticity and streamline visualizations of a submerging sphere with asymmetric wetting at 1.5 m/s.

FIG. 11.

Vorticity and streamline visualizations of a submerging sphere with asymmetric wetting at 1.5 m/s.

Close modal
FIG. 12.

Comparison between experimental results and numerical study. (a) A series of snapshots of experimental study. Reproduced from Truscott, T. T. and Techet, A. H., Phys. Fluids 21(12), 121703 (2009) with the permission of AIP Publishing. (b) Top-view of numerical simulation [correspond to the first two snapshots in (a)].

FIG. 12.

Comparison between experimental results and numerical study. (a) A series of snapshots of experimental study. Reproduced from Truscott, T. T. and Techet, A. H., Phys. Fluids 21(12), 121703 (2009) with the permission of AIP Publishing. (b) Top-view of numerical simulation [correspond to the first two snapshots in (a)].

Close modal

In this study, 3D multi-phase numerical simulations using the WCSPH method are performed for the water entry of spherical projectiles with a wide range of surface wettability conditions. Through detailed analyses, the process of water entry is carefully examined from a hydrodynamic point of view, and it is confirmed that the results are quantitatively and qualitatively consistent with those of previous studies.

The normalized-density based governing equations are solved to treat multi-phase flows with high-density ratios, and the pairwise force is adopted to handle the contact line dynamics. The DBC is adopted for fluid–solid interaction. The accuracy of the PFSPH on modeling wetting characteristics is proved through several numerical verification cases. The motion of the triple line is tracked and the flow characteristics are analyzed in detail. The influence of the wettability of the sphere and the capillary number on the cavity formation is investigated, and the shapes of the cavities are similar to experimental observation. In addition, the critical capillary number where the triple line starts to be unstable and induce the cavity formation is consistent with theoretical prediction and experimental result.

The two-way fluid-rigid body coupling is implemented in this paper to simulate the anisotropic flow behavior from the water entry of a spherical object with anisotropic wettability. The asymmetrical flow characteristics and the shape of the cavity are similar to experimental results.

Many previous numerical analyses of the water entry problem have been conducted, but the 3D quantitative analysis using Lagrangian-based numerical analyses is insufficient until recently. This study confirmed that the multiphase WCSPH method and the pairwise force based on the microscopic origin can accurately simulate the contact-line dynamics and the overall water entry process. Based on these results, it is expected that the methods of this study can be used for improving the understanding of the water entry phenomena and resultant flow characteristics. In addition, an extended study at the large Reynolds number and various wettability conditions will be possible without difficulties.

See the supplementary material videos for more detailed visualizations of our work.

This work was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIT) (Grant No. 2019R1F1A1060159).

The authors have no conflicts to disclose.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

1.
Adami
,
S.
,
Hu
,
X. Y.
, and
Adams
,
N. A.
, “
A new surface-tension formulation for multi-phase SPH using a reproducing divergence approximation
,”
J. Comput. Phys.
229
(
13
),
5011
5021
(
2010
).
2.
Arai
,
E.
,
Tartakovsky
,
A.
,
Holt
,
R. G.
,
Grace
,
S.
, and
Ryan
,
E.
, “
Comparison of surface tension generation methods in smoothed particle hydrodynamics for dynamic systems
,”
Comput. Fluids
203
,
104540
(
2020
).
3.
Aristoff
,
J. M.
and
Bush
,
J. W. M.
, “
Water entry of small hydrophobic spheres
,”
J. Fluid Mech.
619
,
45
78
(
2009
).
4.
Aristoff
,
J. M.
,
Truscott
,
T. T.
,
Techet
,
A. H.
, and
Bush
,
J. W. M.
, “
The water entry of decelerating spheres
,”
Phys. Fluids
22
(
3
),
032102
(
2010
).
5.
Bandara
,
U. C.
,
Tartakovsky
,
A. M.
,
Oostrom
,
M.
,
Palmer
,
B. J.
,
Grate
,
J.
, and
Zhang
,
C.
, “
Smoothed particle hydrodynamics pore-scale simulations of unstable immiscible flow in porous media
,”
Adv. Water Resour.
62
,
356
369
(
2013
).
6.
Batchelor
,
C. K.
and
Batchelor
,
G. K.
,
An Introduction to Fluid Dynamics
(
Cambridge University Press
,
1967
).
7.
Bierbrauer
,
F.
,
Bollada
,
P. C.
, and
Phillips
,
T. N.
, “
A consistent reflected image particle approach to the treatment of boundary conditions in smoothed particle hydrodynamics
,”
Comput. Methods Appl. Mech. Eng.
198
(
41–44
),
3400
3410
(
2009
).
8.
Bonn
,
D.
,
Eggers
,
J.
,
Indekeu
,
J.
,
Meunier
,
J.
, and
Rolley
,
E.
, “
Wetting and spreading
,”
Rev. Mod. Phys.
81
(
2
),
739
(
2009
).
9.
Brackbill
,
J. U.
,
Kothe
,
D. B.
, and
Zemach
,
C.
, “
A continuum method for modeling surface tension
,”
J. Comput. Phys.
100
(
2
),
335
354
(
1992
).
10.
Crespo
,
A. J. C.
,
Gómez-Gesteira
,
M.
, and
Dalrymple
,
R. A.
, “
Boundary conditions generated by dynamic particles in SPH methods
,”
Comput., Mater. Continua
5
(
3
),
173
184
(
2007
).
11.
Ding
,
H.
,
Chen
,
B.-Q.
,
Liu
,
H.-R.
,
Zhang
,
C.-Y.
,
Gao
,
P.
, and
Lu
,
X.-Y.
, “
On the contact-line pinning in cavity formation during solid–liquid impact
,”
J. Fluid Mech.
783
,
504
525
(
2015
).
12.
Do-Quang
,
M.
and
Amberg
,
G.
, “
The splash of a solid sphere impacting on a liquid surface: Numerical simulation of the influence of wetting
,”
Phys. Fluids
21
(
2
),
022102
(
2009
).
13.
Duez
,
C.
,
Ybert
,
C.
,
Clanet
,
C.
, and
Bocquet
,
L.
, “
Making a splash with water repellency
,”
Nat. Phys.
3
(
3
),
180
183
(
2007
).
14.
Fourey
,
G.
,
Oger
,
G.
,
Touzé
,
D. L.
, and
Alessandrini
,
B.
, “
Violent fluid-structure interaction simulations using a coupled SPH/FEM method
,”
IOP Conf. Ser.: Mater. Sci. Eng.
10
(
1
),
012041
(
2010
).
15.
Gilbarg
,
D.
and
Anderson
,
R. A.
, “
Influence of atmospheric pressure on the phenomena accompanying the entry of spheres into water
,”
J. Appl. Phys.
19
(
2
),
127
139
(
1948
).
16.
Gingold
,
R. A.
and
Monaghan
,
J. J.
, “
Smoothed particle hydrodynamics: Theory and application to non-spherical stars
,”
Mon. Not. R. Astron. Soc.
181
(
3
),
375
389
(
1977
).
17.
Gong
,
K.
,
Shao
,
S.
,
Liu
,
H.
,
Wang
,
B.
, and
Tan
,
S.-K.
, “
Two-phase SPH simulation of fluid–structure interactions
,”
J. Fluids Struct.
65
,
155
179
(
2016
).
18.
Hardy
,
R. J.
, “
Formulas for determining local properties in molecular‐dynamics simulations: Shock waves
,”
J. Chem. Phys.
76
(
1
),
622
628
(
1982
).
19.
Howard
,
A. A.
and
Tartakovsky
,
A. M.
, “
Non-local model for surface tension in fluid-fluid simulations
,”
J. Comput. Phys
421
,
109732
(
2020
).
20.
Hu
,
X. Y.
and
Adams
,
N. A.
, “
A multi-phase SPH method for macroscopic and mesoscopic flows
,”
J. Comput. Phys.
213
(
2
),
844
861
(
2006
).
21.
Jo
,
Y. B.
,
Park
,
S.-H.
,
Choi
,
H. Y.
,
Jung
,
H.-W.
,
Kim
,
Y.-J.
, and
Kim
,
E. S.
, “
SOPHIA: Development of Lagrangian-based CFD code for nuclear thermal-hydraulics and safety applications
,”
Ann. Nucl. Energy
124
,
132
149
(
2019
).
22.
Kim
,
N.
and
Park
,
H.
, “
Water entry of rounded cylindrical bodies with different aspect ratios and surface conditions
,”
J. Fluid Mech.
863
,
757
788
(
2019
).
23.
Kintea
,
D. M.
,
Breitenbach
,
J.
,
Thammanna Gurumurthy
,
V.
,
Roisman
,
I. V.
, and
Tropea
,
C.
, “
On the influence of surface tension during the impact of particles on a liquid-gaseous interface
,”
Phys. Fluids
28
(
1
),
012108
(
2016
).
24.
Kordilla
,
J.
,
Tartakovsky
,
A. M.
, and
Geyer
,
T.
, “
A smoothed particle hydrodynamics model for droplet and film flow on smooth and rough fracture surfaces
,”
Adv. Water Resour.
59
,
1
14
(
2013
).
25.
Krastev
,
V. K.
,
Facci
,
A. L.
, and
Ubertini
,
S.
, “
Asymmetric water impact of a two dimensional wedge: A systematic numerical study with transition to ventilating flow conditions
,”
Ocean Eng.
147
,
386
398
(
2018
).
26.
Li
,
D.
,
Zhang
,
J.
,
Zhang
,
M.
,
Huang
,
B.
,
Ma
,
X.
, and
Wang
,
G.
, “
Experimental study on water entry of spheres with different surface wettability
,”
Ocean Eng.
187
,
106123
(
2019
).
27.
Lind
,
S. J.
,
Stansby
,
P. K.
,
Rogers
,
B. D.
, and
Lloyd
,
P. M.
, “
Numerical predictions of water–air wave slam using incompressible–compressible smoothed particle hydrodynamics
,”
Appl. Ocean Res.
49
,
57
71
(
2015
).
28.
Liu
,
M. B.
and
Liu
,
G. R.
, “
Smoothed particle hydrodynamics (SPH): An overview and recent developments
,”
Arch. Comput. Methods Eng.
17
(
1
),
25
76
(
2010
).
29.
Marston
,
J. O.
,
Truscott
,
T. T.
,
Speirs
,
N. B.
,
Mansoor
,
M. M.
, and
Thoroddsen
,
S. T.
, “
Crown sealing and buckling instability during water entry of spheres
,”
J. Fluid Mech.
794
,
506
529
(
2016
).
30.
May
,
A.
and
Woodhull
,
J. C.
, “
The virtual mass of a sphere entering water vertically
,”
J. Appl. Phys.
21
(
12
),
1285
1289
(
1950
).
31.
Monaghan
,
J. J.
, “
Smoothed particle hydrodynamics
,”
Annu. Rev. Astron. Astrophys.
30
(
1
),
543
574
(
1992
).
32.
Monaghan
,
J. J.
and
Kajtar
,
J. B.
, “
SPH particle boundary forces for arbitrary boundaries
,”
Comput. Phys. Commun.
180
(
10
),
1811
1820
(
2009
).
33.
Morris
,
J. P.
, “
Simulating surface tension with smoothed particle hydrodynamics
,”
Int. J. Numer. Methods Fluids
33
(
3
),
333
353
(
2000
).
34.
Morris
,
J. P.
,
Fox
,
P. J.
, and
Zhu
,
Y.
, “
Modeling low Reynolds number incompressible flows using SPH
,”
J. Comput. Phys.
136
(
1
),
214
226
(
1997
).
35.
Nair
,
P.
and
Pöschel
,
T.
, “
Dynamic capillary phenomena using incompressible SPH
,”
Chem. Eng. Sci.
176
,
192
204
(
2018
).
36.
Oger
,
G.
,
Doring
,
M.
,
Alessandrini
,
B.
, and
Ferrant
,
P.
, “
Two-dimensional SPH simulations of wedge water entries
,”
J. Comput. Phys.
213
(
2
),
803
822
(
2006
).
37.
Park
,
S.-H.
,
Jo
,
Y. B.
,
Ahn
,
Y.
,
Choi
,
H. Y.
,
Choi
,
T. S.
,
Park
,
S.-S.
,
Yoo
,
H. S.
,
Kim
,
J. W.
, and
Kim
,
E. S.
, “
Development of multi-GPU-based smoothed particle hydrodynamics code for nuclear thermal hydraulics and safety: Potential and challenges
,”
Front. Energy Res.
8
,
86
(
2020
).
38.
Ren
,
W.
and
E
,
W.
, “
Boundary conditions for the moving contact line problem
,”
Phys. Fluids
19
(
2
),
022101
(
2007
).
39.
Rowlinson
,
J. S.
and
Widom
,
B.
,
Molecular Theory of Capillarity
(
Courier Corporation
,
Mineola
,
2002
).
40.
Shigorina
,
E.
,
Kordilla
,
J.
, and
Tartakovsky
,
A. M.
, “
Smoothed particle hydrodynamics study of the roughness effect on contact angle and droplet flow
,”
Phys. Rev. E
96
(
3
),
033115
(
2017
).
41.
Skillen
,
A.
,
Lind
,
S.
,
Stansby
,
P. K.
, and
Rogers
,
B. D.
, “
Incompressible smoothed particle hydrodynamics (SPH) with reduced temporal noise and generalised Fickian smoothing applied to body–water slam and efficient wave–body interaction
,”
Comput. Methods Appl. Mech. Eng.
265
,
163
173
(
2013
).
42.
Sun
,
P.
,
Zhang
,
A.-M.
,
Marrone
,
S.
, and
Ming
,
F.
, “
An accurate and efficient SPH modeling of the water entry of circular cylinders
,”
Appl. Ocean Res.
72
,
60
75
(
2018
).
43.
Tartakovsky
,
A.
and
Meakin
,
P.
, “
Modeling of surface tension and contact angles with smoothed particle hydrodynamics
,”
Phys. Rev. E
72
(
2
),
026301
(
2005
).
44.
Tartakovsky
,
A. M.
and
Panchenko
,
A.
, “
Pairwise force smoothed particle hydrodynamics model for multiphase flow: Surface tension and contact line dynamics
,”
J. Comput. Phys.
305
,
1119
1146
(
2016
).
45.
Techet
,
A. H.
and
Truscott
,
T. T.
, “
Water entry of spinning hydrophobic and hydrophilic spheres
,”
J. Fluids Struct.
27
(
5–6
),
716
726
(
2011
).
46.
Thoroddsen
,
S. T.
,
Etoh
,
T. G.
,
Takehara
,
K.
, and
Takano
,
Y.
, “
Impact jetting by a solid sphere
,”
J. Fluid Mech.
499
,
139
148
(
2004
).
47.
Truscott
,
T. T.
,
Epps
,
B. P.
, and
Belden
,
J.
, “
Water entry of projectiles
,”
Annu. Rev. Fluid Mech.
46
,
355
378
(
2014
).
48.
Truscott
,
T. T.
and
Techet
,
A. H.
, “
A spin on cavity formation during water entry of hydrophobic and hydrophilic spheres
,”
Phys. Fluids
21
(
12
),
121703
(
2009
).
49.
Wei
,
Z.
and
Hu
,
C.
, “
An experimental study on water entry of horizontal cylinders
,”
J. Mar. Sci. Technol.
19
(
3
),
338
350
(
2014
).
50.
Wei
,
Z.
and
Hu
,
C.
, “
Experimental study on water entry of circular cylinders with inclined angles
,”
J. Mar. Sci. Technol.
20
(
4
),
722
738
(
2015
).
51.
Worthington
,
A. M.
and
Cole
,
R. S.
, “
V. Impact with a liquid surface, studied by the aid of instantaneous photography
,”
Philos. Trans. R. Soc. London, Ser. A
189
,
137
148
(
1897
).
52.
Xu
,
G. D.
,
Duan
,
W. Y.
, and
Wu
,
G. X.
, “
Numerical simulation of oblique water entry of an asymmetrical wedge
,”
Ocean Eng.
35
(
16
),
1597
1603
(
2008
).
53.
Ye
,
T.
,
Pan
,
D.
,
Huang
,
C.
, and
Liu
,
M.
, “
Smoothed particle hydrodynamics (SPH) for complex fluid flows: Recent developments in methodology and applications
,”
Phys. Fluids
31
(
1
),
011301
(
2019
).
54.
Zhang
,
A.-M.
,
Sun
,
P.-N.
,
Ming
,
F.-R.
, and
Colagrossi
,
A.
, “
Smoothed particle hydrodynamics and its applications in fluid-structure interactions
,”
J. Hydrodyn.
29
(
2
),
187
216
(
2017
).

Supplementary Material