Impact of droplets onto surfactant-laden thin liquid films

Abstract

We study the effect of insoluble surfactants on the impact of surfactant-free droplets onto surfactant-laden thin liquid films via a fully three-dimensional direct numerical simulation approach that employs a hybrid interface-tracking/level-set method, and by taking into account surfactant-induced Marangoni stresses due to gradients in interfacial surfactant concentration. Our numerical predictions for the temporal evolution of the surfactant-free crown are validated against the experimental work by Che & Matar (Langmuir, vol. 33, 2017, pp. 12140–12148). We focus on the ‘crown-splash regime’, and we observe that the crown dynamics evolves through various stages: from the growth of linear modes (through a Rayleigh–Plateau instability) to the development of nonlinearities leading to primary and secondary breakup events (through droplet shedding modulated by an end-pinching mechanism). We show that the addition of surfactants does not affect the wave selection via the Rayleigh–Plateau instability. However, the presence of surfactants plays a key role in the late stages of the dynamics as soon as the ligaments are driven out from the rim. Surfactant-induced Marangoni stresses delay the end-pinching mechanisms to result in longer ligaments prior to their capillary singularity. Our results indicate that Marangoni stresses bridge the gap between adjacent protrusions promoting the adjacent protrusions' collision and the merging of ligaments. Finally, we demonstrate that the addition of surfactants leads to surface rigidification and consequently to the retardation of the flow dynamics.

1. Introduction

Droplets impacting on liquid surfaces are observed in a broad range of natural phenomena and industrial applications. Consequently, this area of study drives significant scientific interest; for comprehensive reviews please see e.g. Yarin (2006), Josserand & Thoroddsen (2016) and Cheng, Sun & Gordillo (2022). The complex topological features of droplet impact have attracted the fluid mechanics community since the first ground-breaking experiments from Worthington (1908) and Edgerton (1977) over a century ago. However, it was not until the advent of high-speed imaging and high-resolution computational tools that it became possible to capture the early dynamics of the broad range of scales involved in such phenomena (Thoroddsen 2002; Thoraval et al. 2012; Thoroddsen et al. 2012; Zhang et al. 2012; Agbaglah & Deegan 2014; Josserand, Ray & Zaleski 2016; Fudge, Cimpeanu & Castrejón-Pita 2021).

The dynamics of droplet impact on liquids is extremely diverse and is determined by a wide number of factors including the impact speed and the fluid properties, as well as whether the impact occurs on either a deep pool or a thin layer of fluid (Deegan, Brunet & Eggers 2007; Josserand et al. 2016; Che & Matar 2017). At high impact speeds, the latter situation results in the formation of an upward Worthington jet (Gekle & Gordillo 2010), whereas the former case gives rise to the well-known milkdrop coronet, which is the subject of this paper (Thoroddsen 2002; Deegan et al. 2007; Josserand et al. 2016). Instants after impact on to a thin liquid layer, inertia prompts the formation of an ejecta sheet. After some time, its edge forms a capillary-driven cylindrical-shaped rim that evolves through various stages: first, a combination of Rayleigh–Taylor and Rayleigh–Plateau instabilities grow on the rim (see Zhang et al. 2010; Agbaglah, Josserand & Zaleski 2013), then nonlinearities develop, leading to breakup (through end pinching, see Wang & Bourouiba 2018), and the formation of a cascade of droplet sizes.

Under most natural or industrial conditions, streams are contaminated with surfactants (i.e. surface-active agents), whose concentration variations lead to surface tension gradients, which in turn result in the formation of Marangoni stresses and surface viscosities (see Manikantan & Squires 2020). Previous studies have demonstrated the role of surfactants in capillary singularities (see for instance, Ananthakrishnan & Yeung 1994; Craster, Matar & Papageorgiou 2002; Liao et al. 2004; Kamat et al. 2018; Constante-Amores et al. 2020); here, we study the role of surfactant-induced flows in the late stages of splashing. Che & Matar (2017) demonstrated the role of surfactants on the impact dynamics onto thin liquid films for low and moderate Weber and high Reynolds numbers, e.g. $2.6< We<532$ and $3000< Re<8654$, respectively ($Re=\rho _lUD/\mu _l$ and ${We=\rho _l U^2 D/\sigma }$, where $D$, $\rho _l$, $\mu _l$ $\sigma$ and $U$ stand for the droplet diameter, and its characteristic density, viscosity, surface tension and velocity, respectively). In this past work, the authors showed that the presence of surfactants affects the post-impact dynamics such as the propagation of capillary waves, the growth of the crown and the generation of secondary droplets. However, Marangoni stresses, surfactant concentration dynamics and other important insights into the impact phenomena were not investigated.

The vast majority of past works have focused on the physical understanding of the early dynamics of the ejected sheet using axisymmetric simulations, which is a valid assumption in the limiting case of $t\ll D/U$ (where $t$ stands for time) – see for example Josserand & Zaleski (2003), Josserand et al. (2016) and Agbaglah & Deegan (2014). In contrast, the study of the crown-splash regime requires the use of full three-dimensional (3-D) numerical simulations in order to accommodate the natural occurrence of possible symmetry-breaking events, such as the growth of transverse instabilities in the rim, or rim breakup into droplets. However, 3-D simulations remain a numerical challenge (Gueyffier & Zaleski 1998) and have not often been reported in the literature (Liang et al. 2019; Reijers et al. 2019; Chen et al. 2020; Xavier et al. 2020).

The present study aims to unravel the role of surfactant-induced Marangoni stresses on the interfacial dynamics during the impact of a surfactant-free droplet on a surfactant-laden thin layer of the same liquid in the crown-splash regime. Here, we perform 3-D calculations of splashing in the presence of surfactants by using a hybrid front-tracking method to track the interface location. Section 2 presents the governing dimensionless parameters, the problem geometry, initial and boundary conditions and the numerical validation. Section 3 provides a discussion of the results and concluding remarks are given in § 4.

2. Problem formulation and numerical method

