Analysis of high-speed drop impact onto deep liquid pool

Abstract

The present work is devoted to the analysis of drop impact on a deep liquid pool, focusing on the high-energy splashing regimes caused by large raindrops at high velocities. Such cases are characterized by short time scales and complex mechanisms, thus they have received very little attention until now. The BASILISK open-source solver is used to perform three-dimensional direct numerical simulations. The capabilities of octree adaptive mesh refinement techniques enable capturing of the small-scale features of the flow, while the volume of fluid approach combined with a balanced-force surface-tension calculation is applied to advect the volume fraction of the liquids and reconstruct the interfaces. The numerical results compare well with experimental visualizations: both the evolution of crown and cavity, the emanation of ligaments, the formation of bubble canopy and the growth of a downward-moving spiral jet that pierces through the cavity bottom, are correctly reproduced. Reliable quantitative agreements are also obtained regarding the time evolution of rim positions, cavity dimensions and droplet distributions through an observation window. Furthermore, simulation gives access to various aspects of the internal flows, which allows us to better explain the observed physical phenomena. Details of the early-time dynamics of bubble ring entrapment and splashing performance, the formation/collapse of bubble canopy and the spreading of drop liquid are discussed. The statistics of droplet size show the bimodal distribution in time, corroborating distinct primary mechanisms of droplet production at different stages.

1. Introduction

The impact of raindrops on a deep liquid pool has been extensively studied since the initial works of Worthington (1883, 1908). Various behaviours were obtained, depending primarily on miscible or immiscible characteristics of the two liquids (Lhuissier et al. 2013; Castillo-Orozco et al. 2016), but also on their difference of density (Thomson & Newall 1886; Manzello & Yang 2002; Lherm et al. 2022; Villermaux 2022), the depth of the receiving volume (Macklin & Hobbs 1969; Wang & Chen 2000; Fedorchenko & Wang 2004) and the angle of the impact (Zhbankova & Kolpakov 1990; Okawa, Shiraishi & Mori 2008; Gielen et al. 2017; Liu 2018), to mention only a few parameters.

In the specific case of identical liquids of density $\rho _l$, with a $90^\circ$ impact of the drop falling in the air on a deep volume of target liquid, very different phenomena are still observed, when the liquid viscosity $\mu _l$ and surface tension $\sigma$ with air, the diameter $d$ and speed $U_0$ of the impact drop are varied. Schotland (1960) has first reported the primary effect of Weber number $We=\rho _l {U_0}^{2}d/\sigma$ on the physics of impact, i.e. the ratio of the drop kinetic energy to the energy required to deform the target liquid surface.

A few years later, Engel (1966, 1967) provided a detailed description of various behaviours obtained from the impact of waterdrops on water pools. Based on high-speed pictures and using white particles in the target liquid and red ink in the impact drop, the creation of a large cavity in the flat free surface and the rise of the target liquid in a cylindrical shape were illustrated. A bubble-thin cylindrical sheet of liquid is erected at the upper edge of this crown, which eventually necks in and closes in a bubble dome in the most energetic cases. After that, the cavity shallows and a jet forms at the cavity floor, which flows through the centre of the cavity in case it is open or merges with a downward jet coming from the top of the bubble dome in case the crown is closed. The author suggested that the cavity may vibrate at its maximal expansion if all the drop kinetic energy is not yet transformed into potential energy, and she also assumed that: (i) the maximum possible bubble height is equal to the cavity diameter; (ii) the pressure evolution below the cavity floor should explain the formation of the upward jet; (iii) the liquid of the initial drop is carried by this jet and will eventually form a secondary drop at the centre bottom of the cavity. In the following decades, this analysis has been progressively completed by additional experiments with increasing visualization capabilities (van de Sande, Smith & van Oord 1974; Rodriguez & Mesler 1985; Pumphrey, Crum & Bjo/Rno/ 1989; Rein 1993, 1996; Leng 2001; Manzello & Yang 2002; Yarin 2006; Fedorchenko & Wang 2004; Berberović et al. 2009; Bisighini et al. 2010; Murphy et al. 2015), leading to the characteristic laws of cavity growth and an improved understanding of the phenomena mentioned hereabove.

In the most recent works, attention has been especially focused on the very early stage of impact, where intriguing mechanisms of jet formation and bubble entrapment are observed. Immediately after impact, the contact line between the drop and the receiving liquid, also called the ‘neck region’, moves radially at high speed. For impact with sufficient initial energy, a thin ejecta sheet, composed primarily of the liquid coming from the target pool, is shot out horizontally from the base of the contact surface (Weiss & Yarin 1999; Thoroddsen 2002; Josserand & Zaleski 2003). This ejecta stretches radially and may break up into ligaments and droplets, resulting in the splash of fine sprays (Deegan, Brunet & Eggers 2007; Thoroddsen et al. 2011). For slightly higher Reynolds numbers $Re=\rho _l U_0d/\mu _l$, the presence of distinct two jets, namely the initial ejecta and the late-emerging lamella, were experimentally identified by Zhang et al. (2012) using fast X-ray phase contrast imaging, and the clear distinction between generations of droplets was consequently observed by these authors in a reduced pressure environment. Such interactions of jets were also analysed by Agbaglah et al. (2015) who used a combination of X-ray imaging and axisymmetric simulations. At even higher $Re$, more complicated scenarios of irregular splash were suggested in figure 2(c) of Thoroddsen (2002).

One important output is that this multiplicity of splash may be strongly related to the instabilities observed in the neck region. Indeed, the base of ejecta may become unstable under certain impact conditions, shedding one sign (Agbaglah et al. 2015) or alternate signs (Castrejón-Pita, Castrejón-Pita & Hutchings 2012; Thoraval et al. 2012) of vortex in the liquid below the neck region. A strong interplay between this vortex-shedding course and the inception/breakup of the early liquid sheet was suggested by several recent studies (Thoraval et al. 2013; Agbaglah et al. 2015; Li et al. 2018), but the mechanisms are still not elucidated. Bottom visualizations of drop impact on thin liquid films (Thoraval et al. 2013) have revealed non-axisymmetric behaviours in the process of vortex shedding and ejecta formation, showing that these phenomena break axisymmetry at fine scales and may be linked to three-dimensional mechanisms. Local streamwise vortices, generated at the sharp corners of ejecta at the very early stage of impact, were mentioned by these authors to possibly explain these non-axisymmetric effects. Using ultra-fast imaging, an early azimuthal instability was also detected by Li et al. (2018), which consists of small azimuthal waves that grow at the edge of the outer contact line before the inception of vortex shedding. This clearly poses a new challenge to numerical studies, which mostly assume axisymmetry (Thoraval et al. 2012; Agbaglah et al. 2015). The authors pointed out the need for three-dimensional simulations with sufficient resolutions at the finest scale to explore the vortex structure in the neck region and the interactions between the two corners of the base of the ejecta sheet, which may be a primary mechanism of this instability.

As for the later stages, various physical processes can be obtained, depending on the specific range of impact parameters. An exhaustive review of previous works devoted to the impact of drops on deep pools of the same low-viscosity liquid has been conducted by Murphy et al. (2015), as shown in figure 1. Five different regimes were identified by the authors, depending on the values of the Froude $Fr={U_0}^2/gd$ and Weber numbers, where g is the gravitational acceleration. The two directions of variations of impact speed $U_0$ and drop diameter $d$ are plotted in figure 1, as well as the specific case of raindrops of different sizes falling at terminal speed (solid black line labelled ‘Raindrop TS’). As can be seen in figure 1, if the size and speed of the impact drop increase (moving typically from bottom left to top right on the chart), the following successive regimes can be obtained: (i) $\times$ C&VR, coalescence of the slow drop with the liquid volume, generating a vortex ring that moves downwards as the drop sinks; (ii) $\bigcirc$ S&TJ and $\Delta$ RE, development of a cavity together with a surface wave, forming an upward thin jet from the cavity floor with occasional bubble entrapment; (iii) $\square$ C&TJ, growth of a crown rim that ejects secondary droplets and a deeper crater that rebounds a thicker jet; (iv) $\lozenge$ BC, a large cavity is obtained for the highest speeds and the largest diameters, around which extends vertically a very thin liquid crown that continuously ejects liquid ligaments and secondary droplets. The elevation and radius of the crown keep increasing and its upper part eventually necks in, enclosing a large volume of air.

Figure 1. Various splashing behaviours induced by drop impact on a miscible liquid pool reported in the literature, reproduced from Murphy et al. (2015) with the authorization of the authors: $\times$ C&VR, coalescence and vortex ring; $\Delta$ RE, regular bubble entrainment; $\bigcirc$ S&TJ, swell and thin jet; $\square$ C&TJ, crown and thick jet; $\lozenge$ BC, bubble canopy. Solid black line (Raindrop TS), 0.4–5.8 mm raindrops falling at terminal speed (Gunn & Kinzer 1949); solid grey line (Drop TS ($C_d=1$)), constant $Fr$ for drops falling at terminal speed with an assumed drag coefficient $C_d=1$; dashed line, onset of the bubble canopy regime at $We=2000$. The $d-U_0$ axes indicate the directions of increasing drop diameter and drop speed respectively. Dotted lines, constant drop diameter of $d=0.5$ and 10 mm; dash-dot lines, constant drop speeds of $U_0=0.5$ and $10\ {\rm m}\ {\rm s}^{-1}$.

This last type of behaviour ($\lozenge$ BC in figure 1) is the most energetic regime and only received very limited attention during the past twenty years (Pan et al. 2008; Bisighini et al. 2010; Murphy et al. 2015; Sochan et al. 2018). Such high-energy regimes are usually characterized by shorter time scales and increased complexity of the phenomena involved in the splashing, which makes both experimental and numerical approaches more challenging, and may explain why it has received less attention before. However, this is the case likely to yield the most abundant physical processes and produce the greatest number of aerosol droplets, which is essential for many fundamental techniques (Liang & Mudawar 2016; Breitenbach, Roisman & Tropea 2018; Benther et al. 2021; Lohse 2022) and environmental issues (Murphy et al. 2015; Castillo-Orozco et al. 2016; Bourouiba 2021; Dasouqi, Yeom & Murphy 2021). Indeed, it was estimated that such cases could produce at least 2000 very fine droplets from its initial intrusion (Murphy et al. 2015) and more than 900 droplets larger than $100\ \mathrm {\mu }{\rm m}$ from the later crown expansion (Blanchard & Woodcock 1957), but the details of the production of these tiny droplets and their statistics remain to be found, which motivated the use of numerical simulations in this paper.

