Linear instability of shear thinning pressure driven channel flow

TL;DR

This paper investigates the linear instability of shear thinning viscoelastic channel flows using three constitutive models, revealing that instability occurs when shear thinning exceeds a critical level and is driven by normal stress variations.

Contribution

It provides a comparative theoretical analysis of flow instability mechanisms across three widely used viscoelastic models, highlighting the role of shear thinning and normal stress gradients.

Findings

Instability occurs when shear thinning exceeds a critical logarithmic slope.

Critical pressure drop depends on shear thinning degree and normal stress derivative.

Destabilization involves a quasi-interface in some models, but the mechanism's universality remains uncertain.

Abstract

We study theoretically pressure driven planar channel flow of shear thinning viscoelastic fluids. Combining linear stability analysis and full nonlinear simulation, we study the instability of an initially one-dimensional base state to the growth of two-dimensional perturbations with wavevector in the flow direction. We do so within three widely used constitutive models: the microscopically motivated Rolie-Poly model, and the phenomenological Johnson-Segalman and White-Metzner models. In each model, we find instability when the degree of shear thinning exceeds some level characterised by the logarithmic slope of the flow curve at its shallowest point. Specifically, we find instability for $n<n^*$, with $n^*\approx 0.21, 0.11$ and $0.30$ in the Rolie-Poly, Johnson-Segalman and White-Metzner models respectively. Within each model, we show that the critical pressure drop for the onset of instability obeys a criterion expressed in terms of this degree of shear thinning, $n$, together with the derivative of the first normal stress with respect to shear stress. Both shear thinning and rapid variations in first normal stress across the channel are therefore key ingredients driving the instability. In the Rolie-Poly and Johnson-Segalman models, the underlying mechanism appears to involve the destabilisation of a quasi-interface that exists in each half of the channel, across which the normal stress varies rapidly. In the White-Metzner model, no such quasi-interface exists, but the criterion for instability nonetheless appears to follow the same form as in the Rolie-Poly and Johnson-Segalman models. This presents an outstanding puzzle concerning any possibly generic nature of the instability mechanism. We finally make some briefly comments on the Giesekus model, which is rather different in its predictions from the other three.

Equations

$$\mathbf{\Sigma }=\mathbf{\sigma }+2\eta \mathbf{D}-p\mathbf{I}.$$

$$\nabla \cdot \mathbf{v}=0.$$

$$\nabla \cdot \mathbf{\Sigma }=0.$$

$${\partial }_t\mathbf{W}+\mathbf{v}\cdot \nabla \mathbf{W}=\mathbf{K}\cdot \mathbf{W}+\mathbf{W}\cdot {\mathbf{K}}^T-\frac{1}{{\tau }_d}\left(\mathbf{W}-\mathbf{I}\right)-\frac{2(1-A)}{{\tau }_R}\left[\mathbf{W}+\beta A^{-2\delta }\left(\mathbf{W}-\mathbf{I}\right)\right].$$

$${\partial }_t\mathbf{W}+\mathbf{v}\cdot \nabla \mathbf{W}=\mathbf{K}\cdot \mathbf{W}+\mathbf{W}\cdot {\mathbf{K}}^T-\frac{1}{\tau }\left(\mathbf{W}-\mathbf{I}\right)-\frac{2}{3}\mathrm{Tr}(\mathbf{K}\cdot \mathbf{W})\left[\mathbf{W}+\beta (\mathbf{W}-\mathbf{I})\right].$$

$${\partial }_t\mathbf{W}+\mathbf{v}\cdot \nabla \mathbf{W}=(\mathbf{\Omega }\cdot \mathbf{W}-\mathbf{W}\cdot \mathbf{\Omega })+a(\mathbf{D}\cdot \mathbf{W}+\mathbf{W}\cdot \mathbf{D})-\frac{1}{\tau }\left(\mathbf{W}-\mathbf{I}\right),$$

$${\partial }_t\mathbf{W}+\mathbf{v}\cdot \nabla \mathbf{W}=\mathbf{K}\cdot \mathbf{W}+\mathbf{W}\cdot {\mathbf{K}}^T-\frac{1}{{\tau }^n{\dot{\gamma }}^{n-1}}\left(\mathbf{W}-\mathbf{I}\right).$$

$${\partial }_t\mathbf{W}+\mathbf{v}\cdot \nabla \mathbf{W}=\mathbf{K}\cdot \mathbf{W}+\mathbf{W}\cdot {\mathbf{K}}^T-\frac{1}{\tau }\left(\mathbf{W}-\mathbf{I}\right)-\frac{\alpha }{\tau }{\left(\mathbf{W}-\mathbf{I}\right)}^2,$$

$$\mathbf{\Sigma }=\mathbf{\sigma }+2\eta \mathbf{D}-p\mathbf{I},$$

$$\nabla \cdot \mathbf{\Sigma }=0.$$

$$\nabla \cdot \mathbf{v}=0.$$

$$\mathbf{\sigma }=\mathbf{\sigma }(\mathbf{W}),$$

$${\partial }_t\mathbf{W}+\mathbf{v}\cdot \nabla \mathbf{W}=F(\mathbf{W},\nabla \mathbf{v},\xi ).$$

$$\mathbf{v}(x,y,t)$$

$$\mathbf{W}(x,y,t)$$

$$\mathbf{\sigma }(x,y,t)$$

$$\mathbf{\Sigma }(x,y,t)$$

$$p(x,y,t)$$

$${\partial }_t{\delta \mathbf{W}}_q+iq{u}_0\delta {\mathbf{W}}_q+\delta v_q{\partial }_y{\mathbf{W}}_0=F({\mathbf{W}}_0+\delta {\mathbf{W}}_q,\nabla {\mathbf{v}}_0+e^{-iqx}\nabla (\delta {\mathbf{v}}_q{e}^{iqx}),\xi ),$$

$$\eta {\nabla }^2\delta {\mathbf{v}}_q+\nabla \cdot \delta {\mathbf{\Sigma }}_q-\nabla \delta p_q=0.$$

$$\frac{d{N}_1}{dy}=\frac{d\Sigma }{dy}\frac{d{N}_1}{d\Sigma }=P^{\prime }\frac{d{N}_1}{d\Sigma }=P^{\prime }\frac{d{N}_1}{d\dot{\gamma }}/\frac{d\Sigma }{d\dot{\gamma }}.$$

$$\frac{d{N}_1}{dy}{|}_{\max }>h\left(\frac{d\log \Sigma }{d\log \dot{\gamma }}{|}_{\min }\right),$$

$$P^{\prime \ast }\frac{d{N}_1}{d\Sigma }{|}_{\max }>h\left(\frac{d\log \Sigma }{d\log \dot{\gamma }}{|}_{\min }\right),$$

$$h(x)={\alpha }_0+{\alpha }_1/({\alpha }_2-x),$$

$${\Sigma }^{\ast }\sim \frac{G}{\sqrt{1-a^2}}.$$

$$\frac{d{N}_1}{dy}{|}_{\max }=P^{\prime \ast }\frac{d{N}_1}{d\Sigma }{|}_{\max }>h\left(\frac{d\log \Sigma }{d\log \dot{\gamma }}{|}_{\min }\right),$$

$$h(x)={\alpha }_0+{\alpha }_1/({\alpha }_2-x),$$

$$P^{\prime }{L}_y-2{\Sigma }^{\ast }>0.$$

$$\frac{d{N}_1}{dy}{|}_{\max }=P^{\prime \ast }\frac{d{N}_1}{d\Sigma }{|}_{\max }>h\left(\frac{d\log \Sigma }{d\log \dot{\gamma }}{|}_{\min }\right),$$

$$h(x)={\alpha }_0+{\alpha }_1/({\alpha }_2-x),$$

Full text

I Introduction

Viscoelastic fluids commonly show flow instabilities in which an initially simple flow state gives way to a more complicated one Larson (1992). Such instabilities arise even in the creeping flow regime of low Reynolds number, with the fluid inertia playing little or no role. They stem instead from the nonlinear way in which the viscoelastic dynamics interacts with an imposed flow. Instability typically sets in beyond a critical value of the imposed flow rate scaled by the viscoelastic relaxation time, as expressed by the dimensionless Weissenberg number. For flows with curved streamlines, such instabilities are widespread. They include the viscoelastic Taylor-Couette instability Larson et al. (1990), instability of viscoelastic cone-plate flow Olagunju (1995), and viscoelastic turbulence in plate-plate flow Schiamberg et al. (2006). A criterion for instability based on the Weissenberg number and the curvature of the flow streamlines was put forward in Ref. Pakdel and McKinley (1996).

