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),

(2.14a–c)\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,

(2.15a–c)\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.

Figure 2.

The radial velocity field, $u$, as a function of polar angle, $\theta$, for plastic (solid line) and viscous (dashed line) converging flow in a wedge, with $\alpha ={\rm \pi} /4$ and $r=1$ (see (2.14a–c) and (2.15a–c)).

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

(3.1)\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

(3.2)\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

(3.3)\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

(3.4)\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

(3.5)\begin{equation} \left(r,\theta\right)=\left(r,\theta_Y(r)\varTheta\right), \end{equation}

noting that this gives partial derivatives:

(3.6a,b)\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

(3.7a–d)\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

(3.8)\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

(3.9)\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

(3.10)\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

(3.11)\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}

(3.12)\begin{gather}u_0={-}\frac{A}{r\left(c-\cos\left(2\psi_0\right)\right)}, \end{gather}

(3.13)\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

(3.14)\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,

(3.15)\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

(3.16)\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,

(3.17)\begin{equation} \delta\zeta=\alpha-\epsilon\varPhi(r)-\theta, \end{equation}

and so, as before, we have altered partial derivatives

(3.18a,b)\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

(3.19)\begin{equation} u_0={-}\frac{A}{rc}-\frac{2A}{rc\sqrt{c}}\sqrt{\delta\zeta}+\ldots, \end{equation}

thus we define

(3.20)\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

(3.21)\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

(3.22)\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}

(3.23)\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

(3.24)\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$,

(3.25)\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

(3.26)\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

(3.27)\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

(3.28)\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

(3.29)\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

(3.30)\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

(3.31)\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

(3.32)\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

(3.33)\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

(3.34)\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

(3.35)\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$,

(3.36)\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

(3.37)\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

(3.38)\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

(3.39)\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

(3.40)\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,

(3.41)\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

(3.42)\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}

(3.43)\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}

(3.44)\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

(3.45a–c)\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)

(3.46)\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}

(3.47)\begin{gather} -\hat{u}_1+\frac{1}{\alpha}\hat{v}_1'=\frac{\sqrt{A}\varTheta}{\alpha c}\widehat{u_0}', \end{gather}

(3.48)\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 )$

(3.49a,b)\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

(3.50)\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)

(3.51)\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

(3.52)\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

(3.53)\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

(3.54)\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

(3.55)\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}

(3.56)\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

(3.57)\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$

(3.58)\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}

(3.59)\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.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$.

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

(6.15)\begin{equation} \theta=(\alpha-\epsilon\varPhi)\varTheta, \end{equation}

and velocity expansions

(6.16a,b)\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

(6.17a,b)\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}

(6.18a,b)\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

(6.19a,b)\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

(6.20)\begin{align} &\frac{1}{r}L_r(r^2 u)+L_\varTheta v +v\cot(\alpha\varTheta-\epsilon\varPhi\varTheta)=0, \end{align}

(6.21)\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}

(6.22)\begin{gather} 2\tau_{r\theta}L_r u= \tau_{rr}(L_\varTheta u+L_r v - v), \end{gather}

(6.23)\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

(6.24a,b)\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

(6.25a,b)\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}

(6.26)\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

(6.27a,b)\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

(6.28)\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

(6.29a,b)\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

(6.30a–d)\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$:

(6.31)\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}

(6.32)\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}

(6.33)\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}

(6.34)\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

(6.35a,b)\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

(6.36)\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

(6.37)\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

(6.38)\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

(6.39a–c)\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

(6.40a,b)\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

(6.41)\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

(6.42)\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,

(6.43)\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

(6.44)\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

(6.45a–c)\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

(6.46)\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$.