Lubrication-mediated rebounds off fluid baths

Abstract

We present herein the derivation of a lubrication-mediated (LM) quasi-potential model for droplet rebounds off deep liquid baths, assuming the presence of a persistent dynamic air layer which acts as a lubricating pressure transfer. We then present numerical simulations of the LM model for axisymmetric rebounds of solid spheres and compare quantitatively to current results in the literature, including experimental data in the low-speed impact regime. In this regime the LM model has the advantage of being far more computationally tractable than direct numerical simulation (DNS) and is also able to provide detailed behaviour within the micro-metric thin lubrication region. The LM system has an interesting mathematical structure, with the lubrication layer providing a free-boundary elliptic problem mediating the drop and bath free-boundary evolutionary equations.

1. Introduction

The role of the entrained air layer present throughout droplet rebounds has been well studied. Since the work of Worthington (1881), where it was first questioned, it has been well established that rebounds are only able to occur in the presence of a lubricating layer of air preventing the two liquids contacting. The assumption that within the near-contact region, the air can behave like a lubrication layer (Hicks & Purvis 2010) has directed a rich field of literature of impacts on other substrates such as porous media (Hicks & Purvis 2017). In the case of the impactor being solid, the same assumption can be made in the pre-impact time; that the layer of air helps delay contact (Moore 2021). Recent studies of drops that bounce exist in several configurations, for example off solid hydrophobic surfaces (Kolinski, Mahadevan & Rubinstein 2014), and off liquid films (Tang et al. 2019).

There are fewer studies of drops rebounding off deep baths, particularly in the regime where the underlying fluid is inviscid (or nearly so) and capillary waves are long lived after the droplet impact. These models gained interest in part because of the work of Couder on indefinitely bouncing drops on vibrating baths (Couder et al. 2005) and its fascinating quantum-mechanical analogies (Bush 2015, and references therein). The studies of Galeano-Rios, Milewski & Vanden-Broeck (2017, 2019) and Alventosa, Cimpeanu & Harris (2023) imposed a so-called kinematic match (KM) between the solid sphere or drop and the bath to avoid modelling the air layer and obtained accurate results when compared with Faraday pilot-wave studies of bouncing and walking droplets. In Galeano-Rios et al. (2021) the KM model was validated against experiments and direct numerical simulation (DNS) of the Navier–Stokes equations of single rebounds of solid hydrophobic spheres for low-Weber-number impacts.

Current reduced models for solid or droplet rebounds on deep baths do not capture the dynamics of the air layer which is present throughout the drop-bath dynamics. This paper builds on the work on the impact of two-dimensional drops in Phillips, Cimpeanu & Milewski (2024), where the lubrication air layer is considered as an additional fluid region coupling the drop deformation and the wave equations for the bath. We extend that model to three-dimensions in generality (§ 2), and compare its results to the KM, DNS and experiments for solid hydrophobic sphere impacts on water (§ 3).

2. Model

Consider an almost spherical droplet falling vertically downwards through the air, towards a deep liquid bath initially at rest. Let the subscripts $\alpha \in \{a,b,d\}$ denote the air, bath and droplet, respectively, and introduce notation $\varOmega _\alpha$ for each fluid domain, $\boldsymbol {u}_\alpha$ for the fluid velocity, with the density, dynamic viscosity, kinematic viscosity and pressure of the fluids given as $\rho _\alpha,\mu _\alpha,\nu _\alpha,p_\alpha$, respectively. Finally we use $\beta \in \{b,d\}$ to denote the bath–air and droplet–air boundaries $\partial \varOmega _\beta$, defined by a level set of $F_\beta$ and the surface tension coefficient $\sigma _\beta$ at the air–liquid interfaces. In each fluid region, we begin from the Navier–Stokes equations, and at the interfaces the usual kinematic condition, stress continuity and velocity continuity boundary conditions $\boldsymbol {u}_a = \boldsymbol {u}_\beta$ are taken. Hence,

(2.1)$$\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}_\alpha = 0, \quad \boldsymbol{x} \in \varOmega_\alpha, \end{gather}$$

(2.2)$$\begin{gather}\rho_\alpha ( \partial_t \boldsymbol{u}_\alpha + (\boldsymbol{u}_\alpha \boldsymbol{\cdot} \boldsymbol{\nabla}) \boldsymbol{u}_\alpha) ={-} \boldsymbol{\nabla} p_\alpha + \mu_\alpha \nabla^2 \boldsymbol{u}_\alpha - \rho_\alpha g \boldsymbol{e}_z, \quad \boldsymbol{x} \in \varOmega_\alpha, \end{gather}$$

