Dynamics of a dry-rebounding drop: observations, simulations, and modeling

3.1 Finite volume method and volume of fluid method

The drop was simulated using a two-phase solver, interFoam [19] of OpenFOAM (version 3.0.x) [15, 20], a CFD toolkit software that can be used and exploited under the GNU General Public License (GPL) [16]. The interFoam solver is based on the finite volume method (FVM). In the FVM, the domain of calculation is divided into finite-volume cells which are referred to as control volumes, and physical values (e.g., velocity and pressure) are assigned to the centroid or faces of each cell.

A type of VOF method is used in interFoam to model two-phase flow and to track the free surface. In the present simulation, the two phases of water and air were considered. Note that in the two-phase flow, the volume fraction of liquid, α_l = α, determines the volume fraction of gas, α_g = 1 – α.

3.2 Interface capturing

An efficient method is required to simulate a multiphase flow and capture a sharp interface between the two immiscible phases. VOF methods have a problem with respect to the diffusive interface between two phases. In VOF methods, the volume fraction of each phase is tracked through every control volume. The volume fraction of each phase is expressed by a scalar function, which is referred to as a volume function or a color function. To reproduce the interface between immiscible phases, the volume function needs to keep a steep gradient at the interface. However, the steep gradient readily dissipates because VOF methods solve a momentum equation for a mixture of immiscible phases. Therefore, a special treatment is needed for the interface of volume functions. The relative velocity, U_r, is used to compress the interface between the two phases. Weller [21] proposed a relative velocity between two phases, U_r, as follows:

where U is the velocity field, and C_α is a coefficient set to 1 in the present simulation. This method has proven to be reliable in maintaining a sharp interface [21].

3.3 Surface tension force

Surface tension force is calculated using the continuum surface force (CSF) model [22]:

where σ is the surface tension coefficient and κ is the curvature of the interface between the liquid and gas. κ is given by

in which n̂ is the gradient vector at the face, which is given by

where δ_n is a stabilization factor depending on the volume of grid cells. The typical value of δ_n in our simulation is 1.0 × 10^–5m^–1.

3.4 Velocity-pressure coupling

The momentum equation is given by

where ρ is the mixture density, p is the pressure, g is the gravity vector, h is the position vector in the vertical direction, (g·h)∇ρ is the buoyancy force, and τ is the deviatoric stress. The interFoam solver uses the PIMPLE method, which is a combined velocity-pressure coupling algorithm of the SIMPLE (Semi-Implicit Method for Pressure-Linked Equations) and PISO (Pressure Implicit with Splitting of Operator) algorithm [23]. The PIMPLE algorithm is summarized as the following routine.

1. Momentum prediction: Predict the velocity field using the momentum equation.

2. Pressure solution: Solve the pressure equation and correct flux.

3. Explicit velocity correction: Correct the velocity field with the solved pressure field.

The routine is repeated for certain number of times, which was two times in the present simulation.

3.5 Computation and post-processing

A diameter given by the average diameter of 51 experimental drops was adopted as the initial diameter of the numerical drop, D₀ = 3.7 mm. The initial velocity of the drop was determined using the conservation of mechanical energy:

where g is the gravitational acceleration, and x_c,0 is the initial height of the centroid of the drop. The viscosities of water and air were set to 1.0 × 10^–3 Pa · s and 1.84 × 10^–5 Pa · s, respectively. The field of the initial volume fraction of water was set to α = 1.0 at the interior of the drop and α = 0.0 at the outside of the drop. The surface tension coefficient σ, between water and air was set to 0.07 N · m^–1.