Multiphase flow calculations of drop impact usually consist in laminar simulations, as most previous studies were focused on impact conditions at $Re<10\,000$, and all motions resulting from the impact are characterized by very short characteristic times, so there is no transition to turbulence. The challenge of such computations is thus mostly related to the prediction of multiphase structures. Early numerical works (Harlow & Shannon 1967; Ferreira, Larock & Singer 1985) based on marker and cell and SOLution Algorithm-Volume Of Fluid (SOLA-VOF) techniques of interface tracking could provide some first predictions of the main features of drop impact but could not resolve the small-scale mechanisms. After that, Oguz & Prosperetti (1989, 1990) included the surface-tension effects and Morton, Rudman & Jong-Leng (2000) solved the full Navier–Stokes equations, opening the way to modern simulations that are based either on interface tracking (level-set, front tracking) or interface reconstruction using a transport equation for the volume fraction of one component (VOF method).

Recent simulations of liquid drops impinging onto liquid surfaces are mostly based on the latter category of models. Substantial axisymmetric simulations of the problem (Morton et al. 2000; Josserand & Zaleski 2003; Berberović et al. 2009; Thoraval et al. 2012; Ervik et al. 2014; Agbaglah et al. 2015; Ray, Biswas & Sharma 2015; Deka et al. 2017; Fudge, Cimpeanu & Castrejón-Pita 2021; Alventosa, Cimpeanu & Harris 2023; Fudge et al. 2023) can be found in the literature, focusing on those phenomena that are assumed axisymmetric. For three-dimensional (3-D) calculations, whose primary objective is to capture the non-axisymmetric mechanisms involved in the splashing process, efforts were focused on dynamic refinement techniques (Popinet 2003; Nikolopoulos, Theodorakakos & Bergeles 2007), in order to resolve multiple interfaces resulting from the splash, while ensuring a reasonable grid size. As a workaround, some simulations that only calculated half or a quarter of the 3-D impact have been performed to glimpse into certain 3-D phenomena (Rieber & Frohn 1999; Brambilla & Guardone 2015; Cheng & Lou 2015; Guo & Lian 2017). However, full 3-D direct numerical simulations of drop impact remains very challenging due to the wide range of scales involved as well as the substantial computational resources required (Wu et al. 2021). Only a limited number of full 3-D simulations (Cheng & Lou 2015; Chen et al. 2020; Constante-Amores et al. 2023) have been conducted under the minimum resolution of $96\sim 200$ cells per drop diameter, focusing mainly on the ‘crown splash’ regime at low-energy impact conditions ($Re<2000$ and $We<400$), so it may not directly apply to the highly energetic impact investigated in the present study. For impact at $Re>6000$, complex jet and splashing dynamics has been confirmed by high-speed camera (Thoroddsen 2002; Thoroddsen et al. 2011; Castrejón-Pita et al. 2012; Thoraval et al. 2013; Li et al. 2018), and it has been shown that such intricate phenomena can only be captured under resolutions of at least 1000 cells per drop diameter in previous axisymmetric simulations (Thoraval et al. 2012, 2013), which further complicates the numerical investigations of high-speed drop impact in full 3-D configurations.

In this study, we perform high-resolution direct numerical simulations of a large drop impacting a deep pool of the same liquid at high velocity, which is representative of raindrops at the ocean surface. This configuration has been experimentally studied by Murphy et al. (2015) in detail, primarily as a reference case for a study focused on the influence of oil slicks and oil dispersants on the impact. In their experiments, high-speed videos enabled characterizing the time evolution of the external shape of the cavity and the upper rim of the crown, while microscopic holography provided some statistics of droplets in a specific observation window, which can be used for the validation of our present numerical strategies (§ 3). In the current numerical investigation, further attention is given to the early-time splashing dynamics near the neck region, the internal mechanisms of crater evolution and the distribution of secondary droplets over a larger domain, to explore additional flow details of the high-speed splashing process in various directions.

The paper is organized as follows: the numerical methods and problem statement are described hereafter in § 2, the detailed validation of numerical strategies is conducted in § 3 by comparing the numerical results with valuable experimental data provided by Murphy et al. (2015), the analysis of several mechanisms involved in the splash is presented in § 4 and the statistics of airborne droplets are finally analysed in § 5.

2. Numerical methods and flow configurations

2.1. Main features of the solver

The drop impact of the gas–liquid system is considered as incompressible flow, and it solves a system of the mass balance equation, the momentum balance equation and the advection of one-fluid formulation, called hereafter the colour function

(2.1)$$\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{U}=0, \end{gather}$$

(2.2)$$\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{\mathsf{D}})+\rho\boldsymbol{a}+\sigma\kappa\delta_s\boldsymbol{n} , \end{gather}$$

(2.3)$$\begin{gather}\frac{\partial C}{\partial t}+\boldsymbol{\nabla}\boldsymbol{\cdot}(C\boldsymbol{U})=C\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{U}, \end{gather}$$

where $\boldsymbol {U}$ is the velocity vector, $\rho$ is the density, $P$ is the pressure, $\mu$ is the viscosity, $\boldsymbol{\mathsf{D}}=[\boldsymbol {\nabla }\boldsymbol {U}+(\boldsymbol {\nabla }\boldsymbol {U})^{\rm T}]/2$ is the deformation tensor and $\boldsymbol {a}$ is the body force along the impact direction. The last term in (2.2) is the surface-tension force, with $\kappa$ the curvature of the interface, $\sigma$ the surface-tension coefficient, which is taken as constant in the present study, and the $\boldsymbol {n}$ the unit vector normal to the interface. This term is zero everywhere but at the gas/liquid interface, as controlled by the Dirac function $\delta _s$. The colour function $C$ is transported by (2.3). For incompressible flow, the right-hand term in (2.3) is zero due to (2.1).

The BASILISK framework is a flow solver developed by Popinet (2015) and available for use and development under a free software GPL license (Popinet & collaborators 2013–2023). It solves the time-dependent compressible/incompressible variable-density Euler, Stokes or Navier–Stokes equations with second-order space and time accuracy. The momentum-conserving VOF approach implemented by Fuster & Popinet (2018) is employed here to simulate the problem of gas–liquid two-phase flow. The colour function (2.3) is solved based on a volume fraction advection scheme proposed by Weymouth & Yue (2010), which exhibits complete mass conservation and makes it ideal for highly energetic free-surface flows. The interface then can be represented using the piecewise linear interface capturing VOF method (Rudman 1998). For (2.2), the Crank–Nicholson discretization of the viscous terms is second-order accurate in time and unconditionally stable, while the convective terms are computed using the Bell–Colella–Glaz second-order unsplit upwind scheme (Bell, Colella & Glaz 1989), which is stable for Courant–Friedrichs–Lewy (CFL) numbers smaller than one. The calculation of surface tension is one of the most challenging steps of the process, since no continuous definition of the interface is available in the VOF approach. Here, the balanced-force surface-tension calculation is used (Francois et al. 2006), which is based on the continuum-surface-force approach originally proposed by Brackbill, Kothe & Zemach (1992). In addition, a second-order accurate calculation of the curvature is performed, using the height-function technique developed by Popinet (2009). More detailed descriptions of the numerical schemes can be found in Fuster & Popinet (2018), Popinet (2018), Pairetti et al. (2020) and Zhang, Popinet & Ling (2020).

Cubic finite volumes organized hierarchically in octree are used for space discretization. The octree structure has been developed initially for image processing (Samet 1990) and later applied to computational fluid dynamics (CFD) (Khokhlov 1998) and multiphase flows (Popinet 2003). The basic organization is the following, all details can be found in Popinet (2003, 2015): when a cell is refined, it is divided into 8 cubic cells, whose edges are half the ones of the parent cell. The base of the tree is the root cell and the last cells with no child are the leaf cells. The cell level is the number of times it is refined, compared with the root cell (level 0). To avoid too much complexity in the gradient and flux calculations, the levels of direct and diagonally neighbouring cells are constrained and cannot differ by more than one and all cells directly neighbouring a mixed cell must be at the same level.

BASILISK employs the adaptive mesh refinement (AMR) (van Hooft et al. 2018) to adaptively refine/coarsen the grids based on the wavelet-estimated numerical error of the local dynamics, which makes it especially appropriate for the present application where multiple interfaces and numerous droplets and bubbles are expected. The resolution will be adapted every time step according to the estimated discretization error of the spatially discretized fields (volume fraction, velocity, curvature, etc.). The mesh will be refined as long as the wavelet-estimated error exceeds the given threshold, which eventually leads to a multi-level spatial resolution from the minimum refinement level $L_{min}$ to the maximum refinement level $L_{max}$ over the entire domain.

The BASILISK code (or its ancestor GERRIS) has already been applied to the study of drop impact onto liquid surfaces (Thoraval et al. 2012; Agbaglah et al. 2015; Fudge et al. 2021; Wu et al. 2021; Sanjay et al. 2023; Fudge et al. 2023), and the great superiority of parallel capacity and computational efficiency of the solver were discussed in Wu et al. (2021). Furthermore, the capabilities of the BASILISK solver have been extensively validated on various problems of multi-phase flows in Popinet & collaborators (2013–2023).

2.2. Initial flow configuration

The configuration of drop impact investigated in the present study mimics the one (control case) studied by Murphy et al. (2015) using high-speed video. In the experiments, the drop falls in a $15.2\times 15.2\ {\rm cm}^{2}$ tank filled with seawater to a depth of 8 cm. The measured horizontal and vertical diameters of the drop just before impact are $d_h=4.3$ mm and $d_v=3.8$ mm, respectively, which results in an effective drop diameter $d=(d_v d_h^2)^{1/3}=4.1$ mm. The ratio between the width of the tank $L$ and the drop diameter $d$ is therefore defined as $L=38d$. The speed before impact is measured as $U_0=7.2\ {\rm m}\ {\rm s}^{-1}$, which is approximately 81 % of the drop terminal speed. The properties of seawater are: density $\rho _l=1018.3\ {\rm kg}\ {\rm m}^{-3}$, dynamic viscosity $\mu _l=0.001\ {\rm N}\ {\rm s}\ {\rm m}^{-2}$ and surface tension $\sigma =0.073\ {\rm N}\ {\rm m}^{-1}$.

In our present numerical simulation, the computational domain is reduced down to a cube with a side length $L=16d$, which represents approximately 1/14 of the water tank volume in the experiments, as shown in figure 2(a). It has been found to be the best compromise to avoid any effect of the boundary conditions on the splash, while decreasing as much as possible the dimensions of the domain. The free surface is located at the mid-distance between the bottom and the top, thus the depth of the pool is $H=8d$, to give enough space to the aerosolized droplets without any interaction with the boundaries. The free outflow boundary condition is imposed on the top of the domain, while the default slip boundary condition (symmetry) is applied for the four sidewalls and the bottom. A zoomed-in view of the initial flow set-ups in the vicinity of the drop is depicted in figure 2(b). The shape and impact velocity of the downward-moving drop are initialized to be the same as those in the corresponding experimental set-ups. The initial gap between drop and pool is $\delta =0.1d$, allowing the observation of air sheet entrainment near the contact line. The liquid in the drop and the tank are the same (seawater), which is assumed to have almost no effect on the splashing behaviours, as the properties of the freshwater drop and the target seawater in the experiments are very close. The density and viscosity ratios between seawater and air are $\rho _l/\rho _g=783$ and $\mu _l/\mu _g=56$, which leads to a system of drop-pool impact with dimensionless numbers $Re=30\,060$, $We=2964$ and $Fr=1290$.