Considerable work has also been devoted to the question of whether creeping viscoelastic flows that have straight streamlines (in the initially simple flow state) are also unstable. For fluids that have a constant viscosity as a function of flow rate (ie, fluids that don’t shear thin, or thicken), but that are nonetheless still viscoelastic, a one-dimensional (1D) base state has been shown to be linearly stable against the growth of 2D perturbations with wavevector in the flow direction, $x$, in both boundary driven planar Couette Wilson et al. (1999) and pressure driven Poiseuille flow Ho and Denn (1977), within the Oldroyd B model. However, a nonlinear analysis within the same model suggests a possible sub-critical (finite amplitude) instability at high enough Weissenberg number, given perturbations of large enough amplitude Morozov and van Saarloos (2005). Experiments testing this prediction can be found in Refs. Qin and Arratia (2017); Pan et al. (2013).

The present paper concerns the stability properties of shear thinning fluids (ie, fluids that show a decrease in viscosity with increasing shear rate), in a commonly studied flow geometry with straight streamlines: that of pressure driven planar channel flow. An earlier stability analysis performed within the shear thinning White-Metzner model predicted an initially 1D base state to be linearly unstable to the growth of 2D perturbations with wavevector in the flow direction, for fluids that shear thin strongly enough, at large enough imposed pressure drops Wilson and Rallison (1999); Wilson and Loridan (2015); Castillo and Wilson (2017). Instability has also been predicted in the shear thinning Giesekus and Phan-Thien Tanner models Grillet et al. (2002) and also, very recently, in a model of thixoelastoviscoplastic flow Castillo and Wilson (2018). Instability of shear thinning polymeric solutions in pressure driven flow has recently been confirmed experimentally Picaut et al. (2017); Bodiguel et al. (2015); Poole (2016); Bonn et al. (2011). Instabilities have also long been known to arise during extrusion of polymeric fluids Denn (2001), although often with wall slip involved: a phenomenon that we ignore here.

Motivated by these works, the present study set out to address several outstanding questions concerning the stability properties of shear thinning fluids in pressure driven flows. What level of shear thinning is in general required for instability to arise? Is such an instability generic across all fluids and constitutive models that show (at least) this level of shear thinning? Is the required level of shear thinning independent of all other constitutive properties of the fluid, or does it also depend on (for example) normal stresses? For any fluid showing at least the required minimal level of shear thinning, what is the critical pressure drop required for instability? Can we uncover a criterion for this critical pressure drop in terms of easily measurable rheological signatures, thereby providing a practical guide for when an experimentalist can expect instability? Is the mechanism of instability the same in all constitutive models? And is instability predicted by constitutive models that are built on a microscopic understanding of the underlying polymeric relaxation processes, or only in the phenomenological models studied to date? As will be seen in what follows, some of these questions will prove to have relatively straightforward answers; others remain unresolved by this study.

To address these questions, we have performed an extensive linear stability analysis together with full nonlinear numerical simulations within three widely used constitutive models: the microscopically motivated Rolie-Poly model Likhtman and Graham (2003), and the phenomenoligical Johnson-Segalman Johnson and Segalman (1977), and White-Metzner White and Metzner (1963) models. We shall also make some comments about the Giesekus model Giesekus (1982), which is rather different in its predictions from the other three, in our concluding discussion Sec. VII. As noted above, Giesekus Grillet et al. (2002) and White-Metzner Wilson and Rallison (1999) have already been shown to predict instability of pressure driven shear thinning channel flow. The purpose of studying them again here is to try to elucidate the degree to which the predictions of the different models can – or cannot – be understood within a common framework.

All of the models considered, except White-Metzner, are capable of strong enough shear thinning to give shear banding. We restrict ourselves here to parameter regimes in which thinning is less pronounced, with no banding.

We find an initially 1D base state, in which the flow varies only in the flow-gradient direction, $y$, to indeed be linearly unstable to the onset of 2D perturbations with wavevector in the flow direction, $x$, in all four models. (In Giesekus and White-Metzner this finding simply confirms the earlier predictions; we believe the result to be new in Rolie-Poly and Johnson-Segalman.) Within Rolie-Poly, Johnson-Segalman and White-Metzner, we further perform a comprehensive study of the phase behaviour of this instability as a function of the model parameters and the imposed pressure drop. We thereby calculate the minimal degree of shear thinning needed for instability, within each model, expressed in terms of the logarithmic slope of the flow curve of shear stress as a function of shear rate $\Sigma(\dot{\gamma})$. We find instability below a critical value, $n^{*}$, of the logarithmic slope of the flow curve at its shallowest point, $n=d\log\Sigma/d\log\dot{\gamma}|_{\rm min}$ with $n^{*}\approx 0.21$ in Rolie-Poly, $n^{*}\approx 0.11$ in Johnson-Segalman, and (as calculated previously in Ref. Wilson and Rallison (1999)) $n^{*}\approx 0.30$ in White-Metzner. In the Giesekus model, in contrast, $n^{*}$ depends strongly on the solvent viscosity, and is furthermore much lower than in the other models, varying (with solvent viscosity) in the range [math] to $0.04$. The values in the Rolie-Poly and White-Metzner models are broadly consistent with instability having been observed for the shear thinning exponents of $n=0.21$ and $n=0.19$ in the experimental studies of Refs. Bodiguel et al. (2015) and Poole (2016) respectively. Only the White-Metzner model appears consistent with the presence of instability for the value of $n=0.29$ in the experiments of Ref. Picaut et al. (2017).

Within each of the Rolie-Poly, Johnson-Segalman and White-Metzner models, we show that the critical pressure drop for the onset of instability (expressed in units of the fluid modulus and the channel width) follows, to a good or reasonable level of approximation, a scaling function expressed in terms of (i) the logarithmic slope of the flow curve at its shallowest point (i.e., at the point of minimum slope), and (ii) the maximum derivative of the the first normal stress with respect to the shear stress. In each of these three models, the scaling function is characterised by three dimensionless fitting parameters. Once these fitting parameters have been determined (for any model), a reasonable level of scaling collapse of the critical pressure drop at onset is achieved for the data obtained across all values of the model parameters, for that model.

Our results thereby show that both strong shear thinning and large variations in normal stress across the channel tend to predispose a flow to instability. The criterion just described is only a partial success, however. In particular, the values of the three fitting parameters are different from model to model. We are therefore unable to provide a criterion that is truly universal. Ideally, these fitting parameters would be recast in terms of any additional relevant dimensionless constitutive properties, leading to a re-expressed set of fitting parameters that have the same values across all models. We have been unable to achieve that in this work. It remains an open future goal, which may however be unattainable.

In the Rolie-Poly and Johnson-Segalman models, we give evidence suggesting that the mechanism of instability involves the destabilisation of a quasi-interface at some location in each half of the channel, across which the first normal stress varies rapidly as a function of position, and on either side of which the fluid flows with a different local viscosity. In this way, the instability closely resembles that arising at the interface between layered fluids Wilson and Rallison (1997). Indeed, an extreme form of shear thinning arises when the underlying constitutive curve $\Sigma(\dot{\gamma})$ (for states of homogeneous shear flow) is non-monotonic. The flow state as predicted by a 1D calculation then comprises bands of differing shear rates and normal stresses. This 1D banded state is known to be unstable to the formation of two dimensional perturbations with wavevector in the flow direction, with the instability driven by a jump in first normal stress across the interface between the bands Fielding (2005); Fielding and Olmsted (2006); Fielding and Wilson (2010). The instability found here in the Johnson-Segalman and Rolie-Poly models therefore resembles that of an interface between shear bands, but now in a regime where the flow is no-longer truly banded.

In the White-Metzner model, no such quasi-interface exists, and the eigenfunction characterising the instability is spread diffusely across the channel. An unresolved question is why the same form of criterion for instability onset (albeit with different values of the three fitting parameters) appears to hold in the White-Metzner model as in the Rolie-Poly and Johnson-Segalman models.

The paper is structured as follows. In Sec. II we introduce the constitutive models that we shall study. In Sec. III we define the flow geometry to be considered, then discuss units and parameter values in Sec. IV. In Sec. V we outline our calculation methods, which include linear stability and nonlinear simulations. We then present our results in Sec. VI, for the Rolie-Poly, Johnson-Segalman and White-Metzner models in turn. Finally we present our conclusions and discuss outstanding challenges in Sec. VII.

II Models

We assume the total stress $\bm{{\Sigma}}(\bm{{r}},t)$ in a fluid element at position $\bm{{r}}$ at time $t$ to comprise a viscoelastic contribution $\bm{{\sigma}}(\bm{{r}},t)$ from the entangled polymer chains, a Newtonian contribution of viscosity $\eta$, and an isotropic pressure, with $\bm{{I}}$ the unit tensor:

[TABLE]

The Newtonian part, $2\eta\bm{{D}}$, may arise from a true solvent, or model any fast polymeric relaxation modes that are not ascribed to the viscoelastic stress. Here $\bm{{D}}=\frac{1}{2}(\bm{{K}}+\bm{{K}}^{T})$ is the symmetric part of the velocity gradient tensor, $K_{\alpha\beta}=\partial_{\beta}v_{\alpha}$, where $\bm{{v}}(\bm{{r}},t)$ is the fluid velocity field. Below we shall also use the antisymmetric tensor $\bm{{\Omega}}=\frac{1}{2}(\bm{{K}}-\bm{{K}}^{T})$. The isotropic pressure $p(\bm{{r}},t)$ acts to ensure that the flow remains incompressible:

[TABLE]

Throughout we shall work in the limit of creeping flow, in which the condition of force balance states that the total stress tensor must be divergence free:

[TABLE]

Accordingly, any instabilities that we report stem not from inertial effects, but are of purely viscoelastic origin.

We write the viscoelastic stress in terms of a constant elastic modulus $G$ and a molecular conformation tensor $\bm{{W}}(\bm{{r}},t)$ that describes the deformation of the polymer chains relative to an isotropic undeformed state, with $\bm{{\sigma}}=G\bm{{W}}$ in the Rolie-poly model and $\bm{{\sigma}}=G\left[\bm{{W}}-\bm{{I}}\right]$ in all the other models that we shall consider. The dynamics of $\bm{{W}}(\bm{{r}},t)$ in flow is then prescribed by a viscoelastic constitutive equation Larson (2013). We shall consider four shear thinning constitutive models that are widely used across the rheology literature: the Rolie-Poly model Likhtman and Graham (2003), the Johnson-Segalman model Johnson and Segalman (1977), the Giesekus model Giesekus (1982), and the White-Metzner model White and Metzner (1963). (We map out the phase behaviour of three of these in detail, and comment only briefly on the Giesekus model.)

Three of these models admit such strong shear thinning that they predict shear banding in some regions of their parameter space. However, all the results presented below are for parameter regimes in which shear thinning is more moderate, with the constitutive curve of shear stress as a function of shear rate, $\Sigma_{xy}(\dot{\gamma})$, being monotonic, precluding shear banding. Therefore, in contrast to studies of instabilities between shear bands Fielding (2005), we do not include a stress diffusion term in our constitutive equations.

II.1 Rolie-Poly model

The flow behaviour of an entangled polymeric fluid can be modelled at a microscopic level by considering the molecular dynamics of the constituent polymer chains. The basic notion is that any test chain of interest has its motion laterally constrained by entanglements with other chains. In a mean field approach, the constraining influence of these entanglements is represented by an effective tube to which that test chain is (initially) confined Doi and Edwards (1988). The GLAMM model Graham et al. (2003) provides a stochastic equation for the dynamics of such a test chain in its tube. This model is however computationally prohibitive to work with in fluid dynamical simulations. We shall therefore work with a simpler approximation, in which GLAMM is projected onto a single mode description Likhtman and Graham (2003). The resulting rolie-poly (RP) model prescribes the dynamics of the polymer conformation tensor as follows:

[TABLE]

In this equation, $\tau_{d}$ is the characteristic timescale on which a test chain escapes its tube of constraints via a process of 1D curvilinear diffusion along its own length, known as reptation. The Rouse time $\tau_{R}$ is the much shorter timescale on which the degree of chain stretch relaxes, by the lowest mode of curvilinear breathing of a chain within its tube. The chain stretch $A=\sqrt{3/T\,}$, where $T=\text{tr}\,\bm{{\bm{{W}}}}$. The ratio $\tau_{d}/\tau_{R}$ of these two relaxation times increases with the number of entanglements per chain $Z$. We work throughout with highly entangled polymers and use the simpler, non-stretching version of the Rolie-Poly model, which follows from Eqn. 4 in the limit $\tau_{R}\to 0$:

[TABLE]

We thereby assume that flow rates, and the rates associated with the onset of any instability, remain modest compared with $1/\tau_{R}$. Note that we now denote $\tau_{d}$ simply by $\tau$. The parameter $\beta$ describes the level of “convective constraint release” (CCR) Ianniruberto and Marrucci (2014), in which the relaxation of polymer chain stretch also relaxes entanglement points, thereby also allowing relaxation of tube orientation. It has a range $0\leq\beta\leq 1$, with no current consensus on its value within that range. For $\beta<1$ the model is capable of capturing shear banding, at low values of the solvent viscosity $\eta$. In what follows, we consider only values of $\beta,\eta$ that are outside the banding regime.

II.2 Johnson-Segalman model

Alongside the Rolie-Roly model, which was derived by considering the molecular relaxation processes of reptation, chain stretch relaxation and CCR as just discussed, we shall also consider some constitutive models that are instead phenomenological in origin. We start with the Johnson Segalman model Johnson and Segalman (1977). This is closely related to the Oldroyd B model Spiess (1987), which considers the dynamics of an ensemble of dumbbells in a solvent. Each dumbbell comprises two beads connected by a spring, and is taken to model a single polymer molecule. Each bead is subject to Brownian motion, Stokes drag, and the spring force. In the Oldroyd B model, each dumbbell is assumed to deform in flow in an affine way, following the background solvent. In contrast, in the Johnson-Segalman model each dumbbell is instead taken to slip relative to the solvent in a non-affine way. This phenomenon is described by a slip parameter $a$, which obeys $|a|\leq 1$. The molecular conformation tensor then obeys:

[TABLE]

with a relaxation time $\tau$. For $|a|<1$, the model captures shear thinning. (For $\eta<0.125$, it further captures shear banding. We set $\eta>0.125$ in all of what follows.) For $a=1$, the Oldroyd B model is recovered, with no shear thinning. We note that the Johnson-Segalman model predicts unphysical time-dependent behaviour in shear startup at high strain rates. This pathology appears not to cause problems with the stability calculation here, however: even at high strain rates, the base state about which we linearise is stationary.

II.3 White-Metzner model

The White-Metzner model captures shear thinning by invoking a relaxation time that depends on the frame-invariant strain rate, $\dot{\gamma}\equiv\sqrt{\bm{{D}}:\bm{{D}}}$, as follows:

[TABLE]

The power law index $0<n\leq 1$, giving shear thinning for $n<1$, and recovering the Oldroyd B model for $n=1$.

Note that the White-Mezner model gives power law fluid behaviour even in the limit of $\dot{\gamma}\tau\to 0$, lacking any terminal Newtonian regime of weak shear. As such, it should be treated with some caution in modelling the experimental polymeric solutions that inspired this work, at least in the limit of low strain rates (which are always obtained near the centre of the channel, even at high pressure drops).

II.4 Giesekus model

Another phenomenological model commonly used to describe concentrated polymeric solutions or melts is due to Giesekus Giesekus (1982). It considers an anisotropic drag on polymer chains that are oriented due to flow, and models this via an anisotropy parameter $\alpha$, which must lie in the range $0\leq\alpha\leq 1$. The molecular conformation evolves according to

[TABLE]

with relaxation time $\tau$. For $\alpha>0.5$ the model is capable of capturing shear banding, at low enough values of $\eta$. We consider only $\alpha,\eta$ values that are outwith the banding regime in what follows. For $\alpha=0$, the Oldroyd B model is recovered.

III Flow Geometry

We consider a slab of fluid confined between flat parallel plates at $y=-L_{y}/2$ and $y=+L_{y}/2$. The flow is driven along the channel in the flow direction of positive $\hat{\bm{{x}}}$ by means of a negative pressure gradient per unit length along $x$, of magnitude $P^{\prime}$. The flow is assumed to remain translationally invariant in the vorticity direction $\hat{\bm{{z}}}$, and with no component of velocity in $\hat{\bm{{z}}}$. (These assumptions are in accordance with Squire’s theorem, but should be checked in future in fully 3D studies, in particular with respect to the potentially destabilising influence of second normal stress differences Lockett (1969). ) Boundary conditions of no-slip and no-permeation are assumed at the plates at $y=\pm L_{y}/2$. In the flow direction the channel has length $L_{x}$, with periodic boundary conditions.

Note that we have chosen to use $P^{\prime}$ as the parameter controlling the strength of flow in this study. An alternative choice, commonly used in other studies Grillet et al. (2002); Wilson and Rallison (1999), would be to instead characterise the strength of flow by a Weissenberg number, defined by the shear rate at some point across the channel (or the velocity on the centreline divided by the channel width), multiplied by the fluid’s viscoelastic relaxation time.

IV Units and parameter values

We report all results below in units of length in which the channel width $L_{y}=1$; units of time in which the basic polymeric relaxation timescale $\tau=1$; and units of mass in which the polymeric shear modulus $G=1$.

In these units, there remain just four parameters to explore. For each of the constitutive models that we consider, the first is the single parameter pertaining to the viscoelastic dynamics: the convective constraint release parameter $\beta$ in the Rolie-Poly model, the slip parameter $a$ in the Johnson-Segalman model, the power law index $n$ in the White-Metzner model, and the anisotropy parameter $\alpha$ in the Giesekus model. In any place where we need a single symbol to denote this viscoelastic parameter in a model-independent way, we use $\xi$. The second parameter is the Newtonian solvent viscosity $\eta$, which we take to be always small compared to the scale of the polymer viscosity $G\tau=1$. The third is the channel length $L_{x}$. In our linear stability analysis we report dispersion relations as a continuous function of the magnitude of the wavevector $q\hat{\bm{{x}}}$. However, it is important to realise that this is actually quantized as $q=n\pi/L_{x}$ for any channel of finite length $L_{x}$. In our nonlinear simulations, $L_{x}$ is an explicit parameter. The final parameter is the pressure drop, $P^{\prime}$, characterising the strength of the imposed flow.