The boundary conditions used for the calculation are shown in Table 1. The contact angle between water and air on the bottom patch on which the drop impacted was set to 180° (a perfectly hydrophobic surface), which means that the gradient of α on the bottom boundary is determined as the negative normal vector of the boundary patch. The schematic diagram of the CFD setup is shown in Figure 3. The calculation domain was in the shape of a cube, 1.5 cm on a side (Figure 3(a)). A spherical drop of 3.8 mm diameter was placed at a height of 5 mm above the bottom boundary (Figure 3(a)). To improve efficiency of the calculation, the resolution of the mesh is uniform and finest at the interior of the central rectangular column covering the drop, while it becomes coarser toward the exterior, where the drop will never enter (Figures 3(b) and 3(c)). Number of grid cells above 1.0 cm from the bottom patch was reduced to 30 % of the finest region, and number of grid cells outside a square with a side length 1.0 cm and a center of the drop was reduced to 50 % of the finest region. From the finest grid cells towards boundaries except the bottom, the length of edge of grid cells are expanded with the expansion ratio of smallest length to largest length of the edge length of grid cells. The expansion ratios are 2.81 towards the left, the right, the front, the back patches of the drop and 6.26 towards the upper patch. Two resolutions of the mesh, 5 and 10 grid cells mm^–1 at the finest part of the mesh, were used to validate the effect of the resolution.

Figure 3

Schematic diagram of the CFD setup. (a) System size and initial setup of the volume fraction of water, α. The grid mesh from (b) top and bottom views, and (c) side views. The central cubic domain has the finest and unity resolution, otherwise domain has reduced resolution with expanding cell size outward.

The simulations were performed on a computer equipped with an Intel® Core™ i7-3960X CPU and with 32GB RAM. The simulation results were rendered as movies using ParaView [24]. The interface between water and air was determined by thresholding the volume fraction of water at α = 0.5 to capture the center of transitional region between water and air. Rendered movies were processed using the image processing pipeline that was also used to process the experimental results.

4.1 Assessment of the numerical result

Here, the numerical results are assessed by comparison with the experimental results.

Before assessment of the results from a physical perspective, the effect of the mesh design was validated by evaluating the dependency on the mesh resolution. No particular differences were observed in the results for the two different resolutions, which indicates that the mesh resolution has no significant effect on the result. To inspect the numerical result with the fine resolution, results calculated with the finer mesh (10 cells / mm) were used for further analysis.

Figure 4 shows sequential images of the experimental and simulated drops. The experimental and simulated drops were comparable in that each drop exhibited a stable rebound. The sequence of the deformation (spreading after first impact, forming a disk-shape, shrinking, making a head at the center of the disk, lift-off, shaking of the shape while in the air, and impacting again) was also reproduced in the calculation. For high Weber numbers (We ≥ 15), the experimental result fluctuated, possibly due to asymmetrical expansion and contraction, while the numerical result was stable and had symmetrical expansion and contraction. Deformation for the experimental drop was so sensitive that no symmetrical deformation could be achieved.

Figure 4

Sequential images of experimental drops (upper rows) and simulated drops (lower rows) captured every 5 ms for We = 7 (top), We = 15 (middle), and We = 23 (bottom).

Time evolutions of the height and width of the drop are expected to provide vibrational patterns of the deformation process. Time evolutions of relative diameter D/D₀ for both experimental and simulated drops with different Weber numbers were compared (Figure 5), and the top, middle, and bottom heights of the drops, x₁, x_c, and x₂, respectively (Figure 2 and Figure 6), show that both sets of results have the same vibrational patterns, although for high Weber numbers, the time spans between the first and second expansions and between the first and second impacts for the numerically simulated drops were slightly wider than those for the experimental drops. The time series for the horizontal diameter of the drop during the impact approximately represents how much kinetic energy is converted to surface energy (Figure 5). The time series for the middle height can be considered to represent approximately the potential energy of the drop. Thus, as shown in Figure 6, the time series for the potential energy of the drop for both the experiments and the simulations can be considered to be in agreement.

Figure 5

Time series of width for experimental and numerically simulated drops with (A) We = 7, (B) We = 15, and (C) We = 23. The first expansion and contraction (0-20 ms) have good agreement.

Figure 6

Time series of top height, middle height, and bottom height for experimental and numerically simulated drops with (A) We = 7, (B) We = 15, and (C) We = 23. The difference between the top and bottom heights represents deformation in the vertical direction and the middle height represents the approximate potential energy of the drop.

The dissipated energy during the rebound is very difficult to determine because both the velocity and the surface area of the drop are unknown [25]. One effective way to experimentally estimate the dissipated mechanical energy is to calculate the ratio of the maximum height after the first impact to the initial height, as a ratio of mechanical energy at the maximum height to the initial mechanical energy:

by assuming that the potential energy is equal to the mechanical energy when the drop is at the highest position.

The ratio of mechanical energy at the maximum height to the initial mechanical energy for each Weber number is shown in Figure 7. Both the numerical and experimental results showed a decrease with an increase of the Weber number. The energy loss for the experimental result with high Weber numbers is considered to fluctuate due to asymmetrical deformation during the rebound (Figure 4).

Figure 7

Weber number and the ratio of the mechanical energy of the drop at the maximum height after the impact to the initial mechanical energy. Both the numerical and experimental results showed a decrease in the ratio with an increase of the Weber number.

Through the assessment performed here, the numerical result is considered to be reasonably reliable with respect to the deformation and dissipated energy.

4.2 Quantitation of the energy conversion

Kinetic energy and potential energy were calculated using the following respective equations:

where V is the volume and Ω is the entire domain for the calculation. Under the condition that the width of the interface between water and air is asymptotically limited to zero, the integral over the interface can be reformulated by a volume with the gradient of the volume fraction, ∇α [22]. Thus, the surface energy can be calculated using

The sum of the kinetic and potential energies is the mechanical energy:

In this system, the pressure and volume are considered to be constant, and the energy of interest is the sum of the mechanical and surface energies:

Figure 8 shows the time evolution of the energies calculated from the numerical results. At the impact (t = 0 ms), the mechanical energy begins to decrease rapidly and the surface energy simultaneously begins to increase. When the surface energy reaches a maximum (t ≈ 8 ms), the kinetic energy has a local minimum. After reaching the maximum surface energy, the mechanical energy begins to increase while the surface energy decreases. After takeoff of the drop (t ≈ 15 ms), as evident for high Weber numbers, the conversion between the mechanical energy and surface energy still continues, which is considered to be caused by vibration of the drop in the air. Interestingly, the changes of these energies cancel each other out and are considered to be conserved in the form of the sum of the mechanical and surface energies.

Figure 8

Kinetic, potential, and surface energies of drops as a function of t, the time after impact, with (A) We = 7, (B) We = 15, and (C) We = 23. The dashed line represents the transition point from regime I to regime II.

4.3 Imaginary damped spring model

A poorly elastic shock of a Leidenfrost drop has been modeled by an imaginary spring [4, 6], which is a linear spring model with two mass points that represent the mass of the drop at both ends of the spring. Here, we extend the imaginary spring model by adding a damping term:

where x₁ and x₂ are the heights of the bottom and top of the spring above the plate respectively, and

is the strain of the spring, m is the mass of the drop, D₀ is the initial vertical length of the drop, k is the stiffness of the spring, c is the damping coefficient of the spring (c ≥ 0), and F is the external force loaded at the bottom of the spring. Figure 2 shows a schematic diagram of the imaginary damped spring model. Note that by combining Eqs. (14a) and (14b), the momentum equation for the centroid of the spring, xc=12(x1+x2), can be represented as

which plots the free-fall and bounce-back of the spring.

Let us define the regime in which the drop is in contact with the vapor film over the plate as regime I. In regime I, the height of the bottom of the spring is considered to be fixed (x₂ = 0); therefore,

and

where ϵ_I is the strain in regime I. Equations (14a) and (14b) then become:

where

By solving Eq. (19a), we obtain

where t_I is the time after the impact, A_I is the initial amplitude of the oscillation, ζI=cI2mkI is the damping ratio, ωI=kIm is the undamped angular frequency of the spring, ωd,I=1−ζI2ωI is the under-damped harmonic oscillator, and ψ_I is the phase at the impact.

The time span from the lift-up to the next impact of the drop is defined as regime II. In regime II, the bottom height of the spring is no longer fixed (x₂ ≥ 0) and there is no external force loaded on the bottom mass point (F = 0). Differentiation of Eq. (15) gives

where ϵ_II is the strain in regime II. The combination of Eqs. (14a), (14b), and (22) gives

where

By solving Eq. (23), we obtain

