2 Gradient regularization of $\unicode[STIX]{x1D707}(I)$

As the method we employ applies to general models of viscoplasticity, we offer a derivation of constitutive equation and momentum balance based on the notion of a dissipation potential and the associated hyperstresses and variational principle. Readers not interested in these details can skip immediately to the momentum balance presented below in (2.11).

Consider a strictly dissipative material endowed with frame-indifferent dissipation potential (Edelen Reference Edelen1972, Reference Edelen2005; Goddard Reference Goddard2014; Saramito Reference Saramito2016) depending on the first two spatial gradients of the material velocity field $\boldsymbol{v}(t,\boldsymbol{x})$ ,

(2.1) $$\begin{eqnarray}\unicode[STIX]{x1D713}(\unicode[STIX]{x1D63F},\unicode[STIX]{x1D735}\unicode[STIX]{x1D735}\boldsymbol{v})=\unicode[STIX]{x1D713}(\unicode[STIX]{x1D64C}\unicode[STIX]{x1D63F}\unicode[STIX]{x1D64C}^{\text{T}},\unicode[STIX]{x1D64C}\unicode[STIX]{x1D735}\unicode[STIX]{x1D64C}\unicode[STIX]{x1D735}\boldsymbol{v}\unicode[STIX]{x1D64C}^{\text{T}}),\quad \text{with }\unicode[STIX]{x1D63F}={\textstyle \frac{1}{2}}(\unicode[STIX]{x1D735}\boldsymbol{v}+(\unicode[STIX]{x1D735}\boldsymbol{v})^{\text{T}}),\end{eqnarray}$$

where $\unicode[STIX]{x1D63F}$ is the deformation rate and $\unicode[STIX]{x1D64C}=\unicode[STIX]{x1D64C}(t)$ is an arbitrary time-dependent but spatially independent orthogonal tensor.

A simple special case is represented by the vdW–C–H equation. (The terminology ‘van der Waals–Cahn–Hilliard’ is more accurate historically (Rowlinson Reference Rowlinson1979) than the oft-used ‘Cahn–Hilliard’. We recall that Korteweg proposed a more general stress arising from a Helmholtz free energy contribution of the form $\unicode[STIX]{x1D70C}_{,i}\unicode[STIX]{x1D70C}_{,j}\unicode[STIX]{x1D70C}_{,ij}$ .) Here,

(2.2) $$\begin{eqnarray}\unicode[STIX]{x1D713}(\unicode[STIX]{x1D63F},\unicode[STIX]{x1D735}\unicode[STIX]{x1D735}\boldsymbol{v})=\unicode[STIX]{x1D713}_{0}(\unicode[STIX]{x1D63F})+\unicode[STIX]{x1D712}|\unicode[STIX]{x1D735}\unicode[STIX]{x1D735}\boldsymbol{v}|^{2},\end{eqnarray}$$

where we take $\unicode[STIX]{x1D712}$ to be a positive constant, and the term in the higher-order velocity gradient $|\unicode[STIX]{x1D735}\unicode[STIX]{x1D735}\boldsymbol{v}|^{2}=(\unicode[STIX]{x1D735}\unicode[STIX]{x1D735}\boldsymbol{v})\boldsymbol{\cdot }(\unicode[STIX]{x1D735}\unicode[STIX]{x1D735}\boldsymbol{v})=v_{i,jk}v_{i,jk}$ serves to provide a local regularization whenever $\unicode[STIX]{x1D713}_{0}$ becomes locally non-convex. We recall that the gradient term in the standard scalar version of the vdW–C–H equation involves fluid density in place of $\unicode[STIX]{x1D735}\boldsymbol{v}$ (Cahn & Hilliard Reference Cahn and Hilliard1958; Gurtin Reference Gurtin1996), providing an energy penalty on large gradients. This yields well-posed field equations in the case of non-convex free energy with phase transition and serves to define the equilibrium surface tension in the limit $\unicode[STIX]{x1D712}\rightarrow 0$ .

In the following, we deal with strongly dissipative materials, i.e. strictly dissipative materials devoid of gyroscopic or ‘powerless’ stress (Goddard Reference Goddard2014), for which the partial derivatives of the dissipation potential $\unicode[STIX]{x1D713}(\unicode[STIX]{x1D63F},\unicode[STIX]{x1D735}\unicode[STIX]{x1D735}\boldsymbol{v})$ yield work-conjugate (symmetric Cauchy) stress $\unicode[STIX]{x1D64F}^{(1)}$ and hyperstress $\unicode[STIX]{x1D64F}^{(2)}$ respectively, according to

(2.3a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D64F}^{(1)}=\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D63F}}=\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D63F}}\quad \text{and}\quad \unicode[STIX]{x1D64F}^{(2)}=\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}(\unicode[STIX]{x1D735}\unicode[STIX]{x1D735}\boldsymbol{v})}=2\unicode[STIX]{x1D712}\unicode[STIX]{x1D735}\unicode[STIX]{x1D735}\boldsymbol{v},\end{eqnarray}$$

and the volumetric rate of dissipation is given by

(2.4) $$\begin{eqnarray}\mathfrak{D}=\unicode[STIX]{x1D64F}^{(1)}\boldsymbol{\cdot }\unicode[STIX]{x1D63F}+\unicode[STIX]{x1D64F}^{(2)}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D735}\boldsymbol{v}=\unicode[STIX]{x1D61B}_{ij}^{(1)}\unicode[STIX]{x1D60B}_{ij}+\unicode[STIX]{x1D61B}_{ijk}^{(2)}v_{k,ij}.\end{eqnarray}$$

The hyperstress $\unicode[STIX]{x1D64F}^{(2)}=[\unicode[STIX]{x1D61B}_{ijk}^{(2)}]$ represents a generalized ‘pinch’ $\boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D64F}^{(2)}=[n_{i}\unicode[STIX]{x1D61B}_{ijk}^{(2)}]$ acting on a material plane with unit normal $\boldsymbol{n}$ , which can be represented by a force dipole consisting of equal and opposite forces. The symmetric part $[n_{i}\unicode[STIX]{x1D61B}_{i(jk)}^{(2)}]$ involves forces acting along the line of centres of their points of application, representing a normal pinch or symmetric ‘stresslet’, whereas the antisymmetric part $[n_{i}\unicode[STIX]{x1D61B}_{i[jk]}^{(2)}]$ involves forces acting perpendicular to their line of centres, representing a ‘torque’ or ‘rotlet’. It should be noted that the gradient of vorticity enters into the mechanical power and that most of the concepts of kinematics and stress carry over from the various works on the elasticity of solids, where the strain energy (Helmholtz free energy) rather than the dissipation potential is involved. (Whereas the theory for elastic solids or fluids leads to elastic surface tension, the hyper-dissipative model can apparently yield a surface viscosity, which as far as we know would constitute a novel continuum approach to the subject.)

Thus, upon enslaving microstructural to continuum kinematics as in the classic work of Mindlin (Reference Mindlin1964), achieved by taking his relative displacement gradient $\unicode[STIX]{x1D6FE}_{ij}=0$ , one obtains the linear-elastic analogue of the present work, with his strains $\unicode[STIX]{x1D716}_{ij}$ replacing our strain rates $\unicode[STIX]{x1D60B}_{ij}$ . Our quasi-static equation of equilibrium (2.7) follows from that of Mindlin upon relaxing the assumption of incompressibility and combining his equations (4.1), effectively eliminating his ‘relative stress’ $\unicode[STIX]{x1D70E}_{ij}$ . Mindlin’s work also shows that, within a linear isotropic gradient model, one may anticipate further quadratic terms in $\unicode[STIX]{x1D735}\unicode[STIX]{x1D735}\boldsymbol{v}$ beyond that adopted in the simpler one-parameter model (2.2) of the present study.

The Hamiltonian momentum balances for the hyper-elastic system of Mindlin do not apply to hyper-dissipative systems, except in the limit of quasi-static (i.e. inertialess) motion. Hence, further analysis is required to obtain the relevant balances for the latter, and we begin with the variational principle leading to the quasi-static balance. We note that similar methods have been adopted in past works of Hill (Reference Hill1956) and Leonov (Reference Leonov1988), methods that are made rigorous mathematically by the later works of Edelen, e.g. Edelen (Reference Edelen1972, Reference Edelen2005), which are highlighted and simplified in the survey by Goddard (Reference Goddard2014).

It should be noted that the variational derivative of the functional of $\boldsymbol{v}(\boldsymbol{x},t)$ representing the global dissipation potential,

(2.5) $$\begin{eqnarray}\unicode[STIX]{x1D6F9}[\boldsymbol{v}]=\int _{V}[\unicode[STIX]{x1D713}(\unicode[STIX]{x1D63F},\unicode[STIX]{x1D735}\unicode[STIX]{x1D735}\boldsymbol{v})-p\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{v}]\,\text{d}V,\end{eqnarray}$$

subject to incompressibility $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{v}=0$ in the spatial domain $V$ , is given by

(2.6) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FF}\unicode[STIX]{x1D6F9}[\boldsymbol{v}] & = & \displaystyle \int _{V}\left[\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D63F}}\boldsymbol{\cdot }\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D735}\boldsymbol{v}+\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}(\unicode[STIX]{x1D735}\unicode[STIX]{x1D735}\boldsymbol{v})}\boldsymbol{\cdot }\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D735}\unicode[STIX]{x1D735}\boldsymbol{v}-p\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D6FF}\boldsymbol{v}\right]\,\text{d}V\nonumber\\ \displaystyle & = & \displaystyle \int _{V}[\unicode[STIX]{x1D61B}_{ij}^{(1)}\unicode[STIX]{x1D6FF}v_{i,j}+\unicode[STIX]{x1D61B}_{ijk}^{(2)}\unicode[STIX]{x1D6FF}v_{i,jk}-p\unicode[STIX]{x1D6FF}v_{j,j}]\,\text{d}V\nonumber\\ \displaystyle & = & \displaystyle \int _{V} [(\unicode[STIX]{x1D61B}_{ij}^{(1)}\unicode[STIX]{x1D6FF}v_{i})_{,j}-\unicode[STIX]{x1D61B}_{ij,j}^{(1)}\unicode[STIX]{x1D6FF}v_{i}+(\unicode[STIX]{x1D61B}_{ijk}^{(2)}\unicode[STIX]{x1D6FF}v_{i,j})_{,k}\nonumber\\ \displaystyle & & \displaystyle -\,(\unicode[STIX]{x1D61B}_{ijk,k}^{(2)}\unicode[STIX]{x1D6FF}v_{i})_{,j}+\unicode[STIX]{x1D61B}_{ijk,jk}^{(2)}\unicode[STIX]{x1D6FF}v_{i}-(p\unicode[STIX]{x1D6FF}v_{j})_{,j}+p_{,j}\unicode[STIX]{x1D6FF}v_{j}\! ]\,\text{d}V\nonumber\\ \displaystyle & = & \displaystyle \int _{\unicode[STIX]{x2202}V}\boldsymbol{n}\boldsymbol{\cdot }[(\unicode[STIX]{x1D64F}^{(1)}-\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D64F}^{(2)}-p\unicode[STIX]{x1D644})\boldsymbol{\cdot }\unicode[STIX]{x1D6FF}\boldsymbol{v}+\unicode[STIX]{x1D64F}^{(2)}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FF}\boldsymbol{v}]\,\text{d}S\nonumber\\ \displaystyle & & \displaystyle -\,\int _{V}[\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D64F}^{(1)}-\unicode[STIX]{x1D735}p-(\unicode[STIX]{x1D735}\unicode[STIX]{x1D735})\boldsymbol{\cdot }\unicode[STIX]{x1D64F}^{(2)}]\boldsymbol{\cdot }\unicode[STIX]{x1D6FF}\boldsymbol{v}\,\text{d}V,\end{eqnarray}$$

where pressure $p$ plays its usual role as Lagrange multiplier and use has been made of the divergence theorem for integration by parts. Noting that $\unicode[STIX]{x1D6FF}\boldsymbol{v}=0$ on $\unicode[STIX]{x2202}V$ implies that the surface-tangential gradient of $\unicode[STIX]{x1D6FF}\boldsymbol{v}$ vanishes on $\unicode[STIX]{x2202}V$ , this establishes the following variational theorem.

The latter is given by (2.6) as

(2.7) $$\begin{eqnarray}\displaystyle & & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D64F}^{(1)}-\unicode[STIX]{x1D735}p-(\unicode[STIX]{x1D735}\unicode[STIX]{x1D735})\boldsymbol{\cdot }\unicode[STIX]{x1D64F}^{(2)}=0,\quad \text{i.e. }\unicode[STIX]{x1D61B}_{ij,j}^{(1)}-p_{,i}-\unicode[STIX]{x1D61B}_{jki,jk}^{(2)}=0,\nonumber\\ \displaystyle & & \displaystyle \quad \text{or }\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D64F}-\unicode[STIX]{x1D735}p=0,\quad \text{where }\unicode[STIX]{x1D64F}:=\unicode[STIX]{x1D64F}^{(1)}-\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D64F}^{(2)},\quad \text{i.e. }\unicode[STIX]{x1D61B}_{ij}=\unicode[STIX]{x1D61B}_{ij}^{(1)}-\unicode[STIX]{x1D61B}_{jki,k}^{(2)},\qquad\end{eqnarray}$$

where $\unicode[STIX]{x1D64F}$ obviously serves as the effective stress tensor. Despite the ostensible reduction to a single stress tensor, we should re-emphasize that the hyperstress $\unicode[STIX]{x1D64F}^{(2)}$ can give rise to singular surface stresses balancing discontinuities in the Cauchy stress $\unicode[STIX]{x1D64F}^{(1)}$ , as mentioned below in our analysis of shear bands. We further note that the presence of such effects at the nominal free surface of thin avalanching layers could invalidate theories based on variants of the $\unicode[STIX]{x1D707}(I)$ model, as already suggested by certain strongly non-local models (Henann & Kamrin Reference Henann and Kamrin2014).

