1. Introduction
The ubiquitous phenomenon of drop impact plays a vital role in numerous natural and technological processes, influencing a diverse array of occurrences from surface erosion and soil stability to the performance of advanced manufacturing techniques. In the natural realm, raindrops influence the response of plant surfaces and the transmission of seeds and pathogens (Kim et al. Reference Kim, Park, Gruszewski, Schmale and Jung2019; Bourouiba Reference Bourouiba2021; Roth-Nebelsick et al. Reference Roth-Nebelsick, Konrad, Ebner, Miranda, Thielen and Nebelsick2022), the spreading of petroleum contaminants (Murphy et al. Reference Murphy, Li, d’Albignac, Morra and Katz2015; Castillo-Orozco et al. Reference Castillo-Orozco, Davanlou, Choudhury and Kumar2016) and soil and structure erosion (Assouline Reference Assouline2004; Zhou et al. Reference Zhou, Li, Chen, Xu, Hui and Zhang2009; Fernández-Raga et al. Reference Fernández-Raga, Palencia, Keesstra, Jordán, Fraile, Angulo-Martínez and Cerdà2017; Lardier et al. Reference Lardier, Roudier, Clothier and Willmott2019). In the realm of industrial applications, it plays a crucial role in processes such as spray injecting and cooling (Breitenbach, Roisman & Tropea Reference Breitenbach, Roisman and Tropea2018; Wang et al. Reference Wang, Guo, Xiong and Wang2020; Benther et al. Reference Benther, Pelaez-Restrepo, Stanley and Rosengarten2021), inkjet printing (Lohse Reference Lohse2022), agriculture irrigation (Gart et al. Reference Gart, Mates, Megaridis and Jung2015; Gilet & Bourouiba Reference Gilet and Bourouiba2015) and internal combustion engines (Panão & Moreira Reference Panão and Moreira2005; Moreira, Moita & Panao Reference Moreira, Moita and Panao2010). Scientific research aimed at uncovering the internal mechanisms and fundamental physics of drop impact thus holds immense potential in facilitating various technical applications where this process is involved, but remains challenging due to increased complexities of multiscale physics and interfacial flows (Deegan, Brunet & Eggers Reference Deegan, Brunet and Eggers2007).
Since the pioneering works of Worthington (Reference Worthington1883, Reference Worthington1895) more than a century ago, different scenarios of drop impact, including solid target (Josserand & Thoroddsen Reference Josserand and Thoroddsen2016; Jian et al. Reference Jian, Josserand, Popinet, Ray and Zaleski2018), thin and shallow liquid film (Mukherjee & Abraham Reference Mukherjee and Abraham2007; Coppola et al. Reference Coppola, Rocco and de Luca2011; Sykes et al. Reference Sykes, Cimpeanu, Fudge, Castrejón-Pita and Castrejón-Pita2023), deep liquid pool (Murphy et al. Reference Murphy, Li, d’Albignac, Morra and Katz2015; Ray, Biswas & Sharma Reference Ray, Biswas and Sharma2015; Michon, Josserand & Séon Reference Michon, Josserand and Séon2017), miscible or immiscible liquids (Lhuissier et al. Reference Lhuissier, Sun, Prosperetti and Lohse2013; Castillo-Orozco et al. Reference Castillo-Orozco, Davanlou, Choudhury and Kumar2016; Lherm et al. Reference Lherm, Deguen, Alboussière and Landeau2022) and normal or oblique angle (Leneweit et al. Reference Leneweit, Koehler, Roesner and Schäfer2005; Gielen et al. Reference Gielen, Sleutel, Benschop, Riepen, Voronina, Visser, Lohse, Snoeijer, Versluis and Gelderblom2017; Liu Reference Liu2018), have been categorised and documented in the literature (Rein Reference Rein1993; Yarin Reference Yarin2006; Liang & Mudawar Reference Liang and Mudawar2016). For normal (perpendicular) drop impact onto the deep volume of the same liquid, three dimensionless parameters are usually used to characterise the splashing dynamics, namely the Froude number
$\textit{Fr}=U_0/\sqrt {gd}$
, the Reynolds number
${\textit{Re}}=\rho _l U_0d/\mu _l$
and the Weber number
${\textit{We}}=\rho _lU_0^2 d/\sigma$
, where
$U_0$
is the impact speed,
$d$
the drop diameter,
$g$
the gravitational acceleration,
$\rho _l$
the liquid density,
$\mu _l$
the liquid dynamic viscosity, and
$\sigma$
the liquid surface tension.
In the regime of splashing, significant progress has been made in delineating its early-stage morphological behaviours over the past decades. As the drop touches the pool, a thin liquid ejecta is shot out horizontally from the connection of two masses due to localised pressure buildup (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; Josserand, Ray & Zaleski Reference Josserand, Ray and Zaleski2016) and expands rapidly outwards, resulting in a series of rich splashing processes afterwards. In some cases (
$2000\lt {\textit{Re}}\lt 6000$
), the base of the ejecta may oscillate, shedding a vortex that later lifts a ‘lamella’ (Zhang et al. Reference Zhang, Toole, Fezzaa and Deegan2012; Agbaglah et al. Reference Agbaglah, Thoraval, Thoroddsen, Zhang, Fezzaa and Deegan2015; Sykes et al. Reference Sykes, Cimpeanu, Fudge, Castrejón-Pita and Castrejón-Pita2023) or propelling the ejecta to bend and collide with the surfaces of drop or pool (Thoroddsen et al. Reference Thoroddsen, Thoraval, Takehara and Etoh2011; Agbaglah et al. Reference Agbaglah, Thoraval, Thoroddsen, Zhang, Fezzaa and Deegan2015; Tian et al. Reference Tian, Aljedaani, Alghamdi and Thoroddsen2024), potentially entrapping a large toroidal bubble in the substrate liquid (Thoroddsen Reference Thoroddsen2002; Deegan et al. Reference Deegan, Brunet and Eggers2007). For the highest range of
${\textit{Re}}$
, the most unstable splashing mode occurs and the ‘emerging-rupturing-ejecta’ is observed (Thoroddsen Reference Thoroddsen2002), where the ejecta interacts strongly with the free surfaces, entrapping a row of toroidal bubbles inside the liquid tank (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). Despite the improved understanding of the initial jet multiplicity owing to the enhanced capabilities in both experimental and numerical approaches, the intricate interplays between jets and free surfaces remain fascinating, and the critical conditions that determine the transitions between different splashing behaviours are still not well elucidated.
After the initial stage, a typical drop-impingement-induced structure, composed of a submerging cavity surrounded by a thin-walled liquid crown, is formed. It is well established in the literature that the overall morphology during this period depends greatly on the impact parameters. When the drop impinges at very low energy, a thick tongue (swell) is generated along the drop–pool connection and the drop coalesces with the receiving bulk (Chapman & Critchlow Reference Chapman and Critchlow1967; Rodriguez & Mesler Reference Rodriguez and Mesler1985, Reference Rodriguez and Mesler1988; Mohamed-Kassim & Longmire Reference Mohamed-Kassim and Longmire2004). For slightly higher impact energy, a flared liquid crown rises up cylindrically along the periphery of a semispherical cavity (Fedorchenko & Wang Reference Fedorchenko and Wang2004; Ray et al. Reference Ray, Biswas and Sharma2015). Further increasing the impact energy, the crown is more developed, and liquid ligaments are stretched out of its upper rim, producing and sending out secondary droplets into the air (Okawa, Shiraishi & Mori Reference Okawa, Shiraishi and Mori2006; Deegan et al. Reference Deegan, Brunet and Eggers2007; Guildenbecher et al. Reference Guildenbecher, Engvall, Gao, Grasser, Reu and Chen2014; Wu, Wang & Zhao Reference Wu, Wang and Zhao2020). For these abovementioned cases, after the crown (swell) reaches its maximum position, capillary ripples start to propagate on the interior wall of the crater and subsequently lead to the collapse of the structure, resulting in an axial upward jet rising out of the cavity bottom, known as the ‘Worthington jet’ (Ray et al. Reference Ray, Biswas and Sharma2015; Michon et al. Reference Michon, Josserand and Séon2017).
In the most energetic regime of drop impact, the crown undergoes further elevation and its upper edge necks in before retraction, encapsulating a large gas volume within it (Worthington Reference Worthington1883), which is qualitatively different from the phenomena already mentioned. Upon surface closure, an upward jet and a downward jet protrude from the closure point due to a large axial pressure gradient. The formation of similar vertical jets is also a subject of intense study in other physical processes, such as rapid boiling (Marston & Li Reference Marston and Li2019) and cavitation (Saade et al. Reference Saade, Jalaal, Prosperetti and Lohse2021). The downward jet then moves towards the cavity bottom and eventually mixes with the upward-rising jet coming from the cavity floor (Engel Reference Engel1966). During the final stages, a thin-filmed bubble forms over the pool surface, ultimately breaking up due to combined effects of surface tension and external disturbances (Motzkus, Gensdarmes & Géhin Reference Motzkus, Gensdarmes and Géhin2009; Murphy et al. Reference Murphy, Li, d’Albignac, Morra and Katz2015). The importance of understanding such physical processes under high-energy impact conditions was evidenced by Murphy et al. (Reference Murphy, Li, d’Albignac, Morra and Katz2015) in the context of post-oil-spill accidents on the ocean surface, where additional spraying sources of marine aerosols (Dasouqi, Yeom & Murphy Reference Dasouqi, Yeom and Murphy2021) could be introduced by the impact of raindrops or use of oil dispersant, which could therefore cause adverse effects on the environment, economy and public health.
For drop impact on thin liquid layers, previous studies have provided some quantitative investigations of the critical conditions for large-bubble encapsulation at high-impact velocities (Pan et al. Reference Pan, Cheng, Chou and Wang2008; Pan & Hung Reference Pan and Hung2010; Geppert et al. Reference Geppert, Chatzianagnostou, Meister, Gomaa, Lamanna and Weigand2016; Aljedaani et al. Reference Aljedaani, Wang, Jetly and Thoroddsen2018; Ribeiro et al. Reference Ribeiro, Silva and Panão2023a , Reference Ribeiro, Panão, Barata and Silvab ), focusing primarily on the effects of liquid-film thickness and underlying material on the splash. In such scenarios, the growth of the submerging cavity is greatly inhibited by the finite liquid film, and the interactions between the enlarging cavity and the underlying wall are crucial factors influencing the closure event (Pan et al. Reference Pan, Cheng, Chou and Wang2008). However, in the case of deep receiving pools, the cavity undergoes rapid expansion predominantly in the vertical direction, independent of the presence of walls. As a result, the cavity takes on a hemispherical shape (Engel Reference Engel1966, Reference Engel1967; 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). This distinct crater morphology therefore gives rise to different characteristics and dynamics in the process of bubble-canopy (BC) formation compared with scenarios impacting thin liquid films.
Notably, while certain studies have mentioned the phenomenon of BC formation in deep receiving pools within the context of underwater noise (Medwin et al. Reference Medwin, Nystuen, Jacobus, Ostwald and Snyder1992; Prosperetti & Oguz Reference Prosperetti and Oguz1993), crater expansion (Worthington Reference Worthington1883; Engel Reference Engel1966, Reference Engel1967; Bisighini et al. Reference Bisighini, Cossali, Tropea and Roisman2010), marine aerosols (Murphy et al. Reference Murphy, Li, d’Albignac, Morra and Katz2015) and Rayleigh–Taylor instability (Lherm et al. Reference Lherm, Deguen, Alboussière and Landeau2021), limited attention has been devoted to investigating the underlying determinants of this specific phenomenon. The influence of liquid viscosity and surface tension on the characteristic behaviours of the closure event has been initiated (Murphy et al. Reference Murphy, Li, d’Albignac, Morra and Katz2015; Sochan et al. Reference Sochan, Beczek, Mazur, Ryżak and Bieganowski2018), yet the critical transitional criteria are unsettled and the mechanisms of surface closure remain to be understood. Therefore, in this paper, we study the dynamics of drop impact on the deep pool of the same liquid using direct numerical simulations. Various splashing behaviours, spanning from gentle coalescence to the most energetic splash of BC, are investigated by systematically varying the dimensionless Reynolds and Weber numbers. A comprehensive regime map of drop splash that incorporates both early-time jet dynamics and BC formation, covering a broad parameter space of
$1000\leqslant {\textit{Re}} \leqslant 30000$
and
$100\leqslant We \leqslant 3000$
, is therefore concluded.
The paper is organised as follows. In § 2, we describe the governing equations, numerical approaches and general numerical set-ups. In § 3, detailed comparisons between numerical and experimental results are conducted. In § 4, we present numerical investigations of splashing dynamics under varied impact conditions. The paper ends with the conclusions in § 5.
2. Numerical methods
2.1. Governing equations
To simulate the drop-impact gas–liquid system, we solve the two-phase incompressible Navier–Stokes equations with respective densities (
$\rho _g$
,
$\rho _l$
) and viscosities (
$\mu _g$
,
$\mu _l$
) (subscript
$g$
and
$l$
denote the properties in the gas and in the liquid, respectively):
\begin{align}&\qquad\qquad\qquad\qquad\qquad \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{U}=0, \end{align}
\begin{align}& \rho \left (\frac {\partial \boldsymbol{U}}{\partial t}+(\boldsymbol{U}\boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{U}\right )=-\boldsymbol{\nabla }P+\boldsymbol{\nabla }\boldsymbol{\cdot }(2\mu \unicode{x1D63F}\kern1pt)+\rho \boldsymbol{a}+\sigma \kappa \delta _s\boldsymbol{n}, \\[9pt] \nonumber \end{align}
where
$\boldsymbol{U}$
is the velocity vector,
$t$
is time,
$P$
is the pressure,
$\unicode{x1D63F}=[\boldsymbol{\nabla }\boldsymbol{U}+(\boldsymbol {\nabla }\boldsymbol{U})^T]/2$
is the deformation tensor,
$\boldsymbol{a}$
is the body force along the impact direction,
$\delta _s$
is the Dirac function indicating the interface, and
$\sigma$
,
$\kappa$
,
$\boldsymbol{n}$
are the surface tension, curvature and normal vector of the interface, respectively. The density
$\rho$
and viscosity
$\mu$
take constant values in each phase.
We choose the reference values of drop diameter
$d$
, initial impact speed
$U_0$
, liquid density
$\rho _l$
and gravitational acceleration
$g$
to non-dimensionalise the Navier–Stokes equations. Therefore, the dimensionless variables scaled by the characteristic quantities are
\begin{equation} \hat {\boldsymbol{X}}=\frac {\boldsymbol{X}}{d},\quad \hat {\boldsymbol{U}}=\frac {\boldsymbol{U}}{U_0},\quad \hat {t}=\frac {tU_0}{d},\quad \hat {P}=\frac {P}{\rho _l U_0^2},\quad \hat {\boldsymbol{a}}=\frac {\boldsymbol{a}}{g},\quad \hat {\kappa }=\kappa d,\quad \hat {\delta _s}=\delta _sd. \end{equation}
where
$\boldsymbol{X}$
is the three-dimensional coordinates. We substitute these dimensionless variables into the governing equations. Then (2.2) can be rewritten as
\begin{equation} \frac {\partial \hat {\boldsymbol{U}}}{\partial \hat {t}}+(\hat {\boldsymbol{U}}\boldsymbol{\cdot }\boldsymbol{\nabla})\hat{\boldsymbol{U}}=-\boldsymbol{\nabla }\hat {P}+\frac {1}{\textit{Re}}\boldsymbol{\nabla }\boldsymbol{\cdot }(2\unicode{x1D63F}\kern1pt)+\frac {1}{\textit{Fr}^2}\hat {\boldsymbol{a}}+\frac {1}{\textit{We}}\hat {\kappa }\hat {\delta _s}\boldsymbol{n}. \end{equation}
This final form incorporates three dimensionless numbers:
\begin{equation} {\textit{Re}}=\frac {\rho _l dU_0}{\mu _l},\quad We=\frac {\rho _l dU_0^2}{\sigma },\quad \textit{Fr}=\frac {U_0}{\sqrt {gd}}, \qquad \end{equation}
where the Reynolds number (
${\textit{Re}}$
), Weber number (
${\textit{We}}$
) and Froude number (
$\textit{Fr}$
) represent ratios of inertial forces to viscous, surface tension and gravitational forces, respectively. These dimensionless numbers are crucial for understanding the dominant physical mechanisms governing fluid flows and classifying various splashing behaviours.
2.2. Numerical solver
We use the open-source Basilisk C (Popinet & collaborators 2013–2025) to solve the two-phase flow system. The momentum conserving volume-of-fluid (VOF) numerical scheme is employed here (Fuster & Popinet Reference Fuster and Popinet2018). The gas–liquid interface is tracked and reconstructed using the VOF method (Scardovelli & Zaleski Reference Scardovelli and Zaleski1999), and the volume fraction
$f(\boldsymbol{X},t)$
of the first fluid (liquid) is advected and described as follows:
\begin{equation} \frac {\partial f}{\partial t}+\boldsymbol{\nabla }\boldsymbol{\cdot }(f\boldsymbol{U})=0,\quad f(\boldsymbol{X},t) = \left \{\!\! \begin{array}{rcl} f=1 & \mbox{for} & {\rm liquid} \\ 0\lt f\lt 1 & \mbox{for} & {\rm interface} \\ f=0 & \mbox{for} & {\rm gas}. \end{array}\right . \end{equation}
The local density and viscosity in the computational domain then can be computed by the arithmetic means using
$\rho (f)=f\rho _l+(1-f)\rho _g$
and
$\mu (f)=f\mu _l+(1-f)\mu _g$
. The surface tension is calculated using the continuum surface force (CSF) model described in Francois et al. (Reference Francois, Cummins, Dendy, Kothe, Sicilian and Williams2006). The tree-based structure of spatial discretisation and adaptive mesh refinement (AMR) techniques (Van Hooft et al. Reference van Hooft, Popinet, van Heerwaarden, van der Linden, de Roode and van de Wiel2018) are applied to adaptively follow the smallest structures of the flow, thus concentrating the computational resources on the most concerning regions. Starting from a ‘parent’ cell at level 0, four ‘children’ cells are generated at level 1 in two dimensions. Various levels of local refinement can be present within the spatial domain at the same time. With the maximum level of refinement
$L_{\textit{max}}=n$
, the finest resolution in each dimension is reduced to
$1/2^n$
of the initial computational domain. The grid is refined or coarsened adaptively based on the wavelet-estimated spatially discretised field errors using AMR. The Basilisk solver has been validated extensively on the problems with multiphase complex flows as can be found in Popinet & collaborators (2013–2025).
All the numerical results presented in the main text of this paper were performed on 128 cores using the computational resources at the Advanced Research Computing (ARC) at Virginia Tech.
2.3. Axisymmetric flow configurations
In the literature, the shape of the impact drop before contact was widely assumed as a sphere (Berberović et al. Reference Berberović, van, Nils, Jakirlić, Roisman and Tropea2009; Thoraval et al. Reference Thoraval, Takehara, Etoh, Popinet, Ray, Josserand, Zaleski and Thoroddsen2012; Agbaglah et al. Reference Agbaglah, Thoraval, Thoroddsen, Zhang, Fezzaa and Deegan2015; Ray et al. Reference Ray, Biswas and Sharma2015; Fudge et al. Reference Fudge, Cimpeanu and Castrejón-Pita2021; Wu et al. Reference Wu, Zhang, Xiao and Ni2021), but this may not be the fact in many situations (Soni et al. Reference Soni, Kirar, Kolhe and Sahu2020). In real experiments, right after leaving the nozzle of the injector, the drop will fall freely through the surrounding gas to accelerate to a certain impact velocity. During this displacement, significant deformation might be observed depending on various factors (liquid properties, drop size, temperature and surrounding gas) (Volkov, Kuznetsov & Strizhak Reference Volkov, Kuznetsov and Strizhak2015), which would accordingly alter the physical dynamics of the impact (Scolan & Korobkin Reference Scolan and Korobkin2001; Thoraval, Li & Thoroddsen Reference Thoraval, Li and Thoroddsen2016). Particularly, the shape versus height oscillation was tested and modelled by Thoraval et al. (Reference Thoraval, Takehara, Etoh and Thoroddsen2013), who highlighted the relevance of pinch-off oscillations to the dynamics of the impact. In our current study, we define the drop as an oblate shape based on the experimental measurements of Murphy et al. (Reference Murphy, Li, d’Albignac, Morra and Katz2015). The initial horizontal and vertical diameters of the drop are given as
$d_h=4.3$
mm and
$d_v=3.8$
mm, thus resulting in an effective drop diameter calculated as
$d=(d_vd_h^2)^{1/3}=4.1$
mm.
Axisymmetric configurations for the simulation of drop impact on the same liquid pool. The blue represents the heavier fluid (e.g. water) and the white represents the lighter one (e.g. air). The dimensions shown here are not to scale.
In this study, simulations are performed in axisymmetric configurations as shown in figure 1. The computational domain is a square (cylinder in axisymmetric coordinates) and the side length of the domain is
$L=20d$
. The free surface of the pool is located in the middle of the axial direction and the depth of the pool is
$H=10d$
, providing sufficient space for crater development while minimising boundary influences (see discussion in Appendix B). The initial gap between drop and pool is
$\delta =0.1d$
, which allows taking the effects of air sheet entrapment into account. The free outflow boundary conditions are imposed on the top of the domain, while the default symmetry boundary conditions are applied for the other three boundaries. At the initial time, the mesh around the drop and the free surface of the pool is refined at the maximum refinement level, and the mesh is gradually coarsened away from the interfaces, which promises a good geometric description of the drop and free surface at the beginning. During the simulations, grids are adapted based on the wavelet-estimated errors of the volume fraction (tolerance
$f_{\textit{err}}=1\textrm e-4$
) and velocity (tolerance
$u_{\textit{err}}=1\textrm e-2$
) fields using AMR. The effect of gravitational force is taken into consideration in this study.
As discussed in Appendix A, our preliminary mesh study has found that a minimum grid resolution of at least 1638 cells per drop diameter, corresponding to a maximum refinement level at
$L_{\textit{max}}$
= 15, is necessary for capturing the ‘irregular splash’ regime induced by the most energetic drop impact at the highest
${\textit{Re}}$
and
${\textit{We}}$
range. The ‘primary splashing features’ are approximately grid-converged at
$L_{\textit{max}}=15$
, 16 and 17. Hence, in the main text of this manuscript, all the numerical results are obtained using the maximum refinement level
$L_{\textit{max}}=15$
, for the purpose of parametric investigation through a broad range of parameter space. The minimum level of refinement is fixed to
$L_{\textit{min}}$
= 7 during calculations.
It should be noted that this work is aimed to capture the primary features of high-speed drop impact that are considered symmetric. For those who are interested in analysing the non-symmetric phenomena such as neck instabilities, crown fragmentation, droplet statistics and motion of central jet, full three-dimensional simulations are strongly suggested (Wang et al. Reference Wang, Liu, Bayeul-Lainé, Murphy, Katz and Coutier-Delgosha2023). All the data hereafter are scaled by the initial drop diameter
$d$
and the impact speed
$U_0$
, thus the time is scaled as
$\hat {t}=tU_0/d$
.
3. Comparisons with experiment
This section numerically reproduces one particular case of a water drop impacting a deep water pool. Detailed comparisons between numerical results and available experimental data from Murphy et al. (Reference Murphy, Li, d’Albignac, Morra and Katz2015) are performed to validate the employed numerical schemes in the present study. In this calculation, artificial seawater is used for both the drop and pool, since the properties of the liquids in the drop and pool are almost the same, and air for the surrounding gas. The properties of seawater are density
$\rho _l=1018.3$
kg m
$^{-3}$
, dynamic viscosity
$\mu _l=1.0\times 10^{-3}$
N s m
$^{-2}$
and surface tension
$\sigma =73.0\times 10^{-3}$
N m
$^{-1}$
. The drop is initialised with
$d_h=4.3$
mm,
$d_v=3.8$
mm and
$U_0=7.2$
m s
$^{-1}$
, as same as the experimental measurements in Murphy et al. (Reference Murphy, Li, d’Albignac, Morra and Katz2015). This leads to a drop impact set-up at
${\textit{Re}}=30\,060$
and
${\textit{We}}=2964$
, corresponding to the highest parameters investigated in this study.
3.1. Morphology and kinematics
Typical sequence of splashing events for a water drop impacting a deep water pool. The drop for experiment and simulation has the same parameters:
$d_h=4.3$
mm,
$d_v=3.8$
mm and
$U_0=7.2$
m s−1 (
${\textit{Re}}=30\,060$
,
${\textit{We}}=2964$
). The experiment and simulation are presented with the same scale. The scale bar is
$2.0d$
long, where
$d=4.1$
mm. The images in (a) are extracted from one of the replicated counterparts of the control case in Murphy et al. (Reference Murphy, Li, d’Albignac, Morra and Katz2015), and the snapshots in (b) are produced by our axisymmetric simulation. The green asterisks indicate the tracked positions of the upper rim of the crown in the experiment. The red and blue colours in the simulation represent the fluids from the drop and the pool respectively. See supplementary movie 1 available at https://doi.org/10.1017/jfm.2025.10543.
Figure 2 shows a typical sequence of the event produced by high-speed raindrops impacting the ocean surface. It can be seen that the simulation successfully reproduced most of the distinctive physical features observed by the high-speed camera, including the formation of an initial outward-expanding crown around a disk-like cavity (
$\hat {t}=1.8$
), the transition of the orientation of the crown rim from horizontal-outward to vertical-upward along with the growth of a hemispherical cavity (
$\hat {t}=4.5$
), the inward bending of the upper portion of the crown (
$\hat {t}=9.0$
), the formation of BC and the protrusion of vertical jets from the closure point due to localised pressure buildup (
$\hat {t}=12.6$
), and finally the impingement of the downward central jet to the cavity bottom (
$\hat {t}=24.75$
). An overall good qualitative comparison is obtained. It is also worth noting that an elliptic pattern can be observed in the experimental images at certain instances (e.g.
$\hat {t}=5.4, 12.6, 21.6$
). This feature is caused by lensing effects, where the curvature of the crown acts as a lens, refracting light and producing the observed pattern. As the crown evolves dynamically, changes in the local surface curvature can alter the position of the lensing effect, leading to its shifting appearance in different frames.
Quantitative comparisons of axisymmetric numerical results with experimental measurements of Murphy et al. (Reference Murphy, Li, d’Albignac, Morra and Katz2015). (a) Sketch of the tracked quantities: crown neck position (
$\widetilde {R_r}$
,
$\widetilde {R_h}$
), cavity dimensions (
$\widetilde {C_r}$
,
$\widetilde {C_d}$
). The radial and axial coordinates are rescaled by
$d$
as
$\widetilde{R}$
and
$\widetilde{Z}$
. (b) Trajectory of the crown neck position. The first point is at
$\hat {t}=1.8$
. The time delay in experimental data between points is
$\hat {t}=1.8$
, and the time delay in numerical results is
$\hat {t}=0.45$
. (c) Time evolution of the crown neck radius (
$\widetilde {R_r}$
). (d) Time evolution of the crown neck height (
$\widetilde {R_h}$
). (e) Time evolution of the cavity width (
$\widetilde {C_r}$
). (f) Time evolution of the cavity depth (
$\widetilde {C_d}$
). The vertical dotted lines indicate the timing of the crown closure in axisymmetric simulation. The error bars indicate the standard deviation in experimental data.
To further verify our current numerical strategies, figure 3 presents a comparison between various quantities extracted from axisymmetric simulation and experiment. In the experiments of Murphy et al. (Reference Murphy, Li, d’Albignac, Morra and Katz2015), the crown rim position is defined as the outer edge of the crown neck where ligaments are formed, as indicated by green asterisks in figure 2. In our axisymmetric simulation, the trajectory of the crown neck position (
$\widetilde {R_r}$
,
$\widetilde {R_h}$
) is assessed using the diagrams depicted in figure 3(a). During the crown outward-expanding phase, the crown neck position corresponds to the highest point on the crown interface, as sketched at
$\hat {t}=3.6$
in figure 3(a). In the inward-bending surface closure phase, we track the point on the crown surface closest to the axis of symmetry, as sketched at
$\hat {t}=9.45$
in figure 3(a). The cavity depth (
$\widetilde {C_d}$
) and width (
$\widetilde {C_r}$
) are monitored by tracking the lowest position of the cavity and the cavity’s horizontal radius along the initial quiescent free surface of the pool, respectively.
As shown in figure 3(b), the numerical results for crown neck kinematics show a generally similar trend to the experimental measurements. However, the advancement of the impact in axisymmetric simulation evolves more rapidly compared with real experiments, leading to the accelerated growth of the upper crown and earlier occurrence of bubble closure, as clearly shown in the temporal evolution of
$\widetilde {R_r}$
and
$\widetilde {R_h}$
in figures 3(c) and 3(d). This discrepancy arises from the absence of strong three-dimensional phenomena along the thickened crown rim, as discussed in detail in § 3.2. The discrepancy in crown evolution could accordingly alter the development of the subsurface cavity. Figures 3(e) and 3(f) plot the temporal evolution of the cavity dimensions (
$\widetilde {C_r}$
,
$\widetilde {C_d}$
), revealing general good agreement between numerical and experimental results before closure. For highly energetic impacts like the present case, the cavity remains in its expanding phase after the closure of the upper crown (
$\hat {t}\gt 12.6$
in simulation), thereby drawing gas into the subsurface cavity by pulling the bubble ceiling downwards. As a result, the cavity expansion along the free surface of the receiving pool is significantly impeded as illustrated in figure 3(e). This suggests that an earlier surface closure event would hinder the development of cavity. Interestingly, in the experimental study of Mansoor et al. (Reference Mansoor, Marston, Vakarelski and Thoroddsen2014), further extended cavities were achieved and examined by artificially eliminating this surface closure process for the impacts of superhydrophobic spheres onto the deep liquid pool. The depth of the cavity, as shown in figure 3(f), except for minor disturbances caused by the impingement of the downward-moving jet, follows a similar trend in both the numerical and experimental results.
3.2. Limitations of axisymmetric simulation
In the aforementioned comparisons, a noticeable phase difference in the timing of crown growth and closure event is observed between our axisymmetric simulation and experimental observations. Specifically, the simulation exhibits faster growth and closure of the upper crown compared with the experiments. This discrepancy can be attributed to the inherent limitations of axisymmetric approaches. In real-world phenomena, during high-speed drop impact, crown development is determined by a three-dimensional rim that preserves strong transverse instabilities followed by corrugation formation and continuous droplet shedding (Roisman, Horvat & Tropea Reference Roisman, Horvat and Tropea2006; Deegan et al. Reference Deegan, Brunet and Eggers2007; Murphy et al. Reference Murphy, Li, d’Albignac, Morra and Katz2015). This complex dynamics gives rise to intricate flow patterns and additional interactions with the surrounding environment within the upper crown, which plays a crucial role in delaying the development of the structure (Wang et al. Reference Wang, Liu, Bayeul-Lainé, Murphy, Katz and Coutier-Delgosha2023). However, in axisymmetric simulations, the upper part of the crown evolves into a very thin and highly sensitive interface, serving as the determining factor for the crown dynamics and closure event. The absence of these three-dimensional effects in axisymmetric simulations eventually leads to an underestimation of the time required for the crown to develop and enclose, compared with the experimental observations.
Furthermore, the calculation of the downward-moving jet protruding from the closure point affects the subsequent dynamics. In reality, this jet arises from the rapid local pressure buildup due to the convergence of crown rims from all directions. Previous experiments have clearly demonstrated the non-axisymmetric footprint of this central jet and the off-centre point where it pierces the cavity (Bisighini et al. Reference Bisighini, Cossali, Tropea and Roisman2010; Murphy et al. Reference Murphy, Li, d’Albignac, Morra and Katz2015), indicating that the jet kinematics is greatly influenced by three-dimensional effects. However, this physical process cannot be accurately captured by an axisymmetric configuration, as the jet is only composed of flows along the symmetry plane. As a consequence, the jet in axisymmetric simulation is confined to attach to the symmetry axis and rapidly penetrates the bottom of the cavity floor (
$\hat {t}=24.75$
in figure 2), thereby further accelerating the impact process.
While acknowledging the limitations and simplifications made in the model, axisymmetric simulations remain an effective tool for studying drop-impact phenomena and it has provided valuable insights into the understanding of splashing dynamics, as evidenced by previous studies (Thoraval et al. Reference Thoraval, Takehara, Etoh, Popinet, Ray, Josserand, Zaleski and Thoroddsen2012; Agbaglah et al. Reference Agbaglah, Thoraval, Thoroddsen, Zhang, Fezzaa and Deegan2015; Josserand et al. Reference Josserand and Thoroddsen2016; Thoraval et al. Reference Thoraval, Li and Thoroddsen2016; Deka et al. Reference Deka, Ray, Biswas, Dalal, Tsai and Wang2017; Jian et al. Reference Jian, Channa, Kherbeche, Chizari, Thoroddsen and Thoraval2020; Fudge et al. Reference Fudge, Cimpeanu and Castrejón-Pita2021, Reference Fudge, Cimpeanu, Antkowiak, Castrejón-Pita and Castrejón-Pita2023; Sykes et al. Reference Sykes, Cimpeanu, Fudge, Castrejón-Pita and Castrejón-Pita2023; Osama, Deegan & Agbaglah Reference Osama, Deegan and Agbaglah2024). It is important to emphasise that the main objective of this study is to perform a comprehensive study of drop-impact dynamics across a wide range of parameter space, focusing on those phenomena that exhibit approximate axisymmetry. The validation conducted in § 3.1 has demonstrated its ability to reproduce key splashing features in the process of high-speed drop impact, which is sufficient to support the systematic study conducted in the following sections.
4. Results and discussion
The splashing dynamics of drop impact, as described in § 2.1, depends primarily on several dimensionless parameters. In this parametric study, we maintain the values of
$\rho _l/\rho _g=783$
,
$\mu _l/\mu _g=56$
,
$d_h=4.3$
mm,
$d_v=3.8$
mm and
$U_0=7.2$
m s
$^{-1}$
constant. The effect of dimensionless Reynolds and Weber numbers are explored for
$1000\leqslant {\textit{Re}}\leqslant 30\,000$
and
$100\leqslant We\leqslant 3000$
by systematically varying the liquid viscosity
$\mu _l$
and surface tension
$\sigma$
. Through this variation, we investigate and analyse the diverse splashing behaviours exhibited during high-speed drop–pool impingement.
Regime map of drop impact in terms of dimensionless Reynolds (
${\textit{Re}}$
) and Weber (
${\textit{We}}$
) numbers, showing transitions between various splashing behaviours identified in the present investigation. The filled symbols represent BC formation, whereas the hollow symbols indicate that the crown does not enclose. The solid black lines determine the limit values of BC formation at different ranges of
${\textit{Re}}$
fitted based on the simple relation
$K={\textit{We}}\sqrt {{\textit{Re}}}$
:
$K1=4.9\times 10^4$
,
$K2=9.9\times 10^4$
and
$K3=14.7\times 10^4$
. The solid red lines
$A1$
and
$A2$
indicate the approximate limits between different jet behaviours. Insets show representative impact cases for each regime, with arrows pointing from the impact condition to the corresponding case:
coalescence
$\rightarrow$
${\textit{Re}}=2000$
,
${\textit{We}}=100$
;
ejecta
$\rightarrow$
${\textit{Re}}=2000$
,
${\textit{We}}=400$
;
$\lozenge$
ejecta & lamella
$\rightarrow$
${\textit{Re}}=4000$
,
${\textit{We}}=600$
;
$\bigtriangleup$
bumping & roll jet
$\rightarrow$
${\textit{Re}}=3500$
and
${\textit{We}}=700$
;
$\bigtriangledown$
bumping
$\rightarrow$
${\textit{Re}}=6500$
and
${\textit{We}}=700$
;
$\square$
irregular
$\rightarrow$
${\textit{Re}}=9500$
and
${\textit{We}}=700$
.
Figure 4 presents a comprehensive regime map illustrating various distinctive drop-impact dynamics observed in our simulations. In this map, we encompass both the early-time jet behaviours and overall evolution of the crater structure, which, to the best of the authors’ knowledge, should be first reported herein. Based on our axisymmetric simulations, a wide variety of repeatable behaviours of early jets (see figure 5) and modes of crater evolution (see figure 9) have been identified, underscoring the rich and complex nature of splashing phenomenon. The presented map is the outcome of conducting more than 800 simulations.
4.1. Early-time jet behaviours
Representative cases calculated under different impact conditions, showing various repeatable jet behaviours and vorticity structures.
$\omega$
is the vorticity value, where the red and blue colours represent clockwise and counterclockwise rotation, respectively. (a) Coalescence followed by a downward-moving vortex separation (
${\textit{Re}}=2000$
,
${\textit{We}}=100$
). (b) Combination of ejecta and lamella without vortex separation (
${\textit{Re}}=2500$
,
${\textit{We}}=300$
). (c) Vortex separation from the upper corner of the ejecta leads to a separated lamella (
${\textit{Re}}=4000$
,
${\textit{We}}=400$
). (d) Earlier one-sign vortex shedding without separated lamella (
${\textit{Re}}=6500$
,
${\textit{We}}=400$
). (e) Reconnection between ejecta and drop surface entraps a large toroidal air bubble, leading to a secondary roll jet (
${\textit{Re}}=3500$
,
${\textit{We}}=700$
). (f) Reconnection between ejecta and drop surface without roll jet, where the entrapped bubble sinks rapidly inside the pool (
${\textit{Re}}=6500$
,
${\textit{We}}=700$
). (g) Strong interactions between ejecta and free surfaces, resulting in repeated toroidal bubbles and vortex separation (
${\textit{Re}}=9500$
,
${\textit{We}}=700$
). See supplementary movie 2.
This section focuses on the early-time jet dynamics in the vicinity of the neck region formed between the downward-impinging drop and the receiving pool. This comprises several key splashing regimes depicted in figure 4, including
coalescence,
ejecta formation,
$\lozenge$
ejecta with a late-emerging lamella,
$\bigtriangledown$
bumping between ejecta and drop interface,
$\bigtriangleup$
bumping with secondary roll jet formation, and
$\square$
irregular splash with strong interfacial interactions. Figure 5 illustrates seven representative cases calculated under different impact conditions, showing various repeatable jet shapes and vorticity structures captured by our simulations (see supplementary movie 2).
At
${\textit{We}}\lt 150$
, the very high surface tension generally prevents the emergence of the ejecta jet and the drop coalesces gently with the substrate fluid (
coalescence regime in figure 4). As the drop sinks and spreads, a large vortex is subsequently detached from the upper corner of the connection and travels horizontally underneath the free surface, generating a capillary swell that propagates on the interface (see figure 5
a for a case with
${\textit{Re}}=2000$
,
${\textit{We}}=100$
). A similar ‘no-jet’ boundary at
${\textit{We}}\approx 200$
was also observed in previous experimental (Zhang et al. Reference Zhang, Toole, Fezzaa and Deegan2012) and numerical (Agbaglah et al. Reference Agbaglah, Thoraval, Thoroddsen, Zhang, Fezzaa and Deegan2015) investigations. For
${\textit{Re}}\lt 3000$
and
${\textit{We}}\gt 150$
, the formation of thin ejecta starts to be observed in the neck of the connection between two liquid masses (
ejecta regime in figure 4), as illustrated in figure 5(b) for a case with
${\textit{Re}}=2500$
and
${\textit{We}}=300$
.
Increasing the Reynolds number to
${\textit{Re}}\approx 2000{\sim}9000$
, intriguing jet behaviours are categorised. As shown in figure 5(c) (
${\textit{Re}}=4000$
,
${\textit{We}}=400$
), a vorticity structure of one sign is generated at the upper corner of the ejecta (
$\hat {t}=0.1656$
). However, compared with figure 5(b) (
${\textit{Re}}=2500$
,
${\textit{We}}=300$
) at the same time, the vorticity structure is stronger here and the base of the upper corner is no longer stable, thus transitioning the local dynamics of the ejecta root to a vortex-shedding regime. The detached vortex rings then propagate underneath the free surface in a similar mode to figure 5(a), pushing the ejecta root to climb up along the side of the impact drop. As a consequence, a thicker jet is gradually lifted by the resurfacing vortex due to relatively lower surface tension. The emergence of these two distinct jets has been previously investigated by both experimental (Zhang et al. Reference Zhang, Toole, Fezzaa and Deegan2012; Agbaglah et al. Reference Agbaglah, Thoraval, Thoroddsen, Zhang, Fezzaa and Deegan2015; Sykes et al. Reference Sykes, Cimpeanu, Fudge, Castrejón-Pita and Castrejón-Pita2023) and numerical (Agbaglah et al. Reference Agbaglah, Thoraval, Thoroddsen, Zhang, Fezzaa and Deegan2015; Sykes et al. Reference Sykes, Cimpeanu, Fudge, Castrejón-Pita and Castrejón-Pita2023) works. Following the definition of these authors, the early-emerging jet is hereafter named ‘ejecta’ and the late-emerging jet is named ‘lamella’. In figure 4, this phenomenon is grouped as the
$\lozenge$
ejecta & lamella regime. As for the lifetime of the initial ejecta, it could be entirely absorbed by the impact drop at relatively low
${\textit{We}}$
(Zhang et al. Reference Zhang, Toole, Fezzaa and Deegan2012) or torn into microdroplets by its outward-stretching leading edge (
$\hat {t}=0.333$
and 0.387 in figure 5
c). Figure 5(d) shows the jet dynamics at a slightly higher
${\textit{Re}}$
(
${\textit{Re}}=6500$
,
${\textit{We}}=400$
), where a ‘vortex separation’ event occurs earlier than its counterpart in figure 5(c) (
${\textit{Re}}=4000$
,
${\textit{We}}=400$
) and the formation of the lamella cannot be clearly observed. Instead, the neck could move rapidly along the radial direction in its early evolution, leaving the vortex ring behind the base of the ejecta.
For relatively higher Weber numbers (
${\textit{We}}\gt 700$
) at
${\textit{Re}}\approx 3000{\sim}9000$
, reconnection between the radial-stretching ejecta and the downward-impinging drop occurs as shown in figures 5(e) and 5(f). The uprising ejecta bends and reconnects with the drop surface at early times before vortex shedding, entrapping and shedding a toroidal bubble at its upper corner. Comparing these two cases, it can be found that the reconnection between ejecta and drop at higher
${\textit{Re}}$
(figure 5
f) occurs earlier than the case at lower
${\textit{Re}}$
(figure 5
e). As the angle between the impinging drop and the initial quiescent pool increases with time, the entrapped bubble should become smaller with increasing
${\textit{Re}}$
. Indeed, at higher
${\textit{Re}}$
(figure 5
f), a new neck is established immediately upon the detachment of the toroidal bubble and expands outwards rapidly because of the sharper surface geometry at a small angle between drop and pool, leaving the vortex structure rotating and sinking along the drop–pool interface. In contrast, the bubble is detached at a larger angle at lower
${\textit{Re}}$
(figure 5
e), where the new point of connection propagates primarily in the vertical direction along the drop surface, and the bubble travels quickly to the outer side of the ejecta base (
$\hat {t}=0.333$
) due to pronounced pressure gradient across its curved surface. Accordingly, the motion of the entrapped toroidal bubble leads to different eventual jet dynamics regimes in figure 4: (i)
$\bigtriangleup$
bumping & roll jet: the toroidal bubble travels to the outer side of the ejecta base, the bubble will resurface and interact with the free surface of the pool, generating a secondary roll jet, and the original ejecta will be completely disintegrated; (ii)
$\bigtriangledown$
bumping: the bubble sinks quickly enough and stays at the inner side of the ejecta base, the remaining part of the ejecta will keep extending and the original ejecta will not be replaced. A rough limit between these two types of phenomenon is given at
${\textit{Re}}\approx 4000$
in figure 4.
In the phenomenon mentioned earlier, the distinction from one to two jets in drop impact has brought attention to a previously overlooked aspect of the phenomenon. According to Agbaglah et al. (Reference Agbaglah, Thoraval, Thoroddsen, Zhang, Fezzaa and Deegan2015), the emergence of a second jet is determined by whether the vortex detaches from the upper corner of the ejecta base in the parameter range
${\textit{Re}}\leqslant 5000$
and
${\textit{We}}\leqslant 900$
. These authors have found that the capillary number
$\textit{Ca}={\textit{We}}/{\textit{Re}}=\mu _lU_0/\sigma$
might be the most appropriate parameter to characterise this transition, and a transitional boundary at
$\textit{Ca}=0.2$
was concluded to represent the approximate limit of the vorticity-shedding regime. In our current study, a more realistic boundary (red line
$A1$
in figure 4) is introduced to differentiate between the one-jet regime (
ejecta) and the two-jet regime (
$\lozenge$
ejecta & lamella and
$\bigtriangleup$
bumping & roll jet) based on our numerical results, which looks quite close to the
$\textit{Ca}=0.2$
line, further verifying the validity of the previously established transitional limit.
In the highest range of Reynolds numbers (
${\textit{Re}}\gt 9000$
), the phenomenon of continuous vortex shedding becomes prominent (
$\square$
irregular splash in figure 4). This behaviour is exemplified in figure 5(g) for a case at
${\textit{Re}}=9500$
and
${\textit{We}}=700$
, where the vortex-shedding event occurs much earlier compared with other cases. This early onset of vortex shedding leads to oscillations of the ejecta base immediately after contact. Consequently, the thin ejecta is propelled to bend up and down, resulting in alternating collisions with the surfaces of drop and pool. In this process, toroidal bubbles of alternate signs are entrapped by the newly established neck points and only the leading parts of the ejecta survive and move outwards. As the impact angle increases over time, the bubble entrapment event takes longer to develop, causing the bubble to grow in size until the ejecta is no longer in contact with the interfaces. Eventually, a von-Kármán-type street is formed, characterised by the regular shedding of vortices along the drop–pool interface (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).
Regime map of drop impact in terms of dimensionless Reynolds (
${\textit{Re}}$
) and splashing (
$K={\textit{We}}\sqrt {{\textit{Re}}}$
) numbers. The numerical classifications are shown using the same colour convention as in figure 4. The experimental results reported in Thoraval et al. (Reference Thoraval, Takehara, Etoh, Popinet, Ray, Josserand, Zaleski and Thoroddsen2012) are plotted with different symbols:
, smooth ejecta sheet;
$\bigtriangleup$
, quartering;
$\bigtriangledown$
, bumping;
$\lozenge$
, protrusion;
$\square$
, irregular splashing. Overlap of distinctive splashing regimes can be found between numerical simulations and laboratory experiments. The wavy region indicates the area not explored in the present numerical simulations.
In some previous studies, a secondary dimensionless parameter
$K$
(splashing number), defined as
$K={\textit{We}}\sqrt {{\textit{Re}}}$
, is also often used to characterise the drop-impact phenomena (Deegan et al. Reference Deegan, Brunet and Eggers2007; Thoroddsen et al. Reference Thoroddsen, Thoraval, Takehara and Etoh2011; Thoraval et al. Reference Thoraval, Takehara, Etoh, Popinet, Ray, Josserand, Zaleski and Thoroddsen2012). Figure 6 reorganises the current numerical outcomes in terms of the
${\textit{Re}}$
and
$K$
parameters, and the experimental data from Thoraval et al. (Reference Thoraval, Takehara, Etoh, Popinet, Ray, Josserand, Zaleski and Thoroddsen2012) are incorporated into the plot using different symbols for reference. Despite some deviations between our numerical set-up and their experiments, superpositions of distinctive splashing behaviours can be clearly observed in these two datasets. Both the
ejecta (low
${\textit{Re}}$
) and
$\square$
irregular splash (high
${\textit{Re}}$
) regimes are found in good alignment between numerical simulations and laboratory experiments. At intermediate
${\textit{Re}}$
, the upward protrusion of the ejecta at lower
$K$
in the experiment is well captured in the present vortex detaching and resurfacing category (
$\lozenge$
ejecta & lamella), while the quartering and bumping classes at higher
$K$
fall within our toroidal bubble entrapment regimes (
$\bigtriangleup$
bumping & roll jet and
$\bigtriangledown$
bumping). These results underscore the concordance between our numerical regime map and existing experimental data, further enhancing the reliability of this comprehensive characterisation of drop-impact behaviour.
4.2. Effect of drop shape on transition boundary
Effect of drop shape variation on the key transition boundaries between different splashing regimes. Three drop shapes are considered: blue, prolate (
$\alpha = 1.17$
); red, spherical (
$\alpha = 1.0$
); black, oblate (
$\alpha = 0.88$
). The numerical classifications follow the same symbol conventions as in figure 4. The combined symbol
$\bigtriangleup$
$ \raise1.5pt\hbox{$\bigtriangledown$}$
represents a unified ‘bumping’ regime, encompassing both the
$\bigtriangleup$
bumping & roll jet and
$\bigtriangledown$
bumping regimes, to indicate the reconnection between the ejecta sheet and the drop surface.
As discussed in § 2.3, another key parameter influencing splashing dynamics is the drop shape at the moment of impact. The effect of drop shape on jet formation and bubble entrapment has been highlighted in some previous studies (Thoraval et al. Reference Thoraval, Takehara, Etoh and Thoroddsen2013, Reference Thoraval, Li and Thoroddsen2016; Li et al. Reference Li, Thoraval, Marston and Thoroddsen2018), indicating its potential to shift regime boundaries in splashing behaviour. Following the approach of Thoraval et al. (Reference Thoraval, Takehara, Etoh and Thoroddsen2013), the effective drop diameter is defined as
$d = (d_v d_h^2)^{1/3}$
, and the drop shape is characterised by the aspect ratio
$\alpha = d_v / d_h$
, where
$d_v$
and
$d_h$
denote the vertical and horizontal diameters, respectively. Figure 7 presents additional simulations conducted with varying aspect ratios
$\alpha$
for
$1000 \leqslant {\textit{Re}} \leqslant 10\,000$
and
$200 \leqslant We \leqslant 1500$
, spanning the key regime boundaries while keeping the effective drop diameter constant. Three drop shapes were examined: prolate (
$\alpha = 1.17$
), spherical (
$\alpha = 1.00$
) and oblate (
$\alpha = 0.88$
). These simulations suggest that drop shape plays a subtle yet non-negligible role in influencing the transition boundaries. At high Reynolds numbers (
${\textit{Re}} \gt 6000$
), a systematic rightward shift of the
$\square$
irregular splash transition boundary is observed as
$\alpha$
increases from 0.88 (black) to 1.00 (red). For the prolate drop (
$\alpha = 1.17$
, blue), this trend continues, further suppressing the occurrence of
$\square$
irregular splash in the lower-right region of the regime map and expanding the
$\bigtriangleup \raise1.5pt\hbox{$\bigtriangledown$}$
bumping behaviour toward lower
${\textit{We}}$
. Figure 8(c) presents an example of high-Reynolds-number impact dynamics at
${\textit{Re}} = 10000$
and
${\textit{We}} = 800$
for drops with varying aspect ratios. As
$\alpha$
decreases, the drop becomes flatter, reducing the inclination angle between the drop surface and the quiescent pool. This flattening enhances the neck curvature, increasing the likelihood of rapidly outward-expanding ejecta root instability. As a result, the neck of the connection becomes increasingly destabilised and the ejecta are more prone to swing, leading to intricate neck dynamics (e.g. von-Kármán-type vortex street) and frequent reconnection with surrounding free surfaces (e.g. chains of bubble rings). Similar behaviours were also reported by Thoraval et al. (Reference Thoraval, Takehara, Etoh and Thoroddsen2013) for impacts at
${\textit{Re}} = 12\,900$
, where numerical simulations, supported by experimental observations, demonstrated that oblate drops, especially those with flatter bases, tend to produce more complex neck dynamics and pronounced bubble-ring entrapment during the early stages of impact.
Early-time impact dynamics for drops with varying aspect ratios at
$\alpha =1.17, 1.00, 0.88$
, while maintaining a constant effective drop diameter. (a) Gas–liquid interface evolution over time (left to right) at
${\textit{Re}}=2500$
,
${\textit{We}}=800$
. (b) Neck dynamics at
$\hat {t} = 0.2376$
for
${\textit{Re}}=2500$
,
${\textit{We}}=800$
. (c) Neck dynamics at
$\hat {t}=0.171$
for
${\textit{Re}}=10\,000$
,
${\textit{We}}=800$
. In the top panels of (b) and (c), red and blue denote liquid originating from the drop and the pool, respectively. In the bottom panels of (b) and (c), red and blue indicate clockwise and counterclockwise vorticity in the zoomed regions. Insets show the initial drop shape (not to scale).
On the other hand, in the lower-Reynolds-number range, we observe a leftward shift of both the
$\lozenge$
ejecta & lamella transition and the
$\bigtriangleup \raise1.5pt\hbox{$\bigtriangledown$}$
bumping boundary as
$\alpha$
increases (see figure 7). This trend indicates that a more prolate drop shape promotes the formation of a stronger vortex at the upper corner of the ejecta base, thereby facilitating both vortex shedding (associated with
$\lozenge$
ejecta & lamella) and reconnection events (associated with
$\bigtriangleup \raise1.5pt\hbox{$\bigtriangledown$}$
bumping). At low
${\textit{Re}}$
, the ejecta typically emerges thicker and advances more slowly compared with high-
${\textit{Re}}$
cases due to reduced inertial forces and greater viscous damping. Under such conditions, the redistribution of momentum caused by drop shape becomes more subtle, especially near the ejecta root. Prolate drops focus more of their impact momentum along the vertical axis, thus encouraging the base of the ejecta to turn upward more quickly. More critically, the elongated geometry of a prolate drop enhances the local curvature at the upper corner of the ejecta base as it rises. This sharper curvature brings the ejecta sheet into closer proximity with the descending drop surface, creating favourable conditions for the formation of a strong vortex at the ejecta base. Once generated, this vortex can either detach from the upper corner or exert sufficient force to drive the ejecta base into reconnection with the drop surface, especially at higher Weber numbers, resulting in bubble entrapment as illustrated in figure 8(a). In contrast, oblate drops, with their flatter impact bottom, distribute momentum more radially. This lateral momentum spreading accelerates the horizontal expansion of the ejecta base while delaying its vertical deflection. Although the root of the ejecta eventually transitions from horizontal to vertical growth due to the intrinsic dynamics of sheet propagation, this progression leads to distinct behaviours near the regime transition boundaries for oblate drops. As illustrated in figures 8(a) and 8(b), when the drop shape shifts from prolate to spherical at lower
${\textit{Re}}$
and higher
${\textit{We}}$
, the reconnection point between the ejecta and the drop surface moves outward in the radial direction, resulting in a larger entrapped toroidal bubble. For an oblate drop (
$\alpha = 0.88$
), the reduced vertical momentum and delayed deflection can even prevent the ejecta sheet from making contact with the descending drop surface altogether, thereby suppressing reconnection and shifting the
$\bigtriangleup \raise1.5pt\hbox{$\bigtriangledown$}$
bumping transition boundary further to the right.
These shape-dependent effects, though subtle, meaningfully influence the boundaries between splashing regimes, underscoring the importance of accounting for drop geometry in impact analyses and related applications. Given the difficulty in precisely controlling drop shape in experiments and limited understanding of its role in altering impact phenomenon and transitional dynamics, these numerical insights could provide a valuable foundation for future investigations and furthermore offer guidance for interpreting complex splashing behaviours and for optimising operating conditions in applications such as inkjet printing and additive manufacturing, where drop deformation is often inevitable and impactful.
4.3. From corolla to BC
Direct numerical simulation snapshots for
${\textit{Re}}=18\,000$
, illustrating different drop-impact dynamics and flow field under varied Weber numbers: (a)
${\textit{We}}=200$
, thick tongue with central Worthington jet; (b)
${\textit{We}}=800$
, crown with Worthington jet; (c)
${\textit{We}}=1300$
, semiclosed dome with Worthington jet; (d)
${\textit{We}}=2000$
, BC. For each snapshot, the left part shows the magnitude of axial velocity
$\widetilde {U_z}$
and the right part shows the magnitude of radial velocity
$\widetilde {U_r}$
, where the purple colour means positive and the green means negative. See supplementary movies 3–6.
Another interesting aspect of the drop impact process is the overall evolution of the crater-shaped structure. Figure 9 presents typical examples of different modes of crater evolution for a constant Reynolds number (
${\textit{Re}}=18\,000$
) observed in our simulations, with variations in the Weber number (see supplementary movies 3–6). The left part of each snapshot shows the magnitude of axial velocity
$\widetilde {U_z}$
and the right part shows the magnitude of radial velocity
$\widetilde {U_r}$
, where the purple colour means positive and the green means negative. In figure 9(a), where the liquid has the highest surface tension, the growth of the crown is greatly suppressed, resulting in the formation of a thick tongue (swell) and a U-shaped cavity (
$\hat {t}=5.4$
). Due to the substantially low inertial force to the surface force, the structure reaches its maximum position quickly after impact and undergoes contraction caused by the concentric propagation of capillary waves along the crater wall (
$\hat {t}=18$
). At higher
${\textit{We}}$
(figure 9
b), the crown is more developed and the structure maintains a longer expanding period before collapsing. Further increasing
${\textit{We}}$
, as shown in figure 9(c), a thinner liquid crown is captured and its upper parts bend at higher elevations. In the second frame of figure 9(c) (
$\hat {t}=45$
), where the cavity just reaches its maximum depth, it can be clearly observed that both
$\widetilde {U_z}$
and
$\widetilde {U_r}$
are negative in the upper crown, indicating that the crown is collapsing towards the impact axis. Notably, in this case, cavity contraction initiates from the shallowing of the cavity bottom rather than the propagation of capillary waves coming from the earlier collapsing crown. For the highest range of
${\textit{We}}$
, the crown ultimately necks in as illustrated in figure 9(d), encapsulating a large bubble and forming a central ‘bullet jet’, the same as the case experimentally validated in § 3.
Effect of Weber number on the main crater of drop impact at
${\textit{Re}}=18\,000$
. (a) Gas–liquid interface shape at
$\hat {t}=4.5$
and 15.3. See supplementary movie 7. (b) Time evolution of the axial velocity of the crown tip, obtained by differentiating a smoothed spline fit to its tracked position. (c) Time evolution of the radial velocity of the crown tip, obtained by differentiating a smoothed spline fit to its tracked position. (d) Trajectories of the crown rim. The square marks the start point at
$\hat {t}=0.9$
and the arrow points at the direction of the motion. (e) Time evolution of the cavity width (dashed lines) and depth (solid lines).
Figure 10(a) illustrates the effect of Weber number on the evolution of the main crater. Upon varying the liquid surface tension, a clear pattern emerges: lower surface tension (higher
${\textit{We}}$
) leads to a faster, higher and thinner crown (see supplementary movie 7). Under high
${\textit{We}}$
, the upper crown surges and expands rapidly, experiencing significant effects from the surrounding airflow (refer to figure 12
b and discussion in § 4.4), causing it to bend towards the impact axis. In contrast, surface tension minimally affects the cavity during its initial stage (
$\hat {t}=4.5$
), although its impact becomes more pronounced once the crown collapses (
$\hat {t}=15.3$
) (see also figure 10
e). This phenomenon is particularly evident in cases with smaller Weber numbers (
${\textit{We}}=200$
and 400), where propagation of capillary waves at later times reduces the dimensions of the cavity, leading to the collapse of the structure.
Figures 10(b) and 10(c) show the time evolution of the axial (
$\widetilde{V_z}$
) and radial (
$\widetilde{V_r}$
) velocities of the crown tip, obtained by differentiating a smoothed spline fit to its tracked position. The higher the
${\textit{We}}$
, the larger the
$\widetilde {V_z}$
. At low
${\textit{We}}$
, the stronger surface tension significantly hinders the axial rise of the crown. This can be observed by the earlier time point when
$\widetilde {V_z}$
approaches zero, where the crown reaches its maximum height and retracts shortly after impact. Conversely, the
$\widetilde {V_r}$
curves corresponding to higher
${\textit{We}}$
values (
${\textit{We}}=1300$
, 2000 and 2700) decrease more rapidly and subsequently become negative, signifying the inward-bending activities of its thinner upper crown. Figure 10(d) depicts the trajectories of the crown rim for different Weber numbers. For relatively low
${\textit{We}}$
, the crown grows similarly, retracting vertically after reaching its maximum height while expanding radially (see figure 9
a,b). However, for high
${\textit{We}}$
, as showcased in figure 9(d), the crown rises faster and grows thinner, which allows its upper part to warp towards the impact axis before reaching the maximum height, leading to the eventual BC closure. Interestingly, for a case with intermediate Weber number at
${\textit{We}}=1300$
(see figure 9
c), the crown, despite its inward-bending tendency, fails to enclose. Instead, it attains a maximum height and subsequently undergoes inward collapse. The negative values for both
$\widetilde {V_z}$
(red curve in figure 10
b) and
$\widetilde {V_r}$
(red curve in figure 10
c) in the later stage signify that the rim continues to slowly approach the
$z$
-axis as it descends, where it rejoins the shallowing cavity before complete merger. In figure 10(e), the temporal evolution of the dimensionless cavity width (dashed lines
$\widetilde {C_r}$
) and depth (solid lines
$\widetilde {C_d}$
) is presented. For smaller
${\textit{We}}$
, the effect of surface tension is more significant, preventing the further expansion of the cavity. Here we see the less smooth curves for lower
${\textit{We}}$
cases showing more pronounced effects of surface wave propagation on the cavity evolution.
The progressive effects of
${\textit{We}}$
on the crown evolution reveal the presence of a distinct threshold that marks the initiation of BC phenomenon. Pan et al. (Reference Pan, Cheng, Chou and Wang2008) first reported a non-monotonic transitional boundary in terms of
${\textit{We}}$
and dimensionless film thickness
$\widetilde {h}=H/d$
for the onset of bubble encapsulation during water drop impacting thin water films. They observed a decrease in the transitional
${\textit{We}}$
with an increase in
$\widetilde {h}$
when the film thickness is smaller than the crater size. For impact on deep pools (
$\widetilde {h}\geqslant 2$
), a constant characteristic boundary
${\textit{We}}\approx 2570$
was identified for the BC phenomenon. By compiling previous experimental data in the literature, Murphy et al. (Reference Murphy, Li, d’Albignac, Morra and Katz2015) concluded the occurrence of BC at
${\textit{We}}\gt 2000$
for low-viscosity liquids. A similar value of
${\textit{We}}\approx 2450$
was later obtained by Sochan et al. (Reference Sochan, Beczek, Mazur, Ryżak and Bieganowski2018) in their study on splashing shapes in different liquid systems. In addition, these authors observed that the transitional
${\textit{We}}$
limit in high-viscosity liquid may become significantly larger than in low-viscosity liquid, highlighting the influential role of liquid viscosity in the closure phenomenon. However, due to the limited range of examined liquids (petrol, water and diesel fuel), detailed effects and dependence of transitional
${\textit{We}}$
on
${\textit{Re}}$
were not able to be fully elucidated.
In our current numerical investigations, as shown in figure 4, BC formation is observed across all calculated ranges of
${\textit{Re}}$
when
${\textit{We}}$
surpasses a certain threshold, signifying the dominance of inertial forces over the surface tension forces. The Weber number are explored for
$100\leqslant We\leqslant 3000$
, corresponding to surface tension values in the range of
$\sigma \approx 72.14{\sim}2164.33\times 10^{-3}$
N m
$^{-1}$
in realistic experimental conditions. Generally, the BC starts to be found at
${\textit{We}}\gt 900$
. However, the onset of the ultimate liquid crown depends primarily on the final status of the leading jet, and the occurrence of BC shows great dependence on the types of early-time jet dynamics due to liquid viscosity variations, so they do not adhere to a single scaling relationship. Therefore, as indicated in figures 4 and 6, by employing the simple relation
$K={\textit{We}}\sqrt {{\textit{Re}}}$
, three transitional boundaries that indicate the onset of BC can be determined by fitting our numerical results. Line
$K1=4.9\times 10^4$
separates the BC phenomenon at low
${\textit{Re}}$
within the
ejecta regime (figure 5
b), line
$K2=9.9\times 10^4$
indicates the onset of BC during
$\bigtriangledown$
bumping regime (figure 5
f), and finally line
$K3=14.7\times 10^4$
signifies the occurrence of canopy for the most energetic
$\square$
irregular splash regime (figure 5
g).
Effect of Reynolds number on the main crater of drop impact at
${\textit{We}}=2000$
within the
$\bigcirc$
ejecta early-time splashing regime. (a) Gas–liquid interface shape at
$\hat {t}=7.65$
. (b) Time evolution of the radius of the upper crown (
$\widetilde {R_r}$
).
Figure 11(a) shows the effect of varying Reynolds number
${\textit{Re}}$
on the impact phenomenon. Here we specifically focus on the
ejecta regime (
${\textit{Re}}\lt 3000$
) without involving complex early-time splashing interactions between the ejecta and the free surfaces. At higher
${\textit{Re}}$
, the crown rises faster, and its upper part subsequently bends inward toward the impact symmetry, resulting in an earlier surface closure as shown in figure 11(b). In contrast, at lower values of
${\textit{Re}}$
, viscous forces increasingly counteract liquid inertia, resulting in reduced jet rising velocities. The enhanced viscous effect also helps the crown resist the suction effect due to gas rushing (see discussion in the following § 4.4), delaying the bending of the upper portion. As a result, at lower
${\textit{Re}}$
(with higher viscosity), the crown takes longer to converge toward the symmetry axis, rises higher, and traps a larger gas bubble volume. This behaviour potentially explains the dependence of the bubble closure
${\textit{We}}$
threshold on
${\textit{Re}}$
(corresponding to line
$K1$
in figure 4), corroborating prior experimental findings presented by Murphy et al. (Reference Murphy, Li, d’Albignac, Morra and Katz2015) and Sochan et al. (Reference Sochan, Beczek, Mazur, Ryżak and Bieganowski2018) as referenced earlier.
It is noteworthy that the viscosity ratio between the gas and liquid phases is kept constant at
$\rho _l/\rho _g=56$
in the present simulations, meaning that variations in the Reynolds number
${\textit{Re}}$
alter both the liquid and gas viscosities simultaneously. This differs from typical experimental conditions in the air, where the gas viscosity remains fixed while only the liquid viscosity changes. In our study, the Reynolds number range is
$1000\leqslant {\textit{Re}}\leqslant 30\,000$
, corresponding to liquid viscosities in the range of
$\mu _l\approx 1.00\times 10^{-3}\sim 30.1\times 10^{-3}$
N s m
$^{-2}$
and gas viscosities of
$\mu _g\approx 17.9\times 10^{-6}\sim 536.8\times 10^{-6}$
N s m
$^{-2}$
. Under typical impact conditions in the air (
$\mu _g=18.1\times 10^{-6}$
N s m
$^{-2}$
), the viscosity ratio would be in the range
$\rho _l/\rho _g=55{\sim}1660$
, making the ambient gas in our simulations effectively more viscous at lower
${\textit{Re}}$
. As a result, the increased gas viscosity introduces additional damping effects that could influence crown expansion, closure dynamics and bubble entrapment. Specifically, higher gas viscosity may increase the drag resistance from the surrounding gas and suppress the growth of the crown, leading to a more confined structure compared with realistic air-based experiments. In three dimensions, the enhanced gas viscosity may dampen instabilities along the crown rim, leading to a smoother, less fragmented transverse rim structure. Additionally, the increased viscous resistance in the gas phase can potentially delay crown closure, eventually resulting in larger entrapped bubbles than those observed experimentally. While the present approach ensures numerical consistency and stability across different
${\textit{Re}}$
cases, it is important to recognise that these effects may introduce discrepancies when directly comparing with experiments conducted in the air. Nevertheless, the overall trends observed in our simulations remain valuable for understanding the fundamental mechanisms governing crown formation and closure under varying viscosity conditions. Future studies could explore cases with a fixed gas viscosity to provide a more direct comparison with experimental results.
4.4. Mechanism of surface closure
Surface-sealing phenomena are commonly observed in processes such as water entry of projectiles (Marston et al. Reference Marston, Truscott, Speirs, Mansoor and Thoroddsen2016b ; Kroeze, Rivas & Quetzeri-Santiago Reference Kroeze, Rivas and Quetzeri-Santiago2023) and oscillating bubble explosion near an interface (Li et al. Reference Li, Zhang, Wang, Li and Liu2019; Rosselló et al. Reference Rosselló, Reese and Ohl2022), where a splash curtain is forcefully ejected outwards. While several forces, including inertia, surface tension and ambient pressure, have been proposed as driving factors for surface closure (Clanet Reference Clanet2007; Aristoff & Bush Reference Aristoff and Bush2009; Marston et al. Reference Marston, Mansoor, Thoroddsen and Truscott2016a , Reference Marston, Truscott, Speirs, Mansoor and Thoroddsenb ; Eshraghi, Jung & Vlachos Reference Eshraghi, Jung and Vlachos2020; Wang et al. Reference Wang, Wang, Du, Wang, Wang and Chenguang2022), the intricate interplay among these forces and the dominant mechanism remain a subject of controversy.
Details of airflow information involved in (a)
${\textit{Re}}=18\,000$
and
${\textit{We}}=800$
, open corolla and (b)
${\textit{Re}}=18\,000$
and
${\textit{We}}=2000$
, BC formation. The left part shows pressure amplitude and the right part shows vorticity field and flow streamlines. The boundary between coloured and uncoloured regions (white) indicates the gas–liquid interface. The yellow arrow indicates vortex separation due to air rushing, and the green arrow points at the corresponding low-pressure region. The vorticity field is normalised by
$d/U_0$
and the pressure field is scaled by the initial dynamic pressure of the impact drop
$P_0=(\rho _lU_0^2)/2$
.
To understand the mechanism of surface closure during drop impact, we also need to turn our attention to the intriguing airflow dynamics around the crown. Figure 12 presents detailed airflow information during drop impact, providing numerical insight into the surface closure mechanism. As the liquid crown rises from the surface, external gas rushes into the expanding cavity, causing a reduction in internal pressure (
$\hat {t}=1.8$
). Meanwhile, the crown’s outward expansion is met with increasing hydrodynamic resistance from the liquid pool, which restricts further lateral displacement. The resulting force balance redirects the crown rim upward, leading to a transition from horizontal to vertical motion (
$\hat {t}=7.2$
). As the liquid film continues to expand under the action of inertia, a large vortex (yellow arrow) is generated due to gas rushing and subsequently detaches from the crown rim. This, in turn, results in a low-pressure region within the gradually growing vortex, as indicated by the green arrow at
$\hat {t}=12.6$
. For higher
${\textit{We}}$
, as shown in figure 12(b), the substantial pressure difference between the interior and exterior sides of the crown further compels its thinner film to bend inwards towards the symmetry plane (
$\hat {t}=18$
), leading to the eventual bubble closure event (
$\hat {t}=25.2$
). During this process, one could expect that the interplay of inertia, surface tension and aerodynamic pressure would significantly influence the behaviour of the thin-walled crown. For a case with lower
${\textit{We}}$
in figure 12(a), qualitatively similar airflow profiles emerge during its early stages. However, the relatively higher surface tension mitigates the vertical extension of splashes while encouraging the accumulation of liquid at its upper rim, giving rise to a thicker and shorter crown. In these instances, the pressure difference is less dominant in drawing the crown bending, although a slight tendency can be still observed. Instead of surface closure, this crown reaches its maximum position first and immediately undergoes a contraction phase (see also figure 9
b), leading to the collapse of the structure.
(a) Representative snapshot of velocity field and vorticity structure during high-speed drop impact (
${\textit{Re}}=18\,000$
,
${\textit{We}}=2000$
). The insert focuses on the flow details near the crown rim. The purple arrow points at the local airflow vorticity peak (
$\varOmega$
). (b) Time evolution of the airflow vorticity peak
$\varOmega$
recorded near the crown rim. The pressure profiles along the axial (symmetry) line at different time instants are plotted for cases (c)
${\textit{We}}=800$
and (d)
${\textit{We}}=2000$
. The dimensionless time
$\hat {t}$
is marked along each curve. The red curve indicates the locations of rim height
$\widetilde {R_h}$
.
Figure 13(a) shows a representative frame captured during high-speed drop impact, illustrating the inward-bending movement of the free-standing liquid film. As the crown expands, a large vorticity patch is generated around the rim interface due to gas shearing. Figure 13(b) displays the recorded time evolution of the airflow positive vorticity peak (purple arrow in figure 13
a) for the two distinct cases presented in figure 12. Initially, the vorticity peak magnitude appears significantly stronger for the case with reduced surface tension at
${\textit{We}}=2000$
. As the structure rises further, its rim decelerates (see figure 10
b), and the vorticity is gradually mitigated. Notably, starting from
$\hat {t}\approx 10$
, the upper crown undergoes rapid inward shrinkage and a high-velocity field is formed near the surrounding gas, resulting in an inflection point for the vorticity peak with a minimum value before surface sealing. In contrast, for the case with higher surface tension at
${\textit{We}}=800$
, the rim starts to retract after reaching a maximum position, and the vorticity weakens and vanishes as the structure collapses. In figures 13(c) and 13(d), we plot their pressure profiles along the axial line, with the position of the rim height
$\widetilde {R_h}$
indicated by red curves. Over time, a rapid drop in aerodynamic pressure can be caused due to surface shrinking before surface closure, while the pressure exhibits a gradual increase in the case with higher surface tension.
These results demonstrate the significant role of aerodynamics in the surface-sealing process during drop impact. Further endeavours should focus on advancing our understanding of this phenomenon under varied liquid properties and ambient environments, unravelling the dominant mechanisms governing the transition from liquid curtain expansion to surface sealing. Particularly, for further delineating the transitional boundary between non-closure and closure, combined parameters that integrate both hydrodynamic and aerodynamic effects will be crucial. This aspect remains to be explored in future work.
5. Conclusion
In this work, we have conducted a thorough exploration of high-speed drop impact dynamics through a series of axisymmetric simulations. Through careful comparisons with existing experimental data, our simulation successfully reproduced most of the phenomena captured by camera, showing its ability and effectiveness to reproduce key splashing features. Our meticulous observations have led to the identification of various repeatable early-time splashing patterns within the jets, vortex and overall crater structures, revealing the intricate nature of high-speed splash. We show that, under increased
${\textit{We}}$
, the crown exhibits accelerated growth, reaches greater heights and becomes thinner. This phenomenon facilitates the bending of the crown’s upper portion toward the impact axis under further reduced surface tension, leading to the ultimate crown closure event for the highest
${\textit{We}}$
range. This observation confirmed the existence of a corolla–canopy transition driven by the combined effects of surface tension and aerodynamic pressure, highlighting the subtle interplay of forces at play. We also demonstrated the crucial role of the delicate balance between inertia, surface tension and aerodynamic forces in governing the surface-sealing process on the thin-walled crown. The comprehensive investigation culminated in the development of a regime map for high-speed drop impact splashes, which encompasses both the early jet dynamics and the overall crater evolution, thereby providing a profound and detailed understanding of this intriguing phenomenon.
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2025.10543.
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.
Funding
This work is supported by the scholarship from the China Scholarship Council (CSC) under grant CSC no. 201908320462.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Effect of spatial resolution
In this appendix, we examine the effects of grid resolution on the numerical outcomes of drop impact in axisymmetric configurations. The numerical set-ups and impact conditions employed in this analysis are identical to the validated case outlined in § 3 (
${\textit{Re}}=30\,060$
,
${\textit{We}}=2964$
). Here we vary the maximum mesh refinement level, denoted as
$L_{\textit{max}}$
, from 11 to 17 by utilising the AMR algorithms implemented in the Basilisk code.
A.1. Early-time splashing phenomenon
Early-time splashing behaviours captured under different spatial resolutions in axisymmetric configurations (
${\textit{Re}}=30\,060$
,
${\textit{We}}=2964$
). The snapshots of the VOF two-phase field are demonstrated for (a)
$L_{\textit{max}}=12$
, (b)
$L_{\textit{max}}=13$
, (c)
$L_{\textit{max}}=14$
, (d)
$L_{\textit{max}}=15$
, (e)
$L_{\textit{max}}=16$
and (f)
$L_{\textit{max}}=17$
. Higher maximum refinement levels (d, e, f) show ‘irregular splash’ and alternate bubble entrapment near the neck region, whereas lower levels (a,b) can only capture the emergence of a smooth ejecta. The results indicate that (c)
$L_{\textit{max}}=14$
is the resolution limit where the large angle interaction between the uprising ejecta and the downward-moving drop is first captured (red arrow).
Figure 14 shows the evolution of the VOF two-phase field with time simulated under different spatial resolutions. Two distinct splashing modes are observed when comparing the dynamics near the neck region across different resolutions. At finer resolutions (
$L_{\textit{max}}= 15, 16, 17$
), strong interactions occur between the radial-stretching jet and the surfaces of the drop and pool. In contrast, coarser resolutions (
$L_{\textit{max}}= 12, 13$
) capture only a rising ejecta without secondary jet impingement, due to relatively larger cell size. As for
$L_{\textit{max}}=14$
, the ejecta rises ‘smoothly’ in its early evolution and impinges subsequently to the surface of the downward-moving drop (red arrow at
$\hat {t}=0.225$
in figure 14
c), indicating that this is the resolution limit at which the reconnection activities between the ejecta and the drop can be first captured.
The observed approximate convergence of the primary splashing features with increasing resolution suggests that
$L_{\textit{max}}=15$
(1638 cells per drop diameter) is necessary to capture such high-energy phenomena of ‘irregular splash’ together with toroidal bubble entrapment. As
$L_{\textit{max}}$
increases, the initial ejecta appears earlier and thinner, resulting in even earlier occurrence of alternate bubble entrapment along the drop–pool boundaries, as shown in figure 14(f). Similar patterns of alterations in the splashing morphology with varying maximum refinement levels have been also observed in our previous three-dimensional simulations in Wang et al. (Reference Wang, Liu, Bayeul-Lainé, Murphy, Katz and Coutier-Delgosha2023).
It is noteworthy that physical processes of irregular splash and bubble entrapment have been also numerically observed by Thoraval et al. (Reference Thoraval, Takehara, Etoh, Popinet, Ray, Josserand, Zaleski and Thoroddsen2012) and later been confirmed by Castrejón-Pita et al. (Reference Castrejón-Pita, Castrejón-Pita and Hutchings2012) using advanced experimental techniques. In the numerical study of Thoraval et al. (Reference Thoraval, Takehara, Etoh, Popinet, Ray, Josserand, Zaleski and Thoroddsen2012, Reference Thoraval, Takehara, Etoh and Thoroddsen2013), they found that irregular splash and bubble (ring) entrapment occur at relatively high Reynolds numbers (
${\textit{Re}}\gt 12\,000$
). Despite the smaller range of their study in comparison with our present case, their results suggest that more complex jet dynamics and bubble entrapment can be even expected at higher
${\textit{Re}}$
. Furthermore, using the Gerris code (ancestor of the current employed Basilisk code), these authors concluded that a minimum cell size of more than 1000 times smaller than the drop diameter is essential to capture the jet impingement activities (bumping), which agrees well with the resolution limit at
$L_{\textit{max}}=14$
(820 cells per diameter) observed in this investigation (figure 14
c). However, for observing the intricate irregular splash and bubble entrapment occurring at the highest range of
${\textit{Re}}$
and
${\textit{We}}$
numbers, finer spatial resolution is suggested by both our (figure 14
d,e,f) and their (their supplemental material) convergence studies.
A.2. Crater shape and closure event
Comparisons of numerical results under different maximum mesh refinement levels in axisymmetric configurations. Panels (a) and (b) compare the shapes of air–water interface at
$\hat {t}=0.693$
and
$\hat {t}=1.701$
. The black arrow points at the breakup of the crown tip with the highest resolution. (c) Closure time of the upper part of the crown
$\hat {t}_c$
. (d) Volume of the entrapped large bubble
$V_b$
at the moment of canopy enclosure, where
$V_0$
is the initial volume of the impact drop.
Figures 15(a) and 15(b) compare the shapes of the air–water interface at
$\hat {t}=0.693$
and
$\hat {t}=1.701$
calculated at different maximum refinement levels. At low resolutions (
$L_{\textit{max}}=12, 13$
), the ejecta rises smoothly from the neck, resulting in an upward-expanding crown rim. At
$L_{\textit{max}}=14$
, the evolution of the crown is significantly altered by the entrapped toroidal air bubble (see also figure 14
c). However, for higher resolutions (
$L_{\textit{max}}=15, 16$
and 17), the expansion of the crown rim is upwards and outwards, which is consistent with previous experimental observations in Murphy et al. (Reference Murphy, Li, d’Albignac, Morra and Katz2015), as discussed in § 3. It should be noted that the thickness of the early jet can be very thin on its tip and the liquid layer might be represented by just several cells, which makes the shape of the liquid layer close to the tip very sensitive to any changes in numerical parameters. This effect can cause small variations at the tip of the crown (black arrow in figure 15
a), but the main shape of the crown is well converged at
$L_{\textit{max}}=15, 16$
and 17.
Figure 15(c) shows the closure time of the upper part of the crown
$\hat {t}_c$
and figure 15(d) shows the volume of the entrapped large bubble
$V_b$
at the moment of canopy enclosure, calculated under different maximum refinement levels. Large variations of the closure event are observed when spatial resolution is insufficient. For cases at
$L_{\textit{max}}=15, 16$
and 17, the results are very close to each other.
A.3. Energetics
Now we discuss the numerical convergence with respect to energetics. The kinetic energy
$E_k$
and the gravitational potential energy
$E_g$
are calculated as follows:
\begin{equation} E_k=\frac {1}{2}\displaystyle \int _{V}\rho \|\boldsymbol{U}\|^2{\rm d}V,\quad E_g= \displaystyle \int _{V}\rho \boldsymbol{g}z{\rm d}V-E_{g0}. \end{equation}
where
$g$
is the gravitational acceleration and
$Z$
is the vertical distance to the bottom of the computational domain. The integrals are evaluated over the entire liquid phase. Here,
$E_{g0}$
is the gravitational potential energy of the liquid phase at
$\hat {t}=0$
. The total mechanical energy of the liquid phase is therefore calculated as
$E_m=E_k+E_g$
. Figure 16(a) plots the time evolution of
$E_k$
,
$E_g$
and
$E_m$
calculated under different maximum refinement levels. The results at
$L_{\textit{max}}=15, 16$
and 17 are well converged, whereas the results at lower levels show large variations.
Time evolution of energy aspects calculated under different maximum mesh refinement levels in axisymmetric configurations. (a) Time evolution of kinetic energy
$E_k$
(dotted), gravitational potential energy
$E_g$
(dashed) and the total mechanical energy
$E_m=E_k+E_g$
(solid). In the vertical axis, various energy aspect
$E$
is rescaled by the initial kinetic energy
$E_0$
in the domain. (b) Time evolution of instantaneous dissipation rate
$\epsilon$
collected during simulations.
Figure 16(b) shows the time evolution of the instantaneous dissipation rate calculated directly using the form
\begin{equation} \epsilon =\displaystyle \int _{V}\mu \frac {\partial u_i}{\partial x_j}\frac {\partial u_j}{\partial x_i}{\rm d}V. \end{equation}
where
$u_i$
and
$u_j$
are the velocity components,
$x_i$
and
$x_j$
are the 2D spacial components. Comparing all the curves, we see different dissipating patterns between ‘smooth ejecta’ and ‘irregular splash’. For
$L_{\textit{max}}=15, 16$
and 17 cases, the data exhibit the same qualitative features and magnitudes of early-time dissipation rate, reflecting the rapid transfer of kinetic energy during the chaotic process of the irregular splash. The calculated dissipation rate is generally higher as
$L_{\textit{max}}$
increases, which is reasonable since Basilisk is not energy-conserving by construction and the direct calculation of dissipation rate based on the deformation tensor is greatly dependent on grid resolutions (see also discussion in Mostert & Deike (Reference Mostert and Deike2020)).
Appendix B. Effect of computational domain size
Effect of computational domain size on the cavity depth evolution for drop impact at
${\textit{Re}}=30\,060$
and
${\textit{We}}=2964$
.
Figure 17 shows the effect of computational domain size on the cavity depth evolution for a case at
${\textit{Re}}=30\,060$
and
${\textit{We}}=2964$
. In this comparison, the actual domain size
$L=20d$
employed in the main text and a further increased domain
$L=40d$
are examined. While the domain size is doubled, the maximum refinement level is also adjusted from
$L_{\textit{max}}=15$
for
$L=20d$
to
$L_{\textit{max}}=16$
for
$L=40d$
, in order to maintain the same minimum resolved scale. The two curves are almost identical, except for minor deviations at later stages after surface closure. These slight discrepancies appear in thin flow structures such as crown tip, central jet and surface capillary waves due to the dynamic mesh refinement using AMR, which does not interfere with the overall cavity expansion. The curves are similar to those in figure 3(f), capturing the same physics and showing no significant difference between the two domain sizes. In both simulations, the maximum cavity depth is reached at
$\hat {t}\approx 40$
, corresponding to
$33\,\%$
and
$16.5\,\%$
of the pool depth for
$L=20d$
and
$L=40d$
, respectively. This analysis confirms that the chosen domain size of
$L=20d$
is sufficient for accurately modelling the physical dynamics of the expanding cavity in the present study.
Open access
