1. Introduction
Drop impacts are common in everyday life, industry and nature. Examples include drops splashing on a kitchen sink, precision inkjet printing (Lohse Reference Lohse2022) and raindrops striking leaves and animal fur (Roth-Nebelsick et al. Reference Roth-Nebelsick, Konrad, Ebner, Miranda, Thielen and Nebelsick2022). Since the seminal studies of Worthington (Worthington Reference Worthington1877, Reference Worthington1908) and Edgerton’s iconic images (Edgerton & Killian Reference Edgerton and Killian1954), researchers have strived to uncover the physics of impact phenomena, which include drop deposition, splashing and rebound (Rein Reference Rein1993; Yarin Reference Yarin2006). Early studies primarily focused on splashing, but interest in drop rebound has surged, inspired in part by raindrops bouncing off superhydrophobic leaves (Chen et al. Reference Chen, Xiao, Chan, Lee and Li2011) and driven by advances in surface preparation techniques that can render surfaces hydrophobic and enable rebound (Richard & Quéré Reference Richard and Quéré2000; Okumura et al. Reference Okumura, Chevy, Richard, Quéré and Clanet2003). Despite the impact of these innovations on numerous everyday items, such as windshields and eyewear, our understanding of the physics governing rebound is mostly limited to millimetre-sized drops striking surfaces at moderate and high speeds. This paucity motivates a detailed study of low-speed rebound of submillimetric drops, which are typical in commercial spray nozzles (Makhnenko et al. Reference Makhnenko, Alonzi, Fredericks, Colby and Dutcher2021; Chen & Ashgriz Reference Chen and Ashgriz2022), naturally emerge from ligament breakup (Villermaux Reference Villermaux2020) and can result from splashing of millimetric drops (Yarin Reference Yarin2006). Recent attempts to fill this knowledge gap have been afflicted with contact line friction, which increasingly alters the rebound physics as the impact velocity decreases (Wang et al. Reference Wang, Zhao, Sun, Mehrizi, Lin, Guo and Chen2022; Thenarianto et al. Reference Thenarianto, Koh, Lin, Jokinen and Daniel2023). To address this knowledge gap, we perform experiments on a perfectly non-wetting substrate – a thin viscous film – using submillimetre-sized drops and identify the relevant physics governing rebound using reduced-order modelling, direct numerical simulation (DNS) and energetic arguments.
Over a century ago, Lord Rayleigh observed drops rebounding upon impact (Rayleigh Reference Rayleigh1879) and attributed this phenomenon to an intervening air layer separating the drops and preventing coalescence (Rayleigh Reference Raleigh1899). Since then, experiments conducted under modified atmospheric conditions and corresponding theoretical studies have confirmed his hypothesis (Willis & Orme Reference Willis and Orme2000; Bach, Koch & Gopinath Reference Bach, Koch and Gopinath2004; Chubynsky et al. Reference Chubynsky, Belousov, Lockerby and Sprittles2020). Meanwhile, advancements in experimental and computational techniques have enabled measurement of the minimum thickness of the air layer – of the order of 100 nm – below which van der Waals forces bridge the liquid–air interfaces and capillary forces drive coalescence (Couder et al. Reference Couder, Fort, Gautier and Boudaoud2005a ; Kolinski et al. Reference Kolinski, Rubinstein, Mandre, Brenner, Weitz and Mahadevan2012; de Ruiter et al. Reference de Ruiter, Oh, van den Ende and Mugele2012; Kavehpour Reference Kavehpour2015; Lewin-Jones, Lockerby & Sprittles Reference Lewin-Jones, Lockerby and Sprittles2024; Sprittles Reference Sprittles2024). Thus, drops only rebound without contact from smooth surfaces. Although there are limited solid surfaces smooth enough to promote contactless rebound, liquid–gas interfaces are atomically smooth and ideal for rebound. For example, recent results for the rebound of a drop on a bath of the same fluid showed that the deformation of the drop and bath are dynamically coupled (Alventosa, Cimpeanu & Harris Reference Alventosa, Cimpeanu and Harris2023; Phillips & Milewski Reference Phillips and Milewski2024), and that vibrating the bath can lead to long-lasting bouncing and wave-propelled drop motion (Couder et al. Reference Couder, Protiere, Fort and Boudaoud2005b ; Moláček & Bush Reference Moláček and Bush2013). This dynamic coupling is further evidenced by the complex bubble patterns observed when an impacting drop merges with a bath (Saylor & Bounds Reference Saylor and Bounds2012). If the bath height is reduced to a thin film, rebound can become substrate-independent (Hao et al. Reference Hao2015; Sanjay et al. Reference Sanjay, Lakshman, Chantelot, Snoeijer and Lohse2023b ). Although drops can rebound from heated substrates (Biance et al. Reference Biance, Chevy, Clanet, Lagubeau and Quéré2006), dry ice (Li et al. Reference Li, Liu, Wang and Chen2024) and when the drop is coated with hydrophobic powder (Aussillous & Quéré Reference Aussillous and Quéré2006), these methods inject energy into the system or modify the drop composition, potentially altering the rebound dynamics and obfuscating the fundamental physics.
A drop of silicone oil (
$R = 0.31$
mm) rebounding from a glass slide coated with a thin, viscous film (
$\varGamma / Oh_{\!f}^{1/3} \approx 10^{-2}$
). Here,
$\textit{We}=1.03$
,
$Oh=0.02$
and
$Bo=0.04$
. Subsequent images are separated by 1.44 ms and shifted 0.4 mm to the right.
Rebound dynamics is governed by a competition of inertial, viscous, capillary and gravitational forces. A common set of governing dimensionless parameters describing the problem are
\begin{equation} \begin{array}{ccc} \quad \textit{We} = \dfrac {\rho \, V^2 \, R}{\sigma }, & \quad Bo = \dfrac {\rho \,g \, R^2}{\sigma }, & \quad Oh = \dfrac {\mu }{\sqrt {\rho \,\sigma \, R}}. \quad \end{array} \end{equation}
The Weber number
$\textit{We}$
compares inertial and surface tension forces, Bond number
$Bo$
compares gravitational and surface tension forces and Ohnesorge number
$Oh$
is the ratio of the inertio-capillary to the viscous time scale. Here,
$R$
is the undeformed drop radius,
$V$
is the impact velocity,
$\rho$
is the liquid density,
$\mu$
is the dynamic viscosity,
$\sigma$
is the surface tension coefficient and
$g$
is the gravitational acceleration. Although early studies emphasised
$\textit{We}$
as the primary parameter controlling rebound (Foote Reference Foote1975; Anders, Roth & Frohn Reference Anders, Roth and Frohn1993; Richard & Quéré Reference Richard and Quéré2000; Richard, Clanet & Quéré Reference Richard, Clanet and Quéré2002; Okumura et al. Reference Okumura, Chevy, Richard, Quéré and Clanet2003), recent work has shown that increasing
$Bo$
and
$Oh$
can suppress rebound for
$\textit{We}\gt 1$
(Jha et al. Reference Jha, Chantelot, Clanet and Quéré2020; Sanjay et al. Reference Sanjay, Chantelot and Lohse2023a
). For
$\textit{We}\lt 1$
, the role of
$Bo$
and
$Oh$
is less explored, and attempts to isolate the role of
$\textit{We}$
in rebound have been confounded by the increasing influence of contact line friction as
$\textit{We}$
decreases (Wang et al. Reference Wang, Zhao, Sun, Mehrizi, Lin, Guo and Chen2022; Thenarianto et al. Reference Thenarianto, Koh, Lin, Jokinen and Daniel2023). Impacting a drop on a substrate coated with a viscous film avoids contact line formation, as discussed by Gilet et al. (Reference Gilet, Terwagne, Vandewalle and Dorbolo2008) and Gilet & Bush (Reference Gilet and Bush2012). This technique produces substrate-independent rebound when the film is sufficiently thin and viscous (Hao et al. Reference Hao2015). To determine the criteria where drop rebound is substrate-independent, Sanjay et al. (Reference Sanjay, Lakshman, Chantelot, Snoeijer and Lohse2023b
) analysed the competing influences of the drop Ohnesorge number
$Oh$
and film Ohnesorge number
$Oh_{\!f}=\mu _{\!f} / ( \rho \sigma R )^{1/2}$
, and found that substrate-independent rebound occurs for
$Oh_{\!f} \gtrsim 0.1$
when
$\varGamma \lt \textrm {B} \, Oh_{\!f}^{1/3}$
, where
$\varGamma = h_{\!f}/R$
and
$\textrm {B} = \textrm {B} ( Oh )$
. Since
$\textrm {B} \approx 0.1$
when
$Oh \approx 0.1$
, and
$\textrm {B}$
increases with decreasing
$Oh$
, a conservative condition for substrate-independent bouncing is
$\varGamma /Oh_{\!f}^{1/3} \lt 0.1$
. Notably, this result is derived for
$\textit{We} \gt 1$
and becomes less restrictive for
$\textit{We} \lt 1$
. By following these guidelines, a rigid substrate can be prepared that allows a drop to rebound without contact line formation or dynamic influence from the viscous film. An example of such a rebound is imaged in figure 1.
When a drop impacts a rigid substrate, it first spreads until reaching a maximum equatorial radius
$R_m$
and maximum contact radius
$R_c$
. When
$\textit{We} \gg 1$
, the dimensionless maximum equatorial radius
$\beta = R_{m}/R$
and maximum contact radius
$\beta _{c}= R_{c}/R$
are approximately equal, and two competing scaling laws have been proposed to describe their evolution with
$\textit{We}$
. Clanet et al. (Reference Clanet, Béguin, Richard and Quéré2004) proposed
$\beta \sim \textit{We}^{1/4}$
based on momentum conservation and assuming the drop deforms into a cylinder with negligible height during impact. Alternatively, the conservation of energy between the kinetic energy of the drop at impact and the surface energy at the maximum deformation of the drop predicts
$\beta \sim \textit{We}^{1/2}$
(Chandra & Avedisian Reference Chandra and Avedisian1991; Bennett & Poulikakos Reference Bennett and Poulikakos1993; Okumura et al. Reference Okumura, Chevy, Richard, Quéré and Clanet2003). Although additional arguments and empirical evidence favour the scaling
$\textit{We}^{1/2}$
(Eggers et al. Reference Eggers, Fontelos, Josserand and Zaleski2010; Villermaux & Bossa Reference Villermaux and Bossa2011; Laan et al. Reference Laan, de Bruin, Bartolo, Josserand and Bonn2014; Josserand & Thoroddsen Reference Josserand and Thoroddsen2016; Lee et al. Reference Lee, Laan, de Bruin, Skantzaris, Shahidzadeh, Derome, Carmeliet and Bonn2016), extensive experimental results support the scaling
$\textit{We}^{1/4}$
over a wide range of
$\textit{We}$
(Bartolo, Josserand & Bonn Reference Bartolo, Josserand and Bonn2005; Biance et al. Reference Biance, Chevy, Clanet, Lagubeau and Quéré2006; Tsai et al. Reference Tsai, Hendrix, Dijkstra, Shui and Lohse2011; García-Geijo et al. Reference García-Geijo, Riboux and Gordillo2020; Zhan et al. Reference Zhan, Lu, Liu, Wang, Lv and Liu2021; Li et al. Reference Li, Liu, Wang and Chen2024). When viscous forces are important, these scalings break down (Rein Reference Rein1993), and describing the crossover from the capillary to viscous regime requires consideration of finite
$Oh$
effects (or
$Re=\sqrt {We}/Oh=\rho V R / \mu$
) (Roisman Reference Roisman2009; Eggers et al. Reference Eggers, Fontelos, Josserand and Zaleski2010; Laan et al. Reference Laan, de Bruin, Bartolo, Josserand and Bonn2014; Lee et al. Reference Lee, Laan, de Bruin, Skantzaris, Shahidzadeh, Derome, Carmeliet and Bonn2016; Sanjay & Lohse Reference Sanjay and Lohse2025). For
$\textit{We}\lt 1$
, the deformed drop resembles an oblate ellipsoid rather than a wide cylinder and
$\beta \not \approx \beta _{c}$
. In this case, Richard & Quéré (Reference Richard and Quéré2000) predicted
$\beta -1 = (5/48 )^{1/2} \, \textit{We}^{1/2}$
by assuming the initial kinetic energy of the drop was converted to surface energy but did not test their prediction across a range of
$\textit{We}$
. Recently, efforts to describe and relate
$\beta$
and
$\beta _{c}$
during low-
$\textit{We}$
rebound have used drops with non-negligible
$Bo$
and relied on fits to empirical relations (Liu, Liu & Chen Reference Liu, Liu and Chen2023; Li et al. Reference Li, Liu, Wang and Chen2024). Thus, a comprehensive and validated understanding of spreading for low
$\textit{We}$
is still lacking.
Following spreading, surface tension drives the drop to retract and can lead to rebound. When
$Oh$
is negligible, the maximum retraction rate is independent of impact velocity, and the retraction time
$t_{r} \sim t_{\sigma }$
, where
$t_{\sigma }= ( \rho R^{3} / \sigma )^{1/2}$
is the intertio-capillary time (Bartolo et al. Reference Bartolo, Josserand and Bonn2005). This result is valid for
$\textit{We}\gt 1$
, and since inertia also governs the spreading time
$t_s$
, the interaction time between the substrate and drop
$t_{c}=t_{s}+t_{r}$
, referred to as the contact time henceforth, scales as
$t_c \sim A \, t_{\sigma }$
(Richard et al. Reference Richard, Clanet and Quéré2002), where
$A=f (Oh )$
(Jha et al. Reference Jha, Chantelot, Clanet and Quéré2020). Alternatively, when
$\textit{We}\lt 1$
, the non-dimensional contact time
$t_c/t_{\sigma }$
increases as
$\textit{We}$
decreases (Foote Reference Foote1975; Okumura et al. Reference Okumura, Chevy, Richard, Quéré and Clanet2003). This transition in
$\textit{We}$
dependence signals a shift in the underlying physics. This shift also alters the
$\textit{We}$
dependence of the coefficient of restitution
$\varepsilon$
, a measure of the elasticity of the rebound. For
$\textit{We}\gt 1$
, Anders et al. (Reference Anders, Roth and Frohn1993) noted that
$\varepsilon$
is
$\textit{We}$
dependent and Biance et al. (Reference Biance, Chevy, Clanet, Lagubeau and Quéré2006) showed that
$\varepsilon \sim \textit{We}^{-1/2}$
, which has been verified on various rigid substrates (Wang et al. Reference Wang, Zhao, Sun, Mehrizi, Lin, Guo and Chen2022; Thenarianto et al. Reference Thenarianto, Koh, Lin, Jokinen and Daniel2023). Additionally, Jha et al. (Reference Jha, Chantelot, Clanet and Quéré2020) showed that increasing
$Oh$
decreases the coefficient of restitution
$\varepsilon$
. For
$\textit{We}\lt 1$
, Gilet & Bush (Reference Gilet and Bush2012) suggested
$\textit{We}$
-independent rebound such that
$\varepsilon = f ( Oh )$
. Recent efforts to study
$\varepsilon$
for low-
$\textit{We}$
rebounds observed a peak in
$\varepsilon$
as
$\textit{We}$
decreases, followed by a decrease to
$\varepsilon = 0$
(no rebound) at finite
$\textit{We}$
(Wang et al. Reference Wang, Zhao, Sun, Mehrizi, Lin, Guo and Chen2022; Thenarianto et al. Reference Thenarianto, Koh, Lin, Jokinen and Daniel2023). The suggested functional form of this decrease
$\varepsilon \sim \textit{We}^{1/2}$
is attributed to contact line friction, highlighting the need for contact-line-free conditions to determine the
$\varepsilon$
–
$\textit{We}$
relationship at low
$\textit{We}$
.
These unique trends in the rebound metrics at low
$\textit{We}$
have inspired several models. Okumura et al. (Reference Okumura, Chevy, Richard, Quéré and Clanet2003) derived a simple model for the contact time
$t_{c}$
and maximum equatorial deformation
$\beta$
by treating the drop as a spring–mass system, with spring constant scaling with the surface tension coefficient. Jha et al. (Reference Jha, Chantelot, Clanet and Quéré2020) expanded this model to include damping to describe the rebound dynamics of viscous drops. Moláček & Bush (Reference Moláček and Bush2012) proposed a quasi-static model for impacts on non-wetting surfaces, approximating drop shapes as static during impact, giving analytical predictions for
$t_{c}$
and
$\varepsilon$
in various limits. Concurrently, Chevy et al. (Reference Chevy, Chepelianskii, Quéré and Raphaël2012) derived a quasi-static model for droplet deformation, highlighting the nonlinear elastic behaviour of weakly deformed drops. More recently, Galeano-Rios, Milewski & Vanden-Broeck (Reference Galeano-Rios, Milewski and Vanden-Broeck2017) introduced a kinematic match (KM) method to derive a fully predictive model for the rebound of a rigid hydrophobic sphere on a liquid bath, which has been extensively validated through experiments and DNS (Galeano-Rios et al. Reference Galeano-Rios, Milewski and Vanden-Broeck2017, Reference Galeano-Rios, Milewski and Vanden-Broeck2019, Reference Galeano-Rios, Cimpeanu, Bauman, MacEwen, Milewski and Harris2021). The KM method imposes natural kinematic and geometric constraints on the motion of interacting surfaces, ensuring continuity at the contact interface. This approach is particularly well suited for rebound problems as it provides a robust framework to capture the interplay between fluid dynamics and interfaces, making it highly versatile for a wide range of impact scenarios in mechanics. For example, Agüero et al. (Reference Agüero, Alventosa, Harris and Galeano-Rios2022) applied the KM method to non-Hertzian impact of a solid sphere on an elastic membrane and observed good agreement with experiment. These successful applications of the KM method and its ability to efficiently simulate rebound phenomena with minimal assumptions motivate its extension to drop rebound on rigid substrates.
Herein, we present experiments and DNS results of drop rebound at low
$\textit{We}$
, and develop a robust reduced-order model using the KM method. The experimental apparatus, testing protocol and image processing technique are detailed in
$\S$
2. In
$\S$
3, we present the formulation of our reduced-order model and describe our DNS framework. In
$\S$
4, we compare our experiments and DNS results with our model predictions for the pertinent metrics describing drop rebound (
$\S$
4.1), spreading and retraction dynamics (
$\S$
4.2) and maximum drop deformation and spreading (
$\S$
4.3). In
$\S$
5, we synthesise our findings into new knowledge about drop rebound and contextualise it within the existing literature. Finally, we summarise our key findings in
$\S$
6 and identify a number of open questions inspired by this work.
2. Experiment
2.1. Experimental set-up
We performed experiments using a three-dimensionally (3D) printed drop generator mounted on a motorised linear stage to produce submillimetre-sized silicone oil drops with less than 1 % variability (Harris, Liu & Bush Reference Harris, Liu and Bush2015; Ionkin & Harris Reference Ionkin and Harris2018), as illustrated in figure 2(
$a$
). The latest design, employed for this work, is documented at https://github.com/harrislab-brown/DropGen. Among other improvements, this design allows for precise height control via a motorised linear state and is adapted to use commercially available nozzles from fused deposition modelling 3D printers. The drop generator expels a thin liquid jet from a nozzle by deforming a piezoelectric disk and then quickly relaxes it, retracting the thin jet while surface tension drives pinch-off of a single drop. By timing drop pinch-off to occur as the jet retracts, the drop is directed upward and ‘floats’ momentarily, allowing interfacial oscillations to decay considerably. Here, we report only experiments where the ratio of the drop diameter in the horizontal and vertical directions is between 0.96 and 1.04 just before impact to minimise the effect of residual drop oscillations on the rebound metrics (Biance et al. Reference Biance, Chevy, Clanet, Lagubeau and Quéré2006; Yun & Kim Reference Yun and Kim2019). Supplementary movie 1 available at https://doi.org/10.1017/jfm.2025.10589 demonstrates drop generation. The radius of the drop
$R$
and its velocity at impact
$V$
were varied by using nozzles with different inner diameters and changing the distance between the nozzle and substrate, respectively. The working fluids were silicone oils (Clearco Products) with slightly varying density
$\rho$
and surface tension
$\sigma$
, but a wide range of dynamic viscosity
$\mu$
, which allowed us to explore the role of drop viscosity in rebound. The ranges of the dimensional and non-dimensional parameters explored are summarised in table 1.
Parameter ranges explored in our experiments.
(
$a$
) A rendering of the experimental set-up. The inset depicts the drop generation process with the evolution of the liquid–air interface (left half) illustrated with coloured lines increasing with time from yellow (light) to blue (dark). (
$b$
) Schematic of an impacting drop and measured variables. (
$c$
) Height of the centre of mass
$h$
against time
$t$
for a rebounding drop (
$R=0.2$
mm), where
$\textit{We}=0.25$
,
$Oh=0.03$
and
$Bo=0.02$
. Dashed and dot-dashed lines are parabolic fits to the drop position during free-fall and rebound, respectively. Grey dashed lines indicate the initial contact time
$t=0$
and detachment time
$t=t_c$
.
An atomically smooth, non-wetting substrate was prepared by coating a standard glass slide with a thin film of viscous silicone oil (Gilet & Bush Reference Gilet and Bush2012; Hao et al. Reference Hao2015; Sanjay et al. Reference Sanjay, Lakshman, Chantelot, Snoeijer and Lohse2023b
). To achieve this, we first coated the slide with a thick layer of viscous silicone oil of
$\mu =975$
g
$\,$
cm
$^{-1}$
s
$^{-1}$
. Then, the slide was passed beneath a fixed knife edge to thin and level the film. The mean film thickness
$h_{\!f}$
was measured by weighing the slide, and was
$h_{\!f} = 231 \pm 55$
μm across all experiments. This thickness corresponds to a dimensionless film thickness
$\varGamma = h_{\!f} / R \approx 10^{0}$
and film mobility
$\varGamma / Oh_{\!f}^{1/3} \approx 10^{-2}$
that is well below the
$\varGamma / Oh_{\!f}^{1/3} \leqslant 10^{-1}$
criterion discussed earlier and thus confidently within the regime of substrate-independent drop bouncing (Sanjay et al. Reference Sanjay, Lakshman, Chantelot, Snoeijer and Lohse2023b
). After the slide was prepared it was placed in a closed environment for at least 12 h to allow non-uniformities in the film to relax under gravity. During testing, the substrate was mounted perpendicular to the direction of gravitational acceleration to ensure the air layer between the impacting drop and substrate did not rupture prematurely due to asymmetry (Lo, Liu & Xu Reference Lo, Liu and Xu2017).
2.2. Experimental protocol
To begin a set of experiments, a reservoir connected to the drop generator is filled with silicone oil, and the system is manually purged to remove any trapped air (Harris et al. Reference Harris, Liu and Bush2015). Next, the pulse width controlling the piezoelectric disk actuation is tuned to produce a single drop with minimal surface oscillations at impact. These settings vary based on the working fluid and nozzle. Next, the nozzle is lowered until it rests just above the substrate, and the system height is set to zero. The nozzle is raised to the desired initial drop height, typically 1 or 2 mm, and a droplet is generated. The drop falls under the action of gravity and its motion is recorded from release until bouncing ceases. After coming to rest, an external pulse of air is applied and sweeps the drop from the substrate before it can coalesce with the film, ensuring the surface remains pristine for subsequent experiments. In the few cases where the drop merged with the viscous film, the slide was translated so the next impact occurred sufficiently far from the merge site. The nozzle height is raised incrementally, and three to eight drops are recorded at each height. The experiments conclude once impact stresses become strong enough to pierce the intervening air layer (Hao et al. Reference Hao2015), which occurs at
$\textit{We} \approx 10$
for all
$Oh$
, or drop retraction leads to droplet ejection for low-
$Oh$
drops (Richard & Quéré Reference Richard and Quéré2000; Bartolo, Josserand & Bonn Reference Bartolo, Josserand and Bonn2006), as shown in supplementary movie 2.
The silhouette of the drop was illuminated by an LED backlight and recorded using a Phantom Miro LC311 high-speed camera with a Laowa 25 mm ultra-macro lens. The frame rate ranged from 9000 to 20 000 frames per second. The pixel size in the recordings was calibrated using a microscope calibration slide. The spatial resolution in our experiments varied between 0.003 and 0.008 mm depending on the size of the drop – smaller drops required higher magnification – but, in general, was approximately
$0.02 R$
. Recordings of each experiment were processed using MATLAB to identify the drop silhouette using Canny edge detection and calculated its centre of mass by assuming the drop shape is a surface of revolution about the vertical axis passing through its centroid, and the drop has uniform density. Time was defined such that
$t = 0$
when the drop ‘contacts’ the substrate and
$t=t_{c}$
when it detaches, where the drop and substrate are considered in ‘contact’ when a portion of the drop and substrate occupy the same pixel. As noted previously, no actual contact occurs between the droplet and substrate due to the microscopic mediating air layer. During rebound, the vertical location of the centre of mass
$h$
and top of the drop
$H$
are tracked, as well as the equatorial radius
$r_m$
and contact radius
$r_c$
of the drop, as defined in figure 2(
$b$
). We define
$H$
as the maximum vertical distance between the substrate and any point on the surface of the drop. Consequently, when drop deformation is small, the region at the north pole of the drop is convex upward and
$H$
coincides with the drop height along the vertical axis passing through is centroid. When drop deformation is large, the location along the surface of the drop where
$H$
is measured shifts toward its rim.
We calculate the undisturbed drop radius
$R$
by measuring the average projected area of the drop in the 30 frames preceding impact, and then solve for
$R$
assuming the projected area is a circle. The centre of mass
$h$
during these frames is fitted to a parabola and the velocity of the drop at impact is
$V = \textrm{d}h/\textrm{d}t(t=0)$
. Similarly, the rebound velocity is
$V_r=\textrm{d}h/\textrm{d}t(t=t_{c})$
, where
$h$
is fitted to a parabola for the 30 frames following detachment. For experiments where less than 30 frames separated detachment and subsequent impact, all intervening frames were used to calculate
$V_r$
. Figure 2(
$c$
) plots the centre of mass of a drop
$h$
against time
$t$
, before (blue markers), during (black markers) and after (yellow markers) impact for
$\textit{We}=0.25$
,
$Oh=0.03$
and
$Bo=0.02$
. The vertical dashed lines indicate the initial contact time
$t=0$
and detachment time
$t=t_c$
of the drop, and the dashed and dot-dashed black curves are parabolic fits for
$h$
during free-fall and rebound, respectively. Snapshots of the drop shape during free-fall, contact and rebound are inset.
The elasticity of rebound can be characterised by the coefficient of restitution
$\varepsilon$
:
\begin{equation} \varepsilon ^{2} =\frac {E_{{r}}}{E_{{i}}}= \frac {g(h(t_{c}) - R)+0.5 V_{r}^{2}}{0.5 V^2}, \end{equation}
which compares the energy of the drop at rebound
$E_{ {r}}$
and impact
$E_{ {i}}$
. Although
$E_{ {i}}$
is purely kinetic,
$E_{ {r}}$
has kinetic and gravitational potential contributions owing to drop deformation at detachment. In most prior studies the potential energy has been ignored and
$\varepsilon = V_{r} / V$
(Richard & Quéré Reference Richard and Quéré2000; Biance et al. Reference Biance, Chevy, Clanet, Lagubeau and Quéré2006; De Ruiter et al. Reference De Ruiter, Lagraauw, Van Den Ende and Mugele2015; Hao et al. Reference Hao2015; Jha et al. Reference Jha, Chantelot, Clanet and Quéré2020). Ignoring potential energy implies a spherical drop at detachment (Richard & Quéré Reference Richard and Quéré2000), and is not always negligible when the kinetic energy is small for
$\textit{We}\lt 1$
. A drop can be elongated at detachment when
$\textit{We}\lt 1$
, as shown in the images inset in figure 2(
$c$
). Other authors defined
$E_{ {r}}$
as the gravitational potential energy of the drop at its maximum height after rebound, but this definition ignores possible viscous loss due to the surrounding air during free flight (Wang et al. Reference Wang, Zhao, Sun, Mehrizi, Lin, Guo and Chen2022; Thenarianto et al. Reference Thenarianto, Koh, Lin, Jokinen and Daniel2023). Our definition does not make any of these assumptions, providing energetically consistent and accurate measurements at low
$\textit{We}$
.
The uncertainty of each measurement was determined by considering the spatial and temporal resolution of our experimental set-up. The temporal resolution is related to the time between successive frames
$\Delta t = 0.5 \, (\textrm {fps} )^{-1}$
and the spatial resolution is calculated using the width of a pixel
$L_p$
, as
$\Delta L = 0.5 \, ( L_{p} ) ^{-1}$
. To estimate the uncertainty in the velocity, we first calculated the residuals between the measured centre of mass height
$h$
and the fitted parabola. The standard deviation of these residuals
$\sigma _{{res}}$
served as a measure of the uncertainty in the position data. This uncertainty is propagated to the velocity estimate by assuming that the uncertainty in the velocity is equal to the standard deviation of the residuals, scaled by the curvature of the fitted parabola,
$\delta V \approx 0.5 \, \sigma _{ {res}} \left | 2 a \right |$
, where
$a$
is the quadratic coefficient of the fitted parabola. Uncertainty propagation was performed using Mathematica, where first-order Taylor expansion was used for uncertainty propagation and uncertainties were assumed uncorrelated.
3. Modelling and simulations
3.1. Kinematic match model formulation
We consider a drop of incompressible fluid of constant density
$\rho$
, dynamic viscosity
$\mu$
and surface tension coefficient
$\sigma$
impacting normally on a perfectly hydrophobic (non-wetting), horizontal solid substrate. Before impact, the droplet is spherical with radius
$R$
and the surface of the droplet is stationary relative to its centre of mass. During impact the drop is ‘pressed’ to the substrate rather than in true physical ‘contact’ with it; however, we use the latter in our discussion for consistency and connection with the prior literature. The criteria for contact and its implication for the model results are discussed in Appendix A.
We introduce spherical coordinates
$(\zeta ',\theta ,\phi )$
with
$\theta \in [0,\pi ]$
and
$\phi \in [0,2\pi ]$
(where
$'$
indicates a dimensional variable) to describe the free surface of the drop. Here,
$\theta$
is measured relative to the direction of gravity, and the origin is at the centre of mass of the drop. Due to the assumed axisymmetry, the azimuthal angle
$\phi$
can be ignored, consistent with experiments for all cases considered here. Additionally, we define the associated Cartesian system of coordinates
$(x',z')$
, as shown in figure 3(
$a$
) for an impacting drop. This formulation allows the smooth drop shape to be discretised, as shown in figure 3(
$b$
), and permits numerical approximations that enable efficient simulation of drop rebound across a wide range of
$\textit{We}$
.
(
$a$
) Schematic of an axisymmetric drop impacting a rigid substrate, where
$\zeta ' ( \theta , t )$
is the dimensional interface,
$h'( t )$
the height of the centre of mass of the drop and
$S ( t )$
the contact surface. (
$b$
) Angular discretisation of the right hemisphere of an axisymmetric drop during impact. Filled circles represent mesh points where the drop and substrate touch, while open circles denote points along the liquid–gas interface. (
$c$
) An example of shape reconstruction through Legendre synthesis for modes
$l=2{-}4$
.
3.1.1. Drop motion
The drop is subject to gravitational acceleration
$g \hat {\boldsymbol{z}}$
, and the substrate is at a distance
$h'$
from its centre of mass. At time
$t'=0$
,
$h'=R$
and the centre of mass of the drop moves towards the substrate with velocity
$V$
. The fluid motion within the drop is governed by the Navier–Stokes equations for incompressible flow, with standard kinematic and dynamic boundary conditions. The normal impact on the substrate is modelled by an evolving pressure distribution over the contact area
$S ( t )$
, where the drop surface and substrate coincide in the model.
Following Alventosa et al. (Reference Alventosa, Cimpeanu and Harris2023), we linearise the governing equations and represent the surface deformation of the drop
$\eta '$
as a sum of axisymmetric spherical harmonic modes:
\begin{equation} \eta ' := \zeta '-R = \sum \limits _{l=2}^\infty {\mathscr{A}^{\prime}_lP_l(\cos (\theta ))}. \end{equation}
Here,
$P_l$
are Legendre polynomials and
$\mathscr{A}^{\prime}_l$
denotes the time-dependent amplitude of the
$l{\text{th}}$
mode. The
$l=0$
mode is excluded due to incompressibility. The
$l=1$
mode is also excluded in the present formulation due to the choice of the drop’s centre of mass as the origin, since the
$l=1$
mode shifts the centre of mass. Figure 3(
$c$
) demonstrates the synthesis of Legendre polynomials for modes
$l=2{-}4$
into a sample deformed drop shape. Similarly, the contact pressure
$p'$
exerted by the substrate on the drop is expressed as
\begin{equation} p' = \sum \limits _{l=0}^\infty {\mathscr{B}^{\prime}_lP_l(\cos (\theta ))}, \end{equation}
where
$\mathscr{B}^{\prime}_l$
is the time-dependent amplitude of the
$l{\text{th}}$
Legendre mode.
The evolution of the drop interface is obtained by applying linearised kinematic conditions on the surface in the weakly damped and small-amplitude regime (Lamb Reference Lamb1924; Miller & Scriven Reference Miller and Scriven1968; Phillips & Milewski Reference Phillips and Milewski2024):
\begin{align}&\qquad\qquad\qquad\qquad\qquad\quad \dot {\mathscr{A}}_l^{\prime} = \mathscr{U}_l^{\prime}, \quad \text{ for }\ l = 2,3,\ldots , \infty ; \end{align}
\begin{align}& \dot {\mathscr{U}}_l^{\prime} = -\frac {\sigma l \left (l+2\right )\left (l-1\right )}{\rho R^3}\mathscr{A}_l^{\prime} -2(2l+1)(l-1)\frac {\mu }{\rho R^2}\mathscr{U}_l^{\prime} -\frac {l}{\rho R}\mathscr{B}_l^{\prime}, \quad \text{ for } l = 2,\ldots ,\infty , \end{align}
where
$\mathscr{U}_l^{\prime}$
is the amplitude of the
$l{\text{th}}$
Legendre mode of the surface’s velocity,
$\dot {\zeta '}$
. Equations (3.3) and (3.4) are subject to
\begin{align}& \mathscr{A}_l^{\prime}(t=0) = 0, \quad \text{ for }\ l = 2,\ldots ,\infty ; \end{align}
\begin{align}& \mathscr{U}_l^{\prime}(t=0) = 0, \quad \text{ for }\ l = 2,\ldots ,\infty . \end{align}
The height of the centre of mass is naturally governed by Newton’s second law of motion:
\begin{align}&\qquad\quad \dot {h}' = v', \end{align}
\begin{align}& \dot {v}' = -g + \frac {1}{m}\frac {4\pi R^2}{3}\mathscr{B}_1^{\prime}, \end{align}
subject to
\begin{align}&\,\, h'(t=0) = R, \end{align}
\begin{align}& v'(t=0) = -V, \end{align}
where
$m=( {4\pi }/{3})\rho R^3$
is the mass of the drop and
$g$
is gravitational acceleration.
3.1.2. The kinematic match
The motion of the drop’s surface must satisfy additional compatibility conditions. Namely, when the substrate exerts pressure on the drop, a circular contact area
$\theta \leqslant \theta _c$
forms, which must lie on the
$z'=0$
plane. That is,
\begin{equation} \sum \limits _{l=2}^\infty {\mathscr{A}_l^{\prime}P_l(\cos (\theta ))} = \frac {h'}{\cos (\theta )}-R, \qquad \text{ for } \theta \leqslant \theta _c. \end{equation}
Outside the contact area, the drop surface must be strictly above the solid surface:
\begin{equation} \sum \limits _{l=2}^\infty {\mathscr{A}_l^{\prime}P_l(\cos (\theta ))} \gt \frac {h'}{\cos (\theta )}-R, \qquad \text{ for } \theta \gt \theta _c, \end{equation}
and the contact pressure there must be null, that is,
\begin{equation} \sum \limits _{l=0}^\infty {\mathscr{B}_l^{\prime}P_l(\cos (\theta ))} = 0, \qquad \text{ for } \theta \gt \theta _c. \end{equation}
Finally, at the boundary of the contact area, the drop surface must be tangent to the substrate:
\begin{equation} \frac {\partial }{\partial \theta }\left [\left (R+\sum \limits _{l=2}^\infty {\mathscr{A}_l^{\prime}P_l(\cos (\theta ))}\right )\cos (\theta )\right ]_{\theta =\theta _c} = 0, \end{equation}
which captures the non-wetting nature of the surface. Our formulation results in a fully self-contained model without fitting parameters, which we can now non-dimensionlise.
3.1.3. Non-dimensionalisation
The characteristic length, mass and time are taken to be the undeformed droplet radius
$R$
, mass scale
$\rho R^{3}$
and the inertio-capillary time
$t_{\sigma }= ( ( \rho R^{3} ) / \sigma )^{1/2}$
, respectively. We also define the Bond number
$Bo$
, Weber number
$\textit{We}$
and Ohnesorge number
$Oh$
as in (1.1). This process leads to the dimensionless equations:
\begin{align}&\qquad\qquad\qquad\qquad\qquad\qquad \dot {\mathscr{A}}_l = \mathscr{U}_l, \quad \text{ for } l = 2,\ldots , \infty ; \end{align}
\begin{align}& \dot {\mathscr{U}}_l = -l \left (l+2\right )\left (l-1\right )\mathscr{A}_l -2(2l+1)(l-1) Oh \, \mathscr{U}_l -l \mathscr{B}_l, \qquad \text{ for } l = 2,\ldots , \infty ; \end{align}
\begin{align}&\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad \dot {h} = v; \end{align}
\begin{align}&\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad \dot {v} = - Bo + \mathscr{B}_1; \end{align}
which are subject to
\begin{align}&\qquad \mathscr{A}_l(t=0) = 0, \quad \text{ for } l = 2,\ldots , \infty ; \end{align}
\begin{align}&\qquad \mathscr{U}_l(t=0) = 0, \quad \text{ for } l = 2,\ldots , \infty ; \end{align}
\begin{align}&\qquad\qquad\qquad\quad\,\, h(t=0) = 1; \end{align}
\begin{align}&\qquad\qquad\qquad v(t=0) =- \sqrt {We}; \end{align}
\begin{align}& \sum \limits _{l=2}^\infty {\mathscr{A}_lP_l(\cos (\theta ))} = \frac {h}{\cos (\theta )}-1, \quad \text{ for } \theta \leqslant \theta _c; \end{align}
\begin{align}& \sum \limits _{l=2}^\infty {\mathscr{A}_lP_l(\cos (\theta ))} \gt \frac {h}{\cos (\theta )}-1, \quad \text{ for } \theta \gt \theta _c; \end{align}
\begin{align}&\qquad \sum \limits _{l=0}^\infty {\mathscr{B}_lP_l(\cos (\theta ))} = 0, \quad \text{ for } \theta \gt \theta _c; \end{align}
\begin{align}& \frac {\partial }{\partial \theta }\left [\left (1+\sum \limits _{l=2}^\infty {\mathscr{A}_lP_l(\cos (\theta ))}\right )\cos (\theta )\right ]_{\theta =\theta _c} = 0. \end{align}
We solve this system of equations by first applying numerical approximations, detailed in Appendix B. Briefly, these include truncating the expansions of the surface shape, velocity and pressure distribution with a spectral mesh consisting of Legendre polynomials up to degree
$L=90$
and approximating time derivatives using a first-order implicit Euler method. All simulations used to generate the results in the remaining figures of this paper adopt this value of
$L$
. A schematic of the angular mesh is shown in figure 3(
$b$
). The result of these approximations is a discrete optimisation problem, which we solve using an adaptive timestep and error minimisation scheme, as described in Appendix C. This approach accurately solves the system at reduced computation cost relative to DNS, with the longest simulations completing in
$\mathcal{O}(10)$
CPU minutes and a typical simulation time of less than 1 CPU minute.
3.2. Direct numerical simulation
We develop a computational framework based on the Basilisk open-source code structure (Popinet Reference Popinet2009, Reference Popinet2015) as a comparative benchmark for the mathematical model described in § 3.1, and to interrogate quantities difficult to visualise experimentally. The fully nonlinear implementation is used to resolve the time-dependent Navier–Stokes equations in both liquid and gas phases, with physical properties of the latter taken to be the values for air at room temperature, i.e.
$\rho _g = 0.00121$
g
$ \,$
cm
$^{-3}$
and
$\mu _g = 0.000181$
g
$ \,$
cm
$^{-1}$
$ \,$
s
$^{-1}$
, in an attempt to replicate the experimental conditions. Basilisk employs the volume-of-fluid method (Hirt & Nichols Reference Hirt and Nichols1981; Scardovelli & Zaleski Reference Scardovelli and Zaleski1999) for the implicit representation of the interface describing the transition between the different fluid regions, with a so-called colour function
$f$
taking the value
$1$
inside the liquid and
$0$
inside the gas regions and values
$0 \lt f \lt 1$
representing computational cells cut by a fluid–fluid interface. The adaptive mesh refinement and parallelisation capabilities of the implementation make it a particularly good candidate for this multiscale impact problem. The code base for bouncing problems has been constructed and systematically validated for several years in our extended group, covering fluid–structure interaction (Galeano-Rios et al. Reference Galeano-Rios, Cimpeanu, Bauman, MacEwen, Milewski and Harris2021) and liquid–liquid impact (Alventosa et al. Reference Alventosa, Cimpeanu and Harris2023), with comprehensive numerical campaigns in the wider community (Sanjay et al. Reference Sanjay, Chantelot and Lohse2023a
,
Reference Sanjay, Lakshman, Chantelot, Snoeijer and Lohseb
) providing key insights into related regimes.
(
$a$
) Visualisation of the computational domain comprising liquid (blue) and gas (white) regions, and underlying adaptive grid cell structure in the DNS, for a typical
$\textit{We}=0.253$
impact taken at the point of maximal radius deformation
$t=0.746$
ms. (
$b$
) Dimensionless distance above the substrate
$\bar {z} = z / R$
against time
$\bar {t} = t / t_{\sigma }$
for the top, bottom and centre of mass of a drop for
$\textit{We}= 0.253$
,
$Bo = 0.02$
and
$Oh = 0.03$
. Black solid lines represent experiment, while blue solid lines and yellow dashed lines represent the KM model and DNS, respectively. (
$c$
) Overlay of the predicted shape from the KM model (solid blue) and DNS (dashed yellow) with experimental images corresponding to (
$b$
). The drop radius is
$R=0.203$
mm and the time interval between images is 0.375 ms.
The axisymmetric computational domain set up in cylindrical coordinates measures
$8R$
across both (
$r$
and
$z$
) spatial dimensions. To ensure sufficient resources are dedicated to the delicate bouncing dynamics, we use a strong grid cell refinement level with a minimal grid cell size of approximately
$1/20$
μm, verified to be robust across our parameter regime of interest. This stringent resolution specification was particularly relevant at the lowest Weber numbers studied herein, where the implementation requires most resilience in view of increased numerical stiffness. Adaptivity was tailored around the evolving location of the drop interface and changes in magnitude of the velocity components in the flow, as illustrated in figure 4(
$a$
). Special care was taken during runtime via stringent projection solver tolerances and increased number of internal iterations set for the multigrid cycle to ensure a high degree of accuracy. Furthermore, a zero Dirichlet boundary condition for the colour function
$f$
was prescribed at the level of the solid surface to prevent touchdown and contact, with all tests resulting in the drop rebounding or resting on the surface. We employ a criterion based on the drop interface crossing a threshold level of
$0.02R$
above the substrate to define the initialisation and end of contact as discussed in Appendix A. This criterion enables us to build a consistent comparative framework with the model and experiment. Under these specifications, a typical run is characterised by
$\mathcal{O}(10^5)$
grid cells and
$\mathcal{O}(10{-}100)$
CPU hours of runtime, executed across multiple cores on local high-performance computing architectures. Mesh-independent results in our target metrics have been observed under the conditions above and reported in the sections to follow. All associated code, including set-up scripts and post-processing functionality, are provided in the form of a GitHub repository at https://github.com/rcsc-group/BouncingDropletsSolidSurface.
3.3. Validation
Beyond verification and resolution studies across our parameter space of interest, we also validate our theoretical approaches with experimental data. A typical case illustrating a comparison between the KM model and DNS is the rebound of a drop with radius
$R=0.203$
mm for
$\textit{We}=0.253$
,
$Bo=0.02$
and
$Oh=0.03$
. Figure 4(
$b$
) plots the normalised distance from the substrate
$\bar {z}=z/R$
of the top, bottom and centre of mass of the drop against time
$\bar {t}=t/t_{\sigma }$
. The experimental trajectory is shown in black and aligns well with predictions from the KM model (blue solid line) and DNS (yellow dashed line) for all three locations. Figure 4(
$c$
) shows overlays of the predicted drop shape from the KM model (solid blue line) and DNS (yellow dashed line) onto experimental images, demonstrating strong visual agreement. Supplementary movie 3 shows this rebound event with overlays. Additionally, we note that the accuracy of the KM model improves on a relative basis as
$\textit{We}$
decreases, as demonstrated in supplementary movie 4 for a similar drop at
$\textit{We}=0.023$
. This low-
$\textit{We}$
regime is of high consistency with the underlying simplifying assumptions, one in which resources can also be deployed very efficiently given smaller interface deformations. By contrast, we find that the DNS framework struggles for sufficiently low values of
$\textit{We}$
(which we identified to be of
$\mathcal{O}(10^{-2})$
or lower) without significant customisation or computational redesign, as both numerical stiffness and parasitic capillary currents become progressively more difficult to contain accurately. This makes the cross-validation and interplay with the KM model even more constructive as our investigation reaches distinguished limits of our target parameter regime.
4. Results
We performed experiments over a wide range of
$\textit{We}$
,
$Oh$
and
$Bo$
, as highlighted in supplementary movies 5, 6 and 7. In this section, we compare our results with DNS data and the KM model. We first investigate the metrics defining rebound, the coefficient of restitution
$\varepsilon$
and contact time
$t_{c}$
. Next, we analyse the evolution of the drop shape and time scales associated with spreading and retraction. Last, we show how the maximum deformation of the drop changes over our parameter space and compare with the KM model and energetic arguments.
4.1. Rebound metrics
(
$a$
) Coefficient of restitution
$\varepsilon$
against Weber number
$\textit{We}$
for drops with varying
$Oh$
. Circles denote experiments and stars represent DNS results for
$Bo=0.02$
, a typical value in experiments. Solid and dashed lines are the KM model for
$Bo=0.02$
and
$Bo=0$
, where the black dashed line shows the inertio-capillary limit (
$Oh=0$
;
$Bo = 0$
). (
$b$
) Dimensionless contact time
$\bar {t}_c = t_{c} / t_{\sigma }$
against
$\textit{We}$
for experiments, DNS and model predictions from (
$a$
). The colour bar maps
$Oh$
to colour on a logarithmic scale.
The coefficient of restitution
$\varepsilon$
quantifies the energy recovered during rebound, and behaves differently at low and high
$\textit{We}$
. Figure 5(
$a$
) shows
$\varepsilon$
against
$\textit{We}$
, for various
$Oh$
spanning two orders of magnitude (as denoted by colour). Here, circle markers represent experiments and star markers represent DNS results for
$Bo=0.02$
, the most common value in our low-
$\textit{We}$
experiments. The solid and dashed lines are predictions from the KM model for
$Bo=0.02$
and
$Bo=0$
, respectively. The intercept of a curve with the
$x$
axis signifies the cessation of rebound. The KM model predictions for the experimental trends at low
$\textit{We}$
agree reasonably well with experiment and produce robust predictions even at very low
$\textit{We}$
where the DNS struggles. At high
$\textit{We}$
, the DNS predictions agree better with experiment, where the assumptions of the kinematic model are less justified. The complementary modelling approaches allow us to confidently span the full region accessed by experiment, and beyond.
For
$\textit{We} \gtrapprox 0.5$
,
$\varepsilon$
decreases as
$\textit{We}$
increases for each
$Oh$
in agreement with prior literature (Biance et al. Reference Biance, Chevy, Clanet, Lagubeau and Quéré2006; Jha et al. Reference Jha, Chantelot, Clanet and Quéré2020). In addition, the black dashed line shows the prediction of the KM model in the strict inertio-capillary limit (
$Oh = 0$
;
$Bo = 0$
). Even in the complete absence of gravity and viscosity, our model predicts a similar roll-off at high
$\textit{We}$
due to energy transfer to droplet oscillations not accounted for in the definition of
$\varepsilon$
. Curiously, this high-
$\textit{We}$
roll-off is notably absent for the analogous case of droplet–bath rebounds, where for low
$Oh$
and
$Bo$
, the coefficient of restitution levels off at approximately
$0.3$
(Zhao, Brunsvold & Munkejord Reference Zhao, Brunsvold and Munkejord2011; Wu et al. Reference Wu, Hao, Lu, Xu, Hu and Floryan2020; Alventosa et al. Reference Alventosa, Cimpeanu and Harris2023).
Below
$\textit{We} \approx 0.5$
,
$\varepsilon$
becomes nearly
$\textit{We}$
-independent across roughly two orders of magnitude when
$Oh$
is small. As
$Oh$
increases, the influence of
$\textit{We}$
on
$\varepsilon$
increases. As
$\textit{We}$
decreases,
$Bo$
effects cause
$\varepsilon \rightarrow 0$
at a critical Weber number
$\textit{We}_c$
that increases with
$Oh$
. For the experimental data presented,
$\varepsilon = 0$
represents the highest-
$\textit{We}$
experiment where rebound was not detected (for each viscosity). In related work for low-
$\textit{We}$
rebound on hydrophobic surfaces, a similar roll-off leading to bounce cessation was attributed to contact line friction (Wang et al. Reference Wang, Zhao, Sun, Mehrizi, Lin, Guo and Chen2022; Thenarianto et al. Reference Thenarianto, Koh, Lin, Jokinen and Daniel2023), wholly absent in the current experiments and models. To investigate the physical origin of this threshold in our fully non-wetting scenario, we performed additional KM simulations in the absence of gravity (i.e.
$Bo=0$
, shown as the coloured dashed curves). For the parameters considered here,
$\textit{We}_c$
is eliminated in the absence of gravity, effectively extending bouncing to lower
$\textit{We}$
. Gravity shifts the static equilibrium height of the drop’s centre of mass below one radius above the surface, representing a higher barrier for complete rebound when energy is lost to viscosity or transferred to droplet oscillations during contact. In the low-
$\textit{We}$
regime, the KM model predicts nearly perfect elastic collisions (
$\varepsilon \approx 1$
) in the inertio-capillary limit (
$Oh=0$
;
$Bo=0$
).
The contact time
$t_{c}$
follows the simple scaling
$t_{c} \sim t_{\sigma }$
for moderate and high
$\textit{We}$
(Richard et al. Reference Richard, Clanet and Quéré2002) but depends on
$\textit{We}$
and
$Oh$
at low
$\textit{We}$
. Figure 5(
$b$
) shows the non-dimensional contact time
$\bar {t}_{c} = t_{c} / t_{\sigma }$
against
$\textit{We}$
for drops of varying
$Oh$
. The results correspond with those in figure 5(
$a$
); thus, plot markers, line types and colours hold the same meaning. For all
$Oh$
, the contact time approaches a constant as
$\textit{We}$
increases, in agreement with prior literature (Okumura et al. Reference Okumura, Chevy, Richard, Quéré and Clanet2003). As
$\textit{We}$
decreases,
$t_{c}$
increases at an
$Oh$
-dependent rate; the rate of increase is greater for larger
$Oh$
. Thus, the effects of viscosity on contact time are most influential at low
$\textit{We}$
and may alter the utility of drops in a number of applications involving heat transfer (Shiri & Bird Reference Shiri and Bird2017), for instance. The KM model agrees well here for all
$\textit{We}$
and
$Oh$
, and the inertio-capillary limit sets a lower bound on the contact time and an upper bound on the coefficient of restitution for each
$\textit{We}$
, similar to the case for droplet–bath rebound (Alventosa et al. Reference Alventosa, Cimpeanu and Harris2023).
4.2. Spreading and retraction
(
$a$
) Dimensionless equatorial radius
$\bar {r}_{m} = r_{m} / R$
against dimensionless time
$\bar {t}=t / t_{\sigma }$
for experiments (markers), DNS (dotted lines) and the KM model (solid lines) for
$\textit{We} = 2.269$
(light yellow),
$\textit{We} = 0.253$
(green) and
$\textit{We} = 0.023$
(dark blue), where
$Bo=0.02$
and
$Oh=0.03$
for all curves. (
$b$
) Dimensionless contact radius
$\bar {r}_{c} = r_{c} / R$
against dimensionless time
$\bar {t}=t / t_{\sigma }$
, where line, marker type and colour are the same as (
$a$
). The insets show experimental results for drops of similar
$\textit{We}$
(
$2.5{-}3.5$
) and
$Bo$
(
$0.16{-}0.24$
) for
$Oh=$
0.03, 0.08, 0.29 and 0.7, where marker darkness indicates higher
$Oh$
.
The shift in rebound behaviour between low and high
$\textit{We}$
revealed in
$\S$
4.1 coincides with a shift in the spreading and retraction dynamics. Figure 6(
$a$
) shows the non-dimensional equatorial radius
$\bar {r}_{m} = r_{m} / R$
against non-dimensional time
$\bar {t} = t / t_{\sigma }$
during contact for
$\textit{We}=0.023$
(dark blue),
$\textit{We}=0.253$
(green) and
$\textit{We}=2.269$
(light yellow), with
$Bo=0.02$
and
$Oh = 0.03$
. Circles represent experimental measurements, while dotted and solid lines denote DNS and the KM model results, respectively. In all cases, inertia drives spreading, causing
$\bar {r}_{m}$
to increase after an initial period of localised deformation where
$\bar {r}_{m} \approx 1$
. As
$\textit{We}$
increases, the drop expands quicker and
$\bar {r}_{m}$
reaches a higher maximum. After the drop reaches its maximum extent, surface tension drives retraction, decreasing
$\bar {r}_{m}$
. For
$\textit{We} = 0.023$
, retraction mirrors spreading, and
$\bar {r}_{m}$
evolves symmetrically, consistent with quasi-static-like behaviour (Moláček & Bush Reference Moláček and Bush2012; Chevy et al. Reference Chevy, Chepelianskii, Quéré and Raphaël2012). When
$\textit{We} = 2.269$
, the two phases are dissimilar, and
$\bar {r}_{m}$
overshoots
$\bar {r}_{m} = 1$
during retraction, leading to apparent underdamped oscillations. The inset plot shows the evolution of
$\bar {r}_{m}$
for drops of similar
$\textit{We}$
and
$Bo$
but varying
$Oh$
, where darker marker shade indicates higher
$Oh$
. Increasing viscous dissipation can eliminate the appearance of underdamped behaviour; however, it is not sufficient to recover the symmetry between spreading and retraction. Instead, higher
$Oh$
shortens the spreading phase, reduces the maximum
$\bar {r}_{m}$
and lengthens the retraction phase, thereby increasing asymmetry in the spreading and retraction process.
The increased asymmetry between spreading and retraction as
$\textit{We}$
increases is also observed in the non-dimensional contact radius
$\bar {r}_{c} = r_{c} / R$
. Figure 6(
$b$
) shows
$\bar {r}_{c} = r_{c} / R$
against non-dimensional time
$\bar {t} = t / t_{\sigma }$
for
$\textit{We}=0.023$
(dark blue),
$\textit{We}=0.253$
(green) and
$\textit{We}=2.269$
(light yellow). The data and model predictions correspond with those in figure 6(
$a$
); thus, plot markers, line types and colours hold the same meaning. Unlike
$\bar {r}_{m}$
, there is no delay in the response of
$\bar {r}_{c}$
to impact as
$\bar {t} = 0$
is defined as the moment the drop and substrate ‘contact’. As
$\textit{We}$
increases, the maximum
$\bar {r}_c$
increases, while the time to reach it decreases. For
$\textit{We} = 0.023$
, the spreading and retraction of
$\bar {r}_c$
are symmetric, whereas at
$\textit{We} = 2.269$
, they are asymmetric, qualitatively resembling the behaviour of
$\bar {r}_m$
during rebound. The inset shows the evolution of
$\bar {r}_c$
for drops of similar
$Bo$
and
$\textit{We}$
, but increasing
$Oh$
, as indicated by darker-shaded markers. Similar to
$\bar {r}_m$
, increasing
$Oh$
shortens the spreading phase by reducing the maximum
$\bar {r}_c$
, and extending the retraction phase. The asymmetry between spreading and retraction motivates a detailed look at their relationship across all experiments.
The spreading time
$t_s$
is the time from when the drop first contacts the substrate until
$r_c$
reaches its maximum value
$R_c$
, and is proportional to the retraction time
$t_r$
at low
$\textit{We}$
. In figure 7 we show the relative spreading time
$t_{s} / t_{c}$
, where the contact time is the sum of the spreading and retraction time
$t_{c} = t_{s} + t_{r}$
, against
$\textit{We}$
for all
$Oh$
. Solid lines represent predictions from the KM model for
$Bo=0.02$
, with the black dashed curve representing the inertio-capillary limit (
$Bo=0$
;
$Oh=0$
). When
$\textit{We} \gtrapprox 0.5$
, drops spread faster than they retract, and increasing
$\textit{We}$
increases the difference between
$t_s$
and
$t_r$
. However, the time the drop spends spreading relative to retracting is nearly
$\textit{We}$
-independent for
$\textit{We} \lessapprox 0.5$
, until the drop ceases to bounce. This trend holds for all viscosities, with
$Oh$
limiting the maximum of
$t_{s} / t_{c}$
. Although this maximum approaches 0.5 as
$\textit{We} \rightarrow 0$
in the inertio-capillary limit, it surpasses it and peaks near
$\textit{We} = 0.04$
. The inset in figure 7 shows the evolution of the dimensionless contact radius
$\bar {r}_c$
near its maximum, which is denoted by a black circle. Oscillations due to capillary waves cause an apparent jumpiness in
$t_s$
at low
$Oh$
. In particular, as
$\textit{We}$
decreases from
$\textit{We} = 0.05$
to
$\textit{We}= 0.04$
, the maximum contact radius suddenly switches between its two largest peaks, producing the apparent sudden rise in
$t_s$
. These results reveal that the symmetry of the spreading and retraction process coincides with nearly elastic rebound metrics that are independent of
$\textit{We}$
in the inertio-capillary limit. In contrast, the emergence of asymmetry coincides with the onset of
$\textit{We}$
-dependent trends for
$\varepsilon$
while increasing
$Oh$
similarly decreases
$t_{s}/t_{c}$
and
$\varepsilon$
across all
$\textit{We}$
. As in the prior rebound measurements, the KM model performs best for low
$\textit{We}$
whereas DNS excels at higher
$\textit{We}$
for the spreading and retraction metrics presented throughout this section.
The relative spreading time
$t_{s} / t_{c}$
against Weber number
$\textit{We}$
for drops with varying
$Oh$
. Circles represent experiment and stars denote DNS for
$Bo=0.02$
. Solid lines are predictions from the KM model for
$Bo = 0.02$
, and the black dashed line shows the inertio-capillary limit (
$Oh = 0$
;
$Bo = 0$
). The inset shows the predicted dimensionless contact radius
$\bar {r}_c$
against dimensionless time
$\bar {t} = t / t_{\sigma }$
for
$\textit{We}=0.04$
and
$\textit{We}=0.05$
, where
$Bo = 0$
and
$Oh = 0$
. The black circles show when
$\bar {r}_c$
is maximum.
4.3. Maximum spreading and deformation
We next examine how drop deformation depends on the governing dimensionless parameters and relate deformation and spreading. To do this, we define the maximum equatorial deformation
$\beta - 1$
, where
$\beta = R_{m} / R$
is the non-dimensional maximum equatorial radius and
$R_{m}=\textrm {max}(r_{m}(t))$
. Similarly, we define maximum vertical deformation
$1 - \alpha$
, where
$\alpha = H_{\textit{min}} / 2R$
and
$H_{\textit{min}} = \textrm {min}(H(t))$
. Lastly, we define the maximum spreading radius as
$\beta _{c} = R_{c} / R$
, where
$R_{c} = \textrm {max}(r_{c}(t))$
.
4.3.1. Energy argument
In addition to the KM model, we can compare the drop spreading and deformation with a simpler energetic argument by assuming the drop principally deforms in the
$l=2$
mode in the low-
$\textit{We}$
regime. Under this assumption, the instantaneous drop shape is fully defined by the coefficient
$\mathscr{A}_2^{\prime }$
, with
$r_{m}=R-({1}/{2})\mathscr{A}_2^{\prime }$
and
$h=2R+2\mathscr{A}_2^{\prime }$
.
The anomalous surface energy associated with drop deformations of the
$l$
th mode (for
$l\geqslant 2$
) can be expressed as
\begin{equation} E_l=2 \pi \sigma \frac {(l+2)(l-1)}{2 l+1} \mathscr{A}_l^{\prime 2}, \end{equation}
with the surface energy in the
$l=2$
mode being
$E_2=({8}/{5}) \pi \sigma \mathscr{A}_{2}^{\prime 2}$
(Moláček & Bush Reference Moláček and Bush2012). Neglecting gravity (
$Bo^{2} \ll We$
) and viscous dissipation (
$Oh\ll 1$
), and assuming the initial kinetic energy of the droplet (
$E_V=({2}/{3})\pi \rho R^3 V^2$
) is fully converted into surface energy of the
$l=2$
mode (i.e.
$E_V=E_2$
), one finds
\begin{equation} \beta -1=\sqrt {\frac {5}{48}} \textit{We}^{1 / 2} \end{equation}
and
\begin{equation} 1-\alpha =\sqrt {\frac {5}{12}} \textit{We}^{1 / 2}, \end{equation}
using the relations
$\mathscr{A}_{2}^{\prime }=2R (1-\beta )$
and
$\mathscr{A}_{2}^{\prime }=R (\alpha -1 )$
, respectively, and assuming
$\mathscr{A}_{2}^{\prime }\leqslant 0$
(i.e.
$R_m\geqslant R$
,
$H_{\textit{min}}\leqslant 2R$
) in the deformed state. The result for
$\beta - 1$
(4.2) is identical to that derived by Richard & Quéré (Reference Richard and Quéré2000) using similar arguments. This derivation also suggests a geometric relationship between
$\beta$
and
$\alpha$
:
\begin{equation} 1-\alpha =2(\beta -1). \end{equation}
To connect this deformation to the contact radius
$r_c$
, we appeal to a quasi-static argument expected to be valid for
$\textit{We}\ll 1$
(Chevy et al. Reference Chevy, Chepelianskii, Quéré and Raphaël2012; Moláček & Bush Reference Moláček and Bush2012). The force associated with the creation of additional surface area in the
$l=2$
mode can be expressed as the derivative of the surface energy:
$F_2=- ({\textrm{d}E_2}/{\textrm{d}\mathscr{A}_2^{\prime }})=-({16}/{5})\pi \sigma \mathscr{A}_2^{\prime }$
. Following Chevy et al. (Reference Chevy, Chepelianskii, Quéré and Raphaël2012), we assume small deformations such that the capillary pressure inside the drop is approximately
$2\sigma /R$
, and thus the net vertical force associated with the ‘contact’ region is
$F_{\sigma }=({2\sigma }/{R})\pi r_c^2$
. Balancing these forces, one arrives at the relationship
${r_c^2}/{R^2}=-({8}/{5}) ({\mathscr{A}_2^{\prime }}/{R})$
. This result can be combined with the prior relations to arrive at the prediction
\begin{equation} \beta _c^2=\frac {16}{5}(\beta -1) \end{equation}
or
\begin{equation} \beta _c=\left (\frac {16}{15}\right )^{1/4} \textit{We}^{1 / 4}. \end{equation}
Note that (4.2)–(4.6) have limited predictive power compared with the KM model, which captures viscous effects and the complete evolution of the drop shape. Nevertheless, the favourable comparison between the simple energy argument and KM model elucidates the dominant mode of deformation and motivates the parametric scalings.
4.3.2. Comparisons
Figure 8(
$a$
) shows the maximum equatorial deformation
$\beta - 1$
against Weber number
$\textit{We}$
for all experiments (circles) and DNS (stars). Predictions from the KM model are shown as solid lines for
$Bo=0.02$
, with the black dashed curve showing the inertio-capillary limit (
$Oh=0$
;
$Bo=0$
). The maximum deformation
$\beta - 1$
monotonically increases with
$\textit{We}$
and is well captured by the KM model. The red dot-dashed curve represents (4.2), which gives an excellent prediction for our low-viscosity (
$Oh\lt 0.05$
) data with no fitting parameters. This agreement persists across four orders of magnitude of
$\textit{We}$
, as shown more clearly in the inset with our low-viscosity results and data from prior literature (Okumura et al. Reference Okumura, Chevy, Richard, Quéré and Clanet2003; Clanet et al. Reference Clanet, Béguin, Richard and Quéré2004; Biance et al. Reference Biance, Chevy, Clanet, Lagubeau and Quéré2006; Laan et al. Reference Laan, de Bruin, Bartolo, Josserand and Bonn2014). Beyond
$\textit{We} \approx 10^{2}$
, (4.2) no longer captures the trend of the amalgamated data, which follow a trend more closely resembling
$\beta - 1 \sim \textit{We}^{1/4}$
.
(
$a$
) Maximum equatorial deformation
$\beta - 1$
(
$\beta = R_{m} / R$
) and (
$b$
) maximum vertical deformation
$1-\alpha$
(
$\alpha = H_{\textit{min}} / 2 R$
) plotted against Weber number
$\textit{We}$
. Circles represent experiments and stars denote DNS for
$Bo=0.02$
. Solid curves show the KM model with
$Bo = 0.02$
, and the inertio-capillary limit (
$Oh=0$
;
$Bo=0$
) is represented by black dashed curves. The red dot-dashed curves are predictions from our simplified model (§ 4.3.1). The inset in (
$a$
) shows our low-viscosity data
$Oh \lt 0.05$
(coloured circles) with prior results (diamond (Okumura et al. Reference Okumura, Chevy, Richard, Quéré and Clanet2003), pentagon (Clanet et al. Reference Clanet, Béguin, Richard and Quéré2004), square (Biance et al. Reference Biance, Chevy, Clanet, Lagubeau and Quéré2006), triangle (Laan et al. Reference Laan, de Bruin, Bartolo, Josserand and Bonn2014)). The red dot-dashed curve represents (4.2), and the black solid line with slope 1/4 is used to guide the eye. The inset in (
$b$
) plots
$1-\alpha$
against
$\beta - 1$
for all experiments and the red dot-dashed curve represents (4.4). (
$c$
) Dimensionless spreading area
$\beta _{c}^{2}$
against
$\beta - 1$
for all experiments and DNS results, with (4.5) shown as a red dot-dashed curve. The colour bar maps
$Oh$
to colour on a logarithmic scale.
The maximum vertical deformation
$1-\alpha$
can also be accurately predicted by our KM model and energy arguments. Figure 8(
$b$
) shows
$1- \alpha$
against Weber number
$\textit{We}$
for all experiments (circles) and DNS (stars). The KM model is shown as solid lines for
$Bo=0.02$
, and the black dashed line denotes the inertio-capillary limit (
$Oh=0$
;
$Bo=0$
). The red dot-dashed line represents (4.3). For low
$Oh$
and
$\textit{We}$
, the data follow a power-law trend captured by the KM model and energy argument, but flattens when
$\textit{We} \gtrapprox 0.5$
due to self-occlusion of the drop surface, a result of our definition for
$H$
, and because
$1-\alpha$
must eventually saturate to unity. The KM model and DNS results capture this feature for all
$Oh$
, closely following experimental results. As
$Oh$
increases,
$1-\alpha$
decreases and follows a weaker power-law trend. Additionally, self-occlusion is delayed to larger
$\textit{We}$
or disappears altogether for our largest
$Oh$
. The inset in figure 8(
$b$
) shows
$1-\alpha$
against
$\beta -1$
for all experiments and DNS. The data follow the linear trend predicted by (4.4) for
$\beta -1 \lessapprox 0.3$
. At higher deformations, the data continue to follow a single trend, but necessarily deviate from the energy argument since the deformations are no longer small and
$1-\alpha$
must saturate to unity.
The maximum spreading radius
$\beta _c$
is a key parameter in determining the heat transfer and total interaction between a drop and substrate in many applications. When
$\textit{We} \gt 1$
, an impacting drop flattens into a pancake-like shape, such that
$\beta _{c} \approx \beta - 1$
, making
$\beta - 1$
a convenient and practical proxy for
$\beta _c$
. This is invalid at low
$\textit{We}$
, and necessitates a revised relationship between
$\beta _c$
and
$\beta$
. From our energy argument, (4.5) predicts a linear relationship between the maximum drop deformation and maximum spreading area
$\beta _c^{2}$
. Figure 8(
$c$
) plots
$\beta _c^{2}$
against
$\beta - 1$
for all experiments (circles) and DNS (stars), with (4.5) overlaid as a red dot-dashed curve. Strikingly, this prediction is accurate across all dimensionless parameters explored, and provides a simple relationship connecting drop spreading and deformation in the low-
$\textit{We}$
regime.
5. Discussion
Previous attempts to quantify rebound behaviour for
$\textit{We}\lt 1$
used superhydrophobic substrates, where contact line friction altered and ultimately suppressed rebound (Wang et al. Reference Wang, Zhao, Sun, Mehrizi, Lin, Guo and Chen2022; Thenarianto et al. Reference Thenarianto, Koh, Lin, Jokinen and Daniel2023). We overcame this limitation by using a thin viscous film as the substrate, which allowed rebound at as low as
$\textit{We}=0.0008$
. Additionally, we provide the first systematic study of the role of
$Oh$
for
$\textit{We}\lt 1$
rebounds. Our main finding is a transition to unique trends in the rebound metrics – the coefficient of restitution
$\varepsilon$
and contact time
$t_{c}$
– at low
$\textit{We}$
. These trends were shown to coincide with a change in the spreading and retraction dynamics, which become symmetric in the inertio-capillary limit, leading to nearly elastic rebound. We complemented these findings with reduced-order models formulated to describe low-
$\textit{We}$
rebound, and show that they provide accurate predictions for the rebound metrics, drop deformation and drop spreading.
The coefficient of restitution
$\varepsilon$
and contact time
$t_{c}$
have different
$\textit{We}$
dependencies at low and high
$\textit{We}$
. When
$\textit{We} \gtrapprox 0.5$
, the coefficient of restitution
$\varepsilon$
decreased as
$\textit{We}$
increased, as observed in prior studies (Biance et al. Reference Biance, Chevy, Clanet, Lagubeau and Quéré2006; Jha et al. Reference Jha, Chantelot, Clanet and Quéré2020). When
$\textit{We} \lessapprox 0.5$
,
$\varepsilon$
reached a maximum that depends upon
$Oh$
and
$Bo$
, and becomes
$\textit{We}$
-independent in the inertio-capillary limit. A gravity-induced roll-off in
$\varepsilon$
leading to
$\varepsilon \rightarrow 0$
is the natural route to rebound suppression when
$Bo$
is finite, with
$Oh$
increasing the critical Weber number
$\textit{We}_c$
below which rebound stops. These results differ from the peak and subsequent roll-off in
$\varepsilon$
as
$\textit{We}$
decreases reported for
$\textit{We}\lt 1$
rebounds on superhydrophobic substrates (Richard & Quéré Reference Richard and Quéré2000; Wang et al. Reference Wang, Zhao, Sun, Mehrizi, Lin, Guo and Chen2022; Thenarianto et al. Reference Thenarianto, Koh, Lin, Jokinen and Daniel2023). Although these prior experiments were conducted in low-
$Oh$
conditions, and in some cases used low-
$Bo$
drops, the increased effect of contact line friction at low
$\textit{We}$
was shown to be the primary mechanism for rebound suppression. Consequently, our results show that submillimetre drops, like those from coughing (Dbouk & Drikakis Reference Dbouk and Drikakis2020) or agricultural sprays (Kooij et al. Reference Kooij, Sijs, Denn, Villermaux and Bonn2018; Makhnenko et al. Reference Makhnenko, Alonzi, Fredericks, Colby and Dutcher2021; Chen & Ashgriz Reference Chen and Ashgriz2022), may continue bouncing longer than previously thought when the substrate is non-wetting – an important insight into the spread of their contents. The contact time between the drop and substrate is
$\textit{We}$
-dependent at low
$\textit{We}$
(Foote Reference Foote1975; Okumura et al. Reference Okumura, Chevy, Richard, Quéré and Clanet2003), and our results show that it is also
$Oh$
-dependent in this regime. Increasing
$Oh$
causes
$t_{c}$
to increase at a higher rate as
$\textit{We}$
decreases such that viscosity is most influential on
$t_{c}$
in the low-
$\textit{We}$
regime. Since
$t_{c}$
is a critical parameter governing drop–substrate interaction in various applications (Shiri & Bird Reference Shiri and Bird2017; Díaz et al. Reference Díaz, Garcia-Gonzalez, Bista, Weber, Butt, Stetten and Kappl2022), this finding has immediate practical implications, and motivates the study of low-
$\textit{We}$
rebound using drops with shear-dependent viscosity in future work.
In § 4.2, we compared the drop shape during spreading and retraction for three values of
$\textit{We}$
spanning the transition from low
$\textit{We}$
to high
$\textit{We}$
, and revealed a nearly symmetric spreading and retraction process for low-
$Oh$
, low-
$\textit{We}$
rebound. As
$\textit{We}$
increased, the spread time decreased relative to the retraction time, increasing asymmetry between spreading and retraction and aligning with the changes in rebound behaviour. This change corresponds to a loss of energy during rebound when
$\textit{We}\gt 1$
and was investigated through a careful energy budget analysis by Gilet & Bush (Reference Gilet and Bush2012). Their results showed that when
$\textit{We}\lt 1$
, the initial translational energy was efficiently converted to surface energy during the spreading phase. When
$\textit{We}\gt 1$
, this conversion was inefficient, with an appreciable loss of energy during spreading. This suggests that the time symmetry is broken in the spreading phase at high
$\textit{We}$
. Additionally, Bartolo et al. (Reference Bartolo, Josserand and Bonn2005) identified a limit on the retraction rate of a drop that is substrate-independent for low
$Oh$
. Although this limit does not explain the onset of asymmetry, which requires considering the energy lost during spreading, it helps rationalise how the spreading-to-retraction time ratio evolves as
$\textit{We}$
increases when
$\textit{We}\gt 1$
. In our experiments, this ratio follows the scaling
$t_{s}/t_{r} \sim \textit{We}^{-0.26 \pm 0.01}$
in the low-
$Oh$
limit. We can rationalise this observation by considering the inertio-capillary limit, where the time needed for a drop to spread
$t_{s} \sim R_{c} / V$
, where
$R_{c} \sim R \textit{We}^{-1/4}$
from our energy argument (4.6), and the maximum retraction rate
$t_{r}^{-1} \sim t_{\sigma }^{-1}$
(Bartolo et al. Reference Bartolo, Josserand and Bonn2005). This analysis gives
$t_{s} / t_{r} \sim \textit{We}^{-1/4}$
, in good agreement with our results. We therefore interpret the growing asymmetry in rebound dynamics at
$\textit{We} \gt 1$
as the result of two effects: a loss of energy during spreading that breaks time symmetry early in the rebound process (Gilet & Bush Reference Gilet and Bush2012), and a limited retraction rate that further skews the rebound time scale (Bartolo et al. Reference Bartolo, Josserand and Bonn2005). Recent studies have also noted the near symmetry between spreading and retraction times at low
$\textit{We}$
(Lin et al. Reference Lin, Shi, Zhou, Chen and Li2024). Our results confirm and extend this observation, by demonstrating that time symmetry holds in the inertio-capillary regime and is progressively lost as either
$\textit{We}$
or
$Oh$
increases.
Lastly, we examined the drop deformation and spreading in § 4.3, and showed that the KM model captures drop deformation well across a wide range of
$Oh$
, and agrees with a simple energy model in the low-
$Oh$
limit, indicating that deformation is primarily in the
$l=2$
mode. The literature on the maximum horizontal deformation is vast and two competing scaling arguments have gained support: one based on momentum conservation
$\beta \sim \textit{We}^{1/4}$
(Clanet et al. Reference Clanet, Béguin, Richard and Quéré2004; Bartolo et al. Reference Bartolo, Josserand and Bonn2005; Biance et al. Reference Biance, Chevy, Clanet, Lagubeau and Quéré2006; Tsai et al. Reference Tsai, Hendrix, Dijkstra, Shui and Lohse2011; García-Geijo et al. Reference García-Geijo, Riboux and Gordillo2020) and the other arising from energy conservation
$\beta \sim \textit{We}^{1/2}$
(Chandra & Avedisian Reference Chandra and Avedisian1991; Bennett & Poulikakos Reference Bennett and Poulikakos1993; Okumura et al. Reference Okumura, Chevy, Richard, Quéré and Clanet2003; Eggers et al. Reference Eggers, Fontelos, Josserand and Zaleski2010; Villermaux & Bossa Reference Villermaux and Bossa2011; Laan et al. Reference Laan, de Bruin, Bartolo, Josserand and Bonn2014; Lee et al. Reference Lee, Laan, de Bruin, Skantzaris, Shahidzadeh, Derome, Carmeliet and Bonn2016). Despite extensive efforts, the correct scaling remains debated, in part due to most studies using the maximum equatorial radius
$\beta = R_{m} /R$
to describe horizontal deformation. However, as discussed in Josserand & Thoroddsen (Reference Josserand and Thoroddsen2016),
$\beta$
ignores the initial drop radius and is inadequate to discern scaling arguments at low
$\textit{We}$
. In this work we used the maximum equatorial deformation
$\beta -1$
, and energy arguments predict
$\beta -1 = (5/48 )^{1/2} \, \textit{We}^{1/2}$
, as in Richard & Quéré (Reference Richard and Quéré2000). This prediction is in excellent agreement with our low-
$Oh$
data and requires no fitting parameters. Furthermore, we combined our low-
$Oh$
data with those from previous studies and validated (4.2) from
$\textit{We} = 10^{-4}$
to
$\textit{We} = 10^{2}$
, above which the data begin following a trend closer to
$ \beta - 1 \sim \textit{We}^{1/4}$
. As
$Oh$
increases, (4.2) is less accurate; however, the KM model continues to provide accurate predictions. The KM model also captures the maximum vertical deformation
$1-\alpha$
across all
$Oh$
, and our energy argument equation (4.3) provides a simple and accurate prediction when
$Oh$
and
$Bo$
effects are negligible. In addition to relating vertical and horizontal deformation, our energy argument allows us to relate maximum deformation and spreading. For example, (4.5) predicts a linear relationship between the horizontal deformation and maximum spread area
$\beta _c^2$
, which holds for all
$Oh$
and
$Bo$
tested. This relationship provides a useful tool for determining the contact area of a bouncing drop from its deformation, which is critical for applications such as heat transfer (Shiri & Bird Reference Shiri and Bird2017) or predicting the net charge acquired by a rebounding water drop (Díaz et al. Reference Díaz, Garcia-Gonzalez, Bista, Weber, Butt, Stetten and Kappl2022).
6. Conclusions
We studied a drop rebounding from a non-wetting, rigid substrate at low
$\textit{We}$
using experiments, DNS and reduced-order modelling. We revealed the rebound behaviour of drops at low
$\textit{We}$
for a range of
$Oh$
, showed how increasing
$Oh$
or
$\textit{We}$
breaks the symmetry between the spreading and retraction phases and elucidated the critical role of
$Oh$
and
$Bo$
in suppressing rebound. The maximum drop deformation and spreading were connected through energy arguments, and closely followed the predictions of the KM model irrespective of
$Oh$
. These findings reveal new physical insight into drop rebound and provide practical tools for predicting drop deformation on non-wetting substrates.
The interplay between our theoretical approaches, the KM model and the dedicated DNS set-up enabled us to cross-validate these tools and establish the range of validity of the reduced-order model with a significant level of confidence, with the experimental dataset acting as foundation where comparisons were possible. Thus, we were suitably equipped to use the more efficient methodology in the relevant section of the parameter space and relevant distinguished limits thereof, while also being able to rely on DNS for larger values of
$\textit{We}$
where the model assumptions no longer hold. In general, DNS is preferable for
$\textit{We} \gt 1$
, while the reduced-order model offers high-accuracy predictions with substantially lower computational cost for
$\textit{We} \lt 1$
, but continues to yield reasonable estimates up to
$\textit{We} \approx 10$
. Furthermore, our work also presents the first application of the KM method to the case of a deformable impactor rebounding off a rigid substrate. The quality of the resulting predictions highlights its versatility and naturally invites further applications for it; e.g. its use in the case of a collision of two deformable objects, and in combination with moving meshes and variational methods.
When contact line formation is avoided, the coefficient of restitution
$\varepsilon$
is
$\textit{We}$
-independent at low
$\textit{We}$
in the inertio-capillary limit. Finite-
$Bo$
effects suppress bouncing while increasing
$Oh$
decreases
$\varepsilon$
and increases the critical Weber number
$\textit{We}_{c}$
where rebound ceases. In contrast, the contact time
$t_{c}$
increases as
$\textit{We}$
decreases at a rate that increases with
$Oh$
. The shift in how
$\varepsilon$
and
$t_{c}$
change with
$\textit{We}$
in the low-
$\textit{We}$
and high-
$\textit{We}$
regimes coincides with altered spreading and retraction dynamics. In particular, the nearly elastic rebound observed in the inertio-capillary limit has a nearly symmetric spreading and retraction phase, but increasing
$\textit{We}$
above unity and increasing
$Oh$
decreases
$\varepsilon$
and increases asymmetry in the spreading and retraction process. Using energy arguments, we derive simple equations that accurately predict the maximum vertical and horizontal deformation, and the maximum spreading for low-
$Oh$
impacts with no fitting parameters. The prediction for the maximum horizontal deformation (4.2) reveals that the horizontal deformation of an impacting drop scales as
$\textit{We}^{1/2}$
, which when compared with experimental data holds up to
$\textit{We}\approx 10^{2}$
. Additionally, this simple model accurately relates
$\beta -1$
to the maximum vertical deformation
$1-\alpha$
via (4.4) and maximum contact area
$\beta _{c}^{2}$
via (4.5). For all
$Oh$
, the KM model provides excellent predictions for the maximum spreading and deformation for low-
$\textit{We}$
rebound where DNS struggles due to numerical stiffness and parasitic capillary currents.
Although we focus here on the normal impact of a drop on a rigid surface, the ubiquity of drop impacts makes our findings relevant to many situations. For example, low-
$\textit{We}$
rebound governs the final stage of drop dispersal in many applications involving low-
$Bo$
drops, such as agricultural sprays (Kooij et al. Reference Kooij, Sijs, Denn, Villermaux and Bonn2018; Makhnenko et al. Reference Makhnenko, Alonzi, Fredericks, Colby and Dutcher2021; Chen & Ashgriz Reference Chen and Ashgriz2022) and cough ejecta (Dbouk & Drikakis Reference Dbouk and Drikakis2020; Bourouiba Reference Bourouiba2021), where understanding the spread of drop contents is essential. In both examples, the drops typically impact non-wetting surfaces, as plant leaves are often naturally superhydrophobic (Barthlott & Neinhuis Reference Barthlott and Neinhuis1997), while medical surfaces are frequently rendered superhydrophobic to reduce microbial adhesion and biofluid contamination (Falde et al. Reference Falde, Yohe, Colson and Grinstaff2016). Future work is needed to extend our results to more complex practical and natural settings, which often involve non-Newtonian drops and complaint substrates that have a pronounced effect on drop impact at high
$\textit{We}$
(Bartolo et al. Reference Bartolo, Boudaoud, Narcy and Bonn2007; Basso & Bostwick Reference Basso and Bostwick2020) but remain unexplored for drop rebound at low Weber number. Furthermore, our focus herein has been strictly on the liquid dynamics whereas the coupled air layer dynamics is also a source of rich physics, and has become of increasing interest to modellers (Phillips & Milewski Reference Phillips and Milewski2024; Sprittles Reference Sprittles2024).
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2025.10589.
Funding
The authors gratefully acknowledge the financial support of the National Science Foundation (NSF CBET-2123371) and the Engineering and Physical Sciences Research Council (EPSRC EP/W016036/1). C.T.G. gratefully acknowledges the Hope Street Fellowship for financial support.
Declaration of interests
The authors report no conflict of interest.
Data availability
The KM code and data used in each figure are available at https://github.com/harrislab-brown/LowWeberDropRebound. The DNS code infrastructure, including pre- and post-processing functionality and a short descriptive tutorial, can be found at https://github.com/rcsc-group/BouncingDropletsSolidSurface.
Appendix A. Defining contact
Determining the substrate location is crucial for evaluating drop–substrate interactions. However, this is typically constrained by the spatial resolution of the measurement technique. In this study, the presence of an intervening air layer further complicates identifying the substrate, as it prevents direct contact with the drop. In our experiments the spatial resolution related to the pixel width in our images was approximately
$0.02R$
, which can lead to discrepancies when compared with KM model and simulation results if not considered. Figure 9(
$a$
) shows a drop impacting a substrate with a magnified view of the contact boundary. Increasing the substrate location by the spatial resolution of
$0.02R$
(grey dashed line) shifts the true contact point (open circle) to the perceived contact point (filled circle), expanding the contact radius by an amount
$b$
and necessitating careful criteria for faithfully comparing experimental and modelling results. The non-wetting nature of the surface makes this process particularly sensitive due to the tangency with the substrate.
(
$a$
) Magnified view near a drop’s contact point during impact. The open circle indicates the contact point with the substrate, while the filled circle marks the contact point along a hypothetical substrate at a height corresponding to the experimental spatial resolution
$z=0.02R$
. This threshold increases the contact time and expands the contact radius
$r_c$
by
$b$
. (
$b$
) Dimensionless contact time
$\bar {t}_{c} = t_{c} / t_{\sigma }$
and (
$c$
) dimensionless maximum contact radius
$\beta _{c} = R_{c} / R$
against Weber number for the KM model evaluated at the hypothetical substrate
$0.02 R$
(grey) and true substrate
$0.00 R$
(black) for
$Oh = 0.03$
. Red curves represent predictions from the model of Chevy et al. (Reference Chevy, Chepelianskii, Quéré and Raphaël2012). Solid lines show predictions for
$Bo=0.019$
and dashed lines show results obtained for
$Bo=0$
.
To provide the most accurate predictions between the KM model, DNS and experiments, the former two measure all relevant variables using a hypothetical substrate located
$0.02R$
above the true substrate. This criterion also allows DNS measurements to be made above the intervening air layer, which would otherwise need to be accounted for during post-processing. The parameters most affected by this criterion are the contact time,
$t_{c}$
, and contact radius,
$r_c$
, which both increase as a result. To show the magnitude of this effect, we plot the dimensionless contact time
$t_c / t_{\sigma }$
against
$\textit{We}$
for the KM model with (grey lines) and without (black lines) the
$0.02R$
criterion in figure 9(
$b$
), where solid lines represent
$Bo = 0.019$
and dashed lines
$Bo = 0$
. The influence of the threshold increases as
$\textit{We}$
decreases. Figure 9(
$c$
) shows the KM model for the dimensionless maximum contact radius
$\beta _{c} = R_{c} / R$
versus
$\textit{We}$
for the same threshold values. Again, the threshold effect increases as
$\textit{We}$
decreases; however, unlike contact time, noticeable differences in the predictions persist even for
$\textit{We}\gt 1$
.
By employing a discretised angular mesh for the polar coordinate, only a finite number of contact radii are numerically allowed. As a consequence, the
$0.00R$
criterion predicts piecewise constant functions for the contact radius
$r_c$
, while the
$0.02R$
criterion smooths this discretised behaviour and results in smoother curves, as seen in figure 9(
$c$
).
We also compare these predictions with results obtained using the quasi-static model of Chevy et al. (Reference Chevy, Chepelianskii, Quéré and Raphaël2012). Their work predicts contact time
$t_{c}$
, but additional features about the shape, such as the contact radius
$r_c$
, are also extractable by reproducing their calculations. In figure 9(
$b,c$
), the red curves show predictions from the model of Chevy et al. (Reference Chevy, Chepelianskii, Quéré and Raphaël2012) with the solid line for
$Bo=0.019$
and dashed line for
$Bo=0$
. While their modelling approach differs, their low-
$\textit{We}$
predictions align very well with ours, but the KM model is able to make predictions above
$\textit{We} \approx 10^{-1}$
while also enforcing the natural geometric constraints in the contact region at all points.
Appendix B. Numerical approximations
B.1. The spectral method
We define
$A_l$
,
$U_l$
,
$B_l$
,
$\textit {h}$
and
$\textit {v}$
as approximations of
$\mathscr{A}_l$
,
$\mathscr{U}_l$
,
$\mathscr{B}_l$
,
$h$
and
$v$
, respectively. These represent the solutions to the system equations obtained by truncating the expansions of the surface shape, velocity and pressure distribution in spherical harmonics up to the
$L\textrm{th}$
mode. That is,
\begin{align}&\qquad\qquad\qquad\qquad\qquad \dot {A}_l = U_l, \quad \text{ for } l = 2,\ldots ,L; \end{align}
\begin{align}& \dot {U}_l = -l \left (l+2\right )\left (l-1\right )A_l -2(2l+1)(l-1)\, Oh \, {U}_l -l \, B_l, \qquad \text{ for } l = 2,\ldots , L; \end{align}
\begin{align}&\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad \dot {\textit {h}} = \textit {v}; \end{align}
\begin{align}&\qquad\qquad\qquad\qquad\qquad\qquad\quad \dot {\textit {v}} = - Bo + B_1; \\[9pt] \nonumber \end{align}
which are subject to
\begin{align}&\qquad A_l(t=0) = 0, \quad \text{ for }\ l = 2,\ldots , L; \end{align}
\begin{align}&\qquad U_l(t=0) = 0, \quad \text{ for }\ l = 2,\ldots , L; \end{align}
\begin{align}&\qquad\qquad\qquad \textit {h}(t=0) = 1; \end{align}
\begin{align}&\qquad\qquad\quad \textit {v}(t=0) = -\sqrt {We}; \end{align}
\begin{align}& \sum \limits _{l=2}^{L} {A_lP_l\left (\cos (\theta )\right )} = \frac {\textit {h}}{\cos (\theta )}-1, \quad \text{ for } \theta \leqslant \theta _c; \end{align}
\begin{align}& \sum \limits _{l=2}^L {A_lP_l(\cos (\theta ))} \gt \frac {\textit {h}}{\cos (\theta )}-1, \quad \, \text{ for } \theta \gt \theta _c; \end{align}
\begin{align}&\qquad \sum \limits _{l=0}^{L} {B_lP_l(\cos (\theta ))} = 0, \quad \text{ for }\ \theta \gt \theta _c; \end{align}
\begin{align}& \frac {\partial }{\partial \theta }\left [\left (1+\sum \limits _{l=2}^L {A_lP_l(\cos (\theta ))}\right )\cos (\theta )\right ]_{\theta =\theta _c} = 0. \end{align}
For the present work, calculations were performed using
$L=90$
spherical modes. Convergence was verified by testing resolutions of 120, confirming that all quantities of interest varied by less than 10 % across the full dimensionless parameter range – representing the worst-case deviation – with an average relative error below 5 %.
B.2. Finite differences in time
We approximate the time derivatives by defining
$A_l^k \approx A_l(t=t^k) \text{, } U_l^k \approx U_l(t=t^k)\text{, } B_l^k \approx B_l(t=t^k) \text{, } h^k \approx \textit {h}(t=t^k) \text{, } v^k \approx \textit {v}(t=t^k)$
and applying the implicit Euler method to (B1):
\begin{equation} A^{k+1}_l -\delta _t^k U^{k+1}_l = A^k_l \qquad \text{ for } l = 2, \ldots L, \end{equation}
where
$\delta ^k_t := t^{k+1}-t^k$
. Similarly, from (B2), we have
\begin{equation} \begin{aligned} \delta ^k_t \, l \, (l+2)(l-1) \, A^{k+1}_l + \big [1 + \delta ^k_t \, 2(2l+1)(l-1) \, Oh \big ] \,{U}_l^{k+1} + \delta ^k_t \, l \, B_l^{k+1} =U^k_l,& \\ \text{ for } l = 2, \ldots L,& \end{aligned} \end{equation}
from (B3), we have
\begin{equation} - \delta _t^k \, v^{k+1} + h^{k+1} = h^{k} \end{equation}
and from (B4), we have
\begin{equation} -\delta _t^k \, B_1^{k+1} + v^{k+1} = v^k - \delta _t^k \, Bo; \end{equation}
which are subject to
\begin{align}& \sum \limits _{l=2}^{L} {A_l^{k+1}P_l\left (\cos (\theta )\right )} - \frac {h^{k+1}}{\cos (\theta )} = -1, \quad \text{ for }\ \theta \leqslant \theta _c; \end{align}
\begin{align}& \cos (\theta )\left (1+\sum \limits _{l=2}^{L} {A_l^{k+1}P_l\left (\cos (\theta )\right )}\right ) \gt h^{k+1}, \quad \text{ for }\ \theta \gt \theta _c; \end{align}
\begin{align}&\qquad\qquad \sum \limits _{l=0}^{L} {B_l^{k+1}P_l(\cos (\theta ))} = 0, \quad \text{ for }\ \theta \gt \theta _c; \end{align}
\begin{align}&\qquad \frac {\partial }{\partial \theta }\left [\left (1+\sum \limits _{l=2}^L {A_l^{k+1}P_l(\cos (\theta ))}\right )\cos (\theta )\right ]_{\theta =\theta _c} = 0. \end{align}
B.2.1. Angular discretisation
The spectral mesh is determined by the maximum order of harmonic modes
$L$
in our truncated expansion, using Legendre polynomials of degree
$0$
to
$L$
. The angular mesh is naturally defined as the union of
$\theta = 0$
(the south pole) and the zeros of
$P_L(\cos (\theta ))$
, which is canonical for the Legendre polynomial basis (Boyd Reference Boyd2000). A schematic of this angular mesh is shown in figure 3(
$c$
).
The candidate integer
$q$
defines the number mesh points from the south pole to the boundary of the contact area, where the pressure is allowed to be non-zero. The boundary angle of the contact area is taken as the midpoint between
$\theta _q$
and
$\theta _{q+1}$
. This results in the discrete version of the constraints imposed by (B17)–(B20):
\begin{align}&\qquad \sum \limits _{l=2}^{L} {A_l^{k+1, q}P_l\left (\cos (\theta _i)\right )} - \frac {h^{k+1, q}}{\cos (\theta _i)} = -1, \quad \text{ for }\ i \leqslant q; \end{align}
\begin{align}& \cos (\theta _i)\left (1+\sum \limits _{l=2}^{L} {A_l^{k+1, q}P_l\left (\cos (\theta _i)\right )}\right ) \gt h^{k+1, q}, \quad \text{ for }\ i \gt q; \end{align}
\begin{align}&\qquad\qquad\quad \sum \limits _{l=0}^{L} {B_l^{k+1, q}P_l(\cos (\theta _i))} = 0, \quad \text{ for }\ i \gt q; \end{align}
and
\begin{equation} \frac {\partial }{\partial \theta }\left [\left (1+\sum \limits _{l=2}^L {A_l^{k+1,q}P_l(\cos (\theta ))}\right )\cos (\theta )\right ]_{\left (\theta _q+\theta _{q+1}\right )/2} = 0. \end{equation}
Here,
$A_l^{k+1,q}$
and
$B_l^{k+1,q}$
are the amplitudes of the surface deformation and pressure distribution, respectively, obtained by solving the discrete forms of (B21), (B23), (B25)–(B28), assuming exactly
$q$
mesh points in the contact area:
\begin{align}&\qquad\qquad\qquad\qquad A^{k+1, q}_l -\delta _t^k \, U^{k+1, q}_l = A^k_l, \quad \text{ for } l = 2, \ldots L, \end{align}
\begin{align}& \delta ^k_t\, l \, (l+2)(l-1) \, A^{k+1, q}_l + \big [1 + \delta ^k_t \, 2(2l+1)(l-1) \, Oh \, \big] \, {U}_l^{k+1, q} + \delta ^k_t \, l \, B_l^{k+1, q} =U^k_l,\nonumber \\&\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad \text{ for } l = 2, \ldots L, \end{align}
\begin{align} - \delta _t^k \, v^{k+1, q} + h^{k+1, q} = h^{k} \end{align}
and
\begin{equation} -\delta _t^k \, B_1^{k+1,q} + v^{k+1, q} = v^k - \delta _t^k \, Bo, \end{equation}
where
$q$
in the superscript denotes solutions calculated assuming
$q$
contact points.
B.3. A discrete optimisation problem
For a given number of contact mesh points
$q$
, we can form a closed, linear system of equations with (B23), (B21) and (B25)–(B28). We use (B22) to define
\begin{equation} e^1_q = \begin{cases} 0,\ &\cos (\theta _i)\left (1+\sum \limits _{l=2}^{L} {A_l^{k+1, q}P_l\left (\cos (\theta _i)\right )}\right ) \gt h^{k+1, q}, \ \text{ for all }\ i \gt q; \\ \infty ,\ &\text{otherwise}. \end{cases} \end{equation}
This discrete version of condition (B22) ensures the drop does not touch or penetrate the substrate outside the contact area.
For each
$q$
, the values of
$A_l^{k+1,q}$
are substituted into the second-order approximation of (B24), and the absolute value of the resulting residual is computed as
\begin{align} e^2_q &= \left | \left [1+ \sum \limits _{l=2}^{L} {A_l^{k+1, q}P_l\left (\cos (\theta _q)\right )} \right ]\cos \left (\theta _q\right ) \right . \nonumber\\&\quad \left .-\, \left [1+ \sum \limits _{l=2}^{L} {A_l^{k+1, q}P_l\left (\cos (\theta _{q+1})\right )} \right ]\cos \left (\theta _{q+1}\right ) \right |. \end{align}
Finally, we define
\begin{equation} e_q = \max \big \{e^1_q,e^2_q\big \}. \end{equation}
The value of
$q$
that minimises
$e_q$
defines the contact mesh points at time
$k+1$
; that is,
\begin{equation} m^{k+1} := \text{argmin}_q e_q, \end{equation}
and, consequently,
\begin{align} & A_l^{k+1} := A_l^{k+1,m^{k+1}}, \ h^{k+1} := h^{k+1,m^{k+1}}, \ v^{k+1} := v^{k+1,m^{k+1}}, \end{align}
\begin{align} & B_l^{k+1} := B_l^{k+1,m^{k+1}}\quad \text{and}\quad U_l^{k+1} := U_l^{k+1,m^{k+1}}. \end{align}
B.4. The system matrices
The discretised problem involves solving a linear system of equations at each sampling time for each candidate contact area radius, parameterised by
$q$
. These systems can be expressed in matrix form as
\begin{equation} M^{k+1,q} \, x^{k+1,q} = b^{k+1,q}, \end{equation}
where
\begin{align} & x^{k+1,q} = \Big [\! A_2^{k+1,q}\!, \ldots , A_L^{k+1,q}\!, U_2^{k+1,q}\!, \ldots , U_L^{k+1,q}\!, B_0^{k+1,q}\!, \ldots , B_L^{k+1,q}\!, h^{k+1,q}\!, v^{k+1,q} \!\Big ]^T \!, \end{align}
\begin{align} &\quad b^{k+1,q} = \left [A_2^k,\ldots ,A_L^k,U_2^k,\ldots ,U_L^k,\underbrace {-1,\ldots ,-1}_{q \text{ times}},\underbrace {0,\ldots ,0}_{L+1-q \text{ times}},h^k,v^k-\delta _t^{k} \, Bo \right ]^T \end{align}
and
\begin{equation} \textbf{M}^{k+1,q} = \left [ \begin{array}{ccccc} I_{L-1} & -\delta _t^kI_{L-1} & 0& 0 & 0 \\ \delta _t^k C & D & \delta _t^k E & 0 & 0 \\ F^q & 0 & 0 & G^q & 0 \\ 0 & 0 & H^q & 0 & 0 \\ 0 & 0 & 0 & 1 & -\delta _t^k \\ 0 & 0 & \delta _t^k K & 0 & 1 \end{array} \right ]. \end{equation}
Here,
$I_{m}$
is the identity matrix of size
$m$
, and from (B26) we have
\begin{align} &\quad\,\,\,\, C_{i,j} = j \, (j+1)(j+3) \, \delta _{i,j}, \end{align}
\begin{align} & D_{i,j} = \big (1 + \delta ^k_t \, 2 \, Oh \,j \, (2j+3)\big) \, \delta _{i,j} \end{align}
for all
$i, j = 1, 2, \ldots , L-1$
, where
$\delta _{i, j}$
is the Kronecker delta. Also,
\begin{equation} E = \left [ \begin{array}{cc|ccccc} 0 & 0 & 2 & 0 &\ldots &\ldots &0\\ 0 & 0 & 0 & 3 &0 &\ldots &0\\ \vdots & \vdots & \vdots & 0 &\ddots & &\vdots \\ \vdots & \vdots & \vdots & \vdots & &(L-1) &0\\ 0 & 0 & 0 & 0 &0 &\ldots &L\\ \end{array} \right ], \end{equation}
where the first two zero columns represent modes
$0$
and
$1$
, absent in (B26).
Similarly, from (B21),
\begin{equation} F_{i,j}^q = P_{j+1}(\cos (\theta _{i-1})), \quad \text{ for } i= 1,2,\ldots ,\quad \text{and } j = 1,2,\ldots ,L-1; \end{equation}
and
\begin{equation} G_i^q = \left [ \frac {-1}{\cos (\theta _i)} \right ], \quad \text{ for } i= 1,2,\ldots ,q. \end{equation}
From (B23) we have
\begin{equation} H_{i,j}^q = P_{j-1}(\cos (\theta _{i + q-1})), \qquad \text{for } i= 1, 2,\ldots ,L- q +1\quad \text{and } j = 1,2,\ldots ,L+1; \end{equation}
and from (B28) we obtain
\begin{equation} K = \left [\underbrace { \begin{array}{ccccc} 0 & -1 & 0 & \ldots & 0 \end{array}}_{L+1 \text{ elements}} \right ]. \end{equation}
In summary, the discrete problem is reduced to finding the
$q$
that satisfies (B32). Instead of solving system (B35) for every
$q$
, we solve it for
$q = m^k-2, m^k-1, m^k, m^k+1, m^k+2$
. This approach significantly reduces computational costs without decreasing accuracy, enabling typical runs to complete within minutes.
Appendix C. Computational implementation
C.1. The adaptive time step
We set a maximum dimensionless step size of
$2\pi / (S_L ( L(L+2)(L+1) )^{1/2} )$
, where
$S_L = 16$
, chosen as a fraction
$1/S_L$
of the oscillation period of the fastest harmonic mode. The adaptive time step refines the initial time mesh by adding discrete times as needed. The step size is dynamically adjusted to ensure the contact area evolution is resolved with the highest accuracy permitted by the angular mesh. Specifically, if the model predicts that the number of contact mesh points
$m^{k}$
changes by more than one in a single step, the step is rejected, halved and candidate solutions are recalculated. This approach relies on the assumption that the contact area evolves continuously during impact.
The code iteratively refines the time step as needed to ensure that the motion of the contact area’s boundary is resolved to the accuracy of the angular mesh. Figure 10 illustrates two consecutive time-step refinements between
$t^{k+3}$
and
$t^{k+4}$
, following the rejection of the initial calculation
$t^{k+4}$
. This results in the insertion of
$t^{k+3^*}$
, which, if also rejected, leads to the insertion of
$t^{k+3^{**}}$
. While not shown in the figure, the discrete time sequence is renumbered to include these refinements.
Time line of regularly spaced intervals (black) of spacing
$\delta _{t}^{k}$
. Blue and red lines represent the time after one and two refinement steps between
$t^{k+3}$
and
$t^{k+4}$
, respectively, where
$\delta _{t}^{k}$
is halved at each.
Normalised error
$\bar {E} = 2 \, \textrm {tan}^{-1}(E) / \pi$
versus candidate contact point
$q$
for a drop with 10 contact points. An infinite error
$E = \infty$
(dashed line) is assigned for
$q \lt 10$
using (B29). For
$q \geqslant 10$
, the error is calculated using (B30) and increases with
$q$
, showing a minimum at
$q=10$
. The inset shows a schematic of the accepted solution. (
$b$
) Relative change in the coefficient of restitution
$\tilde {\varepsilon } = \textrm {log} ( |\varepsilon _L-\varepsilon _{L/2}|/\varepsilon _{L/2} )$
against number of modes log(
$L$
). The dashed line is a power-law fit to the data with slope
$m = -1.516$
.
C.2. Finding local minima of the error
Assuming the contact area boundary moves continuously, we can reduce the computational cost of solving the discrete optimisation problem in (B32), which would otherwise require evaluating all possible
$q$
from
$0$
to
$L$
. The error in the angle between the drop and substrate (the quantity to minimise) exhibits clear monotonic behaviour around the minimum. Figure 11(
$a$
) shows the normalised error
$\bar {E} = 2 \, \textrm {tan}^{-1}(E) / \pi$
plotted against
$q$
for a case where the solution indicates that
$q=10$
contact points is the optimal solution, i.e.
$m^{k+1}=10$
. Our normalisation of the error allows us to display both finite and infinite errors in the same plot, with
$\bar {E}=1$
for
$E = \infty$
.
Flowchart of computational implementation. The return variables are the coefficient of restitution
$\varepsilon$
, minimum height of the centre of mass
$h_{{min}}$
, maximum equatorial radius
$R_{{m}}$
, maximum spreading radius
$R_c$
and contact time
$t_{c}$
.
These results indicate that to determine
$m^{k+1}$
, the number of contact points at time
$t^{k+1}$
, it is sufficient to test a subset of five values:
$q = m^{k}-2$
to
$q = m^{k}+2$
. We accept only a minimum of
$e_q$
within one mesh point of
$m^k$
, and these five points provide enough information to identify any local minimum among the three acceptable points. If the minimum occurs at
$q = m^{k}-2$
or
$q = m^{k}+2$
, the solution is rejected, and the time step is halved; otherwise, the solution is accepted.
In practice, we first check if
$q = m^k-1$
is a local minimum (unless
$m^k=0$
). To do this, we compute the errors for
$q = m^{k}-2$
to
$q = m^{k}$
. If
$q = m^{k}-1$
is the minimum, we accept it as the correct solution. If
$q = m^{k}-2$
has the least error, the solution is rejected and the time step is refined. Otherwise, if
$q = m^{k}$
has the least error, we check
$q = m^{k}+1$
. If
$q = m^{k}$
is confirmed as a local minimum, the solution is accepted. If not, we compute the error for
$q = m^{k}+2$
and either accept
$q = m^{k}+1$
as the solution (if it is a local minimum) or reject it and refine the time step before recalculating.
We verify whether any
$q$
causes the drop interface to intersect the substrate outside the region where pressure is allowed to be non-zero. This check is performed for all points in the angular mesh with
$\theta _q\lt \theta _i\lt \pi /2$
. If an intersection occurs, the error in the discrete optimisation is set to
$\infty$
. To confirm convergence in our results, we plot the relative change in the coefficient of restitution
$\tilde {\varepsilon } = \textrm {log} ( |\varepsilon _L-\varepsilon _{L/2}|/\varepsilon _{L/2} )$
for a typical rebound against number of modes log(
$L$
). A power-law trend with slope
$m = - 1.516$
describes the convergence our results and motivates our choice of
$L=90$
throughout this work. A flowchart for this procedure is presented in figure 12.
Open access
