3.1. Splashing morphology

Figure 9.

Illustration of the different splashing behaviours for droplets with varying aspect ratios at (a) $ \textit{We} = 117$ , (b) $ \textit{We} = 467$ and (c) $ \textit{We} = 729$ . Representative time instances are selected to highlight the crown–lamella evolution in spreading, splashing type-1 and splashing type-2 regimes. Here, $h = 2R_{\textit{eq}}$ , and the fluid system corresponds to a water–air system.

We begin the presentation of our results by demonstrating the splashing morphologies exhibited by prolate ( $A_r = 0.75$ ), spherical ( $A_r = 1$ ) and oblate ( $A_r = 1.5$ ) drops upon impacting a liquid film. Figure 9 $(a{-}c)$ presents representative cases illustrating the distinct splashing behaviours observed for $A_r = 0.75$ , 1.0 and 1.5 at $ \textit{We} = 117$ ( $\tau = 9.52$ ), $ \textit{We} = 467$ ( $\tau = 28.57$ ) and $ \textit{We} = 729$ ( $\tau = 59.5$ ), respectively. It can be seen that, at $ \textit{We}=117$ , the spreading regime features a radially expanding crown with minimal rim undulations. As the aspect ratio increases, particularly in oblate drops, these undulations become more prominent as shown in figure 9( $a$ ). At $ \textit{We} = 467$ and $ \textit{We} = 729$ , two distinct splashing behaviours, splashing type-1 and splashing type-2, are observed, as shown in figures 9( $b$ ) and 9( $c$ ), respectively. These splashing modes are governed by different underlying mechanisms, which are discussed below.

Figure 10.

Formation of holes due to crown sheet rupture leading to unstable splashing at $\tau = 35.2$ for $A_r = 2.0$ and $ \textit{We} = 1050$ obtained using two mesh resolutions: ( $a$ ) $\Delta x = 0.0149R_{\textit{eq}}$ and ( $b$ ) $\Delta x = 0.03R_{\textit{eq}}$ . Panels ( $c$ ) and ( $d$ ) depict the dynamics before and after hole rupture at $\tau = 22.61$ and $\tau = 25$ , respectively, for $A_r = 2.0$ and $ \textit{We} = 729$ . Similarly, panels ( $e$ ) and ( $f$ ) show the crown dynamics before and after hole collapse at $\tau = 30.28$ and $\tau = 33.14$ , respectively, for $A_r = 2.0$ and $ \textit{We} = 1050$ . Here, the fluid system considered is a water–air system.

A key feature common to all splashing impacts is the formation of holes in the crown lamella just beneath the crown rim, where the film is thinnest. These ruptures, referred to as holes, are observed across all aspect ratios; however, they become more pronounced for $A_r \gt 1$ . Figure 10( $a$ ) illustrates this phenomenon for $A_r = 2.0$ and $ \textit{We} = 1050$ . The holes either collapse (figure 10 $e{,}f$ ) or result in finger splitting (figure 10 $c{,}d$ ), introducing variability in finger quantification and potential droplet detachment due to the momentum imparted onto the fingers by these ruptured holes. To confirm that these features are physical and not numerical artefacts, we have conducted additional simulations with significantly finer mesh resolution (figure 10 $b$ ). Cossali et al. (Reference Cossali, Marengo, Coghe and Zhdanov2004) reported that crown lamella thickness is non-uniform, typically ranging from $0.07R_{\textit{eq}}$ to $0.12R_{\textit{eq}}$ in experiments across various $ \textit{We}$ and film thicknesses. Accordingly, the minimum grid size was selected to be well below this range, with AMR applied sufficiently in advance of hole formation. The consistent appearance of holes at the same spatial locations and times in coarse and refined meshes confirms that these ruptures are physical phenomena and not artefacts of numerical resolution. These holes are likely triggered by undulations beneath the crown rim, which cause localized rupture in the thinning lamella. In contrast, other studies (Marston et al. Reference Marston, Truscott, Speirs, Mansoor and Thoroddsen2016; Burzynski & Bansmer Reference Burzynski and Bansmer2018) have attributed similar ruptures to the bursting of air bubbles entrapped just before impact and travelling through the crown lamella; however, those rupture sites typically appear farther from the rim. Notably, our simulations show no evidence of such bubble entrapment prior to impact. Upon rupture, these holes generate jets that carry substantial momentum. The intensity of this jetting governs the breakup behaviour. When the momentum released during hole collapse is moderate, the jets thicken and merge into the rim without further detachment. In contrast, when the released momentum is high, the jets elongate and eventually pinch off to form droplets. Thus, hole rupture not only initiates fingers but also provides the momentum that drives droplet ejection, making hole-induced detachment a common feature of the splashing regime. Additional breakup may occur if other instabilities act on the rim. The distinction between these instabilities is discussed next in figure 11.

Figure 11.

Temporal evolution of the crown rim morphologies for a drop with $A_r = 1.5$ , where panels (a) and (c) correspond to the splashing type-1 ( $ \textit{We} = 467$ ) and panels (b) and (d) represent the splashing type-2 ( $ \textit{We} = 729$ ).

Figure 11 compares the two splashing regimes that emerge from the interplay of hole-driven droplet detachment and rim deceleration dynamics. In both splashing regimes, i.e. splashing type 1 and type 2, the collapse of holes beneath the rim produces fingers that can detach droplets. This hole-induced mechanism is therefore common to both type 1 and type 2 splashing. The difference lies in the behaviour of the crown rim itself. As the crown expands radially, the rim experiences momentum loss, causing the rim of the crown to decelerate. This deceleration makes the rim susceptible to RT instability. In the splashing type-1 regime ( $ \textit{We} = 467$ , figure 11 $a$ ), the rim deceleration is weak; the RT growth rate remains low, and perturbations generated along the rim decay without causing further breakup. In contrast, in the splashing type-2 regime ( $ \textit{We} = 729$ , figure 11 $b$ ), the rim undergoes much stronger deceleration, amplifying RT perturbations and producing pronounced rim undulations that evolve into elongated fingers. These fingers elongate vertically and detach near their tips, forming secondary droplets. Thus, while both regimes exhibit hole-induced jet ejection, only type 2 shows rim deceleration-driven detachment.

