Instability and dynamics of thin viscoelastic liquid films (original) (raw)

Instability and dynamics of thin viscoelastic liquid films

G. Tomar 1{ }^{1}, V. Shankar 2,a{ }^{2, a}, S.K. Shukla 2{ }^{2}, A. Sharma 2,b{ }^{2, b}, and G. Biswas 1,c{ }^{1, \mathrm{c}}
1{ }^{1} Department of Mechanical Engineering, Indian Institute of Technology, Kanpur 208016, India
2{ }^{2} Department of Chemical Engineering, Indian Institute of Technology, Kanpur 208016, India

Received 18 March 2006 and Received in final form 22 April 2006 /
Published online: 20 June 2006 - © EDP Sciences / Società Italiana di Fisica / Springer-Verlag 2006

Abstract

The instability, rupture, and subsequent growth of holes in a thin Jeffreys-type viscoelastic film under the influence of long-range van der Waals force are investigated using both linear stability analysis and nonlinear numerical solutions. The linear stability analysis of full governing equations valid for arbitrary wave numbers shows that although fluid rheology does not influence the dominant length scale of the instability, it significantly affects the growth rate. It is shown that neglect of inertia and solvent dynamics results in a nonphysical singularity in the growth rate beyond a critical value of relaxation time. We further carry out numerical simulations of a set of long-wave, nonlinear differential equations (also derived in Rauscher et al., Eur. Phys. J. E 17, 373 (2005)) governing the evolution of the free surface. The nonlinear simulations, in their domain of validity, confirm the results of the linear analysis. Interestingly, results from nonlinear simulations further show that both for Newtonian and viscoelastic liquids, the shape and the dewetting dynamics of a hole are identical when examined in terms of a rescaled time which depends on rheological parameters. Thus, viscoelasticity of Jeffreys type merely accelerates the growth rate, without however affecting the important morphological characteristics.

PAC5. 68.15.+e Liquid thin films - 47.50.-d Non-Newtonian fluid flows - 47.15.gm Thin film flows 47.55.dr Interactions with surfaces

1 Introduction

Liquid films used in fundamental research for studying the physics due to confinement as well as in technological applications are getting progressively thinner ( <100 nm<100 \mathrm{~nm} ), where the role of intermolecular interactions becomes important. These highly confined films on low-energy substrates are inherently unstable due to the long-range interactions, such as van der Waals and electrostatic forces, and spontaneously break up by the amplification of the fluctuations in the free surface. The amplification is largest for the most unstable (i.e. fastest growing) wave, and thus the deformation of the free surface occurs on a particular length scale, which defines the lateral spacing of the resulting pattern. Similarly, the dynamics or the time scale of the pattern formation also is controlled by the most unstable mode. The instability, dynamics and morphology of rheologically simple, purely viscous (Newtonian) thin liquid films are now relatively well understood [1-16]. A recent review article by Thiele [17] summarizes the current status and discusses open questions in the subject of dewetting characteristics of thin Newtonian films. Most applications as well the fundamental experimental stud-

[1]ies of thin-film stability and dewetting, however, employ polymeric liquids such as melts of polystyrene, PMMA, etc. [18-20]. Also, technological applications ranging from optoelectronic coatings, polymer MEMS, sensors to food processing involve non-Newtonian, viscoelastic liquid films [21,22][21,22]. Yet, there is only a limited theoretical understanding of the instability and dynamics of thin polymeric films, especially during the early stages of the amplification of surface instability leading to its break-up. Indeed, experimental studies on thin polymer films are frequently interpreted within the framework of the stability and dynamics of purely viscous Newtonian films. There are some theoretical studies [23-25] on viscoelastic liquid films related to the long-time evolution of instability which commences after the break-up of the film leading to dewetting by the growth of holes. However, the early stages of the instability until the film break-up in viscoelastic films are only recently beginning to be addressed [26-28].

We present here a detailed account of the linear stability characteristics and results from nonlinear simulations of the initial surface instability as well as the later stages of dewetting and hole growth in a thin viscoelastic liquid film. In particular, we address the questions related to the length and time scales of the instability in a viscoelastic liquid film, and compare our results with the behavior of Newtonian films. It is well known that polymer melts as well as polymeric solutions

a{ }^{a} e-mail: vshankar@iitk.ac.in
b{ }^{\mathrm{b}} e-mail: ashutos@iitk.ac.in
c{ }^{\mathrm{c}} e-mail: gtm@iitk.ac.in ↩︎

display viscoelastic properties [29,30][29,30]. The simplest of the constitutive models that capture the essential aspects of viscoelastic liquids are the Maxwell and Jeffreys class of models, which adequately describe the elastic nature of polymer melts and solutions, respectively. A very recent study [25] of viscoelastic dewetting by the growth of holes has employed the Jeffreys model for studying the shape of a rim that surrounds a growing hole using linear stability analysis. Another recent study [28] also employed a Maxwell-type model to understand the linear instability of a confined polymer film subjected to electric fields. There have been a few other studies that have used viscoelastic constitutive relations in related areas of thinfilm flows. For example, in reference [31], an Oldroyd-B constitutive relation (which is nothing but the Jeffreys model with a convected time derivative; see below) was used to model the viscoelastic film in order to study the effect of viscoelasticity on surfactant transport and fluid flow. In reference [27] the liquid-gas interface was treated to be of viscoelastic nature, which was modeled using a Maxwell model, whereas the Newtonian constitutive relation was used to model the bulk of the dewetting film. The interfacial elasticity in such situations was shown to improve the thin-film stability. In references [32,33] the effect of a linear viscoelastic solid on the dewetting characteristics of a Newtonian liquid film above it was analyzed. It was concluded, through a linear stability analysis followed by nonlinear simulations, that the deformable viscoelastic solid substrate has a destabilizing effect on the liquid film and thus enhances the instability.

In an important earlier work, Safran and Klein [26] considered the linear stability of a generic linear viscoelastic model in the long-wave limit, and showed that the long-wave instability of the thin film is absent when the material has a zero-frequency shear modulus or elasticity. In the case of polymer melts and solutions considered in our study, there is no zero-frequency (permanent) elasticity because only liquid polymeric films (without permanent crosslinks) are considered. That the instability in a viscoelastic film is indeed of long-wave nature is not guaranteed in the first place because in solid viscoelastic films (with a finite zero-frequency elastic modulus) [34-37] the instability is in fact at short waves. Moreover, the effect of nonlinearities on the length and time scales in viscoelastic liquid films has not been previously investigated, and needs to be explored. We also address both of the above issues in this paper.

The remainder of this paper is structured as follows. The governing differential equations describing the dynamics of the thin viscoelastic liquid film have been formulated using the Jeffreys model for the constitutive relation in Section 2. A detailed account of the linear stability characteristics of the thin-film system is presented in Section 3, where we demonstrate some intriguing consequences of taking the special limit of a liquid polymer melt: we show that in the noninertial limit (Sect. 3.2), the growth rate shows a singularity when the relaxation time of the polymer increases beyond a certain critical value; we also show how this singularity can be removed by considering either

Fig. 1. A typical Type-I variation in ΔG\Delta G with thickness of the film [11]. In Type-I systems intermolecular forces are attractive for films of all thicknesses except for a very short-range repulsion at l0=0.137 nml_{0}=0.137 \mathrm{~nm}.
inertial effects in the fluid (Sect. 3.3) or by considering a polymer solution (polymer + solvent) in the limit of small solvent concentration (Sect. 3.4). Using our detailed discussion of the linear dispersion relation, we also establish the parameter regimes in which the instability remains “long wave” (i.e. the most unstable wavelength is large compared to the unperturbed film thickness). In this domain of validity, nonlinear evolution equations for the freesurface position are derived using a standard lubricationlike approximation [10], and the resulting nonlinear coupled equations describing the dynamics and morphological evolution of the thin film are solved numerically in Section 4. A comparison is made between the results obtained for viscoelastic liquids with those for Newtonian liquids in order to discern the effect of viscoelasticity on the dynamics and morphological evolution of the film.

2 Formulation

The equations governing fluid motion in the thin film are, respectively, the mass and momentum conservation equations. The fluid is incompressible, and so the mass conservation equation is simply given by ∇⋅v=0\nabla \cdot \mathbf{v}=0. The momentum conservation equation for the fluid in the thin film, after neglecting gravity and evaporation/condensation effects, can be written as

ρDvDt=−∇Π+∇⋅τ\rho \frac{D \mathbf{v}}{D t}=-\nabla \Pi+\nabla \cdot \boldsymbol{\tau}

where ρ\rho is the density of the fluid, Π\Pi is the total pressure (p+ϕ),p(p+\phi), p is the hydrodynamic pressure in the film and ϕ\phi is the excess pressure in the film due to long-range van der Waals interaction (negative of the disjoining pressure), ΠDt\frac{\Pi}{D t} is the substantial derivative, ∇(∂x,∂y,∂z)\nabla\left(\partial_{x}, \partial_{y}, \partial_{z}\right) is the gradient operator, and τ\boldsymbol{\tau} is the stress tensor in the viscoelastic liquid, whose constitutive equation is discussed a little later. The excess pressure ϕ\phi is related to the excess free energy (ΔG)(\Delta G) per unit area of the film as ϕ(h)=∂ΔG∂h\phi(h)=\frac{\partial \Delta G}{\partial h}, where ΔG=A12πh2+Bh2,A,B>0\Delta G=\frac{A}{12 \pi h^{2}}+\frac{B}{h^{2}}, A, B>0. Figure 1 shows the variation of the interaction energy with the thickness of

the thin film. The behavior obtained using the above expression for the potential is for a purely attractive force field, and is referred to as “Type-I” behavior [11], wherein the film will only partially wet the surface of the substrate. The repulsive term in the energy expression (B/hq)\left(B / h^{q}\right) is due to short-range Born repulsion [6,7][6,7]. The minimum of the free energy occurs at a distance l0l_{0}, where

ϕ=∂ΔG∂h=0 and ΔG(l0)−ΔG(∞)=SLW<0\phi=\frac{\partial \Delta G}{\partial h}=0 \quad \text { and } \quad \Delta G\left(l_{0}\right)-\Delta G(\infty)=S^{\mathrm{LW}}<0

SLWS^{\mathrm{LW}} is the Lifshitz-van der Waals (LW) component of the spreading coefficient of the film material on a surface [10,15][10,15].

