2. Methodology

Consider the head-on collision of two initially spherical liquid drops of radius $R$ in a cylindrical coordinate system $(r,z,\theta )$ with a plane of symmetry at $z=0$ and axisymmetry assumed, so that only one drop needs to be simulated in the $(r,z)$-plane. The drops initially move at speed $U$ towards the $z=0$ plane, so that their relative speed of impact is $V=2U$. The base of the upper drop starts from height $z=0.1R$, which is far enough for the lubrication pressure to be negligible (this is discussed in Appendix B). Gravity is negligible during the collision process and experiments never consider any dependence on it.

The Navier–Stokes equations are solved for the velocity, $\boldsymbol {u}=(u_r,u_z)$, and pressure, $p$, in the liquid drop

(2.1)$$\begin{gather} \rho_l \frac{\mathrm{D}\boldsymbol{u}}{\mathrm{D}t}=\boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{\mathsf{T}}, \end{gather}$$

(2.2)$$\begin{gather}\boldsymbol{\mathsf{T}}={-} p\boldsymbol{\mathsf{I}}+ \mu_l(\boldsymbol{\nabla}\boldsymbol{u}+\boldsymbol{\nabla}\boldsymbol{u}^\top), \end{gather}$$

(2.3)$$\begin{gather}0=\boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{u}. \end{gather}$$

Here, the liquid viscosity is $\mu _l$ and the stress tensor in the drop is denoted by $\boldsymbol{\mathsf{T}}$. Initially, the drop is a sphere of radius $R$, centre at $(r,z)=(0,1.1R)$, with uniform velocity $\boldsymbol {u}(t=0)=-U\boldsymbol {e}_z$.

At $r=0$, we apply the symmetry boundary conditions $u_r=0$ and ${\partial u_z}/{\partial r}=0$. On the free surface of the drop, a stress balance boundary condition is applied

(2.4)\begin{equation} \boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\mathsf{T}}=\boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\mathsf{T}}_g-\gamma(\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{n})\boldsymbol{n}+\frac{A}{6{\rm \pi} h^3}\boldsymbol{n}. \end{equation}

Here, $\boldsymbol{\mathsf{T}}_g$ is the stress tensor in the gas and $A$ is the Hamaker constant. The entire gas film thickness is $h(r,t)$, so that, due to the symmetry at $z=0$, the upper free surface is then located at $z=h/2$.

The disjoining pressure term in (2.4), which contains the Hamaker constant, accounts for the nanoscale van der Waals attraction between the drops (Pan et al. Reference Pan, Law and Zhou2008). Without this, in this model the drops always bounce (Sprittles Reference Sprittles2024). The Hamaker constant, $A$, depends on the liquid considered but can be measured or computed independently of drop-collision experiments (i.e. we do not ‘fit it’). We have used Hamaker constants for each liquid as calculated by Huang & Pan (Reference Huang and Pan2021) using Lifshitz theory (Israelachvili Reference Israelachvili2011).

The liquid's dynamics is coupled to a lubrication equation that solves for the gas pressure, $p_g(r,t)$, following the method introduced in Chubynsky et al. (Reference Chubynsky, Belousov, Lockerby and Sprittles2020). The lubrication model is composed of an equation enforcing incompressibility of an evolving gas film and one from the momentum equation for the volume flux, $Q_r$, giving

(2.5)$$\begin{gather} 0=\frac{\partial h}{\partial t}+\frac{1}{r}\frac{\partial (rQ_r)}{\partial r}, \end{gather}$$

(2.6)$$\begin{gather}Q_r={-}\frac{\varDelta_Ph^3}{12\mu_{g}}\frac{\partial p_g}{\partial r}+hu^d_r. \end{gather}$$

Here, gas viscosity is denoted by $\mu _g$, gas-kinetic factor by $\varDelta _P(h(r,t))$ (a function of $h$, discussed below) and $u^d_r=u_r(r,z=h(r,t)/2,t)$ is the component of fluid velocity (on the drop surface) in the radial direction. The lubrication flow has two components: a Poiseuille flow due to the pressure gradient, and a plug flow due to the horizontal motion of the drop surfaces.

The lubrication equation is solved on the lower surface of the drop, and provides normal and tangential stresses that are applied as boundary conditions to the Navier–Stokes equations in (2.4)

(2.7)\begin{equation} \boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\mathsf{T}}_g={-}p_g\boldsymbol{n}+ \frac{h}{2}\frac{\partial p_g}{\partial r}\boldsymbol{e}_r. \end{equation}

The second term on the right-hand side of (2.7) is the tangential stress, and comes from the Poiseuille flow component of the lubrication flow.

For (2.5) and (2.6), at $r=0$, we apply the symmetry boundary condition ${\partial p_g}/{\partial r}=0$, and at the edge of the film the gas pressure, $p_g$, is fixed to a reference gas pressure $p_0$. As described in Chubynsky et al. (Reference Chubynsky, Belousov, Lockerby and Sprittles2020), and confirmed by our own tests, the exact position where this edge is set does not affect the results: away from the thin gas film, as $h$ becomes large, ${\partial p_g}/{\partial r}$ will tend to 0, and $p_g$ will tend to $p_0$ before reaching the exact location of the boundary condition. We set it based on where the gradient ${\partial h}/{\partial r}$ first exceeds a threshold of $10$.

The lubrication model assumes that the gas is incompressible. Analysis of our results shows this is acceptable, with excess pressures usually at most 10 % of atmospheric pressure (this is only violated for the smallest drops at the highest $\textit {We}$, which can peak at 25 %). This is in agreement with the estimates for excess pressure in Hicks & Purvis (Reference Hicks and Purvis2011); future analyses could include compressibility as formulated in Mandre et al. (Reference Mandre, Mani and Brenner2009) and Hicks & Purvis (Reference Hicks and Purvis2013).