Figure 12.

$(a)$ Formation of a canopy is observed for $A_r = 0.5$ , $ \textit{We} = 1050$ . Similar crown structures have been reported previously $(b)$ by Worthington (Reference Worthington1908) for a water-diluted milk droplet impacting a water surface at $ \textit{We} \approx 1430$ , and $(c)$ by Dandekar et al. (Reference Dandekar, Shen, Naar and Bourouiba2025) in an air–water system for a droplet impacting a water pool at $ \textit{We} = 850$ .

Another striking phenomenon observed during drop impact is canopy formation, which occurs when the vertical growth rate of the crown surpasses its radial expansion. This leads to the emergence of a hollow, inward-falling canopy structure, as shown in figure 12. At higher impact velocities, the crown rises rapidly, with the rim directed upward rather than outward. While RT instabilities along the rim continue to initiate cusp formation and subsequent finger development, the breakup of these fingers into secondary droplets continues to be governed by RP instability. Unlike the classical outward-expanding crown, the inward collapse in the canopy regime significantly alters the fragmentation dynamics. Canopy formation was first documented by Worthington (Reference Worthington1908) and has since been reported by several researchers, including Wang et al. (Reference Wang, Liu, Bayeul L.A., Murphy, Katz and Coutier2023) and Dandekar et al. (Reference Dandekar, Shen, Naar and Bourouiba2025), primarily in the context of spherical droplet impacts on deep liquid pools, where the unconfined vertical space permits substantial crown rise. In contrast, our simulations reveal canopy formation even on a liquid film when the impacting drop is prolate. The slender geometry of prolate drops results in a smaller initial contact area and narrower crown base diameter ( $D_{crown}$ ), favouring rapid vertical expansion over radial expansion. Similar behaviour was noted by Thoraval, Li & Thoroddsen (Reference Thoraval, Li and Thoroddsen2016) and Deka et al. (Reference Deka, Ray, Biswas, Dalal, Tsai and Wang2017), who showed that prolate drops are more prone to canopy formation than spherical or oblate drops, particularly in horizontally confined geometries at lower Weber numbers ( $ \textit{We} \approx 30$ –100). At these low values of $ \textit{We}$ , the resulting canopy exhibits minimal rim undulations and lacks pronounced cusp development, leading to what is often termed a spreading canopy. Eshraghi, Jung & Vlachos (Reference Eshraghi, Jung and Vlachos2020) further showed that canopy (splash-curtain) closure is strongly influenced not only by surface tension but also by the pressure difference between the air and the cavity, and they linked this closure to subsequent jet formation. Their study concluded that the key parameter governing canopy formation is the velocity of the air rushing into the cavity, rather than the impact velocity itself.

Figure 13.

Regime map demarcating the coalescence dynamics and splashing morphology of a non-spherical drop impacting a liquid pool for various values of aspect ratio ( $A_r$ ) and Weber number ( $ \textit{We}$ ), with $h = 2.0$ . The regime map also includes results from prior studies for reference and comparison. The parameters and configurations considered in previous studies are summarized as follows. Wang, Chen & Hwang (Reference Wang, Chen and Hwang2002): $h = 2.0$ , $A_r=1$ ; Cossali et al. (Reference Cossali, Coghe and Marengo1997): $h = 1.4$ , $A_r=1$ ; Wang & Chen (Reference Wang and Chen2000): $h = 1.0$ , $A_r=1$ . Note that Deka et al. (Reference Deka, Ray, Biswas, Dalal, Tsai and Wang2017) and Anirudh et al. (Reference Anirudh, Behera and Sahu2024) examined the impact of non-spherical water droplets on a liquid–air interface, and Murphy et al. (Reference Murphy, Li, d’Albignac, Morra and Katz2015) investigated the impact of spherical droplets on an oil–water film. Panels (b)–(e) present magnified views for $A_r = 1$ .

The regime map depicted in figure 13 delineates the coalescence dynamics and splashing morphology of a non-spherical drop impacting a liquid film as a function of the aspect ratio ( $A_r$ ) and Weber number ( $ \textit{We}$ ). The film thickness is fixed at $h = 2.0$ in our simulations. We perform simulations for discrete values of aspect ratio $A_r = 0.2$ , 0.35, 0.5, 0.75, 1.0, 1.5 and 2.0 across Weber numbers $ \textit{We} = 117$ , 262, 467, 729 and 1050. To ensure complete coverage of the regime map, additional simulations are carried out for $A_r = 1.5$ and 2.0 at $ \textit{We} = 502$ ; $A_r = 1.0$ and 1.5 at $ \textit{We} = 527$ ; $A_r = 0.75$ and 1.0 at $ \textit{We} = 591$ ; and $A_r = 0.5$ and 0.75 at $ \textit{We} = 1232$ and 1429. The regime map (figure 13) also incorporates results from previous studies for reference and comparison. Among these, while Deka et al. (Reference Deka, Ray, Biswas, Dalal, Tsai and Wang2017) and Anirudh et al. (Reference Anirudh, Behera and Sahu2024) investigated the impact of non-spherical water drops on a liquid–air interface, the other earlier studies focused exclusively on spherical drops. Notably, Deka et al. (Reference Deka, Ray, Biswas, Dalal, Tsai and Wang2017) examined the dynamics of a non-spherical drop impacting a restricted pool and reported that canopy structures form more readily in prolate drops than in oblate drops, an observation that is similar to that observed in the present study. In addition, Cossali et al. (Reference Cossali, Coghe and Marengo1997), Wang & Chen (Reference Wang and Chen2000) and Wang et al. (Reference Wang, Chen and Hwang2002) explored the influence of film thickness ( $h$ ) on splashing dynamics, albeit for spherical drops ( $A_r = 1$ ). To maintain visual clarity, magnified views of the regime map are shown in figure 13( $b{-}e$ ) for $A_r = 1.0$ .

