Fragmentation of colliding liquid rims

Abstract

We present direct numerical simulations of the splashing process between two cylindrical liquid rims. This belongs to a class of impact and collision problems with a wide range of applications in science and engineering, and motivated here by splashing of breaking ocean waves. Interfacial perturbations with a truncated white noise frequency profile are introduced to the rims before their collision, whose subsequent morphological development is simulated by solving the two-phase incompressible Navier–Stokes equation with the adaptive mesh refinement technique, within the Basilisk software environment. We first derive analytical solutions predicting the unsteady interfacial and velocity profiles of the expanding sheet forming between the two rims, and develop scaling laws for the evolution of the lamella rim under capillary deceleration. We then analyse the formation and growth of transverse ligaments ejected from the lamella rims, which we find to originate from the initial corrugated geometry of the perturbed rim surface. Novel scaling models are proposed for predicting the decay of the ligament number density due to the ongoing ligament merging phenomenon, and found to agree well with the numerical results presented here. The role of the mechanism in breaking waves is discussed further and necessary next steps in the problem are identified.

1. Introduction

Liquid atomisation is a class of challenging multiphase problems (Obenauf & Sojka 2021) featuring a large separation of scales and various interacting physical mechanisms, and is of significance to numerous fields of application including meteorology (Villermaux & Bossa 2009), ink-jet printing (Castrejón-Pita et al. 2015, 2021; Lohse 2022), internal combustion engines (Yarin 2006) and pharmaceutical manufacturing (Mehta et al. 2017). Within the context of air–sea interactions, liquid fragmentation is primarily associated with wave-breaking events, and gives rise to ocean sprays. These spray drops are then transported within the atmospheric boundary layer while exchanging with the latter moisture, momentum and heat during their lifetime; thus leaving their impact on both global and regional climates (Lhuissier & Villermaux 2012; Veron 2015; Deike 2022). Atomisation involves topological changes of the liquid bulk driven by external forces, typically followed by formation of corrugated ligaments subject to capillary breakup, and ends with a number of fragments featuring a broad size distribution, the knowledge of which is crucial for various areas of applications listed above (Villermaux 2007).

Among various types of liquid atomisation problems, the impact of liquid droplets has received much scholarly attention for its ubiquitous presence, rich dynamics and vast range of applications (Yarin 2006; Villermaux 2007; Cheng, Sun & Gordillo 2022) since the pioneering experimental study of Worthington (1877). The impact process features a competition between inertial and capillary forces, which, together with the characteristics of the impacting object and the surrounding gas phase, shapes the final outcome of the original droplet: bouncing, spreading or splashing. Empowered by the rapid pace of sensor developments and increasing computational capacities, past works have elucidated considerable details about the ephemeral kinematic and morphological development of drops impacting with various types of surfaces (Cheng et al. 2022), including liquid films (Thoraval et al. 2013), deep pools (Agbaglah et al. 2015; Wang et al. 2023), smooth solid surfaces (Wildeman et al. 2016; Cimpeanu & Papageorgiou 2018), rough solid surfaces with friction (García-Geijo et al. 2021) or an identical droplet (He, Xia & Zhang 2019; He, Yue & Zhang 2022). Some recent works have also probed the dynamic properties of drop impact including the distribution of impact force and stresses, providing an alternative approach to investigating the impact dynamics at early times (Cheng et al. 2022). Specifically, there have been a series of recent experimental works studying the high-speed impact of droplets with a surface of comparable size as a canonical unsteady fragmentation problem (Wang & Bourouiba 2017, 2018, 2021, 2022; Wang et al. 2018), providing valuable insights into various aspects of the impact process in the limit of large impact Weber number $We$, including the self-similar evolution of the liquid-phase thickness and velocity profile (Wang & Bourouiba 2017), the dynamics of the rim bordering the expanding liquid sheet (Wang et al. 2018), the growth of liquid ligaments and the detachment of liquid drops from their tips (Wang & Bourouiba 2018, 2021) and the partition of mass, momentum and energy during the entire collision process (Wang & Bourouiba 2022).

However, the impact-induced fragmentation of liquid bulks featuring non-spherical initial shapes is also attested to, which has received considerably less attention and remains less understood compared with drop impact problems (Liu et al. 2021). Among these is the collision of liquid rims, which has seen some recent investigations experimentally (Néel, Lhuissier & Villermaux 2020) and numerically (Agbaglah 2021) and is also the focus of the present work. In the two works cited above, an initially intact liquid film is perforated to form small holes, which then expand under surface tension and develop bordering rims travelling at the Taylor–Culick velocity (Taylor 1959; Culick 1960). Neighbouring film holes merge with one another when their bordering rims collide, oscillate and break up into small fragments, a process commonly observed during the rupture of films, which may be induced by rapid radial expansion of liquid shells (Vledouts et al. 2016) or the inflation of a liquid drop interacting with a surrounding airflow (Jackiw & Ashgriz 2022; Ling & Mahmood 2023; Tang, Adcock & Mostert 2023). Agbaglah (2021) placed the two holes immediately adjacent to each other so that the two liquid rims collide at low $We$ values and the fused liquid bridge is found to pinch off under oscillation and form only a few small droplets. Néel et al. (2020) were able to investigate the rim collision phenomena by varying the initial distance between the two perforation sites, thus varying the impact $We$ value within the range of $50 \leq We \leq 200$. Apart from the primary capillary breakup mechanism of fused liquid bridges as discussed by Agbaglah (2021), Néel et al. (2020) also identified a critical $We$ value of 66 beyond which the decelerating transverse lamella rims develop elongating ligaments under the Rayleigh–Taylor (RT) instability, which then produce many secondary fine drops under capillary instabilities featuring a skewed size distribution. Additionally, there have been some early theoretical analyses on the capillary-driven coalescence of two liquid cylinders (Hopper 1993a,b; Eggers, Lister & Stone 1999); however, these are conducted at the creeping flow limit with negligible liquid bulk velocity, and thus may not be directly applicable to the current problem featuring finite collision speeds.

Apart from its presence during the rupture of thin films, the collision of liquid rims is also of specific interest due to its strong resemblance to the secondary wave splashing phenomena observed in ocean wave-breaking events (Kiger & Duncan 2012; Mostert et al. 2022; Erinin et al. 2023a). Namely, two consecutive well-defined splashing phases have been identified (Kiger & Duncan 2012) during the lifetime of a deep-water plunging breaker. The first is the ‘forward splashing’ mechanism occurring right upon the reconnection of the overturning wave front and the sea surface below, which according to Mostert et al. (2022) only produces small amounts of droplets. After this, as shown in figure 1(d), the wave bulk catches up with the decelerated splash-up generated from the initial impact at a relative speed $u_r$ close to the phase speed of the unbroken wave $U_g$. At the indentation region between these two structures, the shape of the wave bulk and the splash-up can be approximated as two cylinders, whose cross-sections feature radii of curvature $r_b$ and $r_s$, which are typically different. The rapid closure of the indentation region leads to the ‘secondary splashing’ phenomenon, which is characterised by a wall of vertically projected small droplets along the transverse direction as reproduced in figure 1(a–c), accompanied by air entrainment within the former indentation region (Kiger & Duncan 2012). Under experimental conditions, fragments generated via this mechanism comprise approximately one third of the total number of droplets produced over the entire wave-breaking process (Erinin et al. 2023a); Mostert et al. (2022) found that this splashing mechanism produces many fragments and can be curbed by strong surface tension (or small Bond numbers $Bo$). While Wang et al. (2016) noted the connection between the corrugated surface of the splash-up and the vortical structures beneath the wave surface, to our knowledge no existing study has analysed the physical mechanism governing the formation of these fragments, and their contribution to the overall droplet distribution associated with wave breaking remains unknown (Andreas et al. 1995; Kiger & Duncan 2012; Veron 2015). Furthermore, this splashing mechanism is found to produce many fragments close to the minimum grid size of Mostert et al. (2022); together with the highly transient nature of wave breaking and the presence of other fragmentation mechanisms, this indicates considerable difficulty in investigating the secondary splashing phenomena directly within the context of wave breaking. While we do not yet reproduce the fragment statistics seen in the studies above, the present work serves as a first step towards understanding the more complex wave splashing phenomena by retaining the major generation mechanism of splash fragments while leaving out many complicating factors, including the size difference between the wave bulk and the splash-up evolving with time, and internal turbulent flow within the liquid phase as a similar approach taken by Gao, Deane & Shen (2021) reveals the connection between the bubble size distributions of destabilising air cylinders and air cavities entrained by plunging breakers.

Figure 1. (a–c) Wave splashing observed in previous numerical (a,b) and experimental (c) studies, adapted from Wang, Yang & Stern (2016) (a), Mostert, Popinet & Deike (2022) (b) and Erinin et al. (2023a) (c), respectively. (d) Sketch showing the ensemble-averaged breaking wave profile taken from Erinin et al. (2023b) after the initial moment of impact in the breaking wave, at the moment of secondary splashing, where dashed lines indicate our simplification of the problem as the collision of two cylindrical rims with radii $r_b$ and $r_s$. In this study we consider the basic case $r_s = r_b$.

In this work, we conduct a comprehensive investigation of the collision between two liquid cylinders with identical size, covering the entire deformation and fragmentation period. The direct comparison and establishment of connections between the rim collision results and the statistics of the secondary wave splashing phenomenon are left for future work, together with the role of gravity and the difference between the sizes of the wave bulk ($r_b$) and the splash-up ($r_s$) which complicate the early-time rim dynamics. Two-phase numerical simulations are conducted to derive detailed flow field information during this highly transient collision process. Our study is structured as follows. We first present in § 2.1 the problem configuration and the parameter space of the current work, and then introduce the numerical method in § 2.2. After providing an overview of the rim collision phenomena in § 3, we quantitatively analyse the development of each part of the expanding liquid bulk successively following a spatial order, namely the kinematics of the spreading liquid sheet (§ 4.1) and its bordering rim (§ 4.2), the growth and merging of transverse ligaments topping the rim (§ 5) and the statistics of fragments shed from the ligaments (§ 6). We conclude the study in § 7 with some remarks on future work.

