2.1 Variational tools

The energy equation and two variational principles are useful in calculating the solution and studying limiting flows. Firstly, let us consider the space of admissible velocity solutions:

(2.10) $$\begin{eqnarray}\displaystyle {\mathcal{V}}=\{\boldsymbol{v}\boldsymbol{ : }\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{v}=0~\text{in }\unicode[STIX]{x1D6FA}\setminus \bar{X};~\boldsymbol{v}=0~\text{on }\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}\setminus \unicode[STIX]{x2202}X;~\boldsymbol{v}\rightarrow \boldsymbol{V}=V\boldsymbol{e}_{y}~\text{on }\unicode[STIX]{x2202}X\}, & & \displaystyle\end{eqnarray}$$

where $\boldsymbol{V}$ is an unspecified constant particle test velocity. For arbitrary $\boldsymbol{v},\boldsymbol{w}\in {\mathcal{V}}$ , the viscous dissipation, plastic dissipation and buoyancy work functionals, respectively, are

(2.11a ) $$\begin{eqnarray}\displaystyle & \displaystyle a(\boldsymbol{v},\boldsymbol{w})=\int _{\unicode[STIX]{x1D6FA}\setminus \bar{X}}\dot{\unicode[STIX]{x1D738}}(\boldsymbol{v})\boldsymbol{ : }\dot{\unicode[STIX]{x1D738}}(\boldsymbol{w})\,\text{d}A, & \displaystyle\end{eqnarray}$$

(2.11b ) $$\begin{eqnarray}\displaystyle & \displaystyle j(\boldsymbol{v})=\int _{\unicode[STIX]{x1D6FA}\setminus \bar{X}}\Vert \dot{\unicode[STIX]{x1D738}}(\boldsymbol{v})\Vert \,\text{d}A, & \displaystyle\end{eqnarray}$$

(2.11c ) $$\begin{eqnarray}\displaystyle & \displaystyle L(\boldsymbol{v})=\unicode[STIX]{x03C0}V=\unicode[STIX]{x03C0}\boldsymbol{V}\boldsymbol{\cdot }\boldsymbol{e}_{y}=-\int _{\unicode[STIX]{x1D6FA}\setminus \bar{X}}\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{e}_{y}\,\text{d}A. & \displaystyle\end{eqnarray}$$

$\unicode[STIX]{x1D6FA}\setminus \bar{X}$

(2.12) $$\begin{eqnarray}\displaystyle a(\boldsymbol{u},\boldsymbol{u})+Y\,j(\boldsymbol{u}) & = & \displaystyle \frac{\unicode[STIX]{x1D70C}}{1-\unicode[STIX]{x1D70C}}\int _{\unicode[STIX]{x1D6FA}\setminus \bar{X}}\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{e}_{y}\,\text{d}A+\int _{\unicode[STIX]{x2202}X}\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D748}\boldsymbol{\cdot }\boldsymbol{n}\,\text{d}s\nonumber\\ \displaystyle & = & \displaystyle -\frac{\unicode[STIX]{x1D70C}}{1-\unicode[STIX]{x1D70C}}\int _{\bar{X}}\boldsymbol{U}\boldsymbol{\cdot }\boldsymbol{e}_{y}\,\text{d}A+\int _{\unicode[STIX]{x2202}X}\boldsymbol{U}\boldsymbol{\cdot }\unicode[STIX]{x1D748}\boldsymbol{\cdot }\boldsymbol{n}\,\text{d}s\nonumber\\ \displaystyle & = & \displaystyle -\frac{\unicode[STIX]{x03C0}\unicode[STIX]{x1D70C}}{1-\unicode[STIX]{x1D70C}}\boldsymbol{U}\boldsymbol{\cdot }\boldsymbol{e}_{y}+\boldsymbol{U}\boldsymbol{\cdot }\int _{\unicode[STIX]{x2202}X}\unicode[STIX]{x1D748}\boldsymbol{\cdot }\boldsymbol{n}\,\text{d}s=\unicode[STIX]{x03C0}U=L(\boldsymbol{u}).\end{eqnarray}$$

Below, (2.12) is used to define the critical yield number.

The velocity minimization principle states that the solution $\boldsymbol{u}$ minimizes the functional ${\mathcal{J}}(\boldsymbol{v})$ over $\boldsymbol{v}\in {\mathcal{V}}$ :

(2.13) $$\begin{eqnarray}{\mathcal{J}}(\boldsymbol{v})=\frac{1}{2}a(\boldsymbol{v},\boldsymbol{v})+Y\,j(\boldsymbol{v})-\frac{\unicode[STIX]{x1D70C}}{1-\unicode[STIX]{x1D70C}}\int _{\unicode[STIX]{x1D6FA}\setminus \bar{X}}\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{e}_{y}\,\text{d}A-\int _{\unicode[STIX]{x2202}X}\unicode[STIX]{x1D70E}_{ij}(\boldsymbol{u})n_{j}V_{i}\,\text{d}s.\end{eqnarray}$$

Note that the traction of the solution at the particle surface (contained within the last term on the right-hand side) satisfies (2.8). Since $\boldsymbol{V}$ is constant, (2.9) applies to the pairing $(\boldsymbol{v},\boldsymbol{V})$ , and on using (2.8) we obtain

(2.14) $$\begin{eqnarray}\displaystyle -\frac{\unicode[STIX]{x1D70C}}{1-\unicode[STIX]{x1D70C}}\int _{\unicode[STIX]{x1D6FA}\setminus \bar{X}}\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{e}_{y}\,\text{d}A-\int _{\unicode[STIX]{x2202}X}\unicode[STIX]{x1D70E}_{ij}(\boldsymbol{u})n_{j}V_{i}\,\text{d}s=\frac{\unicode[STIX]{x03C0}\unicode[STIX]{x1D70C}V}{1-\unicode[STIX]{x1D70C}}-\boldsymbol{V}\boldsymbol{\cdot }\int _{\unicode[STIX]{x2202}X}\unicode[STIX]{x1D748}(\boldsymbol{u})\boldsymbol{\cdot }\boldsymbol{n}\,\text{d}s=-L(\boldsymbol{v}), & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

so that equivalently ${\mathcal{J}}(\boldsymbol{v})=(1/2)a(\boldsymbol{v},\boldsymbol{v})+Y\,j(\boldsymbol{v})-L(\boldsymbol{v})$ .

For the second variational principle, consider the space of admissible stresses ${\mathcal{S}}$ :

(2.15) $$\begin{eqnarray}\displaystyle {\mathcal{S}}=\left\{\unicode[STIX]{x1D748}^{\prime }=\unicode[STIX]{x1D749}^{\prime }-p^{\prime }\boldsymbol{I}\boldsymbol{ : }\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D748}^{\prime }+\frac{\unicode[STIX]{x1D70C}}{1-\unicode[STIX]{x1D70C}}\boldsymbol{e}_{y}=0~\text{in }\unicode[STIX]{x1D6FA}\setminus \bar{X};~\int _{\unicode[STIX]{x2202}X}\unicode[STIX]{x1D748}^{\prime }\boldsymbol{\cdot }\boldsymbol{n}\,\text{d}s=\frac{\unicode[STIX]{x03C0}}{1-\unicode[STIX]{x1D70C}}\boldsymbol{e}_{y}\right\}, & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

i.e. symmetric stress tensors satisfying the Stokes equations and the force balance (2.8). The solution $\unicode[STIX]{x1D748}$ maximizes over $\unicode[STIX]{x1D748}^{\prime }\in {\mathcal{S}}$ the functional

