Edge fracture instability in sheared complex fluids: onset criterion and possible mitigation strategy

TL;DR

This paper provides a comprehensive theoretical analysis of edge fracture instability in sheared complex fluids, deriving an exact onset criterion, mapping phase diagrams, and proposing mitigation strategies based on nonlinear simulations and linear stability calculations.

Contribution

It introduces an analytical expression for the onset of edge fracture and validates it with nonlinear simulations, offering new insights into instability mechanisms and mitigation methods.

Findings

Edge fracture threshold is robust against wetting angle variations.

Analytical onset criterion matches nonlinear simulation results.

Proposed mitigation strategies could improve rheological measurements.

Abstract

We perform a detailed theoretical study of the edge fracture instability, which commonly destabilises the fluid-air interface during strong shear flows of entangled polymeric fluids, leading to unreliable rheological measurements. By means of direct nonlinear simulations, we map out phase diagrams showing the degree of edge fracture in the plane of the surface tension of the fluid-air interface and the imposed shear rate, within the Giesekus and Johnson-Segalman models, for different values of the nonlinear constitutive parameters that determine the dependencies on shear rate of the shear and normal stresses. The threshold for the onset of edge fracture is shown to be relatively robust against variations in the wetting angle where the fluid-air interface meets the hard walls of the flow cell, whereas the nonlinear dynamics depend strongly on wetting angle. We perform a linear stability calculation to derive an exact analytical expression for the onset of edge fracture, expressed in terms of the shear-rate derivative of the second normal stress difference, the shear-rate derivative of the shear stress (sometimes called the tangent viscosity), the jump in shear stress across the interface between the fluid and the outside air, the surface tension of that interface, and the rheometer gap size. Full agreement between our analytical calculation and nonlinear simulations is demonstrated. We also elucidate in detail the mechanism of edge fracture, and finally suggest a new way in which it might be mitigated in experimental practice. Some of the results in this paper were first announced in an earlier letter. The present manuscript provides additional simulation results, calculational details of the linear stability analysis, and more detailed discussion of the significance and limitations of our findings.

Tables

Table 1. Table 1: Parameters and their dimensions in terms of modulus [ G ] delimited-[] 𝐺 [G] , length [ L ] delimited-[] 𝐿 [L] and time [ T ] delimited-[] 𝑇 [T] ; values used in our numerical simulations; and notes concerning each parameter. The first three parameters listed specify our choice of units. The second five are the key physical parameters to be varied in our study (four within each constitutive model); the set from L z subscript 𝐿 𝑧 L_{z} to M 𝑀 M do not affect the key physics, provided each is set to an appropriately large or small value; the final set are numerical parameters and are converged to their appropriate large or small limit. Abbreviations: Johnson-Segalman (JS) and Giesekus (Gk).

parameter	description	dimension	value	notes
$L_y$	channel width	$\left[L\right]$	1.0	unit of length
$G$	polymer modulus	$\left[G\right]$	1.0	unit of stress
$\tau$	polymer relaxation time	$\left[T\right]$	1.0	unit of time
$\dot{\gamma }$	applied shear-rate	${\left[T\right]}^{-1}$	${10}^{-1}\to {10}^2$	important quantity to be varied
$\theta$	equilibrium contact angle	$\left[1\right]$	${60}^{\circ }\to {120}^{\circ }$	important quantity to be varied
$\Gamma$	surface tension	$\left[G\right]\left[L\right]$	$0.0\to 1.0$	important quantity to be varied
$a$	slip parameter (JS)	$\left[1\right]$	$0.0\to 0.6$ (usually 0.3)	sets dependencies ${\Sigma }_{xy}\left(\dot{\gamma }\right),N_2\left(\dot{\gamma }\right)$
$\alpha$	anisotropy parameter (Gk)	$\left[1\right]$	$0.1\to 1.0$ (usually 0.4)	sets dependencies ${\Sigma }_{xy}\left(\dot{\gamma }\right),N_2\left(\dot{\gamma }\right)$
$L_z$	channel length	$\left[L\right]$	$10.0$	large aspect ratio $L_z/L_y$
$\Lambda$	midpoint sample length	$\left[L\right]$	$7.0$	large enough air gap $\left(L_z-\Lambda \right)/L_y$
${\eta }_a$	air viscosity	$\left[G\right]\left[T\right]$	$0.01$	small air viscosity ${\eta }_a/G\tau$
${\eta }_s$	solvent viscosity (JS)	$\left[G\right]\left[T\right]$	0.15	small viscosity ratio ${\eta }_s/G\tau$
${\eta }_s$	solvent viscosity (Gk)	$\left[G\right]\left[T\right]$	0.01 ($\alpha \le 0.5$); 0.15 ($\alpha >0.5$)	small viscosity ratio ${\eta }_s/G\tau$
$\ell$	polymer microscopic length	$\left[L\right]$	0.01	small microscopic length $l/L_y$
${\ell }_{\mu }$	air-polymer interface width	$\left[L\right]$	0.01 (unless stated otherwise)	small microscopic length ${\ell }_{\mu }/L_y$
$M$	molecular mobility	${\left[L\right]}^2{\left[G\right]}^{-1}{\left[T\right]}^{-1}$	0.0001	rapid phase equilibration
$N_y$	numerical mesh	$\left[1\right]$	large	converge until no dependence
$N_z$	numerical mesh	$\left[1\right]$	large	converge until no dependence
$Dt$	numerical timestep	$\left[1\right]$	small	converge until no dependence

Equations

$$\nabla .\mathbf{T}=0.$$

$${\eta }_{\mathrm{s}}{\nabla }^2\mathbf{v}+\nabla .\mathbf{\Sigma }-\nabla p=0,$$

$${\eta }_{\mathrm{a}}{\nabla }^2\mathbf{v}-\nabla p=0$$

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

$${\partial }_t\mathbf{\Sigma }+\mathbf{v}.\nabla \mathbf{\Sigma }=\left(\mathbf{\Sigma }\cdot \mathbf{\Omega }-\mathbf{\Omega }\cdot \mathbf{\Sigma }\right)+a\left(\mathbf{D}\cdot \mathbf{\Sigma }+\mathbf{\Sigma }\cdot \mathbf{D}\right)+2G\mathbf{D}-\frac{1}{\tau }\mathbf{\Sigma }+\frac{{\ell }^2}{\tau }{\nabla }^2\mathbf{\Sigma }.$$

$${\Sigma }_{xy}$$

$$N_2\equiv {\Sigma }_{yy}-{\Sigma }_{zz}$$

$${\partial }_t\mathbf{\Sigma }+\mathbf{v}.\nabla \mathbf{\Sigma }=\left(\mathbf{\Sigma }\cdot \mathbf{\Omega }-\mathbf{\Omega }\cdot \mathbf{\Sigma }\right)+\left(\mathbf{D}\cdot \mathbf{\Sigma }+\mathbf{\Sigma }\cdot \mathbf{D}\right)+2G\mathbf{D}-\frac{1}{\tau }\mathbf{\Sigma }-\frac{\alpha }{\tau }\mathbf{\Sigma }\cdot \mathbf{\Sigma }+\frac{{\ell }^2}{\tau }{\nabla }^2\mathbf{\Sigma }.$$

$${\Sigma }_{xy}$$

$$N_2$$

$${\Lambda }^2=\frac{1}{8\alpha (1-\alpha ){\dot{\gamma }}^2}\left[\sqrt{1+16\alpha (1-\alpha ){\dot{\gamma }}^2}-1\right].$$

$${\partial }_t\varphi +\mathbf{v}.\nabla \varphi$$

$$\mu$$

$$\Gamma =\frac{2\sqrt{2}}{3}{G}_{\mu }{\ell }_{\mu }.$$

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

$$\mathbf{n}\cdot \nabla \mu$$

$$\mathbf{n}\cdot \nabla \varphi$$

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

$$\nabla \cdot \mathbf{T}=0,$$

$$\mathbf{T}=2{\eta }_{\mathrm{s}}\mathbf{D}+\mathbf{\Sigma }-p\mathbf{I}.$$

$${\partial }_t\mathbf{\Sigma }+\mathbf{v}.\nabla \mathbf{\Sigma }=\left(\mathbf{\Sigma }\cdot \mathbf{\Omega }-\mathbf{\Omega }\cdot \mathbf{\Sigma }\right)+a\left(\mathbf{D}\cdot \mathbf{\Sigma }+\mathbf{\Sigma }\cdot \mathbf{D}\right)+2G\mathbf{D}-\frac{1}{\tau }\mathbf{\Sigma }-\frac{\alpha }{\tau G}\mathbf{\Sigma }\cdot \mathbf{\Sigma },$$

