On the role of variable surface viscosity in free viscous films

On the role of variable surface viscosity in free viscous films

Anjishnu Choudhury1, Venkatesh Kumar Paidi2, Sreeram K. Kalpathy2, Current address: Department of Chemical Engineering and Materials Science, University of Minnesota, Minneapolis, MN, USAEmail address for correspondence: sreeram@iitm.ac.in    Harish N. Dixit1 Email address for correspondence: hdixit@iith.ac.in
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. Surface tension varies linearly with surfactant concentration and two phenomenological models are used for variation of surface viscosity with concentration. In the first model suited for dilute concentration limit, surface viscosity varies linearly with surfactant concentration, while in the second model, surface viscosity varies nonlinearly with concentration and diverges at a critical concentration which is termed as the ‘jamming’ limit. Linear stability analysis is carried out with both surface viscosity models and 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 is a sufficient criteria for instability of the system, where is a normalized initial surface particle concentration, is the disjoining pressure number and 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 consistent with experimental observations. It is shown that rupture times can be arbitrarily increased by tuning initial surfactant concentration on the interface offering a fluid dynamical route to stabilization of thin films. 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.

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 () 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.

Figure 1: Schematic of the problem geometry. A free-film either arises in a typical soap bubble in air or between two bubbles or droplets in an emulsion. Bubbles deform as they approach each other forming a flat free film between them.

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 (dimensional), and the liquid has a bulk viscosity and density represented by and (both dimensional) respectively. The free surface is at a height . Further, we let denote the dimensional mean surface tension of the liquid film when the surface concentration of surfactants is maintained at at a fixed reference concentration . Table 1 gives some realistic dimensional values for the various physical parameters used in our study.

Parameter Definition Estimate units
Mean film thickness
Liquid viscosity
Liquid density
Hamaker constant
Surface tension of clean film
Mean surfactant concentration
Surface diffusion coefficient
Surface viscosity
Table 1: Estimates of relevant physical parameters

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 () axis, and is thus considered the most unstable mode. The dimensional momentum balance in the - and directions, and the continuity equation may be written as:

(2.0)
(2.0)
(2.0)

Here, and denote the velocity components in and directions, while is the fluid pressure. We have included a disjoining pressure term in our model, considering films to be thinner than , where van der Waals forces are prevalent (Edwards & Oron 1995; Berg 2010). In the present work, we set where 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):

(2.0)

Here is surface diffusivity, treated as a constant, is the surface gradient operator, is the surface velocity field, and n is the unit vector normal to the free surface. The boundary conditions at the free surface (= (, )) 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:

(2.0)
(2.0)

where, t is the unit vector tangent to the free surface, is the stress tensor and the dilatational and shear components of surface viscosities will hereafter be written in terms of a single parameter, 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 (). These are written as:

(2.0)
(2.0)

2.2 Scalings and non-dimensionalization

The scalings for the different variables are explained next, where symbols with the tilde () decoration denote dimensional versions. The spatial coordinates and are scaled as:

(2.0)

where is the characteristic length along (for e.g., the wavelength of a typical interfacial perturbation). We further assume that , 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:

(2.0)

