2.1. Experimental set-up

A series of droplet impact experiments were conducted utilizing two working fluids: deionized water and silicone oil with viscosity of $5$ cSt. Ranges of experimental parameters are summarized in table 1. A drop-on-demand generator is used to reliably produce droplets with a maximum variation in the diameter of less than $1\,\%$ (Ionkin & Harris Reference Ionkin and Harris2018). This device, along with a schematic of the experimental set-up, is shown in figure 2(a). The drop generator is entirely three-dimensionally printed (3D-printed), with the exception of a small piezoelectric disk, hardware and connective tubing. The deformation of the piezoelectric disk due to an applied voltage pulse acts to expel fluid through a small nozzle. As the fluid exits the nozzle, the piezoelectric disk relaxes, initiating pinch off of the droplets. The droplets then fall under the action of gravity towards the bath. Two visualizations of droplet impact and rebound are shown in figure 2(b,c). The drop generator is mounted on a 3D-printed translation stage, allowing for repeatable changes to the impacting velocity via height increases of the droplet generator. Directly underneath the drop generator is a 3D-printed fluid bath. The bath is $70$ mm in width and length, and $50$ mm deep. The impact location was $25$ mm from the front wall of the bath. This impact location allowed for consistent focus above and below the free surface yet was still sufficiently far from the front panel that the waves created during impact do not have time to reflect and interact with the droplet during contact. For the water experiments, the front and rear walls are constructed using polystyrene that has an equilibrium contact angle of $87.4^{\circ }$ (Ellison & Zisman Reference Ellison and Zisman1954). Being close to $90^{\circ }$, this creates a negligibly small meniscus that allows for detailed photography of the impact from the side (Galeano-Rios et al. Reference Galeano-Rios, Cimpeanu, Bauman, MacEwen, Milewski and Harris2021). For the silicone oil experiments, we use a shorter bath window panel, constructed of extruded acrylic and a thin transparent plastic sheet. The bath was brim-filled to the height of the acrylic window panel, such that the contact line was pinned with angle held at approximately $90^{\circ }$. The drops are imaged using a high-speed camera (Phantom Miro LC 311) and illuminated by a Phlox LED-W back light. Movie data is taken at 10 000 frames-per-second (f.p.s.) with an exposure time of $99.6\,\mathrm {\mu }\mathrm {s}$.

Table 1.

Relevant parameters and their range of values in our experimental study.

Figure 2.

(a) A rendering of the experimental set-up. (b) Experimental montage of impact of a deionized water droplet on a bath of the same fluid. Images are spaced 0.7 ms apart. (c) Spatiotemporal diagram of a deionized water droplet bouncing. The image is constructed by taking a single-pixel-wide stripe of the raw movie footage, and plotting time along the $x$-axis. Panels (b,c) correspond to an impact of deionized water on a bath of the same fluid with ${We}=0.7$, ${Bo} = 0.017$ and ${Oh} = 0.006$. (d) Schematic of the problem.

2.2. Experimental procedure

Care must be taken to ensure that both the fluid interface and the fluid in the reservoir of the drop generator are contaminant-free, as dust or surfactants can modify the physics involved. Prior to each experiment, the bath and tubing are thoroughly cleaned with an isopropyl alcohol solution, flushed with deionized water, and then left to dry in a fume hood with particulate filtering for 30 min. The drop generator and nozzle are cleansed with an ethanol solution, and then flushed with deionized water for five minutes. Gloves are worn at all times to minimize contamination. The drop generator is controlled by an Arduino Uno board, with a simple H-bridge circuit initiating the voltage switching of a DC power supply (Ionkin & Harris Reference Ionkin and Harris2018). The fluid within the bath is periodically flushed, approximately after every 15 droplet impacts to reduce surface contamination (Kou & Saylor Reference Kou and Saylor2008). Overflow from the flushing is caught by a small lip in the bottom of the bath, which is then drained to a waste container. There are two syringes connected to the bath, which allow for fine adjustments of the equilibrium bath depth after flushing.