$$\mathbf{\Sigma }=2\mathbf{D}G\tau +2(\mathbf{D}\cdot \mathbf{\Omega }-\mathbf{\Omega }\cdot \mathbf{D})G{\tau }^2+2(a-\alpha )(\mathbf{D}\cdot \mathbf{D}+\mathbf{D}\cdot \mathbf{D})G{\tau }^2+h.o.t.$$

$$\mathbf{n}\cdot \mathbf{T}+\Gamma \mathbf{n}{\nabla }_{\mathrm{int}}\cdot \mathbf{n}=0,$$

$${\partial }_{t}h=v_z.$$

$$\mathbf{v}={\mathbf{v}}_0+\tilde{\mathbf{v}}=\left(\begin{array}{c}{\dot{\gamma }}_{0}y\\ 0\\ 0\end{array}\right)+\left(\begin{array}{c}{\tilde{v}}_x\\ {\partial }_z\tilde{\psi }\\ -{\partial }_y\tilde{\psi }\end{array}\right),$$

$$\nabla \mathbf{v}=\nabla {\mathbf{v}}_0+\nabla \tilde{\mathbf{v}}=\left(\begin{array}{ccc}0& 0& 0\\ {\dot{\gamma }}_0& 0& 0\\ 0& 0& 0\end{array}\right)+\left(\begin{array}{ccc}0& 0& 0\\ {\partial }_y{\tilde{v}}_x& {\partial }_y{\partial }_z\tilde{\psi }& -{\partial }_y^2\tilde{\psi }\\ {\partial }_z{\tilde{v}}_x& {\partial }_z^2\tilde{\psi }& -{\partial }_y{\partial }_z\tilde{\psi }\end{array}\right).$$

$$\mathbf{\Sigma }={\mathbf{\Sigma }}_0+\tilde{\mathbf{\Sigma }}=\left(\begin{array}{ccc}{\Sigma }_{xx0}({\dot{\gamma }}_0)& {\Sigma }_{xy0}({\dot{\gamma }}_0)& 0\\ {\Sigma }_{xy0}({\dot{\gamma }}_0)& {\Sigma }_{yy0}({\dot{\gamma }}_0)& 0\\ 0& 0& {\Sigma }_{zz0}({\dot{\gamma }}_0)\end{array}\right)+\left(\begin{array}{ccc}{\tilde{\Sigma }}_{xx}& {\tilde{\Sigma }}_{xy}& {\tilde{\Sigma }}_{xz}\\ {\tilde{\Sigma }}_{xy}& {\tilde{\Sigma }}_{yy}& {\tilde{\Sigma }}_{yz}\\ {\tilde{\Sigma }}_{xz}& {\tilde{\Sigma }}_{yz}& {\tilde{\Sigma }}_{zz}\end{array}\right),$$

$$0$$

Full text

I Introduction

Measurements of a fluid’s shear rheology are commonly performed in a torsional flow device, using either a cone-and-plate or cylindrical Couette flow cell. The former comprises a cone rotating relative to a stationary plate (or vice versa). The latter comprises two coaxial cylinders, with the inner cylinder rotating relative to the outer (or vice versa). See the schematic Fig. 1. In each case, the speed of rotation normalised by the gap size (or angle) sets the imposed shear rate, $\dot{\gamma}$. The resulting torque gives the shear stress response, $\sigma$. In a state of steady flow, the shear stress $\sigma$ as a function of shear rate gives the flow curve, $\sigma(\dot{\gamma})$, which is key to characterising a material’s rheological response. Other common tests of a fluid’s shear rheology include oscillatory shear, shear startup, step stress, and stress relaxation tests.

Rheological measurements performed in the regime of linear viscoelastic response are typically well controlled and highly reproducible. In contrast, the measurement of stronger, nonlinear flows are much more challenging. This is particularly true in highly viscoelastic materials such as entangled polymer melts and concentrated polymer solutions, as well as in concentrated suspensions. Attempts to measure the steady state flow curve or transient flow behaviour at high strain rate (or stress) are often beset by flow instabilities that can lead to unreliable data. Many such instabilities depend not only on the bulk rheology of the material in question, but also on the geometry of the flow device used, including the boundary conditions where the fluid sample meets the hard walls of the flow cell, and/or the outside air.

In both of the flow devices sketched in Fig. 1, the fluid (shown in blue) has an interface with the outside air (shown in white). In entangled polymers or concentrated suspensions, this free surface is highly susceptible to destabilising when the material is strongly sheared, particularly in a cone-and-plate or plate-plate device. Above a critical imposed shear rate, the sample edge will deform into a more complicated edge profile, forming an indentation that invades the fluid bulk. In violent cases, some portion of the sample can even be ejected from the measurement region.

This phenomenon, which is known as ‘edge-fracture’, renders accurate measurements of strong flows extremely difficult. Indeed, to quote Snijkers Snijkers and Vlassopoulos (2011), “the effects of edge fracture on rheological measurements performed with standard rotational rheometers using cone-and-plate or plate-plate geometries are catastrophic and render the torque (in a strain-controlled rheometer) or rotation speed (in a stress-controlled rheometer) as measures of bulk rheological properties of the fluid, virtually useless”. Jensen described the phenomenon as “the limiting factor in rotational rheometry” Jensen et al. (2008). Detailed experimental studies of edge fracture can be found in Refs. Lee et al. (1992); Inn et al. (2005); Sui and McKenna (2007); Schweizer and Stöckli (2008); Mattes et al. (2008); Jensen et al. (2008); Dai et al. (2013). Anecdotal reports further pervade the literature. In recent decades, strategies to mitigate edge fracture have been developed, involving specialised guard-ring Mall-Gleissle et al. (2002) or cone-partitioned-plate devices Schweizer (2003); Schweizer et al. (2004); Snijkers and Vlassopoulos (2011); Meissner et al. (1989); Schweizer and Schmidheiny (2013); Costanzo et al. (2018).

Despite the crucial importance of edge fracture to experimental shear rheology, the phenomenon has remained poorly understood theoretically, until recently. In insightful early works, Tanner and coworkers Tanner and Keentok (1983); Keentok and Xue (1999) predicted that edge fracture should arise for a critical magnitude of the second normal stress difference, $|N_{2}(\dot{\gamma})|>\Gamma/R$, where $\Gamma$ is the surface tension of the fluid-air interface and $R$ is some pre-assumed geometrical lengthscale. Taken as a scaling argument, this prediction shows remarkable insight. Indeed, a careful experimental study later confirmed the key role of $N_{2}$ (rather than the first normal stress difference, $N_{1}$) in driving edge fracture Lee et al. (1992); Keentok and Xue (1999). However, Tanner’s prediction will prove only partly consistent with our findings below. In particular, it fails to incorporate the role of the shear stress, alongside the second normal stress, in driving edge fracture.

In a recent Letter Hemingway et al. (2017), we performed a detailed theoretical study of edge fracture, combining linear stability analysis with full nonlinear simulations. Outcomes of this work that we hope will be useful to the experimental rheology community include: (i) a new criterion for the onset of edge fracture that now recognises the importance, alongside $N_{2}$, of the slope of the flow curve of shear stress as a function of shear rate, and the jump in shear stress between the fluid sample and outside air; (ii) a full mechanistic understanding of the edge fracture instability; and (iii) a new suggestion for how edge fracture might potentially be mitigated in experimental practice.

The present paper aims to provide a fuller discussion of the results originally announced in Ref. Hemingway et al. (2017). In particular, we set out here for the first time the details of the analytic linear stability calculation by which we arrived at the new onset criterion. We also provide new simulation results, beyond those in Ref. Hemingway et al. (2017), delineating more fully the dependence of the onset threshold on the nonlinear parameters of the two constitutive models considered (which set the scaling of the shear and normal stresses with strain rate), and demonstrating the role of the wetting angle in the mode of edge fracture. For definiteness we cast the discussion mostly in the language of entangled polymers. However, we note that the constitutive models used here are not restricted to such materials alone, and our results may indeed apply more generally to other classes of fluids. It would be particularly interesting in future work to consider the edge fracture instability in non-Brownian suspensions and in viscoelastic yield stress materials.

Besides edge fracture, another important instability that routinely confounds attempts to measure strong flows of complex fluids is that of wall slip Hatzikiriakos (2015, 2012): the layer of fluid immediately adjacent to the hard wall of the flow cell shows an apparent slip relative to the wall itself. We shall ignore wall slip in what follows, suppressing the phenomenon upfront by assuming conditions of no-slip. This is potentially a major shortcoming of our approach, which should be addressed in future studies.