Now, one can extend (2.7) to include gravitation or other fixed body forces $\boldsymbol{g}$ by replacing $\unicode[STIX]{x1D713}$ with $\unicode[STIX]{x1D713}-\unicode[STIX]{x1D719}\unicode[STIX]{x1D70C}_{s}\boldsymbol{g}\boldsymbol{\cdot }\boldsymbol{v}$ in the first term of (2.6). Then, by a further appeal to d’Alembert’s principle, one can then replace $\boldsymbol{g}$ with the total acceleration $\boldsymbol{g}-d_{t}\boldsymbol{v}$ in the resulting equation of equilibrium to obtain the complete linear momentum balance,

(2.8a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D70C}_{s}\unicode[STIX]{x1D719}d_{t}\boldsymbol{v}=-\unicode[STIX]{x1D735}p+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D64F}+\unicode[STIX]{x1D70C}_{s}\unicode[STIX]{x1D719}\boldsymbol{g}\quad \text{and}\quad \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{v}=0,\quad \text{with }d_{t}=\unicode[STIX]{x2202}_{t}+\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735},\end{eqnarray}$$

a relation that no longer follows directly from the above extremum principle. Here, $\unicode[STIX]{x1D70C}_{s}$ is the constant solid density of the grains, $p$ is the pressure, $\boldsymbol{g}$ is the gravitational acceleration and $d_{t}$ is the material (or ‘substantial’) time derivative. Thus, with reinterpretation of the stress tensor $\unicode[STIX]{x1D64F}$ , the linear momentum balance retains the same form as for a simple (‘non-polar’) material.

2.1 Application to the $\unicode[STIX]{x1D707}(I)$ model

In the model of Jop et al. (Reference Jop, Forterre and Pouliquen2006) (cf. Ref. 1) for granular rheology, the Cauchy stress $\unicode[STIX]{x1D64F}^{(1)}$ is given by a rate-dependent version of Drucker–Prager plasticity with rate-dependent friction coefficient $\unicode[STIX]{x1D707}(I)$ ,

(2.9a ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D64F}^{(1)}=\unicode[STIX]{x2202}_{\unicode[STIX]{x1D63F}}\unicode[STIX]{x1D713}_{0}=\unicode[STIX]{x1D707}(I)p\unicode[STIX]{x1D640},\quad \text{with }\unicode[STIX]{x1D707}(I)=\unicode[STIX]{x1D707}_{0}+\frac{(\unicode[STIX]{x1D707}_{\infty }-\unicode[STIX]{x1D707}_{0})}{(I+I_{\ast })}I,\quad \text{and} & \displaystyle\end{eqnarray}$$

(2.9b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D713}_{0}=\frac{p}{\unicode[STIX]{x1D703}}\left[\unicode[STIX]{x1D707}_{\infty }I+(\unicode[STIX]{x1D707}_{0}-\unicode[STIX]{x1D707}_{\infty })I_{\ast }\ln \left(\frac{I+I_{\ast }}{I_{\ast }}\right)\right]. & \displaystyle\end{eqnarray}$$

Here, the normalized form or ‘director’ $\unicode[STIX]{x1D640}$ , the Euclidean norm of the strain rate $|\unicode[STIX]{x1D63F}|$ and the inertial number are defined respectively by

(2.10a-c ) $$\begin{eqnarray}\unicode[STIX]{x1D640}:=\frac{\unicode[STIX]{x1D63F}}{|\unicode[STIX]{x1D63F}|},\quad |\unicode[STIX]{x1D63F}|:=\sqrt{\text{tr}(\unicode[STIX]{x1D63F}^{2})}\quad \text{and}\quad I=\unicode[STIX]{x1D703}|\unicode[STIX]{x1D63F}|,\quad \text{with }\unicode[STIX]{x1D703}=d\sqrt{2\unicode[STIX]{x1D70C}_{s}/p},\end{eqnarray}$$

where $I_{\ast }$ is an empirical constant (denoted by $\unicode[STIX]{x1D610}_{0}$ in several previous works), $d$ is a representative grain diameter, $p$ is the local pressure, $\unicode[STIX]{x1D703}$ plays the role of an inertial relaxation time and the quantities $\unicode[STIX]{x1D707}_{\infty }\geqslant \unicode[STIX]{x1D707}_{0}$ represent limiting friction coefficients. One obtains the standard version adapted to simple shear upon replacing the Euclidean norm by the norm $\Vert \unicode[STIX]{x1D63F}\Vert =|\unicode[STIX]{x1D63F}|/\sqrt{2}$ employed in Jop et al. (Reference Jop, Forterre and Pouliquen2006) and Ref. 1.

We note that $\unicode[STIX]{x1D713}_{0}$ is marginally convex in the special case of constant $\unicode[STIX]{x1D707}$ (Goddard Reference Goddard2014), which we believe accounts for the general tendency to material instability in perfectly plastic models. The potential multiplicity of solutions will be made more evident by the model of shear banding presented below in § 4.

The momentum balance (2.8) can now be recast as

(2.11) $$\begin{eqnarray}\unicode[STIX]{x1D70C}_{s}\unicode[STIX]{x1D719}d_{t}\boldsymbol{v}=-\unicode[STIX]{x1D735}p+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D707}(I)p\unicode[STIX]{x1D640})-2\unicode[STIX]{x1D712}\unicode[STIX]{x1D6FB}^{4}\boldsymbol{v}+\unicode[STIX]{x1D70C}_{s}\unicode[STIX]{x1D719}\boldsymbol{g}.\end{eqnarray}$$

The third term on the right-hand side of (2.11) arises from the hyperstress $\unicode[STIX]{x1D64F}^{(2)}$ , and the momentum balance of Ref. 1 is obtained by taking $\unicode[STIX]{x1D712}\equiv 0$ . The preceding term represents the standard Cauchy stress and, confirming the analysis of Ref. 1, can be expanded to yield

(2.12) $$\begin{eqnarray}\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D707}(I)p\unicode[STIX]{x1D640})=\frac{\unicode[STIX]{x1D707}(I)(2-\stackrel{\hspace{2.22198pt}\prime }{\unicode[STIX]{x1D707}})}{2}\unicode[STIX]{x1D640}\unicode[STIX]{x1D735}p+(\stackrel{\hspace{2.22198pt}\prime }{\unicode[STIX]{x1D707}}-1)\frac{\unicode[STIX]{x1D707}(I)p}{|\unicode[STIX]{x1D63F}|}[(\unicode[STIX]{x1D640}\unicode[STIX]{x1D735})\unicode[STIX]{x1D735}\boldsymbol{v}]\boldsymbol{\cdot }\unicode[STIX]{x1D640}+\frac{\unicode[STIX]{x1D707}(I)p}{2|\unicode[STIX]{x1D63F}|}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{v}.\end{eqnarray}$$

Then, the momentum balance in (2.8) becomes

(2.13) $$\begin{eqnarray}\unicode[STIX]{x1D70C}_{s}\unicode[STIX]{x1D719}d_{t}\boldsymbol{v}=-\unicode[STIX]{x1D649}\unicode[STIX]{x1D735}p+\frac{(\stackrel{\hspace{2.22198pt}\prime }{\unicode[STIX]{x1D707}}-1)\unicode[STIX]{x1D707}p}{|\unicode[STIX]{x1D63F}|}[(\unicode[STIX]{x1D640}\unicode[STIX]{x1D735})\unicode[STIX]{x1D735}\boldsymbol{v}]\boldsymbol{\cdot }\unicode[STIX]{x1D640}+\frac{\unicode[STIX]{x1D707}p}{2|\unicode[STIX]{x1D63F}|}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{v}-2\unicode[STIX]{x1D712}\unicode[STIX]{x1D6FB}^{4}\boldsymbol{v}+\unicode[STIX]{x1D70C}_{s}\unicode[STIX]{x1D719}\boldsymbol{g},\end{eqnarray}$$

where we employ the notation

(2.14a-c ) $$\begin{eqnarray}\unicode[STIX]{x1D649}=\unicode[STIX]{x1D644}-\frac{(2-\stackrel{\hspace{2.22198pt}\prime }{\unicode[STIX]{x1D707}})\unicode[STIX]{x1D707}}{2}\unicode[STIX]{x1D640},\quad \stackrel{\hspace{2.22198pt}\prime }{\unicode[STIX]{x1D707}}=\frac{I}{\unicode[STIX]{x1D707}(I)}\frac{\text{d}\unicode[STIX]{x1D707}(I)}{\text{d}I}=\frac{\text{d}\log \unicode[STIX]{x1D707}(I)}{\text{d}\log I}\quad \text{and}\quad \stackrel{\hspace{2.22198pt}\prime \prime }{\unicode[STIX]{x1D707}}=\frac{d^{2}\unicode[STIX]{x1D707}(I)}{\text{d}I^{2}}\frac{I^{2}}{\unicode[STIX]{x1D707}(I)}\end{eqnarray}$$

here and below.

3 Linear-stability analysis

In the usual way, the perturbed velocity $\boldsymbol{v}$ and pressure are written as

(3.1a,b ) $$\begin{eqnarray}\boldsymbol{v}=\boldsymbol{v}^{(0)}+\boldsymbol{v}^{(1)}\quad \text{and}\quad p=p^{(0)}+p^{(1)},\end{eqnarray}$$

where the superscripts $(0)$ and $(1)$ denote the base state and the perturbation respectively. The perturbed friction coefficients are given by

(3.2a ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D707}=\unicode[STIX]{x1D707}^{(0)}+\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}{\unicode[STIX]{x2202}p}\right)^{(0)}p^{(1)}+\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D63F}}\right)^{(0)}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{v}^{(1)}\quad \text{and} & \displaystyle\end{eqnarray}$$

(3.2b ) $$\begin{eqnarray}\displaystyle & \displaystyle \stackrel{\hspace{2.22198pt}\prime }{\unicode[STIX]{x1D707}}\,=\,\stackrel{\hspace{2.22198pt}\prime }{\unicode[STIX]{x1D707}}^{(0)}+\left(\frac{\unicode[STIX]{x2202}\stackrel{\hspace{2.22198pt}\prime }{\unicode[STIX]{x1D707}}}{\unicode[STIX]{x2202}p}\right)^{(0)}p^{(1)}+\left(\frac{\unicode[STIX]{x2202}\stackrel{\hspace{2.22198pt}\prime }{\unicode[STIX]{x1D707}}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D63F}}\right)^{(0)}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{v}^{(1)}. & \displaystyle\end{eqnarray}$$

(3.3) $$\begin{eqnarray}\displaystyle & & \displaystyle \unicode[STIX]{x1D70C}_{s}\unicode[STIX]{x1D719}(d_{t}^{(0)}\boldsymbol{v}^{(1)}+\boldsymbol{v}^{(1)}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{v}^{(0)})=-\unicode[STIX]{x1D649}^{(0)}\unicode[STIX]{x1D735}p^{(1)}\nonumber\\ \displaystyle & & \displaystyle \quad +\,(\stackrel{\hspace{2.22198pt}\prime }{\unicode[STIX]{x1D707}}^{(0)}-1)\frac{\unicode[STIX]{x1D707}^{(0)}p^{(0)}}{|\unicode[STIX]{x1D63F}^{(0)}|}(\unicode[STIX]{x1D640}^{(0)}\unicode[STIX]{x1D735})(\unicode[STIX]{x1D640}^{(0)}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{v}^{(1)})+\frac{\unicode[STIX]{x1D707}^{(0)}p^{(0)}}{2|\unicode[STIX]{x1D63F}^{(0)}|}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{v}^{(1)}-2\unicode[STIX]{x1D712}\unicode[STIX]{x1D6FB}^{4}\boldsymbol{v}^{(1)}\nonumber\\ \displaystyle & & \displaystyle \quad +\,\left[\frac{\unicode[STIX]{x1D707}^{(0)}(\stackrel{\hspace{2.22198pt}\prime \prime }{\unicode[STIX]{x1D707}}^{(0)}-\stackrel{\hspace{2.22198pt}\prime }{\unicode[STIX]{x1D707}}^{(0)})}{4p^{(0)}}p^{(1)}+\frac{\unicode[STIX]{x1D707}^{(0)}}{2|\unicode[STIX]{x1D63F}^{(0)}|}(\stackrel{\hspace{2.22198pt}\prime }{\unicode[STIX]{x1D707}}^{(0)}-\stackrel{\hspace{2.22198pt}\prime \prime }{\unicode[STIX]{x1D707}}^{(0)})\unicode[STIX]{x1D640}^{(0)}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{v}^{(1)}\right]\unicode[STIX]{x1D640}^{(0)}\unicode[STIX]{x1D735}p^{(0)},\nonumber\\ \displaystyle & & \displaystyle \hspace{30.00005pt}\text{where }\unicode[STIX]{x1D649}^{(0)}=\unicode[STIX]{x1D644}-\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D640}^{(0)},\quad \text{with }\unicode[STIX]{x1D6FC}=\left[\frac{(2-\stackrel{\hspace{2.22198pt}\prime }{\unicode[STIX]{x1D707}})}{2}\unicode[STIX]{x1D707}\right]^{(0)},\end{eqnarray}$$

and the equations of Ref. 1 are obtained by taking $\unicode[STIX]{x1D735}p^{(0)}=\mathbf{0}$ and $\unicode[STIX]{x1D712}=0$ .

With $|\unicode[STIX]{x1D63F}^{(0)}|^{-1}$ as time scale and $d$ as length scale, we henceforth adopt the following non-dimensional variables:

(3.4) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \bar{\boldsymbol{x}}=\frac{\boldsymbol{x}}{d},\quad \bar{t}=\surd 2|\unicode[STIX]{x1D63F}^{(0)}|t,\quad \bar{\boldsymbol{v}}=\frac{\boldsymbol{v}}{\surd 2|\unicode[STIX]{x1D63F}^{(0)}|d},\\ \displaystyle \bar{p}=\frac{p}{2\unicode[STIX]{x1D70C}_{s}d^{2}|\unicode[STIX]{x1D63F}^{(0)}|^{2}},\quad \bar{\unicode[STIX]{x1D712}}=\frac{\unicode[STIX]{x1D712}}{\surd 2\unicode[STIX]{x1D70C}_{s}d^{4}|\unicode[STIX]{x1D63F}^{(0)}|},\quad \bar{g}=\frac{g}{2d|\unicode[STIX]{x1D63F}^{(0)}|^{2}}.\end{array}\right\}\end{eqnarray}$$

Various factors of $\surd 2$ , included here to simplify the results presented below for simple shear, can be eliminated by substituting the norm $\Vert \unicode[STIX]{x1D63F}\Vert =|\unicode[STIX]{x1D63F}|/\surd 2$ employed in Ref. 1 and mentioned above. For later reference, we note that inspection of the combination of terms $\unicode[STIX]{x1D6FB}^{2}\boldsymbol{v}$ and $\unicode[STIX]{x1D6FB}^{4}\boldsymbol{v}$ in (2.13) indicates that $\bar{\unicode[STIX]{x1D712}}$ is proportional to the square of a ratio of microscopic to macroscopic length scales.

Employing the above non-dimensional variables and dropping the superimposed bars for simplicity, we can write (3.3) in component form as

(3.5) $$\begin{eqnarray}\displaystyle & & \displaystyle \unicode[STIX]{x1D719}[(d_{t}^{(0)}\boldsymbol{v}^{(1)})_{i}+v_{i,j}^{(0)}v_{j}^{(1)}]\nonumber\\ \displaystyle & & \displaystyle \quad =-\unicode[STIX]{x1D615}_{ij}^{(0)}p_{,j}^{(1)}-2\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D60C}_{ij}^{(0)}\unicode[STIX]{x1D60C}_{kl}^{(0)}v_{k,jl}^{(1)}+\unicode[STIX]{x1D6FE}v_{i,jj}^{(1)}-2\unicode[STIX]{x1D712}v_{i,jjll}^{(1)},\text{and }v_{j,j}^{(1)}=0,\nonumber\\ \displaystyle & & \displaystyle \qquad \text{with }\unicode[STIX]{x1D615}_{ij}^{(0)}=\unicode[STIX]{x1D6FF}_{ij}-\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D60C}_{ij}^{(0)},\unicode[STIX]{x1D6FD}=1-\stackrel{\hspace{2.22198pt}\prime }{\unicode[STIX]{x1D707}}^{(0)}\text{ and }\unicode[STIX]{x1D6FE}=\frac{\unicode[STIX]{x1D707}^{(0)}p^{(0)}}{\surd 2},\end{eqnarray}$$

for the case of uniform base pressure $p^{(0)}$ .

The parameters $\stackrel{\hspace{2.22198pt}\prime }{\unicode[STIX]{x1D707}}^{(0)},\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD}$ and $\unicode[STIX]{x1D6FE}$ of this study are related to the parameters $\unicode[STIX]{x1D708},q,r$ and $\unicode[STIX]{x1D702}^{(0)}$ of Ref. 1 respectively by

(3.6a-d ) $$\begin{eqnarray}\stackrel{\hspace{2.22198pt}\prime }{\unicode[STIX]{x1D707}}^{(0)}=\unicode[STIX]{x1D708},\quad \unicode[STIX]{x1D6FC}=q=\unicode[STIX]{x1D707}^{(0)},\quad \unicode[STIX]{x1D6FD}=r=1-\unicode[STIX]{x1D708}\quad \text{and}\quad \unicode[STIX]{x1D6FE}=\unicode[STIX]{x1D702}^{(0)}.\end{eqnarray}$$

We employ $\unicode[STIX]{x1D640}$ for the normalized deformation rate denoted by $\unicode[STIX]{x1D63C}$ in Ref. 1, with different norms as discussed above, and the vector $\boldsymbol{k}$ for the wavenumber in lieu of $\unicode[STIX]{x1D743}$ of Ref. 1, which is defined below as $\imath \boldsymbol{k}$ .

It should be noted that the quantities $\unicode[STIX]{x1D6FE}/\unicode[STIX]{x1D719}$ and $\unicode[STIX]{x1D712}/\unicode[STIX]{x1D719}$ represent ratios of inertial to dissipative forces, replacing the inverse of the Reynolds number for Newtonian fluids. Since these quantities are $O(1)$ in the following analysis, inertia can be assumed to be relatively unimportant. Hence, we are led to characterize the instability found in Ref. 1 and analysed below as material instability, of a type that is already manifest in the quasi-static equations obtained by taking $\unicode[STIX]{x1D719}\equiv 0$ on the left-hand side of (2.11) or (3.5). As discussed below, this type of instability can be attributed to the loss of generalized ellipticity of the highest-order differential operators on the right-hand side of (3.5), as indicated, e.g., by the classical analysis of Browder (Browder Reference Browder1961; Brezis & Browder Reference Brezis and Browder1998). The first of these papers indicates clearly the connection to the transient instability found in Ref. 1. This type of linear instability, with initially large transient growth rates, may of course be more relevant to the numerical simulation of nonlinear instability than the asymptotic linear stability established analytically by the following analysis. In this respect, the transient instability bears a certain resemblance to that arising from the non-normality of linear-stability operators (Trefethen & Embree Reference Trefethen and Embree2005) which can trigger nonlinear effects.

We now consider a base state $\boldsymbol{v}^{(0)}=\unicode[STIX]{x1D647}^{(0)}\boldsymbol{x}$ with spatially uniform velocity gradient $\unicode[STIX]{x1D647}^{(0)}=(\unicode[STIX]{x1D735}\boldsymbol{v}^{(0)})^{\text{T}}$ . Then, with the Fourier-space representation $\hat{f}(\boldsymbol{k})$ of fields $f(\boldsymbol{x})$ and the duality $\boldsymbol{x}\leftrightarrow \imath \hat{\unicode[STIX]{x1D735}}$ and $\unicode[STIX]{x1D735}\leftrightarrow \imath \boldsymbol{k}$ , (3.5) can be written as

(3.7) $$\begin{eqnarray}\displaystyle & & \displaystyle \hat{d}_{t}^{(0)}\hat{\boldsymbol{v}}^{(1)}=\left[\unicode[STIX]{x1D648}-\frac{1}{\unicode[STIX]{x1D719}}(\unicode[STIX]{x1D6FE}k^{2}+2\unicode[STIX]{x1D712}k^{4})\unicode[STIX]{x1D644}\right]\hat{\boldsymbol{v}}^{(1)}+\boldsymbol{m}\hat{p}^{(1)},\quad \text{where }\nonumber\\ \displaystyle & & \displaystyle \quad \unicode[STIX]{x1D648}=\frac{2\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FE}}{\unicode[STIX]{x1D719}}(\unicode[STIX]{x1D640}^{(0)}\boldsymbol{k})\otimes (\unicode[STIX]{x1D640}^{(0)}\boldsymbol{k})-\unicode[STIX]{x1D647}^{(0)}\text{ and }\boldsymbol{m}=-\frac{\imath }{\unicode[STIX]{x1D719}}\unicode[STIX]{x1D649}^{(0)}\boldsymbol{k},\nonumber\\ \displaystyle & & \displaystyle \hspace{60.00009pt}\text{with }k=|\boldsymbol{k}|=\sqrt{k_{i}k_{i}}\text{ and }\hat{d}_{t}^{(0)}=\unicode[STIX]{x2202}_{t}-(\unicode[STIX]{x1D647}^{(0)\text{T}}\boldsymbol{k})\boldsymbol{\cdot }\hat{\unicode[STIX]{x1D735}},\end{eqnarray}$$

where $\unicode[STIX]{x1D644}$ denotes the unit tensor.

It is further worth noting that the term in $\unicode[STIX]{x1D647}^{(0)}$ appearing in the expression for $\unicode[STIX]{x1D648}$ represents convective distortion of the Fourier mode $\hat{\boldsymbol{v}}$ by the base flow, whereas that appearing in the convected derivative, defined by the final term of (3.7), represents wavevector stretching (Thomson Reference Thomson1887; Goddard Reference Goddard2003). We shall show below that the latter term is crucial to the asymptotic stability for $t\rightarrow \infty$ .

Taking the scalar product with $\boldsymbol{k}$ of (3.7) and invoking incompressibility, $\boldsymbol{k}\boldsymbol{\cdot }\hat{\boldsymbol{v}}^{(1)}=0$ , and its consequence,

(3.8) $$\begin{eqnarray}\boldsymbol{k}\boldsymbol{\cdot }(\hat{d}_{t}^{(0)}\hat{\boldsymbol{v}}^{(1)})=\boldsymbol{k}\boldsymbol{\cdot }\unicode[STIX]{x1D647}^{(0)}\hat{\boldsymbol{v}}^{(1)},\end{eqnarray}$$

we obtain $\hat{p}^{(1)}$ as one component of the oblique projection $\boldsymbol{m}\otimes \boldsymbol{k}/\boldsymbol{m}\boldsymbol{\cdot }\boldsymbol{k}$ ,

(3.9) $$\begin{eqnarray}\hat{p}^{(1)}=\frac{1}{\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{m}}\boldsymbol{k}\boldsymbol{\cdot }(\unicode[STIX]{x1D647}^{(0)}-\unicode[STIX]{x1D648})\hat{\boldsymbol{v}}^{(1)}.\end{eqnarray}$$

Substitution of (3.9) into (3.7) then gives

(3.10) $$\begin{eqnarray}\displaystyle & & \displaystyle \hat{d}_{t}^{(0)}\hat{\boldsymbol{v}}^{(1)}=\unicode[STIX]{x1D63C}\hat{\boldsymbol{v}}^{(1)},\quad \text{where}\nonumber\\ \displaystyle & & \displaystyle \quad \unicode[STIX]{x1D63C}=\left(\unicode[STIX]{x1D644}-\frac{(\unicode[STIX]{x1D649}^{(0)}\boldsymbol{k})\otimes \boldsymbol{k}}{\boldsymbol{k}\boldsymbol{\cdot }\unicode[STIX]{x1D649}^{(0)}\boldsymbol{k}}\right)(\unicode[STIX]{x1D648}-\unicode[STIX]{x1D647}^{(0)})+\unicode[STIX]{x1D647}^{(0)}-\frac{1}{\unicode[STIX]{x1D719}}(\unicode[STIX]{x1D6FE}k^{2}+2\unicode[STIX]{x1D712}k^{4})\unicode[STIX]{x1D644}.\end{eqnarray}$$

The first term in round brackets on the right-hand side of (3.10) represents another oblique projection, orthogonal to $\boldsymbol{k}$ , and is independent of $k=|\boldsymbol{k}|$ . It should be noted that $\boldsymbol{k}\boldsymbol{\cdot }\unicode[STIX]{x1D649}^{(0)}\boldsymbol{k}=k^{2}-\unicode[STIX]{x1D6FC}(\boldsymbol{k}\boldsymbol{\cdot }\unicode[STIX]{x1D640}^{(0)}\boldsymbol{k})$ is a positive-definite quadratic form in $\boldsymbol{k}$ since $\unicode[STIX]{x1D6FC}<1$ and $|\unicode[STIX]{x1D640}^{(0)}|=1$ (cf. Ref. 1).

Transformation of (3.10) from the coordinate $\boldsymbol{k}$ to the dual material coordinate $\unicode[STIX]{x1D73F}$ gives the linear-stability equation as the canonical ordinary differential equation (ODE),

(3.11) $$\begin{eqnarray}\displaystyle & & \displaystyle \frac{\text{d}\hat{\unicode[STIX]{x1D74A}}}{\text{d}t}=\unicode[STIX]{x1D63C}(t,\unicode[STIX]{x1D73F})\hat{\unicode[STIX]{x1D74A}},\quad \text{with }\hat{\unicode[STIX]{x1D74A}}(t,\unicode[STIX]{x1D73F})=\hat{\boldsymbol{v}}^{(1)}(t,\boldsymbol{k}),\frac{\text{d}\hat{\unicode[STIX]{x1D74A}}}{\text{d}t}=\left(\frac{\unicode[STIX]{x2202}\hat{\unicode[STIX]{x1D74A}}}{\unicode[STIX]{x2202}t}\right)_{\unicode[STIX]{x1D73F}},\nonumber\\ \displaystyle & & \displaystyle \quad \text{with }\boldsymbol{k}=(\boldsymbol{F}^{(0)})^{-\text{T}}\unicode[STIX]{x1D73F},\text{and }d_{t}\boldsymbol{F}^{(0)}=\unicode[STIX]{x1D647}^{(0)}\boldsymbol{F}^{(0)},\boldsymbol{F}^{(0)}=\unicode[STIX]{x1D644}\;@\;t=0,\end{eqnarray}$$

where $\unicode[STIX]{x1D73F}$ represents the initial wavevector $\boldsymbol{k}$ at $t=0$ and the superscript $-T$ represents inverse transpose. This provides the linear-stability theory for general homogeneous shearing, similar to that treated elsewhere (Goddard Reference Goddard2003), where it is identified as a time-dependent (or ‘non-autonomous’) stability problem.