(2.3)$$\begin{gather}D_t F_\beta = 0, \quad \boldsymbol{x} \in\partial\varOmega_\beta, \end{gather}$$

(2.4)$$\begin{gather}p_\beta \boldsymbol{n}_\beta + \mu_\beta \tau_\beta \boldsymbol{n}_\beta = p_a\boldsymbol{n}_\beta + \mu_a \tau_a \boldsymbol{n}_\beta - \sigma_\beta \kappa_\beta, \quad \boldsymbol{x} \in\partial\varOmega_\beta, \end{gather}$$

where $\tau _\alpha = (\boldsymbol {\nabla } \boldsymbol {u}_\alpha + \boldsymbol {\nabla } \boldsymbol {u}_\alpha ^T)$ is the stress tensor for each fluid, $\kappa _\beta$ denotes curvature of the free surface with outward unit normal vector $\boldsymbol {n}_\beta$, gravity $g$ acting downwards in the $-\boldsymbol {e}_z$ direction, and $D_t$ is the material derivative. These equations are supplemented by imposing decay of velocities in the far field or appropriate boundary conditions in the case of finite domains. We proceed with several modelling approximations which will reduce the system to a tractable state.

2.1. Linear quasi-potential model for the liquid bath

We assume viscous effects are small except in the lubrication layer, and follow methodology from Dias, Dyachenko & Zakharov (2008) which we briefly outline below. Take a Helmholtz decomposition of the liquid as a small perturbation from potential flow $\boldsymbol {u}_b = \boldsymbol {\nabla } \phi _b + \boldsymbol {\nabla } \times \boldsymbol {\psi _b}$, and use tangential stress conditions and a boundary-layer argument to eliminate the $\boldsymbol {\psi }_b$ terms, expressing the rotational part of the flow in terms of $\phi _b$ and $\eta _b$. A more detailed explanation of this argument can be found in Phillips et al. (2024). The free surface of the bath is given by $z=\eta _b(\boldsymbol {x},t)$ where $z=0$ is its undisturbed free surface. The governing linear system then becomes

(2.5)$$\begin{gather} 0= \Delta_H \phi_b + \partial_z^2 \phi_b, \quad z\le0, \end{gather}$$

(2.6)$$\begin{gather}\partial_t \phi_b ={-}\frac{1}{\rho_b}p_a - g\eta_b + \frac{\sigma_b}{\rho_b}\Delta_H \eta_b + 2\nu_b \Delta_H\phi_b, \quad z=0, \end{gather}$$

(2.7)$$\begin{gather}\partial_t \eta_b = \partial_z \phi_b + 2\nu_b \Delta_H \eta_b, \quad z=0, \end{gather}$$

where $\Delta _H$ is used to denote the horizontal Laplacian, the curvature is small and of the form $\kappa _b = -\boldsymbol {\nabla }_H \eta _b$, and velocity continuity $\boldsymbol {u}_b=\boldsymbol {u}_a$ is imposed at the interface $\eta _b$ between the bath and air layer. The terms proportional to $\nu _b$ arise from $\boldsymbol {\psi _b}$ and are the corrections to the pressure and vertical velocity, respectively.

2.2. Air-layer lubrication approximation

Away from the droplet–bath interaction, the air is assumed to be atmospheric pressure and have negligible effect on the system, as we must introduce the domain in which the lubrication effects may be important. The lower part of the droplet is described by $z=\eta _d(\boldsymbol {x}_H,t) = \boldsymbol {X}(t)-S(\boldsymbol {x}_H,t)$, composed of the shape of the lower part of the droplet $S$, and the vertical height of its centre of mass $\boldsymbol {X}$. The domain of $S$ in $\boldsymbol {x}_H$ is chosen such that $S$ is single valued and with bounded gradients. The dynamics of the air layer are described through a lubrication approximation as in Phillips et al. (2024). Consider a lubrication region $\varOmega _a = \varOmega _L\times [\eta _b(r,\theta,t),\eta _d(r,\theta,t)]$ where its footprint is described in cylindrical coordinates as

(2.8)\begin{equation} \varOmega_{L} = \{(r,\theta) : r\in [0, r^*(\theta,t)], \theta \in [0,2{\rm \pi}) \}, \end{equation}

where $r^*$ is the edge of the lubrication region. Denoting $h=\eta _d-\eta _b$, the edge of the lubrication layer is given by a criterion $h(r^*,\theta,t) = \varepsilon$. It is assumed that $\varOmega _L$ is a subset of the domain of $S$. Assuming that the layer between the droplet and the bath is thin, the leading-order balance of the Navier–Stokes equations results in the lubrication balance where $\partial _z p_a = 0$ with

