1. Introduction
The impact of raindrops on a deep liquid pool has been extensively studied since the initial works of Worthington (Reference Worthington1883, Reference Worthington1908). Various behaviours were obtained, depending primarily on miscible or immiscible characteristics of the two liquids (Lhuissier et al. Reference Lhuissier, Sun, Prosperetti and Lohse2013; Castillo-Orozco et al. Reference Castillo-Orozco, Davanlou, Choudhury and Kumar2016), but also on their difference of density (Thomson & Newall Reference Thomson and Newall1886; Manzello & Yang Reference Manzello and Yang2002; Lherm et al. Reference Lherm, Deguen, Alboussière and Landeau2022; Villermaux Reference Villermaux2022), the depth of the receiving volume (Macklin & Hobbs Reference Macklin and Hobbs1969; Wang & Chen Reference Wang and Chen2000; Fedorchenko & Wang Reference Fedorchenko and Wang2004) and the angle of the impact (Zhbankova & Kolpakov Reference Zhbankova and Kolpakov1990; Okawa, Shiraishi & Mori Reference Okawa, Shiraishi and Mori2008; Gielen et al. Reference Gielen, Sleutel, Benschop, Riepen, Voronina, Visser, Lohse, Snoeijer, Versluis and Gelderblom2017; Liu Reference Liu2018), 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 (Reference Schotland1960) 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 (Reference Engel1966, Reference Engel1967) 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 Reference van de Sande, Smith and van Oord1974; Rodriguez & Mesler Reference Rodriguez and Mesler1985; Pumphrey, Crum & Bjo/Rno/ Reference Pumphrey, Crum and Bjo/Rno/1989; Rein Reference Rein1993, Reference Rein1996; Leng Reference Leng2001; Manzello & Yang Reference Manzello and Yang2002; Yarin Reference Yarin2006; Fedorchenko & Wang Reference Fedorchenko and Wang2004; Berberović et al. Reference Berberović, van Hinsberg, Jakirlić, Roisman and Tropea2009; Bisighini et al. Reference Bisighini, Cossali, Tropea and Roisman2010; Murphy et al. Reference Murphy, Li, d'Albignac, Morra and Katz2015), 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 Reference Weiss and Yarin1999; Thoroddsen Reference Thoroddsen2002; Josserand & Zaleski Reference Josserand and Zaleski2003). This ejecta stretches radially and may break up into ligaments and droplets, resulting in the splash of fine sprays (Deegan, Brunet & Eggers Reference Deegan, Brunet and Eggers2007; Thoroddsen et al. Reference Thoroddsen, Thoraval, Takehara and Etoh2011). 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. (Reference Zhang, Toole, Fezzaa and Deegan2012) 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. (Reference Agbaglah, Thoraval, Thoroddsen, Zhang, Fezzaa and Deegan2015) 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 (Reference Thoroddsen2002).
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. Reference Agbaglah, Thoraval, Thoroddsen, Zhang, Fezzaa and Deegan2015) or alternate signs (Castrejón-Pita, Castrejón-Pita & Hutchings Reference Castrejón-Pita, Castrejón-Pita and Hutchings2012; Thoraval et al. Reference Thoraval, Takehara, Etoh, Popinet, Ray, Josserand, Zaleski and Thoroddsen2012) 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. Reference Thoraval, Takehara, Etoh and Thoroddsen2013; Agbaglah et al. Reference Agbaglah, Thoraval, Thoroddsen, Zhang, Fezzaa and Deegan2015; Li et al. Reference Li, Thoraval, Marston and Thoroddsen2018), but the mechanisms are still not elucidated. Bottom visualizations of drop impact on thin liquid films (Thoraval et al. Reference Thoraval, Takehara, Etoh and Thoroddsen2013) 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. (Reference Li, Thoraval, Marston and Thoroddsen2018), 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. Reference Thoraval, Takehara, Etoh, Popinet, Ray, Josserand, Zaleski and Thoroddsen2012; Agbaglah et al. Reference Agbaglah, Thoraval, Thoroddsen, Zhang, Fezzaa and Deegan2015). 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. (Reference Murphy, Li, d'Albignac, Morra and Katz2015), 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.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20231213144154107-0623:S0022112023007012:S0022112023007012_fig1.png?pub-status=live)
Figure 1. Various splashing behaviours induced by drop impact on a miscible liquid pool reported in the literature, reproduced from Murphy et al. (Reference Murphy, Li, d'Albignac, Morra and Katz2015) 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 Reference Gunn and Kinzer1949); 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. Reference Pan, Cheng, Chou and Wang2008; Bisighini et al. Reference Bisighini, Cossali, Tropea and Roisman2010; Murphy et al. Reference Murphy, Li, d'Albignac, Morra and Katz2015; Sochan et al. Reference Sochan, Beczek, Mazur, Ryżak and Bieganowski2018). 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 Reference Liang and Mudawar2016; Breitenbach, Roisman & Tropea Reference Breitenbach, Roisman and Tropea2018; Benther et al. Reference Benther, Pelaez-Restrepo, Stanley and Rosengarten2021; Lohse Reference Lohse2022) and environmental issues (Murphy et al. Reference Murphy, Li, d'Albignac, Morra and Katz2015; Castillo-Orozco et al. Reference Castillo-Orozco, Davanlou, Choudhury and Kumar2016; Bourouiba Reference Bourouiba2021; Dasouqi, Yeom & Murphy Reference Dasouqi, Yeom and Murphy2021). Indeed, it was estimated that such cases could produce at least 2000 very fine droplets from its initial intrusion (Murphy et al. Reference Murphy, Li, d'Albignac, Morra and Katz2015) and more than 900 droplets larger than
$100\ \mathrm {\mu }{\rm m}$ from the later crown expansion (Blanchard & Woodcock Reference Blanchard and Woodcock1957), 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 Reference Harlow and Shannon1967; Ferreira, Larock & Singer Reference Ferreira, Larock and Singer1985) 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 (Reference Oguz and Prosperetti1989, Reference Oguz and Prosperetti1990) included the surface-tension effects and Morton, Rudman & Jong-Leng (Reference Morton, Rudman and Jong-Leng2000) 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. Reference Morton, Rudman and Jong-Leng2000; Josserand & Zaleski Reference Josserand and Zaleski2003; Berberović et al. Reference Berberović, van Hinsberg, Jakirlić, Roisman and Tropea2009; Thoraval et al. Reference Thoraval, Takehara, Etoh, Popinet, Ray, Josserand, Zaleski and Thoroddsen2012; Ervik et al. Reference Ervik, Hellesø, Munkejord and Müller2014; Agbaglah et al. Reference Agbaglah, Thoraval, Thoroddsen, Zhang, Fezzaa and Deegan2015; Ray, Biswas & Sharma Reference Ray, Biswas and Sharma2015; Deka et al. Reference Deka, Ray, Biswas, Dalal, Tsai and Wang2017; Fudge, Cimpeanu & Castrejón-Pita Reference Fudge, Cimpeanu and Castrejón-Pita2021; Alventosa, Cimpeanu & Harris Reference Alventosa, Cimpeanu and Harris2023; Fudge et al. Reference Fudge, Cimpeanu, Antkowiak, Castrejón-Pita and Castrejón-Pita2023) 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 Reference Popinet2003; Nikolopoulos, Theodorakakos & Bergeles Reference Nikolopoulos, Theodorakakos and Bergeles2007), 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 Reference Rieber and Frohn1999; Brambilla & Guardone Reference Brambilla and Guardone2015; Cheng & Lou Reference Cheng and Lou2015; Guo & Lian Reference Guo and Lian2017). 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. Reference Wu, Zhang, Xiao and Ni2021). Only a limited number of full 3-D simulations (Cheng & Lou Reference Cheng and Lou2015; Chen et al. Reference Chen, Shu, Wang and Yang2020; Constante-Amores et al. Reference Constante-Amores, Kahouadji, Shin, Chergui, Juric, Castrejón-Pita, Matar and Castrejón-Pita2023) 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 Reference Thoroddsen2002; Thoroddsen et al. Reference Thoroddsen, Thoraval, Takehara and Etoh2011; Castrejón-Pita et al. Reference Castrejón-Pita, Castrejón-Pita and Hutchings2012; Thoraval et al. Reference Thoraval, Takehara, Etoh and Thoroddsen2013; Li et al. Reference Li, Thoraval, Marston and Thoroddsen2018), 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. Reference Thoraval, Takehara, Etoh, Popinet, Ray, Josserand, Zaleski and Thoroddsen2012, Reference Thoraval, Takehara, Etoh and Thoroddsen2013), 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. (Reference Murphy, Li, d'Albignac, Morra and Katz2015) 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. (Reference Murphy, Li, d'Albignac, Morra and Katz2015), 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
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20231213144154107-0623:S0022112023007012:S0022112023007012_eqn1.png?pub-status=live)
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20231213144154107-0623:S0022112023007012:S0022112023007012_eqn2.png?pub-status=live)
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20231213144154107-0623:S0022112023007012:S0022112023007012_eqn3.png?pub-status=live)
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 (Reference Popinet2015) and available for use and development under a free software GPL license (Popinet & collaborators Reference Popinet2013–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 (Reference Fuster and Popinet2018) 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 (Reference Weymouth and Yue2010), 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 Reference Rudman1998). 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 Reference Bell, Colella and Glaz1989), 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. Reference Francois, Cummins, Dendy, Kothe, Sicilian and Williams2006), which is based on the continuum-surface-force approach originally proposed by Brackbill, Kothe & Zemach (Reference Brackbill, Kothe and Zemach1992). In addition, a second-order accurate calculation of the curvature is performed, using the height-function technique developed by Popinet (Reference Popinet2009). More detailed descriptions of the numerical schemes can be found in Fuster & Popinet (Reference Fuster and Popinet2018), Popinet (Reference Popinet2018), Pairetti et al. (Reference Pairetti, Damian, Nigro, Popinet and Zaleski2020) and Zhang, Popinet & Ling (Reference Zhang, Popinet and Ling2020).
Cubic finite volumes organized hierarchically in octree are used for space discretization. The octree structure has been developed initially for image processing (Samet Reference Samet1990) and later applied to computational fluid dynamics (CFD) (Khokhlov Reference Khokhlov1998) and multiphase flows (Popinet Reference Popinet2003). The basic organization is the following, all details can be found in Popinet (Reference Popinet2003, Reference Popinet2015): 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. Reference van Hooft, Popinet, van Heerwaarden, van der Linden, de Roode and van de Wiel2018) 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. Reference Thoraval, Takehara, Etoh, Popinet, Ray, Josserand, Zaleski and Thoroddsen2012; Agbaglah et al. Reference Agbaglah, Thoraval, Thoroddsen, Zhang, Fezzaa and Deegan2015; Fudge et al. Reference Fudge, Cimpeanu and Castrejón-Pita2021; Wu et al. Reference Wu, Zhang, Xiao and Ni2021; Sanjay et al. Reference Sanjay, Lakshman, Chantelot, Snoeijer and Lohse2023; Fudge et al. Reference Fudge, Cimpeanu, Antkowiak, Castrejón-Pita and Castrejón-Pita2023), and the great superiority of parallel capacity and computational efficiency of the solver were discussed in Wu et al. (Reference Wu, Zhang, Xiao and Ni2021). Furthermore, the capabilities of the BASILISK solver have been extensively validated on various problems of multi-phase flows in Popinet & collaborators Reference Popinet(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. (Reference Murphy, Li, d'Albignac, Morra and Katz2015) 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$.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20231213144154107-0623:S0022112023007012:S0022112023007012_fig2.png?pub-status=live)
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. (Reference Murphy, Li, d'Albignac, Morra and Katz2015).
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 Reference Engel1966; Berberović et al. Reference Berberović, van Hinsberg, Jakirlić, Roisman and Tropea2009; Bisighini et al. Reference Bisighini, Cossali, Tropea and Roisman2010; Ray et al. Reference Ray, Biswas and Sharma2015), 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.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20231213144154107-0623:S0022112023007012:S0022112023007012_fig3.png?pub-status=live)
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 (Reference Thoroddsen2002) (their figure 2c), also known as ‘prompt splash’ (Deegan et al. Reference Deegan, Brunet and Eggers2007). 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.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20231213144154107-0623:S0022112023007012:S0022112023007012_fig4.png?pub-status=live)
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. (Reference Bisighini, Cossali, Tropea and Roisman2010). 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 Reference Engel1967; van de Sande et al. Reference van de Sande, Smith and van Oord1974; Prosperetti & Oguz Reference Prosperetti and Oguz1993; Leng Reference Leng2001; Berberović et al. Reference Berberović, van Hinsberg, Jakirlić, Roisman and Tropea2009; Bisighini et al. Reference Bisighini, Cossali, Tropea and Roisman2010; Jain et al. Reference Jain, Jalaal, Lohse and van der Meer2019). 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. Reference Bisighini, Cossali, Tropea and Roisman2010) 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. (Reference Murphy, Li, d'Albignac, Morra and Katz2015), 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. (Reference Murphy, Li, d'Albignac, Morra and Katz2015). 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. Reference Murphy, Li, d'Albignac, Morra and Katz2015). 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.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20231213144154107-0623:S0022112023007012:S0022112023007012_fig5.png?pub-status=live)
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. (Reference Murphy, Li, d'Albignac, Morra and Katz2015).
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. Reference Thoroddsen, Thoraval, Takehara and Etoh2012; Tran et al. Reference Tran, de Maleprade, Sun and Lohse2013) or contracts into several bubbles along the central line (Thoroddsen, Etoh & Takehara Reference Thoroddsen, Etoh and Takehara2003; Jian et al. Reference Jian, Channa, Kherbeche, Chizari, Thoroddsen and Thoraval2020). 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 Reference Weiss and Yarin1999; Thoroddsen Reference Thoroddsen2002; Josserand & Zaleski Reference Josserand and Zaleski2003; Howison et al. Reference Howison, Ockendon, Oliver, Purvis and Smith2005; Deegan et al. Reference Deegan, Brunet and Eggers2007; Marcotte et al. Reference Marcotte, Michon, Séon and Josserand2019). 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. Reference Thoraval, Takehara, Etoh and Thoroddsen2013; Li et al. Reference Li, Thoraval, Marston and Thoroddsen2018), 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 Reference Thoroddsen2002; Castrejón-Pita et al. Reference Castrejón-Pita, Castrejón-Pita and Hutchings2012; Thoraval et al. Reference Thoraval, Takehara, Etoh, Popinet, Ray, Josserand, Zaleski and Thoroddsen2012). 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.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20231213144154107-0623:S0022112023007012:S0022112023007012_fig6.png?pub-status=live)
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. (Reference Hendrix, Bouwhuis, van der Meer, Lohse and Snoeijer2016) 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 (Reference Hicks and Purvis2010) and Hicks et al. (Reference Hicks, Ermanyuk, Gavrilov and Purvis2012), where half of the impact speed
$U_0/2$ was suggested later by Tran et al. (Reference Tran, de Maleprade, Sun and Lohse2013) and Hendrix et al. (Reference Hendrix, Bouwhuis, van der Meer, Lohse and Snoeijer2016) 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 (Reference Scolan and Korobkin2001) and Thoraval et al. (Reference Thoraval, Takehara, Etoh and Thoroddsen2013).
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20231213144154107-0623:S0022112023007012:S0022112023007012_fig7.png?pub-status=live)
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 (Reference Weiss and Yarin1999). 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 (Reference Oguz and Prosperetti1989) 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 (Reference Eggers, Lister and Stone1999) 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 (Reference Duchemin, Eggers and Josserand2003) 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 (Reference Chandrasekhar2013), 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 Reference Rosenthal1962; de Hoog & Lekkerkerker Reference de Hoog and Lekkerkerker2001; Ashmore & Stone Reference Ashmore and Stone2004; Eggers & Villermaux Reference Eggers and Villermaux2008). Similar effects have also been reported in the experiments of Thoraval et al. (Reference Thoraval, Takehara, Etoh and Thoroddsen2013), 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 Reference Rioboo, Marengo and Tropea2002; Josserand & Zaleski Reference Josserand and Zaleski2003), but it violates the continuity equation. By incorporating Wagner's theory (Wagner Reference Wagner1932), the radial motion of the neck for drop impact on a solid surface has been deduced as
$r_n=\sqrt {3\tau }$ (Riboux & Gordillo Reference Riboux and Gordillo2014; Philippi, Lagrée & Antkowiak Reference Philippi, Lagrée and Antkowiak2016). 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 Reference Thoroddsen2002; Josserand & Zaleski Reference Josserand and Zaleski2003) and air entrapment (Mani, Mandre & Brenner Reference Mani, Mandre and Brenner2010; Philippi et al. Reference Philippi, Lagrée and Antkowiak2016) 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 (Reference Josserand, Ray and Zaleski2016), 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. Reference Jian, Josserand, Popinet, Ray and Zaleski2018; Li et al. Reference Li, Thoraval, Marston and Thoroddsen2018), 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).
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20231213144154107-0623:S0022112023007012:S0022112023007012_fig8.png?pub-status=live)
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. (Reference Thoraval, Takehara, Etoh, Popinet, Ray, Josserand, Zaleski and Thoroddsen2012).
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. (Reference Li, Thoraval, Marston and Thoroddsen2018) (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 (Reference Josserand and Zaleski2003) 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. Reference Thoraval, Takehara, Etoh, Popinet, Ray, Josserand, Zaleski and Thoroddsen2012) 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).
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20231213144154107-0623:S0022112023007012:S0022112023007012_fig9.png?pub-status=live)
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. (Reference Thoraval, Takehara, Etoh, Popinet, Ray, Josserand, Zaleski and Thoroddsen2012). 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. (Reference Thoraval, Takehara, Etoh, Popinet, Ray, Josserand, Zaleski and Thoroddsen2012). 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.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20231213144154107-0623:S0022112023007012:S0022112023007012_fig10.png?pub-status=live)
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 Reference Pumphrey and Elmore1990). 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. Reference Zou, Ji, Yuan, Ren, Ruan and Fu2012; Wang, Kuan & Tsai Reference Wang, Kuan and Tsai2013; Thoraval, Li & Thoroddsen Reference Thoraval, Li and Thoroddsen2016; Deka et al. Reference Deka, Ray, Biswas, Dalal, Tsai and Wang2017). 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 Reference Worthington1908; Engel Reference Engel1966; van de Sande et al. Reference van de Sande, Smith and van Oord1974; Bisighini et al. Reference Bisighini, Cossali, Tropea and Roisman2010; Murphy et al. Reference Murphy, Li, d'Albignac, Morra and Katz2015; Lherm et al. Reference Lherm, Deguen, Alboussière and Landeau2021), 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 Reference Prosperetti, Crum and Pumphrey1989; Prosperetti & Oguz Reference Prosperetti and Oguz1993) and the additional source of airborne droplets caused by liquid film rupture (Resch, Darrozes & Afeti Reference Resch, Darrozes and Afeti1986; Afeti & Resch Reference Afeti and Resch1990; Resch & Afeti Reference Resch and Afeti1991). 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. Reference Bisighini, Cossali, Tropea and Roisman2010; Murphy et al. Reference Murphy, Li, d'Albignac, Morra and Katz2015; Lherm et al. Reference Lherm, Deguen, Alboussière and Landeau2021). High-speed jet ejections out of a liquid interface are commonly observed in many other physical processes, such as bubble bursting (Boulton-Stone & Blake Reference Boulton-Stone and Blake1993; Thoroddsen et al. Reference Thoroddsen, Takehara, Etoh and Ohl2009; Berny et al. Reference Berny, Deike, Popinet and Séon2022), Faraday waves (Hogrefe et al. Reference Hogrefe, Peffley, Goodridge, Shi, Hentschel and Lathrop1998; Zeff et al. Reference Zeff, Kleber, Fineberg and Lathrop2000) and cavity collapse induced by the process of solid/liquid object impact onto fluid target (Worthington & Cole Reference Worthington and Cole1897, Reference Worthington and Cole1900; Gekle & Gordillo Reference Gekle and Gordillo2010; Ray et al. Reference Ray, Biswas and Sharma2015; Jamali et al. Reference Jamali, Rostamijavanani, Nouri and Navidbakhsh2020; Kim et al. Reference Kim, Lee, Bose, Kim and Lee2021). 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 Reference Gekle and Gordillo2010). Surface tension plays an essential role in the jet formation process.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20231213144154107-0623:S0022112023007012:S0022112023007012_fig11.png?pub-status=live)
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. Reference Ray, Biswas and Sharma2015). 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 (Reference Engel1966) 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.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20231213144154107-0623:S0022112023007012:S0022112023007012_fig12.png?pub-status=live)
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. (Reference Murphy, Li, d'Albignac, Morra and Katz2015) (their figure 3). The production of such tiny droplets from the receding films has been studied by Lhuissier & Villermaux (Reference Lhuissier and Villermaux2012) and Dasouqi et al. (Reference Dasouqi, Yeom and Murphy2021).
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 Reference Engel1966; Thoroddsen Reference Thoroddsen2002; Bisighini & Cossali Reference Bisighini and Cossali2011). The tracer field is then advected by following equation:
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20231213144154107-0623:S0022112023007012:S0022112023007012_eqn4.png?pub-status=live)
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. (Reference Lhuissier, Sun, Prosperetti and Lohse2013). 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.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20231213144154107-0623:S0022112023007012:S0022112023007012_fig13.png?pub-status=live)
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 Reference Bisighini and Cossali2011). 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.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20231213144154107-0623:S0022112023007012:S0022112023007012_fig14.png?pub-status=live)
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. (Reference Berberović, van Hinsberg, Jakirlić, Roisman and Tropea2009). The solid line shows the theoretical estimation of the penetration depth proposed by Bisighini et al. (Reference Bisighini, Cossali, Tropea and Roisman2010). 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 Reference Fedorchenko and Wang2004). 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. Reference Birkhoff, MacDougall, Pugh and Taylor1948; Yarin, Rubin & Roisman Reference Yarin, Rubin and Roisman1995) and has been previously applied for the analytical study of crater evolution (Fedorchenko & Wang Reference Fedorchenko and Wang2004; Berberović et al. Reference Berberović, van Hinsberg, Jakirlić, Roisman and Tropea2009). 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. (Reference Berberović, van Hinsberg, Jakirlić, Roisman and Tropea2009) 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. (Reference Bisighini, Cossali, Tropea and Roisman2010) 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
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20231213144154107-0623:S0022112023007012:S0022112023007012_eqn5.png?pub-status=live)
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20231213144154107-0623:S0022112023007012:S0022112023007012_eqn6.png?pub-status=live)
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. (Reference Bisighini, Cossali, Tropea and Roisman2010), 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 Reference Okawa, Shiraishi and Mori2006; Guildenbecher et al. Reference Guildenbecher, Engvall, Gao, Grasser, Reu and Chen2014; Li, Zhang & Liu Reference Li, Zhang and Liu2019; Wu, Wang & Zhao Reference Wu, Wang and Zhao2020) and oblique (Okawa et al. Reference Okawa, Shiraishi and Mori2008; Liu Reference Liu2018) 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. Reference Thoroddsen, Thoraval, Takehara and Etoh2011; Thoraval et al. Reference Thoraval, Takehara, Etoh, Popinet, Ray, Josserand, Zaleski and Thoroddsen2012; Zhang et al. Reference Zhang, Toole, Fezzaa and Deegan2012; Thoraval et al. Reference Thoraval, Takehara, Etoh and Thoroddsen2013), suggesting that the creation of these tiny droplets might be associated with more complicated ‘irregular’ breakup of liquid sheet. According to Deegan et al. (Reference Deegan, Brunet and Eggers2007), 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. Reference Deegan, Brunet and Eggers2007; Marcotte et al. Reference Marcotte, Michon, Séon and Josserand2019). 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.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20231213144154107-0623:S0022112023007012:S0022112023007012_fig15.png?pub-status=live)
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.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20231213144154107-0623:S0022112023007012:S0022112023007012_fig16.png?pub-status=live)
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. (Reference Pairetti, Damian, Nigro, Popinet and Zaleski2020) showed that the most frequent droplet diameter during primary atomization is always near
$2\varDelta$ for all three proposed mesh resolutions; Mostert, Popinet & Deike (Reference Mostert, Popinet and Deike2022) 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.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20231213144154107-0623:S0022112023007012:S0022112023007012_fig17.png?pub-status=live)
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. (Reference Murphy, Li, d'Albignac, Morra and Katz2015), 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 Reference Afeti and Resch1990; Villermaux & Bossa Reference Villermaux and Bossa2009, Reference Villermaux and Bossa2011; Villermaux, Pistre & Lhuissier Reference Villermaux, Pistre and Lhuissier2013).
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.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20231213144154107-0623:S0022112023007012:S0022112023007012_fig18.png?pub-status=live)
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. (Reference Murphy, Li, d'Albignac, Morra and Katz2015), 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. (Reference Murphy, Li, d'Albignac, Morra and Katz2015), 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 Reference Oguz and Prosperetti1989) and the entrapment of isolated bubbles/bubble rings due to unstable oscillations of ejecta base (Weiss & Yarin Reference Weiss and Yarin1999). 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. Reference Thoraval, Takehara, Etoh and Thoroddsen2013; Li et al. Reference Li, Thoraval, Marston and Thoroddsen2018) 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.).
Supplementary movie
Supplementary movie is available at https://doi.org/10.1017/jfm.2023.701.
Acknowledgements
We appreciate the beneficial discussions and help from the BASILISK community. Simulations were performed using the computational resources of Advanced Research Computing at Virginia Tech. We thank three anonymous reviewers for their valuable comments which significantly improved the manuscript.
Funding
This work is supported by a scholarship from the China Scholarship Council (CSC) under the grant CSC no. 201908320462.
Declaration of interests
The authors report no conflict of interest.
Code availability
The codes to reproduce the results presented in this article are available in BASILISK sandbox (http://basilisk.fr/sandbox/HuiW/).
Appendix A. Effect of spatial resolution
A.1. Early-time splashes, droplets and bubbles
In this section, we discuss the effect of minimum spatial resolution on the early-time splashing dynamics of high-speed drop impact. Preliminary tests have been carried out based on the AMR algorithm implemented in the BASILISK code. For cases at $L_{max}=12$ (
$\varDelta =15.6\ \mathrm {\mu }{\rm m}$) and 13 (
$\varDelta =7.8\ \mathrm {\mu }{\rm m}$), the simulations are performed using the same maximum refinement levels throughout the entire computational domain. However, increasing the maximum level beyond this point is currently prohibitively expensive. As a workaround, we employ the mesh refinement strategy shown in figure 2(c) for increased resolutions. For higher maximum levels, we only apply
$L_{max}=14$ (
$\varDelta =3.9\ \mathrm {\mu }{\rm m}$) and 15 (
$\varDelta =1.95\ \mathrm {\mu }{\rm m}$) near the neck region (green area in figure 2c), while the rest of the domain use
$L_{max}=13$. This allows calculating the
$L_{max}=14$, 15 cases to
$t=200\ \mathrm {\mu }{\rm s}$ and
$t=150\ \mathrm {\mu }{\rm s}$ respectively, and gives access to the development of early-time splashing phenomena at finer resolutions. After the initial stage, the extra refinement layer is removed and
$L_{max}=13$ is used for the entire domain to allow the crown to develop.
Figure 19 shows different early-time splashing phenomena captured under varied spatial resolutions. At $L_{max}\leq 12$, a thin ejecta sheet rises ‘smoothly’ from the contact line and detaches rings of secondary droplets nearly axisymmetrically from its rim (figure 19a). At
$L_{max}=13$, the growth of the ejecta sheet remains ‘smooth’ at its early evolution, but it becomes thinner and rises faster with slightly higher spatial resolution. As the impact advances, the uprising ejecta eventually encounters and reconnects with the surface of the downward-moving drop (
$t=170\ \mathrm {\mu }{\rm s}$ in figure 19b), establishing a new contact line that moves downwards along the entire circumference. As a consequence, the liquid sheet above the new connection neck is completely torn into small droplets and an inward-oriented crown is formed hereafter as shown in the last frame of figure 19(b).
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20231213144154107-0623:S0022112023007012:S0022112023007012_fig19.png?pub-status=live)
Figure 19. Effect of maximum mesh refinement level on the early-time splashing behaviours of high-speed drop impact captured at (a) $L_{max}=12$, (b)
$L_{max}=13$, (c)
$L_{max}=14$ and (d)
$L_{max}=15$. From left to right, the corresponding times are 20, 90, 170 and
$360\ \mathrm {\mu }{\rm s}$ after contact. (e) Comparison of air–water interfaces at
$t=50\ \mathrm {\mu }{\rm s}$ captured by calculation at
$L_{max}=12$ (left), high-speed camera (middle) and calculation at
$L_{max}=14$ (right).
For cases at $L_{max}=14$, 15, very similar splashing behaviours are obtained as shown in figures 19(c) and 19(d). Immediately after impact (
$t=20\ \mathrm {\mu }{\rm s}$), a great number of ‘randomly’ distributed microdroplets are sprayed underneath the drop and an ejecta sheet is not observed. Subsequently, a highly disturbed liquid sheet emerges and interacts strongly with drop and pool, entrapping rings of microbubbles in the neck (
$t=90\ \mathrm {\mu }{\rm s}$). An outward-expanding crown is eventually formed above the free surface (
$t=360\ \mathrm {\mu }{\rm s}$). Referring the splashing classifications in Thoroddsen et al. (Reference Thoroddsen, Thoraval, Takehara and Etoh2011) and Thoraval et al. (Reference Thoraval, Takehara, Etoh, Popinet, Ray, Josserand, Zaleski and Thoroddsen2012) based on the dimensionless Ohnesorge (
$Oh=\mu _l/\rho _l\sigma d$) and splash (
$K=We\sqrt {Re}$) numbers, an irregular splash is always observed for their most energetic cases. Furthermore, Thoraval et al. (Reference Thoraval, Takehara, Etoh and Thoroddsen2013) have also found that the entrapment of bubble rings/arcs in the neck region usually occurs at
$Re>12\,000$ for water drops impacting on water pools. As the parameters in our present case are much higher than their range of study, thinner liquid sheet and more complex interfacial deformation/breakup can be even expected.
Figure 19(e) compares the shape of the air–water interface captured at $L_{max}=12$ (left), high-speed camera (middle) and at
$L_{max}=14$ (right)
$50\ \mathrm {\mu }{\rm s}$ after impact. The emerge-ruptured liquid sheet together with a great number of ‘irregularly’ distributed tiny droplets calculated at
$L_{max}=14$ are consistent with the experimental observations, while a relatively smooth ejecta is captured when spatial resolution is insufficient (
$L_{max}=12$ and 13).
Figure 20 shows size distributions of droplets and bubbles at $t=20$ and
$90\ \mathrm {\mu }{\rm s}$ produced by early-time ‘prompt splash’, corresponding to the first and second columns in figures 19(a)–19(d). As can be seen,
$L_{max}=12$, 13 cases produce very few droplets and bubbles at the very early time of impact and their data are not comparable to higher-resolution cases. For cases at
$L_{max}=14$, 15, very similar profiles are obtained for both droplet and bubble distributions, corroborating its intricate early-time flow dynamics near the neck region. It is noteworthy that, the most frequent diameter is always found in the
$2\varDelta \sim 4\varDelta$ range for both droplets and bubbles, and the data are approximately grid converged for diameters greater than
${\sim }4\varDelta$.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20231213144154107-0623:S0022112023007012:S0022112023007012_fig20.png?pub-status=live)
Figure 20. Effect of maximum mesh refinement level on the early-time size distributions of droplets and bubbles. The droplet distributions are shown at (a) $t=20\ \mathrm {\mu }{\rm s}$ (b)
$t=90\ \mathrm {\mu }{\rm s}$. The bubble distributions are shown at (c)
$t=20\ \mathrm {\mu }{\rm s}$, (d)
$t=90\ \mathrm {\mu }{\rm s}$. The vertical dotted lines show
$S_d,S_b=2\varDelta$ at different maximum refinement levels.
A.2. Energetics
In this section, we discuss the energetics in the process of high-speed drop impact. The kinetic energy $E_k$, gravitational potential energy
$E_g$ and surface potential energy
$E_s$ of the liquid phase are calculated as follows:
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20231213144154107-0623:S0022112023007012:S0022112023007012_eqn7.png?pub-status=live)
The integrals are computed over the entire volume and surface of the liquid phase in the domain. The time point when the first contact occurs between drop and pool is defined as $t=0$ ms. Here,
$E_{g0}$ and
$E_{s0}$ are the gravitational and surface potential energies of the liquid phase at
$t=0$ ms. Figure 21(a) shows the time evolution of
$E_k$,
$E_g$,
$E_s$ and their sum
$E_L=E_k+E_g+E_s$ during 1 ms after impact, normalized by the initial kinetic energy of the impact drop
$E_0$. For each term, the calculated results are very close to each other between the two higher maximum refinement levels at
$L_{max}=14$ and 15, whereas the data at
$L_{max}=13$ show a large discrepancy with them.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20231213144154107-0623:S0022112023007012:S0022112023007012_fig21.png?pub-status=live)
Figure 21. (a) Time evolution of energy aspects calculated at different resolutions. From bottom to top on the right for each resolution: gravitational potential energy $E_g$, surface potential energy
$E_s$, kinetic energy
$E_k$ and total mechanical energy
$E_L=E_k+E_g+E_s$ in the liquid phase. (b) Energy budget for the case calculated at
$L_{max}=15$ (figure 19d). From bottom to top on the right:
$E_g$,
$E_s$,
$E_k$,
$E_L$, the total mechanical energy including the kinetic and gravitational potential energies in the gas
$E_M=E_L+E_A$, and finally the total energy including the viscous energy dissipation in both gas and liquid phases
$E_T=E_M+E_d$. The vertical dotted line shows the time point at
$t=0.15$ ms, where the extra refinement layer is removed. The horizontal dotted line shows
$E_T/E_0=0.95$.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20231213144154107-0623:S0022112023007012:S0022112023007012_fig22.png?pub-status=live)
Figure 22. Effect of mesh resolution on droplet statistics in the process of crown fragmentation. Calculations are performed using the same input ‘restart’ file saved from the simulation of figure 19(c) at $t=0.43$ ms, and all the previously generated tiny droplets are removed at the first time step of continuations. (a,b) Show the air–water interfaces of the initial ‘restart’ file before and after removing tiny droplets. (c–e) Show temporal contours of droplet size distributions calculated at
$L_{max}=12$, 13 and 14 respectively. Here, the lower dashed line shows
$S_d=2\varDelta$ and the upper dashed line shows
$S_d=4\varDelta$. (f) Time-averaged droplet size distribution over the time window
$t\in [0.43, 1.17]$ ms. The vertical dotted lines show
$S_d=2\varDelta$ at different resolutions.
Figure 21(b) plots the temporal evolution of the energy budget for the case at maximum level $L_{max}=15$. The kinetic and gravitational potential energies of the gas phase are integrated over the gas volume using the same equations like (A1a–c), and the sum of the total mechanical energy in the gas is represented as
$E_A$. The total mechanical energy in the system is donated as
$E_M=E_L+E_A$. In addition, the energy dissipation due to viscosity in the gas–liquid system can be calculated directly from the deformation tensor as
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20231213144154107-0623:S0022112023007012:S0022112023007012_eqn8.png?pub-status=live)
Therefore, the total energy is given by $E_T=E_M+E_d$. In figure 21(b), it is evident that the total energy budget
$E_T$ is well conserved in the system after the formation of the thin-walled crown at
$t>0.2$ ms. However, during the early ‘prompt splash’ stage (
$t<0.2$ ms), there is a discrepancy where approximately 5 % of the total budget appears to be missing. The appearance of such an error, however, is reasonable regarding the rapid and violent nature of the phenomena involved in such energetic impact conditions. This suggests that resolution at
$L_{max}=15$ remains not sufficient to fully capture the chaotic small-scale features near the neck region, leading to energy losses due to unresolved turbulent dissipation. As the crown arises, the flow settles and the gradients become smoother, the dissipation decreases, leading to improved energy conservation at the later time of impact (
$t>0.2$ ms).
A.3. Sensitivity analysis of droplet statistics
In order to further quantify the sensitivity of droplet statistics to changes in mesh resolution, additional tests have been conducted by re-starting the simulations under different maximum refinement levels based on the same initial input ‘restart’ file. The chosen initial ‘restart’ is the saved simulation result from the case in figure 19(c) ($L_{max}=14$) at
$t=0.43$ ms, where new droplets are primarily produced by the mechanism of crown sheet fragmentation as shown in figure 22(a). In addition, we label and remove all the small droplets in the computational domain at the first time step of the continuations as shown in figure 22(b), to focus solely on the difference of droplet statistics caused by mesh resolutions.
Figures 22(c)–22(e) show contours of droplet size distribution with time for three test cases calculated at $L_{max}=13$, 14 and 15. The following features can be clearly observed by comparing these contours: (i) all cases predict two separated ribbons of droplet size distribution over the entire simulations; (ii) the primary size peaks (large size) are in the
$80\sim 250\ \mathrm {\mu }{\rm m}$ range. The distributions of large droplets share very similar qualitative features, but are not pointwise converged; (iii) secondary size peaks (small size) are captured at all refinement levels for droplet sizes in the range
$2\varDelta \sim 4\varDelta$, indicating that the production of small droplets is greatly affected by minimum cell size. Note that similar effects in small-sized droplets are also found at the early-time splash stage as shown in figure 20.
Figure 22(f) shows time-averaged droplet size distribution over the time window $t\in [0.43, 1.17]$ ms for these three test cases. In this comparison, the data are well converged for large droplets at the primary peak. A great number of small droplets, representing more than
$25\,\%\sim 30\,\%$ of the total droplet count, are present in the range
$2\varDelta \sim 4\varDelta$. As before, this secondary peak shifts towards smaller sizes when
$L_{max}$ is increased.
The aforementioned observations (figures 20 and 22) demonstrate that mesh convergence is only achieved when the droplet size is significantly larger than the cell size and the presence of secondary peaks in numerical results is attributed to the limitation of minimum cell size, suggesting that at least 4 cells per droplet diameter are required for numerical convergence.