The paper is structured as follows. In Sec. II we introduce the flow geometry to be studied. Sec. III sets out the models and methods that we shall use. The results of our nonlinear simulations are discussed in Sec. IV, followed in Sec. V by a linear stability analysis to understand the phase boundary for the onset of edge fracture, as obtained in our nonlinear simulations. This linear analysis also permits a detailed understanding of the mechanism of instability. In Sec. VI, we suggest a way in which our findings indicate a possible practical route to mitigating edge fracture experimentally. Finally, Sec. VII provides conclusions and perspectives for future work.

II Flow geometry

As noted above, measurements of a fluid’s shear rheology are often performed in a rotational flow device, commonly using either a cone-and-plate or cylindrical Couette flow cell, as shown schematically in Fig. 1. For the former device, we assume a small cone angle and large cell radius. In this limit, the curvature of the streamlines becomes negligible and the flow is well approximated by a planar slab of fluid sheared between flat plates, as sketched in Fig. 1 (top right). (In experimental practice, however, we note that very small cone angles are difficult to obtain because of the limited stiffness of rheometer frames; while large samples are forbidden by the limited normal force capacity of the transducers and the volume of sample available.) For the latter device, we consider the limit $\delta R/R\ll 1$, in which the gap between the cylinders, $\delta R$, is small compared with the inner cylinder radius, $R$. In this limit, the flow is again well approximated by a planar slab.

Throughout what follows, therefore, we shall consider a planar slab of fluid in contact with the outside air, as sketched again in Fig. 2. The sample is sheared at some rate $\dot{\gamma}$ by moving the top boundary at some speed $\dot{\gamma}L_{y}$ in the flow direction $\hat{\mathbf{x}}$, into the page. We consider always positive shear rates, $\dot{\gamma}>0$. (For most fluids, the second normal stress $N_{2}(\dot{\gamma})\sim-\dot{\gamma}^{2}$ at low shear rates. For negative flow rates we would need to replace $d|N_{2}|/d\dot{\gamma}$, which appears several times below, by $d|N_{2}|/d|\dot{\gamma}|$.)

We denote the flow-gradient direction by $\hat{\mathbf{y}}$, shown vertically in Fig. 2. The edges of the sample in the vorticity direction $\hat{\mathbf{z}}$ (horizontal) are in contact with the air. Only the left half of the box is shown in Fig. 2: an equivalent interface exists in the right half. We assume that the flow remains translationally invariant in the flow direction, $\hat{\mathbf{x}}$, setting $\partial_{x}$ of all quantities equal to zero, and performing two-dimensional (2D) calculations in the gradient-vorticity plane $y-z$. (All simulation snapshots below will therefore show only this $y-z$ plane.) The velocity vector and stress tensor are themselves however 3D objects. This assumption of only 2D variations should be checked in future, fully 3D studies.

The sample length in the vorticity direction $z$ at the cell mid-height $y=0$ (before the fluid is sheared) will be denoted $\Lambda$. The simulation box has length $L_{z}$ with periodic boundary conditions in the vorticity direction. In the flow-gradient direction, $y$, we consider two different kinds of boundary condition. The first models hard walls at $y=\pm L_{y}/2$, with boundary conditions of no slip or permeation at these walls. The second considers the simplified theoretical geometry of Lees-Edwards biperiodic shear. In that case, all quantities repeat periodically across shear-mapped points on the boundaries of box copies stacked in the $y$ direction, but with adjacent copies moving relative to each other at velocity $\dot{\gamma}L_{y}\hat{\mathbf{x}}$. The threshold for the onset of edge fracture obtained in our numerical simulations will prove to be in excellent agreement between these two cases. Our reason for invoking the simplified biperiodic geometry is that it will allow progress in analytical calculations, which would otherwise be prohibitively complicated.

III Models

III.1 Force balance and incompressibility

We assume that inertia is negligible, and work in the creeping flow limit of zero Reynolds number. In this limit, the condition of force balance requires the total stress in any element of fluid (or air), $\mathbf{{T}}(\mathbf{{r}},t)$, to obey:

[TABLE]

As usual, $\mathbf{{r}}$ denotes position and $t$ time.

Inside the fluid, we assume $\mathbf{{T}}$ to comprise an isotropic contribution described by a pressure $p(\mathbf{{r}},t)$, a Newtonian contribution of viscosity $\eta_{\rm s}$, and a viscoelastic contribution $\mathbf{{\Sigma}}(\mathbf{{r}},t)$ stemming from the complex fluid microstructure (polymer chains, etc.). The condition of force balance inside the fluid then reads:

[TABLE]

in which $\mathbf{{v}}(\mathbf{{r}},t)$ is the fluid velocity.

The air outside the fluid lacks any viscoelastic component, $\mathbf{{\Sigma}}=0$, and has a lower viscosity, $\eta_{\rm a}\ll\eta_{\rm s}$, giving the force balance condition:

[TABLE]

The pressure field $p(\mathbf{{r}},t)$ is determined by assuming the flow to be everywhere incompressible, with the flow velocity $\mathbf{{v}}(\mathbf{{r}},t)$ obeying

[TABLE]

Note that in computing everywhere in the fluid and air the flow rate, our simulation naturally captures the changes shear rate that will be present in the vicinity of any edge-fracturing disturbances in the fluid-air interface.

III.2 Constitutive models

The dynamics of the viscoelastic stress $\mathbf{{\Sigma}}$ in flow is determined by a viscoelastic constitutive equation. In what follows, we shall study two constitutive models that are widely used across the rheological literature. Our aim in studying two different models is to establish the degree to which any predictions concerning the onset of edge fracture are generic across constitutive models, or whether they instead depend on model details. Importantly, indeed, we shall show that the predictions of these two models for edge fracture depend on their respective parameters $a$ or $\alpha$ (defined below) only via the appearance of those quantities in the shear stress $T_{xy}$ and second normal stress difference $N_{2}=T_{yy}-T_{zz}$. In this way, the key physics will prove robust to choice of model.

We note, however, that both models studied here have a single viscoelastic relaxation time. This is an approximation that should be relaxed in future studies, by adopting multi-mode models.

III.2.1 Johnson-Segalman model

The first model that we consider is the Johnson-Segalman model Johnson Jr and Segalman (1977), in which the viscoelastic stress evolves according to

[TABLE]

Here $\mathbf{{D}}=\tfrac{1}{2}(\nabla\mathbf{{v}}+\nabla\mathbf{{v}}^{T})$ is the symmetric part of the rate of strain tensor $\nabla\mathbf{{v}}_{\alpha\beta}=\partial_{\alpha}v_{\beta}$. Its antisymmetric counterpart $\mathbf{{\Omega}}=\tfrac{1}{2}(\nabla\mathbf{{v}}-\nabla\mathbf{{v}}^{T})$ is the vorticity tensor. The parameter $a$ describes a relative slip of the viscoelastic component compared with the deformation of the background solvent. It must lie in the range $-1\leq a\leq 1$. The simpler Oldroyd B model is recovered at $a=1$.

In a stationary homogeneous simple shear flow, the viscoelastic shear stress and second normal stress difference obey the following functions of the imposed shear rate $\dot{\gamma}$:

[TABLE]

For values of the parameters $|a|<1$ and $\eta_{\rm s}<1/8$, the total shear stress $T_{xy}=\Sigma_{xy}+\eta_{\rm s}\dot{\gamma}$ is then a non-monotonic function of the imposed shear rate. In this regime, coexisting bands of differing shear rates can form at a common value of the total shear stress: a phenomenon known as shear banding. In this work, we consider only flows that are not shear banded, and so confine ourselves to values of the solvent viscosity $\eta_{\rm s}>1/8$. The second normal stress is negative, scaling as $-\dot{\gamma}^{2}$ at low shear rates, before saturating to a negative constant at high shear rates.

The spatial gradient terms prefactored by $\ell$ in Eqn. 5 are included in our simulations but do not affect the predictions for edge fracture. (They would be important if the flow were shear banded Lu et al. (2000), as considered in Refs. Skorski and Olmsted (2011); Hemingway and Fielding (2018).) The microscopic lengthscale $\ell$ is small compared with any bulk lengthscales.

III.2.2 Giesekus model

The second constitutive model that we consider is the Giesekus model Giesekus (1982), in which the viscoelastic stress evolves according to:

[TABLE]

Here $\alpha$ is an anisotropy parameter, which captures an increased rate of stress relaxation in any regime where the polymer chains are more strongly aligned. The Oldroyd B model is recovered at $\alpha=0$. In a state of stationary homogeneous simple shear flow, the viscoelastic shear and second normal stresses obey Giesekus (1982)

[TABLE]

in which

[TABLE]