We collect experimental data for the top and bottom of the droplet during free flight. During contact, we track the height of the top of the droplet and that of the centre of the deformed free surface. Since the air layer that separates the droplet and the interface is negligibly thin relative to the scale of the droplet, we assume that this point is also effectively the location of the droplet's south pole. Just after the drop rebounds off the surface, the axisymmetric surface wave created by the impact is partially in the line of sight of the camera, and obscures the bottom of the droplet for a brief period during take-off. These data points have been omitted from the bottom trajectory when reported. The raw movie data are postprocessed using a custom Canny edge detection software implemented in MATLAB 2021b, which quantifies the droplet trajectory, and then computes impact parameters and bouncing metrics (Galeano-Rios et al. Reference Galeano-Rios, Cimpeanu, Bauman, MacEwen, Milewski and Harris2021).

There are several metrics of interest in our study, which we define in what follows. The maximum penetration depth, $\delta$, of a bounce is defined as the position of the bottom of the droplet at the lowest point in the trajectory (relative to the undisturbed interface height). In our experiment, the contact time, $t_c$, is defined as the time duration from which the top of the droplet crosses the height $z=2R$ to the time the top of the droplet returns to that height. Due to the nature of visualization set-up it was impossible to determine precisely when the droplets lost physical contact with the fluid; however, this always occurred before the top of the drop returned to the level $z=2R$. Each bounce was also characterized by its coefficient of restitution, $\alpha$, which is defined here as the negative of the normal exit velocity, $V_e$, divided by the normal impact velocity, $V_0$. This parameter ranges between 0 and 1, and is related to the momentum exchange during impact.

In order to determine the contact time ($t_c$) and coefficient of restitution ($\alpha$), a parabola was fitted using a least-squares method to the incoming and outgoing trajectories, separately, with at least 30 data points prior to impact and at least 40 data points following rebound. The analytical form of the parabolic fit was then used to extrapolate the time at which the sphere crosses the still air–water interface height (which corresponds to a root of the parabolic function). The derivative of the parabolic fit function was then computed analytically and its value evaluated at these times in order to calculate the impact speed, $V_0$, and exit speed, $V_e$ (Galeano-Rios et al. Reference Galeano-Rios, Cimpeanu, Bauman, MacEwen, Milewski and Harris2021).

3.1. Bath interface model

The present work models the bath interface dynamics using a linearized, quasipotential flow model following the work of Galeano-Rios et al. (Reference Galeano-Rios, Milewski and Vanden-Broeck2017). For the problem of a droplet impacting on a free interface, the Navier–Stokes equations govern the flow generated by the bath–droplet interaction. Assuming the flow to be incompressible, isothermal and Newtonian, we can define the fluid velocity vector $\boldsymbol {u} = [u,v,w]^{\rm T}= \boldsymbol {\nabla } \phi + \boldsymbol {\nabla } \times \varPsi$ and the bath interface shape $\eta =\eta (r,\varTheta,t)$, Here, $\phi$ is the scalar potential and $\varPsi$ is the vector stream function. We then linearize the governing equations and boundary conditions about the undisturbed free surface $z=0$. Utilizing the arguments presented in Galeano-Rios et al. (Reference Galeano-Rios, Milewski and Vanden-Broeck2017) and Dias et al. (Reference Dias, Dyachenko and Zakharov2008), we can recast the governing equations to be

(3.1)$$\begin{gather} \nabla^2 \phi {+ \partial_z^2 \phi} = 0,\quad z \leq 0 , \end{gather}$$

(3.2)$$\begin{gather}\partial_t \eta = \partial_z \phi + 2 \nu {\nabla^2 \eta} ,\quad z = 0 , \end{gather}$$

(3.3)$$\begin{gather}\partial_t \phi ={-} g\eta -2\nu \partial_z^2 \phi + \frac{\sigma}{\rho}\kappa - \frac{p_s}{\rho} ,\quad z = 0 , \end{gather}$$

(3.4)$$\begin{gather}\partial_z \phi = 0,\quad {\rm at}\ z=h_0 . \end{gather}$$