where t_II is the time after the lift-up, A_II is the amplitude of the oscillation, ζII=cII2mkII is the damping ratio, ω_II = kIIm is the undamped angular frequency of the spring, ωId,I=1−ζII2ωII is the under-damped harmonic oscillator, and ψ_II is the phase at lift-off.

The coefficients were obtained according to the description given in Appendix A. The damping coefficient for regime I, c_I, was determined to be 0.7 × 10^–3 kgs^–1 using Eq. (A.5) with the result for We = 7 and was reasonably assigned for all Weber numbers in this study, while that for regime II, c_II, was determined to be half the value of c_I. This difference of the damping coefficient indicates that the mechanism for energy loss is different between regimes I and II. The stiffness k, determined by Eq. (A.6), tends to decrease with an increase of the Weber number.

The sum of the kinetic, potential, and elastic energies as the surface energy of the spring model can be calculated for each regime:

where E_ms,I and E_ms,II are sums of the mechanical energy and the surface energy in regime I and regime II, respectively. These energies are shown in Figure 10 with We = 7, 15, 23.

The overall energy loss rate in regime I is expressed as

and the energy loss rate over 1 cycle of oscillation in regime II,

corresponds well for both the simulation and the spring model (Figure 11).

While the vertical strain and energy loss rates of the drop were well explained by the spring model (Figures 9 and 11), the second decrease of the spring model lagged that of the simulated drop (Figure 10). This time lag indicates that true damping factor has an other period than the damping term of the spring model.

Figure 9

Time series of drop’s vertical strain for the experiment, simulation, and the spring model with We = 7.

Figure 10

Time series for the sum of kinetic, potential, and surface energies from the simulation and spring model results. Each dashed line represents the transitional time from regime I to regime II. Major energy losses were observed at two moments, 2 ms and 12 ms after the impact as indicated by the “*” marks. The dash-dotted line shows t = D₀/U_impact, which predicts the start time of the second energy decay of the drop.

Figure 11

Energy loss rates over regime I λ_I, and the energy loss rate over 1 cycle of oscillation in regime II λ_II, with the simulated drop and the spring model.

4.4 Spherical harmonic analysis

A water drop on an oscillating plate [26, 27, 28] or on a superheated plate [29] is known to show spherical harmonic oscillation due to the surface tension force. We will compare the deformations of the drop on the impact with the spherical harmonic oscillation.

An oscillating drop can be represented as a linear combination of spherical harmonic functions as [30, 31, 32]

where R is the radius of an unperturbed sphere-shape drop, A_n (t) = a_n cos(ω_nt + ψ_n) are intensities of the modes, and Ynm are Laplace’s spherical harmonics with orders n = 0, 1, 2, ⋯, and degrees m = –n, –(n – 1), …, n – 1, n.

The deformation of the drop in this study is radially symmetric, thus we only consider Laplace’s spherical harmonics with 0 degree, Yn0. Furthermore, we assume that the deformation of the impacting drop is aligned in vertical or radial axes, so oscillation modes can be narrowed down to n = 1 and n = 2. Therefore Eq. (30) is simplified as

Both of Y10andY20 modes are radially symmetric, but Y20 is horizontally symmetric while Y10 is horizontally antisymmetric, as shown in sequential images of R+A1Y10 and R+A2Y20 in Figure 12(a). Using this difference, harmonic phase shifts of these two modes, ψ₁ and ψ₂ are determined by discriminating the symmetry of deformations of the drop.

Figure 12

A sequential analysis of the deformation of the impacting drop with shared time axis, with We = 7. (a) Spherical harmonic deformations represented by R+A1Y10+A2Y20,R+A1Y10,andR+A2Y20 with R = 1.85 mm, a₁ = 2 mm, and a₂ = 2 mm, together with that of drops of the experiment (Exp.) and the CFD simulation (Sim.). (b) Intensities of Y10andY20,A₁ and A₂, respectively, are plotted as functions of the time after the first impact, t. (c) Sum of mechanical and surface tension of simulated drop is plotted as a funcion of t. Dashed lines in (b) indicates moments of A₁ = 0, at which exponentially energy decays starts as shown in (c).