For time-dependent stability, we recall that the eigenvalues of $\unicode[STIX]{x1D63C}(t,\unicode[STIX]{x1D73F})$ in (3.11) serve mainly to determine local stability at time $t$ and fixed $\unicode[STIX]{x1D73F}$ , with logarithmic growth of spectral energy $|\hat{\unicode[STIX]{x1D74A}}|^{2}(\unicode[STIX]{x1D73F},t)=\hat{\unicode[STIX]{x1D74A}}^{\ast }\boldsymbol{\cdot }\hat{\unicode[STIX]{x1D74A}}$ given by the first equation of (3.11) as twice the Rayleigh quotient, namely

(3.12) $$\begin{eqnarray}\frac{\text{d}}{\text{d}t}\ln |\hat{\unicode[STIX]{x1D74A}}|^{2}=2\frac{\hat{\unicode[STIX]{x1D74A}}^{\ast }\boldsymbol{\cdot }\unicode[STIX]{x1D63C}\hat{\unicode[STIX]{x1D74A}}}{|\hat{\unicode[STIX]{x1D74A}}|^{2}},\end{eqnarray}$$

where $^{\ast }$ denotes complex conjugate. Whenever $\unicode[STIX]{x1D63C}$ has real eigenvalues, we may take $\hat{\unicode[STIX]{x1D74A}}^{\ast }=\hat{\unicode[STIX]{x1D74A}}$ , and the maximum of the right-hand side over all $\hat{\unicode[STIX]{x1D74A}}$ equals $2\unicode[STIX]{x1D706}(t)$ , where $\unicode[STIX]{x1D706}(t)$ is the largest eigenvalue. Otherwise, the complex eigenvalues of $\unicode[STIX]{x1D63C}$ represent stationary points in the complex plane (Didwania & Goddard Reference Didwania and Goddard1993). In the former case, it follows that the greatest eigenvalue $\unicode[STIX]{x1D706}(t)$ implies asymptotic stability in the sense of energy if $\unicode[STIX]{x1D706}(t)<0$ for $t\rightarrow \infty$ .

3.1 Stability of steady simple shear

Following Ref. 1, we consider a homogeneous steady simple shearing as the base state, with non-dimensional versions of the canonical forms given by $\boldsymbol{v}^{(0)}=(x_{2},0,0)$ ,

(3.13a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D647}^{(0)}=\left[\begin{array}{@{}ccc@{}}0 & 1 & 0\\ 0 & 0 & 0\\ 0 & 0 & 0\end{array}\right]\quad \text{and}\quad \unicode[STIX]{x1D640}^{(0)}=\frac{1}{\sqrt{2}}\left[\begin{array}{@{}ccc@{}}0 & 1 & 0\\ 1 & 0 & 0\\ 0 & 0 & 0\end{array}\right],\quad \text{with }|\unicode[STIX]{x1D63F}^{(0)}|=1/\surd 2.\end{eqnarray}$$

Hence, $(\boldsymbol{F}^{(0)})^{-\text{T}}=(\unicode[STIX]{x1D644}-\unicode[STIX]{x1D647}^{(0)\text{T}}t)$ , and the wavevector $\boldsymbol{k}$ is given in terms of the initial wavenumber $\unicode[STIX]{x1D73F}$ by (3.11) as follows. (See equation (40) in the remarkable paper of Thomson (Thomson Reference Thomson1887, §§ 32–39). His exact solution for the perturbation of a simple shear flow governed by the Navier–Stokes equations illustrates the shortcomings of stability analyses based on initial growth rates, as pointed out in several previous works, e.g. Alam & Nott (Reference Alam and Nott1997).)

(3.14) $$\begin{eqnarray}\boldsymbol{k}(\unicode[STIX]{x1D73F},t)=(k_{1},k_{2},k_{3})=(\unicode[STIX]{x1D705}_{1},\unicode[STIX]{x1D705}_{2}-\unicode[STIX]{x1D705}_{1}t,\unicode[STIX]{x1D705}_{3}),\end{eqnarray}$$

with wavevector stretching restricted to the component $k_{2}$ . Appendix B presents a Squire’s-type theorem showing that planar perturbations with $k_{3}=0$ are the least stable.

For planar perturbations, the components of the tensors appearing in (3.10) can be displayed in the 2D format,

(3.15) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \unicode[STIX]{x1D644}-\frac{\unicode[STIX]{x1D649}^{(0)}\boldsymbol{k}\otimes \boldsymbol{k}}{\boldsymbol{k}\boldsymbol{\cdot }\unicode[STIX]{x1D649}^{(0)}\boldsymbol{k}}=\frac{1}{\unicode[STIX]{x1D6F7}_{1}}\left[\begin{array}{@{}cc@{}}\displaystyle \unicode[STIX]{x1D6F7}_{1}+k_{1}\left(\frac{\unicode[STIX]{x1D6FC}}{\sqrt{2}}k_{2}-k_{1}\right) & \displaystyle k_{2}\left(\frac{\unicode[STIX]{x1D6FC}}{\sqrt{2}}k_{2}-k_{1}\right)\\ \displaystyle k_{1}\left(\frac{\unicode[STIX]{x1D6FC}}{\sqrt{2}}k_{1}-k_{2}\right) & \displaystyle \unicode[STIX]{x1D6F7}_{1}+k_{2}\left(\frac{\unicode[STIX]{x1D6FC}}{\sqrt{2}}k_{1}-k_{2}\right)\end{array}\right],\\ \displaystyle \unicode[STIX]{x1D648}-\unicode[STIX]{x1D647}^{(0)}=\frac{1}{\unicode[STIX]{x1D719}}\left[\begin{array}{@{}cc@{}}\displaystyle \unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FE}k_{2}^{2} & \displaystyle \unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FE}k_{1}k_{2}-2\unicode[STIX]{x1D719}\\ \displaystyle \unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FE}k_{1}k_{2} & \displaystyle \unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FE}k_{1}^{2}\end{array}\right],\end{array}\right\}\end{eqnarray}$$

where $\unicode[STIX]{x1D6F7}_{1}=\boldsymbol{k}\boldsymbol{\cdot }\unicode[STIX]{x1D649}^{(0)}\boldsymbol{k}=k^{2}-\sqrt{2}\unicode[STIX]{x1D6FC}k_{1}k_{2}$ and $\unicode[STIX]{x1D6F7}_{2}=\unicode[STIX]{x1D6FE}k^{2}+2\unicode[STIX]{x1D712}k^{4}$ . Hence,

(3.16) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \unicode[STIX]{x1D608}_{11}=\frac{1}{\unicode[STIX]{x1D719}\unicode[STIX]{x1D6F7}_{1}}[-\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FE}k_{2}^{2}(k_{1}^{2}-k_{2}^{2})-\unicode[STIX]{x1D6F7}_{1}\unicode[STIX]{x1D6F7}_{2}],\\ \displaystyle \unicode[STIX]{x1D608}_{12}=\frac{1}{\unicode[STIX]{x1D719}\unicode[STIX]{x1D6F7}_{1}}[-\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FE}k_{1}k_{2}+\unicode[STIX]{x1D719}](k_{1}^{2}-k_{2}^{2}),\\ \displaystyle \unicode[STIX]{x1D608}_{21}=\frac{1}{\unicode[STIX]{x1D719}\unicode[STIX]{x1D6F7}_{1}}[\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FE}k_{1}k_{2}(k_{1}^{2}-k_{2}^{2})],\\ \displaystyle \unicode[STIX]{x1D608}_{22}=\frac{1}{\unicode[STIX]{x1D719}\unicode[STIX]{x1D6F7}_{1}}[\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FE}k_{1}^{2}(k_{1}^{2}-k_{2}^{2})-\unicode[STIX]{x1D6F7}_{1}\unicode[STIX]{x1D6F7}_{2}+2\unicode[STIX]{x1D719}k_{1}(k_{2}-\unicode[STIX]{x1D6FC}k_{1}/\surd 2)],\end{array}\right\}\end{eqnarray}$$

where the components of $\unicode[STIX]{x1D63C}$ without convection are obtained by setting $\unicode[STIX]{x1D719}=0$ in the numerator of the expressions for $\unicode[STIX]{x1D608}_{12}$ and $\unicode[STIX]{x1D608}_{22}$ .

The determinantal condition $\det (\unicode[STIX]{x1D63C}-\unicode[STIX]{x1D706}\unicode[STIX]{x1D644})=0$ gives the eigenvalue $\unicode[STIX]{x1D706}$ of $\unicode[STIX]{x1D63C}$ with largest real part as

(3.17) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D706} & = & \displaystyle \frac{\unicode[STIX]{x1D608}_{11}+\unicode[STIX]{x1D608}_{22}}{2}+\left[\left(\frac{\unicode[STIX]{x1D608}_{11}-\unicode[STIX]{x1D608}_{22}}{2}\right)^{2}+\unicode[STIX]{x1D608}_{12}\unicode[STIX]{x1D608}_{21}\right]^{1/2}\nonumber\\ \displaystyle & = & \displaystyle \frac{1}{2\unicode[STIX]{x1D719}\unicode[STIX]{x1D6F7}_{1}}[\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FE}(k_{1}^{2}-k_{2}^{2})^{2}-2\unicode[STIX]{x1D6F7}_{1}\unicode[STIX]{x1D6F7}_{2}+2\unicode[STIX]{x1D719}k_{1}(k_{2}-\unicode[STIX]{x1D6FC}k_{1}/\surd 2)+\unicode[STIX]{x1D6F7}_{3}^{1/2}],\end{eqnarray}$$

where

(3.18) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6F7}_{3} & = & \displaystyle \unicode[STIX]{x1D6FD}^{2}\unicode[STIX]{x1D6FE}^{2}(k_{1}^{2}-k_{2}^{2})^{4}+2\unicode[STIX]{x1D719}^{2}k_{1}^{2}(\unicode[STIX]{x1D6FC}k_{1}-\surd 2k_{2})^{2}\nonumber\\ \displaystyle & & \displaystyle -\,2\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D719}k_{1}^{2}(k_{1}^{2}-k_{2}^{2})(\surd 2\unicode[STIX]{x1D6FC}k_{1}^{2}-4k_{1}k_{2}+\surd 2\unicode[STIX]{x1D6FC}k_{2}^{2}).\end{eqnarray}$$

Negative values of $\unicode[STIX]{x1D6F7}_{3}$ represent oscillatory behaviour, analogous to the ‘flutter’ instabilities exhibited by certain elasto-plastic models, the subject of a comprehensive review by Bigoni (Reference Bigoni1995).

When convection is neglected, we have

(3.19a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D6F7}_{3}=\unicode[STIX]{x1D6FD}^{2}\unicode[STIX]{x1D6FE}^{2}(k_{1}^{2}-k_{2}^{2})^{4}\quad \text{and}\quad \unicode[STIX]{x1D706}=\frac{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FE}(k_{1}^{2}-k_{2}^{2})^{2}-\unicode[STIX]{x1D6F7}_{1}\unicode[STIX]{x1D6F7}_{2}}{\unicode[STIX]{x1D719}\unicode[STIX]{x1D6F7}_{1}}.\end{eqnarray}$$

Since $\unicode[STIX]{x1D6F7}_{3}$ is non-negative in this case, it follows that there is no oscillatory behaviour without convection. Moreover, since $\unicode[STIX]{x1D6F7}_{3}$ does not depend on $\unicode[STIX]{x1D712}$ , the oscillatory frequencies arising from convection do not depend on $\unicode[STIX]{x1D712}$ . However, the amplitude of oscillation does depends on $\unicode[STIX]{x1D712}$ through the real part of $\unicode[STIX]{x1D706}$ , as discussed further below.

Figure 1.

Stability diagram for (a) without and (b) with convection, for $\unicode[STIX]{x1D712}=0$ , with $I=0.0001$ , $\unicode[STIX]{x1D719}=0.5$ and $p^{(0)}=1$ . The solid curves represent neutral stability.

These relations can be put in polar form, based on the polar representation

(3.20a,b ) $$\begin{eqnarray}k=|\boldsymbol{k}|=\sqrt{k_{1}^{2}+k_{2}^{2}}\quad \text{and}\quad \unicode[STIX]{x1D717}=\text{tan}^{-1}(\unicode[STIX]{x1D705}_{2}/\unicode[STIX]{x1D705}_{1}),\end{eqnarray}$$

where $k$ is the magnitude of the in-plane wavevector and $\unicode[STIX]{x1D717}$ is its angle with the direction of flow. With $c=\cos \unicode[STIX]{x1D717}$ and $s=\sin \unicode[STIX]{x1D717}$ , the eigenvalue from (3.17) can be expressed as

(3.21) $$\begin{eqnarray}\unicode[STIX]{x1D706}=\frac{1}{2\unicode[STIX]{x1D719}\tilde{\unicode[STIX]{x1D6F7}}_{1}}[\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FE}k^{2}(c^{2}-s^{2})^{2}-2\tilde{\unicode[STIX]{x1D6F7}}_{1}\unicode[STIX]{x1D6F7}_{2}+2\unicode[STIX]{x1D719}c(s-\unicode[STIX]{x1D6FC}c/\surd 2)+\tilde{\unicode[STIX]{x1D6F7}}_{3}^{1/2}],\end{eqnarray}$$

where

(3.22) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle \tilde{\unicode[STIX]{x1D6F7}}_{3}=\unicode[STIX]{x1D6F7}_{3}/k^{4}=\unicode[STIX]{x1D6FD}^{2}\unicode[STIX]{x1D6FE}^{2}(c^{2}-s^{2})^{4}k^{4}-2\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D719}c^{2}(c^{2}-s^{2})(\surd 2\unicode[STIX]{x1D6FC}-4cs)k^{2}\\ \displaystyle +\,2\unicode[STIX]{x1D719}^{2}c^{2}(\unicode[STIX]{x1D6FC}c-\surd 2s)^{2}\quad \text{and}\quad \tilde{\unicode[STIX]{x1D6F7}}_{1}=\unicode[STIX]{x1D6F7}_{1}/k^{2}=1-\surd 2\unicode[STIX]{x1D6FC}cs,\end{array}\right\} & & \displaystyle\end{eqnarray}$$