2. Formulation and methodology

2.1. Problem description

The geometrical configuration for the rim collision problem is shown in figure 2(a), where two infinitely long cylindrical liquid rims with diameter $d_0$, density $\rho _l$ and viscosity $\mu _l$ are aligned along the $x$ axis, and set to travel along the $z$ axis with uniform velocities of opposite signs and the same magnitude $U_0$. The liquid cylinders are surrounded by an inert gas phase with density $\rho _g$ and viscosity $\mu _g$, and the liquid–gas interface is characterised by a surface tension coefficient $\sigma$. It is noted that gravitational effects have been neglected in the current set-up; and differing from the configurations of Néel et al. (2020) and Agbaglah (2021), there is no interstitial film connecting the two approaching cylinders. Consequently, four non-dimensional controlling parameters can be written for this problem

(2.1a–d)\begin{equation} We \equiv \frac{\rho_l (2U_0)^2 d_0}{\sigma},\quad Oh \equiv \frac{\mu_l}{\sqrt{\rho_l d_0 \sigma}},\quad \rho^* \equiv \frac{\rho_l}{\rho_g},\quad \mu^* = \frac{\mu_l}{\mu_g}, \end{equation}

where $We$ and $Oh$ are, respectively, the Weber and Ohnesorge numbers comparing inertial and viscous effects with capillary forces, and $\rho ^*$ and $\mu ^*$ are, respectively, the density and viscosity ratios of the liquid and gas phases. In this work, $Oh$ is set as 0.01 in most of the simulations, whereas its influence on the fragment statistics is briefly discussed in Appendix A. Here, $\rho ^*$ and $\mu ^*$ are set as 830 and 55, respectively, which are typical values for an air–water system (Pairetti et al. 2018).

Figure 2. (a) Sketch showing the configuration of the liquid rim collision problem; (b) ensemble-averaged power density spectrum of the white noise signal for generating initial interface perturbations on the cylindrical rims.

We set the width $D$ of the cubic simulation domain as $10d_0$ to allow enough space for the morphological evolution of the coalesced liquid structure. Utilising the symmetry of the splashing phenomena about the $xz$-plane, we only model the merging of the upper halves of the two liquid cylinders to save computational resources. A symmetric boundary condition is therefore applied at the bottom, and an outflow boundary condition is imposed on the top boundary so that fragments produced from the collision can leave the domain from there at late time; the other boundary conditions are set as periodic.

To investigate the sensitivity of the fragmentation process to the initial conditions (Liu & Bothe 2016; Berny et al. 2022), and also taking into account the surface corrugation of the splash-up in wave-breaking events (Kiger & Duncan 2012) which still has not been quantified according to our knowledge, we introduce a random transverse perturbation within a certain wavelength range on the cross-sectional radius of the two cylindrical rims (Pal et al. 2021), in the form of filtered white noise signals characterised by the following two parameters:

(2.2a,b)\begin{equation} \varepsilon_0 \equiv \frac{2\epsilon_0}{d_0},\quad N_{max}, \end{equation}

where $\varepsilon _0$ is the non-dimensionalised characteristic amplitude of perturbation, and $N_{max}$ defines the highest wavenumber among the spectrum of the perturbation signal. The filtered white noise signal is the default type of initial interface perturbation we impose on the rims, as Zhang et al. (2010) did for analysing the linear stability of the crown splash. Occasionally, we also impose single-wavelength sinusoidal perturbations, or a superposition of sinusoidal perturbations with wavelengths $\lambda = D/8,\ D/16$ and $D/32$ for comparison, which will be explicitly denoted by ’Sing.’ and ’Sup.’ when reported. In the case of single-wavelength perturbations, $N_{max}$ corresponds to the perturbation wavenumber.

As discussed above, $We$, $\varepsilon _0$ and $N_{max}$ constitute the parameter space for the present study. Among these, $We$ is varied between 60 and 280 where the coalesced liquid bulk expands vertically to form a lamella, and $\varepsilon _0$ is set as 0.02, 0.04 and 0.06, within the limit of small radial perturbations; $N_{max}$ varies between 15 and 80, whose influence will be discussed in detail in § 5.2.

2.2. Numerical method

The open-source scientific computation toolbox Basilisk (Popinet 2019) is used in this work to solve the two-phase nonlinear, incompressible, variable-density Navier–Stokes equations. A second-order accurate discretisation is applied in both space and time, and a geometric volume-of-fluid (VOF) method in a momentum-conserving formulation is used to maintain a sharp representation of the liquid–gas interface while minimising the parasitic currents induced by surface tension (Popinet 2018; Tang et al. 2021). Capillary effects are modelled as source terms in the Navier–Stokes equations using an adaptation of Brackbill's method (Brackbill, Kothe & Zemach 1992; Popinet 2009), which calculates the interface curvature by taking the finite-difference discretisation of the derivatives of interface height functions (Popinet 2009). The octree-based adaptive mesh refinement scheme based on the estimation of local discretisation errors of the VOF function $f$ and flow velocity $\boldsymbol {u}$ is adopted so as to reduce the computational cost at high resolution levels $L$, which is defined using the minimum grid size

(2.3)\begin{equation} \varDelta = \frac{D}{2^L}. \end{equation}

The results presented in the main body of this work are obtained from three-dimensional simulations at $L = 10$, at which the late-time evolutions of the interface profile and liquid-phase energetics have reached grid independence. The numerical convergence of fragment statistics is also established for fragments with diameter larger than $4\varDelta _{10}$, which is discussed in detail in Appendix A. Results from some two-dimensional simulations are also presented for comparison with the three-dimensional rim dynamics (§ 4.2), and to investigate lamella foot formation at very early times (Appendix B).

Since the present study focuses on the liquid-phase morphological development after the rims begin to coalesce, the distance between the symmetric axes of the two cylindrical rims is set at initialisation as $\Delta z_c = 0.95d_0$, so that the slightly overlapped rims form a line of contact. The interfacial perturbations on the rims are introduced as follows. Firstly, discrete white noise signals with unit variance $\sigma ^2$ are produced using the random number generator provided in Basilisk. Figure 2(b) shows the ensemble-averaged power density spectra of the white noise signals generated in Basilisk at different numbers of realisations $N_{samp}$. It is observed that the power density $P(f)$ is close to the theoretical value of $\sigma ^2$ at all frequencies $f$, matching the requirement of frequency independence for white noise signals. Next, we apply a low-pass filter to these signals so that only the lowest $N_{max}$ wavelengths are preserved, and the filtered signal is normalised so that its standard deviation becomes $\epsilon _0 = 0.5\varepsilon d_0$. The normalised signals $\eta$ are then mapped onto the transverse radius profile $R(x)$ of the liquid rims, in the form of $R(x) = R_0 + \eta (x)$. The two liquid cylinders being positioned close to each other ensures that they coalesce immediately when the simulation starts and there is not sufficient time for capillary effects to smooth out the perturbations, or amplify them to trigger the Rayleigh–Plateau (RP) instability (Pal et al. 2021) so that the cylinders break up prematurely.

3. Overview of rim splashing

Here, we qualitatively describe the rim splashing process as observed in the simulations before analysing the detailed dynamics at each stage in the following sections. Figure 3 presents the isometric view of the splashing phenomenon at $We = 200$, $\varepsilon = 0.06$ and $N_{max} = 25$, whereas figure 4 shows the side view for a few different $(We,\varepsilon )$ configurations with $N_{max} = 25$. Tiny air bubbles are entrained within the indentation space between the two rims following the impact (Thoraval et al. 2013; Josserand, Ray & Zaleski 2016; Erinin et al. 2023a), which in our case have no known effects on the subsequent development of liquid-phase flow fields. The rims coalesce rapidly, causing a dramatic increase in the local liquid-phase pressure (Néel et al. 2020). A very large pressure gradient arises at the surface ‘indentation’ between the two cylinders since the air-phase pressure remains unchanged, leading to the vertical acceleration of fluids near the surface (Longuet-Higgins 2001). The magnitude of the vertical velocity near the contact line can be a few times larger than the initial horizontal velocity $U_0$, causing the contact line to advance upwards rapidly with fluid particles nearby converging to it, while the high-pressure region follows it closely rather than remaining on the axis of symmetry (Philippi, Lagrée & Antkowiak 2016). A thin transverse liquid lamella (Néel et al. 2020) appears, as shown in figure 3(b), and the location of the lamella foot agrees well with (B12) developed in Appendix B using potential flow theory (Riboux & Gordillo 2014). The liquid lamella is also observed in the experimental study of Goswami & Hardalupas (2023) where two neighbouring drops impact with a solid surface simultaneously and expand into contact. Within the parameter space of the present study, this transverse lamella continues to expand vertically under the pressure difference between the two phases and consumes the two impacting cylinders, evolving into a thin film aligned with the $xy$-plane; in the meantime it is continuously decelerated by capillary force, forming a thick rim at its upper border.

Figure 3. Isometric snapshots showing the liquid sheet expansion process at $We = 200,\ \varepsilon = 0.06$ and $N_{max} = 25$. From left to right: $t/\tau _{cap} = 0.029$, 0.113 and 0.454.