In this study, the viscoelastic liquid is modeled using the Jeffreys linear model [29]. The Jeffreys model is a generalization of the simple Maxwell model, with the inclusion of (Newtonian) solvent contribution to the stress tensor. The simple Maxwell model describes a viscoelastic fluid with a single relaxation time: for deformations occurring at time scales shorter than the relaxation time, the viscoelastic fluid behaves like a purely elastic solid, while for time scales longer than the relaxation time, the viscoelastic liquid behaves like a simple Newtonian fluid. In reality, however, a polymeric liquid (which the simple constitutive models purport to describe) has a spectrum of relaxation times [30], unlike the single relaxation time assumed in the simple Maxwell (and Jeffreys) models. Thus, the relaxation time occurring in the constitutive relation must be interpreted as the longest relaxation time exhibited by a real polymeric liquid. There is one further simplification that is inherent in simple models like Maxwell and Jeffreys ones, namely the use of a simple time derivative of the stress tensor instead of the more rigorous (albeit much more complicated because of nonlinearities) “upper-convected” time derivative. These “upper-convected” terms, by which the stress tensor gets advected and rotated by vorticity, can become significant when the shear rate in the flow becomes large [29]. However, in this study, we are interested in studying the spontaneous dewetting of an initially quiescent liquid film due to long-range intermolecular forces, and we neglect the upper-convected terms in the time derivative. Importantly, however, the linear stability analysis is about a quiescent base state, and therefore the upper-convected terms do not make any contributions to the linearized stability equations. This is because the upper-convected terms are quadratic in the stresses and velocity gradients, and when expanded about a quiescent base state, do not make any linear contributions to the analysis. Thus, the results obtained from the present linear stability analysis are exact even with the inclusion of the upper-convected terms. The results obtained from the nonlinear simulations without the upper-convected terms, however, are expected to be qualitatively accurate. A similar approach has been taken in the recent study by Rauscher et al. [25], where the authors studied the decay of a capillary ridge in a viscoelastic liquid. When the ordinary time derivatives are replaced by the upper-convected derivatives, the Maxwell and Jeffreys models are referred to as the “upper-convected Maxwell” model and “Oldroyd-B” model, respectively.

Fig. 2. Schematic of a thin film (<100 nm)(<100 \mathrm{~nm}) in which instability is engendered due to intermolecular forces.

The constitutive relation under the Jeffreys model [29] is given by

τ+λ1∂tτ=η((∇v+∇vT)+λ2∂t(∇v+∇vT))\boldsymbol{\tau}+\lambda_{1} \partial_{t} \boldsymbol{\tau}=\eta\left(\left(\nabla \mathbf{v}+\nabla \mathbf{v}^{T}\right)+\lambda_{2} \partial_{t}\left(\nabla \mathbf{v}+\nabla \mathbf{v}^{T}\right)\right)

where τ\boldsymbol{\tau} is the dimensional stress tensor. λ1\lambda_{1} and λ2\lambda_{2} are the relaxation and retardation time constants, respectively, η(=ηp+ηs)\eta\left(=\eta_{p}+\eta_{s}\right) is the sum of viscosities due to polymer and the solvent contributions, and λ2\lambda_{2} is related to λ1,ηp\lambda_{1}, \eta_{p} and ηs\eta_{s} by the relation

λ2=λ1ηs(ηp+ηs)\lambda_{2}=\lambda_{1} \frac{\eta_{s}}{\left(\eta_{p}+\eta_{s}\right)}

Defining ηr=ηs/η\eta_{r}=\eta_{s} / \eta, we observe that if ηr=0\eta_{r}=0 (or, equivalently, λ2=0\lambda_{2}=0 ), we obtain the case of a pure polymer melt and the Jeffreys model is then equivalent to the simple Maxwell model used for describing the viscoelastic behavior of polymer melts. Another limiting behavior of the model can be obtained by setting ηr=1\eta_{r}=1, which makes the model equivalent to that of a purely Newtonian fluid. Thus, the Jeffreys model can describe the dynamical behavior ranging from viscous Newtonian liquids to polymer melts.

We consider a thin ( <100 nm<100 \mathrm{~nm} ) initially flat film (shown schematically in Fig. 2) of a polymeric liquid (e.g. polydimethyl siloxane, polystyrene) resting on a solid substrate (such that SLW<0S^{\mathrm{LW}}<0 ). In what follows, we consider the planar (i.e. two-dimensional) flow approximation and neglect variations in the zz-direction. Boundary conditions for equation (1) are as follows. The kinematic condition at the free surface of the film in the two-dimensional case can be written as

∂h∂t+vxs∂h∂x=vys\frac{\partial h}{\partial t}+v_{x_{s}} \frac{\partial h}{\partial x}=v_{y_{s}}

where vxsv_{x_{s}} and vysv_{y_{s}} are the xx and yy components of the velocity at the free surface of the film. The no-slip boundary condition is assumed at the liquid-solid interface (y=0)(y=0)

vx=0,vy=0v_{x}=0, \quad v_{y}=0

At the free surface, the zero shear stress condition applies:

τxy=0 at y=h(x)\tau_{x y}=0 \quad \text { at } \quad y=h(x)

Normal stress balance at the free surface y=h(x)y=h(x) yields

−Π+τyy=−ϕ(h)−γκ-\Pi+\tau_{y y}=-\phi(h)-\gamma \kappa

where κ\kappa is the curvature at the interface, and γ\gamma is the surface tension coefficient at the free surface of the film.

3 Linear stability analysis

We first carry out a linear stability analysis of the thinfilm system subjected to disturbances with arbitrary wavelengths. In order to perform a linear stability analysis (LSA) of the system we consider a base state where the liquid film is quiescent and the base state parameters are given by Π0=A/6πh03,h=h0,vx=0,vy=0\Pi_{0}=A / 6 \pi h_{0}^{3}, h=h_{0}, v_{x}=0, v_{y}=0. Small fluctuations in the form of Fourier modes are imposed onto the above base state:

Π=Π0+Π~(y)eikxesth=h0+h~eikxestvx=v~x(y)eikxestvy=v~y(y)eikxest\begin{aligned} \Pi & =\Pi_{0}+\tilde{\Pi}(y) e^{i k x} e^{s t} \\ h & =h_{0}+\tilde{h} e^{i k x} e^{s t} \\ v_{x} & =\tilde{v}_{x}(y) e^{i k x} e^{s t} \\ v_{y} & =\tilde{v}_{y}(y) e^{i k x} e^{s t} \end{aligned}

where kk is an arbitrary wave number of perturbation and ss is the associated growth rate coefficient. In the linear stability analysis, since the fluctuations are assumed to be small, all the free-surface boundary conditions are linearized about the unperturbed free surface at y=h0y=h_{0}. Linearizing the excess pressure ϕ\phi about the base height h0h_{0} we get

ϕ(h)=A6πh03−A(h−h0)2πh04\phi(h)=\frac{A}{6 \pi h_{0}^{3}}-\frac{A\left(h-h_{0}\right)}{2 \pi h_{0}^{4}}

Linearized constitutive equations using the Jeffreys model are obtained from equation (3) as

τ~xx=2η[λ1ηrs+1][λ1s+1](ik)v~xτ~xy=η[λ1ηrs+1][λ1s+1][v~x′+(ik)v~y]τ~yy=η[λ1ηrs+1][λ1s+1](2v~y′)\begin{aligned} & \tilde{\tau}_{x x}=2 \eta \frac{\left[\lambda_{1} \eta_{r} s+1\right]}{\left[\lambda_{1} s+1\right]}(i k) \tilde{v}_{x} \\ & \tilde{\tau}_{x y}=\eta \frac{\left[\lambda_{1} \eta_{r} s+1\right]}{\left[\lambda_{1} s+1\right]}\left[\tilde{v}_{x}^{\prime}+(i k) \tilde{v}_{y}\right] \\ & \tilde{\tau}_{y y}=\eta \frac{\left[\lambda_{1} \eta_{r} s+1\right]}{\left[\lambda_{1} s+1\right]}\left(2 \tilde{v}_{y}^{\prime}\right) \end{aligned}

After substituting for τ~\tilde{\boldsymbol{\tau}} in the momentum equation (1) (where primes denote differentiation with respect to yy ), we obtain the following system of equations that describe the dynamics of small fluctuations in the liquid film: xx-momentum equation:

−Π~k2+η[λ1ηrs+1][λ1s+1][v~y′′′−k2v~y′]=ρsv~y′-\tilde{\Pi} k^{2}+\eta \frac{\left[\lambda_{1} \eta_{r} s+1\right]}{\left[\lambda_{1} s+1\right]}\left[\tilde{v}_{y}^{\prime \prime \prime}-k^{2} \tilde{v}_{y}^{\prime}\right]=\rho s \tilde{v}_{y}^{\prime}

yy-momentum equation:

−Π~′+η[λ1ηrs+1][λ1s+1][v~y′′−k2v~y]=ρsv~y-\tilde{\Pi}^{\prime}+\eta \frac{\left[\lambda_{1} \eta_{r} s+1\right]}{\left[\lambda_{1} s+1\right]}\left[\tilde{v}_{y}^{\prime \prime}-k^{2} \tilde{v}_{y}\right]=\rho s \tilde{v}_{y}

Linearized kinematic condition at y=h0y=h_{0} :

h~s=v~y\tilde{h} s=\tilde{v}_{y}

Linearized shear free boundary condition at y=h0y=h_{0} :

v~x′+ikv~y=0\tilde{v}_{x}^{\prime}+i k \tilde{v}_{y}=0

No-slip boundary condition at the bottom surface y=0y=0 :

v~x=0,v~y=0\tilde{v}_{x}=0, \quad \tilde{v}_{y}=0

Linearized normal stress balance at y=h0y=h_{0} :

−Π~+τ~yy=Ah~2πh04−γk2h~-\tilde{\Pi}+\tilde{\tau}_{y y}=\frac{A \tilde{h}}{2 \pi h_{0}^{4}}-\gamma k^{2} \tilde{h}

Using the continuity equation, v~x=0\tilde{v}_{x}=0 can be written as v~y=0\tilde{v}_{y}=0. Differentiating the xx-momentum equation (Eq. (17)) with respect to yy and subtracting k2k^{2} times the yy-momentum equation (Eq. (18)) we get

v~yIV−(k2+q2)v~yII+k2q2v~y=0\tilde{v}_{y}^{I V}-\left(k^{2}+q^{2}\right) \tilde{v}_{y}^{I I}+k^{2} q^{2} \tilde{v}_{y}=0

where q=k1+ρs(1+sλ1)/((1+sλ1ηr)ηk2)q=k \sqrt{1+\rho s\left(1+s \lambda_{1}\right) /\left(\left(1+s \lambda_{1} \eta_{r}\right) \eta k^{2}\right)}.
A nondimensional measure of the inertial stresses in the film can be obtained by balancing the dimensional inertial and viscous stresses in equation (1), which yields a nondimensional group U=ρA/(6πh0η2)U=\rho A /\left(6 \pi h_{0} \eta^{2}\right) representing the strength of the inertial terms. For typical values of A=10−20 J,ρ=103 kg/m3,h0=10−8 mA=10^{-20} \mathrm{~J}, \rho=10^{3} \mathrm{~kg} / \mathrm{m}^{3}, h_{0}=10^{-8} \mathrm{~m} and η=0.1 Pas\eta=0.1 \mathrm{~Pas} for thin polymeric films, the value of the parameter UU is 5.3×10−95.3 \times 10^{-9}. Because of the smallness of this parameter, we are, at first sight, tempted to conclude that inertial stresses are indeed very small compared to viscous stresses in the film, and hence we first explore the consequences of completely neglecting the inertial terms in our analysis.

3.1 Zero-inertia limit

Neglecting the inertial terms, the above biharmonic equation reduces to

v~yIV−2k2v~yII+k4v~y=0\tilde{v}_{y}^{I V}-2 k^{2} \tilde{v}_{y}^{I I}+k^{4} \tilde{v}_{y}=0

The general solution for the above biharmonic differential equation (Eq. (24)) is given by