High-resolution simulations were performed by solving the two-phase incompressible Navier–Stokes equations with surface tension in a 3-D Cartesian domain $\boldsymbol {x} = (x, y, z)$ (see figure 1a). The interface between the gas and liquid is described by a hybrid front-tracking/level-set method, where (insoluble) surfactant transport is resolved at the interface (Shin et al. 2018). Here, and in what follows, all variables are made dimensionless (represented by tildes) using

(2.1a–f)\begin{equation} \tilde{\boldsymbol{x}}=\frac{\boldsymbol{x}}{D_o}, \quad \tilde{t}=\frac{t}{t_{r}}, \quad \tilde{\boldsymbol{u}}=\frac{\boldsymbol{u}} {U}, \quad \tilde{p}=\frac{p}{\rho_{{l}} U^2}, \quad \tilde{\sigma}=\frac{\sigma}{\sigma_s}, \quad \tilde{\varGamma}=\frac{\varGamma}{\varGamma_\infty}, \end{equation}

where $t$, $\boldsymbol {u}$ and $p$ stand for time, velocity and pressure, respectively. The physical parameters correspond to the liquid density $\rho _l$, liquid viscosity, $\mu _l$, surface tension, $\sigma$, surfactant-free surface tension, $\sigma _s$, $U$ is the drop impact velocity and $D_o$ its initial diameter; hence, the characteristic time scale is $t_{r}= D_o/U$. The interfacial surfactant concentration, $\varGamma$, is scaled by the saturation interfacial concentration, $\varGamma _{\infty }$. As a result of this scaling, the dimensionless equations read

(2.2)\begin{gather} \boldsymbol{\nabla} \boldsymbol{{\cdot}} \tilde{\boldsymbol{u}}=0, \end{gather}

(2.3)\begin{gather} \tilde{\rho} \left(\frac{\partial \tilde{\boldsymbol{u}}}{\partial \tilde{t}}+\tilde{\boldsymbol{u}} \boldsymbol{{\cdot}}\boldsymbol{\nabla} \tilde{\boldsymbol{u}}\right) + \boldsymbol{\nabla} \tilde{p} ={-}\frac{\mathbf{i}_z }{Fr^2}~ + \frac{1}{Re} \boldsymbol{\nabla}\boldsymbol{{\cdot}} [ \tilde{\mu} (\boldsymbol{\nabla} \tilde{\boldsymbol{u}} +\boldsymbol{\nabla} \tilde{\boldsymbol{u}}^T) ] \nonumber\\ \hspace{163pt}+\,\frac{1}{\rm We} \int_{\tilde{A}\tilde{(t)}} (\tilde{\sigma} \tilde{\kappa} \boldsymbol{n} + \nabla_s \tilde{\sigma}) \delta (\tilde{\boldsymbol{x}} - \tilde{\boldsymbol{x}}_{f} )\,{\rm d}\tilde{A}, \end{gather}

(2.4)\begin{gather}\frac{\partial \tilde{\varGamma}}{\partial \tilde{t}}+\nabla_s \boldsymbol{{\cdot}} (\tilde{\varGamma}\tilde{\boldsymbol{u}}_{t})=\frac{1}{Pe_s} \nabla^2_s \tilde{\varGamma}, \end{gather}

where the density and viscosity are given by $\tilde {\rho }=\rho _g/\rho _{l} + (1 -\rho _g/\rho _{l}) \mathcal {H}(\tilde {\boldsymbol {x}},\tilde {t})$ and $\tilde {\mu }=\mu _g/\mu _{l}+ (1 -\mu _g/\mu _{l}) \mathcal {H}( \tilde {\boldsymbol {x}},\tilde {t})$ wherein $\mathcal {H}( \tilde {\boldsymbol {x}},\tilde {t})$ represents a Heaviside function, which is zero in the gas phase and unity in the liquid phase, while the subscripts ‘$g$’ and ‘$l$’ designate the gas and liquid phases, respectively, and $\tilde {\boldsymbol {u}}_{t}= ( \tilde {\boldsymbol {u}}_{s} \boldsymbol {{\cdot }} \boldsymbol {t}) \boldsymbol {t}$ is the tangential velocity at the interface in which $\tilde {\boldsymbol {u}}_{s}$ represents the interfacial velocity and $\kappa$ the curvature. The interfacial gradient is given by $\nabla _s=({\boldsymbol {I}}-\boldsymbol {n}\boldsymbol {n})\boldsymbol {{\cdot }} \boldsymbol {\nabla }$ wherein $\boldsymbol {I}$ is the identity tensor and $\boldsymbol {n}$ is the outward-pointing unit normal. In addition, $\delta$ is a Dirac delta function, equal to unity at the interface and zero otherwise, and $\tilde {A} (\tilde {t})$ is the time-dependent interface area. As justified by the final paragraph of Stone (1990), in our frame of reference $\boldsymbol {u} \boldsymbol {{\cdot }} \boldsymbol {n}=0$, which gives rise to 2.4. The dimensionless groups that appear in the governing equations are defined as

(2.5)\begin{equation} \left.\begin{array}{c} Re =\dfrac{\rho_l U D_o}{\mu_l}, \quad We =\dfrac{\rho_l U^2 D_o}{\sigma_s}, \quad Fr =\dfrac{ U }{\sqrt{g D_o}}, \\ Pe_s=\dfrac{ U D_o}{\mathcal{D}_s},\quad \beta_s= \dfrac{\Re \mathcal{T} \varGamma_\infty}{\sigma_s}, \end{array}\right\} \end{equation}

where $Re$, $We$, $Fr$ and $Pe_s$ stand for the Reynolds, Weber, Froude and (interfacial) Péclet numbers. Gravity is negligible during the impact as indicated by the Froude number value, i.e. $Fr \sim {O}(10^2)$ (similar assumptions were made by Deegan et al. 2007). The parameter $\beta _s$ is the surfactant elasticity number that is a measure of the sensitivity of the surface tension, $\sigma$, to the interface surfactant concentration, $\varGamma$. Here, $\Re$ is the ideal gas constant value ($\Re = 8.314$ J K$^{-1}$ mol$^{-1}$), $\mathcal {T}$ denotes temperature and $\mathcal {D}_s$ stands for the interfacial diffusion coefficient.