In figure 13, we observe that at $ \textit{We} = 117$ , droplets of all aspect ratios undergo complete spreading, merging smoothly with the liquid film without generating rim instabilities or secondary structures. As $ \textit{We}$ increases, secondary droplets begin to detach from fingers through hole rupture (splashing type-1) for $A_r\gt 0.5$ . At $ \textit{We} = 467$ , drops with $A_r = 0.75$ , 1.0, 1.5 and 2.0 display splashing type-1, while those with $A_r = 0.35$ and 0.5 continue to spread. At this Weber number, the drop with $A_r = 0.2$ transitions to the canopy regime. At $ \textit{We} = 729$ , splashing type-2 emerges for $A_r \geqslant 0.75$ , whereas $A_r = 0.5$ remains in splashing type-1, and $A_r = 0.35$ and 0.2 show canopy formation. Additional simulations with smaller $ \textit{We}$ increments in the range $467 \lt We \lt 729$ reveal that higher aspect ratios promote stronger finger growth, eventually leading to splashing type-2. At $ \textit{We} = 1050$ , canopy structures appear for $A_r = 0.5$ , 0.35 and 0.2. The canopy region increases as we increase the $ \textit{We}$ . At higher Weber numbers, canopy formation persists at $ \textit{We} = 1232$ for $A_r = 0.5$ but is absent for $A_r = 0.75$ . In contrast, at $ \textit{We} = 1429$ , both $A_r = 0.5$ and 0.75 exhibit canopy behaviour. Khan et al. (Reference Khan, Jin and Yang2024) experimentally showed that oblate drops are more susceptible to splash-induced fragmentation when impacting a water film, consistent with our observation that splashing becomes more pronounced at larger $A_r$ . The dotted lines in figure 13 represent notional boundaries separating the different regimes, highlighting the coupled influence of drop shape and impact inertia on the resulting dynamics.

Figure 14.

Temporal evolution of the splashing dynamics for droplets with aspect ratios $A_r = 2.0$ (a), 1.0 (b) and 0.5 (c), impacting a liquid film of height $h = 2$ at $ \textit{We} = 729$ .

To gain physical insight into the observation that splashing becomes more prominent for drops with higher aspect ratios ( $A_r$ ), figure 14 presents the temporal evolution of splashing dynamics for drops with $A_r = 2.0$ , 1.0 and 0.5 at $ \textit{We} = 729$ . At $\tau = 0$ , the drop contacts the liquid film, initiating radial crown expansion. It can be seen that due to morphological differences, the oblate drop ( $A_r = 2.0$ ) forms a broader crown than the prolate drop ( $A_r = 0.5$ ), as evident at $\tau = 4.76$ in figure 14. As the crown expands, fluid accumulates at the rim, leading to thickening. The rim acts like a cylindrical jet and becomes unstable due to the influence of surface tension, triggering undulations via the RP instability. The narrower crown of the prolate drop decelerates more slowly than the oblate and spherical drops, resulting in less fluid accumulation at the rim and yielding a thinner structure. This, in turn, affects the RP instability by producing shorter azimuthal wavelengths for prolate drops, as compared with their oblate and spherical counterparts. By $\tau = 17.86$ , the rim undulations become pronounced (mostly fully developed). Subsequently, as rim deceleration intensifies, the RT instability amplifies these undulations into fingers, which are clearly visible in figure. 14. Since rim deceleration is slower for prolate drops, the RT instability grows more slowly, delaying finger evolution. This pattern holds across all prolate drops ( $A_r \lt 1.0$ ), which consistently exhibit delayed fingering compared with spherical and oblate drops. For $A_r = 0.5$ , finger formation occurs much later than shown and is not captured within the current visualization timeframe. As seen in figure 14, at $\tau = 47.62$ , fingers are well developed in oblate and spherical drops, while those in prolate drops remain subdued, indicating reduced instability growth.

The temporal evolution of splashing reveals morphological differences across various non-spherical drop shapes. However, it does not provide a quantitative understanding of the associated mechanisms responsible for node formation and finger growth. To address this, we employ linear stability analysis to identify the instabilities that drive node initiation and development induced by crown geometry. The governing parameters for the stability analysis, such as rim radius, rim deceleration and crown diameter, are accurately extracted from our numerical simulations and incorporated into the linear stability framework to predict the number of nodes that evolve into fingers and their corresponding growth rates. The details of the linear stability analysis are presented in the next section.

3.2. Linear stability analysis

The linear stability analysis examines the behaviour of an infinitesimally small disturbance imposed on the rim of the crown. Figure 15(a–c) illustrates the typical crown shape obtained from numerical simulation, the unperturbed (base state) and perturbed (base state plus disturbance) cross-sections of the rim, respectively. To examine the base state and its stability characteristics, we employ a Cartesian coordinate system $(x, y)$ , with the rim centreline located at $y = Y_0(t)$ , as shown in figure 15 $(b)$ . The time-dependent cross-sectional area of the rim is given by $a_0 = \pi {r_0}^2$ , where $r_0$ is the cross-sectional radius of the rim. The vertical rise velocity of the rim can be represented as $V_0 = {\rm d}Y_0/{\rm d}t$ .

Figure 15.