(2.9)$$\begin{gather} \boldsymbol{\nabla}_H \boldsymbol{\cdot} \boldsymbol{u}_a^H + \partial_z w_a = 0, \quad \boldsymbol{x}\in \varOmega_{a}, \end{gather}$$

(2.10)$$\begin{gather}\mu_a \partial_z^2 \boldsymbol{u}_a^H = \boldsymbol{\nabla}_H p_a, \quad \boldsymbol{x}\in \varOmega_{a}, \end{gather}$$

and velocity continuity holds at both interfaces: $\boldsymbol {u}_a = \boldsymbol {u}_d$ at $z =\eta _d$ and $\boldsymbol {u}_a = \boldsymbol {\nabla }_H \phi$ at $z =\eta _b$. Here $\boldsymbol {\nabla }_H$ is the horizontal gradient operator, and $\boldsymbol {u}_a^H = (u_r,u_\theta )$ is the horizontal velocity of the air layer, such that $\boldsymbol {u}_a = (\boldsymbol {u}_a^H,w_a)$. As expected, the pressure $p_a = P(r,\theta,t)$ is independent of $z$ throughout the air layer. Further, we can solve the equation for $\boldsymbol {u}_a^H$ in terms of $P$:

(2.11)\begin{equation} \boldsymbol{u}_a^H = \frac{\boldsymbol{\nabla}_H P}{2\mu_a}(z-\eta_b)(z-\eta_d) + \frac{\boldsymbol{\nabla}_H \phi_b|_{\eta_b}}{\eta_b- \eta_d}(z-\eta_d) + \frac{\boldsymbol{u}^H_d|_{\eta_d}}{\eta_d - \eta_b}(z-\eta_b). \end{equation}

To obtain the thin-film equation, we integrate vertically the conservation of mass (incompressibility) equation resulting in

(2.12) \begin{equation} \boldsymbol{\nabla}_H \boldsymbol{\cdot} \int_{\eta_b} ^{\eta_d} \boldsymbol{u}_a^H \,\text{d}z +\Big(\boldsymbol{u}_a^H \Big|_{\eta_b} \boldsymbol{\cdot} \boldsymbol{\nabla}_H\eta_b -\boldsymbol{u}_a^H \Big|_{\eta_d} \boldsymbol{\cdot} \boldsymbol{\nabla}_H\eta_d \Big) + (w_a|_{\eta_d} - w_a|_{\eta_b})= 0. \end{equation}

All terms after the integral can be expressed exactly as $\partial _t\eta _d- \partial _t \eta _b = \partial _t h$ from the kinematic conditions at both interfaces. Finally, evaluating the integral in the first term using (2.11) to obtain the form of the flux

(2.13)\begin{equation} \boldsymbol{Q}= \int_{\eta_b}^{\eta_d} \boldsymbol{u}_a^H \,\text{d}z ={-}\frac{1}{12\mu_a}h^3 \boldsymbol{\nabla}_H P + \frac{1}{2}h (\boldsymbol{\nabla}_H\phi_b + \boldsymbol{u}_d^H) \end{equation}

results in the thin-film equation

(2.14)\begin{equation} \partial_t h + \boldsymbol{\nabla}_H \boldsymbol{\cdot} \boldsymbol{Q} = 0, \quad \boldsymbol{x_H} \in \varOmega_{L}. \end{equation}

As we shall see below, this should be considered a nonlinear free-boundary elliptic equation for $P$ where $-h_t$ is considered as a forcing term, with the boundary condition $P=0$ at the evolving boundary curve $\partial \varOmega _{L}$.

2.3. Droplet model

2.3.1. Droplet deformation

Linear damped droplet oscillations are well understood (Lamb 1924) and the results have been used in a variety of contexts, e.g. Aalilija, Gandin & Hachem (2020). In this section we write equivalent equations for these oscillations using the same quasi-potential theory used above for the bath. Consider a capillary-scale droplet of liquid which when unperturbed is a sphere of radius $R_0$. Describing the problem using spherical polar coordinates with azimuthal angle $\theta$ and polar angle $\varphi$, we can take a Helmholtz decomposition for the liquid velocity $\boldsymbol {u}_d = \boldsymbol {\nabla } \phi _d + \boldsymbol {\nabla } \times \boldsymbol {\psi }_d$. Using the same argument argument as above (Phillips et al. 2024) we obtain similar governing equations to those of the bath (2.5)–(2.7):

