1. Introduction
The impact of a liquid droplet with another droplet or onto a solid surface are important basic processes that occur in many applications and in nature. In plants, pathogens are transported by raindrops, which, upon impact with infected leaves, produce an air vortex ring and disperse plant pathogens (e.g. spores of fungi) away from the host plant (Kim et al. Reference Kim, Park, Gruszewski, Schmale III and Jung2019). In applications involving spray devices, the outcome of droplet collisions plays an important role in the development of the number density and the size distribution of the droplets in dense sprays, which in turn substantially affect the subsequent vaporization and combustion of the spray. There are several collision outcomes in binary droplet collision, namely coalescence, bouncing, reflexive separation, stretched separation and splatter (Pan, Chou & Tseng Reference Pan, Chou and Tseng2009; Tang, Zhang & Law Reference Tang, Zhang and Law2012). In applications such as drop-on-demand technologies, droplet bouncing is unfavoured as small quantities of fluid need to be delivered to specific positions on a surface (Zaugg & Wagner Reference Zaugg and Wagner2003; Poozesh et al. Reference Poozesh, Saito, Akafuah and Graña-Otero2016). In agriculture, pesticide droplets during spraying may collide with the previously sprayed droplets, dew or rainwater on leaves. In these collisions, the colliding droplet either stays attached or bounces away and causes contamination of soil (Bergeron et al. Reference Bergeron, Bonn, Martin and Vovelle2000; Gilliom et al. Reference Gilliom2006; Massinon & Lebeau Reference Massinon and Lebeau2012).
The earliest work on the bouncing of a water droplet from a pool of water was carried out by Reynolds (Reference Reynolds1881). Reynolds observed that a water droplet can hover over a pool before sinking, owing to the presence of the air layer between them. A detailed investigation of various aspects of droplet impact phenomena, such as splashing, was experimentally performed by Worthington (Worthington Reference Worthington1908), which sparked the interest of other researchers in this field. Further advancement in this field was achieved by Schotland (Schotland Reference Schotland1960), who performed experiments on the coalescence of water droplets at a gas–water interface. Schotland suggested that the outcome of the collisions was predominantly controlled by fluid and gas inertia and depended on the Weber number (defined as 
$We = \rho U^{2}R/ \gamma$ for droplets with equal size, where 
$\rho$, 
$U$, 
$R$ and 
$\gamma$ are density, relative impact velocity, droplet radius and surface tension, respectively), which measures the relative importance of the droplet inertia compared with the surface tension. Droplets with sufficiently high inertia coalesced whereas lower-Weber-number droplets bounced. Droplet bouncing on an undisturbed interface at variable impacting angles was first studied in detail by Jayaratne & Mason (Reference Jayaratne and Mason1964) experimentally. They established a relationship between the drop radius, impact speed and impact angle at which the bouncing–coalescence threshold occurs between droplets. Building on the work of Gopinath & Koch (Reference Gopinath and Koch2001), Bach, Koch & Gopinath (Reference Bach, Koch and Gopinath2004) studied the droplet impact of small aerosol droplets (
$R \leq 50$ 
$\mathrm {\mu }$m) impacting a fluid bath. They developed a rarefied gas model to describe the dynamics of the gas layer separating the droplet and bath, and used an inviscid potential flow model to describe the transfer of energy from the droplet to the bath during impact. It was found that the criterion for drop bouncing is more sensitive to the gas mean-free path and gas viscosity than to the Weber number itself. More recently, droplets bouncing repeatedly on a vertically oscillated bath have received considerable interest (Couder et al. Reference Couder, Protiere, Fort and Boudaoud2005; Bush & Oza Reference Bush and Oza2020).
Most of the earliest studies were focused on the impact of droplets with a pool of fluid due to the experimental ease and retention of most of the physics applicable in related applications. With rapid development in high-speed digital photography, studies related to droplet collisions have become popular. Binary droplet collisions have been extensively studied in the past. Early experimental studies conducted with water at atmospheric pressure observed two distinct collision outcomes: permanent droplet coalescence and droplet separation after temporary coalescence (Adam, Lindblad & Hendricks Reference Adam, Lindblad and Hendricks1968; Ashgriz & Poo Reference Ashgriz and Poo1990). It was established that the collision outcomes are mainly influenced by the collision Weber number and the impact parameter, B, which measures the deviation of droplet colliding trajectory from that of an exact head-on collision. This study was further extended to alkane droplets by Jiang, Umemura & Law (Reference Jiang, Umemura and Law1992). They found that a low 
$We$ (
$O(1)$) results in droplet bouncing, while a moderately high 
$We$ (
$O(10)$ and above) results in separation with satellite droplet(s) formation (referred to as the reflex separation regime). It was established that, for the droplets to coalesce, the air layer between them must be drained out and the interfaces of droplets must come very close to each other (e.g. of the order of 100 nm), while bouncing happens because the air layer's pressure rises to a point where the droplets are pushed away before the air can drain out entirely. Subsequently, the well-known collision nomogram was proposed in the 
$We-B$ parameter space by Qian & Law (Reference Qian and Law1997). They observed stretched droplet separation at high B values and argued that the different collision behaviour between water and hydrocarbon droplets was due to the differences in their viscosity and surface tensions. They also found that increasing ambient pressure promotes droplet bouncing. More complex behaviours have also been observed for droplet collisions at high 
$We$ numbers or with droplets of different liquids. For example, droplet collision at 
$We > 200$ can result in droplet splattering (Pan et al. Reference Pan, Chou and Tseng2009), and collisions with dissimilar liquid droplets favoured a reflex separation tendency over stretched droplet separation, as compared with droplets of the same liquid (Chen Reference Chen2007).
Droplet size disparity (i.e. droplets with different diameters) adds substantial complexity to droplet collision processes and its effect has not been explored much. Tang et al. (Reference Tang, Zhang and Law2012) studied the dynamics during head-on collision of unequal-size droplets. They proposed a collision regimes diagram in a 
$We - \varDelta$ parameter space, where 
$We = \rho U^{2}R_{f}/ \gamma$ for unequal-size droplets (
$R_{f}$ is the radius of the smaller droplet) and 
$\varDelta$ is defined as the ratio of the large droplet diameter to the small droplet diameter. Their experiment was limited to 
$\varDelta < 3$. They corroborated that droplet separation is suppressed in unequal-size droplet collisions due to enhanced viscous dissipation through the internal motion that helps to stabilize the coalesced droplet. More general collision scenarios that included the oblique collision of unequal-sized liquid droplets were investigated by Chaitanya, Sahu & Biswas (Reference Chaitanya, Sahu and Biswas2021). They reported that the asymmetric flow, caused by both the unequal sizes of droplets and the oblique collision, affects the collision outcomes. The aforementioned studies extensively investigated binary droplet collision, without the presence of a solid surface.
A few recent studies have investigated the dynamics of a falling droplet impacting a sessile droplet, providing some useful insights into the droplet-to-droplet impact on solid surfaces. Impact of a sessile droplet on a hydrophobic substrate was studied by Wang, Feng & Zhao (Reference Wang, Feng and Zhao2008), where they observed four different regimes for the impact outcome: coalescence, complete rebounding, coalescence accompanied by conglutination and partial rebound after conglutination. A similar study on a hydrophilic surface was performed by Wakefield, Tilger & Oehlschlaeger (Reference Wakefield, Tilger and Oehlschlaeger2016). Following this, drop-on-drop impacts on superhydrophobic surfaces were experimentally examined by Damak & Varanasi (Reference Damak and Varanasi2018). They observed that a droplet impacting a sessile droplet on a solid surface (after coalescence) behaved similarly to a single droplet hitting the solid surface during the receding phase. Ramírez-Soto et al. (Reference Ramírez-Soto, Sanjay, Lohse, Pham and Vollmer2020) used both experiments and DNS simulations to study the collision dynamics of an oil droplet impacting a sessile droplet sitting on a superamphiphobic surface. Thereafter, Abouelsoud & Bai (Reference Abouelsoud and Bai2021) investigated the effect of the wettability on the impact outcomes by using different solid surfaces such as smooth glass, aluminum, copper, Teflon and coated glass. They observed complete coalescence for droplet impact on a hydrophobic surface (e.g. coated glass that has a contact angle of 
$119^{\circ }$), whereas both coalescence and bouncing were seen for hydrophilic surfaces that depend on Weber number. Very recently, Ray et al. (Reference Ray, Chi, Zhang and Cheng2023a) experimentally studied the collision of unequal-size droplets on a hydrophilic surface, and observed droplet coalescence, bouncing and partial coalesences after the impact. The experiments were performed at a large size ratio (i.e. 
$\varDelta = O(10)$), as compared with the experiments (i.e. 
$\varDelta = 1.4$) from Abouelsoud & Bai (Reference Abouelsoud and Bai2021). Although the droplet size ratio adopted by Abouelsoud & Bai (Reference Abouelsoud and Bai2021) and Ray et al. (Reference Ray, Chi, Zhang and Cheng2023a) differs considerably, the experiments seemingly corroborated a universal behaviour during head-on bouncing at different size ratios on wetting surfaces, where the restitution coefficient (defined as the ratio of bouncing velocity to impacting velocity of the falling droplet, and denoted as 
$\varepsilon$) is insensitive to the droplet size ratio at various Weber numbers (e.g. from 2–4.6). However, the comparison was made at only two greatly different size ratios without including the intermediate size ratios.
Therefore, we aim to investigate the head-on bouncing of two droplets on a wetting surface with a wider range of droplet size ratios using direct numerical simulation (DNS). The specific objectives of this study are: (1) to analyse the head-on bouncing dynamics of a falling droplet with a droplet on a hydrophilic surface with varying size ratios, (2) to determine the dependency of the coefficient of restitution on droplet size ratio and whether there exists a universal behaviour, and (3) to prove the universal behaviour if it exits, by using scaling analysis.
The rest of this paper is organized as follows. Section 2 summarizes the numerical methodology adopted in this study. Section 3 discusses the numerical and scaling analysis results, as well as validation against experiments. This is followed by a summary of this study in § 4.
2. Numerical methodology and specifications
2.1. Numerical method
The flow is considered incompressible and the fluids are Newtonian. The continuity and Navier–Stokes equations are expressed as
\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}=0, \end{gather}
\begin{gather}\rho\left[ \frac{\partial \boldsymbol{u}}{\partial t} + \boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u} \right]={-}\boldsymbol{\nabla} p+\boldsymbol{\nabla} \boldsymbol{\cdot} (2\mu \boldsymbol{D}) + \gamma\kappa \boldsymbol{n} {\delta}_{s}, \end{gather}
where 
$\boldsymbol {u}$ is the velocity vector, 
$\rho$ is the density, 
$p$ is the pressure, 
$\mu$ is the dynamic viscosity and 
$\boldsymbol {D}$ is the deformation tensor defined as 
$(\boldsymbol {\nabla } \boldsymbol {u} + \boldsymbol {\nabla } \boldsymbol {u}^{T} )/2$. The last term on the right-hand side of (2.2) is the surface tension term that is modelled using the continuum-surface-force model developed by Brackbill, Kothe & Zemach (Reference Brackbill, Kothe and Zemach1992), where 
$\delta _{s}$ is the Dirac delta function, 
$\gamma$ is the surface tension coefficient, 
$\kappa$ is the local curvature and the unit vector 
$\boldsymbol {n}$ is normal to the local interface. Equations (2.1) and (2.2) are solved by using the standard fractional-step projection method.
A conventional geometric volume of fluid (VOF) with height-function curvature estimation is employed to track the interface between the liquid and gas phases. To resolve the interface, a three-dimensional (3-D) octree adaptive mesh refinement is used. The cells in the region of interest including the droplets and region near the droplets are refined to maximum refinement level (i.e. 
$2^{12}$ cells) using a wavelet-based error method (Van Hooft et al. Reference Van Hooft, Popinet, Van Heerwaarden, Van der Linden, De Roode and Van de Wiel2018). This region is specified by a cylindrical region for each simulation (e.g. for the case of 
$\varDelta = 10.0$, the radius of the cylindrical region is 
$R_{s}$ plus 
$2R_{f}$, and the height is taken as the height of the sessile droplet plus 
$3R_{f}$). Outside this region the cells are not refined and kept at a minimum level of refinement (i.e. 
$2^{8}$ cells). A two-volume fractions method is used to track the interface of each droplet. The volume fraction 
$c_{i}$ (
$i = 1, 2$) is governed by the advection equation, which is expressed as
\begin{equation} \frac{\partial c_{i}}{\partial x}+\boldsymbol{\nabla} \boldsymbol{\cdot} ( c_{i} \boldsymbol{u} ) = 0, \end{equation}
where 
$c_{i} = 0$ for the gas phase, 
$c_{i} = 1$ for the liquid phase and 
$0 < c_{i} < 1$ for the gas–liquid interface. The density and dynamic viscosity are calculated using the volume fractions as 
$\rho = ( c_{1} + c_{2} )\rho _{l} + ( 1-c_{1} - c_{2} )\rho _{g}$ and 
$\mu = ( c_{1} + c_{2} )\mu _{l} + ( 1-c_{1} - c_{2} )\mu _{g}$, where subscripts 
$l$ and 
$g$ denote liquid and gas phases, respectively.
The mesh resolution close to the interface has a significant impact on the traditional VOF technique, and the accurate modelling of droplet bouncing necessitates an extremely fine-grained mesh at a significant computational expense. The two-volume fractions method independently describes the two droplets with no explicit conditions imposed for avoiding the overlap of interfaces. The maximum excess of the combined volume fractions, defined as the maximum of 
$(c_{1} + c_{2} -1)$, was used to verify that the interfaces belonging to the two droplets cannot contact and that there is no overlap between the different VOF functions. It is noted that this method would always enforce droplet bouncing for any droplet collision process and, thus, is applied in the present study for droplet bouncing cases that have been verified by experiments reported in previous studies (Coyajee & Boersma Reference Coyajee and Boersma2009; He, Xia & Zhang Reference He, Xia and Zhang2019).
The DNS of droplet collisions on a hydrophilic surface are performed using the open-source software Basilisk C (Popinet Reference Popinet2013). The free contact lines (FCL) condition has been used on solid surfaces to characterize the wettability of the surface by previous studies (Sakakeeny & Ling Reference Sakakeeny and Ling2020) and the same approach is used in this work. The FCL condition represents the motions of the contact line with a fixed contact angle. The height-function method developed by Afkhami & Bussmann (Reference Afkhami and Bussmann2008) is employed to specify the contact angle.
To accurately model the air film between two droplets, an extremely fine grid size is necessary due to the small size (around 100 nm) of the air film (Jiang et al. Reference Jiang, Umemura and Law1992). This significantly increases computational costs, and the effects of rarefied gas dynamics become significant. To overcome these difficulties, a 2-VOF approach has been proposed and used by several researchers (Ramírez-Soto et al. Reference Ramírez-Soto, Sanjay, Lohse, Pham and Vollmer2020; Alventosa, Cimpeanu & Harris Reference Alventosa, Cimpeanu and Harris2023) to study droplet bouncing. By avoiding overlap between the two-volume fraction fields, a situation similar to an incompletely drained air film with a large air film gap can be simulated, ensuring the occurrence of bouncing (see Appendix). It is assumed that the cases studied in this paper belong to a bouncing regime.
2.2. Numerical specifications
In the simulation all variables are made non-dimensional. The governing equations are non-dimensionalized using the radius of the falling droplet (
$R_{f}$), the capillary time scale (
$t_{c}=( \rho _{l}R_{f}^{3}/\gamma )^{1/2}$) and the capillary velocity (
$V_{\gamma }=R_{f}/t_{c}=( \gamma /\rho _{l}R_{f}^{3} )^{1/2}$). Then the dimensionless impact velocity can be expressed as 
$U/V_{\gamma }=( We )^{1/2}$. The relevant dimensionless numbers used in this work are the Weber number 
$We = \rho U^{2} R_{f}/\gamma$, the Ohnesorge number 
$Oh = \mu /( \rho \gamma R_{f} )^{1/2}$, the Bond number 
$Bo = \rho g R_{f}^{2}/\gamma$, the Reynolds number 
$Re = \rho U R_{f}/\mu$ and size ratio 
$\varDelta = R_{s}/R_{f}$. It should be noted that the spherical cap radius of the sessile droplets has been employed to define 
$\varDelta$. The present problem for droplets bouncing on a solid surface is controlled by three parameters (
$We$, 
$\varDelta$ and 
$\theta$) with fixed gas–liquid density, viscosity ratios and ambient pressure. The falling droplet diameter is also fixed at 
$300\ {\rm \mu}{\rm m}$, based on the previously reported experiments from the authors’ group (Ray et al. Reference Ray, Chi, Zhang and Cheng2023a). The 3-D computational domain is illustrated in figure 1. All the simulations performed in this work are three dimensional. The domain is a cube with an edge length 
$L = 8D_{s}$. Thus, the computation domain is big enough to avoid the influence of boundaries. The bottom face is a solid wall with no-slip boundary conditions. A fixed value of the contact angle with the no-slip condition is specified at the wall. Free outflow boundary conditions are specified at all other surfaces. A sessile droplet with the diameter 
$D_{s}$ and a falling droplet with diameter 
$D_{f}$ and velocity 
$U$ in the downward direction are initialized using the two-volume fractions method. The finest grid points in the computational domain correspond to 
$2^{12}$ uniform grids in each direction. Thus, the minimum cell size is 0.02 times the initial falling droplet radius, which in terms of dimensional length is equal to 3 
$\mathrm {\mu }$m. The number of grid points inside the falling droplet corresponds to 100 uniformly spaced grid points. In the present study, water droplet collisions at one atmospheric pressure are considered, again, same as the experiments (Ray et al. Reference Ray, Chi, Zhang and Cheng2023a). The time is non-dimensionalized with the inertial-capillary time scale (i.e. 
$t_{c}=( \rho _{l}R_{f}^{3}/\gamma )^{1/2}$, where 
$R_{f}$ is the radius of the falling droplet). By setting the CFL (Courant–Friedrichs–Lewy) number to 0.5, the time steps in the simulations are determined, resulting in a dimensionless time step of 
$4 \times 10^{-4}$. In a typical simulation, there are 
$8 \times 10^{5}$ grid points in the entire domain, and it takes about 200 h of physical time to perform each simulation on a single processor of 3.2 GHz kernel with one core and 64 GB memory.
Figure 1. The computational domain used for 3-D simulation of droplets collision on a hydrophilic surface.
2.3. Energy budget calculations
To get further insights into the dynamics of the bouncing droplets on a hydrophilic surface, the interchange of energy among different forms is investigated. The energy budget followed in this paper is similar to the one employed by Wildeman et al. (Reference Wildeman, Visser, Sun and Lohse2016). The total energy is the sum of the mechanical energy (
$E_{m}$), surface energy (
$E_{\gamma }$) and viscous dissipation (
$E_{\mu }$) and is expressed as
\begin{equation} E = E_{m}+E_{\gamma}+E_{\mu}, \end{equation}
where 
$E_{m} = E_{k} + E_{g}$, 
$E_{k}$ is the kinetic energy and 
$E_{g}$ is the gravitational potential energy. Here 
$E_{k}$ includes the translational kinetic energy and the oscillation and rotational energies. The energy terms can be found as
\begin{gather} E_{k}=\int_{\varOmega}\frac{1}{2} \rho( \boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{u} ) \,{\rm d}\varOmega, \end{gather}
\begin{gather}E_{\gamma}=\int_{A} \gamma\, {\rm d}A + \int_{A'} \gamma ( 1 - \cos\theta )\,{\rm d}A', \end{gather}
\begin{gather}E_{\mu}=\int_{0}^{t} \left( \int_{\varOmega} \phi \,{\rm d}\varOmega \right) {\rm d}t, \end{gather}
\begin{gather}E_{g}=\int_{0}^{t} \rho \boldsymbol{g} z \,{\rm d}t, \end{gather}
where 
$\boldsymbol {g}$, 
${\rm d}\varOmega$, 
${\rm d}A$ and 
${\rm d}A'$ are acceleration due to gravity, the differential volume, surface area at the liquid–gas interface and surface area at the liquid–solid surface, respectively, and 
$\phi$ is the viscous dissipation rate, expressed as
\begin{align} \phi& = 2 \mu\left[ \left( \frac{\partial u}{\partial x}\right)^{2} + \left( \frac{\partial v}{\partial y}\right)^{2} + \left( \frac{\partial w}{\partial z}\right)^{2} \right]\nonumber\\ &\quad + \mu \left[ \left( \frac{\partial u}{\partial y} + \frac{\partial v}{\partial x}\right)^{2} + \left( \frac{\partial v}{\partial z} + \frac{\partial w}{\partial y}\right)^{2} + \left( \frac{\partial w}{\partial x} + \frac{\partial u}{\partial z}\right)^{2} \right], \end{align}
where 
$u$, 
$v$ and 
$w$ are the components of velocity in the 
$X$, 
$Y$ and 
$Z$ coordinate directions, respectively.
3. Results and discussion
3.1. Comparison with experiments
To validate the present numerical set-up, the head-on collision of droplets on a solid surface at different size ratios is simulated and the simulation results are compared against the experimental results from Ramírez-Soto et al. (Reference Ramírez-Soto, Sanjay, Lohse, Pham and Vollmer2020) and Ray et al. (Reference Ray, Chi, Zhang and Cheng2023a). Two cases with greatly different size ratios and surface wettability are considered, in order to demonstrate the robustness of the DNS model.
In the first case we considered the collision of low-surface tension droplets with a small size ratio on a low-wetting surface. The experiment was performed by Ramírez-Soto et al. (Reference Ramírez-Soto, Sanjay, Lohse, Pham and Vollmer2020) who studied the impact of an oil droplet on a sessile oil droplet with equal volume on a low-wetting surface with a contact angle (
$\theta$) of 164.0
$^\circ$. The initial diameter of the falling droplet is about 2.05 mm, with a 
$We$ number of 1.30. The Bond number (
$Bo$) and Ohnesorge number (
$Oh$) for this case are 0.308 and 0.0216, respectively. The impact process is illustrated in eight different instances during the collision, as shown in figure 2(a) where the moment of initial contact between the falling droplet and the sessile droplet is referred to as 
$t = 0$ ms. As can be seen from figure 2(a), upon the falling droplet impacting the sessile droplet, both droplets deform and spread radially. During this stage (e.g. 
$t = 0\unicode{x2013}5.5$ ms), the kinetic energy of the falling droplets is transferred to the surface energies of both deformed droplets. After attaining the maximum horizontal diameter (e.g. 
$t = 8.0$ ms), the flow inside the droplets reverses and the receding phase begins. During this process, the sessile droplet transfers energy back to the falling droplet, increasing the kinetic energy of the falling droplet. Finally (e.g. 
$t = 25.0$ ms), the falling droplet bounces off while the sessile droplet stays on the solid surface and oscillates. Eventually, the oscillating sessile droplet loses the excess surface energy and attains the equilibrium shape (this stage is not shown in figure 2a). During the collision, part of the initial impact energy is converted to the surface energy of the sessile droplet and another part is lost as viscous dissipation energy inside the droplets and in the surrounding. It is clear from figure 2(a) that these characteristics and the corresponding timings are excellently reproduced by the simulation results. There is a small shift in the trajectory of the falling droplet at 
$t = 25.0$ ms. This shift occurs as the collision is not perfectly head-on in the experiments (i.e. the impact parameter 
$B$ is small 
$O( 10^{-2})$ but not zero). The restitution coefficient and bouncing time calculated from the experiment are 0.615 and 4.40, respectively, while their corresponding values calculated from DNS are 0.629 and 4.35, respectively. Small discrepancies in the restitution coefficient between experimental and numerical values could arise from the presence of surface contamination on the droplets. Past research has demonstrated that even minute amounts of contaminants can influence the dynamics of air film drainage (Klaseboer et al. Reference Klaseboer, Chevaillier, Gourdon and Masbernat2000) and impact interfacial mobility (Vakarelski, Yang & Thoroddsen Reference Vakarelski, Yang and Thoroddsen2022), ultimately affecting the bouncing behaviour.
Figure 2. Snapshots showing the collision of droplets on solid surfaces at different size ratios. (a) Bouncing of an oil droplet impacting a sessile oil droplet of equal volume on a soot-templated glass slide 
$(We, \varDelta, \theta ) = (1.3, 1.08, 164.0)$. (b) Bouncing for water droplets with a greatly different size ratio on a polyvinyl chloride surface 
$(We, \varDelta, \theta ) = (4.6, 7.2, 72.6)$. The images are taken from experimental videos from Ramírez-Soto et al. (Reference Ramírez-Soto, Sanjay, Lohse, Pham and Vollmer2020) and Ray et al. (Reference Ray, Chi, Zhang and Cheng2023a), and the simulation results are shown as lines. The falling droplet shape is demonstrated by an orange colour line while that of the sessile droplet by a light blue colour line. Plots (c,d) show the trajectories of the falling droplet (subpanel (i)) and the sessile droplet (subpanel (ii)) for the case corresponding to (a) and (b), respectively. The dashed line shows experimental droplet trajectories.
In the second case a falling water droplet with a diameter of 0.32 mm impacting a sessile water droplet with a much greater diameter (i.e. 2.29 mm) on a polyvinyl chloride surface in ambient air is considered. The experiment is taken from a previous study from the authors (Ray et al. Reference Ray, Chi, Zhang and Cheng2023a). The dimensionless numbers were 
$We = 4.6$, 
$Oh = 0.0092$, 
$Bo = 0.0034$. The density ratio and viscosity ratio between the surrounding and the droplets were 
$1\times 10^{-3}$ and 
$3\times 10^{-3}$, respectively. The modelling results and the experimental data are summarized in figure 2(b), where the modelling results are overlayed on the experimental images at eight different timings during the collision process, similar to the layout in figure 2(a).
Immediately seen from figure 2(b) is the correspondence between the modelling and experimental results, particularly with respect to the deformation of the falling droplet and the capillary wave near the collision point. It is also obvious from DNS results that a cavity is formed in the sessile droplet near the contact region upon the impact of the falling droplet, which first grows, then reaches a maximum depth and diminishes with time afterwards. The restitution coefficient and the bouncing time (
$t_{b}$) (defined as the time duration from the beginning of contact between the falling droplet and the sessile droplet to the first separation of the droplets) are 0.331 and 1.025 ms, respectively, in the experiment and 0.334 and 1.050 ms, respectively, from the simulation, which agree well with each other.
A direct comparison between the trajectory derived from the simulation results and the experimental data are shown in figure 2(c,d), which corresponds to figure 2(a,b), respectively. The simulated trajectories of the top, centre of mass and bottom positions for both the falling droplet and the sessile droplet are reported. The distances are measured from the solid surface to the top, centre and bottom of the droplets along the 
$z$ axis (vertically upward direction). However, a portion of the falling droplet's bottom position trajectory is not shown due to its entry into the cavity within the sessile droplet. The experimental images do not provide clear visibility of the droplets interface in this region. With comparisons in figure 2 for the two extreme scenarios, it has been demonstrated that the DNS model can adequately capture the dynamics during head-on bouncing of unequal-size droplets at different size ratios, from which the restitution coefficient and bouncing time of droplets can be accurately determined.
3.2. Head-on bouncing dynamics of unequal-size droplets
To get insight into the effect of size ratio on bouncing behaviour, we performed DNS for a wide range of size ratios. It is well known that bouncing is observed for a certain range of 
$We$ for equal-size droplet collision in air (Qian & Law Reference Qian and Law1997). In this study, the ranges of parameters leading to bouncing are selected based on experimental data from studies by Ramírez-Soto et al. (Reference Ramírez-Soto, Sanjay, Lohse, Pham and Vollmer2020), Abouelsoud & Bai (Reference Abouelsoud and Bai2021) and Ray et al. (Reference Ray, Chi, Zhang and Cheng2023a), as well as observations from similar phenomena reported by Qian & Law (Reference Qian and Law1997) and Alventosa et al. (Reference Alventosa, Cimpeanu and Harris2023). The ranges of parameters used in DNS are summarized in table 1.
Table 1. Value of relevant parameter ranges used in this work and their comparison with the previous experiments. Note that the experimental conditions (reported by Ramírez-Soto et al. Reference Ramírez-Soto, Sanjay, Lohse, Pham and Vollmer2020; Ray et al. Reference Ray, Chi, Zhang and Cheng2023a) specified here are corresponding to bouncing cases.
Simulations are performed for a wide range of 
$We$ and size ratios so that universal behaviour can be captured. The restitution coefficient (
$\varepsilon$) is calculated for wider 
$We-\varDelta$ parameter space. The size ratio 
$\varDelta$ is varied from 1.4 to 10.0 while 
$We$ is varied in the range of 2.0 to 4.5. The fluid considered in the present study is water and the falling droplet diameter is chosen as 0.30 mm. The contact angle at the solid surface can affect the bouncing behaviour that can further sophisticate the problem and, therefore, it is kept constant at 
$75^{\circ }$. These conditions are similar to the conditions in which the experiments were performed by Ray et al. (Reference Ray, Chi, Zhang and Cheng2023a) and showed bouncing. The value of density, surface tension and dynamic viscosity of water are 998 kg m
$^{-3}$, 72.9 mN m
$^{-1}$ and 
$1.0\times 10^{-3}$ Pa s, respectively, taken at 
$22\,^{\circ }{\rm C}$. Thus, the calculated values of 
$Oh = 1.0\times 10^{-2}$ and 
$Bo = 3.0\times 10^{-3}$, which are kept fixed.
Figure 3(a) (surface plot) shows the variation of 
$\varepsilon$ with 
$\varDelta$ and 
$We$. It can be seen that the values of 
$\varepsilon$ decline first and then attain almost constant value with increasing size ratio 
$\varDelta$ at a fixed 
$We$. This trend is similar to other 
$We$ considered here. The values of 
$\varepsilon$ are less than 0.7 and, thus, the collisions are inelastic. This inelastic behaviour indicates that a part of the impact energy of the falling droplet is lost during the bouncing process. The energy exchange during the bouncing process is analysed in detail in the subsequent discussion. Also, it is interesting to note that the value of 
$\varepsilon$ at a large size ratio is constant in a wide parameter space and this behaviour is explored further. The restitution coefficient as a function of 
$We$ and 
$\varDelta$ is plotted in figure 3(b). It can be observed that 
$\varepsilon$ remains insensitive to 
$We$ when 
$\varDelta$ is large (indicated by the dark red colour). The comparison between the simulation and experimental values demonstrates a satisfactory level of quantitative agreement.
Figure 3. (a) Surface plot showing the variation of 
$\varepsilon$ with 
$\varDelta$ and 
$We$. (b) Two-dimensional plot showing 
$\varepsilon$ as a function of 
$\varDelta$ and 
$We$. The numerical results are shown as a coloured surface and coloured lines in the plot, while the experimental data (obtained from Ray et al. Reference Ray, Chi, Zhang and Cheng2023a and Abouelsoud & Bai Reference Abouelsoud and Bai2021) are indicated by markers.
To understand the mechanism governing the constant behaviour of 
$\varepsilon$ in head-on bouncing of unequal-size droplets, further analysis is performed at a specific 
$We$, i.e. 
$We = 4.0$. The variation of 
$\varepsilon$ as a function of 
$\varDelta$ at 
$We = 4.0$ is first shown in figure 4(a). It can be seen that 
$\varepsilon$ decreases with increasing 
$\varDelta$ at a lower value of 
$\varDelta$ (e.g. 
$\varDelta < 3$), while becoming insensitive to 
$\varDelta$ with further increase in 
$\varDelta$. The restitution coefficient 
$\varepsilon$ reaches an asymptotic value of 0.293 for large 
$\varDelta$, which is identical to the measured value of 
$\varepsilon$ for droplet impact with a pool under similar conditions as reported by Alventosa et al. (Reference Alventosa, Cimpeanu and Harris2023). Figure 4(b) shows the variation of the dimensionless bouncing time (
$t_{b}$) with respect to 
$\varDelta$. It can be seen that the bouncing time is constant and is independent of 
$\varDelta$ for 
$\varDelta > 7.0$, where bouncing time fluctuates due to a different characteristic response time of the falling droplet and the sessile droplet that is proportional to their natural oscillation time, 
$t_{osc} = ( \rho R^{3}/\gamma )^{1/2}$ (where 
$R$ is the radius), at lower 
$\varDelta$ (
$< 7$). Based on the trends observed in figure 4, the entire 
$\varDelta$ space can be divided into three regimes. Regime I is specified by 
$\varDelta < 3$ (
$\varDelta$ close to 1), where 
$\varepsilon$ is high and changes with 
$\varDelta$. In this regime the falling droplet can recover a larger portion of its initial kinetic energy and, thus, a high 
$\epsilon$ value. Both the droplets have similar mass and oscillate similarly in the first mode (i.e. expansion and compression). As 
$\varDelta$ increases further (e.g. 
$3 < \varDelta <7$), regime II is entered where 
$\varepsilon$ remains constant but the bouncing time varies. The sessile droplet size is bigger than the falling droplet, with the mass ratio of the sessile droplet to the falling droplet at approximately 
$O (10 )$. Resonance can occur in this regime that reduces the bouncing time significantly. For 
$\varDelta > 7$, regime III is present where both 
$\varepsilon$ and bouncing time are insensitive to 
$\varDelta$. Within this regime, the sessile droplet mass is 
$O( 10^{3} )$ times heavier than the falling droplet mass, and capillary wave propagation on the sessile droplet surface becomes dominant.
Figure 4. (a) Coefficient of restitution, and (b) dimensionless bouncing time for a water droplet bouncing from a sessile water droplet on a hydrophilic surface as a function of droplet size ratio, 
$\varDelta$, (
$We = 4.0$, 
$Oh = 1.0 \times 10^{-2}$, 
$Bo = 3.0 \times 10^{-3}$, 
$\theta = 75^\circ$). Horizontal dashed line in subplot (a) corresponds to 
$\varepsilon = 0.293$ and in subplot (b) corresponds to 
$t/t_c = 4.6$. In (b), numerical data are indicated by circular filled markers and experimental data (Ray et al. Reference Ray, Chi, Zhang and Cheng2023a) are denoted as star markers.
Numerical simulations can provide detailed information on the velocity field and energy exchange between the droplets that are further analysed. Figure 5 shows the viscous dissipation rate, the velocity field and the contour plot of velocity magnitude at two different size ratios, namely 
$\varDelta = 2.0$ and 10.0. In figure 5(a) it can be seen that upon impact, the falling droplet induces flow inside the sessile droplet during its downward motion. Also, both droplets are compressed and an increase in the overall surface area of the droplets is noticed. A cavity is seen at the apex of the falling droplet (e.g. 
$t/t_{c} = 0.9$), which is quite dissimilar from a pancake shape (or flat D shape) observed during the impact of a droplet onto a solid surface (Clanet et al. Reference Clanet, Béguin, Richard and Quéré2004). This cavity creates a low-velocity field nearby as indicated by the light blue colour velocity contour. The low velocity creates a low-pressure region, which pulls the liquid in the falling droplet in the upward direction. Thus, large velocity vectors are present in the air near the apex region at 
$t/t_{c} = 1.25$. Although the velocity magnitudes are larger in the region, the kinetic energy of the air is still small due to the low density of air, which is 
$O (10^{-3} )$ times smaller than the droplet density. Also, the viscous dissipation in the air is low as the dynamic viscosity of air is much lower than that of water. It is interesting to note that in this instant (i.e. 
$t/t_{c} = 1.25$), the flow inside the falling droplet near the apex region is upward whereas the flow near the bottom region is still downward. This opposing motion causes the falling droplet to expand along the 
$Z$ axis. At 
$t/t_{c} = 2.40$, the flow in the falling droplet is completely reversed and the velocity vector is upward in trend. As the flow velocity of the bottom part of the falling droplet is greater than the top part of the falling droplet, the bottom of the falling droplet is forced to expand to both sides, as seen at 
$t/t_{c} = 3.0$. At the same time, the sessile droplet retracts, and the water flows from both sides to the centre of the sessile droplet, forcing the sessile droplet to stretch upward from its apex. At 
$t/t_{c} > 4.4$, the upward flow in the sessile droplet is suppressed and further reversed, while the falling droplet keeps travelling upward, eventually leading to the separation of the two droplets (e.g. 
$t/t_{c}= 5.0$). The viscous dissipation rate is prominent near the impact location and near to both droplet surfaces. It is obvious that considerable deformations exist in the sessile droplet at bouncing.
Figure 5. Snapshots of DNS with viscous dissipation rate, velocity vectors and contours of velocity magnitude shown: (a) 
$\varDelta = 2.0$ and (b) 
$\varDelta = 10.0$ (
$We = 4.0$, 
$Oh = 1.0 \times 10^{-2}$, 
$Bo = 3.0 \times 10^{-3}$, 
$\theta = 75^\circ$).
When the droplet size ratio increases from 2.0 to 10.0, as seen in figure 5(b), the overall dynamics and flow patterns are similar to those observed in figure 5(a). Nevertheless, noticeable differences can be observed between the two cases. Specifically, (1) the cavity formed near the apex of the falling droplet does not create a stronger pull, and thus, the velocity field (near the apex of the falling droplet) is not too large (i.e. 
$t/t_{c} = 1.25$); (2) a cavity with larger depth is formed in the sessile droplet (i.e. 
$t/t_{c} = 1. 80$); (3) the deformation of the sessile droplet occurs at a localized region near the collision point; and (4) the shape of the sessile droplet surface near the apex just after the onset of bouncing (i.e. 
$t/t_{c} = 4.70$) is similar to that before impact (i.e. 
$t/t_{c} = 0.00$), whereas the shape is quite dissimilar in the case with a smaller droplet size ratio.
Figure 6 shows the variation of energies with respect to time for the falling and sessile droplets during the collision and bouncing, which are determined by energy budget calculations as described in § 2.3. Five different energies are determined (namely kinetic energy of the falling droplet, kinetic energy of the sessile droplet, change in the surface energy of the falling droplet, change in the surface energy of the sessile droplet and total viscous dissipation energy) at the two size ratios as adopted in figure 5. All these energy terms are non-dimensionalized with 
$\rho R_{f}^{3} V_{\gamma }^{2}$. For instance, the impact kinetic energy of the falling droplet (
$KE$) is 
$( 1/2 ) \rho U^{2} ( 4{\rm \pi} R_{f}^{3}/3 )$, which when divided by 
$\rho R_{f}^{3} V_{\gamma }^{2}$ yields the dimensionless kinetic energy 
$E_{k}=( 1/2 ) \rho U^{2} ( 4{\rm \pi} R_{f}^{3}/3 )/\rho R_{f}^{3} V_{\gamma }^{2} =( 2{\rm \pi} /3 ) U^{2}/V_{\gamma }^{2} = ( 2{\rm \pi} /3 ) U^{2}/( \gamma /\rho R_{f} ) =( 2{\rm \pi} /3 ) U^{2}\rho R_{f} U^{2}/\gamma = ( 2{\rm \pi} /3 ) We$ (here the definition of capillary velocity and the Weber number are used). Thus, for 
$We = 4.0$, 
$E_{k} = ( 2{\rm \pi} /3 ) \times 4\approx 8$ (as seen in figure 6). Let us consider head-on impact at 
$\varDelta = 2.0$ (figure 6a). First, the falling droplet impacts the sessile droplet and the falling droplet decelerates as its motion is opposed and it stretches radially. The kinetic energy of the falling droplet decreases with time and it is converted into the surface energies of both droplets. At 
$t/t_{c} = 0.9$, the falling droplet has minimum kinetic energy as flow reversal occurs, as can be seen from figure 5(a). Maximum deformation of the falling droplet occurs at this time instant (
$t/t_{c} = 0.9$). Afterward the kinetic energy of the falling droplet recovers gradually with some oscillation. The surface energy of the falling droplet increases with time and attains a maximum value at a maximum deformation time instant, as can be seen from figure 5(a). Thereafter, the surface energy oscillates between the maximum and minimum values due to the harmonic oscillation of the falling droplet. It can be noticed that the change in the surface energy of the falling droplet is quite small at the bouncing time. The kinetic energy of the sessile droplet reaches a maximum value at the maximum deformation of the falling droplet as the internal flow is the strongest at around that time, as can be seen from figure 5(a). It reduces to nearly zero at 
$t/t_{c} = 2.0$ as a result of maximum deformation of the sessile droplet, and reaches another peak at 
$t/t_{c} = 3.0$ due to retraction motion enhancing the internal flow, as shown in figure 5(a). The surface energy of the sessile droplet keeps increasing until 
$t/t_{c} = 2.0$, as this is the timing when it reaches maximum deformation, as illustrated in figure 5(a). It decreases during retraction up until 
$t/t_{c} = 3.2$, and increases again afterward as the droplet is being pulled upward. It is clear that a part of the energy from the falling droplet is transferred to the sessile droplet during the collision, either in the form of surface energy or kinetic energy (internal flow). The total viscous dissipation energy increases with time. From this energy budget calculation, it can be seen that the surface energy of the sessile droplet and the kinetic energy of the falling droplet are the dominant energy term at the bouncing time.
Figure 6. Variations of different energies with non-dimensional time. Subscript ‘
$s$’ denotes the sessile droplet and subscript ‘
$f$’ denotes the falling droplets. The dashed grey vertical line represents time of maximum deformation of the falling droplet (
$t_m$) and the solid vertical line denotes the time of bouncing (
$t_b$). Here TE denotes total energy.
The energy distribution before the maximum deformation of the falling droplet is very similar between the two droplet size ratios, as also can be seen from the deformation and internal flow patterns of the falling droplets in figure 5(a,b). However, the energy distribution and oscillation patterns are dramatically different between the large size ratio (figure 6b) and the small size ratio (figure 6a) at 
$t/t_{c} > 0.85$. Again, these correlate with figure 5(b), where capillary waves are generated on the sessile droplet and move on its surface and affect the internal flow motion. For a large size ratio, the deformation and flow inside the sessile droplet occur mostly in a region near the impact and the flow intensity inside is much weaker compared with the case of a small size ratio. Also, the kinetic energy and the change in the surface energies of the falling droplet are smaller at a large size ratio. It is further noticed that the energy distribution near bouncing remains almost unchanged. It can be noted that the gas-phase energies are negligible due to significantly lower values of the density and viscosity of air (
${\sim } O ( 10^{-3} )$ smaller) as compared with water. Also, the change in potential energies of both the failing droplet and the sessile droplet are negligible (<0.3 % of total energy) and, thus, neglected (Ray, Han & Cheng Reference Ray, Han and Cheng2023b).
Although the energy distribution with time provided great insight into the bouncing of unequal droplets, we are particularly interested in the energy distribution at the bouncing time that determines the values of the restitution coefficient. Thus, the variation of the energy distribution at the bouncing time with size ratio is further analysed. For the sake of simplicity and ease of comparison, all the energies are represented as a fraction of the initial impact energy. The important energy terms, that is, total viscous dissipation, and the surface and kinetic energies of the sessile droplet are only plotted as a function of 
$\varDelta$ at bouncing time, as shown in figure 7. The energies of the falling droplet are not shown as the change in its surface energy is negligible and its kinetic energy can be deduced by energy balance, as can be seen from figure 6. It can be seen that the ratio of total viscous dissipation energy to the impact energy, 
$E_{\mu }/E_{0}$, first decreases, then increases followed by further decline and then remains almost constant with increasing 
$\varDelta$. The fractional surface energy of the sessile droplet, 
$\Delta E_{S,s}/E_{0}$, shows an oscillating pattern followed by extremely gradual increments with rising 
$\varDelta$. The fractional kinetic energy of the sessile droplet first increases, then attains an almost constant value with enhancing 
$\varDelta$. The change in the surface energy of the falling droplet is negligible during bouncing, as indicate by near-zero values of 
$\Delta E_{S,f}/E_{0}$ for different size ratios. The kinetic energy of the falling droplet decreases and eventually reaches a nearly constant value as 
$\varDelta$ increases. These energy trends indicate that the fractional energies of the sessile droplet and viscous dissipation attained a nearly constant value at large size ratios.
Figure 7. (a) Variation of the ratio of viscous dissipation energy at bouncing to initial impact energy with 
$\varDelta$. (b) Variation of the ratio of the surface energy of the sessile droplet at bouncing to initial impact energy with 
$\varDelta$. (c) Variation of the ratio of the kinetic energy of the sessile droplet at bouncing to initial impact energy with 
$\varDelta$. (d) Variation of the fractional surface energy of the falling droplet with 
$\varDelta$. (e) Variation of the fractional kinetic energy of the falling droplet with 
$\varDelta$ (
$We = 4.0$).
From the foregoing discussion, it is apparent that energy distribution and oscillation are governed by the magnitude of deformation. As such, the deformation of the sessile droplet is further analysed quantitatively. In order to do so, we introduce the parameter of crater depth and wave amplitude, denoted as 
$\delta$ and 
$A$, respectively. The normalized crater depth and normalized wave amplitude are indicated as 
$\tilde {\delta }$ and 
$\tilde {A}$, respectively, both of which are normalized with respect to 
$R_{f}$. To find the crater depth (
$\delta$), the height of the sessile droplet (denoted by 
$h( t )$) is tracked with time, then it is subtracted from the initial (also undisturbed) height of the sessile (denoted by 
$H$) (see figure 8a). The evolution of wave amplitude with angular coordinate 
$\alpha$ (as defined in figure 8a) and time are shown in figure 8. Here, 
$\alpha = 0^\circ$ denotes the left-hand vertex whereas 
$\alpha = 180^\circ$ represents the right vertex, and 
$\alpha = 90^\circ$ is the apex of the sessile droplet.
Figure 8. Deformation of the sessile droplet. (a) Schematic diagram showing capillary waves on the sessile droplet. The dashed brown line represents an undeformed droplet whereas the solid black line shows a deformed sessile droplet during the bouncing process. (b) Variation of normalized wave amplitude with angle 
$\alpha$ at different time instants. (c) Evolution of normalized crater depth at the 
$Z$ axis with time.
The spatial variation of normalized capillary wave amplitude (
$\tilde {A}$) at different time instants for 
$\varDelta = 10.0$ is shown in figure 8(b). By definition, wave amplitude is zero everywhere at 
$t/t_{c} = 0$. As the falling droplet impacts the sessile droplet, a crater is formed at the apex (
$\alpha = 90^\circ$) that grows with time until it reaches the maximum depth (i.e. 
$t/t_{c} = 2.0$) and then decreases with time. Since the complete shape of the sessile droplet is well defined for a known crater depth at the axis of symmetry (
$Z$ axis), the temporal evolution of the crater depth at the 
$Z$ axis is plotted for two size ratios, 8.0 and 10.0, as illustrated in figure 8(c). It is seen that the cavity depth increases monotonically, reaches a maximum depth (
$\tilde {\delta }_{m}$, where subscript 
$m$ denotes maximum value), and eventually declines. Also, the maximum depth of the crater for 
$\varDelta = 10.0$ is slightly greater than for 
$\varDelta = 8.0$ (i.e. 
$t/t_{c} = 2.0$). The decline is qualitatively and quantitatively similar for both size ratios (i.e. 
$t/t_{c} > 2.0$).
The dimensionless maximum crater depth (
$\tilde {\delta }_{m}$) is plotted as a function of 
$\varDelta$ at a fixed 
$We$ and as a function of 
$We$ at a fixed 
$\varDelta$ in figure 9(a,b). The surface energy of the sessile droplet is related to 
$\tilde {\delta }_{m}$ (Jayaratne & Mason Reference Jayaratne and Mason1964). It is observed that the dimensionless maximum crater depth remains constant with increasing 
$\varDelta$ at a fixed 
$We$, while it decreases linearly with increasing 
$We$ at a fixed 
$\varDelta$. This information is utilized in the next section for scaling analysis.
Figure 9. (a) Variation of the dimensionless maximum crater depth with 
$\varDelta$ at a fixed 
$We$. (b) Variation of the dimensionless maximum crater depth with 
$We$ at a fixed 
$\varDelta$.
Figure 10(a) depicts the variation of normalized crater velocity with time. The crater velocity is defined as the rate of change of the crater depth with time (i.e. 
$u_{ct}= {{\rm d} \tilde {\delta }}/{{\rm d}\tilde {t}}$, where 
$\tilde {t} = t/t_{c}$), and it is obtained by numerically differentiating the crater depth over time. It can be observed that the crater depth velocity initially remains constant, then decreases to a minimum and subsequently sharply increases followed by a gradual increase. Notably, the initial normalized crater velocity is interestingly close to 0.5 (until 
$t/t_{c} = 1.0$). This velocity matches the velocity of the air layer between a water droplet and a pool of the same liquid during droplet impact, as reported by Tran et al. (Reference Tran, de Maleprade, Sun and Lohse2013). This suggests that at these size ratios, the capillary wave on the sessile droplet resembles the waves generated in a pool. In Tran et al. (Reference Tran, de Maleprade, Sun and Lohse2013) it is observed that a thin air film is present between the droplet and the pool, which persists even when the droplet enters inside the pool and creates a crater. This air film exists during the downward motion and break-up, when the crater nearly reaches its maximum depth. Similarly, in the case of droplets bouncing on a solid surface with large size ratios, as studied in this work, an air film exists between the droplets throughout the entire bouncing process. Considering the similarity between these two processes, a comparison is made between the crater velocities associated with them. Figure 10(b) shows the variation of the normalized initial crater velocity (
$u_{ct-0}/U$) with 
$\varDelta$ at a fixed 
$We$. It is observed that the normalized initial crater velocity remains almost constant with increasing 
$\varDelta$.
Figure 10. (a) Variation of the normalized crater velocity with time. (b) Variation of the normalized initial crater velocity with 
$\varDelta$.
In order to gain a more comprehensive understanding of the bouncing behaviour exhibited by different liquids, characterized by their respective 
$Oh$ numbers, further simulations are performed. The variation of the restitution coefficient with 
$Oh$ number, ranging from 
$10^{-5}$ to 
$0.5$, are depicted in figure 11(a) for two different size ratios, considering negligible gravity (
$Bo = 0$). At low 
$Oh$ values, referred to as the inertia regime, the restitution coefficient remains almost constant, with values similar to those shown in figure 3. As 
$Oh$ increases to intermediate values (i.e. 
$O ( 10^{-1} )$), the restitution coefficient decreases monotonically due to increased energy loss caused by viscosity. This region is denoted as the viscous regime. Snapshots of numerical simulations illustrating these three regimes are provided in figure 11(b). At 
$Oh = 10^{-5}$, both the falling droplet and the sessile droplet exhibit significant deformation, similar to the scenario depicted in figure 2(b). As 
$Oh$ increases from 
$10^{-5}$ to 
$0.1$, the deformation in both the droplets decreases. To further quantify the differences between these regimes, the trajectories of the falling droplet (at three locations – the top, centre of mass and bottom of each droplet) and dimensionless crater depth are plotted as a function of time in figure 11(c,d), respectively. It is observed that a shallow crater depth is formed at higher 
$Oh$ numbers. Moreover, in the viscous regime the retraction of the crater occurs slowly, requiring a longer time to reach an equilibrium position.
Figure 11. (a) The restitution coefficient as a function of the 
$Oh$ number. (b) Snapshot of simulations at different 
$Oh$ number (
$Oh = 10^{-5}$ for subpanel (i) and 0.1 for (ii)). (c) The trajectory of the falling droplet (subpanel (i)) and the dimensionless crater depth (subpanel (ii)) at different 
$Oh$ number.
3.3. Scaling analysis
In this subsection we use scaling arguments to rationalize the independence of the coefficient of restitution at large size ratios (in regime III). As revealed in the prior subsection, at a fixed 
$We$, the coefficient of restitution remains constant at large size ratios (figure 4a) and is determined by the energy distribution. From energy calculations as shown in figure 7, we already know that energy distribution near bouncing for a large size ratio does not change much with increasing size ratio. We propose the scaling for different energy terms and use the energy balance to determine the scaling of the restitution coefficient.
To begin with, we write the energy balance at the beginning of impact and the occurrence of bouncing:
\begin{equation} E_{0} = E_{K}+\Delta E_{S}+\Delta E_{g}+ E_{\mu}. \end{equation}
Here 
$E_{0}$ is the initial impact energy, 
$E_{K}$ is the kinetic energy, 
$\Delta E_{S}$ is the change in surface energy, 
$\Delta E_{g}$ is the change in potential energy and 
$E_{\mu }$ is the viscous dissipation until the bouncing time. The subscripts 
$0$, 
$K$, 
$S$, 
$g$ and 
$\mu$ stand for initial, kinetic, surface, gravitational and viscous, respectively. The initial impact energy is the kinetic energy of the falling droplet and is expressed as
\begin{equation} E_{0} = \tfrac{1}{2} {\rm \pi} \rho R_{f}^{3} U^{2}\sim \rho R_{f}^{3} U^{2}, \end{equation}
where 
$U$ is the falling droplet impact velocity at the beginning of collision. Each other term on the right-hand side of (3.1) has a contribution from both the falling droplet and the sessile droplet. Thus, 
$E_{K} = E_{K,f} + E_{K,s}$, where subscript 
$f$ and 
$s$ denote falling droplet and sessile droplet, respectively.
Let 
$V_{b}$ be the velocity of the falling droplet at bouncing time. Then the kinetic energy at bouncing (i.e. 
$t = t_{c}$) can be expressed as
\begin{equation} E_{0} = \tfrac{1}{2} {\rm \pi} \rho R_{f}^{3} V_{b}^{2}\sim \rho R_{f}^{3} V_{b}^{2}. \end{equation} To determine 
$E_{K,s}$, we assume that the falling droplet induced a spherical symmetric radial flow in the sessile droplet. A crater in the shape of a spherical cap of radius of 
$\delta _{m}$ is present on the top of the sessile droplet. Consider a thin liquid shell at radial distance 
$r$ and thickness 
$dr$ as shown in figure 12. The centre of this shell is at the apex of the sessile droplet at the beginning of the impact (denoted as point 
$O$ in figure 12). The velocity of this thin shell (
$V$) is obtained by balancing the momentum of the falling droplet with the momentum of liquid inside the shell in the sessile droplet, i.e. 
$\frac {4}{3} {\rm \pi} \rho R_{f}^{3} U \approx \frac {4}{3} {\rm \pi} \rho r^{3} V$. Thus, 
$V \approx U R_{f}^{3}/ r^{3}$. This approach is similar to the approach used by Tang et al. (Reference Tang, Saha, Law and Sun2018) for droplet bouncing from a pool. The key distinction is that while Tang et al. (Reference Tang, Saha, Law and Sun2018) uses an upper limit of 
$\infty$ (infinity), our study employs 
$R_{s}$ as the upper limit. To find the kinetic energy of the sessile droplet, the kinetic energy of the liquid shell is integrated from 
$\delta _{m}$ to 
$R_{s}$ and can be expressed as
\begin{align} E_{K,s} &= \int_{\delta_{m}}^{R_{s}} \rho( 2{\rm \pi} r^{2} ) V^{2} \,{\rm d}r = \int_{\delta_{m}}^{R_{s}} \rho( 2{\rm \pi} r^{2} ) \left( \frac{U R_{f}^{3}}{r^{3}} \right)^{2} {\rm d}r \nonumber\\ &= 2{\rm \pi} \rho U^{2} R_{f}^{6}\int_{\delta_{m}}^{R_{s}} \frac{{\rm d}r}{r^{4}} = 2{\rm \pi} \rho U^{2} R_{f}^{6}\left| \frac{r^{{-}4+1}}{-4+1} \right|_{\delta_{m}}^{R_{s}} \nonumber\\ & = 2{\rm \pi} \rho U^{2} R_{f}^{6}\left( \frac{1}{\delta_{m}^{3}}- \frac{1}{R_{s}^{3}} \right). \end{align}
By using 
$\delta _{m}\sim R_{f}$ (see figure 9a) and imposing the condition 
$R_{s} > R_{f}$ for a large size ratio, 
$1/\delta _{m}^{3} \gg 1/R_{s}^{3}$. Now, (3.4) can be written as
\begin{equation} E_{K,s}\approx \frac{2}{3}{\rm \pi} \rho U^{2} \frac{R_{f}^{6}}{\delta_{m}^{3}} \approx \frac{2}{3}{\rm \pi} \rho U^{2} \frac{R_{f}^{6}}{R_{f}^{3}} = \frac{2}{3}{\rm \pi} \rho U^{2} R_{f}^{3}. \end{equation}
It should be noted that the surface area of the liquid shell will change with increasing 
$r$, and thus, (
$2{\rm \pi} r^{2}$) in (3.4) can be replaced with (
$a2 {\rm \pi} r^{2}$), where 
$a$ is a constant such that 
$a < 1$. Equation (3.5) becomes
\begin{equation} E_{K,s}\approx \tfrac{2}{3}{\rm \pi} a \rho U^{2} R_{f}^{3} \sim \rho U^{2} R_{f}^{3}. \end{equation}
It should be noted that (3.6) provides the kinetic energy of the sessile droplet at bouncing time (i.e. 
$t =t_{c}$). Combining (3.3) and (3.6), the total kinetic energy at bouncing is expressed as
\begin{equation} E_{K} = E_{K,f} + E_{K,s} \sim \rho U^{2} R_{f}^{3} +\rho U^{2} R_{f}^{3} \sim \rho U^{2} R_{f}^{3}. \end{equation}Thus, the fractional kinetic energy is
\begin{equation} \frac{E_{K}}{E_{0}} \sim \rho U^{2} R_{f}^{3}/ \rho U^{2} R_{f}^{3} \sim O(1), \end{equation}
(3.6) agrees with the trend observed in figure 7(c). For low Bond numbers, 
$\Delta E_{g}$ can be neglected.
Figure 12. Schematic diagram showing the liquid shell inside the sessile droplet.
At a large size ratio, the change in surface energy of the sessile droplet is much greater than the falling droplet at bouncing, i.e. 
$\Delta E_{S,f} \ll \Delta E_{S,s}$ (in fact, 
$\Delta E_{S,f}$ is close to zero, see figure 7d), then 
$\Delta E_{S} = \Delta E_{S,f} + \Delta E_{S,s} \sim \Delta E_{S,s}$. Following the work of Jayaratne & Mason (Reference Jayaratne and Mason1964), the surface energy of the sessile droplet is considered as a proposal to the square of the maximum crater depth times surface tension, i.e.
\begin{equation} \Delta E_{S} \sim \Delta E_{S,s} = C{\rm \pi} \delta_{m}^{2}\gamma \sim \delta_{m}^{2}\gamma, \end{equation}
where 
$C$ is a constant and 
$\delta _{m}$ is the maximum depth of the crater formed at the sessile droplet surface. Therefore,
\begin{equation} \frac{\Delta E_{S,s}}{E_{0}} \sim \delta_{m}^{2}\gamma / \rho R_{f}^{3} U^{2} = \tilde{\delta}_{m}^{2}/We. \end{equation}
From (3.10), for a fixed 
$We$, we obtain that 
$\Delta E_{S,s}$ weakly depends on 
$\varDelta$ as 
$\delta _{m}$ slightly changes with 
$\varDelta$, which agrees with figure 9(a). When 
$We$ is decreased, then 
$\delta _{m}$ also decreases (see figure 9b) such that 
$\tilde{\delta}_{m}^{2}/We$ remains constant in (3.10). This is also observed from figure 7(b).
Due to viscosity, energy is lost inside the falling droplet, the sessile droplet and the surrounding air. Following Alventosa et al. (Reference Alventosa, Cimpeanu and Harris2023), the rate of viscous energy dissipation (per unit volume) inside a droplet can be estimated to scale as 
$\mu ( \boldsymbol {\nabla } \boldsymbol {u} )^{2} \sim \mu U^{2}/R^{2}$ (where 
$R$ is the droplet radius). For a characteristic fluid volume 
$R^{3}$, the viscous energy dissipation rate scales as 
$\mu U^{2}R$. As shown in the prior sections (figure 4b), the bouncing time 
$t_{b} \sim t_{c}$. Thus, the viscous dissipation inside the falling droplet can be expressed as
\begin{equation} E_{\mu, f} \sim \mu U^{2} R_{f} t_{b} = \mu U^{2} R_{f} \sqrt{\rho R_{f}^{3}/\gamma}. \end{equation}The fractional viscous energy loss inside the falling droplet is
\begin{equation} \frac{E_{\mu, f}}{E_{0}} \sim \frac{\mu U^{2} R_{f} \sqrt{\rho R_{f}^{3}/\gamma}}{\rho R_{f}^{3} U^{2}} = Oh. \end{equation}To find the viscous dissipation energy inside the sessile droplet, the thin liquid shell is considered again. By integrating the viscous dissipation in the liquid shell, the viscous dissipation energy inside the sessile droplet is expressed as
\begin{align} E_{\mu,s} &= \int_{0}^{t_{c}}\int_{\delta_{m}}^{R_{s}}\mu\left( \frac{V}{r} \right)^{2} ( a2{\rm \pi} r^{2} ) \,{\rm d}r\, {\rm d}t= \int_{0}^{t_{c}}\int_{\delta_{m}}^{R_{s}}\mu\left( \frac{U R_{f}^{3}}{r^{4}} \right)^{2} ( a2{\rm \pi} r^{2} ) \,{\rm d}r \,{\rm d}t \nonumber\\ & = \int_{0}^{t_{c}}\int_{\delta_{m}}^{R_{s}}( \mu a2{\rm \pi} U^{2} R_{f}^{6} ) \frac{{\rm d}r}{r^{6}} \,{\rm d}t \nonumber\\ & = \int_{0}^{t_{c}} \left( \frac{2}{5}{\rm \pi} \mu a U^{2} R_{f}^{6} \right)\left( \frac{1}{\delta_{m}^{5}} - \frac{1}{R_{s}^{5}} \right){\rm d}t. \end{align}
In (3.13) the variable 
$a$ (such that 
$a < 1$) represents the change in the surface area of a spherical liquid shell as the radius 
$r$ increases. Again, by using 
$\delta _{m} \sim R_{f}$ (see figure 9a) and imposing the condition 
$R_{s}> R_{f}$ for a large size ratio, we observe that 
$1/\delta _{m}^{5}\gg 1/R_{s}^{5}$. Now, (3.13) can be written as
\begin{align} E_{\mu,s} &= \int_0^{t_c}\left( \frac{2}{5} {\rm \pi} \mu a U^2 R_f^6 \right)\left( \frac{1}{\delta_m^5} \right){\rm d}t \approx \int_0^{t_c}\left( \frac{2}{5} {\rm \pi} \mu a U^2 R_f^6 \right)\left( \frac{1}{R_f^5} \right){\rm d}t \nonumber\\ & \approx \frac{2}{5} {\rm \pi} \mu a U^2 R_f t_c \sim \mu U^2 R_f \sqrt{\rho R_f^3 \gamma}. \end{align}The fraction of viscous dissipation for the sessile droplet can be written as
\begin{equation} \frac{E_{\mu,s}}{E_0} \sim \frac{\mu V_s^2 R_i \sqrt{\rho R_f^3 / \gamma}}{\rho R_f^3 U^2} \sim \frac{\mu U^2 R_f \sqrt{\rho R_f^3 / \gamma}}{\rho R_f^3 U^2} = Oh. \end{equation}As the dynamic viscosity of air is much small than water, the viscous dissipation loss in the surrounding air is neglected. Therefore, the total fraction of viscous energy loss is given as
\begin{equation} \frac{E_{\mu}}{E_0} = \frac{E_{\mu,f}+E_{\mu,s}}{E_0} \sim Oh. \end{equation}
For a fixed 
$Oh$, (3.16) indicates that 
$E_{\mu }/E_0$ is insensitive to the size ratio, which agrees with figure 7(a).
Now we replace each energy term in the energy balance, (3.1), with their corresponding scales. Thus, (3.1) can be rewritten after using (3.6), (3.10) and (3.16) as
\begin{equation} 1=\frac{E_K}{E_0}+ \frac{\Delta E_S}{E_0} + \frac{E_{\mu}}{E_0} \sim \frac{V_b^2}{U^2}+O(1) +\frac{\tilde{\delta}_m^2}{We} + Oh. \end{equation}
By definition, 
$V_b/U$ is the restitution coefficient 
$\varepsilon$. After rearranging (3.17), we get the restitution coefficient as
\begin{equation} \varepsilon \sim \left[ O(1) - \frac{\tilde{\delta}_m^2}{We} - Oh\right]^{1/2}. \end{equation}
From (3.18), it can be seen that for a fixed 
$We$ and 
$Oh$, 
$\varepsilon$ is very weak depending on 
$\varDelta$ through. This proves that the restitution coefficient is constant at large size ratios. Furthermore, it can be noted from (3.18) that 
$\varepsilon$ decreases nonlinearly with 
$Oh$ as observed in figure 11(a). This further validates the proposed scaling analysis.
4. Conclusion and outlook
In this paper we studied the head-on collision of two unequal-size droplets on a wetting surface to a wide range of size ratios using the DNS technique. The robustness of the numerical model is established by comparing the DNS of two extreme cases with the experimental data from previous studies. Direct numerical simulation accurately predicts the droplet's shape, the coefficient of restitution and bouncing time. Applying this validated model to more size ratios at different Weber numbers and Ohnesorge numbers, we observed a universal behaviour in head-on bouncing of unequal-size droplets on a wetting surface.
Bouncing is characterized by the restitution coefficient that is related to the fraction of energy lost during bouncing. It is observed that the restitution coefficient decreases and then attains a constant value with an increasing size ratio. In the inertial regime, the restitution coefficient is close to 0.3 and becomes insensitive to size ratios for a size ratio greater than 7.0. Moreover, the bouncing occurs at 4.55 times of inertial-capillary time scale (
$t_{c}$) for large size ratios.
In order to understand the governing mechanism, we looked at the detailed flow patterns and dynamics during collision at a representative Weber number. We observed that the flow and bouncing characteristics are dramatically different at small and large size ratios. To understand the energy distribution associated with the different collision dynamics and flow patterns, we performed energy budget calculations. Results correlate well with the flow patterns and collision dynamics, which indicate a more stable energy distribution at large size ratios near bouncing time, which is mostly contributed by the sessile droplet. Since energy distribution is governed by the deformation characteristics, we quantified the deformation of the sessile droplet by crater depth, and again observed consistency in crater depth velocity. The velocity ratio (crater depth velocity to impact velocity) is equal to 0.5 that is similar to the ratio reported in the bouncing of a droplet on a pool. At large size ratios, a crater is formed due to the impact of the falling droplet and generates capillary waves on the sessile droplet surface. These capillary waves store the energy and create a flow inside the sessile droplet and, thus, dictate the kinetic energy and viscous energy dissipation of the sessile droplet. The cause of the universal behaviour (constant restitution coefficient) is that above a certain size ratio, the deformation characteristics of the sessile droplet become nearly identical at different size ratios. A similar behaviour has been observed for a droplet bouncing on a pool in earlier studies. Two distinct regimes are identified for different Ohnesorge numbers. Finally, we performed a scaling analysis, and rationalized the independence of the coefficient of restitution in the inertial regime and its decrease with increasing Ohnesorge number for large size ratios. Scaling for different energy forms is proposed and scaling for the restitution coefficient is determined using energy balance.
Most of the previous studies on droplet impact on a droplet resting on a solid surface, including the present study, are limited to the similar liquid. The impact of droplets of dissimilar liquids on solid surfaces can be performed in the future. The numerical model developed in this work has the potential to analyse the rich subject of capillary rebounds in nature and industries.
Funding
This work was supported by grants from the Research Grants Council of the Hong Kong Special Administrative Region, China (PolyU P0039589 and PolyU P0046985 for ECS project funded in 2023/24 Exercise).
Declaration of interests
The authors report no conflict of interest.
Appendix
To ascertain whether droplet interfaces, obtained using the 2-VOF approach, overlap, we observed the minimum gap (i.e. 
$h_c$) (representing the air layer) between a falling droplet and a sessile droplet. In figure 13 the evolution of the normalized minimum gap (normalized with 
$R_f$) between the droplets is shown over dimensionless time. As the falling droplet approaches the sessile droplet, the minimum gap sharply decreases. After a collision, the gap remains nearly constant and it increases rapidly upon bouncing. Notably, the normalized minimum gap between the droplets during bouncing is approximately of the order of 
$\Delta x/2$ (where 
$\Delta x$ varies from 0.01 to 0.04 that is the minimum normalized cell size) during collision, consistent with findings from previous studies (He et al. Reference He, Xia and Zhang2019; Ramírez-Soto et al. Reference Ramírez-Soto, Sanjay, Lohse, Pham and Vollmer2020). It can be seen that the minimum gap in the simulations scales with 
$\Delta x/2$ even across different mesh resolutions. Since the gap remains consistently positive, it indicates that the interfaces of the droplets do not overlap. While decreasing the cell size (
$\Delta x$) in simulations can reduce the air gap between droplets, the gap size cannot be reduced below a value of the order of 
$\Delta x/2$, as demonstrated in previous studies (He et al. Reference He, Xia and Zhang2019; Ramírez-Soto et al. Reference Ramírez-Soto, Sanjay, Lohse, Pham and Vollmer2020). Therefore, the minimum gap value in the simulation can be quite larger than its value in the experiment. As observed in the validation section (§ 3.1), the simulated bouncing behaviour closely resembles the experimentally observed behaviour. Therefore, it is argued that the difference in the minimum gap value does not significantly impact the simulation results. To determine the normalized minimum gap between droplets, the interface locations of both the falling and sessile droplets are identified where the volume fraction value is 0.5. Subsequently, the distance of each point on the sessile droplet's interface is measured from every point on the interface of the falling droplet. From these distances, the minimum value is extracted, representing the minimum gap between the droplets.
Figure 13. The variation of the normalized minimum gap (
$\tilde {h}_c$) between the falling droplet and the sessile droplet with dimensionless time at different cell sizes (
$We = 4.0$, 
$\varDelta = 10$).
Open access