$(a)$ A typical crown obtained from numerical simulations. $(b)$ The unperturbed base state of the rim, where $t_{\kern-1.5pt f}$ denotes the crown lamella thickness and $V_{\kern-1.5pt f}$ is the velocity of the crown lamella. $(c)$ The perturbed state of the rim (base state plus disturbance). Here, $a$ represents the cross-sectional area of the cylindrical structure resting on a thin liquid sheet, and $V$ is the rim velocity. The subscript 0 indicates quantities associated with the unperturbed base state, whereas perturbed quantities used in the linear stability analysis are denoted without subscripts.

3.2.1. Base state

We model the unperturbed crown rim as a cylindrical structure resting at the edge of the sheet formed due to capillary forces (Yarin & Weiss Reference Yarin and Weiss1995; Roisman et al. Reference Roisman, Horvat and Tropea2006). The rim is treated as a quasi-one-dimensional liquid structure that grows due to mass influx from the lamella. Accordingly, the evolution of its cross-sectional area is governed by the mass conservation equation at $y=Y_0$ and can be expressed as

(3.1) \begin{align} \frac {{\rm d}a_0}{{\rm d}t} = t_{\kern-1.5pt f}(V_{\kern-1.5pt f} - V_0), \end{align}

where $t_{\kern-1.5pt f}$ denotes the crown lamella thickness, $V_{\kern-1.5pt f}$ is the velocity within the crown lamella at $y=Y_0$ , and $V_0$ represents the upward velocity of the rim. Furthermore, the rim dynamics are governed by the balance between inertia, capillary pressure and the momentum flux from the crown lamella, leading to the momentum balance equation

(3.2) \begin{align} a_0 \frac {{\rm d}{V_0}}{{\rm d} t} = t_{\kern-1.5pt f}(V_{\kern-1.5pt f} - V_0)^2 - \frac {2\sigma }{\rho }, \end{align}

where $\sigma$ is the surface tension and $\rho$ is the liquid density.

3.2.2. Perturbation equations

To perform the linear stability analysis, we perturb the base state with an infinitesimally small transverse disturbance. The evolution of the centreline of the resulting perturbed rim, $Y(x,t)$ , is then analysed. This imposed transverse disturbance induces a velocity component $u$ in the $x$ -direction within the rim, which is much smaller than the vertical velocity $V$ . As a result, the total velocity field is expressed as

(3.3) \begin{align} v = u \boldsymbol{e}_x + V \boldsymbol{e}_y, \end{align}

where $\boldsymbol{e}_x$ and $\boldsymbol{e}_y$ are the unit vectors in the $x$ and $y$ directions, respectively. In the linear stability framework, the wavelength of the long-wave perturbation $(\lambda )$ satisfies

(3.4) \begin{align} | Y - Y_0 | = \epsilon \ll \lambda \quad \text{and} \quad r_0 \ll \lambda , \end{align}

indicating that the perturbation amplitude $\epsilon$ is much smaller than its wavelength. This also ensures that

(3.5) \begin{align} \left | \frac {\partial Y}{\partial x} \right | \ll 1. \end{align}

Under this approximation, the $x$ -component of the velocity $u$ , the velocity gradient $\partial V / \partial x$ and the gradient of the angular velocity of the rim cross-section $\partial \varOmega / \partial x$ are assumed to be sufficiently small, rendering all nonlinear convective terms in the momentum equation second-order in magnitude, and they are neglected in the linear stability analysis (Entov & Yarin Reference Entov and Yarin1984). Here, $\varOmega = {\partial V/\partial x}$ . The resulting linearized perturbation equations are given by

(3.6a) \begin{align} \frac {\partial a}{\partial t} + {a_0} \frac {\partial u}{\partial x} - {t_{\kern-1.5pt f}}( V_{\kern-1.5pt f} - V) = 0, \\[-28pt] \nonumber \end{align}

(3.6b) \begin{align} \rho {a_0} \frac {\partial u}{\partial t} - \frac {\partial P}{\partial x} - \left [ - \rho {t_{\kern-1.5pt f}} (V_{\kern-1.5pt f} -V) u + 2 \sigma \frac {\partial Y}{\partial x} \right ] = 0, \\[-28pt] \nonumber \end{align}

(3.6c) \begin{align} \rho a \frac {\partial V}{\partial t} - P \kappa - \frac {\partial Q}{\partial x} - \big [ \rho {t_{\kern-1.5pt f}} (V_{\kern-1.5pt f} - V)^2 - 2\sigma \big ] = 0, \\[-28pt] \nonumber \end{align}

(3.6d) \begin{align} \rho \frac {\partial L}{\partial t} - \frac {\partial M}{\partial x} - Q - {m_S} = 0. \\[0pt] \nonumber \end{align}

Here, $\kappa$ denotes the curvature of the centreline of the rim, $Q$ is the shear force in the $y$ -direction acting on the cross-section of the rim and $m_S$ represents the distributed moment due to external forces. The magnitude of the longitudinal tensile force ( $P$ ), which accounts for both the capillary pressure within the rim and the contribution from surface tension, is given by

(3.7) \begin{align} P = \pi \sigma r + \sigma {a_0} \frac {\partial ^2 r}{\partial x^2}. \end{align}

The angular momentum of the rim per unit length in the plane ( $L$ ) is given by

(3.8) \begin{align} L = I_0 \varOmega , \quad \text{where} \quad I_0 = \frac {\pi {r_0}^4}{4}. \end{align}

The internal moment resulting from stresses within the cross-section ( $M$ ), produced by the pressure gradient in the $y$ -direction in the rim cross-section due to its acceleration, is given by

(3.9) \begin{align} M = -\rho I_0 \frac {{\rm d} V}{{\rm d} t}. \end{align}

The moment due to the flow from the crown lamella into the rim is given by