The shear stress $\Sigma_{xy}$ is then a non-monotonic function of imposed shear rate $\dot{\gamma}$ for values of the anisotropy parameter $\alpha>1/2$ Giesekus (1982). We shall therefore mostly restrict ourselves to values of $\alpha<1/2$, again in order to avoid shear-banding. Where we do consider a value of $\alpha>1/2$, we take a large enough solvent viscosity such that the total shear stress $T_{xy}=\Sigma_{xy}+\eta\dot{\gamma}$ is monotonic, again avoiding banding. The second normal stress is negative, as in the Johnson-Segalman model: scaling as $-\dot{\gamma}^{2}$ at low shear rates, before saturating to a negative constant at high shear rates.

It is worth remarking that both the Johnson-Segalman and Giesekus models are highly phenomenological. Indeed, although both are widely used in the rheological literature, each has notable pathologies. For example, the Johnson-Segalman model predicts unphysically large and sustained oscillations in the stress in shear startup at high shear rates. However, the spirit of this paper is to derive an instability criterion for edge fracture that is independent of the particular choice of constitutive model, at least in the limit of low strain rates, where both models reduce to a second order fluid as just described. We defer to future work a simulation of more microscopically faithful constitutive models in this geometry: noting, however, that that one such candidate - the Rolie-poly model Likhtman and Graham (2003) - is excluded by in fact predicting zero second normal stress difference. It would be interesting in future work also to consider the recent tube-based modelling approach of Ref. Costanzo et al. (2018).

III.3 Fluid-air coexistence

In our numerical simulations, we model the coexistence of fluid and air using a phase field approach with an order parameter $\phi(\mathbf{{r}},t)$, which obeys Cahn-Hilliard dynamics Bray (2002)

[TABLE]

Here $M$ is the molecular mobility, which we assume to be constant. The chemical potential

[TABLE]

in which $G_{\mu}$ sets the scale of the free energy of demixing per unit volume. The functional dependence of this free energy on $\phi$ captures the coexistence of a fluid phase, in which the order parameter $\phi=1$, with an air phase, in which $\phi=-1$. The elastic modulus $G$ and relaxation time $\tau$ that appear in the constitutive equations 5 and 7 are then made functions of $\phi$, with $G(\phi=1)=1$, $\tau(\phi=1)=1$ and $G(\phi=-1)=0$, with $\tau(\phi=-1)=0.002$. In this way, viscoelastic stresses arise only in the fluid phase.

The fluid and air bulk phases are separated by a slightly diffuse interface of thickness $\ell_{\mu}$, set by the prefactor to the gradient term in the free energy of Eqn. 11. This interfacial thickness is small compared with any bulk lengthscales. Gradients in $\mu$ contribute an additional source term of the form $-\phi\nabla\mu$ to the force balance condition. This is important in the vicinity of the fluid-air interface, where it confers forces arising from the surface tension of the interface. The surface tension that emerges out of these Cahn-Hilliard dynamics is given by:

[TABLE]

In having a slightly diffuse fluid-air interface, our simulations properly capture any motion of the contact line that arises where the fluid-air interface meets the hard walls of the flow cell Kusumaatmaja et al. (2016). We thereby avoid the contact line singularity that would arise if the interface were perfectly sharp.

In our analytical calculations, which are performed in the biperiodic geometry without walls, we shall instead assume the interface between the fluid and air to be infinitely sharp, although still with an equivalent surface tension $\Gamma$. (This sharpness is unimportant in the absence of hard walls, and therefore of any contact line, at least in the early stages of any instability, before bubbles form.) We shall further assume in our analytical calculations that the viscosity of the outside air is negligible compared to that of the fluid. Below we shall demonstrate full agreement between our numerical simulations and analytical calculations, in the biperiodic flow cell, in the physically relevant limit in which the interface is indeed thin compared with any bulk scales, and in which the air viscosity is small compared with that of the fluid.

III.4 Boundary conditions

We adopt periodic boundary conditions in the vorticity direction $z$. Where the fluid meets the hard walls of the flow cell in the $y$ direction, we assume boundary conditions of no slip and no permeation for the fluid velocity $\mathbf{{v}}$, and zero-gradient for the viscoelastic stress $\mathbf{{\Sigma}}(\mathbf{{r}},t)$:

[TABLE]

Here $\mathbf{n}$ is the outward unit vector normal to the wall.

For the phase field that captures the fluid-air coexistence, as just described, the boundary conditions are Yue et al. (2010); Dong (2012)

[TABLE]

In unsheared equilibrium, the contact angle at which the fluid-air interface meets the hard walls of the flow cell is given by $\theta$. A value $\theta=90^{\circ}$ gives a vertical equilibrium interface. A value $\theta>90^{\circ}$ gives an interface convex into the air. A value $\theta<90^{\circ}$ gives a concave interface.

The simplified biperiodic Lees-Edwards geometry has no walls and its equilibrium interface is always vertical, mimicking the case $\theta=90^{\circ}$ with walls. The boundary conditions for this case were discussed in Sec. II above.

III.5 Initial conditions

As the initial condition for our simulations/calculations in shear, we take a state of air-fluid coexistence that has first been equilibrated in the absence of shear. At the end of the equilibration dynamics, just before shear is applied, a small perturbation is added to the interface’s position $h(y)$ along the $z$ direction, $h(y)\to h(y)+10^{-8}\cos(n\pi y/L_{y})$, to trigger edge fracture. We take $n=1$ with walls and $n=2$ in the biperiodic geometry. The viscoelastic stress is initialised so as to lie on the stationary homogeneous constitutive curve, for the shear rate in question.

III.6 Nonlinear simulation method

Each numerical timestep comprises two separate substeps. In the first substep, we enforce the force balance condition in order to update the fluid velocity field $\mathbf{v}$ at fixed phase field $\phi$ and polymer stress $\mathbf{\Sigma}$. We do so using a streamfunction formulation to ensure incompressible flow. In the second substep, we update the phase field and viscoelastic stress, with the velocity field fixed.

The advective terms are implemented using a third order upwinding scheme Pozrikidis (2011). The spatially local terms in the viscoelastic constitutive equation are updated using an explicit Euler scheme Press et al. (1992). To implement the spatially diffusive terms, we use a 2D Fourier spectral method in the Lees-Edward biperiodic geometry. With walls present, we use Fourier modes in the periodic direction $z$, and finite differencing Press et al. (1992) in the flow gradient direction $y$, with wetting conditions implemented using the method described in Ref. Dong (2012).

As noted above, our simulations explicitly model the coexistence of air and polymeric fluid, using a phase field approach: the phase field $\phi=1$ inside the polymeric fluid and $\phi=-1$ in the outside air. In addition to the viscoelastic component in the polymeric phase, we also have a background Newtonian viscosity equal to that of the solvent $\eta_{s}$ inside the polymeric fluid, and a lower value $\eta_{a}$ in the air.

To achieve this spatially dependent Newtonian viscosity, we simulate everywhere a Newtonian fluid with the air viscosity $\eta_{a}$, and further introduce an additional (and nominally viscoelastic) stress tensor ${\boldsymbol{\sigma}}_{s}$ which obeys Oldroyd-B dynamics with a small enough relaxation time $\tau_{s}\ll 1$ so as to be essentially Newtonian, and with a (spatially dependent) elastic modulus $G_{s}(\phi=-1)=0,G_{s}(\phi=1)=(\eta_{s}-\eta_{a})/\tau_{s}$. This approximates a Newtonian solvent with viscosity $\eta_{s}=\eta_{a}+G_{s}(\phi=1)\tau_{s}$ in the polymer phase only and a fluid with viscosity $\eta_{a}$ in the air phase.

III.7 Units and parameter values

The parameters contained in the model equations, geometry and boundary conditions just described are summarised in Table. 1, along with values to be used in our numerics. As can be seen, the four important quantities to be varied in our study are the dimensionless shear rate, $\dot{\gamma}\tau$, the dimensionless surface tension of the air-fluid interface, $\Gamma/GL_{y}$, the (obviously dimensionless) equilibrium contact angle $\theta$, and the dimensionless nonlinear parameter, $a$ or $\alpha$, that set the way in which the shear and normal stresses depend on shear rate in either constitutive model. The other parameters do not affect the key physics, provided each takes a large or small value appropriate to the physical limit of interest. For example, we set microscopic lengths small compared to macroscopic lengths. We further set the viscosities of the solvent and air small compared to that of the viscoelastic component. Even though the separation of lengthscales (or viscosities) that can feasibly achieved in our simulations is less than that in physical reality, the results are unaffected by this limitation, to within negligible corrections. Again for numerical convenience, the value of air viscosity that we use is much larger than that of actual air, but we have checked that the results obtained do not change upon halving the value used.

IV Nonlinear simulation results

