Role of surfactant-induced Marangoni stresses in retracting liquid sheets

In this work, we study the effect of insoluble surfactants on the three-dimensional rim-driven retraction dynamics of thin water sheets in air. We employ an interface-tracking/level-set method to ensure the full coupling between the surfactant-induced Marangoni-stresses, interfacial diffusion, and inertia. Our findings are contrasted with the (Newtonian) dynamics of a liquid sheet edge finding that the surfactant concentration can prevent, or delay, the breakup of the rim. Our simulations use the fastest growing Rayleigh-Plateau instability to drive droplet detachment from the fluid sheet (rim). The results of this work unravel the significant role of Marangoni stresses in the retracting sheet dynamics at large elasticity numbers. We study the sensitivity of the dynamics to the elasticity number and the rigidification of the interface.


Introduction
The capillary break-up of liquid films into droplets plays a major role in various unsteady fluid fragmentation phenomena, such as atomisation or droplet splashing (Villermaux 2020). It is therefore unsurprising that sheet retraction dynamics have received significant interest since the initial experimental observations by Dupre (1867) and Rayleigh (1879). Capillary film retraction follows the puncturing of a static liquid film, where a capillary-induced flow drives the opening of the interface (a hole) at a constant speed, i.e. the Taylor-Culick velocity (Taylor 1959;Culick 1960). Early pictures by Rayleigh (1891) and Ranz (1959) show that the opening of the sheet results in the formation of a rim with roughly cylindrical caps, as the film moves away from the puncture. The desire to understand the fundamental mechanism underlying the development of the hole-driven expansion has led to numerous studies in this field, see for instance Debrégeas et al. (1995); Keller et al. (1995); Brenner & Gueyffier (1999); Fullana & Zaleski (1999); Sünderhauf et al. (2002); Savva & Bush (2009) ;Roisman et al. (2006); Krechetnikov (2010); Lhuissier & Villermaux (2011); Villermaux & Bossa (2011); Gordillo et al. (2011); Agbaglah et al. (2013); Agbaglah (2021). Debrégeas et al. (1995) and Savva & Bush (2009) concluded that the retracting dynamics of a liquid sheet are governed by a balance between inertial, viscous, and surface-tension forces. Accordingly, the Ohnesorge number, Oh, (e.g., ratio of viscous to capillary forces) becomes the most appropriate parameter to parametrise the flow dynamics. Different regimes have been identified that depend on Oh: (i) in the Oh < 1 regime, which is characterised by a nearly inviscid flow, the dynamics are driven by surface tension leading to the formation of capillary waves ahead of the roughly cylindrical rim; (ii) in the Oh > 10 regime viscous forces dominate the dynamics causing the suppression of the rim; and (iii) an intermediate regime bridges the two previous regimes.
The groundbreaking experiments of McEntee & Mysels (1969) proposed that a surface instability along the rim drives the formation of ligaments, and eventually, the ejection of droplets. Several studies have suggested that different physical mechanisms prompt the detachment of these droplets: a Rayleigh-Plateau (RP) instability (discussed below), a nonlinear amplification mechanism (Yarin & Weiss 1995), a Richtmyer-Meshkov instability, and a competition/collaboration between the Richtmyer-Meshkov and Rayleigh-Taylor (RT) instabilities (Krechetnikov 2010). However, Bremond & Villermaux (2006) (for the fragmentation of the lamella of colliding jets), Rieber & Frohn (1999) and Zhang et al. (2010) (for the fragmentation of the sheet-crown from impacts of drops onto a thin layer of liquid) concluded that the RP instability is responsible for the onset of drop detachment. These findings contrast with the previous study by Fullana & Zaleski (1999), that suggest that the growth of the rim during the retraction prevents the onset of a RP instability. Finally, Agbaglah et al. (2013) (for a numerical sheet retraction in a frame of reference moving at the Taylor-Culick velocity) and  (for drop impact onto a solid surface) showed that both RT and RP instabilities are responsible for droplet generation from the rim. The former is important at early times during the deceleration of the ejected sheet, whereas the latter is predominant at longer times once the rim has already formed. According to recent studies by , 2021, the ligament dynamics can result in end-pinching, break up into satellite droplets or recombination into a single fluid volume.
McEntee & Mysels (1969) used surfactants in their experimental work to study the bursting of liquid films without directly endorsing the effects of Marangoni-induced flows. In contrast, recent studies have demonstrated the crucial role of surfactants on the dynamics of the capillary singularity (Ambravaneswaran et al. 2000;Timmermans & Lister 2002;Craster et al. 2002;Liao et al. 2006;Kamat et al. 2018). Constante-Amores et al. (2020) and Kamat et al. (2020) showed that surfactant-induced Marangoni-stresses inhibit end-pinching for retracting liquid threads via the suppression of stagnation points, which in turn leads to flow reversal near the vicinity of the neck.
Three-dimensional numerical simulations of retracting liquid sheets are scarce in the literature due to the need for large aspect ratios to induce both the growth of instabilities at the rim, and the development of the nonlinear flow dynamics. In this work, the role of surfactant-induced Marangoni stresses on the retraction of thin liquid sheets and the detachment of droplets is studied in a three-dimensional, nonlinear framework. This paper is structured as follows: Section 2 introduces the numerical method, governing dimensionless parameters, problem configuration, and validation; Section 3 provides the results from the simulations, and concluding remarks are given in Section 4.