(3.10) \begin{align} m_S = -\rho {t_{\kern-1.5pt f}} {r_0} (V_{\kern-1.5pt f} - V_0) u. \end{align}

3.2.3. Normal mode approximation

We investigate the temporal stability of the base state to infinitesimal two-dimensional disturbances by employing a normal mode approximation. Here, each flow variable is expressed as the sum of the base state and a two-dimensional perturbation, which is further decomposed into an amplitude and a wave-like component as follows:

(3.11) \begin{align} (Y,r,u,Q)(x,y,t) = [Y_0, {r_0}, 0, 0]+ \left (\hat {\epsilon }, \hat {\eta }, \hat {u}, \hat {q} \right ) \exp (\omega _c t + \textrm {i} \alpha x), \end{align}

where, the hats denote the perturbation quantities, $\alpha = 2 \pi /\lambda$ is the (real) wavenumber and $\omega _c$ is the complex frequency, with $\omega = \text{Real}(\omega _c)$ representing the growth rate of the disturbance. For disturbances of the form given in (3.11), $\omega \gt 0$ corresponds to a temporally unstable mode, $\omega \lt 0$ indicates a stable mode and $\omega = 0$ denotes a neutrally stable disturbance. Note that the quantities $r_0$ , $a_0$ and $I_0$ of the rim vary with time; however, their temporal evolution is assumed to be sufficiently slow compared with that of the disturbances. Hence, a quasisteady-state approach is adopted by perturbing the base state at a frozen time. Following the standard procedure of linear stability analysis, i.e. substituting (3.11) into the governing equations, subtracting the base state equations and subsequently linearizing the resultant equations, we obtain the governing dispersion equation for the linear stability analysis. In dimensionless form, the dispersion relation for the growth rate of the disturbance can be written as

(3.12) \begin{align} & 2 \big(2 + {\bar {\alpha }}^2 \big)\pi ^2 {\bar {\omega }}^4 + 4 \big(3 + {\bar {\alpha }}^2 \big)\pi {\bar {W}}_0 {\bar {\omega }}^3 \nonumber \\[2pt]& + \left [ 2\big(4 + {\bar {\alpha }}^2\big) {\bar {W}}_0^2 - 2\pi {\bar {t}_f} \big(2 + {\bar {\alpha }}^2\big) {{\rm d} {\bar {V}_0} \over {\rm d}t} + \pi ^2 {\bar {\alpha }}^2 \big(2 + {\bar {\alpha }}^2 + {\bar {\alpha }}^4\big) \right ] {\bar {\omega }}^2 \nonumber \\[2pt]& + \left [ 2 {\bar {\alpha }}^2 \left (4 - {\bar {t}_f} - {\bar {t}_f} {d {\bar {V}_0} \over dt} \right ) - 4 {\bar {t}_f} {{\rm d} {\bar {V}_0} \over {\rm d}t} + \pi {\bar {\alpha }}^6 + {\bar {\alpha }}^4 (3\pi + 2 {\bar {t}_f}) \right ] {\bar {W}}_0 {\bar {\omega }} \nonumber \\[2pt]& - 2\pi ^2 {\bar {\alpha }}^4 \big(1 - {\bar {\alpha }}^2\big) + 4\pi {\bar {\alpha }}^2 \big(2 + {\bar {\alpha }}^2\big){{\rm d} {\bar {V}_0} \over {\rm d}t} = 0. \end{align}

Here,

(3.13) \begin{align} \omega & = \sqrt { \sigma \over \rho {r_0}^3} {\bar {\omega }}, \quad \alpha = {{\bar {\alpha }} \over r_0}, \quad t_{\kern-1.5pt f} = {\bar {t}_f} r_0, \quad {\dot {V}_0} ={\sigma \over \rho {r_0}^2} {\bar {\dot {V}}}\quad \textrm { and }\nonumber \\[2pt] {\bar {W}}_0 & = \rho t_{\kern-1.5pt f} (V_{\kern-1.5pt f} - V_0) / \sqrt {\rho \sigma r_0}, \end{align}

wherein ${\dot {V}_0} = {{\rm d} V_0 / {\rm d} t}$ is the acceleration of the rim. The ‘bar’ superscript denotes dimensionless quantities, which are omitted in the subsequent analysis for clarity. The current stability analysis adopts the same methodology as that described in Entov & Yarin (Reference Entov and Yarin1984) and Roisman et al. (Reference Roisman, Horvat and Tropea2006).