Here, $p_s(r,\varTheta,t)$ is the contact pressure, $g$ is the gravitational acceleration, $\kappa (r,\varTheta,t)=\nabla ^2 \eta$ is twice the linearized mean curvature of the interface, $h_0$ is the depth of the undisturbed bath and $\rho$, $\sigma$, $\nu = \mu /\rho$ are the fluid density, surface tension and kinematic viscosity, respectively. In this notation, $\partial _{( )}$ denotes partial differentiation with respect to the variable given in the parenthesis and $\nabla ^2 = \partial _{r}^2 + (1/r)\partial _r + (1/r^2)\partial _{\varTheta }^2$. The tangential stress boundary conditions are automatically satisfied in these approximations. As detailed in Galeano-Rios et al. (Reference Galeano-Rios, Milewski and Vanden-Broeck2017), this leading-order theory is valid in the weakly viscous limit when $Re=\sqrt {We}/Oh \gg 1$. A similar bath model was used by Blanchette (Reference Blanchette2016, Reference Blanchette2017), although the viscous correction term was not included in the dynamic boundary condition (3.3).

We assume that the impact occurs in a bath of some viscous fluid which is subject to two boundary conditions, $\partial _n \phi = 0$ on the walls of the bath and $\partial _z \phi =0$ on the bottom of the bath, where $n$ is the outward facing normal of the walls of the bath (Benjamin & Ursell Reference Benjamin and Ursell1954). The former condition implies that $\partial _n \eta =0$ on the walls of the container. These conditions correspond physically to a bath where the working fluid maintains a constant contact angle of $90^{\circ }$ at the walls, and has no-flux boundary conditions along the walls and bottom. Applying these conditions to the governing system of equations, an orthogonal function expansion for the unknowns $\eta$, $\phi$ and their derivatives can be explicitly written.

The orthogonal basis functions $S_{j,m}(r,\varTheta )$ ultimately must satisfy

(3.5)\begin{equation} (\nabla^2 + k_{j,m}^2) S_{j,m} = 0 , \end{equation}

inside of any bounding curve $K$ that exists on the border of the bath and the boundary condition of $\partial _n S_{j,m} = 0$ on $K$. Here $k_{j,m}$ are the eigenvalues of the system, and depend on the choice of the physical domain of the problem. If the boundary curve $K$ is a circle, the functions become $S_{j,m} = {\rm J}_j(k_{j,m} r) \cos {j \varTheta }$. Here ${\rm J}_j$ are Bessel functions of the first kind and $k_{j,m}$ are the solutions to ${\rm J}'_j(k_{j,m}b) = 0$, where $b$ is the bath radius. If we further assume the problem to be axisymmetric, then we can choose $j=0$, and write $S_{0,m} = {\rm J}_0(k_{0,m} r).$ For convenience, we define $k_{0,m}=k_m$ henceforth. We then can express the free surface elevation as

(3.6)\begin{equation} \eta(r,t) = \sum_{m=0}^{\infty} a_m(t) {\rm J}_0(k_{m} r) . \end{equation}

We then rewrite all of the unknowns of the axisymmetric bath problem, using (3.1) and (3.2), as a function of the time varying amplitude coefficients, $a_m(t)$:

(3.7)$$\begin{gather} \eta(r,t) = \sum_{m=0}^{\infty} a_m(t) {\rm J}_0(k_m r) , \end{gather}$$

(3.8)$$\begin{gather}\kappa(r,t)=\nabla^2 \eta ={-}\sum_{m=0}^{\infty} k_m^2 a_m {\rm J}_0(k_m r) , \end{gather}$$

(3.9)$$\begin{gather}\phi(r,z,t) = \sum_{m=0}^{\infty} \left( \frac{{\rm d} a_m}{{\rm d}t} + 2 \nu k_m^2 a_m \right) {\rm J}_0(k_m r) \frac{\cosh{k_m(h_0+z)}}{k_m\sinh{k_m h_0}} . \end{gather}$$

In order to arrive at the final equations of motion for the free surface, we take the decompositions ((3.7)–(3.9)) and substitute them into the dynamic boundary condition (3.3). Rearranging, we find

(3.10)\begin{align} \sum_{m=0}^{\infty} \left[\frac{{\rm d}^2 a_m}{{\rm d}t^2} + 4 \nu k_m^2 \frac{{\rm d}a_m}{{\rm d}t} \,{+} \left(\frac{\sigma k_m^2}{\rho} + g\right)k_m\tanh{(k_m h_0)}a_m \right ] {\rm J}_0(k_mr) ={-}\frac{p_s}{\rho}k_m\tanh{(k_m h_0)}. \end{align}