V Calculation Methods

The flow models introduced in Sec. II above all have the same basic structure, which we recap as follows. The total stress

[TABLE]

In creeping flow, force balance requires:

[TABLE]

The condition of incompressibility gives:

[TABLE]

The viscoelastic stress

[TABLE]

in which the molecular conformation tensor $\bm{{W}}$ evolves according to a viscoelastic constitutive model. Above, we specified several different constitutive models. However all have the same general form:

[TABLE]

The parameter $\xi=\beta$ in the Rolie-Poly model, $\xi=a$ in the Johnson-Segalman model, $\xi=n$ in White-Metzner, and $\xi=\alpha$ in Giesekus. These equations must be solved subject to a constant imposed pressure gradient $-P^{\prime}\hat{\bm{{x}}}$ per unit length in $x$.

For any values of the model parameters $\xi,\eta$ and imposed pressure gradient $P^{\prime}$, we first calculate the one-dimensional (1D) stationary solution to Eqns. 8 - 12 that depends only on the flow gradient direction $y$. (We do so by numerically evolving the model equations to steady state, allowing spatial variations only in $y$.) We denote this initial base state $\bm{{v}}_{0}(y)$, $\bm{{\bm{{W}}}}_{0}(y)$, $\bm{{\sigma}}_{0}(y)$, $\bm{{\Sigma}}_{0}(y)$, $p_{0}(y)$. (It also depends on $\xi,\eta$ and $P^{\prime}$, though it would be cumbersome to write these dependencies explicitly.) The basic question that we now address is whether, for any $\xi,\eta,P^{\prime}$ and channel length $L_{x}$, this 1D base state is stable or unstable with respect to the growth of 2D perturbations that depend on the flow direction $x$ as well as $y$.

We address this question first by performing a linear stability analysis, in which we add to the base state $\bm{{v}}_{0}(y)$, $\bm{{\bm{{W}}}}_{0}(y)$, $\bm{{\sigma}}_{0}(y)$, $\bm{{\Sigma}}_{0}(y)$, $p_{0}(y)$ small amplitude perturbations as follows:

[TABLE]

Substituting these into the governing Eqns. 8 - 12, expanding in powers of the amplitude of the perturbations, and retaining only terms of first order in that amplitude, gives a set of linearised equations governing the dynamics of the perturbations. These equations are valid in the linear regime in which the perturbations remain small. In this linear regime, the $q-$modes are all independent of each other. The main quantities of interest are then, for each value of $q$, the eigenvalue $\omega_{q}$ with the largest real part, and the associated eigenfunction $\delta\bm{{v}}_{q}(y),\delta\bm{{\bm{{W}}}}_{q}(y),\delta\bm{{\sigma}}_{q}(y),\delta\bm{{\Sigma}}_{q}(y),\delta p_{q}(y)$. The sign of $\Re\omega_{q}$ then determines whether the associated perturbation, with a functional form prescribed by the eigenfunction, grows or decays. A positive value $\Re\omega_{q}>0$ signifies linear instability to the growth of a two-dimensional perturbation $\propto\exp(iqx)$. In contrast, if $\Re\omega_{q}<0$ for all modes, the base state is linearly stable.

The ‘fully linearised’ recipe of the previous paragraph gives an exact but cumbersome method for calculating the eigenmodes of the instability, as just described. In calculational practice, however, we used a different method, which is equivalent to that of the previous paragraph in the early time linear regime, as long as the perturbations remain small. Specifically, we write a ‘partly linearised’ viscoelastic constitutive equation, at any wavevector $q$, as follows:

[TABLE]

and the force balance equation as

[TABLE]

In Eqn. 14, $u_{0}$ is the base state speed, defined by $\bm{{v}}_{0}=u_{0}\hat{\bm{{x}}}$. $\delta v_{q}$ corresponds to the magnitude of the y-components of the perturbation velocity $\delta\bm{{v}}_{q}=\delta u_{q}\hat{\bm{{x}}}+\delta v_{q}\hat{\bm{{y}}}$. Note that in Eqn. 14, we have linearised advective terms of the viscoelastic constitutive equation exactly, but left the rest of the constitutive equation written in nonlinear form. Into these partly linearised equations we substitute the time-independent base state $\bm{{W}}_{0}(y)$ and $\bm{{v}}_{0}(y)$ as calculated above, and initialise spatially random (in y) perturbations $\delta\bm{{W}}_{q}$, $\delta\bm{{v}}_{q}$, etc, with a small amplitude. We then time-step these equations and observe the growth (or decay) of the perturbations in time. The slope of the logarithm of the amplitude of the perturbations as a function of time then gives the real part of the eigenvalue $\Re\omega_{q}$. (The imaginary part of the eigenvalue instead merely represents the advection rate of the perturbations Castillo and Wilson (2018) and is not reported here.) The flow pattern that emerges then gives the associated eigenfunction. We have checked this method against earlier calculations that used the fully linearised method, for the case of shear banded flow Fielding and Wilson (2010).

An important aim will be to relate these stability properties to the functional form of the 1D base state profiles of shear rate, $\dot{\gamma}(y)$, shear stress, $\Sigma(y)$, and first normal stress, $N_{1}(y)$; and via these to the functional form of the 0D flow curves $\Sigma(\dot{\gamma})$ and $N_{1}(\dot{\gamma})$, computed (or measured) for a state of homogeneous shear flow. In this way, we seek to give a practical guide to when an experimentalist might expect a given pressure driven flow to be stable or unstable. Note that we denote $\Sigma_{xy}=\Sigma$ and drop the 0 subscript from the base state for ease of notation. We adopt the usual definition of the first normal stress, $N_{1}=\Sigma_{xx}-\Sigma_{yy}$.

The linear analysis just described is valid as long as the perturbations remain small. To study the dynamics once the perturbations have grown (in the unstable regime) to attain a finite amplitude, we perform a full nonlinear simulation in the $x-y$ plane, using methods as in Ref. Fielding and Wilson (2010). We have checked that the early-time dynamics of this full 2D simulation, starting from a 1D base state subject to a small perturbation, agrees with the linear calculation outlined above.

The numerical calculations introduce two further parameters: the numerical timestep $dt$ and the number of numerical gridpoints $N_{y}$ (and $N_{x}$ in the fully 2D simulations). We give values for these in figure captions below, but have also checked that the results are robust to further decreases in $dt$ and increases in $N_{y},N_{x}$.

VI Results

VI.1 Rolie-Poly model

VI.1.1 0D flow curves

We start by discussing the 0D flow curves of the Rolie-Pole model, computed within the assumption of a homogeneous shear flow. (An example would be homogeneous shear between flat parallel plates that translate relative to each other, known as planar Couette flow.) The shear stress as a function of shear rate, $\Sigma(\dot{\gamma})$, is plotted in Fig. 1a). For strain rates roughly in the range $\dot{\gamma}=1-10$, this shows a regime of strong shear thinning in which the stress shows a quasi-plateau at $\Sigma^{*}\approx 0.75$, in which the stress rises only slowly with strain rate. In contrast, the first normal stress $N_{1}(\dot{\gamma})$ rises steadily to an eventual asymptote as $\dot{\gamma}\to\infty$.

VI.1.2 1D base state

We now consider the way in which the 0D flow curves just discussed, computed for a homogeneous shear flow, relate to the 1D base state profiles in a pressure driven channel flow, in which the shear rate and shear stress now vary across the channel. In many places the discussion in this section will refer to the right half of the channel, $y>0$, for definiteness. Corresponding statements apply to the left half, after taking into account the antisymmetry of $\dot{\gamma}(y)$ and symmetry of $N_{1}(y)$ about the channel centreline, $y=0$.

We start by noting that force balance requires the total shear stress to vary across the channel as $\Sigma=P^{\prime}y$, for a 1D flow. Knowledge of the imposed pressure gradient $P^{\prime}$ therefore immediately gives the shear stress as a function of the coordinate $y$ across the channel: starting from the middle of the channel at $y=0$ and moving out towards the right hand wall at $y=L_{y}/2=1/2$, the shear stress $\Sigma$ that forms the vertical axis of Fig. 1a) increases linearly from $\Sigma=0$ at $y=0$ to $\Sigma=P^{\prime}/2$ at $y=1/2$. The corresponding shear rate $\dot{\gamma}(y)$ at any $y$ can therefore be read off as the value of $\dot{\gamma}$ in Fig. 1a) corresponding to the given $\Sigma$, at that $y$. In other words, inverting the flow curve $\Sigma(\dot{\gamma})$ of Fig. 1a) to give $\dot{\gamma}=\dot{\gamma}(\Sigma=P^{\prime}y)$ immediately gives the base state shear rate profile $\dot{\gamma}(y,P^{\prime})$, for any $P^{\prime}$. A family of such curves for several values of $P^{\prime}$ is shown in Fig.  2a).