Figure 1. Schematic representation of the flow configuration, and validation of the numerical procedure: (a) initial configuration, highlighting the computational domain of size (not to scale) $(8D_o \times 8D_o \times 4D_o)$ in a 3-D Cartesian domain, $\boldsymbol {x} = (x, y, z)$, with a resolution of $768^2\times 384$; (b) validation of the numerical results (solid lines) against the experimental data (dashed lines) for the crown rim diameter temporal evolution of a surfactant-free case reported in Che & Matar (2017): effect of varying the dimensionless film depth, $h=H/D_o$, while keeping the other parameters constant ($Re=7514$, $We = 249$ and $Fr = 13.65$).

The nonlinear Langmuir equation is used to describe $\sigma$ in terms of $\varGamma$, this is

(2.6)\begin{equation} \tilde{\sigma}=1 + \beta_s \ln{ (1 -\tilde{\varGamma})}. \end{equation}

The Marangoni stress, $\tilde {\tau }$, is expressed as a function of $\tilde {\varGamma }$ as

(2.7)\begin{equation} \frac{1}{\rm We} \nabla_s \tilde{\sigma}\boldsymbol{{\cdot}} {\boldsymbol{t}} \equiv \frac{\tilde{\tau}}{\rm We} ={-} \frac{\rm Ma}{(1-\tilde{\varGamma})} \nabla_s\tilde{\varGamma}\boldsymbol{{\cdot}} {\boldsymbol{t}}, \end{equation}

where $Ma=\beta _s/We=Re \mathcal {T} \varGamma _\infty /\rho _l U^2 D_o$ is the Marangoni parameter and $\boldsymbol {t}$ is the unit tangent to the interface. Finally, a dimensionless thickness of the liquid layer, $h=H/D_o$ is defined as the ratio between the liquid film thickness and the droplet diameter. Tildes are dropped henceforth.

2.1. Numerical set-up, validation and parameters

Figure 1(a) highlights the geometry and the initial and boundary conditions of the problem, which closely follow previous work by Josserand et al. (2016), Agbaglah & Deegan (2014) and Fudge et al. (2021). A drop of initial size $D_o$ with velocity $U$ impacts a uniform layer of thickness $H$ of the same liquid. The size of the computational domain is $8D_0 \times 8D_0 \times 4D_0$, which is sufficiently large to avoid any effects of artificial reflections from the boundaries. The centre of the drop is located at a small distance above the pool surface (e.g. an initial separation of $0.05D_0$). A no-slip boundary and no-penetration conditions are assumed for the bottom wall of the domain, whereas a no-penetration boundary condition is prescribed for the top and lateral domain boundaries (similar to previous work by Batchvarov et al. 2021).

Figure 1(b) illustrates the ability of the numerical technique to reproduce the temporal evolution of the crown experimentally obtained and reported by Che & Matar (2017). After impact, inertia drives the rapid ejection of a vertical sheet from the pool, then capillarity prompts the formation of a rim and then a crown. We observe excellent agreement for $h=0.22$, while a small offset is observed for the other two film thicknesses. The discrepancy could be attributed to the experimental conditions not considered here, such as the surface waves produced during drop formation or air-induced flows during droplet fall. Nonetheless, the agreement between the numerical method and experimental results is visible. Additional validation cases are found in Appendix A. The numerical technique has also been validated for drop–interface coalescence, which can be considered as a limiting case of drop impact onto a pool, i.e. $U=0$ (see Constante-Amores et al. 2021a for more details). The nonlinear interfacial dynamics has been validated for the capillary breakup of liquid threads (see Constante-Amores et al. 2020, 2021b for more details). The numerical validation of the surfactant equations has been previously presented in Shin et al. (2018).

In terms of mesh resolution, we have ensured that our numerical simulations are mesh independent, and as a consequence, the numerical findings do not change by decreasing the cell size for a resolution of $768^2 \times 384$ (i.e. $D_0/\Delta \boldsymbol {x}=96$ cells, see Appendix B). Additionally, liquid volume and surfactant mass conservation are met under errors of less than $10^{-1}\%$ and $10^{-2}\%$, respectively. Extensive mesh studies for surface-tension-driven phenomena using the same computational method can be found elsewhere, e.g. Batchvarov et al. (2021, 2020) and Constante-Amores et al. (2020, 2021b, 2022, 2023).

The dimensionless values for the investigated phenomenon are consistent with experimentally realisable systems. We assume a water–air system, where the density and viscosity ratios are $\rho _g / \rho _l = 1.2 \times 10^{-3}$ and $\mu _g / \mu _l = 0.018$, respectively. The targeted flow conditions for the surfactant-free base case are based on Deegan et al. (2007), who built a phase diagram in the $We-Re$ space for $h=0.2$. Our study focuses on the ‘crown-splash’ regime (i.e. $We>500$ and $Re <1100$), as we aim to stay away from the ‘microdroplet-splash’ regime (i.e. $We>500$ and $Re <1000$). The latter regime is characterised by the formation of a rapid ejecta sheet which undergoes capillary singularities to form microdroplets, requiring an extremely large resolution to capture droplet formation from the ejecta sheet and from the crown. Thus, we have selected $Re=1000$, $We=800$ and $h=0.2$ as the targeted surfactant-free case.

