Retraction of thin films coated by insoluble surfactants

Abstract

We investigate the retraction of a circular thin film coated with insoluble surfactants, which is punctured at its centre. We assume that the surface pressure of the liquid–gas interface is related to the number density of surfactants through a linear equation of state, which is characterized by a single parameter: the Gibbs dilation modulus. To solve the governing equations and track the deformation of the domain, we use the finite element method with an arbitrary Lagrangian–Eulerian approach where the film surface is sharp. Our simulations show that the surface elasticity introduced by the surfactants slows down the retraction and introduces oscillations at early times. In agreement with previous experiments and theoretical analysis, we find that the presence of surfactants introduces perturbations of the film thickness over progressively larger distances as the surface elasticity increases. The surface perturbations travel faster than the retracting edge of the film at a speed proportional to the Gibbs modulus. For large values of the Gibbs modulus, the interface behaviour approaches that of an incompressible two-dimensional solid. Our analysis sheds light on the effect of insoluble surfactants on the opening of a circular hole in a thin film and can be extended to investigate the onset of surface cracks and fractures.

1. Introduction

Bubble ruptures are important in natural processes (e.g. formation of aerosol droplets and mass transfer between the sea and the atmosphere) (Eggers & Villermaux 2008; Lhuissier & Villermaux 2009) and in industrial processes (e.g. separation technologies such as devolatilization and distillation, destabilizing foams) (Yarin et al. 1999). The reported observations of bubble ruptures in the scientific literature show the retraction of thin films with velocities and shapes that depend on viscosity, surface tension and geometrical parameters (i.e. bubble radius and thickness). In particular, the velocity of the retracting film immediately after the bubble rupture has been studied extensively (Debrégeas, De Gennes & Brochard-Wyart 1998; Bird et al. 2010).

In 1867, Duprè (1867) studied the kinetics of growth of a hole in a film and concluded that a rolled-up rim collects all the material of the disappearing film and the rest of the film remains undisturbed. Hence, assuming that all the surface energy was transformed into kinetic energy, he derived from the energy balance the retracting velocity of the rim, $V =\sqrt { \phi \gamma /\rho H}$, where $\gamma$ is the surface tension of the liquid–gas interface, $\rho$ is the density of the liquid, $H$ is the film thickness, and $\phi$ is a factor for which Duprè derived a value of 4. In 1959 and 1960, Taylor (1959) and Culick (1960) revisited the work of Duprè and independently derived a value $\phi =2$, meaning that not all the surface energy is transformed into kinetic energy. Culick (1960) also brought experimental evidence re-interpreting the experiments from Ranz (1959) that were showing a deviation from the Duprès predictions of about 10 %.

The effect of the liquid viscosity on the shape and the dynamics of planar and circular retracting films was studied more recently both experimentally (Debrégeas, Martin & Brochard-Wyart 1995) and theoretically (Schulkes 1996; Brenner & Gueyffier 1999; Song & Tryggvason 1999; Notz & Basaran 2004; Savva & Bush 2009). These works showed that increasing the viscosity of the liquid flattens the retracting film, prevents the accumulation of the liquid into the rim, and changes the duration of the initial transient acceleration, but it does not change the final speed attained by the film. More recent work has considered also the impact of viscoelasticity (Tammaro et al. 2018; Villone, Hulsen & Maffettone 2019), viscoplasticity (Deka, Pierson & Soares 2019) and shear-thinning viscosity (Kamat, Anthony & Basaran 2020a) on the retraction of planar and axisymmetric films.

The liquid–gas interfacial tension of the retracting film is influenced strongly by the type and the concentration of the surface-active species (surfactants, proteins, polymers, colloids, etc.), and the introduction of a surfactant for tuning the interfacial properties should be accompanied by consideration of its physicochemical properties, including its bulk and interfacial diffusivities, adsorption/desorption kinetics at interfaces, surface viscosity and elasticity (Fuller & Vermant 2012). Typically, small-molecule surfactants like sodium dodecyl sulfate (SDS) display a measurable elasticity under dilational deformations. This behaviour is introduced by changes in surface pressure in response to changes in local concentration, and it is normally referred to as Gibbs elasticity (Bhamla et al. 2014). It can be understood as the consequence of an equation of state relating the surface pressure to the surfactant concentration. On the other hand, proteins, polymers, nanoparticles or other big molecules adsorbed at a fluid–gas interface may experience strong intermolecular interactions, which lead to the emergence of additional shear and dilational surface elasticity that is not obtained from an equation of state (Sagis 2011; Pepicelli et al. 2017). In fact, the surface elastic stresses in these monolayers can display deviatoric components (Gu & Botto 2016) and typically exceed greatly the surface pressure associated with the Gibbs modulus (Frostad et al. 2016).

McEntee & Mysels (1969) performed a systematic experimental study of the retraction speed of films with different thicknesses and displaying ‘mobile’ and ‘rigid’ liquid–gas interfaces. In the case of the ‘mobile film’ (i.e. a solution of water and small-molecule surfactants), they found a modest agreement with Culick's prediction for thicknesses bigger than 1000 Å. For film thinner than 1000 Å, the deviation became increasingly larger until the retraction speed reached around 4 % of the expected value. For the first time, they observed the formation of an aureole: a film with an increased thickness that extends beyond the retracting rim on the hole border. In the case of the so-called ‘rigid films’, a solution that forms two-dimensional crystal structures at the interface, the retracting velocity was always smaller than that predicted by the Taylor–Culick formula, $V =\sqrt {2 \gamma /\rho H}$. In this case, the aureole is not observed, but a buckling of the film and out-of-plane deformations are reported. The buckling is a consequence of the instantaneous local surface tension becoming negative, which means that the surface pressure of the adsorbed layer has exceeded the surface tension of the clean liquid–gas interface (Milner, Joanny & Pincus 1989). Similar observations were reported by Petit, Le Merrer & Biance (2015), who found that the rupture of circular films with increasing concentrations of surfactant displays an increasing number of ‘cracks’ (similar to the buckling observed by McEntee & Mysels 1969) and a decreased retraction speed. More recently, folding of the liquid–air interface was also observed along the diagonal of square-shaped soap films that are retracting (Habibi & Krechetnikov 2021). In this case, the formation of folds is attributed to the convergence of the interfacial compression wave along the diagonals. In agreement with the earlier works, the authors observed a slower retraction speed than that predicted by the Taylor–Culick formula. Finally, even more complicated dynamics is observed when insoluble micron-sized particles are used to coat the interface of the retracting liquid film (Timounay, Lorenceau & Rouyer 2015).