The dispersion relation (3.12) contains two unknown parameters, namely the rim radius ( $r_0$ ) and the rim acceleration ( $\dot {V}_0$ ). Previous studies on the linear stability of node formation during splashing have proposed various approximations for these parameters, each with inherent limitations. Zhang et al. (Reference Zhang, Brunet, Eggers and Deegan2010) estimated $r_0$ by fitting circles to the inner and outer edges of the crown rim. Although this approach is effective during the initial stages of crown formation, it fails to represent the rim geometry accurately at later times when significant deformations occur. Agbaglah & Deegan (Reference Agbaglah and Deegan2014) used axisymmetric numerical simulations to estimate $r_0$ and $\dot {V}_0$ . However, since splashing is fundamentally an asymmetric phenomenon, the axisymmetric assumption reduces the accuracy and applicability of their results. Additionally, Agbaglah et al. (Reference Agbaglah, Josserand and Zaleski2013) conducted simulations on a crown segment within a computational domain spanning two perturbation wavelengths, introducing localized disturbances to investigate finger growth and droplet detachment. This approach does not capture the global temporal evolution of the crown or related phenomena such as hole formation, thus neglecting crucial aspects of rim dynamics. Using potential flow theory, Bisighini et al. (Reference Bisighini, Cossali, Tropea and Roisman2010) analytically modelled the evolution of a crater as a growing and translating sphere within a liquid pool. Their approach can, in principle, be extended to estimate crown dynamics, including $r_0$ and $\dot {V}_0$ , using conservation laws, provided the crown thickness ( $t_{\kern-1.5pt f}$ ) is known. The present study focuses on impacts involving non-spherical drops, where the drop shape deviates significantly from a sphere. Extending the potential flow formulation to accommodate a growing and moving ellipsoid introduces substantial mathematical and computational challenges, making such an approach impractical for our purposes. Thus, in contrast to the previous approaches, we extract $r_0$ and $\dot {V}_0$ directly from numerical simulation results at the exact time when the number of nodes is to be determined. The rim radius is calculated by sampling data at different locations around the rim, as illustrated in figure 16, and averaging the values to account for spatial variations. The rim acceleration is obtained by tracking the vertical position of the rim ( $Y_0$ ) over time and applying a second-order central difference scheme to compute its second derivative. To minimize errors, we restrict the analysis to rim regions that have not yet developed fingers, as finger formation distorts the accuracy of both $r_0$ and $\dot {V}_0$ . It is to be noted that the geometric approximations like those used by Zhang et al. (Reference Zhang, Brunet, Eggers and Deegan2010), which include the entire crown structure, including fingers, are prone to such inaccuracies. However, this approach allows the effects of hole rupture or collapse, the events that dynamically alter the rim structure, to be incorporated when extracting $r_0$ and $\dot {V}_0$ . Following the convention of Roisman et al. (Reference Roisman, Horvat and Tropea2006), the crown thickness ( $t_{\kern-1.5pt f}$ ) is assumed to be $0.5r_0$ . Using these extracted parameters from the numerical simulations, we solve (3.12) to determine the growth rate ( $\omega$ ) as a function of wavenumber ( $\alpha$ ). From (3.12), the dispersion relation depends solely on the rim deceleration, $\dot {V_0}$ . The critical growth rate, $\omega _{\textit{cr}}$ , shows the strongest sensitivity to $\dot {V_0}$ , whereas the corresponding wavenumber, $\alpha _{\textit{cr}}$ , is less sensitive. The rim deceleration obtained from the simulations is substituted into the dispersion relation to construct the dispersion curves, from which the maximum growth rate, $\omega _{\textit{cr}}$ , and the associated wavenumber, $\alpha _{\textit{cr}}$ , are determined. Using these quantities together with the rim radius $(r_0)$ and the crown diameter $(D_{crown})$ , the number of fingers is estimated as $n = (\alpha _{\textit{cr}}D_{crown}) /(2 r_0)$ .

Figure 16.

( $a$ ) Schematic showing the locations of slices used to extract the rim radius ( $r_0$ ) and rim height ( $Y_0$ ). ( $b$ ) Illustration of the procedure for extracting $r_0$ and $Y_0$ from one of the slices shown in panel ( $a$ ).

3.3. Competition between RT and RP instabilities

To further highlight the effect of aspect ratio, a general comparison of the input parameters ( $r_0$ , $\dot {V_0}$ ) with aspect ratio ( $A_r$ ) is made at three different time steps and is shown in figure 17. Additionally, figures 18(a) and (b) illustrate the variations of $\omega_{cr}$ and $\alpha_{cr}$ with $A_r$ at different times during the crown formation stage. The parameters $r_0$ and $\dot {V_0}$ exhibit similar trends with $A_r$ at different time steps. As shown in figure 14, the rim deceleration $\dot {V_0}$ increases with $A_r$ , such that prolate drops experience a smaller deceleration than oblate ones, consistent with figure 17 $(a)$ . The rim radius $r_0$ is governed by two opposing effects: stronger deceleration promotes greater fluid accumulation and increases $r_0$ , whereas a larger crown circumference associated with higher $A_r$ distributes the accumulated liquid over a wider perimeter, reducing $r_0$ . The combined influence of these effects yields the trend in figure 17 $(b)$ , where spherical drops exhibit the largest rim radius. These variations in $r_0$ and $\dot {V_0}$ directly influence the instability behaviour discussed earlier. The growth rate is primarily controlled by the rim deceleration, which increases with $\dot {V_0}$ , while $r_0$ determines the dominant wavenumber and hence the number of nodes that develop along the rim.

Figure 17.

( $a$ ) Variation of rim deceleration $(\dot {V}_0)$ and ( $b$ ) rim radius $(r_0)$ with aspect ratio $(A_r)$ for $ \textit{We} = 729$ at different time instances $\tau = 42.9$ , 47.6 and 52.4.

Figure 18.

(a) Variation of the critical growth rate $(\omega _{\textit{cr}})$ and (b) critical wavenumber $(\alpha _{\textit{cr}})$ with aspect ratio $(A_r)$ for $ \textit{We} = 729$ at different time instances $\tau = 42.9$ , 47.6 and 52.4.

There has been considerable debate regarding the primary mechanism responsible for node formation at the rim of a crown during splashing. Some studies attribute it to RP instability (Zhang et al. Reference Zhang, Brunet, Eggers and Deegan2010), others to RT instability (Krechetnikov & Homsy Reference Krechetnikov and Homsy2009), while some propose a combination of both these instability mechanisms (Roisman et al. Reference Roisman, Horvat and Tropea2006). In this section, we aim to identify the dominant instability mechanism governing node formation for the system considered in the present study. The RP instability manifests itself as symmetric perturbations along the rim axis, driven by variations in the axial velocity ( $u$ ) and the radius of the rim ( $r$ ). In contrast, perturbations in $Q$ and $Y$ lead to a bending of the rim centreline (Roisman et al. Reference Roisman, Horvat and Tropea2006). To isolate the RP instability, we introduce perturbations $\hat {\eta }$ and $\hat {u}$ , while setting $\hat {\epsilon } = 0$ and $\hat {q} = 0$ . Under these assumptions, the dispersion relation from (3.12) simplifies to