Immediately after the impact (t = 0 ms), the drop is in a spherical shape, thus the phase of mode n = 2 must be π2or−π2. After impacting on the plate, the drop spreads radially and forms into a disk-shape, which corresponds to A₂ < 0, therefore ψ2=π2.

The determination of the harmonic phase of the Y₁ mode is more difficult than the Y₂ mode. At the moment of the maximum width of spreading drop, the drop is in a disk-shape which is horizontally symmetric, therefore A₁ = 0. From this moment to the lift-off, the drop forms an antisymmetric matryoshka-shape, which indicates A₁ < 0. After the lift-off, the drop forms into a vertical peanut-shape (e.g. 25 ms) which is horizontally symmetric thus A₁ = 0. By extrapolating from these conditions, A₁ must be a positive value at t = 0. By looking carefully at energy decay curves in Figure 10, we found that the second energy decay starts slightly earlier at increasing Weber number. Let us consider a free end reflection of the drop as a pulse over the vertical axis on the bottom plate as a free end. We envisage that after the impact the drop receives the reflection of its impact velocity and induces a vertical uplift (i.e. A₁ > 0), until the pulse passes over the end. Therefore, a length of the time during which the drop passes over its diameter with the impact velocity, D₀/U_impact, is considered to characterize the phase of the Y10 mode. As shown by the dash-dotted line in Figure 10, t = D₀/U_impact predicts beginning time of the second decay well.

The deformation of Eq. (31) with determined phases ψ₁ and ψ₂ together with experimental and simulated drops is shown in Figure 12(a). The deformation sequence of spherical shape, disk-shape, matrioshka-shape, and peanut-shape, is reproduced by Eq. (31).

A₁ and A₂ as functions of time after the impact are shown in Figure 12(b). Actually, the imaginary spring model represents the Y20 mode: A₂ has same vibration mode with the drop’s vertical strains, ϵ, which is shown in Figure 9. This is because the vertical strain of a Y10 mode is always zero and thus height of the shape represented by Eq. (31) only depends on the Y20 mode.

One of the most important insights from the spherical harmonic analysis is the existence of an another vibrational mode, Y10 in the deformation of the impacting drop, which is not considered in the spring model.

4.5 Energy loss upon the impact

Time series of A₁ and A₂ were compared with that of the total energy E_ms of the simulated drop with We = 7 (Figures 12(b) and (c)). We found that the cycle of the Y10 mode is synchronous with that of the repetitive energy decay. At the time spans when Y10 increases its amplitude (dA1dt>0), the energy starts an exponential decay. Figure 10 shows that the damped imaginary spring model predicts lagged second decay behind the simulated drop. As mentioned in section 4.4, the spring model represents the Y20 mode. Considering that the cycle of the energy decay depends on the Y10 mode and the spring model does not consider the Y10 mode, the mismatch of the decay timing can be explained.

When the drop impacts at the bottom plate, the bottom part of the drop is forced to stop and the vertical velocity is forced to be zero suddenly. This sudden change of the velocity induces the pressure surge (i.e., stagnation pressure). Figure 13(a) shows that a pressure surge occurs at the beginning of the impact. The amount of the pressure surge is roughly estimated as ρD₀U_impact/Δt ∼ 200 Pa, where Δt ≈ D₀/U_impact ∼ 10 ms is the time span to stop the free-fall motion of the drop. This phenomenon is similar to a water hammer with slow valve closure (slower than sound propagation), in which a pressure surge occurs when a fluid in motion is forced to stop. The pressure surge accompanies a pressure-gradient force, which is expressed by the term –∇p in the Eq. (6).

Figure 13

Pressure (left) and magnitude of the velocity (right) inside the drop, from a cross-sectional lateral view at the drop center (We = 7). (a) At the beginning of the impact (t = 2 ms), a pressure surge occurs at the bottom part of the drop, where magnitude of the velocity is nearly zero. and (b) At the retraction from the disk-shaped drop (t = 12 ms), the pressure-induced motion suppresses the retracting motion.