A large body of theoretical and numerical investigation focused on the effects of soluble and insoluble surfactants and surface viscosity on the stability of jets (Martínez-Calvo & Sevilla 2018; Constante-Amores et al. 2020; Martínez-Calvo et al. 2020; Wee et al. 2021) and breakup of liquid threads (McGough & Basaran 2006; Kamat et al. 2018; Martínez-Calvo & Sevilla 2020; Wee et al. 2020). These works show that interfacial viscosity and Marangoni stresses slow down the thinning of liquid filaments and can completely prevent their rupturing. On the other hand, computational and theoretical analysis of the dynamics of a hole nucleated on a thin film has relied mostly on thin-film approximations of planar and circular geometries (Brenner & Gueyffier 1999; Savva & Bush 2009; Munro et al. 2015). The effect of Gibbs elasticity introduced by the surface-active species was included in simple models by Frankel & Mysels (1969) for a planar geometry and by Petit et al. (2015) for a circular geometry. These simplified models explained the mechanism of formation of the aureole and allowed computation of the retraction speed of the film. Lu, Campana & Corvalan (2018) studied the effect of surfactants on the contraction of a small hole in a thin film, showing that surfactants slow down the closing of a small pore. We are not aware of detailed numerical simulations of the retraction of liquid films or of films displaying more solid-like interfacial rheological responses such as yield stress or stress relaxation (Jaensson, Anderson & Vermant 2021). Going beyond simple interfacial constitutive laws considering surface viscosity (Scriven 1960) is especially important because interfacial viscoelasticity can have a profound role in the bubble bursting dynamics and can lead to mechanical instabilities of the liquid–gas interface, which then dominate the dynamics of the bursting film (Petit et al. 2015; Habibi & Krechetnikov 2021; Tammaro et al. 2021).

In this paper, we present detailed finite element simulations of the retraction of a circular thin film of liquid coated with insoluble surfactants. We use an arbitrary Lagrangian–Eulerian (ALE) formulation to track the sharp air–liquid interface, and we study the impact of the Gibbs elasticity introduced by the surfactants on the retraction speed and on the shape of the film. This work represents a first step into the investigation of the onset of cracks and fractures on the surface of retracting circular holes.

2. Problem formulation

We consider a film of thickness $H$ made of a Newtonian liquid of density $\rho$ and viscosity $\mu$ that is surrounded by air. At time $t=0$, a circular hole of initial radius $r_0(0)=10 H$ is nucleated at the centre of the film and grows in time as the thin film retracts progressively in the radial direction. The choice $r_0(0) \gg H$ guarantees that the hole does not close due to the azimuthal curvature (Lu et al. 2018). A schematic description of this system is displayed in figure 1. The initial length of the film is denoted by $L$, and to avoid finite-size effects, we assume $L \gg H$. At time $t=0$, the film is assumed to close through a semicircle of radius $H/2$; see figure 1. We consider a cylindrical coordinate system with the origin at the centre of the hole, and assume that the problem is symmetric around the $z$-axis; see figure 1. The motion of the liquid in the film is determined by the Navier–Stokes equations (e.g. Aris 2012)

(2.1)\begin{equation} \rho \left( \frac{\partial \boldsymbol{v}}{\partial t} + \boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{v} \right)= \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{T} , \end{equation}

where $\boldsymbol {v}$ is the velocity of the liquid, and the stress tensor $\boldsymbol {T}$ follows the Newtonian constitutive law

(2.2)\begin{equation} \boldsymbol{T}={-} p \boldsymbol{I} + \mu \left(\boldsymbol{\nabla} \boldsymbol{v} + \boldsymbol{\nabla} \boldsymbol{v}{}^{\rm T} \right) , \end{equation}

with $p$ the pressure in the liquid, which enforces the incompressibility of the liquid:

(2.3)\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{v} =0. \end{equation}

The pressure in the outside air is assumed to be zero. At time $t=0$, the fluid is considered to be quiescent.

Figure 1. Schematic description of the retraction of a thin film coated with insoluble surfactants. At time $t=0$, a hole of radius $r_0(0)$ is nucleated at the centre of the film. The film has thickness $H$ and length $L \gg H$. The problem is axisymmetric around the $z$-axis. The fluid is Newtonian and incompressible, and it retracts towards the right edge. The initial surface number density of surfactants is $c_{s,0}$.

The air–liquid interface of the film is coated with insoluble surfactants whose local surface number density is denoted by $c_s$. Following previous works (Lu et al. 2018; Kamat et al. 2020b), at the initial time, the number density of surface-active molecules is homogeneous and equal to $c_{s,0}$. The interfacial tension $\gamma$ between the liquid and the air is related to the local number density $\gamma (c_s)$ through an equation of state that can be derived directly from thermodynamic arguments. This behaviour has been termed Gibbs elasticity. Here, we assume that the surface number density of surfactants is much smaller than its value at maximum packing, and that the air–liquid surface tension decreases linearly with $c_s$, until it reaches zero (Petit et al. 2015):

(2.4)\begin{equation} \gamma (c_s) = \max\left(0, \gamma_0+\left.\frac{{\rm d} \gamma}{{\rm d} c_s}\right|_{c_{s,0}} c_s \right). \end{equation}