(2.15)$$\begin{gather} \Delta \phi_d = 0,\quad r< R_0, \end{gather}$$

(2.16)$$\begin{gather}\partial_t \phi_d ={-}\frac{1}{\rho_d}p_a - \frac{\sigma_d}{\rho_d} \kappa_d - 2\nu_d \partial_r^2 \phi_d,\quad r=R_0, \end{gather}$$

(2.17)$$\begin{gather}\partial_t \zeta_d = \partial_r \phi_d + 2\nu_d\left(\frac{1}{r^2}\Delta_s \zeta_d - \frac{1}{r^3}\int \Delta_s \phi_d \,\text{d}t\right), r=R_0, \end{gather}$$

where surface of the droplet is given by $r=R_0+\zeta _d(\theta,\varphi )$, and $\kappa _d$ is the mean curvature

(2.18a,b)\begin{equation} \kappa_d ={-} \left( \frac{2}{R_0^2} \zeta_d + \frac{1}{R_0^2} \Delta_s \zeta_d \right), \quad \Delta = \frac{1}{r^2}\partial_r (r^2 \partial_r) + \frac{1}{r^2}\Delta_s, \end{equation}

with $\Delta _s$ the Laplace–Beltrami operator (i.e. the surface Laplacian). The equations are valid in the small-oscillation limit and have been truncated to contain terms that are linear in the Ohnesorge number (see Appendix A for a discussion of the non-dimensional equations). One may proceed with an eigenfunction decomposition of $\zeta _d$ and $\phi _d$ as

(2.19a,b)\begin{equation} \phi_d = \sum_{l=2}^\infty \sum_{m={-}l}^l a_{lm}(t) r^l Y_l^m, \quad \zeta_d = \sum_{l=2}^\infty \sum_{m={-}l}^l c_{lm}(t) Y_l^m, \end{equation}

where the $a_{lm}(t), c_{lm}(t)$ are coefficients of the $Y_l^m (\theta,\varphi ) = {\rm e}^{{\rm i}m\varphi }P_l^m(\cos \theta )$ spherical harmonics that satisfy

(2.20)\begin{equation} {-}\Delta_s Y_l^m = (l)(l+1)Y_l^m. \end{equation}

The $P_l^m$ are the associated Legendre polynomials of degree $l$ and order $m$. The terms $l=0,1$ have been omitted as $l=0$ corresponds to a dilation of the sphere and $l=1$ corresponds to translational and rotational symmetries which need to be considered separately as they will be affected by the external forces due to the air layer. Substitution into the boundary conditions (2.16)–(2.17) leads to the differential equation (where leading order contributions have been determined with respect to the Ohnesorge number, Oh, and terms of order $Oh^2$ have been omitted)

(2.21)\begin{equation} \ddot{c}_{l,m} + 2\lambda_{l} \dot{c}_{l,m} + \omega_{l}^2c_{l,m} ={-}\frac{l\hat{p}_{l,m}}{\rho R_0}, \end{equation}

where $\hat {p}_{l,m}$ is the $Y_{l,m}$ coefficient of the lubrication air pressure $p_a$. The coefficients of the velocity potential are given to leading order by $a_{l,m}={\dot {c}_{l,m}}/{(lR_0^{l-1}})$. The damping coefficients and modal frequencies are given by

(2.22a,b)\begin{equation} \lambda_l = (2l+1)(l-1)\frac{\mu}{\rho R_0^2}, \quad \omega_l^2 = l(l-1)(l+2) \frac{\sigma}{\rho R_0^3}, \end{equation}

as seen in the literature (Lamb 1924; Aalilija et al. 2020; Alventosa et al. 2023). A similar equation may be obtained for the coefficients $a_{lm}(t)$; however, these are unnecessary for the lubrication-mediated (LM) model and have been omitted for brevity.

2.3.2. Equations for the centre of mass

The location of the droplet is determined through updating the position of the centre of mass, given by $\boldsymbol {X}(t)$. We shall disregard drag effects from the air, instead only considering contributions from the lubrication layer, then the equation of motion is

(2.23)\begin{equation} m\ddot{\boldsymbol{X}} = F ={-}\int_{\varOmega_L} p_a \, \boldsymbol{n_d} \,\text{d}a - mg \, \hat{\boldsymbol{z}}, \end{equation}