Problem formulation and numerical method
Numerical simulations were performed by solving the transient two-phase Navier-Stokes equations in a three-dimensional Cartesian domain x = (x, y, z) (see figure 1). The interface is captured via a hybrid front-tracking/level-set method; and a convectivediffusion equation is solved for the surfactant transport along the deforming interface (Shin et al. 2018). All variables are made dimensionless following: where, t, u, and p stand for time, velocity, and pressure, respectively. The physical parameters correspond to the liquid density ρ l , viscosity, µ l , surface tension, σ, surfactantfree surface tension, σ s , and u T C = 2σ s /(ρh 0 ), is the well-known Taylor-Culick speed; here h 0 stands for the sheet thickness. In this study we focus on the low-Oh regime; hence, the characteristic time scale is given by the capillary time, t inv = h 0 /u T C = ρ l h 3 0 /(2σ s ). As seen, the interfacial surfactant concentration, Γ , is also made dimensionless through the saturation interfacial concentration, Γ ∞ . As a result of the scaling in equation (2.1), the dimensionless forms of the governing equations for the flow and the surfactant transport are respectively expressed as which correspond to the equations of mass and momentum conservation and the convective-diffusion equations for the interfacial concentration, respectively. Here, the densityρ and viscosityμ are defined byρ = ρ g /ρ l + (1 − ρ g /ρ l ) H x,t and µ = µ g /µ l + (1 − µ g /µ l ) H x,t wherein H x,t represents a smoothed Heaviside function. In this work, H x,t is zero in the gas phase and unity in the liquid. The subscript g designates the gas phase,ũ t = (ũ s · t) t stands for the velocity vector tangential to the interface in whichũ s corresponds to the interfacial velocity, κ stands for the interfacial curvature obtained from the Lagrangian interface structure, and ∇ s = (I − nn) · ∇ stands for the surface gradient operator wherein I is the identity tensor and n is the outward-pointing unit normal to the interface. Finally,x f is the parameterisation of the interfaceÃ(t), and δ represents a Dirac delta function that is non-zero whenx =x f only. The dimensionless groups that govern the dynamics are defined as where Oh and P e s stand for the Ohnesorge and (interfacial) Peclet numbers, respectively, while β s is the surfactant elasticity number which provides a measure of the sensitivity of σ to Γ ; here, is the ideal gas constant value ( = 8.314 J K −1 mol −1 ), T denotes temperature and D s stands for the diffusion coefficient. The non-linear Langmuir equation expresses σ in terms of Γ , i.e.σ = 1 + β s ln (1 −Γ ), and the Marangoni stress, τ , is given as a function ofΓ asτ = ∇ sσ · t = −β s /(1 −Γ )∇ sΓ · t. Tildes are dropped henceforth.
We present a brief summary of the numerical method used in this study. The Navier Stokes equations are solved by a finite volume method on a staggered grid (Harlow & Welch (1965)). The computational domain is discretised by a fixed regular grid (i.e. Eulerian grid) and the spatial derivatives are approximated by standard centred difference discretisation, except for the non-linear term, which makes use of a secondorder essentially non-oscillatory (ENO) scheme, Sussman et al. (1994). The interface is tracked explicitly by an additional Lagrangian grid by using the Front-Tracking method, and the interface is reconstructed by using the Level Contour Reconstruction Method Shin & Juric (2002. Communication between the Eulerian and Lagrangian grid (for the transfer of the geometric information of the interface) is done by using the discrete delta function and the Immersed Boundary Method of Peskin (1977). The advection of the Lagrangian interface is done by integrating dx f /dt = V with a second-order Runge-Kutta method, where V stands for the interfacial velocity which has been calculated by interpolation from the Eulerian velocity. The code is parallelised through an algebraic domain-decomposition technique and the communication between subdomains for data exchange is managed by the Message Passing Interface (MPI) protocol. More details of the numerical method used to solve the above equations is described in detail by Shin et al. (2018) 2.1. Numerical setup and validation Figure 1a illustrates the geometry used in this study, we have initialised the sheet film dynamics using a linearly unstable configuration as proposed by Agbaglah et al. (2013). Thus, we have considered a thin liquid sheet of thickness h 0 , and initial length L 0 , connected to a cylindrical rim of radius a 0 . The initial shape of the rim is given by a(ω, ε) = a 0 [1 + ε cos (ωz)], where ε and ω stand for the perturbation amplitude and the growth rate, respectively. The dispersion relation for the Rayleigh-Plateau instability suggests that the largest growth rate, ω, is obtained for a wavenumber k ∼ 0.6970, which corresponds to λ max = 2π/k = 9.016a 0 (where λ max is the most unstable wavelength). All simulations were run using an initial perturbation amplitude of ε = 0.25 with k ∼ 0.6970, unless stated otherwise in the text.
Additionally to the dimensionless flow parameters, the flow dynamics are also controlled by the ratio of the film thickness to the radius of the rim, i.e. e = h 0 /a 0 . In this study, we consider the zero acceleration limit, and e = 0.2, following the work by Agbaglah et al. (2013) where their stability analysis and dispersion relationship depended on the initial acceleration. As the liquid sheet retracts, the thickness of the rim grows over time as the rim engulfs the sheet, e decreases, and the rim deceleration vanishes resulting in the rim dynamics being fully driven by the RP instability. Therefore, the flow dynamics depend on two competing time scales: the time scale given by the RP instability in the (spanwise) z-direction (T RP ), and the time scales of the rim growth (T ret ). The time scale of the action of surface tension and droplet detachment is given by the capillary time, T RP = ρ l a 3 0 /σ s . On the other hand, mass conservation results in the scaling a(da/dt) ∼ hu T C , leading to T ret = a/(da/dt) ∼ ρ l a 4 0 /(2σ s h 0 ); for T RP << T ret , then, the liquid sheet must satisfy h 0 /a 0 << 1, as a 2 0 ∼ h 0 L 0 , which leads to the next condition 4 h 0 /L 0 << 1 that is accomplished at very high aspect ratios (see Mirjalili et al. (2018) for more details). Finally, we justify the initial sinusoidal perturbation by looking into the breakup time T B as a function of T B = f (Oh, ε) (see expression 5 from Driessen et al. (2013)); for The computational domain corresponds to 6.75λ max × λ max × λ max (i.e., 304.29h 0 × 45.08h 0 × 45.08h 0 ) where the x-coordinate is aligned to the length of the sheet. The domain is sufficiently large in the streamwise direction to avoid the effect of artificial reflections from the boundary of this direction. The simulations are initialised with fluids at rest in the absence of gravity. The periodic boundary condition is imposed at all the variables in the (spanwise) z-direction of the domain. A no-penetration boundary condition is prescribed for the bottom and top domain. The domain is discretised with a regular cubic mesh of size ∆x = h 0 /6, ensuring that our numerical predictions are mesh independent; our results do not vary with decreasing grid size. Note that similar mesh sizes have been used in previously studies, Agbaglah (2021) Figure 1b contrasts our numerical results against the expected Taylor-Culick retraction velocity for Oh = 0.1, = 0 and e = 2. As seen, the constant retracting velocity is ∼ 8% smaller than the predicted Taylor-Culick velocity. This difference is consistent with observation by Agbaglah (2021) and Song & Tryggvason (1999) and are explained by the large density and viscosity ratios. The retracting dynamics of a liquid thread and non-linear dynamics before pinch-off have been validated previously by Constante-Amores et al. (2020). We refer to Shin et al. (2018) for the accuracy of the surfactant equations.

Results
We start the discussion of the results by presenting a phenomenological picture of the interfacial dynamics in a phase diagram in a β s -Oh space. These dynamics range from a nearly inviscid sheet (Oh ∼ 10 −3 ) to typical soap films (Oh ∼ 10 −1 ). Figure 2 shows the regime map in terms of the interfacial dynamics predicted from our numerical simulations. Our results show that two distinct regimes can be identified based on the outcome of the retracting dynamics. At low β s and Oh, surfactant-driven Marangoni stresses do not promote the reopening of the adjacent sheet connected to the rim; thus, capillary waves result in the formation of a hole behind the rim. This hole grows in the spanwise direction and separates the rim from the sheet. In fact, these dynamics have previously been observed by Mirjalili et al. (2018) and Constante-Amores et al. (2021b). The rest of the β s -Oh space is dominated by surfactant-induced Marangoni stresses which reopen the adjacent sheet to promote droplet detachment via RP instability. Below, we focus on this regime, and provide an extensive explanation of the interfacial dynamics.
The spatio-temporal interfacial dynamics for the surfactant-free case, with Oh = 0.0833, are found in Figures 3a-e. As seen, the initial rim instability grows to develop nonlinearities that eventually lead to a capillary singularity (break up). At early stages of the simulation, the fluid is pulled from the rim to form an elongated ligament (see figure 3c). Over time, the ligament undergoes end-pinching, where capillarity overcomes viscosity to form a bulbous end (see figure 3d) terminating in the break up of a droplet (see figure 3e). The interfacial shape prior to the capillary singularity closely resembles experimental observations by Zhang et al. (2010) and , 2021. Additionally, we observe surface waves resulting from the development of capillary instabilities at the upstream of the rim.
Attention is now turned to the effect of surface-active agents by varying the surfactant elasticity parameter, β s , with Oh = 0.0833, P e s = 100 and Γ = Γ ∞ /2. For β s = 0.1, t = 49.49 t = 63.63 t = 106.06 t = 141.42 t = 189.50  the spatio-temporal interfacial dynamics resemble that of a surfactant-free case, from the formation of the ligament to the eventual end-pinching: the dynamics results in droplet shedding (see figure 3f-j). In this case, surfactant is accumulated at the base of the cylindrical rim as it gradually increases its size. Once the ligament is formed, τ induces a flow towards the ligament (see figure 3h) resulting in a longer ligament prior to pinch-off. However, as surfactant-induced Marangoni-stresses are not sufficiently strong to evacuate Γ , the dynamics results in modest surfactant concentration gradients, and, in agreement with Constante-Amores et al. (2020), the ligaments breaks up. By increasing β s , τ act to reduce Γ -gradients by inducing a flow from the rim towards the sheet, and from the rim towards the ligament (regions of low surfactant concentration), resulting in a severe reduction of the rim radius (displayed in figure 3). The surfactant-driven flow leads to the reopening of the fluid sheet adjacent to the retracting rim thereby increasing its film thickness (see figure 4i); the Marangoni-induced flow also acts to suppress the capillary waves ahead of the rim (see figure 3k-o and 3j-n for β s = 0.3 and β s = 0.5, respectively). We now turn our attention to the two-dimensional interfacial shape represented by Γ and spanwise Marangoni stresses, the tangential component of the surface velocity u tz , the streamwise surface velocity, u tx (where the velocity is given in the reference frame of u x at z = 0) presented in figure 4 at t = 141.42 and t = 189.50, respectively. As observed, the surfactant-driven flow, on the spanwise direction, results on the retardation of the dynamics brought by surfactant-induced interfacial rigidification (see the dampening of u tx in the frame of difference in figure 4b,f: as β s > 0, u tx < 0). We observe that surfactant concentration gradients at t = 141.42 give rise to Marangoni stresses which drive a flow from the horizontal part of the rim towards the base of the ligament, and from the tip of the ligament to its base (for high β s , the ligament tip has a nearly uniform distribution of Γ ). The combined effect of Marangoni stresses is to bridge the surfactant gradient between the tip of the ligament and the rim (see figure 4h). This effect is sufficiently strong to promote the development of a thinner bulbous end, and cause a surfactant-driven escape from end-pinching to a longer ligament (in agreement with Constante-Amores et al. (2020) and Kamat et al. (2020)). As expected, in a latter stage, the nearly uniform distribution of surfactant concentration results in the elimination of Marangoni stresses (see figure 4h). Figures 4c and g show that that u tz > 0 upstream and u tz < 0 downstream of the neck, i.e. an stagnation point is found in between. As the neck narrows, capillary-driven flow causes further neck thinning and the creation of a velocity maximum at u tz , as illustrated in figure 4c,g which is indicative of singularity formation. The u tz profile for the surfactant-free case is characterized by the presence of a large velocity maximum and two stagnation points, with the neck located in between. By close inspection of the u tz profile of the surfactant-laden cases it becomes clear that only one stagnation point is found near the neck for high β s . The surfactant-induced Marangoni stresses result in the suppression of one of the stagnation points-in agreement with Constante-Amores et al. (2020, 2021c. Figure 4j,k shows the vorticity field, i.e. ω = × u, in terms of the interface location for the surfactant-free and surfactant-laden cases for t = 141.42 and t = 189.50 (the location is given in a reference frame where the the rim position is found at z = 0). For all cases, a primary vortex (with a positive sign) is formed at the front region of the bulbous edge and is shed into the ambient fluid. For the surfactant-free case, we observe a negative vortex-ring near the bulbous neck resulting from a change of curvature; this vortex ring grows over time ending in breakup. For the surfactant-laden case at t = 141.42, τ induce vorticity from the bulbous tip to the neck and from the rim to the ligament. This flow causes the neck to reopen; after escaping the singularity the flow through the neck triggers the formation of a jet towards the centre of the bulbous region (these findings are consistent with Constante-Amores et al. (2020)). Additionally, the escape from pinch-off causes the length of the ligament to grow over time in the (x)direction, forming a secondary bulbous and a secondary neck as well as a shift in the axial curvature signs (see figure 4k).
To provide conclusive evidence that Marangoni stresses drive the reopening of the neck, we run an additional simulation in which Marangoni stresses, or surface tension gradients, are deactivated, i.e. ∇ s σ = 0, to isolate the effects of a reduction in surface tension from those rising from Marangoni stresses. The Marangoni-suppressed case resembles the surfactant-free dynamics culminating in end-pinching (see Supplemental animation). Additionally, we performed an additional surfactant-free simulation at the "effective" Ohnesorge number, obtained by using the surface tension achieved by the surfactant, i.e. Oh ef f = µ/ √ ρh 0 σ 0 = 0.103 (for β s = 0.5 and Γ 0 = 0.5) finding similar results. Consequently we claim Marangoni stresses are responsible of the change of interfacial dynamics discussed above.
Finally, figure 5 shows the effect of surfactant dynamics on the filament's tip location, L, the kinetic energy, E k = V (ρu 2 )/2dV, and the length of the ligament (defined from the tip of the bulbous edge to its rim). Here, the kinetic energy has been normalised by the surface energy E s = A 0 σ s , where A 0 is the initial surface of the sheet. The evidence for a surfactant-driven retardation of the dynamics can be seen by inspection of E k and the streamwise location of the tip, which reveals that increasing β s monotonically decreases the overall value of both E k and L. Thus, the interfacial flow dynamics are retarded by Marangoni stresses resulting in a reduction in the retraction velocity. In figure 5c, we report the length of the ligament until pinch-off, the addition of surfactant increases the ligament length prior to its breakup becoming particularly pronounced at high β s .

Concluding remarks
Results from three-dimensional numerical simulations of the retracting dynamics of thin liquid sheets in the presence of insoluble surfactants were presented. The liquid properties in the simulations were chosen to be consistent with a realistic air/water system, and the retracting velocity of the rim formed at the end is in good agreement with the Taylor-Culick speed. For the surfactant-laden cases, we have demonstrated that surfactant-induced Marangoni stresses drive a flow from the high surfactant concentration regions to the low concentration ones, reducing the kinetic energy, affecting the location of the bulbous tip and suppressing end-pinching. With increasing elasticity number, the Marangoni stresses play a major role in the interfacial dynamics with the progressive elimination of the capillary wave structures upstream of the rim, where a complete interfacial rigidification is observed for large elasticity numbers. Additionally, Marangoni stresses drive flow towards the sheet resulting in the reopening of the film thickness adjacent to the rim. Future research avenues are related to examine the role of solubility on the flow dynamics.