v~y=C1e−ky+C2ye−ky+C3eky+C4yeky\tilde{v}_{y}=C_{1} e^{-k y}+C_{2} y e^{-k y}+C_{3} e^{k y}+C_{4} y e^{k y}

where {C1,C2,C3,C4}\left\{C_{1}, C_{2}, C_{3}, C_{4}\right\} are constants to be fixed by boundary and interface conditions. Solving for constants using equations (19-21) we obtain

C1=[h~scosh⁡(h0k)][(sinh⁡(2h0k)−2h0k)]C2=−[h~ks(cosh⁡(h0k)+sinh⁡(h0k))][sinh⁡(2h0k)−2h0k]C3=[h~scosh⁡(h0k)][sinh⁡(2h0k)−2h0k]C4=−[h~ks(cosh⁡(h0k)−sinh⁡(h0k))][sinh⁡(2h0k)−2h0k]\begin{aligned} & C_{1}=\frac{[\tilde{h} s \cosh \left(h_{0} k\right)]}{\left[\left(\sinh \left(2 h_{0} k\right)-2 h_{0} k\right)\right]} \\ & C_{2}=-\frac{[\tilde{h} k s\left(\cosh \left(h_{0} k\right)+\sinh \left(h_{0} k\right)\right)]}{\left[\sinh \left(2 h_{0} k\right)-2 h_{0} k\right]} \\ & C_{3}=\frac{[\tilde{h} s \cosh \left(h_{0} k\right)]}{\left[\sinh \left(2 h_{0} k\right)-2 h_{0} k\right]} \\ & C_{4}=-\frac{[\tilde{h} k s\left(\cosh \left(h_{0} k\right)-\sinh \left(h_{0} k\right)\right)]}{\left[\sinh \left(2 h_{0} k\right)-2 h_{0} k\right]} \end{aligned}

The pressure can be obtained by substituting the above solution for v~y\tilde{v}_{y} in the following equation (from linearized xx-momentum equation (17)):

Π~=η(λ1ηrs+1)(λ1s+1)[v~y′′′−k2v~y′]k2\tilde{\Pi}=\frac{\eta\left(\lambda_{1} \eta_{r} s+1\right)}{\left(\lambda_{1} s+1\right)} \frac{\left[\tilde{v}_{y}^{\prime \prime \prime}-k^{2} \tilde{v}_{y}^{\prime}\right]}{k^{2}}

Fig. 3. Variation in the most dominant wave number (kmh0)\left(k_{m} h_{0}\right) with γˉ\bar{\gamma}.

Substituting the expression obtained for pressure above into the normal stress balance equation (Eq. (22)) we obtain the complete dispersion relation for the growth rate ss :

(8kηηrλ1cosh⁡2(h0k))s2+s[8kηcosh⁡2(h0k)−4k3h0×γλ1+A2πh03λ1(4h0k−2sinh⁡(2h0k))+2k2γλ1×sinh⁡(2h0k)]+A2πh03(4h0k−2sinh⁡(h0k))+2k2γ(sinh⁡(2h0k)−2h0k)=0\begin{aligned} & \left(8 k \eta \eta_{r} \lambda_{1} \cosh ^{2}\left(h_{0} k\right)\right) s^{2}+s\left[8 k \eta \cosh ^{2}\left(h_{0} k\right)-4 k^{3} h_{0}\right. \\ & \times \gamma \lambda_{1}+\frac{A}{2 \pi h_{0}^{3}} \lambda_{1}\left(4 h_{0} k-2 \sinh \left(2 h_{0} k\right)\right)+2 k^{2} \gamma \lambda_{1} \\ & \left.\times \sinh \left(2 h_{0} k\right)\right]+\frac{A}{2 \pi h_{0}^{3}}\left(4 h_{0} k-2 \sinh \left(h_{0} k\right)\right) \\ & +2 k^{2} \gamma\left(\sinh \left(2 h_{0} k\right)-2 h_{0} k\right)=0 \end{aligned}

The above-obtained full dispersion relation can be written in terms of nondimensional groups defined as

γˉ=(2πh02γ/A)W=1γˉλ1A6πh03ηsˉ=sγˉ6πh03ηλ1A\begin{aligned} \bar{\gamma} & =\left(2 \pi h_{0}^{2} \gamma / A\right) \\ W & =\frac{1}{\bar{\gamma}} \frac{\lambda_{1} A}{6 \pi h_{0}^{3} \eta} \\ \bar{s} & =s \bar{\gamma} \frac{6 \pi h_{0}^{3} \eta}{\lambda_{1} A} \end{aligned}

where γˉ\bar{\gamma} is a nondimensional energetic parameter representing the ratio of surface tension and intermolecular forces, WW is the “Weissenberg number” [29] which is a nondimensional relaxation time of the viscoelastic liquid, and sˉ\bar{s} is the nondimensional growth rate.

The full dispersion relation can now be written as

8kh0γˉ1/2cosh⁡2(kh0)Wηrsˉ2+[8kh0γˉ1/2cosh⁡2(kh0)−3Wγˉ(1−(kh0)2γˉ)(2γˉ1/2sinh⁡(2kh0)−4kh0γˉ1/2)]sˉ−3γˉ(1−(kh0)2γˉ)(2γˉ1/2sinh⁡(2kh0)−4kh0γˉ1/2)=0\begin{aligned} & 8 k h_{0} \bar{\gamma}^{1 / 2} \cosh ^{2}\left(k h_{0}\right) W \eta_{r} \bar{s}^{2}+\left[8 k h_{0} \bar{\gamma}^{1 / 2} \cosh ^{2}\left(k h_{0}\right)\right. \\ & \left.-3 W \bar{\gamma}\left(1-\left(k h_{0}\right)^{2} \bar{\gamma}\right)\left(2 \bar{\gamma}^{1 / 2} \sinh \left(2 k h_{0}\right)-4 k h_{0} \bar{\gamma}^{1 / 2}\right)\right] \bar{s} \\ & -3 \bar{\gamma}\left(1-\left(k h_{0}\right)^{2} \bar{\gamma}\right)\left(2 \bar{\gamma}^{1 / 2} \sinh \left(2 k h_{0}\right)-4 k h_{0} \bar{\gamma}^{1 / 2}\right)=0 \end{aligned}

Figure 3 shows the variation in kmh0k_{m} h_{0} ( kmk_{m} is the dimensional, most dominant wave number corresponding to the maximum growth rate) with γˉ\bar{\gamma}, where kmk_{m} is obtained by solving equation (31) for maximum sˉ\bar{s}, by setting dsˉ/dk=0\mathrm{d} \bar{s} / \mathrm{d} k=0. Because the parameters ηr\eta_{r} and WW appear
only in the coefficients of sˉ\bar{s} and sˉ2\bar{s}^{2} (in Eq. (31)), the critical and the most dominant wave numbers are therefore independent of ηr\eta_{r} and WW. Thus, we can conclude that the most dominant length scale of the problem is governed only by the thermodynamics of the system (γˉ)(\bar{\gamma}) and is independent of the rheological parameters η,ηr\eta, \eta_{r} and WW. For typical values of the Hamakar constant A=10−20 JA=10^{-20} \mathrm{~J}, γ=10−2 J/m2\gamma=10^{-2} \mathrm{~J} / \mathrm{m}^{2} and h0=10−8 mh_{0}=10^{-8} \mathrm{~m}, the value of γˉ\bar{\gamma} is estimated to be ∼(2π102)∼103\sim\left(2 \pi 10^{2}\right) \sim 10^{3}. In the regime γˉ>103\bar{\gamma}>10^{3}, kmh0k_{m} h_{0} is small (<10−2)\left(<10^{-2}\right) and decreases as 0.707/(γˉ1/2)0.707 /\left(\bar{\gamma}^{1 / 2}\right) with an increase in γˉ\bar{\gamma}. Therefore, we can define the length scale in the xx-direction at which patterns will form as L=h0γˉ1/2=[2πγh03A]1/2L=h_{0} \bar{\gamma}^{1 / 2}=\left[\frac{2 \pi \gamma h_{0}^{3}}{A}\right]^{1 / 2}, which is clearly large compared to h0h_{0} for typical values of γˉ\bar{\gamma} for polymeric thin films.

It is appropriate, therefore, to define a small parameter ϵ=h0/L\epsilon=h_{0} / L (i.e. γˉ=ϵ−2\bar{\gamma}=\epsilon^{-2} ). Corresponding to the long length scale LL, there is also a slow time scale proportional to ϵ−4\epsilon^{-4}, and consequently the dimensional growth rate ss scales as ϵ4\epsilon^{4}. The scaled growth rate, sˉ\bar{s}, however, remains O(1)O(1) as ϵ≪1\epsilon \ll 1. We choose a renormalized nondimensional wave number K=kLK=k L. Now, writing equation (27) in terms of ϵ\epsilon and KK and expanding the expression in series for ϵ→0\epsilon \rightarrow 0, we obtain, to leading order

sˉ2Wηr+sˉ(1−W(1−K2)K2)−(1−K2)K2=0\bar{s}^{2} W \eta_{r}+\bar{s}\left(1-W\left(1-K^{2}\right) K^{2}\right)-\left(1-K^{2}\right) K^{2}=0

Henceforth, to avoid a clumsy notation, sˉ\bar{s} would be written simply as ss. From the dispersion relation (Eq. (32)), the two solutions to the growth rate are obtained as

s1=−(1−W(K2−K4))2ηrW+(1−W(K2−K4))2+4ηrW(K2−K4)2ηrWs2=−(1−W(K2−K4))2ηrW−(1−W(K2−K4))2+4ηrW(K2−K4)2ηrW\begin{aligned} s_{1} & =\frac{-\left(1-W\left(K^{2}-K^{4}\right)\right)}{2 \eta_{r} W} \\ & +\frac{\sqrt{\left(1-W\left(K^{2}-K^{4}\right)\right)^{2}+4 \eta_{r} W\left(K^{2}-K^{4}\right)}}{2 \eta_{r} W} \\ s_{2} & =\frac{-\left(1-W\left(K^{2}-K^{4}\right)\right)}{2 \eta_{r} W} \\ & -\frac{\sqrt{\left(1-W\left(K^{2}-K^{4}\right)\right)^{2}+4 \eta_{r} W\left(K^{2}-K^{4}\right)}}{2 \eta_{r} W} \end{aligned}

We observe that no traveling solutions exist (i.e. ss is always purely real) for any values of ηr\eta_{r} and WW, unlike the result reported in [28] for the instability of a viscoelastic liquid film under the influence of an applied electric field. One of the two solutions to ss is always negative (i.e. decaying perturbations) irrespective of the values taken by ηr\eta_{r} and WW. The other solution becomes unstable in the range of wave numbers 0<K<10<K<1, and there are no wave numbers beyond K=1K=1 for which the system is unstable. For shorter wavelengths (K≫1)(K \gg 1), the growth rate asymptotically converges to −1/W-1 / W irrespective of the values taken by ηr\eta_{r}, except for ηr=1\eta_{r}=1. Maximizing s1s_{1} obtained from equation (33), we obtain an expression for the maximum growth rate as

smax⁡=−4+W+16+8(2ηr−1)W+W28ηrWs_{\max }=\frac{-4+W+\sqrt{16+8\left(2 \eta_{r}-1\right) W+W^{2}}}{8 \eta_{r} W}