In the same way, combining the shear and normal stress flow curves $\Sigma(\dot{\gamma})$ and $N_{1}(\dot{\gamma})$ into a parametric plot $N_{1}(\Sigma)$, as shown in the inset to Fig. 1b), immediately gives the normal stress $N_{1}(\Sigma=P^{\prime}y)$ as a function of position across the channel. A family of these base state normal stress curves $N_{1}(y)$ for several values of $P^{\prime}$ is shown in Fig. 2b). A key feature in each curve (apart from at the lowest value of $P^{\prime}$) is the existence of a fairly well localised region for which $N_{1}(y)$ increases rapidly with $y$. This stems directly from the strong shear thinning seen in Fig. 1a). As the shear stress increases with $y$ across the channel, one will at the position $y^{*}=\Sigma^{*}/P^{\prime}$ corresponding to the stress quasi-plateau at $\Sigma=\Sigma^{*}\approx 0.75$ in Fig. 1 reach the shear thinning regime in which the shear rate increases strongly with $\Sigma$, and so with $y$. Correspondingly, the normal stress $N_{1}(\dot{\gamma})$ also increases strongly with $y$. This steep rise is only seen for pressure drops above about $P^{\prime}\approx 1.5$, however: for lower pressures, the stress $\Sigma=P^{\prime}y$ does not attain the shear thinning regime $\Sigma^{*}\approx 0.75$ before the wall is reached. In this way, the stress at the wall, $\Sigma_{\rm wall}=P^{\prime}/2$, determines how far up the axis of Fig. 1a) is explored.

The discussion of the previous paragraph can be expressed mathematically as:

[TABLE]

From this, we indeed see that the regime of large $dN_{1}/dy$ in the base state profile $N_{1}(y)$ in Fig. 2b) corresponds to that of strong shear thinning (small $d\Sigma/d\dot{\gamma}$) in the flow curve of Fig. 1a), for any pressure drop large enough that this shear thinning regime is indeed explored, within the window of stresses $-P^{\prime}/2<\Sigma<P^{\prime}/2$ inside the channel.

VI.1.3 Instability to 2D perturbations

We now report the results of our linear stability analysis for the dynamics of 2D perturbations to the 1D base states just discussed. Recall that we defined above $\omega_{q}$ to be the eigenvalue with the largest real part, at any wavevector $q$. The dispersion relation $\Re\omega_{q}(q)$ is shown in Fig. 3. A value $\Re\omega_{q}>0$ signifies that the 1D base state is unstable to the growth of a 2D perturbation $\propto\exp(iqx)$. For the pair of parameter values $\beta,\eta$ in Fig. 3, we indeed find instability for imposed pressure drops exceeding a critical threshold value, $P^{\prime}>P^{\prime*}(\beta,\eta)\approx 7.8$.

The corresponding first normal stress part of the eigenfunction, $\delta N_{1}(y)$ is shown in the inset of Fig. 3, for a wavevector $q=4.5$. We have checked that peaks of the eigenfunction are in essentially the same location as the derivative of the base state, $dN_{1}/dy$, and with the same properties of symmetry as $y\to-y$. This indicates that the mode of instability corresponds to a displacement of the base state. This displacement is furthermore in the same direction in both halves of the channel, corresponding to a ‘sinuous’ mode. (A ‘varicose’ mode would instead correspond to equal and opposite displacements in the two halves of the channel. Note that ‘sinuous’ vs ‘varicose’ is often defined in terms of the symmetry of the perturbation to the streamfunction, but we have checked that this accords with that of the perturbation to $N_{1}$.) This will be substantiated by the results of our nonlinear simulations in Fig. 7 below. In closely resembling the derivative of the base state, the eigenfunction appears strongly localised in the region in which the first normal stress $N_{1}(y)$ changes rapidly with $y$.

For any pair of parameter values $\beta,\eta$, there either exists a critical value $P^{\prime*}(\beta,\eta)$ of the pressure drop above which instability first arises, as in Fig. 3, or the flow is stable for all the values of imposed $P^{\prime}$ that we have explored. By performing linear stability calculations over the full range of values of the convective constraint release parameter, $0<\beta<1$, for several values of the solvent viscosity $\eta$, we find the family of curves of neutral stability $P^{\prime*}(\beta,\eta)$ shown by the solid curves in Fig. 4. For any fixed $\beta$, the flow is always stable at low $P^{\prime}$. (This is to be expected: the flow becomes Newtonian in this limit.) For the lowest value of solvent viscosity shown, $\eta=0.025$, increasing $P^{\prime}$ gives instability onset at a critical $P^{\prime*}(\beta,\eta=0.025)$ for all values of $\beta$. For larger values of $\eta$, there exists at low $\beta$ a window of pressure drops at which instability arises, with re-entrant stability at very high $P^{\prime}$; for larger $\beta$, the flow is stable at all $P^{\prime}$.

The kinks apparent in some of the solid curves in Fig. 4 arise not from numerical difficulties, but because the most unstable mode switches at the location of these kinks, with a corresponding switch of the dominant wavevector $q^{*}$.

VI.1.4 Criterion for instability onset

Within the Rolie-Poly model, we have seen that an initially 1D base state flow is linearly unstable to the onset of 2D perturbations over a wide range of pressure drops and model parameter values. Our aim now is to uncover a criterion that relates the pressure drop at which the instability first sets in (the lower boundaries in Fig. 4) to the properties of the underlying 1D base state, and thence to the properties of the 0D flow curves, thereby providing experimentalists with a practical guide for when any given pressure driven flow should become unstable.

As noted above, the eigenfunction associated with the instability is strongly localised in the region where $N_{1}(y)$ changes steeply as a function of position across the channel. Accordingly, we postulate that a high value of $dN_{1}/dy$, maximised across the channel coordinate $y$, tends to predispose the system to instability. Let us denote this by $dN_{1}/dy|_{\rm max}$.

We further recall from Eqn. 16 that $dN_{1}/dy$ increases with imposed pressure drop $P^{\prime}$. Indeed, even a non-shear thinning model such as Oldroyd B would have values of $dN_{1}/dy$ that increase indefinitely with $dP^{\prime}$ in a channel flow. We know, however, that channel flow of the Oldroyd B model is linearly stable to 2D perturbations Wilson et al. (1999); Ho and Denn (1977). Accordingly, a high value of $dN_{1}/dy$ alone cannot be a sufficient condition for linear instability. We therefore additionally postulate the presence of shear thinning to be a necessary condition, such that the region of high $dN_{1}/dy$ is focused into some localised window across the channel, as discussed above for the curves in Fig. 2b). We take as a measure of shear thinning the logarithmic slope of the flow curve, minimised across the channel. For any pressure drop exceeding $2\Sigma^{*}$, where $\Sigma^{*}$ is the stress at the quasi-plateau in Fig. 1a), this simply equals the value of $d\log\Sigma/d\log\dot{\gamma}$ minimised on the flow curve, which we accordingly denote $d\log\Sigma/d\log\dot{\gamma}|_{\rm min}$. We assume finally that the condition for instability is

[TABLE]

where $h$ is some as yet unknown increasing adimensional function of its argument. (For the inequality to be dimensionally consistent then of course requires a prefactor $G/L_{y}$ on the right hand side, which is however equal to one in our chosen units.)

Recalling that $dN_{1}/dy=P^{\prime}dN_{1}/d\Sigma$, this condition for instability can equivalently be rewritten as a criterion on the value of the imposed pressure:

[TABLE]

in which $dN_{1}/d\Sigma|_{\rm max}$ is the maximum value of $dN_{1}/d\Sigma$ at any point across the channel $y$. For pressure drops such that the shear stress at the wall exceeds the quasi-plateau value in Fig. 1a), this simply corresponds to $dN_{1}/d\Sigma$ maximised over the parametric curves of $N_{1}(\Sigma)$ in the inset to Fig. 1b). We have confirmed, further, that this corresponds closely to the value of $dN_{1}/d\Sigma$ as calculated at the shear rate corresponding to the point of minimum logarithmic slope of the constitutive curve. Accordingly, we now use $dN_{1}/d\Sigma|_{\rm max}$ to denote this value.

To test this criterion, we plot in Fig. 5 the left hand side of 18 as a function of the argument of $h$ on the right hand side. To do so, for each parameter pairing $\beta,\eta$ explored in Fig. 4 we read off the critical pressure drop from that Fig. 4, and compute $d\log\Sigma/d\log\dot{\gamma}$ minimised over the flow curve of Fig. 1a) for that $\beta,\eta$, together with $dN_{1}/d\Sigma|_{\rm max}$ as just defined (corresponding closely to the value of this quantity maximised over the parametric normal stress curve inset in Fig. 1b). The results are collected in Fig. 5 into a set of curves, with each curve containing all the data for a single value of $\beta$. We obtain reasonable data collapse onto the function

[TABLE]

with the same set of fitting parameters $\alpha_{0}=25.81$, $\alpha_{1}=0.332$ and $\alpha_{2}=0.205$ for all $\beta,\eta$.