Figure 4. Snapshots showing the liquid sheet expansion process at $We = 60,\ \varepsilon = 0.06$ (a–c), $We = 200,\ \varepsilon = 0.06$ (d–f) and $We = 60, \varepsilon = 0.02$ (g–i). From left to right: $t/\tau _{cap} = 0.91$, 1.82 and 2.73. For all three cases, $N_{max} = 25$.

Without the initial interface perturbation, the liquid rim would remain intact and mostly smooth during the entire simulation period up to $t/\tau _{cap} = 2.73$, where $\tau _{cap} \equiv \sqrt {\rho _l d_0^3/8\sigma }$ is the capillary time scale corresponding to the initial rim diameter. When perturbation is introduced, we observe the amplification of transverse perturbation waves on the rim as shown in figure 4. Depending on the value of $We$, these waves will either slowly increase in amplitude and reach nonlinear development, whose growth rate increases with the initial perturbation amplitude $\varepsilon _0$, as shown in the first and last row of figure 4, or generate slender ligaments which continue to elongate on the top of ‘cusp’ structures (Gordillo, Lhuissier & Villermaux 2014), as shown in figure 3 and figure 4(d–f). It is noted that, for all three cases presented in figure 4, the number density of transverse ligaments (or the characteristic perturbation wavenumber when no such ligaments form) decreases over time rather than remaining a constant; a phenomenon which we will analyse in more detail in § 5. In the meantime, lamella expansion gradually slows down as the vertical position of its bordering rim reaches saturation, and begins to slightly retract at low $We$ values, as can be seen in figure 4(b,c). As is the paradigm of droplet formation in many previous fragmentation studies (Villermaux 2007), here, the transverse ligaments act as the direct and only source of fine drops as the latter intermittently detach from the ligament tips, whose statistics will be presented in § 6. Under the restoring effects of the capillary force, most of the fragments undergo prolate–oblate shape oscillations after their detachment, and a small portion of them may cross paths and merge to form larger droplets during their flight (Tang et al. 2023). Note that in our simulations all fragments keep moving upwards before crossing over the top boundary and leaving the simulation domain; while in wave splashing scenarios without wind forcing they will ultimately fall back to the sea surface under gravity and be destroyed (Mostert et al. 2022).

4. Liquid lamella expansion

In this section, we study the dynamics of the liquid lamellae consisting of the expanding sheet (§ 4.1) and its bordering rim (§ 4.2). To simplify the problem, which features random initial perturbation, here, we do not consider the formation of liquid ligaments. Instead, we average the transverse cross-section of the coalesced liquid bulk along the $x$ axis following Wang & Bourouiba (2017) so that the lamella expansion process can be described in a quasi-one-dimensional manner, characterised by a one-dimensional velocity and thickness profiles $u_y(\kern0.7pt y,t)$ and $h(\kern0.7pt y,t)$.

4.1. Liquid sheet kinematics

The unsteady evolution of the expanding sheet profile is of both fundamental and practical importance for impact problems, especially for predicting their maximum spread radius (Yarin 2006; Josserand & Thoroddsen 2016; Wang & Bourouiba 2017), and the first step towards modelling the profile evolution of expanding sheets is to understand the fluid motion within them.

We first derive the velocity profile within liquid lamella sheets in Cartesian coordinates for $We \gg 1$ and $Oh \ll 1$. At early times, the vertically expanding sheet is bounded at its lower end by the line of collision between the two cylinders, which we call the lamella foot. The trajectory of this foot is given by $y_n = 2\sqrt {U_0t/R_0}$, per the analysis given in Appendix B and following Gordillo, Riboux & Quintero (2019). At later times, this lamella foot becomes increasingly indistinct as the colliding cylinders merge.

We now proceed to discuss the kinematics of the fluid within the vertical lamella sheet itself, subject to these considerations, following the general analytical strategy of Wang & Bourouiba (2017). After neglecting viscous, compressibility and capillary effects, the quasi-one-dimensional momentum equation in the vertical direction can be written as

(4.1)\begin{equation} \frac{\partial u_y}{\partial t} + u_y \frac{\partial u_y}{\partial y} = 0, \end{equation}

which may be presented alternatively in the Lagrangian form along the characteristics ${{\rm d}y}/{\rm d} t = u_y$

(4.2)\begin{equation} \frac{{\rm D} u_y}{{\rm D} t} = 0, \end{equation}

where ${\rm D}/{\rm D} t$ is the total derivative. Equation (4.2) suggests that the velocity of a fluid particle remains unchanged within the liquid sheet as it travels vertically upwards. We can then integrate along the characteristics and obtain the motion of fluid particles within a Lagrangian frame following it

(4.3)\begin{equation} y = u_y t + \xi, \end{equation}

where $u_y$ and $\xi$ are the initial vertical velocity and position of fluid particles. Following Yarin & Weiss (1995), Wang & Bourouiba (2017) assume that the initial velocity is proportional to the initial position, namely $u_y = \beta \xi$, similar to the velocity field of a stagnation-point flow. Consequently, (4.3) may be rewritten in the Eulerian reference frame as

(4.4)\begin{equation} u_y(\kern0.7pt y,t) = \frac{y}{t + \dfrac{1}{\beta}}, \end{equation}

where $1/\beta$ corresponds to the time needed to set up the post-collisional velocity profile within the liquid bulk, which according to Wang & Bourouiba (2017) is of the same order as the collision time scale $d/U_0$. As sheet expansion further progresses so that $t \gg 1/\beta$, (4.4) asymptotes to

(4.5)\begin{equation} u_y(\kern0.7pt y,t) = \frac{y}{t}. \end{equation}

We measure the vertical velocity profile within the expanding sheet from our numerical simulations at different times for $We = 120$, and first plot them in the inset of figure 5(a). Consistent with our assumption, the vertical velocity $u_y$ is observed to scale linearly with the vertical position $y$, with the slope decreasing over time. In the main plot we rescale the horizontal axis by $U_0 t$, after which the velocity profiles at different times all collapse onto a single straight line $y=x$ for $y \leq 2U_0 t$, agreeing with the theoretical model (4.5). For $We = 120$, the collision time scale is $d/U_0 \approx 0.52 \tau _{cap}$, while figure 5(a) indicates that the velocity profile at $t/\tau _{cap} = 0.73$ is already very close to the analytical solution (4.5), which suggests that, in fact, the initialisation time scale $1/\beta \ll d/U_0$, and model (4.5) is applicable to the early collision stage where $t \approx d/U_0$. Note that, as time elapses, the vertical velocity deviates from (4.5) at progressively smaller values of $y/U_0 t$, where the fluid parcels move away from the sheet towards the bordering rim. Indeed, in figure 5(b) we also plot the vertical velocity normalised by $y/t$ at different times, $We$ values and perturbation waveforms. All results presented therein feature a range where the normalised velocity $u_y t/y = 1$, suggesting that (4.5) remains valid at different initial configurations within our current parameter space.

Figure 5. (a) Liquid sheet velocity profile at $We = 120$, scaled according to (4.5); (b) verification of (4.5) at different values of $We$, time and perturbation waveforms. ‘Sing’. indicates that the initial perturbation we impose features a single wavenumber $N_{max}$, while ‘Sup’. denotes a combination of sinusoidal perturbations with wavelengths $\lambda = D/8,\ D/16$ and $D/32$.

Having established the liquid velocity profile within the lamella sheet, we can further solve for its thickness profile utilising the continuity equation, which can be written as follows for a thin sheet expanding in the $y$ direction:

(4.6)\begin{equation} \frac{\partial h}{\partial t} + \frac{\partial (u_y h)}{\partial y} = 0, \end{equation}

combined with (4.5), this yields

(4.7)\begin{equation} t \frac{\partial h}{\partial t} + y\frac{\partial h}{\partial y} + h = 0. \end{equation}

This can be solved by separation of variables in the form of $h(\kern0.7pt y,t) = f(t)g(\kern0.7pt y)$ to obtain the evolution of sheet thickness $h$, which is written in a non-dimensional formulation as

(4.8)\begin{equation} \frac{hU_0t}{R_0^2} = f \left(\frac{y}{U_0 t} \right), \end{equation}

note that this result differs from the axisymmetric configuration where $h U_0^2 t^2 / R_0^3$ evolves self-similarly without an explicit time dependence (Wang & Bourouiba 2017; Gordillo et al. 2019).

We plot the lamella thickness profiles for $We = 120$ in figure 6(a), where it can be seen that the ‘bulge’ centred around $y=0$ gradually flattens into an extended thin liquid sheet, pushing the bordering rim further along the vertical direction. The bordering rim is connected to the sheet via a neck, reminiscent of capillary waves upstream of inviscid liquid rims receding at the Taylor–Culick velocity (Savva & Bush 2009). Figure 6(b) further shows the profiles rescaled by (4.8). The ‘bulk’ region where $y/R_0 \leq \sqrt {2 U_0 t/R_0}$ initially retains its cylindrical shape during the initialisation period $t \sim 1/\beta$, and only comes to agreement with (4.8) when the lamella foot disappears and the bulk can no longer be decisively told apart from the lamella. The non-dimensionalised lamella profile in figure 6(b) is found to be well described by the following functional form $f(x)$:

(4.9)\begin{equation} f(x) = 0.5081 \,{\rm e}^{{-}x} + 2.782 \,{\rm e}^{{-}2x}. \end{equation}

Figure 6. (a) Liquid sheet profiles at $We = 120$; (b) comparison between interface profiles non-dimensionalised according to (4.8) and the exponential fit (4.9).