Fig. 4. Variation in ratio of the most dominant wave numbers obtained using the full dispersion relation (Kmfull )\left(K_{m}^{\text {full }}\right) and the one obtained using the long-wave dispersion relation (Kmlong )\left(K_{m}^{\text {long }}\right) with γˉ\bar{\gamma}.

Figure 4 shows the variation in the ratio Kmfull /Kmlong K_{m}^{\text {full }} / K_{m}^{\text {long }} with γˉ\bar{\gamma}, where Kmfull K_{m}^{\text {full }} is the most dominant wave number obtained using the full dispersion relation (Eq. (31)) and Kmlong K_{m}^{\text {long }} is the most dominant wave number obtained using equation (32). We observe that beyond a critical γˉ\bar{\gamma} (∼102)\left(\sim 10^{2}\right), the ratio Kmfull /Kmlong →1K_{m}^{\text {full }} / K_{m}^{\text {long }} \rightarrow 1. Therefore, for typical values of γˉ∼103\bar{\gamma} \sim 10^{3} (estimated a little earlier) the long-wave assumption is indeed justified. The maximum growth rate always corresponds to the nondimensional wave number K=1/2K=1 / \sqrt{2}. Note, importantly, that this wave number does not depend on ηr\eta_{r} and WW. Therefore, viscoelasticity is predicted to have no effect on the wavelength of the fastest growing mode in the linear theory. The signature of this length scale is also reflected in the nonlinear evolution, where the domain size of the evolving patterns is intimately related to the wavelength of the fastest growing mode. The time of rupture can be estimated from the linear stability analysis as the time it would take for the surface fluctuation to reach the bottom solid surface, if it were to grow at a rate predicted by the linear theory:

Tr=1smax⁡log⁡e(1ζ)T_{r}=\frac{1}{s_{\max }} \log _{e}\left(\frac{1}{\zeta}\right)

3.2 Special case of a polymer melt (ηr=0)\left(\eta_{r}=0\right)

We now examine the expression for the growth rate for the special case of a polymeric melt (ηr=0)\left(\eta_{r}=0\right). The dispersion relation (Eq. (32)) reduces for this case to

s=(K2−K4)1−W(K2−K4)s=\frac{\left(K^{2}-K^{4}\right)}{1-W\left(K^{2}-K^{4}\right)}

Figure 5 shows the variation in the growth rate with wave number KK for different values of WW. For W=0.1W=0.1, the

Fig. 5. Variation of the growth rate with wave number KK for different values of WW, for ηr=0\eta_{r}=0.

Fig. 6. Variation of the growth rate (from the full dispersion relation) with wave number KK for different values of WW, for ηr=0\eta_{r}=0 and γˉ=103\bar{\gamma}=10^{3}.
maximum growth rate is small and for larger values of KK it converges to s=−10s=-10 (not shown in the figure). For W=4W=4, the growth rate diverges at K=1/2K=1 / \sqrt{2}. Beyond W=4W=4, the growth rate becomes unbounded for two wave numbers between which there is a zone of negative growth rate (shown for W=5W=5 and W=10W=10 in Fig. 5). This singular behavior of the growth rate beyond W=4W=4 is also observed in growth rates obtained from the full dispersion relation (Eq. (31)) for ηr=0\eta_{r}=0 and γˉ=103\bar{\gamma}=10^{3} (Fig. 6).

The fastest growing normal mode for a pure polymeric melt is K=1/2K=1 / \sqrt{2} for W≤4W \leq 4 which is again exactly the same as that obtained for Newtonian liquids. For W<4W<4, the maximum growth rate is given by s=1/(4−W)s=1 /(4-W). For W≥4W \geq 4, the growth rate becomes unbounded at two wave numbers, one on each side of K=1/2K=1 / \sqrt{2}. The

two wave numbers are given by K1=1+1−4/W2K_{1}=\sqrt{\frac{1+\sqrt{1-4 / W}}{2}} and K2=1−1−4/W2K_{2}=\sqrt{\frac{1-\sqrt{1-4 / W}}{2}}. Between these two wave numbers, equation (37) gives a negative growth rate. The maximum ss (which is negative) between K1K_{1} and K2K_{2} approaches zero for larger values of WW (Fig. 3). For W>4W>4, the range of wave numbers between K1K_{1} and K2K_{2} becomes a stable zone, the expression s=1/(4−W)s=1 /(4-W) gives the maximum of the negative growth rates (at K=1/2K=1 / \sqrt{2} ) in this stable zone. The gap between these two wave numbers (a zone of stable wave numbers) increases with an increase in WW and K1K_{1}, K2K_{2} approach 0 and 1 , respectively, in the limit of large WW. A very similar prediction has been reported recently by Wu and Chou [28] for the instability of a viscoelastic film under the influence of applied electric fields.

This nonphysical singularity for W≥4W \geq 4, discussed in more detail below, points to the incomplete description of the dynamics in this noninertial regime. Indeed, we show below that the inclusion of fluid inertia (in the absence of solvent) or the inclusion of solvent dynamics (in the absence of inertia) both remove the singularity in the growth rate. In other words, the limit of small inertia is not exactly the same as the limit of zero inertia in a polymer melt. However, as shown below, the noninertial description remains accurate for W<4W<4 in the absence of solvent dynamics.

3.3 Effect of fluid inertia for ηr=0\eta_{\mathrm{r}}=0

We now examine whether the nonphysical unbounded growth rate for W≥4W \geq 4 in case of polymer melts obtained in the previous subsection is a consequence of the neglected inertial terms in our foregoing analysis. The inertial terms were neglected because the nondimensional prefactor UU (multiplying the inertial term in the momentum equation) was estimated to be very small (∼10−3)\left(\sim 10^{-3}\right), and we implicitly assumed that the growth rate will remain finite, as a result of which the inertial term should remain much smaller than the viscous terms. However, when the growth rate obtained from such an analysis diverges beyond a critical value of WW, the product of a very small UU and a very large growth rate becomes finite, and this is a signal that inertial terms must be included in our analysis. Consequently, we now carry out a linear stability analysis for the case of a polymer melt (ηr=0)\left(\eta_{\mathrm{r}}=0\right), but with the inclusion of inertial terms in the fluid momentum equations. Substituting the solution of the biharmonic equation (Eq. (23)), with the boundary conditions (Eqs. (18-21)), in the normal balance condition (Eq. (22)) we obtain the full dispersion relation as F(s,k,W,U,γˉ)=0F(s, k, W, U, \bar{\gamma})=0 (see App. A). Figure 7 shows the variation in the growth rate with the nondimensionalized wave number KK for W=3.8,4.0W=3.8,4.0 and 6.0 . We observe that the singularity that was present in the ss vs. kk curves in the noninertial case for W>4W>4 is completely absent, and the ss vs. kk curves exhibit the usual behavior with the inclusion of fluid inertia at all wave numbers. This implies that the inertial effects in the present case are an example of “singular perturbation”, where a small parameter ( UU in the present discussion) multiplies a physical variable (growth

Fig. 7. Variation of the growth rate with wave number KK with the inclusion of fluid inertia: W=3.8,W=4.0W=3.8, W=4.0 and W=6.0W=6.0. U=10−6,γˉ=103U=10^{-6}, \bar{\gamma}=10^{3} and ηr=0\eta_{\mathrm{r}}=0 (pure polymer melt).

Fig. 8. Effect of fluid inertia on the variation of the maximum growth rate with W:U=0W: U=0 and U=10−6U=10^{-6} in the case of a polymer melt (ηr=0,γˉ=103)\left(\eta_{\mathrm{r}}=0, \bar{\gamma}=10^{3}\right).
rate in this case) which becomes large such that the product of the two quantities cannot be ignored. The maximum value of ss is still at K=1/2K=1 / \sqrt{2}, upon the inclusion of inertia. While the growth rate remains finite with increase in WW, it nonetheless increases very rapidly when WW is increased beyond 4. For instance, the maximum growth rate obtained for W=6.0W=6.0 is four orders of magnitude larger than for the case W=3.8W=3.8. The maximum growth rate for W=4.0W=4.0 is ∼250\sim 250. The dominant wave number is Km=1/2K_{m}=1 / \sqrt{2} for all values of WW with U≠0U \neq 0 in contrast to the case with U=0U=0 where, for W>4W>4 there exists two wave numbers K1K_{1} and K2K_{2} (around K=0.707K=0.707 ) for which the growth rate diverges. Between these two wave numbers, we showed earlier that there is a negative growth rate which is nonphysical for U=0U=0. For the case with U≠0U \neq 0

Fig. 9. Effect of fluid inertia on the variation of the maximum growth rate with WW in the case of a polymer melt ( ηr=0\eta_{r}=0, γ˙=103)\left.\dot{\gamma}=10^{3}\right).
there exist two roots for W>4W>4 in the range (K1,K2)\left(K_{1}, K_{2}\right), one negative and another positive value corresponding to large growth rates.

Figure 8 compares the maximum growth rates obtained by the dispersion relation including the inertial effects with the ones obtained under the noninertial approximation. Both the results match closely (ratio ∼1\sim 1 ) below W=4W=4 but beyond it the maximum growth rate predicted using the noninertial approximation diverges, and the ratio of the two growth rates falls to zero for W>4W>4. Figure 9 shows the variation in the maximum growth rate smax s_{\text {max }} with WW. The curve corresponding to U=0U=0 shows a divergence at W=4W=4, whereas, even for a small value of U=10−6U=10^{-6} the divergence is removed and the growth rate corresponds to a very large but finite value of smax s_{\text {max }}. For W<4W<4 the curves corresponding to U=10−6U=10^{-6} closely follows the one for U=0U=0 but above W≥4W \geq 4 it bifurcates and asymptotically converges to a large value for large WW. Therefore, we conclude that in the case of polymer melts above a critical value of WW, the growth rate becomes large and inertial terms have to be included in the analysis, whereas for W<4W<4 the assumption of neglecting inertia is valid.

It is instructive at this point to examine the physical reasons behind the divergence of the growth rate in the noninertial limit, and the removal of divergence upon the inclusion of inertia. First it is useful to inquire why the divergence happens precisely at W=4W=4 in the noninertial case. To answer this, we first note that for purely viscous liquid films, the maximum value of the nondimensional growth rate smax =1/4s_{\text {max }}=1 / 4. This can be written as smax =smax ∗T=1/4s_{\text {max }}=s_{\text {max }}^{*} T=1 / 4, where smax ∗s_{\text {max }}^{*} is the dimensional growth rate, and TT is the long time scale used to nondimensionalize s∗s^{*}. Also note that WW is defined as λ1/T\lambda_{1} / T. Without loss of generality, we can absorb the factor of 4 in the denominator in the definition of TT, thus resulting in smax⁡∗(4T)=1s_{\max }^{*}(4 T)=1. If we redefine our long time scale
as 4T4 T, then our definition of the Weissenberg number W~≡λ1/(4T)\tilde{W} \equiv \lambda_{1} /(4 T), and for this definition, the divergence in the growth rate will happen at W~≥1\tilde{W} \geq 1. Qualitatively, when W~>1\tilde{W}>1, the relaxation time in the liquid is larger than the viscous flow time, and in this regime, the deformation in the polymeric liquid occurs at time scales shorter than the relaxation time. This implies that when W~>1\tilde{W}>1, the polymeric liquid behaves more like an elastic solid. It is well known [38] that in a purely elastic solid, in the absence of inertia, the response of the material to applied forcing is instantaneous, and this solid-like behavior of the liquid film for W~>1\tilde{W}>1 manifests itself as a divergent growth rate (i.e. instantaneous growth of fluctuations in the material).