Each wave mode in the bath is described by a forced, damped harmonic oscillator equation.

3.2. Droplet interface model

Additionally, we wish to recover a similar set of equations that describe the gravity-capillary waves present in the droplet, and then couple these equations to the motion of the bath. The full derivation of the droplet oscillation model can be found throughout prior works (Lamb Reference Lamb1924; Tsamopoulos & Brown Reference Tsamopoulos and Brown1983; Courty, Lagubeau & Tixier Reference Courty, Lagubeau and Tixier2006; Chevy et al. Reference Chevy, Chepelianskii, Quéré and Raphaël2012; Balla, Tripathi & Sahu Reference Balla, Tripathi and Sahu2019) and is briefly summarized below.

We begin by utilizing spherical harmonics to decompose the droplet radius in a spherical domain,

(3.11)\begin{equation} \xi(\theta,t) = R + \sum_{l=1}^{\infty} \sum_{n={-}l}^l \beta_l^n(t) Y_l^n(\cos{\theta},\phi) , \end{equation}

with $Y_l^n = P_l^n(\cos {\theta }) e^{i n \phi }$. Due to the axisymmetry of the problem, we set $n=0$ and the spherical harmonics reduce to associated Legendre polynomials, $P_l^n (\cos (\theta ))$. For convenience, we write $\beta _l^n=\beta _l$ henceforth. We assume that the velocity potential takes the same form as the decomposition of the interface (Lamb Reference Lamb1924; Balla et al. Reference Balla, Tripathi and Sahu2019). We then turn to an energy conservation equation of the form

(3.12)\begin{equation} \frac{{\rm d} T}{{\rm d}t} = \dot{W} - \epsilon , \end{equation}

with $T = K + G$, $\dot {W}$ and $\epsilon$, as the total energy (sum of the kinetic energy $K$ and potential energy $G$) of the drop, the rate of work done on the droplet interface, and the viscous dissipation, respectively. We can express $T$, $\dot {W}$ and $\epsilon$ using the decomposition (3.11), and substitute these expressions into the conservation of energy equation. Then, utilizing the linearized kinematic boundary condition yields a set of forced, damped harmonic oscillators that describe the amplitude of each individual spherical mode,

(3.13)\begin{align} \sum_{l=1}^{\infty} \left[ \frac{{\rm d}^2 \beta_l}{{\rm d}t^2} + 2 \alpha_{l} \frac{{\rm d} \beta_l}{{\rm d}t} + \omega_{l}^2 \beta_l \right ] = \sum_{l=1}^{\infty} \left[ -\frac{(2l+1)l}{2\rho R} \int_0^{\rm \pi} p_s(\theta) \sin{\theta} P_l(\cos{\theta}) \,{\rm d}\theta + g \delta_{1l} \right ] . \end{align}

We drop the sums, and arrive at the result

(3.14)\begin{equation} \frac{{\rm d}^2 \beta_l}{{\rm d}t^2} + 2 \alpha_{l} \frac{{\rm d} \beta_l}{{\rm d}t} + \omega_{l}^2 \beta_l ={-}\frac{(2l+1)l}{2\rho R} \int_0^{\rm \pi} p_s(\theta) \sin{\theta} P_l(\cos{\theta}) \,{\rm d}\theta + g \delta_{1l}, \end{equation}

with

(3.15)$$\begin{gather} \alpha_{l} = (2l+1)(l-1)\frac{\mu}{\rho R^2} , \end{gather}$$

(3.16)$$\begin{gather}\omega_{l}^2 = l(l-1)(l+2) \frac{\sigma}{\rho R^3} \end{gather}$$

and $\delta _{1l}$ is the Kronecker delta function. This model is valid in the weakly viscous limit, when $Oh\ll 1$. An extension of the free droplet model to arbitrary ${Oh}$ can be found in other prior works (Miller & Scriven Reference Miller and Scriven1968; Moláček & Bush Reference Moláček and Bush2012; Chandrasekhar Reference Chandrasekhar2013).