In (2.6), the term $-({h^3}/{12\mu _{g}})({\partial p_g}/{\partial r})$ is the flow rate for plane Poiseuille flow predicted by the Navier–Stokes equations with no-slip boundary conditions. The rarefied-gas adjustment to this term, derived by Chubynsky et al. (Reference Chubynsky, Belousov, Lockerby and Sprittles2020), is the gas-kinetic factor, plotted in figure 2

(2.8)\begin{equation} \varDelta_P(\textit{Kn})=1+6.88\textit{Kn} +(6\textit{Kn}/{\rm \pi})\ln(1+2.76\textit{Kn}+0.127\textit{Kn}^2), \end{equation}

which is a function of the (local) Knudsen number, $\textit {Kn}=\lambda /h(r,t)$, where $\lambda$ is the mean free path of the gas.

As $\textit {Kn}\to 0$, (2.8) has the asymptotic form (to first order in $\textit {Kn}$) of the flow-rate correction found when, instead of no slip, the Navier slip boundary condition is applied: $\varDelta _P\sim 1+6.88\textit {Kn}$ (Lauga & Stone Reference Lauga and Stone2003). As $\textit {Kn}\to \infty$, it has the correct asymptotic behaviour approaching free molecular flow: $\varDelta _P\sim (6\textit {Kn}/{\rm \pi} )\ln (4\textit {Kn}^2/{\rm \pi} ^2)$ (Cercignani et al. Reference Cercignani, Lampis and Lorenzani2004). The asymptotic expressions are shown in figure 2. The combined expression is derived using asymptotic matching and is fitted to within 2 % of solutions of the Boltzmann equation from Cercignani et al. (Reference Cercignani, Lampis and Lorenzani2004) (also shown in figure 2), and so accurately describes the flow rate of Poiseuille flow across the full range of Knudsen numbers. Here, diffuse scattering of molecules from the boundaries is assumed, which is reasonable for most systems. This expression assumes constant temperature, and the simulations in Cercignani et al. (Reference Cercignani, Lampis and Lorenzani2004) use the linearised Bhatnagar–Gross–Krook operator, which is reasonable for monatomic gases or those with small molecules, such as nitrogen.

Li (Reference Li2016) used a similar but less well-matched expression than (2.8), formulated as an effective viscosity, $\mu _{g,P}(\textit {Kn})=\mu _g/\varDelta _P$; however, they used this to solve the Navier–Stokes equations in the gas film, rather than exploiting the lubrication approach, as in Chubynsky et al. (Reference Chubynsky, Belousov, Lockerby and Sprittles2020) and here.

For an ideal gas, the mean free path is related to the gas pressure $p$, collisional cross-sectional area $d^2$ and temperature $T$, by

(2.9)\begin{equation} \lambda=\frac{k_BT}{\sqrt{2}{\rm \pi} d^2p}. \end{equation}

Additionally, assuming that the excess pressure in the gas film is small compared with the ambient pressure, in calculating the local $\textit {Kn}$ we can use the mean free path at the ambient pressure $p_0$, and assume constant temperature: i.e. we assume that the mean free path is constant along the film and inversely proportional to the ambient pressure. We denote the ratio of the ambient pressure to a reference pressure as $\mathcal {P}=p_0/p_{atm}$ (i.e. the ambient pressure relative to standard atmospheric pressure $p_{atm}=101$,325 Pa). Then for each gas we can use known values for the mean free path at atmospheric pressure $\lambda _{atm}^{gas}$, which differ due to the different collisional cross-sectional area $d^2$, to find

(2.10)\begin{equation} \lambda=\lambda^{gas}_{atm}/\mathcal{P}. \end{equation}

This means that both the gas species and the ambient pressure are captured in our model through the choice of $\lambda$. Values used are $\lambda _{atm}^{air}=69$ nm, $\lambda _{atm}^{nitrogen}=70$ nm and $\lambda _{atm}^{helium}=180$ nm.

Our open-source code (available here, DOI: 10.5281/zenodo.10927009) is implemented in the open-source C++ finite element package oomph-lib (Heil et al. Reference Heil, Hazel and Matharu2022). We solve the Navier–Stokes equations (2.1)–(2.3) in the liquid drop on a triangulated mesh, with the gas manifesting itself through the drop's boundary conditions. An implicit, adaptive time stepper is used. We use oomph-lib's pseudo-solid arbitrary-Lagrangian–Eulerian mesh, where the boundary of the mesh is moved with the normal velocity of the liquid at the drop surface, and the interior elements are moved by treating them as an elastic solid. The time derivatives are adjusted to account for the motion of the mesh (Cairncross et al. Reference Cairncross, Schunk, Baer, Rao and Sackinger2000). Periodically the domain is remeshed using an adaptive method (Zienkiewicz & Zhu Reference Zienkiewicz and Zhu1987) that refines the size of the mesh elements as needed to minimise the inter-element stress jump. Here, our focus is purely on predicting the transitions between bouncing and merging. Capturing the post-contact dynamics is tricky but achievable in the current computational framework – it is essentially a routine computational fluid dynamics calculation, involving no additional physics, but is beyond the scope of this work.

A comparison of our computed drop profiles with experiments in Pan et al. (Reference Pan, Law and Zhou2008) are shown in figure 1, although the real power of our model is in being able to predict the transitions between bouncing and merging, rather than just capture the macroscopic motion of the drops.

Figure 1.