Now that this function $h(x)$ is known, we can finally reconstruct the critical pressure drop for any pair of parameter values $\beta,\eta$ by computing $h(x)/y$ with $x=d\log\Sigma/d\log\dot{\gamma}_{\rm min}$ (for that $\beta,\eta$) and $y=dN_{1}/d\Sigma|_{\rm max}$ (for that $\beta,\eta$). This gives the dashed lines in Fig. 4. These indeed fit the numerical data quite well at low values of the imposed pressure, thereby capturing the initial onset of instability. The fits however fail to capture the regime of re-entrant stability, and so depart from the numerical data around the ‘nose’ in that data.

The form of the function $h$ reveals that instability can never be attained for a value of the logarithmic derivative of the flow curve, minimised across the flow curve, $d\log\Sigma/d\log\dot{\gamma}|_{\rm min}$, that exceeds a value roughly equal to $0.2$. This sets a basic condition on how flat the logarithmically plotted flow curve must be to admit instability: were the flow curve to be approximated by a power law $\Sigma\sim\dot{\gamma}^{n}$ over the quasi-plateau region, instability requires $n<0.21$, in the Rolie-Poly model. This is broadly consistent with instability having been observed for the shear thinning exponents of $n=0.21$ and $n=0.19$ in the experimental studies of Refs. Bodiguel et al. (2015) and Poole (2016) respectively, but apparently at odds with instability having been observed for $n=0.30$ in the experiments of Ref. Picaut et al. (2017).

Indeed, in Fig. 6 we represent the value of $d\log\Sigma/d\log\dot{\gamma}|_{\rm min}$ by dotted contour lines in the plane of the model parameters $\eta$ and $\beta$. We also show by a cross the value of solvent viscosity below which we find instability, at any pressure drop, for any value of $\beta$. The regime of instability is indeed found to be that for which the minimum logarithmic derivative is less than about $0.21$.

The thick black line in Fig. 6 shows the contour $n=0$, below which the constitutive curve of stress against strain rate is non-monotonic, admitting shear banding. Previous studies Fielding (2005); Fielding and Wilson (2010); Fielding and Olmsted (2006) in the Johnson-Segalman model showed an initially 1D shear banded state to be linearly unstable to the onset of 2D perturbations in which the interface between the bands destabilises. While we do not explore this region in detail here, we have confirmed that the interface is indeed unstable in the Rolie-Poly model, increasing the generality of the previous finding in the Johnson-Segalman model. (To extend our study into the banding regime, diffusive terms were added to the right hand side of the Rolie-Poly constitutive equation Lu et al. (2000).) In this way, the instability of a strongly shear thinning – but not shear banding – fluid that is the main focus of this work can potentially be understood as an instability of the quasi-interface in the region of steeply increasing $N_{1}(y)$ in Fig. 2b).The addition of a stress diffusion term outside the banding regime would have the effect of reducing the steepness of variation of $N_{1}(y)$ within the channel, thereby potentially having a slightly stabilising effect on the overall flow Wilson and Rallison (1997).

VI.1.5 Non-linear dynamics

The linear analysis just described is valid in the regime where the 2D perturbations to the 1D base state remain small. To study the dynamics once the perturbations have grown (in the unstable regime) to attain a finite amplitude, we performed fully nonlinear 2D simulations. Sample snapshots of the ultimate state are shown in Fig. 7. The left panel confirms the basic mode of instability to be one of displacement of a quasi-interface between regions of low and high first normal stress. The right panel shows the accompanying difference between the velocity fields of the final 2D and initial 1D velocity fields, revealing a system-spanning vortex (which is modest in amplitude compared with the velocities in the underlying base state). We note that the amplitude of velocity perturbation seen here in our nonlinear simulations is smaller compared with the base flow, than in the experiments of Ref. Bodiguel et al. (2015); Poole (2016).

VI.2 Johnson-Segalman model

VI.2.1 0D flow curves

The flow curves of the Johnson-Segalman model, describing states of homogeneous shear flow, are shown in Fig. 8. These have the same basic form as for the Rolie-Poly model in Fig. 1, with a stress quasi-plateau associated with a regime of strong shear thinning. An important difference between the two models, however, is that whereas the stress $\Sigma^{*}$ at the quasi-plateau in the Rolie-Poly model is more or less independent of the CCR parameter, $\beta$, in the Johnson-Segalman model it depends strongly on the slip parameter $a$, scaling as:

[TABLE]

This diverges in the limit $a\to 1$, in which the Johnson-Segalman model becomes the Oldroyd B model.

A second difference between the Rolie-Poly and Johnson-Segalman models is that the minimum logarithmic slope of the flow curve, $d\log\Sigma/d\log\dot{\gamma}|_{\rm min}$, is independent of the slip parameter $a$ in the Johnson-Segalman model, whereas it depends strongly on the CCR parameter $\beta$ in the Rolie-Poly model. (Recall Fig. 6.)

VI.2.2 1D base states

The 1D base states $\dot{\gamma}(y)$ and $N_{1}(y)$ of the Johnson-Segalman model are shown in Fig. 9. These can be related to the 0D flow curves of Fig. 8 in the same way as for the Rolie-Poly model. In particular, the regime of strong shear thinning in the flow curve $\Sigma(\dot{\gamma})$ gives rise to the region of rapid variation of shear rate $\dot{\gamma}(y)$ and first normal stress $N_{1}(y)$ as a function of the position $y$ across the channel.

VI.2.3 Instability to 2D perturbations

We now consider whether the 1D base states just discussed are stable or unstable to the growth of 2D perturbations, for any set of model parameter values $a,\eta$, channel length $L_{x}$ and imposed pressure drop. For a particular pair of values of $a,\eta$, we show in Fig. 10 the dispersion relation of the real part of the most unstable eigenvalue as a function of wavevector $q$ (which recall must be quantised as $n\pi/L_{x}$ for any channel of finite length $L_{x}$.) As in the Rolie-Poly model, we find onset of instability above a critical pressure drop $P^{\prime*}(\eta,a)$. The dispersion relations have very similar form to those in the Rolie-Poly model. (Recall Fig. 3.) So too does the associated first normal stress part of the eigenfunction $\delta N_{1}(y)$, as shown in the inset to Fig. 10. Indeed, it is heavily localised in the region of steep gradient in $N_{1}(y)$ of the underlying base state. The instability therefore again appears to correspond to a displacement of the quasi-interface formed from the region of steep variation in $N_{1}(y)$ across the channel in the underlying base state. This will be confirmed by the results of our nonlinear simulations below.

This is consistent with earlier studies of the Johnson-Segalman model in the shear banding regime, $\eta<0.125$, in which a 1D shear banding flow was shown to be unstable to 2D undulations along the interface between the bands Fielding (2005); Fielding and Wilson (2010); Fielding and Olmsted (2006), as noted above. However, we do not consider this shear banding regime here, confining ourselves instead to values of the solvent viscosity $\eta>0.125$ for which the constitutive curve is monotonic.

By performing linear stability calculations over a range of values of the slip parameter, $0<a<1$, for several values of the solvent viscosity $\eta$, we find a family of curves of neutral stability $P^{\prime*}(a,\eta)$ in Fig. 11. For any fixed $a$, the flow is always stable at low $P^{\prime}$. (This is to be expected: the flow becomes Newtonian in this limit.) For the lowest four values of solvent viscosity explored in Fig. 11, we find a critical value of pressure drop $P^{\prime*}$ above which instability arises, at all values of $a$. For the highest value of viscosity explored, there is a restricted window of values of $a$, above about $0.7$, that admit instability.

We also explored much higher values of imposed pressure drops, up to $P^{\prime*}=60.0$, for the single value of solvent viscosity $\eta=0.16$. The results (not shown) reveal a regime of re-entrant stability above a pressure drop of, for example, about $35$ for $a=0.0$, and about $55$ for $a=0.8$. We do not explore this further here: such pressure drops may in any case be unattainable in practice.

VI.2.4 Criterion for onset of instability

We now seek to determine whether the criterion (18) developed above in the Rolie-Poly model, relating the pressure drop at which the first instability sets in to the properties of the underlying 1D base state, and thence to the the 0D flow curves, also applies in the Johnson-Segalman model. Recall that this criterion reads:

[TABLE]

with

[TABLE]

in which $d\log\Sigma/d\log\dot{\gamma}_{\rm min}$ refers (for each $a,\eta$) to the minimum logarithmic derivative of the flow curve in Fig. 8a), and $dN_{1}/d\Sigma|_{\rm max}$ refers (for each $a,\eta$) to the maximum derivative of the parametrically plotted normal stress curve in the inset to Fig. 8b). (In fact, it refers to the value of this derivative computed at a shear rate corresponding to the point of minimuum logarithmic slope of the flow curve. This coincides very nearly with the maximum value of $dN_{1}/d\Sigma$.)