in which the trigonometric relations $c^{2}-s^{2}=\cos 2\unicode[STIX]{x1D717}$ and $2cs=\sin 2\unicode[STIX]{x1D717}$ also apply.

3.2 Transient instability

Here, we consider the modifications of the stability analysis of Ref. 1 by the inclusion of convection and the vdW–C–H gradient terms. Since that analysis is strictly valid only for the growth rates inferred from the initial state, the wavenumber $\boldsymbol{k}$ is to be interpreted here as its initial value $\unicode[STIX]{x1D73F}$ . We shall employ the term ‘transient instability’ (or ‘stability’) to denote positive (or negative) growth rates based on these initial values. The parameters employed in this study are the values proposed by Jop et al. (Reference Jop, Forterre and Pouliquen2005), namely $\unicode[STIX]{x1D707}_{0}=0.383$ , $\unicode[STIX]{x1D707}_{\infty }=0.643$ and $I_{\ast }=0.279$ , unless otherwise specified.

Some contours of the initial growth rate, represented by the real part of the eigenvalue $\unicode[STIX]{x1D706}$ from (3.18), are shown in the stability diagram of figure 1, where panels (a) and (b) represent stability with and without convection respectively ( $\unicode[STIX]{x1D647}^{(0)}\equiv \mathbf{0}$ in all terms except $I$ ), for $I=0.0001$ and $\unicode[STIX]{x1D712}=0$ . The solid curves represent neutral stability, $\unicode[STIX]{x1D706}=0$ , with lighter zones representing unstable regions. Figure 1(a) is identical to the result of Ref. 1, whereas figure 1(b) shows that convection causes strong distortion of the neutral stability curve, eliminating parts of two unstable branches. It should be noted that convection has a more pronounced effect in the regions of small $k$ , as it represents a contribution of order zero in $k$ to (3.21), compared with terms of higher order in $k$ arising from the dissipative stresses.

Figure 2 illustrates the effect of $\unicode[STIX]{x1D712}$ , where it is obvious that the vdW–C–H model provides a wavenumber cutoff, shrinking unstable zones accordingly. By choosing very large $\unicode[STIX]{x1D712}$ one can, not surprisingly, eliminate the unstable zones almost completely.

Figure 2.

Stability without convection for (a) $\unicode[STIX]{x1D712}=0$ and (b) $\unicode[STIX]{x1D712}=0.005$ , with $I=0.0001$ , $\unicode[STIX]{x1D719}=0.5$ and $p^{(0)}=1$ . The solid curves represent neutral stability.

Figure 3.

Stability with convection for increasing $\unicode[STIX]{x1D712}$ values: (a) $0$ , (b) $0.001$ , (c) $0.1$ and (d) $1$ , with $I=0.0001$ , $\unicode[STIX]{x1D719}=0.5$ and $p^{(0)}=1$ .

Figure 4.

Figure 3 in polar variables.

Figure 3 shows the combined effects of convection and wavenumber cutoff for various values of $\unicode[STIX]{x1D712}$ , with $I=0.0001$ . It is once again clear that the region of instability is greatly reduced by increasing $\unicode[STIX]{x1D712}$ . This is perhaps made clearer by the plots in terms of polar variables in figure 4. It should be noted that the growth rates in $0\leqslant \unicode[STIX]{x1D717}<\unicode[STIX]{x03C0}$ are also repeated in $\unicode[STIX]{x03C0}\leqslant \unicode[STIX]{x1D717}<2\unicode[STIX]{x03C0}$ when convection is included.

As discussed above, oscillatory behaviour with $\text{Im}(\unicode[STIX]{x1D706})\neq 0$ arises for negative $\unicode[STIX]{x1D6F7}_{3}$ , and, writing $X=k^{2}>0$ , we see from (3.21) that Re $(\unicode[STIX]{x1D706})$ and $\unicode[STIX]{x1D6F7}_{3}$ are given by separate quadratic polynomials in $X$ . Then, in order to have undamped oscillation, we must require the simultaneous conditions

(3.23) $$\begin{eqnarray}\unicode[STIX]{x1D6F7}_{3}<0\quad \text{and}\quad \text{Re}(\unicode[STIX]{x1D706})\geqslant 0,\quad \text{for }X>0,\end{eqnarray}$$

on these two quadratics. While we have not been able to provide a straightforward and rigorous algebraic proof of the impossibility of (3.23), a detailed numerical investigation for various values of $k$ at closely spaced increments of $\unicode[STIX]{x1D717}$ in $(0,2\unicode[STIX]{x03C0})$ and for various values of $\unicode[STIX]{x1D712}$ we failed to achieve the condition. Hence, we are led to conjecture that oscillatory solutions are always damped.

Figure 5 shows the corresponding imaginary part of the eigenvalue $\unicode[STIX]{x1D706}$ for a particular set of parameter values, where the finger-like regions bounded by the dashed lines contain the non-zero imaginary values of $\unicode[STIX]{x1D706}$ . Also shown are two contours of $\text{Re}(\unicode[STIX]{x1D706})$ , the neutral stability curve and a curve representing damping. In line with the above conjecture, the oscillatory behaviour is damped.

Following Ref. 1, we present in figure 6 a stability diagram in the $I$ – $\unicode[STIX]{x0394}\unicode[STIX]{x1D707}$ plane to show the effect of $\unicode[STIX]{x1D712}$ with and without convection, for $\boldsymbol{k}=(0.1,0.8)$ (or $k=0.806$ , $\unicode[STIX]{x1D717}=82.9^{\circ }$ ). The curves shown there are loci of neutral stability, and the regions inside the closed curves or beneath the open curves are the stable regions. The curve $\unicode[STIX]{x1D712}=0$ in figure 6(a) is a magnified version of that given in Ref. 1, indicating that the perfectly plastic case of constant $\unicode[STIX]{x1D707}$ , with $\unicode[STIX]{x0394}\unicode[STIX]{x1D707}=0$ , is always unstable. However, this instability is apparently removed at a critical $\unicode[STIX]{x1D712}$ value between 0.001 and 0.002, where the closed curve opens up to form an unbounded stable region. Interestingly, this critical $\unicode[STIX]{x1D712}$ value is nearly 10 times larger when convection is included, as illustrated by figure 6(b).

In summary, the preceding results indicate that convection distorts the regions of transient instability, completely eliminating instability for some directions of $\boldsymbol{k}$ while leading to short-wavelength instability for others, in the absence of cutoff by the gradient terms in $\unicode[STIX]{x1D712}$ .

Figure 5.

Oscillatory behaviour resulting from convection, for $\unicode[STIX]{x1D712}=0$ , $I=0.0001$ , $p^{(0)}=1$ and $\unicode[STIX]{x1D719}=0.5$ . The oscillation frequency is non-zero within the finger-like regions bounded by dashed curves and equal to zero upon and outside these curves. The solid and dotted curves are contours of the real part of $\unicode[STIX]{x1D706}$ , the former representing neutral stability.

Figure 6.

Stability in the $I$ – $\unicode[STIX]{x0394}\unicode[STIX]{x1D707}$ plane, (a) without and (b) with convection, showing the effect of $\unicode[STIX]{x1D712}$ , where $\boldsymbol{k}=(0.1,0.8)$ ( $k=0.806,\unicode[STIX]{x1D717}=82.9^{\circ }$ ),  $p^{(0)}=1$ and $\unicode[STIX]{x1D719}=0.5$ .

3.3 Asymptotic behaviour

According to (3.14), for large $t$ and $\unicode[STIX]{x1D705}_{1}\neq 0$ , we have $|k_{2}|\gg |k_{1}|$ , resulting in the following asymptotic forms:

(3.24) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \unicode[STIX]{x1D608}_{11}\sim k_{2}^{2}[(\unicode[STIX]{x1D6FD}-1)\unicode[STIX]{x1D6FE}-2\unicode[STIX]{x1D712}k_{2}^{2}]/\unicode[STIX]{x1D719},\quad \unicode[STIX]{x1D608}_{12}\sim -\unicode[STIX]{x1D608}_{12}\sim \unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FE}k_{1}k_{2}/\unicode[STIX]{x1D719},\\ \displaystyle \unicode[STIX]{x1D608}_{22}\sim -k_{2}^{2}[\unicode[STIX]{x1D6FE}+2\unicode[STIX]{x1D712}k_{2}^{2}]/\unicode[STIX]{x1D719}\quad \text{and}\quad \unicode[STIX]{x1D706}\sim \unicode[STIX]{x1D608}_{11},\end{array}\right\}\end{eqnarray}$$

up to terms of relative magnitude $O(|k_{2}|^{-1})$ . Since $\unicode[STIX]{x1D6FD}-1=-\stackrel{\hspace{2.22198pt}\prime }{\unicode[STIX]{x1D707}}\leqslant 0$ , for $\unicode[STIX]{x1D707}_{\infty }\geqslant \unicode[STIX]{x1D707}_{0},$ it follows from (3.24) that the local growth rate for $\unicode[STIX]{x1D712}\geqslant 0$ is never positive for large $t$ whenever $\unicode[STIX]{x1D705}_{1}\neq 0$ . In the exceptional case $\unicode[STIX]{x1D705}_{1}=0$ , where $\boldsymbol{k}(t)\equiv \unicode[STIX]{x1D73F}=(0,\unicode[STIX]{x1D705}_{2})$ , it is easy to show by means of (3.16)–(3.18) that $\unicode[STIX]{x1D706}$ is never positive for $\unicode[STIX]{x1D712}\geqslant 0$ . Therefore, we conclude that any transient instability is eventually killed off by wavevector stretching, irrespective of vdW–C–H regularization. The resulting asymptotic stability is illustrated by a numerical solution of (3.11). For the incompressible planar flows considered here, this is facilitated by the following closed-form solution for the stream function.

3.3.1 Scalar solution to (3.11)

For planar flow, the condition of incompressibility can be expressed in the usual way in terms of a stream function $\unicode[STIX]{x1D713}(x_{1},x_{2})$ as

(3.25a,b ) $$\begin{eqnarray}v_{1}=\unicode[STIX]{x2202}_{x_{2}}\unicode[STIX]{x1D713},\quad v_{2}=-\unicode[STIX]{x2202}_{x_{1}}\unicode[STIX]{x1D713},\end{eqnarray}$$

which in terms of Fourier transforms becomes

(3.26) $$\begin{eqnarray}\displaystyle & & \displaystyle \hat{v}_{1}=\imath k_{2}\hat{\unicode[STIX]{x1D713}},\quad \hat{v}_{2}=-\imath k_{1}\hat{\unicode[STIX]{x1D713}}\quad \text{or}\quad \hat{\boldsymbol{v}}=\imath \boldsymbol{k}^{\bot }\hat{\unicode[STIX]{x1D713}},\nonumber\\ \displaystyle & & \displaystyle \quad \text{where }\boldsymbol{k}^{\bot }=(k_{2},-k_{1})=\unicode[STIX]{x1D64C}\boldsymbol{k},\unicode[STIX]{x1D64C}=\left[\begin{array}{@{}cc@{}}0 & 1\\ -1 & 0\end{array}\right],\end{eqnarray}$$

$\unicode[STIX]{x1D64C}$ is orthogonal and $\boldsymbol{k}^{\bot }$ is identical to the flipped wavevector introduced in Ref. 1. It is easy to see that the preceding linear relations also apply to perturbations. Hence, after making use of $\text{d}\boldsymbol{k}/\text{d}t=-\unicode[STIX]{x1D647}^{(0)\text{T}}\boldsymbol{k}$ , the ODE can be reduced to the form

(3.27) $$\begin{eqnarray}\boldsymbol{k}\frac{\text{d}\hat{\unicode[STIX]{x1D713}}}{\text{d}t}=\unicode[STIX]{x1D63D}\boldsymbol{k}\hat{\unicode[STIX]{x1D713}},\quad \text{where }\unicode[STIX]{x1D63D}=\unicode[STIX]{x1D64C}^{\text{T}}\unicode[STIX]{x1D63C}\unicode[STIX]{x1D64C}+\unicode[STIX]{x1D647}^{(0)\text{T}}.\end{eqnarray}$$

Now, $\boldsymbol{k}$ and $\boldsymbol{k}^{\bot }$ serve as orthogonal bases for 2D vectors and it is easy to show that $\boldsymbol{k}^{\bot }\boldsymbol{\cdot }\unicode[STIX]{x1D63D}\boldsymbol{k}=0$ , i.e. the right-hand side of (3.27) has zero component in direction $\boldsymbol{k}^{\bot }$ . By Fredholm’s theorem, applied implicitly in Ref. 1, this establishes that $\boldsymbol{k}$ is in the range of $\unicode[STIX]{x1D63D}$ and, hence, is an eigenvector of $\unicode[STIX]{x1D63D}$ . Thus, (3.27) involves only components in the direction $\boldsymbol{k}$ , given by orthogonal projection as

(3.28) $$\begin{eqnarray}\frac{\text{d}\hat{\unicode[STIX]{x1D713}}}{\text{d}t}=\unicode[STIX]{x1D6EC}(t,\unicode[STIX]{x1D73F})\hat{\unicode[STIX]{x1D713}},\quad \text{where }\unicode[STIX]{x1D6EC}=\frac{\boldsymbol{k}\boldsymbol{\cdot }\unicode[STIX]{x1D63D}\boldsymbol{k}}{k^{2}}.\end{eqnarray}$$

Here,