We now present the results of our nonlinear simulations. The basic phenomenon of edge fracture is exemplified by the five snapshots shown in Fig. 3 at some long time $t$ after the switch-on of shear at time $t=0$, for the Giesekus model simulated between hard walls. At any given imposed strain rate, an air-fluid interface with a high surface tension is undisturbed by the flow and retains its equilibrium shape, as in snapshot (a). (We shall denote such undisturbed states by black crosses in the phase diagram Fig. 4, which we discuss below.) For intermediate values of the interfacial surface tension, the interface partially fractures, displacing in the $z$ direction by a distance $O(L_{y})$ that is set by the gap between the rheometer plates in the $y$ direction, before finally settling to a new steady state shape, different from its unsheared equilibrium one. See snapshot (b) in Fig. 3. (We denote such partially fractured states by circles in Fig. 4.) For low values of the interfacial surface tension, the interface fully fractures, displacing in the $z$ direction a distance $O(\Lambda)$ set by the sample width in that direction. See snapshots (c)-(e) in Fig. 3. (We denote such fully fractured states by yellow closed triangles in Fig. 4.) In such cases, the system never attains a new steady state. Instead, the sample may (for example) completely de-wet either wall, and/or air bubbles may invade the fluid. Which of these happens in practice depends on the wetting angle in a way that we shall investigate further below.

In Fig. 4, we collect into a phase diagram the results of many such simulations across a full range of values of the surface tension and imposed shear rate, for both the Johnson-Segalman model (left panel) and the Giesekus model (right panel). The symbol shapes (crosses, circles, triangles) are as described in the previous paragraph. The state-points corresponding to the snapshots of Fig. 3 for the Giesekus model are indicated by letters (a)-(e) in the right panel. As can be seen, broadly the same phase behaviour arises in both models, with stability at high values of interfacial surface tension, and instability for low surface tension. At any fixed (low) value of surface tension, we find stability for weak flows $\dot{\gamma}\to 0$, as expected, followed by a window of instability for intermediate values of the shear rate, and finally a regime of re-entrant stability at high strain rates. We shall discuss this re-entrance in more detail below. It arises from a saturation (in both constitutive models) of the growth of $|N_{2}|$ with shear rate $\dot{\gamma}$, which may not be physically realistic. Even if the re-entrant regime does in principle exist in real fluids, it may be hard to access in experimental practice: one would either need to sweep the shear rate upward through the regime of shear rates where edge fracture does arise, or alternatively perform a shear startup at a high shear rate, which is in itself likely to be unstable. From an experimental viewpoint, it is likely that polymers melts mostly lie in the unstable regime of Fig. 4, at all but the lowest strain rates; whereas some polymer solutions, with their lower modulus $G$, may be in the stable regime.

The colorscale of the symbol fill in Fig. 4, in parameter regimes where (either partial or full) edge fracture arises, denotes the degree to which the interface has deformed at long times after the switch-on of flow at time $t=0$. This is defined as follows. We denote by $h(y,t)$ the location along $z$ of the interface at height $y$ across the gap at time $t$. In the sheared state at any time $t$, ignoring any bubbles that form, we subtract the leftmost position of the interface from the rightmost position to give $D(t)={\rm max}(h(y,t))-{\rm min}(h(y,t))$. From this we subtract the value of $D(t=0)$ in the initial unsheared state. (For an equilibrium contact angle $\theta=90^{\circ}$, $D(0)=0$.) The difference between these, normalised by the gap width $L_{y}$, gives the degree of interfacial deformation in shear, $d(t)=[D(t)-D(0)]/L_{y}$. The values of this quantity shown in yellow in the triangles of Fig. 4 should be interpreted as lower bounds, because no steady state is reached as $t\to\infty$ in such cases.

The black solid line in Fig. 4 shows the threshold between a fully stable interface (states shown by crosses) and an interface that partially fractures (states shown by circles). This is calculated by measuring, in the simulation performed at each value of $\Gamma,\dot{\gamma}$, the weight in the Fourier modes $\exp(iqy)$ of the interfacial displacement $h(y,t)-h(0,t)$ as a function of time $t$. In the unstable regime (below the black line), this (initially) grows exponentially for some modes; in the stable regime (above the black line), it decays exponentially for all modes, from our initially small perturbation $O(10^{-8})$. Plotting the growth (or decay) constant (ie, the eigenvalue) of the most dangerous mode as a function of surface tension $\Gamma$ at any fixed $\dot{\gamma}$ enables us to find the value of $\Gamma(\dot{\gamma})$ at which the interface is neutrally stability, with an eigenvalue of zero for the most dangerous mode. This gives the threshold shown by the black solid line. The red dashed line is the analytical prediction given by the linear stability calculation of Sec. V below, which is performed in the limit of low strain rate.

We emphasise that Fig. 4 makes quantitative predictions for the onset of edge fracture as a function of shear rate and surface tension that are testable experimentally, and we hope that this will motivate future experimental work along these lines.

So far, we have explored the phase behaviour across a wide range of values of surface tension and imposed shear rate, within the Johnson-Segalman model for one fixed value of the slip parameter $a$, and within the Giesekus model for one fixed value of the anisotropy parameter $\alpha$. We now explore the dependence of this phase behaviour on the value of $a$ (in Johnson-Segalman) and $\alpha$ (in Giesekus). This is shown in Fig. 5. In each case, the neutral stability curve extracted from the nonlinear simulations in the way described in the previous paragraph is shown as a solid curve, and the prediction of the linear stability analysis of Sec. V is shown by the dashed curve. The results for the Johnson-Segalman model are given in the left panel. As can be seen, the unstable regime becomes more limited for large values of $a$. This is consistent with a scenario in which the edge fracture instability is driven by the second normal stress $N_{2}$: the Johnson-Segalman model recovers the Oldroyd B model in the limit $a\to 1$, with zero second normal stress $N_{2}(\dot{\gamma},a=1)=0$ (recall Eqn. 6), but non-zero first normal stress. The results for the Giesekus model are shown in the right panel. The unstable regime again becomes more limited for small values of $\alpha$, consistent with the Oldroyd B model being recovered for $\alpha=0$.

The simulation results shown so far have all been for a value of the equilibrium contact angle $\theta=90^{\circ}$, for which the interface is initially flat. To check for robustness with respect to this choice of boundary condition, we now compare the threshold for the onset of edge fracture for different values of the equilibrium contact angle, for a channel with hard walls. Recall that for the particular case of $\theta=90^{\circ}$, we calculated the threshold by identifying the value of surface tension (at any fixed imposed strain rate) at which the eigenvalue of the most dangerous mode crosses zero. A method of determining the threshold that applies more easily across values of $\theta\neq 90^{\circ}$ is instead to identify the value of surface tension (at any fixed imposed shear rate) below which the degree of interfacial deformation, $d$, first exceeds $0.1$. (Good agreement is obtained between these two methods at $\theta=90^{\circ}$; the eigenvalue method is only used at $\theta=90^{\circ}$.) As seen in Fig. 6, the threshold for the onset of edge fracture is indeed robust (to within small corrections) against variations of the equilibrium contact angle.

Despite this robustness in the onset threshold, the spatial mode of instability does vary with the equilibrium contact angle, particularly for values of the surface tension deep inside the unstable regime. This can be seen in the snapshots of Fig. 7. For values of the contact angle $\theta<90^{\circ}$, for which the fluid tends to wet the walls, the air first invades the fluid at a central location across the channel width $y$, leading to the formation of bubbles within the sample’s bulk. For values $\theta>90^{\circ}$, the air instead invades along the channel walls, often leading the fluid eventually to de-wet the walls entirely. (The air from the left side of the cell shown in Fig. 7 finally joins up with that from the right side, not shown.)

We note that larger angles, such as $\theta=120^{\circ}$, correspond to an experimental case where a polymer sample is squeezed between a cone and plate without trimming. In contrast, $\theta=90^{\circ}$ corresponds to the ideal case where the surface of the sample is part of a sphere. Therefore, results for experimental papers on melts are likely to fall in between these two cases. Indeed, the results in Fig. 7 corroborate the experimental observations of Fig. 3 of Schweizer and Stöckli (2008), in which a higher shear layer was seen near the walls for an initially convex sample profile, whereas a higher shear layer was seen mid-gap for an initially concave sample profile.

All the simulation results presented so far have been for a fluid sheared between hard walls. We finally compare those results to the corresponding results obtained in a biperiodic flow cell. As seen in Fig. 6, good agreement is obtained between these two cases, particularly at low strain rates $\dot{\gamma}\to 0$. Note that this agreement is obtained after rescaling the box size $L_{y}$ by a factor $2$ (as discussed further in the caption of Fig. 6) to allow for the (relatively trivial) fact that the minimum possible wavevector $q_{\rm min}=\pi/L_{y}$ with hard walls, whereas $q_{\rm min}=2\pi/L_{y}$ in the biperiodic geometry. This will enable us to perform a linear stability analysis for the onset of edge fracture in the simpler geometry of the biperiodic flow cell, for which the eigenfunctions are harmonic functions of $y$. It is to this linear analysis that we now turn.