Plotting the left hand side of 21 as a function of the argument of $h$ on the right hand side, we find a reasonable fit, now with $\alpha_{0}=60.0,\alpha_{1}=1.15$ and $\alpha_{2}=0.11$. This gives the fitted onset boundaries shown by dashed lines in Fig. 11, for values of $a$ less than around $0.6$. As can be seen, it performs reasonably well for the lowest values of solvent viscosity explored, where the pressure onset values are modest. It performs less well, however, for higher solvent viscosities, for which the pressure onset values are much larger. This could be due to the fact that the quasi-interfacial region of steep $N_{1}(y)$ becomes close (on the scale of its own thickness) to the centre of the channel at such large pressure drops: a feature not taken into account in our analysis. However, we ignore this complication, on the grounds that such pressure drops may be difficult to obtain in practice in any case.

The value $\alpha_{2}=0.11$ obtained in the above fitting shows that the value of the logarithmic slope of the flow curve in the quasi-plateau thinning regime should not exceed $n=0.11$ in any fluid for instability to be observed. This sets a basic requirement on how steep the (shallowest part of the) flow curve must be to avoid instability, as explored further in Fig. 12. The value $0.11$ in the Johnson-Segalman model shows that this model predicts, overall, greater stability than the Rolie-Poly model: a fluid has to shear thin more strongly to be predicted unstable by Johnson-Segalman than by Rolie-Poly.

Whereas the procedure just discussed proved successful in the Rolie-Poly model for all values of the CCR parameter $\beta$, in the Johnson-Segalman model it is restricted to values of the slip parameter $a$ that are not too close to $1.0$. (We used only $0<a<0.6$ in the fitting just described.) For values of $a$ close to $1.0$, the height of the stress quasi-plateau in the flow curves diverges, with the Johnson-Segalman model tending to the Oldroyd B model in the limit $a\to 1$: recall Fig. 8a). This brings an additional requirement for instability: that the imposed pressure drop is large enough that the range of shear stresses within the channel, $0<|\Sigma|<P^{\prime}L_{y}/2$ encompasses the value of the stress quasi-plateau $\Sigma^{*}\sim G/\sqrt{1-a^{2}}$. This is necessary to ensure that the region in which $N_{1}(y)$ varies rapidly with $y$ is indeed present within the channel. Accordingly, we update our criterion to now be a double requirement: both that  21 is satisfied, and that

[TABLE]

It is this double criterion that is represented by the dashed lines in Fig. 11, with condition 21 being more stringent for lower values of $a$, and criterion 23 more stringent for higher values of $a$. (In the Rolie-Poly model, we did not need to state the additional condition 23 explicitly, because it is automatically satisfied in all regimes where  18 also holds.)

VI.2.5 Non-linear dynamics

To study the dynamics once the perturbations have grown to attain a finite amplitude, we performed fully nonlinear 2D simulations. Sample snapshots of the ultimate state are shown in Fig. 13. The left panel confirms the basic mode of instability to be one of displacement of a quasi-interface between regions of low and high first normal stress, as in the Rolie-Poly model. (Compare Fig. 7.) The right panel shows the accompanying difference between the velocity fields of the final 2D and initial 1D velocity fields, again revealing a system-spanning vortex.

VI.3 White-Metzner model

VI.3.1 0D flow curves

The flow curves of the White-Metzner model can be calculated analytically. In our usual units in which the modulus $G=1$ and the basic timescale $\tau=1$, one has $\Sigma=\dot{\gamma}|\dot{\gamma}|^{n-1}+\eta\dot{\gamma}$ and $N_{1}=\dot{\gamma}^{2n}$. An important point to note here is that, although the flow curve does display shear thinning in this model, it does so with a constant logarithmic slope, $d\log\Sigma/d\log\dot{\gamma}=n$ (at least up to small corrections arising from the solvent viscosity). Accordingly, the region of shear thinning is not focused into a narrow stress window in which the flow curve shows a quasi-plateau, as for the Rolie-Poly model in Fig. 1 and the Johnson-Segalman model in Fig. 8.

VI.3.2 1D base states

Combining the form of the White-Metzner model’s flow curve $\Sigma=\dot{\gamma}|\dot{\gamma}|^{n-1}$ (in which we ignore small corrections in due to the solvent viscosity) with the fact that the shear stress varies as function of position across the channel as $\Sigma=P^{\prime}y$, enables us to compute the shear rate as a function of position across the channel in the 1D base state: $\dot{\gamma}=\operatorname{sign}(y)(P^{\prime}|y|)^{1/n}$. Combining this with the form of the normal stress as a function of strain rate, $N_{1}=\dot{\gamma}^{2n}$, then gives the normal stress as a function of position across the channel, $N_{1}=(P^{\prime}y)^{2}$ (as long as $\eta\dot{\gamma}$ remains small). An important consequence of the absence of any quasi-plateau in the White-Metzner model’s flow curve is that there is no quasi-interface across which the first normal stress varies rapidly as a function of position across the channel, as there is for the Rolie-Poly model in Fig. 2 and for the Johnson-Segalman model in Fig. 9. In consequence, we expect the basic physics of any instability to differ from the interfacial mode found for Rolie-Poly and Johnson-Segalman.

VI.3.3 Instability to 2D perturbations

Instability has previously been reported in the White-Metzner model by Wilson et al. Wilson and Rallison (1999); Castillo and Wilson (2017); Wilson and Loridan (2015). The purpose of our studying the same model again in this work is to examine whether the results can be cast into the same basic format as for the Rolie-Poly and Johnson-Segalman models above, and thereby to try to establish the extent to which any instability of shear thinning channel flow is generic across constitutive models.

We therefore now consider whether the 1D base states just discussed are stable or unstable to the onset of 2D perturbations, for any set of model parameter values $n,\eta$, channel length $L_{x}$ and imposed pressure drop. For a particular pair of values of $n,\eta$, we show in Fig. 14 the dispersion relation of the real part of the most unstable eigenvalue as a function of wavevector $q$ (which recall must be quantised as $n\pi/L_{x}$ for any channel of finite length $L_{x}$). As in the Rolie-Poly and Johnson-Segalman models, we find onset of instability above a critical pressure drop $P^{\prime*}(\eta,n)$. The dispersion relations have broadly similar form to those for the Rolie-Poly and Johnson-Segalman models in Figs. 3 and 10.

The associated first normal stress part of the eigenfunction, $\delta N_{1}(y)$, is shown in the inset to Fig. 14. A crucial difference between the White-Metzner model and the other two models is that the eigenfunction is diffusely spread across the channel in White-Metzner, rather than being heavily localised, as in Rolie-Poly and Johnson-Segalman. This is consistent with the absence of any steep gradient in $N_{1}(y)$ in the underlying base state. The instability in White-Metzner therefore does not appear to correspond to a displacement of any quasi-interface formed from a region of steep variation in $N_{1}(y)$ across the channel in the underlying base state.

By performing linear stability calculations over a range of values of the shear thinning parameter, $0<n<1$, for several values of the solvent viscosity in the physically relevant regime $\eta<G\tau=1$, we find a family of curves of neutral stability $P^{\prime*}(n,\eta)$ in Fig. 15.

VI.3.4 Criterion for onset of instability

In the Rolie-Poly and Johnson-Segalman models, we constructed a criterion for the critical pressure drop needed to observe instability. Recall that this reads:

[TABLE]

with

[TABLE]

in which $d\log\Sigma/d\log\dot{\gamma}|_{\rm min}$ refers (for each $\xi,\eta$) to the minimum logarithmic derivative of the flow curve, and $dN_{1}/d\Sigma|_{\rm max}$ refers (for each $\xi,\eta$) to the maximum derivative of the parametrically plotted normal stress as a function of shear stress. In the White-Metzner model, $d\log\Sigma/d\log\dot{\gamma}=n$ is constant across the flow curve (to within small corrections set by the solvent viscosity). Accordingly, we drop the subscript “min” from $d\log\Sigma/d\log\dot{\gamma}|_{\rm min}$. As noted above, this absence of any quasi-plateau in the flow-curve means that there is no strongly localised region of high gradient in $dN_{1}/dy$ forming a quasi-interface at some location across the channel. Instead, $dN_{1}/d\Sigma$ is now maximum close to the channel wall, and it is at this location that we accordingly now calculate $dN_{1}/d\Sigma|_{\rm max}$. This derivative therefore cannot be calculated from a one-off inspection of the flow curve in White-Metnzer, but must involve explicit knowledge of the stress at the wall of the channel.

The criterion just discussed was built on arguments concerning the existence of a quasi-interface in $N_{1}(y)$ in the underlying base state of 1D channel flow, in the Rolie-Poly and Johnson-Segalman models. This stemmed in turn from the quasi-plateau in the models’ basic flow curves $\Sigma(\dot{\gamma})$ for homogeneous shear flow. In view of the absence of any quasi-plateau (and corresponding quasi-interface) in the White-Metzner model, there is no reason, a priori, to expect the criterion still to apply. Nonetheless, we explore whether it does by plotting $dN_{1}/d\Sigma|_{\rm max}$ as a function of $d\log(\Sigma)/d\log(\dot{\gamma})=n$ in Fig. 16. As can be seen, a good fit to the function $h$ is obtained, with $\alpha_{0}=1.21,\alpha_{1}=0.55$ and $\alpha_{2}=0.31$. This then enables us to reconstuct a fit to the critical pressure drop as a function of $\eta,n$, as shown by the dashed line in Fig. 15, giving reasonable agreement with the numerical data. Finally, we show in Fig. 17 contours of the logarithmic slope of the flow curve (minimised over the flow curve) in the plane of $\eta,n$. Also shown by crosses are the values of $\eta,n$ that mark the boundary between the region of stability (for all values of $P^{\prime}$) at high $n$ from that of instability (for some values of $P^{\prime}$) at low $n$.