where $\boldsymbol {n_d}$ is the outward normal of the droplet, $m = 4\rho {\rm \pi}R_0^3/3$ is its mass, $g$ is the gravitational constant and $\hat {\boldsymbol {z}}$ is the unit vector in the vertical direction. The boundary of the droplet is now given by $\boldsymbol {X} + (R_0+\zeta _d(\theta,\varphi )) \boldsymbol {\hat {r}}$, where $\boldsymbol {\hat {r}}$ is a radial unit vector relative to the position $\boldsymbol {X}$.

We note that we have disregarded global rotational motion of the droplet (forced by both asymmetric pressure and shear stresses in the lubrication layer) which would significantly add to the modelling complexity, and which will not be relevant for the vertical (axisymmetric) impacts we shall consider next. We also expect these effects to be small for quasi-normal impacts (e.g. Galeano-Rios et al. 2019).

2.4. The lubrication-mediated droplet-bath system

We can now propose a closed system of equations for the evolution of the thin film, liquid bath and droplet. Within the bath and the droplet one may reduce the system further using the framework of Dirichlet to Neumann (DtN) maps. Consider the ‘trace’ of the potentials on the boundary of the domains in which Laplace's equation are solved:

(2.24a,b)\begin{equation} \varPhi_b(\boldsymbol{x_H},t) = \phi_b|_{z=0}, \quad \varPhi_d(\theta,\varphi,t) = \phi_d|_{r=R_0}. \end{equation}

The DtN map, denoted $D$, is the linear map that provides the normal derivative of the potential given its value on the boundary. In particular we denote

(2.25a,b)\begin{equation} (\partial_z \phi_b)|_{z=0} = D_b \varPhi_b, \quad (\partial_r \phi_d)|_{r=R_0} = D_d\varPhi_d. \end{equation}

We shall see that these maps are easily expressed in terms of eigenfunction expansions. Making use of DtN maps, we now have the evolution system

(2.26)$$\begin{gather} \partial_t \varPhi_b ={-}\frac{1}{\rho_b} p_a -g\eta_b+ \frac{\sigma_b}{\rho_b}\Delta_H \eta_b + 2\nu_b \Delta_H\varPhi_b , \end{gather}$$

(2.27)$$\begin{gather}\partial_t \eta_b = D_b \varPhi_b + 2\nu_b \Delta_H \eta_b \doteq F_b, \end{gather}$$

(2.28)$$\begin{gather}\partial_t \varPhi_d ={-}\frac{1}{\rho_d} p_a - \frac{\sigma_d}{\rho_d} \kappa_d - 2\nu_d \partial_r^2 \varPhi_d, \end{gather}$$

(2.29)$$\begin{gather}\partial_t \zeta_d = D_d \varPhi_d + 2\nu_d\left(\frac{1}{r^2}\Delta_s \zeta_d - \frac{1}{r^3} \Delta_s D_d^{{-}1} \zeta_d \right) \doteq F_d , \end{gather}$$

(2.30)$$\begin{gather}m\ddot{\boldsymbol{X}} = F ={-}\int_{\varOmega_L} p_a \, \boldsymbol{n_d} \,\text{d}a - mg \, \hat{\boldsymbol{z}}, \end{gather}$$

(2.31)$$\begin{gather}\boldsymbol{\nabla}_H \boldsymbol{\cdot} \boldsymbol{Q} = W+F_d-F_b, \quad \boldsymbol{x_H} \in \varOmega_{L}. \end{gather}$$

In the surface potential equations, the pressure $p_a$ is non-zero only in $\varOmega _L$ and $W={{\rm d}\kern0.7pt\boldsymbol {X}}/{{\rm d}t}$ is the vertical velocity of the droplet's centre of mass. In the kinematic condition for the droplet, we used the leading-order balance $\partial _t \zeta _d = D_d \varPhi _d$ to express the right-hand side locally in time. While $D_d$ was defined using the velocity potential, its inverse acts on $\zeta _d$ here. This is admissible since it is an invertible linear operator on functions with mean zero, and we may restrict the inverse map to have mean zero also. Equation (2.31), where $\boldsymbol {Q}$ is given by (2.13), is an inhomogeneous elliptic equation for the pressure with the boundary condition that $p_a=0$ on the free boundary $\partial \varOmega _L$.

3. Axisymmetric solid impacts and numerical results

We now consider the simplest case in which this framework can be used and compared with prior simulations and experiments: the rebound of solid (hydrophobic) spheres (Galeano-Rios et al. 2021). This considerably simplifies the system, reducing the dimension and eliminating droplet oscillations. Considering the axisymmetric domain for the bath of finite radius $r=L$, and imposing Neumann boundary conditions there for both $\eta _b$ and $\varPhi _b$, we expand these in terms of Bessel functions of the first kind,