Lastly, we note that, although the liquid lamella expands vertically to form a thin film, the latter does not suffer from spontaneous perforation during our simulation period, although this is observed in figure 14 of Vledouts et al. (2016) and marks the onset of bag film fragmentation in droplet aerobreakup problems (Ling & Mahmood 2023; Tang et al. 2023), where the liquid films are subject to radial accelerations. Neither does the lamella film experience destabilisation under shear force arising from its interaction with the surrounding gas phase (Villermaux & Bossa 2011; Riboux & Gordillo 2015), which may arise at larger $We$ values. We therefore do not take special measures to artificially stabilise (Liu & Bothe 2016) or perforate (Chirco et al. 2022) the expanding lamellae in this study.

4.2. Bordering rim evolution

As is discussed in § 4.1, while the expanding lamella sheet abides by the velocity profile (4.5) and the self-similar thickness profile (4.8), these two models break down for the bordering rim, which demands separate scaling laws to describe its kinematics. Similar rim structures are ubiquitous in impact problems and act as the crucial link between the expanding sheet and shedding droplets (Wang et al. 2018), and understanding their motion lays the foundation for further theoretical analysis of ligament merging, which we will perform in § 5.2.

In a quasi-one-dimensional framework, the kinematics of the bordering rim can be characterised by two parameters, namely its average vertical position $y_{rim}$ and diameter $b_{rim}$. We first present their evolution at different $We$ values in figures 7(a) and 7(b), respectively; where time and length are scaled using $U_0/R_0$ and $R_0$. It is seen that both $y_{rim}$ and $b_{rim}$ first increase with time and eventually saturate; and $y_{rim}$ and $b_{rim}$ show different dependences on $We$, as the former increases and the latter decreases with $We$, although these $We$-dependencies have become very subtle by $We = 200$. Further, the evolution of $y_{rim}$ for $We = 200,\ \varepsilon = 0.06$ and $We = 200,\ \varepsilon = 0.04$ in figure 7(a) is virtually the same, which shows that the dynamic behaviour of the rim is largely independent of the initial perturbation amplitude $\varepsilon$.

Figure 7. (a,b) The evolution of the vertical position $y_{rim}$ (a) and the rim thickness $b_{rim}$ (b) over time, compared with solutions of (4.10)–(4.12) at corresponding $We$ values (solid lines). Early-time measurements from two two-dimensional simulations with $We = 80$ and 160 are also included. (c,d) Results in (a,b) rescaled using (4.13a,b).

We now seek to further compare our numerical results obtained in figures 7(a) and 7(b) with available theoretical models. Following Gordillo et al. (2019), the following mass- and momentum-conservation equations can be proposed in the non-dimensionalised form to describe the dynamic behaviour of the advancing lamella rim in the inviscid limit

(4.10)$$\begin{gather} \frac{\rm \pi}{8} \frac{{\rm d} b_{rim}^2}{{\rm d} t} = [u_y(\kern0.7pt y_{rim}, t) - v_{rim}]h(\kern0.7pt y_{rim}, t), \end{gather}$$

(4.11)$$\begin{gather}\frac{{\rm d} y_{rim}}{{\rm d} t} = v_{rim}, \end{gather}$$

(4.12)$$\begin{gather}\frac{\rm \pi}{8} \frac{{\rm d}}{{\rm d} t} (b_{rim}^2 v_{rim}) = u_y(\kern0.7pt y_{rim}, t){[u_y(\kern0.7pt y_{rim}, t) - v_{rim}]} h(\kern0.7pt y_{rim}, t) - \frac{8}{We}, \end{gather}$$

where the cylinder radius $R_0$, initial impact velocity $U_0$ and their quotient $R_0/U_0$ are chosen as the reference scales for length, velocity and time. The numerical solutions of the ordinary differential equation system (4.10)–(4.12) at various $We$ values, with initial conditions as defined in Appendix B, are presented in figures 7(a) and 7(b), respectively, as transparent solid lines. Since most of our measurements from three-dimensional simulations are taken when $U_0 t/R_0 > 1$, we also present results of two-dimensional simulations at $We = 80$ and 160 using solid dots, which extend to much earlier times. An excellent agreement between the predictions of (4.10)–(4.12) and the two-dimensional results is observed. Three-dimensional measurements of $y_{rim}$ and $b_{rim}$ are found to saturate earlier and become smaller than their two-dimensional counterparts and theoretical predictions as time elapses. These observations can be explained by taking into account the formation of ligaments on the lamella rim. This mechanism arises due to highly nonlinear transverse perturbations imposed on the liquid rim, which are not present in two-dimensional simulations. At sufficiently large times, the growth of ligaments causes a substantial mass flux away from the liquid rim, which becomes more prominent at higher $We$ values, as is shown in figure 4. Other factors may also play a role, e.g. the initial overlapping between the two rims in three-dimensional simulations.

In figures 7(c) and 7(d) we present the evolution of $y_{rim}$ and $b_{rim}$ again, where time is now non-dimensionalised using the capillary time scale $\tau _{cap}$. The growth of both $y_{rim}$ and $b_{rim}$ in three-dimensional simulations is consistent with a power law of $\sqrt {t/\tau _{cap}}$ up to $t/\tau _{cap} \approx 1.4$, after which deviation from this power law is observed. It is also found that prefactors of $\sqrt {We}$ (for $y_{rim}$) and ${We}^{-1/4}$ (for $b_{rim}$) can collapse the evolution data reasonably well, especially for $b_{rim}$ as a single master curve is clearly observed in figure 7(d), leading to the following scaling arguments:

(4.13a,b)\begin{equation} \frac{y_{rim}}{R_0} \propto \sqrt{We} \sqrt{\frac{t}{\tau_{cap}}},\quad \frac{b_{rim}}{R_0} \propto {We}^{{-}1/4} \sqrt{\frac{t}{\tau_{cap}}}, \end{equation}

which we will use for further theoretical analysis of the ligament merging phenomenon in § 5.2, while a complete determination for (4.13a,b) valid for very late time remains for future work. The scaling of $We$ in (4.13a,b) agrees with the experimental results of Villermaux & Bossa (2011) and Wang et al. (2018) for drops impacting a small surface; although the dependence on time is different, as their models aim at describing the entire sheet expansion–retraction motion. Nonetheless, similar empirical $\sqrt {t}$ scalings have been proposed by Mongruel et al. (2009), Thoroddsen, Takehara & Etoh (2012) and Visser et al. (2015) for quantifying the radial position of the expanding lamellae in the inertial regime, indicating that this simple form of time dependence can still describe the rim kinematics reasonably well between its formation and the onset of the retraction motion. Note that (4.9) implies the evolution of the liquid momentum carried by the rim $p_{rim} \equiv {\rm \pi}\rho _l D b_{rim}^2 \dot {y}_{rim}$ is independent of the values of $We$, even though its average vertical velocity $\dot {y}_{rim}$ and volume $\varOmega _{rim} \equiv {\rm \pi}b_{rim}^2 D$ do depend on $We$. This contrasts with the axisymmetric results of Wang & Bourouiba (2022) where the rim volume remains independent of $We$ due to their axisymmetric configuration.

5. Transverse liquid ligaments

5.1. Formation and growth

In the scenario considered by Wang & Bourouiba (2021), liquid ligaments grow slowly out of the corrugated bordering rim along its azimuthal direction, which they ascribed to a combination of local geometry, pulling effects of inertial force associated with rim deceleration and the global liquid-phase mass conservation. For this study, and given our perturbation profile, we find the transverse ligaments form very early for $We \geq 120$, nearly at the same time when the lamella is born out of the indentation region between the two cylinders, as presented in the simulation snapshots of figure 8. A closer look at figure 8(b) reveals that the ligaments are produced preferentially from the concave regions along the two perturbed cylinders, suggesting that the ligament formation process is closely associated with the initial rim perturbation waveform. Indeed, a previous investigation by Gordillo, Onuki & Tagawa (2020) suggests that the ligaments originate from the non-uniform initial distribution of normal interface velocity, which is in turn determined by the upstream liquid velocity and the curved initial cavity profile at the moment of impact. Note also that the nascent ligaments do not always project exactly vertically along the vertical $y$-axis; rather, they can display complex twisting and surface oscillation motions, and may grow in an oblique direction. At the same time, the liquid sheet beneath its bordering rim also features a wavy surface not perfectly aligned with the $xy$-plane. These are most likely due to the difference in the initial interface perturbations imposed on the two colliding rims, which gives rise to oscillatory motions under capillary effects.

Figure 8. Snapshots taken from a simulation case at $We = 200,\ \varepsilon = 0.06$ and $N_{max} = 25$ showing ligaments generated from the ‘indentation’ region between two colliding rims. From left to right: $t/\tau _{cap} = 0.045,\ 0.091$ and 0.136.

We now move on to discuss the subsequent growth of the height of these ligaments after their formation, which is shown in figure 9(a) for rim collision at $We = 120$ and 160, where we present the height evolution of three individual ligaments at each $We$. The shedding of the first fine drop from these ligaments is characterised by a kink around $t = 0.4 \tau _{cap}$. While this initial pinch-off may happen at even earlier times as $We$ increases, as shown in figure 7 of Wang & Bourouiba (2018); this is still much later than the onset of the micro-splashing phenomena investigated by Thoroddsen et al. (2012), which happens at $t/\tau _{cap} = O(10^{-4})$ (see e.g. their figure 5b). Before this first pinch-off event, the height $h_{lig}$ of different ligaments increases following a similar trend, and the growth rate at $We = 160$ is higher than that at $We = 120$. Given that the initial phase of growth of these ligaments is governed primarily by inertio-capillary effects, we compare it in figure 9(a) with the self-similar power law of $h \propto t^{2/3}$, proposed by Lai, Eggers & Deike (2018) in their investigation of inertio-capillary-dominated collapse of small surface bubbles. It is found that the height increase of the majority of transverse ligaments (except ligament 1 at $We = 120$) agrees better with the linear growth model. This is most likely because the formation of fast jets observed by Lai et al. (2018) is preceded by focusing of interfacial capillary waves at the bottom of the bubble cavity while the liquid bulk remains quiescent whereas, here, the bulk velocity plays a vital role in driving the closure of rim indentation, and may thus modify the rate of jet growth. In addition, a recent study by Gordillo & Blanco-Rodríguez (2023) suggests that the exponent of power laws dictating the evolution of jet radius and speed is dependent on the initial geometry of the collapsing air cavity, which may also account for the difference between our results and those of Lai et al. (2018) since the concave regions on our perturbed cylindrical surfaces do not feature a uniform radius of curvature.