At the point $\gamma (c_s^{*})=0$, we considered the surface to be buckled and to support no tension (Marmottant et al. 2005). This can be expressed using the piecewise function that is shown in figure 2. In (2.4), $\gamma _0$ is shorthand notation for $\gamma (c_{s,0})$ and indicates the surface tension at the initial time $t=0$. The constant $({\rm d} \gamma /{\rm d} c_s)|_{c_{s,0}}$ is the linear dilation module of the interface and is negative, which means that increasing the concentration of surfactants decreases the surface tension of the air–liquid interface. The coefficient $({{\rm d} \gamma }/{{\rm d} c_s})|_{c_{s,0}}$ is commonly referred to as the Gibbs modulus. Strictly speaking, the surface tension given by (2.4) is not differentiable at the point where $\gamma (c_{s}^{*})=0$, and this could be problematic when solving the nonlinear system obtained from the discretization. However, in our simulations, we find that $c_{s}$ is always smaller than $c_{s}^{*}$, and no special treatment is required to guarantee that the surface tension is well behaved. This model has been used previously by Petit et al. (2015), and it represents a simplification of the true dependence of the interfacial tension on the number density of surfactants, which is nonlinear in the general case. In Appendix B, we show that considering a logarithmic equation of state yields results that are similar qualitatively to those obtained with the linear equation of state given by (2.4).

Figure 2. Plot of the surface tension normalized with $\gamma _0$ as a function of the area density of surfactants rescaled with $c_{s,0}$. The slope of the curve is equal to minus the elastocapillary number defined as $E = -({c_{s,0}}/{\gamma _0})({{\rm d} \gamma }/{{\rm d} c_s})|_{c_{s,0}}$. The interface cannot attain negative tension.

The balance of stresses at the interface is given by

(2.5)\begin{equation} \boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{T} ={-}\gamma (c_s)\, \boldsymbol{n}\, \boldsymbol{\nabla}_s \boldsymbol{\cdot} \boldsymbol{n}+\boldsymbol{\nabla}_s \gamma (c_s), \end{equation}

with $\boldsymbol {n}$ the unit vector normal to the film interface pointing into the air; see figure 1. We have neglected the viscous stresses due to the airflow, which could lead to flapping instabilities (Lhuissier & Villermaux 2009). In this work, we also neglect the effect of surface viscosity and viscoelasticity (Fuller & Vermant 2012; Jaensson et al. 2021), which will be the object of future studies. We assume a no-stress boundary condition at the end of the film, $r=r_0+L$:

(2.6)\begin{equation} \boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{T} = \boldsymbol{0}. \end{equation}

However, the choice of a different boundary condition has a negligible effect on the retraction of the film because the domain is very large, $L\gg H$ and $L \gg r_0$.

In typical experimental conditions, once the hole has been nucleated, the retraction of the film occurs over a time scale of a few milliseconds (Savva & Bush 2009; Petit et al. 2015). This is much faster than the rate of mass transfer of surfactants between the bulk and the surface, which can thus be neglected (Lin, Frostad & Fuller 2018). It follows that the evolution of the number density of surfactants is described by

(2.7)\begin{equation} \frac{\partial c_s}{\partial t} + \boldsymbol{\nabla}_s \boldsymbol{\cdot} \left( c_s \boldsymbol{v}_s \right)+ c_s ( \boldsymbol{\nabla}_s \boldsymbol{\cdot}\boldsymbol{n} ) (\boldsymbol{v} \boldsymbol{\cdot}\boldsymbol{n}) = D_s\, \nabla_s^{2} c_s , \end{equation}

where $\boldsymbol {v}_s$ is the projection of the velocity on the plane tangent to the fluid–air interface,

(2.8)\begin{equation} \boldsymbol{v}_s = ( \boldsymbol{I}- \boldsymbol{n}\boldsymbol{n}) \boldsymbol{\cdot}\boldsymbol{v} , \end{equation}

and $D_s$ is the diffusion coefficient of the surfactants along the surface. Similarly, $\boldsymbol {\nabla }_s$ denotes the gradient operator tangential to the interface,

(2.9)\begin{equation} \boldsymbol{\nabla}_s= ( \boldsymbol{I}- \boldsymbol{n}\boldsymbol{n}) \boldsymbol{\cdot} \boldsymbol{\nabla} . \end{equation}

The second term on the left-hand side of (2.7) describes the advection of surfactants along the interface. The third term on the left-hand side describes the range of change of number density due to local dilation or compression of the interface.

The governing equations are made dimensionless using the thickness of the film $H$ as the characteristic length, the speed $\sqrt {\gamma _0 / \rho H}$ as the characteristic velocity, the initial surface tension $\gamma _0$ as the characteristic tension, and the initial surface number density $c_{s,0}$ as the characteristic number density. The characteristic stresses and the characteristic time are then given by $\sqrt {\mu ^{2} \gamma _0 / \rho H^{3}}$ and $\sqrt {\rho H^{3} / \gamma _0 }$, respectively. We find that the problem is determined by three dimensionless numbers. The Ohnesorge number,

(2.10)\begin{equation} Oh = \sqrt{\frac{\mu^{2} }{\rho \gamma_0 H}} , \end{equation}

compares viscous forces to inertial and surface tension forces. The surface tension forces are evaluated using the value of the surface tension at the initial number density of surfactants, $\gamma _0$. The elastocapillary number,

(2.11)\begin{equation} E = \left.-\frac{c_{s,0}}{\gamma_0}\,\frac{{\rm d} \gamma}{{\rm d} c_s}\right|_{c_{s,0}} , \end{equation}

compares the Gibbs elastic modulus to the initial surface tension of the film. For $E \gg 1$, the liquid–air interface is very elastic, and small changes in $c_s$ yield large changes in surface tension. Conversely, for $E \ll 1$, the surfactants have little effect on the surface tension. In the limit $E \rightarrow 0$, the case of a clean interface is recovered. The elastocapillary number can be related to the surfactant strength parameter (Leal 2007; Kamat et al. 2020b), in the limit of surfactant concentration being much smaller than its maximum value. Finally, the surface Péclet number,

(2.12)\begin{equation} Pe_s =\sqrt{\frac{\gamma_0 H}{\rho D_s}} , \end{equation}

compares the rate of transport of surfactants by advection to the rate of transport by diffusion along the surface. Owing to the fast time scale of the hole opening process, the Péclet number is typically very large. For typical surfactants adsorbed on the interface of a thin film of water, we have $Pe_s \approx 10^{5}$. In this work, we fix $Pe_s=1000$ and study the effects of $E$ and $Oh$ on the retraction of the film. We verified that increasing $Pe_s$ did not change the results presented.