3.3. Pressure forcing during impact

There is still an additional unknown in the bath mode (3.10) and drop mode (3.14) equations: the applied pressure distribution $p_s(r,t)$. This is generally a function of the properties of the fluid medium, the impacting speed of the object, the shape of the impacting object and the motion of the gas that surrounds the fluid. However, in this work, we will assume that the viscosity of the ambient gas is small relative to the fluid bath such that the flow within the small air film is negligible, and the pressure acts solely to apply an upwards hydrodynamic force on the droplet. In non-dimensional terms, we can construct two additional restrictions for our model, $\rho _g / \rho \ll 1$ and ${Oh}_g = \mu _g / \sqrt {\rho \sigma R} \ll {Oh}$ following prior work (Moláček & Bush Reference Moláček and Bush2012). In the current experimental and DNS work, $\rho _g/ \rho \sim O(10^{-3})$ and ${Oh}_g \sim O(10^{-4})$, thus the influence of these two additional parameters is indeed negligible. We have verified this for our typical experimental parameters through DNS, and find that both a four-fold increase and decrease in the ambient density and viscosity from air properties at standard temperature and pressure (STP) produces negligible changes to the trajectory, instantaneous shape of the droplet, and free surface shape throughout the interaction of the droplet and bath (Appendix A).

The radial extent of the pressure distribution is generally unknown for impact problems, and constitutes an additional problem that we must solve. In the present work, we assume that this unknown pressure distribution takes the form

(3.17)\begin{equation} p_s(r,t) = \frac{F(t)}{{\rm \pi} R^2} H_r(r/r_c(t)), \end{equation}

where $F$ is the instantaneous magnitude of the contact force, evaluated at $r=0$ and $H_r$ is an assumed spatial profile of the pressure in the contact region. For this distribution, we can use a function that resembles the true shape of the pressure distribution during contact. The contact region, $A_c$, will be assumed to a simply connected disk, following Galeano-Rios et al. (Reference Galeano-Rios, Milewski and Vanden-Broeck2017) and Korobkin (Reference Korobkin1995). This allows us to write a single unknown $r_c(t)$ to fully describe the temporal evolution of this region of contact. Blanchette (Reference Blanchette2016, Reference Blanchette2017) used a fixed parabolic pressure shape function

(3.18)\begin{equation} H_r(r) = \begin{cases} C \left( 1 - \left(\dfrac{r}{R}\right)^2 \right), & r\leq R ,\\ 0, & r > R . \end{cases} \end{equation}

Here, $C$ is the constant magnitude of the pressure at $r=0$ and $R$ is the undeformed radius of the droplet. Blanchette (Reference Blanchette2016, Reference Blanchette2017) chose the value of the magnitude $C$ such that $\int _0^b p_s r \,{\rm d}r= {\rm \pi}R^2$. Thus, the pressure acting on the bath interface in the respective models had a constant pressure shape function $H_r$ for all times during contact. However, simulation results from Galeano-Rios et al. (Reference Galeano-Rios, Milewski and Vanden-Broeck2017, Reference Galeano-Rios, Cimpeanu, Bauman, MacEwen, Milewski and Harris2021) show that the shape of the pressure distribution at the surface of a fluid bath due to an impacting, non-wetting sphere is flatter and more similar to a top-hat function for most times, and that the spatial extent of the distribution changes continuously with time during impact. Additionally, for a droplet impacting on a solid surface, the pressure in the air film has been inferred by de Ruiter et al. (Reference de Ruiter, Lagraauw, Van Den Ende and Mugele2015). The air film thickness during a bounce was measured using interferometry and the pressure estimated using a lubrication model. The film pressure in both the impacting and rebounding regimes is approximately uniform, with deviations from uniformity only near the edge of the film. In related work, the impact pressure between a droplet and a wettable solid substrate has been studied extensively by Mandre, Mani & Brenner (Reference Mandre, Mani and Brenner2009) and Mani, Mandre & Brenner (Reference Mani, Mandre and Brenner2010), and their results indicate that the impact pressure increases sharply near the contact line, likely a consequence of the decreased air film thickness in that region. For the case of droplets bouncing on a deep pool, Tang et al. (Reference Tang, Saha, Law and Sun2019) measured the air film thickness and found the film thickness to be significantly more uniform in both impacting and rebounding stages for ${We}$ values similar to those explored in the present work, presumably as a result of the deformability of the substrate and impactor. Our predictions from DNS (presented and discussed in § 4) similarly suggest a more uniform air film thickness for the present problem, and a nearly uniform pressure profile during all stages of rebound.

