2.1. Main features of the solver

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

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

(2.2)$$\begin{gather}\rho\left(\frac{\partial \boldsymbol{U}}{\partial t}+(\boldsymbol{U}\boldsymbol{\cdot}\boldsymbol{\nabla})\boldsymbol{U}\right)=-\boldsymbol{\nabla} P+\boldsymbol{\nabla}\boldsymbol{\cdot}(2\mu\boldsymbol{\mathsf{D}})+\rho\boldsymbol{a}+\sigma\kappa\delta_s\boldsymbol{n} , \end{gather}$$

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

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

The BASILISK framework is a flow solver developed by Popinet (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$.

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.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.

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.

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.

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.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.

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).

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).

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).

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.

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.

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.

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:

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

Figure 13 shows the distribution of drop liquid at different stages of impact. Shortly after contact ($t=0.5$ ms), the main part of the drop is still concentrated at the bottom of the ‘bulk cavity’ and radially spreading thin film extends from its edges. For the next few milliseconds, the drop flattens and expands rapidly into a thin liquid layer coated along the interior surface of the pool. Meanwhile, liquid threads are emitted from the rim of the drop liquid and then break up into small droplets on its tips ($t=3$ ms), which is comparable to the formation of azimuthal destabilization on the retracting flattened drop edge at lower impact energy reported in Lhuissier et al. (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.

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.

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

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

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

where $\alpha$ and $\zeta$ donate the dimensionless crater radius and axial coordinate of the centre of the sphere and the dimensionless penetration depth is expressed as $\alpha +\zeta$. As explained by Bisighini et al. (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.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.

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

5.2. Droplet statistics

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

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

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

The sensitivity of droplet statistics to changes in mesh resolution during sheet fragmentation has been discussed in Appendices A.1 and A.3, showing that large droplets are approximately numerically converged but small droplets close to the minimum cell size ($S_d\approx 2\varDelta \sim 4\varDelta$) are generally overpredicted. Similar effects have also been captured by other high-resolution simulations involving droplet production using BASILISK code: the mesh study in Pairetti et al. (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.

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.

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.