Figure 9. (a) The evolution of liquid ligament length measured at different $We$ values, compared with the $t^{2/3}$ scaling law of Lai et al. (2018) and a linear growth model. (b) Vertical component of liquid velocity $u_y$ measured within liquid sheets and ligaments, showing the ballistic region within the ligament proposed by Gekle & Gordillo (2010). ‘Sup’ denotes that the initial rim perturbation is a superposition of sinusoidal signals with wavelengths $\lambda = D/8,\ D/16$ and $D/32$.

To better understand the fluid motion within the growing ligaments, we measure the liquid-phase vertical velocity $u_y$ from the bottom of the expanding sheet up to the tip of a single ligament, and plot it in figure 9(b) as a function of the vertical coordinate $y$. Two linear scaling regimes are observed; the first one is well described by $u_y = y/t$, which corresponds back to the velocity profile (4.5) we established for the expanding sheet. As the fluid particles move higher up and away from the sheet, its velocity first increases and then abruptly decreases as it enters the neck and rim region, respectively; a similar abrupt deceleration is also noted by Wang & Bourouiba (2022) for drop collision problems. Interestingly, the combination of the neck and rim causes a constant decrease in the vertical velocity of approximately $0.7U_0$ for different $We$ values, times and perturbation waveforms. This pattern is not explicable from the scaling model (4.13a,b) since, according to it, the fluid velocity should be equal to the average rim velocity $\dot {y}_{rim}$ after deceleration, which is always one half of the sheet velocity $y_{rim}/t$ before deceleration. The implication is that the fluid velocity at the ligament root is always faster than the average rim velocity, which Wang & Bourouiba (2021) ascribed to the additional acceleration due to the interface curvature at the rim–ligament junction. After this constant offset at the bordering rim which is most likely a viscous effect (Ghabache, Séon & Antkowiak 2014), $u_y$ grows linearly again with the vertical position $y$ with the same slope as the sheet region. This second linear scaling regime most likely corresponds to the ballistic region (although in the present study gravity is not included) in Worthington jets identified by Gekle & Gordillo (2010), where the liquid particles travel at constant speed upwards before being slowed down once again at the bulb by capillary effects.

While Wang & Bourouiba (2022) were able to demonstrate theoretically that energy dissipation at the bordering rim is responsible for the local rapid deceleration of fluid particles, which we also observed in figure 9(b), the detailed nonlinear dissipation mechanism is not captured by their one-dimensional model. As discussed in § 1, there is still a dissipation deficit of $15\,\%$ not covered by their calculations occurring at early time. To help elucidate this discrepancy, we present contour plots in figure 10 for a simulation case at $We = 200$, showing the distribution of the liquid-phase viscous dissipation rate $\epsilon _d \equiv \mu _l (\partial u_i / \partial x_j) (\partial u_j / \partial x_i)$ within the centre plane $x=0$ for $t/\tau _{cap} \leq 0.36$, covering the early deformation period $t/\tau _{cap} \leq 0.2$ where Wang & Bourouiba (2022) observed the dissipation deficit (see their figure 19). It is found that, when the lamella is born at very early time and its bordering rim has not yet fully developed (figure 10a), there is an extremely high concentration of $\epsilon _d$ located at the lamella foot, agreeing with the simulation results of Wildeman et al. (2016) for drops impacting a smooth surface (see their figure 4) and Fudge et al. (2023) for drops impacting a liquid pool (see their figure 8b). As time elapses, the bordering rim takes shape, with its neck featuring relatively high concentration of $\epsilon _d$, whereas the dissipation at the lamella foot weakens and eventually becomes negligible by $t/\tau _{cap} = 0.364$ (figure 10d), matching the saturation trend shown in figure 19 of Wang & Bourouiba (2022) for the dissipation deficit. This decay pattern of liquid-phase dissipation may be explained as follows. The liquid particles feeding the lamella at early times mostly come from a very narrow boundary straddling the corners on either side of the lamella foot (Riboux & Gordillo 2014), and they generate vorticity (Batchelor 2000; Li et al. 2018) and experience capillary deceleration while traversing the highly curved free surface, leading to very large values of viscous dissipation. Since this early-time deceleration is geometry induced, Fudge et al. (2023) further hypothesised that the magnitude of the velocity gradient remains largely unchanged at different flow configurations so that the early-time dissipation rate is proportional to the liquid viscosity $\mu _l$, although this is not directly verified in the present work. As the liquid sheet extends further upwards, the interface curvature at the lamella foot decreases and capillary deceleration becomes much weaker, hence the decrease in the dissipation rate $\epsilon _d$. Eventually, as the lamella foot is no longer discernible from the flattened liquid bulk and the capillary effects become negligible, the inviscid sheet velocity profile (4.5) will be established.

Figure 10. Contour plots visualising the two-dimensional distribution of instantaneous liquid-phase dissipation rate $\epsilon _d$ within the centre plane $x/D = 0.5$ for $t/\tau _{cap}=0.045$ (a), 0.091 (b), 0.182 (c) and 0.364 (d), where $We = 200$ and $\varepsilon _0 = 0.06$.

These observations of the liquid-phase dissipation evolution therefore suggest that the unidentified dissipation deficit found by Wang & Bourouiba (2022) is most likely linked with the early-time lamella formation, which also features strong viscous effects and, according to Wildeman et al. (2016) and Ó Náraigh & Mairal (2023), this kind of dissipation is a type of ‘general head loss’ imposed by the deformation mode of the impacting object and independent of the detailed impact parameters. The head loss in impact problems dissipates a fixed fraction of the initial kinetic energy via recirculating flows (Wildeman et al. 2016; Ó Náraigh & Mairal 2023), which matches the observations of Wang & Bourouiba (2022). It is also noted that, besides the early-time lamella foot and the late-time rim neck, the flow field elsewhere within the coalesced liquid bulk shown in figure 10 is nearly inviscid, supporting our derivation of the centreline velocity profile within the liquid sheet (4.5) based on inviscid flow assumptions. While we have not fully explored the influence of $Oh$ on the lamella expansion process, it can be expected that, in the inviscid limit where $Oh \ll 1$, its influence will be confined to the narrow ‘viscous’ regions identified in this section, and thus will not significantly affect the deformation of the droplet bulk, or the ligament dynamics and fragmentation mechanisms to be discussed below.

5.2. Ligament merging phenomenon

As is noted in § 3, the number of transverse ligaments on the rim will decrease over time as they merge with their neighbours. This merging process is observed only when the initial perturbation is not monochromatic such that the nascent ligaments are not equidistantly spaced, and it turns out to be crucial in maintaining the growth of ligaments and the continuation of fragment shedding. Figure 11(a–c) shows the development of ligaments formed out of monochromatic perturbation waveforms, and it is observed that, after two rounds of pinching-off events at their tips, the ligaments can no longer sustain their own growth and are re-absorbed back into the underlying liquid rim whereas ligaments formed out of filtered white noise perturbations at the same values of $We, \varepsilon$ and $N_{max}$ merge with their neighbours and survive end-pinching, continuing to grow and shed fragments until the end of simulation, as observed in figure 11(d–f).

Figure 11. Snapshots at $We = 120,\ \varepsilon = 0.06$ and $N_{max} = 25$ showing ligament evolution from monochromatic initial perturbations (a–c) and filtered white-noise perturbations (d–f). Re-absorption of ligaments back into the rim is observed for the monochromatic perturbation case after two cycles of drop shedding. From left to right: $t/\tau _{cap} = 0.2, 0.4$ and 0.6.

The ligament merging process is shown in more detail in the snapshots of figure 12, where the merging ligaments are observed to be located on rim ‘cusp’ structures extruding from the liquid sheet beneath, an indication of the non-uniform incoming mass distribution along the bordering rim (Gordillo et al. 2014; Wang & Bourouiba 2018). When ligament merging occurs, the roots of two neighbouring ligaments approach each other; liquid is then drawn upwards from the underlying cusp in between the two ligaments, causing them to coalesce into a thicker and more corrugated ligament; while the length of the fused ligament does not differ significantly from those of its parents. Ligament merging is thus capable of delaying the next end-pinching event since the fused ligament takes longer to be stretched, thus allowing incoming fluid from the rim to sustain their growth and evade re-absorption into the rim when the depleted ligament becomes too short, as observed in figure 11. Note that the monotonically decreasing trend of ligament numbers observed by us and Wang & Bourouiba (2021) differs from the early experimental work of Thoroddsen & Sakakibara (1998) on drop impact, where they also observed splitting of liquid fingers aside from their merging behaviour, so that the total finger number remains approximately constant.

Figure 12. Snapshots taken from a simulation case at $We = 160,\ \varepsilon = 0.04$ and $N_{max} = 25$ showing ligaments merging on the corrugated rim bordering the expanding sheet, while shedding fragments via the end-pinching mechanism. From left to right: $t/\tau_{\rm cap} = 0.73, \, 1.09, \, 1.45$.