The scalings in (2.2) 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 Marangoni stresses), and would be more apt for ultrathin films () (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:

(2.0)

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 :

(2.0)

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 (). The actual choice of will have implications on the model used for variable surface viscosity effects which will be discussed later. The surface tension is scaled with respect to its clean interface value (i.e. in the limit whereas the surface viscosity, , is scaled with , the choice of which depends on the type of variation used for (see Appendix B),

(2.0)

Several dimensionless characteristic numbers appear upon rewriting 2.1-2.1 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 (non-dimensional surface tension). Further, the Marangoni number M and the surface viscosity gradient parameter, , represent the variation of surface tension and surface viscosity with surfactant concentration. A linear variation with concentration is assumed for surface tension:

(2.0)

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

(2.0)
(2.0)

We use the linear model (2.2) for an interface with dilute surfactants, thus the reference concentration, has a value in the dilute limit. For the nonlinear model (2.2), the jammed state concentration value, of the system would intuitively be the suitable non-dimensionalizing factor for 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 is related to the surface viscosity gradient parameter by the relation,

(2.0)

where and 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 , in accordance with earlier studies on surface viscosity by Scriven (1960).

Parameter Definition Order of magnitude estimates
Aspect ratio
Boussinesq number
Disjoining pressure number
Marangoni number
Surface viscosity gradient number
Peclet number
Reynolds number
Table 2: Estimates of non-dimensional system parameters evaluated using characteristic velocity scale, and dimensional estimates from Table 1.

The disjoining pressure number, in Table 2, defined as has been deliberately scaled to contain , 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 quantity, a weak Marangoni limit is defined by rescaling it as , where the new O(1) parameter is . Then, Marangoni effects would drop out in the lubrication limit when terms are ignored. Following the definition of in (2.2), the weak Marangoni limit would imply, at leading order,

(2.0)

For the surface viscosity, 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, may be rescaled such that it is retained at leading order, and the boundary conditions are self-consistent. However, as (aspect ratio) can be made arbitrarily small, presenting the surface viscosity terms at 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 (2.1-2.1) reduce to the following when only the zeroth (leading) order terms are considered, denoted with the superscript

(2.0)
(2.0)
(2.0)

Using the symmetry conditions (2.1) at gives:

(2.0)
(2.0)

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

(2.0)

Integrating (2.3) and comparing it with (2.3), and substituting (2.3-2.3), we obtain:

(2.0)

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

(2.0)

This is also consistent with (2.3). The kinematic condition (2.1) and the surfactant transport equation (2.1) provide two nonlinear evolution equations for and , with 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 (2.3) and (2.3) (refer Appendix A.2, A.3). To these, we substitute (2.2) or (2.2) to account for variable viscosity. The final system of three nonlinear evolution equations read:

(2.0)
(2.0)
(2.0)

or,

(2.0)

Here, are as defined in Table 2. The superscript (0) has been omitted in (2.3 - 2.3) and the remainder of the paper for easy readability. Equations (2.3 - 2.3) reduce to those derived by DeWit et al. (1994) in the limit , and to those of Matar (2002) when we set , in the absence of surfactant solubility.

3 Linear Stability Analysis

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

(3.0)
(3.0)
(3.0)

Here is the growth rate and 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:

(3.0)
(3.0)

where, and are the non-dimensional base states for surfactant concentration in the linear model and NPM respectively, as has already been discussed in deriving (2.2). 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 and in (3) and (3) respectively.

3.1 Parametric studies

Figure 2: Dispersion curves for linear and nonlinear models. Left panels are for the linear model (equation 3) with and the right panels are for NPM (equation 3) with for varying , and . (a,b) Varying with fixed; (c,d) varying Bo with fixed; (e,f) varying M with fixed. Other parameters are fixed at . As is clearly evident, cut-off wavenumber, , is insensitive to changes in and and depends solely on . Figure 2 shows that growth rate is highly sensitive to surface viscosity effects in NPM and has been analysed in detail in section §3.2.
Figure 3: Dispersion curves for linear and nonlinear models. Left panels are for the linear model (equation 3) with and the right panels are for NPM (equation 3) with . (a,b) Varying Re with fixed; (c,d) varying and in the linear and NPM respectively with fixed. Other parameters are fixed at .

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, and a cut-off wavenumber, ,

(3.0)

above which the film is stable. Fig. 2 and 2 show the linear stability characteristics with variation in the parameter for the two models. We set for the dispersion curves of NPM. Intuitively, a decrease in 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 . 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: Re and Pe 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 Re 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 Re. As regards to Pe, a higher surface diffusivity of the surfactants (lower values of Pe) 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 ) which is very different from that of the linear model. This is discussed in §3.2 in detail. Both the models follow (3.0). Equations (3) and (3) become identical when the parameters are chosen to satisfy (2.2) in the limit of .

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 . 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 . This limit of (3) is distinctly different from that of the linear model in (3) which is valid in the limit of small . 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, . As , . Equation (3) can then be rewritten in terms of :

(3.0)

where the coefficients to are functions of and can be obtained by comparing (3.0) with (3):

(3.0)
(3.0)
(3.0)
(3.0)
(3.0)

In the limit , Equation (3.0) is a standard singular perturbation problem for growth rate . Fig. 4 shows plots of the sole positive root of (3.0) for different values of . It is apparent that the growth rate scales with and singular root appears in the region . We therefore define a rescaled growth parameter , so that equation (3.0) now becomes

(3.0)

A standard expansion in in the form results in two regular roots in the outer region, i.e. . Focussing on the root which is unstable, we obtain

(3.0)
(3.0)

To obtain the singular root of the inner region which has a boundary layer-like structure, we define a rescaled variable , where . The solution to the inner region in terms of this rescaled variable at leading order is:

(3.0)

where . Figure 4 shows a comparison of inner solution, (dot-dashed line), outer solution, (dashed line) and the numerical root, (solid line) obtained from (3.0). Examining (3.0), it can be shown that the necessary condition for stable growth rates () is . Hence we arrive at a necessary condition for obtaining stable modes in the space :

(3.0)

We can therefore define a critical concentration, . Any value of should result in inherently stable films. Also note that since is, by definition, bounded by: , we can hypothesize that no amount of surfactants can stabilize the system if the product of the system parameters, . Examining the system parameters in real physical systems, we find that 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.

Figure 4: (a) Unstable root for dispersion relation (3.0) with varying . Growth rate scales with and curve exhibits a sharp gradient near . (b) Comparison of perturbation solutions showing inner, , (dot-dashed) and outer, , (dashed) lines along with full dispersion relation from equation (3.0) compared with numerical roots (solid). Other parameters are fixed at

It is easy to derive a simple relationship between the amplitude of tangential velocity at the interface, and height perturbation, using the kinematic condition (2.3):

(3.0)

The above relationship suggests that the tangential velocity at the interface, , has a phase-difference with the interface height, since 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 in the jamming limit. This dependence can be written in the simple form

(3.0)

where we have used the definition of and functional dependence on other parameters, , is not reproduced for simplicity. Substituting into (3.0), we get

(3.0)

As , it is clear that , 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 (2.3), (2.3) 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 (2.3) - (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, (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:

(4.0)
(4.0)
(4.0)

where is the fastest growing wavelength, determined from linear stability analysis. The perturbation amplitude for surfactant concentration, is set using:

(4.0)

which is obtained from linear stability analysis. Here, we use and for the linear and nonlinear models respectively. Further, we set , with being the maximum value of the growth rate predicted by the linear theory, at . The time steps used in the initial stages are of the order and reduces to at the latter stages of our numerical simulations. The spatial domain is discretised into 256 nodal points () 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 and use adaptive time-stepping with lower end values of the order . We assume the film to have ruptured when the minimum height of the film, .

4.2 ‘Dilute limit’ analysis

4.2.1 Evolution and parametric studies

Figure 5: Spatiotemporal evolution of interface height (a), tangential velocity (b) and surfactant concentration (c) for the linear model obtained by solving equations (2.3)-(2.3) at various times . Initial condition at is shown with a dotted curve and rupture profile shown with dashed curve occurs at . The other parameters used are . (d) Comparison of numerically obtained growth rate and linear theory predictions for the same set of parameters. Here is the difference between the maximum and minimum values of at any time (see Dey et al. (2019)).

The nonlinear dynamics using the linear model (2.2) 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 , (b) concentration , (c) velocity at the interface, for a typical case: . The profiles have been plotted on spatial coordinates over a domain length . Given that the boundary conditions are periodic, successive crests/troughs in the profiles are separated by , 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 and indicate the possibility of self-similar solutions close to rupture.

Figure 6: Nonlinear results from the linear model obtained by solving (2.3)-(2.3). (a) Effect of varying Bo on height profiles at the time of rupture with fixed. corresponds to the case of zero surface viscosity and exhibits a cusp-like solution. Increasing causes the interface to flatten near rupture. (b) Evolution of minimum film thickness () with time and comparison of rupture time predictions from linear stability analysis (solid) and nonlinear simulations (circles) for varying Bo as used in (a). (c) Effect of varying on height profiles at the time of rupture with . As increases, the film profile near rupture approaches a cusp-like solution. (d) Evolution of minimum film thickness () with time and comparison of rupture time predictions from linear stability analysis (solid) and nonlinear simulations (squares) for varying . Other parameter values are fixed at .

Figure 6 shows the effect of surface viscosity (Bo) and its linear concentration dependence (represented by ) 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 produces a more cusp-like rupture, as shown in figure 6. This is understandable since surface viscosity from (2.2) at rupture location reduces to simply becoming since near rupture. Therefore increasing reduces the local value of surface viscosity near rupture location. It has to be noted that the maximum value of cannot exceed unity as a result of positivity of . The upper limit of causes the non-dimensional surface viscosity to become identical to the non-dimensional surface concentration, i.e. . Since at the rupture location, 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 in figure 6. This profile matches very well with the case of zero surface viscosity, 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 as seen in Fig. 6.

4.2.2 Self-similar solutions

Figure 7: Results from self-similar analysis. (a) Plot of (circles), (diamonds), and (squares) vs plotted in the log-log scale for reveals power-law dependence. Relevant scaling exponents are extracted from the slopes of these profiles (see text for details). In the limit of small , the slopes assume the values of -1.5 (dot-dashed line), -0.3 (dotted line) and 1 (dashed line). (b) Self-similar film profiles at the rupture location close to the break-up time for . (c) Variation of surface viscosity with time at rupture location for linear model with (solid) and constant surface viscosity (dashed), . The value of is also shown in the figure and it denotes the value of reference surface viscosity on a clean interface. At rupture location, surfactant concentration goes to zero causing to approach . Note that surface viscosity at rupture location for is lower than that for a case with constant surface viscosity . (d) Self-similar film profiles at the rupture location close to the break-up time for the case of constant surface viscosity, . Interface profiles in (b) exhibit a higher value of curvature than in (d) owing to lower surface viscosity as shown in (c). Other parameters are .

Next, we investigate the presence of self-similar behaviour in free film profiles in the vicinity of rupture location. Equations (2.3)-(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 , and in the form:

(4.0)

where is the similarity variable and , and being the spatial location and time of rupture, respectively. Substituting (4.2.2) into (2.3)-(2.3) and subsequent simplification yields:

(4.0)
(4.0)
(4.0)

To obtain scalings followed by , and , we follow the approach used by Matar (2002) and Vaynblat et al. (2001). These scalings allow us to extract the exponents without knowledge of rupture time a priori. The natural logarithm of the curvature , velocity and surfactant concentration on the surface at are plotted against the natural logarithm of film thickness at . Three algebraic equations in terms of can be generated by comparing the slopes of , and , 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 in all the terms of equations (4.0)-(4.2.2) to vanish. An algebraic equation can be written by setting the powers of 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 . Most of the solutions can be discarded on physical grounds. For example, since during rupture as evident in figure 5, we require and . Similarly, tangential velocity, 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 in (4.2.2). We can also argue on similar grounds that in the definition of similarity variable . Of the many solutions that are possible for , 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 dependence in all the terms since the equations do not admit perfect self-similarity as was noted in earlier works. The terms which retain dependence in (4.0)-(4.2.2) become vanishingly small in the limit of . 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 -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 . 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 . We consider four test cases for the study as governed by equations (2.3)-(2.3): (i) which corresponds to the case of zero surface viscosity as reported by Vaynblat et al. (2001); (ii) where surface viscosity and concentration assume equivalence; (iii) which corresponds to the case of variable surface viscosity and is a canonical case for ; and (iv) which corresponds to the case of constant surface viscosity as reported by Matar (2002).

Figure 8: Comparison of self-similarity scaling exponents with various surface viscosity models. For reference, surface viscosity in the linear model assumes the form and allowing for comparison with existing literature for various values of and .

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

Results for are shown in figures 7 and 7. Following the procedure as explained above, the following set of exponents are found: . These exponents values are different when compared to the constant viscosity case studied by Matar (2002). In our case, this occurs for and yields the exponents 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 , surface viscosity all along the interface reduces to the expression, . 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: 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 . 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 , since 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 . Our self-similarity analysis for the variable surface viscosity model results in scalings that lie between the regimes (zero surface viscosity) and (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 and (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 . The interface develops towards a more cusp-like profile (closer to profiles for ) 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 , surface viscosity effects are suppressed completely at the rupture point and a cusp-like profile of a clean interface, (or of a particle-laden interface with no surface viscosity, ) is recovered. This also means that however large the initial surface viscosity we choose, we can get a cusp-like profile if we have . 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 (4.0-4.2.2) for various cases. With (case (i)) and (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 (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 and for case (iv), we get . 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, , in equation (4.0) 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

Figure 9: Equation (2.2) plotted for three different values of exponent, . The higher the value of the exponent, faster the divergence of surface viscosity with concentration. In the dilute limit, the three models approach the linear model, albeit with different slopes.
Figure 10: Evolution of film profiles using NPM obtained by solving equations (2.3), (2.3) and (2.3). (a) Height profiles (solid) and surfactant profile (dot-dash) at the instance of rupture for . The dashed lines represent the initial perturbation. (b) Height profiles for three different values of surfactant concentration ranging from dilute (), intermediate () and jammed limit () occurring at three different rupture times, respectively. The interface profile is remarkably different in the three cases, so is the rupture time. (c) Corresponding Tangential velocities at rupture, , for varying , plotted just before rupture. (d) Maximum tangential velocity (dot-dash) and rupture time, (solid) plotted for various values of . Other parameters are fixed at .

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 (2.2). Figure 9 illustrates the variation of for NPM (2.2) for various values of . Higher values of causes to diverge at lower values of . To compare film profiles for both the models in dilute limit, the values of and are chosen so as to satisfy the relation (2.2). 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 and 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 , an intermediate limit at and a near-jammed state limit at