(3.14) \begin{align} 2\omega ^2 + W_0\omega + \alpha ^2 \big(\alpha ^2 - 1\big) = 0. \end{align}

In contrast, for analysing the RT instability, we introduce perturbations $\hat {\epsilon }$ and $\hat {q}$ , while keeping $\hat {\eta } = 0$ and $\hat {u} = 0$ . This is because the RT instability develops perpendicular to the rim axis, in the direction of its deceleration, and is triggered solely by perturbations in $Y$ and $Q$ (Roisman et al. Reference Roisman, Horvat and Tropea2006). Under these conditions, the dispersion relation derived from (3.12) simplifies to

(3.15) \begin{align} \pi \big(2+\alpha ^2 \big)\omega ^2 + W_0\big(4+\alpha ^2 \big)\omega + 2\pi \alpha ^2 = 0. \end{align}

Figure 19.

Dispersion curves ( $\omega$ versus $\alpha$ ) of the RP instability (dashed red line) and the combined RP + RT instability (solid black line) for $A_r=1.0$ , $ \textit{We}=1050$ at $\tau = 71.43$ .

The dispersion curve obtained from the simplified RP instability relation (3.14) is shown in figure 19. It exhibits a parabolic shape, with $\omega \gt 0$ over a finite wavenumber range, indicating linear instability. The curve exhibits a distinct ‘most dangerous’ mode, corresponding to the wavenumber ( $\alpha _{\textit{cr}}$ ) at which the growth rate is maximal ( $\omega _{\textit{cr}}$ ), and a ‘cutoff’ mode, beyond which disturbances decay ( $\omega \lt 0$ ). In contrast, the RT instability relation (3.15) yields no positive real roots for $\omega$ , implying that this mode is stable and does not independently lead to node formation. This observation is consistent with the findings of Roisman et al. (Reference Roisman, Horvat and Tropea2006). Hence, the formation of nodes is primarily attributed to the RP instability. However, when both RP and RT instabilities are considered by solving the full dispersion relation (3.12), the resulting growth rate is higher than that from the RP mode alone (see figure 19). Notably, this trend is further supported by the data in table 1, which shows that the most unstable wavelength remains nearly unchanged between the two cases when compared with the change in the growth rate. This behaviour is consistent across a range of aspect ratios ( $A_r$ ), as presented in table 1, which provides a quantitative comparison of $\omega _{\textit{cr}}$ and $\alpha _{\textit{cr}}$ at different instants for both instability modes. These results suggest that RP instability determines the dominant wavelength and number of fingers, while RT instability significantly alters the growth rate of the perturbations. Thus, RP instability governs the spatial pattern of node formation, while RT instability contributes primarily to their amplification. In addition to the planar crown assumption used in the linear analysis presented in this section, we have also performed a curved-rim stability analysis to assess the effect of crown curvature on the dispersion characteristics. Our results show no significant change in either the maximum growth rate ( $\omega _{\textit{cr}}$ ) or the corresponding wavenumber ( $\alpha _{\textit{cr}}$ ). A detailed discussion of this analysis has been included in Appendix A.

Table 1.

Comparison of the maximum growth rate ( $\omega _{\textit{cr}}$ ) and the critical wavenumber ( $\alpha _{\textit{cr}}$ ) obtained from (3.14) (RP mode) and (3.12) (combined RP and RT modes) for Weber numbers ( $ \textit{We} = 729$ ) for various aspect ratios ( $A_r$ ) at different time instants. Here, $D_{crown}$ denotes the crown diameter, ${\mathcal A} = (\omega _{\textit{cr}, \textit{RP}+\textit{RT}}-\omega _{\textit{cr}, \textit{RP}})/\omega _{\textit{cr}, \textit{RP}} ) \times 100$ and ${\mathcal B} = (\alpha _{\textit{cr}, \textit{RP}+\textit{RT}}-\alpha _{\textit{cr}, \textit{RP}})/\alpha _{\textit{cr}, \textit{RP}} ) \times 100$ . The analysis is conducted for a water–air fluid system.

Figure 20.

$(a)$ Dispersion curves for $A_r = 2.0$ at $ \textit{We} = 467$ , 729 and 1050. The number of nodes corresponding to the most unstable mode for $(b)$ $ \textit{We} = 467$ at $\tau = 25.7$ , $(c)$ $ \textit{We} = 729$ at $\tau = 36.9$ and $(d)$ $ \textit{We} = 1050$ at $\tau = 57.1$ . The identified nodes are highlighted in red and numbered accordingly. The fluid system considered here is a water–air system.

3.4. Parametric study

Figure 20( $a$ ) shows the dispersion curves ( $\omega$ versus $\alpha$ ) for various Weber numbers ( $ \textit{We}$ ) obtained from linear stability analysis. The respective dispersion curves indicate the values of $(\alpha _{\textit{cr}}, \omega _{\textit{cr}})$ corresponding to each Weber number. It can be seen that increasing $ \textit{We}$ increases the growth rate of the most unstable mode ( $\omega _{\textit{cr}}$ ), indicating that higher impact velocities intensify the instability of the crown rim. A higher $\omega _{\textit{cr}}$ implies that small perturbations on the rim grow more rapidly, thereby increasing the likelihood of node formation. This observation is consistent with the regime boundaries shown in figure 13, where higher $ \textit{We}$ values correspond to more frequent splashing events. These unstable nodes can evolve into fingers, eventually detaching as secondary droplets, which is a distinctive feature of the splashing process. Figure 20( $b{-}d$ ) show the snapshots obtained from the numerical simulations highlighting the number of nodes for $ \textit{We} = 467$ , 729 and 1050, respectively. The nodes are circled in red and numbered. The number of nodes observed for different Weber numbers agrees well with the theoretical prediction from the linear stability analysis. Interestingly, we observe that the number of nodes does not increase monotonically with $ \textit{We}$ , despite the higher growth rates predicted by the linear stability analysis. This suggests that while a higher $ \textit{We}$ increases the growth rate of the most unstable mode (thereby accelerating crown development), the node formation is affected by other factors, such as the rim radius ( $r_0$ ) and local crown irregularities, including hole ruptures, which can disrupt regular node formation.