To solve the set of equations (2.1)–(2.7), we use the finite element method with an ALE formulation to track the interface deformation (Villone et al. 2014, 2019). Briefly, the velocity of the mesh, $\boldsymbol {v}_{m}$, at any point in the computational domain is obtained by solving a Laplace equation ${\nabla }^{2}\boldsymbol {v}_{m}=\boldsymbol {0}$ with Dirichlet boundary conditions $\boldsymbol {v}_{m}=\boldsymbol {v} \boldsymbol {\cdot } \boldsymbol {n}\boldsymbol {n}$ on all the boundaries. For more details on the mesh deformation algorithm, we refer to the work of Villone et al. (2014). An example of a highly deformed mesh obtained with our method is shown as a supplementary figure available at https://doi.org/10.1017/jfm.2022.412. This method captures correctly the location of the film surface, $h(r,t)$, even when it is not a single-valued function of $r$ (Notz & Basaran 2004). We use a triangular mesh that is more refined near the hole and becomes coarser away from it. We use quadratic interpolation for $\boldsymbol {v}$ and $c_s$, and linear interpolation for $p$. To avoid finite-size effects, the length of the computational domain is chosen to be always larger than $L=1400 H$. We do not enforce the symmetry with respect to the line $z=0$; instead, we simulate the entire liquid film. The time evolution of the fields is computed using an adaptive second-order implicit Euler method. The nonlinear system of equations resulting from the discretization is solved using the Newton method with relative tolerance $10^{-4}$. This value of the tolerance was sufficient to give convergent and accurate numerical results. In our simulations we fixed $H= \rho = \gamma _0= c_{s,0}=1$ and $D_s=0.001$. By choosing the above values, the dimensionless Taylor–Culick velocity is equal to $\sqrt {2}$. We control the Ohnesorge number and the elastocapillary number $E$ by varying the viscosity $\mu$ and the linear dilation elastic modulus $({{\rm d} \gamma }/{{\rm d} c_s})|_{c_{s,0}}$. The numerical results are presented in dimensionless form.

3. Validation of the code

To validate the ALE finite element method used here, we perform simulations with no surfactants, $E=0$, and compare the results with previous results by Savva & Bush (2009), who studied the retraction of a thin film of viscous liquid using the lubrication approximation. We consider an initial hole with $r_0(0)=50$ and $Oh=1$. In figure 3, we report the numerical results compared with the data extracted from figures 16 and 17 of Savva & Bush (2009). In figure 3(a), we show the evolution of the retraction speed $\dot {r}_0$, and in figure 3(b), we report the height profile of the film $h(r,10)$ at time $t=10$. The numerical simulations agree with the results of Savva & Bush (2009), with the small discrepancy between the two datasets being probably due to errors in manually extracting data from their plots. In Appendix A, we report additional validation of the numerical method in the case of a surfactant-laden interface, $E>0$.

Figure 3. Validation of the numerical method against the results of Savva & Bush (2009). (a) Retraction speed of the film scaled by the Taylor–Culick velocity, $\sqrt {2}$. (b) Film thickness profile at $t=10$.

4. Numerical results

We now use our numerical method to study the retraction of a liquid film with insoluble surfactants. All the simulations presented in this section are performed considering quiescent initial conditions. We consider values of the elastocapillary number $E$ that are comparable to those found in the experiments of Petit et al. (2015) and Mitrinova et al. (2013).

4.1. Retraction speed with insoluble surfactants

We begin our analysis by investigating the effects of surfactants on the retraction speed of the film. In figure 4, we report the evolution of the retraction speed $\dot {r}_0$ for different elastocapillary numbers and for two Ohnesorge numbers. It is apparent that in the case $E \neq 0$, the curves display a qualitatively different behaviour from the case of a surface with zero elasticity, $E=0$. The surface elasticity introduced by the surfactants generates oscillations at early times, which are then damped progressively. Simulations performed at $Oh=1$ show that by increasing the importance of the liquid viscosity, the oscillations of the retraction speed are damped.

Figure 4. Retraction speed of the film, $\dot {r}_0$, normalized with the Taylor–Culick speed, plotted as a function of time for (a) $Oh=0.1$ and (b) $Oh=1$. The surface elasticity introduced by the surfactants generates oscillations at early times and reduces the final retraction speed.

To investigate the origin of the oscillations of the retraction speed, we plot the surfactant concentration at different times during the transient in figure 5. For all the cases investigated, there is an initial fast compression regime that is dominated by the inertial acceleration of the liquid. During this fast retraction regime, the speed is insensitive to the presence of surfactants, which is demonstrated in figure 4 by the collapse of the speed computed at different $E$ onto the speed computed at $E=0$. Figure 5 shows that the initial retraction quickly increases the surfactant concentration, but it always remains below the value $1+E^{-1}$ at which the surface tension becomes zero; see figure 2. After the initial inertial acceleration, the decrease of surface tension at the rim slows down the retraction speed. The fast increase of surface concentration at the retracting rim also generates a steep gradient of surfactant concentration. In all the cases shown in figure 5, the surface concentration first increases and then decreases, achieving a maximum value near the retracting tip. Since the Marangoni stresses are proportional to the surface concentration gradient, they are tangential to the surface but act in an opposite direction before and after the concentration maximum. This means that in principle, they could contribute positively or negatively to the retraction speed. However, we show that the Marangoni stresses before and after the maximum do not contribute equally to the retraction speed. Indeed, at early times, the film is only slightly deformed and the tangential vector is mostly vertical in the region $r-r_0<1$ and essentially horizontal for $r-r_0>1$. In the region $r-r_0>1$, the surfactant concentration is decreasing steeply and the tangent vector is aligned with the $r$-axis, resulting in Marangoni stresses directed along the motion of the retracting edge. Conversely, the Marangoni stresses in the region $r-r_0<1$ are directed mostly along the $z$-axis, and they do not contribute to the retraction speed. In summary, the sudden initial inertial compression of the surface generates a local accumulation of surfactants that slows down the retracting rim. At longer times, the Marangoni stresses drive viscous flows that increase the velocity of the rim, thus explaining the oscillations observed in figure 4. Finally, we find that the surface elasticity does not significantly change the length of the transient regime, which appears to be comparable to the case of a clean surface, $E=0$.