Figure 2. Initial numerical configurations of 3-D simulation. (a) Overall view of the computational domain and the initial mesh structure at a plane across the centre of the impact drop ($z=0$). (b) Closeup view of initial flows around the impact drop. (c) Mesh refinement strategy at the initial stage ($S1$). A higher maximum level of refinement at $L_{max}=14$ is imposed near the neck region (green area) to capture the early-time splashing behaviours, while $L_{max}=13$ is employed for the rest of the domain.

As illustrated in figure 2(a), the initial mesh configuration around the drop is generated based on the AMR algorithm using the estimated discretization errors of the volume fraction ($\,f_{Err}=1e-6$) and the velocity ($u_{Err}=1e-4$) fields, which promises a rounded geometric description of the initial drop and lowers the RAM requirement of initialization. The mesh is coarsened gradually down to the given minimum level of refinement ($L_{min}$) away from the drop interface. The region around the free surface of the pool is initially refined at $L_{max}=11$ to avoid any divergence issue. Once the simulation starts, the mesh will be redistributed adaptively based on AMR, using the volume fraction field with tolerance $f_{Err}=1e-4$ and the velocity field with tolerance $u_{Err}=1e-1$ as adaption criteria. Additionally, we remove those droplets that approach the boundaries of the computational domain, since these tiny droplets have few effects on the evolution of the main impact but are expensive to track. In real experiments, they are considered to evaporate at one point and this is not our focus in this paper. Using the initial contact centre as the reference point, any droplet whose centroid lies outside of the region of a semi-sphere with diameter $D_{rm}=15d$ will be removed from the simulation. The effect of gravity is taken into consideration in this study.

2.3. Minimum spatial resolution

Direct numerical simulation consists in resolving all scales of the flow, from the smallest relevant ones (here, the smallest water droplets ejected in the air or air bubbles entrained in the water) up to the large scales of the problem (here, the volume of water tank that receives the initial impact drop). Ideally, it means that the grid resolution of the air/water domain should be fine enough to capture all air/water interfaces created at all steps of the splashing process.

We have carried out extensive tests to explore the effect of minimum spatial resolution on the splashing behaviours of drop impact (see Appendix A), and it has been found that the early dynamics of splash and air entrapment in the neck region affects significantly the subsequent phenomena at high $Re$ and $We$ impact conditions. Our preliminary tests suggest that a grid resolution of at least 1024 cells per drop diameter is essential to capture the converged primary features of ‘prompt splash’ associated with early-time droplets/bubbles productions (see Appendix A.1 for details), but remains insufficient to fully resolve the early-time small-scale features in the neck region based on the energy aspect of view (see Appendix A.2 for details), highlighting the increased challenge of 3-D numerical study for such an intricate physical process. This will need to apply at least a maximum refinement level at $L_{max}=14$, corresponding to an equivalent uniform grid of more than 4.3 trillion $({(2^{14})}^3)$ cells. Note that the more the maximum refinement level is increased, the more the time steps are reduced to comply with the CFL condition, which means that the calculation CPU time is not proportional to the number of cells. Although the total number of cells decreases significantly by employing the AMR algorithm, it is still far too expensive to perform a long-term simulation in full 3-D configurations at this level.

Therefore, based on the primary physical characteristics of the impact at different times, we have divided the simulation into three consecutive stages, namely the early-time splashing stage at $t\leq 0.2$ ms ($S1$), the crown formation stage at $0.2< t<4$ ms ($S2$) and the bubble canopy stage at $t\geq 4$ ms ($S3$). At stage $S1$, a very thin liquid sheet emerges on the neck region and interacts strongly with drop and pool, rupturing and producing numerous very fine droplets and bubbles, thus the primary objective here is to capture the flow dynamics near the contact region in finer scales. Figure 2(c) shows the mesh refinement strategy employed at the initial stage. A higher maximum level $L_{max}=14$ ($\varDelta =3.9\ \mathrm {\mu }{\rm m}$) is enforced in the vicinity of the impact centre to solve the neck dynamics at smaller scales (green area), while $L_{max}=13$ ($\varDelta =7.8\ \mathrm {\mu }{\rm m}$) is employed for the rest of the domain to capture the general dynamics. At stage $S2$, a coherent liquid sheet is developed above the extra refinement layer and secondary droplets are generated incessantly from the top of the crown. Thus, we remove the extra refinement layer here and use $L_{max}=13$ for the entire computational domain to capture the dynamics of crown fragmentation and droplet shedding. At stage $S3$, the droplets pinched from the top of the crown are much fewer and generally larger, compared with the ones from previous stages. Therefore, we restart the simulation with $L_{max}=12$ ($\varDelta =15.6\ \mathrm {\mu }{\rm m}$) over the whole computational domain, which allows capturing of the main physical dynamics of crown and cavity until the end of the simulation ($t=48$ ms). This mesh refinement strategy ensures that the simulation is doable in a full 3-D configuration and can be accomplished in a ‘reasonable’ time, considering the available computational resources for the moment.

All the numerical results presented in the main text of this paper were performed on 1024 cores for 33.5 days, which consumes more than $8.21 \times {10}^5$ CPU-hours in total, using the computational resources on Advanced Research Computing at Virginia Tech.

3. Comparisons with the experiments

In this section, we validate the numerical scheme applied in the present study by conducting qualitative and quantitative comparisons between simulated results and experimental data provided by Murphy et al. (2015).

3.1. Morphology

Figure 3(a) shows the global evolution of air–water interfaces generated by high-speed drop–pool impact during 48 ms after contact. The side views of numerical results (bottom) are compared with experimental high-speed images (top) in time. The drop is initialized above the surface of the pool and the time point of contact is defined as $t=0$ ms. Once the simulation starts, the drop will fall downwards to hit the receiving pool driven by the combination of initial impact speed and gravitational force. Directly after impact ($t=1$ ms), a cylindrical wave around a flat-bottomed disk-like cavity is produced due to the initial violent penetration, and substantial liquid ligaments emanate almost horizontally from the thickened rim of the crown, spraying a large number of micro-droplets into the air. In the following few milliseconds, the drop keeps expanding and eventually spreads into a thin layer of liquid that is distributed along the interior surface of the cavity, stretching the cavity into a typical hemispherical shape that has been widely reported in the literature (Engel 1966; Berberović et al. 2009; Bisighini et al. 2010; Ray et al. 2015), as seen at $t=3$ ms. By $t=7$ ms, the upper rim of the crown has started to proceed inwards and the orientation of ligaments has transitioned from almost horizontal to vertical, which eventually leads to the upcoming closure event on the upper part. At the time when the crown necks in ($t\approx 12$ ms), flows from the sidewalls of the cylindrical wave meet along the impact axis and a large volume of air is encapsulated, generating a central liquid jet that moves spirally from the merging point ($t=18$ ms). The downward-moving jet keeps growing and eventually pierces the cavity floor ($t=37$ ms), disturbing the retraction of the cavity bottom. Finally, the rebound of the cavity bottom produces a broad upward jet that merges with the previous central column of fluid, leaving several air bubbles inside the liquid tank, evidenced at $t=48$ ms.

Figure 3. Qualitative comparisons between numerical and experimental results. (a) Overall evolution of air–water interfaces during 48 ms after impact. From left to right, the experiment shows $-1$, 1, 3, 7, 12, 18, 41 and 52 ms after impact, and the simulation shows $-0.07$, 1, 3, 7, 12, 18, 37 and 48 ms after impact. The red stars indicate the tracked positions of the upper rim of the crown. The scale bar is 10 mm long. (b) Closeup view of early-time splashing behaviours during $450\ \mathrm {\mu }{\rm s}$ after impact. From left to right, the experiment shows 49, 148, 246, 345 and $443\ \mathrm {\mu }{\rm s}$ after first contact, and the simulation shows 50, 150, 250, 350 and $450\ \mathrm {\mu }{\rm s}$ after first contact. The scale bar is 1 mm long. The qualitative comparisons show that the simulation successfully reproduced all the distinctive features observed in the experiments. See also supplementary movie available at https://doi.org/10.1017/jfm.2023.701.

A close-up view of the early-time splashing behaviours near the neck region during $450\ \mathrm {\mu }{\rm s}$ after contact is provided and compared with high-speed experimental holograms in figure 3(b). As shown in the first frame ($t\approx 50\ \mathrm {\mu }{\rm s}$), a cloud of very fine droplets is scattered underneath the drop immediately after contact and no coherent ejecta sheet is observed, which looks very similar to the irregular splash of microdroplets for the case at $Re=29\,000$ and $We=1800$ in Thoroddsen (2002) (their figure 2c), also known as ‘prompt splash’ (Deegan et al. 2007). After the initial irregular sprays, more target liquid is gradually pushed out by the impinging drop, forming a thin-walled crown that continuously sheds secondary droplets from its rim. The morphological behaviours of the early-time splash captured by numerical simulation are quite consistent with experimental observations.

These qualitative comparisons show that our simulation reproduced all the distinctive features observed in the experiments. The correct prediction of the early-time splashing behaviours, the transition of the ligaments’ orientation and the exact time when the upper rim of the crown necks in and the central spiral jet pierces the cavity bottom are especially convincing. A general good agreement to the experiments is obtained.

3.2. Kinematics of crown and cavity

The kinematic behaviours of the crown and cavity are manually tracked from experiments and simulations, and compared quantitatively in figure 4. Using the impact centre as a reference point, the time evolution of the upper rim of the crown (marked by red stars in figure 3), its radial distance (figure 4a) and height (figure 4b), are measured during 24 ms after impact. The height is defined as the distance between the initial quiescent free surface of the pool and the position where ligaments are formed. The trajectory of the upper rim of the crown is then plotted in figure 4(c). For the first few milliseconds, the crown expands rapidly along the horizontal radial direction and reaches its maximum radius very soon ($t\approx 3$ ms). After the maximum horizontal position, the upper rim starts to proceed inwards under the effect of surface tension. Meanwhile, the leading edge of the crown rises continuously in the vertical direction and approaches its maximum height right before it necks in ($t\approx 12$ ms). After the closure of the upper part, this point vibrates slightly near the impact axis along with the shrinking of the entrapped large air bubble.

Figure 4. Analysis of quantitative data measured with respect to the initial impact centre. (a) Evolution of the crown radius. (b) Evolution of the crown height. (c) Trajectory of the upper rim of the crown. (d) Evolution of the cavity radius. (e) Evolution of the cavity depth. (f) Evolution of the cavity volume. The black dashed line in (e) shows the theoretical prediction of penetration depth using the proposed model in Bisighini et al. (2010). The error bars indicate the standard deviation in experimental data.

