3.1. Conservation laws; prototypical viscoplastic rheology

Consider the flow of a two-dimensional, incompressible, complex fluid, describing the geometry by a Cartesian coordinate system $({\hat {x}},{\hat {y}})$ . The fluid velocity is $\boldsymbol{{\hat {u}}} = ({\hat {u}},{\hat {v}})$ . Conservation of mass and momentum demand that

(3.1) \begin{equation} \boldsymbol{\nabla \cdot {\hat {u}}} = 0, \end{equation}

and

(3.2) \begin{equation} \boldsymbol{0} = \boldsymbol{ \nabla \cdot } \boldsymbol{{\hat {\tau }}} + \hat {\boldsymbol{f}} - \boldsymbol{\nabla } {\hat {p}}, \end{equation}

where $\hat {p}$ is the pressure and $\hat {\boldsymbol{f}}$ represents any body force like gravity (we use the hat notation to denote dimensional variables; this decoration is removed using suitable scalings, to arrive at a dimensionless formulation, as outlined presently).

One of the simplest and most popular constitutive models for a yield-stress fluid of the form of (2.2) is the Herschel-Bulkley law, with tensorial formulation

(3.3) \begin{equation} \hat {\boldsymbol{\tau }} = \left (\frac {{\tau _{_{_P}}}}{\hat {\dot \gamma }} + K \hat {\dot \gamma }^{n-1} \right ) \hat {\dot {\boldsymbol{\gamma }}} \qquad \textrm {for} \qquad \hat \tau \gt {\tau _{_{_P}}}. \end{equation}

The plastic viscosity, ${\mu }(\hat {{\dot \gamma }})=K\hat {{\dot \gamma }}^{n-1}$ takes the same form as the power-law fluid, with a consistency $K$ and power-law index $n$ . The yield stress is $\tau _{_{_P}}$ ; if $\hat \tau \lt {\tau _{_{_P}}}$ , this threshold for yield is not breached and the deformation rates must all vanish, translating to $\hat {\dot \gamma }_{ij}=0$ . An important special case of (3.3) is when $n=1$ and the rate-dependent part of the stress takes the form of a constant viscous stress. This is the Bingham fluid, which is popular because of its simplicity, not its realism. Fits to data obtained in rheometers more typically generate power-law exponents $n\lt 1$ , although the power-law form itself is usually only a rougher representation of real fluid behaviour (Balmforth et al. Reference Balmforth, Frigaard and Ovarlez2014; Bonn et al. Reference Bonn, Denn, Berthier, Divoux and Manneville2017; Coussot Reference Coussot2018; Frigaard Reference Frigaard2019a ).

Importantly, when ${\hat {\tau }}\lt {\tau _{_{_P}}}$ , no relation is imposed between the stresses and strain rates. In two dimensions, the two force-balance equations (3.2) are then insufficient to determine the two independent components of $\hat {\boldsymbol \tau }$ and the pressure $\hat {p}$ . That is, the stress state is formally indeterminate over any plugged region, which is an awkward feature of the problem that will recur often in our discussion.

3.2. Scaling and dimensionless formulation

It is helpful to reformulate the equations in dimensionless form and identify key dimensionless groups. We have already identified the Bingham number (2.3) as the ratio of viscous and yield-stress scales; for the Herschel–Bulkley law, given length and velocity scales $\mathcal{L}$ and $\mathcal{U}$ , the Bingham number is

(3.4) \begin{equation} \textit{Bi} = \frac {{\tau _{_{_P}}}}{K\left ({\mathcal{U}}/{\mathcal{L}}\right )^n}, \end{equation}

with $\textit{Bi} \gg 1$ denoting the plastic limit of the problem. Alternatively, when the velocity scale is not imposed, but a body force like gravity is present, a stress scale $\mathcal{P}$ can be defined (such as $\mathcal{P} = \rho g \mathcal{L}$ ). In this circumstance, one can define a velocity scale

(3.5) \begin{equation} \mathcal{U} = \mathcal{L} \left ( \frac {\mathcal{P}}{K} \right )^{\frac 1n} \quad \left (\equiv \mathcal{L}^{1+\frac 1n} \left ( \frac {\rho g}{K} \right )^{\frac 1n} \right ), \end{equation}

where $\rho$ and $g$ denote density and gravitational acceleration and the Bingham number in (3.4) becomes $\textit{Bi}={\tau _{_{_P}}}/\mathcal{P}$ (or ${\tau _{_{_P}}}/(\rho g \mathcal{L})$ ). This parameter must remain order one if the applied stresses are to exceed the yield stress and force motion. Strictly speaking, $\textit{Bi}$ now represents the ratio of yield stress to the imposed stress, and (for gravity) is sometimes referred to as the ‘Oldroyd number’ instead. The plastic limit is then characterised by the special value of $\textit{Bi}$ for which the fluid speed approaches zero. That is, the plastic limit corresponds to the threshold for motion.

For two-dimensional flow without body forces, we write the dimensionless governing equations using either a Cartesian coordinate system $(x,y)$ or general curvilinear coordinates $(s,\eta )$ built from the arc length along some prescribed curve and its normal:

(3.6) \begin{equation} \begin{aligned} &\frac {\partial {u}}{\partial {x}}+\frac {\partial {v}}{\partial {y}} = 0, \cr &\frac {\partial {{\tau _{_{\textit{XX}}}}}}{\partial {x}} + \frac {\partial {{\tau _{_{\textit{XY}}}}}}{\partial {y}} = \frac {\partial {p}}{\partial {x}},\cr &\frac {\partial {{\tau _{_{\textit{XY}}}}}}{\partial {x}} - \frac {\partial {{\tau _{_{\textit{XX}}}}}}{\partial {y}} = \frac {\partial {p}}{\partial {y}}, \end{aligned} \qquad \qquad \begin{aligned} &\frac {\partial {U}}{\partial {s}} + (1-\kappa \eta )\frac {\partial {V}}{\partial {\eta }} - \kappa V = 0, \cr &\frac {\partial {{\tau _{_{\textit{SS}}}}}}{\partial {s}} + (1-\kappa \eta )\frac {\partial {{\tau _{_{SN}}}}}{\partial {\eta }} - 2\kappa {\tau _{_{SN}}} = \frac {\partial {p}}{\partial {s}},\cr &\frac {\partial {{\tau _{_{SN}}}}}{\partial {s}} - (1-\kappa \eta )\frac {\partial {{\tau _{_{\textit{SS}}}}}}{\partial {\eta }} + 2\kappa {\tau _{_{\textit{SS}}}} = \frac {\partial {p}}{\partial {\eta }}, \end{aligned} \end{equation}

where

(3.7) \begin{equation} \kappa = \frac {\partial {\theta }}{\partial {s}}, \end{equation}

is the curvature of the prescribed curve. Here and throughout this work, capital-letter subscripts are used to identify tensor components. The velocity field in Cartesian coordinates is $(u,v)$ ; the curvilinear relative is $(U,V)$ . The strain-rate components are given by

(3.8) \begin{equation} \begin{aligned} &{\dot \gamma _{_{\textit{XX}}}} = 2\frac {\partial {u}}{\partial {x}},\cr &{\dot \gamma _{_{\textit{XY}}}} = \frac {\partial {u}}{\partial {y}}+\frac {\partial {v}}{\partial {x}}, \end{aligned} \qquad \qquad \begin{aligned} &{\dot \gamma _{_{\textit{SS}}}} = \frac {2}{1-\kappa \eta }\left (\frac {\partial {U}}{\partial {s}}-\kappa V\right ),\cr &{\dot \gamma _{_{SN}}} = \frac {1}{1-\kappa \eta }\left (\frac {\partial {V}}{\partial {s}}+\kappa U\right )+ \frac {\partial {U}}{\partial {\eta }}, \end{aligned} \end{equation}

and the Herschel–Bulkley constitutive law is

(3.9) \begin{equation} \boldsymbol{\tau } = \left ({\dot \gamma }^{n-1} +\frac {\textit{Bi}}{{\dot \gamma }}\right ) \dot {\boldsymbol{\gamma }} \qquad \textrm {for} \qquad \tau \gt \textit{Bi}. \end{equation}

The nonlinear rheology in (3.9) presents some challenges for computing numerical solutions to (3.6)–(3.9), primarily because of the yield condition. A convenient, alternative approach that circumvents such challenges is provided by ‘regularising’ the constitutive law: (3.9) is replaced by a smoothed version with a well-defined and finite relationship between stress and strain rate everywhere. More precise techniques also exist, most notably the augmented Lagrangian approach, which provides an iterative solution to the exact rheology directly. Some details of both of these schemes, along with other numerical considerations and observations, are presented in Appendix A.

Alternatively, as in any fluid mechanical problem, asymptotic approaches based on an extreme difference in length scales provide a useful tool for interrogating viscoplastic flows (Balmforth Reference Balmforth2019). Shallow flows can be attacked with a generalisation of viscous lubrication theory, although the approach has been incorrectly denigrated in the past because of the fallacy of the so-called lubrication paradox (Balmforth & Craster Reference Balmforth and Craster1999; Hewitt & Balmforth Reference Hewitt and Balmforth2012). More generally, and as illustrated earlier, viscous deformation becomes restricted to narrow boundary layers in the plastic limit, which are amenable to boundary-layer analysis (Oldroyd Reference Oldroyd1947b ; Balmforth et al. Reference Balmforth, Craster, Hewitt, Hormozi and Maleki2017). Such asymptotic approaches can become technically involved; the interested reader can find further details in Appendix C.

3.3. Extremum principles and minimum dissipation

The (dimensionless) governing equations can be cast into a variational form that identifies extremum principles satisfied by the velocity and stress fields (see Prager (Reference Prager1954) and, for example, Frigaard (Reference Frigaard2019b )). We quote these principles in settings for which we ignore any body forces (i.e. gravity) and the velocity field is prescribed on all the boundaries. The principles therefore apply to the steady flow problems considered in § 4 and § 5, and provide a useful constraint associated with the net dissipation rate on the solutions in the plastic limit.

The first principle is that, of all the divergence-free velocity fields $\boldsymbol{v}$ that match up with the velocities along the boundaries, the actual velocity field $\boldsymbol{u}$ minimises the functional

(3.10) \begin{equation} \Phi [\boldsymbol{v}] = \int \int _{\mathcal{D}} \left \{ \textit{Bi} + \frac {[{\dot \gamma }(\boldsymbol{v})]^{n}}{n+1}\right \} {\dot \gamma }(\boldsymbol{v}) \; \textrm {d} x\textrm {d} y, \end{equation}

where $\mathcal{D}$ denotes the domain occupied by the fluid and ${\dot \gamma }(\boldsymbol{v})$ denotes the strain-rate invariant built from $\boldsymbol{v}$ . In this minimisation, the trial velocity fields are not linked through the constitutive law to any stress field.

The second principle is that, of all the possible stress fields $\tilde {\boldsymbol{\sigma }}$ that satisfy the force-balance equations, independently of the constitutive law, the actual stress $\boldsymbol{\sigma }$ maximises

(3.11) \begin{equation} \Psi [\tilde {\boldsymbol{\sigma }}] = \int _{\partial \mathcal{D}} \boldsymbol{u \cdot } \tilde {\boldsymbol{\sigma }} \boldsymbol{\cdot } \boldsymbol{\hat {n}} \; \textrm {d} s - \frac {n}{n+1} \int \int _{\mathcal{D}} \left \{ \textrm {Max}\left [0, \tau (\tilde {\boldsymbol{\sigma }}) - \textit{Bi} \right ] \right \}^{1+\frac 1n} \; \textrm {d} x\textrm {d} y, \end{equation}

where $\partial \mathcal{D}$ denotes the boundary of the fluid domain, with unit outward normal $\boldsymbol{\hat {n}}$ and infinitesimal arc length $\textrm {d} s$ , and $\tau (\tilde {\boldsymbol{\sigma }})$ denotes the second invariant of the trial stress field.

Necessarily, when the extrema are achieved, we also have that $\Phi [\boldsymbol{u}]=\Psi [ \boldsymbol{\sigma }]$ , which corresponds to the integral constraint

(3.12) \begin{equation} \int _{\partial \mathcal{D}} \boldsymbol{u \cdot } \boldsymbol{\sigma } \boldsymbol{\cdot } \boldsymbol{\hat {n}} \; \textrm {d} s = \int \int _{\mathcal{D}} \tau {\dot \gamma } \; \textrm {d} x\textrm {d} y \equiv \mathcal{E}, \end{equation}