We will use a simple polynomial that resembles a smoothed top hat in this work, with

(3.19)\begin{equation} H_r(r/r_c(t)) = \begin{cases} C \left( 1 - \left(\dfrac{r}{r_c(t)}\right)^6 \right), & r\leq r_c(t), \\ 0, & r > r_c(t) . \end{cases} \end{equation}

In order to remain consistent with our linearization, we do not allow $r_c$ to exceed $R$. Requiring that the integral of the pressure over the contact area is $F(t)$, we find

(3.20)\begin{equation} \int_0^{b} H_r(r/r_c) r \, {\rm d}r = \frac{R^2}{2} , \end{equation}

which sets the constant $C$ in the pressure shape function $H_r(r/r_c)$. Our bath model relies on the decomposition of the fluid motion into a linear superposition of infinitely many waves with wavenumbers $k_m$. Therefore, we apply a similar decomposition to this pressure function $p_s=\sum _{m=0}^{\infty } \,{\rm d}_m {\rm J}_0(k_m r)$. Since we are working in a cylindrical domain, we will choose the zeroth-order Bessel function of the first kind as our orthogonal basis function, and thus the $d_m$ are the Fourier–Bessel coefficients of the function $H_r$:

(3.21)\begin{equation} d_m = \frac{2}{(b {\rm J}_1(k_m))^2}\int_0^b H_r(r) r {\rm J}_0(k_m r) \,{\rm d}r, \end{equation}

with the domain extending from $r=[0,b]$. The reconstruction of the top-hat function in Fourier–Bessel space converges too slowly to be of practical use (Storey Reference Storey1968), also noted by Blanchette (Reference Blanchette2016), and as such we use a polynomial expression that resembles a smoothed top hat. Additionally, we tested higher-order polynomials (corresponding to a larger flat region), and found increasingly poor convergence behaviour, similar to that of the top hat (see Appendix B for a case study on the sensitivity of the results to the choice of shape function). The ultimate choice of shape function used here thus represents a practical compromise.

Substituting in the definition of the pressure (3.17) into (3.10), performing the Fourier–Bessel decomposition, we find

(3.22)\begin{align} \frac{{\rm d}^2 a_m}{{\rm d}t^2} + 4 \nu k_m^2 \frac{{\rm d}a_m}{{\rm d}t} + \left(\frac{\sigma k_m^2}{\rho} + g\right)k_m\tanh{(k_m h_0)}a_m ={-}\frac{F}{\rho {\rm \pi}R^2} d_{m} k_m\tanh{(k_m h_0)} , \end{align}

which govern the evolution of bath wave modes $m$. Similarly, substituting the pressure (3.17) into (3.14), we write

(3.23)\begin{equation} \frac{{\rm d}^2 \beta_l}{{\rm d}t^2} + 2 \alpha_{l,0} \frac{{\rm d} \beta_l}{{\rm d}t} + \omega_{l,0}^2 \beta_l ={-}\frac{F}{2 {\rm \pi}\rho R^3} c_{l} (2l+1) l + g \delta_{1l} . \end{equation}

The coefficients $c_{l}$ result from the mode decomposition of the projection of the pressure into spherical space,

(3.24)\begin{equation} c_{l} = \int_0^{\rm \pi} H_r(\theta) \sin(\theta) P_l(\cos{\theta}) \,{\rm d}\theta , \end{equation}

which naturally arise in the derivation of (3.14). Additionally, the definition of the pressure (3.17) reduces the governing equation of the $l=1$ ‘translational’ mode to

(3.25)\begin{equation} \frac{{\rm d}^2 \beta_1}{{\rm d}t^2} ={-}\frac{F}{m} + g , \end{equation}