The evolution of the submerging cavity has been thoroughly discussed in experiments, simulations and theories, and substantial attention has been particularly given to the estimation of the cavity geometric dimensions in previous studies (Engel 1967; van de Sande et al. 1974; Prosperetti & Oguz 1993; Leng 2001; Berberović et al. 2009; Bisighini et al. 2010; Jain et al. 2019). Figure 4(d) demonstrates the temporal variation of the horizontal radius at the intersection edge between the cavity and the initial free surface. The cavity grows rapidly on the horizontal plane at the beginning due to the initial impinging momentum and the outward-expanding speed decreases gradually. After the closure of the upper crown, the increase of the cavity radius slows down and is nearly linear. Starting from $t\approx 3$ ms, the horizontal radius at the bottom of the crown becomes wider than at the top, which may presumably provide an overall momentum towards the central axis, thus propelling more liquid from the receiving pool to the sidewalls of the crater. Figure 4(e) plots the time evolution of the cavity depth measured from the lowest point of the cavity bottom, showing good agreement between numerical, experimental and theoretical (Bisighini et al. 2010) results. The cavity keeps expanding in depth and reaches its maximum position 24 ms after impact, which takes almost double the time of the maximum crown position ($t\approx 12$ ms). Finally, the volume of the cavity is calculated as half of an ellipsoid as same as in Murphy et al. (2015), as shown in figure 4(f).

Figure 4 confirms this conclusion: both the trajectory of the upper rim of the crown and the geometric dimensions of the submerging cavity are found in very good agreement with experimental measurements. A very reliable quantitative agreement is obtained between simulation and experiment.

3.3. Droplets in observation window

For further validating the present numerical strategies, the statistics of secondary droplets in a small field of view are extracted from the simulation and compared with experimental measurements. The specific location of the observation window is shown in figure 5(a), where a $10\times 10\ {\rm {mm}}^{2}$ square field of view is placed 13 mm above the free surface of the pool, which is similar to the experimental set-up of the high-speed camera in Murphy et al. (2015). In their experiments, the holographic frames for droplet analysis were selected at the time when the first upward-rising droplet exits the top of the observation window ($t\approx 3\sim 4$ ms), and droplet statistics were ensemble averaged over all replicates (more than 25). For our numerical results here, the time-averaged droplet statistics are obtained using 100 time slices over the time window $t\in$ [3, 4] ms. Figure 5(c) shows the distribution of secondary droplets in the vertical direction within the observation window. Similar to the experimental measurements, most droplets are still concentrated at a lower position, which can be also clearly seen in the close-up view of the observation window at $t=3.5$ ms in figure 5(b). Figure 5(d) plots the time-averaged droplet size distribution. Compared with the bimodal distribution in experiments, an inapparent secondary peak can be still found from the numerical results under relatively wider bins (doubling the size of bins used in Murphy et al. 2015). Two primary plateaux centred around $50\ \mathrm {\mu }{\rm m}$ and $225\ \mathrm {\mu }{\rm m}$ are observed, which is in good agreement with experimental data. The successful reproduction of droplet statistics in a specific observation window is highly encouraging, as it gives confidence for further comprehensive and in-depth analysis over a broader spatial domain and distinct temporal stages in the following sections.

Figure 5. Comparisons of droplet statistics between numerical and experimental data captured in a specific field of view during $3\sim 4$ ms after impact. (a) Overall schematic view of the relative position of the observation window. (b) Closeup view of secondary droplets in the observation window at $t=3.5$ ms. (c) Vertical distribution of secondary droplets. (d) Size distribution of secondary droplets. The numerical droplet statistics presented in (c,d) are time-averaged data using 100 time slices over the time window $t\in [3, 4]$ ms. The experimental data are ensemble averaged using more than 25 replicates as originally presented in figures 17 and 18 in Murphy et al. (2015).

4. Overall dynamics of splashing

4.1. Early-time dynamics

The very early splashing dynamics that occurs shortly after contact (${<}100\ \mathrm {\mu }{\rm s}$) in the process of drop impact onto a liquid surface has been widely discussed during the past twenty years. Conventionally, for low impact velocities, an air disc is entrapped under the centre of the impact drop by lubrication pressure and it later breaks up into chains of micro-bubbles (Thoroddsen et al. 2012; Tran et al. 2013) or contracts into several bubbles along the central line (Thoroddsen, Etoh & Takehara 2003; Jian et al. 2020). As the contact line expands radially, sheet-like liquid ejecta is sent out nearly axisymmetrically from the outer edge of the neck and possibly emits rings of tiny droplets from its rim for sufficiently large Reynolds numbers (Weiss & Yarin 1999; Thoroddsen 2002; Josserand & Zaleski 2003; Howison et al. 2005; Deegan et al. 2007; Marcotte et al. 2019). However, with increased impact velocity ($Re>7000$), azimuthal undulations were found to grow at the base of the ejecta and the entrapment of bubble rings was observed near the neck region (Thoraval et al. 2013; Li et al. 2018), which thereby breaks the axisymmetry of the motions in fine scales. For even higher impact velocities, irregular splash and intricate vortex-shedding progress were experimentally observed (Thoroddsen 2002; Castrejón-Pita et al. 2012; Thoraval et al. 2012). Nevertheless, the study of such a complex flow structure remains very challenging for both experimental and numerical approaches as it mainly occurs in microscopic length within the time scale of several microseconds after contact.

Figure 6 shows the simulated early-time dynamics in the neck region between drop and pool from the bottom view. Here, it is focused on the stage before the outer edge of the neck $R_n$ overtakes the outline of the impact drop $R_d$, where $R_d=d_h/2$. The first frame and the last frame reach $25\,\%$ and $89\,\%$ of the drop size respectively. The air–water interfaces are coloured by the volume fraction of the passive tracer (see § 4.4) and the opacity of the interfaces is set to 0.5, which enables us to visualize the interfacial dynamics for both inside and outside regions of the neck.

Figure 6. Early-time neck dynamics near the contact region induced by high-speed drop impact onto deep liquid pool observed from the bottom view. From left to right and top to bottom, the first three frames are shown 4, 6 and $8 \mathrm {\mu }{\rm s}$ after first contact, where the ‘nearly axisymmetric’ bubble rings are entrapped from the neck of the connection. The black arrow points at the central air disc. The yellow arrow indicates the early-time entrapment of an air void. The red arrows show the formation of a new bubble ring. The last four images show a smaller magnification 10, 15, 32 and $50\ \mathrm {\mu }{\rm s}$ after impact. The outer edge is the downward-moving drop, the inner edge is the contact line of the neck and the central irregular disc is the entrapped air pocket. Azimuthal instabilities and liquid ejecta are developed along the neck. The green arrow indicates the entrapment of bubble arcs due to the oscillations of the base of the early fingers/ejecta. The orange arrow indicates the entrapment of bubble ring due to sheet impingement. The outer line of the neck has not reached the size of the impact drop here ($R_n< R_d$). The scale bar is $500\ \mathrm {\mu }{\rm m}$ long.

Upon contact, a central disc of air (black arrow at $t=4\ \mathrm {\mu }{\rm s}$ in figures 6 and 7a) is entrapped due to air pressure build-up between the pool and the drop, therefore cushioning the impact. For drop impact on the deep pool of the same liquid, Hendrix et al. (2016) showed that both the upper (drop bottom) and the lower (pool) interfaces react similarly to the local pressure increase in the air layer, and they deform nearly symmetrically. This leads to a doubled total volume of air pocket, in contrast to impact cases such as a rigid drop impact onto a pool or a liquid drop on a solid substrate, where only one deformable interface exists. The universal scaling law for the volume of the entrapped air pocket was given as $V_p/V_d\sim St^{-4/3}$, where $V_p$ is the air pocket volume, $V_d$ is the drop volume and $St=\rho _lR_dU_0/\mu _g$ is the Stokes number showing the competing effects of the viscous force in the air layer and the inertial force of the drop. The initial radius of the entrapped air pocket can be theoretically predicted using $R_p/R_d=3.8(4/St)^{1/3}$, as proposed by Hicks & Purvis (2010) and Hicks et al. (2012), where half of the impact speed $U_0/2$ was suggested later by Tran et al. (2013) and Hendrix et al. (2016) to evaluate the $St$ for liquid–liquid impact. For the present case, the estimated air pocket radius should be $\sim$0.17 mm, while the measured initial radius of the air pocket from numerical results is around $\sim$0.41 mm. One possible explanation for the slightly larger numerical air disc is the inadequate spatial resolution near the neck region, which results in the partially resolved early neck dynamics, as discussed in Appendix A.2. It is worth noting that the bottom curvature of the drop at the moment of impact may also alter the entrapment of air pocket as discussed in Scolan & Korobkin (2001) and Thoraval et al. (2013).

Figure 7. Flow field and vorticity structure in the vicinity of the neck region on the vertical slice at $z=0$ (see the dashed line at $t=10\ \mathrm {\mu }{\rm s}$ in figure 6). The red and blue colours represent counterclockwise and clockwise rotation respectively. (a) Entrapment of air disc and bubble rings $4\ \mathrm {\mu }{\rm s}$ after first contact. The black arrow points at the central air disc, and the yellow arrow shows the entrapment of axisymmetric air void in the neck, as indicated in the first frame of figure 6. The scale bar is $50\ \mathrm {\mu }{\rm m}$ long. (b) Formation of azimuthal fingers from the neck at $t=12\ \mathrm {\mu }{\rm s}$. Secondary droplets are emitted from its tips. Vortex shedding of the alternate signs from the base of the early fingers/ejecta generates a von Kármán-type structure along the drop/pool boundary, with occasional air bubbles/bubble arcs entrapment as indicated by the green arrow in figure 6. The scale bar is $100\ \mathrm {\mu }{\rm m}$ long. (c) Collision between the ejecta and the downward-moving drop leads to the entrapment of a large air ring at $t=73\ \mathrm {\mu }{\rm s}$. The orange arrow indicates the entrapment of bubble ring due to sheet impingement (see also the last frame of figure 6). The scale bar is $500\ \mathrm {\mu }{\rm m}$ long.