Figure 5. Profile of surfactant concentration at the surface of the film computed at different times during the initial transient regime. The initial fast acceleration of the film quickly increases the concentration of surfactants near the retracting edge, $r -r_0 \approx 0$. A strong gradient of surfactant concentration is developed during this initial phase.

In figure 6, we plot the steady-state retraction speed, normalized with the Taylor–Culick speed, as a function of $E$. In the limit of a clean surface with no surfactants, $E \rightarrow 0$, our results recover the Taylor–Culick speed. As $E$ is increased, the speed of the retracting film decreases. As discussed in Savva & Bush (2009), we find that the steady-state retraction speed is not sensitive to the initial shape of the rim; see Appendix C. Similarly, the same steady-state retraction speed is achieved if a concentration of surfactants smaller than in the rest of the film is considered in the circular rim as the initial condition; see Appendix C. In the same plot, we compare our simulations with the results of Petit et al. (2015), who used a simplified lubrication model to calculate the retraction speed. Our results agree with their simplified model at large elastocapillary numbers, but deviate from them at moderate and small $E$, because some of the hypotheses considered in their simplification break down at $E \approx 1$. Note that, in contrast to the case of clean films, circular and planar surfactant-laden films attain different steady-state retraction speeds (Petit et al. 2015).

Figure 6. Steady-state retraction speed of the film, normalized with the Taylor–Culick speed and plotted as a function of the elastocapillary number. The Ohnesorge number has no effect on the final retraction speed. At large $E$, our results agree with the simplified model of Petit et al. (2015).

In figure 6, it is apparent that numerical simulations at different Ohnesorge numbers demonstrate that viscosity does not play any role in the steady-state retraction speed, but it changes only the transient values. As already discussed by Savva & Bush (2009) for a thin film without surfactants, viscous stresses are internal forces and do not contribute to the rate of change of the total momentum of the film, $\boldsymbol {P}$. The argument by Savva & Bush (2009) can be extended to the case of a surfactant-laden film by considering an integral balance of momentum:

(4.1)\begin{equation} \frac{{\rm d}}{{\rm d}t} \boldsymbol{P} = \frac{{\rm d}}{{\rm d}t} \int_{V(t)} \rho \boldsymbol{v} \,{\rm d}V= \int_{S(t)} \boldsymbol{T} \boldsymbol{\cdot} \boldsymbol{n} \,{\rm d}S+\int_{S_{bound}} \boldsymbol{T} \boldsymbol{\cdot} \boldsymbol{n} \,{\rm d}S , \end{equation}

where $V(t)$ and $S(t)$ are the volume and the surface of the retracting film, and $S_{bound}$ is the boundary positioned at the end of the film, $r_0+L$. By substituting the boundary conditions (2.5) and (2.6), we notice that the last surface integral in the right-hand side vanishes, and the first surface integral in the right-hand side is rewritten as

(4.2)\begin{equation} \frac{{\rm d}}{{\rm d}t} \int_{V(t)} \rho \boldsymbol{v} \,{\rm d}V= \int_{S(t)} \left( -\gamma (c_s)\, \boldsymbol{n}\,\boldsymbol{\nabla}_s \boldsymbol{\cdot} \boldsymbol{n}+\boldsymbol{\nabla}_s \gamma (c_s) \right) {\rm d}S . \end{equation}

Equation (4.2) shows that, similarly to the case of a clean film, the liquid viscosity does not contribute to the rate of change of the film momentum when surfactants coat its surface. To confirm the insight obtained from equation (4.2), we plot in figure 7 the radial component of the total momentum of the film $P_r$. The curves at different $Oh$ values collapse on each other, thus confirming the conclusions drawn from (4.2).

Figure 7. Radial component of the total momentum of the film plotted as a function of time, for different values of $Oh$ and $E$.

To investigate the mechanism responsible for the slowdown of the steady-state retraction in the case $E>0$, we consider the balance of energy in the film. To do so, we take the scalar product of the momentum balance, given by (2.1), with $\boldsymbol {v}$, and then integrate over the film volume using the Reynolds transport theorem to obtain

(4.3)\begin{equation} \frac{{\rm d}}{{\rm d}t} \int_{V(t)} \frac{\rho}{2}\,\boldsymbol{v}^{2} \,{\rm d}V = \int_{V(t)} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{T} \boldsymbol{\cdot} \boldsymbol{v} \, {\rm d}V . \end{equation}

Using integration by parts in the integral on the right-hand side, and using the boundary conditions (2.5) and (2.6), we obtain

(4.4)\begin{align} \frac{{\rm d}}{{\rm d}t} \int_{V(t)} \frac{\rho}{2}\,\boldsymbol{v}^{2} \,{\rm d}V +\int_{V(t)} \boldsymbol{T} \boldsymbol{:} \boldsymbol{\nabla} \boldsymbol{v} \,{\rm d}V ={-} \int_{S(t)} \gamma (c_s)\, \boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{n}\,\boldsymbol{\nabla}_s \boldsymbol{\cdot} \boldsymbol{n} \,{\rm d}S +\int_{S(t)} \boldsymbol{\nabla}_s \gamma (c_s) \boldsymbol{\cdot} \boldsymbol{v} \, {\rm d}S . \end{align}

The two integrals on the left-hand side represent the rate of change of kinetic energy and the rate of dissipated energy, respectively. The two integrals on the right-hand side represent the power input that drives the retracting process. The integrand in the first integral on the right-hand side is the product of the surface tension times the rate of change of the surface area due to the normal dilation. Since the surface area of the film decreases in time and the surface tension is always positive, this integral is always positive and introduces energy in the system. The second integral represents the effect of the Marangoni stresses. The sign of this term can change in time during the transient, but at steady state, the product $\boldsymbol {\nabla }_s \gamma (c_s) \boldsymbol {\cdot } \boldsymbol {v}$ is positive because surface flow is directed from regions of low surface tension to regions of large surface tension.