V Linear stability analysis

We now perform a linear stability analysis to derive a criterion for the onset of edge fracture. This calculation is performed in the limit of low strain rates within the biperiodic flow geometry. We now take the viscosity of the air outside the fluid to be zero (it was small but non-zero in our nonlinear simulations), and assume the fluid-air interface to be infinitely thin (it had a slightly diffuse width in our nonlinear simulations), although still with a finite surface tension $\Gamma$. Readers who do not wish to follow the details can skip directly to the result in Eqn. 37.

V.1 Governing equations

We first recall for convenience the governing equations.

V.1.1 Force balance in the fluid bulk

Inside the fluid bulk, the flow obeys the incompressibility condition

[TABLE]

together with the force balance condition

[TABLE]

in which the total stress

[TABLE]

V.1.2 Bulk viscoelastic constitutive equation

The viscoelastic stress obeys a constitutive equation

[TABLE]

with $-1<a<1$ and $\alpha=0$ in the Johnson-Segalman model, and $a=1$ and $0<\alpha\leq 1$ in the Giesekus model. Making the twin assumptions that both the time derivative and the advective derivative on the left hand side of Eqn. 19 can be neglected (we shall return below to justify these assumptions), we find that for slow flows, $\dot{\gamma}\tau<1$, the solution to Eqn. 19 can be written as an expansion in increasing powers of the flow rate tensors $\mathbf{{D}},\mathbf{{\Omega}}$:

[TABLE]

We ignore the higher order terms (h.o.t.) in what follows. In this case, the Johnson-Segalman and Giesekus models both reduce to this so-called ‘second-order fluid’. Note that it was identified in Ref. Porteous and Denn (1972) that the second order fluid contains an unphysical mode of instability because of the inherent instability of the rest state and incompatibility with the second law of thermodynamics if the disturbance timescale is not large compared to the characteristic relaxation timescale of the fluid. Because we are using the second order fluid in our analytical work to calculate the boundary of instability onset, where the disturbance timescale is by definition infinity, our use of the second order fluid is deemed safe in this regard. We further note that the analytical results obtained within this second order fluid calculation agree quantitatively with our numerical simulations of the Johnson-Segalman and Giesekus models at low strain rates (and qualitatively at high strain rates).

V.1.3 Force balance across fluid-air interface

The condition of force balance across the interface between the fluid and outside air is given by

[TABLE]

where $\mathbf{{n}}$ is the interface normal into the fluid and $\nabla_{\rm int}$ is the interfacial gradient operator. (In the diffuse interface case, this condition of interfacial force balance emerges naturally from having a term of the form $\phi\nabla\mu$ in the bulk force balance equation. In the sharp interface limit that we consider in this analytical calculation, the interfacial condition provided by Eqn. 21 must be imposed separately.)

V.1.4 Motion of fluid interface

Finally, the position $h(y)$ along the $z$ axis of the interface at any height $y$ across the flow cell moves with the $z-$component of fluid velocity as

[TABLE]

V.2 Base state and small perturbation

To perform a linear stability analysis, we represent the state of the system by an underlying homogeneous base state, denoted by subscript [math]. This base state is stationary, in the sense of being a time-independent solution of the model equations, but not (necessarily) stable against the onset of edge fracture. It corresponds to a state of stationary shear in an initially unfractured sample, in which the interface is flat and the flow uniform. To this, we add a small perturbation, denoted by over-tildes, which represents the precursor of edge fracture.

Accordingly, in the fluid bulk we write the velocity field

[TABLE]

in which $\tilde{v}_{x}=\tilde{v}_{x}(y,z,t)$ and $\tilde{\psi}=\tilde{\psi}(y,z,t)$. In restricting the velocity components to depend only on $y$ and $z$, and in writing the $y$ and $z$ velocity components as the $y-z$ curl of a streamfunction $\psi$, we automatically ensure that the incompressibility condition of Eqn. 16 is obeyed.

The strain rate tensor follows from Eqn. 23 as:

[TABLE]

We similarly write the viscoelastic stress tensor as the sum of a homogeneous base state plus small perturbations:

[TABLE]

in which $\tilde{\Sigma}_{ij}=\tilde{\Sigma}_{ij}(y,z,t)$. The total stress tensor $\mathbf{{T}}=\mathbf{{T}}_{0}+\tilde{\mathbf{{T}}}$ follows an analogous componentwise format, which we do not write out.

Our strategy now is to substitute Eqns. 23 to 25 into the governing Eqns. 16 to 21, and expand to first order in the amplitude of the perturbations. This will result in a set of linearised equations that govern the dynamics of the perturbations, the solution of which will allow us to delineate the regime in which the perturbations grow to give edge fracture.

V.3 Linearised equations of motion

V.3.1 Linearised force balance condition in the fluid bulk

The conditions of force balance and incompressibility, Eqns. 16 to 18, are linear in the quantities they contain. Therefore, the linearised force balance condition is simply $\nabla.\tilde{\mathbf{{T}}}=0$. This can be written componentwise as:

[TABLE]

The first of these is the $x$ component of force balance. The second is the curl in the $y-z$ plane of the $y$ and $z$ components of force balance. Writing the perturbed total stress as the sum of the perturbed Newtonian stress, the perturbed viscoelastic stress, and the perturbed pressure, $\tilde{\mathbf{{T}}}=2\eta\tilde{\mathbf{{D}}}+\tilde{\mathbf{{\Sigma}}}-\tilde{p}\,\mathbf{{I}}$, we can exactly rewrite Eqns. 26 as:

[TABLE]

in which the 2D Laplacian $\nabla\equiv\partial_{y}^{2}+\partial_{z}^{2}$. (The pressure does not appear in Eqn. 27 because $\partial_{x}\tilde{p}=0$ in the first equation and the $yz$ curl of $\nabla\tilde{p}$ is zero in the second.)

Substituting into Eqn. 27 the componentwise forms of $\tilde{\Sigma}_{ij}$ that will be set out in the next subsection, we find finally the linearised bulk equations:

[TABLE]

V.3.2 Linearised bulk viscoelastic constitutive equation

For parameter values in the vicinity of the threshold of onset of instability in the phase diagram (as in Fig. 4, for example), edge fracture develops only slowly. This allows us to neglect the term $\partial_{t}\mathbf{{\Sigma}}$ in the viscoelastic constitutive equation. At the level of terms that are linear in the perturbation, the advective term, $\mathbf{{v}}\cdot\nabla\mathbf{{\Sigma}}$, is also negligible, because the perturbed quantities vary only in the $y-z$ plane, which is orthogonal to the velocity of the base flow. For imposed shear rates $\dot{\gamma}_{0}\tau\ll 1$, this enables us to work with the form of the viscoelastic constitutive equation given by Eqn. 20. Expanding this to linear order in the amplitude of the perturbations, we get

[TABLE]

Here we use the notation that $\mathbf{{A}}^{\rm S}=\tfrac{1}{2}(\mathbf{{A}}+\mathbf{{A}}^{\rm T})$, for any tensor $\mathbf{{A}}$. Componentwise, Eqn. 29 can be written as:

[TABLE]

V.3.3 Linearised force balance across air-fluid interface

In similar spirit, we write the position along the $z$ axis of the interface at any location $y$ across the gap as $h_{0}+\tilde{h}(y,t)$. We further choose the origin of the $z$ axis to lie at the location of the (unperturbed) interface, so that $h_{0}=0$, with the positive $z$ direction defined so as to have fluid for $z>0$ and air for $z<0$. To first order, the interface normal $\hat{\mathbf{{n}}}$ in Eqn. 21 is then written $\hat{\mathbf{{n}}}=\hat{\mathbf{{z}}}-\partial_{y}\tilde{h}\hat{\mathbf{{y}}}$. In linearised form, the interfacial force balance condition, Eqn. 21, is then written componentwise as:

[TABLE]

Here we have denoted by $N_{2}$ the second normal stress difference in the base state, $\Sigma_{yy0}-\Sigma_{zz0}$, and by $\Delta\sigma$ the shear stress in the base state, $T_{xy0}=\Sigma_{xy0}+\eta_{\rm s}\dot{\gamma}_{0}$. We use the notation $\Delta\sigma$ rather than simply $\sigma$, because the relevant quantity is the jump in shear stress between the fluid and outside air. (In this linear calculation we are however assuming the stress in the outside air to be zero, so $\Delta\sigma=\sigma$. $N_{2}$ is always zero in the outside air, with no $\Delta$ notation needed in that case.) The notation $\tilde{T}_{ij}|_{z=0^{+}}$ denotes the $ij$-th component of the perturbation to the polymer stress immediately next to the interface, just inside the polymer phase.