Subsequently, the outer line connecting the drop and pool expands rapidly in the radial direction, entrapping bubble rings along the drop/target boundary. In the second panel of figure 6, it can be seen that a new round of air cylinder (red arrow) is entrapped at the neck and is pinched off shortly into the bulk within the time scale of ${<}0.5\ \mathrm {\mu }{\rm s}$. By $t=8\ \mathrm {\mu }{\rm s}$, up to 10 bubble rings are present within the target volume and most of them are entrapped axisymmetrically as concentric circles. Note that the neck of connection between drop and pool remains rather smooth at this time ($t<10\ \mathrm {\mu }{\rm s}$) and no valid liquid ejecta is observed around it, implying that the fundamental mechanism of bubble ring entrapment here differs from the jet-induced air encapsulation predicted by Weiss & Yarin (1999). Figure 7 shows the air–water interface overlaid by the vorticity field near the neck region on the vertical intersection across the drop centre ($z=0$). As indicated in figures 6 and 7(a), an air void (yellow arrow) is rolled up along the entire circumference of the neck and later entrapped axisymmetrically inside the liquid phase. This is reminiscent of the entrapment of toroidal bubbles numerically captured by Oguz & Prosperetti (1989) for two drops collision, where the basic question ‘Does the contact line between two approaching surfaces move outwards fast enough to prevent further contact after the initial one?’ was discussed. At the moment of contact, a very small liquid bridge of radius $r_n$ with high curvature ‘meniscus’ connects two liquid masses and a thin outer air sheet retracts outwards rapidly driven by surface tension and local pressure gradient. A self-similar repeated reconnection of this air gap was predicted at the very early time of contact, thus enclosing a number of tiny toroidal bubbles in the liquid phase. Focusing on the viscous regime of $Re_l\ll 1$, $Re_l=\sigma R_n/(\rho _l\nu _l^2)$, the analytical analysis of Eggers, Lister & Stone (1999) drew the scaling laws for the radius of the entrapped toroidal bubble $r_b\propto {r_n}^{3/2}$ and the width of the thin air gap connecting the bubble $r_g\propto r_n^2$. Further investigation of Duchemin, Eggers & Josserand (2003) for inviscid drop coalescence explained that the occurrence of each pinch-off event ($i_{th}$) depends only on the local dynamics of air gap width $r_g^i$, where $r_g^i=(r_n^i)^2$. For each self-similar pinch-off succession, the distance between the tip of the meniscus and the reconnection point was given as $r_c^i=10(r_n^i)^2$ and the time interval can be estimated as $t_c^i=7.6(r_n^i)^3$. The author emphasized that the reconnection ceases when $r_n>0.05$, thus $r_c^i\leq 0.025$ and $t_c^i\leq 0.00095$, where the space and time coordinates are scaled by $R_d$ and $\sqrt {\rho _l{R_d}^3/\sigma }$. The evolution of the bubble rings was, however, not able to be predicted by the analytical model because of their highly non-circular shape and 3-D rotations and stretching. Based on this theoretical model, the distances and time intervals of the neighbouring bubble rings for the present case can be estimated ${\leq }53.75\ \mathrm {\mu }{\rm m}$ and ${\leq }11.18\ \mathrm {\mu }{\rm s}$. Measurements from our simulation show that the maximum distance and time interval between bubble rings are around ${\sim }31\ \mathrm {\mu }{\rm m}$ and ${\sim }0.5\ \mathrm {\mu }{\rm s}$, which are of the same order as analytical estimates. It should be noted that the existence of the central air disc as well as the variation of drop bottom curvature at the moment of impact may alter the local dynamics, but the phenomena should be qualitatively similar.

As the drop sinks, these air rings are stretched longitudinally while rotating and eventually break up into necklaces of micro-bubbles due to surface-tension Rayleigh instability as shown at $t=50\ \mathrm {\mu }{\rm s}$ in figure 6. According to Chandrasekhar (2013), for a stationary hollow gas cylinder of radius $r_b$ in the liquid, the most unstable wavelength is estimated as $\lambda _m=2{\rm \pi} r_b/0.484$, and the characteristic time of breakup is given as $t_\sigma =1.22\sqrt {\rho _l{r_b}^3/\sigma }$. In the second frame of figure 6, the two most recent bubble rings are measured to be entrapped with diameters ${\sim }16\ \mathrm {\mu }{\rm m}$, which would result in an estimated lifespan of $t_\sigma =3.3\ \mathrm {\mu }{\rm s}$ based on the seawater properties. However, our simulation shows that these bubble rings can maintain their tubular shapes for over ${\sim }29\ \mathrm {\mu }{\rm s}$ before eventually breaking, surpassing a duration of more than $8t_\sigma$. This discrepancy can be attributed to the stabilizing influence of rotation around the air cylinder, which effectively suppresses the growth rate of disturbances, thereby delaying the breakup of bubble rings (Rosenthal 1962; de Hoog & Lekkerkerker 2001; Ashmore & Stone 2004; Eggers & Villermaux 2008). Similar effects have also been reported in the experiments of Thoraval et al. (2013), who measured the breakup duration of bubble rings in the range $4\sim 12t_\sigma$ based on the properties of methanol.

During these early times of impact, the majority of the liquid remains unperturbed while the connection of two liquid masses propagates outwards instantaneously. Theoretically, the spreading law of $r_n=\sqrt {2\tau }$ can be derived from the truncated sphere approximation based solely on the geometric considerations (Rioboo, Marengo & Tropea 2002; Josserand & Zaleski 2003), but it violates the continuity equation. By incorporating Wagner's theory (Wagner 1932), the radial motion of the neck for drop impact on a solid surface has been deduced as $r_n=\sqrt {3\tau }$ (Riboux & Gordillo 2014; Philippi, Lagrée & Antkowiak 2016). However, the contact line dynamics during drop impact is more intricate in reality and involves potential interplays of various factors. Notably, both liquid viscosity (Thoroddsen 2002; Josserand & Zaleski 2003) and air entrapment (Mani, Mandre & Brenner 2010; Philippi et al. 2016) have been identified as important factors shaping the local dynamics under certain impact conditions, thus cushioning and delaying the contact. A novel dimensionless jet number, $J=St^2Re^3$ ($St=\mu _g/\rho _l\,{\rm d}U_0$), thus has been introduced by Josserand, Ray & Zaleski (2016), in order to characterize the interplay between the effects of liquid viscosity and gas cushioning on the impact dynamics. It has provided valuable insights into two asymptotic regimes of phenomena: (i) at small $J$ ($J\ll 1$), the air cushioning is found insignificant and the speed of the jet is mostly selected by the viscous balance of the liquid, resulting in a sequence of toroidal bubbles that are usually entrapped before the emergence of a liquid jet; (ii) at large $J$ ($J\gg 1$), the cushioning of air becomes more dominant, leading to a weaker jet velocity and a higher likelihood of simultaneous formation of the liquid jet right after air pocket entrapment. For the present studied case, the dimensionless numbers $Re=30\,060$ and $St=5.9\times 10^{-7}$ result in a value of $J=9.4$, which is close to one, indicating that liquid viscous and air cushioning could have comparable effects on the early neck dynamics. Therefore, here, we employ the form $r_n=R_n/R_d=C\sqrt {3(T-T_0)U_0/R_d}$ to describe the instantaneous motion of the neck (Jian et al. 2018; Li et al. 2018), where the origin of the time scale $T$ ($T=0$) is taken as the virtual time when the free-falling drop would first contact the free surface of the pool, assuming that both the interfaces of the drop and pool do not deform during this initial contact, and $T_0$ is an adjustable offset to represent the time delay of contact caused by viscous forces and air cushioning. In this simulation, the drop is initialized above the receiving bulk at the height $\delta =0.1d$ (see § 2.2), thus the simulation starts at $T\approx -55\ \mathrm {\mu } {\rm s}$. The coefficients $C=1.22$ and $T_0=12.7\ \mathrm {\mu }{\rm s}$ are therefore obtained by fitting the numerical measurements, and the overall good match between simulated neck motion and analytical prediction is plotted in figure 8(a).

Figure 8. Early-time dynamic behaviours of the neck region. (a) Evolution of the neck radial position $R_n$. The solid line shows the theoretical estimate using the form $r_n=R_n/R_d=C\sqrt {3(T-T_0)U_0/R_d}$, where $C=1.22$ and $T_0=12.7\ \mathrm {\mu }{\rm s}$ are obtained by fitting the numerical measurements. (b) Evolution of the ejecta angle $\theta$ measured from the vertical central slices in figure 7 (two sides), using the definition sketch proposed by Thoraval et al. (2012).

Starting from $t\approx 10\ \mathrm {\mu }{\rm s}$, an azimuthal instability is observed along the neck. Here, the diameter of the outer contact line reaches only $38\,\%$ of the drop size. Figure 9 shows the closeup view of the azimuthal undulations at the neck of connection between drop and pool at $t=10\ \mathrm {\mu }{\rm s}$, where the entire periphery of the neck is visible. Irregular undulations are captured on both sides of the neck, namely the outer liquid ejecta and the inner air sheet. Liquid fingers are initiated arbitrarily from the outer side and some of them break up immediately on their tips (white dashed circle), producing the very first generation of secondary droplets (see also figure 7b). At some locations where fingers are not found (or long wavelength), smooth ejecta and air sheets are present. According to Li et al. (2018) (their figure 8), the breakup of the ejecta from a liquid sheet into triangular teeth can be attributed to the localized disturbances within the liquid. These disturbances trigger vortex-shedding progress along the drop/pool boundary, leading to the destabilization of the ejecta base at various locations. Josserand & Zaleski (2003) estimated that the thickness of the initial ejecta should be of the order of a viscous length scale $e_j\approx \sqrt {\mu _lt/\rho _l}$, which for our present case would be ${\sim }3.15\ \mathrm {\mu }{\rm m}$. The measured ejecta thickness from the two rectangular regions in figure 9 is approximately ${\sim }5.5\ \mathrm {\mu }{\rm m}$, which is very close to the theoretical prediction. We have measured that the characteristic wavelength and amplitude of the initial undulations are around ${\sim }30\ \mathrm {\mu }{\rm m}$ and ${\sim }21\ \mathrm {\mu }{\rm m}$, which are generally longer than the thickness of the ejecta base. As the neck spreads radially, alternate vortices are shed from the base of the ejecta, thus forming a similar structure of von Kármán-type street (Thoraval et al. 2012) along the drop/target contact as shown in 7(b). Meanwhile, the oscillations of the base of the ejecta pull the local air sheet on the corner of it, entrapping occasionally some isolated bubbles/bubble arcs on both/alternate sides of the base (green arrow in figure 6).

Figure 9. Irregular azimuthal undulations on the neck region between the drop and the pool $10\ \mathrm {\mu }{\rm s}$ after first contact. The central irregular plate is the entrapped air disc. The white circle indicates the early-time breakups of ejecta fingers. The scale bar is $100\ \mathrm {\mu }{\rm m}$ long.

After the initial small amplitude oscillations combined with the ‘complete’ disintegration of the ejecta jets, a coherent liquid sheet eventually emerges along the contact line until it impacts the drop surface (see $t=32\ \mathrm {\mu }{\rm s}$ in figure 6), entrapping a larger bubble ring along the neck, as evidenced at $t=50\ \mathrm {\mu }{\rm s}$ in figures 6 and 7(c). Figure 8(b) plots the time evolution of the ejecta angle $\theta$ measured from vertical central slices in figure 7, based on the definition sketch proposed by Thoraval et al. (2012). The simulation shows that $\theta$ has an overall linear increase at the beginning and drops sharply at the moment of ‘sheet impingement’, similar to the numerical measurements for cases at $Re\approx 6000$ in Thoraval et al. (2012). Furthermore, at higher $Re$, our simulation demonstrates that the angle of the ejecta base undergoes multiple reversals during the very early stages of the ejecta jet formation, indicating repeated interactions between the ejecta jets and the interfaces of the drop and pool. Once the new neck of connection is established, the fast-moving rim will stretch the ejecta sheet and tear it immediately into multiple liquid ‘tori’, which later break up and produce a large number of similar-sized microdroplets, as illustrated in figure 7(c) and discussed later in § 5.