Equation (4.4) states that the power input generated by the Marangoni stresses and by the reduction of surface energy is partly transformed into kinetic energy and partly dissipated. In the case of a clean surface, $E=0$, there are no Marangoni stresses, $\boldsymbol {\nabla }_s \gamma (c_s)=0$, and the surface tension is constant. In the case of a planar film without surfactants, Savva & Bush (2009) showed that the surface energy is divided equally into kinetic and dissipated energy. We can use (4.4) to analyse the effect of surfactants on the rate of change of energy in the retracting film. Their presence has two opposite effects: (i) it reduces power input due to surface energy by decreasing the value of the surface tension; and (ii) it generates Marangoni stresses that can increase the power input. Besides changing the total power input that drives the retraction of the film, the presence of surfactants also impacts the rate of change of kinetic energy and the rate of change of dissipated energy. To understand how the presence of surfactants changes the different terms appearing in (4.4), in figure 8 we plot the rate of change of kinetic energy ${\rm d} E_{K}/{\rm d}t$, the dissipated energy $\mathcal {D}$, and the power input ${\rm d} E_\gamma /{\rm d}t$, defined as

(4.5)$$\begin{gather} \frac{{\rm d} E_{K}}{{\rm d}t}=\frac{{\rm d}}{{\rm d}t} \int_{V(t)} \frac{\rho}{2}\,\boldsymbol{v}^{2} \,{\rm d}V, \end{gather}$$

(4.6)$$\begin{gather} \mathcal{D} =\int_{V(t)} \boldsymbol{T} \boldsymbol{:} \boldsymbol{\nabla} \boldsymbol{v} \,{\rm d}V, \end{gather}$$

(4.7)$$\begin{gather} \frac{{\rm d} E_\gamma}{{\rm d}t}={-} \int_{S(t)} \gamma (c_s)\, \boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{n}\, \boldsymbol{\nabla}_s\boldsymbol{\cdot} \boldsymbol{n} \,{\rm d}S +\int_{S(t)} \boldsymbol{\nabla}_s \gamma (c_s) \boldsymbol{\cdot} \boldsymbol{v} \, {\rm d}S . \end{gather}$$

Figure 8(a) reveals that increasing the elastocapillary number reduces the overall power input, ${\rm d} E_\gamma /{\rm d}t$, compared to the case of a clean interface, $E=0$. This means that the increased contribution of the Marangoni stresses on the right-hand side of (4.4) is not sufficient to make up for the reduction of surface tension due to the increased surfactant concentration. While this leads to a reduction of the power input that drives the retraction process, figure 8(b) shows that the surface elasticity also reduces the rate of energy dissipated in the film, $\mathcal {D}$. As a consequence, more energy is converted into kinetic energy in the case $E>0$ than in the case $E=0$. However, the reduction of the dissipated energy is not sufficient to make up for the reduction of the power input due to the surfactants. Indeed, figure 8(c) shows that in the case $E>0$, the rate of change of the total kinetic energy of the film is smaller than in the case $E=0$. These findings explain the mechanism behind the slowing down of the steady-state retracting speed observed in figure 6. Finally, figure 8(d) shows that the numerical error in the computation of (4.4) is smaller than 1 % in the cases $E=10$ and $E=100$.

Figure 8. Plot of the different terms appearing in the balance of energy, given by (4.4), for the case $Oh=0.1$. The plots in (a) represent the power input due to the surface tension and the Marangoni stresses. The curves in (b) represent the rate of energy dissipation. In (c), we report the rate of change of the total kinetic energy of the film. (d) Sum of the terms of (4.4) normalized with ${{\rm d} E_\gamma }/{{\rm d}t}$ shows that the deviation is always smaller than 1 % in the cases $E=10$ and $E=100$.

4.2. Surface elasticity effects on the film thickness

We now turn our attention to the effects of surfactant-induced elasticity on the shape of the retracting film that is quantified using $h(r)$, shown schematically in figure 1. In figure 9, we report the thickness of the film at different times, for $Oh=1$. In figure 9(a), we report the case of a clean surface, $E=0$. In this case, the fluid is collected into the rim of the retracting front, and the film thickness $h(r)$ is disturbed a few dimensionless distances from the hole, reaching its far-field value right after the rim neck. As a consequence, the volume of the rim grows progressively as the film retracts and the hole grows, and the thickness of the film near the edge appears to grow continuously in time. Conversely, in the case of a surfactant-laden interface with $E=10$, figure 9(b) shows that the rim of the retracting film appears to stop growing at longer times and to reach a much smaller thickness. Another difference from the case $E=0$ is that the disturbance of the film profile $h(r)$ decays on a much larger scale; see the inset in figure 9(b). This result agrees with the experimental observations (McEntee & Mysels 1969; Liang, Chan & Choi 1996) and theoretical predictions (Frankel & Mysels 1969; Petit et al. 2015) that the presence of surfactants drives the formation of an aureole of film with increased thickness ahead of the retracting rim. As demonstrated in figure 5, the compression of the interface near the retracting hole increases the concentration of surfactants and generates a gradient of surface tension. As a consequence, the fluid entrained by the retracting film is no longer collected into the rim but is pushed down the film by the surface tension gradient (Frankel & Mysels 1969), thus increasing progressively its thickness over a longer length scale.

Figure 9. The thickness of the film, $h(r)$, at different times in the case $Oh=1$. (a) The case of a clean interface, $E=0$. (b) The case $E=10$.

The difference between the shape of a retracting clean interface, $E=0$, and a retracting surfactant-laden interface, $E \neq 0$, is illustrated further in figures 10 and 11, where we show the shape of the film and the radial component of the velocity at different times. In figure 10, the retracting clean interface collects the fluid in a rim that grows progressively in time, and the radial velocity decays to zero very close to the rim neck, leaving the rest of the film unperturbed. Conversely, figure 11 shows that the rim of an interface with $E=10$ does not grow in time. Instead, the radial velocity decays over a much longer length and perturbs a much larger portion of the film. Such a long-ranged velocity field pushes the fluid away from the retracting edge, which stops the growth of the rim. This phenomenon is a consequence of the accumulation of surfactants at the surface of the retracting edge. In turn, this generates a gradient of surfactant concentration (see figure 5), which drives the long-ranged fluid flow displayed in figure 11.

Figure 10. Radial velocity field inside a retracting film with no surface elasticity, $E=0$, and moderate viscous effects, $Oh=0.1$. Most of the liquid is collected inside the rim.