Figure 21.

Evolution of the front view of the crown for $(a)$ $A_r = 0.75$ , $(b)$ $A_r = 1.0$ , $(c)$ $A_r = 1.5$ and $(d)$ $A_r = 2.0$ , for $ \textit{We} = 729$ at $\tau = 24$ , 36, 48. Here, $x$ and $y$ are dimensionless horizontal and vertical distances, respectively. The fluid system considered here is a water–air system.

Building on the insights from figure 20, which illustrated the relationship between the growth rate ( $\omega$ ), Weber number ( $ \textit{We}$ ) and the number of nodes ( $n$ ) along the crown rim, figure 21 presents the temporal evolution of crown morphology for different aspect ratios ( $A_r$ ) at $ \textit{We} = 729$ . This figure further highlights the influence of drop shape on crown dynamics by showing front-view snapshots of the crown at various non-dimensional times. Figure 21 $(a{-}d)$ illustrate how the crown shape evolves for aspect ratios $A_r = 0.75$ , 1.0, 1.5 and 2.0, respectively. At early times ( $\tau = 24$ , red dashed line), the crown forms a relatively smooth structure with small perturbations. As the time progresses to $\tau = 36$ (black solid line) and further to $\tau = 48$ (blue dash–dotted line), the crown develops more pronounced undulations. Additionally, the temporal difference in height reached by the rim of the crown is significantly higher for prolate drops compared with oblate drops. This suggests that oblate drops experience higher rim deceleration ( $\dot {V_0}$ ) over a longer duration compared with prolate drops. This amplifies the RT instability, increasing the growth rate of nodes and promoting subsequent finger formation. This trend is consistent with the observation in figure 13, which shows that an increase in $A_r$ makes the drop more prone to splashing.

Figure 22.

Comparison of the predicted number of nodes from the linear stability analysis (blue dash–dot line) with the observed number of nodes from numerical simulations (black solid line) for $ \textit{We} = 729$ at $(a)$ $A_r = 0.75$ , $(b)$ $A_r = 1.0$ and $(c)$ $A_r = 1.5$ .

It is important to note that while linear stability analysis provides a crucial estimate for the onset of splashing by predicting the number of nodes formed on the crown rim using the critical wavenumber ( $\alpha _{\textit{cr}}$ ), not all of these nodes evolve into fully developed fingers. Nonlinear effects, such as merging and splitting of nodes induced by holes in the crown (figure 10), play a crucial role in governing the subsequent evolution of the rim structure. Additionally, as shown in figure 21, while some of the initially formed nodes evolve into fingers, others diminish over time and do not develop into visible fingers. This damping effect results from the complex interplay between the local momentum distribution along the rim and the instability growth rate. Specifically, insufficient fluid accumulation in certain rim regions suppresses node growth, preventing them from evolving into fingers.

Figure 22 illustrates the temporal evolution of the nodes predicted by linear stability analysis for droplets with aspect ratios $A_r = 0.75$ , 1.0 and 1.5 at $ \textit{We}=729$ . It is observed that the number of nodes is not constant but decreases over time. This trend is clearly visible in figure 22, where a steady decline in node count can be observed, even as the crown diameter $D_{crown}$ continues to grow. The observed reduction in node count is attributed to the accumulation of fluid from the lamella into the rim, which increases the rim radius $(r_0)$ over time, as also seen in figure 17( $b$ ). As $r_0$ increases, the most unstable wavelength grows longer, resulting in fewer nodes forming around the rim. Notably, the number of nodes predicted by linear theory represents an upper bound, with not all predicted nodes necessarily evolving into visible fingers. The conversion of nodes into fingers depends on rim deceleration, which governs the growth rate of the RT instability. Additionally, the local redistribution of fluid within the rim further modulates this transformation, underscoring the complex interaction between rim geometry and the growth of instability. For the prolate case ( $A_r = 0.75$ ), the comparison in figure 22 shows that the predicted node count remains higher than the value observed in the simulations. This behaviour differs slightly from the spherical and oblate cases, where the predicted count is lower than the simulated value. This can be explained by the combined effect of the rim radius and crown diameter. Since the prolate drop has the smallest rim radius ( $r_0$ ), it exhibits a higher critical wavenumber ( $\alpha _{\textit{cr}}$ ), leading to a slight overprediction of the number of nodes in the linear stability analysis. At the same time, the crown diameter ( $D_{crown}$ ) is also the smallest, which leads to closer spacing between neighbouring fingers. As the crown evolves, these closely spaced fingers merge more readily, resulting in a smaller number of distinct nodes in the simulation. Also in prolate drops, weaker deceleration results in reduced fluid accumulation at the rim and a minor increase in $r_0$ over time (as observed in figure 17 $b$ ). Consequently, the most unstable wavelength remains nearly constant, and the number of predicted nodes exhibits minimal temporal variation. In contrast, oblate droplets experience more substantial deceleration, which drives more fluid into the rim, significantly increasing $r_0$ and reducing the number of nodes over time. This behaviour is consistently captured in both theoretical predictions and simulations. Initially, the oblate drop, despite having a higher rim radius compared with the prolate drop, has a higher number of nodes due to its increased rim circumference. However, due to the rapid increase in rim radius driven by high deceleration, the node count in the oblate case eventually decays much faster than that of the prolate and spherical drops.