4.2. Formation of bubble canopy

Figure 10 shows the internal flow field of the liquid phase at different stages, overlapped by velocity vector and pressure fields at the cross-section. The grid-based arrow is oriented by velocity field and its value is represented by colour. Within the first millisecond, the drop has not yet fully sprawled out and most of the momentum is still concentrated in the impact drop, thus generating a broad high-pressure area along the drop/pool boundary. Above the initial free surface, an outward-expanding liquid crown arises along the periphery of the cavity. As the drop moves downwards, more liquid is therefore pushed away from the pool and transported to the uprising crown. At high impact velocities, thin liquid ligaments are generated along the rim of the crown and break up on its tips by instability, producing moderate and large-scale secondary droplets (see § 5). As discussed above in § 3.2, the crown reaches a maximum horizontal radial position soon after impact ($t\approx 3$ ms) and its rim thickens and bends towards the impact axis under the effect of surface tension, while it simultaneously rises in the vertical direction. The outline of the cavity expands rapidly on the surface of the pool and overtakes the crown expansion at around $t\approx 4$ ms, which also facilitates generating an inward-directed momentum on the crown top to enclose the upper part. During the period of crown expansion, the maximum velocity is reached on the top where the liquid film is thinnest. Right before the closure ($t\approx 12$ ms), it is visible that the velocity on the crown rim points almost horizontally inwards. Flows from all directions on the liquid wall meet and interact at the instant of closing, generating a sharp pressure rise around the point of closure, which later protrudes a downward-moving jet evidenced at $t=16$ ms. Meanwhile, ligaments above the dome tangle and merge, possibly shedding several large-scale droplets on the top of the dome, as shown at $t=16$ ms.

Figure 10. Internal flows of high-speed drop impact overlapped by velocity and pressure fields. From left to right and top to bottom, the corresponding times are 1, 3, 11, 16, 24, 32, 38 and 48 ms after impact. The velocity magnitudes are scaled by the drop impact speed $U_0$. The pressure field is scaled by the initial dynamic pressure of the impact drop $P_0=(\rho _l{U_0}^2)/2$. The scale bar is 10 mm long.

It can be found in the literature that large bubble entrapment owing to drop impact onto liquid pool occurs mainly in two types of mechanisms. Firstly, under low impact energy, a vortex-induced roll jet may be formed and grows into a thick liquid tongue, which later collapses near the pool surface to entrap an air bubble (Pumphrey & Elmore 1990). The shape of the drop at the time of impact is crucial for such mechanism and it occurs almost exclusively for prolate-shaped drops (Zou et al. 2012; Wang, Kuan & Tsai 2013; Thoraval, Li & Thoroddsen 2016; Deka et al. 2017). Secondly, with sufficiently high impact energy, the thin-walled liquid crown rises up higher above the pool due to violent collision and its rim bends towards the impact axis while rising vertically, enveloping the air bubble above the target pool. The formation of such complex flow structures resulting from high-speed drop–pool collisions has been observed from time to time by experiments during the past century (Worthington 1908; Engel 1966; van de Sande et al. 1974; Bisighini et al. 2010; Murphy et al. 2015; Lherm et al. 2021), but is probably firstly solved in three dimensions by high-resolution numerical simulations in the present work. Interestingly, although driven by different mechanisms, similarities are still observed between them, such as the generation of the central jets, the reconnection of the downward-moving jet as well as the bursting of the final ‘floating bubble’. Compared with the low-energy counterparts, the rather intense interfacial deformation at high impact velocities surely introduces some new phenomena, which could accordingly influence the subsequent physical process like the generation of underwater noise induced by air bubble entrainment (Prosperetti, Crum & Pumphrey 1989; Prosperetti & Oguz 1993) and the additional source of airborne droplets caused by liquid film rupture (Resch, Darrozes & Afeti 1986; Afeti & Resch 1990; Resch & Afeti 1991). Further investigation and analysis could be initiated in the future to compare and contrast these seemingly similar dynamic patterns.

Figure 11 shows the formation of the central jet from the top of the dome. The jet is formed side biased as the flows on the upper part of the crown are not perfectly axisymmetric and they do not arrive at the merger point simultaneously, which is also consistent with one of the few available experimental observations in literature (Bisighini et al. 2010; Murphy et al. 2015; Lherm et al. 2021). High-speed jet ejections out of a liquid interface are commonly observed in many other physical processes, such as bubble bursting (Boulton-Stone & Blake 1993; Thoroddsen et al. 2009; Berny et al. 2022), Faraday waves (Hogrefe et al. 1998; Zeff et al. 2000) and cavity collapse induced by the process of solid/liquid object impact onto fluid target (Worthington & Cole 1897, 1900; Gekle & Gordillo 2010; Ray et al. 2015; Jamali et al. 2020; Kim et al. 2021). In general, all these jets are ejected as a consequence of a very large axial pressure gradient created at the jet base, which therefore can be further classified based on the way the large overpressure is created and the length scale at which pressure variations take place (Gekle & Gordillo 2010). Surface tension plays an essential role in the jet formation process.

Figure 11. Formation of the central spiral jet inside the bubble canopy. (a) Dynamics of the liquid jet at $t=24$ ms. The scale bar is 4 mm long. (b) Jet motions observed from the bottom view, showing 18, 23, 24 and 38 ms after impact. The scale bar is 1 mm long. The interface is contoured by velocity field.

4.3. Cavity contraction

Now we focus on the contraction of the structure. For drop impact at low and moderate velocities, the liquid crown collapses and generates capillary waves that move radially towards the cavity bottom along the interior of the crater, after it reaches the maximum expansion. The capillary waves then meet concentrically at the cavity bottom, forming a classic upward Worthington jet that may break up on its tips (Ray et al. 2015). At higher impact velocities, qualitatively different phenomena are observed. As shown in figure 10, during the cavity expansion stage, the direction of the velocity field is outward and upward around the enlarging cavity and the cavity continues to expand even after the closure of the upper part. The cavity depth reaches its maximum value at $t\approx 24$ ms as the cavity radius extends continuously along its edge (see also figure 4d,e), which is very different from the measurements for an even more energetic case impacting two times faster ($U_0=15.98\ {\rm m}\ {\rm s}^{-1}$) in Engel (1966) where the crown and cavity arrive at their maximum positions at approximately the same time. At the moment that the penetration depth reaches its maximum value, the flow around the cavity bottom should have reached its minimum value of essentially zero and shortly be redistributed by inertia force, as seen at $t=24$ ms in figure 10. At approximately the same time ($t\approx 25$ ms), the downward spiral jet arrives at the cavity bottom and penetrates deeply into the bulk, creating a subcavity that moves also spirally, which may potentially accelerate the re-establishment of the velocity field. By $t\approx 32$ ms, it can be observed that the velocity field has been completely reversed in the pool and a new circulation has been established. The flow around the cavity in the pool has transitioned from outward–upward expansion to inward–downward contraction. Such a change of the flow direction would result in a pressure buildup around the lower part of the cavity, thus pushing the cavity to shallow. Air bubbles are entrapped in the bulk in this process as can be seen at $t=32$, 38 and 48 ms in figures 10 and 12(a). In the last frame of figure 10, an upward jet eventually rises from the cavity floor and merges with the previous downward jet.

Figure 12. Kinematic behaviours of the entrapped large bubble. (a) Successive positions of the vertical slices for the entrapped large bubble. The time interval between each curve is 4 ms. (b) Time evolution of the vertical centroid position (left axis) and the vertical speed (right axis) of the entrapped large bubble. The bubble sinks at the expansion stage and then starts to shallow from its bottom due to the concentric axial pressure, which eventually leads to a floating air bubble above the pool surface.

Figure 12(a) draws the successive shapes of the entrapped large bubble and (b) plots the temporal motion of the bubble centroid in the vertical direction, showing the shallowing steps of the cavity. The bubble expands to a maximum position and then contracts from its bottom, eventually generating a toroidal air bubble floating at the top of the pool.

As for the final collapsing stage, the central column thickens and merges with the outer bubble wall, creating a horseshoe-shaped bubble that transforms later into a hemisphere. For the next long period of time (more than 300 ms), this toroidal bubble will be stretched and thinned under the action of surface tension and eventually ruptures due to instabilities. Subsequently, the film cap recedes along the periphery of the ‘ruptured hole’ from one side, scattering fine-scale liquid droplets from the ‘receding rim’ to the air, as shown in the experiments by Murphy et al. (2015) (their figure 3). The production of such tiny droplets from the receding films has been studied by Lhuissier & Villermaux (2012) and Dasouqi et al. (2021).

4.4. Transportation of drop liquid

In the simulation, a passive tracer field $f_p$ is initially added to the impact drop for the purpose of following the trace of the drop liquid. Experimentally speaking, this can be achieved by adding colours/dye to the liquids (Engel 1966; Thoroddsen 2002; Bisighini & Cossali 2011). The tracer field is then advected by following equation:

(4.1)\begin{equation} \frac{\partial f_p}{\partial t}+\boldsymbol{U_f}\boldsymbol{\cdot}\boldsymbol{\nabla} f_p=0. \end{equation}

Figure 13 shows the distribution of drop liquid at different stages of impact. Shortly after contact ($t=0.5$ ms), the main part of the drop is still concentrated at the bottom of the ‘bulk cavity’ and radially spreading thin film extends from its edges. For the next few milliseconds, the drop flattens and expands rapidly into a thin liquid layer coated along the interior surface of the pool. Meanwhile, liquid threads are emitted from the rim of the drop liquid and then break up into small droplets on its tips ($t=3$ ms), which is comparable to the formation of azimuthal destabilization on the retracting flattened drop edge at lower impact energy reported in Lhuissier et al. (2013). By the time of closing, it can be observed that these ‘drop threads’ meet at the closure point ($t=14$ ms) and are transported backwards to the bulk by the downward-moving jet ($t=20$ ms). At the receding stage, the thin drop film starts to propagate towards the cavity bottom due to the surface-tension capillary waves and is later penetrated deeply into the pool when the central spiral jet impinges the cavity bottom, seen at $t=32$ and 48 ms.

Figure 13. Transportation of the passive tracer for drop liquid. From left to right and top to bottom, the corresponding times are 0.5, 3, 14, 20, 32 and 48 ms after impact. The scale bar is 10 mm long. Azimuthal destabilization is captured at the edge of the drop film, which therefore produces secondary droplets from its tips. At the shallowing stage, the drop liquid recedes to the cavity bottom and is later penetrated and mixed inside the target pool by the downward-moving jet.

