Numerical Investigation of Water Entry of Half Hydrophilic and Half Hydrophobic Spheres
A numerical simulation to investigate the water entry of half-half sphere which is hydrophobic on one hemisphere and hydrophilic on the other is performed. Particular attention is given to the simulation method based on solving the Navier-Stokes equations coupled with VOF (volume of fluid) method and CSF (continuum surface force) method. Numerical results predicted experimental results, validating the suitability of the numerical approach to simulate the water entry problem of sphere under different wetting conditions. Numerical results show that the water entry of the half-half sphere creates an asymmetric cavity and “cardioid” splash, causing the sphere to travel laterally from the hydrophobic side to the hydrophilic side. Further investigations show that the density ratio and mismatch of asymmetric in wetting condition affect the trajectory, velocity, and acceleration of the half-half sphere during water entry. In addition, the total hydrodynamic force coefficient is investigated as a result of the forces acting on the sphere during water entry dictated by the cavity formation.
The water entry problem of a solid object impacting on a liquid surface has challenged researchers for centuries. This study is directly relevant to many different applications, such as ship slamming , skipping stones , water walking insects , and landing of seaplane . The water entry of sphere, the most canonical object, presents many unforeseen mechanisms which are important for researchers to study the physics during water entry.
The first systematic study of water entry problem, such as splashes and cavity formation, was published more than a century ago by Worthington , who used single-flash photography to investigate the process of solid objects impacting on liquid surfaces. He observed that the water entry of spheres coated with soot formed a distinct subsurface cavity whereas an already wetted sphere did not.
May  also did research on the effect of surface condition of a sphere on its water entry cavity. He described that “At 30 ft/sec it was found that a cavity was always produced but that the cavity was generally smaller when the sphere was clean.”
In Worthington’s study, the sphere is described as “wetted sphere” and “sphere coated with soot”; in May’s study, the sphere is described as “clean sphere” and “handled sphere.” Essentially, “wetted sphere” and “clean sphere” in their studies, to put it more precise, are hydrophilic sphere; “sphere coated with soot” and “handled sphere” are hydrophobic sphere.
Generally, the surface condition of the sphere, for example, static wetting angle, can dramatically affect the water entry phenomenon. Hydrophobic sphere (such as the sphere coated with soot) typically forms cavities, whereas hydrophilic sphere (such as an already wetted sphere) does not. Duez et al.  presented a theoretical model, arguing that the condition to create a water entry cavity is that the impact velocity must be above a critical magnitude, and the critical velocity is found to be dependent on the sphere’s static contact angle.
The theoretical work of water entry problem can be traced back to von Kàrmàn  and Wagner , who investigated the phenomenon related to the impact of objects on liquid surfaces. They considered the solid body entering a fluid at large Reynolds and Weber numbers, where viscosity and surface tension were neglected. In the recent years, there has been a growing interest in the physics of impact, particularly due to the development of computer. Greenhow and Moyo  simulated the free surface deformations of initially calm water caused by the impact of a horizontal circular cylinder by using a nonlinear two-dimensional numerical calculation. Lin  developed a two-dimensional numerical model to simulate a moving body in a free surface flow based on the cut-cell technique in a fixed-grid system, and the volume of fluid (VOF) method was applied to track interface. Villanueva and Amberg  investigated the wetting phenomenon of the solid objects in liquid sintering-like flows by numerical simulation. The wetting phenomenon is modeled using the coupled Cahn-Hilliard/Navier-Stokes system, which is called the phase-field method. Do-Quang and Amberg [12, 13] used the same method to simulate the impact of a small solid sphere on a liquid surface, where surface wettability and liquid properties are considered. Recently, Abraham et al.  modeled the forces acting on the sphere during the early stage of water entry by coupling Navier-Stokes equations with the Sheer Stress Transport (SST) model in the software ANSYS CFX.
In addition to surface static wetting angle, transverse spin can affect the trajectory and cavity formation of a sphere after water entry. Truscott and Techet [15–17] showed that transverse spin imparts mismatch surface velocity on either side of the sphere and finally led to asymmetrical air cavity formation around the sphere, which caused curvature to the trajectory of a sphere during water entry.
The ultimate result presented herein is that the asymmetrical cavity formation can also be produced by simply dropping a half-half sphere that is hydrophobic on the hemisphere side and hydrophilic on the other. This phenomenon is the main focus of this paper. In this paper, the water entry of a half-half sphere is investigated by numerical simulation. The numerical and experimental results are compared to validate the numerical method. Then several cases of the water entry by the half-half sphere are illustrated. The trajectories, velocities, accelerations of the sphere, and the forces acting on it during water entry are also investigated.
2. Numerical Model
2.1. Governing Equations
In this study, the volume of fluid (VOF) method is used to track the volume fraction of each of the fluids throughout the domain. This method introduces a variable for each phase: the volume fraction of the phase in the computational cell. If the fluid’s volume fraction in the cell is denoted as , then we have the following:the cell is empty of the fluid; the cell is full of the fluid; the cell contains the interface between the fluid and other fluids.
Based on the local value of , the appropriate properties and variables will be assigned to each control volume within the domain. In each control volume, the volume fractions of all phases sum to unity. In this article, the fluid flow is described by the Navier-Stokes equations for a three-dimensional incompressible two-phase flow. The two phases are water and air, respectively. The water-air interface is defined as .
The governing equations for fluid flow are expressed by (1) and (2). These equations represent conservation of momentum and mass, respectively.where is the density; is the velocity vector; is the dynamic viscosity; is the pressure; is the gravity acceleration; is the surface tension force density.
The density and viscosity appearing in the momentum equations are determined by the presence of the component phases in each control volume. For a two-phase system in this study, the density and viscosity are given by (3) and (4), respectively.where is the volume fraction of liquid phase; are the density and viscosity of the liquid and gas phase, respectively. The detailed magnitude of phase properties is shown in Table 1.
The continuum surface force (CSF) model proposed by Brackbill et al.  has been implemented so that the addition of surface tension to the VOF calculation results in a source term in the momentum equation.
Let be the surface normal, defined as the gradient of , the volume fraction of the phase.
The curvature, , is defined in terms of the divergence of the unit normal, :where
The surface tension can be written in terms of the pressure jump across the surface. The force at the surface can be expressed as a volume force using the divergence theorem. It is this volume force that is the source term which is added to the momentum equation. It has the following form:
In this study, only two phases are presented in a cell, then , , and the surface tension term can be written aswhere is the volume-averaged density computed using (3).
Many numerical researchers intended to simulate wettability, such as Do-Quang’s work [12, 13]; the water entry problem is modeled using the coupled Cahn-Hilliard/Navier-Stokes system, where surface contact angle is used to correct the unit vector normal to the solid surface. Another model deals with the wettability problem presented by Brackbill et al. , called CSF model. This method is more applicable. Ashish Saha and Mitra  used this method to deal with a microfluidic capillary flow. Mirzaii and Passandideh-Fard  used this method in his research to simulate the water entry of sphere.
An option to specify a static contact angle is applied according to the work by Brackbill et al. . The contact angle is used to adjust the surface normal in cells adjacent to the wall.where and are the unit vectors normal and tangential to the wall, respectively.
2.2. Computational Domain
Setting up an appropriate computational domain is a key step in correctly simulating the multiphase problem. In this paper, the CFD computation deals with a water entry problem of a single solid sphere with different wetting properties and kinematics impacting on a free liquid surface, so a 3D computational domain is built. The schematic diagram of the solution domain is illustrated in Figure 1. Spheres are dropped from a given distance above the free surface, and impact occurs at with an approximate impacting speed . For the convenience of comparing with experiment data, the sphere diameter chosen in this study is and the mass ratio of the sphere is , respectively.
Atmospheric pressure, turbulence quantities, and air mass fraction are specified at the top inlet boundary. The right and left boundaries are set up as wall boundary condition. The bottom of the computation domain is the pressure outlet boundary condition, and the static pressure and turbulence quantities are specified. At the surface of the solid sphere, the no-slip wall boundary condition is used. Numerically, the approach to implementing slip depends on the method used to track the interface. VOF methods utilize cell face normal velocities to advect volume fractions . This implies that the methodology includes an “implicit” slip length at no-slip boundary. This condition is fit to many multiphase flows, such as Do-Quang and Amberg [12, 13] and Ashish Saha and Mitra .
To reduce the computational cost, the symmetry-plane boundary condition is used in the simulation. As illustrated in Figure 1(a), the core of the sphere is located at z-x plane. Because the sphere only moves at the z-x plane, there is no displacement in the y-direction. Therefore, the symmetry plane is located just at z-x plane through the sphere core. The computational domain and boundary at the symmetry plane are illustrated in Figure 1(b), the free liquid surface is 20D below the inlet boundary, the sphere is located 10D above the free surface, the coordinate zero point is located in the free liquid surface, and the positive z-direction is vertically downward.
The grid was generated using the commercial software ANSYS ICEM as a preprocessor for mesh generation. The mesh employed could be structured, unstructured, or hybrid. Structured mesh is identified by regular connectivity, while unstructured mesh is irregular connectivity. A hybrid mesh contains a mixture of structured portions and unstructured portions. Compared to the unstructured mesh, structured mesh has better convergence and higher resolution, but much more work in preprocess. In the present research, the structured mesh is employed throughout the whole domain, and O-grid is built around the sphere. An O-grid with in dimension is built around the sphere for clustering, as shown in Figure 2. The mesh law for clustering is “BiGeometric” with ratio 1.1.
A mesh refinement study was performed in which the grid size was gradually increased: coarse mesh, fine mesh, and finer mesh. The coarse mesh case contained 455670 elements. In contrast, the fine mesh case consists of 987600 elements. A further refinement was carried out with a mesh of 1538340 elements. This three-level refinement was applied in the hydrophobic cases (Figures 3 and 4) and also in the hydrophilic cases (Figures 5 and 6). The cavity shape and the sphere depth are compared, and the finer mesh with 1538340 elements was found to be applicable. Calculations were performed on a Linux PC-cluster with 20 processors. The time step convergence was monitored and the simulation was considered to have converged when residuals of all conserved variables fell below . The present simulations required an average of about 30 subiterations to make the solution converge at each time step. The computation took about 48 hours for each water entry with the finer mesh.
(a) Experimental results 
(b) Coarse mesh
(c) Fine mesh
(d) Finer mesh
(a) Experimental results 
(b) Coarse mesh
(c) Fine mesh
(d) Finer mesh
2.3. Numerical Treatment
The numerical simulations were carried out using ANSYS FLUENT software, a finite volume code in dealing with multiphase flow problem. The boundary conditions, initial conditions, user defined functions, and the method of solving each equation were all specified in a user interface. The time-dependent governing equations were discretized by using the finite volume method; the SIMPLE algorithm was used for coupling between pressure and velocity with a green-gauss cell based gradient evaluation option to solve Reynolds-Averaged Navier-Stokes (RANS) equations. All the solution variables were solved via second-order upwind discretization scheme for improving accuracy. The PRESTO (pressure staggering option) scheme was used for the pressure interpolation. Standard k-ε model  was carried out in the software ANSYS FLUENT.
3. Results and Discussions
3.1. Numerical Validation
In order to validate the numerical methodology in the present paper, simulations are made to directly compare with experiment results in Truscott and Techet . The simulation results are plotted together with the experimental results as shown in Figures 7 and 8. As illustrated in Figure 7, the sphere is hydrophilic with a static contact angle = 60°, and its impact speed is = 1.72. In this hydrophilic case, a simple vertical jet ejects from the free surface after impact but a subsurface air cavity does not form, and this numerical result agrees well with the experiment result.
Another simulation was also made to demonstrate the ability of modeling wetting phenomenon in this paper, as illustrated in Figure 8. The sphere is hydrophobic with a static contact angle = 120°, and its impact speed is = 1.72, the same as the case in Figure 7. The numerical results show a subsurface cavity following the sphere and subsequently pinches-off, at some point below the free surface. The cavity size and the pinch-off location from numerical simulation are well agreed with the experimental result.
Comparing the hydrophilic case with hydrophobic one, the spheres are identical in terms of diameter, and kinematic properties, however, are only different in surface wettability that exhibit very different water entry phenomenon: in the hydrophobic case, a classic splash crown and subsurface air cavity form; the air cavity grows and subsequently collapses at pinch-off, whereas no such behavior is observed in the hydrophilic case. Figures 7 and 8 illustrate the fundamental differences between the impact of hydrophilic spheres and hydrophobic spheres, and our numerical results predicted experimental results well validating the trustworthy of the numerical approach to simulate the water entry problem of sphere under considerations of wettability.
3.2. Flow Patterns
More detailed water entry photos of cavity formation are shown in Figure 9. The sphere mass ratio and impact speed = 2.37 m/s. Figure 9(a) shows a hydrophilic sphere ( = 60°) entering water, and no cavity forms, whereas Figure 9(b) shows a hydrophobic sphere ( = 120°) entering water, resulting in a cavity formation case. The essential difference of the water entry cavity between the hydrophilic and hydrophobic sphere is the same as described in the previous section: a big air cavity is created in the hydrophobic case, while no such behavior is observed in the hydrophilic case. If the sphere is a half-half sphere that has a hydrophilic surface on one side and hydrophobic surface on the other, what the water entry phenomenon should be? We contend that it is possible to form a half cavity on one side and no cavity on the other side at a moderate water entry velocity .
(a) Hydrophilic sphere
(b) Hydrophobic sphere
(c) Half-half sphere
Figure 9(c) shows the water entry phenomenon of this half-half sphere. The sphere is hydrophilic on the right half side with static contact angle = 60° and hydrophobic with static contact angle = 120° on the left half side. As illustrated in Figure 9(c), the asymmetric surface wetting condition makes an asymmetric water entry cavity and causes the sphere to travel laterally from the hydrophobic side (left half) to the hydrophilic side (right half).
The streamline and isobaric line around a sphere during water entry are shown in Figures 10 and 11. Sphere diameter and impact speed are all the same as the cases in Figure 9. Because of the symmetry of the flow field, the streamline and isobaric line are shown together in each image. In Figure 10(a), the sphere is hydrophilic, there is no cavity formation, and a simple jet is observed. A complicated flow field with several vortices around the jet in the gas phase can be observed. Figure 10(b) shows the streamline and isobaric line around a hydrophobic sphere during water entry. There is a big vortex appearing in the gas phase above the free surface. Figure 11 shows the streamline and isobaric line around a half-half sphere. The unbalanced momentum associated with the asymmetric cavity formation causes the sphere to travel laterally, and this asymmetric cavity makes the flow field more complex.
(a) Hydrophilic sphere
(b) Hydrophobic sphere
(b) Pressure contour
The most exciting behavior of the half-half sphere during water entry is not only the asymmetric cavity but also the “cardioid” splash as shown in Figure 12. As the sphere impacts the free water surface, water is drown up the right hydrophilic side, creating a fluid wedge and no cavity formation; air is entrained on the left hydrophobic side, creating an air cavity. In the half-half case as illustrated in Figure 12, the initial splash radiates outward asymmetrically at first and begins to curl up along the equator as the sphere descends. A fluid wedge forms apparently from ms after impact. At time around ms, the wedge has almost completely traversed the cavity and is still affecting the splash asymmetry, and a cardioid shape is seen in the top view, with the wedge fully bisecting the cavity. More details of the fluid wedge are shown in Figure 13. The fluid wedge is drawn into the center of the cavity along the equator. The side view in Figure 13 shows the formation of asymmetric cavity in this early stage of water impact. The top view in Figure 13 shows a distinct “cardioid” splash and shows the fluid wedge bisecting the cavity from right side to the left.
It should be noted that, in this paper, the impact velocity is not large with the purpose of highlighting the difference of water entry cavity created by hydrophilic and hydrophobic sphere. When the velocity is too large, the effect of wettability becomes to be little, and no significant difference can be observed on the shape of cavity. So in most of the literatures [6, 12, 13, 15] and also in this paper, the impact velocity is always <5 m/s. When the velocity is small enough, the wettability can dramatically affect the water entry phenomenon, such as Figure 14.
(a) Hydrophilic sphere = 60°
(b) Hydrophobic sphere
(c) Half-half sphere = 120°/60°
In this simulation, a wetting liquid film is observed that develops immediately after impact. Korobkin and Pukhnachov  detailed the thin fluid film that develops during the early stage of water impact and climbs up along the impacting object. This film is evident in Figure 14(a). The film dynamics is seen to strongly differ depending on the wettability of the sphere surface. For the hydrophilic case considered in Figure 14(a), the film is seen to follow the sphere and closes up at the north pole of the sphere, and no cavity is created. On the contrary, for the hydrophobic case shown in Figure 14(b), the film is seen to detach from the sphere before reaching the pole. For the half-half case shown in Figure 14(c), the film is asymmetric in the right and left side of the sphere, and the splash asymmetry is already beginning in the very early stage of water impact. On the left hydrophobic side, the film detaches from the sphere at some place leading to cavity formation; on the right hydrophilic side, the film climbs upon the sphere and meets at the north pole of the sphere, resulting in a noncavity formation case. The film gathers in the north pole of the sphere and creates fluid wedge. This wedge continues to move towards the cavity formation side, bisecting the cavity into two distinct symmetric holes and consequently forms a “cardioid” shape splash.
The most obvious and exciting behavior of the half-half sphere during water entry is not only the asymmetric air cavity and “cardioid” splash but also the curvature in trajectory. The lateral force (x-direction) induced by the unbalance momentum associated with the asymmetric cavity formation moves the sphere from the hydrophobic side to the hydrophilic side. Sphere trajectories for three cases are plotted in Figure 15, for the same impact speed = 2.37 m/s and the same half-half surface condition = 120°/60° but with different density ratio . The z- and x-positons are normalized by the diameter of the sphere. The free water surface corresponds to z/d = 0. Since gravity plays an important role in the water entry of the sphere, the mass ratio of the sphere compared with water is considered.
As illustrated in Figure 15, the lighter sphere tends to have more curvature in its trajectory than the other two heavier spheres. For low density ratio, the inertial forces are diminished compared to hydrodynamic forces, and lateral forces become more significant. For high density ratio, such as , about the density ratio of steel, the inertial forces dominate and the half-half surface condition does not have much effect on the trajectory of the sphere.
The cavity formation for the three half-half spheres with different density ratios is shown in Figure 16. Figure 16 shows a distinct difference in cavity formation for the three spheres despite the identical impact parameters that density ratio has a significant effect on the cavity formation for the half-half sphere cases. The heaviest sphere, with density ratio , descends more quickly and does not have a distinct lateral translation. The lightest sphere, with density ratio , travels an apparent lateral displacement and descends slower than others.
The hollow triangle markers in Figure 15 indicate the position of the sphere when pinch-off occurs. The depth for the sphere with density ratio is deeper than the other two. This trajectory information does not reveal the exact quantity of the velocity and acceleration of the spheres, but qualitative velocity and acceleration differences can be seen in Figures 15 and 16. The lighter sphere experiences a much larger deceleration in both z- and x-direction compared to the other two and thus descends slower and has a more significant lateral travel. Further discussion of this is found in the following section.
In order to determine the hydrodynamic force acting on the half-half sphere during water entry, the displacement, velocity, and acceleration are got directly from simulation, as shown in Figure 17.
As illustrated in Figure 17, the lighter sphere with density ratio = 1.2 loses vertical speed more quickly and descends slower. The heavier sphere with density ratio = 7.8 does not lose speed but accelerates due to its large density ratio and descends more quickly.
In x-direction, the lighter sphere has more lateral travel with a larger lateral speed. The lateral component of the velocity increases with asymmetric cavity formation and is the lowest for the biggest density ratio = 7.8 spheres, as illustrated in Figure 17(d). For the three density ratio cases, = 2.7 and = 7.8 have the similar trend, the velocity of x-direction increases with asymmetric cavity formation. But for the lightest sphere = 1.2, after about t = 45 ms, Ux/ decreases.
From the vertical velocity illustrated in Figure 17(c), we can see that the sphere with = 1.2 loses vertical speed rapidly, and this will cause much smaller cavity formation. Because the lateral translation for the half-half cases is led by the asymmetric cavity formation, when the cavity decreases in volume, the velocity in lateral direction will decrease apparently. For the cases = 2.7 and = 7.8, the asymmetric cavity persists well, so the lateral speed increases steadily. It can be expected that when the cavity collapses thoroughly, the lateral velocity will come to zero finally.
The acceleration of the sphere, received directly from simulation, can be used to determine the forces acting on the sphere along its descent. The acceleration is directly proportional to the force. The z-component of acceleration (Figure 17(e)) shows that each case has a constant acceleration before the moment of impact. When impact occurs, acceleration peak against the direction of the gravity is seen in about = 0~5 ms. When pinch-off occurs (about ms), the acceleration peak appears again with the opposite direction, indicating that a pulse acting on the sphere appeared.
The lateral acceleration of the sphere (Figure 17(f)) shows initial acceleration to the positive x-direction for all cases. Each case has a short window of increasing acceleration initially, which then decreases towards cavity pinch-off. When pinch-off occurs, the lateral acceleration behaves as a significant peak, indicating a pressure pulsation appears in the cavity connected to the sphere.
Accelerations got directly from the numerical simulation were used to calculate the total hydrodynamic forces  acting on the sphere during water entry by Newton’s second law. In z-direction,where is the mass of the sphere, is the gravity acceleration, and is the total hydrodynamic force acting on the sphere that results from surface tension forces, viscous traction forces, and pressure. So stand for component of the total hydrodynamic force in z-direction. It is important to reiterate that we simply concern the total hydrodynamic force and not a drag force in tradition sense.
Regardless of the nature of the hydrodynamic force, the total hydrodynamic force can be inferred from the acceleration of the sphere using (11). Then the total hydrodynamic force coefficient can be calculated from the acceleration data.where is the instantaneous magnitude of the velocity of the sphere, is the density of the water, , and R is the radius of the sphere.
The coefficient of the hydrodynamic force in x-direction is the same as that in z-direction.where is the component of the total hydrodynamic force in x-direction.
The force coefficients and for cases with different density ratio are shown in Figure 18 comparing the three density ratios for the same half-half surface condition: = 120°/60°. As illustrated in Figure 18(a), curves are similar up to t≈10 ms, and the pinch-off time for these cases is almost the same at about t≈45 ms, revealing that the mass ratio does not affect the pinch-off time. The impact forces occurring immediately after impact are independent of mass ratio and surface condition, but only relevant to the impact speed. When pinch-off occurs, as illustrated in Figure 18, the force coefficient curve shows a positive peak revealing that a pulse acting on the sphere appeared resulting from cavity pinch-off.
The numerical simulation results reveal that, by a change in the surface wetting conditions, it is possible to change the water entry behavior when a sphere impacts onto a free liquid surface. In addition, the numerical simulations of dynamic wetting during water entry are scarce in literature. The numerical results show good agreement with the experimental results. Therefore, a trustworthy simulation result is demonstrated even for complex wetting cases, such as the half hydrophilic and half hydrophobic sphere.
This water entry phenomenon of the half-half sphere shows an asymmetric water entry cavity and a “cardioid” splash, thus resulting in a lateral travel from the hydrophobic side to the hydrophilic side. The trajectory, velocity, and acceleration are investigated by several cases with different sphere density ratios. Results show that the lateral travel is affected by the density ratio. The lightest sphere has the most lateral travel whereas the heaviest sphere descends more quickly. Future additional simulations should be done to directly match experimental results and the analysis.
Overall, the water entry of half-half sphere with mismatch in surface wetting condition gives a new perspective on water entry problem by pointing to the unexpected role of surface wettability. This has applications to many hydrodynamics problems. Moreover, consideration should be taken into account when trajectory, hydrodynamic force, and cavity shape are important factors in engineering and military applications.
The authors declare that there is no conflict of interests regarding the publication of this paper.
This project is supported by the Natural Science Foundation of Heilongjiang Province, China (Grant no. A201409), the Special Foundation for Harbin Science and Technology Innovation Talents of China (Grant no. 2013RFLXJ007), and the Fundamental Research Funds for the Central Universities, China (Grant no. HIT. NSRIF. 201159).
M. Faltinsen, Sea Loads on Ships and Offshore Structures, Cambridge University Press, Cambridge, UK, 1990.
T. von Kàrmàn, “The impact on seaplane floats during landing,” NACA Technical Report 321, 1929.View at: Google Scholar
A. M. Worthington, A Study of Splashes, Longmans Green and Company, New York, NY, USA, 1908.
H. Wagner, “Phenomena associated with impacts and sliding on liquid surfaces,” Zeitschrift für Angewandte Mathematik und Mechanik, vol. 12, p. 193, 1932.View at: Google Scholar
B. E. Launder and D. B. Spalding, Lectures in Mathematical Models of Turbulence, Academic Press, London, UK, 1972.