(3.1a,b)\begin{equation} \eta_b(r,t) = \sum_{j=1}^\infty a_j(t) {\rm J}_0(k_j r), \quad \varPhi_b(r,t) = \sum_{j=1}^\infty b_j(t) {\rm J}_0(k_j r), \end{equation}

where the $a_j(t), b_j(t)$ are coefficients of the Bessel functions, and the $k_j$ satisfy ${\rm J}_0^\prime (k_j L) = 0$. We note that in this basis $\phi _b = \sum _{j=1}^\infty b_j {\rm J}_0(k_j r) {\rm e}^{k_j z}$ to satisfy Laplace's equation and the DtN map is expressed simply as

(3.2)\begin{equation} D_b \varPhi = \partial_z\phi_b|_{z=0} = \sum_{j=1}^\infty k_j b_j {\rm J}_0(k_j r). \end{equation}

Hence the DtN operator in this basis corresponds to multiplication by $k_j$. Similarly, the horizontal Laplacian $\Delta _H$ results in a Bessel multiplier $-k_j^2$. The system can now be written for the Bessel coefficients as

(3.3)$$\begin{gather} \dot{a}_j = k_jb_j - 2\nu_l k_j^2 a_j, \end{gather}$$

(3.4)$$\begin{gather}\dot{b}_j ={-}\frac{1}{\rho_b}\hat{p}_j - \frac{\sigma_b}{\rho_b}k_j^2 a_j - 2\nu_l k^2 b_j , \end{gather}$$

(3.5)$$\begin{gather}m\ddot{Z} = {2{\rm \pi}} \int_0^{r^*} p_a\; r\text{d}r - mg, \end{gather}$$

with the elliptic equation for pressure

(3.6)\begin{equation} {-}\frac{1}{r} \partial_r \left(\frac{r}{12\mu_a} h^3 \partial_r p_a + \frac{r}{2}h \partial_r \varPhi_b \right) = \partial_t h = W- (D_b \varPhi_b + 2\nu_l \Delta_H \eta_b), \end{equation}

where $W=\dot {Z}$, $h = Z-\sqrt {R_0^2-r^2} - \eta _b$ for $r\le R_0$, $r^*$ is defined by $h(r^*,t)=\varepsilon$ for $\varepsilon \ll 1$, and where $\hat {p}_j$ are the Bessel coefficients of the pressure provided by

(3.7)\begin{equation} \hat{p}_j = \frac{2}{(LJ_0(k_jL))^2}\int_0^L p_a(r,t) {\rm J}_0(k_j r) r \,\text{d}r. \end{equation}

Similar expressions apply for obtaining $a_j$ from $\eta _b$ and $b_j$ from $\varPhi _b$. The pressure can be integrated from (3.6) using $p_a(r^*,t)=0$ and the symmetry condition $\partial _r p_a(r^*,t)=0$:

(3.8) \begin{equation} p_a ={-} \int_{r^*}^r \Big[ \frac{ 12 \mu _a}{r'h(r',t)^3} \Big(\int _0 ^{r^\prime} r^{\prime \prime} \partial_t h(r'',t) \,\text{d}r^{\prime \prime} \Big) + \frac{ 6 \mu _a}{h(r',t)^2} \partial_r \varPhi_b(r',t) \Big] \,\text{d}r^{\prime} . \end{equation}

3.1. Numerical implementation

We shall compute a solution approximating the system (3.3)–(3.8) by truncating the Bessel expansions at a large $N$ and integrating the resulting differential equations (3.3)–(3.5) with a fourth-order Runge–Kutta scheme. The equation for the pressure (3.8) is calculated using a trapezium rule. Since Bessel expansions are ill-conditioned on regular grids, we use an oversampled non-uniform grid in $r$ with $M$ points distributed on Chebyshev collocation points on $[0,L]$ using $\theta _j = (j-1){\rm \pi} /(M-1)$ with $j = 1\ldots M$ and $r_j = {L(1+\cos (\theta _j))}/{2}$. Hence we compute two matrices: an $M\times N$ matrix which evaluates an $N$-term Bessel series at $r_j$ and an $N\times M$ matrix corresponding to calculating the projection formula (3.7) for the Bessel coefficients of a function given at $r_j$.

3.2. Model results

We first present results of the rebound of two representative spheres of different radii and densities rebounding off a deep-water bath against a sweep of initial downward velocities. The choice of parameters is displayed in table 1 and were chosen to correspond to the data for low-velocity impacts in Galeano-Rios et al. (2021). This will permit the results of the LM model derived within this paper to be compared with the KM model (Galeano-Rios et al. 2017) and DNS of the Navier–Stokes equations. In that study, experiments were also performed using hydrophobically coated solid spheres.