and is simply an expression of the fact that the power supplied through boundaries (the left-hand side) must balance the net dissipation rate ( $\mathcal{E}$ , the second integral). Importantly, in the plastic limit, $\tau {\dot \gamma } \to \textit{Bi}{\dot \gamma }$ and the velocity minimisation corresponds to selecting the trial velocity field that has the least net rate of dissipation. In addition, if a body force is present (cf. § 6), the work performed by that force must also be added to the variational problems and to (3.12) (Prager Reference Prager1954; Frigaard Reference Frigaard2019b ).

3.4. Plasticity theory

For $\textit{Bi}\gg 1$ , the theory of perfect plasticity applies (Hill Reference Hill1950; Prager & Hodge Reference Prager and Hodge1951); key concepts for generating solutions in this limit are provided below. Additional details of slipline analysis, and a sample construction of a relevant slipline pattern, can be found in Appendix B.

3.4.1. Sliplines

In terms of our two-dimensional Cartesian coordinate system, the constitutive law in the ideal plastic limit reduces to

(3.13) \begin{equation} {\tau ^2_{_{\textit{XY}}}}+{\tau ^2_{_{\textit{XX}}}}=\textit{Bi}^2, \end{equation}

which motivates the definition of the new variable $\theta$ in

(3.14) \begin{equation} ({\tau _{_{\textit{XX}}}},{\tau _{_{\textit{XY}}}}) = \frac {\textit{Bi}}{{\dot \gamma }}({\dot \gamma _{_{\textit{XX}}}},{\dot \gamma _{_{\textit{XY}}}}) = \textit{Bi}(-\sin 2\theta, \cos 2\theta ). \end{equation}

We may then manipulate the force-balance equations (3.2) into the hyperbolic pairing

(3.15) \begin{equation} \left (\frac {\partial }{\partial {y}}+\cot \theta \frac {\partial }{\partial {x}}\right )(p + 2 \textit{Bi} \theta ) = 0 \qquad \& \qquad \left (\frac {\partial }{\partial {y}}-\tan \theta \frac {\partial }{\partial {x}}\right )(p - 2 \textit{Bi} \theta ) = 0 . \end{equation}

That is, we have the Riemann invariant $p+2\textit{Bi}\theta$ along the curves with $\textrm {d} y/\textrm {d} x = \tan \theta$ ; the curves with $\textrm {d} y/\textrm {d} x = -\cot \theta$ have $p-2\textit{Bi}\theta = \textrm{constant}$ . These characteristic curves are the ‘sliplines’ of classical plasticity theory (Hill Reference Hill1950; Prager & Hodge Reference Prager and Hodge1951). Evidently, the two sets of curves are orthogonal to one another, and are commonly referred to as the $\alpha$ and $\beta $ -lines. Convention takes the two sliplines to form a local right-handed coordinate system.

The sliplines form a net that spans all the yielded regions (two examples are shown in figure 2). A number of properties of the network follow from (3.15). Two, in particular, are often quoted and termed ‘Hencky’s rules’; we describe these more thoroughly in Appendix B.1. One important consequence of them is that when a slipline of one family (i.e. $\alpha$ or $\beta$ ) is straight, then all the other lines of the same family must also be straight if sliplines from the other slipline family connect them together. This relatively strong geometrical constraint impacts the structure of the yielded regions and provides a pathway to solving the slipline equations analytically in some problems. In particular, it implies that when one slipline is straight, the acceptable slipline fields include square checkerboards and centred circular fans, as illustrated in figure 2(b), or a set of involutes to a curve such as a circle (as we encounter later).

In general, however, the utility of sliplines is limited because rigid plugs can be present, over which (3.13) no longer holds, and any slipline can be taken to be a yield surface bounding these plugs (Hill (Reference Hill1950); Prager & Hodge (Reference Prager and Hodge1951); see also figure 2). Worse, discontinuities can also appear in the stress field. The sliplines do not remain continuous across these features, implying scars can disfigure the slipline pattern. The (total) normal and shear stresses remain continuous across any such discontinuity, but there can be a jump in the tangential normal stress (i.e. ${\sigma _{_{\textit{SS}}}}={\tau _{_{\textit{SS}}}}-p$ , if the curvilinear coordinate system is aligned with the discontinuity). These conditions translate to the requirement that the sliplines meet the discontinuity at equal angles.

Figure 2. Illustrations of the slipline field, and Prandtl’s (a) cycloid and (b) punch solutions. The shaded regions indicate plugs. In (b), the dotted lines show sample streamlines.

Complicating matters further is that, in many problems, boundary conditions are provided on the velocity, not the stress. One cannot then construct the slipline field without a consideration of the velocity field. The components of the velocity field along the sliplines, $(v_\alpha, v_\beta )$ , satisfy the so-called Geiringer equations

(3.16a,b) \begin{equation} \frac {\partial {v_\alpha }}{\partial {s_\alpha }}=v_\beta \frac {\partial {\theta }}{\partial {s_\alpha }} \qquad \text{and} \qquad \frac {\partial {v_\beta }}{\partial {s_\beta }}=-v_\alpha \frac {\partial {\theta }}{\partial {s_\beta }}, \end{equation}

which follow from the constraint of incompressibility (3.1). These equations imply that the plastic flow lines (i.e. the characteristics of the velocity field) coincide with the sliplines. Therefore, if the velocity field is directed along one set of sliplines, these curves also correspond to streamlines and the velocity along them remains fixed. Again, a complication is that the ideal plastic velocity field need not remain continuous. In fact, as implied by the name, a slipline can be taken as the locus of a jump in tangential velocity. In other words, velocity discontinuities can be present, but they are aligned with sliplines. Figure 2(b) also plots a selection of streamlines, corresponding to the slipline field shown. Here, fluid enters the domain at an angle to the $45^\circ$ checkerboard of sliplines. Once fluid reaches the fan, however, slip occurs, allowing the velocity to re-align with the circular arcs of the fan. Thereafter, one of the slipline families corresponds to streamlines.

Although the slipline field can be used to directly construct a range of perfectly plastic solutions, the difficulties implicit in the construction can make such an approach unwieldy and it is often not possible to proceed analytically. While techniques have been developed to ease numerical constructions of the slipline field (Collins Reference Collins1982), our perspective here is different. With the advent of augmented Lagrangian numerical techniques, one can solve two-dimensional problems of viscoplastic steady flow with relative ease, even close to the plastic limit. This circumvents any need to construct numerical slipline fields directly. Instead, with a numerical viscoplastic solution in hand, the slipline analysis provides a convenient tool to diagnose the solution structure, pointing to the existence of any analytical special cases and providing a means to verify the fidelity of computations. We illustrate these ideas later, using some specific flow configurations.

3.4.2. Exact solutions and constructions

Two classical exact solutions to the slipline problem date back to Prandtl and are illustrated in figure 2. The solution illustrated in figure 2(b) consists of a patchwork of circular fans and triangular checkerboards. The rightmost sliplines (a $\beta $ -line for the bottom half, and an $\alpha $ -line for the top) provide a yield surface. This construction was offered by Prandtl as a slipline solution for plastic indentation, following along the lines of traditional developments for the indentation of a punch in linear elasticity. Prandtl’s solution (rotated by ninety degrees) is commonly used in soil mechanics for the failure of a foundation above a cohesive soil.

A second solution offered by Prandtl applies to the compression of a plastic layer, or a squeeze flow. In this case the sliplines are cycloids: if $\theta =\theta (y)$ and $p=\textit{Bi}[ -\Gamma x + \Phi (y)]$ , then the force-balance equations reduce to

(3.17) \begin{equation} \frac{\partial}{\partial y} (\cos 2\theta) = -\Gamma \quad \& \quad \frac{\partial}{\partial y} ( \Phi + \sin 2 \theta) = 0, \end{equation}

If $\theta =0$ at $y=y_0$ and $\theta = (1/4)\pi$ at $y=y_a$ , then the $\alpha $ -lines follow from

(3.18) \begin{equation} \cos 2\theta = \frac {y-y_a}{y_0-y_a}, \quad \frac {\textrm {d} y}{\textrm {d} x} = \tan \theta = \sqrt {\frac {y_0-y}{y+y_0-2y_a}}, \quad \Gamma = -(y_0-y_a)^{-1}. \end{equation}

Hence

(3.19) \begin{equation} 2\tan ^{-1} \sqrt {\frac {2}{v}-1} - \sqrt {v(2-v)} + 1 - \frac {1}{2}\pi = \frac {x_0-x}{y_0-y_a}, \quad v=\frac {y_0-y}{y_0-y_a}, \end{equation}

if $x=x_0$ for $y=y_a$ , which is the equation for a cycloid. Note that this curve is scale free, in the sense that any scaling of $y$ and $x$ leaves the equation unchanged (as long as $\Gamma$ can be scaled suitably too). The resulting slipline field is actually what is illustrated in figure 2(a).

Note that Prandtl’s two solutions illustrate an important feature of the plasticity problem: the squeezing walls generating the cycloid solutions in figure 2(a) are assumed to be ‘fully rough’, a condition demands that the fluid slides along the wall with a local shear stress that equals the yield stress. For a viscoplastic fluid, such slip prompts the appearance of a thin wall boundary layer over which the tangential velocity is brought back to zero. In other words, fully rough boundary conditions on the walls translate to no slip, and the slipline angle is either $\theta =0$ or $\theta =(1/2)\pi$ . Conversely, the surface on the left in figure 2(b) is free, implying that the shear stress must vanish, and $\theta =\pm (1/4)\pi$ . Such a condition on the slipline angle also applies along lines of symmetry, such as the midline between the squeezing plates in figure 2(b). The $y$ -axis can further be taken as another line of symmetry at the centre of the squeeze flow, if the shaded region is assumed to plug up with yield surfaces along the bordering sliplines. The correspondence between no-slip or free-slip boundary conditions and the fixing of the slipline angle recurs in the examples we present below.

More generally, the slipline equations (when expressed in terms of the two local radii of curvature, or some other geometrical variables; Appendix B.1) and Geiringer’s equations reduce to the telegraph equation under a change of independent variables (from arc length to an angle parameterisation; Hill (Reference Hill1950)). Classical solutions to that problem then allow one to write down exact solutions to the plasticity problem in terms of single integrals if information is provided along a known slipline. Even if such information is not provided, one can still use this solution as the basis of a numerical scheme (Dewhurst & Collins Reference Dewhurst and Collins1973; Collins Reference Collins1982). This strategy significantly expanded the library of slipline solutions available in the plasticity literature (e.g. Collins (Reference Collins1970, Reference Collins1982)), a resource that largely remains untapped in analyses of viscoplastic fluid flows near the plastic limit. In fact, this technique leads to semi-analytical solutions, one of which is closely related to the slipline problem discussed in § 4 (Hill Reference Hill1950; Chakrabarty Reference Chakrabarty1979). Other exact solutions or constructions follow from searching for self-similar solutions (e.g. Taylor-West & Hogg (Reference Taylor-West and Hogg2021, Reference Taylor-West and Hogg2023)).

4.1. Narrow domains; breaking the incoming plug

A suite of numerical solutions for the jet problem are displayed in figure 4. Here, the Bingham number is fixed close to the plastic limit ( $\textit{Bi}=2048$ ), the cross-stream height is set at $\ell _y= 3/2$ and the downstream length $\ell _x$ is varied. For small domain lengths $\ell _x$ , flow enters and exits the domain as a moving plug of almost uniform width bordered by two thin shear layers (figure 4 c). This localisation of fluid deformation by the yield stress, and the creation of extensive blockages, is rather different from traditional Stokes flow, for which fluid deforms throughout the domain (and is equivalent to the localisation of viscoplastic flow around the disk in figure 1). The shear layers have a structure that can be predicted by Oldroyd’s boundary-layer theory (see figure 5 a,c), as outlined in § C.1.1. That analysis indicates that the net dissipation $\mathcal{E}$ over the upper half of the domain is approximately $\ell _x\textit{Bi}$ (cf. (C9)) and the area of the corresponding yielded region is given by (C8); both diagnostics are plotted in figure 4(a,b) for the full series of computations.

Figure 4. Viscoplastic jets for $(\textit{Bi},\ell _y)=(2048, {3}/{2})$ and varying $\ell _x$ . (a) Net dissipation rate and (b) area of the yielded regions for $y\gt 0$ , plotted against $\ell _x$ . (c)–(i) Density maps of $\log _{10}{\dot \gamma }$ on the $(x,y)$ -plane for the values of $\ell _x$ indicated. For each $\ell _x$ , two solutions are shown: for the upper (filled blue stars in (a,b)), no-slip conditions are applied at $y=\pm \ell _y$ ; for the lower (open red squares in (a,b)), free-slip conditions are used instead. The insets in (a,b) show transitional, free-slip solutions at $\ell _x=1.275$ and $\ell _x=1.9$ . Open (blue) stars in (a)–(b) represent solutions with $\ell _y\gg 3/2$ , for which the yielded regions do not touch the top and bottom boundaries. The grey lines in (a,b) show the predictions in (C8)–(C9). The solid blue and red lines in (a) show predictions for $\mathcal{E}$ based on the slipline fields in figure 6(a). The vertical dashed lines at $\ell _x=({\pi }/{4})+(1/2)$ and $\ell _x\approx 1.9$ indicate the flow pattern transitions arising when one of the plugs breaks; that at $\ell _x=3.824$ indicates where the slipline solution predicts that the yielded region collides with the walls at $y=\pm \ell _y$ .