(3.29) $$\begin{eqnarray}\unicode[STIX]{x1D6EC}=\frac{1}{\unicode[STIX]{x1D719}k^{2}\unicode[STIX]{x1D6F7}_{1}}\{\unicode[STIX]{x1D719}[k_{1}k_{2}(k^{2}+\unicode[STIX]{x1D6F7}_{1})-\surd 2\unicode[STIX]{x1D6FC}k_{1}^{4}]+\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FE}k^{2}(k_{1}^{2}-k_{2}^{2})^{2}-\unicode[STIX]{x1D6F7}_{1}\unicode[STIX]{x1D6F7}_{2}k^{2}\}\end{eqnarray}$$

is the eigenvalue for eigenvector $\boldsymbol{k}$ of $\unicode[STIX]{x1D63D}$ .

It follows that the solution to the linear-stability problem for constant $\unicode[STIX]{x1D647}^{(0)}$ is given by

(3.30) $$\begin{eqnarray}\hat{\boldsymbol{v}}^{(1)}=\hat{\unicode[STIX]{x1D713}}_{0}\unicode[STIX]{x1D64C}\exp \left\{\int _{0}^{t}\unicode[STIX]{x1D6EC}(t^{\prime },\unicode[STIX]{x1D73F})\,\text{d}t^{\prime }\right\}\boldsymbol{k}=\hat{\unicode[STIX]{x1D713}}_{0}\unicode[STIX]{x1D64C}\exp \left\{\unicode[STIX]{x1D644}\int _{0}^{t}\unicode[STIX]{x1D6EC}(t^{\prime },\unicode[STIX]{x1D73F})\,\text{d}t^{\prime }-\unicode[STIX]{x1D647}^{(0)\text{T}}t\right\}\unicode[STIX]{x1D73F},\end{eqnarray}$$

which gives

(3.31) $$\begin{eqnarray}|\hat{\boldsymbol{v}}^{(1)}|=|\hat{\unicode[STIX]{x1D713}}_{0}|k\exp \left\{\int _{0}^{t}\unicode[STIX]{x1D6EC}(t^{\prime },\unicode[STIX]{x1D73F})\,\text{d}t^{\prime }\right\},\end{eqnarray}$$

where, in all of the above formulae,

(3.32) $$\begin{eqnarray}\boldsymbol{k}(t,\unicode[STIX]{x1D73F})=\exp \{-\unicode[STIX]{x1D647}^{(0)\text{T}}t\}\unicode[STIX]{x1D73F},\quad \text{with }k=[\unicode[STIX]{x1D73F}\boldsymbol{\cdot }\exp \{-\unicode[STIX]{x1D647}^{(0)}t\}\exp \{-\unicode[STIX]{x1D647}^{(0)\text{T}}t\}\unicode[STIX]{x1D73F}]^{1/2}.\end{eqnarray}$$

It is easy to show that (3.31) gives the same result as that obtained by substitution of $\hat{\boldsymbol{v}}^{(1)}=\unicode[STIX]{x1D713}\boldsymbol{k}^{\bot }$ into (3.12). Moreover, when convection is neglected, one finds that $\boldsymbol{k}\boldsymbol{\cdot }\unicode[STIX]{x1D63D}\boldsymbol{k}=\boldsymbol{k}^{\bot }\boldsymbol{\cdot }\unicode[STIX]{x1D63C}\,\boldsymbol{k}^{\bot }$ , and, hence, that $\boldsymbol{k}^{\bot }$ is an eigenvector of $\unicode[STIX]{x1D63C}$ with eigenvalue $\unicode[STIX]{x1D6EC}$ , as already found in Ref. 1.

Making use of the asymptotic form (3.24), we find for large $k$ with $|k_{2}|\gg |k_{1}|$ that

(3.33) $$\begin{eqnarray}\unicode[STIX]{x1D6EC}\sim \frac{k_{2}^{2}}{\unicode[STIX]{x1D719}\tilde{\unicode[STIX]{x1D6F7}}_{1}}[(\unicode[STIX]{x1D6FD}-\tilde{\unicode[STIX]{x1D6F7}}_{1})\unicode[STIX]{x1D6FE}-2\tilde{\unicode[STIX]{x1D6F7}}_{1}\unicode[STIX]{x1D712}k_{2}^{2}]\sim \unicode[STIX]{x1D706},\quad \text{with }\tilde{\unicode[STIX]{x1D6F7}}_{1}=\unicode[STIX]{x1D6F7}_{1}/k^{2}\sim \unicode[STIX]{x1D6F7}_{1}/k_{2}^{2},\text{for }t\rightarrow \infty ,\end{eqnarray}$$

where $\unicode[STIX]{x1D706}$ is the eigenvalue of $\unicode[STIX]{x1D63C}$ with largest real part, which confirms the asymptotic behaviour inferred previously from (3.24).

Figure 7.

The magnitude $|\hat{v}(\unicode[STIX]{x1D73F},t)|$ of various Fourier modes versus $t$ for $\unicode[STIX]{x1D712}=10^{-6}$ and initial conditions (a) $\unicode[STIX]{x1D73F}=(-0.5,0.866)$ , (b) $\unicode[STIX]{x1D73F}=(-0.707,-0.707)$ , (c) $\unicode[STIX]{x1D73F}=(-3.54,-3.54)$ and (d) $\unicode[STIX]{x1D73F}=(0.259,0.966)$ .

Figure 7 presents semi-log plots of (3.31), with $|\hat{\unicode[STIX]{x1D713}}_{0}|=1/\unicode[STIX]{x1D705}$ and, hence, $|\hat{\boldsymbol{v}}^{(1)}|=1$ at $t=0$ , as obtained from the integration of (3.28) by means of the finite-difference approximation (FDA) and the MATLAB^® numerical integrator ‘ode45’ for two values $\unicode[STIX]{x1D705}=1$ and $\unicode[STIX]{x1D705}=5$ , and parameter values

(3.34a-f ) $$\begin{eqnarray}I^{(0)}=0.001,\quad I_{\ast }=0.2790,\quad \unicode[STIX]{x1D707}_{0}=0.3830,\quad \unicode[STIX]{x1D707}_{\infty }=0.6430,\quad \unicode[STIX]{x1D719}=0.5,\quad p^{(0)}=1,\end{eqnarray}$$

with cessation of integration for values of $|\hat{\unicode[STIX]{x1D74A}}|$ less than 0.001. Figures 7(b) and 7(c) serve to illustrate the effect of $\unicode[STIX]{x1D705}=|\unicode[STIX]{x1D73F}|$ for $\unicode[STIX]{x1D717}=45^{\circ }$ . These figures are not changed drastically by taking $\unicode[STIX]{x1D712}=1$ , whereas the transient instability in figure 7(d) is completely eliminated.

As pointed out above in the paragraph following (3.5), material instability can also be characterized as the loss of generalized ellipticity in the quasi-static (inertialess) field equations. The latter can be obtained by omitting all inertia terms from (3.11) and (3.27), which gives the following expression for $\unicode[STIX]{x1D6EC}$ :

(3.35) $$\begin{eqnarray}\unicode[STIX]{x1D6EC}=\frac{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FE}(k_{1}^{2}-k_{2}^{2})^{2}}{\unicode[STIX]{x1D6F7}_{1}(\boldsymbol{k})}-(\unicode[STIX]{x1D6FE}k^{2}+2\unicode[STIX]{x1D712}k^{4})=-\frac{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FE}(\unicode[STIX]{x1D709}_{1}^{2}-\unicode[STIX]{x1D709}_{2}^{2})^{2}}{\unicode[STIX]{x1D6F7}_{1}(\unicode[STIX]{x1D743})}+\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D709}^{2}-2\unicode[STIX]{x1D712}\unicode[STIX]{x1D709}^{4},\end{eqnarray}$$

where $\unicode[STIX]{x1D743}=\imath \boldsymbol{k}$ , the so-called symbol, is the algebraic representation of the differential operator $\unicode[STIX]{x1D735}$ (Renardy & Rogers Reference Renardy and Rogers2006). Hence, the positive-definite quadratic form $\unicode[STIX]{x1D6F7}_{1}(\unicode[STIX]{x1D743})$ in real-valued $\unicode[STIX]{x1D743}$ represents an anisotropic Laplacian. The latter can be reduced to the usual Laplacian on an appropriate coordinate system, with $\unicode[STIX]{x1D6F7}_{1}^{-1}$ representing the corresponding Green’s function.

Now, generalized ellipticity may be defined as the requirement that the differential operator represented by $-\unicode[STIX]{x1D6EC}$ be positive definite in real-valued $\unicode[STIX]{x1D709}$ (Browder Reference Browder1961; Brezis & Browder Reference Brezis and Browder1998; Renardy & Rogers Reference Renardy and Rogers2006), which in the present context requires that $\unicode[STIX]{x1D712}>0$ . Therefore, the vdW–C–H regularization is essential for material stability according to this criterion, and it seems worthwhile to illustrate how this regularization might serve to impart a diffuse length scale to shear bands.

4 Model of a steady shear band and its nonlinear stability

Based on the preceding analysis, with large $k_{2}$ representing dominant gradients in the direction $y=x_{2}$ , we assume a single stabilized shear band localized at $y=0$ in an initially unperturbed shear flow with $\boldsymbol{g}=0,\unicode[STIX]{x1D735}p=0,$ and $\boldsymbol{v}=(y,0)$ will take the form of a stable fully developed flow in the $x=x_{1}$ direction with $\boldsymbol{v}=(u(y),0)$ and with the partial derivatives $(\unicode[STIX]{x2202}_{t},\unicode[STIX]{x2202}_{x})$ of all quantities vanishing. Hence, after one integration with respect to $y$ , the momentum balance (2.11) reduces to the ODE

(4.1) $$\begin{eqnarray}\displaystyle u^{\prime \prime \prime } & = & \displaystyle f(u^{\prime })=\frac{\unicode[STIX]{x1D707}(I)-\unicode[STIX]{x1D707}(I^{(0)})}{2\surd 2\unicode[STIX]{x1D712}I^{(0)2}}\,\text{sgn}(u^{\prime })\nonumber\\ \displaystyle & = & \displaystyle \frac{(\unicode[STIX]{x1D707}_{\infty }-\unicode[STIX]{x1D707}_{0})I_{\ast }(|u^{\prime }|-1)\,\text{sgn}(u^{\prime })}{2\surd 2\unicode[STIX]{x1D712}I^{(0)}(I^{(0)}+I_{\ast })(I^{(0)}|u^{\prime }|+I_{\ast })},\nonumber\\ \displaystyle & & \displaystyle \quad \text{with }I=I^{(0)}|u^{\prime }|\text{and }I=I^{(0)}|u^{\prime }|\text{and }I^{(0)}=1/\surd p,\end{eqnarray}$$

where primes denote differentiation with respect to $y$ and we employ the non-dimensional variables (3.4). We further invoke the condition $u^{\prime }\rightarrow 1$ for $|y|\rightarrow \infty$ , corresponding to the unperturbed flow. Hence, assuming that $u(y)$ is an odd function of $y$ , we may focus on the half-space $y>0$ with further conditions $u(0)=0$ .

As discussed below, the form of the ODE (4.1) indicates that the limit $\unicode[STIX]{x1D712}\rightarrow 0$ leads to a well-known singular perturbation for small $\unicode[STIX]{x1D712}$ , in which one may generally neglect the vdW–C–H terms in the ‘outer region’ lying outside a thin shear band of thickness $O(\unicode[STIX]{x1D712}^{1/2})$ , which we recall represents the ratio of microscopic to macroscopic length scales.

To pursue an analytical treatment, we note that the function $f(u^{\prime })$ in (4.1) can be written as the derivative $\text{d}\tilde{\unicode[STIX]{x1D713}}_{0}/\text{d}u^{\prime }$ of a modified form of the dissipation potential $\unicode[STIX]{x1D713}_{0}$ ,

(4.2) $$\begin{eqnarray}\tilde{\unicode[STIX]{x1D713}}_{0}(u^{\prime })=\frac{(\unicode[STIX]{x1D707}_{\infty }-\unicode[STIX]{x1D707}_{0})I_{\ast }}{2\surd 2\unicode[STIX]{x1D712}I^{(0)3}}\left[\frac{|u^{\prime }|}{1+I_{\ast }/I^{(0)}}-\ln \left(\frac{I^{(0)}|u^{\prime }|}{I_{\ast }}+1\right)\right].\end{eqnarray}$$

Then, by a standard method, (4.1) can be integrated twice to yield $u^{\prime }$ as an implicit function of $y$ ,

(4.3) $$\begin{eqnarray}y=Y(u^{\prime },u_{0}^{\prime })=\int _{u^{\prime }}^{u_{0}^{\prime }}\frac{\text{d}w}{\sqrt{2[\tilde{\unicode[STIX]{x1D713}}_{0}(w)-\tilde{\unicode[STIX]{x1D713}}_{0}(1)]}}.\end{eqnarray}$$

It should be noted that (4.1)–(4.2) imply that $\tilde{\unicode[STIX]{x1D713}}_{0}(w)-\tilde{\unicode[STIX]{x1D713}}_{0}(1)\rightarrow 0$ quadratically and $y\rightarrow \infty$ logarithmically in $|w-1|$ as $w\rightarrow 1$ .

It is an easy matter to convert the preceding relation to the second implicit form

(4.4) $$\begin{eqnarray}\displaystyle & & \displaystyle u=y+\int _{u^{\prime }}^{u_{0}^{\prime }}\frac{(w-1)\,\text{d}w}{\sqrt{2[\tilde{\unicode[STIX]{x1D713}}_{0}(w)-\tilde{\unicode[STIX]{x1D713}}_{0}(1)]}}\sim y+\unicode[STIX]{x1D6E5}_{0},\quad \text{for }|u^{\prime }|\rightarrow 1,\nonumber\\ \displaystyle & & \displaystyle \hspace{16.00008pt}\text{where }\unicode[STIX]{x1D6E5}_{0}=\int _{1}^{u_{0}^{\prime }}\frac{(w-1)\,\text{d}w}{\sqrt{2[\tilde{\unicode[STIX]{x1D713}}_{0}(w)-\tilde{\unicode[STIX]{x1D713}}_{0}(1)]}},\end{eqnarray}$$