Different from the gas/liquid interface that can be recorded directly by the camera, the ‘virtual’ drop/pool interface is usually invisible if the fluids in the drop and the receiving pool are the same, but the estimation of the kinematic behaviours of this layer is vital for theoretical studies of crater evolution (Bisighini & Cossali 2011). As indicated in figure 14(a), the boundaries between different components can be easily differentiated here by the isosurface $f_p=0.5$. Figure 14(b) plots the temporal variations for the positions of the upper point $T_a$, lower point $T_b$ as well as the thickness of the drop tracer $T_\delta$ along the impact axis. As expected, the upper point $T_a$ moves rapidly adjoining the air, while the penetration speed of the lower point $T_b$ is greatly decelerated by the reacting flows from the target liquid, thus causing the decrease in drop thickness. It can be seen that the drop deforms significantly and its thickness decreases noticeably during an initial dimensionless time period $\tau =tU_0/d\leq 2$. For the later stages, $T_\delta$ reaches a plateau and will not experience much difference throughout the expansion stage. Slight fluctuations of $T_\delta$ can be anticipated at the contraction stage when the drop liquid recedes into the cavity bottom.

Figure 14. (a) Sketch of the drop penetration. Boundaries between different fluid components are differentiated by isosurface $f_p=0.5$. The positions of the upper point $T_a$, lower point $T_b$ and the thickness of the drop tracer $T_\delta$ along the vertical axis of symmetry are tracked. (b) Time variations of $T_a$, $T_b$ and $T_\delta$ along the axial direction. The dashed line shows the asymptotic solution proposed by Berberović et al. (2009). The solid line shows the theoretical estimation of the penetration depth proposed by Bisighini et al. (2010). The dimensionless time $tU_0/d=2$ is indicated by the vertical dotted line.

Following previous investigations in the literature, a critical dimensionless time $tU_0/d\approx 2$ is usually used to subdivide the evolution of drop impact into two phases (Fedorchenko & Wang 2004). During the first phase ($tU_0/d\leq 2$), the drop deforms and extends above the bulk cavity, and the interfaces of air/drop and drop/pool are clearcut. The penetration speed of the drop/pool interface during this time can be approximated as half of the impact speed $U_p\approx U_0/2$, which is well known from the penetration mechanics (Birkhoff et al. 1948; Yarin, Rubin & Roisman 1995) and has been previously applied for the analytical study of crater evolution (Fedorchenko & Wang 2004; Berberović et al. 2009). At times $tU_0/d>2$, the flow effects in the thin drop layer are negligibly small and the cavity expansion thus can be approximated as the shape of the drop/pool interface. Berberović et al. (2009) developed a theoretical approach to estimate the penetration depth of drop impact $T_d$ based only on the linear momentum balance of the liquid around the cavity and gave an asymptotic solution as $T_d=2^{-4/5}(5t-6)^{2/5}$ for $tU_0/d>2$ at high $Fr$, $We$ and $Re$ numbers. The predicted results using the above asymptotic formula are shown in figure 14(b) (dashed line) and a reasonable agreement is found with the simulated results. Since the effects of surface tension, viscosity and gravity are neglected in this theory, it cannot predict accurately the later stages of impact where the deformation of cavity shape becomes significant due to gravity and capillary waves. Bisighini et al. (2010) proposed a theoretical model based on the potential flow theory that accounted for the effects of inertia, gravity, viscosity and surface tension for sufficiently high Reynolds and Weber numbers, using the combination of the sphere expansion and its translation along impact axis. A system of ordinary differential equations is obtained and numerically solved using initial conditions

(4.2)$$\begin{gather} \ddot{\alpha}=-\frac{3}{2}\frac{\dot{\alpha}^2}{\alpha}-\frac{2}{\alpha^2We}-\frac{1}{Fr}\frac{\zeta}{\alpha}+\frac{7}{4}\frac{\dot{\zeta}^2}{\alpha}-\frac{4\dot{\alpha}}{\alpha^2Re}, \end{gather}$$

(4.3)$$\begin{gather}\ddot{\zeta}=-3\frac{\dot{\alpha}\dot{\zeta}}{\alpha}-\frac{9}{2}\frac{\dot{\zeta}^2}{\alpha}-\frac{2}{Fr}-\frac{12\dot{\zeta}}{\alpha^2Re}, \end{gather}$$

where $\alpha$ and $\zeta$ donate the dimensionless crater radius and axial coordinate of the centre of the sphere and the dimensionless penetration depth is expressed as $\alpha +\zeta$. As explained by Bisighini et al. (2010), the initial conditions can be obtained from the initial phase ($tU_0/d\leq 2$) using the forms $\dot {\alpha }\approx 0.17$, $\alpha \approx \alpha +0.17\tau$, $\dot {\zeta }\approx 0.27$, $\zeta \approx -\alpha _0+0.17\tau$ and the dimensionless width of the cavity can be estimated using the geometrical conditions $W=2\sqrt {\alpha ^2-\zeta ^2}\approx 2\sqrt {(\alpha _0+0.17\tau )^2-(0.27\tau -\alpha _0)^2}$. By fitting the simulated bulk cavity width using the least-mean-square method, the constant is obtained as $\alpha _0=0.79$ in the present case, thus the initial conditions for (4.2) and (4.3) are $\alpha (2)=1.13$, $\dot {\alpha }(2)=0.17$, $\zeta (2)=-0.25$ and $\dot {\zeta }(2)=0.27$. As shown in the solid line of figure 14(b), the temporal variation of the predicted depth of the drop/pool interface agrees very well with our numerical results during the expansion stage. As for the retraction stage, discrepancies between the theoretical prediction and the simulation/experiment become more pronounced as demonstrated in figure 4(e), which can be explained as the fact that the shape of the cavity does not follow the spherical expansion anymore due to the impingement of central spiral jet and the propagation of capillary waves, so the model is invalid for the retraction phase.

5. Airborne droplets

The production of secondary droplets and their distribution in the process of normal (perpendicular) (Okawa, Shiraishi & Mori 2006; Guildenbecher et al. 2014; Li, Zhang & Liu 2019; Wu, Wang & Zhao 2020) and oblique (Okawa et al. 2008; Liu 2018) drop impact on the liquid surface has been discussed in some previous studies, focusing mostly on the fairly ‘regular’ splashing regimes under relatively low $Re$ and $We$. At high-energy conditions, complex splashing processes have been experimentally observed (Thoroddsen et al. 2011; Thoraval et al. 2012; Zhang et al. 2012; Thoraval et al. 2013), suggesting that the creation of these tiny droplets might be associated with more complicated ‘irregular’ breakup of liquid sheet. According to Deegan et al. (2007), at least three sources of secondary droplets can be identified during the impact: (i) prompt instability of ejecta sheet occurring immediately after contact, which produces very small droplets, (ii) early-time rim instability of ejecta sheet that produces medium-sized droplets and (iii) rim instability of crown that produces large droplets from ligaments. However, these different mechanisms are typically interdependent and the earlier ones influence the later ones, which therefore further complicates the characterization of airborne droplets. The understanding of the governing mechanisms of droplet ejection and their population remain insufficient, and the statistical data on droplet production are very limited in the literature.

5.1. Source of secondary droplet

Figure 15 qualitatively illustrates some primary mechanisms of droplet production during the impact captured by simulation. The left panel shows stages and locations where droplets are generated and the right panel demonstrates the mesh structures on the overlaid section. As shown in figure 15(a), a great number of microdroplets are sprayed immediately after contact due to the very early breakups of ejecta film, so-called the ‘prompt splash’ (Deegan et al. 2007; Marcotte et al. 2019). After the intricate early splash, a ‘coherent’ liquid sheet rises up around the contact line and forms the cylindrical crown, emanating thin liquid ligaments on its rims from a fairly regular distance. A sustained droplet ejection, so-called the ‘crown splash’, is therefore observed. The size of droplets produced at this stage is greatly influenced by the thickness of ligaments and generally increases with time. Figure 15(b) shows the crown fragmentation at $t=0.65$ ms. From the right panel, it can be clearly seen that the droplets produced at the present time instant are sufficiently larger than the minimum cell scale. Figure 15(c) shows the production of very small droplets from ‘secondary impact’ caused by previous generations. When the first-born droplets fall back and impinge on the pool, it may produce another splash, which is generally partially resolved. The radii of these smallest droplets are represented approximately by the smallest size of the cell (right panel). Besides the above sources, very small child drops can be also ejected from bubble bursting, which is usually not fully resolved as same as figure 15(c). Fewer larger droplets are also observed later from the downward-moving central jet after the closure of the upper crown.

Figure 15. Different mechanisms of droplet production at different stages of impact. The images are shown under different magnifications: (a) $t=0.03$ ms, the ‘prompt splash’ that occurs at the very early time of impact near the neck region due to irregular rupture/breakup of ejecta; (b) $t=0.65$ ms, the sustained ‘crown splash’ due to the breakup of thin ligaments on the top of the crown rim; (c) partially resolved tiny droplets near the pool surface produced by secondary impact and bubble bursting. The left panel shows the locations of splashes and the right panel demonstrates the mesh structures nearby.

5.2. Droplet statistics

We now discuss the statistics of droplets. Here, we only analyse the droplet information captured with higher resolutions ($L_{max}=13$ and 14) for the first 4 ms after impact. The time variation of the total number of droplets is plotted in figure 16(a). It should be noted that a short-time disturbance is presented at $t\approx 0.2$ ms on the curve, which might be explained as the loss of smallest droplets within the extra refinement layer at $L_{max}=14$ (see § 2.3), as it happens at approximately the same time when many droplets reach this height and the extra refinement layer is removed. Shortly after contact, a large number of very fine droplets are scattered from the emerge-ruptured ejecta (corresponding to figure 15a), which is reflected in figure 16(a) for the sharp increase of droplet count at the beginning. The maximum droplet count is found 0.65 ms after impact with around 4340 droplet population, and the dynamics around this time is qualitatively shown in figure 15(b). A change of the decreasing slope can be found at $t\approx 1.7$ ms in figure 16(a), which indicates the timing when lots of droplets start to exit the field of view and are removed from the computational domain.

Figure 16. (a) Temporal evolution of the total number of secondary droplets. The first vertical dotted line indicates the time point at $t=0.2$ ms and the second dotted line indicates the maximum droplet count at $t\approx 0.65$ ms. (b) Temporal evolution of the total mass of secondary droplets $M_d$, scaled by the mass of the impact drop $M_0$.

Figure 16(b) plots the time evolution of the mass ratio of the total secondary droplets ($M_d$) to the impact drop ($M_0$). An overall increasing trend of the mass transfer from pool to air is found. The number of droplets starts to decrease after the peak ($t\approx 0.65$ ms) but the total mass of droplets keeps increasing, which also reveals the fact that the droplets produced at the later ‘crown splash’ stage are in much larger scales than at the early ‘prompt splash’.