We have examined the case of adding surfactant to the thin liquid film but not the impacting droplet as it is less challenging to control the characteristics of surfactant-free drops experimentally. We have considered insoluble surfactants whose critical micelle concentration (CMC), i.e. $\varGamma _\infty \sim {O}(10^{-6})$ mol m$^{-2}$ for NBD- PC (1-palmitoyl-2-12-[(7-nitro-2-1,3-benzoxadiazol-4-yl)amino]dodecanoyl-sn-glycero-3-phosphocholine); thus, we have explored the range of $0.1 <\beta _s<0.9$ which corresponds to ${O}(10^{-7})<\textrm {CMC}<{O}(10^{-6})$ mol m$^{-2}$, for typical values of $\sigma _s$. We have set $Pe_s=10^2$ following Batchvarov et al. (2020) and Constante-Amores et al. (2020) who demonstrated that the interfacial dynamics is weakly dependent on $Pe_s$ beyond this value.

For a water droplet of a typical size of $D_o \sim {O}(10^{-3})$ m, the impact time scale $T_{imp}=t_r=D_o/U \sim {O}(10^{-4}\unicode{x2013}10^{-3})$ s; whereas the Marangoni time scale $T_{\tau }$ can be estimated by a balance between the Marangoni and the viscous stresses, resulting in $T_{\tau } \sim \mu _l D^2 / (h \Delta \sigma ) \sim {O}(10^{-4}\unicode{x2013}10^{-3})$ s (see Che & Matar 2017 for more details). Therefore, surfactant-induced Marangoni stresses will play a crucial role in the flow physics. A discussion of the results is presented next.

3. Results

First, we discuss the early-time dynamics of the surfactant-free impact at $Re=1000$ and $We=800$. Figure 2(a) shows the 3-D representation of the interface at early times, i.e. to show the formation of transverse instabilities growing at the rim. As seen in this figure, instants after impact, we observe the formation of an axisymmetric thin layer of ejecta fluid which draws fluid from the droplet eventually giving rise to a Peregrine sheet, in agreement with Deegan et al. (2007) and Josserand et al. (2016). As time evolves, the Peregrine sheet grows upwards and expands radially while its film thickness reduces; here, the radial coordinate $r$ is defined at the initial point of impact. A closer inspection of the interfacial shape of the ejecta sheet at $t=0.5$ demonstrates that its orientation is nearly horizontal, i.e. parallel to the pool. We also see surface tension forces inducing the formation of a hemispherical tip at the edge of the ejecta, which leads to a pressure gradient between the tip and the adjacent sheet, which results in fluid motion from the tip towards the sheet. This behaviour leads to the formation of a sheet perpendicular to the pool at longer times; this is in agreement with the experimental observations by Zhang et al. (2010) and Deegan et al. (2007). Figure 2(b) shows the temporal evolution of the rim size for the surfactant-free case at the crossing of the $x$–$z$ plane ($y=4.90$). The peak observed in the range $t=3\unicode{x2013}5$ is due to the growth of corrugations in the rim. As the rim elongates to form the crown, we observe the growth of an azimuthal instability on the surface. Our results indicate that these naturally occurring rim instabilities select wavelengths that are consistent with the most unstable Rayleigh–Plateau (RP) instability. We have measured the average distance ($\lambda$) between the ‘undulations’ or crests seen along the rim. We then normalised the result by the most unstable RP instability wavelength, given by $\lambda ^*= 2 {\rm \pi}/k= 9.016 a$, where $k$ and $a$ are the wavenumber of the fastest-growing instability mode and the rim radius, respectively (Zhang et al. 2010; Agbaglah et al. 2013). In figure 2(c) we plot the selected wavelength as a function of the rim radius; for the cases studied here, $\lambda /\lambda ^* = 1.0 \pm 0.1$. These results are in agreement with the experimental results from Zhang et al. (2010) and the theoretical and computational work from Agbaglah et al. (2013) and Constante-Amores et al. (2022). A closer inspection of the top image of figure 2(c) shows that the interfacial corrugations are not always evenly distributed along the rim, as observed in experiments (Thoraval et al. 2012). In particular, we observe large wavelengths, which are the result of interactions between adjacent instability modes, which is consistent with experimental observations.

Figure 2. Wavelength selection in the crown-splash regime. Panel (a) presents the early interfacial dynamics through a 3-D representation of the predicted interface location for $Re=1000$ and $We=800$, at $t=(0.5, 1.0, 1.5)$ corresponding to columns one to three, respectively. Panel (b) plots the temporal evolution of the rim radius. Panel (c) shows the selection of the wavelength $\lambda$ normalised with the theoretical RP instability, $\lambda ^*$ as a function of the local rim radius, $a$, at early times of the simulation for $t<5$. A 3-D representation of the crown for the dimensionless time $t=7.5$ is also presented, highlighting the predicted wavelength between crests.

Figure 3(a–d) shows a 3-D representation of the temporal evolution of the interface, from the time instabilities become visible to the time droplets breakup by end pinching. At these longer times ($t>5$), the rim is weakly affected by the dynamics occurring at the impacting point, and interfacial RP instabilities trigger the development of ligaments (see figure 3a). With increasing time, liquid ligaments grow perpendicular to the rim, resulting in a local pressure-gradient-induced flow from the ligament to the tip, triggering the formation of a capillary-induced blob and a neck (displayed in figure 3b). The neck thins and stretches, due to a capillary-driven flow, eventually resulting in the capillary pinch-off of a drop (the well-known ‘end-pinching’ mechanism proposed by Stone & Leal 1989), see figure 3(c). We note that the interfacial dynamics closely resembles experimental observations by Zhang et al. (2010), Josserand & Zaleski (2003) and Deegan et al. (2007), as well as the dynamics of retracting sheets reported by Wang & Bourouiba (2018) and Constante-Amores et al. (2022).

Figure 3. Effect of $\beta _s$ on the drop impact dynamics for insoluble surfactants. Spatio-temporal evolution of the 3-D interface shape for surfactant-free cases, (a–d), and surfactant-laden cases for $\beta _s = (0.1, 0.5, 0.7)$ corresponding to panels (e–h), (i–l) and (m–p), respectively. Here, the dimensionless parameters are $Re=1000$, $We=800$ and $h=0.2$. For the surfactant-laden cases, $Pe_s=100$ and $\varGamma =\varGamma _\infty /2$, the colour indicates the value of $\varGamma$ and the legend is shown in panel (e).