which gives $u^{\prime }$ implicitly as a function of $u$ and $y$ , and shows that $u_{0}^{\prime }=u^{\prime }(0)$ is a free parameter that must be specified to complete the solution. The quantity $\unicode[STIX]{x1D6E5}_{0}$ , the intercept at $y=0$ of the asymptotic solution for $y\rightarrow \infty$ , represents one-half of the apparent slip on a shear band of zero thickness at $y=0$ as seen in the far field. It can be calculated by numerical quadrature for given $u_{0}^{\prime }$ and provides quite a useful criterion for convergence of the numerical solution considered next. The quadrature was performed in the present study by means of the MATLAB^® function ‘integral.m’.

It should be noted that under the rescaling

(4.5a-d ) $$\begin{eqnarray}\bar{y}=y/I^{(0)}\surd \unicode[STIX]{x1D712},\quad \bar{u}=u/I^{(0)}\surd \unicode[STIX]{x1D712},\quad \bar{\unicode[STIX]{x1D713}}=I^{(0)2}\unicode[STIX]{x1D712}\tilde{\unicode[STIX]{x1D713}}_{0},\quad \bar{\unicode[STIX]{x1D6E5}}_{0}=\unicode[STIX]{x1D6E5}_{0}/I^{(0)}\surd \unicode[STIX]{x1D712}\end{eqnarray}$$

(where overbars are not to be confused with prior usage to denote the present non-dimensional variables), $u^{\prime }$ is invariant and the relation (4.4) gives, upon division by $I^{(0)}\surd \unicode[STIX]{x1D712}$ , a description of the inner layer discussed above in terms of $\bar{y}$ , $\bar{u}$ , $\bar{\unicode[STIX]{x1D713}}$ and $\bar{\unicode[STIX]{x1D6E5}}_{0}$ .

Since numerical methods are generally necessary to treat (4.3) or (4.4), it is more efficacious to integrate (4.1) numerically. Upon replacement of $u$ and $y$ by the scaled variables $\bar{u}$ and $\bar{y}$ defined in (4.5), the ODE (4.1) maintains the same form, with

(4.6) $$\begin{eqnarray}\unicode[STIX]{x1D712}\rightarrow 1,\quad I^{(0)}\rightarrow 1,\quad I_{\ast }\rightarrow I_{\ast }/I^{(0)}.\end{eqnarray}$$

Moreover, one can show that the replacement of $\unicode[STIX]{x1D712}$ by $\unicode[STIX]{x1D712}\unicode[STIX]{x1D707}_{0}$ in (4.5) converts the parameter $\unicode[STIX]{x0394}\unicode[STIX]{x1D707}=\unicode[STIX]{x1D707}_{\infty }-\unicode[STIX]{x1D707}_{0}$ to $\unicode[STIX]{x0394}\unicode[STIX]{x1D707}/\unicode[STIX]{x1D707}_{0}$ , with $\unicode[STIX]{x1D707}_{0}\rightarrow 1$ elsewhere. Since the same transformation applies to the full momentum balance, it serves to justify the general validity of the stability plot of Barker et al. (Reference Barker, Schaeffer, Bohorquez and Gray2015), represented by figure 6 in the present paper. We do not adopt this additional scaling in the present paper, since $\unicode[STIX]{x1D707}_{0}$ is not varied.

As pointed out above, the quantity $u_{0}^{\prime }=u^{\prime }(0)$ is an unknown which must in principle be specified by the complete solution to the full momentum balance. We recall that such a numerical solution was carried out in Ref. 1 for a non-homogeneous shear flow for $\unicode[STIX]{x1D712}=0$ . While the authors found instability, it is not clear what role nonlinear effects may play and, in line with comments in the introduction, whether their numerics serve to rule out asymptotically stable states with sharp shear bands.

It would take us well beyond the scope of the present work to pursue numerical solutions of the full momentum balance. Instead, we provide below a one-dimensional analysis of the nonlinear stability of steady shear bands against normal sinusoidal perturbations as a function of $u^{\prime }(0)$ . For illustrative purposes, we consider here the assumption that, up to factors of order unity, $u^{\prime }(0)\sim k_{2,max}$ , where $k_{2,max}$ is the component $k_{2}$ of the wavenumber in the vicinity of the maximum growth rate according to the linear-stability analysis. As indicated in the preceding paragraphs, a reasonable approximation for the transient growth rate is provided by the eigenvalue $\unicode[STIX]{x1D706}$ given by (3.17), whenever real. Hence, we have made use of the MATLAB^® program ‘fminunc’, with the quasi-Newton option, to determine numerically the unconstrained minimum in $-\unicode[STIX]{x1D706}(\boldsymbol{k})$ . Thus, for parameter values $\unicode[STIX]{x1D707}_{0}=0.383,\unicode[STIX]{x1D707}_{\infty }=0.643$ , $I_{\ast }=0.279$ , $I^{(0)}=0.001$ , $\unicode[STIX]{x1D712}=10^{-6}$ , we find $\boldsymbol{k}_{max}\approx \pm (2,28)$ .

Figure 8 shows the corresponding curve of $u/I^{(0)}\surd \unicode[STIX]{x1D712}$ versus $y/I^{(0)}\surd \unicode[STIX]{x1D712}$ obtained for $u^{\prime }(0)=28$ by means of the MATLAB^® numerical integrator ‘bvp4c.m’, a finite-difference code applicable to two-point boundary-value problems. Also shown as a dashed line in figure 8 is the asymptote for large $y/I^{(0)}\surd \unicode[STIX]{x1D712}$ , whose intercept at $y=0$ represents the quantity $\unicode[STIX]{x1D6E5}_{0}$ defined in (4.4). We note that the integration is sensitive to the choice of the initial FDA mesh, as defined by the number of mesh points $N$ and the FDA step-size $\unicode[STIX]{x0394}y$ , and the comparison of the intercept calculated from $u(y_{f})-y_{f}$ , where $y_{f}=N\unicode[STIX]{x0394}y$ , against the exact value $\unicode[STIX]{x1D6E5}_{0}$ calculated numerically from the last equation in (4.4).

Figure 8.

Velocity profile in a model steady-state shear band with $u_{0}^{\prime }=28$ , $I^{(0)}/I_{\ast }=0.036$ , $\unicode[STIX]{x1D707}_{0}=0.383$ and $\unicode[STIX]{x1D707}_{\infty }=0.643$ .

We note that in the case of a constant friction coefficient $\unicode[STIX]{x1D707}_{0}=\unicode[STIX]{x1D707}_{\infty }$ the one-dimensional form of (2.11) becomes

(4.7) $$\begin{eqnarray}2\unicode[STIX]{x1D712}u^{\prime \prime \prime \prime }-\frac{\unicode[STIX]{x1D707}_{0}p}{\surd 2}(\text{sgn}(u^{\prime }))^{\prime }=0,\end{eqnarray}$$

which formally involves a Dirac delta at points of transition from $u^{\prime }=0$ to $u^{\prime }\neq 0$ . It is clear that the term in $\unicode[STIX]{x1D712}$ represents a singular surface traction which balances the jump in frictional stress at such points. In any case, (4.7) can be integrated immediately to give the odd solution in $y$ with representation appropriate to the infinite half-space $y>0$ ,

(4.8) $$\begin{eqnarray}u=\left\{\begin{array}{@{}ll@{}}c_{0}, & \text{for }y>\unicode[STIX]{x1D6FF},\\ p\unicode[STIX]{x1D707}_{0}y^{3}/12\surd 2\unicode[STIX]{x1D712}+u_{0}^{\prime }y, & \text{for }0\leqslant y<\unicode[STIX]{x1D6FF}.\end{array}\right.\end{eqnarray}$$

The continuity condition $u=c_{0}$ at $y=\unicode[STIX]{x1D6FF}$ provides one condition on the three quantities $c_{0},u_{0}^{\prime }$ and $\unicode[STIX]{x1D6FF}$ . For $u_{0}^{\prime }>0$ , it follows that there is a discontinuity $p\unicode[STIX]{x1D707}_{0}\unicode[STIX]{x1D6FF}^{2}/4\surd 2\unicode[STIX]{x1D712}+u_{0}^{\prime }$ in $u^{\prime }$ at $y=\unicode[STIX]{x1D6FF}$ , leading to the Dirac delta anticipated above in (4.7). At any rate, upon specification of any two quantities, the solution consists of a rigid unyielded region $y>\unicode[STIX]{x1D6FF}$ riding on a shear band in $0\leqslant y\leqslant \unicode[STIX]{x1D6FF}$ . It appears that other solutions given by piecewise cubics in $y$ are possible for finite regions $0\leqslant y\leqslant L$ .

Based on the preceding analysis, one can envisage solutions with any number of diffuse shear bands interspersed with more gently sheared regions, which not only reflects multiplicity of solution but is no doubt related to the mesh-size dependence in various numerical treatments of the underlying field equations (Belytschko et al. Reference Belytschko, Chiang and Plaskacz1994).

While the above steady shear bands represent one possible solution of the steady field equations, it remains to show that they represent stable points in the space of all steady solutions.

4.1 Nonlinear stability of shear bands against parallel shearing

Without attempting a full two-dimensional analysis, we consider the nonlinear stability governed by the one-dimensional form of (2.11) for $u(t,y)$ ,

(4.9) $$\begin{eqnarray}\unicode[STIX]{x2202}_{t}u=\unicode[STIX]{x2202}_{y}\left[h(\unicode[STIX]{x2202}_{y}u)-2\unicode[STIX]{x2202}_{y}^{3}u\right],\quad \text{with }h(\unicode[STIX]{x2202}_{y}u)=\frac{1}{\surd 2}\unicode[STIX]{x1D707}(I)\,\text{sgn}(\unicode[STIX]{x2202}_{y}u)\text{ and }I=I^{(0)}|\unicode[STIX]{x2202}_{y}u|,\end{eqnarray}$$

where $u,y$ denote the variables $\bar{u},\bar{y}$ defined in (4.5) and $t$ denotes the scaled time

(4.10) $$\begin{eqnarray}\bar{t}=t/\unicode[STIX]{x1D719}\unicode[STIX]{x1D610}^{(0)4}\unicode[STIX]{x1D712}.\end{eqnarray}$$

It is obvious that this scaling of variables restricts all parameter dependence to the function $\unicode[STIX]{x1D707}(I)$ and that the partial differential equation (4.9) reduces to the rescaled ODE (4.1) at steady state.

To investigate the stability of the steady-state velocity $u=u_{s}(y)$ , we consider the special case of an initial condition on (4.9) given by a perturbation of $u_{s}$ . After considering several possibilities, we have settled on what appeared to be a particularly unstable perturbation represented by the sinusoidal form

(4.11) $$\begin{eqnarray}u(0,y)=u_{0}(y)=u_{s}(y)+A\sin k_{2}y,\end{eqnarray}$$

with constant amplitude $A$ and spatial frequency $k_{2}$ . We employ a standard FDA with spatial discretization on nodal points $y_{i},i=1,\ldots ,N$ , with $y_{1}=0$ and with representative spacing $\unicode[STIX]{x0394}y_{i}$ . We then employ the method of lines (MOL) (Schiesser Reference Schiesser2012) to solve the ODEs resulting from the discretization of (4.9) numerically by means of the MATLAB^® stiff integrator ode15s.m. The details are summarized below in appendix C, and, as pointed out there, we employ a certain number $M$ of invariant ‘ghost nodes’ at each end of the $y$ interval at which $u(t,y_{i})=u_{s}(y_{i})$ .

Figure 9.

Stable shear band with $k_{2}=27.8$ for $u_{0}^{\prime }=28.0$ , $I^{(0)}/I_{\ast }=0.036$ , $\unicode[STIX]{x1D707}_{0}=0.3830$ and $\unicode[STIX]{x1D707}_{\infty }=0.6430$ . The solid curves represent the initial steady-state and sinusoidally perturbed shear band with amplitude $A=300$ , while the points ● represent the final shear band.

Figure 10.

Maximally stable $u_{0}^{\prime }$ versus $k_{2}$ for 72 different combinations $I^{(0)}/I_{\ast }=0.1,1,10$ , $\unicode[STIX]{x0394}\unicode[STIX]{x1D707}=0,0.1,0.2,0.4$ and $M=10$ ghost nodes. The number of mesh points ranged from 700 to 1000, with FDA step-size $\unicode[STIX]{x0394}y=0.1$ , and initial guesses of the form $k_{2}=u_{0}^{\prime }=20,30,40,50$ .

Stability of the initial steady-state profile $u_{s}(y)$ is measured by the mean square departure $\unicode[STIX]{x1D716}=\sum _{i}|u(t_{max},y_{i})-u_{s}(y_{i})|^{2}$ at the time $t=t_{max}$ required for convergence of the unsteady solution. The most stable shear band for various values of the parameters $\unicode[STIX]{x1D707}_{0},\unicode[STIX]{x1D707}_{\infty },I^{(0)}$ can be determined by fixing one member of the pair $k_{2},u_{0}^{\prime }$ , where $u_{0}^{\prime }=u_{s}^{\prime }(0)$ , and minimizing $\unicode[STIX]{x1D716}$ with respect to the other one.

However, a few initial calculations revealed that $k_{2,opt}\hspace{2.22198pt}\dot{=}\hspace{2.22198pt}u_{0}^{\prime }$ , so that it is somewhat more efficient to minimize the quantity $\unicode[STIX]{x1D716}$ with respect to both members of the pair. This is accomplished by means of the MATLAB^® program fminunc.m, with initial guess $k_{2}=u_{0}^{\prime }$ , for various values of $u_{0}^{\prime }$ . Thus, with initial guess $u_{0}^{\prime }=k_{2}=28$ , one finds $u_{0}^{\prime }=28$ and $k_{2}=27.8$ as the optimal or maximally stable pair, for which the original, perturbed and final shear bands are shown in figure 9. In this calculation, the $y$ interval is divided into $N$ nodes, with $M\ll N$ ghost nodes at each end of the interval. As indicated in figure 9, the sinusoidal perturbation is suppressed at the ghost nodes in order to maintain the correct boundary conditions on the unsteady solution.

Contrary to our initial hope of establishing a unique point $k_{2},u_{0}^{\prime }$ , we find instead a fairly definite locus of optimal points in the $k_{2}$ – $u_{0}^{\prime }$ plane with $k_{2,opt}\dot{=}u_{0,opt}^{\prime }$ . This is summarized by the scatter plot in figure 10 of 72 different runs with diverse values of the various parameters. The least-squares straight line shown there is given by $u_{0}^{\prime }=0.998k_{2}$ , the root mean square fractional error of which is approximately 1 %. Remarkably, the results appear to be almost independent of the sinusoidal perturbation amplitude $A$ , although $A$ does have some influence on the convergence of the unsteady solution.

For completeness, we have also investigated the stability of the homogeneous shear field with above sinusoidal perturbation, i.e. the stability of the state with $u_{s}=y$ in (4.11). With $u=0$ imposed at $y=0$ as the only constraint, we find unbounded growth near $y=0$ without achieving a profile that resembles our steady shear band. Although one might conclude that our steady shear band does not represent a general point of attraction in the space of steady-state solutions, and is perhaps attainable only through 2D effects, we are inclined to attribute the form of the instability to our finite-difference implementation of the method of lines. While it would be interesting to employ more robust spectral methods, this would take us well beyond the scope of the present paper.

5 Pure shear

To illustrate the effects of base-state shearing, we consider a base state defined by $\boldsymbol{v}^{(0)}=(x_{1},-x_{2})/2$ , with

(5.1) $$\begin{eqnarray}\unicode[STIX]{x1D63F}^{(0)}\equiv \unicode[STIX]{x1D647}^{(0)}=\frac{1}{\surd 2}\unicode[STIX]{x1D640}^{(0)}=\frac{1}{2}\left[\begin{array}{@{}cc@{}}1 & 0\\ 0 & -1\end{array}\right].\end{eqnarray}$$

Hence, for planar disturbances, the relation for $\boldsymbol{k}(t,\unicode[STIX]{x1D73F})$ in (3.11) becomes

(5.2) $$\begin{eqnarray}\boldsymbol{k}=\exp (-t\unicode[STIX]{x1D647}^{(0)\text{T}})\unicode[STIX]{x1D73F},\quad \text{with }k_{1}=\unicode[STIX]{x1D705}_{1}\unicode[STIX]{x1D705}_{2}/k_{2}\text{ and }k_{2}=\unicode[STIX]{x1D705}_{2}\exp (t/2).\end{eqnarray}$$

Then, by the methods employed above, we find the eigenvalue $\unicode[STIX]{x1D706}$ of $\unicode[STIX]{x1D63C}$ with largest real part to be

(5.3) $$\begin{eqnarray}\displaystyle & & \displaystyle \unicode[STIX]{x1D706}=\frac{1}{2\unicode[STIX]{x1D719}\unicode[STIX]{x1D6F7}_{1}}\left[4\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FE}k_{1}^{2}k_{2}^{2}-2\unicode[STIX]{x1D6F7}_{1}\unicode[STIX]{x1D6F7}_{2}+\unicode[STIX]{x1D719}\left(k_{1}^{2}-k_{2}^{2}-\frac{\unicode[STIX]{x1D6FC}}{\surd 2}k^{2}\right)+\unicode[STIX]{x1D6F7}_{3}^{1/2}\right],\nonumber\\ \displaystyle & & \displaystyle \quad \text{where }\unicode[STIX]{x1D6F7}_{1}=\left(1-\frac{\unicode[STIX]{x1D6FC}}{\surd 2}\right)k_{1}^{2}+\left(1+\frac{\unicode[STIX]{x1D6FC}}{\surd 2}\right)k_{2}^{2},\nonumber\\ \displaystyle & & \displaystyle \quad \text{and }\unicode[STIX]{x1D6F7}_{3}=2k_{1}^{2}k_{2}^{2}[8\unicode[STIX]{x1D6FD}^{2}\unicode[STIX]{x1D6FE}^{2}k_{1}^{2}k_{2}^{2}+4\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D719}(k_{1}^{2}-k_{2}^{2}-\unicode[STIX]{x1D6FC}k^{2}/\surd 2)+\unicode[STIX]{x1D719}^{2}(\unicode[STIX]{x1D6FC}^{2}-2)].\qquad\end{eqnarray}$$

