Numerical investigation of air mediated droplet bouncing on flat surfaces

A liquid droplet can bounce off a flat substrate independent of surface wettability if the impact occurs at low velocities, i.e., We of less than seven. In this case, the droplet spreads on a sub-micrometer air layer and rebounds subsequently without any direct contact with the surface. We have numerically investigated the process of air layer formation beneath the droplet. The numerical simulations are validated using experimental results available in the literature based on morphology of the droplet interface and thickness of the air layer. Numerical results revealed that the formation of a high pressure zone at the center of impact deforms the droplet to a kink shape at the moment of impact. The deformation leads to displacement of high pressure zone from center to kink edge of the droplet interface. Further investigation of pressure and velocity of air beneath the droplet divulged that high pressure region at the kink edge suppresses air flow at the inner region while accelerating flow at the outer region. In addition, it is demonstrated that fluid flow at the kink edge where droplet interface has the minimum distance from the substrate resembles Couette flow. It is demonstrated that the deformation of droplet along with displacement of high pressure region from the center to kink edge are responsible for stabilizing the air layer beneath the droplet and consequently spreading and receding of droplet over a thin air cushion. © 2017 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/). [http://dx.doi.org/10.1063/1.4993837]


I. INTRODUCTION
Droplet impingement on wetted or dry surfaces is ubiquitous in industrial and natural processes which attracts attention of researchers during the last century.The phenomenon encompasses various physical mechanisms which are responsible for the wide variety of outcomes.The history of systemic investigation of the droplet impact goes to the nineteenth century when Worthington 1 studied droplet impact various surfaces.Over the last decades numerous experimental, numerical, and analytical studies have been performed to investigate unknown aspects of the phenomenon.
Considering normal impingement of a spherical droplet on a flat surface with initial diameter of D 0 and impact velocity of U 0 , the following non-dimensional groups are governing the phenomenon, where ρ l , µ l , and σ are representative of the liquid density, viscosity, and surface tension, respectively.We, Re, Oh, and Bo denote the Weber, Reynolds, Ohnesorge, and Bond numbers respectively, and τ indicates dimensionless time.
Despite the large body of literature, this phenomenon is far away from being thoroughly understood, and many of its fascinating aspects are being reported by researchers.One of the interesting applications, reported by Couder et al., 2 is air mediated droplet bouncing on a vibrating bath of the same fluid.Guided and piloted by its own wave field, the droplet exhibits behavior previously thought to be exclusive to the world of quantum mechanics (Bush 3 ).This phenomenon, along with others such as multiple bouncing of droplet on a bath before undergoing partial coalescence (Charles and Mason, 4 Chen, Mandre, and Feng, 5 Bach, Koch, and Gopinath 6 ), and the droplet bouncing a soap film (Gilet and Bush 7 ) suggest remarkable effect of surrounding gas on behavior of the droplet.
In the case of droplet impact on solid surfaces, formation of bubble at the center of droplet is reported by many researchers which are comprehensively reviewed by Josserand and Thoroddsen. 8handra and Avedisian 9 reported that as the droplet approaches the substrate, pressure rise at the center of impact forms a dimple at the bottom of the droplet which leads to air entrapment and bubble formation.Similar air entrainment in the case of droplet impact on water pool is also reported by Pumphrey and Elmore, 10 who performed experimental studies and postulated some plausible explanations of the phenomenon.Numerical investigation of droplet impingement on a solid surface, performed by Mehdi-Nejad, Mostaghimi, and Chandra, 11 provided some details about the pressure in the gas phase and deformation of droplet before wetting the surface.Their numerical results also confirmed bubble entrainment at the center of droplet.
Possibility of droplet spreading on thin air layer without any direct contact between droplet and solid surface is postulated for the first time by Mandre, Mani, and Brenner. 12This hypothesis, further examined by Kolinski, Mahadevan, and Rubinstein 13 who demonstrated that droplets rebound from superhydrophilic surfaces without ever touching the solid, because of presence of nanometer-scale air film layer.Detailed experimental investigations of this phenomenon, carried out by de Ruiter et al., 14 revealed that droplets with low impact velocities will bounce off the surface without any contact with solid.They used dual wavelength reflection interference microcopy method (de Ruiter, Mugele, and van den Ende 15 ) to measure thickness of the air layer and acquired detailed information about behavior of air underneath the droplet during the entire bouncing process.Based on the observations provided by de Ruiter, van den Ende, and Mugele, 16 a kink is formed at the edge of the droplet which is spreading outward while the center of the dimple remains stationary.For higher impact velocities, film thickness decreases until liquid-solid contact is formed and droplet wets the surface based on wettability characteristics of the surface.
The aim of the present study is to investigate the formation of air cushion and subsequent droplet bouncing on a flat substrate regardless of the surface wettability.This study provides a detailed quantitative investigation of pressure and velocity fields in the thin air layer underneath the droplet using numerical simulations.We claim that simultaneous droplet deformation and pressure pick displacement beneath the droplet are responsible for stabilization of air layer.Moreover, it is demonstrated that droplet deformation leads to formation of three distinct regions which are characterized by air flow behavior.

II. MATHEMATICAL FORMULATION
The flow is assumed to be laminar, incompressible, and axisymmetric.Besides, since fluid flow is isothermal, density and viscosity of each phase is considered to be constant.Momentum and continuity equations for incompressible laminar flows presenting both phases in the whole domain can be written as follows, where P, ρ, and µ, represent pressure, density and dynamic viscosity of the fluid.The last term in momentum equation, f σ , corresponds to the surface tension effect.
In order to take into account properties of liquid and gas phases, a so-called volume of fluid method (Hirt and Nichols 17 ) has been used.In this method two phases are considered as one effective fluid throughout the whole domain.Properties of the effective fluid are evaluated using weighted average of each phase properties.For instance, density and viscosity of the effective fluid are calculated as follows, where α is volume fraction function and is defined as It can be shown that the transport equation of the indicator function, α, can be written as following, Since surface tension is applied only at the interface of liquid and gas, f σ should be expressed as a proper relation which is non-zero at transitional area and zero elsewhere in the domain.A common approach for dealing with this situation is the Continuum Surface Force (CSF) model proposed by Brackbill, Kothe, and Zemach 18 which expresses curvature of the interface as follows, where α is the volume fraction function.The term for surface tension could be written as following, where σ is surface tension of the liquid.

III. COMPUTATIONAL PROCEDURE
Small spatial and temporal scales encountered in air cushioning phenomenon makes it necessary to use appropriate numerical methods for solving governing equations.In this study Gerris Flow Solver (Popinet 19 ) which is an open-source code for solving multiphase flows has been used.The advantage of this code lies in the utilization of quadtree hierarchical data structure for spatial FIG. 1.A Schematic of computational domain.
discretization which makes it possible to achieve high resolutions at the interface.Besides, a generalized piecewise linear interface calculation scheme (Scardovelli and Zaleski 20 ) has been used to prevent interface smearing due to the numerical diffusion.Moreover, in order to increase accuracy of curvature calculations, height function method (Torrey et al. 21) is utilized.In this study projection method of Chorin 22 is utilized to solve Navier-Stokes equations.Besides, discretization in time is performed using a second order staggered discretization of volume fraction and pressure (Popinet 23 ).For stability considerations the time step is adapted in each time step to satisfy CFL stability conditions (CFL number is chosen to be 0.5).
As illustrated in Fig. 1, simulations were carried out on an axisymmetric domain with a size of 2D 0 in both horizontal and vertical directions.Neumann boundary condition is utilized for top and right boundaries of the domain.At the bottom of the domain no-slip boundary condition is applied for velocity, and in order to make sure that droplet bouncing is not a result of hydrophobicity of the substrate, a very low contact angle (θ = 9 o ) is used for bottom edge to mimic a super-hydrophilic surface, as used in the experiments of de Ruiter et al. 14 .
An extremely thin air layer trapped underneath the droplet makes it necessary to have enough refined grids in the vicinity of the substrate (Fig. 2).Various heights for the refined region has been examined in order to ensure independency of simulations from size of the refined region which lead to the selection of h = 0.01D 0 , because further increment of h did not affect the results.Furthermore, grid independency test is performed to verify that simulations are not dependent on grid sizes (Fig. 3).At least 11 levels of grid refinement is required, because coarser grids are unable to capture the air layer beneath the droplet.Droplet profiles for various levels of refinement (Fig. 3) reveals that there is no significant change at the interface of droplet by further refining the grids.Consequently, all simulations in the rest of the paper are performed using 14 level of grid refinement.

IV. RESULTS AND DISCUSSIONS
Numerical results are used in this section to investigate air layer trapped underneath the droplet which leads to droplet bouncing on flat surfaces.In the first step, numerical simulations are validated against experimental results reported by de Ruiter et al. 14  predicted interface and presence of air layer beneath the droplet during spreading and receding phases of impact confirms that there is no direct contact between droplet and the surface.

A. Mechanism of air layer formation
Non-wetting regime of impact which is characterized by super hydrophobic like behavior of droplet impact occurs at lower impact We numbers (We < 7.0).In this regime, droplet takes sombrerolike dimple shape at initial stages of impact and spreads on a persisting air layer underneath the droplet (Fig. 6).During droplet spreading, the stationary trapped air volume increases until the droplet reaches to its maximum spreading diameter.During the retraction phase where the kink at the edge of droplet moves inward, the volume of entrained air decreases until the droplet detaches from the surface.It is important to mention that unlike the cases such as bouncing of droplet on vibrated bath or Leidenfrost effect of droplets over heated surfaces, there is no external force to assist the thin air layer formation.
Formation of the air film is initiated by creation of a high pressure region at the center of impact synchronized with a favorable deformation of droplet interface is responsible for persistent air layer formation at the bottom of droplet.Consequently the droplet deforms to take a kink shape and pushes air to the lower pressure regions.Due to the droplet deformation, high pressure zone moves towards the kink edge (Fig. 7), where the droplet interface has minimum distance from the substrate.The displacement of high pressure zone pushes out air flow (Fig. 8) and creates a stationary air bucket near the center of impact.The process of pushing air from the high pressure zone and subsequent air squeezing, continues until droplet reaches to its maximum spreading diameter.As can be seen in Fig. 7  momentum and capillary forces beneath the droplet.In numerical simulations with low We numbers only one kink forms at the droplet interface.However, by increasing the We number and consequently increasing the fluid's inertial two or multiple kinks may occur at the droplet interface.
Maximum pressure at the impact center increases by increasing the droplet impact velocity (We > 7.0).Higher degrees of droplet deformation, which is a result of pressure increase, reduces the minimum distance of droplet interface at the kink edge.If the minimum distance of the droplet profile become less than the molecular mean free path, air layer collapses and the droplet wets the substrate as shown in Fig. 9.The wetting regime which occurs in most of droplet impacts is characterized by direct contact of droplet and the substrate.In this regime, behavior of droplet during advancing and receding depends on substrate characteristics.Deformation of droplet and the air pocket entrained at the center of the droplet forms a bubble which is repeatedly reported in the literature (e.g.Chandra and Avedisian, 9 Mehdi-Nejad, Mostaghimi, and Chandra 11 ).

B. Characteristics of the air layer beneath the droplet
Numerical simulation of air cushioning phenomenon resolving the thin air layer beneath the droplet makes it possible to investigate the characteristics of air flow using velocity and pressure fields.Droplet interface during spreading phase for We = 1.5 and Oh = 0.0028 is illustrated in Fig. 10, which also includes pressure and velocity plots along a horizontal line (red line) in the air AIP Advances 7, 095003 (2017) FIG. 9. Droplet-air interface at early stages of impact, occurred in wetting regime (We = 9.0, Oh = 0.0028).layer beneath the droplet.The maximum pressure occurs at the kink edge where the air layer has minimum thickness.Higher pressure at this location, causes different flow scenarios at the left and right sides of the kink.The left side, air flow is being suppressed.On the other hand, the pressure difference produces high velocity air flow at the right side of the kink.As can be seen in Fig. 10 fluid flow in the air layer starts to accelerate near the high pressure zone which is located at the position where droplet interface has minimum distance from the solid surface.The location of maximum pressure (i.e. kink) varies as the droplet spreads.
In order to further examine the thin air layer flow behavior, the air velocity profile along a vertical line at a fixed location at r/D 0 = 0.06 is plotted in Fig. 11.The presented times are carefully chosen to demonstrate the moving kink before approaching and consequently passing the vertical line.Velocity profile along the vertical line over which the pressure pick is moving demonstrates a behavior similar to Couette flow with positive, negative and without pressure gradient.
The vertical red line at time τ = 0.016 is placed at the right side of pressure pick.The velocity profile represents a Couette flow profile with negative pressure gradient (pressure decreases in the direction of upper motion).In this situation, velocity of the droplet interface, and pressure difference caused by creation of the high pressure zone at the kink edge are responsible for inducing velocity in the gas field.When the kink edge is situated at r/D 0 = 0.6, velocity profile associates with Couette flow in the absence of pressure gradient which means that moving droplet interface is the only driving force inducing velocity in the gas phase.By further movement of kink edge to a position where the kink and consequently the pressure pick pass the location of r/D 0 = 0.6, the pressure gradient becomes positive (pressure increases in the direction of moving interface).
Having considered various behaviors of fluid flow, the gas layer beneath the droplet can be divided into three distinct regions namely inner, kink, and outer regions (Fig. 12).In the inner region, positive pressure gradient either suppresses or creates reverse velocity field.As a result, stationary FIG.12.A schematic of different regions in thin air layer.FIG. 13.Pressure and velocity fields along vertical lines situated in inner, kink and outer regions (We = 1.5, Oh = 0.0028, τ = 0.017).gas entrained at the center of droplet is mainly a consequence of high pressure at edge of the kink.In the kink region which is situated under the kink edge of the droplet, there is a Couette flow without pressure gradient and a linear velocity field occurs at this region.At the outer region, pressure gradient is the main driving force for induction velocity in the gas layer.
Velocity profile in the gas field is plotted along vertical lines situated at different regions in the air layer in Fig. 13.It can be seen that velocity at the outer region is mainly induced by pressure difference between maximum pressure at the kink region and the outer region.Magnitude of velocities in this area is much larger than those in the other regions because of greater pressure differences.Linear velocity profile occurs at the kink edge where there is no pressure gradient.Induced reverse flow and suppression of velocity in the inner region are consequences of pressure difference between kink and inner regions.

V. CONCLUSIONS
Air mediated bouncing of droplet independent of the substrate wettability is investigated using numerical simulations.A quadtree based adaptive mesh refinement is used in order to have extremely fine grid sizes at the interface and bottom of the droplet.Simulation results revealed that there is a persistent air layer underneath the droplet during spreading and receding phases and droplet does not have any direct contact with the solid.Pressure rise at the center of impact deforms the droplet to take a dimple shape and consequently high pressure zone shifts from center to the kink edge of the droplet interface.
Impacting droplet may either spread on a air layer or wet the surface depending on the impact We number.In the non-wetting regime which occurs at We numbers less than 7.0, droplet spreads on a thin film of air without direct contact with the solid.In this regime droplet spreads, recedes and rebounds.On the other hand, in the wetting regime the distance of the kink shape decreases to values less than molecular mean free path and droplet comes into contact with the solid surface.Bubble formation at the center of impinging droplet is a consequence of the deformation of droplet and wetting regime.
During droplet spreading on air layer in bouncing regime, air layer beneath the droplet can be divided into three distinct regions namely inner, kink, and outer regions.In the outer region the pressure difference is the main driving force inducing velocity in the gas, while in the kink region movement of droplet interface is the main driving force inducing velocity in the gas layer.The inner region is characterized by suppression of velocity field which is a result of high pressure zone at the kink edge.

FIG. 10 .
FIG. 10.Illustration of air layer behavior by presenting a) droplet-air interface, b) pressure field, and c) velocity field along horizontal red line (We = 1.5, Oh = 0.0028).

FIG. 11 .
FIG. 11.Illustration of air layer behavior by presenting a) droplet-air interface and b) velocity field along a vertical red line at r/D 0 = 0.6 (We = 1.5, Oh = 0.0028).