Next, we proceed to study the role of surfactants in the flow dynamics; as illustrated by figure 3(e–p) where we observe the droplet impact at various elasticity parameters $\beta _s$. As seen in this figure, regardless of the value of $\beta _s$, large surfactant concentration gradients are found on the outer interface of the ejecta sheet. This result supports the hypothesis that the ejecta sheet is generated at the interface of the surfactant-laden pool. On the other hand, the inner interface of the ejecta sheet arises from the surfactant-free drop, i.e. $\varGamma$ vanishes on the inner surface, see figure 3(e). This observation, that the fluid sources of the inner and outer surfaces of the ejecta sheet are the pool and the droplet, respectively, is in good agreement with the work by Josserand et al. (2016). A close inspection of the distribution of $\varGamma$ in the outer sheet shows that $\varGamma$-values are maximum at the base of the ejecta, and their values reduce upwards (see for example panel (e) of figure 3). This is the result of an increase of the interfacial area, that leads to a dilution of surfactant at the interface, thus the surfactant concentration (surface tension) decreases (increases). In addition, the interfacial expansion results in large convective effects that transport surfactant towards the rim; see figure 3(e), in which larger values of $\varGamma$ are found in the rim than in its adjacent sheet.

The increase (decrease) of $\varGamma$ ($\sigma$) implies that a large surface deformation is required to satisfy the normal stress balance at the interface. Similar to the surfactant-free case, capillary-induced flow results in the formation of a rim on the sheet edge, whose thickening leads to the onset of destabilisation. This is in agreement with previous studies from Asaki, Thiessen & Marston (1995), who reported a surfactant-driven dampening effect on the capillary waves.

In addition, we expect higher (lower) values of $\varGamma$ at the rim wave crests (instability waves) as they belong to a radially converging (diverging) region, which has the effect of locally increasing (reducing) $\varGamma$. The gradients in $\varGamma$ result in surfactant-induced Marangoni stress flows from the lower-$\sigma$ radially converging to the higher-$\sigma$ radially diverging regions. This slows down the interfacial dynamics, and the formation of instabilities at the rim, in agreement with the idealised case presented by Constante-Amores et al. (2022). By increasing $\beta _s$ the surfactant distribution along the interface is enhanced, as displayed in figures 3(i–l) and 3(j–m) for $\beta _s=0.5$ and $\beta _s=0.7$, respectively.

Surfactant-induced Marangoni stresses only delay the development of the ligament from the rim; these eventually break up via the end-pinching mechanism to form droplets, as shown in figure 3(g,k,o). These results offer an explanation for the experimental observations of Che & Matar (2017). Figure 4 displays the selected instability wavelength in terms of $\beta _s$ and $a_0$. Indeed, a close inspection indicates that the surfactant does not seem to greatly affect the selection of the wavelength; the predicted $\lambda$ is consistent with the most unstable RP instability at a ratio of $\lambda /\lambda ^* = 1.0 \pm 0.1$. Consequently, according to our results, the crown dynamics is mostly driven by the RP instability, although longer ligaments developed prior to break up in the presence of surfactants, as observed in figure 3(p). This indicates that Marangoni stresses promote the retardation from end pinching, which, under certain circumstances may even lead to the neck re-opening, further delaying the splash, a phenomenon previously demonstrated by Kamat et al. (2020) and Constante-Amores et al. (2020, 2022).

Figure 4. Effect of the elasticity parameter, $\beta _s$, on the selection of wavelength of the undulations normalised with the theoretical RP instability, $\lambda ^*$, as a function of the local rim radius, $a$.

Two-dimensional projections, in the $x$–$z$ plane ($y=4)$, of the interfacial shape, $\varGamma$, $\tau$ and the radial component of the interfacial velocity $u_{tz}$ in terms of $\beta _s$ and a fixed value of $Pe_s$, are displayed in figures 5 and 6 corresponding to $t=5$ and $t=20$, respectively. As described above, at early times, i.e. $t=5$, the drop impact results in a non-uniform interfacial surfactant distribution, with high surface concentrations at the base of the ejected sheet. As seen, the interfaces of the surfactant-free and the $\beta _s=0.1$ surfactant-laden cases are practically indistinguishable. However, as the surfactant concentration increases, $\varGamma$ distributions trigger a surfactant-driven dynamics from the sheet base to its edge (see figure 5b). As $\beta _s$ increases, long ejecta sheets are promoted, resulting in the reduction of the rim and sheet thickness due to mass conservation (displayed in the magnified panel of figure 5a). The increase in $\beta _s$ also enhances the $\varGamma$ distribution and increases the value of $\tau$ (see the increase of $\tau$ with $\beta _s$ in figure 5c). Furthermore, figure 5(c) highlights the presence of large $\tau$ peaks on both sides of the rim due to the accumulation of $\varGamma$, and sudden drops at the base of the rim, suggesting a fluid transport from the rim to the sheet. In addition, a larger $\tau$ maximum is observed from the side of the inner sheet as the inner sheet is devoid of $\varGamma$, leading to an uneven flow motion inside the sheet. The variation of the streamwise interfacial tangential velocity $u_{t_{z}}$, along the arc length, is shown in figure 5(d). This figure, at the early stages of the dynamics, exhibits $u_t>0$ over the majority of the domain due to the dominance of the vertical radial expansion of the ejecta sheet. A strong variation of $u_{t_{z}}$ is observed around the rim with $u_{t_{z}}$ tending to zero at a sufficiently large distance away from the impactregion.

Figure 5. Effect of the elasticity parameter, $\beta _s$, on the flow dynamics at $t=5$. Two-dimensional projections of the interface, $\varGamma$, $\tau$ and $u_{tz}$ in the $x$–$z$ plane ($y=4$) are shown in (a–d), respectively. In panel (a), a magnified view of the ejecta sheet is also presented. Note that the abscissa in (a) corresponds to the $x$ coordinate, and in (b–d) to the arc length, $s$. The arc length $s$ corresponds to the $x$–$z$ plane ($y = 4$) intersecting the interface; $s$ has been normalised on the full extent of $s$ associated with the length of the impact region in each case. The diamond shapes in panels (a,b) indicate the location of the crown. All parameters remain unchanged from figure 3. Surf., surfactant.