When convection is neglected, we find that

(5.4a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D6F7}_{3}=16\unicode[STIX]{x1D6FD}^{2}\unicode[STIX]{x1D6FE}^{2}k_{1}^{4}k_{2}^{4}\quad \text{and}\quad \unicode[STIX]{x1D706}=\frac{4\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FE}k_{1}^{2}k_{2}^{2}-\unicode[STIX]{x1D6F7}_{1}\unicode[STIX]{x1D6F7}_{2}}{\unicode[STIX]{x1D719}\unicode[STIX]{x1D6F7}_{1}}.\end{eqnarray}$$

Figure 11 shows the effects of $\unicode[STIX]{x1D712}$ and convection on stability, and figure 12 illustrates oscillatory behaviour. It should be noted that the unstable regions in both figures are rotated by approximately $45^{\circ }$ relative to the corresponding regions for simple shear, which is relevant to the associated shear-band models.

Oscillatory behaviour occurs when $\unicode[STIX]{x1D6F7}_{3}<0$ in (5.3), or in the region

(5.5) $$\begin{eqnarray}8\unicode[STIX]{x1D6FD}^{2}\unicode[STIX]{x1D6FE}^{2}k_{1}^{2}k_{2}^{2}+4\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D719}(k_{1}^{2}-k_{2}^{2}-\unicode[STIX]{x1D6FC}k^{2}/\surd 2)+\unicode[STIX]{x1D719}^{2}(\unicode[STIX]{x1D6FC}^{2}-2)<0,\end{eqnarray}$$

the boundaries of which are readily found to be

(5.6) $$\begin{eqnarray}k_{1}=\pm \left[\frac{\unicode[STIX]{x1D719}}{2\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FE}}\left(1+\frac{\unicode[STIX]{x1D6FC}}{\surd 2}\right)\right]^{1/2},\end{eqnarray}$$

upon solving for $k_{2}^{2}$ in terms of $k_{1}^{2}$ . These are shown in the stability diagram of figure 12, where, in contrast to simple shear, there now exist both stable and unstable oscillations.

Since $k_{1}k_{2}=\unicode[STIX]{x1D705}_{1}\unicode[STIX]{x1D705}_{2}$ , a constant, from (5.2), the asymptotic behaviour with or without convection is stable,

(5.7) $$\begin{eqnarray}\unicode[STIX]{x1D706}\sim -\frac{1}{\unicode[STIX]{x1D719}}\unicode[STIX]{x1D6F7}_{2}\sim -\frac{1}{\unicode[STIX]{x1D719}}(\unicode[STIX]{x1D6FE}k_{2}^{2}+2\unicode[STIX]{x1D712}k_{2}^{4}),\quad \text{for }t\rightarrow \infty ,\end{eqnarray}$$

with relative error $O(k_{2}^{-2})=O(\text{e}^{t})$ . Clearly, we have stability even for $\unicode[STIX]{x1D712}=0$ .

Figure 11.

Stability for pure shear without convection, for $\unicode[STIX]{x1D712}$ values of (a) $0$ and (b) $0.001$ , and with convection, for $\unicode[STIX]{x1D712}$ values of (c) $4\times 10^{-5}$ and (d) $0.001$ , for $I=0.0001$ , $\unicode[STIX]{x1D719}=0.5$ and $p^{(0)}=1$ .

Figure 12.

Oscillatory behaviour for pure shear resulting from convection, for $\unicode[STIX]{x1D712}=0$ , $I=0.0001$ , $p^{(0)}=1$ and $\unicode[STIX]{x1D719}=0.5$ . The oscillation frequency is non-zero within the central regions bounded by dash-dotted vertical lines. The solid and dotted curves are contours of the real part of $\unicode[STIX]{x1D706}$ , with the former representing neutral stability.

Once again, one may derive a model for the shear band resulting from material instability. Thus, making use once more of the MATLAB^® program ‘fminunc’ to determine the unconstrained minimum in $-\unicode[STIX]{x1D706}(\boldsymbol{k})$ numerically, for parameter values $\unicode[STIX]{x1D707}_{0}=0.383,\unicode[STIX]{x1D707}_{\infty }=0.643,I_{\ast }=0.279,I^{(0)}=0.001,\unicode[STIX]{x1D712}=10^{-6}$ , we find $\boldsymbol{k}_{max}\approx \pm (26.0,22.5)$ . This corresponds roughly to the classical ${\approx}45^{\circ }$ shear-band orientation found in previous experiments and numerical simulations. Figure 13 shows the growth of the Fourier mode, and it is seen from panel (b) that the transient growth can become quite large before being quenched by a combination of wavevector stretching and damping by higher gradients.

The above result leads us to express (2.11), with $\boldsymbol{g}=0$ and $\unicode[STIX]{x1D735}p$ , on a coordinate system oriented at $45^{\circ }$ to that involved in (5.1), reducing it to the form (4.1) and (5.1) to the form (3.13). Therefore, the analysis of the shear band reduces essentially to that of § 4, with appropriate modification of the boundary condition on $u^{\prime }(0)$ . On this new coordinate system, the above value of $\boldsymbol{k}_{max}$ becomes ${\approx}\pm (2,34)$ as opposed to the value $\pm (2,28)$ found in § 4 for simple shear. Given the rough agreement of these values and in light of the non-unique solution to our shear-band model, it hardly seems worthwhile to repeat the calculation leading to figure 8 with an only slightly modified boundary condition on $u^{\prime }(0)$ .

Figure 13.

The magnitude $|\hat{\unicode[STIX]{x1D74A}}(\unicode[STIX]{x1D73F},t)|$ of initially stable and unstable Fourier modes versus $t$ for planar pure shear for $\unicode[STIX]{x1D712}=10^{-6}$ and initial conditions (a) $\unicode[STIX]{x1D73F}=(-0.5,-0.866)$ and (b) $\unicode[STIX]{x1D73F}=(1,0.0001)$ .

6 Concluding remarks

Our major findings are adequately summarized in the abstract and the introduction. To briefly recapitulate the most important of these, we find for both planar simple shear and pure shear that the linear instability of the $\unicode[STIX]{x1D707}(I)$ model identified by Barker et al. (Reference Barker, Schaeffer, Bohorquez and Gray2015) is modified through convection by the base flow, giving way to long-time stability induced by Kelvin wavevector stretching. The addition of gradient effects via the vdW–C–H model provides a wavenumber cutoff that serves to stabilize the dynamical equations over the entire time domain, to regularize the quasi-static field equations and to assign a diffuse length scale to eventual shear bands.

We find that steady shear bands are stable against steady parallel sinusoidal shear fields, provided that the normal velocity gradient of the shear band is very nearly equal to the wavenumber of the sinusoidal perturbation. To obtain a unique shear band, it is necessary to assume some preferred wavenumber, e.g. the most unstable one according to linear theory. A challenge for future work is to elucidate shear-band formation in the presence of more complex spatial perturbations to homogeneous shearing.

In summary, we conclude that the (Hadamard) short-wavelength instability found by Barker et al. (Reference Barker, Schaeffer, Bohorquez and Gray2015) is connected to the loss of ellipticity in the quasi-static field equations. Although transient and eventually quenched by wavevector stretching according to linear theory, the instability is doubtless problematical for numerical simulation and it may also trigger nonlinear instabilities. The addition of higher-gradient effects like those considered in the present work should go a long way towards alleviating such problems.

As a practical matter, the present work may suggest a stratagem for numerical simulation of materially unstable viscoplasticity, in which a transient form of the vdW–C–H or similar regularization is invoked to minimize transient instabilities, with the possibility of describing the structure of shear bands on longer time scales if so desired. Our simple model of a shear band based on the vdW–C–H model suggests that spatial boundary conditions, required in principle by any higher-gradient model, can be chosen somewhat arbitrarily, as the higher gradients tend to have a spatially localized domain of influence.

As additional future work, it would be worthwhile to apply the current theory to other homogeneous shear flows, such as axisymmetric straining of the kind that arises in the standard quasi-static tests of soil mechanics or in more rapid converging-hopper flows. While somewhat distinct from the issue of material instability, an investigation should be carried out of the coupling of pressure gradients in the base flow to the perturbed momentum balance, an effect neglected in Ref. 1 and the present study.

As a deeper theoretical issue, it would be interesting to investigate the possible relation between the weakly non-local model of the present study and the fully non-local variants of the $\unicode[STIX]{x1D707}(I)$ model proposed by Pouliquen & Forterre (Reference Pouliquen and Forterre2009) and more recently by Kamrin and co-workers (see, e.g., Henann & Kamrin Reference Henann and Kamrin2014). We recall that the latter model can be tied to a Ginzburg–Landau formalism, which shares a certain kinship to the vdW–C–H model of the present study (Gurtin Reference Gurtin1996).