As is evident from Figs. 15 to  17, the results are largely independent of the solvent viscosity $\eta$ in the White-Metzner model. This should not be surprising: the solvent contribution to the flow curve stress, $\Sigma=\dot{\gamma}|\dot{\gamma}|^{n-1}+\eta\dot{\gamma}$, is much smaller relative to the viscoelastic one, for realistic values of $\eta<G\tau=1$, at least in the window of stresses explored across the channel in the vicinity of instability onset. In this way, the collapse obtained with $\eta$ in Figs. 15 and 16 is relatively trivial. The more interesting feature of these results is the fact that the instability threshold appears to follow the same functional dependency on $n$, albeit with different values of the fitting parameters $\alpha_{0},\alpha_{1},\alpha_{2}$.

VII Conclusions

We have combined linear stability analysis with full numerical simulation to study pressure driven channel flow within three widely used constitutive models of shear thinning viscoelastic flow: the microscopically motivated Rolie-Poly model, and the phenomenological Johnson-Segalman and White-Metzner models. We have shown an initially 1D base state to be linearly unstable to the onset of 2D perturbations with wavevector in the flow direction, in all three models, for a high enough degree of shear thinning. That we indeed find instability across several constitutive models suggests that the instability may be generic across shear thinning polymeric fluids. Within each model, we have calculated the minimal degree of shear thinning needed to observe instability, finding instability below a critical value $n^{*}$ of the logarithmic slope of the flow curve at its shallowest point, $n=d\log\Sigma/d\log\dot{\gamma}|_{\rm min}$, with $n^{*}\approx 0.2$ in Rolie-Poly, $n^{*}\approx 0.11$ in Johnson-Segalman, and $n^{*}\approx 0.3$ in White-Metzner.

Within each of the Rolie-Poly, Johnson-Segalman and White-Metzner models, we have shown the critical adimensional pressure drop at the onset of instability to follow, to a reasonable level of approximation, a scaling function expressed in terms of (i) this degree of shear thinning, $n$, and (ii) the maximum derivative of the the first normal stress with respect to the shear stress. In each model, the scaling function is of the same functional form. However, it is characterised by three fitting parameters, $\alpha_{0},\alpha_{1},\alpha_{2}$, the values of which differ across the three models. In this way, the onset criterion obtained here is clearly only a partial success. An open challenge for future work is to recast these fitting parameters in terms of any additional relevant dimensionless constitutive properties, leading to a re-expressed set of fitting parameters that have the same values across all models, thereby providing a truly universal criterion for the onset of instability. This goal might be unattainable, however.

In the Rolie-Poly and Johnson-Segalman models, we have given evidence suggesting that the mechanism of instability involves the presence of a quasi-interface at some location in each half of the channel, across which the first normal stress and shear rate both vary steeply. In this way, the instability appears similar to that arising at the interface between layered fluids or shear bands. In the White-Metzner model, no such quasi-interface exists. An unresolved puzzle is why the same criterion for instability (albeit with different values of the fitting parameters) appears to hold in the White-Metzner model as in the Rolie-Poly and Johnson-Segalman models.

Recently, the original predictions of instability within the White-Metzner model Wilson and Rallison (1999) were extended by including an additional stress contribution that matches the shear thinning of White-Metzner, but that responds instantaneously to the imposed flow, and lacks any normal stress contribution Castillo and Wilson (2017). As the contribution from this term increased relative to the viscoelastic contribution from White-Metzner, the level of instability decreased and stability was eventually restored. The overall scaling of the modulus and relaxation time of the White-Metzner model were also respectively decreased and increased with the increasing additional stress contribution. The decrease of the modulus, in particular, decreases the normal stress. This work appears consistent with our finding, of increasing levels of stability with decreasing levels of normal stress. In future work, it would be interesting to see whether the predictions of Castillo and Wilson (2017) indeed accord with our criterion.

As things currently stand, then, two possible mechanisms have been identified for instability: a jump in $N_{1}$ across an interface between two fluid layers Wilson and Rallison (1997), and an instability of shear thinning flow, even if $N_{1}$ varies smoothly across the channel Wilson and Rallison (1999). This raises the question of whether, in the various models considered here, the instabilities are related primarily to rapid variations in $N_{1}$ or in $\dot{\gamma}$ across the channel. The answer to this remains unclear at present, but it is worth collecting the evidence for the spatial nature of the eigenfunctions in the different models. In the Rolie-Poly and Johnson-Segalman models, both the viscosity and $N_{1}$ vary rapidly in the quasi-interface region, and the eigenfunction is localised here. In the White-Metzner model Wilson and Rallison (1999), the viscosity and $N_{1}$ both vary as power laws across the channel, and the eigenfunction appears delocalised right across the channel. In a very recent study Castillo and Wilson (2018) of a model in which $N_{1}$ varies smoothly across the channel but the viscosity varies relatively rapidly near the walls, the eigenfunction was found to relatively localised near the walls. Refs Miller and Rallison (2007a, b) considered layered viscoelastic fluids with matched viscosities, but different relaxation times and so normal stresses. The interface was found to be unstable, even with the addition of a surface tension. As things stand, therefore, it remains unclear whether instabilities of the kind studied in this work are generically driven primarily by variations across the channel in the base state viscosity, in the base state normal stress, or even in some dynamical property not captured by the stationary underlying base state.

We now make some comments about the Giesekus model, for which instability was previously found in Ref. Grillet et al. (2002), and for which we have performed here a more modest range of numerical explorations. The stability behaviour of this model is rather different from that of the other three models. In particular, the maximum logarithmic slope $n^{*}$ of the flow curve (minimised across the flow curve) that permits instability is much smaller in this model, and furthermore depends strongly on the solvent viscosity. For values of $\eta$ that we explored in the range [math] to $0.015$, $n^{*}$ appears to vary in the range [math] - $0.04$. Therefore, a fluid must be much more highly shear thinning to be predicted unstable by Giesekus than by the other three models. Furthermore, for any set of values of $\alpha,\eta$, the window of pressure drops for which instability arises (if at all) is very much narrower than in the other models, with re-entrant stability following at a pressure drop not much greater than that for instability onset. Taken overall, then, the Giesekus model appears to predict a much greater level stability in channel flow. It is probably just a coincidence, although worth remarking, that the same model also predicted much greater levels of stability against the formation of transient bands in shear startup (compared with the Rolie-Poly model) in Ref. Moorcroft and Fielding (2014).

An aspect of the Giesekus model noted in Ref.Giesekus (1982) is that for a non-linear relaxation parameter $\alpha=1$, the shear and normal stress constitutive curves exactly match those of the co-rotational Maxwell model, which is equivalent to Johnson Segalman for $a=0$. In this way, the 1D base state of the Giesekus model at $\alpha=1$ is the same as that of the Johnson-Segalman model at $a=0$. Puzzlingly, the former appears stable, within Giesekus dynamics (data not shown), while the latter appears unstable, within the Johnson-Segalman model (recall Fig. 11). This observation again suggests that some feature of constitutive dynamics not represented by the stationary base state plays a role in determining instability.

All the constitutive models used in this study are designed to describe shear thinning viscoelastic fluids. Given that they predict rather different stability properties from each other (quantitatively for the Rolie-Poly, Johnson-Segalman and White-Metzner model and even qualitatively for the Giesekus model), this suggests that prediction of stability properties could be used as a helpful constraint in model building, beyond the fitting to homogeneous rheological functions as is more usually carried out.

Finally, it is worth remarking that we have considered here a time-independent and x-independent 1D base state, about which 2D perturbations may or may not grow. In practice, one may instead be interested in whether a 1D state of fluid flowing into and out of a channel of a finite length becomes unstable: in particular, in whether the instability has time to develop for a channel of any given length. To answer that question fully, one would have to simulate that more dynamic situation. A reasonable estimate may, however, be gained from our results by comparing the timescale $\omega^{-1}$ for instability to develop with the residence time in the channel, $L_{x}/V$, where $V$ is the velocity at the channel midpoint.

Acknowledgements — We thank Mike Graham, Hugo Castillo-Sanchez and Helen Wilson for useful discussions. We also acknowledge the SOFI CDT, Durham University, EPSRC (EP/L015536/1), the European Research Council under the European Union Seventh Framework Programme (FP7/2007-2013)/ ERC grant agreement number 279365, and Schlumberger Cambridge Research for funding.