Figure 6. Effect of the elasticity parameter, $\beta _s$, on the flow dynamics at $t=20$. Two-dimensional projections of the interface, $\varGamma$, $\tau$ and $u_{tz}$ in the $x$–$z$ plane ($y=4$) are shown in (a–d), respectively. Note that the abscissa in (a) corresponds to the $x$ coordinate, and in (b–d) to the arc length, $s$. The arc length $s$ corresponds to the $x$–$z$ plane ($y = 4$) intersecting the interface; $s$ has been normalised on the full extent of $s$ associated with the length of the impact region in each case. The diamond shapes in panels (a,b) indicate the location of the crown. All parameters remain unchanged from figure 3.

The influence of $\beta _s$ at long times, $t=20$, is examined in figure 6. Qualitatively, the surfactant effects are particularly evident as the dynamics has been slowed down by the surface rigidification brought about by the Marangoni stresses, as can be seen in figure 6(d) (the overall magnitude of $u_{t_z}$ is smaller with increasing $\beta _s$). As $\beta _s$ decreases ($\sigma$ increases), the surfactant is highly convected as $\tau$ tends to zero, lowering the interfacial tension and consequently increasing surface deformation. As $\beta _s$ increases, the gradient in $\varGamma$ ($\tau$) decreases (increases) near the rim, and the Marangoni stresses predominantly retard the flow and ‘rigidify’ the interface.

As can be observed in the results presented in figure 6(c), an increase in $\beta _s$ results in the strengthening of $\tau$. Figure 6(d) provides conclusive evidence of the surface rigidification as most of $u_{tz}$ is characterised by a negative value due to the decay of the dynamics. The variation of $u_{tz}$ decreases as $\beta _s$ increases, strengthening the Marangoni stresses arising from the surfactant-induced surface tension gradients. Finally, according to the foregoing results, we conclude that the interfacial dynamics for surfactant-laden cases in the ‘crown regime’ resembles its surfactant-free counterpart except for the retardation due to the surface rigidification brought about by the presence of surfactant-induced Marangoni stresses.

In figure 7 we plot the component of the velocity variation in the sheet-normal direction, $u_z$, within the sheet, in a frame of reference moving with the inner sheet at $t=5$. For the surfactant-free case, the liquid within the sheet follows a parabolic profile into the rim; this effectively corresponds to a no-slip condition with $u_z$ pinned to the interface, due to the absence of surfactants. In contrast, the retardation brought about by the presence of surfactants via Marangoni stresses results in shear flow for the surfactant-laden case. The higher surfactant concentration in the outer rim wall results in $\tau$ triggering the deceleration of the fluid entering into the bulbous end near the outer wall.

Figure 7. Profiles of the velocity component, $u_z$, in the sheet-normal direction inside the sheet for the surfactant-free and surfactant-laden $(\beta _s=0.5)$ cases at $t=5$. Here, $u_z$ has been calculated in an inner-sheet frame of reference along the cross-stream direction, $x$, for the axial location along the rim neck. Note that the axial distance has been normalised with the average sheet thickness, i.e. $x/h_s$ (thus $x/h_s=0$ and $x/h_s=1$ correspond to the inner and outer sheets, respectively). All parameters remain unchanged from figure 3.

Figure 8(a) presents the late stages of a surfactant-laden ligament before the end-pinching mechanism. In the absence of surfactants, the competition between viscosity and capillarity drives the formation of and dynamics of the bulbous edge on the ligament tip. We proceed to calculate the relative importance of viscous to surface tension forces via the local Ohnesorge number for the ligament (i.e. $Oh=\mu _l / \sqrt {\rho _l \sigma _s R}$, where $R$ stands for the average radius of the ligament) to explain the retardation phenomenon observed in figure 8(a). The resulting local $Oh \sim 0.07$ (e.g. low $Oh$-regime), and subsequently the local flow dynamics affecting the ligament (i.e. surfactant-driven retardation from end pinching) can be explained by the mechanism proposed by Constante-Amores et al. (2020), Kamat et al. (2020) and Hoepffner & Paré (2013). Figures 8(b,c) and 8(d,e) show $\varGamma$ and $\tau$, and the tangential velocity components $u_t$ for the framed panels in figure 8(a), respectively. Figure 8(b,d) demonstrates that, near the neck (marked with a blue diamond symbol) $\tau$ results in $u_{tx}>0$, preventing the formation of a second stagnation point required to trigger the capillary singularity. As shown by Constante-Amores et al. (2020), two stagnation points sandwiched between the neck boundaries are needed to induce end pinching, as can be observed in 8(d,e) (‘SP’ stands for stagnation points). Surface tangential stresses prevent capillary break-up for a longer time than in the surfactant-free case, leading to longer ligaments and an increase in surface area. This results in the dilution (increase) of surfactant (surface tension) and the eventual reduction in the strength of surfactant-induced Marangoni stresses thus initiating end pinching.

Figure 8. (a) Spatio-temporal evolution of a surfactant-laden ligament and its retardation from pinch-off driven by surfactant-induced Marangoni flow when $\beta _s=0.5$. (b,c) Represent $\varGamma$ and $\tau$ as a function of the arc length $s$ for the framed panels of (a), i.e. $t=13.49$ and $t=15.49$, respectively; (d,e) represent the tangential velocity $u_t$ as a function of the arc length $s$. The axial location in (a) has been normalised over the distance of the ligament on the axial direction; the abscissa direction has been normalised with respect to its value at $x=0$ for each single panel. The arc lengths have been normalised over the full extent of $s$. The diamond markers represent the location of the neck; ‘SP’ in (e) indicates the location of the stagnation points. All parameters remain unchanged from figure 3.