which clearly governs the droplet centre of mass motion $\beta _1=-z_{cm}$. Evidence from the simulations of Galeano-Rios et al. (Reference Galeano-Rios, Milewski and Vanden-Broeck2017) indicates that impact trajectory is very sensitive to the instantaneous size of the contact area. Utilizing a constant pressed area for the pressure, as done in Blanchette (Reference Blanchette2016), does not produce results that compare well with experiment (Appendix B), particularly for cases of small ${We}$. Appendix B details how the choice of this pressure shape function modifies the predicted trajectory of the droplet for typical experimental parameters. The trajectory is largely insensitive to the choice of pressure shape function, but incorporating a time-dependent contact radius is essential for agreement. The method for determining both $F(t)$ and $r_c(t)$ are discussed in the next section.

3.4. Modelling contact

The contact force $F(t)$ is determined through the use of a ‘1-Point’ kinematic match (1PKM) condition. Essentially, we enforce contact only at a single point; the centre of our axisymmetric domain. Thus, the additional constraint can be written as

(3.26)\begin{equation} \eta(r=0,t) = z_{cm}(t) - \xi(\theta=0,t) = \sum_{m=0}^{\infty} a_m(t) = z_{cm}(t) -\left(R + \sum_{l=2}^{\infty}\beta_l(t)\right) . \end{equation}