Figure 11. Radial velocity field inside a retracting film with a surfactant-laden surface displaying $E=10$ and moderate viscous effects, $Oh=0.1$. The rim does not grow in time, and the fluid entrained by the retracting edge is pushed far from the rim by gradients of surface tension.

We quantify the maximum height attained by the retracting film as a function of time, and plot it in figure 12. In the case of a clean interface, $E=0$, the maximum height grows in time as more liquid is collected into the rim. Conversely, in the case $E>0$, after the initial transient growth, the maximum height of the film plateaus and reaches a steady-state value. In this regime, the fluid is no longer collected into the rim as the hole progressively opens. By comparing figures 12(a) and 12(b), we see that increasing the Ohnesorge number changes the transient and the steady-state value of the maximum height attained by the film. To confirm that the gradients of surface tension push fluid away from the retracting rim and stop the growth of the film thickness, we report the radial component of the velocity along the surface in figure 13(a) for the case $Oh=1$ and $E=10$. The plot shows that the velocity decays slowly with the distance from the retracting edge because it is driven by the long-ranged Marangoni stresses. In figure 13(b), we show that the long-ranged velocity field shown in figure 13(a) is such that at long times, $\partial v_{r}/\partial r + v_r/r$ is zero over a large portion of the film. It follows from the conservation of mass that $\partial v_{z}/\partial z \approx 0$ along most of the film, which then retracts essentially without changing its thickness. There is a small region near the edge of the film where the radial velocity is changing in space and the interface is deforming. This region is where the long-time change of the maximum thickness reported in figures 12(a,b) takes place. Qualitatively similar results are obtained for different values of $Oh$ and $E$. The finding that most of the film retracts without changing its thickness is compatible with the reduction of dissipated energy compared to the case of a clean interface as shown in figure 8.

Figure 12. Maximum height of the retracting film is plotted as a function of time for (a) $Oh=0.1$ and (b) $Oh=1$. For $E>1$, the maximum height reaches a steady-state value and the rim of the retracting film stops growing.

Figure 13. (a) Radial component of the velocity evaluated at the surface of the film at different times. The flow field decays slowly with the distance from the retracting edge. (b) Local spatial change of radial velocity $\partial v_{r}/\partial r + v_r/r$. There is a large region of the film where two consecutive elements of fluid move at the same speed, which means that the film thickness does not change. The dimensionless numbers used in these plots are $Oh=1$ and $E=10$.

In figure 14, we plot the steady-state value of the maximum rim height as a function of $E$, for three Ohnesorge numbers. The maximum height decreases as $E$ increases, and it tends towards an Ohnesorge-dependent asymptotic value at large $E$. This means that, as the interface becomes more elastic, the volume of fluid collected into the rim of the retracting film decreases. To understand why the maximum height attained at steady state depends on the Ohnesorge number, we investigate the concentration of surfactants near the retracting edge. Indeed, as shown in figures 10–14, this region is where the film attains its maximum height. In figures 15(a) and 15(b), we report the surfactant concentration profile at different times in the case $E=10$, for two values of the Ohnesorge number. It is apparent that there is a small region near the edge where the surfactant concentration increases and then decreases. This non-monotonic behaviour is more pronounced in the case $Oh=1$ than in the case $Oh=0.1$. This change of slope of the surfactant concentration profile near the edge implies that the radial component of the Marangoni stresses, $-E \, \boldsymbol {\nabla }_s c_s$, changes sign. Figure 16 shows that the change of sign of the Marangoni stresses generates a recirculating field inside the retracting rim. Such a flow field tends to reduce the maximum height of the film. Therefore, since the gradient of surface concentration is larger in the case $Oh=1$ than in the case $Oh=0.1$, the maximum thickness of the film is smaller.

Figure 14. Maximum height of the retracting film at steady state plotted as a function of the elastocapillary number $E$. The maximum height decreases with $E$ and seems to reach an $Oh$-dependent asymptotic value.

Figure 15. Profiles of surfactant concentration near the retracting edge in the case $E=10$, for (a) $Oh=1$, and (b) $Oh=0.1$. The surface concentration first increases and then decreases, leading to Marangoni stresses that are directed oppositely.

Figure 16. The $z$-component of the velocity field for the case $Oh=1$, $E=10$ and $t=10$. The velocity attains both positive and negative values due to a weak recirculation inside the retracting edge, which is a consequence of the change of sign of the Marangoni stresses.

The results shown so far demonstrate that the presence of surfactants and the gradients of surface tension change significantly the shape of the retracting film. Their effect is to push the liquid entrained by the rim down the film, leading to a long-ranged deformation of the film compared to the case of a clean interface. We now investigate the speed at which the disturbance of the film height interface propagates during the opening of the hole. The propagation speed is computed by looking at the largest radial position of the interface, $r^{*}$, that has $|h(r^{*})-0.5| >\epsilon$, with $\epsilon$ being a numerical tolerance set to $\epsilon =10^{-8}$. By tracking the evolution of $r^{*}$ in time, we compute the propagation speed of the thickness disturbance of the film. An example of the evolution of $r^{*}$ is shown in figure 17. After an initial transient, $r^{*}$ grows linearly in time, and we compute the propagation speed from the slope of a linear fit. By repeating this procedure for different $E$, we construct the plot of the propagation speed as a function of the elastocapillary number, shown in figure 18. It is apparent that the propagation speed follows two regimes. At small values of $E$, the propagation speed approaches asymptotically the Taylor–Culick speed. This means that the disturbance of the film thickness travels slightly faster than the retraction speed of the film. Conversely, for large values of $E$, the propagation speed grows as $1.75 \, E^{1/2}$. The scaling of the propagation speed at large $E$ agrees with the simplified model used by Petit et al. (2015), but we find a slightly larger coefficient $1.75$ rather than $\sqrt {2}$. Simulations performed at different Ohnesorge numbers show that the propagation speed of the disturbance is insensitive to the viscosity of the film. At large $E$, the interface is behaving effectively as a compressible two-dimensional solid. Indeed, the propagation speed of the height perturbation is proportional to the two-dimensional dilation modulus, $({{\rm d} \gamma }/{{\rm d} c_s})|_{c_{s,0}}$, which is similar to the case of the speed of sound for an elastic solid (Liang et al. 1996). As $E$ is increased, the elasticity of the interface becomes progressively more important approaching the behaviour of an incompressible interface. In the limit $E \rightarrow \infty$, the interface is effectively incompressible, and the perturbation of the film reaches the edges instantaneously.