(2.16) $$\begin{eqnarray}{\mathcal{K}}(\unicode[STIX]{x1D748}^{\prime })=\int _{\unicode[STIX]{x2202}X}\unicode[STIX]{x1D70E}_{ij}^{\prime }n_{j}U_{i}\,\text{d}s-\frac{1}{2}\int _{\unicode[STIX]{x1D6FA}\setminus \bar{X}}(\Vert \unicode[STIX]{x1D749}^{\prime }\Vert -Y)_{+}^{2}\,\text{d}A,\end{eqnarray}$$

where $(f)_{+}$ denotes the positive part of $f$ . Furthermore, for any admissible velocity and stress we have

(2.17) $$\begin{eqnarray}{\mathcal{K}}(\unicode[STIX]{x1D748}^{\prime })\leqslant {\mathcal{K}}(\unicode[STIX]{x1D748})={\mathcal{J}}(\boldsymbol{u})\leqslant {\mathcal{J}}(\boldsymbol{v}).\end{eqnarray}$$

These results date back to Prager (Reference Prager1954). A general derivation is explained at length in Huilgol (Reference Huilgol2015).

2.2 Yielding of flows

An important feature of yield-stress fluids is their ability to hold particles statically stable under certain circumstances. From the physical perspective, as discussed in our scaling, the density difference between particle and fluid leads to shear stresses of order $\hat{\unicode[STIX]{x1D70F}}_{b}=(\hat{\unicode[STIX]{x1D70C}}_{p}-\hat{\unicode[STIX]{x1D70C}}_{f}){\hat{g}}\hat{L}$ . It is physically intuitive that, if the yield stress $\hat{\unicode[STIX]{x1D70F}}_{Y}$ is sufficiently large (implicitly meaning $\hat{\unicode[STIX]{x1D70F}}_{Y}\gg (\hat{\unicode[STIX]{x1D70C}}_{p}-\hat{\unicode[STIX]{x1D70C}}_{f}){\hat{g}}\hat{L}$ ), then the yield stress will not be exceeded anywhere in the flow. Since the velocity is zero on $\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}$ there would then be no motion in the fluid and hence no motion of the particle. The ratio $\hat{\unicode[STIX]{x1D70F}}_{Y}/\hat{\unicode[STIX]{x1D70F}}_{b}$ is the yield number and the above argument predicts stationary flow for sufficiently large $Y$ . We can verify this intuition analytically, using the tools of the previous section (see e.g. Putz & Frigaard Reference Putz and Frigaard2010).

Secondly, we might re-examine the maximization (2.16). Since each admissible $\unicode[STIX]{x1D748}^{\prime }$ satisfies (2.8), we may also write

(2.18) $$\begin{eqnarray}{\mathcal{K}}(\unicode[STIX]{x1D748}^{\prime })=\frac{1}{1-\unicode[STIX]{x1D70C}}L(\boldsymbol{u})-\frac{1}{2}\int _{\unicode[STIX]{x1D6FA}\setminus \bar{X}}(\Vert \unicode[STIX]{x1D749}^{\prime }\Vert -Y)_{+}^{2}\,\text{d}A.\end{eqnarray}$$

The first term is independent of $\unicode[STIX]{x1D748}^{\prime }\in {\mathcal{S}}$ and the second term is negative definite. Any admissible stress field with $\Vert \unicode[STIX]{x1D749}^{\prime }\Vert \leqslant Y$ will therefore maximize ${\mathcal{K}}$ . Equally, if any such $\unicode[STIX]{x1D748}^{\prime }\in {\mathcal{S}}$ can be found, then it follows that the true stress must also have $\Vert \unicode[STIX]{x1D749}\Vert \leqslant Y$ . This in turn implies that the solution stress does not exceed the yield stress, hence that the strain rate is zero everywhere, and consequently that $\boldsymbol{u}=0$ . Via the equivalence (2.17), we see that ${\mathcal{K}}(\unicode[STIX]{x1D748})={\mathcal{J}}(\boldsymbol{u})=0$ , in this limit. Thus, the condition for the particle not to yield the fluid is independent of the precise type of fluid (Bingham, Casson, Herschel–Bulkley, Robertson–Stiff, etc.): any viscoplastic fluid remains preferentially static.

The energy equation can be used to define the minimal limiting yield stress for static flow, via a critical value of $Y$ :

(2.19) $$\begin{eqnarray}Y_{c}=\sup _{\boldsymbol{v}\in {\mathcal{V}},\,\boldsymbol{v}\not =0}\frac{L(\boldsymbol{v})}{j(\boldsymbol{v})},\end{eqnarray}$$

with $\boldsymbol{u}=\boldsymbol{U}=0$ if $Y\geqslant Y_{c}$ (Putz & Frigaard Reference Putz and Frigaard2010). The variational form of the velocity minimization can then be used to show that as $Y\rightarrow Y_{c}^{-}$ the viscous dissipation of the solution $\boldsymbol{u}$ will decay at least one order faster than the plastic dissipation. More precisely, as $Y\rightarrow Y_{c}^{-}$ ,

(2.20a,b ) $$\begin{eqnarray}\displaystyle a(\boldsymbol{u},\boldsymbol{u})=O([Y_{c}-Y]^{2}),\quad j(\boldsymbol{u})=O(Y_{c}-Y), & & \displaystyle\end{eqnarray}$$

(2.21) $$\begin{eqnarray}\displaystyle \frac{a(\boldsymbol{u},\boldsymbol{u})}{j(\boldsymbol{u})}=O(Y_{c}-Y). & & \displaystyle\end{eqnarray}$$

Returning to (2.12), we see that the plastic dissipation must balance with the buoyancy work as $Y\rightarrow Y_{c}^{-}$ , i.e.

(2.22) $$\begin{eqnarray}Y\,j(\boldsymbol{u})\sim L(\boldsymbol{u})=\unicode[STIX]{x03C0}U,\end{eqnarray}$$

giving us directly the decay of the settling velocity as $Y\rightarrow Y_{c}^{-}$ .

2.3 Resistance formulation

The formulation [M] adopted is as a mobility problem. A different formulation is as a resistance problem [R], where we suppose that a constant particle velocity $\hat{\boldsymbol{U}}_{p}$ is imposed and the force on the particle $\hat{\boldsymbol{F}}$ is measured. The particle velocity scales the velocity and is used to define a viscous stress scale $\hat{\unicode[STIX]{x1D70F}}_{p}=\hat{\unicode[STIX]{x1D707}}\hat{U} _{p}/\hat{L}$ , with which the Stokes equations are scaled:

(2.23a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D748}^{\ast }=0,\quad \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}^{\ast }=0\quad \text{in }\unicode[STIX]{x1D6FA}\setminus \bar{X}.\end{eqnarray}$$

Lengths are scaled as previously. The ‘asterisk’ is used to distinguish variables from the [M]-scaled variables, where needed. No body forces are present in the fluid domain. The constitutive laws are again (2.4) with $Y$ replaced by $B$ , which is the Bingham number: $B=\hat{\unicode[STIX]{x1D70F}}_{Y}\hat{L}/\hat{\unicode[STIX]{x1D707}}\hat{U} _{p}$ . The motion is driven only by boundary conditions on the particle, i.e.  $\boldsymbol{u}^{\ast }=\boldsymbol{e}_{y}$ on $\unicode[STIX]{x2202}X$ .