In many cases of fluid dynamics, a pressure-gradient force is a driving force of a flow (e.g., a channel flow). However, in this case, the pressure-gradient force consequently dumps the motion inside the drop. The pressure-induced force associated with the stagnation pressure is an important factor in the Drop Deformation and Breakup (DDB) model [33, 34] introduced by Ibrahim et al., which successfully predicts the deformation of spray drops.

The pressure-gradient force generated at the impact induces an upward flow inside the drop in the time span of D₀/U_impact from the impact. The upward flow is subsequently reflected at the end of the drop due to the surface tension. Therefore, the flow induced by the pressure-gradient force generates a vertical vibrational motion. The Y10 mode found in our spherical harmonics analysis represents this vertical vibrational motion.

Meanwhile after the impact, a major part of the free-fall motion of the drop is converted to radially spreading motion, which is also subsequently reflected at the end of the drop due to the surface tension. This motion forms a vibrational motion represented by the Y20 mode. Therefore, the free-fall motion of the drop before the impact is converted to two motions after the impact, the pressure-induced motion and inertial motion, represented by Y10 and Y20 modes, respectively.

To investigate the breakdown of the force acting to dump inside the drop, viscous force, pressure-gradient force, and total external force (viscous force, pressure-gradient force, surface tension force, and gravitational force) were summed inside the drop of the CFD simulation for each time step. Specifically, we calculated an α-weighted summation of each force over grid cells with a negative value of inner product of a force and velocity, as expressed in ∑_i,fi·Ui<0α_if_i, where i is an index of grid cells, f_i is a force at i-th cell, U_i is velocity at i-th cell, and α_i is the volume fraction of water phase at i-th cell. Considering the deformation of the bounce of the 3D drop where spreading and shrinking in vertical and radial directions, these forces were split into vertical and radial components. The radial force component shows major forces within the drop in a disk-shape, while the vertical component shows ones within the drop in a cylinder-shape. Figure 14 shows that amongst the forces damping the motion of the drop, the pressure-gradient force dominates the external forces (right-hand side of Eq. (6)) and the viscosity effect is fairly small. The small impact of the viscosity on the drop deformation was also reported by Renardy et al. [11].

Figure 14

The viscous force ∇ · τ, pressure-gradient force –∇p, and the external forces (right-hand side of Eq. (6)) damping the motion of the drop for the vertical and radial directions were obtained from the simulated results (We = 7). The dashed line represents the transition point from regime I to regime II. The pressure-gradient force dominates the external force and the contribution of the viscous force is fairly small. At the two major energy losses (*), the external forces damp the velocity, the former against the vertical direction and the latter against the radial direction.

The pressure-induced motion resists the free-fall inertial motion at the beginning of the impact (t = 2 ms), and then resists the inertial motion when the drop is retracting from the disk-shape (t = 12 ms). Figure 14 shows that the pressure-gradient force dominantly resists the motion of the drop at these two timings. At the retraction of the drop, the direction of the pressure-induced motion is inherently downward, however, due to the disk-shaped drop as a flow field and the existence of the bottom plate, it is forced to advance radially (Figure 13(b)). As the result, the retracting inertial motion of the drop is dumped by this radially spreading pressure-induced motion.

Most part of the pressure-induced motion decays during the two resistances to the inertial motion, however, still remains with a small intensity after the lift-off, causing a small energy loss starting at t = 25 ms as shown in Figures 10 and 12(c).

Note that we have considered just the first dry-rebound. Biance et al. [4] have shown restitution coefficients of successive dry-rebounds of a drop with diameter of 1 mm. They reported that the restitution coefficient e is relatively low at the first impact (e ∼ We^–1/2, called as poorly elastic shocks) but very close to 1 after multiple bounces (called as quasi-elastic shocks). With lower Weber number, the energy loss tends to be small as shown in Figure 7. Thus, one reason of the small energy loss after multiple bounce is the low Weber number. Biance et al. also reported that in the quasi-elastic shocks the vibration of drop’s diameter is in phase with the flight of the drop. The other reason of the small energy loss after multiple bounce is considered that the pressure surge disappears due to the inertial motion synchronized with the bounce of the drop. The synchronized Y20 mode, which suppresses the impact velocity, avoids the sudden stop at the bottom and generating the pressure surge.