Figure 17. Plot of the evolution of the largest radial position where the film thickness is perturbed from its initial value. The position of the thickness perturbation $r^{*}$ grows in time and reaches a final speed. The linear fit is used to extract its propagation speed. The slope of the Taylor–Culick velocity is reported for comparison.

Figure 18. Plot of the propagation speed of the perturbation of the film thickness as a function of the elastocapillary number. At large $E$, the propagation speed grows as $1.75 \, E^{1/2}$.

4.3. Low Ohnesorge simulations and damping of capillary waves

As the Ohnesorge number decreases, the viscosity becomes less important and capillary waves are formed during the retraction of the film. For very small values of $Oh$, the capillary waves locally reduce the thickness to a small fraction of its initial value, which can lead to the breakup of the retracting film and filaments (Eggers 1993; Notz & Basaran 2004; Castrejón-Pita, Castrejón-Pita & Hutchings 2012; Anthony et al. 2019) and the formation of satellite droplets and bubbles (Bird et al. 2010). Previous studies on the retraction of jets showed that soluble surfactants have a stabilizing effect against the breakup of the film by damping the amplitude of the capillary waves (Constante-Amores et al. 2020). In this subsection, we discuss the effect of surfactants on the capillary waves formed during the hole opening process, and the fact that the Marangoni stresses can prevent the pinch-off of almost inviscid retracting filaments (Kamat et al. 2020b).

In figure 19, we report the results of simulations performed for $Oh=0.01$. In the absence of elasticity, $E=0$, figure 19(a) shows that capillary waves are propagated in front of the retracting rim, and their amplitude appears to grow in time. The growth of the capillary waves can lead eventually to the formation of a very thin neck and the subsequent breakup of the film into droplets. Conversely, figures 19(b,d) show that the presence of surfactants damps the propagation of capillary waves along the film, and hinders their amplitude. By increasing the elastocapillary number, the elasticity introduced by the surfactants becomes progressively more important than the surface tension, which damps the growth of the capillary waves. The damping of the capillary waves is likely an effect of the vorticity developed due to the gradient of surfactants, as discussed in Kamat et al. (2020b).

Figure 19. Film thickness profiles at different times for $Oh=0.01$ and (a) $E=0$, (b) $E=1$, (c) $E=5$, and (d) $E=10$. The surface elasticity introduced by the surfactants reduces the amplitude of the capillary waves.

5. Conclusions

In this work, we performed a detailed numerical study of the retraction of a thin axisymmetric film of liquid with insoluble surfactants adsorbed onto the air–liquid interface. The presence of surfactants changes the surface tension locally, leading to an elastic response of the interface upon dilation or compression. We used a sharp-interface model, and we tracked its deformation using an arbitrary Lagrangian–Eulerian (ALE) approach. To model the elasticity of the surfactant-laden interface, we used a simplified equation of state that relates linearly the surface tension to the number density of surfactants until zero tension is reached, at which point we assume that the liquid–gas interface bears no tension.

We find that the retraction speed of the hole undergoes oscillations at early times, which are caused by the elasticity of the surfactant-coated interface. These oscillations are caused by the fast initial compression of the interface, which generates a local accumulation of surfactant near the hole edge that slows down the retraction. After the slowdown, the speed increases again due to the formation of a steep gradient of concentration of surfactant and large Marangoni stresses that promote the opening of the hole. The steep gradients of surfactants are then smoothed progressively over longer distances by the Marangoni flow. At longer times, the retraction speed reaches a steady-state value that decreases as the importance of surface elasticity over surface tension increases. Our findings are in agreement with previous experimental observations (Petit et al. 2015; Habibi & Krechetnikov 2021) and theoretical considerations (Frankel & Mysels 1969). The balance of energy in the film reveals that the slowdown of the steady-state speed of the retracting edge is caused ultimately by the decrease of the surface tension due to the accumulation of surfactants.

The numerical simulations show that the profile of the retracting film is changed significantly with respect to the case of a clean interface. In the case of a clean interface, the fluid entrained by the retracting film is collected in the rim, which grows in time as the hole opens progressively. In the case of a surfactant-coated interface, instead the rim reaches a steady-state size after an initial transient. The steady-state size of the rim decreases as the elasticity of the interface increases. In this case, the gradients of surface tension prevent the growth of the rim by pushing the liquid entrained during the retraction down the film length. As a consequence of mass conservation, the liquid entrained by the hole-opening process is distributed over much larger distances, and the film thickness is perturbed over distances much larger than the rim size. The perturbation to the film thickness has been termed an aureole in previous works. We find that the disturbance of the film thickness propagates ahead of the rim with a speed that depends on the elastocapillary number. At large values of the elastocapillary number, the disturbance propagates with a speed proportional to the square root of the Gibbs modulus, thus behaving effectively as a compressible two-dimensional solid. Finally, we find that the presence of surfactants damps the capillary waves that are formed for small values of the Ohnesorge number.

The present work represents a first step towards understanding the effects of a more complicated and solid-like interfacial rheology on the retraction of thin films. Previous studies showed that high-speed deformations of structured interfaces possessing solid-like behaviour and stress relaxation may lead to qualitatively different behaviours compared to their quasistatic counterpart (Petit et al. 2015; Pitois, Buisson & Chateau 2015; Poulichet & Garbin 2015; Huerre, De Corato & Garbin 2018; Garbin 2019; Tammaro et al. 2021). In future works, we will investigate the retraction of fluid with interfaces that have complex behaviour, including pressure-dependent surface viscosity, concentration-dependent bending modulus, shear elasticity, stress relaxation and surface yield stress (Fuller & Vermant 2012; Beltramo et al. 2017; Jaensson et al. 2021). We believe that including some of these features in our simulations could shed light on the mechanical instabilities observed in the retraction of films with structured interfaces (Petit et al. 2015; Tammaro et al. 2021).