Of course, with the inclusion of inertia (Chapt. 3, [38]), an incompressible purely elastic solid admits shear waves with speed proportional to G/ρ\sqrt{G / \rho} (transverse speed of sound), where GG is the elastic modulus and ρ\rho is the density of the material. The shear wave propagation in a purely elastic material involves a time scale that is dictated by the shear wave speed, and this makes the response to applied forcings noninstantaneous, albeit very rapid. What happens in a polymeric liquid film with the presence of inertia is very similar for W~>1\tilde{W}>1, where the time scale for deformations is not governed by the viscous time TT, but by the propagation of shear waves in the viscoelastic liquid.

3.4 Effect of nonzero solvent viscosity (ηr)\left(\eta_{r}\right), with zero inertia

We now return to the dispersion relation in the absence of inertia, equation (35), and examine the growth rate in the limit ηr→0\eta_{r} \rightarrow 0. We recall that ηr\eta_{r} is the ratio of solvent to total viscosity of the polymer solution, and the limit ηr→0\eta_{r} \rightarrow 0 corresponds to the limiting case of a polymer melt. We find that in the limit ηr→0\eta_{r} \rightarrow 0, we obtain the maximum growth rate as s=1/(4−W)s=1 /(4-W) for W<4W<4 in agreement with the result for ηr=0\eta_{r}=0. However, for W=4W=4, the maximum growth rate is s=1/(4ηr)s=1 /\left(4 \sqrt{\eta_{r}}\right) which diverges as ηr−1/2\eta_{r}^{-1 / 2} as ηr\eta_{r} approaches zero (i.e. the limiting case of a polymer melt). For W>4W>4, the growth rate for ηr→0\eta_{r} \rightarrow 0 from equation (35) can be written as s=(W−4)/4ηrW+1/(W−4)s=(W-4) / 4 \eta_{r} W+1 /(W-4), which diverges as ηr−1\eta_{r}^{-1} as ηr\eta_{r} approaches zero. Beyond W=4W=4, the growth rate diverges at a faster rate than for W=4W=4. Thus, it emerges from our analysis that, in the absence of fluid inertia, the case of a polymer melt with ηr=0\eta_{r}=0 is a special singular limit, in the sense that even for very small (but finite) values of ηr\eta_{r}, the maximum value of the growth rate is a finite value for W=4W=4, but for ηr\eta_{r} identically equal to zero, we find a divergence of the maximum value of the growth rate. The physical reason for the removal of divergence by even a small amount of solvent is that the presence of solvent allows for an additional route for dissipation and flow.

Setting ηr=1\eta_{r}=1, we obtain the dispersion relation for a Newtonian fluid in the Jeffreys model multiplied by (1+(1+ Ws)W s) due to temporal derivatives (1+λ1∂t)\left(1+\lambda_{1} \partial_{t}\right) of the stress in the constitutive relation:

(s−(K2−K4))(1+Ws)=0\left(s-\left(K^{2}-K^{4}\right)\right)(1+W s)=0

Fig. 10. Variation of the growth rate with wave number KK for different values of ηr\eta_{r} with W=3W=3. The dashed line corresponds to s=−1/3s=-1 / 3. The growth rate asymptotically converges for large KK to s=−1/Ws=-1 / W for all ηr<1\eta_{r}<1.

Of the two growth rates obtained from the above dispersion relation, the one of importance is s=K2−K4s=K^{2}-K^{4} as the other corresponds to a negative constant (s=−1/W)(s=-1 / W). From the constitutive equation for Newtonian fluids we get the dispersion relation as s=K2−K4s=K^{2}-K^{4} which can be obtained by setting ηr=1\eta_{r}=1 and W=0W=0 in the Jeffreys model. The maximum growth rate for a Newtonian fluid is equal to 1/41 / 4 and it corresponds to the wave number K=1/2[10,15]K=1 / \sqrt{2}[10,15], and this wave number is the same as that obtained for polymer solutions (ηr≠1)\left(\eta_{r} \neq 1\right) in the present study. Figure 10 shows the variation in ss with KK for different values of ηr\eta_{r} (for W=3W=3 ). With an increase in ηr\eta_{r} the growth rate decreases due to increased viscous behavior of the system which opposes the instability. For ηr=1\eta_{r}=1, the growth rate does not converge to −1/W-1 / W as for other values of ηr\eta_{r}. Shorter wavelengths decay as (K2−K4)\left(K^{2}-K^{4}\right) similar to Newtonian fluids. Figure 11 shows a variation in ss with KK for different values of WW for ηr=0.5\eta_{r}=0.5. With an increase in WW, we see an increase in the maximum growth rate thus implying faster dynamics due to increased elastic behavior. For larger values of KK, the growth rate asymptotically converges to −1/W-1 / W. For W→∞W \rightarrow \infty, the growth rate variation with wave number is given by s=(K2−K4)/ηrs=\left(K^{2}-K^{4}\right) / \eta_{r}. The maximum growth rate corresponds to K=1/2K=1 / \sqrt{2} and is given by 1/4ηr1 / 4 \eta_{r}. Figure 12 shows the variation in the maximum growth rate with WW for different values of ηr\eta_{r}. The maximum growth rate corresponding to a particular WW decreases with an increase in ηr\eta_{r}. The maximum growth rate increases with an increase in WW, except for ηr=1\eta_{r}=1 (Newtonian fluid) for which the growth is independent of W(smax⁡=0.25)W\left(s_{\max }=0.25\right).

Thus, the key conclusions that emerge from our linear stability analysis are the following: i) The most dominant length scale in the linear analysis is unaffected by elastic effects in a viscoelastic fluid, and remains the same as in a Newtonian fluid. ii) The maximum growth rate increases

Fig. 11. Variation of the growth rate with wave number KK for different values of W.ηrW . \eta_{r} is 0.5 for all the cases.

Fig. 12. Variation in the maximum growth rate with the Weissenberg number WW for different ηr\eta_{r}.
with an increase in elastic effects in the fluid (signified by the parameter WW ), and for the special case of a noninertial polymer melt, it diverges when the nondimensional relaxation time W≥4W \geq 4. However, an inclusion of the inertial terms (proportional to UU ) in the analysis leads to finite but large growth rates for W≥4W \geq 4. The dominant wave number in case of polymer melts remains the same as for Newtonian films. Thus, based on our linear stability results, we expect that the elasticity in a polymeric liquid can change the dynamics of the instability, but should not have a significant effect on the length scales. iii) A careful examination of the growth rate for ηr→0\eta_{r} \rightarrow 0 in the absence of inertia indicated that even for very small (but finite) values of ηr\eta_{r}, there is no divergence of the maximum growth rate for W≥4W \geq 4, but the growth rate becomes large as ηr\eta_{r} approaches zero.

In the following section, thin-film equations are derived under the long-wave assumption and lubrication approximation in order to study the nonlinear evolution of the instability. While deriving the nonlinear equations, we consider explicitly the noninertial limit, but we also include the solvent contribution to the stress. Based on our above linear stability results, we expect our long-wave, nonlinear evolution equations to be valid for all values of WW if ηr>0\eta_{r}>0, and for W<4W<4 if ηr\eta_{r} is identically equal to zero.

4 Numerical solution of long-wave, nonlinear evolution equations

The nonlinear evolution equation for the height field h(x,t)h(x, t) of the free surface can be obtained from the governing equations and constitutive relations for a thin viscoelastic liquid film using standard methods [10] (also see App. B for a brief outline of our derivation). This nonlinear model is identical to that derived by Rauscher et al. [25], who linearized the nonlinear equations for understanding the stability of the shape of the rim that surrounds a growing hole. However, the issues related to the length and time scales of the surface instability as well as the numerical solutions of the nonlinear equations have not been previously addressed. Further, results from our linear stability of the complete governing equations, discussed above, also delineate the conditions for the applicability of the long-wave, inertialess, nonlinear model.

The governing partial differential equations in nondimensional parameters H(X,T),VˉXs(X,T)H(X, T), \bar{V}_{X_{s}}(X, T) and I(X,T)I(X, T), where X=x/L,Y=y/h0,T=t/(3e4(ηh0/γ))X=x / L, Y=y / h_{0}, T=t /\left(3 e^{4}\left(\eta h_{0} / \gamma\right)\right), H=h/h0,VˉXs(X,T)=ϵ3vxs/(γ/η)H=h / h_{0}, \bar{V}_{X_{s}}(X, T)=\epsilon^{3} v_{x_{s}} /(\gamma / \eta) and I=I= ∫0H(X,T)VˉX(Y)dY\int_{0}^{H(X, T)} \bar{V}_{X}(Y) \mathrm{d} Y obtained (in App. B), are

∂H∂T+3∂I(X,T)∂X=03ηrW∂I(X,T)∂T−3ηrWVˉXs∂H∂T+3I(X,T)=−W∂∂T(∂Ωˉ∂XH3)+3W2∂H∂T∂Ωˉ∂XH2−∂Ωˉ∂XH3ηrW∂VˉXs∂T+VˉXs=−W2∂∂T(∂Ωˉ∂XH2)−12∂Ωˉ∂XH2\begin{gathered} \frac{\partial H}{\partial T}+3 \frac{\partial I(X, T)}{\partial X}=0 \\ 3 \eta_{r} W \frac{\partial I(X, T)}{\partial T}-3 \eta_{r} W \bar{V}_{X_{s}} \frac{\partial H}{\partial T}+3 I(X, T)= \\ -W \frac{\partial}{\partial T}\left(\frac{\partial \bar{\Omega}}{\partial X} H^{3}\right)+\frac{3 W}{2} \frac{\partial H}{\partial T} \frac{\partial \bar{\Omega}}{\partial X} H^{2}-\frac{\partial \bar{\Omega}}{\partial X} H^{3} \\ \eta_{r} W \frac{\partial \bar{V}_{X_{s}}}{\partial T}+\bar{V}_{X_{s}}=-\frac{W}{2} \frac{\partial}{\partial T}\left(\frac{\partial \bar{\Omega}}{\partial X} H^{2}\right)-\frac{1}{2} \frac{\partial \bar{\Omega}}{\partial X} H^{2} \end{gathered}

The total pressure Ωˉ\bar{\Omega} including the short-range Born repulsion term is given by

Ωˉ=−∂2H∂X2+13H3(1−(lH)6)\bar{\Omega}=-\frac{\partial^{2} H}{\partial X^{2}}+\frac{1}{3 H^{3}}\left(1-\left(\frac{l}{H}\right)^{6}\right)

The first term in equation (42) is due to interfacial tension at the top surface, and the second term is due to the Lifshitz component of van der Waals interaction and a shortrange Born repulsion term. The parameter l=l0/h0l=l_{0} / h_{0}, where l0l_{0} is the cut-off distance at which the excess free energy has a minimum [10,11][10,11]. The dimensional long length scale in the xx-direction is h0(2πh02γ/A)1/2h_{0}\left(2 \pi h_{0}^{2} \gamma / A\right)^{1 / 2} and the viscous
time scale is 3(2πh02γ/A)2(ηh0/γ)3\left(2 \pi h_{0}^{2} \gamma / A\right)^{2}\left(\eta h_{0} / \gamma\right). We now perform a linear stability analysis under the long-wave approximation using equations (39-41) obtained above. The equations are linearized by perturbing the base state with a nondimensional wave number KK as