To the knowledge of the authors, a scaling law for the ligament numbers accounting for their merging dynamics is not yet available. The early work of Marmanis & Thoroddsen (1996) found that the number of liquid fingers at the maximum spread radius for high-speed drop impacts scales with $Re^{3/4}$ ($Re$ being the impact Reynolds number), without accounting for their evolution over time, whereas the recent $We^{3/8}$ scaling model proposed by Wang & Bourouiba (2021) is based on the understanding that the ligaments are formed from a subset of rim corrugations arising from a combined RT and RP instability; a physical mechanism which is most likely not yet active within our parameter space as we find the formation of ligaments more closely linked with the initial interface perturbation geometry. The early theoretical analysis of Yarin & Weiss (1995) predicted that perturbed free rims will spontaneously develop ‘cusp’ structures due to nonlinearity, where two neighbouring rim sections impinge and give rise to free jets. This physical picture closely resembles our present observations, although Yarin & Weiss (1995) did not proceed to develop scaling models for the splashing fragments. We therefore seek to derive a new ligament number scaling model accounting for the merging dynamics in this section, and compare it with the simulation results.

As the first step towards quantifying the ligament dynamics, we seek to determine the evolution of the average ligament width $w_{lig}$ over time by inspecting the evolution of fragment diameter $d_{frag}$, which can be easily reconstructed using the droplet-tracking algorithm of Chan et al. (2021). The connection between these two quantities is established in figure 13(a), where we plot the ratio between $d_{frag}$ and $w_{lig}$ at the instant of pinch-off. The ligament diameters are measured at the cross-section corresponding to one half of the total ligament length. Most of the measured data are found to scatter between 1 and 2, centred around 1.4, which is close to the average value of 1.5 as found by Wang & Bourouiba (2018). These are below the theoretical value of 1.89 as predicted by the RP instability (Gordillo & Gekle 2010), indicating that end-pinching is indeed the dominant fragment production mechanism for the present study, and that the diameter of fragments $d_{frag}$ remains in proportion to the width of their parent ligaments $w_{lig}$.

Figure 13. (a) Measurement of the ratio between the diameter $d_i$ of the detaching fragment and the width $w_i$ of its originating ligament at different ejection times. (b) The evolution of the fragment diameter $d_{frag}$ of ejected fragments. The results in (b) have been ensemble-averaged across three realisations for each $(We, \varepsilon )$ configuration, and rescaled by $We^{-0.75}$ in the main plot.

The inset of figure 13(b) shows the diameters of the fragments $d_{frag}$ vs their time of formation, where the data for each configuration of $(We, \varepsilon )$ have been averaged across three individual realisations to reduce the range of scatter. It is found that larger fragments are generally produced at later times and smaller $We$ values. The scatter in data increases over time, which is most likely due to the increase in the diameter difference and the surface corrugation of remaining ligaments as they merge with one other. The main plot suggests that the evolution of individual fragment diameter in figure 13(b) roughly scales with $We^{-3/4} \sqrt {t/\tau _{cap}}$, while noting that our fitted prefactor ${We}^{-3/4}$ is not conclusive and remains to be further validated at higher $We$ values. Since the fragment size remains in proportion to the width of the parent ligament according to figure 13(a), it can be inferred that

(5.1)\begin{equation} w_{lig} \propto d_{frag} \propto \sqrt{\frac{t}{\tau_{cap}}}. \end{equation}

Gordillo et al. (2014) and Wang & Bourouiba (2018) proposed the following drift velocity of ligaments $u_{drift}$ on top of a liquid cusp to characterise their migration:

(5.2)\begin{equation} u_{drift} = [u_s(\kern0.7pt y_{rim}) - u_{rim}] \sin \alpha, \end{equation}

where $u_s(\kern0.7pt y_{rim}) = y_{rim}/t$ is the velocity at which liquids enters the rim from the expanding liquid sheet, calculated from the self-similar expanding sheet velocity profile (4.5), and $u_{rim}$ is the vertical rim velocity determined by differentiating the scaling law for $y_{rim}$ in (4.13a,b). It is noted that both $u_s(\kern0.7pt y_{rim})$ and $u_{rim}$ scale as $\sqrt {We}\ t^{-1/2}$, and $\alpha$ is the angle between the local rim and the horizontal plane, as shown in figure 14. Overall, (5.2) suggests that the migration of ligaments is driven by the tangential projection of the net incoming fluid velocity along the corrugated rim.

Figure 14. Main plot: sketch showing the quantities defined in § 5.2 for developing the ligament merging model (5.6). Inset: sketch showing the local geometry of the junction region at the ligament base.

As $u_{drift}$ drives the ligament migration and causes the ligament number density $N_{lig}$ to decrease, the average transverse ligament spacing $1/N_{lig}$ on the rim consequently increases. Therefore

(5.3)\begin{equation} u_{drift} \propto \frac{{\rm d}}{{\rm d} t} \left(\frac{1}{N_{lig}} \right). \end{equation}

The averaged value of $\sin \alpha$ can be determined by inspecting the geometry of the junction region between the ejected ligaments and the lamella rim, which is shown in the inset of figure 14. Making use of the law of cosines

(5.4)\begin{equation} \frac{w_{lig}}{2b_{rim}} = \cos \left(\frac{\rm \pi}{2} - \alpha\right) = \sin \alpha. \end{equation}

This, combined with (5.1) and (4.13a,b), yields

(5.5)\begin{equation} \sin \alpha = C(We), \end{equation}

suggesting that the rim slope $\alpha$ remains unchanged over time and depends only upon the impact Weber number. Taking into account that $\tan \alpha \approx \varepsilon _{rim} N_{lig}$, $\alpha$ remaining constant also indicates that as the average ligament spacing $1/N_{lig}$ becomes larger over time owing to the ongoing merging of ligaments, the rim corrugation $\varepsilon _{rim}$ also increases proportionally to maintain a constant local slope.

By further incorporating (5.3) and (5.2), the following model predicting the evolution of the ligament number density $N_{lig}$ can be derived:

(5.6)\begin{equation} N_{lig} \propto {\left(\frac{t}{\tau_{cap}} \right)}^{{-}1/2}. \end{equation}

In the main plots of figure 15 we show the decay of the ligament number density $N_{lig}$ at different values of $N_{max}$ and $We$. Figure 15(a) shows that while increasing $N_{max}$ leads to the formation of more ligaments at early time; $N_{lig}$ appears to reach saturation and does not increase proportionally with $N_{max}$ when $t/\tau _{cap} = 0.18$, and it is likely that viscous damping effects are particularly strong when the lamella is ejected at the early impact stage, which may smooth out short-wavelength components in the initial perturbation spectra (see our figure 10 and relevant discussions). Larger $N_{max}$ also leads to faster decay of $N_{lig}$, as expected, since the average distance between neighbouring ligaments becomes smaller, and for all values of $N_{max}$, figure 15(a) suggests that the evolution of $N_{lig}$ becomes largely similar for $t/\tau _{cap} > 1.5$, where the remaining ligaments take much longer to migrate and merge. Figure 15(b) shows that the decay of $N_{lig}$ does not appear to depend strongly on $We$, and thus lends some support to the prediction of (5.6) that it is $We$-independent, although ensemble averaging would be needed to ascertain this due to the randomness in the initial perturbation waveform. The insets of figure 15 show the evolution of $N_{lig}$ again in log–log axes, which is also compared with model (5.6). The inset of figure 15(a) indicates that the measured results collapse well and agree with (5.6) for $N_{max} \geq 35$ throughout the period of measurement, whereas at $N_{max} = 25$ the decay of $N_{lig}$ is initially slower, and only matches the prediction of $t^{-1/2}$ for $t/\tau _{cap} \geq 1$. The measurement of $N_{lig}$ extends to earlier times in figure 15(b), and its inset suggests that the evolution of $N_{lig}$ at different $We$ is initially slower and matches model (5.6) only after $t/\tau _{cap} \geq 0.2$. This slower decay at early times observed in both insets arises most likely because the rim slope $\alpha$ takes a finite period of time to develop before reaching the steady-state value given by (5.5). Overall, these results indicate that model (5.6) offers a good working description of the ligament merging phenomenon occurring within our parameter space. However, the approximations we make for its derivation restricts its application to the scenario when $N_{lig}$ is large and the slope $\alpha$ of the corrugated rim has reached the steady-state value predicted by (5.5).

Figure 15. Evolution of the ligament number density $N_{lig}$ at different values of $N_{max}$ with $We = 120$ (a) and different $We$ with $N_{max} = 60$ (b). The insets compare the evolution of $N_{lig}$ with our model (5.6).

6. Droplet generation and characteristics

Since the expanding sheet remains intact during our simulation period, fragments are formed solely through the breakup of liquid ligaments originating from the bordering rim. The majority of fragments are produced via the end-pinching mechanism, a process that coincides with the ligament merging phenomena and is already visible in the snapshots presented in figures 4 and 12. Namely, the ligament tips are decelerated by capillary force and produce surface corrugations, which in turn create a local pressure gradient within the ligament neck that drains the liquid towards the tip and triggers pinch-off. After this, the enlarged ligament tip detaches as a primary drop (Gordillo & Gekle 2010), whose diameter is proportional to the parent ligament width as shown in figure 13(a). Occasionally, satellite drops are formed from the small amount of remnant liquid within the neck before it can be fully reabsorbed into the ligament after end-pinching, as also shown in figure 5(b) of Wang & Bourouiba (2018). In contrast to the primary drops, the production of these satellite drops is sensitive to initial liquid-phase velocity perturbations, thus introducing randomness to the jet breakup process. However, they do not affect the size and velocity of the subsequent primary drops ejected, as shown by Berny et al. (2022) in the instance of jet-droplet production by bursting of surface bubbles. Possibly due to the difference in their generation mechanisms, our ligaments can grow much longer than those Wang & Bourouiba (2018) observed in their experiments, and at late times, some particularly long ligaments may break up due to the RP instability and shed multiple primary fragments at a time, as can be seen in the rightmost ligament in figures 4(e) and 2(f) of Liu et al. (2022). This appears inconsistent with the conclusion of Wang & Bourouiba (2021) that a ligament can only pinch off to produce one primary drop at a time under the chaotic dripping regime.

