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.

## I. INTRODUCTION

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* = Δ*ρgL*^{2}/*σ*, 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 $\theta 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.

## II. NUMERICAL MODELING

### A. Fundamentals of SPH

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)

where $r\u21c0$ denotes the point vector, $fri\u21c0$ is a function at the position of center particle *i*, and subscript *j* represents the neighboring particles of the center particle *i*. *V*_{j}, *m*_{j}, and *ρ*_{j} are the volume, mass, and density of the *j*th particle, respectively. $Wij|ri\u21c0\u2212rj\u21c0|,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).

### B. Governing equations

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),

where t is time and *P*, *ν*, *g*, and *f*_{s} 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

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),

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),

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)

where *c*_{0} 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.

### C. Surface tension

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),

where *σ* is the surface tension and *κ*_{i} and $n\u21c0i$ 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),

where *σ*_{αβ} is the surface tension between phases *α* and *β*, $T\tau \u21c0$ and $Tn\u21c0$ 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),

where *s*_{ij} 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 *s*_{ij}. 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),

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).

### D. Fluid solid interaction

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

where $v\u21c0$, $v\u21c00$, $w\u21c0$, and $w\u21c00$ are the velocity, initial velocity, angular velocity, and initial angular velocity, *f*_{i} denotes a total force acting on the *i*th solid particle from neighboring fluid particles, and $I\u2194$ and $\tau \u21c0$ are the momentum of inertia tensor and torque vector, respectively.

## III. NUMERICAL INVESTIGATION

### A. Validations of wetting characteristics modeling

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

where *R*_{1} and *R*_{2} 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.

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 (R_{0} = 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.

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.

### B. Numerical model

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/m^{3} and 0.005 Pa s, respectively. The gas density and the viscosity are set to be 1 kg/m^{3} 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.

## IV. RESULTS AND DISCUSSION

### A. Water entry of a hydrophilic sphere and a hydrophobic sphere

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).

### B. Effect of wettability on the contact-line dynamics

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).

### C. Effect of wettability and capillary number on the impact behavior

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.

### D. Water entry of a sphere with anisotropic wettability

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).

## V. CONCLUSION

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.

## SUPPLEMENTARY MATERIAL

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

## ACKNOWLEDGMENTS

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

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

## DATA AVAILABILITY

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