H=1+ζe(iKX+sT)I=I~e(iKX+sT)VˉX,I=VˉX,Ie(iKX+sT)\begin{aligned} H & =1+\zeta e^{(i K X+s T)} \\ I & =\tilde{I} e^{(i K X+s T)} \\ \bar{V}_{X, I} & =\bar{V}_{X, I} e^{(i K X+s T)} \end{aligned}

where ss is the nondimensional growth rate and ζ\zeta is the nondimensional amplitude of the perturbation in height. Substituting the parameters into governing equations, linearizing and eliminating I~\tilde{I} and VˉXs\bar{V}_{X_{s}}, equations (39-41) can be reduced to obtain the dispersion relation as (in a similar manner to the previous subsection)

ηrWs2+(1−W(K2−K4))s−(K2−K4)=0\eta_{r} W s^{2}+\left(1-W\left(K^{2}-K^{4}\right)\right) s-\left(K^{2}-K^{4}\right)=0

The parameters I~\tilde{I} and VˉXs\bar{V}_{X_{s}} can be obtained as

I~=sζe−iπ/23KVˉXs=ζ(K3−K)2(Ws+1)(Wηrs+1)e−iπ/2\begin{gathered} \tilde{I}=\frac{s \zeta e^{-i \pi / 2}}{3 K} \\ \bar{V}_{X_{s}}=\frac{\zeta\left(K^{3}-K\right)}{2} \frac{(W s+1)}{\left(W \eta_{r} s+1\right)} e^{-i \pi / 2} \end{gathered}

The dispersion relation so obtained from the long-wave analysis agrees exactly with the result obtained from the full dispersion relation in Section 3, thereby providing a validation of our derivation of long-wave evolution equations.

We next proceed to the discussion of results obtained by the numerical solution of the nonlinear evolution equations. The nondimensionalized coupled partial differential equations in H(X,T),VˉXs(X,T)H(X, T), \bar{V}_{X_{s}}(X, T) and I(X,T)I(X, T) (Eqs. (39-41)) are first discretized in space (X)(X) and for every discretized point in space, a set of coupled ordinary differential equations (ODEs) in time are obtained.

Fig. 13. Morphology evolution with a random initial perturbation on a simulation domain of 3Λ(ηr=0.5,W=4.0)3 \Lambda\left(\eta_{r}=0.5, W=4.0\right).

The ODEs thus obtained are solved with periodic boundary condition over an interval sufficiently large compared to the nondimensional wavelength of the fastest growing mode according to the linear theory (given by 22π2 \sqrt{2} \pi ). The initial conditions are chosen as small-amplitude (ζ=0.01)(\zeta=0.01) sinusoidal or random perturbations. We use Gear’s algorithm which is most suited for solving coupled, stiff ODEs to solve these set of ODEs. Four hundred grid points in a scale of 2Λ2 \Lambda were found to be sufficient for simulations, where Λ\Lambda is the most dangerous wavelength in the linear stability analysis. Figure 13 shows the evolution of instability, starting with a small-amplitude random perturbation of the free interface (ζ=0.01,ηr=0.5\left(\zeta=0.01, \eta_{r}=0.5\right. and W=4.0)\left.W=4.0\right), until the formation of a dry patch, i.e. the first time when the film surface touches the solid substrate. A domain of 3Λ3 \Lambda has been chosen for simulations when random initial perturbation are used. The initially random disturbances reorganize themselves on the length scale of the dominant wavelength (22π)(2 \sqrt{2} \pi). Thus, the length scale of domains obtained in nonlinear simulations is in agreement with predictions of the linear theory.

Figures 14(a) and (b) show the complete nonlinear evolution of the film from the formation and growth of holes leading to the formation of isolated droplets. The film is initially perturbed with a sinusoidal perturbation with an amplitude of 0.01 . The wavelength of the sinusoidal perturbation is taken as the dominant wavelength from the LSA. Simulations have been carried out over a length of two wavelengths with periodic boundary conditions. We observe that the evolution is symmetric, i.e. both the waves evolve together in exactly the same way. We have also verified that if the simulation is carried out in a domain size less than the critical wavelength (predicted by the linear theory) i.e. beyond which growth coefficient is negative, the perturbations do not grow, and after some time all the perturbations die out (data not shown). The initially perturbed interface (curve 1) in Figure 14(a) becomes unstable and the amplitude of perturbation increases with two crests and two troughs each of which will later correspond, respectively, to the formation of a drop and a hole. A rim forms rapidly around the expanding hole (curve 4 in Fig. 14(a)) and the rim height increases as the hole grows. Curves 5 and 6 show that as the rim starts growing it is influenced by the neighboring rims and both adjacent rims grow and coalesce. Curve 7 shows the valley formation just before the rim merger. Finally, the merger of rims takes place and they form a single drop. Dewetting after hole formation is much faster than the film break-up. After rupture, the time taken for the rim formation, collision and finally the formation of a single drop is much smaller. Curve 8 shows the final configuration of the film where the scales in XX and YY directions are different ( XX being much larger than YY ). From the above discussion, there are three distinct phases of evolution. In the first phase, the interface is unstable and the growth of instability ruptures the film. In the second phase, the dry spot formed grows in size and a rim forms around each hole. The rims grow in height and the adjacent rims later merge leading

Fig. 14. Evolution of the film with an initial sinusoidal perturbation (a) before rupture and (b) after rupture on a 2Λ2 \Lambda domain (ηr=0.5,W=0.001\left(\eta_{r}=0.5, W=0.001\right. and W=1)\left.W=1\right).
to the formation of a single drop. The width of the rim is observed to be less than the dominant wavelength, about one-half of the wavelength. In the third phase, the drops further grow in height in order to attain an equilibrium contact angle. It is useful to define a “rupture time” TrT_{r} as the time at which the film surface touches the bottom surface (H=l0)\left(H=l_{0}\right) and another time scale TeT_{e} is the time after which not more than 1%1 \% change in morphology is observed in simulations. With an increase in WW, there is a decrease in the rupture time (Fig. 15), thus showing that the elasticity enhances the instability along with the

Fig. 15. Variation in rupture time (relative to that of a Newtonian liquid) with WW for ηr=0.01,ηr=0.5\eta_{r}=0.01, \eta_{r}=0.5 and ηr=1.0\eta_{r}=1.0; Tr=13.1588T_{r}=13.1588 for a Newtonian fluid.

Fig. 16. Variation in TeT_{e} (relative to that of a Newtonian liquid) with WW for ηr=0.01,ηr=0.5\eta_{r}=0.01, \eta_{r}=0.5 and ηr=1.0;Te=13.2689\eta_{r}=1.0 ; T_{e}=13.2689 for a Newtonian fluid.
intermolecular interactions which engender the instability in thin films. Also, as ηr\eta_{r} is increased we observe an increase in the rupture time, thus showing the stabilizing effect of fluid viscosity. Compared to Newtonian liquids, the nondimensional rupture time is always smaller for viscoelastic solutions. The rupture time asymptotically converges to a constant value for large values of WW : Tr=6.597T_{r}=6.597 for ηr=0.5\eta_{r}=0.5 and W=1000W=1000 and Tr=6.589T_{r}=6.589 for W=3000W=3000 (not shown in the figure). For these cases the (nondimensional) time of rupture as predicted by the linear stability analysis is 9.2288 and 9.2165 , respectively. A similar trend is observed for the corresponding TeT_{e} relative to the TeT_{e} for a Newtonian fluid (Fig. 16).

Figure 17 shows a comparison between the results from the simulation for W=2W=2 and W=8W=8 with a low value of ηr=0.01\eta_{r}=0.01. For W=2W=2, the growth of instability is much slower than that for W=8W=8. At nondimensional T=T=

Fig. 17. Growth of instability for W=2W=2 and W=8W=8 with ηr=0.01\eta_{r}=0.01 (close to a polymer melt). Curves with symbols are data points from simulations, while curves without symbols are results from linear stability analysis.

Fig. 18. Comparison between growth of instability for W=W= 0.1,W=20.1, W=2 and W=6.0W=6.0 with ηr=0.5\eta_{r}=0.5. A comparison with linear stability results is also shown along with nonlinearsimulation results
0.4 , when the film with W=2W=2 is just starting to evolve ((1−Hmin⁡)=0.015)\left(\left(1-H_{\min }\right)=0.015\right), the one with W=8W=8 is near to its rupture ((1−Hmin⁡)=0.64)\left(\left(1-H_{\min }\right)=0.64\right). For polymeric solutions with very low solvent component, growth rates are quite high and are highly sensitive to WW. In this figure, continuous curves without symbols correspond to prediction from the linear stability analysis.

Figure 18 shows a comparison between LSA predictions and nonlinear simulations for three different WW with ηr=0.5\eta_{r}=0.5. Nonlinear simulations agree with linear stability predictions in the early times but after some time, the growth rate increases further and deviates from LSA results, as would be expected.

Fig. 19. Variation of the rim height with the hole diameter for W=0.001W=0.001 and W=1.0(ηr=0.5)W=1.0\left(\eta_{r}=0.5\right). This variation is for the hindered growth of rims in which two adjoining rims later coalesce to form a single drop. In the present discussion rim heights are considered even after coalescence.

Fig. 20. Variation of the rim height with the hole diameter for ηr=0.5\eta_{r}=0.5 and ηr=1.0(W=1)\eta_{r}=1.0(W=1). This variation is for the hindered growth of rims in which two adjoining rims later coalesce to form a single drop. In the present discussion rim heights are considered even after coalescence.

Figure 19 shows the variation in the rim height with the hole diameter for W=0.001W=0.001 and W=1.0W=1.0 with ηr=0.5\eta_{r}=0.5. For both the cases the curves overlap each other thus showing that the evolution of morphology in both cases is identical. Figure 20 shows similarly the variation of the rim height with the hole diameter for both ηr=0.5\eta_{r}=0.5 and ηr=1\eta_{r}=1 (Newtonian fluid) with W=1.0W=1.0. This suggests that both WW and ηr\eta_{r} have little effect on the details of the dynamical evolution of the dewetting film, indicating that viscoelasticity does not play a major role in the morphological evolution of the dewetting film. It also

Fig. 21. Growth of the hole diameter with rescaled time (T−Tr)/(Te−Tr)\left(T-T_{r}\right) /\left(T_{e}-T_{r}\right), where TrT_{r} is the rupture time and TeT_{e} is the final time, when the system has attained a time-invariant metastable state and beyond which not more than 1%1 \% change in the evolution is observed in simulations. Data for ηr=0.5\eta_{r}=0.5, W=4W=4 and ηr=1.0\eta_{r}=1.0 (Newtonian liquid).
shows the two distinct phases of hole growth. In the first phase, two adjacent rims grow together symmetrically hindering each other’s lateral movement. While growing they form a valley in-between, which as they are pushed against each other grows in height. In the later stages, marked by a sudden increase in slope (Figs. 19 and 20), the rims coalesce and form a single drop. This marks the second phase of hole growth in which the drops ripen towards a metastable state with an equilibrium contact angle. In the first phase, the rim height varies as DmD^{m} with 0<m<10<m<1 (where DD is the hole diameter) while in the second phase it varies linearly. The value of mm is found to be 0.525 by a curve-fitting procedure. Before the coalescence of the rims, each hole corresponds to two rims, whereas after the rim merger each drop corresponds to one growing hole. Due to volume constraints (conservation of mass) a growth of hole corresponds to an increase in the height of the drop.

