1. Introduction
The dynamics of unstable waves developing on a viscous film flowing down an inclined plane is a classical problem in fluid mechanics, often referred to as the Kapitza problem (e.g. Chang Reference Chang1994). The instability features in a number of everyday settings, including water flow down roads and windscreens on rainy days. Analogous instabilities appear on flowing films of complex fluids, and thereby used to rationalize natural phenomena such as mud surges (e.g. Ng & Mei Reference Ng and Mei1994; Coussot Reference Coussot1997). In this latter context, the instability has been used to illustrate some of the novelties introduced by non-Newtonian behaviour, such as inertialess elastic instability (Shaqfeh, Larson & Fredrickson Reference Shaqfeh, Larson and Fredrickson1989), stabilization by a yield stress (Balmforth & Liu Reference Balmforth and Liu2004) and discontinuous shear thickening (Darbois Texier et al. Reference Darbois Texier, Lhuissier, Metzger and Forterre2023; Balmforth Reference Balmforth2025). Worded differently, the wavy instabilities on a falling film of complex fluids offer a probe into material behaviour: the spatiotemporal dynamics of the film provides a different setting than simple rheometric flows, and can be compared with relatively straightforward laboratory experiments. A similar vein has been mined for granular media (Edwards & Gray Reference Edwards and Gray2015). A goal of the current paper is to extend these ideas to thixotropic fluids.
Thixotropic fluids have a time-dependent microstructure that builds up when the fluid is at rest, leading to an increase in effective viscosity, but is reversibly broken down by flow, thereby lowering the fluid’s resistance (Mewis & Wagner Reference Mewis and Wagner2009). This behaviour can lead to pronounced hysteresis, or ‘thixotropic loops’, in rheometers when stresses or strain rates are cycled up and down. A wide range of fluids exhibit thixotropic behaviour, including natural clay suspensions, industrial drilling fluids and cements, printing inks and paints, oils and grease, and food products such as mayonnaise and ketchup (Barnes Reference Barnes1997; Mewis & Wagner Reference Mewis and Wagner2009).
Thixotropic hysteresis can be particularly pronounced when the material behaviour features so-called ‘viscosity bifurcations’ (e.g. Coussot et al. Reference Coussot, Nguyen, Huynh and Bonn2002; Dullaert & Mewis Reference Dullaert and Mewis2006; Møller et al. Reference Møller, Mewis and Bonn2006; Alexandrou, Constantinou & Georgiou Reference Alexandrou, Constantinou and Georgiou2009; Wachs, Vinay & Frigaard Reference Wachs, Vinay and Frigaard2009; Dimitriou & McKinley Reference Dimitriou and McKinley2014). One such transition arises when the fabric of the fluid evolves relatively slowly towards a structured state that is highly viscous or ‘jammed’ (solid-like), but collapses abruptly under shear. Dramatic failure events can then occur at the initiation of flow. A second type of viscosity bifurcation arises once motion is underway and the microstructure broken down, if that structure is able to swiftly recover and jam again below a critical shear rate. This behaviour has been observed in a number of fluids, and implies that, should stresses decline sufficiently to permit shear rates to descend below the critical shear rate, the viscosity abruptly increases and flow becomes arrested.
Hewitt & Balmforth (Reference Hewitt and Balmforth2013) previously explored the dynamics of thixotropic fluids with viscosity bifurcations for inertialess gravity currents (see also Oishi, Martins & Thompson Reference Oishi, Martins and Thompson2017). In this setting, the sudden release and collapse of a reservoir seeds the runout of the current towards a final deposit. During runout, the thixotropic rheology of the material leads to a characteristic two-layer structure to the flow, with lower-lying destructured material conveying along a structured top layer. Under certain conditions, Hewitt & Balmforth (Reference Hewitt and Balmforth2013) observed that the interface between the layers developed what looked like wavy instabilities. This observation was rationalized by drawing an analogy with the instabilities that appear in multilayer, viscously stratified or non-Newtonian flows (e.g. Loewenherz & Lawrence Reference Loewenherz and Lawrence1989; Loewenherz, Lawrence & Weaver Reference Loewenherz, Lawrence and Weaver1989; Chen Reference Chen1993, Reference Chen1995; Charru & Hinch Reference Charru and Hinch2000; Balmforth, Craster & Toniolo Reference Balmforth, Craster and Toniolo2003; Ern, Charru & Luchini Reference Ern, Charru and Luchini2003). The idea that instability can be triggered by the deformation of rheologically induced interfaces has also featured in the discussion of other shear flows of complex fluids (e.g. Renardy Reference Renardy1995; Fielding Reference Fielding2005, Wilson & Fielding Reference Wilson and Fielding2006; Fielding Reference Fielding2007; Nghe et al. Reference Nghe, Fielding, Tabeling and Ajdari2010; Cromer, Cook & McKinley Reference Cromer, Cook and McKinley2011; Nicolas & Morozov Reference Nicolas and Morozov2012; Renardy & Renardy Reference Renardy and Renardy2017; Castillo & Wilson Reference Castillo and Wilson2020).
Nevertheless, in order to attack the gravity current problem, Hewitt & Balmforth (Reference Hewitt and Balmforth2013) significantly simplified the rheological model that they employed. In particular, they assumed that the microstructure adapted instantaneously to local shear rates. As a consequence, wherever stresses prompted a viscosity bifurcation, it was assumed that the material instantaneously switched state. Without spatial diffusion, the border between structured and destructured fluid then becomes a jump discontinuity that evolves as a material interface. The limitations inherent in this approximation (which were necessary to simplify the gravity current model), lead one to question whether interfacial instability can arise in low-Reynolds-number thixotropic flows.
The purpose of the present article is to reconsider the linear stability of inertialess, thin films of thixotropic material, with a view to addressing this question. We therefore return to the original rheological model proposed by Hewitt & Balmforth (Reference Hewitt and Balmforth2013), which aligns with many others used in the literature on complex fluids (Coussot et al. Reference Coussot, Nguyen, Huynh and Bonn2002; Dullaert & Mewis Reference Dullaert and Mewis2006; Møller et al. Reference Møller, Mewis and Bonn2006; Alexandrou et al. Reference Alexandrou, Constantinou and Georgiou2009). We thereby accommodate the time-dependent relaxation of the microstructure and also its spatial diffusion, formulating a model suitable for a thin-film setting (§ 2). Armed with this model, we construct steady base states corresponding to sheet-like flow down an inclined plane (§ 3), and then test their linear stability (§ 4). We find that the instability detected by Hewitt & Balmforth (Reference Hewitt and Balmforth2013) can indeed feature, and interrogate its dynamics more thoroughly, including at the finite amplitudes where instability saturates (§ 5). For completeness, Appendix A presents a summary of the Hewitt–Balmforth reduction and its predictions for the current problem of the wavy instabilities on a falling fluid film.
2. Shallow flow model
2.1. Dimensional formulation
Consider a complex fluid flowing down a plane inclined at angle
$\theta$
. The geometry is described by a Cartesian coordinate system
$({\hat {x}},{\hat {z}})$
, orientated such that the
$x$
-axis points down the slope, see figure 1. The fluid velocity is
$\boldsymbol{{\hat {u}}} = ({\hat {u}},{\hat {w}})$
. The fluid is shallow, with a mean depth
$\mathcal H$
that is much smaller than the characteristic length scale for variations over the plane,
$\mathcal L$
. The slope angle is also taken to be relatively small, leading us to set
\begin{equation} \tan \theta = \varepsilon = \frac {{\mathcal H}}{{\mathcal L}} \ll 1 . \end{equation}
In this setting, we consider waves evolving over spatially periodic domains of length of
$2\pi {\mathcal L}/k$
, where
$k$
is an
$O(1)$
dimensionless wavenumber. The local fluid depth is
${\hat {z}} = \hat {h}({\hat {x}},{\hat {t}}\kern1pt)$
.
Geometry of the thin-film model: a shallow layer of thixotropic fluid, with viscosity dictated by a structure function
$\lambda ({\hat {x}},{\hat {z}},{\hat {t}}\kern1pt)$
, flows down an incline of slope
$\tan \theta =\varepsilon ={\mathcal H}/{\mathcal L}$
, where
$\mathcal H$
is the mean depth.
Assuming that the fluid is incompressible and inertia is negligible, conservation of mass and momentum demand that
\begin{equation} \frac {\partial {\hat {u}}}{\partial {\hat {x}}} + \frac {\partial {\hat {w}}}{\partial {\hat {z}}} = 0, \end{equation}
and
\begin{equation} \textbf {0} = \rho \tilde {\boldsymbol{g}} - \boldsymbol{\nabla }{\hat {p}} + \boldsymbol{\nabla }\boldsymbol{\cdot } \boldsymbol{\hat \tau }. \end{equation}
Here,
$\tilde {\boldsymbol{g}} = (g \sin {\theta }, 0, -g \cos {\theta })$
, with constant gravitational acceleration
$g$
, the pressure is
$\hat {p}$
and the deviatoric stress is
$\boldsymbol{{\hat {\tau }}}=\{{\hat {\tau }}_{\textit{ij}}\}$
, which we relate to the rate of strains by the generalized Newtonian fluid model,
\begin{equation} {\hat {\tau }}_{\textit{ij}} = {\hat {\mu }}(\lambda )\, \hat {\dot {\gamma }}_{\textit{ij}} , \end{equation}
where
\begin{equation} \boldsymbol{\hat {\dot \gamma }} = \big\{ \hat {\dot {\gamma }}_{\textit{ij}}\big\} = \boldsymbol{\nabla }\boldsymbol{\hat u} + (\boldsymbol{\nabla }\boldsymbol{\hat u})^T . \end{equation}
The tensors have second invariants,
$\hat {\dot {\gamma }} = \sqrt {(1/2)\sum_{i,j}\hat {\dot \gamma }_{\textit{ij}}\hat {\dot \gamma }_{\textit{ij}}}$
and
${\hat {\tau }} = \sqrt {(1/2)\sum_{i,j}{\hat {\tau }}_{\textit{ij}}{\hat {\tau }}_{\textit{ij}}}$
. The viscosity
${\hat {\mu }}(\lambda )$
, is controlled by the parameter
$\lambda ({\hat {x}},{\hat {z}},{\hat {t}}\kern1pt)$
, which gauges the degree of microstructural order.
There is no slip at the base of the fluid and we assume that the top surface is stress-free,
\begin{equation} \begin{aligned} \boldsymbol{{\hat {u}}} &= 0 \quad \text{at} \quad {\hat {z}} = 0, \\ \left ({\hat {\tau }}_{\textit{ij}} - {\hat {p}} \delta _{\textit{ij}}\right )n_j &= 0 \quad \text{at} \quad {\hat {z}} = \hat {h}({\hat {x}},{\hat {t}}\kern1pt), \end{aligned} \end{equation}
where
$\boldsymbol{n}$
is the normal to the upper surface
${\hat {z}} = \hat {h}$
, and the kinematic condition there is
\begin{equation} \frac {\partial \hat {h}}{\partial {\hat {t}}} + {\hat {u}} \frac {\partial \hat {h}}{\partial {\hat {x}}} = {\hat {w}} \quad \text{at} \quad {\hat {z}} = \hat {h}({\hat {x}},{\hat {t}}\kern1pt). \end{equation}
The omission of surface tension is potentially problematic in thin films with relatively small length scales, as in the Kapitza problem (e.g. Chang Reference Chang1994). We assume here that the film has larger length scale so that capillary effects are small, our purpose being to interrogate instabilities driven by rheological effects rather than how surface tension might stabilize them.
2.2. Dimensionless leading-order formulation
To remove the dimensions (and hats) from the equations, we introduce the rescalings
\begin{align}&\qquad\qquad\qquad\quad {\hat {t}} = \frac {{\mathcal L}}{{\mathcal U}}t,\quad {\hat {x}} = {\mathcal L} x, \quad ({\hat {z}},\hat {h}) = {\mathcal H}(z,h), \end{align}
\begin{align}& {\hat {u}} = {\mathcal U} u, \quad {\hat {w}}=\varepsilon {\mathcal U} w,\quad {\hat {p}} = \rho g {\mathcal H} p \cos {\theta }, \quad {\hat {\tau }}_{\textit{ij}} = \frac {\mu _{_0} {\mathcal U}}{{\mathcal H}}\tau _{\textit{ij}}, \quad {\hat {\mu }} = \mu _{_0} \mu , \end{align}
where
$\mu _{_0}$
is a characteristic viscosity and the speed scale
\begin{align} {\mathcal U} = \frac {{\mathcal H}^3 \rho g}{{\mathcal L} \mu _{_0}} \cos {\theta } . \end{align}
With these scalings, and omitting terms of
$O(\varepsilon ^2)$
, (2.3) reduces to the lubrication equations
\begin{equation} 0 = 1 - \frac {\partial p}{\partial x} + \frac {\partial \tau _{\textit{xz}}}{\partial z} ,\quad 0 = -1 -\frac {\partial p}{\partial z}. \end{equation}
The dominant component of the rate of strain tensor is the vertical shear
$\dot {\gamma }_{\textit{xz}} = \partial u/\partial z+O(\varepsilon ^2)$
, and to leading order, the stress conditions in (2.6b
) become
\begin{equation} p=\tau _{\textit{xz}} = 0 \qquad \text{at} \qquad z=h(x,t). \end{equation}
The kinematic condition (2.7) is unchanged after scaling (but for the removal of the hats). Equations (2.11) and (2.12) imply that the pressure is hydrostatic,
\begin{equation} p = h - z, \end{equation}
and the shear stress is given by
\begin{align} \tau _{\textit{xz}} &= (h-z)\left ( 1-\frac {\partial h}{\partial x}\right ) = (1-\zeta ) \tau _b, \quad \tau _b = h \left (1-\frac {\partial h}{\partial x}\right )\!, \quad \zeta = \frac {z}{h}, \nonumber\\ &= \mu (\lambda )\frac {\partial u}{\partial z}, \end{align}
where
$\tau _b$
is the basal shear stress.
Finally, mass conservation integrated over the fluid layer implies
\begin{equation} \frac {\partial h}{\partial t} = - \frac {\partial }{\partial x}\int _0^h u \:\mbox{d}{z} = - \frac {\partial }{\partial x}\int _0^h\left (h-z\right )\frac {\partial u}{\partial z}\:\mbox{d}{z} . \end{equation}
Since
\begin{equation} \frac {\partial u}{\partial z} = \frac {\tau _{\textit{xz}}}{\mu (\lambda )} = \frac {(1-\zeta )\tau _b}{\mu (\lambda )}, \end{equation}
the constitutive law can now be plugged into (2.15) to arrive an evolution equation for the fluid thickness. The exercise is complicated by the structure parameter
$\lambda$
, however, which retains an unknown vertical profile, as discussed next.
2.3. Constitutive law
Thixotropic rheology has been described using simple constitutive models (Barnes Reference Barnes1997; Mewis & Wagner Reference Mewis and Wagner2009; Larson & Wei Reference Larson and Wei2019) in which the structure parameter
$\lambda (x,z,t)$
satisfies the (dimensional) evolution equation,
\begin{equation} \frac {\partial \lambda }{\partial {\hat {t}}}+ {\hat {u}}\frac {\partial \lambda }{\partial {\hat {x}}}+ {\hat {w}}\frac {\partial \lambda }{\partial {\hat {z}}} = \frac {1-\lambda }{T} - \alpha \lambda \tilde {\dot {\gamma }} + K {\nabla} ^2\lambda , \end{equation}
and controls the fluid rheology through the viscosity law (2.4). The order parameter lies in a range
$[0,1]$
; for
$\lambda =0$
, the fluid has no effective microstructure, while it is fully structured when
$\lambda =1$
. The evolution equation in (2.17) contains a restructuring term
$(1-\lambda )/T$
that drives
$\lambda$
towards the fully structured state, and a destruction term dictated by the local shear rate
$\tilde {\dot {\gamma }}$
. Restructuring is characterized by a healing time scale
$T$
, and destruction is parameterized by
$\alpha$
. In some other related models, the rate of destructuring is taken to be proportional to a power of
$\tilde {\dot {\gamma }}$
, or a power of the stress invariant
$\tilde \tau$
; we use (2.17) in view of its simplicity. We have further included a diffusive term in (2.17) to account for any spatial structural diffusion, with diffusivity
$K$
, although we will be mostly concerned with the limit in which this physical effect is relatively weak.
For the viscosity law, we follow Hewitt & Balmforth (Reference Hewitt and Balmforth2013) and choose
\begin{equation} \tilde \mu (\lambda ) = \frac {\mu _{_O}}{(1-\lambda )[(1-\lambda )(1-a)+a]} =\mu _{_O}\mu (\varLambda ), \end{equation}
where
$\mu _{_O}$
is the characteristic viscosity and
$a$
is a parameter with
$0\lt a\lt 1$
. The resulting constitutive law is qualitatively similar to others proposed in the literature, in that it produces comparable flow curves in steady, uniform shear (Coussot et al. Reference Coussot, Nguyen, Huynh and Bonn2002; Dullaert & Mewis Reference Dullaert and Mewis2006; Møller et al. Reference Møller, Mewis and Bonn2006; Alexandrou et al. Reference Alexandrou, Constantinou and Georgiou2009). The precise form, however, proves analytically convenient, as discussed below (and is written slightly differently to that in Hewitt & Balmforth (Reference Hewitt and Balmforth2013)).
In the shallow limit, following the scalings of § 2.2, the dimensionless version of the constitutive model is
\begin{align} {\mathcal T} \left ( \lambda _t + u \lambda _x + w \lambda _z \right ) = 1 - \lambda - \varGamma {\dot {\gamma }} \lambda + \kappa \lambda _{\textit{zz}},\nonumber\end{align}
\begin{align}u_z = \tau _{\textit{xz}} A(\lambda ) ,\qquad A(\lambda ) \equiv [\mu (\lambda )]^{-1} = (1-\lambda )[(1-\lambda )(1-a)+a], \end{align}
where have used
$(t,x,z)-$
subscripts as shorthand for partial derivatives (except for the shear stress
$\tau _{\textit{xz}}$
), and the dimensionless groups,
\begin{equation} {\mathcal T} = \frac {{\mathcal U} T}{{\mathcal L}}, \quad \varGamma = \frac {\alpha T{\mathcal U}}{{\mathcal H}} , \quad \kappa = \frac {\textit{KT}}{{\mathcal H}^2}, \end{equation}
represent a characteristic relaxation time, healing rate and diffusivity. The definition of the reciprocal viscosity function
$A(\lambda )$
proves useful in the further analytical developments of the appendices, where its usage also highlights how the results are insensitive to the detailed analytical choice of the viscosity law.
The addition of diffusion in (2.19) demands the inclusion of further boundary conditions on
$\lambda (x,z,t)$
at the base and surface of the film. Assuming that structure is neither created nor destroyed by interfacial interaction, we adopt the no-flux conditions,
\begin{equation} \lambda _z(x,0,t)=\lambda _z(x,h,t)=0 . \end{equation}
2.4. Flow curves
Viscosity bifurcations are commonly exposed and illustrated using the flow curve of a model fluid; i.e. by considering how the shear stress depends on shear rate under conditions of steady, spatially uniform shear. For the current rheological model, the relations (2.19) imply the flow curve,
\begin{equation} \lambda = \frac {1}{1+\varGamma {\dot {\gamma }}} \quad {\rm and} \quad \tau = |\tau _{\textit{xz}}| = \frac {(1+\varGamma {\dot {\gamma }})^2}{\varGamma (a+\varGamma {\dot {\gamma }})} . \end{equation}
Consequently,
$\varGamma$
dictates the dimensionless shear rate at which destructuring becomes effective, whereas
$a$
controls the shape of the flow curve, as illustrated in figure 2. For
$\varGamma {\dot {\gamma }}\gg 1$
, we recover the Newtonian law,
$\tau ={\dot {\gamma }}$
. Conversely, for
$\varGamma {\dot {\gamma }}\ll 1$
, we find
\begin{equation} \tau = \frac {1}{\varGamma a} \left [ 1 + \frac {2a-1}{a}\varGamma {\dot {\gamma }}+\cdots \right]\!. \end{equation}
Unless
$a=0$
, the flow curve therefore intersects the
$\tau \hbox{-}$
axis at a finite stress,
$\tau\!_{_A} = (\varGamma a)^{-1}$
. If
$a\gt (1/2)$
, the flow curve increases monotonically away from this yield point towards higher shear rates, in the manner of a simple (non-thixotropic) yield-stress fluid, see the lower flow curves in figure 2(a). Beginning from fully structured fluid (
$\lambda =1$
and
${\dot {\gamma }}=0$
), the threshold stress above which motion must take place is then
$\tau \gt \tau\!_{_A}$
.
Steady-state flow curves for (a) varying
$a$
(scaling
$\dot {\gamma }$
and
$\tau$
by
$\varGamma$
), and (b)
$(\varGamma ,a)=(8,{1/5})$
. The special values,
$(0,\tau\!_{_A})$
and
$({\dot {\gamma }}_{_C},\tau _{_C})$
, are indicated. The dotted line in (a) shows the locus of
$({\dot {\gamma }}_{_C},\tau _{_C})$
for varying
$a$
; the specific flow curves shown have
$a=0$
,
$0.1$
,
$0.225$
,
$0.4$
,
$0.666$
and 1. The (red) arrows in (b) indicate the path of a hysteretic loop taken on first increasing then decreasing the stress, starting from an initially structured state with
$\varLambda =1$
.
For
$a\lt (1/2)$
, on the other hand, the flow curve becomes non-monotonic, first descending to lower stresses for a range of shear rates. The curve subsequently reaches a minimum at
$({\dot {\gamma }},\tau )=({\dot {\gamma }}_{_C},\tau _{_C})$
, and then ascends to higher stresses. The minimum stress is given by
\begin{align} \tau _{_C} = \frac {4(1-a)}{\varGamma } \qquad \left ( {\dot {\gamma }}_{_C} = \frac {1-2a}{\varGamma }, \quad \lambda _{_C} = \frac {1}{2(1-a)} \right)\!, \end{align}
and implies an alternative stress threshold for motion of
$\tau \gt \tau _{_C}$
. The upper four flow curves in figure 2(a) illustrate this alternative behaviour. As illustrated in figure 2(b), such flow curves imply excursions along hysteretic loops on first increasing, then decreasing the stress, in the conventional manner of thixotropy; sudden jumps in shear rate, or viscosity bifurcations, arise at the stresses
$\tau =\tau\!_{_A}$
and
$\tau =\tau _{_C}$
. The multiplicity of possible states for
$\tau\!_{_A}\lt \tau \lt \tau _{_C}$
, also implies a sensitivity to initial preparation or flow history. The limiting case
$a=0$
, implying
$\tau\!_{_A}\to \infty$
and no finite yield stress for the fully structured state, is similar to the models proposed by Coussot et al. (Reference Coussot, Nguyen, Huynh and Bonn2002) and Møller et al. (Reference Møller, Mewis and Bonn2006).
3. Base states
Steady sheet-like flows are given by
$\tau =1-z$
,
$u=U(z)$
and
$\lambda =\varLambda (z)$
, where
\begin{equation} \begin{aligned} \kappa \varLambda _{\textit{zz}} &= \varGamma \varLambda U_z - (1 - \varLambda ), \cr U_z &= (1-z) A(\varLambda ), \end{aligned} \quad U(0)=\varLambda '(0) = \varLambda '(1) = 0. \end{equation}
Sample solutions to this system are shown in figure 3 for two representative parameter settings, the first case without a viscosity bifurcation, the other case including one.
Base states for (a,b,c)
$a=(3/5)$
and (d,e,f)
$a={1/5}$
, with
$\varGamma =8$
and
$\kappa =10^{-j}$
,
$j=3,4,\ldots ,8$
. In (e–f), the green dot–dash line indicates
$z=Y$
and
$\tau =1-Y$
, with
$Y$
given by (3.6). The (red) dashed lines indicate the diffusionless solution (computed assuming the same
$Y$
in (d–f)); the (red) dotted line in (f) shows the untraced part of the flow curve.
In the limit that
$\kappa \to 0$
, (2.19) or (3.1) imply that the solutions follow stress-controlled sections of the flow curve. Consequently, to achieve flow, the stress must exceed one of the flow thresholds,
$\tau\!_{_A}$
and
$\tau _{_C}$
, somewhere within the layer. Given that
$\tau =1-z$
, the maximum stress is unity and occurs at the base. Therefore, flow requires
$\tau\!_{_A}\lt 1$
or
$\tau _{_C}\lt 1$
. For a monotonic flow curve with
$a\gt (1/2)$
, the relevant criterion is
$\tau\!_{_A}^{-1}=\varGamma a\gt 1$
. Conversely, if the flow curve first bends downwards and
$a\lt (1/2)$
, then the base of the fluid can be destructured and flowing if
$1\gt \tau _{_C}$
. Nevertheless, if
$\tau\!_{_A}\gt 1\gt \tau _{_C}$
, it is also possible that the entire fluid layer is fully structured and plugged up; only when
$\tau\!_{_A}\lt 1$
are we assured that flow must take place. For
$\tau\!_{_A}\gt 1\gt \tau _{_C}$
, we have the potential for two possible states, one flowing and one fully plugged.
Diffusionless, base-flow solutions therefore trace out a locus on the stress-shear-rate plane which begins at the bottom of the
$\tau -$
axis (i.e. the top of the fluid layer). For
$a\gt (1/2)$
and
$1\gt \tau\!_{_A}$
, the locus moves up that axis until it reaches the yield point,
$\tau =\tau\!_{_A}$
, then shifts to finite shear rates along the monotonic flow curve (all the while descending through the fluid layer); the yield point arises at the height
$z=z_{_A}=1-\tau\!_{_A}$
. In other words, the fluid is fully structured above
$z=z_{_A}$
, and is rafted along as a rigid plug by an underlying layer of partially structured material. This leads to profiles for the velocity and shear rate that are continuous and similar to those for a Bingham fluid (Balmforth et al. Reference Balmforth, Craster, Rust and Sassi2006). Figure 3(a–c) illustrate such profiles for
$(\varGamma ,a)=(10,(3/5))$
. Solutions to (3.1) for
$\kappa \gt 0$
are also included in this figure, and demonstrate how the switch at
$z\gt z_{_A}$
becomes smoothed by small diffusion, but otherwise the profiles remain similar. Notably, the upper section of the fluid layer is not fully structured (
$\varLambda \lt 1$
), but is, nevertheless, plug-like with
$\varLambda \approx 1$
and
$U_z\approx 0$
.
The other example in figure 3(e–f) displays a case with a non-monotonic flow curve and a viscosity bifurcation. Here, the locus of the base state on the
$({\dot {\gamma }},\tau )\hbox{-}$
plane first moves up the
$\tau \hbox{-}$
axis, but then switches abruptly over to the right-hand part of the flow curve at a stress level
$\tau =1-Y$
that lies between the two thresholds,
$\tau\!_{_A}$
and
$\tau _{_C}$
. Without diffusion, the jump is discontinuous; for
$\kappa \gt 0$
, the transition becomes smoothed. The abrupt jump leads to a sharper transition from the superficial plug to much less structured material underneath, with repercussions on the velocity and shear-rate profiles.
In principle, any level
$z=Y$
for which
$\tau\!_{_A}\gt 1-Y \gt \tau _{_C}$
could locate a jump from structured to destructured fluid in the diffusionless problem; the level
$Y$
must be selected as an initial condition (cf. Wilson & Fielding Reference Wilson and Fielding2006). The choice shown in figure 3(e–f) by the dashed (red) line corresponds to the particular location that is selected for
$\kappa \to 0$
in the problem with diffusion (solid blue lines). To find this location, we consider (3.1) in more detail: for
$\kappa \ll 1$
, the transition arises over a narrow region with a thickness of
$O(\sqrt \kappa )$
. Setting
$z=Y+\sqrt \kappa \zeta$
, (3.1) then reduces to
\begin{equation} \kappa \varLambda _{\textit{zz}} \equiv \varLambda _{\zeta \zeta } \sim \varPsi (\varLambda ;Y,\varGamma ,\textrm {Bi}) = (1-\varLambda ) \{ \varGamma \varLambda [(1-\varLambda ) (1-a)+a] (1-Y) - 1 \} . \end{equation}
The level
$z=Y$
now follows by imposing the matching conditions,
\begin{equation} \varLambda \to 1 \ \textrm {for} \ \zeta \to \infty \quad \textrm {and} \quad \varLambda \to \varLambda _a \ \textrm {for} \ \zeta \to -\infty , \end{equation}
where
\begin{equation} \varLambda _a = \frac {1}{2(1-a)} \left [ 1 - \sqrt {1 - \frac {4(1-a)}{\varGamma (1-z)}} \right ] \end{equation}
corresponds to the structure function along the right-hand branch of the flow curve. The exercise is equivalent to demanding that
\begin{equation} \int _{\varLambda _a}^{\varLambda _*} \varPsi (\varLambda ;Y,\varGamma ,\textrm {Bi}) \; \textrm {d}\varLambda = 0 , \end{equation}
which implies that
\begin{equation} Y = 1-\frac {9(1-a)}{\varGamma (2-a)(1+a)} , \end{equation}
and successfully locates the transition level in figure 3(e–f).
4. Stability theory
4.1. Linear equations
Infinitesimal perturbations to the base states, with
\begin{equation} h = 1+{\check \eta } e^{{\textit{ikx}}+\sigma t}, \quad u = U(z)+\psi _z(z) e^{{\textit{ikx}}+\sigma t}, \quad \lambda = \varLambda (z)+{\check \lambda }(z) e^{{\textit{ikx}}+\sigma t}, \end{equation}
wavenumber
$k$
and (complex) growth rate
$\sigma =\sigma _r+i\sigma _i$
, satisfy
\begin{align} \psi _{\textit{zz}} = \frac {U_z{\check \tau }}{1-z} + (1-z) A'(\varLambda ) {\check \lambda },\nonumber\end{align}
\begin{align}{\mathcal T} [ (\sigma +ikU) {\check \lambda } - ik \varLambda _z\psi ] = - \left ( 1 + \varGamma U_z \right ) {\check \lambda } - \varGamma \varLambda \psi _{\textit{zz}} + \kappa {\check \lambda }_{\textit{zz}} .\end{align}
Here,
$\psi (z)$
represents a stream function such that
$(u,w)=(U+\psi _ze^{{\textit{ikx}}+\sigma t}, -ik\psi e^{{\textit{ikx}}+\sigma t})$
, and the shear stress perturbation has amplitude
\begin{equation} {\check \tau } = [1 - ik(1-z)]{\check \eta } . \end{equation}
The boundary conditions translate to
\begin{equation} [\sigma + ik U(1)] {\check \eta } = -ik \psi (1), \qquad \psi (0)=\psi _z(0)={\check \lambda }_z(0)={\check \lambda }_z(1)+{\check \eta }\varLambda _{\textit{zz}}(1)=0 \end{equation}
(after a shift of the moving surface position back to
$z=1$
in the last one).
The problem in (4.2)–(4.3) can be solved by discretizing in
$z$
, and then using Chebyshev differentiation matrices to turn the eigenvalue problem into a matrix one. Though limited in spatial resolution, this technique expedites the search for unstable modes. Once found, MATLAB’s inbuilt solver bvp4c can be used to compute those modes with better resolution, and then continue them to different parameter settings. Alternatively, we may take different limits of the problem and reduce the equations further to obtain some more explicit results. That exercise is accomplished in Appendices B–E.
4.2. Results
(a) Growth rates
$\sigma _r$
and (b) scaled phase speeds
$c/U_z(0)=-\sigma _i/kU_z(0)$
, for varying
$k$
, with
$a=({1}/{20})(0,1,2,\ldots ,6)$
(from red to blue) and
$(\varGamma ,\kappa ,{\mathcal T})=(8,10^{-4},1)$
. The most unstable modes for the three lowest values of
$a$
are indicated by stars. The inset compares the growth rates, scaled
$k^2$
, with the predictions of the long-wave analysis of Appendix B (dashed lines). The eigenfunctions,
$\psi _z/{\check \eta }$
and
${\check \lambda }/{\check \eta }$
, of the most unstable mode for
$a=0$
are plotted in (c) and (d). The level
$z=Y$
from (3.6) is shown by the light grey line. Also plotted in (d) is
$\varLambda _z$
for the base state, scaled to have the same maximum amplitude as
$|{\check \lambda }/{\check \eta }|$
(dots).
Figure 4(a,b) displays sample numerical solutions to the linear stability problem, plotting growth rates and phase speeds against
$k$
for various values of
$a$
, with
$(\varGamma ,\kappa ,{\mathcal T})=(5,10^{-4},1)$
. For the higher values of
$a$
, the flow curve is monotonically increasing, and the flow is stable. When
$a$
is decreased sufficiently into the regime in which a viscosity bifurcation occurs, however, a band of unstable wavenumbers appears, extending from
$k=0$
up to a critical value
$k_\textit{crit}$
. The phase speeds of the unstable modes are close to
$U_z(0)$
, a feature that can be uncovered from a long-wave analysis (Appendix B).
The eigenfunction of the most unstable mode over all wavenumber
$k$
for the case with
$a=0$
is plotted in figure 4(c,d). For this mode, an increase in fluid depth enhances downslope flow speed (
$\psi _z$
; figure 4
c), whilst strongly perturbing the structure function over the transition region buffering the destructured lower layer from the overlying plug-like flow (figure 4
d). Indeed,
$\check \lambda$
is opposite in phase to
$\check \eta$
and its profile is approximately dictated by
$\varLambda _z(z)$
, implying that structural perturbations are mostly due to a vertical shift of the transition region. In other words, the instability operates by destructuring fluid underneath the wave crests, which accelerates flow, enabling the wave to entrain more fluid and thereby amplify. A more detailed interrogation of the instability mechanism follows from the long-wave analysis (Appendix B): to leading order in
$k$
, perturbations take the form of waves propagating with phase speed
$c\sim U_z(0)$
. These waves perturb the flow over the destructured lower layer, where quasisteady viscous perturbations turn out to be stabilizing, but the delay introduced by relaxation contributes to instability. More significant are the associated destabilizing vertical shifts of the transition region. All these effects appear at
$O(k^2)$
in the long-wave expansion, and point to an interfacial-type instability as suggested by Hewitt & Balmforth (Reference Hewitt and Balmforth2013).
Growth rate as a density over the
$(k,{\mathcal T})$
plane for the values of
$a$
indicated, with
$(\varGamma ,\kappa )=(8,10^{-4})$
. The red contours show the stability boundary, and the dots indicate the most unstable mode over all wavenumber. In (e), the triangle shows the scaling
${\mathcal T}\propto k^{-2}$
.
Growth rate as a density over the
$(k,a)$
plane for the values of
$\mathcal T$
indicated, all with
$(\varGamma ,\kappa )=(8,10^{-4})$
. The red contours show the stability boundary, and the dots indicate the most unstable mode over all wavenumber. The growth rates are scaled by their maxima; the colour scale can be inferred from (f), which plots the maximum growth rate (achieved at the smallest values of
$a$
) against
$\mathcal T$
. The (green) stars show the specific values of
$\mathcal T$
used in panels (a)–(e), whereas the (red) hexagrams show the additional cases that are also presented in figure 8(a). The triangle shows the scaling
$\sigma _r\propto {\mathcal T}^{-1}$
.
Growth rate as a density over the
$(k,\kappa )$
plane for the values of
$\mathcal T$
indicated, all with
$(a,\varGamma )=({1/5},8)$
. The red contours show the stability boundary, and the dots indicate the most unstable mode over all wavenumber. The triangle in (b) shows the scaling
$\kappa \propto k^4$
. The growth rates are scaled by their maxima; the colour scale can be inferred from (f), which plots the maximum growth rate against
$\mathcal T$
. The stars show the specific values of
$\mathcal T$
used in panels (a)–(e), whereas the filled circles show additional cases that are also presented in figure 8(c). The solid line shows the prediction in the limit of large diffusion
$\kappa \gg 1$
(Appendix E).
Stability boundaries on the (a)
$(k\sqrt {\mathcal T},a)$
, (b)
$(k,{\mathcal T})$
and (c)
$(k\sqrt {\mathcal T},\kappa )$
planes. The panels correspond to figures 5–7, with the stability boundaries shown from red to blue for increasing
$\mathcal T$
,
$a$
and
$\mathcal T$
, respectively. The dashed lines show additional cases, corresponding to
$a=0.1$
for panel (b), or the extra points plotted in figures 6(f) and 7(f) (for panels (b) and (c)). The stars along the left-hand axes show the predictions of the long-wave, small
$\kappa$
analysis of Appendix B. In (c), the triangles for
$\kappa =10^2$
indicate predictions using the large diffusion theory of Appendix E; the triangles along
$\kappa =10^{-8}$
indicate predictions in the zero-diffusion limit (Appendix D).
A wider view of the unstable regime is presented in figures 5–7, which plot growth rates as densities over the
$(k,{\mathcal T})$
,
$(k,a)$
and
$(k,\kappa )$
planes (respectively), each time holding the other parameters fixed. All three figures contain five cases, for varying
$a$
in figure 5, and varying
$\mathcal T$
in figures 6 and 7. In every example, the instability has a long-wave character, with
$k_\textit{crit}=O(10^{-1})$
or less. The stability boundaries extracted from the figures are assembled in figure 8, plotting
$k_\textit{crit}$
against
$\mathcal T$
,
$a$
or
$\kappa$
(including some additional cases to better highlight the effect of varying the other parameters).
A first conclusion that can be drawn from figure 5 is that instability is eventually suppressed for sufficiently small relaxation time
$\mathcal T$
when diffusion is not strong, regardless of the value of
$a$
. Indeed, in the limit
${\mathcal T}\to 0$
and
$\kappa \ll 1$
, one can further reduce the thin-film model to a nonlinear diffusion equation, as summarized in Appendix C. The base states in this limit can then be shown to be stable (at least in the absence of inertia).
More curiously, figures 5 and 6 also demonstrate that instability extends into the regime without a viscosity bifurcation (
$a\gt (1/2)$
), provided the relaxation time
$\mathcal T$
is long enough. This result is also apparent in figures 6 and 8. The rescaling used in the latter figure emphasizes how stability eventually depends on the combination
$k\sqrt {\mathcal T}$
for
${\mathcal T}\gg 1$
. Evidently, when
$\mathcal T$
becomes sufficiently large, the destabilizing effect of viscous relaxation overwhelms stabilization by quasistatic viscous effects, even when there is no viscosity bifurcation and the base state possesses no sharp interface. Instabilities associated with the bulk of the fluid, rather than a sharp interface, have also been observed for other shear flows of complex fluids (Renardy & Renardy Reference Renardy and Renardy2017; Castillo & Wilson Reference Castillo and Wilson2020).
Finally, figure 7 illustrates the impact of varying the diffusion coefficient
$\kappa$
. Earlier figures display results for
$\kappa =10^{-4}$
, which is sufficiently small that diffusion does not have a dramatic effect on the base states (other than broadening the interface-like transition region). For low relaxation time (
${\mathcal T}\lesssim 4$
in figure 7) the impact of diffusion on the linear stability is more subtle: instability becomes restricted to a low-wavenumber band with
$k_\textit{crit} \propto \kappa ^{1/4}$
. Instability therefore becomes suppressed for
$\kappa \to 0$
. By contrast, with higher
$\mathcal T$
, growth rates converge to finite values. The limit of small diffusion is considered in Appendix D; corresponding predictions for the position of the stability boundary are included in figure 8(c).
Figure 7 also highlights how instability persists even with relatively strong diffusion. For
$\kappa \gg 1$
, the structure function becomes uniform across the film, and the evolution equations in (2.19) and (2.15) simplify further, as outlined in Appendix E. The prediction for the large-
$\kappa$
limit of the stability boundary are again shown in figure 8(c). Note that the uniform profile of the structure function implies that, in this limit, instability cannot be the result of interfacial dynamics, but solely due to bulk viscous relaxation.
5. Nonlinear dynamics
To explore the nonlinear dynamics of unstable waves, we solve the thin-film equations (2.15) and (2.19) numerically in periodic-in-
$x$
domains with length
$\ell =2\pi /k$
(some details of the numerical scheme are provided in Appendix F). A sample solution for a case in which a viscosity bifurcation occurs (
$a={1/5}$
) and vertical diffusion is relatively small (
$\kappa =10^{-4}$
) is shown in figure 9. The initial perturbation to the base flow seeds the growth of an unstable wave, which amplifies according to linear theory. The instability subsequently saturates at finite amplitude to generate a steadily propagating nonlinear wave. The final wave travels slightly faster than the original unstable linear mode, as evident from the snapshots in figure 9(c), which plots the wave height profile in the linear wave frame (with phase speed
$c$
). Also shown (figure 9
d), are snapshots of the interface between structured and destructured fluid (identified as the contour where
$\lambda (x,Y,T)=(4/5)$
). In line with the spatial structure of the linear mode, the interface elevates below the rise to the crest of the wave, as fluid locally destructures due to increased shear stresses.
Numerical solution to (2.14), (2.15) and (2.19), for
$(a,\varGamma ,{\mathcal T},\kappa )=({1/5},8,10^2,10^{-4})$
in a periodic domain of length
$20\pi$
. In (a), we present snapshots of the evolving profiles of
$h(x,t)$
and
$\lambda (x,z,t)$
. The times of the snapshots are indicated by the stars in (b), which plots the time series of
$\sqrt {\langle (h-1)^2\rangle }$
; the dashed line shows the expected linear growth. Also shown is another numerical solution with
$a=(3/5)$
. Snapshots of
$h$
and the level at which
$\lambda =(4/5)$
(dashed lines in (a)) are plotted in (c) and (d) for the first solution (with
$a={1/5}$
), in the frame of linearly unstable wave (which travels at phase speed
$c$
). The snapshots are spaced by 100 time units.
The final nonlinear wave is shown in more detail in figure 10, which presents the distribution of the structure function, together with a selection of streamlines (in the frame of the wave). Also shown (in figure 10
b) is the final wave for a numerical solution with
$a=(3/5)$
. Although this second example does not feature a viscosity bifurcation, and the division between structured and destructured fluid is smoother, the dynamics of the instability is similar, with linear growth again saturating to create the nonlinear wave illustrated (see also figure 9
b).
Final nonlinear wave solutions for (a)
$a={1/5}$
and (b)
$=(3/5)$
, with
$(\varGamma ,{\mathcal T},\kappa ,k)=(8,10^2,10^{-4},0.1)$
. Shown are density plots of
$\lambda (x,z,t)$
, with a superposed selection of streamlines (red lines, in the frame of the final wave, which has speed
$c_n$
).
(a) Wave amplitudes
$(h_{max}-h_\textit{min})$
, (b) linear growth rates
$\sigma _r$
and (c) phase speeds
$-\sigma _i/k$
, plotted against
$\mathcal T$
for a suite of computations with
$(\varGamma ,a,\kappa ,k)=(8,{1/5},10^{-4},0.1)$
. The wave profiles for the cases shown by (red) circles are plotted in (d) (shifted horizontally to align their maxima). The squares indicate the computation also shown in figures 9 and 10(a).
A wider suite of computations is shown in figure 11. This suite corresponds to traversing a vertical section at
$k=0.1$
through figure 5(b). Waves are unstable within a band of relaxation times
$\mathcal T$
over this section, and figure 11 shows the corresponding wave amplitudes for the final, saturated states. Those amplitudes decrease continuously to zero at the edges of the band, illustrating how the transition to instability is supercritical. The nonlinear waves travel at speeds that are a per cent or so higher than the linear phase speed
$c$
(also plotted in figure 11
c), which in turn exceeds the mean flow speed
$U_*=0.207$
by a factor of approximately three or more.
6. Discussion
This paper has focussed on the dynamics of instabilities in thin films of thixotropic fluid flowing down an inclined plane. A previously suggested interfacial instability (Hewitt & Balmforth Reference Hewitt and Balmforth2013) has been confirmed to arise, and its details have been interrogated. The instability has a long-wave character and is weak, in the sense that the time scale for growth is rather longer than that for propagation. However, instability arises over a wide choice of the settings of the problem parameters; only when the time scale for microstructural relaxation becomes relatively short is instability suppressed. The instability persists even when there are no viscosity bifurcations in the underlying rheological model, or when there is a relatively large vertical diffusion of the order parameter representing microstructure. Both results are surprising given that the discontinuous interface hypothesized by Hewitt & Balmforth (Reference Hewitt and Balmforth2013) to be key to instability is then completely absent. Instead, the phase lag introduced by the microstructural relaxation primarily drives instability (cf. Renardy & Renardy Reference Renardy and Renardy2017; Castillo & Wilson Reference Castillo and Wilson2020).
An exploration of the nonlinear dynamics (in spatially periodic domains) reveals that instability saturates into steadily propagating nonlinear waves that travel close to linear phase speeds and faster than the mean flow. The transition to instability is supercritical, with wave amplitudes reaching a fraction (20 % or so) of the fluid depth. Interactions and secondary instabilities of the nonlinear waves have not been explored.
One key ingredient of the current exploration is that it is based on a thin-film reduction of the original governing equations. In this reduction, the dominant resistance to fluid motion arises through shear stresses across the depth of the film; the extensional stresses in the film’s plane are omitted from the force balance equations. One awkward consequence is that the two-layer structure of the flow becomes suspicious: the upper layer is assumed to be strongly structured and plug-like, yet is simultaneously able to deform within the plane of the film. This is permitted in the model by letting the viscosity become sufficiently large over the upper layer to render the flow profile plug-like, but not so large that the extensional stresses impact force balance (cf. Balmforth et al. Reference Balmforth, Craster and Toniolo2003). This constraint translates to a thickness-dependent limit on the viscosity: in the dimensionless notation of § 2.2, the terms involving the extensional stresses appear at
$O(\varepsilon ^2)$
in the thin-film limit. Therefore, as long as the viscosity
$\mu$
remains much smaller than
$\varepsilon ^{-2}$
, there is no issue with neglecting those stresses.
Nevertheless, a more careful approach would be to solve the full, Orr–Sommerfeld-like stability problem, relaxing the thin-film approximation whilst simultaneously including a finite Reynolds number. This extension would permit one to identify precisely how the extensional stresses eventually rigidify the upper layer and arrest any deformation of the free surface, whilst gauging the impact of inertia.
Acknowledgements
N.J.B. would like to thank the Isaac Newton Institute for Mathematical Sciences, Cambridge, for support and hospitality during the programme `Anti-diffusive dynamics: from subcellular to astrophysical scales’, where work on this paper was undertaken.
Funding
This work was supported by EPSRC grant no EP/K032208/1.
Declaration of interests
The author reports no conflict of interest
Appendix A. Instantaneous, non-diffusive adjustment
Hewitt & Balmforth (Reference Hewitt and Balmforth2013) consider a (rather drastic) simplification of the thixotropic model of § 2.3 in which there is no diffusion (
$\kappa =0$
) and the fluid adjusts instantaneously to the steady-shear flow curve (
${\mathcal T}=0$
). Without diffusion, the two-layer structure of the base states, in which a superficial structured plug lies above a destructured shear layer (§ 3), continues to apply even when flow becomes unsteady. As long as the stress at the interface remains between the special values
$\tau\!_{_A}$
and
$\tau _{_C}$
, that division is imprinted by the initial state of the fluid and acts as a material surface. In addition, if the microstructure adjusts instantaneously to the flow dynamics (
${\mathcal T}=0$
), the upper layer remains plug-like and fully structured whereas the lower layer lies along the flow curve, with
\begin{equation} \lambda = \frac {1-\sqrt {1-4(1-a)(\varGamma \tau )^{-1}}}{2(1-a)}, \quad \tau =(1-h_x)(h-z). \end{equation}
The thin-film model now devolves to two mass conservation equations; one for the whole film and one for the destructured layer,
\begin{equation} \frac {\partial h}{\partial t} + \frac {\partial }{\partial x} F(h,h_x,Y) = 0, \qquad \frac {\partial Y}{\partial t} + \frac {\partial }{\partial x} G(h,h_x,Y) = 0, \end{equation}
where the flux functions,
\begin{equation} F(h,h_x,Y) = \int _0^h{u\:\mbox{d}{z}} = \left (1 - h_x\right )\int _0^Y \left (h-z\right )^2 \left (1-\lambda \right ) \left [(1 - \lambda )(1-a) + a\right ] \mbox{d}{z} , \end{equation}
and
\begin{equation} G(h,h_x,Y) = \int _0^Y{u\:\mbox{d}{z}} =\left (1 - h_x\right )\int _0^Y \left (h-z\right )\left (Y-z\right ) \left [(1 - \lambda )(1-a) + a\right ] \mbox{d}{z} , \end{equation}
with
$\lambda$
given by (A1), can be evaluated analytically although the explicit formulae are convoluted and not very informative.
Equations (A2) have the uniform equilibrium solution
$h=1$
and
$Y=Y_0$
, with
$0 \lt Y_0 \lt 1-4\varGamma ^{-1}(1-a)$
(ensuring
$\tau _{_C}\lt \tau (z=Y)\lt \tau\!_{_A}$
). Normal-mode perturbations to this base state of the form
$(h,Y)=(1,Y_0)+(\check {h},\check {Y})e^{\sigma t + ikx}$
, with wavenumber
$k$
and growth rate
$\sigma$
, satisfy a quadratic dispersion relation. It follows from this relation that the uniform flow is either stable, or one of the modes is unstable over all wavenumbers, with a growth rate that increases monotonically with
$k$
and converges to a constant for
$k\to \infty$
(see Hewitt & Balmforth Reference Hewitt and Balmforth2013).
Typical solutions of the dispersion relation are shown in figure 12. For the parameter settings chosen,
$(a,\varGamma )=((1/5),5)$
, all of the possible choices for the interface position
$Y_0$
lead to instability. The choice
$Y\approx 0.193$
is most unstable; the base flow selected by diffusion, with
$Y$
given by (3.6), is slightly less unstable. As highlighted by the amplitude ratios
$\check {Y}/\check {h}$
also plotted, the unstable mode perturbs
$Y$
more strongly that
$h$
. In other words, the mode has an interfacial character. By contrast, the damped mode has a relatively strong surface signature.
(a) Growth rate
$\sigma _r=\textrm{Re} (\sigma )$
, (b) phase speed
$c$
and (c) amplitude ratio
$\check {Y}/\check {h}$
for
$Y_0=0.193$
(blue) and (3.6) (red), with
$(a,\varGamma )=((1/5),5)$
. The unstable (stable) modes are plotted by the solid (dashed) lines. Panels (d)–(e) show corresponding plots for the unstable mode, with
$k=10^3$
and varying
$Y_0$
; the stars indicate the two cases in (a) and (b). The dash–dot line in (a) indicates how the growth rates are modified when linear diffusion terms are added to equations (A2), with a diffusivity of
$2\times 10^{-6}$
.
Nonlinear computations with (A2), including linear diffusion terms with equal diffusivities of
$2\times 10^{-6}$
. The periodic domain has length
${2\pi }/{5}$
;
$(\varGamma ,a)=(5,{1/5})$
and
$Y_0=0.193$
. Panels (a i) and (b i) show density plots of
$Y(x,t)$
; panels (a ii) and (b ii) are plots of the initial (dashed) and final (solid) profiles of
$h$
(blue) and
$Y$
(red); the lighter lines show intermediate snapshots (every 100 time units). In each plot, the solution is shown in a frame translating close to that of the final nonlinear waves. In (a), the initial perturbation to
$Y=Y_0$
is
$10^{-4}\sin 10x$
; for (b), that perturbation is a random superposition of the first ten Fourier modes, with mean amplitudes of
$10^{-4}$
. Panel (c) displays Max
$(Y-Y_0)$
against time (solid blue line for the solution in (a); dot–dash red line for that in (b)). The dashed line shows the growth expected for the linear
$k=10$
mode.
The system (A2) can be solved numerically with periodic boundary conditions to explore the nonlinear dynamics of the interfacial instability, starting from a small perturbation to the uniform base flow. To ease this task, we add linear diffusion terms to both equations in (A2) (with equal diffusivities of
$2\times 10^{-6}$
). These diffusion terms control the linear instability and ensure that the nonlinear solutions remain smooth. As illustrated figure 12(a), for
$(\varGamma ,a,Y_0)=(5,{1/5},0.193)$
, a broad, but finite, band of unstable modes then arises, with maximum growth for
$k\approx 10$
.
Figure 13 displays the results of such initial-value computations, starting with two different initial conditions, in domains of length
${2\pi }/{5}$
. In the first example, the initial condition is given by
$(h,Y)=(1,Y_0+10^{-4}\sin 10x)$
, where
$Y_0=0.193$
. For the second example, the initial perturbation to the interface
$z=Y_0$
corresponds to a random superposition of the first 10 Fourier modes, with mean amplitudes of
$10^{-4}$
.
For the case initialized with a
$k=10$
perturbation, the unstable mode develops as predicted by linear theory. In the nonlinear regime, the instability saturates to form two nonlinear, shock-like waves. The second solution displays a richer dynamics: the initial condition seeds the growth of five waves. On reaching finite amplitude, these waves interact, leading to a sequence of coarsening events, somewhat like that seen for roll waves in inertial viscous films (Balmforth & Mandre Reference Balmforth and Mandre2004). Eventually, a single wave remains.
Note that there is an awkward issue hidden in this reduced model: the upper layer is assumed to be fully structured, but it is also able to deform along the plane of the film, despite the implied infinite viscosity. This issue reflects an assumption inherent in the thin-film model, as discussed at the end of § 6 and Appendix D, and can be rationalized in the same manner as the pseudoplugs of thin films of yield-stress fluids (Walton & Bittleston Reference Walton and Bittleston1991; Balmforth & Craster Reference Balmforth and Craster1999).
Appendix B. Small
$\boldsymbol{k}$
After differentiating the base-state equations (3.1) in
$z$
, and some algebraic manipulations, the linear stability equations can be recast in the form
\begin{equation} \begin{aligned} (\psi +{\check \eta } U)_{\textit{zz}} =& - ik(1-z)A{\check \eta } + (1-z) A' ({\check \lambda } + {\check \eta } \varLambda _z) , \cr {\mathcal T} [ (\sigma +ikU) {\check \lambda } - ik \varLambda _z\psi ] =& - \left ( 1 + \varGamma U_z \right ) ({\check \lambda }+{\check \eta }\varLambda _z) - \varGamma \varLambda (\psi +{\check \eta } U)_{\textit{zz}} + \kappa ({\check \lambda } + {\check \eta }\varLambda _z)_{\textit{zz}} ,\nonumber\\[6pt] \end{aligned} \end{equation}
where
$A=A(\varLambda )$
and
$A'=A'(\varLambda )$
are evaluated for the base state. In the limit
$k\ll 1$
, we may then introduce the sequences
\begin{equation} \sigma = k \sigma _1 + k^2 \sigma _2 + \ldots , \quad \psi = \psi _0 + k \psi _1 + \ldots , \quad {\check \lambda } = {\check \lambda }_0 + k {\check \lambda }_1 + \ldots \end{equation}
to find a long-wave solution in which we select a normalization that fixes
$\check \eta$
. Here, we adopt
${\check \eta }=1$
.
When
$\kappa$
is small, the fact that
${\check \lambda }(z)$
and
$\varLambda _z(z)$
satisfy different conditions at
$z=0$
is not significant. The leading-order solution is then
\begin{equation} \sigma _1 = - i U_z(0), \quad {\check \lambda }_0 = - \varLambda _z, \quad \psi _0 = - [U-z U_z(0)]. \end{equation}
At the following order of
$k$
, we then have
\begin{equation} \begin{aligned} \psi _{1zz} =& -i(1-z)A + (1-z) A' {\check \lambda }_1 , \cr i {\mathcal T} (1-z) U_z(0) \varLambda _z =& - \left ( 1 + \varGamma U_z \right ) {\check \lambda }_1 - \varGamma \varLambda \psi _{1zz} + \kappa {\check \lambda }_{1zz} , \end{aligned} \end{equation}
with boundary conditions
\begin{equation} \sigma _2 = -i \psi _1(1), \qquad \psi _1(0)=\psi _{1z}(0)={\check \lambda }_{1z}(0)={\check \lambda }_{1z}(1) =0 . \end{equation}
The growth is therefore determined at
$O(k^2)$
, and contains two contributions when
$\kappa \ll 1$
and
$a\lt 1/2$
(parameter settings we focus on here). The first contribution is from the destructured layer in
$0\lt z\lt Y$
: here, we omit the diffusion term in (B4) and solve algebraically for
$\psi _{1zz}$
. Then,
\begin{equation} \sigma _2^{(D)} \sim -i \int _0^Y (1-z)\psi _{1zz} \textrm {d} z = - \int _0^Y \frac { {\mathcal T} U_z(0) (1-z) A' \varLambda \varLambda _z + A} {1 + \varGamma (1-z) A' \varLambda ^2}(1-z)^2 \textrm {d} z , \end{equation}
in which
$\varLambda \sim \varLambda _a$
from (3.4). Note that
$A'\lt 0$
but
$1 + \varGamma (1-z) A' \varLambda ^2 \gt 0$
. Hence,
$\sigma _2^{(D)}$
contains two contributions: the first is proportional to
$\mathcal T$
and positive; the second is negative. In other words, over the lower-lying destructured layer, quasisteady viscous effects are stabilizing, but relaxation is destabilizing.
The second contribution arises from adjustments over the thin transition region around
$z=Y$
. If we again introduce the inner coordinate
$\zeta =(z-Y)/\sqrt \kappa$
to resolve this region, then we observe that
\begin{equation} {\mathcal L} {\check \lambda }_1 \sim i \kappa ^{-\frac 12}{\mathcal T} (1-Y) U_z(0) \varLambda _\zeta + O(1) , \end{equation}
where
\begin{equation} {\mathcal L} = \frac {\partial ^2 }{\partial \zeta ^2} - 1 - \varGamma (1-z) (A + A'\varLambda ) . \end{equation}
But (3.1) implies that
\begin{equation} {\mathcal L} \varLambda _\zeta = - \sqrt \kappa \varGamma A \varLambda . \end{equation}
Hence, on multiplying (B7) by
$\varLambda _\zeta$
and integrating over the transition region, we find
\begin{equation} i \kappa ^{-\frac 12}{\mathcal T} (1-Y) U_z(0) \int _{-\infty }^\infty \varLambda _\zeta ^2 \textrm {d} \zeta + O(1) \sim -\sqrt \kappa \varGamma \int _{-\infty }^\infty A \varLambda {\check \lambda }_1 \textrm {d} \zeta , \end{equation}
given the matching conditions for
$\zeta \to \pm \infty$
(which indicate that
$\varLambda$
and
${\check \lambda }_1$
and their derivatives in
$z$
are all order one). The balance in (B10) demands that
${\check \lambda }_1=O(\kappa ^{-1})$
over the transition region. Given (B9) and (B7), the solution for
${\check \lambda }_1$
must therefore be some multiple of
$\varLambda _\zeta$
at leading order:
${\check \lambda }_1 \sim - i \alpha \varLambda _\zeta$
. Equation (B10) then implies that
\begin{equation} \alpha \sim \frac { {\mathcal T} (1-Y) U_z(0)\int _{-\infty }^\infty \varLambda _\zeta ^2 \textrm {d} \zeta } {\kappa \varGamma \int _{-\infty }^\infty A \varLambda \varLambda _\zeta \textrm {d} \zeta } . \end{equation}
Moreover, from the first relation in (B4),
\begin{equation} \psi _{1\zeta \zeta } \sim \alpha \kappa (1-Y) A' \varLambda _\zeta , \end{equation}
and so the transition region makes a contribution to the growth rate of
\begin{equation} \sigma _2^{(J)} \sim -i (1-Y) \int _{-\infty }^\infty \psi _{1\zeta \zeta } \frac {\textrm {d}\zeta }{\sqrt \kappa } \sim \sqrt \kappa \alpha (1-Y)^2 A(Y^-), \end{equation}
where
$A(Y^-)$
denotes the reciprocal viscosity function evaluated as
$z\to Y$
from below. Note that
$\sigma _2^{(J)}$
is positive and always contributes towards instability. Moreover,
${\check \lambda }_1 \sim -i \alpha \varLambda _\zeta$
represents a vertical shift of the transition layer due to the translation of wave-like perturbations with speed
$-\textrm{Im} (\sigma _1)=U_z(0)$
. In other words, the motion of the transition region is destabilizing.
Appendix C. Fast relaxation with small diffusion
When
${\mathcal T}\to 0$
, the structure function satisfies
\begin{equation} \kappa \lambda _{\textit{zz}} = (1-\lambda )\{\varGamma \tau (1-\varLambda )[(1-\varLambda )(1-a)+a] - 1\} . \end{equation}
In the further limit
$\kappa \ll 1$
, we then arrive at the structure function,
\begin{equation} \lambda = \left \{\!\! \begin{array}{ll} \varLambda _a, & 0\lt z\lt Y , \\ 1, & Y\lt z\lt h , \end{array} \right . \end{equation}
where
\begin{equation} \varLambda _a(\tau ) = \frac {1}{2(1-a)}\left [ 1 - \sqrt {1 - \frac {4(1-a)}{\varGamma \tau }} \right]\!, \quad Y = h - \frac {9(1-a)}{\varGamma (1-h_x)(2-a)(1+a)}, \end{equation}
which corresponds to a generalization of the solution to (3.2)–(3.6), but applying for the full nonlinear dynamics of the film. Armed with this profile, we may compute the flux
\begin{equation} \int _0^h (h-z) {\dot {\gamma }} dz = \frac {h^2}{\tau _b^2} \int _{{\tau _{_{Y}}}}^{\tau _b} {\hat \tau }^2 {\mathcal A}({\hat \tau }) \textrm {d}{\hat \tau } = h^2 G(\tau _b) , \end{equation}
where
$\tau _b=h(1-h_x)$
is the basal shear stress,
\begin{equation} {\tau _{_{Y}}}=(h-Y)(1-h_x)=\frac {9(1-a)}{\varGamma (2-a)(1+a)}, \end{equation}
and
${\mathcal A}(\tau ) = (1-\varLambda _a)[(1-\varLambda _a)(1+a)+a]$
is the reciprocal of the viscosity function, rewritten as a function of
$\tau$
using (C3). Introducing this flux into the evolution equation
$h_t + Q_x = 0$
, gives the nonlinear diffusion problem,
\begin{equation} h_t + \big[h^2 G(\tau _b) \big]_x = 0 . \end{equation}
Linear stability is dictated by the dispersion relation,
\begin{equation} \sigma = -ik [2G(1)+G'(1)] - k^2 G'(1). \end{equation}
But
\begin{equation} G'(1) \equiv {\tau _{_{Y}}}^2 \varGamma ({\tau _{_{Y}}}) + \int _{{\tau _{_{Y}}}}^1 \tau ^2 \varGamma '(\tau ) \; \textrm {d}\tau , \end{equation}
where
$\varGamma (\tau )\equiv \tau {\mathcal A}(\tau )$
denotes the inverse of the ascending (right-hand) branch of the flow curve in (2.22). Thus, as
${\tau _{_{Y}}}^2 \varGamma ({\tau _{_{Y}}}) \equiv \tau _{_C}^2{\dot {\gamma }}_{_C} \gt 0$
and
$\varGamma '(\tau )\gt 0$
, we conclude
$G'(1)\gt 0$
. The film is therefore stable in this limit. Note that the situation may be different if the base state were to contain sections of the descending (left-hand) branch of the flow curve.
Appendix D. Vanishing diffusion
When
$\kappa \to 0$
and there is a viscosity bifurcation (
$a\lt 1/2$
), the interface
$z=Y(x,t)$
evolves as a material surface. Mass conservation over the plug-like structured layer then demands that
\begin{equation} (h-Y)_t + [(h-Y) {u_{_P}}]_x = 0 , \end{equation}
where
${u_{_P}}\equiv u(x,Y,t)$
is the plug speed.
The base states arising in the
$\kappa \to 0$
limit of the diffusive problem have interface positions given by (3.6), and flow profiles following the flow curve. Base states are given by
$(u,\lambda ,h,Y)=(U,\varLambda ,1,Y_0)$
. With
$\kappa \to 0$
, the linearized structure function equation (4.2) simplifies, reducing to an algebraic relation for
$\check \lambda$
. One must then integrate only the stream function equation to solve the linear problem over
$0\lt z\lt Y_0$
. The overlying plug-like flow can be taken into account by applying a boundary condition at
$z=Y_0$
: a linearization of (D1) indicates
\begin{equation} \check {Y} = {\check \eta } + \frac {ik (1-Y_0) \check {u}_{_P}}{\sigma +ik{U_{_P}}} , \end{equation}
whereas the continuity of the velocity at
$z=Y$
indicates that
\begin{equation} \check {u}_{_P} = \psi _z(Y_0) + U_z(Y_0) \check {Y}. \end{equation}
Rewriting the kinematic condition at the free surface as
\begin{equation} (\sigma +ik{U_{_P}}) {\check \eta } = -ik\psi (1) = -ik[\psi (Y_0) + (1-Y_0) \check {u}_{_P}], \end{equation}
and then eliminating
$\check {Y}$
and
$\check {u}_{_P}$
from (D2)–(D4) furnishes the condition needed at
$z=Y_0$
.
Note that the plug-like layer is again able to deform in this limit, despite its apparently fully structured state. Once more, this reflects the omission of the extensional stresses from the shallow-film model. Therefore, even though
$\kappa \ll 1$
, diffusion is assumed sufficient to prevent
$\lambda$
from truly reaching unity. The viscosity function is thereby prevented from diverging, and simply becomes sufficiently large to suppress vertical shear, but not horizontal extension. In other words, the limit implicitly assumes that
$\varepsilon ^2\mu (\varLambda )\ll 1$
(given that the extensional stresses impact the force balance equations at
$O(\varepsilon ^2)$
).
Appendix E. Strong diffusion
In the limit
$\kappa \gg 1$
, the structure function becomes locally uniform across the film:
$\lambda \sim \overline {\lambda }(x,t)$
. The velocity profile is then parabolic and a reduction of the thin-film model follows after integrating the dimensionless structure–function equation in (2.19) over the film thickness,
\begin{align} h_t + (hU)_x &= 0, \end{align}
\begin{align} {\mathcal T} \left (\overline {\lambda }_t + U \overline {\lambda }_x \right ) &= 1 - \overline {\lambda } - \frac {3\varGamma U \overline {\lambda }}{2h}, \end{align}
\begin{align} U \equiv \frac 1h \int _0^h (h-z)u_z\;\textrm {d} z &= \frac {1}{3} h^2 (1-h_x) A(\overline {\lambda }) . \end{align}
Base states exist for
$\varGamma \gt 8(1-a)$
and are given by
\begin{equation} h = 1, \quad \overline {\lambda } = \varLambda = \frac { 1 - \sqrt {1-8\varGamma ^{-1}(1-a)} }{2(1-a)} , \quad U = U_* = \frac {1}{3} A(\varLambda ) , \end{equation}
corresponding to a point on the right-hand branch of the flow curve, with shear rate
$U_*$
and unit shear stress. Linear perturbations with growth rate
$\sigma$
and wavenumber
$k$
have the dispersion relation,
\begin{equation} {\Big[{\mathcal T}(\sigma +ikU_{*}) + 1 + \frac {3\varGamma }{2}(U_*+\varLambda U_{\lambda })\Big]} {\big[\sigma + ik(3-ik)U_{*}\big]} - \frac {3}{2}ik(1-ik)\varGamma \varLambda U_{*} U_{\lambda } = 0 , \end{equation}
with
\begin{equation} U_\lambda = \frac 13 A'(\varLambda ) = -\frac {1}{3} [2(1-\varLambda )(1-a)+a] . \end{equation}
Solutions to this dispersion relation are illustrated in figure 14. Flow is again unstable for sufficiently long relaxation time, and the stability characteristics are similar to those found for finite
$\kappa$
.
For long waves, the dispersion relation simplifies further: we find
$\sigma \sim -ik c + \sigma _2 k^2$
, where
\begin{equation} c = 3U_* \left [ 1- \frac {\varGamma \varLambda U_\lambda }{2+3\varGamma (U_*+\varLambda U_\lambda )} \right]\!, \qquad \sigma _2 = \frac {2{\mathcal T}(c-U_*)(c-3U_*) - U_*(2+3\varGamma U_*)} {2+3\varGamma (U_*+\varLambda U_\lambda )} . \end{equation}
The
$O(k^2)$
growth rate again exposes the destabilizing effect of relaxation (
$c\lt 3U_*$
and
$2+3\varGamma (U_*+\varLambda U_\lambda )\gt 0$
) and the stabilization of quasisteady viscous effects.
(a) Growth rates and (b) phase speeds as functions of
$k$
in the strong-diffusion limit for
$(a,\varGamma ,{\mathcal T})=({1/5},8,14)$
. The (red) solid lines show the predictions of (E5), the (blue) points show numerical solutions of the full linear stability problem with
$\kappa =10$
. The remaining panels show density plots of growth rate (from (E5)) over the (a)
$(k,{\mathcal T})$
(
$a={1/5}$
), (b)
$(k,a)$
(
${\mathcal T}=14$
) and (c)
$({\mathcal T},a)$
planes (
$k=0.2$
), all with
$\varGamma =8$
. The (red) solid lines in (c), (d) and (e) show the stability boundaries and the dots show the maximum growth over all
$k$
or
$\mathcal T$
.
As shown in figure 15, nonlinear solutions to (E3) in periodic spatial domains initially demonstrate the growth predicted by (E5), but then saturate into steady nonlinear waves. The initial condition used here seeds the growth of the most unstable linear mode (with
$k=0.2$
; see figure 15
a). The nonlinear waves that emerge, as illustrated in figure 15(b,c), display destructuring underneath the rise of the wave crests, and structural recovery over the fall to the following trough. These nonlinear waves appear in a supercritical bifurcation as
$a$
is lowered through the stability threshold.
A strong-diffusion solution to (E3) in a periodic domain of length
$20\pi$
, with
$(a,\varGamma ,{\mathcal T})=( 1/2,8,14)$
. The initial condition is a random superposition of the first 10 Fourier modes with mean amplitudes of
$10^{-5}$
. Shown are (a) a time series of
$\sqrt {\langle (h-1)^2\rangle }$
, and snapshots of (b)
$h(x,t)$
and (c)
$\lambda (x,t)$
(every five time units, with the final profile shown darker). The dashed line in (a) indicates the growth expected for the most unstable mode (with
$k=0.2$
). The snapshots are shown in a moving frame with speed
$V=0.64$
. The final wave profiles from a wider suite of computations with varying
$a$
are shown below, plotting (d) wave amplitude
$\Delta h \equiv h_{max}-h_\textit{min}$
against
$a$
, and (e)
$h$
and (f)
$\lambda$
for the values of
$a$
indicated by the stars in (d). The grey dashed line shows the stability boundary in (d), and the profiles are shifted to align the wave crests (i.e.
$h_{max}$
) in (e) and (f). The dot–dash line in (a), circle in (d) and dashed lines in (e) and (f) show corresponding results from the full thin-film model, with
$\kappa =10$
(plotting
$\langle \lambda \rangle$
in (f)); the most obvious differences arise during the initial transients.
Appendix F. Details of the numerical scheme used to compute nonlinear solutions
We first map the fluid domain to a rectangle by transforming to the new spatial variables
$(x,\zeta )$
, where
$\zeta =z/h(x,t)$
). After discretizing the new variables, spatial derivatives are evaluated using spectral differentiation matrices (supplementing Chebyshev matrices for
$\partial /\partial \zeta$
with matrices for
$\partial /\partial x$
based on the fast Fourier transform (Weideman & Reddy Reference Weideman and Reddy2000)). The discrete spatial grid in
$(x,\zeta )$
is then an
$N\times M$
collection of collocation points, where practically we take
$N=M=128$
or more. The resulting set of coupled ordinary differential equations are solved using MATLAB’s ODE15s. As initial conditions, we adopt
$h(x,0)=1 + A\sin kx$
and
$\lambda (x,\zeta ,0) = \varLambda (\zeta )$
, where
$A$
is typically chosen to be
$5\times 10^{-4}$
and
$\varLambda (z)$
is one of the base-flow solutions of § 3.
To examine the effect of spatial resolution, solutions with different choices for
$N$
and
$M$
(but setting
$N=M$
) were computed and compared. For the same two solutions of figure 9, the final amplitude of
$\sqrt {\langle (h-1)^2\rangle }$
changes by less than a per cent when halving or even quartering the number of grid points; the changes to the final nonlinear wave speed are smaller still.
Note that MATLAB’s ODE15s uses preset relative and absolute error tolerances of
$10^{-3}$
and
$10^{-6}$
, respectively, by default. These were reset to
$10^{-5}$
and
$10^{-9}$
in order to provide additional assurance about the fidelity of the time integration. Nevertheless, the two solutions of figure 9 changed by less than a tenth of a per cent when recomputed with the original tolerances.
Figure 9(b) also compares the two nonlinear computations for
$N=M=128$
with the predictions from linear theory; over the period of exponential growth, the computations match the linear predictions. In more detail, the growth rate extracted from the time series of the nonlinear solution agrees with the linear eigenvalue to better than a per cent. In summary, there was no evidence of any issue with spatial resolution or time integration error.
Open access