The sensitivity of droplet statistics to changes in mesh resolution during sheet fragmentation has been discussed in Appendices A.1 and A.3, showing that large droplets are approximately numerically converged but small droplets close to the minimum cell size ($S_d\approx 2\varDelta \sim 4\varDelta$) are generally overpredicted. Similar effects have also been captured by other high-resolution simulations involving droplet production using BASILISK code: the mesh study in Pairetti et al. (2020) showed that the most frequent droplet diameter during primary atomization is always near $2\varDelta$ for all three proposed mesh resolutions; Mostert, Popinet & Deike (2022) examined two different maximum levels on droplet spray in the action of breaking waves, where rapid increases of droplet count for droplet radii in the $\varDelta \sim 2\varDelta$ range are always observed (figures 8 and 9 in their supplementary materials). These results suggest that 4 cells per droplet diameter are essential to obtain numerical convergence, and droplets less than 4 cells per diameter are not reliable and should be considered as under-resolved structures.

Figure 17(a) plots the contour map of droplet size distribution per bin size $\Delta r$ during $0\sim 4$ ms after impact. The equivalent mean diameter of each droplet $S_d$ is calculated as a sphere using the integrated volume fraction of liquid phase. Droplet sizes at $S_d=2\varDelta$ (lower) and $S_d=4\varDelta$ (upper) are indicated by horizontal red dashed lines. Figure 17(b) shows various instantaneous time slices and time-averaging windows of droplet size distributions. The shift of the small-sized ribbon can be clearly seen between stage $S1$ and stage $S2$ when $L_{max}$ is decreased due to the effect of minimum cell size.

Figure 17. (a) Temporal contour of droplet size distribution during $0\sim 4$ ms of impact. The vertical dotted line shows the time point at $t=0.2$ ms, where $L_{max}$ changes from 14 to 13. The horizontal dashed lines indicate the length scales at $S_d=2\varDelta$ (lower) and $S_d=4\varDelta$ (upper). (b) Droplet size distributions at different time slices in (a). The filled areas show numerical results at $t=0.03$ ms right after contact and at $t=0.2$ ms where crown splash starts, calculated at $L_{max}=14$. Time-averaged droplet size distribution is calculated using 100 time slices for time windows $t\in [0.2, 2]$, $[2, 3]$ and $[3, 4]$ ms. Most $S_d<60\ \mathrm {\mu }{\rm m}$ droplets are produced from the early-time splash at $t<200\ \mathrm {\mu }{\rm s}$. Large droplets ($S_d>100\ \mathrm {\mu }{\rm m}$) are only generated from ligament breakups at the ‘crown splash’ stage. Bi-model distribution of droplet size is found 2 ms after impact in the domain.

Initially, a lot of very small droplets ranging in $S_d\approx 15\sim 30\ \mathrm {\mu }{\rm m}$ are generated immediately upon impact as shown at $t=0.03$ ms in figure 17(b), which make up only around 30 % of the total population, suggesting that the initial smallest droplets should be treated with caution. In the experiments of Murphy et al. (2015), a ring of fine spray containing mostly $S_d\approx 6\sim 19\ \mathrm {\mu }{\rm m}$ droplets was also recorded at the instant of impact ($6.22\ \mathrm {\mu }{\rm m}\ {\rm pixel}^{-1}$), in the same magnitude as our numerical results. However, the existence of smaller droplets is not able to be determined at this moment owing to the limited resolutions for both numerical and experimental approaches. Droplets larger than $100 \ \mathrm {\mu }{\rm m}$ are very few at $t\leq 0.2$ ms, indicating that large droplets are only produced by the fragmentation of rim ligaments at the ‘crown splash’ stage.

After the initial splashing stage, an ‘coherent’ liquid sheet rises cylindrically from the contact line and disseminates secondary droplets from rim ligaments as shown in figure 15(b). During stage $S2$, the newly detached droplets increase gradually in size as ligaments thicken and merge, corroborating the gradually increased top boundary of the contour in figure 17(a). The averaged droplet distributions over the time windows $t\in [0.2, 2]$, $[2, 3]$, and $[3, 4]$ ms are plotted in figure 17(b). Comparing size distribution at $t\in [0.2, 2]$ ms with $t=0.2$ ms, very similar profiles can be found in the $20\sim 60\ \mathrm {\mu }{\rm m}$ range, reflecting that most of them are generated from the initial ‘prompt splash’ or the very early breakup of ‘crown splash’ at $t<0.2$ ms. During $t=0.2\sim 2$ ms, droplets ranging in $80\sim 250\ \mathrm {\mu }{\rm m}$ increase significantly due to the mechanism of crown fragmentation as shown in figure 22(d). Starting from $t\approx 2$ ms, the earlier droplets have moved far from the impact region and subsequently exit the observation domain while the newly pinched larger droplets are still near the crown rim, therefore forming bimodal distributions of droplet size as shown in figure 17. As the impact advances, the upper peak broadens and its centre shifts to larger sizes, revealing the fact that both the escaped and newly generated droplets increase in size. However, for the small-sized peak, it is still hard to determine its centre and amplitude due to the limitation of the minimum cell size. The bimodal size distribution of droplets has been also reported by some other studies involving liquid sheet fragmentation (Afeti & Resch 1990; Villermaux & Bossa 2009, 2011; Villermaux, Pistre & Lhuissier 2013).

Lastly, we would like to focus on those droplets that tend to fall back to the liquid bulk, which is vital for the estimation of mass transfer through the air–sea interface. In practice, those droplets that move along the impact direction (downwards) and whose centroid are located less than $100\ \mathrm {\mu }{\rm m}$ to the free surface of the liquid pool are selected and assumed as the ones that will re-join the pool. Figure 18(a) shows the temporal contour of the size distribution for the ‘re-merging’ droplets and (b) shows its count in time. Initially, a big part of very small droplets vibrate near the free surface and tend to re-join the bulk shortly after generation. At the very early stage of ‘crown splash’, very thin ligaments stretch out from the downward-bending crown and send out tiny droplets towards the target liquid as demonstrated in figure 3(b), which thereby contribute to a second primary peak of the ‘re-merging’ droplets. Despite the fact that a majority of these ‘re-merging’ droplets are produced from the partially resolved structures ($S_d<4\varDelta$), this information still provides valuable insights into the behaviour of secondary droplets and their interactions with the receiving pool, corroborating the chaotic interfacial activities demonstrated in figures 6 and 7. As the impact advances, the direction of the liquid ligaments transitions quickly from horizontal–outward to vertical–upward and the droplets are generally pinched off upwards in larger angles accordingly, thus the number of the ‘re-merging’ droplets becomes insignificant in the domain. Finally, it is worth noting that the presence of wind, not accounted for in the present study, may significantly advect the smallest droplets and delay or prevent them from rejoining the bulk.

Figure 18. Statistics of droplets that tend to re-merge with the liquid bulk. (a) Temporal contour of the droplet size distribution for the ‘re-merging’ droplets. The vertical dotted line shows the time point at $t=0.2$ ms, where $L_{max}$ changes from 14 to 13. The horizontal dashed lines indicate the length scales at $S_d=2\varDelta$ (lower) and $S_d=4\varDelta$ (upper). (b) Evolution of the ‘re-merging’ droplet count with time.

6. Conclusion

In this work, the high-energy splash of drop impact onto a deep volume of the same liquid pool is investigated by performing high-resolution direct numerical simulations in three dimensions. The calculations are conducted in the exact same configurations of the control case previously studied in Murphy et al. (2015), in order to perform detailed comparisons with experimental data and prepare an in-depth analysis of the splashing dynamics. The BASILISK open-source solver, which combines a VOF description of the gas/liquid interfaces and an adaptive octree grid refinement, is employed to simulate the process of drop impact. Qualitative and quantitative comparisons between numerical and experimental results have been conducted in terms of the morphological behaviours of splash, kinematics of crown and cavity as well as the distributions of secondary droplets through a particular field of view, which efficiently validated the present numerical strategies and therefore enabled the discussion of the internal mechanisms for high-speed drop impact afterwards.

Following the experimental observations of Murphy et al. (2015), we have performed a detailed investigation on flow physics and splashing behaviours of high-speed drop impact, serving as an important supplementary study for this issue. Firstly, the very early instantaneous motions in the vicinity of the neck region are discussed under sufficient time resolution. We observe the existence of two different mechanisms of air entrapment in the neck of connection, namely the entrapment of axisymmetric bubble rings driven by high localized pressure at $R_n/R_d<35\,\%$ (Oguz & Prosperetti 1989) and the entrapment of isolated bubbles/bubble rings due to unstable oscillations of ejecta base (Weiss & Yarin 1999). Moreover, the simulation successfully captures the initialization of irregular azimuthal undulations along the outer edge of the neck, which thereby breaks the axisymmetry of the motions in microscopic scales. Thereafter, detailed information on the internal flows, such as velocity, pressure and passive tracer fields, is extracted to explain the corresponding physical phenomena observed in experiments. We show that azimuthal destabilization occurs on the edges of the flattening drop, breaking up on its tips and participating in the production of child droplets. These ‘liquid threads’, emitted from the initial impact drop, grow together with the uprising crown and meet upon closure, and they are eventually transported backwards to the bulk along with the penetration of the central spiral jet. Lastly, we present the statistics of airborne droplets produced by high-energy drop–pool collisions. The results show that a great number of very fine microdroplets are produced immediately after contact by irregular ‘prompt splash’ within the time scale of $t<200\ \mathrm {\mu }{\rm s}$, composing the most populated small and moderate sizes. The earliest tiny droplets are sprayed underneath the drop and vibrate near the pool surface, suggesting the great possibility that a big part of them may return to the liquid bulk shortly after birth. Large droplets ($S_d>100\ \mathrm {\mu }{\rm m}$) are only observed from the fragmentation of rim ligaments during the ‘crown splash’ stage, which therefore forms a gradually increasing but narrowing distribution ribbon, reflecting the merging and thickening activities of ligaments on the rim of the crown. Finally, the bimodal size distribution of secondary droplets is found within the entire domain.

We are aware that the splashing dynamics induced by the most energetic regime of drop impact is far more complicated than has been mentioned here. Further analysis of multi-scale flow physics needs to be conducted under sufficient spatial and temporal resolutions in the future. Details of bubble ring entrapment are still not well understood, as well as its motions and breakups. The physical mechanism that is responsible for the early azimuthal instability remains unknown. Recent experimental observations (Thoraval et al. 2013; Li et al. 2018) have suggested that this early dynamics may be greatly affected by the intricate vortical motions and 3-D instabilities, which certainly added more complexities to the analysis. For drop impact on the same liquid pool, large parameter space of $Re$ and $We$ need to be studied to determine the critical conditions of large bubble entrapment (bubble canopy regime), as well as to explore its influence on the production of secondary droplets and closure event (closure time, closure height, large bubble volume, floating bubble and its burst). The formation and motion of the central spiral jet also need to be analysed in detail and compared with various types of liquid column jets. Statistics of droplets and bubbles need to be collected under various impact conditions to form a more comprehensive database, which could therefore reflect directly the mass/momentum transfer through gas–liquid interfaces as well as enlighten various applications where this process is involved (oil dispersant, spray cooling, metallurgy, disease transmission, etc.).