Table 1. Parameters used in simulations of a water bath and various solid spheres.

For such studies, the coefficient of restitution is usually taken as the negative ratio of the speeds between the time of impact $t_{imp}$, taken to be when the south pole of the sphere first crosses $z=0$ (correspondingly the centre of mass crosses $z=R_0$), and the time of lift-off $t_{lift}$, which is when the sphere leaves the same height on an upwards trajectory. The time spent below $z=0$ is the contact time $t_c = t_{lift}-t_{imp}$. In Galeano-Rios et al. (2017) for low impact velocities, the square coefficient of restitution $\alpha ^2$, which is the ratio of mechanical energies before and after impact (using the $z=0$ crossing as reference height), is used, as it can capture very small rebounds where $\alpha ^2<0$ where the sphere detaches but it's centre of mass does not reach the $z=R_0$ after rebound. Thus

(3.9)\begin{equation} \alpha^2 = \frac{E_{out}}{E_{in}} = \frac{W_{imp}^2}{W_{detach}^2-2g(Z-R_0)}, \end{equation}

where $W_{imp}$ and $W_{detach}$ are the vertical velocities at the moment of impact $t_{imp}$ and separation $t_{detach}$, when $\text {min}(h)>\varepsilon$ is first achieved after impact. In figure 1, $\alpha ^2$, the maximum penetration depth $\delta = -(Z_{min}-R_0)$, and the pressing time $t_p=t_{detach}-t_{imp}$ are displayed for the LM and KM models and DNS. In Galeano-Rios et al. (2017) $t_p$ is favoured over $t_c$ to capture the rebounds which do not return to $z=0$.

Figure 1. Comparison of $\alpha ^2$, $\delta$ and $t_p$ for two parameter regimes, with LM ($\circ$), compared with DNS ($\Diamond$) and KM ($*$) from Galeano-Rios et al. (2021).

Figure 1 demonstrates that the agreement is good between the LM model against both the results from the KM and DNS, with general trends being upheld. We note here that the latter two methods also involve certain approximations. The KM model assumes no lubrication layer (the method matches the position and velocities of the free surface to those of the sphere) and uses a tangential contact between the sphere and the bath throughout the impact. In addition, for the KM model we are comparing with, the evolution equations for the bath are linear and with the same quasi-potential approximation that we use here. The DNS results in Galeano-Rios et al. (2021) were performed on the open source Gerris where the solid drop is modelled by using artificially high viscosity and surface tension coefficients. Gerris uses a volume of fluid formulation and interfaces are approximated using continuous density changes (see Popinet 2009).

We conclude by providing numerical results for a configuration which can be compared with the KM model, DNS and also experimental values, as shown in table 2. The experimental configuration used a rectangular container while our simulation is cylindrical. In table 2 we varied the radius of the bath over comparable values to their length and width dimensions, resulting in a range of values for our simulations.

Table 2. Results for a falling sphere of radius $R_0= 0.83\times 10^{-4}$ m, density $\rho = 1200$ kg m$^{-3}$ and initial velocity $W_0 = -0.35\,{\rm m}\,{\rm s}^{-1}$ onto a deep-water bath. The range in the LM data arises from varying bath radius $L$.

The results for the LM model restitution coefficient $\alpha$ is at the lower end of the experimental range, the depth of impact is somewhat below the experimental range and the impact time is well within the experimental range. In addition to modelling assumptions, one potential source of error is that the physics of the hydrophobic coating used in experiments is not considered in any of the models.