In order to have a correspondence between problems [M] and [R], we align the direction of particle motion with $\boldsymbol{e}_{y}$ , assumed to be a symmetry axis of the particle. The mapping between problems [M] and [R] is captured by

(2.24a-c ) $$\begin{eqnarray}\boldsymbol{u}=\frac{Y}{B}\boldsymbol{u}^{\ast },\quad U=\frac{Y}{B},\quad \unicode[STIX]{x03C0}\frac{B}{Y}=F^{\ast },\end{eqnarray}$$

where the scaled force exerted on the particle has magnitude $F^{\ast }$ and is in the direction  $-\boldsymbol{e}_{y}$ .

One reason for introducing [R] is that it may be more convenient to compute the solution of problem [R] (rather than [M]), and use this to evaluate $Y_{c}$ . In controlling relative numerical errors, the velocities in problem [R] are strictly $O(1)$ . Secondly, as argued by Putz & Frigaard (Reference Putz and Frigaard2010), optimality of $L(\boldsymbol{v})/j(\boldsymbol{v})$ in (2.19) is attained by the solution $\boldsymbol{u}$ to [M], in the limit that $Y\rightarrow Y_{c}^{-}$ . On rescaling $\boldsymbol{u}$ with $U$ , we arrive at the equivalent solution to [R], for $B=Y/U$ . The functional $L(\boldsymbol{v})/j(\boldsymbol{v})$ in (2.19) also rescales and thus $Y_{c}$ may be calculated alternatively as

(2.25) $$\begin{eqnarray}Y_{c}=\lim _{Y\rightarrow Y_{c}^{-}}\frac{L(\boldsymbol{u})}{j(\boldsymbol{u})}=\lim _{Y\rightarrow Y_{c}^{-}}\frac{\unicode[STIX]{x03C0}U}{j(\boldsymbol{u})}=\lim _{B\rightarrow \infty }\frac{\unicode[STIX]{x03C0}}{j(\boldsymbol{u}^{\ast })}.\end{eqnarray}$$

Although finding the solution $\boldsymbol{u}^{\ast }$ , as $B\rightarrow \infty$ , may become numerically difficult, at least $j(\boldsymbol{u}^{\ast })$ is an integral quantity and should be more stable to compute than a pointwise quantity, allowing the asymptotic behaviour at large $B$ to be estimated effectively from a numerical solution.

2.4 The plastic drag coefficient

A different way of evaluating the critical stopping limit, which is useful for comparing between problems, involves the drag coefficient. Conventionally, a drag coefficient is defined by scaling the force on an object with an inertial stress scale and cross-sectional area. Here we consider inertialess flows close to yielding. The plastic drag coefficient instead uses the yield stress in place of an inertial stress. For the 2D flow, in either formulation:

(2.26) $$\begin{eqnarray}\displaystyle C_{D}^{p}=\left\{\begin{array}{@{}l@{}}\left[\displaystyle \frac{\hat{F}}{\hat{\ell }_{\bot }\hat{\unicode[STIX]{x1D70F}}_{Y}}\right]^{[R]}=\left[\displaystyle \frac{F^{\ast }}{\ell _{\bot }B}\right]^{[R]}\quad \text{for problem [R]},\quad \\ \left[\displaystyle \frac{\hat{\unicode[STIX]{x1D70C}}{\hat{g}}\hat{A}_{p}}{\hat{\ell }_{\bot }\hat{\unicode[STIX]{x1D70F}}_{Y}}\right]^{[M]}=\left[\displaystyle \frac{\unicode[STIX]{x03C0}}{\ell _{\bot }Y}\right]^{[M]}\quad \text{for problem [M]},\quad \end{array}\right. & & \displaystyle\end{eqnarray}$$

where $\hat{\ell }_{\bot }=\ell _{\bot }\hat{L}$ is the linear dimension of the particle perpendicular to the flow direction. Other similar problems may be compared by scaling the particle force $\hat{F}$ , as above with $\hat{\ell }_{\bot }\hat{\unicode[STIX]{x1D70F}}_{Y}$ . We see from the above again that $F^{\ast }=\unicode[STIX]{x03C0}B/Y$ . As we increase the yield stress and approach $Y\rightarrow Y_{c}^{-}$ , the plastic drag coefficient approaches a critical value:

(2.27) $$\begin{eqnarray}C_{D,c}^{p}=\frac{\unicode[STIX]{x03C0}}{\ell _{\bot }Y_{c}}=\left[C_{D}^{p}\right]_{Y\rightarrow Y_{c}^{-}}^{[M]}=\left[C_{D}^{p}\right]_{B\rightarrow \infty }^{[R]}.\end{eqnarray}$$

Critical values $C_{D,c}^{p}$ have been reported for various particles with simple geometry. For example, Tokpavi et al. (Reference Tokpavi, Magnin and Jay2008) calculated $C_{D,c}^{p}$ for a circular disc and Nirmalkar et al. (Reference Nirmalkar, Chhabra and Poole2012) calculated $C_{D,c}^{p}$ for a square. Both studies used the resistance formulation [R] and extrapolated to the critical value. Alternatively, Putz & Frigaard (Reference Putz and Frigaard2010) calculated $C_{D,c}^{p}$ for 2D elliptical particles using problem [M].

2.5 The yield limit and perfect plasticity

Under either formulation [M] or [R] we have seen that the viscous dissipation terms apparently play a minor role in the limit of zero flow. If we neglect the viscous terms in the constitutive law, in yielded regions we have

(2.28) $$\begin{eqnarray}\unicode[STIX]{x1D749}\approx B\frac{\dot{\unicode[STIX]{x1D738}}}{\Vert \dot{\unicode[STIX]{x1D738}}\Vert }.\end{eqnarray}$$

This is reminiscent of a rigid–perfectly plastic material (Hill Reference Hill1950), which may be considered as the limit of the Prandtl–Reuss model for perfectly plastic materials, in the case that the elastic modulus is infinitely large. The stress at no point can exceed the yield stress in (2.28). If the stress is strictly lower than the yield stress, the material is modelled as rigid, but formally (2.28) is an approximation in the limit of vanishingly small strain rates.

In a 2D flow, the steady momentum equations plus the von Mises yield criterion form a closed set of hyperbolic equations for the plastic regions of flow (see Chakrabarty Reference Chakrabarty2012). Considering now only problem [R] we have

(2.29a ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\tilde{\unicode[STIX]{x1D749}}-\unicode[STIX]{x1D735}\tilde{p}=0, & \displaystyle\end{eqnarray}$$

(2.29b ) $$\begin{eqnarray}\displaystyle & \displaystyle \Vert \tilde{\unicode[STIX]{x1D749}}\Vert =B. & \displaystyle\end{eqnarray}$$

$\unicode[STIX]{x1D6FC}$

$\unicode[STIX]{x1D6FD}$

$\unicode[STIX]{x1D719}$

$\unicode[STIX]{x1D6FC}$

$x$

(2.30a-c ) $$\begin{eqnarray}\tilde{\unicode[STIX]{x1D70E}}_{xx}=-\tilde{p}-B\sin (2\unicode[STIX]{x1D719}),\quad \tilde{\unicode[STIX]{x1D70E}}_{yy}=-\tilde{p}+B\sin (2\unicode[STIX]{x1D719}),\quad \tilde{\unicode[STIX]{x1D70F}}_{xy}=B\cos (2\unicode[STIX]{x1D719}).\end{eqnarray}$$

The principle stress directions bisect the right angle between the characteristics. Further features of the characteristics are described in Hill (Reference Hill1950) and Chakrabarty (Reference Chakrabarty2012).

For the present discussion, the most interesting fact is that calculating admissible stress and velocity fields from analysis of the characteristics (or indeed any admissible field constructed by any other means) allows us to make lower- and upper-bound estimates of the force on a particle translating at unit speed within the plastic medium; see appendix A. In turn, these lead to estimates of the limiting plastic drag coefficient and hence $Y_{c}$ . Note that, because we consider a 2D problem, the drag forces we estimate below (and elsewhere in the paper) are forces per unit depth of the particle in the third $z$ -direction.

1. (1) Lower bound: as outlined in (A 2) the load calculated by using an admissible stress field is a lower bound for the actual force, i.e.

   (2.31) $$\begin{eqnarray}[\tilde{F}^{p}]_{L}=\int _{\unicode[STIX]{x1D6E4}}\tilde{\unicode[STIX]{x1D748}}\boldsymbol{\cdot }\boldsymbol{n}\,\text{d}S.\end{eqnarray}$$
   Here $\unicode[STIX]{x1D6E4}$ is the boundary of the region $\unicode[STIX]{x1D6E5}$ that consists of the particle plus any plugs attached to the particle surface and translating with the same velocity. Evidently, $\unicode[STIX]{x1D6E4}$ may be replaced by the boundary of the particle. The characteristic network simply provides a convenient method for computing an admissible stress field.

2. (2) Upper bound: as outlined in (A 4) the mechanical work dissipated in any admissible velocity field $\tilde{\boldsymbol{v}}$ is an upper bound for the actual force, i.e.

   (2.32) $$\begin{eqnarray}[\tilde{F}^{p}]_{U}=B\left(\int _{\unicode[STIX]{x1D6FA}\setminus \bar{\unicode[STIX]{x1D6E5}}}\Vert \tilde{\dot{\unicode[STIX]{x1D738}}}^{\prime }\Vert \,\text{d}V+\int _{S_{d}(\tilde{\boldsymbol{v}})}|\unicode[STIX]{x0394}\tilde{\boldsymbol{v}}|\,\text{d}S_{d}(\tilde{\boldsymbol{v}})\right).\end{eqnarray}$$
   Here $S_{d}(\tilde{\boldsymbol{v}})$ contains all curves that contain velocity discontinuities in $\unicode[STIX]{x1D6FA}\setminus \bar{\unicode[STIX]{x1D6E5}}$ and $\unicode[STIX]{x0394}\tilde{\boldsymbol{v}}$ represents the velocity discontinuity vector. We may compute an admissible velocity field in any way. Again, we may replace $\unicode[STIX]{x1D6E5}$ with $X$ , above.

If the lower and upper bounds are equal, the calculated drag force is exact. In most problems there is an uncertainty between the two bounds and this usually stems from the velocity field giving a poor estimate. This can be reduced by finding a better admissible velocity field, but this becomes hard for complex geometries. A final point to acknowledge here is that the admissible stress fields that we calculate are incomplete. In constructing the characteristic meshes, we typically find a bounded region around the particle on which the stress field (for calculating the lower bound) is resolved. Our primary interest is to make comparisons of the perfectly plastic flow in this region with that of the analogous viscoplastic flow, and in particular in the limit of zero flow. Thus, we do not attempt to construct stress fields for (2.29) outside of the characteristic net computed. It is simply assumed that the computed stress field can be coupled with such a stress field in the exterior. In cases where the perfectly plastic and limiting viscoplastic solutions appear to coincide, we comment that the viscoplastic solution has computed a stress field with $\Vert \unicode[STIX]{x1D749}\Vert \leqslant B$ in the unyielded exterior region.

4.1 The limit $\unicode[STIX]{x1D712}\rightarrow 0$

The limit $\unicode[STIX]{x1D712}\rightarrow 0$ gives the greatest agreement between drag coefficients computed from the viscoplastic flow and the characteristics analysis. The structure of the solution that we have seen in our viscoplastic computations consists of (up to) five regions as $\unicode[STIX]{x1D712}\rightarrow 0$ . Referring to one quarter of the flow domain, these regions are as follows.

1. (i) A rigid plug attached to the front (and rear) of the particle, of width ${\sim}\ell _{\bot }/2$ as $\unicode[STIX]{x1D712}\rightarrow 0$ , and forming a triangular region fore–aft of the particle, making an angle $\unicode[STIX]{x03C0}/4$ with the $y$ -axis.

2. (ii) A rotating zone that is centred at the ‘corner’ of the rigid plug, where it contacts the particle. The fluid rotates on circular streamlines about this point and appears to have approximately constant speed in the azimuthal direction. The radius of this region is ${\sim}\ell _{\bot }/\sqrt{2}$ and the azimuthal speed is ${\sim}1/\sqrt{2}$ .

3. (iii) A thin shear layer between the attached rigid plug and the rotating zone. The length of this boundary layer is $\ell _{\bot }/\sqrt{2}$ . Across the layer the normal component of velocity is continuous and the tangential drop from $1/\sqrt{2}$ to $0$ .

4. (iv) An outer boundary layer, surrounding the rotating zone, across which the azimuthal velocity drops from $1/\sqrt{2}$ to $0$ . The circumferential length of this boundary layer is $3\unicode[STIX]{x03C0}\ell _{\bot }/(4\sqrt{2})$ . We have not studied the thickness, but it does not appear to grow as $\unicode[STIX]{x1D712}\rightarrow 0$ .

5. (v) A ‘side region’ that connects the rotating regions fore and aft. This region has width ${\sim}\ell _{\bot }/2$ and height ${\sim}\unicode[STIX]{x1D712}\ell _{||}$ (ellipse) or ${\sim}\ell _{||}$ (rectangle). This region is absent in the case of the diamond shape.

We can estimate the contributions to $j(\boldsymbol{u}^{\ast })$ from each of the above regions as follows: (i) $0$ ; (ii) ${\sim}3\unicode[STIX]{x03C0}\ell _{\bot }/2$ ; (iii) ${\sim}2\ell _{\bot }$ ; (iv) ${\sim}3\unicode[STIX]{x03C0}\ell _{\bot }/2$ ; (v) ${\sim}\unicode[STIX]{x1D712}\ell _{||}/2$ (ellipse) or ${\sim}\ell _{||}/2$ (rectangle). Combining these we find: $C_{D,c}^{p}=\lim _{B\rightarrow \infty }[j(\boldsymbol{u}^{\ast })/\ell _{\bot }]\sim 2+3\unicode[STIX]{x03C0}$ , for all three geometries as $\unicode[STIX]{x1D712}\rightarrow 0$ , which agrees with the characteristics estimates at leading order.

In order to estimate convergence with respect to $\unicode[STIX]{x1D712}$ , we would need to know the next order of terms in each of the above regions, which is not apparent. Figure 17 plots the critical plastic drag coefficient for different $\unicode[STIX]{x1D712}$ , as $\unicode[STIX]{x1D712}\rightarrow 0$ , computed with $B=10^{4}$ . It is rather interesting that the diamond shape shows no variation with $\unicode[STIX]{x1D712}$ for $\unicode[STIX]{x1D712}<1$ . As the diamond has no ‘sides’ and only the tip region, the main role of $\unicode[STIX]{x1D712}$ is to determine the width of the particle, from which we have seen there results an attached fore–aft plug. The unyielded plug is thus indistinguishable from the particle and we see no variation with  $\unicode[STIX]{x1D712}$ .

Figure 17.

Critical plastic drag coefficient for $\unicode[STIX]{x1D712}\rightarrow 0$ , computed with $B=10^{4}$ .

By the same analogy, the parts of the ellipse and rectangle that lie within the fore–aft angular plug regions as $\unicode[STIX]{x1D712}\rightarrow 0$ should not have any influence on the flow convergence with $\unicode[STIX]{x1D712}$ . This suggests that $\unicode[STIX]{x1D712}$ dependence comes only from the side region, denoted (v) above. For the ellipse the surface tangent that is at $\unicode[STIX]{x03C0}/4$ to the $y$ -axis meets the surface at

(4.1) $$\begin{eqnarray}\displaystyle y=\frac{\ell _{||}}{2}\sqrt{\frac{\unicode[STIX]{x1D712}^{2}}{1+\unicode[STIX]{x1D712}^{2}}}\sim \frac{\unicode[STIX]{x1D712}\ell _{||}}{2}. & & \displaystyle\end{eqnarray}$$

The contribution from this side region to $j(\boldsymbol{u}^{\ast })$ is thus ${\sim}\unicode[STIX]{x1D712}\ell _{||}$ and to $C_{D,c}^{p}$ is ${\sim}\unicode[STIX]{x1D712}\ell _{||}/\ell _{\bot }\sim \unicode[STIX]{x1D712}^{2}$ . For the rectangle, we have instead that the contribution to $C_{D,c}^{p}$ is ${\sim}\ell _{||}/\ell _{\bot }\sim \unicode[STIX]{x1D712}$ . We see that these estimates do agree with the (approximately) quadratic and linear variations in figure 17, for ellipse and rectangle, respectively.

Regarding $B$ dependence, generally this must be examined numerically. However, for the diamond we may progress via a rescaling argument as follows. Defining $\check{\boldsymbol{x}}=2\boldsymbol{x}/\ell _{\bot }$ we rescale the lengths while keeping velocities unchanged in problem [R]. If we also rescale $\check{\unicode[STIX]{x1D70F}}_{ij}^{\ast }=(\ell _{\bot }/2)\unicode[STIX]{x1D70F}_{ij}^{\ast }$ and $\check{p}^{\ast }=(\ell _{\bot }/2)p^{\ast }$ , the Stokes equations are satisfied by the rescaled stress tensor, which is characterized by a single parameter,

(4.2) $$\begin{eqnarray}\check{B}=\frac{\ell _{\bot }}{2}B.\end{eqnarray}$$

The re-mapped diamond intersects the $\check{x}$ -axis at $\pm 1$ and the $\check{y}$ -axis at $\pm \unicode[STIX]{x1D712}$ . For any fixed $\unicode[STIX]{x1D712}<1$ the limit $B\rightarrow \infty$ is the same as $\check{B}\rightarrow \infty$ .

We have already seen that as $\check{B}\rightarrow \infty$ the solution will have a rigid plug attached to the diamond that intersects the $\check{y}$ -axis at $\pm 1$ , and then regions (ii)–(iv), as described above. The contributions to $j(\check{\boldsymbol{u}}^{\ast })$ from regions (ii)–(iv) are $3\unicode[STIX]{x03C0}$ , $4$ and $3\unicode[STIX]{x03C0}$ , respectively. However, since now $\check{\boldsymbol{u}}^{\ast }$ depends only on $\check{B}$ , we may postulate convergence,

(4.3) $$\begin{eqnarray}j(\check{\boldsymbol{u}}^{\ast })\sim 4+6\unicode[STIX]{x03C0}+2A\check{B}^{-t}\quad \text{as}~\check{B}\rightarrow \infty ,\end{eqnarray}$$

for constants $A$ and $t>0$ . On rescaling lengths, $j(\boldsymbol{u}^{\ast })=(\ell _{\bot }/2)j(\check{\boldsymbol{u}}^{\ast })\sim \ell _{\bot }(2+3\unicode[STIX]{x03C0}+A\check{B}^{-t}),~\text{as}~\check{B}\rightarrow \infty$ . From this we deduce that

(4.4) $$\begin{eqnarray}C_{D,c}^{p}=\lim _{B\rightarrow \infty }\frac{j(\boldsymbol{u}^{\ast })}{\ell _{\bot }}\sim 2+3\unicode[STIX]{x03C0}+O(B^{-t}).\end{eqnarray}$$

The exponent $t$ can be deduced economically from computation of $\check{\boldsymbol{u}}^{\ast }$ as $\check{B}\rightarrow \infty$ . We have found $t\approx 0.70$ from fitting a power-law expression to our data.

4.2 A general framework

We would like to understand the limit $Y\rightarrow Y_{c}^{-}$ . The rescaling leading to (2.24) clarifies that $Y\rightarrow Y_{c}^{-}$ and $B\rightarrow \infty$ are the same limit, but does not actually determine the relationship between the two parameters. As well as clarifying this relationship, we would like to be able to determine the decay rates. The only general results we have are the decay bounds (2.20) and (2.21) for problem [M], as $Y\rightarrow Y_{c}^{-}$ .

To this end, let us suppose that the particle velocity $U\sim b(\unicode[STIX]{x1D712})(1-Y/Y_{c})^{\unicode[STIX]{x1D708}}$ , for some positive exponent $\unicode[STIX]{x1D708}=\unicode[STIX]{x1D708}(\unicode[STIX]{x1D712})$ . The bound (2.20) implies that $\unicode[STIX]{x1D708}\geqslant 1$ and (2.22) gives us $j(\boldsymbol{u})\sim b(\unicode[STIX]{x1D712})(1-Y/Y_{c})^{\unicode[STIX]{x1D708}}\unicode[STIX]{x03C0}/Y_{c}$ . The rescaling of (2.24) then tells us that

(4.5) $$\begin{eqnarray}B=\frac{Y}{U}\sim \frac{Y_{c}(\unicode[STIX]{x1D712})}{b(\unicode[STIX]{x1D712})}(1-Y/Y_{c})^{-\unicode[STIX]{x1D708}},\end{eqnarray}$$

which couples the limits $Y\rightarrow Y_{c}^{-}$ and $B\rightarrow \infty$ . To find $\unicode[STIX]{x1D708}$ we must estimate the limiting solutions.

Suppose we compute problem [R] at large $B$ and substitute the velocity $\boldsymbol{u}^{\ast }$ into (2.25), i.e. we approximate $Y_{c}$ using the limiting values of $j(\boldsymbol{u}^{\ast })$ . If the contributions to $j(\boldsymbol{u}^{\ast })$ come from the entire flow domain, it seems clear that determination of $\unicode[STIX]{x1D708}$ is best carried out numerically. This occurs primarily for $\unicode[STIX]{x1D712}\sim O(1)$ . The best-studied example is the circular disc. Looking back at figure 4 it becomes evident that $j(\boldsymbol{u}^{\ast })$ will have contributions from at least three regions: the boundary layer close to the particle, the outer shear layer, and the zone of near-constant azimuthal velocity. The plug region converges to a fixed size as $B\rightarrow \infty$ (Tokpavi et al. Reference Tokpavi, Magnin and Jay2008), and so each of these regions continues to contribute at leading order. Tokpavi et al. (Reference Tokpavi, Magnin and Jay2008) find that the boundary layer thickness between the small rigid plug on the disc equator and the particle scales like $B^{-0.52}$ and the convergence of their plastic drag coefficient to its limiting value scales like $B^{-0.55}$ (i.e.  $\unicode[STIX]{x1D708}\approx 1.82$ ). Our results (see figure 8 b) also have $\unicode[STIX]{x1D708}\approx 2$ , consistent with Tokpavi et al. (Reference Tokpavi, Magnin and Jay2008). Putz & Frigaard (Reference Putz and Frigaard2010) calculate the limiting exponents for elliptical particles at $\unicode[STIX]{x1D712}=0.5$ , 1 and 2, as $\unicode[STIX]{x1D708}=1.37\pm 0.23$ ,  $1.20\pm 0.26$ and $1.84\pm 0.25$ , respectively. The convergence rate for $\unicode[STIX]{x1D712}=1$ does not agree with our present computations, nor with those of Tokpavi et al. (Reference Tokpavi, Magnin and Jay2008), although the values for $\unicode[STIX]{x1D708}=0.5$ and $2$ are consistent with our calculations. Evidently, calculating each such $\unicode[STIX]{x1D712}\sim O(1)$ geometry is laborious. The small- $\unicode[STIX]{x1D712}$ limit offers some simplification. The procedure outlined for the diamond, leading to the estimate $\unicode[STIX]{x1D708}\approx 1.43$ for all $\unicode[STIX]{x1D712}<1$ ( $t\approx 0.70$ ), should also be applicable to the rectangle and ellipse for sufficiently small $\unicode[STIX]{x1D712}$ (as the side regions vanish). Large $\unicode[STIX]{x1D712}$ also potentially allows for simplification (see § 4.3 below).

Lastly, note that the separation of geometric dependence from the decay in the solution as $Y\rightarrow Y_{c}^{-}$ , leading to (4.5), is justified by numerous examples. Classical 1D flows such as Poiseuille flow (pipe or planar) have $\unicode[STIX]{x1D708}=2$ . This also seems to be the case for more complex flows along uniform ducts; for example, Frigaard & Scherzer (Reference Frigaard and Scherzer2000) find $\unicode[STIX]{x1D708}=2$ for buoyancy-driven two-fluid exchange flows in circular pipes, with either concentric or stratified interfaces. These flows all have a single shear stress active along the bounding surface or at an interface, which seems to lead to $\unicode[STIX]{x1D708}=2$ . For 2D flows along uneven channels, Roustaei et al. (Reference Roustaei, Chevalier, Talon and Frigaard2016) find both $\unicode[STIX]{x1D708}\approx 2$ and $\unicode[STIX]{x1D708}\approx 1.5$ , depending on the channel shape, with the latter relating to limiting flows with an extensional stress component.

4.3 Long particles

The only way in which we might make general deductions about the limit $Y\rightarrow Y_{c}^{-}$ is if the viscoplastic solutions develop distinct geometric forms, in which predictable regions of the flow make the dominant contribution to $j(\boldsymbol{u}^{\ast })$ . This might be true in the case of boundary layers, which occur in at least two situations. Firstly, for the ellipse and rectangle at large $\unicode[STIX]{x1D712}$ , we have seen that there is a sharp transition from the particle velocity to the velocity of the large plug regions along the sides of the particle, that are in rigid rotation. Secondly, for the diamond the particle transitions into a fan-like structure over a relatively thin layer; see figure 16(a). We consider these two cases below. Note that in both cases the length of the boundary layer is expected to scale with  $\ell _{||}$ .

4.3.1 Ellipse and rectangle

At large $\unicode[STIX]{x1D712}$ , the general solution structure is as illustrated in figure 13(a). Along the sides of the particle there is an inner boundary layer (width $\unicode[STIX]{x1D6FF}_{i}$ ), then the rigid plug and finally an outer boundary layer (width $\unicode[STIX]{x1D6FF}_{o}$ ). The rigid plug has width $\ell _{p}$ and rotates with angular velocity $\unicode[STIX]{x1D714}_{p}$ . While the rigid plug does not contribute to $j(\boldsymbol{u}^{\ast })$ , the outer boundary layer does. Equally, there is another contribution from near the tip of the particle, where the fluid pushed out from in front of the plug and dragged upwards in the boundary layer turns and moves back downwards into the plug and outer boundary layer. Thus, we may write $j(\boldsymbol{u}^{\ast })=j_{i}(\boldsymbol{u}^{\ast })+j_{o}(\boldsymbol{u}^{\ast })+j_{t}(\boldsymbol{u}^{\ast })$ , denoting the contributions from the three regions.

Writing $\boldsymbol{u}^{\ast }=(u^{\ast },v^{\ast })$ , we may evaluate the no-flux condition (2.9) along the $x$ -axis for problem [R]:

(4.6) $$\begin{eqnarray}\displaystyle 0 & = & \displaystyle \frac{\ell _{\bot }}{2}+\int _{\ell _{\bot }/2}^{\ell _{\bot }/2+\unicode[STIX]{x1D6FF}_{i}}v^{\ast }(x,0)\,\text{d}x+\int _{\ell _{\bot }/2+\unicode[STIX]{x1D6FF}_{i}}^{\ell _{\bot }/2+\unicode[STIX]{x1D6FF}_{i}+\ell _{p}}v^{\ast }(x,0)\,\text{d}x+\int _{\ell _{\bot }/2+\unicode[STIX]{x1D6FF}_{i}+\ell _{p}}^{\ell _{\bot }/2+\unicode[STIX]{x1D6FF}_{i}+\ell _{p}+\unicode[STIX]{x1D6FF}_{o}}v^{\ast }(x,0)\,\text{d}x\nonumber\\ \displaystyle & = & \displaystyle O(\unicode[STIX]{x1D712}^{-1/2})+O(\unicode[STIX]{x1D6FF}_{i})-O(\unicode[STIX]{x1D714}_{p}\ell _{p}^{2})-O(\unicode[STIX]{x1D714}_{p}\ell _{p}\unicode[STIX]{x1D6FF}_{o}).\end{eqnarray}$$

In the limit of $\unicode[STIX]{x1D712}\rightarrow \infty$ , our computations of the viscoplastic solution suggest that $\unicode[STIX]{x1D714}_{p}\sim \unicode[STIX]{x1D712}^{-1}$ and $\ell _{p}\sim \unicode[STIX]{x1D712}^{1/2}$ (see figure 18). This suggests that, as $\unicode[STIX]{x1D712}\rightarrow \infty$ , the second and third terms on the right-hand side should balance: $\unicode[STIX]{x1D6FF}_{i}\sim O(1)$ . To leading order, the mass flux in the yielded boundary layer dragged upwards by the plate will turn back within the tip region into the rigidly rotating plug.

Figure 18.

Geometric behaviour at large $\unicode[STIX]{x1D712}$ , computed for the rectangular particle at $B=10^{4}$ using problem [R]: (a) $\unicode[STIX]{x1D714}_{p}$ ; (b) $\ell _{p}$ (symbols from computation and solid lines fitted, with the indicated slopes).

We may also evaluate $j_{i}(\boldsymbol{u}^{\ast })$ , $j_{o}(\boldsymbol{u}^{\ast })$ and $j_{t}(\boldsymbol{u}^{\ast })$ directly. Figure 19 presents the variation in these contributions at large $\unicode[STIX]{x1D712}$ , evaluated numerically at $B=10^{4}$ using problem [R]. The inner boundary layer contribution $j_{i}(\boldsymbol{u}^{\ast })\sim \unicode[STIX]{x1D712}^{0.5}$ is clearly dominant. The outer boundary layer contribution $j_{o}(\boldsymbol{u}^{\ast })$ is approximately constant and the dissipation in the tip region decays as $j_{t}(\boldsymbol{u}^{\ast })\sim \unicode[STIX]{x1D712}^{-0.3}$ .

Figure 19.

Contributions to the plastic dissipation at large $\unicode[STIX]{x1D712}$ , computed for the rectangular particle at $B=10^{4}$ using problem [R]: (a) $j_{i}(\boldsymbol{u}^{\ast })$ ; (b) $j_{o}(\boldsymbol{u}^{\ast })$ ; (c) $j_{t}(\boldsymbol{u}^{\ast })$ (symbols from computation and solid lines fitted, with the indicated slopes).

We now consider fixed large $\unicode[STIX]{x1D712}$ and variations with $(1-Y/Y_{c})$ (equivalently $B$ ), starting with the inner boundary layer. The velocity $\boldsymbol{u}^{\ast }$ decays from ${\sim}1$ to ${\sim}O(1-Y/Y_{c})$ across the width $\unicode[STIX]{x1D6FF}_{i}$ and has a length ${\sim}\ell _{||}$ along the particle. It follows that the contribution to $j(\boldsymbol{u}^{\ast })$ is ${\sim}\ell _{||}$ . This must be compatible with problem [M]. On rescaling of $\boldsymbol{u}^{\ast }$ we find

(4.7) $$\begin{eqnarray}\unicode[STIX]{x03C0}\frac{b(\unicode[STIX]{x1D712})}{Y_{c}(\unicode[STIX]{x1D712})}(1-Y/Y_{c})^{\unicode[STIX]{x1D708}}\sim j(\boldsymbol{u})=Uj(\boldsymbol{u}^{\ast })\sim b(\unicode[STIX]{x1D712})\ell _{||}(1-Y/Y_{c})^{\unicode[STIX]{x1D708}},\end{eqnarray}$$

suggesting that $Y_{c}(\unicode[STIX]{x1D712})\sim 1/\ell _{||}\sim \unicode[STIX]{x1D712}^{-1/2}$ . We might explicitly take the leading-order term as corresponding to the drop in velocity along the full length of the particle,

(4.8) $$\begin{eqnarray}j(\boldsymbol{u}^{\ast })\sim j_{i}(\boldsymbol{u}^{\ast })\sim 2\ell _{||}+O(\unicode[STIX]{x1D712}^{-1/2})+O(1-Y/Y_{c}),\end{eqnarray}$$

accounting for geometric effects and small non-zero velocity at the edge of the boundary layer. This results explicitly in $Y_{c}(\unicode[STIX]{x1D712})\sim \unicode[STIX]{x03C0}/(2\ell _{||})$ at leading order.

This asymptotic behaviour is entirely compatible with our earlier results. The critical plastic drag coefficient is

(4.9) $$\begin{eqnarray}C_{D,c}^{p}=\frac{\unicode[STIX]{x03C0}}{\ell _{\bot }Y_{c}}=\lim _{B\rightarrow \infty }\frac{j(\boldsymbol{u}^{\ast })}{\ell _{\bot }}.\end{eqnarray}$$

In every example computed (and those estimated from the lower-bound characteristics solutions) we have $C_{D,c}^{p}\sim 2\unicode[STIX]{x1D712}$ at large $\unicode[STIX]{x1D712}$ . This coincides precisely with the leading-order contribution $j(\boldsymbol{u}^{\ast })\sim 2\ell _{||}$ and with $Y_{c}(\unicode[STIX]{x1D712})\sim \unicode[STIX]{x03C0}/(2\ell _{||})$ .

The boundary layer contribution to the viscous dissipation $a(\boldsymbol{u}^{\ast },\boldsymbol{u}^{\ast })$ is ${\sim}\ell _{||}/\unicode[STIX]{x1D6FF}_{i}$ , which translates to

(4.10) $$\begin{eqnarray}a(\boldsymbol{u},\boldsymbol{u})=U^{2}a(\boldsymbol{u}^{\ast },\boldsymbol{u}^{\ast })\sim b^{2}(\unicode[STIX]{x1D712})\ell _{||}\frac{(1-Y/Y_{c})^{2\unicode[STIX]{x1D708}}}{\unicode[STIX]{x1D6FF}_{i}}.\end{eqnarray}$$

Our analysis of (4.6) suggests that $\unicode[STIX]{x1D6FF}_{i}\sim (1-Y/Y_{c})^{k}$ , as $\unicode[STIX]{x1D712}\rightarrow \infty$ . The bound (2.20) implies $k\leqslant 2(\unicode[STIX]{x1D708}-1)$ , whereas (2.21) is more restrictive, implying $k\leqslant \unicode[STIX]{x1D708}-1$ (note also $\unicode[STIX]{x1D708}\geqslant 1$ ). For the narrowest possible boundary layer $k=\unicode[STIX]{x1D708}-1$ and for the widest possible $k=0$ . However, numerical evidence is that $k>0$ . In terms of $B$ we have $\unicode[STIX]{x1D6FF}_{i}\sim (Y_{c}/b)^{k/\unicode[STIX]{x1D708}}B^{-k/\unicode[STIX]{x1D708}}$ . As we have determined that $\unicode[STIX]{x1D6FF}_{i}$ is independent of $\unicode[STIX]{x1D712}$ as $\unicode[STIX]{x1D712}\rightarrow \infty$ we deduce that

(4.11) $$\begin{eqnarray}b(\unicode[STIX]{x1D712})\sim Y_{c}(\unicode[STIX]{x1D712})\sim \unicode[STIX]{x1D712}^{-1/2},\quad \text{as}~\unicode[STIX]{x1D712}\rightarrow \infty .\end{eqnarray}$$

We may then deduce that

(4.12a ) $$\begin{eqnarray}\displaystyle & \displaystyle j(\boldsymbol{u})\sim \frac{\ell _{||}Y_{c}}{B}\sim B^{-1}, & \displaystyle\end{eqnarray}$$

(4.12b ) $$\begin{eqnarray}\displaystyle & \displaystyle a(\boldsymbol{u},\boldsymbol{u})\sim \frac{\ell _{||}Y_{c}^{1/\unicode[STIX]{x1D708}}}{b^{1/\unicode[STIX]{x1D708}-2}}B^{-2+k/\unicode[STIX]{x1D708}}\sim \unicode[STIX]{x1D712}^{-1/2}B^{-2+k/\unicode[STIX]{x1D708}}, & \displaystyle\end{eqnarray}$$

(4.12c ) $$\begin{eqnarray}\displaystyle & \displaystyle B\sim (1-Y/Y_{c})^{-\unicode[STIX]{x1D708}}. & \displaystyle\end{eqnarray}$$

Figure 20 shows the values $\unicode[STIX]{x1D708}$ calculated by fitting our computations, all at $\unicode[STIX]{x1D712}=100$ , for the three geometries. We observe different convergence rates $\unicode[STIX]{x1D708}$ , with the smaller exponents coming from the rectangle.

So far this is non-conclusive. We may try to estimate the exponents $\unicode[STIX]{x1D708}$ and $k$ directly from physical grounds or simpler assumptions. A simple balance of the buoyancy force with a viscous shear stress in the boundary layer suggests $\unicode[STIX]{x1D6FF}_{i}\sim \sqrt{\ell _{||}Y_{c}}B^{-1/2}$ , i.e.  $k=\unicode[STIX]{x1D708}/2$ . For the flat-plate limit we might follow the boundary layer analysis of Oldroyd (Reference Oldroyd1947), who derived $k/\unicode[STIX]{x1D708}=1/3$ , meaning here $\unicode[STIX]{x1D708}=3/2$ for the narrowest boundary layer. Recently, however, Balmforth et al. (Reference Balmforth, Craster, Hewitt, Hormozi and Maleki2016) have shown that Oldroyd’s analysis and scaling are erroneous for the flat plate, but does apply to thin free shear layers. Following their analysis, the boundary layer close to a flat plate has $k/\unicode[STIX]{x1D708}=1/2$ , to which our results should converge as $\unicode[STIX]{x1D712}\rightarrow \infty$ . It may be, however, that for $\unicode[STIX]{x1D712}=100$ we are still far from the flat-plate limit $\unicode[STIX]{x1D712}\rightarrow \infty$ .

Figure 20.

Convergence at large $B$ for the three geometries at $\unicode[STIX]{x1D712}=100$ : symbols are computed and lines are power-law fits with exponent $m$ (the slope of the fitted lines in log–log). This gives $\unicode[STIX]{x1D708}\approx 1.61$ , 1.32 and 1.60 for ellipse, rectangle and diamond, respectively.

Let us now examine the next largest contribution to $Y_{c}$ , i.e.  $j_{o}(\boldsymbol{u}^{\ast })$ . As well as the $\unicode[STIX]{x1D712}$ dependence considered earlier, both $\unicode[STIX]{x1D714}_{p}$ and $\ell _{p}$ may vary with $B$ , say

(4.13a,b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D714}_{p}\propto B^{-q},\quad \ell _{p}\propto B^{s}. & & \displaystyle\end{eqnarray}$$

To maintain the leading-order balance in (4.6) we require $q=2s+k/\unicode[STIX]{x1D708}$ . We may now estimate the contribution of the outer boundary layer in problem [R] to the plastic dissipation, say $j_{o}(\boldsymbol{u}^{\ast })$ . Across the outer boundary layer the velocity drops from ${\sim}\unicode[STIX]{x1D714}_{p}\ell _{p}$ to zero, along the length ${\sim}\ell _{||}$ , i.e.

(4.14) $$\begin{eqnarray}j_{o}(\boldsymbol{u}^{\ast })\sim \unicode[STIX]{x1D714}_{p}\ell _{p}\ell _{||}\sim \unicode[STIX]{x1D712}^{0}B^{-s-k/\unicode[STIX]{x1D708}}.\end{eqnarray}$$

In terms of the limit $B\rightarrow \infty$ the dependence of $\ell _{p}$ observed computationally appears very weak, so that we might infer $s\approx 0$ , i.e. the extent of the yielded region in the limit of no flow depends on $\unicode[STIX]{x1D712}$ but not significantly on $B$ . We conclude that $j_{o}(\boldsymbol{u}^{\ast })/j(\boldsymbol{u}^{\ast })\rightarrow 0$ , both as $B\rightarrow \infty$ and as $\unicode[STIX]{x1D712}\rightarrow \infty$ .

Finally, it is interesting to reflect physically on why for large $\unicode[STIX]{x1D712}$ the viscoplastic model and perfectly plastic model differ, in terms of outer yield surface and velocity. In the absence of viscosity in the rigid–plastic model, the plate does not drag the material as it is moving and hence no rigid rotation needs to be formed. This may also explain why in the viscoplastic solution we get a significantly smaller plastic drag coefficient for the rectangle than the upper bound based on the proposed velocity fields in § 3.3.2. The dissipation in the tip region is in fact smaller in the plastic solution than that in the viscoplastic solution (decreasing like $\unicode[STIX]{x1D712}^{-1/2}$ compared to $\unicode[STIX]{x1D712}^{-0.3}$ ). However, the upper-bound plastic solution has leading-order term ${\sim}(1+\sqrt{2})(\sqrt{\unicode[STIX]{x03C0}}/2)\unicode[STIX]{x1D712}^{1/2}$ , whereas the viscoplastic solution has leading-order contribution $(\sqrt{\unicode[STIX]{x03C0}}/2)\unicode[STIX]{x1D712}^{1/2}$ .

4.3.2 Diamond

The flow around the long diamond ( $\unicode[STIX]{x1D712}\gg 1$ ) shares many characteristics with the small- $\unicode[STIX]{x1D712}$ flow of § 4.1, but is simpler in having only three regions contributing to $j(\boldsymbol{u}^{\ast })$ . Referring to one quarter of the flow domain, these three regions are as follows.

1. (i) A rotating zone that is centred at the ‘corner’ of the diamond on the $x$ -axis. The fluid rotates on circular streamlines about this point and appears to have approximately constant speed in the azimuthal direction. The radius of this region is ${\sim}\ell _{||}/(2\cos \unicode[STIX]{x1D6FC})$ and the azimuthal speed is ${\sim}\text{sin}\unicode[STIX]{x1D6FC}$ . Here we have $\tan \unicode[STIX]{x1D6FC}=\unicode[STIX]{x1D712}^{-1}$ , i.e.  $\unicode[STIX]{x1D6FC}<\unicode[STIX]{x03C0}/4$ is the angle that the top of the diamond makes with the $y$ -axis.

2. (ii) A thin shear layer between the particle and the rotating zone. The length of this boundary layer is ${\sim}\ell _{||}/(2\cos \unicode[STIX]{x1D6FC})$ . Across this layer the normal component of velocity is continuous and the tangential velocity drops from $\cos \unicode[STIX]{x1D6FC}$ to $0$ .

3. (iii) An outer boundary layer, surrounding the rotating zone, across which the azimuthal velocity drops from ${\sim}\text{sin}\unicode[STIX]{x1D6FC}$ to $0$ . The circumferential length of this boundary layer is $\ell _{||}(\unicode[STIX]{x03C0}/2+\unicode[STIX]{x1D6FC})/(2\cos \unicode[STIX]{x1D6FC})$ .

We can estimate the contributions to $j(\boldsymbol{u}^{\ast })$ from each of the above regions as follows: (i) ${\sim}2\ell _{||}(\unicode[STIX]{x03C0}/2+\unicode[STIX]{x1D6FC})\tan \unicode[STIX]{x1D6FC}$ ; (ii)  ${\sim}2\ell _{||}$ ; (iii)  ${\sim}2\ell _{||}(\unicode[STIX]{x03C0}/2+\unicode[STIX]{x1D6FC})\tan \unicode[STIX]{x1D6FC}$ . Combining these we find

(4.15) $$\begin{eqnarray}C_{D,c}^{p}=\lim _{B\rightarrow \infty }[j(\boldsymbol{u}^{\ast })/\ell _{\bot }]\sim \frac{2\ell _{||}[1+(\unicode[STIX]{x03C0}+2\unicode[STIX]{x1D6FC})\tan \unicode[STIX]{x1D6FC}]}{\ell _{\bot }}\sim 2\unicode[STIX]{x1D712}+4(\unicode[STIX]{x03C0}-\tan ^{-1}\unicode[STIX]{x1D712}),\end{eqnarray}$$

as $\unicode[STIX]{x1D712}\rightarrow \infty$ . Note that the leading-order contribution has come from the boundary layer (ii) as $\unicode[STIX]{x1D712}\rightarrow \infty$ . This expression agrees with that of the characteristics analysis.