This additional constraint allows us to determine the unknown contact force $F(t)$. Contact between the droplet and the bath ends when the magnitude of the contact force as predicted by the kinematic match becomes negative. We note that this 1PKM model is a significant simplification of the full kinematic match successfully used to study related impact problems (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 full kinematic match predicts the evolution of the contact area and the contact pressure distribution (without requiring an assumption for $H_r$) by imposing natural geometric and kinematic constraints, essentially considering additional equations to solve at each time step. The algorithm requires iteration at each time step, and the minimization of a tangency boundary condition is used to determine the correct contact area and pressure shape.

Lastly we turn to the unknown contact radius $r_c(t)$. By not restricting the deformations of the bath and droplet interface with the use of additional tangency and distributed kinematic match conditions, we find that the results of our simulation consistently produce interfacial shapes that cross each other. The amount of overlap between the two interfaces is generally small, for the typical experimental parameters in our study the maximum overlap is less than $0.05R$. However, we can use the predictions from both interfacial models to determine the exact location where the two interfaces cross and separate, and use this as the instantaneous radius of contact $r_c(t)$. Thus, at each time, contact between the bath and drop is ensured at both $r=0$ and $r=r_c(t)$. In order to enforce contact within the entirety of the contact region, a full kinematic match would be required – this circumvents the need for any assumptions on the pressure profile shape, but is substantially more computationally expensive. Our contact radius criterion is similar to that of the numerical model presented in prior work on droplet rebound from solid substrates (Moláček & Bush Reference Moláček and Bush2012). While this method is unphysical, it yields accurate predictions for the contact radius as compared with DNS as demonstrated in § 5.

3.5. Summary

Choosing a time scale of $t_{\sigma }=\sqrt {\rho R^3/\sigma }$ (with $\tau = t/t_{\sigma }$), a length scale of $R$, and a force scale of $2 {\rm \pi}\sigma R$ (with $f = F / 2 {\rm \pi}\sigma R )$, we recast the governing equations in non-dimensional form as

(3.27)$$\begin{gather} \eta(r,\tau) = \sum_{m=0}^{M} a_m(\tau) {\rm J}_0(k_m r) , \end{gather}$$

(3.28)$$\begin{gather}\xi(\theta,\tau) = 1 + \sum_{l=2}^{L}\beta_l(\tau) P_l(\cos{\theta}) , \end{gather}$$

(3.29)$$\begin{gather}\frac{{\rm d}^2 {a_m}}{{\rm d}\tau^2} + 4 {Oh} {k_m}^2 \frac{{\rm d} {a_m}}{{\rm d}\tau} + ({k_m}^2 + {Bo}){k_m}\tanh{({k_m} {h_0})}{a_m} ={-}2{f} d_{m,0} {k_m} \tanh{({k_m} {h_0})} , \end{gather}$$

(3.30)$$\begin{gather}\frac{{\rm d}^2 {\beta_l}}{{\rm d}\tau^2} + 2 {Oh} (2l+1)(l-1) \frac{{\rm d} {\beta_l}}{{\rm d}\tau} + l(l-1)(l+2) {\beta_l} ={-} {f} c_{l} (2l+1)l + {Bo}\delta_{1l}, \end{gather}$$

(3.31)$$\begin{gather}\frac{{\rm d}^2 z_{cm}}{{\rm d} \tau^2} = \frac{3}{2} f - {Bo} , \end{gather}$$

(3.32)$$\begin{gather}\eta(0,\tau) = {z}_{cm}(\tau) - {\xi}(\theta=0,\tau) . \end{gather}$$

Equations (3.29) and (3.30) describe the evolution of the bath and droplet oscillation modes, respectively. Equation (3.31) governs the vertical motion of the droplet's centre of mass. Equation (3.32) couples these equations all together, with $r=0$ as the single point of ‘contact’ enforced between the droplet and the bath and allows for determination of the unknown $f(\tau )$. These equations are solved using standard ODE numerical integration techniques. The shape of the bath and droplet can be reconstructed at any time $t$ via the sums in (3.27) and (3.28), respectively.

The complete model is valid when $Re = \sqrt {We}/Oh \gg 1$ and $Oh\ll 1$. Also, since the model is linearized about the undeformed state, we anticipate it to hold when deformations remain small, further suggesting $Bo\ll 1$ and $We\ll 1$. However, we later demonstrate through direct comparison with experiment and DNS that the model remains predictive even for moderate ${We}$.

3.6. Numerical methods

We solve these equations using a backward Euler method, ensuring a minimum of 100 time steps within the inertio-capillary time $t_{\sigma }$. An implicit method was chosen, following Galeano-Rios et al. (Reference Galeano-Rios, Cimpeanu, Bauman, MacEwen, Milewski and Harris2021), as the instantaneous size of the pressure distribution acting on both the droplet and bath at the next time step is unknown. Treating the pressure explicitly on either the droplet or the bath can lead to non-physical behaviour in the system. We used $M=150$ modes for the bath interface and $L=55$ modes for the droplet interface. These values were determined by running simulations of a ${We} =0.7$, ${Bo} =0.017$, ${Oh} =0.006$ impact and assessing convergence as described in what follows. First, we kept the number of droplet modes fixed at $L=15$ and increased the number of bath modes from $30$ to $500$ in increments of $25$. Then, the simulation was run again, fixing the number of bath modes at $75$ and increasing droplet modes from $15$ to $200$. Sufficient convergence was determined if the maximum absolute value of the difference in centre of mass trajectories during contact changed by less than $1\,\%$. Finally, both the droplet and bath number of modes were increased simultaneously, and convergence was still observed. These values are similar to comparable to those found in Blanchette (Reference Blanchette2016) ($M=200$, using sine functions as the basis functions in a square bath), and Moláček & Bush (Reference Moláček and Bush2012) ($L=150$, but found good accuracy in comparison with experimental data at $L=20$). Additionally, once mode convergence was determined, we decreased the time step of the simulation in increments until time step convergence was similarly reached using the same criterion. Unexpectedly, in a fully converged simulation, there is a time step threshold below which the algorithm results in unstable oscillations in the magnitude of $F(t)$. This time step threshold is typically at least three orders of magnitude smaller than $t_{\sigma }$. This apparent instability deserves awareness and future attention, but does not affect the results presented in the present work. The bath size was set to $b=25 R$ which was determined to be sufficiently large such that reflected waves did not influence the droplet during impact. In order to find the instantaneous contact radius, we take two line segments from the reconstruction of the droplet and bath interface shapes. From these we can write a linear system of equations for four unknowns: the $(r,z)$-pair of the intersection location, the normalized distance from the starting point of the first line segment to the intersection, and the normalized distance from the starting point of the second line segment to the intersection (Schwarz Reference Schwarz2022). We then loop over every line segment to find every intersection. We take the largest of this set as $r_c(t)$. In the reconstructions of the interfaces at each time step, at least 5000 points are used in both $\theta$ and $r$ to ensure that error is minimized. All code associated with the implementation of the model is available at https://github.com/harrislab-brown/BouncingDroplets.