Turning to impact details for this case, figures 2(a) and 3(b), show three distinct phases of rebound. First, a rapid expansion of the air layer, driven by the inertia of the impacting drop; second, a slow contraction of the layer; third, a rapid suction event due to the retraction of the free surface. The expansion of the region lasts $O(1\,\text {ms})$, the contraction lasts $O(10\,\text {ms})$, and the suction event lasts $O(1\,\text {ms})$ and precedes separation. A highlight of the LM model is the ability to capture and quantify the behaviour within the air layer during the rebound. In particular, we can study the evolution of the pressure profile and corresponding height layer within the system, quantities which are challenging to extract even in the DNS setting. In what follows, we have non-dimensionalised all lengths with the radius of the sphere (e.g. $r^*=r/R_0$, and area $A^* = A/R_0^2$), the pressure with the surface tension pressure jump across a spherical interface $P^* = P/(2\sigma R_0^{-1})$, and the force with $F^* = F/(2\sigma R_0)$. Throughout the impact, the height of the air layer retains thickness $O(1\,\mathrm {\mu }{\rm m})$, as expected from the literature (Tang et al. 2019). However, there is some variation as shown in figure 2(a), with a thickness of 0.5–2 $\mathrm {\mu }$m in the bulk of the layer and a narrowing on its boundary to 0.15–0.5 $\mathrm {\mu }$m. In effect, for the bulk of the impact, the air layer forms a quasi-static cushion constricted on its circular boundary. This narrowing effect is enhanced when the lubrication region boundary ‘turns around’ and immediately prior to detachment, where the lubrication layer has shrunk to a small disc at the base of the sphere. Turning to the profile of the pressure within the air layer, shown in figure 3(b), we find a range of 0.9–1.5 $2\sigma R_0^{-1}$ over the vast majority of the region. The maximum pressure is felt at impact and is 7.6 $2\sigma R_0^{-1}$ (indicating an inertially driven pressure) and the minimal pressure being $-$3.8 $2\sigma R_0^{-1}$ at liftoff. A pressure spike occurs where the layer thins at its edge, aligning with results found in Hendrix et al. (2016). Figure 3(a) captures the pressure profiles during expansion and contraction. Small oscillations observed in the pressure profile at two intermediate times, we believe, are numerical artefacts, as they are lessened by increased mesh resolution. Lesser refined meshes had larger oscillations, however, with minimal effect on the overall dynamics (restitution coefficient, impact times, generated waves, etc.) of the system.

Figure 2. (a) Thickness of the air layer, $h(r^*,t)$. (b) Free surface of the bath, $\eta ^*_b(r^*,t)$, shows the sinking depth of the sphere relative to the resultant waves on the free surface.

Figure 3. (a) Lubrication-layer pressure profiles $P^*(r^*,t)$ at selected times. Dark and light curves correspond to expanding and contracting layers, respectively. (b) Colourmap showing near uniform pressure in the interior and its rapid expansion and slow contraction.

Figure 4 shows a temporal slice of pressure, under the south pole of the sphere $P(0,t)$, and the total pressure force acting on the sphere. The pressure clearly has fine-scale behaviour at the initial impact and detachment but these features do not substantially affect the overall force on the system. Most of the work on the drop is due to the quasi-uniform pressure at intermediate times. Figure 2(b) shows the full wavefield during impact.

Figure 4. Comparison between the pressure force $F^*$ felt by the droplet, the pressure $P^*$ along $r=0$ during impact, and the pressed area $A_p^*$. Note initial (inertial) and final (suction) spikes. The consequences of the force on the wider system is well tracked by the lubrication-layer area.

All computations presented were performed in Matlab with ranges $N=2^{11}-2^{12}$, $M=2^8-2^9$. For lower impact velocities this was an over-resolution, whereas the finer grids were needed for stronger impacts. The domain length $L$ varied from $12R_0$ to $16R_0$, the latter ensuring the impact is unaffected by waves reflecting from the boundary, even for longer impacts. The longest simulations took $O(10\,\text {h})$ on a standard personal computer running Matlab, although the use of an adaptive time step would probably reduce this considerably. The only free parameter in the system is $\varepsilon$, the cutoff for lubrication, and this was set to ${\rm 40}\,\mathrm {\mu }{\rm m}$, which corresponds to approximately 5 % the radius of the sphere $R_0$, though values between 10 % and 1 % were also tested and the model displayed insensitivity within this range of $\varepsilon$. The fast decay in pressure in (3.8) as $h$ increases renders the solution insensitive to larger values of the cutoff $\varepsilon$.

4. Conclusions

The LM model herein efficiently resolves certain features of droplet rebounds heretofore unattainable computationally in reduced models. Future work includes computational studies including drop deformation and non-axisymmetric impact, coupled with the continued validation through results compared with the wider literature (e.g. results from studies such as Tang et al. 2018; Alventosa et al. 2023 etc.).

This work follows several (e.g. Galeano-Rios et al. 2017, 2019) which were motivated by Faraday bouncing droplet dynamics (Bush 2015). We expect as future development to use this model with Faraday pilot-wave systems, as an accurate reproduction of the vertical motion of the droplet is thought to be crucial to the observed statistical features resembling quantum mechanics (Durey, Milewski & Wang 2020). The present formulation of a lubrication layer mediating the impact of two bodies of fluid also has potential interest in the mathematics community which has recently rigorously addressed questions of whether fluid surfaces can intersect (Fefferman, Ionescu & Lie 2016).