While it may seem counter intuitive to discuss a jump in the shear stress across an interface, it is important to emphasise that by shear stress we mean $\Sigma_{xy}$, consistent with the main flow direction being $\hat{\mathbf{x}}$ and flow gradient direction being $\hat{\mathbf{y}}$. It is this quantity that jumps across the interface, which has its normal in the vorticity direction $\hat{\mathbf{z}}$. Because of the mutually different directions involved, a jump in $\Sigma_{xy}$ across a (perfectly flat) interface with normal along $\hat{\mathbf{z}}$ does not violate the condition of force balance across the interface.

Inserting into Eqn. 31 the componentwise form of the perturbation to the total stress, $\tilde{T}_{ij}=\tilde{\Sigma}_{ij}+2\eta_{\rm s}\tilde{D}_{ij}-\tilde{p}\delta_{ij}$, and eliminating the perturbation to the pressure $\tilde{p}$ via use of the linearised bulk force balance equation, we find the final linearised condition of force balance across the fluid-air interface:

[TABLE]

V.3.4 Linearised equation of interfacial motion

The linearised form of the equation of interface motion, Eqn. 21, is:

[TABLE]

V.4 Criterion for edge fracture

The linearised bulk equations, Eqns. 28, the linearised condition of interfacial force balance, Eqns. 32 and the linearised equation of interfacial motion, Eqn. 33, can now together be solved to determine whether, for any given interfacial tension $\Gamma$ and imposed flow rate $\dot{\gamma}$, the heterogeneous perturbations grow towards an edge fractured state, or decay to leave a flat interface.

To leading order in $\dot{\gamma}$, and at any wavevector $q$ in the $y$ direction, the solution to the linearised bulk force balance condition, Eqn. 28, has the following normal mode form:

[TABLE]

In these equations, $k=q/\sqrt{1+\beta}$, with $\beta=G^{2}\tau^{4}(1+\alpha-a)(a-\alpha)\dot{\gamma}^{2}/(G\tau+\eta_{\rm s})^{2}$. The constants $A,B$ and $C$ can be determined by imposing the three boundary conditions of componentwise force balance across the fluid-air interface, Eqns. 32, although we do not write out the expressions for $A,B,C$ here.

A colourmap of this analytical solution for the eigenfunction $\tilde{\psi}(y,z)$ and $\tilde{v}_{x}(y,z)$ at a fixed time $t$ is shown in the left two panels of Fig. 8. Excellent agreement is obtained with the form extracted from the early time regime of our full nonlinear simulation, shown in the right two panels. As can be seen, the perturbation due to any interfacial disturbance of wavelength $2\pi/q$ decays on the same lengthscale $O(q^{-1})$ into the bulk.

The corresponding eigenvalue $\omega$ follows by inserting the normal mode solution for the interfacial position:

[TABLE]

together with the solution of Eqn. LABEL:eqn:eigenfunction for the streamfunction into the linearised equation of interface motion, Eqn. 33. This gives

[TABLE]

to within a small correction $O(\eta_{\rm s}/G\tau)$.

In the left panel of Fig. 9, we plot as coloured solid lines the time-evolution of the logarithm of the amplitude of the five modes with the lowest values of $q$ consistent with the boundary conditions, as obtained in the early time regime of our full nonlinear simulations. Each should follow the linear form $\omega(q)t$. Accordingly, we perform a best fits (black dashed lines) to the slopes of these lines. These best-fit slopes are plotted as a function of wavevector in the right panel to construct dispersion relation $\omega(q)$ of growth rate as a function of wavevector. We have plotted this as solid lines in the right panel, from simulations with three different numerical grids and fluid-air interfacial thicknesses. The results of the analytical prediction of Eqn. 36 (which assumes zero interfacial thickness) are shown by the black dashed line for comparison, with good agreement at low wavevectors.

Perturbations grow at any wavevector $q$ if their eigenvalue $\omega(q)>0$. This condition is most readily satisfied for the mode with the lowest wavevector that is consistent with the biperiodic boundary conditions, $q=2\pi/L_{y}$, consistent with the dispersion relation having its maximum value at the lowest $q$ in Fig. 9 (right). Accordingly, our final criterion for an initially flat fluid-air interface to undergo edge fracture is given by

[TABLE]

This criterion (after rescaling by the factor 2 discussed at the end of Sec. IV, to account for the difference between the biperiodic and walled geometries) is marked by the dashed line in Fig. 4, and gives good agreement at low shear rates with the onset of fracture in our full nonlinear numerical simulations. An alternative, equivalent form of this criterion is to write, more simply:

[TABLE]

Appearing on the left side of this onset criterion, written in the form (37), are several rheological quantities pertaining to the base state, as follows. First is the shear stress $\sigma(\dot{\gamma})$, which is shown as a function of shear rate by the solid red line in Fig. 10 (top). (Recall that $\Delta\sigma=\sigma$, within our assumption that the outside air viscosity is zero.) Second, in the denominator, is the derivative of this quantity, $d\sigma/d\dot{\gamma}$, which is shown by the red dashed line in Fig. 10 (middle). This is the local slope of the flow curve, sometimes referred to as the tangent viscosity. Third is the derivative of the amplitude of the second normal stress with respect to strain rate, $\partial_{\dot{\gamma}}|N_{2}(\dot{\gamma})|$, which is shown by the blue dashed line in Fig. 10 (middle). These three quantities all combine to give the quantity plotted in green in Fig. 10 (bottom), which indeed follows the shape of the instability thresholds in Fig. 4. The region of re-entrant stability at high strain rates can therefore now be understood as arising from the decrease of $\partial_{\dot{\gamma}}|N_{2}(\dot{\gamma})|$ at large strain rates, consistent with $|N_{2}(\dot{\gamma})|$ saturating to a constant as $\dot{\gamma}\tau\gg 1$.

It is worth a reminder at this point that the criterion derived in this section applies to the biperiodic flow geometry. In order to apply to the experimentally realisable case of flow between hard walls, the box size must be rescaled by a factor $2$, as described at the end of Sec. IV and in the caption of Fig. 6.

V.5 Mechanism of edge fracture

The results of our analytical calculation allow us to understand the basic mechanism of the edge fracture instability, as follows. Were the interface between the fluid and air to remain perfectly flat, the jump $\Delta\sigma$ in shear stress across it would be consistent with force balance. Imagine now that a small interfacial tilt $\partial_{y}\tilde{h}$ arises, as shown in the second column of Fig. 11. This exposes the jump in shear stress across the interface, potentially disturbing the $x-$component of force balance across the interface, as expressed in the first of Eqns. 31. To recover this $x-$component of interfacial force balance, a counterbalancing perturbation $\tilde{T}_{xz}|_{z=0^{+}}=iqh\Delta\sigma$ is required, as shown in the third column of Fig. 11. To maintain the $x$-component of force balance in the fluid bulk (recall the first of Eqns. 26), a corresponding perturbation $\tilde{T}_{xy}$ is needed, as sketched in the fourth column of Fig. 11. This is achieved via a perturbation $\tilde{\dot{\gamma}}=\partial_{y}\tilde{v}_{x}=qh\Delta\sigma/\sigma^{\prime}(\dot{\gamma})$ in the shear rate (fifth column of Fig. 11). The second normal stress $N_{2}\approx-(1-a+\alpha)\dot{\gamma}^{2}$ in the fluid bulk then suffers a corresponding perturbation $\tilde{T}_{yy}-\tilde{T}_{zz}|_{\rm shear}=-qh\Delta\sigma\,|N_{2}|^{\prime}(\dot{\gamma})/\sigma^{\prime}(\dot{\gamma})$ (sixth column of Fig. 11). This must be counterbalanced by an equal and opposite extensional perturbation: $\tilde{T}_{yy}-\tilde{T}_{zz}|_{\rm ext}=4\partial_{y}\partial_{z}\tilde{\psi}=-4\partial_{z}\tilde{v}_{z}=4q\tilde{v}_{z}$ (seventh column of Fig. 11). This in turn demands a $z$-component of fluid velocity gradient (eight column of Fig. 11) and so of fluid velocity (ninth column). This finally further convects the interface, $\partial\tilde{h}/\partial t=\tilde{v}_{z}=\tfrac{1}{4}\Delta\sigma h\,|N_{2}|^{\prime}(\dot{\gamma})/\sigma^{\prime}(\dot{\gamma})$, enhancing the original interfacial tilt (tenth column of Fig. 11) with a growth rate $\omega=\tfrac{1}{4}\Delta\sigma\,|N_{2}|^{\prime}(\dot{\gamma})/\sigma^{\prime}(\dot{\gamma})$. This is indeed consistent with Eqn. 36 at zero surface tension (which is the limit in which the argument of this paragraph has been constructed). This mechanism resembles that of instabilities between layered viscoelastic fluids as studied (with different interfacial and wavevector orientations) in previous works Hinch et al. (1992); Wilson and Rallison (1997); Nghe et al. (2010).

