1. Introduction
The steady flow of a viscous fluid in a planar wedge is a fundamental problem in fluid mechanics (Drazin & Riley Reference Drazin and Riley2006). First documented by Jeffery (Reference Jeffery1915) and Hamel (Reference Hamel1917), the flow is a rare example of a general solution to the Navier–Stokes flow at all Reynolds numbers. ‘Jeffery–Hamel flow’ has been studied extensively (see, for e.g. Dean Reference Dean1934; Rosenhead Reference Rosenhead1940; Fraenkel Reference Fraenkel1962; Banks, Drazin & Zaturska Reference Banks, Drazin and Zaturska1988). For convex wedges, with half-angle 
$\alpha \leqslant 90^{\circ }$, converging and diverging solutions exist at vanishing Reynolds numbers for all wedges, while at larger Reynolds numbers the flow forms thin boundary layers at the walls. These layers are stable for converging flow but unstable for diverging flow, resulting, above a critical Reynolds number, in separation, and the formation of a central jet and regions of reversed flow at the walls. In essence this occurs because the weak viscous stresses are unable to balance the adverse pressure gradient of the diverging flow.
Steady converging flows of plastic materials also have a long history. Nadai (Reference Nadai1924) derived a solution for the converging flow of a perfectly rigid plastic in a wedge, with a constant friction imposed on the walls (solution presented in English by Hill Reference Hill1950), and Shield (Reference Shield1955) solved the corresponding problem in an axisymmetric conical geometry. Notably, both of these flow solutions exhibit slip and divergent shear rates at the boundaries. Motivated by industrial applications such as wire drawing, polymer processing, extrusion of pastes and the discharge of granular materials from hoppers and silos, more recent studies have examined the converging flows of other non-Newtonian fluids. For example, the slow converging flow of incompressible power-law fluids has been studied by Durban (Reference Durban1986), Brewster et al. (Reference Brewster, Chapman, Fitt and Please1995) and Nagler (Reference Nagler2017), while the viscoelastic problem has been studied by Hull & Pearson (Reference Hull and Pearson1984). The work by Brewster et al. (Reference Brewster, Chapman, Fitt and Please1995) is of particular interest, as they consider a power-law fluid of small shear exponent, for which they construct similarity solutions under the assumption of purely radial flow, and demonstrate the need for a complex boundary layer structure in order to enforce no slip at the walls.
Viscoplastic fluids are a particular class of non-Newtonian fluid which act as a rigid plastic or flow as a viscous fluid, depending on whether the stress is less than or exceeds a critical yield stress, respectively. This behaviour is common for slurries and suspensions, and the viscoplastic model has wide ranging applications in geophysics and industry (Ancey Reference Ancey2007; Balmforth, Frigaard & Ovarlez Reference Balmforth, Frigaard and Ovarlez2014; Frigaard, Paso & de Souza Mendes Reference Frigaard, Paso and de Souza Mendes2017). The non-Newtonian behaviour is characterised by the dimensionless Bingham number, 
${\textit {Bi}}$, which measures the magnitude of the yield stress to the viscous stresses developed by the flow. (The precise definition of 
${\textit {Bi}}$ for the flow analysed in this study is given in (2.13).)
For the slow converging flow of a viscoplastic fluid, Cristescu (Reference Cristescu1975) studied the stress distribution produced by a given kinematically feasible radial velocity profile in a cone, without attempting to solve the full system of partial differential equations required to explicitly determine the coupled velocity and stress fields. Motivated by strip drawing in the regime of low yield stress or fast drawing speed, Sandru & Camenshi (Reference Sandru and Camenshi1979) analysed the converging flow of a viscoplastic material when the Bingham number was relatively small. They imposed ‘friction’ boundary conditions in the form of an imposed ratio between shear and normal stresses at the boundaries. No-slip boundary conditions can be treated as the limiting case of this approach, in which the deviatoric normal stress vanishes and the shear stress is maximised. Although this case was not the focus of Sandru & Camenshi (Reference Sandru and Camenshi1979) and they stop short of explicit evaluation of the velocity fields, we demonstrate that their perturbation solution may be used to calculate the velocity fields and, in particular, emphasise that the solution implies purely radial flow is not possible. Durban (Reference Durban1984, Reference Durban1985, Reference Durban1986) calculated flow profiles and stress fields for both a planar wedge and an axisymmetric conical die, under constant friction boundary conditions, as opposed to no-slip boundary conditions. These results are restricted by the assumption that the wall friction is low compared to the yield stress, as relevant to well lubricated, metal-forming processes at high temperatures, and the flow profiles are assumed to be purely radial and approximately uniform. Later, Alexandrov, Lyamina & Mishuris (Reference Alexandrov, Lyamina and Mishuris2007) studied viscoplastic flow in a wedge under maximum friction boundary conditions, their definition of which reduces to the condition that the deviatoric stress consists only of shear stresses at the wall. Under this boundary condition they find that the radial velocity must be constant at the boundary, which is inconsistent with a slipping regime in which the radial velocity increases towards the vertex of the wedge. They further consider constitutive laws with a saturation stress (where the stress remains bounded as strain rates diverge) and find that the classical assumptions of a radially independent stress state and purely radial flow are inconsistent when shear stresses become comparable to the shear yield stress. They identify the cause of these difficulties as the divergent strain rates that occur at the wall in the perfectly plastic solution. However, they do not attempt to derive velocity or stress fields for either the standard viscoplastic model, or for one with a saturation stress. Finally, Ara et al. (Reference Ara, Khan, Sultan and Ullah2019) recently performed a numerical study of wedge flow of a Bingham fluid including heat transfer, in the presence of a magnetic field. Again, this study assumed a purely radial flow to derive ordinary differential equations for similarity solutions. These solutions are found to be dependent on radius through the authors’ definition of Bingham number, and so the conclusions of this work rely on an implicit assumption that this radial dependence is sufficiently weak to be neglected in the expression of conservation of mass.
In this paper, we present a comprehensive analysis of converging flow of viscoplastic fluid through a wedge or cone, that satisfies no-slip boundary conditions, in both the plastically 
$({\textit {Bi}} \gg 1)$ and viscously 
$({\textit {Bi}}\ll 1)$ dominated regimes. We initially carry out a full analysis for a Bingham fluid in a planar wedge, before extending the theory to a Herschel–Bulkley fluid in a wedge and a Bingham fluid in an axisymmetric conical geometry, both of which share the same analytical framework to the flow of a Bingham fluid through a wedge, but are algebraically more involved. Most importantly, in all cases we find that the viscoplasticity introduces a weak angular flow towards the centre of the wedge or cone, explaining the challenges faced by previous attempts to find purely radial solutions. Furthermore, in the plastically dominated regime, boundary layers are found at the wall, allowing the solution to satisfy no slip. The angular extent of these boundary layers is found to decrease with Bingham number as 
${\textit {Bi}}^{-1/2}$ (or 
${\textit {Bi}}^{-1/(N+1)}$ for a Herschel–Bulkley fluid of flow index 
$N$) and depends also on the radial distance from the apparent apex of the wedge or cone, 
$r$. Direct numerical simulations for a Bingham fluid in a planar wedge are carried out using the finite element method, verifying the analytical results in both regimes.
In the plastic regime, the converging flow field away from the boundaries is given to leading order by the solutions derived by Nadai (Reference Nadai1924) and Shield (Reference Shield1955), with viscous effects becoming significant within relatively thin boundary layers in which the shear rates are relatively large and the shear stresses exceed the yield stress. The velocity field in this regime may be determined using viscoplastic boundary layer techniques, originally developed by Oldroyd (Reference Oldroyd1947). Oldroyd explored a distinguished limit in which viscous shear stresses and plastic extensional stresses are both significant in the streamwise momentum balance, which occurs for dimensionless boundary layer widths of order 
${\textit {Bi}}^{-1/3}$. Viscoplastic boundary layers were re-examined by Piau (Reference Piau2002), who argued for a different distinguished scaling from Oldroyd, in which viscous shear stresses dominate the momentum balance, and boundary layer widths instead scale like 
${\textit {Bi}}^{-1/2}$. Most recently, Balmforth et al. (Reference Balmforth, Craster, Hewitt, Hormozi and Maleki2017) developed a consistent boundary layer theory for viscoplastic flows, clarifying the conditions under which Oldroyd's and Piau's scalings are appropriate. Specifically, they determined that where the boundary layer is adjacent to a rigid wall, the boundary layer structure is consistent with Piau's arguments (although the exact scaling can vary depending on the far field pressure gradient), whereas when bounded only by unyielded or plastically deforming fluid, the appropriate scaling is Oldroyd's. In the problems studied by Balmforth et al. (Reference Balmforth, Craster, Hewitt, Hormozi and Maleki2017), the boundary layer solutions are only explicitly asymptotically matched to rigid plugs. However, in the problem considered in this work, the velocity is varying at the outer edge of the boundary layer due to the converging geometry, and thus our asymptotic construction is rather different. It will be shown that the solution does share some similarities with other viscoplastic flows which exhibit regions of nearly plastic deformation, such as flow in an eccentric annulus (Walton & Bittleston Reference Walton and Bittleston1991), flow around a cylinder (Tokpavi, Magnin & Jay Reference Tokpavi, Magnin and Jay2008) or thin-layer gravitational flow (Balmforth & Craster Reference Balmforth and Craster1999). In the weakly yielded regions, or ‘pseudo-plugs’, present in these flows, the second invariant of the deviatoric stress does not exceed the yield stress at leading order, but does so at higher orders. These regions are typically delineated from fully yielded material by ‘fake yield surfaces’. In accord with the previous results for viscoplastic boundary layers, we show that the width of the layer scales like 
${\textit {Bi}}^{-1/2}$, as expected for flow bounded by rigid walls, but that it is matched, via an intermediate asymptotic layer encompassing a fake yield surface, to a weakly yielded, pseudo-rigid plastic solution in the bulk of the wedge. The existence of these two asymptotic layers is a consequence of the need to both impose no slip and to regularise the divergent shear rates found at the wall in the rigid plastic solution of Nadai (Reference Nadai1924).
We first define the problem and outline the solution for the converging wedge flow of a Bingham fluid (§ 2). We then construct the asymptotic solution for the plastic regime and compare the results to direct numerical simulations (§ 3). The asymptotic solution for the viscous regime is derived and compared to direct simulations in § 4. We focus initially on the Bingham fluid and planar geometry, since this case most clearly demonstrates the asymptotic structure with minimal algebraic complications. We then generalise the problem to a Herschel–Bulkley fluid (§ 5), and to an axisymmetric conical geometry (§ 6), which both follow the same asymptotic structure as the solutions in §§ 3,4, but are algebraically more involved due to rheology (§ 5) and geometry (§ 6). We summarise and conclude our study in § 7. There are also four appendices in which we provide some of the algebraic details required for the derivation of the flow solutions.
2. Problem definition: Bingham fluid in a planar wedge
We analyse steady slow incompressible flow of a viscoplastic fluid through a planar wedge of half-angle 
$\alpha$ (see figure 1). We restrict our focus to 
$\alpha \leqslant {\rm \pi}/2$, for which the geometry represents a converging flow in a wedge, as opposed to a sink flow external to a wedge when 
$\alpha >{\rm \pi} /2$. A volume flux per unit width, 
$Q$ (taken to be positive for flow towards the vertex), is imposed and we assume that the wedge is sufficiently long that end effects may be ignored. The flow field is predominantly in the radial direction and is assumed symmetric about the centre-line of the wedge. Plane polar coordinates, 
$(r,\theta )$, are defined relative to the apparent vertex of the wedge. The key dependent variables are the velocities in the radial and angular directions, denoted by 
$u$ and 
$v$ respectively.
Figure 1. A schematic of the flow geometry.
We initially assume the constitutive law for a Bingham fluid, before generalising to the Herschel–Bulkley model in § 5. Thus, the constitutive law is given by
\begin{equation} \begin{cases} \begin{pmatrix} \tau_{rr} \\ \tau_{r\theta} \end{pmatrix} = \left(\mu+\dfrac{\tau_c}{\dot{\gamma}}\right) \begin{pmatrix} \dot{\gamma}_{rr}\\ \dot{\gamma}_{r\theta} \end{pmatrix} & \text{if } \tau> \tau_c,\\ \dot{\gamma}=0 & \text{otherwise}, \end{cases}\end{equation}
where 
$\mu$ is the Bingham plastic viscosity and 
$\tau _{c}$ is the yield stress. The variables 
$(\tau _{rr},\tau _{r\theta })$ and 
$(\dot {\gamma }_{rr} ,\dot {\gamma }_{r\theta })$ are the components of the deviatoric stress and strain-rate tensors, respectively, 
$\tau$ and 
$\dot {\gamma }$ are the corresponding second invariants given by
\begin{equation} \dot{\gamma}^2=\dot{\gamma}_{rr}^2+\dot{\gamma}_{r\theta}^2 \quad \text{and} \quad \tau^2=\tau_{rr}^2+\tau_{r\theta}^2, \end{equation}and the strain-rate tensor is defined by
\begin{equation} \dot{\boldsymbol{\gamma}}=\left(\boldsymbol{\nabla}\boldsymbol{u}\right)+\left(\boldsymbol{\nabla}\boldsymbol{u}\right)^\textrm{T}, \end{equation}
where 
$\boldsymbol{u}$ is the velocity vector and T is the transpose. Due to the flow convergence we expect that 
$\dot {\gamma }>0$, and hence that the fluid is yielding everywhere and rigid plugs do not form in the flow. As such, we need only consider the 
$\tau >\tau _c$ case of the constitutive law (2.1).
We have the following dimensional quantities in our problem: density, 
$\rho$, viscosity, 
$\mu$, yield-stress, 
$\tau _c$, volume flux per unit width, 
$Q$, and typical radial distance from the vertex, 
$R$. We scale velocities by 
$Q/R$, strain rates by 
$Q/R^2$, stresses and pressures by 
$\mu Q/R^2$ and, unless stated otherwise, all variables will be assumed dimensionless. On the further assumption that the flow is sufficiently slow to ignore inertial effects, the non-dimensional governing equations in polar coordinates are given by
\begin{gather} \frac{1}{r}\frac{\partial \left(ru\right)}{\partial r} + \frac{1}{r}\frac{\partial v}{\partial \theta} = 0, \end{gather}
\begin{gather}\frac{\partial p}{\partial r} = \frac{\partial \tau_{rr}}{\partial r}+\frac{1}{r}\frac{\partial \tau_{r\theta}}{\partial \theta}+\frac{2}{r}\tau_{rr}, \end{gather}
\begin{gather}\frac{1}{r} \frac{\partial p}{\partial \theta} = \frac{\partial \tau_{r\theta}}{\partial r}-\frac{1}{r}\frac{\partial \tau_{rr}}{\partial \theta}+\frac{2}{r}\tau_{r\theta}, \end{gather}which express incompressibility (2.4) and the balance of momentum in the radial (2.5) and angular (2.6) directions. The deviatoric stress and strain-rate tensors are given by
\begin{gather} \begin{pmatrix} \tau_{rr} \\ \tau_{r\theta} \end{pmatrix} = \left(1+\frac{{\textit{Bi}}}{\dot{\gamma}}\right) \begin{pmatrix} \dot{\gamma}_{rr}\\ \dot{\gamma}_{r\theta} \end{pmatrix}, \end{gather}
\begin{gather} \dot{\gamma}_{rr}= 2\frac{\partial u}{\partial r}, \quad \dot{\gamma}_{r\theta}=\frac{1}{r}\frac{\partial u}{\partial\theta}+\frac{\partial v}{\partial r}-\frac{v}{r}, \end{gather}
\begin{gather} \dot{\gamma}=\sqrt{\left(\frac{1}{r}\frac{\partial u}{\partial\theta}+\frac{\partial v}{\partial r}-\frac{v}{r}\right)^2+4\left(\frac{\partial u}{\partial r}\right)^2}. \end{gather}The boundary conditions are
\begin{gather} \left(u,v\right)=\left(0,0\right)\quad \text{at } \theta={\pm}\alpha, \end{gather}
\begin{gather}v=\frac{\partial u}{\partial \theta}=0\quad \text{at } \theta=0, \end{gather}which express no slip at the edges of the wedge (2.10) and symmetry about the mid-line (2.11). There is, additionally, an integral expression for the volume flux per unit width given by
\begin{equation} \int_{-\alpha}^{\alpha}ru\,\textrm{d}\theta ={-}1. \end{equation}The residual dimensionless parameter is given by
\begin{equation} {\textit{Bi}}=\frac{R^2\tau_c}{\mu Q}, \end{equation}
which defines a Bingham number, giving the ratio of plastic to viscous stresses. Note that this Bingham number depends on the scale of the radial distance from the vertex, due to the geometry imposing no other length scale. For given rheological parameters and volume flux, the parameters define a critical dimensional radial distance, 
$r_c=(\mu Q/\tau _c)^{-1/2}$, below which viscous forces dominate, and beyond which plastic forces dominate the flow. We will concern ourselves with the 
${\textit {Bi}}\gg 1$ (plastically dominated) and 
${\textit {Bi}}\ll 1$ (viscously dominated) regimes, which apply sufficiently far from and sufficiently close to the vertex, respectively. We note that we could have chosen to scale radial distances by 
$r_c$, in which case all dimensional parameters are removed from the governing equation. However, our exposition is clearer and easier to compare with previous studies when based upon the dimensionless variables and governing equations given above. We note that, due to the quadratic dependence of 
${\textit {Bi}}$ on 
$R$, we anticipate that our asymptotic expansions will be functions of 
$r^2{\textit {Bi}}$. The absence of inertial terms in (2.4)–(2.6) requires that the Reynolds number, 
$Re=\rho Q/\mu$, is sufficiently small. In practice this requirement varies depending on the regime of interest. For the 
${\textit {Bi}}\gg 1$ regime it will be sufficient that 
$Re=O(1)$, while for the 
${\textit {Bi}}\ll 1$ regime we will require 
$Re\ll {\textit {Bi}}$.
2.1. An outline of the solution
Before performing a detailed analysis of the problem, we outline the key results. In this discussion we assume that the wedge extends from a large dimensional radius, 
$r\gg r_c$, to a small radius 
$r\ll r_c$. Thus, in the far field of the wedge, and away from the walls, the leading-order flow is given by Nadai's solution for converging plastic flow (Nadai Reference Nadai1924; Hill Reference Hill1950),
\begin{equation} u={-}\frac{A}{r\left(c-\cos2\psi\right)}, \quad v=0, \quad \frac{\textrm{d}\psi}{\textrm{d}\theta}=c\sec2\psi-1. \end{equation}
Here, 
$A$ and 
$c$ are constants determined by the flux and boundary conditions, respectively, which are given explicitly later in (3.14) and (3.15). The function 
$\psi$ measures the angle between the directions of the principle stresses and the coordinate axes, and is a function of the polar angle, such that 
$\tau _{r\theta }=\tau _{rr}\tan 2\psi$. Viscous shear-stress-dominated flow will be shown to be confined to a thin boundary layer of angular thickness 
$O((r^2{\textit {Bi}})^{-1/2})$, allowing the solution to satisfy no slip. The variation of the angular extent of the boundary layer with 
$r$ is such that the boundary layer is of a constant Cartesian thickness along the no-slip boundary. The first-order corrections to the plastic flow (2.14a–c) will be shown to also be 
$O((r^2{\textit {Bi}})^{-1/2})$ and are driven by the flow within the boundary layer. A thin intermediate layer is also required between the bulk of the wedge and the boundary layer, for asymptotic matching of the respective solutions in these regions. One interpretation of this boundary layer structure is that the plastically dominated bulk solution, (2.14a–c), suffers from two different inconsistencies at the wall, namely diverging shear rates and slip, which require addressing at different scales. The intermediate layer plays the role of regularising the divergent shear rates, while the boundary layer enforces the no-slip boundary condition.
In contrast, at the outflow of the wedge, the shear rate is much higher due to radial convergence, and viscous forces dominate the flow. The leading-order flow is the classical Stokes solution,
\begin{align} u={-}\frac{\cos(2\theta)-\cos(2\alpha)}{r\left(\sin(2\alpha)-2\alpha\cos(2\alpha)\right)}, \quad v=0, \quad p={-}\frac{2\cos2\theta}{r^2\left(\sin2\alpha-2\alpha\cos2\alpha\right)}+\text{const.}, \end{align}
which defines a purely converging flow for all 
$\alpha \leqslant {\rm \pi}/2$. Here, the first-order corrections due to plasticity are 
$O(r^2{\textit {Bi}})$. Comparison of the radial velocity profiles in the plastic and viscous limits (figure 2), shows that the viscous profile has enhanced flow at the centre of the wedge compared to the plastic flow. A consequence of this is that we expect an angular velocity to be induced towards the centre of the wedge, away from the walls, to satisfy conservation of mass as the fluid transitions from the plastic profile at the inflow to the viscous profile at the outflow. This is indeed borne out and quantified in the analytical solution constructed below.
3. Plastic regime 
$({\textit {Bi}}\gg 1)$
In the regime of large Bingham number, 
${\textit {Bi}}\gg 1$, or equivalently at large distances from the apparent apex of the wedge, the constitutive equation (2.7) has three main regimes for different sizes of the strain rate, 
$\dot {\gamma }$: plastic stresses dominate the momentum balance when the strain rate is 
$O(1)$; viscous stresses dominate in a thin boundary layer where the strain rate is high; and both become significant for an intermediate magnitude of the strain rate.
3.1. The bulk solution
In the bulk of the wedge, away from the boundaries, we expect that 
$\dot {\gamma }=O(1)$. In this case the constitutive equation is dominated by the yield stress and we find that the deviatoric stresses are determined up to 
$O(1)$ by the plastic flow equation
\begin{equation} \tau_{ij}={\textit{Bi}}\frac{\dot{\gamma}_{ij}}{\dot{\gamma}}, \end{equation}which implies that, at leading order, the magnitude of the deviatoric stress is everywhere equal to the yield stress, and that the deviatoric stress aligns with the rate of strain. As for the flow in an eccentric annulus (Walton & Bittleston Reference Walton and Bittleston1991) and gravitationally driven shallow layer flow (Balmforth & Craster Reference Balmforth and Craster1999), while the stress does not exceed the yield stress at leading order, the material is nonetheless undergoing deformation, and so must be yielded. This region is therefore analogous to ‘pseudo-plugs’ found in other flow scenarios (Walton & Bittleston Reference Walton and Bittleston1991; Balmforth & Craster Reference Balmforth and Craster1999).
Following the approach of Nadai (Reference Nadai1924) in the bulk of the flow we introduce a variable, 
$\psi$, such that
\begin{equation} \begin{pmatrix} \tau_{rr} \\ \tau_{r\theta} \end{pmatrix} = {\textit{Bi}} \begin{pmatrix} \cos{2\psi}\\ \sin{2\psi} \end{pmatrix}. \end{equation}
We must have 
$\psi (r,0)=0$ by symmetry, and in the analysis that follows will focus on the region 
$0\leqslant \theta \leqslant \alpha$, with symmetry allowing the automatic construction of the flow in 
$-\alpha \leqslant \theta \leqslant 0$. In the solution of Nadai (Reference Nadai1924), (2.14a–c), the velocity does not vanish on the wall. Hence, to satisfy the no-slip boundary condition, we expect the shear rate, 
$\partial u/\partial \theta$, to become large in a small region near the walls. This leads to a different regime where the viscous stresses are no longer negligible in the balance of momentum. We define the angle, 
$\theta _Y$, at which this change of dynamical balance occurs by
\begin{equation} \theta_Y=\alpha-\epsilon\varPhi(r), \end{equation}
where 
$\epsilon$ is a small parameter depending on 
${\textit {Bi}}$ in a way that is yet to be determined, and we have introduced a potential dependence on 
$r$ through the function 
$\varPhi$. The notation 
$\theta _Y$ alludes to the fact that this is a fake yield surface, delineating weakly yielded fluid from highly sheared, fully yielded fluid adjacent to the boundary.
The governing equation (3.2) can be rearranged to obtain
\begin{equation} \tan\left(2\psi\right)=\frac{\tau_{r\theta}}{\tau_{rr}}=\frac{\dot{\gamma}_{r\theta}}{\dot{\gamma}_{rr}}.\end{equation}
At the wall, 
$\dot {\gamma }_{rr}$ vanishes, and in the thin boundary layer, 
$\dot {\gamma }_{r\theta }$ is large. Hence, both ratios on the right-hand side of (3.4) are large in the boundary layer, and we require 
$\psi \to {\rm \pi}/4$ as 
$\theta \to \theta _Y$ in order for the bulk solution to match to the boundary layer. This dependence of 
$\psi$ on 
$r$ through the boundary condition at 
$\theta =\pm \theta _Y(r)$, leads to apparent difficulties in the asymptotic investigation of the boundary layer due to divergent strain rates which obfuscate the straightforward expansion of the flow variables. To avoid this difficulty, we use strained coordinates, introducing a coordinate transform that transforms the fake yield surfaces, 
$\theta =\pm \theta _{Y}(r)$, to 
$\varTheta =\pm 1$. Specifically we make the transformation of independent variables
\begin{equation} \left(r,\theta\right)=\left(r,\theta_Y(r)\varTheta\right), \end{equation}noting that this gives partial derivatives:
\begin{align} \frac{\partial}{\partial r} \to \frac{\partial}{\partial r} + \frac{\partial\varTheta}{\partial r}\frac{\partial}{\partial\varTheta} = \frac{\partial}{\partial r} + \frac{\epsilon\varPhi'\varTheta}{\alpha-\epsilon\varPhi}\frac{\partial}{\partial\varTheta}, \qquad \frac{\partial}{\partial \theta} \to \frac{\partial\varTheta}{\partial\theta}\frac{\partial}{\partial\varTheta} = \frac{1}{\alpha-\epsilon\varPhi}\frac{\partial}{\partial\varTheta}. \end{align}
In the plastically dominated region we seek a solution that is, to leading order, radial. The boundary layer may drive an azimuthal flow of order 
$\epsilon$, so we search for a solution of the form
\begin{align} u= u_0+\epsilon u_1 + \ldots, \quad v= \epsilon v_1 + \dots, \quad \psi= \psi_0 + \epsilon\psi_1 + \dots, \quad p=p_0+\epsilon p_1 +\dots. \end{align}We have three governing equations. The curl of the momentum balance gives
\begin{align} &\Bigg(\frac{1}{\left(\alpha-\epsilon\varPhi\right)^2}\frac{\partial^2}{\partial\varTheta^2}-\left(r\frac{\partial}{\partial r}+\frac{\epsilon r\varPhi'\varTheta}{\alpha-\epsilon\varPhi}\frac{\partial}{\partial\varTheta}\right)^2-2\left(r\frac{\partial}{\partial r}+\frac{\epsilon r\varPhi'\varTheta}{\alpha-\epsilon\varPhi}\frac{\partial}{\partial\varTheta}\right)\Bigg)\sin2\psi\nonumber\\ &\quad +\frac{2}{\alpha-\epsilon\varPhi}\frac{\partial}{\partial\varTheta}\left(r\frac{\partial}{\partial r}+ \frac{\epsilon r\varPhi'\varTheta}{\alpha-\epsilon\varPhi}\frac{\partial}{\partial\varTheta}+1\right) \cos2\psi=0; \end{align}conservation of mass gives
\begin{equation} \left(\frac{\partial}{\partial r}+\frac{\epsilon\varPhi'\varTheta}{\alpha-\epsilon\varPhi}\frac{\partial}{\partial\varTheta}\right)\left(ru\right)+\frac{1}{\alpha-\epsilon\varPhi}\frac{\partial v}{\partial\varTheta}=0;\end{equation}
and, provided we remain above 
$O({\textit {Bi}}^{-1})$ in (3.7a–d), at which point viscous terms re-enter the constitutive equation, the rigid plastic approximation holds and the orientation of the stress tensor, parameterised through 
$\psi$ (3.2), is given by
\begin{equation} 2\tan2\psi\left(r\frac{\partial u}{\partial r}+\frac{\epsilon r\varPhi'\varTheta}{\alpha-\epsilon\varPhi}\frac{\partial u}{\partial\varTheta}\right)=\frac{1}{\alpha-\epsilon\varPhi}\frac{\partial u}{\partial\varTheta}+r\frac{\partial v}{\partial r} + \frac{\epsilon r\varPhi'\varTheta}{\alpha-\epsilon\varPhi}\frac{\partial v}{\partial \varTheta} -v.\end{equation}3.2. Leading-order bulk flow
At 
$O(1)$ the solution of (3.8)–(3.10) corresponds to Nadai's (Reference Nadai1924) solution up to the slight alteration of the factors of 
$\alpha$ due to the scaling to unit angle. Specifically, we have
\begin{gather} \alpha\varTheta ={-}\psi_0+\frac{c}{\sqrt{c^2-1}}\arctan\left(\sqrt{\frac{c+1}{c-1}}\tan\psi_0\right), \end{gather}
\begin{gather}u_0={-}\frac{A}{r\left(c-\cos\left(2\psi_0\right)\right)}, \end{gather}
\begin{gather}p_0=2{\textit{Bi}}\,c\ln(r)+{\textit{Bi}}\,c\ln(c-\cos2\psi_0)+\text{const.} \end{gather}
The constant, 
$c$, is determined by the boundary condition 
$\psi _0(1)={\rm \pi} /4$, as required to match to a shear-stress-dominated boundary layer (see discussion after (3.4)). This gives the implicit equation
\begin{equation} \alpha+\frac{\rm \pi}{4}=\frac{c}{\sqrt{c^2-1}}\arctan\left(\sqrt{\frac{c+1}{c-1}}\right). \end{equation}
Finally, 
$A$ is determined by the volume-flux condition. To leading order,
\begin{equation} A=\frac{2c(c^2-1)}{4\alpha+{\rm \pi}+2c}. \end{equation}3.3. The intermediate layer
Asymptotic analysis of this leading-order solution (see Appendix A) verifies that the shear rate diverges as we approach the fake yield surface, 
$\theta =\theta _Y$, thus we will require an intermediate layer around 
$\varTheta =1$, to regularise the diverging terms. We proceed following Walton & Bittleston (Reference Walton and Bittleston1991) and Balmforth & Craster (Reference Balmforth and Craster1999), to derive a governing equation that reduces to the appropriate dynamical balances as we move out of the intermediate layer into either the bulk or boundary layer regions. We define this intermediate layer by the further coordinate transformation 
$\delta \zeta =\theta _Y-\theta =(\alpha -\epsilon \varPhi )(1-\varTheta )$, where 
$\zeta$ is our new angular coordinate and 
$\delta$ is another ordering parameter. We re-label velocity components
\begin{equation} \begin{pmatrix} u \\ v \end{pmatrix} = \begin{pmatrix} U \\ \epsilon V \end{pmatrix}, \end{equation}
where 
$U$ and 
$V$ are assumed to be 
$O(1)$. It is assumed, and later verified, that 
$\delta \ll \epsilon \ll 1$, so that the intermediate layer is much narrower than both the bulk region and the boundary layer. The 
$\epsilon$ scaling of the azimuthal velocity is inherited from the bulk solution (3.7a–d).
Note that,
\begin{equation} \delta\zeta=\alpha-\epsilon\varPhi(r)-\theta, \end{equation}and so, as before, we have altered partial derivatives
\begin{equation} \frac{\partial}{\partial r} \to \frac{\partial}{\partial r} -\frac{\epsilon\varPhi'}{\delta}\frac{\partial}{\partial \zeta}, \quad \frac{\partial}{\partial \theta} \to -\frac{1}{\delta} \frac{\partial}{\partial\zeta}. \end{equation} Expansion of the leading-order plastic solution for small 
$\delta$ (see Appendix A), gives
\begin{equation} u_0={-}\frac{A}{rc}-\frac{2A}{rc\sqrt{c}}\sqrt{\delta\zeta}+\ldots, \end{equation}thus we define
\begin{equation} U=U_0(r)+\delta^{1/2}U_1(r,\zeta)+\ldots,\quad \text{with} \ U_0(r)={-}\frac{A}{cr},\end{equation}
and the term 
$U_1(r,\zeta )$ capturing the weak variation of 
$U$ with 
$\zeta$ over the width of the intermediate layer. In addition to the assumption that 
$\delta \ll \epsilon \ll 1$, we will now also work under the assumption that 
$\epsilon ^2\ll \delta$ which ensures that the additional terms from the coordinate transformation do not enter the leading-order balance, and is also required for the second term in (3.19) to be greater than contributions to the velocity from the order-
$\epsilon$ correction terms in the bulk solution.
Under this assumption, verified below, expanding the shear stress gives
\begin{equation} \tau_{r\theta}={-}\frac{1}{r}\delta^{{-}1/2}\partial_\zeta U_1+{\textit{Bi}}-2r^2\delta{\textit{Bi}}\left(\frac{\partial_r U_0}{\partial_\zeta U_1}\right)^2+\ldots,\end{equation}
where we now use short-hand notation for partial derivatives for notational clarity. The second term on the right-hand side of (3.21) is constant and so does not contribute to the momentum balance, hence we have both viscous and plastic terms when 
$\delta ={\textit {Bi}}^{-2/3}$. The leading-order radial and angular expressions of momentum balance are now given by
\begin{gather} {\textit{Bi}}^{{-}1}\partial_r p = \frac{1}{r^2} \partial_{\zeta}\partial_{\zeta}U_1 + 2r\partial_\zeta \left[\left(\frac{\partial_rU_{0}}{\partial_\zeta U_1}\right)^2\right], \end{gather}
\begin{gather}{\textit{Bi}}^{{-}1/3}\partial_\zeta p = 2r\partial_\zeta\left[\frac{\partial_rU_0}{\partial_\zeta U_1}\right], \end{gather}
which, along with the fact that 
$p$ is 
$O({\textit {Bi}})$ from the bulk solution, imply 
$p={\textit {Bi}}\,P(r)=2c{\textit {Bi}}\ln {r}$ to leading order. Thus, the radial pressure gradient is unchanged to leading order by the intermediate layer.
The asymptotic behaviour of the shear stress in the bulk solution is given by
\begin{equation} \tau_{r\theta}={\textit{Bi}}\sin2\psi_0={\textit{Bi}}-2 c{\textit{Bi}}\delta\zeta+\ldots={\textit{Bi}}-2 c{\textit{Bi}}^{1/3}\zeta {\ldots}, \end{equation}
(see Appendix A). Then, integrating (3.22) once with respect to 
$\zeta$ and using (3.24) and (3.21) at 
$O({\textit {Bi}}^{1/3})$ to set the constant of integration, we establish a cubic equation in 
$\partial _\zeta U_1$,
\begin{equation} (\partial_\zeta U_1)^3-r^2\zeta \partial_r P(\partial_\zeta U_{1})^2+2r^3(\partial_rU_{0})^2=0,\end{equation}
in the intermediate layer. Substituting for 
$U_0(r)$ and 
$P(r)$, we find the form
\begin{equation} (\partial_\zeta U_1)^3-2cr\zeta (\partial_\zeta U_{1})^2+\frac{2A^2}{c^2r}=0.\end{equation}For the angular component of the velocity, from the incompressibility equation we have
\begin{equation} \frac{\partial}{\partial r}\left(rU\right)-\frac{\epsilon r\varPhi'}{\delta}\frac{\partial U}{\partial \zeta} -\frac{\epsilon}{\delta}\frac{\partial V}{\partial \zeta}=0, \end{equation}which, on substitution of (3.20), leads to
\begin{equation} V=V_0(r)-\delta^{1/2}r\varPhi'U_1(r,\zeta)+\ldots. \end{equation}
It remains to find 
$V_0$ and 
$\varPhi$ by solving in the boundary layer (§ 3.4).
We deduce some matching information for the boundary layer by considering the behaviour of the cubic equation for 
$\partial _\zeta U_{1}$, (3.26). As 
$\zeta \to -\infty$, in (3.26), the diverging positive second term can only balance with the first term, and so 
$\partial _\zeta U_1\to -\infty$ and we have
\begin{equation} \partial_\zeta U_{1}\sim 2cr\zeta. \end{equation}Hence, as we move into the boundary layer, we integrate (3.29) to find the velocity profile satisfies
\begin{equation} U=U_0(r)+cr{\textit{Bi}}^{{-}1/3}\zeta^2+{\textit{Bi}}^{{-}1/3}W(r)+\ldots, \end{equation}
where 
$W(r)$ is an arbitrary function of integration. Immediately, we identify the new distinguished scaling as 
${\textit {Bi}}^{-1/3}\zeta ^2=O(1)$, or 
$\zeta =O({\textit {Bi}}^{1/6})$ which, when related back in terms of 
$\theta$, gives 
$\theta -\theta _Y=O({\textit {Bi}}^{-1/2})$. Thus we have derived the width of the boundary layer as 
$\epsilon ={\textit {Bi}}^{-1/2}$, and can verify that 
$\epsilon ^2={\textit {Bi}}^{-1}\ll \delta \ll \epsilon$, as was assumed in deriving the leading order stress balance in the intermediate layer above.
3.4. Boundary layer solution
The boundary layer solution is constructed by introducing the scaled coordinate, 
$\epsilon \tilde {\phi } = \alpha -\theta$. From mass conservation we find the natural scalings
\begin{equation} \begin{pmatrix} u \\ v \end{pmatrix} = \begin{pmatrix} \tilde{U} \\ \epsilon \tilde{V} \end{pmatrix}, \end{equation}
where 
$\tilde {U}$ and 
$\tilde {V}$ are assumed 
$O(1)$. Then, expanding the shear stress gives
\begin{equation} \tau_{r\theta} ={-}\frac{1}{\epsilon r}\frac{\partial \tilde{U}}{\partial\tilde{\phi}} +{\textit{Bi}} + O(\epsilon^2{\textit{Bi}}),\end{equation}
where the sign of the 
${\textit {Bi}}$ term on the right-hand side arises due to 
$\partial \tilde {U}/\partial \tilde {\phi }<0$ (since 
$U$ varies from 0, at the wall, to some negative value in the bulk of the flow). Thus, viscous stresses balance the pressure gradient since 
$\epsilon ^3{\textit {Bi}}={\textit {Bi}}^{-1/2}\ll 1$ and so the plastic shear-stress terms are omitted in the radial momentum balance. We find that, as in the intermediate layer, the pressure is constant across the extent of the boundary layer to leading order, and is thus set by the pressure in the bulk and intermediate regions. The radial momentum equation then becomes
\begin{equation} \frac{\partial P}{\partial r} \equiv \frac{2c}{r} = \frac{1}{r^2}\frac{\partial^2 \tilde{U}}{\partial\tilde{\phi}^2}, \end{equation}which integrates to give
\begin{equation} \tilde{U} = cr(\tilde{\phi}^2 - 2\varPhi(r)\tilde{\phi}). \end{equation}
The constants of integration have been chosen so that 
$\partial \tilde {U}/\partial \tilde {\phi }=0$ at 
$\tilde {\phi }=\varPhi$, and 
$\tilde {U}=0$ at 
$\tilde {\phi }=0$.
As 
$\tilde {\phi }\to \varPhi$ we may write
\begin{equation} \tilde{U}={-}cr\varPhi^2+cr(\tilde{\phi}-\varPhi)^2, \end{equation}
which is consistent with the 
$\zeta \to -\infty$ limit, (3.30), and, from (3.20), determines 
$\varPhi$,
\begin{equation} \varPhi(r)=\frac{\sqrt{A}}{cr}. \end{equation}
Thus we have found that the angular width of the boundary layer, 
$\epsilon \varPhi$, does indeed depend on 
$(r^2{\textit {Bi}})^{-1/2}$, as anticipated. The Cartesian width of this boundary layer is 
$\epsilon r\varPhi ={\textit {Bi}}^{-1/2}\sqrt {A}/c$. It is therefore a constant, independent of the radial distance from the apex of the wedge. We can also now calculate the additional shear stress at the wall due to the viscoplastic boundary layer. From (3.32) we have
\begin{equation} \left.\tau_{r\theta}\right\vert_{\theta=\alpha}={\textit{Bi}}-\frac{1}{\epsilon r} \left.\frac{\partial\tilde{U}}{\partial\tilde{\phi}}\right\vert_{\tilde{\phi}=0}+\ldots= {\textit{Bi}}+\sqrt{{\textit{Bi}}}\,\frac{2\sqrt{A}}{r}+\ldots, \end{equation}
demonstrating that the additional shear stress at the wall is 
$O({\textit {Bi}}^{1/2})$ and is proportional to 
$1/r$. Thus, the additional shear force over the extent of the upper wall, from 
$r=r_1$ to 
$r=r_2$, assuming the plastic regime applies throughout this domain, is given by
\begin{equation} 2\sqrt{{\textit{Bi}}}\sqrt{A}\ln\left(\frac{r_2}{r_1}\right). \end{equation}For the polar velocity, using incompressibility, we have
\begin{equation} \frac{\partial \tilde{V}}{\partial\tilde{\phi}}={-}\frac{1}{\epsilon}\frac{\partial \tilde{V}}{\partial\theta} ={-}\frac{\partial v}{\partial\theta}=\frac{\partial}{\partial r}(r\tilde{U})=2cr\tilde{\phi}^2-2\sqrt{A}\tilde{\phi}. \end{equation}
Integrating and using 
$\tilde {V}=0$ at 
$\tilde {\phi }=0$, yields
\begin{equation} \tilde{V}=\frac{cr}{3}(2\tilde{\phi}^3-3\varPhi\tilde{\phi}^2). \end{equation}
In particular, by matching 
$\tilde {V}$ to 
$V_0$ (see (3.28)), we find the negative angular flow,
\begin{equation} V_0(r)={-}\frac{cr\varPhi^3}{3}={-}\frac{A\sqrt{A}}{3c^2r^2}. \end{equation}This velocity away from the boundary is transmitted through the intermediate layer to the bulk region. To determine the profile of this weak angular flow in the body of the wedge requires exploring higher orders in the bulk solution. This is done in § 3.5.
It is interesting to consider the limit 
$\alpha \to 0$ with 
$r\alpha =1$, for which 
$A/c^2\to 1$. In this limit: the leading-order bulk solution reduces to a uniform plug flow; the intermediate layer vanishes; the leading-order solution in the uniform width, 
$O({\textit {Bi}}^{-1/2})$, boundary layer becomes a parabolic profile; and the outflow from the boundary layer, 
$\epsilon V_0$, vanishes. Thus the Bingham–Poiseuille flow in a channel is recovered, as it must be, in this limit.
3.5. Higher orders in the bulk
At 
$O(\epsilon )$, the governing system of equations (3.8)–(3.10) becomes
\begin{align} &\left(\frac{1}{\alpha^2}\frac{\partial^2}{\partial\varTheta^2}-\left(r\frac{\partial}{\partial r}\right)^2-2r\frac{\partial}{\partial r}\right)\left(\psi_1\cos2\psi_0\right)-\frac{2}{\alpha}\frac{\partial}{\partial\varTheta}\left(r\frac{\partial}{\partial r}+1\right)\left(\psi_1\sin2\psi_0\right) \nonumber\\ &\quad ={-}\!\left(\frac{\varPhi}{\alpha^3}\frac{\partial^2}{\partial\varTheta^2}+\frac{\varPhi\varTheta}{2\alpha}\frac{\partial}{\partial\varTheta}\right)\sin2\psi_0+\frac{\varPhi\varTheta}{\alpha^2}\frac{\partial^2}{\partial\varTheta^2}\cos2\psi_0, \end{align}
\begin{gather} \frac{\partial\left(ru_1\right)}{\partial r}+\frac{1}{\alpha}\frac{\partial v_1}{\partial\varTheta}=\frac{\varPhi\varTheta}{\alpha}\frac{\partial u_0}{\partial \varTheta},\end{gather}
\begin{gather} 2r\tan2\psi_0\frac{\partial u_1}{\partial r}-\frac{1}{\alpha}\frac{\partial u_1}{\partial\varTheta}-r\frac{\partial v_1}{\partial r}+v_1-4\psi_1u_0\sec^22\psi_0 \nonumber\\ =2\frac{\varPhi\varTheta}{\alpha}\tan2\psi_0\frac{\partial u_0}{\partial\varTheta}+\frac{\varPhi}{\alpha^2}\frac{\partial u_0}{\partial\varTheta}, \end{gather}
where the identities, 
$r({\textrm {d}\varPhi }/{\textrm {d}r})=-\varPhi$ and 
$r({\partial u_0}/{\partial r})=-u_0$, have been used. The right-hand sides of (3.42)–(3.44) determine the 
$r$ dependence of 
$u_1$, 
$v_1$ and 
$\psi _1$, and so we define
\begin{equation} u_1=\frac{\hat{u}_1(\varTheta)}{r^2}, \quad v_1=\frac{\hat{v}_1(\varTheta)}{r^2}, \quad \text{and} \quad \psi_1=\frac{\hat{\psi}_1(\varTheta)}{r}. \end{equation}
Making these substitutions, along with 
$\varPhi ={\sqrt {A}}/{cr}$, results in the coupled ordinary differential equations (ODEs)
\begin{gather} \frac{1}{\alpha^2}(\hat{\psi}_1\cos2\psi_0)''+\hat{\psi}_1\cos2\psi_0\nonumber\\ = \frac{\sqrt{A}}{\alpha c}\left(\frac{\varTheta}{\alpha}\left(\cos2\psi_0\right)''-\frac{1}{\alpha^2}\left(\sin2\psi_0\right)''-\frac{\varTheta}{2}\left(\sin2\psi_0\right)' \right), \end{gather}
\begin{gather} -\hat{u}_1+\frac{1}{\alpha}\hat{v}_1'=\frac{\sqrt{A}\varTheta}{\alpha c}\widehat{u_0}', \end{gather}
\begin{gather} -4\hat{u}_1\tan2\psi_0-\frac{1}{\alpha}\hat{u}_1'+3\hat{v}_1-4\hat{\psi}_1\hat{u}_0\sec^22\psi_0=\frac{\sqrt{A}}{\alpha c}\left(2\varTheta\tan2\psi_0+\frac{1}{\alpha}\right)\hat{u}_0'. \end{gather}
There are two boundary conditions for 
$\hat {v}_1(\varTheta )$
\begin{equation} \hat{v}_1(0)=0, \quad \hat{v}_1(1)={-}\frac{A\sqrt{A}}{3c^2},\end{equation}which respectively represent the symmetry of the flow at the mid-line and matching to the angular outflow from the boundary layer. Additionally, the symmetry conditions demand that
\begin{equation} \hat{\psi}_1(0)=0.\end{equation}
The final boundary condition arises from matching the deviatoric stresses between the intermediate layer and the bulk solution. Specifically, in the intermediate layer, using the transformed partial derivatives, (3.18a,b), the asymptotic form of the radial velocity, (3.20), and the identity 
$r\varPhi '=-\varPhi$, we have (in the notation of the intermediate layer)
\begin{equation} \tau_{rr}={-}\frac{2r{\textit{Bi}}\,\delta^{1/2}}{\partial_\zeta U_1}\left(\frac{\textrm{d} U_0}{\textrm{d} r}-\frac{\epsilon}{\delta^{1/2}}\varPhi'\partial_\zeta U_1\right)+\ldots= \frac{2U_0}{\partial_\zeta U_1}{\textit{Bi}}^{2/3}-2\varPhi{\textit{Bi}}^{1/2}+\ldots. \end{equation}
By considering the 
$\zeta \to \infty$ limit of the cubic (3.26), we find that 
$\partial _\zeta U_1\to -A/(r\sqrt {c^3\zeta })$, and so
\begin{equation} \tau_{rr}=2\sqrt{c\zeta}{\textit{Bi}}^{2/3}-2\varPhi{\textit{Bi}}^{1/2}+\ldots. \end{equation}Meanwhile, in the bulk solution, we have
\begin{gather} \tau_{rr}={\textit{Bi}}\cos2\psi_0-2\epsilon{\textit{Bi}}\,\psi_1\sin2\psi_0+\ldots=2{\textit{Bi}}\sqrt{\alpha c\left(1-\varTheta\right)}-2{\textit{Bi}}^{1/2}\psi_1+\ldots, \end{gather}
where we have used the expansion of 
$\cos 2\psi _0$ and 
$\sin 2\psi _0$ about 
$\varTheta =1$, given in Appendix A. Finally, using the substitution 
$\delta \zeta =\theta _Y(1-\varTheta )=\alpha (1-\varTheta )+O(\epsilon \delta )$, we see that the leading-order terms match, and that matching at 
$O({\textit {Bi}}^{1/2})$ requires
\begin{equation} \psi_1(r,\varTheta=1)=\varPhi(r)=\frac{\sqrt{A}}{cr} \implies \hat{\psi}_1\left(1\right)=\frac{\sqrt{A}}{c}. \end{equation}
Furthermore, we may determine the behaviour of all the dependent variables in the regime 
$\lvert 1-\varTheta \rvert \ll 1$ (see Appendix B, setting 
$N=1$ for the Bingham model).
The boundary-value problem (3.46)–(3.48) was solved using a shooting method. First we use the local form of the dependent variables (B6)–(B8) to step a small distance, 
$d$, away from the singular point before integrating to 
$\varTheta =0$ where the boundary conditions (3.49a,b) and (3.50) determine the unknown coefficients, 
$E$ and 
$F$, in (B6)–(B8). The dependence of the solutions on 
$d$ was investigated and it was found that profiles were essentially independent of 
$d$ for 
$d<10^{-5}$. Numerical solutions are plotted for a range of values of 
$\alpha$ in figure 3. We note that the radial velocity is enhanced by the effects of the boundary layer, and that an angular velocity away from the boundary is induced, which reaches a maximum at some interior location. Also, the magnitude of the shear stress is promoted throughout the domain. The magnitude of these corrections decreases systematically with increasing wedge half-angle.
Figure 3. The first-order corrections to the velocities, 
$u$ and 
$v$, and the stress orientation function, 
$\psi$, in the plastic regime, as functions of the scaled angular coordinate 
$\varTheta$, for 
$\alpha ={\rm \pi} /6$ (solid), 
${\rm \pi} /4$ (dashed), 
${\rm \pi} /3$ (dash-dotted) and 
${\rm \pi} /2$ (dotted).
Physically, as we deviate from a perfectly rigid plastic, the radial flow is reduced by no slip at the wall, over an increasingly thick boundary layer. Thus, for the same volume flux, the radial flow must be enhanced in the centre of the wedge, and conservation of mass requires an angular flow from the wall, where radial flow is small, towards the centre of the wedge, where the flow is enhanced. The dependence on half-angle arises via the decreasing boundary layer strength with increasing half-angle. This occurs through a competition between the radial pressure gradient at the edge of the boundary layer (which decreases with increasing 
$\alpha$, supporting a lower angular gradient of the shear stress potentially widening the boundary layer) and the magnitude of the non-vanishing velocity at the boundary in Nadai's (Reference Nadai1924) plastic solution (2.14a–c) (which decreases with 
$\alpha$, thus requiring a thinner boundary layer to enforce no slip for the same radial pressure gradient). The balance between these effects is reflected in the 
$\sqrt {A}/c$ factor in the boundary layer width, 
$\epsilon \varPhi$. This term turns out to be a decreasing function of 
$\alpha$, and so the boundary layer and the induced velocity perturbations are weaker for larger half-angles (at the same radial position, 
$r$). The stress behaviour is also non-trivial since the enhanced radial flow has the potential to increase both the normal radial stress, and the shear stress, but cannot increase both, since the magnitude of the deviatoric stress is everywhere equal to the yield stress at this order. However, as we have enhanced shear stress at the walls due to the viscoplastic shear layer, it is natural to expect the shear stress to be enhanced across the domain, as is observed.
3.6. Composite solutions
To produce leading-order radial and angular velocity profiles, it is useful to define a composite solution that is valid across the whole domain. We define such a solution, using the limiting behaviours of the plastic and boundary layer solutions in the vicinity of 
$\theta _Y$, as
\begin{gather} u_{comp}(r,\theta)= \begin{cases} \displaystyle u_0+U+\frac{A}{rc}+\frac{2A\sqrt{\theta_Y-\theta}}{rc\sqrt{c}} & \text{for } \theta<\theta_Y, \\ \displaystyle\tilde{U}+U+\frac{A}{rc}-rc{\textit{Bi}}\left(\theta-\theta_Y\right)^2 & \text{for } \theta\geqslant\theta_Y, \end{cases} \end{gather}
\begin{gather}v_{comp}(r,\theta)= \begin{cases} \displaystyle\epsilon \left(v_1 + V +\frac{A\sqrt{A}}{3r^2c^2}+\frac{2A\sqrt{A}}{r^2c^2\sqrt{c}}\sqrt{\theta_Y-\theta}\right) & \text{for } \theta<\theta_Y, \\ \displaystyle\epsilon \left(\tilde{V}+V+\frac{A\sqrt{A}}{3r^2c^2}-\sqrt{A}{\textit{Bi}}\,\left(\theta-\theta_Y\right)^2\right) & \text{for } \theta\geqslant\theta_Y, \end{cases} \end{gather}
where 
$\theta _Y(r)$ is the location of the fake yield surface defined by (3.3) and (3.36). These profiles (for fixed 
$r$) are continuous at 
$\theta =\theta _Y$, since the velocities here are equal to those for the intermediate layer solution at 
$\zeta =0$, on both branches of the piece-wise definition.
3.7. Numerical simulations
Numerical simulations of the complete governing partial differential equations (2.4)–(2.6), were carried out using the finite element method as implemented by FEniCS (Logg, Mardal & Wells Reference Logg, Mardal and Wells2012; Alnæs et al. Reference Alnæs, Blechta, Hake, Johansson, Kehlet, Logg, Richardson, Ring, Rognes and Wells2015). The domains were meshed with 60 000 triangular elements, with a higher concentration at the wall and near the outflow, where higher strain rates are expected. The fluid is yielded everywhere, so we do not expect the occurrence of plug regions and the associated singular Jacobian matrix of the discretised equations, when applying the Newton method. Nonetheless, the nonlinearity in the equations causes problems when trying to use Newton iterations to converge to the solution from an arbitrary initial state. We found that using an augmented-Lagrangian method, as described by Saramito (Reference Saramito2016), resulted in slow convergence but could be used to produce an effective initial guess for the Newton method, allowing a subsequent fast convergence to a Newton residual of less than 
$10^{-8}$. For the following results 
${\textit {Bi}}$ has been set at 1, but the domains have a radial range of 
$r\in [0.01,500]$, allowing the effective Bingham number, 
$r^2{\textit {Bi}}$, to vary from 
$10^{-4}$ to 250 000, within a single simulation. The low value of this parameter at the outflow means the purely viscous radial solution (2.15a–c) is a good approximation to the velocity at the outflow, and is therefore imposed as a Dirichlet boundary condition here, while the large value at the inflow means the leading-order asymptotic viscoplastic solution is an appropriate boundary condition here.
Firstly, we confirm the predicted behaviour of the boundary layer thickness in the plastic regime. For a given radius, the boundary layer thickness was determined from the numerical radial flow profile by the angle at which
\begin{equation} \frac{\partial u}{\partial \theta}= \left(2{\textit{Bi}}\,\frac{A^2}{c^2r}\right)^{1/3}, \end{equation}
which is the shear rate predicted at the centre of the matching region by setting 
$\zeta =0$ in the equation (3.26). Figure 4 shows the measured boundary layer thickness alongside the predicted value, 
$\epsilon \varPhi (r)$, with 
$\varPhi$ given by (3.36), at a number of radial positions for three different wedge angles, confirming the analytical result. Also shown are example velocity profiles from the numerical solution for 
$\alpha ={\rm \pi} /4$, in a small region near the upper boundary 
$\theta =\alpha$. We note that the boundary layer is of constant width, in a Cartesian sense, due to the 
$1/r$ dependence of 
$\varPhi$, and that this is borne out by the numerical results.
Figure 4. (a) The width of the viscoplastic boundary layer determined analytically (dashed) and numerically (symbols), at a range of radial positions and wedge angles. (b) Profiles of radial velocity determined numerically for 
$\alpha ={\rm \pi} /4$. The quantity 
$r(\alpha -\theta )$ defines a small Cartesian distance from the wall at 
$\theta =\alpha$. The theoretical position of the fake yield surface is shown by a dotted line.
Figure 5 provides comparisons between the numerical and composite asymptotic flow profiles for a choice of wedge angle (
$\alpha ={\rm \pi} /4$) and radius (
$r=100$), which is sufficiently distant from the inflow and outflow boundaries to exhibit the fully developed profile. The numerical simulations accurately reproduce the predicted velocity profiles in both the bulk and boundary layer regions. In particular, the predicted weak negative angular flow is confirmed in these simulations.
Figure 5. The asymptotic (dashed line) and numerical (circles) radial and angular velocities as function of angle, 
$\theta$, for 
$\alpha ={\rm \pi} /4$, 
${\textit {Bi}}=1$ and 
$r=100$. The numerical data points shown are only a small sub-sample of the full numerical solutions, which are of a much higher resolution. (a,b) Plot the radial and angular velocities respectively, across the domain. (c,d) Are close-ups of the boundary layers for (a,b) respectively. The position of the fake yield surface, 
$\theta _Y$, is depicted by a red dotted line in all panels (and is almost indistinguishable from the boundary, 
$\theta ={\rm \pi} /4$, in the upper two panels).
Finally, we note the small discrepancy in panel (b) of figure 5 and explore whether this can be explained by the absence of higher-order terms in the asymptotic expansions. This is confirmed by figure 6, which shows the difference between the leading-order composite solutions and the numerical solutions, measured by the 
$L_2$-norm at fixed 
$r$
\begin{gather} \Delta u(r)=\sqrt{\int_0^\alpha\left[u_{num}(r,\theta)-u_{comp}(r,\theta)\right]^2\,\textrm{d}\theta}, \end{gather}
\begin{gather}\Delta v (r)=\sqrt{\int_0^\alpha\left[v_{num}(r,\theta)-v_{comp}(r,\theta)\right]^2\,\textrm{d}\theta}. \end{gather}
For the velocity in the radial direction, the composite solution neglects terms of 
$O({\textit {Bi}}^{-1/2})$, which results in 
$\Delta u$ proportional to 
${\textit {Bi}}^{-1/2}r^{-2}=r^{-2}$ for fixed 
${\textit {Bi}}$. For the velocity in the angular direction we must consider the order of terms neglected above 
$O({\textit {Bi}}^{-1/2})$. From (3.30), we expect an 
$O({\textit {Bi}}^{-1/3})$ correction to the radial flow in the boundary layer. This in turn would drive a correction of 
$O(\epsilon {\textit {Bi}}^{-1/3})={\textit {Bi}}^{-5/6}$ in the velocity in the angular direction at the top of the boundary layer. Hence, without calculating higher orders explicitly, we expect the second-order correction to the velocity to be 
$O({\textit {Bi}}^{-5/6})$ or, including radial dependence, 
$O(r^{-8/3}{\textit {Bi}}^{-5/6})$, since the expansion parameter is 
$r^2{\textit {Bi}}$ and the leading-order flow has 
$r^{-1}$ dependence. Hence, we expect 
$\Delta v$ to scale like 
$r^{-8/3}$. Figure 6 confirms this relationship for the three wedge angles studied, providing strong evidence for the validity of our asymptotic structure. The deviation from the predicted trends at 
$r\approx 500$ is due to the leading-order profiles being imposed here as boundary conditions, and so the difference is at the level of the implementation precision.
Figure 6. The difference between the leading-order asymptotic and numerical solutions for the velocity fields, evaluated as 
$\Delta u$ and 
$\Delta v$, as a function of radial position, for 
$\alpha ={\rm \pi} /6$, 
${\rm \pi} /4$ and 
${\rm \pi} /3$. The slope markers show the predicted slopes of 
$-2$ (top) and 
$-8/3$ (bottom).
Figure 7 shows a density plot of log strain-rate and radial velocity profiles from the numerical solution for the wedge of half-angle 
$\alpha ={\rm \pi} /3$. These demonstrate how the plastically dominated regime, with a thin boundary layer structure, is valid for large values of 
$r$, but transitions to a different regime, without a boundary layer, towards the apparent apex of the wedge. For sufficiently small values of 
$r$ the viscous stresses are dominant and a new asymptotic solution can be derived. This is detailed in the next section.
Figure 7. (a) A density plot of log strain rate from the 
$\alpha ={\rm \pi} /3$ numerical simulation. (b) The scaled radial velocity, 
$ru$, as a function of angle, 
$\theta$, for 
$\alpha ={\rm \pi} /3$ and 
$r=1$ (dotted), 
$4$ (dashed) and 
$20$ (solid). The 
$\theta$ axis is reversed so that the no-slip boundary is located at the left of the figure.
4. The viscous regime
In the regime of small Bingham number, 
${\textit {Bi}}\ll 1$, or, equivalently, at small distances from the apex of the wedge, the viscoplastic problem is a regular perturbation of the viscous Stokes problem, and so we do not require a boundary layer structure. This regime has been studied previously by Sandru & Camenshi (Reference Sandru and Camenshi1979), so we initially outline this solution (with notation adjusted for consistency) before exploring it in greater detail by carrying out numerical integration not evaluated by Sandru & Camenshi (Reference Sandru and Camenshi1979), and comparing to the results of direct numerical simulations in § 4.1.
In this regime the leading-order flow is given by the viscous solution (2.15a–c), which produces an 
$O({\textit {Bi}})$ term in the stress. Hence the solution takes the form
\begin{equation} u= u_0 + {\textit{Bi}}\, u_1 + \ldots, \quad v= {\textit{Bi}} \,v_1 + \ldots, \quad p= p_0 + {\textit{Bi}}\,p_1 + \dots, \end{equation}
with 
$u_0$ and 
$p_0$ given by (2.15a–c). By conservation of mass and for homogeneity in 
$r$ at 
$O({\textit {Bi}})$ in the expression of conservation of momentum, we define
\begin{equation} u_1=rf'(\theta) \quad \text{and} \quad v_1={-}2rf(\theta), \end{equation}
where 
$f(\theta )$ is to be determined and 
$'$ denotes differentiation with respect to 
$\theta$. Note that this 
$r$ dependence of 
$u_1$ and 
$v_1$ was anticipated, due to the expansion parameter, in fact, being 
$r^2{\textit {Bi}}$ and the leading-order velocity depending upon 
$r^{-1}$. By eliminating pressure from the conservation of momentum equations, we find that 
$f(\theta )$ satisfies the differential equation
\begin{equation} f'''(\theta)+4\,f'(\theta)+H(\theta)=C, \end{equation}
where 
$C$ is a constant of integration and
\begin{equation} H(\theta)=\frac{\textrm{d}}{\textrm{d} \theta}\left(\frac{\sin2\theta}{\sqrt{\sin^22\theta+\left(\cos2\alpha-\cos2\theta\right)^2}}\right)+\frac{2\left(\cos2\theta-\cos2\alpha\right)}{\sqrt{\sin^22\theta+\left(\cos2\alpha-\cos2\theta\right)^2}}. \end{equation}The governing equation (4.3) is subject to four boundary conditions
\begin{equation} f(0)=f''(0)=0, \quad f(\alpha)=f'(\alpha)=0, \end{equation}
which represent symmetry at 
$\theta =0$ and no slip at 
$\theta =\alpha$, which are sufficient to evaluate 
$f(\theta )$ and 
$C$. Having determined 
$C$, the pressure perturbation is given by
\begin{equation} p_1=C\ln(r)+\frac{1+2\cos^22\alpha-3\cos2\alpha\cos2\theta}{\sqrt{\sin^22\theta+\left(\cos2\theta-\cos2\alpha\right)^2}}+\text{const.} \end{equation} We integrate (4.3) numerically by first finding a particular solution from an initial value problem with initial conditions given at 
$\theta =0$, then determining 
$C$ and the coefficients of the complementary solution by the boundary conditions at 
$\alpha$. Specifically, we find a solution of the form
\begin{equation} f(\theta)=f_p(\theta)+\frac{C}{4}\theta+D\sin2\theta, \end{equation}
where 
$f_p$ is the unique solution to the initial value problem
\begin{equation} f_p'''(\theta)+4f_p'(\theta)+H(\theta)=0, \end{equation}
subject to 
$f_p(0)=f_p'(0)=f_p''(0)=0$. Then 
$C$ and 
$D$ are determined by the boundary conditions at 
$\theta =\alpha$
\begin{equation} C= \frac{4(2f_p(\alpha)\cos2\alpha-f_p'(\alpha)\sin2\alpha)}{\sin2\alpha-2\alpha\cos2\alpha} \quad \text{and} \quad D= \frac{\alpha f_p'(\alpha)-f_p(\alpha)}{\sin2\alpha-2\alpha\cos2\alpha}. \end{equation} Velocity profiles for the first-order corrections to the leading-order viscous flow are given in figure 8. Note, again, the negative angular velocity, as anticipated. Also, the radial velocity is reduced at the centre of the wedge and enhanced towards the walls, as required to flatten the radial velocity profile as 
$r$ increases. In this regime the magnitude of the corrections increases with increasing half-angle. This can be explained by considering the magnitude of the deviatoric stress in the leading-order solution,
\begin{equation} \tau\equiv\sqrt{\tau_{rr}^2+\tau_{r\theta}^2}=\frac{2\sqrt{1+\cos^22\alpha-2\cos2\theta\cos2\alpha}}{r^2\left(\sin2\alpha-2\alpha\cos2\alpha\right)}, \end{equation}
which can be shown to decrease with increasing 
$\alpha$. Thus, the fluid is less strongly stressed by the leading-order flow for larger wedge half-angles, and so the plasticity has a larger impact on the flow.
Figure 8. The first-order corrections to the velocities, 
$u$ and 
$v$, in the viscous regime, as functions of the scaled angle, 
$\theta /\alpha$, for a wedge of half-angle 
$\alpha ={\rm \pi} /6$ (solid), 
${\rm \pi} /4$ (dashed), 
${\rm \pi} /3$ (dash-dotted) and 
${\rm \pi} /2$ (dotted).
The additional shear stress at the wall due to the viscoplasticity in this regime is given by
\begin{equation} {\textit{Bi}}\,\left.\tau_{r\theta}^{(1)}\right|_{\theta=\alpha}={\textit{Bi}}+{\textit{Bi}}\left.\left(\frac{1}{r}\frac{\partial u_1}{\partial \theta}+\frac{\partial v_1}{\partial r}-\frac{v_1}{r}\right)\right|_{\theta=\alpha}={\textit{Bi}}(f''(\alpha)+1), \end{equation}
where 
$f''(\alpha )$ is determined numerically. This correction to the shear stress at the wall is plotted as a function of 
$\alpha$ in figure 9(a). An interesting feature of this profile is the minimum attained for an angle of 
$\alpha \approx 32^{\circ }$. The behaviour of the constant 
$C$ is also of particular interest, as the first-order correction to the radial pressure gradient is given from (4.6) by 
${\textit {Bi}}\,C/r$. The behaviour of 
$C=r\partial p_1/\partial r$ as a function of 
$\alpha$ is given in figure 9b), showing that 
$C$ varies from 2.55 (approximately) at 
$\alpha ={\rm \pi} /2$ and diverges as 
$\alpha \to 0$. In particular, we find that the perturbation to the shear stress at the wall, 
$\tau _{r\theta }^{(1)}\sim 3/2$, and 
$C\sim 3/(2\alpha )$ for 
$\alpha \ll 1$ as required for consistency with the plane Poiseuille solution in the limit 
$\alpha \to 0$ with 
$r\alpha =1$ fixed (see Appendix D).
Figure 9. Features of the perturbation solution for the planar converging flow of a viscoplastic in the low Bingham number regime: (a) the additional wall shear stress, 
$\tau _{r\theta }^{(1)}|_{\theta =\alpha }$, and (b) the scaled additional radial pressure gradient, 
$r\partial p_1/\partial r=\mathcal {C}$, as functions of wedge half-angle, 
$\alpha$.
4.1. Numerical simulations
For numerical simulations in the viscous regime, a radial extent of 
$r\in [0.1,10]$ and a range of small values for 
${\textit {Bi}}$ were used. Again the simulations were run in FEniCS, on a mesh of 60 000 elements. Here, we imposed the purely viscous radial solution (2.15a–c) as Dirichlet boundary conditions at both inflow and outflow, and examined the flow profiles at 
$r=1$, which was found to be sufficiently far from these boundaries not to be affected by the approximations made in the boundary conditions. Due to the low values of the Bingham number, all simulations were solved directly by Newton iteration from the purely radial Stokes solution (2.15a–c) converging to residuals of less than 
$10^{-8}$.
Figure 10 shows example velocity profiles for 
$\alpha ={\rm \pi} /4$, and with 
${\textit {Bi}}=0.01$, demonstrating good agreement between asymptotic theory and numerical computation. As with the plastic regime, we measured the difference between the asymptotic profiles and the numerical solutions for 
$\alpha ={\rm \pi} /6$, 
${\rm \pi} /4$ and 
${\rm \pi} /3$, this time at fixed 
$r=1$ and variable 
${\textit {Bi}}$. With the absence of a boundary layer structure, there is no need for composite solutions, and so we measure the differences, 
$\Delta u$ and 
$\Delta v$, between the first-order asymptotic expressions and the numerical solutions. Both differences were found to scale like 
${\textit {Bi}}\,^2$, as expected having neglected second-order terms in the expansions.
Figure 10. The asymptotic (dashed line) and numerical (circles) radial and angular velocities as functions of angle, 
$\theta$, for 
${\textit {Bi}}=0.01$ and 
$r=1$. (a,b) Show the radial and angular profiles, respectively, for a wedge of half-angle 
$\alpha ={\rm \pi} /4$.
5. Herschel–Bulkley flow
In this section we extend our examination of converging flows to Herschel–Bulkley fluids. Importantly, the key features of the boundary layer structure determined for Bingham fluids when 
${\textit {Bi}}\gg 1$ (§ 3) carry over to Herschel–Bulkley fluids and here we generalise our analysis to these materials. The Herschel–Bulkley viscoplastic model is defined by the dimensional constitutive law
\begin{equation} \boldsymbol{\tau} = \left(K\dot{\gamma}^{N-1}+\frac{\tau_c}{\dot{\gamma}}\right) \dot{\boldsymbol{\gamma}}, \end{equation}
where 
$K$ is the consistency, 
$N$ the flow index, and 
$\tau _c$ the yield stress. This model extends the Bingham model, which is reproduced for 
$N=1$ and 
$K=\mu$, by allowing for shear thickening (
$N>1$) and shear thinning (
$N<1$). For this model we have the critical distance
\begin{equation} r_c=\left(\frac{KQ^{N}}{\tau_c}\right)^{{1}/{2N}}, \end{equation}
above which the flow is dominated by plastic stresses. For a typical radial distance, 
$R$, we scale velocities by 
$Q/R$, strain rates by 
$Q/R^2$ and stresses and pressures by 
$K Q^N/R^{2N}$ resulting in the non-dimensional constitutive equation
\begin{equation} \begin{pmatrix} \tau_{rr} \\ \tau_{r\theta} \end{pmatrix} = \left(\dot{\gamma}^{N-1}+\frac{{\textit{Bi}}}{\dot{\gamma}}\right) \begin{pmatrix} 2\dfrac{\partial u}{\partial r} \\ \dfrac{1}{r}\dfrac{\partial u}{\partial\theta}+\dfrac{\partial v}{\partial r}-\dfrac{v}{r} \end{pmatrix}, \end{equation}where the Bingham number is now given by
\begin{equation} {\textit{Bi}}=\frac{\tau_c R^{2N}}{KQ^{N}}. \end{equation}
The critical radial length scale is equivalent to asserting that the asymptotic analysis will break down for 
$r=O({\textit {Bi}}^{-{1}/{2N}})$ in non-dimensional variables.
In the 
${\textit {Bi}}\gg 1$ regime (or equivalently at large distances from the apex of the wedge), the bulk solution is unchanged at leading order since the viscous stresses are neglected here, and hence is given by (3.11)–(3.13). In the intermediate layer, using the same notation as for the Bingham case, albeit with a different small parameter, 
$\delta$, the radial velocity takes the form
\begin{equation} U=U_0(r)+\delta^{1/2}U_1(r,\zeta), \end{equation}and the shear stress is given by
\begin{equation} \tau_{r\theta}={-}\frac{1}{r^N}\delta^{{-}N/2}\left|\partial_\zeta U_1\right|^{N-1}\partial_\zeta U_1+{\textit{Bi}}-2r^2\delta{\textit{Bi}}\left(\frac{\partial_r U_0}{\partial_\zeta U_1}\right)^2+\ldots. \end{equation}
Hence, a balance of viscous and plastic stresses is achieved in an intermediate layer of width 
$\delta$, provided
\begin{equation} \delta^{{-}N/2}=\delta{\textit{Bi}} \implies \delta={\textit{Bi}}^{-{2}/({N+2})}, \end{equation}
and we find that 
$U_1$ satisfies the equation
\begin{equation} \lvert\partial_\zeta U_{1}\rvert^{N+1}\partial_\zeta U_{1}-2cr^N\zeta\lvert\partial_\zeta U_{1}\rvert^2+2r^{N-2}\frac{A^2}{c^2}=0,\end{equation}
where 
$A$ and 
$c$ are defined as previously for the bulk solution and given by (3.14) and (3.15). In the limit 
$\zeta \to -\infty$, 
$\partial _\zeta U_{1}\to -\infty$ and leading-order terms of (5.8) then determine
\begin{equation} \partial_\zeta U_{1}={-}r\left({-}2c\zeta\right)^{1/N}, \end{equation}which, integrated, reveals that the transition into the boundary layer occurs when
\begin{equation} \delta^{1/2}\zeta^{({N+1})/{N}}=O(1), \end{equation}or, equivalently, when
\begin{equation} \lvert \theta-\theta_Y \rvert=O(\delta^{(({N+2})/{2(N+1)})})=O({\textit{Bi}}^{-{1}/({N+1})}). \end{equation}
Thus, the viscoplastic boundary layer is of thickness 
$\epsilon ={\textit {Bi}}^{-{1}/({N+1})}$ (cf. Piau & Debiane Reference Piau and Debiane2004) and, after scaling, we have the boundary layer equation (using the same notation as in the Bingham problem)
\begin{equation} \frac{\partial P}{\partial r}\equiv \frac{2c}{r}=\frac{1}{r^{N+1}}\frac{\partial}{\partial\tilde{\phi}}(\lvert\partial_{\tilde{\phi}}\tilde{U}\rvert^{N-1}\partial_{\tilde{\phi}}\tilde{U}), \end{equation}which integrates to
\begin{equation} \tilde{U}=\frac{N}{N+1}\left(2c\right)^{{1}/{N}} r[(\varPhi-\tilde{\phi})^{({N+1})/{N}}-\varPhi^{({N+1})/{N}}], \end{equation}
and matching to leading order at 
$\tilde {\phi }=\varPhi$ implies that
\begin{equation} \varPhi(r)=\hat{\varPhi}r^{-{2N}/({N+1})}, \end{equation}where
\begin{equation} \hat{\varPhi}=\frac{1}{c}\left(\frac{\left(N+1\right)A}{N2^{1/N}}\right)^{{N}/({N+1})}.\end{equation}The additional shear stress due to the viscoplastic boundary layer is given, to leading order, by
\begin{align} \left.\left(\frac{1}{r^N}\left|\frac{\partial u}{\partial \theta}\right|^{N-1}\frac{\partial u}{\partial \theta}\right)\right|_{\theta=\alpha} &=\left.-\left(\frac{1}{\epsilon^Nr^N}\left|\partial_{\tilde{\phi}} \tilde{U}\right|^{N-1}\partial_{\tilde{\phi}} \tilde{U}\right)\right|_{\tilde{\phi}=0}\nonumber\\ &=2{\textit{Bi}}^{{N}/({N+1})}c\varPhi ={\textit{Bi}}^{{N}/({N+1})}\left(\frac{2\left(N+1\right)A}{Nr^2}\right)^{{N}/({N+1})}. \end{align}Using incompressibility, we find the polar velocity at the top of the boundary layer
\begin{align} v(r,\varTheta=1)=\epsilon \tilde{V}(r,\tilde{\phi}=\varPhi)={-}\epsilon r^{-({3N+1})/({N+1})}\frac{2N^3}{(2N+1)(N+1)^2}\hat{\varPhi}^{({2N+1})/{N}} (2c)^{1/N}. \end{align} We note that the radial dependence of both the boundary layer width and the first-order corrections to the velocity field are consistent with expansions in 
$r^{2N}{\textit {Bi}}$, as anticipated. Unlike for a Bingham fluid, the Cartesian thickness of the boundary layer, 
$r\epsilon \varPhi (r)=\epsilon \hat {\varPhi }r^{(1-N)/(1+N)}$, is dependent on 
$r$ – growing with 
$r$ for a shear-thinning fluid and decreasing for a shear-thickening fluid. This can be understood by considering the magnitude of the strain rate in the bulk, which decreases with 
$r$, resulting in a higher viscosity for a shear-thinning fluid, and effectively reducing the local Bingham number (the ratio of plastic and viscous stresses) relative to the global Bingham number, 
${\textit {Bi}}$. Since the boundary layer thickness increases with decreasing Bingham number, we find a larger boundary layer thickness at larger distances from the apex for a shear-thinning fluid (and the reverse for a shear-thickening fluid).
Provided 
$\epsilon \gg {\textit {Bi}}^{-1}$, as is the case for 
$N>0$, the first-order correction to the bulk solution is found at an order before viscous stresses need to be considered. Hence, we can proceed to solve for these corrections as in § 3.5. The dependent variables in the bulk are given by the following series expansions:
\begin{gather} \psi=\psi_0(\varTheta)+\epsilon r^{-{2N}/({N+1})}\hat{\psi}_1(\varTheta)+\ldots, \end{gather}
\begin{gather}u=\frac{\hat{u}_0(\varTheta)}{r}+\epsilon r^{-({3N+1})/({N+1})} \hat{u}_1(\varTheta)+\ldots, \end{gather}
\begin{gather}v=\epsilon r^{-{(3N+1)}/({N+1})} \hat{v}_1(\varTheta)+\ldots, \end{gather}
and following through the analysis with the new expression for 
$\varPhi$ gives a similar set of ODEs to those for the Bingham case, (3.46)–(3.48). These equations, along with the asymptotic behaviour at 
$\varTheta =1$, are detailed in Appendix B. They are numerically integrated using a shooting technique as in § 3.5.
Example profiles of the first-order corrections, for 
$\alpha ={\rm \pi} /4$ and different values of 
$N$, are given in figure 11. The general dependence of these profiles on 
$N$ can be explained as follows. For a shear-thinning (thickening) fluid, the high shear at the walls results in a lower (higher) effective viscosity, allowing the pressure gradient to support higher (lower) shear rates and a thinner (thicker) boundary layer, resulting in a smaller (larger) adjustment to the leading-order plastic velocities and stresses in the bulk of the wedge. The velocity in the radial direction depends on 
$N$ in a non-trivial way. For a small flow index, for which the constitutive law is almost plastic, the profile is largely independent of 
$\varTheta$, corresponding to a reasonably plastic-like velocity profile. As the index increases the profile becomes more strongly sheared and, above some value of the flow index, the magnitude of the radial flow enhancement also reduces.
Figure 11. First-order corrections to the velocity and stress fields in the bulk as functions of the rescaled angle, 
$\varTheta$, for a Herschel–Bulkley fluid with flow index 
$N$ (shown in legend). A solid line depicts the results for the Bingham model, 
$N=1$. All solutions are for 
$\alpha ={\rm \pi} /4$.
6. Bingham fluid in a cone
In this section we examine converging flow of a Bingham fluid through an axisymmetric cone. The results obtained for the two-dimensional flow through a wedge are generalised to this three-dimensional setting. Importantly, we find that viscoplasticity induces angular velocity away from the boundaries towards the axis of the cone and enhances boundary shear stress, and, when the material is only weakly yielded in the bulk, viscoplastic boundary layers develop. The conical geometry, however, introduces some algebraic complexity to the governing equations and asymptotic solution.
In the axisymmetric conical geometry we employ spherical coordinates, 
$(r,\theta ,\phi )$, where 
$r$ is the distance from the apparent vertex of the cone, 
$\theta \in [0,{\rm \pi} ]$ is now the polar angle measured from the axis of symmetry of the cone, and 
$\phi \in [0,2{\rm \pi} )$ is the azimuthal angle measured from an arbitrary choice of axis perpendicular to the axis of symmetry. The rigid walls are located at 
$\theta =\alpha$, the half-angle of the cone, and we assume independence of 
$\phi$ and no flow in the 
$\phi$ direction. Thus, we define the velocity components as 
$\boldsymbol {u}=(u,v,0)$.
This problem is characterised by the following dimensional parameters (and for clarity of exposition we focus only on the Bingham rheology): density, 
$\rho$, viscosity, 
$\mu$, yield stress, 
$\tau _c$, typical radial distance from the apex, 
$R$, and the volume flux, 
$\mathcal {Q}$. (We note that, for flow in a cone, the volume flux, 
$\mathcal {Q}$, is dimensionally distinct from the volume flux per unit width in a wedge, 
$Q$). The flow variables are rendered dimensionless as follows: velocities are scaled by 
$\mathcal {Q}/R^2$, strain rates by 
$\mathcal {Q}/R^3$ and stresses and pressures by 
$\mu \mathcal {Q}/R^3$. The resulting non-dimensional equations are
\begin{gather} \frac{1}{r^2}\frac{\partial}{\partial r}(r^2u)+\frac{1}{r\sin\theta}\frac{\partial}{\partial\theta}(v\sin\theta)=0, \end{gather}
\begin{gather}\frac{\partial p}{\partial r}=\frac{\partial \tau_{rr}}{\partial r}+ \frac{1}{r}\frac{\partial \tau_{r\theta}}{\partial \theta}+\frac{3}{r}\tau_{rr}+\frac{\cot\theta}{r}\tau_{r\theta}, \end{gather}
\begin{gather}\frac{1}{r}\frac{\partial p}{\partial\theta}=\frac{\partial \tau_{r\theta}}{\partial r}+\frac{1}{r}\frac{\partial\tau_{\theta\theta}}{\partial \theta}+\frac{3}{r}\tau_{r\theta}+\frac{\cot\theta}{r}(\tau_{\theta\theta}-\tau_{\phi\phi}), \end{gather}
\begin{gather}\boldsymbol{\tau}=\left(1+\frac{{\textit{Bi}}}{\dot{\gamma}}\right)\boldsymbol{\dot{\gamma}}, \end{gather}
\begin{gather} \dot{\gamma}_{rr}=2\frac{\partial u}{\partial r}, \quad \dot{\gamma}_{r\theta}=\frac{1}{r}\frac{\partial u}{\partial \theta}+\frac{\partial v}{\partial r}-\frac{v}{r}, \quad \dot{\gamma}_{\theta\theta}=\frac{2}{r}\frac{\partial v}{\partial \theta}+\frac{2u}{r}, \quad \dot{\gamma}_{\phi\phi}=\frac{2v\cot\theta}{r}+\frac{2u}{r}, \end{gather}
\begin{gather} \dot{\gamma}=\sqrt{\frac{1}{2}(\dot{\gamma}_{rr}^2+\dot{\gamma}_{\theta\theta}^2+\dot{\gamma}_{\phi\phi}^2+2\dot{\gamma}_{r\theta}^2)}, \end{gather}representing incompressibility (6.1), conservation of momentum in the radial and polar directions (6.2) and (6.3), the Bingham constitutive law (6.4) and the definition of the non-zero components and second invariant of the strain-rate tensor (6.5a–d) and (6.6). The Bingham number for this problem is given by
\begin{equation} {\textit{Bi}}=\frac{R^3\tau_c}{\mu \mathcal{Q}}. \end{equation}
Again, due to this dependence of the Bingham number on 
$R$, we anticipate the expansions to proceed as functions of 
$r^3{\textit {Bi}}$. Now we have the critical length scale, 
$r_c=(\mu \mathcal {Q}/\tau _c)^{1/3}$. The plastically dominated regime, 
${\textit {Bi}}\gg 1$, corresponds to the dimensional distance from the apex of the cone being significantly larger than 
$r_c$, and the viscously dominated regime, 
${\textit {Bi}}\ll 1$, occurs for radial distances much smaller than 
$r_c$. The Reynolds number is given by 
$Re=\rho \mathcal {Q}/(\mu R)$, and again for the following it will be sufficient that 
$Re=O(1)$ and 
$Re\ll {\textit {Bi}}$, to neglect inertial terms in the 
${\textit {Bi}}\gg 1$ and 
${\textit {Bi}}\ll 1$ regimes, respectively. We note that, unlike the planar problem, the Reynolds number also has radial dependence and so, for given material parameters and volume flux, it will not be possible to neglect inertial terms for arbitrarily small distances from the apex of the cone.
The equations (6.1)–(6.3) are to be solved subject to boundary conditions
\begin{equation} u=v=0\ \text{at } \theta=\alpha, \quad v=\frac{\partial u}{\partial \theta}= 0 \text{ at } \theta=0, \end{equation}corresponding to no slip and axisymmetry, respectively. We additionally have the integral expression for the volume flux,
\begin{equation} 2{\rm \pi}\int_0^\alpha r^2u\sin\theta \,\textrm{d}\theta={-}1. \end{equation}The corresponding problem for a rigid plastic was solved by Shield (Reference Shield1955), assuming purely radial flow and a radially independent deviatoric stress state. The solution corresponds closely to Nadai's solution for a planar wedge, and to make the correspondence clearer we adjust Shield's notation slightly by parameterising the non-zero components of the deviatoric stress via
\begin{align} \tau_{rr}=\frac{2}{\sqrt{3}}{\textit{Bi}}\,\cos2\psi, \quad \tau_{\theta\theta}=\tau_{\phi\phi}={-}\frac{1}{\sqrt{3}}{\textit{Bi}}\cos2\psi, \quad \tau_{r\theta}={\textit{Bi}}\sin2\psi, \quad \psi=\psi(\theta), \end{align}
which ensures that the deviatoric stress tensor is trace free and everywhere has constant second invariant, 
${\textit {Bi}}$. The equality of 
$\tau _{\theta \theta }$ and 
$\tau _{\phi \phi }$ is a consequence of the purely radial flow.
The solution is then given by
\begin{gather} u={-}\frac{\mathcal{A}}{r^2}\exp\left({-}2\sqrt{3}\int_0^\theta\tan2\psi \,\textrm{d}\theta\right), \quad \frac{\partial p}{\partial r}=\frac{2\mathscr{c}{\textit{Bi}}}{r}, \end{gather}
\begin{gather} \psi'(\theta)=\mathscr{c}\sec2\psi-\sqrt{3}-\tfrac{1}{2}\cot\theta\tan2\psi, \end{gather}
where 
$\mathscr {c}$ is determined by the boundary conditions
\begin{equation} \psi(0)=0, \quad\psi(\alpha)=\frac{\rm \pi}{4}, \end{equation}
and 
$\mathcal {A}$ is determined by the flux condition (6.9).
The viscous problem of purely radial converging flow in a cone was first solved by Harrison (Reference Harrison1919), and later explored in great detail by Ackerberg (Reference Ackerberg1965), who found that the inclusion of inertial terms resulted in non-radial solutions occurring at the outflow of the cone. For vanishing Reynolds number, as considered by Harrison (Reference Harrison1919), the problem is more straightforward, with a purely radial solution given by
\begin{align} u=\frac{3}{2{\rm \pi} r^2}\frac{\cos^2\alpha-\cos^2\theta}{\left(1-\cos\alpha\right)^2\left(1+2\cos\alpha\right)}, \quad p=\frac{1}{{\rm \pi} r^3}\frac{1-3\cos^2\theta}{\left(1-\cos\alpha\right)^2\left(1+2\cos\alpha\right)}+\text{const}. \end{align}As with the planar wedge, the viscous velocity profile (6.14a,b) exhibits enhanced velocity in the centre of the wedge, compared to the perfectly plastic velocity profile (6.11a,b). As the fluid flows from large to small distances from the apex, the strain rate increases due to the converging nature of the flow, and the material evolves from plastically dominated to viscously dominated behaviour. Thus, superimposed upon the radial dependence of the flow due to the converging geometry, we expect an additional acceleration at the centre of the cone, and velocity reduction at the outer surface of the cone. Conservation of mass then requires a flow of the fluid in the negative polar direction, towards the centre of the cone and away from the walls, and the quantification of this velocity field is a significant outcome of this analysis.
6.1. The plastic regime
The plastic solution (6.11a,b) and (6.12) does not satisfy no slip on the boundary of the cone, and the strain rate becomes unbounded here. So, for a viscoplastic fluid in the plastic regime, 
${\textit {Bi}}\gg 1$ (or equivalently, sufficiently far from the apex), we construct a viscoplastic boundary layer and intermediate layer, as for the wedge solution. Since the boundary layer and intermediate region are relatively thin, the curvature in the 
$\phi$ direction is negligible and the equations and solutions are identical to the planar geometry in these regions. In particular, we have the same boundary layer thickness, 
$\epsilon ={\textit {Bi}}^{-1/2}$, and intermediate region thickness, 
$\delta ={\textit {Bi}}^{-2/3}$. On the other hand, the 
$r$ dependent terms, 
$\varPhi$, 
$U_0$ and 
$V_0$, which are determined from matching with the bulk solution, differ.
In the bulk, we use the same strained coordinate approach, defining
\begin{equation} \theta=(\alpha-\epsilon\varPhi)\varTheta, \end{equation}and velocity expansions
\begin{equation} u=u_0+\epsilon u_1+\ldots, \quad v=\epsilon v_1+\ldots. \end{equation}
However, we need to introduce a different parameterisation of the deviatoric stress, which allows for a weak flow in the polar direction, and hence a non-zero normal stress difference, 
$\tau _{\theta \theta }-\tau _{\phi \phi }$. One such parameterisation is
\begin{gather} \frac{\tau_{rr}}{{\textit{Bi}}}=\frac{2}{\sqrt{3}}\cos2\psi\cos2\chi, \quad \frac{\tau_{\theta\theta}}{{\textit{Bi}}}={-}\frac{1}{\sqrt{3}}\cos2\psi\cos2\chi+\cos2\psi\sin2\chi, \end{gather}
\begin{gather} \frac{\tau_{\phi\phi}}{{\textit{Bi}}}={-}\frac{1}{\sqrt{3}}\cos2\psi\cos2\chi-\cos2\psi\sin2\chi, \quad \frac{\tau_{r\theta}}{{\textit{Bi}}}=\sin2\psi, \end{gather}
which reduces to the previous parameterisation (6.10a–d) for 
$\chi =0$ and gives all possible symmetric, trace-free stress states satisfying 
$\tau \equiv \sqrt {\tau _{ij}\tau _{ij}/2}={\textit {Bi}}$, and 
$\tau _{r\phi }=\tau _{\theta \phi }=0$. We define asymptotic series for 
$\psi$ and 
$\chi$ by
\begin{equation} \psi=\psi_0(r)+\epsilon\psi_1(r,\theta)+\ldots, \quad \chi = \epsilon\chi_1(r,\theta)+\ldots, \end{equation}
so that, at 
$O(1)$, 
$\chi =0$ and the stress decomposition is precisely that given in (6.10a–d), for the purely radial solution.
With the coordinate transform (6.15), the full governing equations (6.1)–(6.6) can be reduced to
\begin{align} &\frac{1}{r}L_r(r^2 u)+L_\varTheta v +v\cot(\alpha\varTheta-\epsilon\varPhi\varTheta)=0, \end{align}
\begin{align} &(L_\varTheta L_r+3L_\varTheta)\tau_{rr}+(L_\varTheta^2-L_r^2-3L_r)\tau_{r\theta}-L_rL_\varTheta\tau_{\theta\theta} +L_\varTheta(\cot(\alpha\varTheta-\epsilon\varPhi\varTheta)\tau_{r\theta})\nonumber\\ &\quad -L_r(\cot(\alpha\varTheta-\epsilon\varPhi\varTheta)(\tau_{\theta\theta}-\tau_{\phi\phi}))=0, \end{align}
\begin{gather} 2\tau_{r\theta}L_r u= \tau_{rr}(L_\varTheta u+L_r v - v), \end{gather}
\begin{gather} 2\tau_{\theta\theta}L_r u = \tau_{rr}(2L_\varTheta v+2u).\end{gather}
Here, (6.20) is conservation of mass, (6.21) is the curl of the momentum balances and (6.22) and (6.23) are a statement of the proportionality between corresponding components of the deviatoric stress and strain-rate tensors. The linear operators, 
$L_r$ and 
$L_\varTheta$, are given by
\begin{equation} L_r=r\frac{\partial}{\partial r}+\frac{\epsilon r\varPhi'\varTheta}{\alpha}\frac{\partial}{\partial\varTheta}, \quad L_\varTheta= \frac{1}{\alpha-\epsilon\varPhi}\frac{\partial}{\partial\varTheta}=\left(\frac{1}{\alpha}+\frac{\epsilon \varPhi}{\alpha^2}+\ldots\right)\frac{\partial}{\partial\varTheta}. \end{equation} Expansion of the governing equations at 
$O(1)$ gives the equivalent of Shield's solution
\begin{equation} u_0:=\frac{\hat{u}_0}{r^2}={-}\frac{\mathcal{A}}{r^2}\exp\left({-}2\sqrt{3}\int_0^{\varTheta}\alpha\tan2\psi_0 \, \textrm{d}\varTheta\right), \quad \frac{\partial p_0}{\partial r}=\frac{2c{\textit{Bi}}}{r}, \end{equation}
\begin{equation} \psi_0'(\varTheta)=\alpha(\mathscr{c}\sec2\psi_0-\sqrt{3}-\tfrac{1}{2}\cot\alpha\varTheta\tan2\psi_0), \end{equation}
where 
$\mathscr {c}$ is determined by the boundary conditions
\begin{equation} \psi_0(0)=0, \quad\psi_0(1)=\frac{\rm \pi}{4}, \end{equation}
and 
$\mathcal {A}$ is determined by the flux condition
\begin{equation} 2{\rm \pi}\alpha \int_0^1 \hat{u}_0(\varTheta)\sin\alpha\varTheta \,\textrm{d}\varTheta=1.\end{equation}
We define 
$\hat {U}_0=\hat {u}_0(\varTheta =1)$ since, unlike for the planar solution, this does not have a closed-form expression in terms of 
$\mathcal {A}$ and 
$\mathscr {c}$. We note that this solution (6.25a,b)–(6.28) exhibits a non-vanishing radial velocity and a divergent shear rate as 
$\varTheta \to 1$, just as for the flow through a planar wedge. We therefore must introduce both the intermediate and boundary layers in order to enforce the boundary conditions. However, this matching is identical to §§ 3.3 and 3.4 and we deduce that
\begin{equation} \varPhi(r)\equiv\frac{\hat{\varPhi}}{r^{3/2}}=\sqrt{-\frac{\hat{U}_0}{\mathscr{c}r^3}}, \quad V_0(r)={-}\frac{\mathscr{c}r\varPhi^3}{3}={-}\frac{\mathscr{c}\hat{\varPhi}^3}{3}r^{{-}7/2}. \end{equation}
Thus, the boundary layer width scales like 
$(r^3{\textit {Bi}})^{-1/2}$ and the first-order corrections to the velocities scale like 
$r^{-2}(r^3{\textit {Bi}})^{-1/2}$, as expected given the true expansion parameter 
$r^3{\textit {Bi}}$. We may now use these expressions to calculate the effect of the boundary layer in the bulk, and, in particular, the induced angular velocity. To do so we define
\begin{align} u_1=r^{{-}7/2}\hat{u}_1(\varTheta), \quad v_1=r^{{-}7/2}\hat{v}_1(\varTheta), \quad \psi_1=r^{{-}3/2}\hat{\psi}_1(\varTheta), \quad \chi_1=r^{{-}3/2}\hat{\chi}_1(\varTheta), \end{align}
and expand the governing equations (6.20)–(6.23) at order 
$\epsilon$ to obtain four ODEs for 
$\hat {u}_1$, 
$\hat {v}_1$, 
$\hat {\psi }_1$ and 
$\hat {\chi }_1$:
\begin{align} &{-}\frac{3}{2}\hat{u}_1+\frac{1}{\alpha}\hat{v}_1'+\hat{v}_1\cot\alpha\varTheta=\frac{3\hat{\varPhi}\varTheta}{2\alpha}\hat{u}_0', \end{align}
\begin{align}
&\frac{2}{\alpha^2}(\hat{\psi}_1\cos2\psi_0)''+\frac{3}{\alpha}(\hat{\chi}_1\cos2\psi_0)'
-\frac{\sqrt{3}}{\alpha}(\hat{\psi}_1\sin2\psi_0)'+\frac{2}{\alpha}(\hat{\psi}_1\cot\alpha\varTheta
\cos2\psi_0)'\nonumber\\ &\quad
+\left(\frac{9}{2}\hat{\psi}_1+6\hat{\chi}_1\cot\alpha\varTheta\right)\cos2\psi_0
=\frac{\hat{\varPhi}}{\alpha^2}\left(\frac{3\sqrt{3}\varTheta}{2}(\cos2\psi_0)''-\sqrt{3}(\cos2\psi_0)'
\right.\nonumber\\ &\quad
- \left. \frac{2}{\alpha}(\sin2\psi_0)''-\frac{9\alpha\varTheta}{4}(\sin2\psi_0)'-((\cot\alpha\varTheta+\alpha\varTheta\textrm{cosec}^2\alpha\varTheta)\sin2\psi_0)'\right),
\end{align}
\begin{gather} -7\sqrt{3}\hat{u}_1\tan2\psi_0-\frac{2}{\alpha}\hat{u}_1'+9\hat{v}_1-8\sqrt{3}\hat{\psi}_1\hat{u}_0+\frac{4}{\alpha}\hat{\psi}_1\hat{u}_0'\tan2\psi_0\nonumber\\ =\frac{\hat{\varPhi}}{\alpha^2}(3\sqrt{3}\alpha\varTheta\tan2\psi_0+2)\hat{u}_0', \end{gather}
\begin{gather} 3\hat{u}_1-\frac{4}{\alpha}\hat{v}_1'-8\sqrt{3}\hat{\chi}_1\hat{u}_0={-}\frac{3\hat{\varPhi}\varTheta}{\alpha}\hat{u}_0', \end{gather}with boundary conditions
\begin{equation} \hat{\psi}_1=\hat{v}_1=0 \text{ at } \varTheta=0, \quad \hat{v}_1={-}\frac{\mathscr{c}\hat{\varPhi}^3}{3} \text{ at } \varTheta=1, \end{equation}
from symmetry at 
$\varTheta =0$ and matching to the boundary layer solution at 
$\varTheta =1$. As in the planar case, matching of 
$\tau _{rr}$ to the intermediate layer requires
\begin{equation} \hat{\psi}_1=\frac{3\sqrt{3}\hat{\varPhi}}{4}\quad \text{at}\ \varTheta=1. \end{equation}
Similarly, analysis of the normal stress difference, 
$\tau _{\theta \theta }-\tau _{\phi \phi }$, in the intermediate layer gives
\begin{align} \tau_{\theta\theta}-\tau_{\phi\phi} &\equiv\left(1+\frac{{\textit{Bi}}}{\dot{\gamma}}\right)\left(\frac{2}{r}\frac{\partial v}{\partial \theta}-\frac{2v}{r}\cot\theta\right)\nonumber\\ &=\left(-\frac{{\textit{Bi}}}{\dfrac{1}{r\sqrt{\delta}}\dfrac{\partial U_1}{\partial \zeta}}+\ldots\right)\left(\frac{2\epsilon r\varPhi'}{r \sqrt{\delta}}\frac{\partial U_1}{\partial \zeta}+\ldots\right) =3{\textit{Bi}}^{1/2}\varPhi+\ldots, \end{align}
using (3.20), (3.28) and the identity 
$r\varPhi '=-3\varPhi /2$. While in the bulk solution
\begin{equation} \tau_{\theta\theta}-\tau_{\phi\phi}\equiv2{\textit{Bi}}\cos2\psi\sin2\chi=4{\textit{Bi}}^{1/2}\chi_1\sqrt{2\alpha\left(2c-\cot\alpha\right)\left(1-\varTheta\right)}+\ldots,\end{equation}
as 
$\varTheta \to 1$, using the behaviour of 
$\psi _0$ in the neighbourhood of 
$\varTheta =1$ (see Appendix C). Thus, 
$\hat {\chi }_1$ diverges like 
$(1-\varTheta )^{-1/2}$ as 
$\varTheta \to 1$, which is consistent with the leading-order behaviour of 
$\hat {\chi }_1$ from the differential equations (6.31)–(6.34) close to 
$\varTheta =1$. Although the divergence of 
$\hat {\chi }_1$ may seem problematic for the asymptotic order of our expansions, we expect the plastic regime to break down within the intermediate layer, for which 
$1-\varTheta =O(\delta )$. Here, we have 
$\chi =O(\epsilon \delta ^{-1/2})\ll 1$, since 
$\delta \gg \epsilon ^2$, and so 
$\chi$ remains small throughout the region of validity of the plastic equations. The reason for this divergence of the normal stress difference is that the gradient of the polar velocity, 
$\partial v/\partial \theta$, becomes large in the boundary layer, while the polar velocity itself remains small, thus 
$\tau _{\theta \theta }$ becomes significantly larger than 
$\tau _{\phi \phi }$.
The boundary-value problem given by (6.31)–(6.36) was again solved by a shooting method. In this situation we have potential singularities at both ends of the domain, due to the 
$\cot \alpha \varTheta$ terms, so asymptotic analysis was required to step a small distance away from 
$\varTheta =0$ and 
$\varTheta =1$ (see (C4)–(C8) in Appendix C) before numerically shooting towards the centre of the domain, and fixing unknown coefficients by enforcing continuity at some central point. Again, the solutions were found to be essentially identical for sufficiently small choices of the step sizes, and the method converges effectively for any sufficiently central choice of matching point. The profiles given in figure 12 were produced with a step size of 
$d=10^{-8}$ away from the singularities at both ends of the domain, and with matching point at 
$\varTheta =0.5$.
Figure 12. The first-order corrections to the velocities, 
$u$ and 
$v$, and the stress orientation functions, 
$\psi$ and 
$\chi$, for viscoplastic flow in a cone in the plastic regime, as functions of the scaled polar coordinate 
$\varTheta$, for 
$\alpha ={\rm \pi} /6$ (solid), 
${\rm \pi} /4$ (dashed), 
${\rm \pi} /3$ (dash-dotted) and 
${\rm \pi} /2$ (dotted).
From the solutions to the boundary-value problem, we find that the velocity profiles are qualitatively similar to those for the planar wedge (figure 3), with a negative polar velocity out from the boundary layer towards the centre of the cone and an enhancement of the velocity in the radial direction, both of which are increased for decreasing half-angle, due to the boundary layer thickness, 
$\hat {\varPhi }$, increasing for decreasing 
$\alpha$. The shear stress is again enhanced throughout the domain, due to the greater shear stress at the wall, and we find that the normal stress difference becomes significant at the edge of the boundary layer, 
$\varTheta =1$, via the divergence of 
$\hat {\chi }_1$, as discussed above.
6.2. The viscous regime
The viscous regime, 
${\textit {Bi}}\ll 1$, can be tackled by expanding the governing equations around the leading-order Stokes solution. We define the asymptotic series
\begin{equation} u=u_0+{\textit{Bi}}\, u_1+\ldots,\quad v={\textit{Bi}}\, v_1+\ldots, \quad p=p_0+{\textit{Bi}} p_1+\ldots, \end{equation}
where 
$u_0$ and 
$p_0$ are given by (6.14a,b). For homogeneity in 
$r$ in the perturbation to the stress tensor, and using conservation of mass, we can write
\begin{equation} u_1\equiv r\hat{u}_1(\theta)=\frac{r}{\sin\theta}g'(\theta), \quad v_1\equiv r\hat{v}_1(\theta)={-}\frac{3}{\sin\theta}g(\theta), \end{equation}
where 
$g$ is a function to be determined, and 
$'$ represents differentiation with respect to 
$\theta$. From the expression of conservation of momentum in the angular direction at 
$O({\textit {Bi}})$, we find that the angular pressure gradient is independent of 
$r$, and so the radial pressure gradient is independent of 
$\theta$, giving
\begin{equation} r\frac{\partial p_1}{\partial r}=\mathcal{H}(\theta)+\frac{1}{\sin\theta}g'''-\frac{\cos\theta}{\sin^2\theta} g''+\frac{1+6\sin^2\theta}{\sin^3\theta}g'=\mathcal{C}, \end{equation}
where 
$\mathcal {C}$ is an arbitrary constant, to be determined, and 
$\mathcal {H}(\theta )$ is given by
\begin{align} \mathcal{H}(\theta)=\frac{\textrm{d}}{\textrm{d}\theta}\left(\frac{\sin2\theta}{\sqrt{\sin^22\theta+3\left(\cos2\theta-\cos2\alpha\right)^2}}\right)+\frac{2\cos^2\theta+6\cos2\theta-6\cos2\alpha}{\sqrt{\sin^22\theta+3\left(\cos2\theta-\cos2\alpha\right)^2}}. \end{align}Thus we have the third-order ODE,
\begin{equation} g'''-\cot\theta g''+(6+\textrm{cosec}^2\theta)g'=(\mathcal{C}-\mathcal{H}(\theta))\sin\theta. \end{equation}
For 
$\hat {u}_1(0)$ to remain finite, we need 
$g(0)=0$. For symmetry we require 
$g$ to be an even function of 
$\theta$, which ensures 
$u_1$ is even and 
$v_1$ is odd. Finally, we have the no-slip boundary condition at the wall, which gives
\begin{equation} g(\alpha)=g'(\alpha)=0. \end{equation} To solve the ODE numerically we shoot from 
$\theta =d\ll 1$, to avoid the potential singularity at 
$\theta =0$, with initial conditions
\begin{equation} g(d)=\lambda d^2, \quad g'(d)=2\lambda d, \quad g''(d)=2\lambda, \end{equation}
and determine 
$\lambda$ and 
$\mathcal {C}$ by enforcing the two boundary conditions at 
$\theta =\alpha$, (6.44). The solutions are found to be essentially independent of 
$d$ for 
$d<10^{-5}$. Profiles of 
$\hat {u}_1$ and 
$\hat {v}_1$ are given in figure 13, demonstrating the anticipated negative polar velocity away from the boundary and towards the centre of the cone. As in the planar case, the first-order corrections are larger for larger 
$\alpha$, due to the magnitude of the deviatoric stress in the leading-order viscous solution being smaller at larger half-angles, resulting in a less yielded fluid and more significant plasticity effects.
Figure 13. The first-order corrections to the velocities, 
$u$ and 
$v$, for viscoplastic flow in a cone in the viscous regime, as functions of the scaled polar coordinate 
$\theta /\alpha$, for 
$\alpha ={\rm \pi} /6$ (solid), 
${\rm \pi} /4$ (dashed), 
${\rm \pi} /3$ (dash-dotted) and 
${\rm \pi} /2$ (dotted).
The additional shear stress at the wall is given by
\begin{equation} {\textit{Bi}}\,\left.\tau_{r\theta}^{(1)}\right|_{\theta=\alpha}={\textit{Bi}}+{\textit{Bi}}\left.\left(\frac{1}{r}\frac{\partial u_1}{\partial \theta}+\frac{\partial v_1}{\partial r}-\frac{v_1}{r}+1\right)\right|_{\theta=\alpha}={\textit{Bi}}\left(\frac{g''(\alpha)}{\sin\alpha}+1\right), \end{equation}
which is plotted as a function of 
$\alpha$ in figure 14a), demonstrating a similar behaviour to that for the planar wedge, with a minimum attained for a slightly smaller half-angle of 
$\alpha \approx 27^{\circ }$. The constant 
$\mathcal {C}$ is also of interest since the first-order correction to the radial pressure gradient is given by 
${\textit {Bi}}\,\mathcal {C}/r$, thus the behaviour of 
$\mathcal {C}=r\partial p_1/\partial r$ with 
$\alpha$ is given in figure 14b), showing that 
$\mathcal {C}$ varies from 4.20 (approximately) at 
$\alpha ={\rm \pi} /2$ and diverges as 
$\alpha \to 0$. In particular, we find that the perturbation to the shear stress at the wall, 
$\tau _{r\theta }^{(1)}\sim 4/3$, and 
$\mathcal {C}\sim 8/(3\alpha )$ for 
$\alpha \ll 1$ as required for consistency with the pipe Poiseuille solution in the limit 
$\alpha \to 0$ with 
$r\alpha =1$ fixed (see Appendix D).
Figure 14. Features of the perturbation solution for the conical converging flow of a viscoplastic in the low Bingham number regime: (a) the additional wall shear stress, 
$\tau _{r\theta }^{(1)}|_{\theta =\alpha }$, and (b) the scaled additional radial pressure gradient, 
$r\partial p_1/\partial r=\mathcal {C}$, as functions of cone half-angle, 
$\alpha$.
7. Discussion and conclusions
Asymptotic solutions have been found for the converging flow of a viscoplastic fluid through a planar wedge and conical geometry, in the plastically and viscously dominated regimes, and verified for the case of a Bingham fluid in a wedge using direct finite element simulations. A key feature in both regimes is that no purely radial solution is possible. Instead, a weak angular flow is induced away from the boundaries and towards the centre of the domains. This result can be explained as a consequence of mass conservation as the fluid flows from a plastically dominated to viscously dominated velocity profile, experiencing enhanced flow velocity at the centre of the domain and reduced velocity near the walls as a consequence of no-slip conditions.
In the plastically dominated regime (
${\textit {Bi}}\gg 1$) for a Bingham fluid, thin viscoplastic boundary layers of width 
$\epsilon ={\textit {Bi}}^{-1/2}$ are required at the walls, allowing the flow to satisfy no slip. The angular extent of these boundary layers also depends on the radial distance from the apparent apex, 
$r$, scaling as 
$\epsilon r^{-1}$ in a wedge and 
$\epsilon r^{-3/2}$ in a cone. Additionally, an intermediate asymptotic layer is required to regularise divergent strain rates close to the wall. The weak angular flow, induced by the boundary layer, is also 
$O(\epsilon )$ and both the boundary layer and angular flow are weaker for larger half-angles. The shear stress is enhanced compared to the rigid plastic solution, exceeding the yield stress by 
$O(\epsilon {\textit {Bi}})$ at the walls.
For a Herschel–Bulkley fluid of flow index 
$N$, the same features are observed in the plastically dominated regime, with the distinguished scaling given instead by 
$\epsilon ={\textit {Bi}}^{-1/(N+1)}$ and the radial dependence of the boundary layer width given by 
$\epsilon r^{-2N/(N+1)}$ for flow in a wedge. We choose not to explore in detail, in this paper, the problem of a Herschel–Bulkley fluid in an axisymmetric conical geometry in the plastically dominated regime (
${\textit {Bi}}\gg 1$). However, we can deduce by simple scaling arguments, analogous to those in §§ 5, 6, that the boundary layer width scales as 
$(r^{3N}{\textit {Bi}})^{-1/(N+1)}$, where 
$N$ is the flow index. Notably, the boundary layer is of constant Cartesian thickness for the flow of a Bingham material through a wedge and for a Herschel–Bulkley fluid of flow index 
$N=1/2$ through a cone, but not for the other situations.
In the viscously dominated regime (
${\textit {Bi}}\ll 1$), the weak angular flow is 
$O({\textit {Bi}})$ and, in contrast to the plastically dominated regime, is stronger for larger half-angles. This is due to the leading-order viscous flow shearing the fluid less strongly for larger half-angles, resulting in less strongly yielded fluid and a more significant effect of the yield stress. The shear stress at the wall is again enhanced compared to the purely viscous solution, with the excess shear stress being 
$O({\textit {Bi}})$.
The direct finite element simulations, carried out in FEniCS (Alnæs et al. Reference Alnæs, Blechta, Hake, Johansson, Kehlet, Logg, Richardson, Ring, Rognes and Wells2015) using a combination of augmented-Lagrangian and Newton methods, strongly supported the validity of the asymptotic solutions derived in the case of a Bingham fluid in a planar wedge. Boundary layer widths and velocity profiles were found to be in excellent agreement with the theoretical results. The challenges posed by the converging geometry and nonlinear constitutive equation leave room for a more expansive and accurate numerical study. In particular, we were unable to reproduce accurately the predictions for both the plastic and viscous regimes in the same simulations due to the vastly disparate scales involved in the different regimes.
This study has tackled the simplest and most generic scenario for converging flows of viscoplastic materials, demonstrating the emergence of viscoplastic boundary layers when the viscous stresses are weak, and the development of an angular flow. It would be interesting to analyse the motion in related problems, such as those driven by a body force (e.g. gravity) or situations in which there is wall slip or significant inertia. There have been several experimental investigations of extrusions and flows through dies and contractions (for e.g. Wildman et al. Reference Wildman, Blackburn, Benton, McNeil and Parker1999; Jay, Magnin & Piau Reference Jay, Magnin and Piau2002; Rabideau et al. Reference Rabideau, Moucheront, Bertrand, Rodts, Roussel, Lanos and Coussot2010; Luu, Philippe & Chambon Reference Luu, Philippe and Chambon2015), but the authors are unaware of any experimental observations that may be used to validate the predictions of this theory. In particular it would be interesting to detect and quantify experimentally the emergent angular flow predicted by our analysis.
Funding
This work was funded by the Engineering and Physical Sciences Research Council, UK. (EP/R513179/1).
Declaration of interests
The authors report no conflict of interest.
Appendix A
Expansion of the equation for 
$\psi _0$, (3.11), for 
$1-\varTheta \ll 1$ gives
\begin{equation} \psi_0 = \frac{\rm \pi}{4} -\sqrt{\alpha c\left(1-\varTheta\right)} + \frac{2}{3}\alpha\left(1-\varTheta\right)+\ldots. \end{equation}
Substituting this into the definitions of the scaled radial velocity, 
$\hat {u}_0$, (3.12), and the components of deviatoric stress, 
$\tau _{rr}$ and 
$\tau _{r\theta }$, (3.2), gives
\begin{gather} \hat{u}_0={-}\frac{A}{c}-\frac{2A}{c}\sqrt{\frac{\alpha}{c}\left(1-\varTheta\right)}-\frac{8A\alpha}{3c^2}\left(1-\varTheta\right)+\ldots, \end{gather}
\begin{gather}\tau_{rr}\equiv {\textit{Bi}}\cos2\psi_0=2{\textit{Bi}}\sqrt{\alpha c\left(1-\varTheta\right)}-\tfrac{4}{3}{\textit{Bi}}\,\alpha\left(1-\varTheta\right)+\ldots, \end{gather}
\begin{gather}\tau_{r\theta}\equiv{\textit{Bi}}\sin2\psi_0={\textit{Bi}}-2{\textit{Bi}}\,\alpha c\left(1-\varTheta\right)+\ldots. \end{gather}Appendix B
For a Herschel–Bulkley fluid with flow index 
$N$, the generalisation of equations (3.46)–(3.48) is given by
\begin{gather} \frac{1}{\alpha^2}(\hat{\psi}_1\cos2\psi_0)''+\frac{4}{(N+1)^2}\hat{\psi}_1\cos2\psi_0-\frac{2(1-N)}{\alpha(1+N)}(\hat{\psi}_1\sin2\psi_0)'\nonumber\\ \hspace{1.5pc}=\frac{\hat{\varPhi}}{\alpha^2}\left(\frac{2N}{N+1}\varTheta\left(\cos2\psi_0\right)''- \frac{1-N}{N+1}\left(\cos2\psi_0\right)'-\frac{1}{\alpha}\left(\sin2\psi_0\right)''\right.\nonumber\\ \hspace{-9pc} \left. - \frac{N}{N+1}\alpha\varTheta\left(\sin2\psi_0\right)'\right), \end{gather}
\begin{gather} -\frac{2N}{N+1}\hat{u}_1+\frac{1}{\alpha}\hat{v}_1'=\frac{2N}{N+1}\frac{\hat{\varPhi}\varTheta}{\alpha}\hat{u}_0', \end{gather}
\begin{gather} -\frac{6N+2}{N+1}\hat{u}_1\tan2\psi_0-\frac{1}{\alpha}\hat{u}_1'+\frac{4N+2}{N+1}\hat{v}_1-4\hat{\psi}_1\hat{u}_0\sec^22\psi_0\nonumber\\ \quad =\frac{\hat{\varPhi}}{\alpha}\left(\frac{4N}{N+1}\varTheta\tan2\psi_0+\frac{1}{\alpha}\right)\hat{u}_0', \end{gather}with boundary conditions
\begin{equation} \hat{v}_1(0)=\hat{\psi}_0(0)=0, \quad \hat{v}_1(1)={-}\frac{2N^3}{\left(2N+1\right)\left(N+1\right)^2}\hat{\varPhi}^{2+{1/N}} \left(2c\right)^{1/N}, \end{equation}
where 
$\hat {\varPhi }$ is given by (5.15) and 
$c$ and 
$A$ are given by (3.14) and (3.15). When 
$N=1$ we have 
$\hat {\varPhi }=\sqrt {A}/c$ and the equations reduce to (3.46)–(3.48), as required.
The boundary condition for 
$\hat {\psi }_1$ at 
$\varTheta =1$ can be deduced by matching the deviatoric normal stress, 
$\tau _{rr}$, between the bulk and intermediate layer, as detailed for the Bingham case in § 3.5. This gives
\begin{equation} \hat{\psi}_1(1)=\frac{2N}{N+1}\hat{\varPhi}. \end{equation} Using (A1) and (A2), we analyse the limiting behaviour of the equations (B1)–(B3) for 
$\lvert 1-\varTheta \rvert \ll 1$, which yields the following local forms of all three dependent variables:
\begin{gather} \hat{\psi}_1=\frac{2N}{N+1}\hat{\varPhi}+E\sqrt{\frac{\alpha}{c}\left(1-\varTheta\right)}+\ldots, \end{gather}
\begin{gather}\hat{u}_1=F+\left(\frac{6N+2}{N+1}F+\frac{2A}{c^2}\left(\frac{c\hat{\varPhi}}{\alpha}-E\right)\right)\sqrt{\frac{\alpha}{c}\left(1-\varTheta\right)}+\ldots \end{gather}
\begin{gather}\hat{v}_1={-}\frac{2N^3\hat{\varPhi}^{2+{1/N}}}{\left(2N+1\right)\left(N+1\right)^2}\left(2c\right)^{1/N}- \frac{4N}{N+1}\frac{A\hat{\varPhi}}{c}\sqrt{\frac{\alpha}{c}\left(1-\varTheta\right)}+\ldots, \end{gather}
where 
$E$ and 
$F$ are arbitrary constants.
Appendix C
Expansion of the differential equation for 
$\psi _0$, (6.26), at 
$\varTheta =1$ gives the asymptotic series
\begin{equation} \psi_0=\frac{\rm \pi}{4}-\sqrt{\frac{\alpha}{2}\left(2\mathscr{c}-\cot\alpha\right)\left(1-\varTheta\right)}+\frac{2\sqrt{3}}{3}\alpha\left(1-\varTheta\right)+\ldots. \end{equation}Matching the leading-order radial velocity to the intermediate and boundary layers gives
\begin{equation} \hat{u}_0\left(\varTheta=1\right)={-}\mathscr{c}\hat{\varPhi}^2, \end{equation}and expansion of the leading-order solution, (6.25a,b), gives
\begin{equation} \hat{u}_0={-}\mathscr{c}\hat\varPhi^2\left(1+\frac{2\sqrt{6}}{\sqrt{2\mathscr{c}-\cot\alpha}}\sqrt{\alpha\left(1-\varTheta\right)}+\frac{16\alpha}{2 \mathscr{c}-\cot\alpha}\left(1-\varTheta\right)+\ldots\right). \end{equation}Then, analysis of the first-order equations, (6.31)–(6.34), and application of boundary conditions, (6.35a,b) and (6.36), reveals that solutions have the local form
\begin{gather} \hat{\psi}_1=\frac{3\sqrt{3}}4\hat{\varPhi}+\mathcal{E}\sqrt{1-\varTheta}+\ldots, \end{gather}
\begin{gather} \hat{u}_1=\mathcal{F}+\left(\frac{7\sqrt{6\alpha}}{2\sqrt{2 \mathscr{c}-\cot\alpha}}\mathcal{F}+\mathscr{c}\hat{\varPhi}^2\left(\frac{\sqrt{6}\hat{\varPhi}}{\sqrt{\alpha\left(2 \mathscr{c}-\cot\alpha\right)}}-\frac{2\sqrt{3}\mathcal{E}}{\left(2\mathscr{c}-\cot\alpha\right)^{3/2}}\right)\right)\sqrt{1-\varTheta}\nonumber\\ \hspace{-23.5pc} +\ldots, \end{gather}
\begin{gather} \hat{v}_1={-}\frac{\mathscr{c}\hat{\varPhi}^3}{3}\left(1+\frac{9\sqrt{6}}{\sqrt{2\mathscr{c}-\cot\alpha}}\sqrt{\alpha\left(1-\varTheta\right)}+\ldots\right),\end{gather}
where 
$\mathcal {E}$ and 
$\mathcal {F}$ are undetermined constants. We note that 
$\hat {\chi }_1$ is in fact determined algebraically by (6.34), so we require no boundary conditions for 
$\hat {\chi }_1$. To leading order at 
$\varTheta =1$ we have
\begin{equation} \hat{\chi}_1=\frac{3\sqrt{2}\hat{\varPhi}}{8\sqrt{\alpha\left(2\mathscr{c}-\cot\alpha\right)\left(1-\varTheta\right)}}+\ldots, \end{equation}as required for matching of the normal stress difference between the intermediate layer and bulk solution, (6.37) and (6.38).
For 
$\varTheta \ll 1$ we have the leading-order local forms
\begin{align} \hat{\psi}_0=\frac{\alpha(\mathscr{c}-\sqrt{3})}{2}\varTheta+\ldots, \quad \hat{u}_0={-}\mathcal{A}+\ldots, \quad \hat{\psi}_1=\tilde{\mathcal{E}}\varTheta+\ldots, \quad \hat{u}_1=\tilde{\mathcal{F}}+\ldots, \end{align}
where 
$\tilde {\mathcal {E}}$ and 
$\tilde {\mathcal {F}}$ are undetermined constants.
Appendix D
For small Bingham number, 
${\textit {Bi}}\ll 1$, we consider the plane Poiseuille flow of a Bingham fluid in a channel of height 
$2$ and volume flux per unit width 
$-1$, driven by pressure gradient 
$\textrm {d}p/\textrm {d} x=G$. From the volume-flux condition we can write
\begin{align} -\frac{1}{2}&=\int_0^1u\,\textrm{d} y=u_p y_c+\int_{y_c}^1 u\,\textrm{d} y={-}\int_{y_c}^1 y\frac{\textrm{d}u}{\textrm{d} y}\,\textrm{d} y={-}\int_{y_c}^1y\left(Gy-{\textit{Bi}}\right)\textrm{d} y \nonumber\\ &={-}\frac{G}{3}+\frac{1}{2}{\textit{Bi}}+\ldots, \end{align}
where 
$u_p$ is the velocity of the plug which occupies the region 
$-y_c < y < y_c$. In the above we have used integration by parts, made the substitution 
$\textrm {d}u/\textrm {d} y=Gy-{\textit {Bi}}$, using the integral of the streamwise conservation of momentum equation in the yielded region, and finally used the fact that 
$y_c=O({\textit {Bi}})$ in neglecting terms in the final expression. Thus we find
\begin{equation} G=\tfrac{3}{2}+\tfrac{3}{2}{\textit{Bi}}+\ldots, \end{equation}and, by a global force balance, the shear stress at the wall is given by
\begin{equation} \tau_{xy}=G=\tfrac{3}{2}+\tfrac{3}{2}{\textit{Bi}}+\ldots. \end{equation}
Similarly, for the pipe Poiseuille flow of a Bingham fluid along the axis of a cylinder of unit radius, with volume flux 
$-1$ driven by a pressure gradient 
$\textrm {d}p/\textrm {d}z=G$, we can write
\begin{align} -1&=2{\rm \pi}\int_0^1ru\,\textrm{d}r={\rm \pi} r_c^2u_p+2{\rm \pi}\int_{r_c}^1 ru\,\textrm{d}r={-} {\rm \pi}\int_{r_c}^1 r^2\left(\frac{1}{2}Gr-{\textit{Bi}}\right)\textrm{d}r \nonumber\\ &={-}\frac{\rm \pi}{8}G+\frac{\rm \pi}{3}{\textit{Bi}}+\ldots, \end{align}and so
\begin{equation} G=\frac{8}{\rm \pi}+\frac{8}{3}{\textit{Bi}}+\ldots, \end{equation}and, by a global force balance, the shear stress at the wall is given by
\begin{equation} \tau_{rz}=\frac{1}{2}G=\frac{4}{\rm \pi}+\frac{4}{3}{\textit{Bi}}+\ldots. \end{equation}