We first analyse the time evolution of various fragment properties during the rim collision process. Figure 16(a) shows the evolution of the total number density of primary fragments $N_{frag}$, where it is observed to generally increase over time towards saturation, despite infrequent decreases due to coalescence of primary fragments with another fragment or a neighbouring ligament. It is also noted that $N_{frag}$ increases with both $We$ and $\varepsilon$, as larger $We$ and $\varepsilon$ encourages the growth and subsequent pinch-off of liquid ligaments.

Figure 16. The evolution of the total number density $N_{frag}$ (a) and the ejection velocity (b) of primary fragments, compared with the rim velocity $u_{rim}$ derived from solving (4.10)–(4.12) (solid transparent lines) and (4.13a,b) (dash-dotted line).

The decrease of drop production rate $\Delta N_{frag}/\Delta t$ can be explained as follows. The ongoing ligament merging phenomenon causes the ligament number density to decrease, hence fewer fragments can detach at the same time. Meanwhile, merged ligaments become more corrugated and thicker, thus end-pinching events occur at larger ligament widths as time elapses. Consequently, the time interval between successive end-pinching events also becomes longer, since the necking time scale

(6.1)\begin{equation} t_{neck} \equiv 1.13 \sqrt{\frac{\rho_l w_{lig}^3}{\sigma}} \end{equation}

increases with $w_i$ according to the experimental results of Wang & Bourouiba (2018).

Figure 16(b) shows the evolution of the fragment ejection speed $u_{ej}$, the speed of the drop at the moment it detaches from its parent ligament. The fastest drops found in our simulations feature $u_{ej}$ slightly over $3U_0$, which is comparable to the fragment ejection velocity in the prompt splashing phenomena (Burzynski, Roisman & Bansmer 2020). Regardless of the detailed geometrical features of the parent ligaments, when scaled with the initial collision speed $U_0$, the decay of $u_{ej}$ over time does not significantly depend on either $We$ or $\varepsilon _0$. We further compare the measured $u_{ej}$ values with the predictions of both the rim dynamic equations (4.10)–(4.12) and the scaling model (4.13a,b). Both capture the early-time evolution of $u_{ej}$ up to $t/\tau _{cap} \approx 0.5$, after which the decrease of $u_{ej}$ slightly steepens and shows a better agreement with the predictions of (4.10)–(4.12). This is most likely due to the late-time capillary deceleration being not well represented by (4.13a,b). We note that, while the rim velocity predicted by (4.10)–(4.12) is closer to $u_{ej}$, they tend to over-predict the three-dimensional simulation results, as already observed in figure 7, suggesting the existence of a gap between the actual rim velocity and the fragment ejection velocity. This agrees with the earlier results of Liu et al. (2022), where their figure 5 also shows fragment speeds remaining higher than the rim velocity. Similar measurements of the fragment ejection velocity have also been reported by Thoroddsen et al. (2012) for micro-splashing in drop impact problems at much larger values of $We$, which were well explained by the theory of Riboux & Gordillo (2015). However, in Riboux & Gordillo (2015) the lamella rim fragments under RT and capillary instabilities, and therefore the velocity and size of ejected droplets are directly determined by the rim. In Thoroddsen et al. (2012) fragmentation also occurs shortly after the emergence of the lamella sheet, based on which Riboux & Gordillo (2015) modelled the droplet ejection velocity using flow field information at the lamella foot. These differ from our scenario where the ejection of droplets are governed by the end-pinching of ligaments erupting on the lamella rim (Wang & Bourouiba 2018). Development of theoretical models capable of predicting fragment ejection speed in our case thus requires detailed knowledge of ligament growth and the merging dynamics, which is out of the scope of the current work.

Since rim collision is a transient process where both the total number and the individual size of fragments increase over time, as shown in figure 13, it is of interest to determine how the fragment size and velocity distribution functions evolve with time. We first show the fragment size distributions $n(r/R_0)$ in figure 17(a) at different times $t_c/\tau _{cap}$, which are computed by sampling over a time window of $t_c - 0.45\tau _{cap} \leq t \leq t_c + 0.45\tau _{cap}$. To ensure statistical convergence of the data, three individual realisations are computed for each initial configuration $(We, \varepsilon, N_{max})$ and all results presented in figures 17 and 18 have been averaged across these ensembles. It can be seen that, initially, fragments with $r \leq 0.2R_0$ are produced at early time. While the number density $n$ within this range remains largely unchanged over time and does not appear to depend strongly on $r$, this should be treated with caution due to a lack of fully established grid independence of fragment statistics at small sizes, as discussed in Appendix A. As time elapses, the number of larger fragments increases and causes the falling tail of the distribution to move further to the right, indicating that larger drops are fewer and produced later in time. This can be attributed to the ongoing ligament merging process, as it increases the thickness of individual ligaments.

Figure 17. (a) The evolution of time- and ensemble-averaged size distribution function $n(r/R_0)$ of all fragments produced by colliding rims at $We = 200$, $\varepsilon = 0.06$ and $N_{max} = 25$. (b) The fragment size probability distribution function $f(r/R_0)$ compared with the experimental data and model of Néel et al. (2020). (c,d) The influence of $We$ (c) and $\varepsilon$ and $N_{max}$ (d) on the fragment size distribution function, where $\varepsilon = 0.06$ and $N_{max} = 25$ for all simulation results presented in (c), and $We = 200$ for those presented in (d).

Figure 18. (a,b) Ensemble-averaged vertical (a) and in-plane (b) components of fragment ejection velocity $u_y$ and $u_{xz}$ calculated at various $We$ values, with $N_{max} = 25$. (c) Evolution of the fragment velocity distribution over time, obtained at $We = 200, \varepsilon = 0.06$ and $N_{max} = 25$. (d) The probability distribution functions of fragment velocities at different $We$ values.

We now seek to compare our numerical results with the experimental data and the theoretical model of Néel et al. (2020). It is noted therein that the variation in the fragment size distribution for the rim collision problem arises from two sources, namely those of the transverse ligament size and of the fragments produced from the breakup of a single ligament. A linear superposition of these two effects yields the following size distribution function:

(6.2)\begin{equation} p \left(\zeta = \frac{r}{\bar{r}} \right) = \frac{2{(mn)^{{(m+n)}/{2}}}}{\varGamma(m) \varGamma(n)} \zeta^{{(m+n)}/{2} - 1} K_{m-n} (2\sqrt{mn\zeta}), \end{equation}

where $K_{m-n}$ is a modified $(m-n)$th-order Bessel function of the second kind and $m$ and $n$ reflect the roughness of the distribution of ligament widths and corrugation amplitudes of individual ligaments. Néel et al. (2020) fit their experimental data at $We = 193$ using (6.2) with $m=40$, $n=5$, which we can reproduce in figure 17(b) together with our fragment size probability density function at $We = 200$. Note that the fragment sizes were originally normalised by Néel et al. (2020) using the average fragment diameter $\bar {d}$ in their figure 15(b), which is shown in their figure 13(b) to saturate at large $We$ values and remain proportional to their interstitial sheet thickness $h$. The sheet thickness is in turn related to the pre-collisional rim radius $R_0$ via the rim collision Weber number, given in their work as $We = 16 R_0/h$. This allows us to estimate their average fragment diameter as

(6.3)\begin{equation} \bar{d} = \chi h = \frac{16 \chi}{We} R_0. \end{equation}

We find that setting the coefficient $\chi$ to 5 leads to an excellent match between our simulation results with $(We, \varepsilon, N_{max}) = (200, 0.06, 25)$ and the re-normalised data of Néel et al. (2020) for $r \geq 4\varDelta _{10}$. Figure 13(b) of Néel et al. (2020) suggests a $\chi$ value of approximately 25 for their controlled rim production setup, which is larger than our fitted value by a factor of 5. This might be because rim fragmentation has completed in the experiments, and the larger fragments produced at later times increase the value of $\bar {d}$. Our numerical results differ from Néel et al. (2020) for $r \leq 4\varDelta _{10}$, where their model (6.2) exhibits a fall-off not found in our data. This more uniform portion of fragment size distribution for small $r$ values is also found in figure 13(b) of Lhuissier & Villermaux (2012), where they attributed it to the transverse impact between adjacent rim ligaments. The differences between their experimental and our numerical configurations may also play an important role, as their expanding rims feature toroidal shapes and therefore coalesce within a finite period of time. Last but not least, the two liquid rims are connected by an interstitial thin film in the configuration of Néel et al. (2020), whereas in our case the rims come into contact directly.

The dependence of the time-averaged fragment size distribution on the controlling parameters $We, \varepsilon$ and $N_{max}$ is further shown in figures 17(c) and 17(d). It can be seen that the number density $n$ of small fragments with $r \leq 4\varDelta _{10}$ continues to increase with $We$, $\varepsilon$ and $N_{max}$, consistent with our observations in figure 16(a). More specifically, here, the sensitivity of $n$ to $\varepsilon$ and $N_{max}$, as shown in figure 17(d), further supports our understanding that the differences between the numerical and experimental initial conditions causes the difference between the corresponding results. The number density of intermediate fragments with $4\varDelta _{10} \leq r \leq 0.2 R_0$ also increases with $We$ in figure 17(c), however, different from that of smaller fragments, it appears to asymptote to ${(r/R_0)}^{-1/2}$ for $We \geq 240$ or $N_{max} \geq 60$. The tail of the distribution functions consisting of even larger fragments with $r \geq 0.2 R_0$ remains approximately independent of the controlling parameters, and its decaying trend agrees well with a power law of ${(r/R_0)}^{-7/2}$. While this finding differs from that of Néel et al. (2020), where $m$ and $n$ decrease with increasing $We$ and cause the slope of the tail to decrease correspondingly, this may again be due to differences in detail of the initial conditions. Overall, figures 17(c) and 17(d) suggest that, as $We$ increases, the fragment size distribution within $r \geq 4\varDelta$ asymptotes to a regime independent of the controlling parameter, and well described by a power-law decay with a break in slope, which may originate from the insensitivity of the ligament merging and breakup phenomena to the initial perturbation configurations. It is noted that a similar transition between two power-law regimes has been observed in the droplet size distributions associated with wave breaking by Mostert et al. (2022) and Erinin et al. (2023a), although their distributions feature steeper slopes, with the power law transitioning from $r^{-2}$ to $r^{-6}$.