Comparison of our simulations (red outlines) with the experiments of Pan, Law & Zhou (Reference Pan, Law and Zhou2008) (images have been rotated; snapshots in time from left to right, at intervals of approximately 100 $\mathrm {\mu }$s, see Pan et al. (Reference Pan, Law and Zhou2008) for exact timings). Examples of drop bouncing near to the soft transition (decreasing impact speed results in coalescence; (a)) and the hard transition (increasing speed results in coalescence; (b)). Physical parameters are: $\rho =762$ kg m$^{-3}$, $\mu _l=2.128$ mPa s, $\mu _g=18.27$ $\mathrm {\mu }$Pa s, $\gamma =26.5$ mN m$^{-1}$, $A=5.0\times 10^{-20}$ J; above: $R=170.6$ $\mathrm {\mu }$m, $V=0.486$ m s$^{-1}$; below: $R=167.6$ $\mathrm {\mu }$m, $V=0.992$ m s$^{-1}$.

Figure 2.

(a) Schematic showing the initial position of the drop. The dashed blue line shows the position of the second drop, which is included via symmetry. (b) The gas-kinetic factor $\varDelta _P$ as a function of $\textit {Kn}$, showing the simulation data from Cercignani, Lampis & Lorenzani (Reference Cercignani, Lampis and Lorenzani2004) used in Chubynsky et al. (Reference Chubynsky, Belousov, Lockerby and Sprittles2020) in fitting it. The dashed lines show the asymptotic expressions for large and small $\textit {Kn}$.

3. Comparison with experimental results of Huang & Pan (Reference Huang and Pan2021)

Table 1 presents quantitative predictions for $\textit {We}_S$ and $\textit {We}_H$ from our computational model and compares these with experiments from Huang & Pan (Reference Huang and Pan2021). The physical parameters used, taken from Huang & Pan (Reference Huang and Pan2021), are given in the table's caption. Note that we define $\textit {We}$ in terms of the drop radius, so that $\textit {We}=\rho _lV^2 R/\gamma$; most experimental articles including Huang & Pan (Reference Huang and Pan2021) use the drop diameter, so that the values reported in those articles are twice those referenced here. Huang & Pan (Reference Huang and Pan2021) reported $\textit {We}_S$ and $\textit {We}_H$ for both the lower-speed transition from merging to bouncing with increasing $\textit {We}$ (the ‘soft transition’ at $\textit {We}_S$) and the higher-speed transition from bouncing to merging (the ‘hard transition’ at $\textit {We}_H$). For comparison, also listed are values found for $\textit {We}_S$ and $\textit {We}_H$ when the GKE corrections are removed by setting $\varDelta _P=1$.

Table 1.

Values of the critical Weber numbers, $\textit {We}_S$ and $\textit {We}_H$, found in our simulations with and without gas-kinetic effect (GKE) corrections, in comparison with the experimental values found in Huang & Pan (Reference Huang and Pan2021). Entries labelled NB (‘no bounce’) show where no transition between merging and bouncing was found, as the drops always merged. The starred entry shows where no soft transition was observed for the values of $\textit {We}$ used in the experiments. The simulations results are $\pm 5$ of the smallest digit shown (e.g. 0.265 means the value lies in the interval 0.26–0.27; 4.05 means the interval 4.0–4.1).

The physical parameters of the liquids and gas, shown below, are those used in Huang & Pan (Reference Huang and Pan2021). Decane: $\rho =728$ kg m$^{-3}$, $\mu _l=0.82$ mPa s, $\gamma =23.8$ mN m$^{-1}$, $A=5.0\times 10^{-20}$ J; dodecane: $\rho =750$ kg m$^{-3}$, $\mu _l=1.33$ mPa s, $\gamma =24.9$ mN m$^{-1}$, $A=5.0\times 10^{-20}$ J; tetradecane: $\rho =759$ kg m$^{-3}$, $\mu _l=2.05$ mPa s, $\gamma =26.0$ mN m$^{-1}$, $A=5.2\times 10^{-20}$ J; water: $\rho =998$ kg m$^{-3}$, $\mu _l=0.92$ mPa s, $\gamma =72.0$ mN m$^{-1}$, $A=3.7\times 10^{-20}$ J; air: $\mu _g=18.6$ $\mathrm {\mu }$Pa s, $\lambda =69.0$ nm.

What is immediately apparent is that the simulations which exclude the GKE predict transition values that are wildly different from the experimentally observed ones. Specifically, the bouncing window is massively extended, with values for $\textit {We}_S$ often three orders of magnitude too low and $\textit {We}_H$ often double what it should be. Without GKE, regimes in which no bouncing (NB) is seen experimentally cannot be reproduced, making this model particularly inept for small drops, where merging becomes the dominant outcome of a collision. In contrast, simulations with GKE routinely predict the correct magnitudes for the transitions and are able to reproduce NB regimes. We now analyse these results in further detail.

For the $R=300$ $\mathrm {\mu }$m hydrocarbon (decane, dodecane and tetradecane) drops, our simulation results agree very well with the experimental values for $\textit {We}_S$ and $\textit {We}_H$. Note that no soft transition was observed experimentally for $R=300$ $\mathrm {\mu }$m tetradecane; this makes sense, as our computed $\textit {We}_S$ is below the minimum $\textit {We}$ possible in those experiments. Our simulations also agree very well with the hard transition for $R=500$ $\mathrm {\mu }$m water and $R=150$ $\mathrm {\mu }$m dodecane drops. For all of the small drops, both our simulations and the experiments found that NB regime can be observed, in contrast to the simulations without the gas-kinetic factor when a bouncing regime was always found, highlighting again that the GKEs are essential to predicting these bounce–merge transitions.

The agreement between the simulations and the experiments is less good for the intermediately sized drops, and, unlike Huang & Pan (Reference Huang and Pan2021), for $R=150$ $\mathrm {\mu }$m decane and $R=115$ $\mathrm {\mu }$m tetradecane we observe the NB regime. However, at these radii we predict that $\textit {We}_S$ and $\textit {We}_H$ are particularly sensitive to radius, as shown in figure 6 and discussed in § 6. Our values for $R=115$ $\mathrm {\mu }$m tetradecane are very close to the experiments when plotted on the $R-\textit {We}$ regime diagram. The agreement is also less good for $R=150$ $\mathrm {\mu }$m tetradecane, and the soft transitions for $R=150$ $\mathrm {\mu }$m dodecane and $R=500$ $\mathrm {\mu }$m water, and we cannot currently explain the cause of these discrepancies. However, the agreement is still far superior to the simulations without GKE, see table 1.