Figure 21 shows the variation in the hole diameter with a rescaled normalized nondimensional time θ=\theta= (T−Tr)/(Te−Tr)\left(T-T_{r}\right) /\left(T_{e}-T_{r}\right) for ηr=0.5\eta_{r}=0.5 and ηr=1.0\eta_{r}=1.0 (Newtonian liquid). The times TrT_{r} and TeT_{e} for the case with ηr=0.5\eta_{r}=0.5, W=4.0W=4.0 are obtained from the simulations as 9.0202 and 9.0722 , respectively. For the case ηr=1.0,Tr=13.1642\eta_{r}=1.0, T_{r}=13.1642 and Te=13.2689T_{e}=13.2689. The hole diameter is observed to scale as θq\theta^{q}, where q=0.75q=0.75 initially and increases gradually to q=1q=1. Remarkably, the variation of the hole diameter with the rescaled time collapses onto the same curve for both Newtonian liquids and polymer solutions, thus indicating that even the dynamical evolution is identical with respect to the rescaled time. Experiments on growths of holes in thin polymer films also show a variation in the hole diameter as D∝(T−Tr)qD \propto\left(T-T_{r}\right)^{q}, where q∼0.7q \sim 0.7 [19].

A central feature that emerges from our analysis of the Jeffreys model is the insensitivity of the length

scale of instability to the rheological details, even though these are very important for the dynamics. This basic conclusion may also hold for more complicated rheological models (such as the FENE-P model [29]), details of which however need to be investigated. Certainly, the precise morphology and dynamics are likely to change for more complicated rheology. However, whether these effects can be “renormalized” to the Newtonian behavior by an appropriate rescaling of the dynamics remains an open question. It should be noted that the above remarks relate only to viscoelastic fluids that lack a permanent or zero-frequency elastic modulus. Indeed, viscoelastic solids are expected to even exhibit a different length scale of instability due permanent elasticity [34,37]. In particular, for the case of a nearly elastic solid, the length scale of instability shifts to the short-wave regime also independent of the shear modulus.

5 Conclusions

We have studied the surface instability, its length and time scales, and morphology evolution in an unstable, thin, viscoelastic liquid film dewetting under the influence of long-range van der Waals attractive force, using both linear stability analysis and numerical solutions of a set of nonlinear evolution equations. The viscoelastic nature of the liquid was modeled using both Maxwell and Jeffreys class of models, which incorporate elasticity of polymeric liquids through a characteristic relaxation time. We have carried out a comprehensive linear stability analysis valid for arbitrary wave numbers, as well as nonlinear numerical simulations in order to assess the roles of fluid viscoelasticity and nonlinearities in the processes of the growth of instability, film rupture and dewetting by hole growth.

Our results show that the length scale of instability is completely unaltered by the rheological details of viscoelasticity, thus showing that the length scale of patterns is governed entirely by thermodynamic factors - surface tension and the destabilizing interactions, such as van der Waals force, rather than by transport processes that are controlled by the fluid viscoelasticity. The latter, however, critically influence the time scale and the dynamics including the early growth of instability and subsequent dewetting, both of which are accelerated by increasing fluid elasticity. For the Maxwell fluid (which describes a liquid polymer melt) under the zero-inertia limit, the growth rate of the most unstable wave diverges in a window of wave numbers which enlarges when the Weissenberg number WW becomes increasingly larger than a critical value, W=4W=4. A similar divergence in a Maxwell liquid film subjected to electric fields was predicted in reference [28]. We demonstrated that in such cases the inertial terms cannot be ignored. Subsequent analysis including fluid inertia showed that even for extremely small values of a nondimensional group UU characterizing fluid inertia, the singularity of the growth rate disappears. The maximum growth rate corresponds to the scaled wave number K=0.707K=0.707. In the absence of inertia, the aforementioned singularity can also be removed by inclusion of solvent viscosity as in the Jeffreys
model, wherein the solvent dynamics provides an additional mode of dissipation. However, the maximum growth rate remains large and scales as smax⁡∼1/ηr1/2s_{\max } \sim 1 / \eta_{r}^{1 / 2} for W=4W=4 and smax⁡∼1/ηrs_{\max } \sim 1 / \eta_{r} for W>4W>4 in the limit of a polymer melt ηr→0\eta_{r} \rightarrow 0. The ultrafast growth of fluctuations results from a predominantly solid-like elastic response where viscous dissipation becomes ineffective at short time scales and inertial forces become dominant. This fast growth has interesting implications for rapid and precise patterning of a polymer melt film in soft lithography.

Interestingly, results from our nonlinear simulations indicate that some of the most readily observed morphological details such as the rim height sss s. hole diameter are invariant to the rheological details of the liquid, by an appropriate rescaling of the dynamics. Even though the kinetics of hole growth is greatly influenced by fluid elasticity, the scaling of the hole diameter with a rescaled time θ,D∼θq\theta, D \sim \theta^{q}, also remains unaltered with elasticity. Our results thus underscore the precise role of viscoelasticity in the morphology and dynamics of a dewetting polymeric film and help in the interpretation of experimental observations by suggesting which features can be attributed to bulk rheology as opposed to other factors such as wall slippage.

Appendix A.

The complete dispersion relation alluded to in Section 3.3 is given by the determinant

F(s,k,W,U,γˉ)=∣1111−q−kqka31a32a33a34a41a42a43a44∣F(s, k, W, U, \bar{\gamma})=\left|\begin{array}{cccc} 1 & 1 & 1 & 1 \\ -q & -k & q & k \\ a_{31} & a_{32} & a_{33} & a_{34} \\ a_{41} & a_{42} & a_{43} & a_{44} \end{array}\right|

where the variables aija_{i j} are given by

a31=ekh0(k2+q2)h02a32=2k2h02ekh0a33=e(2q+k)h0(k2+q2)h02a34=2e(q+2k)h0k2h02a41=ekh0(P1−P2+P3−qh0P4)a42=eqh0(−P5−P2+P3−kh0P4)a43=e(2q+k)h0(−P1−P2+P3+qh0P4)a44=e(q+2k)h0(P5−P2+P3+kh0P4)\begin{aligned} & a_{31}=e^{k h_{0}}\left(k^{2}+q^{2}\right) h_{0}^{2} \\ & a_{32}=2 k^{2} h_{0}^{2} e^{k h_{0}} \\ & a_{33}=e^{(2 q+k) h_{0}}\left(k^{2}+q^{2}\right) h_{0}^{2} \\ & a_{34}=2 e^{(q+2 k) h_{0}} k^{2} h_{0}^{2} \\ & a_{41}=e^{k h_{0}}\left(P_{1}-P_{2}+P_{3}-q h_{0} P_{4}\right) \\ & a_{42}=e^{q h_{0}}\left(-P_{5}-P_{2}+P_{3}-k h_{0} P_{4}\right) \\ & a_{43}=e^{(2 q+k) h_{0}}\left(-P_{1}-P_{2}+P_{3}+q h_{0} P_{4}\right) \\ & a_{44}=e^{(q+2 k) h_{0}}\left(P_{5}-P_{2}+P_{3}+k h_{0} P_{4}\right) \end{aligned}

Here, P1=γˉqh0(−3k2+q2)h02sˉ,P2=3γˉ2k2h02(1+sˉW)P_{1}=\bar{\gamma} q h_{0}\left(-3 k^{2}+q^{2}\right) h_{0}^{2} \bar{s}, P_{2}=3 \bar{\gamma}^{2} k^{2} h_{0}^{2}(1+\bar{s} W), P3=3γˉ3k4h04(1+sˉW),P4=sˉ2U(1+sˉW)P_{3}=3 \bar{\gamma}^{3} k^{4} h_{0}^{4}(1+\bar{s} W), P_{4}=\bar{s}^{2} U(1+\bar{s} W) and P5=P_{5}= 2γˉk3h03sˉ2 \bar{\gamma} k^{3} h_{0}^{3} \bar{s}.

Appendix B.

In the long-wave limit, a lubrication-type approximation [10] can be employed to derive the evolution equations for the free-surface position. This appendix provides

a brief outline of our derivation. The xx-momentum equation then yields, in the long-wave limit,

∂Π∂x=∂τxy∂y\frac{\partial \Pi}{\partial x}=\frac{\partial \tau_{x y}}{\partial y}

and the yy-momentum equation reduces to

∂Π∂y=0\frac{\partial \Pi}{\partial y}=0

Introducing nondimensional parameters using H=h/h0H=h / h_{0}, T∗=t/(ηh0/γ),V∗=v/(γ/η),τ∗=τ/(γ/h0),Π∗=T^{*}=t /\left(\eta h_{0} / \gamma\right), \mathbf{V}^{*}=\mathbf{v} /(\gamma / \eta), \boldsymbol{\tau}^{*}=\boldsymbol{\tau} /\left(\gamma / h_{0}\right), \Pi^{*}= Π/(γ/h0),y∗=y/h0\Pi /\left(\gamma / h_{0}\right), y^{*}=y / h_{0} and x∗=x/h0x^{*}=x / h_{0}. Henceforth, the ∗∗∗* * * superscript from all nondimensional variables will be dropped in order to keep the notation simple. The constitutive relation for τxy\tau_{x y} in terms of the above-defined nondimensional parameters under the long-wave assumption can be written as

Wi∂τxy∂T+τxy=Wiηr∂∂T(∂Vx∂y)+∂Vx∂yW i \frac{\partial \tau_{x y}}{\partial T}+\tau_{x y}=W i \eta_{r} \frac{\partial}{\partial T}\left(\frac{\partial V_{x}}{\partial y}\right)+\frac{\partial V_{x}}{\partial y}

where WiW i is the Weissenberg number defined as

Wi=(λ1γηh0)W i=\left(\frac{\lambda_{1} \gamma}{\eta h_{0}}\right)

The normal stress balance at h(x)h(x) (Eq. (8)) in terms of nondimensional parameters can be written as

Π=−∂2H∂x2+A2πh0313H3\Pi=-\frac{\partial^{2} H}{\partial x^{2}}+\frac{A}{2 \pi h_{0}^{3}} \frac{1}{3 H^{3}}

and the kinematic equation (Eq. (5)) can be written as

∂H∂T+Vxs∂H∂x=Vys\frac{\partial H}{\partial T}+V_{x_{s}} \frac{\partial H}{\partial x}=V_{y_{s}}

where VxsV_{x_{s}} and VysV_{y_{s}} are nondimensional velocities and xx and yy direction velocities at the free surface.

Equation (B.2) implies that Π\Pi is independent of yy and is a function only of xx and TT, i.e. Π(x,T)\Pi(x, T). The xx momentum equation (Eq. (B.1)) can be integrated to give