Figure 5. (a,c) The viscoplastic shear layer for $\ell _x=1$ (cf. figure 4 c), and (b,d) the wall boundary layer for $\ell _x=6$ (cf. figure 4 i). The vertical dashed lines in (a,b) show the cuts at which the horizontal velocity profiles of (b,d) are drawn. The dashed curves in (a) shows the yield surfaces predicted by shear-layer theory (§ C.1.1). In (c), the velocity profiles (shown by red lines) are collapsed by shifting each by $(1/2)$ , then plotting them against $\zeta =(y-(1/2))/Y$ , where $Y$ is the local half-width of the shear layer. In (d), the profiles are collapsed by scaling $u$ by the plug speed $u_*$ , and then plotting the curves against $\zeta =(y- (3/2)/\eta _*$ , where $\eta _*$ is the local boundary-layer thickness. The blue dots in (c,d) show the predictions $(1/2)-(1/4)\zeta (3-\zeta ^2)$ (§§ C.1.1) and $(2+\zeta )\zeta$ (C.1.2).

Figure 6. The two slipline fields for $\ell _x=3$ (cf. figure 4 g) for the (b) no-slip and (c) free-slip cases. Additional diagnostics are plotted in (c,d). In (a,b), half of the domain shows the numerical solution for $\log _{10}{\dot \gamma }$ along with reconstructions of the curves of constant $p\pm 2\textit{Bi}\theta$ , and the insets display the geometry of the yielded regions and some special points of the slipline construction. For (c), we plot the $x-$ position of point $C$ ( $x_{_C}$ ) and the slipline angle $(\theta _{_E}$ ) at point $E$ against $\ell _x$ for no-slip solutions with $\ell _y=3$ . In (d), the scaled dissipation rate $\mathcal{E}/\textit{Bi}$ , $\theta _{_E}$ and $x_{_C}$ are plotted against $\ell _y$ , for free-slip solutions with $\ell _y=1.5$ . The stars in (c,d) display corresponding results from the numerical solutions (with $x_{_C}$ measured from the right-hand border of the plastic region in (a), then from the centre of the shear layer at $y=\ell _y$ in (b) and $\theta _{_E}$ from the centre of the upper shear layer at $x=(1/4)$ ).

As the domain becomes longer, the stresses exerted on the moving plug increase; eventually, this increase causes the plug to break. A region of almost perfectly plastic deformation then appears near the inlet that allows the jet to adjust its width (figure 4 d,e). This widening reduces the net dissipation rate below that expected for an unbroken plug (figure 4 a), and abruptly increases the area of yielded plug (figure 4 b). Evidently, the widening of the jet is sufficient to reduce the velocity jumps across the shear layers, and the net dissipation rate within them, despite the lengthening of those layers and the additional dissipation incurred over the plastic region.

As illustrated in figure 6(a,b), the perfectly plastic flow can be described by slipline theory: much as for Prandtl’s punch, two centred fans emerge from the edges of the inlet, with a triangular wedge between them. Indeed, the intrusion can be viewed as an indentation problem, with different boundary conditions applying on the left wall. The fans extend outwards and meet along the centreline of the jet; the slipline pattern then continues out to the right, with all the sliplines becoming curved due to the symmetry condition along $y=0$ ( $|\theta |=({1}/{4})\pi$ ). The plastic region terminates to the right along a final slipline that fronts the widened moving plug. Shear layers continue to buffer the jet from above and below. More details of the construction of the slipline field are given in Appendix B. Note that the checkerboard of sliplines filling the wedge at the inlet, in combination with the rigid-body inflow velocity, allows that region to solidify into another example of an serendipitous plug. A different example, in which the inflow is not rigid and the checkerboard remains fully plastic, is presented below in § 4.5.

Figure 6(b,c) also illustrates the slipline field reconstructed directly from the numerical solutions, by plotting the curves of constant $p\pm \textit{Bi} \tan ^{-1}({\tau _{_{\textit{XX}}}}/{\tau _{_{\textit{XY}}}})$ (cf. equations (3.15)). This more practical strategy to determine the sliplines successfully confirms most of the details of the slipline field, but breaks down in the shear layers and plugs, and fails to trace all the sliplines in the fans. An example is which the reconstruction is even poorer will be encountered later in figure 10(b).

From the slipline solution, one can construct the force balance on the widened plug, which controls the size of the plastic region. Appendix B.2 summarises the details of this calculation, which provides predictions for a number of key characteristics of the flow pattern, as plotted in figures 4(a) and 6(b,c). For example, the net dissipation rate can be determined, either by a direct calculation using Geiringer’s equations (3.16), or by exploiting the power balance in (3.12) and calculating the power input along the borders of the slipline pattern.

The analysis of the slipline patterns also predicts the domain length at which the uninterrupted moving plug breaks: the critical value is $\ell _x\sim (1/4)\pi +(1/2)$ , as seen in figure 4(a,b). Because sudden jumps arise in the properties of the flow pattern at this transition (including the area of the yielded region), owing to the breakage of the plug, we refer to this type of bifurcation as ‘first order’, following the terminology of phase-transition theory. As we shall see below, continuous, or ‘second-order’ bifurcations are also possible, that do not involve the breakage of any plug. Note that the transition for the numerical solutions in figure 4 is slightly shifted and ‘blurred’: the finite Bingham number used for the computations adjusts the domain length at which the plug breaks and appears to lead to flow patterns with mixed character over a narrow window surrounding that transition (see the inset in figure 4 a). However, it is hard to tell whether this latter feature is genuine or an artefact of the numerical scheme, the solution having perhaps failed to properly converge (see Appendix A).

In fact, there are actually two suites of computations shown in figure 4: one with no-slip walls at $y=\pm \ell _y$ , the other with free slip there. As long as $\ell _x$ remains less than about $1.9$ , the two sets of solutions are identical. In other words, two different problems have the same solution. This coincidence arises because of the cloaking effect of the plug that intervenes between the flowing region and the wall for $\ell _x \lesssim 1.9$ : if the wall lies inside a plug, that entire region is rigid and so the conditions applying upon the wall, and even its shape are irrelevant. This notion of cloaking dates back at least to Oldroyd (see Balmforth, Craster & Hewitt (Reference Balmforth, Craster and Hewitt2021)), and has been considered in detail for the motion of differently shaped particles encased in viscoplastic fluid (Chaparian & Frigaard Reference Chaparian and Frigaard2017a ,Reference Chaparian and Frigaard b ). The difference between the two cases that emerges for longer domains arises because the plug breaks for free-slip walls at $\ell _x\approx 1.9$ , decloaking the wall (see below in § 4.3).

4.2. Wide domains; no-slip walls

For longer domain lengths, the solutions depend on the boundary condition applied along the top and bottom boundaries (see e.g. figure 4 a,b,f). With a no-slip condition, the jet continues to widen gradually until it meets the walls at $y=\pm \ell _y$ , for domains lengths near $\ell _x\approx 3.6$ . Thereafter, as in figure 4(h), the bordering shear layers narrow to become described by the viscoplastic wall-boundary-layer theory of § C.1.2. That wall layer is predicted to have constant thickness and a parabolic velocity profile. The former feature is well reproduced by the numerical solution, as illustrated in figure 5(b). The parabolic velocity profile is less well reproduced, however; see figure 5(d). The discrepancy here arises not because the boundary-layer flow takes a different form to that predicted, but because of the algorithm used to compute the solutions: the algorithm exploits a Fourier cosine series in $y$ , and the solutions only satisfy a no-slip condition at $y=\ell _y$ after the addition of a numerical ‘sponge layer’ at the top boundary (spanning $\ell _y\lt y\lt 1.1\ell _y$ ). Over the sponge, the yield stress and viscosity are artificially increased to suppress motion. In finer detail, however, this device does not fully incorporate no slip, leaving a velocity profile with noticeable departures from parabolic form. In other words, the boundary-layer theory exposes an artefact in the numerical solution.

When the yielded region reaches the no-slip wall near $\ell _x\approx 3.6$ , the containment of the flow by the no-slip boundaries limits the yielded area and enhances the dissipation rate above what would have been incurred had there been no collision with the wall. Both feature are illustrated in figure 4(a,b), in which $\mathcal{E}$ and the yielded area are compared with corresponding data for solutions with larger $\ell _y$ . The collision of the yielded region with the wall indicates a second transition in flow pattern for the no-slip problem. In this case, however, no plug breaks and the pattern properties change continuously. The transition is therefore second-order.

According to the slipline solution, the widened jet meets the walls when $\ell _x\approx 3.824$ for this specific value of $\ell _y$ . Because the shear layers have finite thickness in the numerical solutions, the transition occurs at the slightly smaller value of $\ell _x\approx 3.6$ . Once the yielded region collides with the wall, its structure becomes fixed (see figure 4 h,i), with the lengthening of the domain in $x$ serving only to impact the length and thickness of the wall layer.

4.3. Wide domains; free-slip walls

For free-slip walls, the solutions develop differently with increasing domain length: at $\ell _x\approx 1.9$ , the yielded region widens abruptly, with the moving plug then filling the entire cross-stream domain (e.g. figure 4 f). The sudden jump in yielded area indicates that the transition is again first-order. In this case, the transition arises because the stationary plug that cloaks the wall for smaller $\ell _x$ breaks suddenly. This de-cloaking event leaves a pattern for $\ell _x\gt 1.9$ containing sliplines that meet the wall at $45^\circ$ (see figure 4 f); the new slipline pattern is shown in more detail in figure 6(b). The inset of figure 4(b) further illustrates how the transition again becomes blurred for the finite- $\textit{Bi}$ computation, with patterns of mixed character appearing once more.

In regard to the minimisation of $\mathcal{E}$ , the widening of the moving plug incurs no additional dissipation at the top and bottom boundaries in view of the free-slip condition that holds there; the plug now freely slides. This means that, once the domain becomes sufficiently wide, the dissipation rate can be kept lower by eliminating the shear layers of the moving plug, despite the broadening of the perfectly plastic flow region, paving the way for the flow transition at $\ell _x\approx 1.9$ . The switch is reflected by a kink in the plot of $\mathcal{E}$ near $\ell _x\approx 1.9$ (figure 4 a), and a corresponding jump in the yielded area (figure 4(b)). Both metrics become independent of $\ell _x$ beyond the transition because the sliding plug shields the extended yielded zone from whatever lies to the right, and can accommodate any increase in domain length (see figure 4 f– 4(i)).

Note that the density plots of the strain rate in figure 4 display substructure beyond the identification of the yielded regions for both no-slip and free-slip top walls: the shear and boundary layers bordering the plastically deforming region are highlighted by elevated levels of $\dot \gamma$ . In addition, enhanced dissipation arises at the junctions between different sections of the slipline pattern due to the smoothing of the associated discontinuities in the plastic velocity field. Altogether, the viscoplastic shear and boundary layers in combination with the structure of the slipline field rationalise the form of the flow pattern.

4.4. Further details of the plugs; stress indeterminacy

Figure 7. The scaled stress component ${\tau _{_{\textit{XX}}}}/\textit{Bi}$ for (a) $\ell _x=1$ , (b) $5/4$ and (c) $ 7/5$ , showing $y\gt 0$ . The upper density plot is computed with the augmented Lagrangian scheme; the lower one exploits a regularisation of the constitutive law (Papanastasiou (A1), with $\beta =10^{4}$ ). The red contours indicate where $\tau =\textit{Bi}$ ; the dashed green contours identify curves of constant ${\dot \gamma }=10^{-3}$ and $10^{-4}$ . Cuts of ${\tau _{_{\textit{XX}}}}/\textit{Bi}$ along $x=(1/3)\ell _x$ and $(2/3)\ell _x$ are plotted to the right of the density maps (upper and lower panels, respectively); cuts along $y=(1/4)$ and $3/4$ are plotted below and above (respectively). In these cuts, the results with the augmented Lagrangian algorithm are plotted as solid blue, and for regularisation with dashed red.

Over the plugs in figure 4, the indeterminacy of the stress state implies that multiple solutions should be possible, even though the precise numerical algorithm employed selects a particular one. This feature is illustrated in figure 7 where computations with the augmented Lagrangian algorithm for $\ell _x=1$ , $5/4$ and $ 7/5$ are compared with corresponding solutions computed using the popular Papanastasiou regularisation (see Appendix A). The two agree over the yielded regions, provided the regularisation parameter $\beta$ is chosen to be sufficiently large. In panels (a,b), the yielded region consists of only the shear layer; for the example in panel (c), the yielded region is larger, the central plug having broken.

Over the plugs the two stress solutions do not agree, because the augmented Lagrangian scheme converges to a stress–velocity solution that does not satisfy the constitutive law of the generalised Newtonian fluid employed by the regularised scheme. However, near the bifurcation at $\ell _x\sim (1/4)\pi +(1/2)$ (panel (b)), the range of possible stress states evidently becomes more restricted. As a consequence, the two stress solutions align more closely over the impending yielded region. At bifurcation, the range of possible stress solutions narrows to the single, unique, yielded solution, which correspond to the mode of failure of the moving plug.

Figure 7 also highlights two other practical issues with the computations: for the augmented Lagrangian scheme, the yield surface predicted by $\tau =\textit{Bi}$ can be somewhat rough; it is arguably better to use a contour of a small value of $\log _{10}{\dot \gamma }$ . Indeed, we have followed this practice everywhere but in figure 7. On the other hand, with a regularisation of the viscosity, strain rates of order $\beta ^{-1}$ are not relevant to a yield-stress fluid and it can be awkward to detect true yield surfaces. For this algorithm, figure 7 demonstrates that it is better to use either $\tau =\textit{Bi}$ or a contour of constant $1\gg {\dot \gamma }\gg \beta ^{-1}$ to locate effective yield surfaces. Even with such choices, however, the scheme fails to capture all of the yield surface in figure 7(b). Other practical issues with regularising the constitutive model are discussed by Frigaard & Nouar (Reference Frigaard and Nouar2005).

4.5. Indentation; the appearance of stress discontinuities

Prandtl’s punch problem in figure 2(b) corresponds to a minor variation on the jet configuration, in which we apply different boundary conditions on the left boundary and take the domain to be sufficiently large that fluid plugs up well before reaching all the other boundaries. In fact, a simple periodic version of Prandtl’s punch follows from adopting free-slip conditions at $y=\pm \ell _y$ and setting $u(0,y)=u_0(y)$ , where $u_0(y)$ is a $2\ell _y-$ periodic function with zero average. Figure 8 presents solutions for such a problem, in which $u_0(y)$ is prescribed as either a square wave or a sine function. For large $\textit{Bi}$ , we arrive at solutions for which we recognise a periodic repetition of the fans of Prandtl’s slipline pattern (panels (a,b), left-hand density plots). Much further from the plastic limit, at lower $\textit{Bi}$ , some of the features remain (panels (a,b), right-hand density plots).

Figure 8. Indentation solutions for $u(0,y)=u_0(y)$ with (a) $u_0=$ sgn $((1/2)-|y|)$ and (b) $u_0=\sin y$ (and symmetry conditions along $y=\pm 1$ and $x=\ell _x$ ). For each case, solutions for $\textit{Bi}=1$ and $2048$ are shown, plotting $\log _{10}{\dot \gamma }$ over the $(x,y)$ -plane (with the same colour map). The green lines show sample streamlines; sliplines are plotted in black. In (c), the scaled stress component ${\tau _{_{\textit{XX}}}}/\textit{Bi}$ is shown for both solutions with $\textit{Bi}=2048$ ; the thicker (red) contour shows the yield surfaces. Cuts of ${\tau _{_{\textit{XX}}}}/\textit{Bi}$ are plotted to the left, taken from the vertical dotted lines on the density maps; the blue lines are for the sinusoid, red for the square wave. (d) A similar plot for the solutions with $\textit{Bi}=1$ .

Note that the velocity field over the triangular checkerboards in Prandtl’s solution is uniform, and so these regions are consistent with plugs. This remains true for the periodic version shown in figure 8(a) which uses the square-wave ‘indentation’, $u_0(y)=$ sgn $((1/2)-|y|)$ . By contrast, when $u_0(y)=\sin y$ (figure 8 b), the velocity field over those regions is sheared and the checkerboards then add to the yielded area. In other words, the checkerboards again provide settings for serendipitous plugs.

Aside from this difference, two solutions in figure 8 evidently possess the same stress solution for $\textit{Bi}\gg 1$ (see figure 8 c). This demonstrates another novel feature of viscoplastic flow near the plastic field: there can be multiple solutions with the same stress solution but different velocity field. This form of non-uniqueness is different from the cloaking effect of a plug, as it arises for yielded flow. The multiplicity is also only a feature of the plastic limit: for lower Bingham number (at $\textit{Bi}=1$ in figure 8 d), viscous effects distinguish the two solutions.

The indentation problem illustrates yet another feature of the plastic limit: Prandtl’s solution demands that $\theta =\pm \frac \pi 4$ , or ${\tau _{_{\textit{XY}}}}=0$ and $|{\tau _{_{\textit{XX}}}}|=\textit{Bi}$ , along the left wall. In fact, as seen in figure 8(c), $\tau _{_{\textit{XX}}}$ jumps from $-\textit{Bi}$ to $\textit{Bi}$ as one passes through the centre of the fans. For the square-wave example, this stress distribution and jump is not significant, as the entire left boundary, barring the centres of the fans, lies within the plugged-up checkerboards. For the sinusoid, however, the velocity field through the checkerboards is inconsistent with rigid-body motion, leaving them yielded. But the velocity boundary condition and continuity demand $v=\partial v/\partial y = \partial u/\partial x=0 \quad {\rm at} \quad x=0$ . If the left-hand boundary were yielded, the last of these conditions implies that ${\tau _{_{\textit{XX}}}}=0$ , in disagreement with Prandtl’s solution. Thus, the stress field must contain a discontinuity along $x=0$ to rescue the slipline solution. Indeed, the cuts shown to the left of figure 8(c) provide numerical evidence that $\tau _{_{\textit{XX}}}$ suddenly jumps from $\pm \textit{Bi}$ in $x\gt 0$ to zero at $x=0$ . As mentioned in § 3.4, such features are acceptable for a perfectly plastic material.

What is more surprising, however, is that the numerical solutions in figure 8 are not computed for the perfectly plastic problem, but with finite $\textit{Bi}$ . One might therefore imagine that the stress discontinuities are not true features of these numerical solutions, but become smoothed in some way by the addition of the viscous stress to the constitutive law. In fact, this is not the case: the solutions for $\textit{Bi}=1$ as well as those for $\textit{Bi}=2048$ appear to share genuine stress discontinuities along $x=0$ . In other words, the viscous stress does not smooth out these feature, in conflict with one’s expectations from Stokes flow (where the strain rates and stresses are continuous for non-singular solutions). A closer inspection at figure 8(c,d) reveals how this is apparently accomplished: a yield surface with ${\dot \gamma }=0$ appears along $x=0$ . Because of this yield line, the stress need no longer converge to ${\tau _{_{\textit{XX}}}}=0$ at $x=0$ for any $\textit{Bi}$ .

Stress discontinuities can also emerge at finite $\textit{Bi}$ within the flow domain. For example, Hewitt & Balmforth (Reference Hewitt and Balmforth2017) considered the viscoplastic version of Taylor’s swimming sheet, which is closely related to the indentation problems in figure 8, except that the left-hand boundary is no longer straight. At both large and order-one Bingham numbers, Hewitt & Balmforth (Reference Hewitt and Balmforth2017) present solutions with stress discontinuities.

5.1. Pure translation

Solutions for flow around a disk in pure translation ( ${u_{W}}=1$ , ${U_{W}}=0$ ) with three different annular gaps, $\unicode{x1D6E5} =\ell -1$ , are shown in figure 10. For a sufficiently wide gap, the fluid plugs up further away from the inner cylinder, rendering the outer cylinder irrelevant and the flow equivalent to that around an isolated disk, as illustrated earlier in figure 1. Randolph & Houlsby (Reference Randolph and Houlsby1984) provided a slipline solution for the perfectly plastic version of this problem which is compared with the numerical solution (at finite $\textit{Bi}$ ) in figure 10(a). The slipline field contains triangular checkerboards attached to the front and back of the disk, with semicircular fans positioned to either side; between the fans and checkerboards, the sliplines are involutes of the disk. As mentioned earlier, the velocity field across the checkerboards is consistent with solid-body motion, rendering these regions into serendipitous plugs in the numerical solution. The cores of the semicircular fans also solidify into residual plugs.

Figure 10. Translating cylinders for varying outer radius $\ell$ (as indicated), showing density maps of $\log _{10}{\dot \gamma }$ along with sample streamlines (green). A selection of sliplines from Randolph & Houlsby’s solution are included in (a) (green lines), and from reconstructions using the numerical stress field in (b,c) (black lines). (d) Net dissipation rate $\mathcal{E}/\textit{Bi}$ and (e) yielded area are plotted against $\unicode{x1D6E5} =\ell -1$ for a wider suite of computations. The (red) dashed lines display the predictions of lubrication theory (§ C.2); the dot-dashed line in (d) shows $\mathcal{E}/\textit{Bi}$ for Randolph & Houlsby’s solution (labelled RH). In (e), two flow patterns transitions are indicated: the values of $\unicode{x1D6E5}$ at which either the yielded region or the residual plug collides with the outer wall. The vertical dashed lines indicate the prediction for the former using Randolph & Houlsby’s solution. $\textit{Bi}=2048$ and $n=1$ .

Randolph & Houlsby’s slipline construction is analytical, and the net dissipation rate $\mathcal{E}$ can be calculated by hand along with the drag $F_{X}$ and torque $T$ on the cylinder

(5.4) \begin{equation} |{F_{X}}| = 4 \left(\pi +2\sqrt {2} \right) = \mathcal{E}, \quad T \equiv \oint \left [ {\tau _{_{R\Phi }}} \right ]_{r=1} \textrm {d} \phi = 0. \end{equation}

Importantly, Randolph & Houlsby also demonstrate that one can find an acceptable stress field over the attached plugs and surrounding stationary plug; as already noted, the stress field is indeterminate over these regions, but to fully complete a slipline solution, one must show that there is an admissible stress over each plug that satisfies the force-balance equations (3.2) and yield condition ( ${\tau ^2_{_{\textit{XX}}}}+{\tau ^2_{_{YY}}}\lt \textit{Bi}^2$ ).

The other two solutions displayed in figure 10 indicate how the flow, and Randolph & Houlsby’s slipline field, becomes distorted when the outer boundary lies too close to the inner cylinder ( $\ell \lt 2+(1/4)\pi \approx 2.79$ for $\textit{Bi}\to \infty$ ). Other features remain, such as the residual plugs and the viscoplastic shear and boundary layers. Note that the slipline patterns reconstructed from the numerical solutions here are rather imprecise (especially in figure 10 b), which illustrates that construction of these characteristics can be delicate, particularly near non-plastic flow features (e.g. shear layers and rigid plugs).

The progression of the flow pattern as once proceeds from the wide-gap (Randolph–Houlsby) solution to the thin-gap limit is shown more quantitatively in figure 10(d,e). Here, net dissipation rates and yielded areas are presented for a wider suite of computations with varying $\ell$ . There are two second-order transitions in flow pattern present in this progression. The first transition arises when the yielded regions meet the shrinking outer cylinder, which occurs for $\ell =2+ {\pi }/{4}$ according to Randolph & Houlsby’s slipline solution. The second arises only when $\textit{Bi}$ is finite, and corresponds to the collision of the outer wall with the residual plugs lying near the top and bottom of the disk. Both transitions are second order because slipline patterns exist to either side of the transitions that continuously connect to one another; no sudden jumps arise in pattern geometry, and no plugs break.

5.2. The narrow-gap limit; annular squeeze flow

When the gap between the cylinders is narrow, $\ell -1=\unicode{x1D6E5} \ll 1$ , lubrication theory can be used to build solutions analytically. We illustrate this construction for the case of a uniformly translating inner disk ( ${u_{W}}=1$ , ${U_{W}}=0$ ); i.e. an annular squeeze-flow problem. Figure 11 presents numerical solutions in this limit, and a comparison with predictions from lubrication analysis, which is summarised in more detail in § C.2. Figure 10(e,g) also includes the corresponding predictions for the dissipation rate and plug area. The former is quoted in (C36), whereas the latter equals the area of the annulus because the plugs all become small in the thin-gap limit (the relevant Bingham number is $B=\unicode{x1D6E5} ^2\textit{Bi}$ , rather than $\textit{Bi}$ , implying flow becomes effectively Newtonian for $\unicode{x1D6E5} \to 0$ at fixed $\textit{Bi}$ ).

The density plots of the strain-rate invariant in panel (a) illustrate how the lubrication flow adopts a distinctive character in which a relatively weakly sheared region occupies the centre of the gap. This plug-like central flow is the thin-gap relative of the nearly plastic flow region in wider annuli (cf. figures 1 and 10), and is buffered from the walls by layers of stronger shear (which are the cousins of the boundary layers in figure 1). This layered flow pattern is the hallmark of the lubrication limit (§ C.2). Because the central region is not fully rigid, common practice is to refer it as a ‘pseudo-plug’, following Walton & Bittleston (Reference Walton and Bittleston1991). Earlier work has sometimes incorrectly identified the pseudo-plug as a true plug, and then falsely claimed an inconsistency in the lubrication analysis because of the weaker shear that remains present (the fallacy of the lubrication paradox).

The solution in figure 11, with $\unicode{x1D6E5} =1/5$ , is not particularly close to the narrow-gap limit. Consequently, there are noticeable discrepancies with the predictions of lubrication theory in the velocity profiles displayed in figure 11(b,c). The numerical solution matches the asymptotic predictions more closely for narrower gaps, as illustrated in figure 12 which presents numerical solutions for $\unicode{x1D6E5} ={1}/{20}$ . In this example, the main disagreements arise around the border between the pseudo-plug and the more strongly sheared layers. In fact, the lubrication velocity profile here can be further improved by smoothing out the splice between the two regions over a narrow intervening transition zone (Walton & Bittleston Reference Walton and Bittleston1991; Putz et al. Reference Putz, Burghelea, Frigaard and Martinez2008; Muravleva Reference Muravleva2015).

Figure 11. Squeeze flow in the narrow gap (of radial width $\unicode{x1D6E5}$ ) between two cylinders, with the inner cylinder moving to the right at unit speed and the outer cylinder stationary. (a–c) Solutions for $\unicode{x1D6E5} = 1/5$ : (a) density plots of $\log _{10}{\dot \gamma }$ for the values of $\textit{Bi}$ indicated in each quadrant; (b) cuts through the angular velocity along the midline $r=1+(1/2)\unicode{x1D6E5}$ (blue) are compared with the pseudo-plug velocity (C34) predicted by lubrication theory (red dashed); (c) cuts from the rotation rate at $\phi =(1/4)\pi$ (blue) are compared against the asymptotic velocity profile (C30) (red dotted). The dashed lines in (a) and the dots in (c) indicate the predicted fake yield surfaces $r=1+Y$ and $r=1+\unicode{x1D6E5} -Y$ from (C33).

Figure 12. The thin-gap squeeze-flow solution for $\unicode{x1D6E5} = {1}/{20}$ ( $n=1$ ). Density plots of $\log _{10}{\dot \gamma }$ are plotted in panels (a–c) for $\textit{Bi}=2048$ : (a) shows the entire solution over $0\lt \phi \lt {1}/{2}\pi$ ; (b) and (c) show magnifications near $\phi =0$ and ${1}/{2}\pi$ , respectively. The black dashed lines show the fake yield surfaces, $r=1+Y$ and $r=1+\unicode{x1D6E5} -Y$ , from (C33). The blue and white dashed lines show Prandtl’s cycloids (3.19), intersecting either $\theta =0$ or the predicted edge of true plug at $\phi ={1}/{2}\pi$ . Panels (d–h) add further solutions with $\textit{Bi}=128$ , $512$ , $2048$ and $8096$ . The plug borders for solutions are plotted (in red) in (d) and (f). In (e), the yield surfaces near $\phi =0$ are scaled, using $[r(\phi )-1]/[r(0)-1]$ (for the lower) or $[1+\unicode{x1D6E5} -r(\phi )]/[1+\unicode{x1D6E5} -r(0)]$ (for the upper), then replotted against $\phi /[r(0)-1]$ or $\phi /[1+\unicode{x1D6E5} -r(0)]$ . The blue dots show Prandtl’s cycloid. In (e), the prediction (C74) of Appendix C.3 for the edges of the plugs at $\phi = {1}/{2}\pi$ are indicated by (blue) stars, and the dashed lines indicate the radial borders given by (C33). In (g) and (h), cuts of the angular velocity along the midline ( $r=1+(1/2) \unicode{x1D6E5}$ ) and of the rotation rate at $\phi =(1/4)\pi$ (blue lines) are compared with the predictions in (C30) and (C34) (red dotted lines). The (red) dots in (h) show the predicted fake yield surfaces from (C33).

More awkwardly, lubrication theory also predicts that the pseudo-plug (the weakly sheared central region) persists all the way around the annulus. The numerical solutions, however, display plugs attached to the walls near $\phi =0$ and $\pi$ , and true plugs embedded within the pseudo-plug near $\phi =\pm (1/2)\pi$ . These features are not captured by leading-order lubrication analysis (§ C.2), primarily because their angular scale is too small. Despite this, their presence can be detected in the details of that analysis.

First, the attached plugs near $\phi =0$ and $\pi$ can be anticipated because the lubrication analysis predicts that the pseudo-plug expands to fill the entire gap if $\sin \phi =0$ (see (C33) and figures 11 a and 12 a). But this further implies that the entire gap then truly plugs up. Evidently, away from the asymptotic limit, such unyielded points broaden into the attached plugs over an angular scale comparable to the gap.

Second, the lubrication analysis fails to recover the true plugs at $\phi =\pm (1/2)\pi$ because it implicitly assumes that the pseudo-plug solution is always weakly sheared. But, at $\phi =\pm (1/2)\pi$ , the angular velocity of the pseudo-plug reaches a maximum (see figures 11 b and 12 g) and the shear in the pseudo-plug locally vanishes. This opens up the possibility of a different asymptotic solution in which a true plug arises at the centre of the gap. One can construct this solution by adapting a method presented by Walton & Bittleston (who considered axial viscoplastic flow down an eccentric annular duct; see also Frigaard & Ryan (Reference Frigaard and Ryan2004); Putz et al. (Reference Putz, Burghelea, Frigaard and Martinez2008); Liu et al. (Reference Liu, Balmforth and Hormozi2019); Balmforth et al. (Reference Balmforth, Craster and Hewitt2021)). This higher-order extension of the lubrication theory is provided in Appendix C.3, and predicts embedded plugs with angular borders that are indicated by stars in figure 12(f).

The higher-order lubrication analysis assumes that the angular length of the plug is relatively small, but not as small as the gap. Consequently, the plug is predicted to end at a single angular position. The numerical solutions, however, display yield surfaces with a definite shape, albeit of the scale of the gap (see figures 11 and 12), similar to the attached plugs. In fact, both yield surfaces closely follow Prandtl cycloidal sliplines, as illustrated in figure 12(b,c).

We rationalise this observation by first noting that on the scale of the narrow gap, the walls appear straight and parallel with a nearly constant pressure gradient, $\partial P / \partial s=-\Gamma \textit{Bi}$ . At $\phi =0$ , we have $Y\to 0$ or $\Gamma =2 \Delta ^{-1}$ , whereas $\Gamma = [(1/2)\unicode{x1D6E5} -Y({1}/{2}\pi )]^{-1}$ near $\phi = {1}/{2}\pi$ . Inserting these two options into (3.19) gives the sliplines also plotted, choosing either $x_0=0$ or the location of the edge of the plug near $\phi =(1/2)\pi$ predicted by (C74).

To reinforce this point, we exploit the scale invariant form of the cycloid noted in § 3.4.2 to collapse all the yield surfaces measured from the numerical solutions near $\phi =0$ : for the lower yield surfaces, we scale $r(\phi )-1$ and $\phi$ by $r(0)-1$ , then scale $1+\unicode{x1D6E5} -r(\phi )$ and $\phi$ by $1+\unicode{x1D6E5} -r(0)$ for the upper yield surface. As shown by figure 12(e), this scaling successfully collapses the observed yield surfaces onto Prandtl’s cycloid.

Note that the floating plugs at $\phi =\pm (1/2)\pi$ correspond to the small-gap limit of the residual plugs at the core of the semi-circular fans of Randolph & Houlsby’s slipline field (figures 1 a,b and 10 a). Indeed, the former also disappear as $\textit{Bi}\to \infty$ ( $Y\to {1}/{2}$ in this limit; see (C74)). On the other hand, the plugs attached to the inner cylinder at $\phi =0$ and $\phi =\pi$ are the remnants of the triangular plugs attached to the front and back of the disk in figure 1 (and are no longer serendipitous). The ambient plug that surrounds and cloaks the outer cylinder in figure 1 becomes compressed into the plugs attached to the outer cylinder at $\phi =0$ and $\phi =\pi$ .

5.3. Translation with rotation

Figure 13. Translating and rotating cylinders for varying rotation speed $\Omega$ (as indicated; translation speed in positive $x$ direction is again scaled to unity). On the left of each plot, we display a density map of $\log _{10}{\dot \gamma }$ for the numerical solution along with sample streamlines. On the right, we show the corresponding slipline solution, with various important points indicated. Plugs are shaded black, and the $\alpha$ ( $\beta$ ) lines are plotted as red (blue) lines.

Returning again to the limit in which the outer cylinder is sufficiently distant that it does not affect the flow, we now allow the inner cylinder to rotate as well as translate (Hewitt & Balmforth Reference Hewitt and Balmforth2018). Solutions in the plastic limit are illustrated in figure 13. Surprisingly, although one expects that rotation breaks sideways symmetry (up-down in the figures), the stress field remains symmetrical and equivalent to that of Randolph & Houlsby as long as $\Omega \lt ({1}/{\sqrt 2})$ (figure 13 a). This feature becomes reflected in the drag and torque, which remain at the values given in (5.4) to leading order; see figure 14. Nevertheless, the velocity field is asymmetrical, as illustrated in figure 13(a). This solution demonstrates how it is possible for an asymmetrical velocity field to be consistent with a symmetrical stress field in the plastic limit, since only the orientation of the strain-rate tensor needs to be symmetrical in (3.14), not the invariant $\dot \gamma$ itself.

Figure 14. (a) Forces $|{F_{X}}|$ and dissipation rates $\mathcal{E}$ , scaled by $\textit{Bi}$ , plotted against $\Omega$ for translating and rotating cylinders. In (b,c) we show the corresponding torques $T/\textit{Bi}$ and yielded areas, respectively. The stars indicate data taken from numerical solutions. The black, red and blue lines show the slipline predictions for the three different flow patterns, which are illustrated in the insets to (b) (density maps of $\log _{10}{\dot \gamma }$ ) and shown in more detail in figure 13 (with the same colour scale). In (a), the red stars show $\mathcal{E}$ , as computed from the velocity field, and the red circles indicate $|{F_{X}}+\Omega T|$ .

When $\Omega$ reaches $({1}/{\sqrt 2})$ , the asymmetrical velocity field is no longer able to remain consistent with the symmetrical stress field. This arises because the (anti-clockwise) tangential velocity along the section $O_+ D_+$ of the disk is $\Omega -\sin \phi$ , $(1/4)\pi \lt \phi \lt (1/2)\pi$ , whereas the plastic region on the other side of the wall boundary layer has zero angular velocity. Thus, between the plastic region and disk there is a tangential slip velocity of $\sin \phi -\Omega$ , which must have the same sign as the local shear stress, ${\tau _{_{R\Phi }}}=-\textit{Bi}$ . This can only be true if $\Omega \lt ({1}/{\sqrt {2}})$ .

For higher rotation rates, the stress field necessarily deviates from the Randolph & Houlbsy form and becomes asymmetrical, as seen in figure 13(b). In particular, the rightmost involutes leaving the front of the cylinder no longer meet at $(x,y)=(\sqrt 2,0)$ . Instead, a gap opens, with the attached plug separated from the surrounding ambient plug by the circular slip surface $C_-C_+$ . This arc is centred at point $E$ , which is the centre of rotation of the combined velocity field of the cylinder, $(1,0)+\Omega (-y,x) = \Omega (y_{\Omega }-y,x)$ with $y_{\Omega }=\Omega ^{-1}$ .

The arc $C_-C_+$ widens as $\Omega$ is increased still further. Eventually, for $\Omega =1$ , the rotation centre $E$ reaches the centre of the fan $O_+$ and the attached plug meets the semi-circular fan on that side of the disk. At this stage, the entire upper fan plugs up, as illustrated in figure 13(c). This third slipline pattern was constructed previously by Hewitt & Balmforth (Reference Hewitt and Balmforth2018). As highlighted by the area of the yielded region plotted in figure 14(c), the first switch in flow pattern at $\Omega =1/\sqrt 2$ is second order because the intersecting slipline fields are compatible. The second switch at $\Omega =1$ is first order, because the rotating plug of the last pattern breaks.

The drag, torque and dissipation rate can be calculated analytically for all of the three slipline patterns in figure 13 (see § B.3). In each case, $\mathcal{E} = |{F_{X}}|+|\Omega T|$ . The predictions are included in figure 14. Note that the torque calculated numerically for solutions with $\textit{Bi}=2048$ does not vanish for $\Omega \lt ({1}/{\sqrt 2})$ because of a residual viscous torque stemming from the asymmetrical viscoplastic boundary and shear layers. For $\textit{Bi}\to \infty$ , however, the first slipline pattern (rather surprisingly) implies that a translating disk can rotate at a rate $0\lt \Omega \lt ({1}/{\sqrt 2})$ without requiring any torque or incurring any additional dissipation beyond that expended for the non-rotating disk.

Strictly speaking, as noted earlier, to establish that the slipline patterns in figure 13 are the exact solutions in the plastic limit one must further demonstrate that there is an acceptable stress field over each of the plugs. Randolph & Houlsby perform this construction for $\Omega =0$ , but none has been provided for the asymmetrical or rotating plugs of figure 13(b,c). Following Hewitt & Balmforth (Reference Hewitt and Balmforth2018) and § B.3, one can establish there are matching velocity and stress fields elsewhere which satisfy force balance, the constitutive law and boundary conditions. Hence, except for the extension of the stress field into the plugs, the extremum principles imply that $\mathcal{E} = |{F_{X}}|+|\Omega T|$ and the solutions are exact. That said, the numerical solutions provide evidence that admissible stress fields can be found for the plugs. Therefore, the slipline patterns do not only provide bounds on the net drag or dissipation rate, as sometimes they are constructed to do, but are the true solutions. Nor is it strictly necessarily to compute both $\mathcal{E}$ and $|{F_{X}}|+|\Omega T|$ independently and verify they are equal (they must be if all is consistent).

5.4. Squirmers

The two examples shown in figure 10(a) and 13(a) repeat the possibility for two solutions to possess the same stress field but different velocity fields in the plastic limit. A third example is shown in figure 15(a,c), which corresponds to an idealised model of a circular microswimmer that ‘squirms’ through fluid by activating a tangential surface motion (Lighthill Reference Lighthill1975). More specifically, we consider domains for which the outer boundary $r=\ell$ is sufficiently far away that it is hidden by the surrounding plug, and, in (5.1), we take ${U_{W}}(\phi )=\phi (\pi -|\phi |)$ for $-\pi \lt \phi \lt \pi$ (sinusoidal surface velocities were considered by Supekar et al. (Reference Supekar, Hewitt and Balmforth2020)). The (dimensionless) translation speed $u_{W}$ is now a parameter of the problem, although for a genuine self-propelled squirmer, $u_{W}$ must be selected by demanding that no net forces are exerted on the inner disk (see below).

Figure 15. Model squirmers with a surface angular velocity ${U_{W}}(\phi )=\phi (\pi -|\phi |)$ for $-\pi \lt \phi \lt \pi$ , and $\textit{Bi}=2048$ ( $n=1$ ). In (a), the (logarithmic) strain rate is plotted for a Randolph & Houlsby-like solution with a translation speed of ${u_{W}}=0.4\textit{Bi}^{-1/2}$ ; a boundary-layer solution with ${u_{W}}=5\textit{Bi}^{-1/2}$ is presented in (b). A magnified view of the first solution in the first quadrant is shown in (c), along with sample streamlines (green). The boundary layer against the squirmer is shown in more detail in (d). A corresponding plot for another boundary-layered solution with ${u_{W}}=0.7\textit{Bi}^{-1/2}$ is presented in (e) (the colour scale is the same as for (d)); the dotted line shows the prediction (C23) for the boundary-layer width. A suite of computations with varying $\textit{Bi}^{1/2}{u_{W}}$ are presented in (f,g). Plotted are (f) the drag force $F_{X}$ and net dissipation rate $\mathcal{E}$ , and (g) the yielded area. In (d), the (red and blue) solid lines show the predictions in (5.5) and (C24); the (blue) dashed line shows (5.4). The self-propelled squirmer, with zero net drag, is indicated by the star.

All three examples in figures 10(a), 13(a) and 15(a,c) have the same (leading-order) slipline pattern, and all experience the same drag force $|{F_{X}}|\approx 4(\pi +2\sqrt 2)\textit{Bi}$ (figure 15 f). One minor difference is that the triangular sections at the front and back of the disk that correspond to checkerboards in the slipline field are no longer fully plugged up for the squirmer. The viscoplastic boundary layer at $r=1$ also persists all the way around the disk. These discrepancies arise because the squirmer’s surface velocity is not consistent with solid-body motion over the checkerboard, unlike in figures 10(a) and 13(a). Instead, the checkerboards contain true plastic deformation, or at least partly so, with the fine structure of the flow pattern tracking the $45^\circ$ sliplines (see figure 15 c). Besides this, the circular fan turns out to mostly have a rigidly rotating velocity field, leading to a persistent serendipitous plug, instead of the residual plug of figure 1.

Much as with the rotating, translating cylinder considered above, the squirmer’s stress solution takes the Randolph–Houlsby form only for as long as the translation speed $u_{W}$ is not too large. A different solution emerges when $u_{W}$ exceeds a critical value of $O(\textit{Bi}^{-1/2})$ , which is much less than the the $O(1)$ tangential surface motion ${U_{W}}(\phi )$ . Because this transition value is so small, the dissipation in the Randolph–Houlsby-type solutions is dominated by the boundary layer against the cylinder (see figure 15 a,c,d). The power balance (5.2) then reduces to

(5.5) \begin{equation} \mathcal{E} \sim \textit{Bi} \oint | {U_{W}} | \; \textrm {d} \phi, \end{equation}

since ${\tau _{_{R\Phi }}}\sim -\textit{Bi}$ sgn $({U_{W}})$ over the boundary layer (§ C.1.2).

Beyond the critical transition in the translation speed, the flow becomes completely confined to the viscoplastic boundary layer around the inner cylinder and there is no extended region of plastic deformation, as illustrated in figure 15(b,e). In this case, the boundary-layer theory of § C.1.2 provides analytical predictions for the boundary-layer width and drag force $F_{X}$ , with the net dissipation rate again given by (5.5). In figure 15(e), the prediction for the boundary-layer thickness (equation (C23)) agrees satisfyingly with the numerical solution, except near $\phi =0$ and $\pi$ , where the boundary-layer approximation locally breaks down.

The net drag and dissipation rate for a suite of computations with varying $u_{W}$ are plotted in figure 15(f) and compared with the predictions stated above. Also displayed is the area of the yielded region, which abruptly jumps to small values when flow localises to the squirmer surface. That is, the transition in flow pattern is first order, precipitated by the breakage of the plug around the localised boundary-layer flow. The transition also arises when the drag for the boundary-layer solution matches that for Randolph & Houlsby’s construction (see figure 15 f), which implies

(5.6) \begin{equation} {u_{W}} = \frac {\frac {1}{3}\textit{Bi}^{-1/2}}{\sqrt {\pi +2\sqrt {2} - 1}} {\left [ \int _0^\pi \frac {|{U_{W}}|^3\; \textrm {d} \phi }{\sin \phi } \right ]}^{\frac {1}{2}}, \end{equation}

from boundary-layer theory (§ C.1.2; Supekar et al. (Reference Supekar, Hewitt and Balmforth2020)). This criterion is consistent with a minimisation of the net dissipation rate (§ 3.3) because, hidden beyond the leading-order term in (5.5), lies a contribution of ${u_{W}} |{F_{X}}|$ from translation.

Note that the force $F_{X}$ must vanish for a self-propelled squirmer. As indicated in figure 15(f), this condition selects the particular locomotion speed ${u_{W}} \approx 1.54\textit{Bi}^{-1/2}$ , shown by the filled star. The analytical prediction from boundary-layer analysis

(5.7) \begin{equation} {u_{W}} \sim \frac {1}{3} \textit{Bi}^{-1/2} \left [ \int _0^\pi \frac {|{U_{W}}|^3\; \textrm {d} \phi }{\sin \phi } \right ]^{\frac 12}, \end{equation}

(§ C.1.2; Supekar et al. (Reference Supekar, Hewitt and Balmforth2020)) is slightly smaller, as seen in figure 15(f).

Figure 16. Model squirmers for higher translation speeds $u_{W}$ ; $(n,\textit{Bi})=(1,2048)$ . On the left, strain-rate maps and sample streamlines are plotted for (a) ${u_{W}}=45\textit{Bi}^{-1/2}\approx 1$ , (b) ${u_{W}}=80\textit{Bi}^{-1/2}\approx 1.77$ and (c) ${u_{W}}=250\textit{Bi}^{- 1/2}\approx 5.52$ . The lower half of the plots show slipline fields, constructed using the argument described in the main text. The thicker (red) lines in (a,b) reconstruct the curve inside the disk that appears to generate the involutes. A wider suite of computations with varying $u_{W}$ is presented on the right, plotting (d) the drag force $F_{X}$ and net dissipation rate $\mathcal{E}$ , and (e) the yielded area. In (d), the (red and blue) dashed lines show (5.4) and $\mathcal{E}=|{F_{X}}|/u_{W}$ . The Randolph–Houlsby (RH) slipline pattern features to the right of the vertical line indicated.

Though not relevant to locomotion, one can consider much larger translation speeds of $O(1)$ , as illustrated in figure 16. Because the boundary-layer width grows with $u_{W}$ , the flow delocalises again in this parameter range, forming a wider region of plastic deformation (panels (a,b)). Eventually, once this region extends sufficiently far from the disk, the pattern transforms into another with the Randolph & Houlsby form (panel (c)). The drag is then approximately given by (5.4), with a dissipation rate of $\mathcal{E}\approx {u_{W}} {F_{X}}$ for large $u_{W}$ .

The nearly plastic solution for ${u_{W}}\lesssim 2.45$ has an interesting structure, as highlighted by the slipline patterns shown in the lower halves of figure 16(a,b). Here, rather than attempt to reconstruct the slipline using the stress field (as in figure 10), we make use of the observation that there are no plugs within the yielded region and the $y$ -axis must correspond to a slipline as it is a line of symmetry. Hencky’s rules then imply all the sliplines of the same family as the $y$ -axis are straight, and Geiringer’s equations demand that the other family correspond to streamlines. We therefore reconstruct the sliplines by finding contours of constant velocity direction and streamlines. This leads to the patterns shown in figure 16(a,b,c), with the latter largely reproducing the Randolph & Houlsby slipline field.

Note that the reconstruction fails once the sliplines cross a shear layer, which happens near the front and back of the disk, where the slipline pattern from the $y$ -axis meets the local checkerboards there. For ${u_{W}}\lesssim 2.45$ , the slipline field is composed of involutes generated by a curve lying inside the disk (see figure 16 a,b). Importantly, there is no boundary layer along the surface of the disk that dictates the slipline angle there (unlike for the Randolph–Houlsby patterns). The geometry of the curve from which the involutes are built must therefore be controlled by the strain-rates at the disk’s surface.

Finally, figure 17 presents solutions for a squirmer with a different surface velocity distribution. In this case, the angular surface velocity is more step-like, with $U_{W}(\phi )=\tanh (5\sin \phi )$ . Unlike in previous examples, the nearly plastic solution shown in figure 17(a) with $u_{w}=0$ escapes the Randolph & Houlsby form, because of the appearance of a stress discontinuity (indicated by the green arrow; the yield surfaces are a little difficult to determine in this example owing to the presence of a wide region held very close to the yield stress). The discontinuity also appears to persist for the solution with much lower $\textit{Bi}$ shown in figure 17(b), breaking into two pieces that each track curved yield lines (black arrows). The nearly vertical one close to $y=\phi =0$ manifests mainly in the component ${\tau _{_{RR}}}=-{\tau _{_{\Phi \Phi }}}$ , in accord with the jump conditions stated in § 3.4.1. Thus, this squirmer offers another example in which stress discontinuities appear, even at finite $\textit{Bi}$ .

Figure 17. Sqirmer solutions for surface velocity, $u_{w}=0$ and $U_{W}(\phi )=\tanh (5\sin \phi )$ , and (a) $\textit{Bi}=2048$ and (b) $\textit{Bi}=1$ . Shown are sample streamlines (green) superposed on density maps of $\log _{10}{\dot \gamma }$ ; in (a), the diagnosed slipline field is indicated by the black lines in $y\lt 0$ . Scaled stress components for the solution with $\textit{Bi}=1$ are shown in (c,d). The plugs (diagnosed by the contour ${\dot \gamma }=10^{-4}$ ) are shaded dark blue with dashed red borders (the yield surfaces). The arrows point to what appear to be stress discontinuities.

5.5. Translation with axial motion

Moving beyond two-dimensional flow, the relatively simple geometry of the annular gap (figure 9) offers a convenient setting in which to explore the impact of adding flow in the third spatial dimension on the dynamics in the nearly plastic limit. In particular, we can consider three-dimensional flow fields $(u,v,w)$ that depend only on the planar coordinates $(x,y)$ . This task was accomplished by Hewitt & Balmforth (Reference Hewitt and Balmforth2018) and Hewitt & Balmforth (Reference Hewitt and Balmforth2022), who considered the flow around a cylinder translating at different angles to its axis, taking the annular gap to be sufficiently wide that the outer boundary is cloaked by a plug (rendering the problem equivalent to that of an infinitely long cylinder in an infinite domain). Further results for this problem are presented in figure 18. When the angle of motion is not closely aligned with the cylinder axis, the in-plane flow field is similar to the Randolph–Houlsby pattern; the boundary layers, shear layers and extended regions of nearly plastic deformation all survive the addition of an order-one out-of-plane velocity component (at large Bingham number). The nearly plastic regions, however, no longer possess an embedded slipline structure. Nevertheless, the structure of the two-dimensional flow pattern broadly carries over to a genuinely three-dimensional flow.

Figure 18. Flow around infinitely long cylinders translating with respect to their axes at angles $\unicode{x1D6E5}$ of (a) $ {7\pi }/{22}$ , (b) $ \pi/ 8$ and (c) $ {\pi }/{200}$ . Shown are density plots of $\log _{10}{\dot \gamma }$ with sample streamlines for the in-plane velocity field (green, $y\gt 0$ ) and contours of constant $w$ (black, $y\lt 0$ ). $\textit{Bi}=2048$ . (d) A magnification of the boundary layer for the solution also shown in (c); panel (e) shows a similar plot for a fourth solution with localised boundary-layer flow at $\unicode{x1D6E5} =5\pi \times 10^{-5}$ . The force components, $F_{X}$ and $F_{Z}$ (solid red and blue), and yielded area (dotted) are shown in (f) as functions of $\unicode{x1D6E5}$ , for a wider set of computations. The filled stars show the angles for the solutions in (a,b,c,e), whereas the open (yellow) stars indicate the limits $({F_{X}},{F_{Z}})=(8\sqrt 2+4\pi, 0)\textit{Bi}$ and $(0,2\pi )\textit{Bi}$ , expected for $\unicode{x1D6E5} =(1/2)\pi$ and $\unicode{x1D6E5} =0$ , respectively. Panel (g) replots the data for the small window of angles (of order $\textit{Bi}^{-1}$ ) over which the drag force reorientates. The dashed lines in (e) and (g) shows the predictions from boundary-layer analysis (Hewitt & Balmforth Reference Hewitt and Balmforth2022) of ${F_{X}}\sim 9\pi \textit{Bi}^2\unicode{x1D6E5}$ and boundary-layer thickness ( $\sqrt {2/\textit{Bi}}$ ).

As the direction of motion becomes closer to axial, the flow pattern loses its similarity with the Randolph–Houlsby slipline field (see figure 18 b). More interestingly, when translation is almost completely aligned to the axial direction, fluid flow becomes primarily axial, and mostly confined to a narrow boundary layer around the cylinder (panel (c)). However, sideways motion still persists, generating a significant transverse force $F_{X}$ . Indeed, it is only when the direction of motion reaches a relatively narrow window of angles close to the axial direction (with a width of $O(\textit{Bi}^{-2/(n+1)})$ ) that the sideways force abruptly drops from $O(\textit{Bi})$ values to zero (see figure 18 f,g). Hewitt & Balmforth (2018, 2022) focused on the impacts of this feature in swimming problems in which relatively long waves along a cylindrical worm or tail drive locomotion through a viscoplastic fluid. In particular, it was shown that because of the abrupt reorientation of the force, locomotion took the form of ‘burrowing’, with each cylindrical segment of the worm or tail moving almost axially (Hewitt Reference Hewitt2024).

The sharp switch in the drag force arises because, in the perfectly plastic problem, axial motion can be accounted for purely by slip over the cylinder surface (which broadens into a viscoplastic boundary layer for finite $\textit{Bi}$ ). But sideways translation, no matter how small, cannot. Thus, an arbitrarily small, but finite sideways translation must always create a region of perfectly plastic deformation away from the cylinder, generating a finite sideways drag. In the viscoplastic problem, this sideways drag can be switched off when the boundary layer against the cylinder consumes all motion, leading to the $O(\textit{Bi}^{-2/(n+1)})$ window of angles (a true boundary-layered solution from within this window is shown in figure 18 e). The switch from a localised boundary layer to extended sideways flow takes the form of a first-order transition in flow pattern in which part of the plug surrounding the cylinder breaks. The transition is visible in figure 18(g) as the kink in $F_{X}$ and sharp rise in yielded area near $\unicode{x1D6E5} =2.5\times 10^{-4}$ .

The perfectly plastic solution that emerges in the limit $\textit{Bi}\to \infty$ for translation angles $1\gg \unicode{x1D6E5} \gg O(\textit{Bi}^{-2/(n+1)})$ is similar to that for the translation of a two-dimensional disk with a free-slip condition on its surface (Martin & Randolph Reference Martin and Randolph2006; Supekar et al. Reference Supekar, Hewitt and Balmforth2020). This arises because most of the axial motion is confined to the viscous boundary layer; the in-plane velocities are much smaller, and not so confined. Viscoplastic boundary layer theory, on the other hand, demands that the stress components, $\tau _{_{R\Phi }}$ and $\tau _{_{RZ}}$ , dominate close to the cylinder. But this further implies that $\tau _{_{RZ}}\gg {\tau _{_{R\Phi }}}$ since $w\gg (u,v)$ . The boundary layer therefore enforces a stress-free condition on the plastic flow further from the cylinder. Nevertheless, the axial velocity does not vanish outside the boundary, and merely becomes as small as $u$ and $v$ . The problem does not therefore correspond precisely to flow around a two-dimensional, stress-free disk, although the flow pattern looks similar (see Supekar et al. (Reference Supekar, Hewitt and Balmforth2020)).

6.1. Failure

Failure problems like the initial collapse of the block in figure 19 are commonplace in geotechnical engineering, in addressing the stability of slopes, embankments and other structures (Terzaghi Reference Terzaghi1943). For cohesive soil, the problem is traditionally dealt with within the framework of ideal plasticity theory, which corresponds to the plastic limit for failure of a block of yield-stress fluid. The problem is also more commonly couched in terms of the limit analysis of plasticity (Prager & Hodge Reference Prager and Hodge1951), wherein one exploits the corresponding extremum principles to § 3.3 to provide bounds on the critical value of $\textit{Bi}$ at which failure occurs.

For a rectangular block, the critical value for failure, $\textit{Bi}_{crit}$ , depends on the aspect ratio $\ell$ (Liu et al. Reference Liu, Balmforth, Hormozi and Hewitt2016), as summarised in figure 20(a). In addition, at $\textit{Bi}=\textit{Bi}_{crit}$ deformation adopts a specific spatial pattern, or ‘failure mode’. Such patterns correspond to the dynamical modes for incipient collapse in elastic systems, or the neutral modes at the onset of linear instability in a viscous flow, with the critical value $\textit{Bi}_{crit}(\ell )$ playing the role of the stability threshold. However, the yield stress ensures that viscoplastic failure is equivalent to a nonlinear eigenproblem (with an indeterminate equilibrium stress), unlike a conventional linear stability analysis. Demanding symmetry along $x=0$ also rules out failure modes that correspond to sideways buckling, divorcing the current formulation from the classical problem of Euler’s elastic column.

Figure 20. Failure modes for the problem sketched in figure 19. (a) The critical Bingham number $\textit{Bi}_{crit}$ above which there is no motion as a function of the aspect ratio $\ell = \mathcal{L}/\mathcal{H}$ , and (b) the heights, $Z_0$ (filled circles) and $Z_\ell$ (open squares), at which failure occurs along $x=0$ and $x=\ell$ , respectively. Four distinct modes are identified in (a,b) and plotted with different colours; samples of each are illustrated in the remaining panels (e–f), which present density plots of $\log _{10}(\dot \gamma )$ along with selected streamlines. The solution at $\ell =1$ , shown further in the inset to panel (a), appears to have mixed character. In (c), a corresponding slipline pattern is also shown (in $x\lt 0$ ); the predictions of slipline theory are shown by (red) solid lines in (a,b). The dashed lines in (a,b) indicates the aspect-ratio-independent theoretical limit for $\ell \gt 1$ (Lyamin & Sloan Reference Lyamin and Sloan2002a ,Reference Lyamin and Sloan b ; Martin Reference Martin2011). The inset to (f) shows a reconstruction of the slipline field. The failure modes here are identified by computing solutions at sequentially larger $\textit{Bi}$ until the maximum strain rate over the block, ${\dot \gamma }_{max}$ , falls below $2\times 10^{-3}$ ; the critical Bingham number is then determined by extrapolating ${\dot \gamma }_{max}$ to zero.

For relatively slender blocks, the column of material fails primarily at its base, below a vertically falling plug. The distinctive flow pattern, or the failure mode, that results is illustrated in figure 20(c), and can be built by slipline theory (Chamberlain et al. Reference Chamberlain, Sader, Landman and White2001). The slipline pattern, also shown in the figure, contains a checkerboard stemming from the free surface at the side, together with a fan emanating from the lower outer corner. The sliplines in both are no longer straight, but curve because of gravity and the impact of hydrostatic pressure on the slipline invariants (see Appendix B.4). Note that a section of the expected plastic region appears to plug up in the computation, a feature permitted by the shear layer along its top border. This leads to a serendipitous plug at the side that rotates out of position.

The slipline solution illustrated in figure 20(c) provides a prediction for $\textit{Bi}_{crit}(\ell )$ that is added to figure 20(a). Other properties of the failure mode can also be established, such as the heights, $Z_0$ and $Z_\ell$ , up to which failure occurs at $x=0$ and $x=\ell$ , respectively (cf. figure 20 c). These alternative measures offer clearer diagnostics of the transitions in flow patterns (or failure mode) that take place as the aspect ratio is varied; see figure 20(b). The first transition arises for $\ell \approx 0.5$ and is first order, involving the breakage of the upper plug along the symmetry line. The breakage eliminates purely vertical collapse (and the slipline field of panel (c)), and instead permits the upper corners of the initial block to rotate out sideways. Widespread plastic deformation takes places underneath the rigid rotation.

The failure mode changes again when the aspect ratio $\ell$ reaches values closer to unity, with a prominent shear layer emerging that pierces the block close to its core. Figure 20 identifies what appear to be four distinct failure modes altogether. In panels (a,b), the four modes are plotted with different colours, and examples of each mode are shown in (c,d,e,f). In these failure computations (which proceed in an iterative manner in order to determine $\textit{Bi}_{crit}$ ), the modes and the transitions between them are less straightforward to characterise because strain rates are necessarily small, which can obscure the geometry of the yield surfaces, plugs and plastic regions.

Once $\ell$ exceeds unity, collapse only arises in the vicinity of the vertical face (figure 20 f); the central core of the block never fails, and the critical value becomes independent of aspect ratio. The solution then corresponds to the failure mode of an infinitely wide block, bounds for which have been provided by Lyamin & Sloan (Reference Lyamin and Sloan2002a ,Reference Lyamin and Sloan b ). This failure mode has $\textit{Bi}_{crit}\approx 0.2648$ , $Z_0=1$ and $Z_\ell \approx (1/2)$ (as seen in figure 20 a,b), and possesses a slipline solution that has been constructed by Martin (Reference Martin2011). This slipline field is rather intricate, however, and Martin’s strategy to build it is more like our slipline reconstruction from the numerical viscoplastic flow computation (cf. the inset in figure 20 f) than an explicit construction using the slipline equations, as in Appendix B. In fact, beyond the relatively narrow limit, the modes of failure and their corresponding slipline patterns are fairly complicated; for no aspect ratio do they take the form of simple arcs of failure between rigid blocks (as is often assumed in the geotechnical engineering literature in order to establish simple bounds; e.g. Terzaghi (Reference Terzaghi1943); Chamberlain et al. (Reference Chamberlain, Horrobin, Landman and Sader2004)). Nevertheless, the failure modes again form the patchwork of plugs, plastic regions and boundary layers that we now expect in the plastic limit.

6.2. Runout; thin-film spreading

After failure, the block collapses into a free-surface gravity current that runs out to lower the gravitational stresses towards the yield stress everywhere. The time-dependent dynamics of the spreading gravity current can be interrogated by numerical simulation (see Liu et al. (Reference Liu, Balmforth, Hormozi and Hewitt2016) and references therein; Valette et al. (Reference Valette, Pereira, Riber, Sardo, Larcher and Hachem2021)), although computations can be time consuming and tricky. An alternative is to assume that flows are relatively shallow ( $\ell \ll 1$ ) and exploit lubrication theory, as has been done for viscous gravity currents (e.g. Huppert Reference Huppert1982). Balmforth et al. (Reference Balmforth, Craster, Rust and Sassi2006) and Appendix C.2.2 review details of the viscoplastic version of this analysis, which furnishes an evolution equation for the local depth of the fluid layer $h(x,t)$ . Related situations in which the analysis has been exploited include bearings and blade coating (Hewitt & Balmforth Reference Hewitt and Balmforth2012; Lister & Hinton Reference Lister and Hinton2022), the flow of lava and mud (Blake Reference Blake1990; Hinton & Hogg Reference Hinton and Hogg2022; Hinton, Hewitt & Hogg Reference Hinton, Hewitt and Hogg2023), the growth of mountain ranges (Ribinskas et al. Reference Ribinskas, Ball, Penney and Neufeld2024; Taylor-West & Hogg Reference Taylor-West and Hogg2024), the viscoplastic linings of lung airways (Shemilt et al. Reference Shemilt, Horsley, Jensen, Thompson and Whitfield2022, Reference Shemilt, Horsley, Jensen, Thompson and Whitfield2023) and other problems involving surface tension (Ross, Wilson & Duffy Reference Ross, Wilson and Duffy2001; Jalaal & Balmforth Reference Jalaal and Balmforth2016; Jalaal, Stoeber & Balmforth Reference Jalaal, Stoeber and Balmforth2021; Ball & Balmforth Reference Ball and Balmforth2024).

Figure 21. A shallow-layer solution to (C39) for the collapse of a block above a horizontal base plate for $\textit{Bi}=0.004$ and $\ell =10$ . The initial shape is given by $h(x,0)=h_\infty +1+\tanh [100(x-(1/2)\ell )]$ , where $h_\infty =10^{-3}$ denotes a pre-wetted film coating the entire base plate (added to ease numerical computations, along with the regularisation parameter, $Y_\infty =10^{-6}$ , introduced in the replacement $Y=$ Max $(Y_\infty, h-\textit{Bi}/|\partial h / \partial x|)$ in (C39)). The geometry of the slump is sketched in (a). Panel (b) shows snapshots of $h(x,t)$ and $Y(x,t)$ ; the final profile of $h(x,t)$ (as given by (6.3)) is shown by the red dots. In (c), scaled snapshots of $Y$ and $u_p$ for $t\gt 0.1$ are plotted against $\xi =x/X(t)$ (with $X(t)$ defined by $h(X,t)=1.01h_\infty$ ); $Y$ is scaled by its maximum value, $u_p$  by twice the value taken where $\xi =(2/3)$ . Panel (d) shows times series of $h_m(t)=h(0,t)$ and $X(t)$ , with the dashed lines show the predictions from solving (C46).

Figure 21 presents a sample numerical solution to the thin-film evolution equation (C39), in which a block of fluid collapses and runs out over the underlying horizontal plane. Note that, because shear stresses vanish along the free surface (in the absence of surface tension), a superficial pseudo-plug must appear during runout, as sketched in figure 21(a). The relatively strongly sheared layer underneath the pseudo-plug has depth $Y(x,t)\equiv h-\textit{Bi}/|\partial h / \partial x|$ ; as the fluid runs out, and surface slopes and stresses decline, this depth descends to $y=0$ . Importantly, runout becomes limited by the yield stress and fluid eventually reaches a slumped state with a sloped surface, unlike in the corresponding viscous problem (in the absence of surface tension).

As noted by Nye (Reference Nye1951) and others, the final state is given by $Y\to 0$ , or

(6.3) \begin{equation} h = \sqrt { 2\textit{Bi} (X-x)}, \end{equation}

where $x=X$ denotes the fluid edge (see Appendix C.2.2). In the example shown in figure 21, the entire block yields and flows on its path to the final state. When the initial block is sufficiently wide, or the yield stress large enough, however, part of the block remains intact, leaving a flat-topped final deposit (Balmforth et al. Reference Balmforth, Craster, Rust and Sassi2006). Here, we consider only fully collapsed states.

The convergence to the final state can be described analytically and follows a relatively long, algebraic dependence on time (see Matson & Hogg (Reference Matson and Hogg2007); Hogg & Matson (Reference Hogg and Matson2009); Appendix C.2.2). For a Bingham fluid the convergence takes the form $t^{-1}$ , as illustrated in figure 21(d). During this passage, $\partial u_{p}/\partial x\gt 0$ , where $u_p$ is the speed of the pseudo-plug, indicating a state of horizontal expansion (cf. figure 21 c).

Note that dambreaks from a step-like initial condition are problematic in lubrication theory because $|\partial h / \partial x|$ becomes arbitrarily large at the step, but should be asymptotically small for the theory to remain valid (which is partly why the step in the initial condition in figure 21 has been smoothed out slightly). Lubrication theory cannot therefore address the initial failure of the block, as in § 6.1. Indeed, the theory predicts that the block always fails (the yield criterion being $Y=h-\textit{Bi}/|\partial h / \partial x|\gt 0$ , which is inevitably satisfied if $|\partial h / \partial x|\to \infty$ ).

6.3. Slumped states with finite slope

Armed with slipline theory, one can attempt to advance beyond the thin-film limit and build the final shape for slumps with finite slope. In this context, Nye provided an elegant construction for gravity-driven plastic flow with a free surface in the context of modelling the shape of the snout of a glacier. This construction applies equally to viscoplastic flow in the plastic limit.

Details of Nye’s construction are provided in Appendix B.4, which begins by artificially providing data along the outer arc of a circular fan, as illustrated in figure 22(a). The opening angle of the fan and its height $h_0$ act as two parameters. The surprising feature of Nye’s construct is that, as one proceeds to shallower depths, the slipline field converges to a special solution regardless of the precise details of the initial slipline (i.e. the fan angle and $h_0$ ). Indeed, another choice for the initial slipline based on Prandtl’s cycloid (see Appendix B.4) leads to the slipline field shown in figure 22(b). This alternative construction converges more quickly to the special solution.

Figure 22. Surface profiles constructed for (a–d) horizontal extension and (e,f) horizontal compression. Nye’s original construction of the slipline field from part of a circular fan is shown in (a); the construction in (b) uses a modified Prandtl cycloid to launch the sliplines. In (c), a magnification near the nose is shown for the slipline field of (b); the shaded region indicates the plug introduced by Nye to avoid the sliplines unphysically turning downwards from $y=0$ . The constructions in (d,e) uses the minimisation scheme outlined in Appendix B.4 to build the slipline field using a section of the free surface for $h_1 \gt (h/Bi) \gt h_2$ with $h_1=12$ and $h_2=2$ . The red dots show surface profiles constructed for the different values of $h_2$ marked by stars ( $h_2=0.505\pi$ , $(5/3)$ , 2, 3 and 4 in (c); $h_2=0.01$ , $(1/3)$ , $(2/3)$ , $(4/3)$ and $(8/3)$ in (e)). The dashed lines show the profile predicted by (higher-order) shallow-layer theory. The solid line in (d) shows the profile from (b). All the profiles are aligned at the surface height $h=8$ indicated by the dot-dashed lines. The overlaid plots in (d,e) show magnifications near the nose.

The profile of the special solution is displayed in figure 22(d). One awkward detail of this solution is that the steepening of the free surface near the fluid edge inevitably leads to the lowest $\alpha $ -line eventually curving downwards rather than upwards into the fluid (see Nye (Reference Nye1951) and Appendix B.4). Nye argued that this breakdown of the special solution could by cured by placing a plug below the distinguished $\alpha $ -line that curves initially upwards from the base before descending and intersecting the fluid edge, as illustrated by the shaded region in figure 22(c).

Unfortunately, Nye’s construction is not relevant to the dambreak problem shown in figure 21: his glacier snout is assumed to be in a state of horizontal compression. By contrast, as illustrated in figure 21(c), shallow viscoplastic dambreaks reach their the final state under horizontal expansion, a result that carries over to finite depths (Liu et al. Reference Liu, Balmforth, Hormozi and Hewitt2016). When the surface boundary conditions are adjusted to account for this difference, Nye’s construction fails to work: the slipline calculation no longer converges, but diverges from any special solution. The mathematical reason underlying the ‘stability’ of Nye’s construction has not previously been established; nor has the ‘instability’ of its extensional counterpart.

To provide the extensional version of Nye’s construction a different procedure is needed. One simple strategy for the task is outlined in Appendix B.4, and leads to the results shown in figure 22(e). This procedure begins with a guess for a section of the free-surface position, builds the resulting slipline field down to the base, then updates the surface position to correct the bottom boundary condition ( $\theta (x,0)=0$ ). The strategy can also be used to build solutions for horizontal compression; figure 22(c) demonstrates that the results agree with Nye’s construction, except near the end of the profile where the $\alpha $ -lines bend downward.

The breakdown of the slipline field near the fluid edge is also inevitable when the fluid is under horizontal extension (Appendix B.4). Awkwardly, once the $\alpha $ -lines begin to curve downwards it no longer becomes possible to extend the slipline field down from the surface to positions where $y=\theta =0$ , and the strategy of Appendix B.4 fails to produce a physical solution. For horizontal extension, the breakdown is more noticeable (see figure 22 e,f), and it does not seem possible to rectify the situation by adding a plug against the base to allow the slipline network to be continued to the fluid edge.

As an analytical alternative to sliplines, one can continue the asymptotics of lubrication theory to higher order to generate more accurate profiles than the leading-order result in (6.3). This exercise was accomplished by Dubash et al. (Reference Dubash, Balmforth, Slim and Cochard2009) and Liu et al. (Reference Liu, Balmforth, Hormozi and Hewitt2016), leading to the predictions

(6.4) \begin{equation} h =\left \{ \begin{array}{l} \sqrt {\dfrac {1}{4}\pi ^2\textit{Bi}^2+2\textit{Bi}(X-x)}-\dfrac {1}{2}\pi \textit{Bi} \cr \sqrt {2\textit{Bi}(X-x)}+\dfrac {1}{2}\pi \textit{Bi} \end{array}\right ., \end{equation}

for horizontal compression or extension, respectively. Note that the profile in extension ends at $x=X$ with a finite depth (indicating an implausible vertical cliff face at the edge), and both profiles reduce to parameterless curves when using the variables, $(x,h)/\textit{Bi}$ . The predictions in (6.4) are added to figure 22(d,e) for comparison with the slipline solutions.

Figure 23. Final profiles from a suite of numerical simulations with varying $\textit{Bi}$ that begin with either a square or triangle of unit area (Liu et al. Reference Liu, Balmforth, Hormozi and Hewitt2016). The insets in (a) show the two sets of final profiles. These are then scaled by $\textit{Bi}$ and replotted in the main panel. A magnification near the front is displayed in (b). The dot-dashed line shows the leading-order asymptotic prediction from (6.3). The dashed line is the improved profile from (6.4) and the solid line shows the construction from figure 22(c). The latter two curves are shifted horizontally to align them with the numerical simulations at the position where $h=10\textit{Bi}$ .

Finally, the slipline results and higher-order asymptotic predictions for horizontal extension are replotted in figure 23 and compared with the final profiles from a suite of numerical simulations with varying $\textit{Bi}$ provided by Liu et al. (Reference Liu, Balmforth, Hormozi and Hewitt2016). These simulations are conducted by adopting either a square or triangular initial shape with the same area (the latter being a right-angle triangle with width $\ell$ at the base, but height 2 at $x=0$ ). Aside from some relics from the initial shape, the final profiles are successfully collapsed on scaling by $\textit{Bi}$ . The scaled profiles match well with (6.4) and the slipline construction over distances exceeding about $4\textit{Bi}$ (i.e. four times the dimensional length scale ${\tau _{_{_P}}}/(\rho g)$ ). However, very close to the front, the latter terminate at finite height, whereas the simulations abruptly descend and bend around to create a small overhang. Evidently, (6.4) adequately describes the final shape of a slump except for the region at the nose, whose detailed structure remains unresolved.