Figure 9 presents the temporal evolution of the neck radius for the surfactant-free and surfactant-laden cases (we cannot report the late stages of the neck due to the lack of snapshots up to its capillary singularity). As seen, at short times, the presence of surfactants at the interface leads to the reopening of the neck ($\textrm {d}r/\textrm {d}t \sim 0$), and subsequently, to a short period in which surfactants induce retardation from end pinching. At a later time, as surfactants are evacuated from the neck, the dynamics eventually results in the interfacial singularity.

Figure 9. Temporal evolution of the neck radius for the surfactant-free and surfactant-laden cases ($\beta _s=0.5$), highlighting the surfactant-driven retardation of the interfacial capillarity.

Next, we focus on drop formation and the splash phenomena: we have identified three different modes of droplet shedding. The first mode corresponds to droplet detaching from the ligament via end pinching. The mean lengths of the ligament $\langle l\rangle$ which undergoes mode-1 of ejection are: $\langle l\rangle =0.775$, $\langle l\rangle =0.864$, $\langle l\rangle =0.992$ and $\langle l\rangle =1.141$, which correspond to the surfactant-free and surfactant-laden cases with $\beta _s=(0.1,0.5,0.7)$, respectively. This agrees with a surfactant-driven retardation of the capillary singularity. In the second mode we observe the initial ejection of filaments that are rapidly recombined into a single filament, see figure 11(a–e). The third mode sees the formation of satellite (smaller) droplets during the stretching of the filament neck. These satellite droplets are not observed for large values of the elasticity parameter, in agreement with Kovalchuk et al. (2018). For the surfactant-free case, figure 10, ligaments merge due to the presence of local cusps arising from a non-uniform distribution of mass per unit length (in agreement with Gordillo, Lhuissier & Villermaux (2014) and Wang & Bourouiba (2018)), or by the rim parameter not being a multiple of the RP instability. In the vicinity of a cusp, the rim is not perpendicular to the flow entering the rim; instead, it takes an angle $\theta$, as shown in figure 10(a). This, oblique flow induces a local drift velocity along the rim. In a reference frame along the rim, the drift velocity can be calculated from the projection of the incoming velocity $u_{in}$ in the longitudinal direction of the rim as $u_{drift}=u_{in}(t) \sin (\theta )$ (see Wang & Bourouiba 2018 for more details). Under the conditions of figure 10, we obtain $u_{drift} \sim 0.106$ m s$^{-1}$ and $\theta \sim 20^{\circ } \pm 3^{\circ }$. In contrast, from simulations, the velocity field yields $u_{in} \sim 0.27$ m s$^{-1}$ which leads to $u_{drift} \sim 0.083$ m s$^{-1}$, and a good agreement with the direct measurement. The collision of the drifting protrusions results in the formation of a corrugated ligament (see figure 10e), which succumbs to end pinching, resulting in drop formation (not shown). We confirm that we have made the same calculation for the other two flow-induced cuspids which give in similar results.

Figure 10. Spatio-temporal evolution of two adjacent rim protrusions which shift in the spanwise direction induced by local cusps, resulting in their collision; $\theta$ denotes the angle made by the left crest with the rim. The time difference between snapshots is $\Delta t=1.0$. The parameters remain unchanged from figure 2.

The presence of surfactant enhances the ligament-merging mode due to surfactant-driven Marangoni stresses that induce a local motion along the rim, and a radial shift of the ligament position along the rim. This is highlighted in figure 11, which presents the temporal evolution of two protrusions that eventually become ligaments, and merge prior to pinch-off.

Figure 11. Spatio-temporal evolution of two adjacent rim protrusions shown panels (a–e) which shift in the spanwise direction due to surfactant-induced Marangoni stresses, resulting in their collision and merging. Panel (f) shows a 3-D reconstruction of the $|\tau |$ profile with respect to the arc length, $s$, across a plane cutting the rim in (a) in half; the arrows represent the direction of $|\tau |$. Here, $\beta _s=0.5$, and all other parameters remain unchanged from figure 3.

Figure 11(a) demonstrates that the highest $\varGamma$ values are found at the protrusions as the fluid motion along the rim to the protrusions acts to accumulate surfactant locally. They are characterised by regions of a radially converging dynamics, and subsequently, increased (reduced) $\varGamma$ ($\sigma$) locally. These results demonstrate that the effect of Marangoni stresses from adjacent protrusions is to bridge their gap, promoting the collision and merging of the ligaments. Figure 11(f) shows a 3-D reconstruction of the $|\tau |$ profile with respect to the arc length, $s$, across a plane cutting the rim seen in 11(a) in half; the arrows represent the direction of $|\tau |$, evidencing the enhancement of the filament merging. To the best of our knowledge, this surfactant-induced mechanism of droplet shedding has not been reported previously. Figure 11 also shows that surfactants are actively convected from the rim to the film. Consequently, the crown is thinner and survives longer in the presence of surfactants, in agreement with the experimental observations reported by Che & Matar (2017).

The total number of ejected droplets for our surfactant-free case is found to be approximately $30$ up to times $t < 22.5$. In contrast, the number of secondary droplets (at $We,Re=800,1000$) as reported by Agbaglah & Deegan (2014), corresponds to 30–32 so our numerical predictions are in good agreement. This also verifies that a RP instability is the primary mechanism for wavelength selection.

