1. Introduction
Droplet impacts are integral to a diverse range of natural and industrial processes (Mohammad Karim Reference Mohammad Karim2023; Hu et al. Reference Hu, Chu, Shan, Wu, Dong and Wang2024). For example, raindrop impacts on liquid and compliant biological surfaces facilitate pathogen transport (Wu et al. Reference Wu, Basu, Kim, Sorrells, Beron-Vera and Jung2024; Xu, Zhang & Lv Reference Xu, Zhang and Lv2024; Shen et al. Reference Shen, Kulkarni, Jamin, Popinet, Zaleski and Bourouiba2025). Here, we present new insights from experiments and modelling of low-Weber-number impacts of submillimetric droplets on liquid baths, a regime governing the terminal bouncing and coalescence of agricultural and biological sprays.
Since the seminal work of Worthington (Reference Worthington1882), researchers have sought to understand droplet impacts on a wide range of surfaces (Rein Reference Rein1993; Josserand & Thoroddsen Reference Josserand and Thoroddsen2016). The majority of these efforts have focused on moderate and high-Weber-number impacts, which can result in spreading, retraction or splashing (Yarin Reference Yarin2006) for impacts on solid surfaces, and also corona splashing and coalescence on liquid surfaces (Thoroddsen, Etoh & Takehara Reference Thoroddsen, Etoh and Takehara2008; Kavehpour Reference Kavehpour2015; Burzynski, Roisman & Bansmer Reference Burzynski, Roisman and Bansmer2020). Notably, whether a droplet coalesces with a bath during impact depends upon whether the thin gas layer separating the liquid bodies can be sufficiently drained during their interaction (Couder et al. Reference Couder, Fort, Gautier and Boudaoud2005; Yarin Reference Yarin2006; Kavehpour Reference Kavehpour2015). Thus, at high Weber numbers, coalescence is observed readily due to rapid gas layer drainage, while at moderate Weber numbers, impacting droplets partially or fully rebound (Burzynski et al. Reference Burzynski, Roisman and Bansmer2020; Alventosa, Cimpeanu & Harris Reference Alventosa, Cimpeanu and Harris2023), provided that the separating gas film does not thin below a critical thickness (Chubynsky et al. Reference Chubynsky, Belousov, Lockerby and Sprittles2020; Sprittles Reference Sprittles2024). Interestingly, while reducing the bath height leads to substrate-independent rebound properties (Sanjay et al. Reference Sanjay, Lakshman, Chantelot, Snoeijer and Lohse2023b ; Gabbard et al. Reference Gabbard, Agüero, Cimpeanu, Kuehr, Silver, Barotta, Galeano-Rios and Harris2025), vertically oscillating the bath can replenish the intervening gas layer and sustain droplet bouncing (Couder et al. Reference Couder, Fort, Gautier and Boudaoud2005).
Droplets bouncing on vibrating liquid baths have attracted significant interest due to their dynamics showing quantum-like behaviour reminiscent of pilot-wave theory (Bush & Oza Reference Bush and Oza2020). More recently, droplets bouncing on baths driven at two superimposed frequencies have been shown to self-propel through interaction with their own wavefields, giving rise to collective behaviours and the emergence of ‘superwalkers’, droplets larger than the maximum stable size on single-frequency baths (Valani, Slim & Simula Reference Valani, Slim and Simula2019). Galeano-Rios, Milewski & Vanden-Broeck (Reference Galeano-Rios, Milewski and Vanden-Broeck2019) derived a model for superwalkers with detailed features for the wave field. Their work treated the bouncing droplet as a rigid sphere and only partially captured the dynamics of superwalkers, suggesting that droplet deformation plays an integral role.
For droplet impact on non-vibrating baths, whether a droplet bounces or coalesces depends on the fluid properties, drop size, impact speed and bath depth (Couder et al. Reference Couder, Fort, Gautier and Boudaoud2005; Yarin Reference Yarin2006; Kavehpour Reference Kavehpour2015; Tang et al. Reference Tang, Saha, Law and Sun2019). For fixed droplet size and fluid properties, the bouncing regime is bounded by two transitions: a low-speed transition from coalescence to bouncing and a high-speed transition from bouncing back to coalescence (Pan & Law Reference Pan and Law2007; Zhao, Brunsvold & Munkejord Reference Zhao, Brunsvold and Munkejord2011). At moderate Weber number, the contact time between the droplet and bath, and the coefficient of restitution, which quantifies the translational energy recovery, are approximately constant, but become dependent on impact speed at low Weber numbers (Jayaratne & Mason Reference Jayaratne and Mason1964; Zhao et al. Reference Zhao, Brunsvold and Munkejord2011; Zou et al. Reference Zou, Wang, Zhang, Fu and Ruan2011; Wu et al. Reference Wu, Hao, Lu, Xu, Hu and Floryan2020; Sanjay et al. Reference Sanjay, Chantelot and Lohse2023a ), consistent with trends observed for droplet impact on non-wetting rigid substrates (Gabbard et al. Reference Gabbard, Agüero, Cimpeanu, Kuehr, Silver, Barotta, Galeano-Rios and Harris2025). Although bath-depth effects can alter rebound metrics (Pan & Law Reference Pan and Law2007; Zou et al. Reference Zou, Wang, Zhang, Fu and Ruan2011) and the early-time dynamics of the air film have been shown to vary with impact conditions (Thoroddsen et al. Reference Thoroddsen, Etoh and Takehara2008), we focus here on impacts on deep baths, which have been shown to yield robust rebound properties (Alventosa et al. Reference Alventosa, Cimpeanu and Harris2023). While these studies capture key rebound scalings at moderate Weber numbers, low-Weber-number impacts and models capable of resolving the droplet and bath deformation remain comparatively sparse.
Fully resolving the dynamics of a deformable impactor is challenging, as numerical multiscale effects typically span three orders of magnitude (Couder et al. Reference Couder, Fort, Gautier and Boudaoud2005), and natural time scales are short (see figure 3b of Galeano-Rios, Milewski & Vanden-Broeck Reference Galeano-Rios, Milewski and Vanden-Broeck2017). Early attempts to address these difficulties relied on asymptotic analysis and linearised formulations of the free surface equations of motion for impacts against inviscid, incompressible fluids (Wagner Reference Wagner1932; Howison, Ockendon & Wilson Reference Howison, Ockendon and Wilson1991; Korobkin Reference Korobkin2004). A weakly viscous framework introduced by Lamb (Reference Lamb1895), and revisited by Dias, Dyachenko & Zakharov (Reference Dias, Dyachenko and Zakharov2008), was combined with a formulation accounting for drop-bath coupling through the kinematic match (KM) model (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). This approach imposes natural geometric and kinematic constraints to solve for the evolution of the contact area between the impacting bodies. The KM model does not take into account the thin gas film, replacing it with an infinitely thin contact surface that transmits the pressure between the fluid bodies. Although providing remarkable agreement with experiments in the weakly viscous, low-to-moderate Weber number regime, existing implementations have been restricted to rigid impactors.
Recently, Alventosa et al. (Reference Alventosa, Cimpeanu and Harris2023) introduced a simplified version of the KM framework that allowed both impacting bodies to deform through orthogonal mode decompositions, enforce contact at a single point and assume a polynomial pressure profile at contact. This 1-point kinematic match (1PKM) model predicts droplet rebound dynamics near the inertio-capillary regime, including trajectories, coefficient of restitution, contact times and maximum deflection, in close agreement with experiments and direct numerical simulation (DNS). Phillips & Milewski (Reference Phillips and Milewski2024) and Phillips, Cimpeanu & Milewski (Reference Phillips, Cimpeanu and Milewski2025) solved directly for the pressure at contact by accounting for the gas film mediating rebound between the droplet and bath. These models are particularly well suited for weakly viscous flows impacted at low speeds, and were validated against direct numerical simulations and rigid-body impacts. Other reduced-order models have incorporated droplet deformation by representing the droplet as a spring–mass system (Blanchette Reference Blanchette2016, Reference Blanchette2017), initially using a single vertical spring coupled to a linear bath response and later an octahedral spring network calibrated to the drop’s natural oscillation modes.
In this manuscript, we develop a KM model from first principles that allows both the droplet and the bath surface to deform upon impact. Previous KM-based approaches treated the droplet as approximately rigid (Galeano-Rios et al. Reference Galeano-Rios, Milewski and Vanden-Broeck2019) and, while successful in predicting integral rebound quantities, such as the coefficient of restitution, contact time and surface pressure profiles, remained limited when droplet deformation became non-negligible. Building on the weakly viscous bath formulation of Galeano-Rios et al. (Reference Galeano-Rios, Cimpeanu, Bauman, MacEwen, Milewski and Harris2021) and Alventosa et al. (Reference Alventosa, Cimpeanu and Harris2023), the present model imposes full kinematic match conditions, allowing droplet deformation through a spectral representation of the surface of the droplet in spherical coordinates using Legendre polynomials, as done by Alventosa et al. (Reference Alventosa, Cimpeanu and Harris2023) and Phillips & Milewski (Reference Phillips and Milewski2024). This new framework yields detailed predictions for the evolution of the pressure field and pressed area without explicitly resolving the intervening gas layer. We also extend experimental measurements to submillimetric drops impacting at velocities approaching the lower rebound limit, reporting contact time and coefficient of restitution. These measurements fill a previously sparse region of parameter space and provide stringent benchmarks for models of droplet–bath impacts.
The manuscript is organised as follows. Section 3 states the governing equations, the simplifications and the boundary conditions. Section 4 presents a detailed description of the discretised system of equations to be solved numerically, complemented by the details of the numerical implementation in Appendices A and C. Section 5 presents the main results and comparisons with previous modelling efforts, including experimental results from the literature (Alventosa et al. Reference Alventosa, Cimpeanu and Harris2023), and new experiments performed in the low-
$W\kern -0.15em e$
regime. Finally, § 6 discusses our main results and contextualises them using available literature.
2. Experimental set-up
We performed experiments using submillimetre-sized silicone oil droplets. The droplets were produced by a piezoelectric droplet generator mounted on a linear stage, allowing precise control over droplet size and release (Harris, Liu & Bush Reference Harris, Liu and Bush2015; Ionkin & Harris Reference Ionkin and Harris2018), as shown in figure 1(
$a$
). Details on the generator design and droplet formation process are available in previous work (Gabbard et al. Reference Gabbard, Agüero, Cimpeanu, Kuehr, Silver, Barotta, Galeano-Rios and Harris2025). The radius of the droplet
$R_{d}$
and its impact velocity
$V_{0}$
were controlled by adjusting the nozzle size and its initial height above the bath. Both the droplet and bath were the same liquid: silicone oil with density
$\rho = \rho _{d}= 870$
kg m
$^{-3}$
, dynamic viscosity
$\mu =\mu _{d} = 0.00174$
g cm
$^{-1}$
s
$^{-1}$
and surface tension
$\sigma = \sigma _{d} = 18.7$
g s
$^{-2}$
. All experiments were performed in air at room temperature.
(
$a$
) A rendering of the experimental set-up. (
$b$
) Height of the centre of mass
$h$
against time
$t$
for a bouncing drop, where
$W\kern -0.15em e _d\,=1.285$
,
$O\kern -0.04em h_d\,=0.0238$
and
${B\kern -0.04em o\,}=0.049$
. The dashed lines are parabolic fits to the incoming and outgoing data and the dotted line denotes the drop radius
$R_{d} = 0.328$
mm. (
$c$
) Select experimental images corresponding to panel (
$b$
). The dashed line indicates the height of the undisturbed free surface at the point of contact.
The experiment begins by filling a clear, rectangular container (68 mm
$\times$
68 mm
$\times$
37 mm) with silicone oil slightly beyond brimful conditions, ensuring unobstructed side-view imaging of the droplet’s contact with the bath. The connecting tubing is purged to eliminate trapped air (Harris et al. Reference Harris, Liu and Bush2015) and the pulse width of the piezoelectric disk is tuned to produce an upward-moving droplet, allowing interfacial oscillations to subside before impact (Gabbard et al. Reference Gabbard, Agüero, Cimpeanu, Kuehr, Silver, Barotta, Galeano-Rios and Harris2025). For the lowest
$W\kern -0.15em e$
experiments, a subsequent droplet impact was required. In these cases, both droplet and bath deformations were observed to decay before the next impact, and any effect of remnant deformations on rebound metrics was smaller than the experimental uncertainty. The nozzle is then positioned approximately 1 mm above the free surface to release a droplet, which falls under gravity and impacts the bath, resulting in rebound, coalescence or transient floating. Herein, we focus on rebounds, which are defined as a rebounding droplet whose south pole rises above the undisturbed interface. For each nozzle height, three to six trials are conducted before systematically adjusting the height to achieve the desired impact velocity. This process is repeated until rebound ceases, which occurs at both sufficiently low and high impact velocities.
The droplet silhouette was captured using a Phantom Miro LC311 high-speed camera with a Laowa 25 mm ultra-macro lens. Most experiments were recorded at 15 000 fps, except for those with the smallest droplets, which required a higher frame rate of 39 000 fps due to their shorter contact times. Spatial calibration was performed using a microscope calibration slide, yielding an average resolution of 0.00469 mm pixel−1. Droplet shape and trajectory were extracted from the recordings using MATLAB. The droplet’s centre of mass was calculated in each frame by assuming a surface of revolution about the vertical axis. The droplet height
$h$
was defined as the vertical distance from the centre of mass to the undisturbed bath surface. Figure 1(
$b$
) plots
$h$
versus time
$t$
for a representative experiment, with blue and gold markers showing the incoming and outgoing trajectories, respectively. Time
$t=0$
marks the first visible contact between the droplet and bath, and the grey region indicates the contact period during which
$h$
cannot be reliably determined. The droplet radius
$R_{d}$
is determined using the projected area of the droplet in the frames prior to contact and assuming the projected shape is a circle. The drop velocity at impact
$V_{0} = V (t=0 )$
was determined by fitting a parabola to
$h (t )$
for
$t\lt 0$
and evaluating its slope at
$t=0$
. The rebound velocity
$V_{r}=V (t=t_{c} )$
was computed by fitting a parabolic fit to
$h (t )$
for
$t\gt t_{c}$
and evaluating its slope at
$t=t_{c}$
, where
$t_{c}$
is the contact time. The uncertainty of each measurement was determined using the spatial and temporal resolution of our experimental set-up, corresponding to the width of a pixel and time between successive frames, respectively. Uncertainty propagation was performed using Mathematica, where first-order Taylor expansion was used for uncertainty propagation and uncertainties were assumed uncorrelated.
Figure 1(
$c$
) shows experimental images at select moments: during free fall (
$t=-1.4$
ms), initial contact with the bath (
$t=0$
ms), maximum vertical compression during impact (
$t=1.3$
ms), liftoff from the surface (
$t=6.8$
ms) and during rebound (
$t=8.9$
ms). Liftoff is defined as the first frame where the droplet’s south pole rises above the undisturbed bath surface. Due to post-impact oscillations, the droplet is typically non-spherical at liftoff. Therefore, in addition to measuring the rebound velocity
$V_{r}$
, we also extract the height of the centre of mass of the drop at rebound
$h_{r}=h (t=t_{c} )$
.
3. Problem formulation
We consider the three-dimensional incompressible flow of a fluid bath of infinite depth and infinite lateral extension, which is initially at rest. The fluid has density
$\rho$
, kinematic viscosity
$\nu$
and surface tension
$\sigma$
. At time
$t=0$
, the undisturbed free surface of the bath is subject to a non-coalescing normal impact by a droplet of density
$\rho _d$
, kinematic viscosity
$\nu _d$
and surface tension
$\sigma _d$
. The interface of the droplet is assumed to be spherical at the moment of first contact, having a radius
$R_d$
and zero relative surface velocity with respect to the centre of mass of the droplet.
We use cylindrical coordinates
$(r,\varphi ,z,t)$
with the
$z$
pointing normally to the undisturbed free surface and away from the bath (i.e. gravity is given by
$-g\hat {\boldsymbol{z}}$
), with
$z = 0$
corresponding to the undisturbed free surface, as illustrated in figure 2. The axial symmetry of the problem around the
$z$
axis allows us to ignore the dependence on
$\varphi$
, as is henceforth done.
3.1. The bath model
Schematic of the problem. An undeformed droplet (thick solid black line) of radius
$R_d$
and fluid properties
$(\sigma _d, \rho _d, \nu _d)$
is depicted above the surface of a fluid bath with properties
$(\sigma , \rho , \nu )$
, shown with a solid grey line. The surface of the droplet is described in a non-inertial spherical reference frame by
$\xi '(t, \theta )$
, whose origin is the centre of mass of the droplet. The height of the bath’s surface, denoted by
$\eta (r, t)$
, is described in a fixed cylindrical reference frame whose origin coincides with the initial point of impact.
We model the bath using a linearised, quasi-potential flow approximation, initially presented by Lamb (Reference Lamb1895) and revisited by Dias et al. (Reference Dias, Dyachenko and Zakharov2008), and follow closely the notation and formulation presented by Galeano-Rios et al. (Reference Galeano-Rios, Milewski and Vanden-Broeck2017). This fluid model assumes that the free-surface elevation can be described by
$z = \eta (r, t)$
and introduces a velocity potential
$\phi (r, z, t)$
, as well as the pressure on the free surface due to the effect of the impacting droplet
$p(r, t)$
. We further assume that the integrity of the separating air layer is maintained throughout the impacts and simplify our model by assuming that this air layer is infinitely thin and serves only to transfer pressure between the drop and the bath. This treatment is justified in the low-Weber-number regime as a limiting case of lubrication theory where the pressure variation across the film is negligible to leading order, an assumption supported by the stable rebounds observed in our experiments. Furthermore, the impact is assumed to be perfectly non-wetting, so that the contact angle where the free surface detaches from the drop is constant and equal to
$180^\circ$
at all times during contact.
We define
$V_0$
to be the impacting velocity of the droplet with the bath at
$t=0$
, when contact is initiated (i.e. the centre of mass of the droplet is at a distance
$R_d$
above the free surface).
Taking
$R_d$
,
$T_d = \sqrt {\rho _d R_d^3 / \sigma _d}$
and
$\rho _d R_d^3$
to be the characteristic length, time and mass, respectively, we adopt the following dimensionless numbers:
\begin{align} Re = \frac {R_d^2}{T_d \nu }, \quad Bo = \frac {gT_d^2}{R_d}, \quad O\kern -0.04em h_d\, = \frac {\nu _d T_d}{R_d^2}, \quad W\kern -0.15em e _d\, = \frac {\rho _d V_0^2 R_d}{\sigma _d}, \quad D\kern -0.04em r\, = \frac {\rho _d}{\rho }, \quad Sr = \frac {\sigma _d}{\sigma }.\\[-10pt]\nonumber \end{align}
that is, the Reynolds, Bond, Ohnesorge and Weber numbers, the density ratio (similar to an Atwood number), and the ratio of the surface tensions of the two fluids. The non-dimensional quasi-potential formulation (Galeano-Rios et al. Reference Galeano-Rios, Milewski and Vanden-Broeck2017) results in the following set of equations:
\begin{align} &\qquad\qquad\partial _t \eta =\frac {2}{\textit{Re}} \varDelta _H \eta + \partial _z \phi , \quad z = 0; \end{align}
\begin{align} \phi _t &= -{B\kern -0.04em o\,} \eta + Sr\, Dr\, \kappa [\eta ] + \frac {2}{\textit{Re}} \varDelta _H \phi - Dr \,p , \quad z = 0; \end{align}
which are subject to
\begin{equation} \eta , \phi , |\boldsymbol{\nabla }\phi | \to 0 \quad \text{as}\; \; r\to \infty \quad \text{for}\;\ z = 0, \end{equation}
and
\begin{align} \eta (r,0) &= 0\quad \forall\; r\geqslant 0, \end{align}
\begin{align} \phi (r,z,0) &= 0\quad \forall\; r\geqslant 0,\; z\leqslant 0. \end{align}
Here,
$\varDelta _H = \partial _{rr} + (1/r)\,\partial _r$
denotes the two-dimensional Laplacian for axially symmetric functions and
$\kappa$
is twice the mean-curvature operator; after linearisation,
$\kappa$
is reduced to
$\varDelta _H$
. The function
$\phi$
is harmonic on its domain (
$\Delta \phi = 0,\ z \leqslant 0$
), thus allowing for the definition of the Dirichlet-to-Neumann operator for the Laplace problem in said domain. This operator is given by
\begin{equation} \partial _z \phi (r, z = 0, t) = N \varPhi (r, t) = \frac {1}{2\pi } \lim _{\epsilon \to 0^+} \int _{\mathbb{R}^2 \setminus B(r, \epsilon )} \frac {\varPhi (r, t) - \varPhi (s, t)}{|r-s|^3} \,{\rm d}A(s) \end{equation}
(see Galeano-Rios et al. Reference Galeano-Rios, Milewski and Vanden-Broeck2017), in which, the correspondence
$\{(x, y, 0) \in \mathbb{R}^3\} \equiv \mathbb{R}^2$
has been made and
$\varPhi (r,t):= \phi (r,z=0,t)$
. Therefore, (3.7), together with axial symmetry, reduces this three-dimensional problem to a one-dimensional configuration.
3.2. The droplet model
The droplet is subject to gravitational forces and to the pressure distribution from the air layer that separates the bath and droplet in the non-coalescing impact. The pressure on the droplet is assumed to be identical to the pressure on the bath, as would be expected in the lubrication limit when the intervening air layer is thin.
3.2.1. The centre of mass of the droplet
The non-dimensional height
$h(t)$
of the centre of mass of the droplet is therefore governed by
\begin{align}\dot {h}(t) &= v(t),\end{align}
\begin{align}\dot {v}(t) &= -\textit{Bo} + \frac {3}{4\pi } \int _{C(t)} p(r,t)\, {\rm d}A,\end{align}
which are subject to
\begin{equation} h(0) = 1 \end{equation}
and
\begin{equation} v(0) = -\sqrt {{\textit{We}}_d\,}, \end{equation}
where
$C(t)$
is the orthogonal projection of the pressed surface
$S(t)$
onto the
$z = 0$
plane (see figure 3) and the same units as with the bath model were used for non-dimensionalisation. Moreover, we assume that
$C(t)$
is simply connected; and, therefore, as follows from the axial symmetry assumptions, is a disk of radius
$r_c(t)$
.
Schematic of deformations during impact. The surface of the bath is shown with thin grey solid lines outside the pressed surface
$S(t)$
and with a thick grey solid line in the pressed surface
$S(t)$
. The droplet interface is shown with a thick black line. The orthogonal projection of
$S(t)$
onto the
$z = 0$
plane is
$C(t)$
, shown with a thick dark grey dashed line. Variables
$h(t)$
,
$\eta (r, t)$
and
$r_c(t)$
correspond to the height of the centre of mass of the droplet, the elevation of the free-surface of the bath and the radius of
$C(t)$
, respectively. The origin of the
$(x',z')$
system of reference is attached to the centre of mass of the droplet. The separation between the droplet and the bath over the pressed surface is introduced solely for the purpose of better visualisation of the two, as droplet and bath interfaces are predicted to fully coincide within
$S(t)$
.
3.2.2. Droplet deformation
To analyse the droplet deformation, we use spherical coordinates
$(\xi ,\theta ,\tau )$
whose origin is at the centre of mass of the droplet. Figure 3 shows that the spherical coordinates are based on the axes
$x'$
and
$z'$
, which were chosen to be respectively parallel to the
$r$
and
$z$
axes from the cylindrical coordinates that are used to describe the fluid of the bath. That is,
\begin{equation} x' = r,\quad z' = z-h(t),\quad \xi ' = \sqrt { r^2+\left (z'\right )^2},\quad \theta = \arccos \left (\frac {z'}{\xi '}\right ). \end{equation}
Moreover, the axial symmetry of the problem implies that we can ignore the dependence on the azimuthal angle
$\tau$
.
We assume that the radial coordinate that describes the shape of the drop
$\xi '$
can be written in terms of a perturbation away from a spherical shape,
$\xi '(\theta ,t) = R_d + \zeta (\theta ,t)$
. We also adopt the hypothesis of small deformations, that is, terms of the order of
$\zeta ^2$
are neglected, and we assume
$\partial _{\theta } \zeta \ll 1$
. We can thus express the surface perturbation as
\begin{equation} \zeta (\theta , t) = \sum \limits _{l = 2}^{\infty } \mathcal{A}_l(t) P_l(\cos (\theta )), \end{equation}
where
$P_l$
is the
$l$
th-order Legendre polynomial. These polynomials form a complete set of orthogonal functions defined on the interval
$[-1, 1]$
. The assumptions of small deformations and fluid incompressibility imply that the term
$l =0$
does not appear in (3.13). The term
$l=1$
is also absent from (3.13) because the centre of mass is fixed at
$z'=0$
.
We also decompose the pressure distribution due to contact, obtaining
\begin{equation} p\left (\theta ,t\right ) = \sum _{l = 0}^{\infty } \mathcal{B}_l(t) P_l\left (\cos (\theta )\right ). \end{equation}
We note that all terms
$l\geqslant 0$
appear in the sum in this case. We recover
$B_l$
from
$p(\theta , t)$
via
\begin{equation} \mathcal{B}_l(t) = \frac {1}{2l + 1} \int _{0}^{\pi } p(\theta , t) P_l(\cos (\theta )) \sin (\theta )\, {\rm d}\theta . \end{equation}
The surface modes satisfy the following set of equations (Lamb Reference Lamb1895; Blanchette Reference Blanchette2016; Alventosa et al. Reference Alventosa, Cimpeanu and Harris2023; Phillips & Milewski Reference Phillips and Milewski2024):
\begin{align}\dot {\mathcal{A}}_l &= \mathcal{U}_l \quad \text{ for } \ l \geqslant 2,\end{align}
\begin{align}R_d^3 \rho _d \dot {\mathcal{U}}_l &= - \sigma _d l(l+2)(l-1) \mathcal{A}_l- 2 R_d \mu _d (2l+1)(l-1)\mathcal{U}_l-R_d^2l\mathcal{B}_l \quad \text{ for }\, \ l \geqslant 2,\end{align}
where
$\mu _d = \rho _d \nu _d$
is the dynamic viscosity of the droplet. Observe that, when the pressure terms are zero in (3.17), the uncoupled surface modes of the drop are damped harmonic oscillators with natural frequencies
$w_l^2 = l(l-1)(l+2) \sigma _d / ( \rho _d R_d^3 )$
(Lamb Reference Lamb1895). This motivates the choice for characteristic time,
$T_d$
, in (3.1), as it has the same order of magnitude as the period of the first modes of oscillation.
Using the same characteristic units as for the bath equations, the dimensionless form of (3.16) and (3.17) is
\begin{align}\dot {A}_l &= U_l \quad \text{ for } \ l \geqslant 2,\end{align}
\begin{align}\dot {U}_l &= - l(l+2)(l-1) A_l-2O\kern -0.04em h_d\, (2l+1)(l-1)U_l-lB_l \quad \text{ for } \ l \geqslant 2,\end{align}
where
$A_l = \mathcal{A}_l/R_d$
,
$U_l=\mathcal{U}_lT_d/R_d$
,
$B_l = \mathcal{B}_l R_d/ \sigma _d$
, which are subject to
\begin{equation} A_l(0) = 0,\ U_l(0) = 0 \quad \text{ for } \ l \geqslant 2. \end{equation}
We note that (3.16) and (3.17) are valid in the weakly viscous limit, which requires the condition
$O\kern -0.04em h_d\, \sqrt {l} \lt 0.03$
(Moláček & Bush Reference Moláček and Bush2012). We satisfy this condition by truncating the Legendre series in a way that preserves the inequality; see table 1.
Parameter ranges explored in this manuscript.
3.3. The contact model
We define
$r_M(t)$
as the smallest distance from the
$z$
axis to the point on the surface of the droplet for which the plane that is tangent to the interface of the droplet is vertical. That is,
$r_M(t)$
is an upper bound for the maximum possible contact radius that a given drop shape can have when in contact with the fluid.
We can then introduce the function
$z_d(r, t)$
such that
$h(t) + z_d(r, t)$
describes the vertical height of the bottom section of the surface of the droplet for all
$r\leqslant r_{M}(t)$
, and we set
$z_d(r, t) = \infty$
for all other radii. Points on the axisymmetric surface of the drop are parametrised by
$\xi (\theta , t) = \xi '(\theta , t)/R_d$
and
$\theta$
; where
$\theta \in [ 0, \pi ]$
and
$t\geqslant 0$
. Thus the function
$z_d$
is given by
\begin{equation} z_d(r,t) = \begin{cases} \displaystyle \min _{\theta }\bigl \{\cos (\theta )\,\xi (\theta ,t)\colon \sin (\theta )\,\xi (\theta ,t)=r\bigr \}, & r \leqslant r_M(t), \\ \infty , & r \gt r_M(t). \end{cases} \end{equation}
Given the above-mentioned definition, we can express our KM constraints as
\begin{align} p(r, t) &= 0, \qquad r \gt r_c(t),\; t \geqslant 0,\end{align}
\begin{align}\eta (r, t) &= h(t) + z_d(r, t), \quad r \leqslant r_c(t),\; t\geqslant 0,\end{align}
\begin{align} \eta (r, t) &\lt h(t) + z_d(r, t), \quad r \gt r_c(t),\; t\geqslant 0\end{align}
and
\begin{equation} \partial _r \eta (r, t) = \partial _r z_d(r, t), \quad \text{at} \ r = r_c(t),\; t \geqslant 0, \end{equation}
where the final constraint requires that the drop and the bath share the same tangent cone at the boundary of the contact area. In other words, the KM imposes that the contact angle must be exactly
$\pi$
.
3.4. Model summary
Equations (3.2), (3.3), (3.8), (3.9), (3.7), (3.18) and (3.19) describe the set of coupled equations that must be solved.
The axisymmetric impact of a deformable droplet against a fluid bath is thus modelled by the solution to
\begin{align}&\qquad\qquad\qquad\quad \partial _t \eta =\frac {2}{\textit{Re}} \varDelta _H \eta + N\varPhi, \quad \forall \;r\geqslant 0 ; \end{align}
\begin{align}&\qquad \partial _t\varPhi = -{B\kern -0.04em o\,} \eta + Sr \, Dr \varDelta _H \eta + \frac {2}{\textit{Re}} \varDelta _H \varPhi - Dr \,p, \quad \forall\; r\geqslant 0;\\[-11pt]\nonumber \end{align}
\begin{align}&\qquad\,\, N\varPhi (r, t) = \frac {1}{2\pi } \lim _{\epsilon \to 0^+} \int _{\mathbb{R}^2 \setminus B(r, \epsilon )} \frac {\varPhi (r, t) - \varPhi (s, t)}{|r-s|^3} \,{\rm d}A(s); \end{align}
\begin{align}&\qquad\qquad\qquad\qquad\qquad\qquad \dot {h} = v; \end{align}
\begin{align}&\qquad\qquad\qquad \dot {v} = -{B\kern -0.04em o\,} + \frac {3}{4\pi } \int _{r\leqslant r_c(t)} p(r,t) \,{\rm d}A;\\[-11pt]\nonumber \end{align}
\begin{align}&\qquad\qquad\qquad\qquad\qquad U_l = \dot {A}_l, \quad \forall\; l\geqslant 2; \end{align}
\begin{align}& \dot {U}_l = - l(l+2)(l-1) A_l-2O\kern -0.04em h_d\, (2l+1)(l-1)U_l-lB_l, \quad \forall \;l\geqslant 2; \end{align}
for all
$t\geqslant 0$
, where
\begin{equation} B_l(t) = \frac {1}{2l + 1} \int _{0}^{\pi } p(\theta , t) P_l(\cos (\theta )) \sin (\theta )\, {\rm d}\theta . \end{equation}
These equations are subject to the velocity potential and the free-surface elevation having zero first derivatives at the axis of symmetry, and to the boundary, initial and matching conditions given by
\begin{align}&\qquad\qquad\quad \eta \to 0, \quad \text{as} \ r\to \infty, \; \forall\; t\geqslant 0; \end{align}
\begin{align}& \phi \to 0,\quad |\boldsymbol{\nabla }\phi | \to 0,\quad \text{as} \ \big(r^2+z^2\big)\to \infty, \;\forall\; t\geqslant 0; \end{align}
\begin{align}&\qquad\qquad\qquad \eta (r,0) = 0, \quad \forall\; r\geqslant 0; \end{align}
\begin{align}&\qquad\qquad \phi (r,z,0) = 0, \quad \forall \;r\geqslant 0, \;z\leqslant 0; \end{align}
\begin{align}&\qquad\qquad\qquad\quad h(0) = 1; \end{align}
\begin{align}&\qquad\qquad\qquad v(0) = -\sqrt {{\textit{We}}_d\,}; \end{align}
\begin{align}&\qquad\qquad\qquad\quad r_c(0)=0; \end{align}
\begin{align}&\qquad\quad\,\, p(r, t) = 0, \quad \forall\; r \gt r_c(t),\; t\geqslant 0; \end{align}
\begin{align}&\quad \eta (r, t) = h(t) + z_d(r), \quad \forall\; r \leqslant r_c(t),\; t\geqslant 0; \end{align}
\begin{align}&\quad\, \eta (r,t) \lt h(t) + z_d(r), \quad \forall\; r \gt r_c(t),\; t\geqslant 0; \end{align}
\begin{align}&\quad \partial _r \eta (r, t) = \partial _r z_d(r, t), \quad \text{at} \ r = r_c(t),\; \forall \;t\geqslant 0; \end{align}
where
\begin{equation} z_d(r, t)= \cos (\theta ) \left (1 + \sum \limits _{l = 2}^{\infty } A_l(t) P_l(\cos (\theta )) \right )\quad\; \forall\; r \leqslant r_c(t). \end{equation}
We note that the system of (3.26) is a closed system, containing a non-linear subproblem which is to determine the contact radius
$r_c(t)$
over which the integral in (3.26e
) is calculated.
4. Discretisation and numerical approximation
We delineate herein the approximations that convert the above-mentioned continuous system into a discrete set of evolution equations. We shall truncate the series of the surface deformation and velocity series at mode
$M+1$
; hence, the summations are taken over
$l=2,\ldots ,M+1$
. Spatially, we discretise
$r$
uniformly in intervals
$\delta _r$
, denoting
$r_i=i \delta _r$
for
$i=0,\ldots K$
. At the domain boundary,
$r_{max} = K \delta _r$
, we set
$\eta (r_{max}, t) = 0$
.
We introduce an adaptive time step
$\delta ^k_t = t^{k+1} - t^{k}$
, with
$t^1 = 0$
, and consider discrete (in time and space if applicable) approximations of all variables of interest, namely
$\eta$
,
$\phi$
,
$p$
,
$h$
,
$v$
,
$r_c$
,
$B_l$
,
$A_l$
and
$U_l$
.
We use superscripts to indicate time indices and subscripts to indicate radial position. We thus define
$\eta ^k_i$
as the approximant to
$\eta (i\delta _r,t^k)$
, and similarly
$\phi ^k_i$
and
$p^k_i$
.
The vectors
$\eta ^k$
,
$\phi ^k$
,
$p^k$
in
$\mathbb{R}^{K+1}$
whose components are given by the values for
$i=0,\ldots ,K$
are represented by the omission of the subscript. We define
$h^k$
as the approximant to
$h(t^k)$
, and similarly for
$v^k$
and
$r_c^k$
.
Finally, we also define
$A^k_l$
as the approximant to
$A_l(t^k)$
, and similarly for
$B_l^k$
and
$U^k_l$
.
The vectors
$A^k$
,
$B^k$
and
$U^k$
in
$\mathbb{R}^{M}$
whose components are given by the values for
$l = 2,\ldots ,M+1$
are represented without subscripts.
Typically, we set
$\delta _r = 1/50$
to guarantee sufficient resolution under the droplet and choose
$r_{max}=100$
(noting these are non-dimensionalised with the drop radius), which is large enough to prevent waves from reaching the boundary of the domain over the course of a rebound to mimic the behaviour of a bath of infinite lateral extension.
The pressed area
$C(t)$
is a circular disk with radius
$r_c(t)$
(see figure 3). This variable is resolved to the accuracy provided by the mesh, whose discrete-time approximation is
$r_c^k = \delta _r(m(t^k) - 1/2)$
if
$m(t^k) \gt 0$
,
$r^k_c = 0$
otherwise, where the function
$m(t^k)$
takes non-negative integer values and represents the number of contact points of the discretised system. The
$1/2$
term added to the expression of
$r^k_c$
reflects the assumption that the boundary of the contact area is found exactly at the mid-point between nodes
$m(t^k)-1$
and
$m(t^k)$
.
To find the value of
$m(t^{k+1})$
in our implicit finite difference scheme, we formulate a different system of equations for each possible number of contact points
$q$
; thus generating candidate solutions parametrised by
$q$
and set
$m(t^{k+1})$
to be equal to the value of
$q$
that best satisfies the restrictions imposed by (3.28). Details of this optimisation problem are discussed in the following section.
Time derivatives are approximated using an implicit, variable-step backward difference formula of order 2 (VS-BDF2) scheme (see Celaya, Aguirrezabala & Chatzipantelidis Reference Celaya, Aguirrezabala and Chatzipantelidis2014; Denner Reference Denner2014). To this end, we introduce the variable
\begin{equation} s^k = \delta _t^k/\delta _t^{k-1}, \end{equation}
and we define
\begin{equation} a^k = \bigl(1 + 2s^k\bigr)/\bigl(1+s^k\bigr), \quad b^k = -\bigl(1 + s^k\bigr) \quad \text{and} \quad c^k = \bigl (s^k\bigr )^2/ (1 + s^k); \end{equation}
and approximate all time derivatives by
\begin{equation} \dot {f}\bigl(t^{k+1}\bigr)\approx \dot {f}^{k+1} := \frac {a^k f^{k+1}+b^k f^{k}+c^k f^{k-1}}{\delta _t^k}. \end{equation}
We use second-order accurate matrix approximations
$\boldsymbol{\unicode{x1D649}}$
and
${\boldsymbol{\varDelta }}_{\boldsymbol{\unicode{x1D643}}}$
for the spatial operators
$N$
, given by (3.7), and
$\varDelta _H = \partial _{rr} + (1/r)\,\partial _r$
. These linear operators in matrix notation have their derivations in Appendix B of Galeano-Rios et al. (Reference Galeano-Rios, Milewski and Vanden-Broeck2017). For instance, the horizontal Laplacian is approximated by
\begin{equation} \left ({\boldsymbol{\varDelta }}_{\boldsymbol{\unicode{x1D643}}} \eta ^{k+1}\right )_i := \left \{ \begin{array}{cc} \dfrac { \eta ^{k+1}_{i-1}-2 \eta ^{k+1}_i+ \eta ^{k+1}_{i+1}}{\delta _r^2} + \dfrac {\eta ^{k+1}_{i+1} - \eta ^{k+1}_{i-1}}{2 (i-1) \delta _r^2}, & i\gt 1; \\[10pt] 4 \dfrac { \eta _{i+1}^{k+1} -\eta _i^{k+1} }{ \left (\delta _{r}\right )^2}, & i=1; \end{array}\right . \end{equation}
where we note that
$\partial _r \eta = \partial _r \phi =0$
has been imposed at the axis of symmetry.
The far-field conditions (3.28a
) are also built into the approximation when
$i = K$
, by setting boundary conditions
$\eta ^{k}_{K+1} =\phi ^{k}_{K+1} = 0$
for all
$k$
.
The integral of the contact pressure that appears on (3.26e
) can be represented by a
$1\times q$
matrix
$\boldsymbol{\unicode{x1D64E}}(q)$
; i.e.
\begin{equation} \boldsymbol{\unicode{x1D64E}}(q) p^{k} \approx \int _{C(t)} p(r,t^k)\,{\rm d}A = 2\pi \int _0^{r_c} p(r, t^k) r \,{\rm d}r. \end{equation}
Similarly, the integral on (3.27) is approximated via another single-row matrix
\begin{equation} B^k_l = \boldsymbol{\unicode{x1D647}}(l) p^k \approx \frac {1}{2l+1} \int _0^\pi p(\theta , t^k) P_l(\cos (\theta )) \boldsymbol{\cdot }\sin (\theta )\, {\rm d}\theta . \end{equation}
Matrix
$\boldsymbol{\unicode{x1D647}}(l)$
transforms a vector of pressure values on the evenly spaced radial mesh of the bath into the coefficients of its (truncated) expansion into axisymmetric spherical harmonics. We highlight that matrix
$\boldsymbol{\unicode{x1D647}}(l)$
will depend on
$A^{k+1}$
, because the elevation angle
$\theta$
that corresponds to a given point in the radial mesh of the bath changes with
$A^{k+1}$
. Figure 4 illustrates this dependence. Since
$\theta$
depends nonlinearly on the mode amplitudes
$A^k$
, it is found using a Newton root-finding method.
Dependence of the angle of the point of application of the pressure on the deformation of the droplet. Panels show the angle
$\theta$
that corresponds to
$r=\delta _r$
when the droplet is (a) undeformed and (b) significantly deformed. The lower height of the centre of mass in panel (b) results in a larger value of
$\theta$
.
Nevertheless, for a given droplet shape, operator
$\boldsymbol{\unicode{x1D647}}(l)$
is linear in
$p^k$
. A description of the matrix representation of the linear operators
$\boldsymbol{\unicode{x1D64E}}(q)$
and
$\boldsymbol{\unicode{x1D647}}(l)$
is given in Appendix A. The equations for the discretised model are given in Appendix B.
5. Results
5.1. Validation
Table 1 lists the full range of physical parameters used across all simulations and experiments. We first verify the predictions of our model by simulating the impact of a water droplet of radius
$R_d = 0.035$
cm onto a water bath at velocity
$V_0 = 38$
cm s−1, corresponding to a characteristic time of
$T_d = 0.77$
ms. This scenario matches the experimental case studied by Alventosa et al. (Reference Alventosa, Cimpeanu and Harris2023), enabling a direct comparison between our model’s predictions, and their experimental and numerical results.
The first row of figure 5 shows the droplet shape at equally spaced times during impact. The droplet initially spreads upon contact, reaching a maximum contact area before surface tension drives retraction and powers rebound. The first and last panels correspond to the initial and final moments of contact between the droplet and the bath. The second row of figure 5 presents the dimensionless pressure distribution along the droplet surface as a function of azimuthal angle in spherical coordinates. This distribution resembles a step function confined to the contact region, with accuracy determined by the number of modes retained in the simulation. The oscillations observed during spreading arise from truncating the pressure-field expansion. The third row of panels in figure 5 shows the amplitude attributed to each mode of the Legendre decomposition for the dimensionless pressure distribution. The magnitude of
$B_l$
generally decreases with
$l$
, consistent with the physical interpretation that higher order oscillation frequencies are less important – an interpretation supported by prior results on droplet rebound on non-wetting substrates (Gabbard et al. Reference Gabbard, Agüero, Cimpeanu, Kuehr, Silver, Barotta, Galeano-Rios and Harris2025). However, pressure localisation at near detachment excites higher modes, as shown in the column that corresponds to
$t/T_d = 3.42$
. These higher-frequency excitations are damped in proportion to viscosity and mode number squared (see (3.26g
)), ensuring higher modes have little effect on the overall dynamics. Here,
$\theta = 0$
corresponds to the droplet’s north pole. If the coordinate system is rotated so that
$\theta = 0$
corresponds to the south pole, all coefficients with odd
$l$
would change sign, reflecting the odd symmetry of
$P_l$
. Because higher-order modes become more prominent during spreading, truncation must be chosen carefully to ensure this phase is properly resolved. A sensitivity analysis justifying the truncation level used herein is provided in Appendix C.
Simulation of a water droplet impacting a water bath at
$V_0 = 38$
cm s−1, corresponding to
$W\kern -0.15em e _d\, = 0.7, {B\kern -0.04em o\,} = 0.017, O\kern -0.04em h_d\,=0.006$
and
$M= 20$
. Columns represent snapshots of the impact. From top to bottom, the rows show the spatial reconstruction of the impact, the dimensionless pressure distribution in azimuthal spherical coordinates at the droplet’s surface and the amplitude of the pressure distribution in Legendre’s decomposition as a function of the mode number.
We go on to compare the predicted trajectories of the north pole, south pole and centre of mass of the droplet against the experimental measurements of Alventosa et al. (Reference Alventosa, Cimpeanu and Harris2023) for the same parameters. Note that the north pole of the droplet was experimentally defined as the highest point of the droplet silhouette as captured by a side-view camera, which corresponds to a unique point along the vertical centreline at low
$W\kern -0.15em e _d\,$
, but can correspond to a circumferential ring of points at higher
$W\kern -0.15em e _d\,$
. The droplet was considered to be in contact with the bath when its north pole was below
$z=2R_d$
; we adopt this criterion to allow for a direct comparison with the prior experimental data, where the exact moment of separation was difficult to determine. The south pole was determined when the droplet was in contact with the surface of the bath and not obscured by the meniscus formed between the bath and its container, as discussed by Alventosa et al. (Reference Alventosa, Cimpeanu and Harris2023). Numerically, these parameters are defined in the same way when post-processing the simulations; the only distinction is that the south pole can be tracked at all times during the simulations.
Figure 6(
$a$
) shows the vertical position of the north pole, south pole and centre of mass of the droplet, as well as the surface of the fluid bath against non-dimensional time. Predictions from the full kinematic match (KM) developed in this manuscript are shown as yellow solid lines, while the DNS simulations, 1-point-kinematic match predictions (1PKM) and experimental measurements from Alventosa et al. (Reference Alventosa, Cimpeanu and Harris2023) are shown as a dotted blue line, dash-dotted purple line and red dashed line, respectively. The full KM model captures the dynamics with close agreement to both prior models and experimental data. In particular, it reproduces the motion of the droplet’s centre of mass with excellent accuracy, an essential feature for reliably predicting rebound characteristics such as the coefficient of restitution.
Comparison between model predictions and measurements for a water droplet impacting a water bath at
$W\kern -0.15em e _d\, = 0.7, {B\kern -0.04em o\,} = 0.017, O\kern -0.04em h_d\,=0.006$
and
$M = 20$
. Results from the full KM model (yellow solid lines) are shown alongside experiments (red dashed lines), DNS simulations (blue dotted lines) and the 1-point KM model predictions (dashed purple lines), both from Alventosa et al. (Reference Alventosa, Cimpeanu and Harris2023). (
$a$
) Vertical positions of the droplet top, centre of mass and bottom, together with the bath surface. (
$b$
) Non-dimensional contact radius. (
$c$
) Non-dimensional maximum width of the droplet, all plotted as functions of non-dimensional time.
Figure 6(
$b$
) compares the predicted evolution of the instantaneous contact radius of the droplet among the full KM model, DNS simulations and 1PKM model. Overall, the methods agree reasonably well with the full KM model attaining a similar maximum contact radius as DNS, and a similar time between first and last contact when compared with the 1PKM model. Two distinct features arise in the full KM predictions. First, as the droplet approaches its maximum spread, small oscillations in the radius of the contact area are visible. These oscillations result from surface waves propagating up and down the droplet, as a consequence of the moving pressure field, and can be clearly seen in the rebound animation in supplementary movie 1 available at https://doi.org/10.1017/jfm.2026.11291. Oscillations of the location of the contact boundary and of the pressure field (see figure 7
a) are not reported in experimental observations or DNS results. In experiments, the pressure field cannot be measured directly and the definition of contact usually requires a cut-off criterion due to the intervening air layer. This criterion determines what is resolved and may prevent oscillations of the type predicted by our model from appearing in the data. DNS results were obtained by volume-of-fluid methods where interfaces are reconstructed from approximate level sets, which may also lead to damping of oscillations. Given the small slopes of the surface of the droplet at these locations, even small differences in the cut-off thickness of the air-layer can absorb the oscillations that our model predicts. We also note that some oscillations in the pressure are seen in the sharp interface lubrication layer model for solid rebounds of Phillips & Milewski (Reference Phillips and Milewski2024). More interestingly, these oscillations are not observed in the 1PKM. This discrepancy is expected, as the 1PKM model excites a more restrictive set of surface modes (namely those used in the prescribed form of the pressure field). The other noticeable difference in the prediction of the full KM is the brief lift-off period just before rebound, around
$t/T_d = 3$
. This effect also results from neglecting the intervening air layer and assuming direct liquid–liquid contact. While this simplification is consistent with previous studies (Alventosa et al. Reference Alventosa, Cimpeanu and Harris2023) and supported by the results presented here, it produces a short interval where the droplet surface rises above the bath by a height that is less than the typical air-film thickness (
${\sim} 100$
nm; Couder et al. Reference Couder, Fort, Gautier and Boudaoud2005).
(
$a$
) Evolution of the dimensionless pressure field along the pressed area in cylindrical coordinates. The red dashed line represents the extent of the numerical contact radius. (
$b$
) Evolution of the dimensionless mode amplitudes
$A_l$
for the first modes of oscillation. Line thicknesses are inversely proportional to the mode number. The shaded area shows the timespan when the droplet and bath are in contact. Both panels correspond to the same impact as in figure 5.
Figure 6(
$c$
) compares the non-dimensional maximum width of the droplet predicted by the full KM model, DNS simulations and the 1PKM model. Unlike in panel (
$b$
), where the contact radius is restricted to discrete values by the mesh, the width in the full KM model is obtained from the spectral representation of the droplet. The Newton method is used to locate the vertical tangent plane at the droplet surface, allowing the width to vary continuously. All three approaches show good agreement, particularly in the maximum width reached during the initial spreading phase. The oscillation periods predicted by the full and 1PKM models are mutually consistent, although both are slightly longer than those from DNS. Since these differences primarily influence the post-rebound droplet shape, they have little effect on accurately capturing rebound metrics, as will be shown in the next section.
A key strength of the full KM model developed here is its ability to directly predict the pressure distribution across the contact area, whereas in the 1PKM model, the spatial form of the distribution must be prescribed (Alventosa et al. Reference Alventosa, Cimpeanu and Harris2023). Figure 7(
$a$
) presents the dimensionless pressure field for the water drop impact studied in figures 5 and 6, with the
$x$
and
$y$
axes denoting time and radial distance from the axis of symmetry, respectively. The highest pressures occur near the boundary of the contact region. Importantly, the spatio-temporal evolution of the contact pressure is highly non-trivial, underscoring the advantage of the full KM model in capturing these dynamics. Figure 7(
$b$
) shows the evolution of the first ten mode amplitudes. Higher order modes are not plotted as their amplitudes were negligible compared with those presented in the figure. The first mode of oscillation
$A_2$
is dominant, having an absolute peak at
$t/T_d \approx 1$
that roughly corresponds to the time of maximum contact area. The average absolute amplitude of a mode
$A_\ell$
is approximately in inverse proportion to its index
$\ell$
.
We compare the radial pressure profiles computed by the KM model against DNS results from Alventosa et al. (Reference Alventosa, Cimpeanu and Harris2023) in figure 8. The profiles are shown for the three dimensionless times analysed by Alventosa et al. (Reference Alventosa, Cimpeanu and Harris2023): start of impact (
$t_1=0$
), maximum spreading (
$t_2=0.6$
) and retraction (
$t_3=2.1$
). The KM model demonstrates good quantitative agreement with the fully resolved DNS, capturing both the pressure magnitude and the sharp drops in pressure near the moving contact line. We highlight that
$t_1 = 0$
was defined as the time at which contact would be starting (single point of contact) by Alventosa et al. (Reference Alventosa, Cimpeanu and Harris2023) for an airless model. For this time, the KM model predicts zero pressure, but the DNS already shows an effect of the droplet presence. The overall comparison supports the case for the use of the KM method, especially when considering the simplifying assumptions introduced here, as well as the gains in computational cost with respect to DNS.
Normalised pressure distribution versus dimensionless cylindrical coordinate
$r/R_d$
along the surface of the bath. Solid lines denote the KM model and points represent DNS results from Alventosa et al. (Reference Alventosa, Cimpeanu and Harris2023). Plots correspond to three (dimensionless) times:
$t_1 = 0.0$
,
$t_2 = 0.6$
and
$t_3 = 2.1$
, which are thickness and colour coded. Vertical dashed lines indicate the contact radii extracted from the DNS. Simulation parameters correspond to
$W\kern -0.15em e _d\, = 0.7$
,
$Bo = 0.017$
and
$O\kern -0.04em h_d\, = 0.006$
as defined in figure 5.
5.2. Rebound metrics
To compare our simulation results with experimental values, we measure the contact time
$t_c$
as the difference between the times when the north pole of the droplet crosses
$z=2R_d$
during descent and rebound; that is, following the experimental convention. Moreover, the coefficient of restitution
$\alpha$
is defined as the square root of the ratio of the rebound energy
$E_r$
to impact energy
$E_0$
:
\begin{equation} \alpha = \sqrt {\frac {E_{r}}{E_{0}}} = \sqrt {\frac {g(h_{r}-R_d) + 0.5V_{r}^2}{0.5 V_{0}^2}}, \end{equation}
where
$h_r = h(t = t_c)$
represents the centre of mass height at rebound, while
$V_0$
and
$V_r$
represent the centre of mass velocities at impact and rebound, respectively. This definition sets the potential energy of the droplet at
$t=0$
to zero and allows for a direct comparison with experiments. We define the maximum penetration depth
$\delta$
as minus the lowest height of the south pole during impact, as measured from
$z=0$
, consistent with the experimental definition used by Alventosa et al. (Reference Alventosa, Cimpeanu and Harris2023).
Non-dimensional parameters as a function of
$W\kern -0.15em e _d\,$
. (
$a$
) Maximum penetration depth, (
$b$
) contact time and (
$c$
) coefficient of restitution for a droplet with the same non-dimensional parameters as in figure 5. Experimental results are red circles, predictions from the full kinematic match model are yellow stars, DNS results are blue dotted lines and predictions from the 1-point kinematic match from Alventosa et al. (Reference Alventosa, Cimpeanu and Harris2023) are purple dash-dotted lines.
We further compare predictions from the full KM model and 1-point KM model derived by Alventosa et al. (Reference Alventosa, Cimpeanu and Harris2023) with new experiments. In our experiments, we generated submillimetre-sized droplets and tracked their motion as they impact a bath of the same liquid. The drop trajectory and shape were determined from high-speed recordings, and we measured the contact time
$t_c$
and coefficient of restitution
$\alpha$
for
$W\kern -0.15em e _d\,$
down to
$10^{-2}$
. Below this threshold, droplets either coalesced with the bath or transiently floated, depending on
$B\kern -0.04em o\,$
. Here, we define either outcome as a non-bounce event, where a coalescing drop merges with the bath during impact or its initial rebound, while a transiently floating drop persists longer on the surface: it never rises during the first rebound and only coalesces during a later oscillation on the bathB.
Figure 10 shows the coefficient of restitution
$\alpha$
, dimensionless contact time
$t_{c} / T_{d}$
and dimensionless maximum penetration depth
$\delta / R_{d}$
versus
$W\kern -0.15em e _d\,$
for a range of
$B\kern -0.04em o\,$
values, indicated by colour on a logarithmic scale. Experiments are denoted by circles, while predictions from the 1-point KM and full KM model are shown as dashed and solid lines, respectively. As most experiments involve low-
$W\kern -0.15em e _d\,$
impacts, the maximum vertical deformation is obscured in experiments and only model predictions are compared. Error bars are associated with experimental uncertainties, as discussed in Appendix B.
(
$a$
) Coefficient of restitution
$\alpha$
, (
$b$
) dimensionless contact time
$t_{c} / T_{d}$
and (
$c$
) dimensionless maximum deflection
$\delta / R_{d}$
versus Weber number
$W\kern -0.15em e _d\,$
for drops impacting a bath of the same liquid. Circles are experimental results, and dashed and solid lines are predictions from the 1-point and full kinematic match models, respectively. The colour bar maps
$B\kern -0.04em o\,$
to colour on a logarithmic scale and
$O\kern -0.04em h_d\,$
ranged from 0.0234 to 0.0351 in the simulations.
In figure 10(
$a$
), three trends in
$\alpha$
are observed. When
$W\kern -0.15em e _d\, \gtrsim 1$
, the coefficient of restitution
$\alpha$
increases with
$W\kern -0.15em e _d\,$
, approaching a
$B\kern -0.04em o\,$
-independent value of
${\approx} 0.3$
, consistent with prior literature (Bach, Koch & Gopinath Reference Bach, Koch and Gopinath2004; Wu et al. Reference Wu, Hao, Lu, Xu, Hu and Floryan2020; Alventosa et al. Reference Alventosa, Cimpeanu and Harris2023). Below
$W\kern -0.15em e _d\, \approx 1$
,
$\alpha$
increases as
$W\kern -0.15em e _d\,$
decreases over a range that widens as
$B\kern -0.04em o\,$
decreases. This increasing trend ceases once
$B\kern -0.04em o\,$
effects become significant, causing
$\alpha$
to roll off to zero, i.e. drops transition from bouncing to transiently floating on the interface. This rolloff mirrors that recently reported for droplets bouncing on a rigid, non-wetting substrate (Gabbard et al. Reference Gabbard, Agüero, Cimpeanu, Kuehr, Silver, Barotta, Galeano-Rios and Harris2025), but here, both the onset of rolloff and the cessation of rebound occur at higher
$W\kern -0.15em e _d\,$
due to the deformable substrate. For the smallest
$B\kern -0.04em o\,$
experiments (dark blue), there is no roll off as bouncing drops coalesce with the bath at low
$W\kern -0.15em e _d\,$
rather than transiently float. As the kinematic match framework does not consider the air layer, and thus cannot capture coalescence, we only simulate the
$B\kern -0.04em o\,$
associated with a bounce-to-float transition in experiment. The 1PKM and full KM models capture all three trends in
$\alpha$
, with exceptional accuracy in predicting the transition from bouncing to floating. Only predictions for the 1PKM are provided at the largest
$Bo$
(yellow) as the high deformation associated with the significantly larger droplet size in this regime prevented full KM simulations from converging.
Figure 10(
$b$
) shows the dimensionless contact time
$t_{c} / T_{d}$
versus
$W\kern -0.15em e _d\,$
. For
$W\kern -0.15em e _d\, \gtrsim 1$
, the contact time is
$W\kern -0.15em e _d\,$
-independent but increases as
$W\kern -0.15em e _d\,$
decreases in the low-
$W\kern -0.15em e _d\,$
limit, similar to drop rebound on a rigid, non-wetting substrate (Gabbard et al. Reference Gabbard, Agüero, Cimpeanu, Kuehr, Silver, Barotta, Galeano-Rios and Harris2025). The contact time begins rising at higher
$W\kern -0.15em e _d\,$
as
$B\kern -0.04em o\,$
increases. The 1PKM and full KM models accurately predict these trends and show quantitative agreement except for the largest
$W\kern -0.15em e _d\,$
. Figure 10(
$c$
) compares predictions for the maximum deflection
$\delta / R_{d}$
for the 1-point and full KM models. The model predictions are nearly identical across all
$W\kern -0.15em e _d\,$
. The initial decrease in deflection transitions to nearly constant deflection in the low-
$W\kern -0.15em e _d\,$
limit, where the asymptotic value depends solely on
$B\kern -0.04em o\,$
.
6. Discussion
This is the first work to implement the kinematic match method to model a deformable impactor colliding with a liquid bath. Our results show that this method accurately reproduces our experimental results as well as experimental and numerical results found in the literature with significantly reduced computational costs compared with the direct numerical simulations of Alventosa et al. (Reference Alventosa, Cimpeanu and Harris2023). A typical run of the full kinematic match takes 4 CPU hours, representing a twelvefold reduction compared with the average DNS run time reported by Alventosa et al. (Reference Alventosa, Cimpeanu and Harris2023).
Our results demonstrate the usefulness and versatility of the kinematic match approach when coupled with linearised fluid models, particularly for problems involving complex contact dynamics in the low-Weber-number limit. This capability is especially relevant for improving models of superwalkers (Valani et al. Reference Valani, Slim and Simula2019) and other bouncing droplet phenomena where droplet deformation is essential. Furthermore, it naturally applies to problems where the liquid bath is replaced with an elastic medium (Agüero et al. Reference Agüero, Alventosa, Harris and Galeano-Rios2022). These findings, combined with the fact that the contact conditions imposed by (3.22)–(3.25) are largely model-agnostic, encourage further exploration of the KM method to many more applications, including its use with nonlinear fluid models, extension to non-axisymmetric cases and use within variational formulations (e.g. for problems that involve more complex geometry). Indeed, provided there is no actual wetting (in the case of a droplet on a solid surface or a solid impactor on a liquid surface), coalescence or droplet breakup, the approach presented here is, in principle, generalisable to almost any impact regime.
A significant insight from the present work is that the nonlinear problem of determining the pressed area between two deformable impactors can be successfully addressed by dividing it into two iterative cycles. The external cycle explores the possible shapes of one impactor, while the internal cycle determines the extent of the pressed area. A remaining avenue for exploration is to reduce the iterative nature of the solution by parametrising the geometry and solving a nonlinear system of equations in which the geometric nonlinearities are absorbed into algebraic nonlinearities of the parametrisation.
Lastly, this study offers new insights into the physics of droplet rebound on a liquid bath, enabling direct comparisons between impacts on rigid and liquid substrates at low Weber number. For small droplets on a liquid bath, the energy recovered during rebound remains nearly constant at moderate
$W\kern -0.15em e _d\,$
, but increases when
$W\kern -0.15em e _d\, \lesssim 1$
for low
$O\kern -0.04em h_d\,$
(Zhao et al. Reference Zhao, Brunsvold and Munkejord2011; Wu et al. Reference Wu, Hao, Lu, Xu, Hu and Floryan2020; Alventosa et al. Reference Alventosa, Cimpeanu and Harris2023). We further show that at even lower
$W\kern -0.15em e _d\,$
,
$\alpha$
exhibits a roll-off that depends critically on
$B\kern -0.04em o\,$
, analogous to behaviour observed for droplets impacting non-wetting substrates (Gabbard et al. Reference Gabbard, Agüero, Cimpeanu, Kuehr, Silver, Barotta, Galeano-Rios and Harris2025). However, our experiments and model reveal that a deformable substrate shifts this rolloff to significantly higher
$W\kern -0.15em e _d\,$
, suggesting that substrate compliance may limit droplet transport via bouncing. The contact time
$t_{c}$
scales with
$T_{d}$
(Zhao et al. Reference Zhao, Brunsvold and Munkejord2011; Alventosa et al. Reference Alventosa, Cimpeanu and Harris2023) similar to droplets bouncing on a rigid substrate (Richard, Clanet & Quéré Reference Richard, Clanet and Quéré2002), but shifted towards longer interaction times. At lower
$W\kern -0.15em e _d\,$
, the contact time increases, with
$t_{c} \rightarrow \infty$
for large
$B\kern -0.04em o\,$
droplets that can transiently float on the bath (Couder et al. Reference Couder, Fort, Gautier and Boudaoud2005), while for low
$B\kern -0.04em o\,$
, contact time progressively increases until it exceeds the drainage time for the intervening air film, leading to coalescence. These results provide the most detailed characterisation to date of submillimetre-sized droplet impacts in the low-
$W\kern -0.15em e _d\,$
regime. The exceptional accuracy of the kinematic match model in this limit makes it a powerful predictive tool for a range of timely applications, including the spread of cough ejecta (Bourouiba Reference Bourouiba2021) and agricultural sprays (Makhnenko et al. Reference Makhnenko, Alonzi, Fredericks, Colby and Dutcher2021) along liquid surfaces.
Supplementary movie
Supplementary movie is available at https://doi.org/10.1017/jfm.2026.11291.
Acknowledgements
The authors gratefully acknowledge the support of the Hope Street Fellowship (C.T.G.).
Funding
This work was supported by the Engineering and Physical Sciences Research Council (EPSRC) under projects EP/P031684/1 (C.A.G-R.) and EP/N018176/1 (C.A.G-R. and P.A.M.); by Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP) under grants Sprint-17/50295-8 and 2023/07076-4 (C.R.); and by the National Science Foundation under grant NSF CBET-2123371 (D.M.H.).
Declaration of interests
The authors report no conflict of interests.
Data availability statement
The data that support the findings of this study are openly available in GitHub at https://github.com/harrislab-brown/km-droplet-onto-bath.
Appendix A. Operators
$\boldsymbol{L}(\ell )$
and
$ \boldsymbol{S}(q)$
A.1. Construction of
$\boldsymbol{L}(\ell )$
In the main text, we introduced the matrix operator
$\boldsymbol{\unicode{x1D647}}(\ell )$
that projects the discrete pressure values at a given time
$p(\theta _i) :=p(\theta _i, t)$
in the droplet onto the Legendre-mode coefficients
$B_\ell$
. In this appendix, we derive
$\boldsymbol{\unicode{x1D647}}(\ell )$
, rendering
$B_\ell$
as explicit linear combinations of the nodal pressures
$p(\theta _i)$
. We begin by recalling the continuous projection
\begin{equation} B_\ell \;=\;\frac {1}{2\ell +1}\int _{0}^{\pi } p(\theta )\,P_\ell \bigl (\cos \theta \bigr )\,\sin \theta \,{\rm d}\theta . \end{equation}
We split the integral to be approximated and use the change of variables formula to get
\begin{equation} B_l = \frac {1}{2l+1} \int _0^\pi p(\theta ) P_l(\cos (\theta )) \boldsymbol{\cdot }\sin (\theta ) {\rm d}\theta = \frac {1}{2l+1} \sum _{i=0}^{q} \int _{\cos (\theta _i)}^{\cos (\theta _{i+1})} \bar {p}(x) P_l(x) {\rm d}x, \end{equation}
where
$\bar {p}(x) = p(\arccos (x))$
. Under the assumption that
$p$
is a piecewise linear function of
$x = \cos (\theta )$
on the original spherical coordinate system, we employ a piecewise linear representation,
\begin{equation} \bar {p}(x) = ax + b = p(\theta _{i+1})\frac {x - \cos (\theta _i)}{\cos (\theta _{i+1}) - \cos (\theta _i)} + p(\theta _{i})\frac {x - \cos (\theta _{i+1})}{\cos (\theta _{i}) - \cos (\theta _{i+1})} , \end{equation}
where
$a$
and
$b$
are chosen so that
$p(\cos (\theta _i))$
is equal to the prescribed query values at
$\cos (g(i \delta _r, z^{k, j})) := \cos (\theta _i)$
, where
$g$
is a function that translates from cylindrical coordinates to spherical coordinates. We can integrate each term in (A2) exactly, resulting in a linear combination of the pressure values at the endpoints of each subinterval:
\begin{equation} \int _{\cos (\theta _i)}^{\cos (\theta _{i+1})} \bar {p}(x) P_l(x) \,{\rm d}x = a \int _{\cos (\theta _i)}^{\cos (\theta _{i+1})} x P_l(x) \,{\rm d}x + b \int _{\cos (\theta _i)}^{\cos (\theta _{i+1})} P_l(x) \,{\rm d}x. \end{equation}
We compute these integrals exactly since the antiderivatives
$F_l(x) = \int P_l(x) \,{\rm d}x$
and
$G_l(x) = \int x P_l(x) \,{\rm d}x$
satisfy the following Bonnet-type recurrence relations:
\begin{align}&\qquad (2l+1) P_l(x) = \frac {\mathrm{d}}{{\rm d}x} \bigl ( P_{l+1}(x) - P_{l-1}(x)) \implies F_l = \frac {P_{l+1} - P_{l-1}}{2l+1}, \end{align}
\begin{align}& (l+1) P_{l+1}(x) = (2l+1) x P_l(x) - l P_{l-1}(x) \implies G_l = \frac {(l+1)F_l + l F_{l-1}}{2l+1} . \end{align}
The estimation of
$B_l$
then becomes
\begin{equation} B_l \approx \frac {1}{2l+1} \sum _{i=0}^{q} \left (a_i G_l\bigg |_{\cos (\theta _i)}^{\cos (\theta _{i+1})} + b_i F_l\bigg |_{\cos (\theta _i)}^{\cos (\theta _{i+1})} \right ) , \end{equation}
where
\begin{align} a_i &= \frac {p(\theta _{i+1}) }{\cos (\theta _{i+1}) - \cos (\theta _i)} + \frac {p(\theta _{i})}{\cos (\theta _{i}) - \cos (\theta _{i+1})} , \end{align}
\begin{align} b_i &= p(\theta _{i+1})\frac { - \cos (\theta _i)}{\cos (\theta _{i+1}) - \cos (\theta _i)} + p(\theta _{i})\frac { - \cos (\theta _{i+1})}{\cos (\theta _{i}) - \cos (\theta _{i+1})} . \end{align}
Approximating the Legendre modes from the pressure values thus requires computing the value of
$P_l(\cos (\theta _i))$
for all
$ l = 2, \ldots M+1$
, for all
$ i = 0, \ldots q$
. This is done in linear time as a function of
$q$
and
$M$
, the only assumption being the linearity of the pressure as a function of
$x = \cos (\theta )$
. It can be noted that recurrence relations exist for antiderivatives of the form
$\int x^n P_l(x) \,{\rm d}x$
too, thus allowing for integration of the pressure assuming piecewise polynomial behaviour of arbitrary order.
A.2. Construction of
$\boldsymbol{S}(q)$
We recall that
$\boldsymbol{\unicode{x1D64E}}(q)$
is defined to approximate the integral of the pressure in cylindrical (bath) coordinates,
\begin{equation} \boldsymbol{\unicode{x1D64E}}(q) p^{k} \approx 2\pi \int _0^{r_c} p(r, t^k) r\, {\rm d}r, \end{equation}
where pressure on the contact area between the drop and the bath has been reduced to a one-dimensional function
$p(\boldsymbol{\cdot }, t^k)$
because of our axisymmetric assumptions, and
$t^k$
is a constant. The area of integration was originally
$C(t^k)$
, a two-dimensional (2-D) flat disk, turned into a one-dimensional (1-D) integral.
We have query points
$p^k_j := p(j \delta _r, t^k)$
for
$j = 0, \ldots , q-1$
, where
$r_c = (q - 1)\delta _r + ( {\delta _r}/{2})$
. We split the integral into subintervals,
\begin{equation} \boldsymbol{\unicode{x1D64E}}(q) p^{k} \approx 2\pi \left ( \int _{(q-1)\delta _r }^{(q-1/2)\delta _r} p(r, t^k) r \,{\rm d}r + \sum _{i=0}^{q-2}\int _{i \delta _r}^{(i+1)\delta _r} p(r, t^k) r \,{\rm d}r\right ). \end{equation}
Under the assumption that the pressure distribution is a linear interpolation of the query points from
$0$
to
$r_c$
, we calculate the integral in each subinterval as a linear combination of the query points. For example,
\begin{equation} \int _0^{\delta _r} \left (p^k_0 + \frac {p^k_1 - p^k_0}{\delta _r} r\right ) r \,{\rm d}r = p^k_0 \frac {\delta _r^2}{2} + \left (\frac {p^k_1 - p^k_0}{\delta _r}\right ) \frac {\delta _r^3}{3} = \delta _r^2 \left (\frac {p^k_0}{6} + \frac {5p^k_1}{6}\right ). \end{equation}
By adding up the contribution of the integral in each subinterval, we obtain the following vector:
\begin{equation} \boldsymbol{\unicode{x1D64E}}(q) = \pi \delta _r^2\left [ \frac {1}{3}, 2, 4, 6, \ldots , 2(q-2), \frac {3}{2}(q-1) - \frac {1}{4}\right ]. \end{equation}
Appendix B. Discretised model
We formulate the system for a fixed unknown value
$q=m(t^{k+1})$
in what follows. Hence, the system of (3.26) takes the following discretised form for all
$k \geqslant 1$
:
\begin{align} {a^k {\eta }^{k+1}+b^k {\eta }^{k}+c^k {\eta }^{k-1}} &= \delta _t^k \left ( \frac {2}{\textit{Re}} {\boldsymbol{\varDelta }}_{\boldsymbol{\unicode{x1D643}}}\eta ^{k+1} + \boldsymbol{\unicode{x1D649}}\phi ^{k+1}\right )\!; \end{align}
\begin{align} {a^k {\phi }^{k+1}+b^k {\phi }^{k}+c^k {\phi }^{k-1}} &= \delta ^k_t \left (\!-\textit{Bo} \eta ^{k+1} + \textit{Sr} \, \textit{Dr} {\boldsymbol{\varDelta }}_{\boldsymbol{\unicode{x1D643}}}\eta ^{k+1}\right . \notag\\ &\quad \left . +\, \frac {2}{\textit{Re}} {\boldsymbol{\varDelta }}_{\boldsymbol{\unicode{x1D643}}}\phi ^{k+1} - Dr \,p^{k+1}\!\right )\! ; \end{align}
\begin{align} {a^k {v}^{k+1}+b^k {v}^{k}+c^k {v}^{k-1}} &= \delta ^k_t \left (-{B\kern -0.04em o\,} + \frac {3}{4\pi } \boldsymbol{\unicode{x1D64E}}(q) p^{k+1}\right )\!; \end{align}
\begin{align} {a^k {h}^{k+1}+b^k {h}^{k}+c^k {h}^{k-1}}{h} &= \delta ^k_t v^{k+1}; \end{align}
\begin{align} {a^k {A_l}^{k+1}+b^k {A_l}^{k}+c^k {A_l}^{k-1}} &= \delta ^k_t U^{k+1}_{l}, \quad \forall\; l = 2, 3, \ldots M+1; \end{align}
\begin{align} {a^k {U_l}^{k+1}+b^k {U_l}^{k}+c^k {U_l}^{k-1}}{U_l}&= \delta ^k_t \left (- l(l+2)(l-1) A^{k+1}_l \right . \notag\\ &\quad \left .-\, 2O\kern -0.04em h_d\, (2l+1)(l-1) U^{k+1}_{l}-lB^{k+1}_l \right ) \forall\; l \notag\\&= 2, 3, \ldots M+1. \end{align}
\begin{align} B^k_l = \boldsymbol{\unicode{x1D647}}(l) p^{k}, \quad \forall\; l &= 2, \ldots , M+1. \end{align}
\begin{align} p^{k}_{i} &= 0, \quad \forall\; i \gt q, \end{align}
subject to discretised initial and boundary conditions, and the discretised approximations of the system of restrictions (3.28). In particular, for all
$k\geqslant 1$
, we have
\begin{align} \eta ^k_i &= h^k + z_d(i \delta _r, t^k), \quad \forall \;i \leqslant q; \end{align}
\begin{align} \eta ^k_i &\lt h^k + z_d(i \delta _r, t^k), \quad \forall\; i \gt q; \end{align}
\begin{align} \frac {\eta ^k_{q+1} - \eta ^k_{q}}{\delta _r} &= \partial _r z_d((q - 1/2) \delta _r, t^k),\quad \text{if } q \gt 0. \end{align}
Before making further substitutions, we recast the system of (B1a )–(B1f ) into two coupled matrix forms
\begin{equation} \boldsymbol{\unicode{x1D64C}}(\delta _t^k) W^{k+1} = F, \quad \boldsymbol{\unicode{x1D64D}}(\delta ^k_t) Y^{k+1} = G; \end{equation}
where
\begin{align} \boldsymbol{\unicode{x1D64C}}(\delta _t^k) &= \left [ \begin{array}{ccccc} a^k\left (\boldsymbol{\unicode{x1D644}}- \dfrac { 2}{\textit{Re}} {\boldsymbol{\varDelta }}_{\boldsymbol{\unicode{x1D643}}}\right ) & -\delta _t^{k}\boldsymbol{\unicode{x1D649}} & 0_{K\times K} & 0_{K\times 1} & 0_{K\times 1} \\[6pt] \delta _t^{k} \left ( \displaystyle Bo \boldsymbol{\unicode{x1D644}} - \textstyle Sr\, \textstyle Dr {\boldsymbol{\varDelta }}_{\boldsymbol{\unicode{x1D643}}}\right ) & a^k \boldsymbol{\unicode{x1D644}} - \frac { \textstyle 2 \delta ^k_t}{ \textstyle Re} {\boldsymbol{\varDelta }}_{\boldsymbol{\unicode{x1D643}}} & \delta _t^{k}Dr \boldsymbol{\unicode{x1D644}} & 0_{K\times 1} & 0_{K\times 1} \\[6pt] 0_{1\times K} & 0_{1\times K} & \frac {\textstyle 3\delta ^k_t}{\textstyle 4\pi } \boldsymbol{\unicode{x1D64E}}(q) & a^k & 0 \\[6pt] 0_{1\times K} & 0_{1\times K} & 0_{1\times K} & -\delta ^k_t & a^k \end{array} \right ]\!,\end{align}
\begin{align}W^{k+1} &= \left [ \begin{array}{ccccc} \eta ^{k+1} & \phi ^{k+1} & p^{k+1} & v^{k+1} & h^{k+1} \end{array} \right ]^{\top },\end{align}
\begin{align}\boldsymbol{\unicode{x1D64D}}(\delta _t^k) &= \left [ \begin{array}{ccc} a^k\boldsymbol{\unicode{x1D644}} & -\delta _t^k\boldsymbol{\unicode{x1D644}} & \boldsymbol{\mathsf{0}} \\[6pt] \delta ^k_t\boldsymbol{\unicode{x1D64B}} & a^k \boldsymbol{\unicode{x1D644}} + \delta ^k_t \boldsymbol{\unicode{x1D64F}} & \delta ^k_t \boldsymbol{\unicode{x1D651}} \end{array} \right ]\!,\end{align}
with
\begin{align}&\qquad P_{i, j} = (j+1)(j+3)j\delta _{\textit{ij}}, \quad \text{ for }\ i, j = 1 \ldots M, \end{align}
\begin{align}&\qquad\quad\,\, T_{i, j} = (2j+3)j\delta _{\textit{ij}}, \quad \text{ for }\ i, j = 1 \ldots M , \end{align}
\begin{align}&\qquad\qquad V_{i, j} = (j+1)\delta _{\textit{ij}}, \quad \text{ for }\ i, j = 1 \ldots M, \end{align}
\begin{align}&\qquad\qquad\,\, Y^{k+1} = \left [ \begin{array}{ccc} A^{k+1} & U^{k+1} & B^{k+1} \end{array} \right ]^{\top },\\[-11pt]\nonumber \end{align}
\begin{align}& F = \left [ \begin{array}{c} - b^k \eta ^{k} - c^k \eta ^{k-1} \\[2pt] -b^k \phi ^{k} - c^k \phi ^{k-1} \\[2pt] -b^k v^k - c^k v^{k-1} \\[2pt] -b^k h^k - c^k h^{k-1} \end{array} \right ]\!, \quad G = \left [ \begin{array}{c} - b^k A^{k} - c^k A^{k-1} \\[2pt] -b^k U^{k} - c^k U^{k-1} \end{array} \right ]\!, \end{align}
where
$\boldsymbol{\unicode{x1D644}}$
and
$\boldsymbol{\mathsf{0}}$
are the identity and zero matrix of size
$M$
, respectively. We note that
$F$
and
$G$
depend on information from previous times only. The linear operators in (B3) represent the linear bath and drop dynamics, respectively. They are coupled via (B2), which describe the free surface height compatibility with respect to the drop deformation during contact. This compatibility introduces nonlinearities in the system.
We note that
$\boldsymbol{R}$
is a rectangular matrix of size
$2M \times 3M$
and
$\boldsymbol{Q}$
is of size
$(2K+2) \times (3K+2)$
, corresponding to
$2K+2$
equations and
$3K+2$
unknowns. Further manipulations are thus needed to make them square matrices.
The same procedure described by Galeano-Rios et al. (Reference Galeano-Rios, Milewski and Vanden-Broeck2017) to reduce the number of unknowns from
$3K+2$
to
$2K+2$
for the case of a rigid sphere is applied to this bath-drop matrix system. As such, we apply relations (B1h
) and (B2a
) to reduce the effective number of unknowns.
This substitution makes the right-hand side of the first equation
$A^{k+1}$
-dependent. We obtain a modified matrix system,
\begin{equation} \boldsymbol{Q}'(\delta _t^k, q) W^{\prime k+1} = F'(A^{k+1}, q) = F'. \end{equation}
Similarly, we recast the other matricial system by incorporating
$B^{k+1}$
into the left-hand side of (B3), therefore having the following formulation for the reduced system:
\begin{equation} \boldsymbol{\unicode{x1D64D}}' (\delta _t^k) Y^{\prime k+1} = G'(\delta _t^k, B^{k+1}); \end{equation}
where
\begin{align}&\qquad\qquad\qquad\qquad\quad \boldsymbol{\unicode{x1D64D}}^{\prime}(\delta _t^k) = \left [ \begin{array}{cc} a^k\boldsymbol{\unicode{x1D644}} & -\delta _t^k\boldsymbol{\unicode{x1D644}} \\[3pt] \delta ^k_t\boldsymbol{\unicode{x1D653}} & a^k \boldsymbol{\unicode{x1D644}} \end{array} \right ]\!, \\[-10pt] \nonumber \end{align}
\begin{align}&\qquad\qquad\quad X_{i, j} = (i+1)(i+3)i \ \delta _{\textit{ij}} \quad \text{ for } i, j = 1 \ldots M, \\[-10pt] \nonumber \end{align}
\begin{align}&\qquad\qquad\qquad\qquad Y^{\prime k+1} = \left [ \begin{array}{cc} A^{k+1} & U^{k+1} \end{array} \right ]^{\top },\\[-10pt]\nonumber \end{align}
\begin{align}& G'(\delta _t^k, B^{k+1}) = \left [ \begin{array}{c} - b^k A^{k} - c^k A^{k-1} \\ -b^k U^{k} - c^k U^{k-1} - \left (\begin{array}{c} 2\delta _t^k B^{k+1}_2 \\[5pt] 3\delta _t^k B^{k+1}_3 \\ \vdots \\ (M+1)\delta _t^k B^{k+1}_{M+1} \end{array} \right ) \end{array} \right ]\!.\end{align}
Equations (B12) and (B13) are therefore analogous to (3.19) and (3.6) from Agüero et al. (Reference Agüero, Alventosa, Harris and Galeano-Rios2022) and Galeano-Rios et al. (Reference Galeano-Rios, Milewski and Vanden-Broeck2017), to render them square matrices.
Further technical details regarding the computational implementation are given in Appendix C. To simplify notation, deformation variables
$B_l$
and
$A_l$
without the
$k$
superscript in subsequent sections will refer to the numerical solutions, unless stated otherwise.
Pseudocode solve_system() to find the coupled solution to (B 3) byiterating on
$B^{k+1,j}$
Appendix C. Computational implementation
C.1. An iteration on pressure
The system of (B12) and (B13) has to be solved together with the nonlinear set of equations given by (B2). We address the interdependence of
$A^{k+1}$
in both (B12) and (B13) by using the iterative procedure outlined in Algorithm1 to compute
$B^{k+1}$
.
Starting from an initial guess
$B^{k+1,1}$
, we iteratively compute
$B^{k+1,j}$
for
$j\geqslant 1$
, stopping when convergence yields the pressure coefficient vector
$B^{k+1}$
, which by construction satisfies both matrix systems at once.
The above-mentioned manipulations transform the problem of finding the solution of the next time step to finding the limit of
$B^{k+1, j}$
in the following iterative subroutine:
\begin{align} B^{k+1,1} &= B^k, \end{align}
\begin{align} Y^{k+1, j} &= \big(\boldsymbol{\unicode{x1D64D}}'(\delta _t^k, q)\big)^{-1}G'(\delta _t^k, B^{k+1,j}), \end{align}
\begin{align} W^{\prime k+1, j} &= \big(\boldsymbol{\unicode{x1D64C}}'(\delta _t^k, q)\big)^{-1}F'(\delta _t^k, A^{k+1, j}), \end{align}
\begin{align} B^{k+1,j+1}_l &= \boldsymbol{\unicode{x1D647}}(l) p^{k+1,j}. \end{align}
The set of (C6) is prescribed under the assumption that the contact area for that time step
$t^k$
contains exactly
$q$
contact points in the discretised mesh, which is itself an unknown to be found using another iterative method described later.
C.2. An iteration on geometry
The contact radius is parametrised by
$q$
and remains to be determined. We note that from the system of restrictions (B2), only (B2b
) and (B2c
) have not been verified in the previous discussion. These two conditions will determine which value for
$q$
is assigned to
$r^{k+1}_c/\delta _r$
.
By following the procedure described in Appendix C.1, each triplet
$(k, \delta _t^k, q)$
determines a set of candidate solutions for the next time step,
${}_q{\eta ^{k+1}}$
,
${}_q{\phi ^{k+1}}$
,
${}_q{p^{k+1}}$
,
${}_q{h^{k+1}}$
,
${}_q{v^{k+1}}$
,
${}_q{r^{k+1}_{c}}$
,
${}_q{B^{k+1}}$
,
${}_q{A^{k+1}}$
and
${}_q{U^{k+1}}$
.
We define the error of this potential solution by the functional
After convergence of the pressure iteration
$B^{k+1,j}\to B^{k+1}$
, condition (C7) is used to eliminate any solution where the drop and bath surfaces intersect. For the remaining candidates, the residual in (B2c
) provides a quantitative measure of the error associated with the parameters
$(\delta _t^k, q)$
. An exact partial differential equation solution would drive this residual to zero by satisfying (B2c
) exactly. We then select the number of contact points at
$t^{k+1}$
that minimises this residual to define the pressed area,
\begin{equation} r^{k+1}_c := \delta _r \boldsymbol{\cdot }\text{argmin}_q \ e(q, \delta _t^k). \end{equation}
Pseudocode to solve the discrete formulation at the next time step
C.3. Numerical implementation considerations
To guarantee a physically continuous evolution of the pressed surface
$S(t)$
, we enforce a smooth update of the discretised contact radius via
\begin{align} \bigl \lvert r^{k+1}_c - r^{k}_c \bigr \rvert \leqslant \delta _r. \end{align}
Whenever this criterion is violated, we halve the time step (
$\delta _t^k \leftarrow \delta _t^k/2$
) and discard all candidate solutions computed with the larger step. The search then restarts under the refined temporal resolution, ensuring that each update of the contact area proceeds without jumps.
The numerical implementation further ensures that surface waves do not reach the boundaries of the computational domain, preserving the validity of the imposed assumptions. Additionally, condition
$|\boldsymbol{\nabla }\eta | \lt 1$
is verified for all solutions at all times, ensuring the free surface slope remains within the linear limit consistent with the linearised curvature approximation.
As already described, time derivatives are discretised with a second-order backward difference formula (BDF2). To bootstrap BDF2, the first increment uses a first-order implicit Euler step, providing the two initial time levels needed for the higher-order scheme.
Integral operators
$\boldsymbol{S}(q)$
and
$\boldsymbol{\unicode{x1D647}}(l)$
used in (B1c
) and (B1g
) are implemented using a trapezoidal rule, as described in Appendix A.
The simulations reported here use a spectral representation of the droplet with
$M=20$
modes (3.13). This truncation level was chosen to ensure the stability of the inner pressure optimisation loop, which becomes unstable for
$M \gtrsim 65$
. To validate this choice, simulations were repeated with
$M=40$
and
$M=60$
modes, ensuring a relative change of at most 5 % in all relevant variables and yielding visually indistinguishable plots.
To sum up, the full algorithm is laid out in Algorithm2. After initialising the solution and taking the first time step, the main loop generates multiple candidates (indexed by
$q$
), evaluates them against the tangent-error metric (C7) and selects the optimal value of contact points to advance. Any overly rapid numerical change in
$r^k_c$
triggers the aforementioned time-step reduction and recomputation.
Open access