Focusing on a particular set of results, figure 3 shows our simulations for the decane $R=300$ $\mathrm {\mu }$m drops for $\textit {We}$ either side of $\textit {We}_S$ and $\textit {We}_H$. The main plot shows the minimum film thickness on a log scale, with videos in the Supplemental Material showing the dynamics. Note that we regularly see film thicknesses less than $100$ nm, compared with the mean free path of air at atmospheric conditions of $69$ nm, justifying the need for gas-kinetic corrections.

Figure 3.

Time evolution of minimum film height at $\textit {We}$ just above and below $\textit {We}_S=0.265$ (blue lines) and $\textit {We}_H=6.95$ (orange lines); for $R=300$ $\mathrm {\mu }$m decane drops. Drop profiles are shown at key stages of the collision for the Weber numbers ($\textit {We}=0.3$ and $\textit {We}=6.5$) in the bouncing regime close to the transition (videos are available as Supplementary Material Movies 1 and 2 available at https://doi.org/10.1017/jfm.2024.722). Below the drop profiles, the lower boundary of the drop is shown, zoomed vertically by a factor of 500, so that the shape of the gas film can be seen. The black dashed line is the critical height (3.1) predicted by Chubynsky et al. (Reference Chubynsky, Belousov, Lockerby and Sprittles2020).

In blue ($\textit {We}=0.3$, shown in the top panels of figure 3) are snapshots of the simulation for $\textit {We}$ just above $\textit {We}_S$, with a corresponding video in Supplementary Material Movie 1. For $\textit {We}$ just below $\textit {We}_S$, contact occurs once the bulk of the drops has starting receding (at the point labelled v in figure 3), as previously found in Bach et al. (Reference Bach, Koch and Gopinath2004). Contact happens at the ‘kink’ at the edge of the film, with a significant gas bubble trapped within.

In orange ($\textit {We}=6.5$, shown in the bottom panels of figure 3) are snapshots of the simulation for $\textit {We}$ just below $\textit {We}_H$, with the video given in Supplementary Material Movie 2. Here, the minimum height is reached at the kink at the edge of the gas film (point III in figure 3), and it is here that contact happens for $\textit {We}$ just above $\textit {We}_H$. The gas film is noticeably affected by the macroscopic behaviour of the drop, such as at point IV in figure 3, where the top of the drop has become concave before rebounding, causing the oscillations in film height seen at point IV. These effects can mostly clearly seen in Movie 2. This film dynamics has not been seen experimentally, unlike for drop-on-solid impacts where interferometry has been able to measure it (de Ruiter, Van Den Ende & Mugele Reference de Ruiter, Van Den Ende and Mugele2015c).

The black line shows the predicted critical film height $h_c$, derived in Chubynsky et al. (Reference Chubynsky, Belousov, Lockerby and Sprittles2020), at which the disjoining pressure becomes dominant, overcoming the surface tension, and pulls the interfaces together (i.e. creates contact). This is based on a linear stability analysis of the lubrication equation, which is complicated by the height-dependent gas-kinetic factor $\varDelta _P$. Notably, the critical film height does not depend on the impact speed or the drop radius, and is calculated by solution to

(3.1)\begin{equation} \frac{h_c^5}{\varDelta_P(\lambda/h_c)}=\frac{A^2}{96 {\rm \pi}^2\sigma_{max}\gamma \mu_g}, \end{equation}

where $\sigma _{max}$ is an estimated growth rate of the fastest growing mode, which we must impose. In an experiment, the instability will be initiated by a thermal fluctuation, whereas in a simulation it will be a small numerical error. We use an estimate for the size of a thermal fluctuation to approximate this, as described in Chubynsky et al. (Reference Chubynsky, Belousov, Lockerby and Sprittles2020). For decane at room temperature, the characteristic size of interfacial thermal fluctuations is $h_{fluc}\sim 0.1$ nm (Sides, Grest & Lacasse Reference Sides, Grest and Lacasse1999). This must grow to cover the gas film of thickness $h\sim 20$ nm in a characteristic time of $t_{c}\sim 0.1$ ms (these are approximations based on the simulations shown in figure 3). This means that, for an initial perturbation which grows exponentially to span the gap during the collision period, we require $\sigma _{max}=\log (h/h_{fluc})/t_c=5\times 10^4$ s$^{-1}$. The resulting estimate for $h_c$ is plotted in figure 3. The theory predicts that, if the minimum film height falls below $h_c$, the film becomes linearly unstable and the disjoining pressure can drive the drops into contact, as shown for $\textit {We}=0.2$, 0.25, 7.0 and 7.5 in figure 3.

4. Influence of gas pressure

Having validated our approach at atmospheric pressures, we now consider whether our model is able to capture how transitions between bouncing and merging vary with gas type and ambient pressure, for fixed temperature – these effects have only been considered experimentally in Qian & Law (Reference Qian and Law1997). There, tetradecane drops with a radius of $100$–$200$ $\mathrm {\mu }$m were considered for nitrogen and helium environments at different pressures. Given the small variation in viscosity for these two gases (differing only by approximately 10 %), we fix $\mu _g$ at an intermediate value of $18.4$ $\mathrm {\mu }$Pa s and run simulations with varied ambient gas pressure, which we recall only appears in the model due to the incorporation of GKEs (imposed via the mean free path) and varied $\textit {We}$ (by changing the impact speed, fixing $R=150$ $\mathrm {\mu }$m).

The inset to figure 4 shows our results for $\textit {We}_S$ and $\textit {We}_H$ as a curve plotted against the ambient gas pressure, with merging on one side of the curve and bouncing on the other. This shows excellent agreement with the experiments of Qian & Law (Reference Qian and Law1997) and Huang & Pan (Reference Huang and Pan2021), for both nitrogen and helium gases. For a large enough fixed pressure the curve becomes multivalued, showing the merge–bounce–merge transitions. At sufficiently low pressure, drops always merge and here this occurs very near atmospheric pressure for nitrogen (and hence air also). Surprisingly, given good agreement elsewhere, one higher-pressure nitrogen experiment from Qian & Law (Reference Qian and Law1997) does not agree with our simulations, and in fact shows $\textit {We}_H$ decreasing with gas pressure. We do not have an explanation for this, and think it is worthy of further experimental analysis: note there are relatively few data points shown in the relevant figure 7(c) in Qian & Law (Reference Qian and Law1997).

Figure 4.

Regime diagram for $R=150$ $\mathrm {\mu }$m tetradecane drops. The Weber number is varied by changing the velocity. The white region is where the drops bounce. The contour shows the relative contact position, $r_c/r_{max}$, with labelled contact modes. The inset shows the experimental and computed transitions for the different gases at the gas pressure relative to the reference atmospheric pressure. Above are the drop profiles at the moment of contact for characteristic contact modes, at the numbered points in the regime diagram. Points 1, 2a, 3a and 3b are all for $\textit {We}=12$ and with $\lambda$ corresponding to $\mathcal {P}= 0.6, 1, 1.9, 2.0$ (for air), respectively. Point 2b is $\textit {We}=0.5$, $\mathcal {P}=1$. Below the drop profiles, the shape of the gas film is shown (zoomed vertically by a factor of 500) so that the position of contact can be seen. The regions of the contour corresponding to each contact mode have been separated by the dashed white lines, and are discussed in §5.

The main section of figure 4 shows $\textit {We}_S$ and $\textit {We}_H$ against the inverse mean free path. The white region is where the drops bounce, and the coloured region where they merge. The gas type and ambient pressure only appear in the model through the mean free path, calculated from the gas pressure using (2.10). The experimental values have also been plotted in terms of their mean free path ($\lambda =70$ nm for nitrogen and $\lambda =180$ nm for helium), using (2.10) to account for both ambient pressure and gas type. A great boost for the model proposed, with GKEs at its core, is the collapse of experimental data onto a single transition curve showing that the dependence on both ambient pressure and gas species can be characterised purely by the mean free path. As expected, as $\mathcal {P}\to \infty$ (or equivalently, $\lambda \to 0$, where the mean free path will be much smaller than the gas film heights), $\textit {We}_S$ and $\textit {We}_H$ tend to those that are found when the gas-kinetic factors are removed.

Previous analyses of drop impact on solids and liquid films have attempted to build regime maps based on the modes of contact (Zhang et al. Reference Zhang, Soori, Rokoni, Kaminski and Sun2021a). Here, by covering all of the relevant parameter space with simulation data (totalling 798 simulations, shown in figure 10), we introduce a new quantitative method for presenting this information and identifying changes in contact mode. This is done by contour plotting relevant measures for the position and time at which contact occurs. For position, in figure 4 the relative contact position is given, i.e. the radial position of contact, $r_c$, divided by the maximum radial extent of the drop $r_{max}$ at the time of contact, see illustration in figure 4. For time, in figure 5, the normalised contact time indicates how far through the bouncing process contact occurs: the ratio of the contact time, $t_c$, to the equivalent time taken for the drops to bounce, $t_b$, when the disjoining pressure is removed (this the time until the drops are again at their initial separation distance). Notably, this denominator, $t_b$, does not depend on the gas pressure, and is plotted in the inset, along with theoretical predictions for large and small $\textit {We}$.

Figure 5.

Regime diagram for $R=150$ $\mathrm {\mu }$m tetradecane drops, as in figure 4, but instead the contour shows normalised contact time, the ratio of the contact time to the equivalent time for the drops to bounce without the disjoining pressure (the time until the drops are again at their initial separation distance). The white region is where the drops bounce, and the contact modes shown in figure 4 are labelled. The bounce time is plotted in red in the inset, where the dashed lines showing asymptotic predictions: for large $\textit {We}$ the limit to the natural frequency of the drop $(t=2.2(\rho R^3/\gamma )^{1/2})$ (Rayleigh Reference Rayleigh1879), and for small $\textit {We}$ the prediction from Gopinath & Koch (Reference Gopinath and Koch2002).

5. Contact modes

The contours in figures 4 and 5 show that there are different contact modes: distinct ways that the drops can come into contact. Characteristic drop and film profiles at the moment of contact are shown at the top of figure 4. Where contact occurs will have a significant effect on how the drops merge beyond the initial contact, and the size of the entrapped gas bubble. There are 3 distinct modes that can be seen, which we call the ‘inner kink’ mode, the ‘outer kink’ mode and the ‘film’ mode. To understand these modes and the boundaries in phase space between them, we can look at the contours in figures 4 and 5 together.

The inner kink mode (with characteristic profile shown at point 1 in figure 4) occurs at high $\textit {We}$ and/or low pressure. Contact occurs at a ‘kink’ near the centre of the drop where the local curvature is high, where the entrapped gas bubble meets with an expanding film (as described by Mandre et al. Reference Mandre, Mani and Brenner2009), with this mode having the smallest relative contact position of all the modes. The normalised contact time is also small, as contact happens before the deformed surface of the drop extends beyond the initial radius of the drop. This is the high-speed mode seen for drop-on-solid impact, e.g. in de Ruiter et al. (Reference de Ruiter, Van Den Ende and Mugele2015c), and the ‘film’ mode, seen in Zhang et al. (Reference Zhang, Soori, Rokoni, Kaminski and Sun2021a).

The outer kink mode (with characteristic profile shown at point 2a in figure 4) occurs to the right of the inner kink mode on the regime diagrams. Here, the drop is still ‘spreading’, and contact now occurs at a kink at the edge of the film, where the film meets the more spherical part of the drop. This mode is equivalent to the lower-speed contact seen for drop-on-solid impact in de Ruiter et al. (Reference de Ruiter, Van Den Ende and Mugele2015c), the ‘kink’ mode seen by Zhang et al. (Reference Zhang, Soori, Rokoni, Kaminski and Sun2021a) and the ’kink’ mode found numerically by Chubynsky et al. (Reference Chubynsky, Belousov, Lockerby and Sprittles2020). It is also the mode of the simulations for $\textit {We}$ above $\textit {We}_H$ in figure 3.

The outer kink mode also extends below the soft transition, with characteristic profile shown at point 2b in figure 4. The normalised contact time is largest here, with contact happening once the drops start rebounding, as shown in figure 3.

There is a clear jump in the relative contact position across the boundary between the outer kink mode and the inner kink mode. In this part of the regime diagram, during the part of the collision when contact occurs there are two local minima in the film: the two ‘kinks’ where contact occurs for the inner and outer kink modes. This can be seen most clearly in the zoomed gas film profile at point 2a in figure 4. The boundary between the kink modes is formed due to a change in the kink which has minimum film height at the moment of contact. A flip in where this minima is can be seen in figure 3: the snapshots of the $\textit {We}=6.5$ simulation show that at point I the minima is relatively near the centre of the drop (at the inner kink), but that by point II, the minima is at the outer kink. For drop-on-solid impacts this two-kink behaviour was observed experimentally and discussed in detail by de Ruiter et al. (Reference de Ruiter, Van Den Ende and Mugele2015c), and contact occurring at both an inner and an outer kink is seen by Li et al. (Reference Li, Langley, Tian, Hicks and Thoroddsen2017).

The film mode has two characteristic contact profiles shown at points 3a and 3b in figure 4. Here, the drop has spread significantly, and contact can occur either in the middle of the gas film (3a) or at the edge (3b). This is the ‘film mode’ at high pressures seen in Chubynsky et al. (Reference Chubynsky, Belousov, Lockerby and Sprittles2020), and appears similar to the contacts seen in impacts of ultra-viscous drops in Langley et al. (Reference Langley, Li and Thoroddsen2017), where an extended air film formed before contact occurred. The exact position of contact depends on where numerical noise triggers the disjoining pressure driven instability, but in reality thermal fluctuations could trigger this film collapse (Fetzer et al. Reference Fetzer, Rauscher, Seemann, Jacobs and Mecke2007; Sprittles et al. Reference Sprittles, Liu, Lockerby and Grafke2023). This is why the film mode region of the contour in figure 4 is ‘speckled’, i.e. non-uniform in colour: each simulation has made contact at a different noise-derived point, which is sensitive to the precise numerical set-up. From a practical perspective, we may expect the position of contact seen in experiments to be stochastic in nature in this region, as is indeed the case for the similar drop–solid experiments in Langley et al. (Reference Langley, Li and Thoroddsen2017).

Across the boundary between the outer kink mode and the film mode, there is a jump in the normalised contact time, most clearly at the moderate $\textit {We}$ where the film mode first appears. This is because when the outer kink contact disappears, the drop spends time spreading before the film mode collapse can occur. This also highlights that having information from both temporal and spatial information for the contact helps define modes in a less ambiguous manner.

For drop-on-solid impact, Chubynsky et al. (Reference Chubynsky, Belousov, Lockerby and Sprittles2020) also observed contact at the very end of the bouncing process, during the collapse of the dimple where gas is trapped between the drops (as shown at points vii and VII in figure 3). This occurred for high-pressure, low-speed impacts, but we have not observed this in our drop-collision simulations. It is yet to be determined if this is a fundamental difference between drop collisions and drop-on-solid impacts, or is simply due to the different parameter regimes of the studies.

6. Dependence on drop size

The majority of studies have focused on how transitions depend on impact speed, and have thought of $\textit {We}$ as a dimensionless measure of this. Here, to systematically study the dependence of $\textit {We}_S$ and $\textit {We}_H$ on drop radius, considered experimentally in Huang & Pan (Reference Huang and Pan2021), we simulate the collision of tetradecane drops at atmospheric pressure (fixed $\lambda$) with varying drop radius and impact speed. Figure 6 shows the resulting regime map in terms of $\textit {We}$ and drop radius, with the contour showing relative contact position as in figure 4. It shows very good agreement with the experimental results of Huang & Pan (Reference Huang and Pan2021), and that for a fixed $\textit {We}$ the drops can bounce or merge depending on $R$.

Figure 6.

Above: regime diagram for head-on collisions of tetradecane drops at atmospheric pressure (and therefore fixed $\lambda$). For a particular radius, $\textit {We}$ is varied by changing only the impact velocity. The contour shows the relative contact as in figure 4. The dashed red line shows $\textit {We}_S$ and $\textit {We}_H$ when GKEs are not considered. Below: the lower section of the above plot, but with a log scale on the y-axis to show the soft transitions more clearly. Also included is the merging–halting transition discussed in § 8.

In figure 6, the boundary between the inner kink (light region) and outer kink modes (dark region) can clearly be seen. The film mode does not appear here, and was only observed for higher ambient pressures than the atmospheric pressure used to create this figure.

The dashed red lines in figure 6 show what happens to the hard transition when GKEs are removed, by setting $\varDelta _P=1$. This curve is dramatically different – for $R=150$ $\mathrm {\mu }$m, it is well over twice the experimental value for $\textit {We}_H$, and increases as the radius decreases, despite the experiments and our primary simulations showing that the bouncing regime vanishes for small drops. Moreover, without GKEs, $\textit {We}_S$ is below $\textit {We}=0.001$, far beneath the soft transitions found experimentally or by the simulations with GKE, as shown in table 1 and the bottom plot of figure 6.

As $R$ becomes large, $\textit {We}_H$ tends to a constant value. For simulations up to 1 mm, $\textit {We}_H$ tends to $10.5\pm 0.5$, and without GKEs, $\textit {We}_H$ tends to $12.5\pm 0.5$. Here contact occurs in the film mode, once the drop has spread significantly. Contact is determined by the stability of the extended gas film. If the height of the film is below the critical height for the van der Waals driven instability, contact will occur. Our results suggest that (for large $R$ at least) the stability of the film depends on the radius only via the Weber number. The height (and stability) also depend on the mean free path, and hence the values of $\textit {We}_H$ with and without GKEs do not converge for large $R$.

The transition boundary from figure 6 is plotted again in figure 7, this time in terms of velocity and drop radius. This reveals that for some values of the velocity (here, around $V=1$ m s$^{-1}$) there is also a merge–bounce–merge transition with increasing radius (see the arrow in figure 7). The experiments of Huang & Pan (Reference Huang and Pan2021) also demonstrate this, but they only have 3 data points (as shown in figure 7), and did not present this in terms of velocity. By systematically varying the radius and presenting the results in this way, we are able to identify this merge–bounce–merge transition with fixed velocity.

Figure 7.

The regime diagram for tetradecane, as in figure 6, but in dimensional variables, showing critical velocities as a function of drop radius, at atmospheric pressure. The blue arrow shows that for fixed velocity, there is an additional merge–bounce–merge transition.

Both the transitions found with fixed radius (in figures 4 and 5) and with fixed $\lambda$ (figure 6) could be plotted against an inverse global Knudsen number $R/\lambda$. However, the curves are similar but not the same, demonstrating that $\textit {We}_S$ and $\textit {We}_H$ are not simply a function of a global Knudsen number of the ratio of $\lambda$ and $R$.

The theory of Chubynsky et al. (Reference Chubynsky, Belousov, Lockerby and Sprittles2020) (the solution of (3.1)) tells us whether or not the film is stable at a given height; however, without being able to predict this height for a given set of parameters it cannot tell us a priori whether or not a given pair of colliding drops will merge or bounce; clearly such a theory would be highly sought after. Some previous analytic models, such as those of those of Mandre et al. (Reference Mandre, Mani and Brenner2009), Huang & Pan (Reference Huang and Pan2021) and Sharma & Dixit (Reference Sharma and Dixit2021), have predicted the minimum film heights, but these are yet to incorporate full GKEs which we have seen is essential. The model of Zhang & Law (Reference Zhang and Law2011) includes a gas-kinetic factor and predicts both $\textit {We}_S$ and $\textit {We}_H$, but it uses a parameter fitted to experimental results and requires solving a nonlinear ordinary differential equation. Bach et al. (Reference Bach, Koch and Gopinath2004) considered $\textit {We}_S$ for small drops with small deformations, predicted a minimum height incorporating Knudsen effects, and used the results of their simulations to fit a critical height, comparing favourably with the experiments of Qian & Law (Reference Qian and Law1997).

An ideal model for predicting the bounce–merge transitions would characterise $\textit {We}_S$ and $\textit {We}_H$ as functions of a rarefaction parameter $\lambda /H$, where $H$ is a characteristic length scale of the gas film height. Even then, however, figures 4 and 5 show us how complex the situation is, with different modes of contact that would have to be predicted by any simplified theory hoping to understand the instability mechanism. Furthermore, such a model should predict both soft and hard transitions, meaning most simple expressions will fail.

García-Geijo, Riboux & Gordillo (Reference García-Geijo, Riboux and Gordillo2024) compared several models for the minimum height of the film during the initial approach of a drop towards a solid, but these theories could easily be adapted to drop–drop collisions. Combined with a theory for the critical height at which the disjoining pressure causes contact, this could be a step towards such a model. However, at present, they only predict the height of the ‘inner kink’, and we have shown that (at least in the parameter regime considered), the bounce–merge transitions are between contact at later times, such as at the ‘outer kink’.

7. Water drops

Similar to the results for tetradecane in § 4, we consider the bounce–merge transitions for water in figure 8. Again, we plot $\textit {We}_S$ and $\textit {We}_H$ against the inverse mean free path, and compare with the experiments of Qian & Law (Reference Qian and Law1997) for different gases and ambient pressures. These experiments varied $\textit {We}$ by varying both the drop radius (between 100 and 200 $\mathrm {\mu }$m) and the impact speed, and the results are reported only in terms of $\textit {We}$. As with the hydrocarbons discussed above, for a fixed value of $\textit {We}$, the drops can bounce or merge depending on the drop radius. We therefore vary $\textit {We}$ by only changing the impact speed, and run two sets of simulations: for $R=100$ $\mathrm {\mu }$m and $R=150$ $\mathrm {\mu }$m.

Figure 8.

Computed transitions for water drops collision of two different radii, compared with the experiments of Qian & Law (Reference Qian and Law1997).

At lower pressure, the experimentally reported transitions lie between the two bounce–merge transition curves found from our model, demonstrating the sensitivity to the drop radius, and the need for experiments to report the radius and velocity independently. The agreement is less good for the higher-pressure (larger $\lambda ^{-1}$) experiments, similar to the high-pressure tetradecane in nitrogen experiment discussed in § 4. Again, there is a large amount of experimental uncertainty in the relevant figures 6(c) and 8(c) of Qian & Law (Reference Qian and Law1997), and this regime is clearly worthy of further experimental study.

8. Very-low-speed impacts

For sufficiently low impact speeds, there is an additional transition that cannot be seen in the figures above. Here, the lubrication pressure is enough to slow the drops until they are stationary. This is shown for tetradecane in figure 9, with halting occurring between $\textit {We}=10^{-6}$ and $\textit {We}=10^{-5}$. A critical Weber number curve for this transition is also shown in the bottom part of figure 6. This was predicted theoretically by Bach et al. (Reference Bach, Koch and Gopinath2004), who suggest a value of $\textit {We}=9.8\times 10^{-5}$ based on a fit for the coefficient of restitution. Note that our simulations include only the lubrication pressure from the gas film, not the Stokes drag from the surrounding air which would cause halting to occur for higher $\textit {We}$ than seen here. The velocities here are far smaller than those investigated experimentally with free drops, and so this transition has not yet been observed. Altogether then, there can be three transitions in the outcome of drop collisions with increasing $\textit {We}$: starting from $\textit {We}=0$, with increasing $\textit {We}$ the drops first fail to merge (‘halting’), then merge, then bounce, then merge again!

Figure 9.

Very-low-speed $R=150$ $\mathrm {\mu }$m tetradecane collisions at atmospheric pressure, showing the transition from merging to halting.

10. Future directions

The results in this article suggest some regimes that would be of particular interest for future experimental investigation. In particular, experiments with controlled drop size at different ambient pressures, looking at: (i) the critical pressure below which the bouncing regime disappears, and (ii) the trend towards constant $\textit {We}_H$ for high ambient pressure.

In this article we consider only head-on collisions. Notably, for off-centre collision some models for $\textit {We}_S$ and $\textit {We}_H$ are fitted through the values for the head-on transitions (Hu et al. Reference Hu, Xia, Li and Wu2017; Al-Dirawi & Bayly Reference Al-Dirawi and Bayly2019). Clearly, however, to capture off-centre collisions correctly, the computational model must become three-dimensional, with the lubrication equations becoming two-dimensional. Here, it is the dramatic increase in computational cost that is the main limiting factor, but our lubrication approach should help to reduce this burden and the extension to three dimensions should be the focus of important future work.

Our computations are currently only valid up to the moment of contact, and it would be necessary to perform ‘numerical surgery’ on the mesh to continue the simulation beyond that point (Anthony et al. Reference Anthony2023). This is certainly possible and will be an avenue for future work. An alternative approach would be to use our results at the moment of contact as an initial condition to a solver better suited to simulations involving changes in topology, but that cannot predict when and where the initial contact will occur.

The model described in this paper can also easily be applied to impacts of drops onto smooth solids, as in Chubynsky et al. (Reference Chubynsky, Belousov, Lockerby and Sprittles2020). Recent experiments have been able to measure the thickness of the gas layer with remarkable accuracy (Kaviani & Kolinski Reference Kaviani and Kolinski2023), and comparisons could easily be made using this model.

Beyond drop–drop collisions and drop-on-solid impacts, there are good prospects of including additional physical processes. Drop impacts onto films could be considered with the addition of a second lubrication equation for the liquid film (Duchemin & Josserand Reference Duchemin and Josserand2020). Thin liquid films have been used in experiments as an analogue for smooth solids, but unexplained discrepancies between the behaviours of solid and films have been identified (Lo et al. Reference Lo, Liu and Xu2017; Pack et al. Reference Pack, Hu, Kim, Zheng, Stone and Sun2017; Zhang et al. Reference Zhang, Soori, Rokoni, Kaminski and Sun2021b), with different contact modes observed. By incorporating a model for a liquid film, our methodology could isolate the differences of a film as compared with a smooth solid.

The impact of a drop onto a liquid bath could also be modelled, using a second meshed domain for the bath solving the Navier–Stokes equations, along with the required adjustments to the lubrication equation. In Tang et al. (Reference Tang, Saha, Law and Sun2019), the dynamics of the trapped gas film was measured in an impact with a bath of thickness similar to the drop radius. Our methodology could investigate these dynamics, identify contact modes and predict the critical Weber numbers for merging to occur. By combining this approach with a liquid thin-film model, a regime diagram of Weber number vs bath thickness could be constructed.

As investigated by Pan et al. (Reference Pan, Law and Zhou2008), the critical Weber numbers are also sensitive to the value of the Hamaker constant. We have used Hamaker constants for each liquid as calculated by Huang & Pan (Reference Huang and Pan2021) using Lifshitz theory (Israelachvili Reference Israelachvili2011), and it is not fitted to any drop-collision results. This constant will be different for each pair of contacting materials, which must be considered carefully when comparing drop-on-solid and drop-on-film impacts.

If a drop is placed above a solid surface which is hotter than the liquid's boiling temperature, it can levitate on a layer of its own evaporating vapour, in the ‘Leidenfrost effect’ (Quéré Reference Quéré2013). The model described in this article has previously been extended to quasi-static drops in Leidenfrost conditions (without GKEs), and used to predict the equilibrium shapes of these drops (Chakraborty, Chubynsky & Sprittles Reference Chakraborty, Chubynsky and Sprittles2022). Transitions between bouncing and contact have been observed for impacts in Leidenfrost conditions (Tran et al. Reference Tran, Staat, Prosperetti, Sun and Lohse2012; Chantelot & Lohse Reference Chantelot and Lohse2021), though at higher Weber numbers than for impacts in isothermal conditions. Leidenfrost drops that are statically levitating have also been shown to start spontaneously ‘trampolining’ (Liu & Tran Reference Liu and Tran2020; Graeber et al. Reference Graeber, Regulagadda, Hodel, Küttel, Landolf, Schutzius and Poulikakos2021). To investigate these dynamic processes, the model could be coupled to heat flow inside the drop, along with a more complex evaporation model. This would allow the thermal Marangoni effect to be included, which has been predicted to have a significant effect on the stability of Leidenfrost drops (van Limbeek et al. Reference van Limbeek, Ramírez-Soto, Prosperetti and Lohse2021).