Next, we turn our attention to the size distribution of the droplets ejected from the impact. Figure 12(a) shows the probability density function (PDF), at $t=20$, scaled by the initial droplet diameter. The sampling is conducted after the impact reaches $t=20$ and over three snapshots separated by time intervals of $\Delta t=0.2$. To get a smoother PDF, more snapshots should be included, but the dynamics is too rapid and due to the extreme computational cost of the current numerical simulation, this task was not undertaken. It is important to note that, as $\beta _s$ increases, a higher droplet size is favoured, as seen in figures 12(b–d) and 12(e–g), where the droplet volume, average surfactant concentration and average velocity for each droplet at $t=15$ and $t=20$, respectively, are shown. The average velocity and surfactant concentration are calculated as $\langle |\boldsymbol {u}|\rangle =\int _{\varOmega }|\boldsymbol {u}| \,\textrm {d} \varOmega / \varOmega$ and $\langle \varGamma \rangle =\int _{\varOmega }\varGamma \,\textrm {d}\varOmega / \varOmega$, respectively, where $\varOmega$ represents the surface area of the droplet. We observe that fewer droplets are produced with increasing $\beta _s$; fewer droplets are produced due to the rigidification brought about by surfactant-driven Marangoni stresses (see figure 12b,e). For the high values of $\beta _s$, we observe some droplets with large volume corresponding to the merging of ligaments. We also observe that small droplets are produced for the surfactant-free and the $\beta _s = 0.1$ surfactant-laden cases, (see figure 12e). These small droplets correspond to satellite droplets, which are characterised by $V, \langle |\boldsymbol {u}|\rangle$ and $\langle \varGamma \rangle \xrightarrow {} 0$, in good agreement with Kovalchuk & Simmons (2018) and Che & Matar (2017). We conclude that the addition of surfactants leads to a larger droplet size distribution for two reasons: the retardation of the dynamics, and the suppression of end pinching through the reopening of the neck driven by a surfactant-induced flow.

Figure 12. Metrics of splashing as a function of the elasticity parameter $\beta _s$. Panel (a) shows the probability density function (PDF) of the droplet size; droplet sizes are made dimensionless by the initial droplet diameter $D_0$. Panels (b–d) and (e–g) show the volume, average surfactant concentration and average velocity for every droplet at $t=15$ and $t=20$, respectively ($d_i$ stands for the enumeration of droplets). All other parameters remain unchanged from figure 3.

Figure 12(c,f) indicates that the average $\langle \varGamma \rangle$ of the droplets is lower than the initial surfactant concentration of the pool (viz. $\varGamma =0.5$). This agrees well with the experimental studies of Blanchard & Syzdek (1972) and Constante-Amores et al. (2021c) on bursting bubbles. Figure 12(c) highlights that $\langle \varGamma \rangle$ increases with decreasing $\beta _s$ as the surfactant has been highly convected towards the rim. Lastly, figure 12(d,g), which shows $\langle |\boldsymbol {u}|\rangle$ at $t=15$ and $t=20$, respectively, demonstrates the effect of rigidification brought about by surfactant-induced Marangoni flow which corresponds to an overall reduction in $\langle |\boldsymbol {u}|\rangle$ for the surfactant-laden droplets.

Lastly, we plot in figure 13 the effect of surfactants on the interfacial area, and kinetic energy defined as $E_k=\rho \int _V \boldsymbol {u}^2/2 \,\textrm {d}V$, normalised by their initial values. Inspection of the surface area plot in figure 13(a) reveals a linear growth rate in interfacial area for all cases at short times; at longer times, a monotonic increase of surface area is observed with increasing $\beta _s$ due to surfactant-induced effects discussed above. In 13(b), we see that the presence of surfactants by increasing $\beta _s$ acts to decrease the overall value of $E_k$ (e.g. rigidification brought about surfactants).

Figure 13. Temporal evolution of the total interfacial area (a), and the kinetic energy (b), normalised by their initial values, for the surfactant-free and surfactant-laden cases. All parameters remain unchanged from figure 3.

4. Concluding remarks

We have presented, for the first time, a detailed analysis of the effect of insoluble surfactants on the dynamics of droplet impact on thin-film liquid layers using high-resolution 3-D simulations. This article focuses on the crown-splash regime, in which the transverse destabilisation of the rim leads to the formation of ligaments and droplets. Thus, we have focused on conditions of high Reynolds and Weber numbers. The nature of the hybrid front-tracking level-set numerical simulations enables the coupling, and analysis, of inertia, capillarity, interfacial diffusion and Marangoni stresses owing to the presence of surfactant-induced surface tension gradients. According to our results, the addition of surfactants does not significantly affect the selection of the wavelength of the transverse rim instability; as the predicted wavelength is consistent with the most unstable RP instability. The instability give rise to the formation of ligaments which eventually end pinch into droplets. In the presence of surfactants, we observe a delay in end pinching, driven by surfactant-induced Marangoni stresses, resulting in longer ligaments. We identify three modes of droplet ejection, and demonstrate that surfactant-induced Marangoni stresses result in the bridging of the spanwise surfactant concentration gradients between adjacent ligaments in the rim, eventually leading to their merging. The presence of surfactant-induced Marangoni stresses also leads to interfacial rigidification, which is observed through a reduction of the surface velocity and the kinetic energy, and a larger interface area.

This research is of importance to many applications as droplet impact onto finite-depth liquid layers is found in a wide range of industrial and daily-life applications, from raindrops impinging onto puddles to screen and inkjet printing. This research focuses on the crown-splashing regime but we acknowledge that the presence of surfactants will play a crucial role in other $We-Re$ regions, e.g. at the microdroplet-splash regime. Although we have only focused on insoluble surfactants, the presence of surfactant solubility will lead to additional richness and complexity, and we anticipate that the addition of soluble surfactants may result in detrimental effects, and they will be the subject of future work. Additionally, surfactant-induced rheological effects could be important and of great relevance to industry as most industrial fluids are formulated with dispersants (surfactants).

Finally, it is well known that the thickness of the film layer affects the interfacial dynamics. Additionally, most of the literature has focused on the limiting case of $t\ll D/U$ in which axisymmetric simulations suffice. The main novelty of this work is the study of the role of the surfactant at longer times after rim formation. At these times, surfactants play a major role in the late stages of droplet detachment, and, to the best of our knowledge, this manuscript is the first computational study on this phenomenon. Consequently, future studies could focus on the surfactant-induced effects of droplet impacting onto various liquid film thicknesses, and the ejecta sheet.