Here, we show that the power-law scaling ${(r/R_0)}^{-1/2}$ we observed in figures 17(c) and 17(d) at small fragment sizes can be derived from the ligament merging dynamics previously established in this work. The ligament merging time scale is defined as the ratio between the average ligament spacing $L_0/N_{lig}$ and the drift velocity $u_{drift}$, which can be evaluated based on (5.2)

(6.4)\begin{equation} \Delta t_{merge} \equiv \frac{L_0}{N_{lig} u_{drift}} \propto N_{lig}^{{-}2}. \end{equation}

The rate of droplet shedding is controlled by the ligament necking process, thus $\Delta t_{shed} = t_{neck}$ as given in (6.1) (Wang & Bourouiba 2018). Consequently, the total number of fragments shed from all ligaments between two consecutive merging events can be estimated as

(6.5)\begin{equation} N_{frag} = N_{lig} \frac{\Delta t_{merge}}{\Delta t_{shed}} \propto N_{lig}^{{-}1} w^{{-}3/2}, \end{equation}

where $w$ is the average width of ligaments.

As we observed in figure 13(a,b), $w \propto r \propto \sqrt {t/\tau _{cap}}$, whereas (5.6) suggests $N_{lig} \propto {(t/\tau _{cap})}^{-1/2}$. Substituting these two scalings into (6.5) leads to

(6.6)\begin{equation} N_{frag} \propto {(t/\tau_{cap})}^{{-}1/4} \propto {(r/R_0)}^{{-}1/2}. \end{equation}

These derivations suggest that the ${(r/R_0)}^{-1/2}$ scaling at small fragment sizes can be explained as a competition between ligament merging and end-pinching, whereas the ${(r/R_0)}^{-7/2}$ scaling found at larger fragment size ranges still awaits further analysis. It is likely that this steeper dependence on $r$ originates from the presence of corrugations on individual ligaments, or the variation of ligament width across the bordering rim (Néel et al. 2020), neither of which has been considered in the derivations above.

Lastly, we discuss the fragment velocity statistics. Figures 18(a) and 18(b) present the ensemble-averaged vertical (figure 18a) and in-plane (figure 18b) components $u_y$ and $u_{xz}$ of the fragment ejection velocity, which are plotted as functions of the fragment radius $r$. Since most of the initial liquid momentum is deflected in the vertical direction over the collision process, $u_y$ remains a few times to a decade larger than $u_{xz}$. Figure 18(a) shows that the $u_y$ values for fragments within the size range of $0.1 \leq r/R_0 \leq 0.4$ collapse reasonably well when scaled by the initial rim velocity $U_0$, and are well predicted by a power-law decay model $u_y \propto {(r/R_0)}^{-1}$. The velocity distribution at larger fragment sizes deviate from this power-law scaling, most likely due to a combination of late-time effects including amplified ligament corrugations, rim deceleration and the onset of RP instability on the ligaments. The in-plane velocity component values $u_{xz}$ measured in figure 18(b) are more scattered compared with figure 18(a), but the distribution can still be roughly described by the same power-law decay model $u_{xz} \propto {(r/R_0)}^{-1}$. These observations further corroborate the fragment size distribution model (6.6) we proposed, since the power-law exponent of $-1$ for $u_y$ can be derived based on our observations in figures 16(b) and (5.1)

(6.7)\begin{equation} u_y \propto {(t/\tau_{cap})}^{{-}1/2} = {[(t/\tau_{cap})^{1/2}]}^{{-}1} \propto {(r/R_0)}^{{-}1}. \end{equation}

As the only source of the liquid in-plane motion is the transverse drifting of ligaments, the in-plane velocity component $u_{xz}$ can be estimated by the drifting velocity $u_{drift}$. Combining (5.2) and (5.6) similarly leads to

(6.8)\begin{equation} u_{xz} \propto u_{drift} \propto N_{lig} \propto {(t/\tau_{cap})}^{{-}1/2} \propto {(r/R_0)}^{{-}1}. \end{equation}

The good agreement between the velocity scaling models derived above and our numerical results once again highlights the importance of the rim ligament merging phenomenon in determining the size and velocity statistics of splashing fragments.

Figures 18(c) and 18(d) show the number density of fragments as a function of their travel speed $u_{frag}$. Figure 18(c) suggests that, while most of the fragments produced at early time feature a skewed velocity distribution peaking at $u_{frag} \approx 2.8U_0$, the maximum fragment speed decreases and more fragments travelling at lower speeds are recorded as time elapses, and the distribution function has developed a plateau by $t/\tau _{cap} = 1.36$. This suggests that, as the ligaments continue to grow in length and break up, the droplets produced come to span uniformly across a large range of travelling speeds, with fewer drops featuring very slow or particularly fast speeds. Figure 18(d) shows the velocity distribution at different $We$ values averaged over the entire collision process. It is found that, as $We$ increases, the velocity distribution becomes broader, and for $We$ beyond 240, the right tail of the velocity distribution appears to reach a $We$-independent regime, similar to our observation in figure 17(c) for the size distribution of fragments, and the fastest speeds recorded is around $u_{frag} \approx 3.2 U_0$. While a direct comparison with breaking wave statistics (Mostert et al. 2022; Erinin et al. 2023a) is out of the scope of the current work, it is noted that the shapes of velocity distribution functions obtained here in figures 18(c) and 18(d) differ from their counterparts in wave-breaking phenomena (see, e.g. figure 16(d) of Mostert et al. (2022) and figure 12 of Erinin et al. 2023a), as the latter are skewed and narrower than our distributions. This may be due to the presence of gravity in wave breaking, which may arrest the ligament fragmentation process and define a short time scale for completing the splashing phenomenon, thereby reducing the total number of fragments. These splash fragments are themselves also decelerated by gravity. In the velocity distribution, this would appear as a narrowing of the distribution, with lower velocities at the peak and the higher velocities represented by a long tail. Other fragmentation mechanisms in the wave-breaking phenomena may also alter the shape of velocity distribution, e.g. the bursting of surface bubbles (Lhuissier & Villermaux 2012; Berny et al. 2021), whose contribution to the production of droplets is known to be significant after the wave splashing phase.

7. Conclusions

We have investigated the collision and subsequent fragmentation of perturbed liquid rims, focusing on the range of $120 \leq We \leq 280$ that allows for the growth and merging of transverse ligaments and the production of fine drops from such ligaments via the pinch-off mechanism. We look into different parts of the post-collisional liquid bulk as it evolves over time, and our key findings are summarised below.

Firstly, following the quasi-one-dimensional approach of Wang & Bourouiba (2017), we derive analytical solutions of the liquid velocity and free surface profiles for the vertically expanding lamella sheet, which are shown to be in good agreement with the numerical results. Capillary effects have been neglected in this model for sheet evolution, but prove significant for the dynamics of the bordering rim. We then compared the growth of the rim position and thickness with the theoretical model of Gordillo et al. (2019), and develop scaling laws collapsing the data.

Secondly, we analyse the behaviour of transverse ligaments on top of the lamella rim. The ligaments produce fragments primarily via the end-pinching mechanism, and when the initial perturbation waveform is polychromatic, they migrate on the rim and merge with each other to form thicker and more corrugated ligaments, thus preventing their absorption into the rim and sustaining the fragmentation process. A novel scaling model is derived for predicting the evolution of ligament number density based on the migration speed model of Wang & Bourouiba (2018).

Lastly, we present the size and velocity statistics associated with the rim collision phenomenon. An excellent agreement between our fragment size distribution and the experimental results of Néel et al. (2020) is found within the range of grid convergence ($r \geq 4\varDelta _{10}$). The fragment size distribution becomes insensitive to the initial configurations when $We$ or $N_{max}$ further increases, which can be described using a power law with a break in slope. A theoretical model is proposed predicting the power-law distribution observed for $r \leq 0.2R_0$. Over time, the fragment speed $u_{frag}$ develops a largely uniform spread over the range of $0 \leq u_{frag} \leq 3.2 U_0$ as their parent ligaments continue to grow vertically and decelerate to form slower drops.

The implications of the present work are manifold. Firstly, it sheds new light on the fluid physics involved in a liquid impact problem that has not received much attention prior to the recent works of Néel et al. (2020) and Agbaglah (2021). Furthermore, the results we obtained are also of reference value for ongoing research works on spherical drop impacts, especially the influence of initial perturbations on the fragmentation process, which may also be present during the early-time prompt splashing phenomena observed in previous experimental studies (Burzynski et al. 2020; Wang et al. 2023). Lastly, this work serves as a stepping stone towards understanding the secondary splashing phenomenon observed in wave-breaking events and the associated fragment statistics (Mostert et al. 2022; Erinin et al. 2023a), and provides the basis for investigating the influence of other physical mechanisms not covered in the present work, e.g. viscosity, gravity and air-phase turbulence effects.