τxy=dΠ dx[y−H(x,t)]\tau_{x y}=\frac{\mathrm{d} \Pi}{\mathrm{~d} x}[y-H(x, t)]

where the constant of integration is fixed using the freesurface boundary condition (Eq. (7)). Substituting τxy\tau_{x y} from the above equation into the constitutive relation (Eq. (B.3)), we get

Wi∂∂T( dΠ dx(y−H(x,T)))+(dΠ dx(y−H(x,T)))=∂∂Y(Wiηr∂Vx∂T+Vx)\begin{aligned} & W i \frac{\partial}{\partial T}\left(\frac{\mathrm{~d} \Pi}{\mathrm{~d} x}(y-H(x, T))\right)+\left(\frac{\mathrm{d} \Pi}{\mathrm{~d} x}(y-H(x, T))\right)= \\ & \quad \frac{\partial}{\partial Y}\left(W i \eta_{r} \frac{\partial V_{x}}{\partial T}+V_{x}\right) \end{aligned}

Further integrating the above equation with respect to yy we get

Wi∂∂T( dΠ dx(y22−yH(x,T)))+(dΠ dx(y22−yH(x,T)))=(Wiηr∂Vx∂T+Vx)\begin{aligned} & W i \frac{\partial}{\partial T}\left(\frac{\mathrm{~d} \Pi}{\mathrm{~d} x}\left(\frac{y^{2}}{2}-y H(x, T)\right)\right) \\ & +\left(\frac{\mathrm{d} \Pi}{\mathrm{~d} x}\left(\frac{y^{2}}{2}-y H(x, T)\right)\right)=\left(W i \eta_{r} \frac{\partial V_{x}}{\partial T}+V_{x}\right) \end{aligned}

Rewriting the normal stress balance equation (Eq. (B.5)) as

Π2πh02γA=−2πh02γA∂2H∂x2+13H3\Pi \frac{2 \pi h_{0}^{2} \gamma}{A}=-\frac{2 \pi h_{0}^{2} \gamma}{A} \frac{\partial^{2} H}{\partial x^{2}}+\frac{1}{3 H^{3}}

and rescaling xx as X=(A/2πh02γ)1/2xX=\left(A / 2 \pi h_{0}^{2} \gamma\right)^{1 / 2} x, we obtain the “long length scale” LL of the problem as

L=(2πh04γA)1/2L=\left(\frac{2 \pi h_{0}^{4} \gamma}{A}\right)^{1 / 2}

The nondimensional pressure at the top surface can now be written as

Πˉ=−∂2H∂X2+13H3\bar{\Pi}=-\frac{\partial^{2} H}{\partial X^{2}}+\frac{1}{3 H^{3}}

where Π=ϵ2Πˉ,X=ϵx,ϵ=h0/L\Pi=\epsilon^{2} \bar{\Pi}, X=\epsilon x, \epsilon=h_{0} / L and Y=yY=y. In terms of these rescaled parameters equation (B.9) can now be written as

Wi∂∂T(ϵ3∂Πˉ∂X(Y22−YH))+ϵ3∂Πˉ∂X(Y22−YH)=ηrWi∂Vx∂T+Vx\begin{aligned} W i \frac{\partial}{\partial T} & \left(\epsilon^{3} \frac{\partial \bar{\Pi}}{\partial X}\left(\frac{Y^{2}}{2}-Y H\right)\right)+\epsilon^{3} \frac{\partial \bar{\Pi}}{\partial X}\left(\frac{Y^{2}}{2}-Y H\right)= \\ & \eta_{r} W i \frac{\partial V_{x}}{\partial T}+V_{x} \end{aligned}

Evidently, from equation (B.13), the nondimensional VxV_{x} can be rescaled as Vx=Vˉxϵ3V_{x}=\bar{V}_{x} \epsilon^{3}. Substituting VxsV_{x_{s}} with the rescaled Vˉxs\bar{V}_{x_{s}} into the kinematic condition we get

∂H∂T+VˉXsϵ4∂H∂X=VVx\frac{\partial H}{\partial T}+\frac{\bar{V}_{X_{s}}}{\epsilon^{4}} \frac{\partial H}{\partial X}=V_{V_{x}}

From the above equation, TT and VyV_{y} can be rescaled as Tˉ=\bar{T}= T/ϵ4T / \epsilon^{4} and VˉY=ϵ4VY\bar{V}_{Y}=\epsilon^{4} V_{Y} and defining I=∫0H(X,Tˉ)VˉX(Y)dYI=\int_{0}^{H(X, \bar{T})} \bar{V}_{X}(Y) \mathrm{d} Y as the total volume flux per unit width passing across a plane at a location XX, we can rewrite the kinematic equation as

∂H∂Tˉ+∂I(X,Tˉ)∂X=0\frac{\partial H}{\partial \bar{T}}+\frac{\partial I(X, \bar{T})}{\partial X}=0

Writing equation (B.13) in terms of the above rescaled parameters and evaluating at Y=H(X,Tˉ)Y=H(X, \bar{T}), we get a differential equation for the interfacial flow velocity (VˉXs)\left(\bar{V}_{X_{s}}\right) :

ηrWˉi∂VˉXs∂Tˉ+VˉXs=−Wˉi2∂∂Tˉ(∂Πˉ∂XH2)−12∂Πˉ∂XH2\eta_{r} \bar{W} i \frac{\partial \bar{V}_{X_{s}}}{\partial \bar{T}}+\bar{V}_{X_{s}}=-\frac{\bar{W} i}{2} \frac{\partial}{\partial \bar{T}}\left(\frac{\partial \bar{\Pi}}{\partial X} H^{2}\right)-\frac{1}{2} \frac{\partial \bar{\Pi}}{\partial X} H^{2}

Integrating equation (B.13) with respect to YY and using Leibnitz’s theorem, we obtain

ηrWˉi∂I(X,T)∂Tˉ−ηrWˉiVˉXs∂h∂T+I(X,T)=−Wˉi3∂∂T(∂Πˉ∂XH3)+Wˉi2∂H∂T∂Πˉ∂XH2−13∂Πˉ∂XH3\begin{aligned} \eta_{r} \bar{W} i \frac{\partial I(X, T)}{\partial \bar{T}} & -\eta_{r} \bar{W} i \bar{V}_{X_{s}} \frac{\partial h}{\partial T}+I(X, T)= \\ & -\frac{\bar{W} i}{3} \frac{\partial}{\partial T}\left(\frac{\partial \bar{\Pi}}{\partial X} H^{3}\right) \\ & +\frac{\bar{W} i}{2} \frac{\partial H}{\partial T} \frac{\partial \bar{\Pi}}{\partial X} H^{2}-\frac{1}{3} \frac{\partial \bar{\Pi}}{\partial X} H^{3} \end{aligned}

Without loss of generality, we can define W=Wˉi/3W=\bar{W} i / 3 and T=Tˉ/3T=\bar{T} / 3 and rewrite the equations above to obtain three nonlinear coupled partial differential equations for H(X,T),VˉXs(X,T)H(X, T), \bar{V}_{X_{s}}(X, T) and I(X,T)I(X, T) (Eqs. (39-41)).

References

E. Ruckenstein, R.K. Jain, Chem. Soc. Faraday Trans. 270, 132 (1974).
R.K. Jain, E. Ruckenstein, J. Colloid Interface Sci. 54, 108 (1976).
M.B. Williams, S.H. Davis, J. Colloid Interface Sci. 90, 220 (1982).
S.G. Bankoff, S.H. Davis, Physicochem. Hydrodyn. 9, 5 (1987).
C. Redon, F. Brochard-Wyart, F. Rondelez, Phys. Rev. Lett. 66, 715 (1991).
V.S. Mitlin, J. Colloid Interface Sci. 156, 491 (1993).
V.S. Mitlin, N.V. Petviashvili, Phys. Lett. A 192, 323 (1994).
F. Brochard-Wyart, P.G. de Gennes, H. Hervert, C. Redon, Langmuir 10, 1566 (1994).
P. Lambooy, K.C. Phelan, O. Haugg, G. Krausch, Phys. Rev. Lett. 76, 1110 (1996).
A. Oron, S.H. Davis, S.G. Bankoff, Rev. Mod. Phys. 69, 931 (1997).
A. Sharma, R. Khanna, J. Colloid Interface Sci. 195, 42 (1997).
T.G. Stange, D.F. Evans, W.A. Hendrickson, Langmuir 13, 4459 (1997).
A. Ghatak, R. Khanna, A. Sharma, J. Colloid Interface Sci. 212, 483 (1999).
C.K. Lin, C.C. Hwang, W.Y. Uen, J. Colloid Interface Sci. 231, 379 (2000).
A. Sharma, Eur. Phys. J. E 12, 397 (2003).
J. Becker, G. Grün, R. Seemann, H. Mantz, K. Jacobs, K.R. Mecke, R. Blossey, Nature Mater. 2, 59 (2003).
U. Thiele, Eur. Phys. J. E 12, 409 (2003).
A. Sharma, G. Reiter, J. Colloid Interface Sci. 178, 383 (1996).
G. Reiter, P. Auroy, L. Auvray, Macromolecules 29, 2150 (1996).
C. Neto, K. Jacobs, R. Seemann, R. Blossey, J. Becker, G. Grün, J. Phys.: Condens. Matter 15, S421 (2003).
A. Sathyagal, G. Narsimhan, Chem. Eng. Commun. 111, 161 (1992).
Y.L. Zhang, O.K. Matar, R.V. Craster, J. Colloid Interface Sci. 262, 130 (2003).
S. Herminghaus, R. Seemann, K. Jacobs, Phys. Rev. Lett. 89, 056101 (2002).
S. Herminghaus, K. Jacobs, R. Seemann, Eur. Phys. J. E 12, 101 (2003).
M. Rauscher, A. Münch, B. Wagner, R. Blossey, Eur. Phys. J. E 17, 373 (2005).
S.A. Safran, J. Klein, J. Phys. II 3, 749 (1993).
G. Narsimhan, Z. Wang, J. Colloid Interface Sci. 291, 296 (2005).
L. Wu, S.Y. Chou, J. Non-Newtonian Fluid Mech. 125, 91 (2004).
R.B. Bird, R.C. Armstrong, O. Hassager, Dynamics of Polymeric Liquids, Vol. 1 (Wiley & Sons, New York, 1987).
R. Larson, Structure and Rheology of Complex Fluids (Oxford, 1998).
Y.L. Zhang, O.K. Matar, R.V. Craster, J. Non-Newtonian Fluid Mech. 105, 53 (2002).
S. Kumar, O.K. Matar, J. Colloid Interface Sci. 273, 581 (2004).
O.K. Matar, V. Gkanis, S. Kumar, J. Colloid Interface Sci. 286, 319 (2005).
A. Ghatak, M.K. Chaudhury, V. Shenoy, A. Sharma, Phys. Rev. Lett. 85, 4329 (2000).
G. Reiter, Phys. Rev. Lett. 87, 186101 (2001).
W. Monch, S. Herminghaus, Europhys. Lett. 52, 525 (2001).
V. Shenoy, A. Sharma, J. Mech. Phys. Solids 50, 1155 (2002).
L.D. Landau, E.M. Lifshitz, Theory of Elasticity, Vol. VII (Butterworth, London, 1995).