On the role of variable surface viscosity in free viscous films

TL;DR

This paper investigates how variable surface viscosity, influenced by surfactant concentration, affects the stability and rupture behavior of free viscous films through theoretical modeling, stability analysis, and nonlinear simulations.

Contribution

It introduces a coupled nonlinear model incorporating variable surface viscosity and tension effects, revealing new stability criteria and rupture dynamics for free films with surfactants.

Findings

Surface viscosity significantly influences film stability.

Jamming surfactant concentration can stabilize films.

Rupture profiles differ markedly with variable surface viscosity.

Abstract

The stability of a thin liquid film bounded by two free surfaces is examined in the presence of insoluble surface active agents. The surface active agents not only cause gradients in surface tension, but could also render surface viscosity to be significant, and variable, as a function of their concentration. A set of three coupled nonlinear evolution equations are derived for the film height, concentration of surface active agents, and the horizontal liquid velocity which governs the dynamics of the free film. Suitable phenomenological models for variation of surface tension and surface viscosity with the concentration of surface active agents are incorporated in the interfacial stress boundary condition. Linear stability analysis reveals the effect of various non-dimensional parameters, specifically the retarding nature of surface viscosity and Marangoni effects on film rupture. An analysis of the `jamming' limit of surfactant concentration reveals that $\Gamma_0^{(nl)}<3 \mathbb{D}/M$ is a sufficient criteria for instability of the system, where $\Gamma_0^{(nl)}$ is a normalized initial surface particle concentration, $\mathbb{D}$ is the disjoining pressure number and $\textit{M}$ is the Marangoni number. Nonlinear simulations suggest film profiles at rupture are also qualitatively different from those reported in earlier studies, revealing that free films in the jamming limit are remarkably stable and their free surfaces behave like immobile interfaces. Furthermore, self-similar exponents extracted from our nonlinear simulations are used to explain topological differences between zero and constant surface viscosity in the vicinity of rupture.

Tables

Table 1. Table 1: Estimates of relevant physical parameters

Parameter	Definition	Estimate	units
$H$	Mean film thickness	${10}^{-9}-{10}^{-7}$	$m$
$\mu$	Liquid viscosity	${10}^{-3}$	$kg{m}^{-1}{s}^{-1}$
$\rho$	Liquid density	${10}^3$	$kg{m}^{-3}$
$A$	Hamaker constant	${10}^{-21}-{10}^{-19}$	$J$
$\widetilde{{\sigma }_0}$	Surface tension of clean film	${10}^{-3}-{10}^{-2}$	$N{m}^{-1}$
$\widetilde{{\Gamma }_0}$	Mean surfactant concentration	${10}^{-6}$	$mol{m}^{-2}$
$D_s$	Surface diffusion coefficient	${10}^{-12}-{10}^{-8}$	$m^2{s}^{-1}$
$\widetilde{{\eta }_0}$	Surface viscosity	${10}^{-12}-{10}^{-6}$	$kg{s}^{-1}$

Table 2. Table 2: Estimates of non-dimensional system parameters evaluated using characteristic velocity scale, and dimensional estimates from Table 1 .

Parameter	Definition	Order of magnitude estimates
$ϵ=\frac{{\tilde{h}}_0}{L}$	Aspect ratio	${10}^{-3}-{10}^{-2}$
$\text{Bo}=\frac{\widetilde{{\eta }_0}}{\mu h_0}$	Boussinesq number	${10}^{-1}-{10}^1$
$𝔻=\frac{A}{6\pi {\tilde{\sigma }}_0ϵL^2}$	Disjoining pressure number	${10}^{-3}-{10}^{-2}$
$\text{M}=\frac{{\Gamma }_0}{{\tilde{\sigma }}_0}\left(\frac{\partial \tilde{\sigma }}{\partial \tilde{\Gamma }}\right)$	Marangoni number	$0-5\times {10}^{-3}$
$\beta =\frac{{\Gamma }_0}{{\tilde{\eta }}_0}\left(\frac{\partial \tilde{\eta }}{\partial \tilde{\Gamma }}\right)$	Surface viscosity gradient number	$0-1$
$\text{Pen}=\frac{A}{6\pi \mu h_0{D}_s}$	Peclet number	${10}^{-1}-{10}^2$
$\text{Rey}=\frac{\rho A}{6\pi {\mu }^2ϵ{\tilde{h}}_0}$	Reynolds number	${10}^{-2}-1$

Equations

$$\rho ({\tilde{u}}_t+\tilde{u}{\tilde{u}}_x+\tilde{v}{\tilde{u}}_z)=-({\tilde{P}}_{\tilde{x}}+{\tilde{\Phi }}_{\tilde{x}})+\mu ({\tilde{u}}_{\tilde{x}\tilde{x}}+{\tilde{u}}_{\tilde{z}\tilde{z}}),$$

$$\rho ({\tilde{u}}_t+\tilde{u}{\tilde{v}}_x+\tilde{v}{\tilde{v}}_z)=-({\tilde{P}}_{\tilde{z}}+{\tilde{\Phi }}_{\tilde{z}})+\mu ({\tilde{v}}_{\tilde{x}\tilde{x}}+{\tilde{v}}_{\tilde{z}\tilde{z}}),$$

$${\tilde{u}}_{\tilde{x}}+{\tilde{v}}_{\tilde{z}}=0.$$

$$\frac{\partial \tilde{\Gamma }}{\partial \tilde{t}}+({\nabla }_s\cdot \text{n})\tilde{\Gamma }(\text{n}\cdot \text{u})+{\nabla }_s\cdot (u_s\tilde{\Gamma })=D_s{\nabla }^2\tilde{\Gamma }.$$

$$-\text{n}\cdot ||\text{mathsfbi}T||\cdot \text{n}=2\kappa \sigma +2\kappa \left(k^s+{\mu }^s\right){\nabla }_s\cdot \text{u},$$

$$-\text{t}\cdot ||\text{mathsfbi}T||\cdot \text{n}=\text{t}\cdot {\nabla }_s\sigma +\left(k^s+{\mu }^s\right)\text{t}\cdot {\nabla }_s{\nabla }_s\cdot \text{u}+\text{t}\cdot {\nabla }_s\left(k^s+{\mu }^s\right){\nabla }_s\cdot \text{u},$$

$${\tilde{h}}_{\tilde{t}}+\tilde{u}{\tilde{h}}_{\tilde{x}}=\tilde{v}(\mathrm{kinematic}),$$

$${\tilde{u}}_{\tilde{z}}=0,and\tilde{v}=0(\mathrm{squeezing}\mathrm{mode}\mathrm{symmetry}).$$

$$z=\frac{\tilde{z}}{H},x=\frac{\tilde{x}}{L},$$

$$u=\frac{\tilde{u}}{A/6\pi \mu ϵL^2},v=\frac{\tilde{v}}{A/6\pi \mu {ϵ}^2{L}^2},P=\frac{\tilde{P}}{A/6\pi ϵL^3}.$$

$$t=\frac{\tilde{t}}{6\pi \mu {ϵ}^3{L}^3/A}.$$

$$\Gamma (x,t)=\frac{\tilde{\Gamma }(\tilde{x},\tilde{t})}{{\tilde{\Gamma }}_0}.$$

$$\sigma (x,t)=\frac{\tilde{\sigma }(\tilde{x},\tilde{t})}{\widetilde{{\sigma }_0}};\eta (x,t)=\frac{\tilde{\eta }(\tilde{x},\tilde{t})}{{\tilde{\eta }}_0}.$$

$$\sigma =1-\text{M}\Gamma (x,t).$$

$$\eta =1+\beta \left(\Gamma (x,t)-1\right),$$

$$\eta =\frac{1}{{[1-\Gamma (x,t)]}^{\alpha }}.$$

$$\frac{{\tilde{\Gamma }}_0}{{\tilde{\Gamma }}_{max}}=\frac{\beta }{\alpha \left(1-\beta \right)}=\frac{{\Gamma }_0^{(nl)}}{{\Gamma }_0^{(l)}}$$

$$\sigma =1,and{\sigma }_x=0.$$

$$u_{zz}^{(0)}=0,$$

$$P_z^{(0)}=v_{zz}^{(0)},$$

$$u_x^{(0)}+v_z^{(0)}=0.$$

$$u^{(0)}=c(x,t),$$

$$v^{(0)}=-c_{x}z$$

$$-P^{(0)}={\mathbb{D}}^{-1}{h}_{xx}+2h_x{u}_z-2v_z.$$

$$P^{(0)}=-2c_x-{\mathbb{D}}^{-1}{h}_{xx}^{(0)}.$$

$$u_z^{(0)}={\mathbb{D}}^{-1}{\sigma }_x=-{ϵ}^2{\mathbb{D}}^{-1}\text{M}{\Gamma }_x^{(0)}=0.$$

$$h_t+{(ch)}_x=0,$$

$${\Gamma }_t+{(c\Gamma )}_x=\frac{{\Gamma }_{xx}}{\text{Pen}},$$

$$\text{Rey}(c_t+c{c}_x)-{\mathbb{D}}^{-1}{h}_{xxx}-\frac{3}{8h^4}{h}_x-4c_{xx}=4\frac{h_x{c}_x}{h}+\text{Bo}(1+\beta (\Gamma -1))\frac{c_{xx}}{h}$$

$$+\text{Bo}\beta \frac{c_x{\Gamma }_x}{h}-{\mathbb{D}}^{-1}\text{M}\frac{{\Gamma }_x}{h}.$$

$$(\mathrm{for}\mathrm{linear}\mathrm{model}\mathrm{using}(\text{ref}eta1))$$

$$\text{Rey}(c_t+c{c}_x)-{\mathbb{D}}^{-1}{h}_{xxx}-\frac{3}{8h^4}{h}_x-4c_{xx}=4\frac{h_x{c}_x}{h}+\text{Bo}\frac{c_{xx}}{(1-\Gamma {)}^{\alpha }h}$$

$$+\alpha \text{Bo}\frac{c_x{\Gamma }_x}{(1-\Gamma {)}^{\alpha +1}h}-{\mathbb{D}}^{-1}\text{M}\frac{{\Gamma }_x}{h}.$$

$$(\mathrm{for}\mathrm{NPM}\mathrm{using}(\text{ref}eta2))$$

$$h(x,t)=\frac{1}{2}+\hat{h}{e}^{(ikx+st)},$$

Full text

1 Introduction

Surface active agents are often used to modify surface rheological properties of dispersed media such as liquid foams, emulsions, and soap films. The interactions of these particle-laden interfaces are of great interest in several industrial applications like foams and detergents, inks (Tadros, 2014), Pickering emulsions (Anjali & Basavaraj, 2018), groundwater treatment (Ye et al., 2012), and soil remediation (Mao et al., 2015). Often surface active agents become important through their inadvertent presence, for instance, as contaminants in coating or printing processes. They typically alter the wetting characteristics and hydrodynamics of a system by giving rise to interfacial stresses and affecting the surface properties of the system. The relevant surface properties of interest are surface tension, Marangoni effects and concentration-dependent surface viscosity. A unified model that can characterize the interplay between these effects is necessary to truly predict the stability of such systems. As subject of the present work, we formulate and solve a mathematical model that addresses some of these interplays in a thin free liquid film, covered with insoluble surface-active agents at its free surface.

The configuration of a free liquid film has direct practical relevance to soap bubbles and cosmetic foams. A soap bubble is a thin spherical shell of water with air at a slightly elevated pressure trapped inside. Drainage of the aqueous layer due to gravity would lead to thinning and van der Waals forces, eventually causes rupture. Similar dynamics are seen in Pickering emulsions (Ramsden, 1903; Pickering, 1907). Though the process is conceptually simple, the physics involved in determining the precise bubble breakup time is complicated. The presence of surfactants create Marangoni stresses at the fluid-liquid interface which will delay the breakup time, (de Gennes, 2001). On the other hand, high pressure inside the bubble can promote rupture, often leading to catastrophic breakup. The idealised model studied here also has great relevance to bubbly suspensions and emulsions, though we do not make a distinction between the two in the current study. The model examined in the present work idealizes the space between two bubbles in an emulsion as a semi-infinite thin film of liquid in the longitudinal direction, with two free surfaces above and below it (See figure 1). The free surfaces are covered with a layer of surface active agents (hereafter simply referred to as ‘surfactants’) that are insoluble in the bulk of the film. Though we use the word ‘surfactants’ to describe any surface active agents, some of the findings of the study is expected to be relevant to interfacial flows with tiny colloidal particles on an interface. We neglect gravitational effects in the current study owing to the very thin films ($O(100nm)$) considered here. In such thin films, long range van der Waals forces are included, and for simplicity we ignore any stochastic effects such as Brownian motion. Our configuration differs from that of (Naire et al., 2001, 2000; Braun et al., 2000) who studied a vertical free film undergoing gravitational drainage, but owning to very small gravitational effects for such thin films, we believe our work may have relevance even to draining films.

Before proceeding further, it is useful to highlight existing literature on modelling stability of thin free liquid films. Several prior investigations (Prevost & Gallez, 1986; Hatziavramidis, 1992; Gallez et al., 1993; Erneux & Davis, 1993; Vaynblat et al., 2001) have modelled the dynamics of thin free films in the framework of lubrication theory, accounting for effects of capillarity, evaporation, condensation, and bulk viscosity. Among these, Hatziavramidis (1992) was the earliest to explore the role of surfactants, which illustrated the conflict between thermally- and surfactant-induced Marangoni flows. In their work, the dynamics of a free film were approximated by those of a wall film with half the thickness of a free film. A detailed presentation of the nonlinear evolution equations using long-wave theory for a clean free film was first given by Erneux & Davis (1993). The concept was later extended to include insoluble surfactants at the surface by DeWit et al. (1994), for which one arrives at a system of three coupled nonlinear evolution equations governing the film dynamics. Furthermore, this work was extended to soluble surfactants by Chomaz (2001), but he ignored the role of surface viscosity. It was not until the work of Edwards & Oron (1995) and later by Naire et al. (2001, 2000) and Braun et al. (2000) that surface viscosity was included in these models. Naire et al. (2001) considered the configuration of a vertically draining film while Edwards & Oron (1995) studied a horizontal film without surfactants and Marangoni effects. A model incorporating Marangoni effects, surface viscosity, and bulk solubility of surfactants in horizontal free films was studied by Matar (2002). However, in all these works, surface viscosity is always assumed to be constant throughout the surface and independent of surfactant concentration.

In the present work, we focus on horizontal non-draining films with surfactant-laden interfaces. We neglect the effect of solubility focussing primarily on surface rheology. In particular, we use a model framework that allows for a variation of surface viscosity as function of the surfactant concentration. Our work primarily stems from a need to describe systems where surfactants or colloidal particles exhibit large concentration variations on the interface. Danov et al. (1999) show that at sufficiently high surfactant concentrations, the interface appears to form a rigid shell accompanied by an enhancement in surface viscosity. Rising spherical bubbles and sedimenting drops covered with surfactants are known to have a drag coefficient different from that of a clean interface. In the limit of Stokes flows, such systems were modelled using the ‘spherical cap’ approximation (Leal, 2007) where a part of the drop surface is assumed to be rigid (immobile) while the remainder of the drop assumed to be devoid of surfactants (Sadhal & Johnson, 1983). Diffusion of surfactants on the interface makes such an approximation unrealistic and one is expected to find a smooth variation in concentration continuously reducing from leeward to windward side of drop/bubble. Thus, in addition to Marangoni effects, we also introduce surface viscosity effects which depend on local concentration of surfactants as demonstrated in prior experiments by Lopez & Hirsa (2000). This allows us to have a continuous variation in mobility on the interface consistent with diffusion effects of surfactants. A nonlinear phenomenological model describing this dependence is incorporated in the governing equations to construct a robust framework.

We adopt a scaling based on viscous forces balancing the inter-molecular van der Waals interactions. A lubrication approximation is taken in the equations owing to the small aspect ratio of the free films. It is found that, variable surface viscosity leads to new self-similar solutions for film evolution in the ‘dilute’-limit when surface active agents are sparsely distributed across the surface and the sensitivity of surface viscosity toward particle concentrations is weak. In addition, we also find that surfactants could potentially render a free film remarkably stable in the ‘jamming’-limit.

We first derive a set of governing equations by non-dimensionalizing the Navier-Stokes equations and interfacial boundary conditions in §2. We explore the linear stability of the system in §3 and apply perturbation techniques to reveal an important criteria in §3.2. A brief discussion of numerical methods and the numerical solutions to the governing equations is provided in §4. This is followed by self similar solutions found in the ‘dilute’ regime in §4.2.2, and important findings in the ‘jamming’ limit in §4.3. Finally, conclusions are made in §5, connecting our results to prior experiments and findings.

2 Problem formulation

2.1 Problem geometry and governing equations

Figure 1 shows the idealized two-dimensional problem setup. A thin film of incompressible liquid, extending infinitely in the lateral direction, is bounded by a passive gas phase above and below it. The free surfaces on both sides are covered with surfactants which are assumed to be insoluble in the bulk of the film. These agents could represent either surfactants in the classical sense, significantly affecting the liquid surface tension, or serve as colloidal particulates which alter the surface rheology. The film has mean thickness $2H$ (dimensional), and the liquid has a bulk viscosity and density represented by $\mu$ and $\rho$ (both dimensional) respectively. The free surface is at a height $\tilde{h}(\tilde{x},\tilde{t})$. Further, we let $\tilde{\sigma}_{0}$ denote the dimensional mean surface tension of the liquid film when the surface concentration of surfactants $\tilde{\Gamma}$ is maintained at at a fixed reference concentration $\tilde{\Gamma}_{0}$. Table 1 gives some realistic dimensional values for the various physical parameters used in our study.

Following prior work (DeWit et al., 1994; Erneux & Davis, 1993; Matar, 2002), we will consider only the squeezing (or varicose) mode which is symmetric about the horizontal ($\tilde{x}$) axis, and is thus considered the most unstable mode. The dimensional momentum balance in the $\tilde{x}$- and $\tilde{z}$ directions, and the continuity equation may be written as:

[TABLE]

[TABLE]

[TABLE]

Here, $\tilde{u}$ and $\tilde{v}$ denote the velocity components in $\tilde{x}$ and $\tilde{z}$ directions, while $\tilde{P}$ is the fluid pressure. We have included a disjoining pressure term $\tilde{\Phi}$ in our model, considering films to be thinner than $\sim 100nm$, where van der Waals forces are prevalent (Edwards & Oron, 1995; Berg, 2010). In the present work, we set $\tilde{\Phi}=A/6\pi(2\tilde{h})^{3}$ where $A$ is the dimensional Hamaker constant representing van der Waals attraction between the two free surfaces separated by the liquid. The attractive van der Waals forces will be the prime destabilizing factor for the film. The dynamics of the surface active species at the free surface neglecting their solubility in the bulk, is governed by a convection-diffusion equation (Matar, 2002):

[TABLE]

Here $D_{s}$ is surface diffusivity, treated as a constant, $\nabla_{s}$ is the surface gradient operator, $u_{s}$ is the surface velocity field, and n is the unit vector normal to the free surface. The boundary conditions at the free surface ($\tilde{z}$= $\tilde{h}$($\tilde{x}$, $\tilde{t}$)) are the normal stress continuity and the shear stress balance. The vector form of these equations in the presence of variable surface tension and surface viscosity effects are taken from the work by Naire et al. (2001), written as:

[TABLE]

[TABLE]

where, t is the unit vector tangent to the free surface, $\mathsfbi{T}$ is the stress tensor and the dilatational $(k^{s})$ and shear $(\mu^{s})$ components of surface viscosities will hereafter be written in terms of a single parameter, $\tilde{\eta}=(k^{s}+\mu^{s})$ in the rest of the paper. A simplified version of the above equations in scalar form have been stated in Appendix A. The other conditions that complete the system are the kinematic condition at the free surface, and the symmetry condition at the center line ($\tilde{z}=0$). These are written as:

[TABLE]

[TABLE]

2.2 Scalings and non-dimensionalization

The scalings for the different variables are explained next, where symbols with the tilde ($\sim$) decoration denote dimensional versions. The spatial coordinates $\tilde{z}$ and $\tilde{x}$ are scaled as:

[TABLE]

where $L$ is the characteristic length along $\tilde{x}$ (for e.g., the wavelength of a typical interfacial perturbation). We further assume that $\epsilon=H/L\ll 1$, i.e. a small aspect ratio film, which would permit us to employ the lubrication approximation. The velocity components and the fluid pressure are scaled as follows:

[TABLE]

The scalings in (10) reflect a balance between viscous stresses and the disjoining pressure gradient caused due to van der Waals forces. We note that our choice of scalings is different from the ones used by Matar (2002) (extensional viscous stresses $\sim$ Marangoni stresses), and would be more apt for ultrathin films ($\leq 100~{}nm$) (Edwards & Oron, 1995; Ting et al., 1984; Berg, 2010). The natural choice for the characteristic time scale is the ratio of the characteristic scales of velocity and length. Thus, the dimensionless time variable is expressed as:

[TABLE]

Since surface properties depend on surfactant concentration, suitable scalings are needed for the concentration, surface tension, and surface viscosity. The surfactant concentration is scaled in terms of the reference concentration $\tilde{\Gamma}_{0}$:

[TABLE]

This reference concentration could, for instance, be an initial equilibrium surface concentration, or a critical concentration of colloidal surface active agents in a jammed-state interface ($\tilde{\Gamma}_{max}$). The actual choice of $\tilde{\Gamma}_{0}$ will have implications on the model used for variable surface viscosity effects which will be discussed later. The surface tension $\tilde{\sigma}$ is scaled with respect to its clean interface value (i.e. in the limit $\tilde{\Gamma}_{0}(x,t)\to 0)$ whereas the surface viscosity, $\tilde{\eta}$, is scaled with $\tilde{\eta_{0}}$, the choice of which depends on the type of variation used for $\tilde{\eta}$ (see Appendix B),

[TABLE]

Several dimensionless characteristic numbers appear upon rewriting 1-8 using the scalings mentioned above. These are summarized in table 2 along with their typical order of magnitude estimates based on physical values given in table 1. Among these, the parameters representing dominant effects are Bo (non-dimensional surface viscosity) and $\mathbb{D}$ (non-dimensional surface tension). Further, the Marangoni number M and the surface viscosity gradient parameter, $\beta$, represent the variation of surface tension and surface viscosity with surfactant concentration. A linear variation with concentration is assumed for surface tension:

[TABLE]

For the surface viscosity, a linear as well as a nonlinear model is used:

[TABLE]

We use the linear model (15) for an interface with dilute surfactants, thus the reference concentration, $\tilde{\Gamma}_{0}$ has a value in the dilute limit. For the nonlinear model (16), the jammed state concentration value, $\tilde{\Gamma}_{max}$ of the system would intuitively be the suitable non-dimensionalizing factor for $\tilde{\Gamma}$ as the surface viscosity is expected to diverge in the jammed state, analogous to bulk viscosity (Maki & Kumar, 2011; Krieger & Dougherty, 1959; Quemada, 1977). Recent studies in polymer blends also reveal such a dependence of surface viscosity on nano-particles concentration (Vandebril et al., 2010). Note that the nonlinear phenomenological model (NPM) simplifies to a linear model in the limit of dilute surfactant concentration. Accordingly, the exponent $\alpha$ is related to the surface viscosity gradient parameter $\beta$ by the relation,

[TABLE]

where $\Gamma^{(l)}_{0}$ and $\Gamma^{(nl)}_{0}$ denote the respective non-dimensional base state values in the linear and nonlinear models. Further details on the validity and derivation of the surface viscosity models are given in Appendix B. We also note that the NPM for surface viscosity reduces to a constant surface viscosity when $\Gamma\to 0$, in accordance with earlier studies on surface viscosity by Scriven (1960).

The disjoining pressure number, $\mathbb{D}$ in Table 2, defined as $\mathbb{D}=A/6\pi\sigma_{0}\epsilon L^{2}$ has been deliberately scaled to contain $\epsilon$, so that capillary effects are retained in the normal stress boundary condition, as also described in Appendix A.1. It represents the ratio of intermolecular van der Waals interactions to surface tension forces. A disjoining pressure interpretation of the capillary effects is apt as it represents the overall wetting characteristics of the thin-film system, as has been done in an earlier work (Alleborn & Raszillier, 2004). With respect to Marangoni and surface viscosity effects, we restrict ourselves to the default regime examined in prior works (DeWit et al., 1994; Matar, 2002) in which both Marangoni effects and surface viscosity only appear in the first order correction of the tangential stress balance condition. Accordingly, if M is an $O(1)$ quantity, a weak Marangoni limit is defined by rescaling it as $\textit{M}=\epsilon^{2}\hat{\textit{M}}$, where the new O(1) parameter is $\hat{\textit{M}}$. Then, Marangoni effects would drop out in the lubrication limit when $O(\epsilon^{2})$ terms are ignored. Following the definition of $\sigma$ in (14), the weak Marangoni limit would imply, at leading order,

[TABLE]

For the surface viscosity, $\eta$ is retained as an O(1) quantity, but it would still appear only in the first (and higher) order corrections to the stress balance conditions as shown by Matar (2002). To explore the strong surface viscosity regime, $\eta$ may be rescaled such that it is retained at leading order, and the boundary conditions are self-consistent. However, as $\epsilon$ (aspect ratio) can be made arbitrarily small, presenting the surface viscosity terms at $O(\epsilon)$ stress balance conditions (as in the present formulation) could still be deemed acceptable despite the diverging nature of NPM.

2.3 Evolution equations

With the scalings presented in §2.2, the governing equations (1-3) reduce to the following when only the zeroth (leading) order terms are considered, denoted with the superscript $(0)$

[TABLE]

Using the symmetry conditions (8) at $z=0$ gives:

[TABLE]

where $c(x,t)$ is still an unknown function independent of $z$. The leading order normal stress balance (see Appendix A.1) is written as:

[TABLE]

Integrating (20) and comparing it with (24), and substituting (22-23), we obtain:

[TABLE]

The leading order tangential stress balance (see Appendix A.2), on the other hand, simplifies to: (using 18)

[TABLE]

This is also consistent with (22). The kinematic condition (7) and the surfactant transport equation (4) provide two nonlinear evolution equations for $h^{(0)}$ and $\Gamma^{(0)}$, with $c(x,t)$ still undetermined. As discussed by (DeWit et al., 1994) and (Matar, 2002), the system can be completed only through a third equation obtained by considering the first order corrections of (19) and (26) (refer Appendix A.2, A.3). To these, we substitute (15) or (16) to account for variable viscosity. The final system of three nonlinear evolution equations read:

[TABLE]

[TABLE]

or,

[TABLE]

Here, $\Rey~{}\&~{}\Pen$ are as defined in Table 2. The superscript (0) has been omitted in (27 - 2.3) and the remainder of the paper for easy readability. Equations (27 - 2.3) reduce to those derived by DeWit et al. (1994) in the limit $\textit{Bo}=0$, and to those of Matar (2002) when we set $\beta=0$, in the absence of surfactant solubility.

3 Linear Stability Analysis

We linearize the system of three nonlinear equations (27), (28), and (2.3) or (2.3) as the case be, by perturbing the dependent variables about a uniform base state as follows:

[TABLE]

Here $s$ is the growth rate and $k$ is the wavenumber of the perturbation. Using standard linear stability analysis, the following dispersion relations are obtained for the linear and nonlinear surface viscosity models:

[TABLE]

where, $\Gamma_{0}^{(l)}$ and $\Gamma_{0}^{(nl)}$ are the non-dimensional base states for surfactant concentration in the linear model and NPM respectively, as has already been discussed in deriving (17). The above dispersion relations are in agreement with earlier studies (Erneux & Davis, 1993; Matar, 2002; DeWit et al., 1994) when constant surface viscosity is considered. The main novelty of the present work is a careful study of variable surface viscosity and these enter through the variables $\beta$ and $\alpha$ in (3) and (3) respectively.

3.1 Parametric studies

The dispersion curves for the linear and NPM models are plotted in Fig. 2 and 3 for a range of parameter values chosen from table 2. All the curves exhibit a fastest growing mode, $k_{max}$ and a cut-off wavenumber, $k_{c}$,

[TABLE]

above which the film is stable. Fig. 2 and 2 show the linear stability characteristics with variation in the parameter $\mathbb{D}$ for the two models. We set $\Gamma_{0}^{(nl)}=0.9$ for the dispersion curves of NPM. Intuitively, a decrease in $\mathbb{D}$ results in greater opposition from surface tension towards film breakup, and a smaller driving force from intermolecular van der Waals attractions. Accordingly, it is seen that the maximum growth rates reduces with decreasing $\mathbb{D}$. Figures 2 - 2 are the linear stability curves for different values of Bo and M for both the models. Both surface viscosity and the solutal Marangoni effect act as film stabilizers, as inferred from the growth rate plots. Note that, varying Bo has drastic effects in case of NPM (fig. 2).

Figures 3-3 are plots of the dispersion relation for the other system parameters: $\Rey$ and $\Pen$ for both the models. The nature of the dispersion curves plotted for these parameters are also intuitive. Physically, it may be understood that a greater value of $\Rey$ implies more inertial forces relative to viscous forces. Though both the parameters act as retarding effects towards growth of instability, the contribution of the former is more significant, leading to smaller growth rates at higher $\Rey$. As regards to $\Pen$, a higher surface diffusivity of the surfactants (lower values of $\Pen$) facilitates higher gradients of surfactant-concentrations to destabilize the film evolution which can be clearly observed in fig. 3. The linear stability curves for NPM, governed by (3) may seem qualitatively similar to those of linear model, but has a distinctive boundary layer structure (for $k\to O\left(1-\Gamma\right)^{\alpha}$) which is very different from that of the linear model. This is discussed in §3.2 in detail. Both the models follow (36). Equations (3) and (3) become identical when the parameters are chosen to satisfy (17) in the limit of $\tilde{\Gamma}\ll\tilde{\Gamma}_{max}$.

3.2 ‘Jammed limit’ analysis

In this section, we show that it is possible to obtain a simple perturbation solution to the dispersion relation (3) in the limit of large surfactant concentrations, i.e. for $\Gamma_{0}^{(nl)}\rightarrow 1$. In this limit, the surfactant concentration approaches the jamming limit and mobility of surfactant molecules reduces to zero. It will be shown later (see §4.3) that the tangential velocity indeed reduces to zero as $\Gamma_{0}^{(nl)}\rightarrow 1$. This limit of (3) is distinctly different from that of the linear model in (3) which is valid in the limit of small $\Gamma_{0}^{(nl)}$. It is shown below that the nature of the cubic equation changes in the jamming limit which is not possible using the linear model. To proceed with the perturbation analysis, it is convenient to define a small parameter, $\delta=\left(1-\Gamma_{0}^{(nl)}\right)^{\alpha}$. As $\Gamma_{0}^{(nl)}\rightarrow 1$, $\delta\rightarrow 0$. Equation (3) can then be rewritten in terms of $\delta$:

[TABLE]

where the coefficients $a_{1}$ to $a_{5}$ are functions of $k$ and can be obtained by comparing (37) with (3):

[TABLE]

In the limit $\delta\rightarrow 0$, Equation (37) is a standard singular perturbation problem for growth rate $s$. Fig. 4 shows plots of the sole positive root of (37) for different values of $\delta$. It is apparent that the growth rate $s$ scales with $\delta$ and singular root appears in the region $k\leq\delta$. We therefore define a rescaled growth parameter $Y=s/\delta$, so that equation (37) now becomes

[TABLE]

A standard expansion in $\delta$ in the form $Y^{(o)}=Y_{0}^{(o)}+\delta Y_{1}^{(o)}+...$ results in two regular roots in the outer region, i.e. $k\gg\delta$. Focussing on the root which is unstable, we obtain

[TABLE]

To obtain the singular root of the inner region $0\leq k\leq\delta$ which has a boundary layer-like structure, we define a rescaled variable $X=k/\delta$, where $X\sim\it{O}(1)$. The solution to the inner region in terms of this rescaled variable at leading order is:

[TABLE]

where $\displaystyle b_{2}=\frac{2Bo}{\Rey}X^{2};~{}~{}b_{3}=\left(\frac{\mathbb{D}^{-1}\textit{M}\Gamma_{0}}{\Rey}-\frac{3}{\Rey}\right)X^{2}$. Figure 4 shows a comparison of inner solution, $Y_{0}^{(i)}$ (dot-dashed line), outer solution, $Y_{0}^{(o)}$ (dashed line) and the numerical root, $Y$ (solid line) obtained from (43). Examining (46), it can be shown that the necessary condition for stable growth rates ($Y_{0}^{(i)}<0$) is $b_{3}>0$. Hence we arrive at a necessary condition for obtaining stable modes in the space $0\leq k\leq\delta$:

[TABLE]

We can therefore define a critical concentration, $\Gamma_{c}=3/\mathbb{D}^{-1}\textit{M}$. Any value of $\Gamma_{0}^{(nl)}>\Gamma_{c}$ should result in inherently stable films. Also note that since $\Gamma_{0}^{(nl)}$ is, by definition, bounded by: $0<\Gamma_{0}^{(nl)}<1$, we can hypothesize that no amount of surfactants can stabilize the system if the product of the system parameters, $\mathbb{D}^{-1}\textit{M}<3$. Examining the system parameters in real physical systems, we find that $\mathbb{D}^{-1}\textit{M}\sim\it{O}\left(10^{-1}\right)$ in real systems suggesting that though foams can be stable for long durations, indefinite stability is not guaranteed. This is consistent with the observation in Kloek et al. (2001) who studied bubble dissolution with complex interfacial and bulk rheologies.

It is easy to derive a simple relationship between the amplitude of tangential velocity at the interface, $\hat{c}$ and height perturbation, $\hat{h}$ using the kinematic condition (27):

[TABLE]

The above relationship suggests that the tangential velocity at the interface, $c(x,t)$, has a phase-difference with the interface height, $h(x,t)$ since $\hat{c}$ is a complex number. Since interface height is symmetric about the rupture location, tangential velocity is anti-symmetric. This result is unsurprising since as the interface height reduces towards rupture, fluid is driven away from the rupture location. An important outcome of the above analysis is that growth rate scales with the small parameter $\delta$ in the jamming limit. This dependence can be written in the simple form

[TABLE]

where we have used the definition of $\delta$ and functional dependence on other parameters, $F$, is not reproduced for simplicity. Substituting $s$ into (48), we get

[TABLE]

As $\Gamma_{0}^{(nl)}\rightarrow 1$, it is clear that $c(x,t)\rightarrow 0$, thus rendering the interface immobile in the jammed limit. Both the anti-symmetric nature of tangential velocity and approach to immobile limit are confirmed through full numerical solutions of equations (27), (28) and (2.3) in section §4.

4 Nonlinear evolution

4.1 Numerical methods

The nonlinear evolution and dynamics of the free-film was studied by numerically solving the set of partial differential equations (27) - (2.3) for a range of parameter values. A Fourier pseudo spectral (FPS) method is used to solve these PDEs, since it is most suited to handle periodic boundary conditions. Discretization of the spatial derivatives in the Fourier domain reduces the set of PDEs to ordinary differential equations (Trefethen, 2000). Time marching is done using an implicit Adams-Moulton method (trapezoidal rule) using MATLAB function fsolve to solve the linearized equations in fourier space. An adaptive time stepping, $\Delta t\propto\kappa^{-1}$(curvature) is also implemented to capture the dynamics at later times, when the evolution proceeds much faster. The initial conditions used in our simulations are as follows:

[TABLE]

where $\lambda=2\pi/k_{max}$ is the fastest growing wavelength, determined from linear stability analysis. The perturbation amplitude for surfactant concentration, $\Delta\hat{\Gamma}$ is set using:

[TABLE]

which is obtained from linear stability analysis. Here, we use $\Gamma_{0}^{(l)}=2$ and $\Gamma_{0}^{(nl)}\in[0,0.95]$ for the linear and nonlinear models respectively. Further, we set $h_{0}=0.5,~{}\Delta\hat{h}=0.01$, with $s_{max}$ being the maximum value of the growth rate predicted by the linear theory, at $k=k_{max}$. The time steps used in the initial stages are of the order $\Delta t\sim O(10^{-1})$ and reduces to $\Delta t\sim O(10^{-3})$ at the latter stages of our numerical simulations. The spatial domain is discretised into 256 nodal points ($\Delta x=1/256$) in all the simulations, ensuring grid-independence of the results. While exploring self-similarity, the grid size and time step values are further refined. We set $\Delta x=1/4096$ and use adaptive time-stepping with lower end values of the order $\Delta t\sim O(10^{-5})$. We assume the film to have ruptured when the minimum height of the film, $h_{min}\leq 10^{-4}$.

4.2 ‘Dilute limit’ analysis

4.2.1 Evolution and parametric studies

The nonlinear dynamics using the linear model (15) for surface viscosity in a film with dilute concentration of surfactants are discussed first. We remind that this is a special limit of NPM as discussed in §2.2. Figure 5 shows spatiotemporal evolution profiles for (a) film height $h$, (b) concentration $\Gamma$, (c) velocity at the interface, $c$ for a typical case: $\mathbb{D}=0.02,~{}\textit{M}=10^{-3},~{}\beta=10^{-3},~{}\textit{Bo}=1,~{}\Pen=20,~{}\Rey=2$. The profiles have been plotted on spatial coordinates over a domain length $\lambda$. Given that the boundary conditions are periodic, successive crests/troughs in the profiles are separated by $\lambda$, consistent with expectations of the spinodal dewetting mechanism. The numerically obtained growth rate is also validated with the predictions of linear theory, as is shown in fig. 5. At short times, the linear theory predictions are in agreement with the nonlinear solutions, but deviate at later times closer to rupture, when the nonlinear terms become significant. Parametric studies on film evolution in the ‘dilute’ regime reveal that qualitative features of evolution follow the trends observed in earlier studies (Matar, 2002; Erneux & Davis, 1993). The rapid formation of sharp ruptures in the profiles of $h$ and $\Gamma$ indicate the possibility of self-similar solutions close to rupture.

Figure 6 shows the effect of surface viscosity (Bo) and its linear concentration dependence (represented by $\beta$) on the interface profile and film evolution kinetics. For smaller values of Bo, the surface viscosity contribution is lesser, and a sharp, pointed rupture is observed (Erneux & Davis, 1993; DeWit et al., 1994). Increasing Bo flattens the cusps formed in the vicinity of rupture. Surface viscosity may be interpreted as a diffusivity of surface velocity or momentum, hence it neutralizes velocity gradients at the surface. Therefore, there is greater advection of the liquid onto either side near the rupture location resulting in a flatter rupture profile. On the other hand, increasing the parameter $\beta$ produces a more cusp-like rupture, as shown in figure 6. This is understandable since surface viscosity from (15) at rupture location reduces to simply becoming $\eta|_{x=x_{r}}\approx 1-\beta$ since $\Gamma\ll 1$ near rupture. Therefore increasing $\beta$ reduces the local value of surface viscosity near rupture location. It has to be noted that the maximum value of $\beta$ cannot exceed unity as a result of positivity of $\eta$. The upper limit of $\beta=1$ causes the non-dimensional surface viscosity to become identical to the non-dimensional surface concentration, i.e. $\eta|_{\beta=1}=\Gamma$. Since $\Gamma\rightarrow 0$ at the rupture location, $\eta$ too vanishes at the rupture location. This causes the film profile to bear resemblance to that of a zero surface viscosity case as is clearly evident with the cusp-like profile at $\beta=1$ in figure 6. This profile matches very well with the case of zero surface viscosity, $\textit{Bo}=0$ as shown in figure 6. This is further discussed in §4.2.2. The rupture times appear to be of the same order of magnitude for a wide range of values of $\beta$ as seen in Fig. 6.

4.2.2 Self-similar solutions

Next, we investigate the presence of self-similar behaviour in free film profiles in the vicinity of rupture location. Equations (27)-(2.3) do not lead to perfect self-similar solutions, but contain a few terms with weak time dependence. We assume solutions for the unknowns $h(x,t)$, $\Gamma(x,t)$ and $c(x,t)$ in the form:

[TABLE]

where $\chi$ is the similarity variable and $\tau=t_{r}-t$, $x_{r}$ and $t_{r}$ being the spatial location and time of rupture, respectively. Substituting (4.2.2) into (27)-(2.3) and subsequent simplification yields:

[TABLE]

To obtain scalings followed by $h$, $\Gamma$ and $c$, we follow the approach used by Matar (2002) and Vaynblat et al. (2001). These scalings allow us to extract the exponents $i,j,k,l$ without knowledge of rupture time a priori. The natural logarithm of the curvature $(h^{min}_{xx})$, velocity $(c_{max})$ and surfactant concentration on the surface $(\Gamma_{min})$ at $x_{r}$ are plotted against the natural logarithm of film thickness $(h_{min})$ at $x_{r}$. Three algebraic equations in terms of $i,j,k,l$ can be generated by comparing the slopes of $(h^{min}_{xx})$, $(c_{max})$ and $(\Gamma_{min})$, shown in figure 7, with relations obtained from equation (4.2.2). A number of choices are possible for the fourth equation. For self-similarity, we require the exponent of $\tau$ in all the terms of equations (56)-(4.2.2) to vanish. An algebraic equation can be written by setting the powers of $\tau$ to zero in the above equations and each such equation can be solved along with the three algebraic equations obtained from numerical simulation. This leads to multiple solutions for $i,j,k,l$. Most of the solutions can be discarded on physical grounds. For example, since $h(x,t),\Gamma(x,t)\rightarrow 0$ during rupture as evident in figure 5, we require $i>0$ and $l>0$. Similarly, tangential velocity, $c(x,t)$ rapidly increases on either side of rupture location since fluid is rapidly driven away from the rupture location as interface height reduces to zero. This requires $k<0$ in (4.2.2). We can also argue on similar grounds that $j>0$ in the definition of similarity variable $\chi$. Of the many solutions that are possible for $i,j,k,l$, a large number of solutions are ruled out by the above constraints on these variables. It has to be noted that the above procedure does not eliminate $\tau$ dependence in all the terms since the equations do not admit perfect self-similarity as was noted in earlier works. The terms which retain $\tau$ dependence in (56)-(4.2.2) become vanishingly small in the limit of $\tau\rightarrow 0$. This ensures that the system evolves towards a self-similar state as one approaches the time of rupture. This analysis gives physical insights into the dominant forces in operation at the onset of rupture, essentially among the terms with vanishing $\tau$-dependence.

Since the primary aim of this work is to shed light on effects of variable surface viscosity, we restrict analysis in this section only to the case of varying $\beta$. For comparison, we also report the case of zero surface viscosity. Recall that surface viscosity varies with concentration in the linear model by the relation $\eta=1+\beta(\Gamma-1)$. We consider four test cases for the study as governed by equations (27)-(2.3): (i) $Bo=0$ which corresponds to the case of zero surface viscosity as reported by Vaynblat et al. (2001); (ii) $\beta=1$ where surface viscosity and concentration assume equivalence; (iii) $\beta=0.5$ which corresponds to the case of variable surface viscosity and is a canonical case for $0<\beta<1$; and (iv) $\beta=0$ which corresponds to the case of constant surface viscosity as reported by Matar (2002).

Results for case (ii) are not shown here but will be constrasted with other cases. Further $\beta$ was varied from $0.1$ to $0.9$ in steps of $0.1$ to study the effect of $\beta$ on similarity behaviour. It is found that scaling exponents remained unaffected with $\beta$ and the reason for the same is explained below.

Results for $\beta=0.5$ are shown in figures 7 and 7. Following the procedure as explained above, the following set of exponents are found: $i=l=6/13,j=4/13,k=-9/13$. These exponents values are different when compared to the constant viscosity case studied by Matar (2002). In our case, this occurs for $\beta=0$ and yields the exponents $i=l=1/2,j=1/4,k=-3/4$ which agree with the values reported in Matar (2002). The exponents with variable surface viscosity are only slightly different from the constant viscosity case, nevertheless they reveal that self-similar behaviour is altered with variable surface viscosity. For case (ii) with $\beta=1$, surface viscosity all along the interface reduces to the expression, $\eta=\Gamma$. Since surfactant concentration at rupture location becomes vanishingly small, surface viscosity too becomes negligible in the vicinity of rupture. The following exponents are found for this case: $i=l=1/3,j=1/2,k=-1/2$ which is identical to the case with zero surface viscosity as obtained in Vaynblat et al. (2001) and corresponds to case (i) in our study with $Bo=0$. It has to be noted that surface viscosity has a finite value away from the rupture point and this leads to deviation of the interface shape from that reported in Vaynblat et al. (2001) in the far-field. A summary of all the cases with relevant scaling exponents in each case is given in figure 8.

It may be intuitive to conclude that if a concentration-dependent surface viscosity is present, the self-similar scalings will revert to those obtained for $Bo=0$, since $\Gamma\rightarrow 0$ near the rupture location. Though the same is suggested by Matar (2002), our analysis shows that this is not entirely true and obvious, for there is a finite surfactant concentration gradient at the rupture location, even though $\Gamma\rightarrow 0$. Our self-similarity analysis for the variable surface viscosity model results in scalings that lie between the regimes $Bo=0$ (zero surface viscosity) and $\beta=0$ (constant surface viscosity). The reason underlying this result is further illustrated by fig. 7 showing the temporal evolution of the surface viscosity at the rupture point for $\beta=0$ and $\beta=0.1$ (concentration driven surface viscosity). Even though the surface viscosity starts with a higher value than the constant viscosity case, the effect decreases near the rupture point for $\beta>0$. The interface develops towards a more cusp-like profile (closer to profiles for $Bo=0$) decreasing the surface viscosity even below than that of the constant surface viscosity case. The same may also be understood by comparing fig. 7 that exhibits a flatter profile with fig. 7 that exhibits a profile with sharper cusps. For the extreme case of $\beta=1$, surface viscosity effects are suppressed completely at the rupture point and a cusp-like profile of a clean interface, $\Gamma=0$ (or of a particle-laden interface with no surface viscosity, $Bo=0$) is recovered. This also means that however large the initial surface viscosity we choose, we can get a cusp-like profile if we have $\beta=1$. Self-similarity is not seen in simulations with NPM in the jammed limit due to the nature of the model.

The scaling exponents discussed above and given in figure 8 also reveal the dominant balance of forces in the vicinity of rupture. This can be determined by substituting the exponents into the similarity equations (56-4.2.2) for various cases. With $Bo=0$ (case (i)) and $\beta=1$ (case (ii)), the dominant balance is between inertia, van der Waals and viscous forces and the interface forms a cusp during rupture. In the case of $0\leq\beta<1$ (cases (iii) and (iv)), the dominant balance is between inertia, van der Waals and surface viscous forces. Viscous forces become sub-dominant as one approaches the rupture point. In all the above cases, capillary and Marangoni forces play negligible role near rupture consistent with observations in earlier studies. The above dominant balance was further verified by comparing the magnitude of each term in the evolution equation (2.3).

It is interesting to note that for cases (iii) and (iv), another set of consistent scaling exponents are also found. For case (iii), these are $i=l=3/4,j=1/2,k=-9/8$ and for case (iv), we get $i=l=1,j=1/2,k=-3/2$. The latter set of exponents were not reported by Matar (2002) perhaps ignored due to the fact that a large value of Peclet number was used in their study. As a result, the algebraic equation arising out of surfactant diffusion term, $1-2j=0$, in equation (57) was ignored. Nevertheless, it is found that the above scaling exponents do not affect the dominant balance of forces near rupture. It is likely that these new exponents will change the self-similar profiles reported in figure 7, but a detailed investigation on this will be considered in a future study.

In the next section, we investigate nonlinear evolution of the interface at high concentrations using nonlinear surface viscosity model.

4.3 ‘Jammed’ limit analysis

The rigidification of the interface in the limit of jammed state concentration of surfactants is perhaps the most interesting aspect of this work. The notion that films at jammed limit are very stable as suggested by linear stability analysis can be visualised better in nonlinear simulations using the nonlinear phenomenological model (NPM) for surface viscosity given by (16). Figure 9 illustrates the variation of $\eta$ for NPM (16) for various values of $\alpha$. Higher values of $\alpha$ causes $\eta$ to diverge at lower values of $\Gamma$. To compare film profiles for both the models in dilute limit, the values of $\alpha,~{}\beta,~{}\tilde{\Gamma}_{0}$ and $\tilde{\Gamma}_{max}$ are chosen so as to satisfy the relation (17). Both the profiles appear to be in agreement in this limit, as shown in Appendix B, figure 12. Having validated the code using NPM, we now study the effect of high surfactant concentration for which NPM is ideally suited. The height profile and the surfactant concentration profile for the particular value $\Gamma_{0}^{(nl)}=0.8$ and $\alpha=2$ of NPM are shown in fig. 10. The film stiffens and tends to be immobile all along the interface except in a narrow region at the centre of the domain where height and concentration profiles are at a minimum. Unlike in earlier cases, the film approaches rupture in a very narrow region. This evolution quickly deviates from a linear evolution (not shown) which mainly follows a smooth amplifying sinusoidal structure. A better way to demonstrate the difference in film shape in NPM is to plot the film height profile for different initial surfactant concentrations. Figure 10 compares film shapes at rupture for three different initial concentrations: the dilute limit at $\Gamma_{0}^{(nl)}=0.1$, an intermediate limit at $\Gamma_{0}^{(nl)}=0.6$ and a near-jammed state limit at $\Gamma_{0}^{(nl)}=0.9$. The profiles in the dilute limit exhibit a large radius of curvature and is similar to what is in earlier sections. In the near-jammed limit, rupture occurs over a very narrow region within the interface. It has to be noted that the interface flattens at the point of rupture but has a very small radius of curvature. The case of $\Gamma_{0}^{(nl)}=1$ is a singular limit and cannot be solved numerically. This is consistent with the perturbation solution carried out in section §3.2. Moreover, rupture time diverges with concentration as shown in figure 10.

A particularly interesting result is the evolution of tangential velocity, $c(x,t)$, on the interface with varying concentration as shown in figure 10. It is clear that the maximum tangential velocities decreases and asymptotes towards zero with increasing concentration, suggesting a transition from a mobile interface for $\Gamma_{0}^{(nl)}=0$ to an immobile interface at jammed state. Such a surface viscosity-driven transition from partially mobile to immobile interfaces have also been reported by Danov et al. (1999). For drainage of a thin liquid film between two gas bubbles, interfacial resistance arising from surface viscosity have been reported. Consequently, the kinetics of film break-up is also found to slow down drastically with increase in concentration, and diverges to infinity as $\Gamma_{0}^{(nl)}\to 1$, suggesting the formation of a solid-like rigid interface which does not deform.

A key result of the paper is summarised in figure 10 which shows the transition of the interface from mobile to immobile limit associated with a divergence of rupture time as surfactant concentration approaches the jamming limit. This result offers a simply way to model fluid-fluid systems with variable surface viscosity effect to describe a smooth transition from free-slip regime to no-slip region such as in the case of a surfactant-laden sedimenting drop (see Leal (2007)). Recent experiments by Mayarani et al. (2018) on an evaporating droplet containing polystyrene and soft micro-gel particles reveal that micro-gel particles form an immobile interface. As a result, mean-squared displacement of adsorbed polystyrene particles reveal a transition from liquid-like to solid-like behaviour. This effect can potentially be modelled using a nonlinear surface viscosity model used in the current study.

5 Summary and Discussions

We have formulated a unified model that includes various interfacial effects due to the presence of surface active agents in a thin free film. Apart from intermolecular forces and Marangoni effect, it is shown that surface viscosity effects play a key role in determining the stability of a thin film. It is well known that Marangoni effects due to surfactant driven surface tension gradients stabilizes a thin film, but this stabilization does not explain the long shelf-life of certain surfactant-stabilized emulsions and Pickering emulsions. Marangoni stabilization only tends to delaying rupture by a small factor. We have shown that concentration-dependent surface viscosity mediated stabilization is a missing element in most earlier studies and offers a theoretical tool to analyse long shelf-life of such emulsions. As a canonical problem, stability of a free film is studied and the theory can be easily extended to bounded films and other interfacial problems. To better illustrate the dependence of surface viscosity on surfactant concentration, two distinct phenomenological models are employed. In the linear model, surface viscosity varies linearly with surfactant concentration and is expected to be suited for dilute limits. The nonlinear model allows us to probe the role of surface viscosity profile in greater detail. In the nonlinear model (NPM), surface viscosity varies nonlinearly with concentration such that it diverges at a critical concentration which is termed as the ‘jammed limit’. The two models are shown to concur in the dilute limit and this lends validity to the nonlinear model at all concentrations.

In the linear model, surface viscosity stabilizes the film with increasing Boussinesq number but does not affect the cut-off wavenumber. On the other hand, NPM drastically alters the stability characteristics of the thin film. The growth rate is found to scale with $\delta=(1-\Gamma_{0}^{(nl)})^{\alpha}$ where $0<\Gamma_{0}^{(nl)}<1$ is the non-dimensional concentration scaled in terms of the ‘jammed limit’ concentration, $\Gamma_{max}$. The parameter $\delta$ can thus be made arbitrarily small near the jamming limit. This leads to a dramatic stabilization of the thin film and significantly delays rupture in the nonlinear simulations. Using standard perturbation techniques, analysis of the dispersion relation yields a simple criteria for film stability. It is found that when $\Gamma_{0}^{(nl)}>3\mathbb{D}/\textit{M}$, instability of the thin film is completely suppressed. This suggests that surfactant concentration can be arbitrarily increased beyond a critical value ($=3\mathbb{D}/\textit{M}$) to achieve permanent stabilization of the film. But in reality, such concentrations are difficult to achieve for realistic parameter values of M and $\mathbb{D}$. A key finding of the present paper is that rupture times can be arbitrarily increased by varying surfactant concentration. This stabilization is achieved via surface viscosity effects which increase nonlinearly with surfactant concentration. Large surface viscosities render the interface immobile which helps counteract the rupture process by slowing down film evolution.

The findings of this work have relevance to many other experimental and theoretical works. Experiments with rising drop/sedimenting bubbles (Leal, 2007) reveal that surfactants enhance the drag force on a drop/bubble. This effect has been modelled theoretically using a spherical cap approximation (Leal, 2007; Sadhal & Johnson, 1983) wherein surfactants are assumed to be in a jammed state on the leeward side of the drop/bubble and are assumed to be absent on the rest of the surface. A no-slip condition is employed on the spherical cap and a stress-free condition on the rest of the surface. NPM employed here leads to a smooth transition from the no-slip to the stress-free condition on the interface and prevents a jump in the interface conditions as was employed by Sadhal & Johnson (1983). Our results are also consistent with Danov et al. (1999) who found that interface mobility reduces with increasing surfactant concentration. This has direct significance to the stabilty of Pickering emulsions (Wasan et al., 1992; Langevin, 2000). There is sufficient anecdotal evidence in literature on Pickering emulsions (Anjali & Basavaraj, 2018) which show that removing adsorbed particles (by changing the pH of the aqueous phase) reduces emulsion stability dramatically. The treatment of the surface active agents as ‘surfactants’ also allows one to look at proteins in food systems as surface active agents (Murray & Ettelaie, 2004). A few studies have found that truly stable foams can only be formed when coagulated proteins form a rigid solid network on the interface (Walstra, 2003; Martin et al., 2002). Our work serves as a mathematical deduction of the various phenomena described above by modelling the surface viscosity.

A number of interesting features remain unexplored. Though nonlinear surface viscosity effects sucessfully explain the transion from mobile to immobile interfaces, they do not account for elasticity found in particle-laden interfaces. Liquid marbles (Aussillous & Quèrè, 2001) and particle-coated bubbles can sustain non-spherical shapes upon compression (Mahadevan et al., 2005). This could be due to anistropic surface tension or elasticity or both. Studies with particle-covered Landau-Levich dip-coating flows (Dixit & Homsy, 2013) reveal that elastic effects play an important role in explaining experimental power-law dependence between coating thickness and withdrawal speed. It is therefore imperative that future studies incorporate a combination of variable surface viscosity and elasticity to better predict characteristics of particle/surfactant-covered interfacial systems, especially when jamming is expected occur.

An earlier version of this work was presented by one of the authors (A.C.) at the AIChE Annual Meeting 2018 and A.C. thanks Dr. Ashutosh Ricchariya (LVPEI, Hyderabad) and Department of Biotechnology, Govt. of India for travel support. H.N.D and S.K. thank Pavan Pujar from NIT Surathkal who carried out preliminary analysis of this work. The authors gratefully acknowledge Prof. James J. Feng (UBC, Vancouver) and Dr. K. Badarinath (IIT Hyderabad) for their critical reading of the manuscript and constructive suggestions. A.C was supported by a doctoral fellowship from the Ministry of Human Resource Development, Govt. of India and H.N.D thanks partial support from an Early Career Research Award (ECR/2015/000086).

Appendix A Interfacial boundary conditions and $x$-momentum equation

A.1 Normal stress balance

For a 1D interface, the normal stress balance assumes the form:

[TABLE]

We resolve the $\nabla_{s}$ curvilinear gradient operator and the tensor products to get a PDE in cartesian coordinates. This is done analytically and validated by the open-source technical computing software, Mathematica. The dimensional equation turns out to be:

[TABLE]

Here the surface dilatational $(k^{s})$ and surface shear $(\mu^{s})$ viscosities are combined into a parameter $\eta$, since they occur in additive pairs in the equations. This is only observed for a 1D interface. Inserting relevant scalings (9-14) as discussed in the main text, we observe the pressure-term appears at leading order followed by the capillary-term at $O(\epsilon)$. Hence, the Disjoining pressure number $(\mathbb{D}^{-1})$ (or $\sigma$) is rescaled to retain capillary effects. We obtain the final simplified non-dimensional PDE, which reads:

[TABLE]

Considering only the leading order terms, we get (24) used in deriving the film evolution equations.

A.2 Tangential stress balance

For a 1D interface, the tangential stress balance assumes the form:

[TABLE]

The dimensional equation turns out to be:

[TABLE]

Inserting relevant scalings and simplifying, we obtain the final non-dimensional PDE:

[TABLE]

The above equation reduces to (26) when the leading order terms are considered. The first order corrections of this equation used to derive (2.3) or (2.3) can also be realized when $O(\epsilon^{2})$ terms are considered. As can be observed clearly, the surface viscosity enters the equations in the first order corrections of tangential stress boundary condition, and has highly nonlinear coupling with the height and various velocity gradient terms.

A.3 $x$-momentum equation

The non-dimensional $x$-momentum equation valid to $O(\epsilon^{2})$ take the form,

[TABLE]

Expanding all terms in powers of $\epsilon^{2}$, the leading order equation becomes

[TABLE]

and the first order equation at $O(\epsilon^{2})$ assumes the form

[TABLE]

where superscripts [math] and $1$ denote leading and first order correction terms respectively in the asymptotic expansion (1). Integrating (67) and using symmetry boundary condition at $z=0$ (see eqn. 8), we obtain an expression for $u^{(1)}_{z}$ which when compared to the first order correction of (A.2), gives the third nonlinear evolution equation (2.3) (or (2.3)).

Appendix B Notes on the validity of scalings for $\tilde{\eta}_{s}$

The linear model adds small corrections to the surface viscosity to exhibit weak relationship with surfactant concentration. We may consider it to be apt when the surfactants are sparsely distributed on the interface (dilute limit):

[TABLE]

where, $\tilde{\eta}_{0}$ is the surface viscosity at a surfactant concentration within dilute limit $\tilde{\Gamma}=\tilde{\Gamma}_{0}$. The definition of $\beta$ as given in Table 2 can also be realised in the above equation. The surface viscosity for NPM is defined as,

[TABLE]

which diverges at $\tilde{\Gamma}\to\tilde{\Gamma}_{max}$ and reduces to a linear model when expanded in the limit $\tilde{\Gamma}\ll\tilde{\Gamma}_{max}$:

[TABLE]

Here, $\tilde{\eta}_{ref}$ is the surface viscosity for a clean interface ($\Gamma\rightarrow 0$). From (B), we also get the surface viscosity for a clean interface as,

[TABLE]

Comparing (B) and (70), we get an equation to relate the model parameters at a dilute concentration. This is visualized in fig. 11. Another necessary condition is satisfied by equating the slopes (or derivatives of (B) and (70)):

[TABLE]

Using (71) in the above relation, we obtain the equivalence condition (17). Note the advantage of defining the linear model for an arbitrary $\tilde{\Gamma}_{0}$ and not for $\tilde{\Gamma}_{0}=0$ (clean interface). This definition renders the linear model valid about any value of reference surfactant concentration, $\tilde{\Gamma}_{0}$ when $\beta$ at that concentration is known. Further, to relate the linear model with NPM, (72) will change accordingly, as now the NPM cannot be expanded in the limit of $\tilde{\Gamma}\ll\tilde{\Gamma}_{max}$.

Having estabilished the criteria for comparing the linear and nonlinear models in the dilute limit, we show that this comparison is robust and yields identical film profiles even in the nonlinear regime as shown in figure 12 which shows the comparison of height profiles from the two models at the instance of rupture.