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. (Reference Néel, Lhuissier and Villermaux2020) and Agbaglah (Reference Agbaglah2021), 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. Reference Pairetti, Popinet, Damián, Nigro and Zaleski2018).

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 Reference Liu and Bothe2016; Berny et al. Reference Berny, Deike, Popinet and Séon2022), and also taking into account the surface corrugation of the splash-up in wave-breaking events (Kiger & Duncan Reference Kiger and Duncan2012) 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. Reference Pal, Crialesi-Esposito, Fuster and Zaleski2021), 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. (Reference Zhang, Brunet, Eggers and Deegan2010) 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 Reference Popinet2019) 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 Reference Popinet2018; Tang et al. Reference Tang, Mostert, Fuster and Deike2021). Capillary effects are modelled as source terms in the Navier–Stokes equations using an adaptation of Brackbill's method (Brackbill, Kothe & Zemach Reference Brackbill, Kothe and Zemach1992; Popinet Reference Popinet2009), which calculates the interface curvature by taking the finite-difference discretisation of the derivatives of interface height functions (Popinet Reference Popinet2009). 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. Reference Pal, Crialesi-Esposito, Fuster and Zaleski2021) so that the cylinders break up prematurely.

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 Reference Yarin2006; Josserand & Thoroddsen Reference Josserand and Thoroddsen2016; Wang & Bourouiba Reference Wang and Bourouiba2017), 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 (Reference Gordillo, Riboux and Quintero2019). 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 (Reference Wang and Bourouiba2017). 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 (Reference Yarin and Weiss1995), Wang & Bourouiba (Reference Wang and Bourouiba2017) 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 (Reference Wang and Bourouiba2017) 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 Reference Wang and Bourouiba2017; Gordillo et al. Reference Gordillo, Riboux and Quintero2019).

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 Reference Savva and Bush2009). 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. (Reference Vledouts, Quinard, Vandenberghe and Villermaux2016) and marks the onset of bag film fragmentation in droplet aerobreakup problems (Ling & Mahmood Reference Ling and Mahmood2023; Tang et al. Reference Tang, Adcock and Mostert2023), 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 Reference Villermaux and Bossa2011; Riboux & Gordillo Reference Riboux and Gordillo2015), which may arise at larger $We$ values. We therefore do not take special measures to artificially stabilise (Liu & Bothe Reference Liu and Bothe2016) or perforate (Chirco et al. Reference Chirco, Maarek, Popinet and Zaleski2022) 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. Reference Wang, Dandekar, Bustos, Poulain and Bourouiba2018), 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. (Reference Gordillo, Riboux and Quintero2019), 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 (Reference Villermaux and Bossa2011) and Wang et al. (Reference Wang, Dandekar, Bustos, Poulain and Bourouiba2018) 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. (Reference Mongruel, Daru, Feuillebois and Tabakova2009), Thoroddsen, Takehara & Etoh (Reference Thoroddsen, Takehara and Etoh2012) and Visser et al. (Reference Visser, Frommhold, Wildeman, Mettin, Lohse and Sun2015) 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 (Reference Wang and Bourouiba2022) where the rim volume remains independent of $We$ due to their axisymmetric configuration.

5.1. Formation and growth

In the scenario considered by Wang & Bourouiba (Reference Wang and Bourouiba2021), 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 (Reference Gordillo, Onuki and Tagawa2020) 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 (Reference Wang and Bourouiba2018); this is still much later than the onset of the micro-splashing phenomena investigated by Thoroddsen et al. (Reference Thoroddsen, Takehara and Etoh2012), 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 (Reference Lai, Eggers and Deike2018) 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. (Reference Lai, Eggers and Deike2018) 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 (Reference Gordillo and Blanco-Rodríguez2023) 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. (Reference Lai, Eggers and Deike2018) 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. (Reference Lai, Eggers and Deike2018) 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 (Reference Gekle and Gordillo2010). ‘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 (Reference Wang and Bourouiba2022) 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 (Reference Wang and Bourouiba2021) 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 Reference Ghabache, Séon and Antkowiak2014), $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 (Reference Gekle and Gordillo2010), where the liquid particles travel at constant speed upwards before being slowed down once again at the bulb by capillary effects.

While Wang & Bourouiba (Reference Wang and Bourouiba2022) 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 (Reference Wang and Bourouiba2022) 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. (Reference Wildeman, Visser, Sun and Lohse2016) for drops impacting a smooth surface (see their figure 4) and Fudge et al. (Reference Fudge, Cimpeanu, Antkowiak, Castrejón-Pita and Castrejón-Pita2023) 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 (Reference Wang and Bourouiba2022) 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 Reference Riboux and Gordillo2014), and they generate vorticity (Batchelor Reference Batchelor2000; Li et al. Reference Li, Thoraval, Marston and Thoroddsen2018) 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. (Reference Fudge, Cimpeanu, Antkowiak, Castrejón-Pita and Castrejón-Pita2023) 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 (Reference Wang and Bourouiba2022) is most likely linked with the early-time lamella formation, which also features strong viscous effects and, according to Wildeman et al. (Reference Wildeman, Visser, Sun and Lohse2016) and Ó Náraigh & Mairal (Reference Ó Náraigh and Mairal2023), 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. Reference Wildeman, Visser, Sun and Lohse2016; Ó Náraigh & Mairal Reference Ó Náraigh and Mairal2023), which matches the observations of Wang & Bourouiba (Reference Wang and Bourouiba2022). 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. Reference Gordillo, Lhuissier and Villermaux2014; Wang & Bourouiba Reference Wang and Bourouiba2018). 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 (Reference Wang and Bourouiba2021) differs from the early experimental work of Thoroddsen & Sakakibara (Reference Thoroddsen and Sakakibara1998) 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 (Reference Marmanis and Thoroddsen1996) 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 (Reference Wang and Bourouiba2021) 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 (Reference Yarin and Weiss1995) 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 (Reference Yarin and Weiss1995) 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. (Reference Chan, Dodd, Johnson and Moin2021). 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 (Reference Wang and Bourouiba2018). These are below the theoretical value of 1.89 as predicted by the RP instability (Gordillo & Gekle Reference Gordillo and Gekle2010), 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. (Reference Gordillo, Lhuissier and Villermaux2014) and Wang & Bourouiba (Reference Wang and Bourouiba2018) 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).