V.6 Comparison with Tanner’s prediction

We now compare our criterion for the onset of edge fracture, written again here for clarity:

[TABLE]

with Tanner’s original prediction:

[TABLE]

In Tanner’s argument, $R$ is the radius of an (artificially) assumed initially semicircular crack in the sample edge. To establish a closer correspondence between Tanner’s criterion and ours, $R$ must now must be replaced by the dominant wavelength of instability $L_{y}$. To within $O(1)$ prefactors, the difference between Tanner’s prediction and ours then lies in replacing

[TABLE]

Given negligible air viscosity, the jump $\Delta\sigma$ in shear stress across the interface between the fluid and air simply equals the shear stress $\sigma$ in the fluid bulk, as noted above. For most fluids, in the limit of small shear rates, $N_{2}\sim-\dot{\gamma}^{2}$ and $\sigma=\eta\dot{\gamma}$. Tanner’s $|N_{2}|$ on the left hand side of of (41) then simply equals our expression on the right hand side. (In non-Brownian hard sphere suspensions, these scalings no longer hold. For example, the second normal stress scales linearly with the shear rate. It would be interesting in future work to study the edge fracture instability in that class of materials.)

The identification by Tanner of the role of $N_{2}$ in driving edge fracture was a remarkable early insight, and indeed Tanner discussed carefully a mechanism for edge fracture based on radial stress balance. However, this fortuitous agreement between the expressions given above should not be over-interpreted. Indeed, whereas our expression can be connected directly to the detailed mechanism of instability by tracing through the argument in the previous subsection, no such detailed mechanistic insight can be gained from Tanner’s prediction. Furthermore, Tanner’s prediction says nothing about the role of the shear stress in driving the instability.

At higher shear rates, the simple power law scalings $\sigma\sim\dot{\gamma}$ and $N_{2}\sim-\dot{\gamma}^{2}$ noted above no longer hold (in general), and our prediction departs from Tanner’s. Indeed (assuming a monotonic $N_{2}(\dot{\gamma})$) Tanner predicts a monotonic phase boundary $\Gamma(\dot{\gamma})$ even for a material in which the derivative of $N_{2}$ with respect to shear rate is a non-monotonic function of shear rate. This is in stark contrast to our non-monotonic one, with its re-entrant stability at high strain rates.

VI Possible mitigation strategy

Finally, and perhaps most importantly, the results of our linear stability calculation suggest a recipe via which edge fracture might be mitigated in experimental practice. We recall, in particular, that the left hand side of our onset criterion (37) contains the term $\Delta\sigma$, which is the jump in shear stress between the fluid and outside medium. In any experiment where the outside medium is air, the shear stress outside the fluid will be negligible and we simply have $\Delta\sigma=\sigma$.

Instead immersing the flow cell in an immiscible Newtonian ‘bathing fluid’ with a viscosity larger than that of air, more closely matched to that of the study-fluid, the jump $\Delta\sigma$ in shear stress between the study and bathing fluids will clearly be reduced. We explore this suggestion in Fig. 12a), showing the thresholds for onset of instability for successively increasing values of $\eta_{\rm a}$, indeed finding greater stability at larger $\eta_{\rm a}$. The dot-dashed lines in Fig. 12 show the results of a linear stability analysis performed as in Sec. V above, but generalised to include a non-zero viscosity of the outside air. Fig. 12b) explores the suppression of instability for different levels of surface tension, as a function of $\eta_{\rm a}$ and shear rate.

Another obvious strategy for mitigating edge fracture would be to try ensure as large an interfacial surface tension as possible, by suitable choice of the (Newtonian) bathing medium.

VII Conclusions

In this work, we have performed a detailed theoretical study of edge fracture, combining direct nonlinear simulations with linear stability analysis for the initial onset of edge fracture, and finding full agreement between these. We have derived a criterion of the onset of edge fracture, and provided a detailed understanding of the physical mechanism that drives the instability. Our results also suggest a new strategy via which edge fracture might potentially be mitigated in experimental practice. We also suggest that, in containing the second normal stress difference, our criterion for the onset of edge fracture may be used as a means to determine that quantity experimentally

Despite these successes, our work contains several notable limitations that should be addressed in future studies. We discuss these finally in turn.

First, we have considered only constitutive models that predict a negative second normal stress. Positive values of the $N_{2}$ were discussed as being potentially stabilising against edge fracture in Ref. Keentok and Xue (1999). It would be interesting in future work to re-do the present calculation in a constitutive model that predicts a positive $N_{2}$.

Second, we have considered only fluids with a finite terminal relaxation time, $\tau$, for which the shear stress scales as $\sigma\sim\dot{\gamma}\tau$ and the second normal stress as $N_{2}\sim-(\dot{\gamma}\tau)^{2}$ in the limit $\dot{\gamma}\tau\ll 1$. It would be interesting in future work to consider the case of non-Brownian suspensions Denn and Morris (2014), in which $N_{2}$ instead scales linearly with shear rate.

Third, we have assumed in the derivation of our criterion a base flow corresponding to a state of steady shear. Edge fracture is also widely seen in transient rheological tests such as shear startup. Future theoretical studies could profitably consider the effects of a time-dependent base state on the onset of edge fracture.

Fourth, we have considered the limit of planar Couette flow, arguing that this geometry provides a good approximation to cylindrical Couette flow in the small gap limit, and to cone-and-plate (or plate-plate) flows for small gap and large device radius. If this assumption is correct, our results should apply equally to cylindrical Couette and cone-and-plate. In view of this, an outstanding puzzle is why edge fracture is discussed much more commonly in the literature for cone-and-plate (or plate-plate) than for cylindrical Couette. One possibility is that the highly viscoelastic materials for which edge fracture occurs can only be studied in cone-and-plate (or plate-plate). Another possibility is that there is an additional feature, not captured in our analysis, that renders cone-and-plate (or plate-plate) much more highly susceptible to edge fracture than cylindrical-Couette. For example, do the hoop stresses present in cone-and-plate (or plate-plate) greatly accelerate the nonlinear dynamics of edge fracture, even while leaving the criterion for initial onset unaltered? Indeed, in plate-plate and cone-plate devices, the crack propagates inwards towards regions with smaller radius of curvature, whereas in cylindrical Couette it propagates orthogonally to the direction in which the curvature changes. A simpler explanation is that for the same size of crack, the fraction of the surface disturbed in curved Couette is smaller than in cone-plate (or plate-plate) due to the larger surface of the measuring fixture.

Fifth, we have ignored bulk instabilities that can lead to complicated secondary flows within the sample. We have thereby assumed that such instabilities are deferred until beyond the onset of edge fracture, and/or that edge fracture and bulk instabilities occur in different classes of materials. It would be interesting in future work to carry out a simulation study that allows a crossover, as a function of some material property, from a situation in which bulk instabilities first arise as a function of increasing strain rate, to one in which edge fracture dominates instead.

Sixth, and related, we have ignored the well known bulk instability in which an initially homogeneous shear flow gives way to the formation of coexisting shear bands. In recent works, it has been demonstrated theoretically that shear bands must lead inevitably to full edge fracture in some parameter regimes Skorski and Olmsted (2011); and conversely that even modest edge disturbances can lead to an apparent bulk shear banding, even in materials that would not show any true shear banding in the absence of edge effects Hemingway and Fielding (2018). Indeed, possible edge fracture in experiments concerning shear banding in entangled polymers has been discussed in Refs. Sui and McKenna (2007); Li and McKenna (2015); Schweizer and Stöckli (2008); Li et al. (2013); Wang et al. (2014); Li et al. (2014); Boukany et al. (2015). Clearly, more remains to be done to understand the interplay between bulk shear banding and surface instabilities such as edge fracture, in both steady and transient flow protocols.

Seventh, we have ignored wall slip, which occurs widely in strongly sheared entangled polymers. Indeed, wall slip remains under studied from a theoretical viewpoint, with a notable paucity of constitutive models for the layer of fluid immediately adjacent to the wall of the flow cell, as compared with those for the bulk fluid. Future theoretical studies should consider the relative dominance of and/or interplay between edge fracture and wall slip.

Finally, our calculations have ignored inertia, working throughout in the creeping flow limit of zero Reynolds number. Inertia may be an importance factor in the non-linear dynamics of edge fracture, and should be considered in future simulation studies.

*Acknowledgements – * The research leading to these results has received funding from the European Research Council under the EU’s 7th Framework Programme (FP7/2007-2013) / ERC grant number 279365. We thank Mike Cates, Gareth McKinley, Peter Olmsted, Roger Tanner and Dimitris Vlassopoulos for enjoyable discussions.