2. General problem definition

We consider a particle in two or three dimensions, and assume a Bingham constitutive law for the surrounding fluid, so that the deviatoric stress tensor $\boldsymbol{\tau }$ satisfies the constitutive equation

(2.1) \begin{equation} \boldsymbol{\tau }=\left (\mu +\frac {\tau _Y}{\dot {\gamma }}\right )\dot {\boldsymbol{\gamma }} \text{ when } \tau \geqslant \tau _Y, \quad \text{and } \dot {\gamma }=0 \text{ otherwise}. \end{equation}

Here, $\tau _Y$ is the yield stress, $\mu$ is the plastic viscosity, $\dot {\boldsymbol{\gamma }}\equiv \boldsymbol{\nabla }\boldsymbol{u}+ (\boldsymbol{\nabla }\boldsymbol{u} )^{\rm T}$ is twice the symmetric strain rate tensor for the velocity $\boldsymbol{u}$ , and $\dot {\gamma }$ and $\tau$ are the second invariants of $\dot {\boldsymbol{\gamma }}$ and $\boldsymbol{\tau }$ , respectively. We assume no-slip on the stationary particle ( $\boldsymbol{u}=\boldsymbol{0}$ ), and impose a uniform velocity $\boldsymbol{U}_0$ of magnitude $U_0$ in the far field. In all cases, we non-dimensionalise lengths by a typical length scale of the particle (e.g. its radius) $a$ , velocities by $U_0$ , and stresses and pressures by the viscous stress scale $\mu U_0/a$ . Relabelling all quantities to their dimensionless counterparts, the constitutive equation becomes

(2.2) \begin{equation} \boldsymbol{\tau }=\left (1+\frac {\textit{Bi}}{\dot {\gamma }}\right )\dot {\boldsymbol{\gamma }} \text{ when } \tau \geqslant Bi, \quad \text{and } \dot {\gamma }=0 \text{ otherwise}. \end{equation}

The dimensionless parameter

(2.3) \begin{equation} Bi = \frac {a\,\tau _Y}{\mu U_0} \end{equation}

is the Bingham number, representing a dimensionless yield stress. In all cases, we investigate the small Bingham number regime, ${\textit{Bi}}\ll 1$ . We can also form a Reynolds number

(2.4) \begin{equation} Re = \frac {\rho U_0 a}{\mu }, \end{equation}

which we assume to be sufficiently small that inertia does not enter the governing equations at the orders that we consider. The neglect of fluid inertia certainly requires that $Re\ll 1$ . However, we additionally require that inertial effects are smaller than those due to the yield stress, and the precise requirement in terms of ${\textit{Bi}}$ and $Re$ will be shown to be geometry-dependent, and will be given below. The key quantity that we wish to determine is the dimensionless force on the particle (or the dimensionless force per unit width in two dimensions) $\boldsymbol{F}$ .

The inertialess governing equations express mass conservation and the balance of forces, and are given by

(2.5) \begin{equation} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}=0\quad \hbox{and}\quad \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\tau }-\boldsymbol{\nabla }p=0, \end{equation}

where $p$ is the pressure. The boundary conditions are no-slip on the particle surface ( $\boldsymbol{u}=\boldsymbol{0}$ ), and $\boldsymbol{u}\to \hat {\boldsymbol{U}}_0$ (a unit vector in the direction of the imposed far-field velocity) as $|\boldsymbol{x}|\to \infty$ . Where we report numerical solutions of this problem below, we will refer to it as the ‘full problem’, and indicate the particular geometry considered, in order to distinguish it from solutions of the leading-order problem in a given asymptotic regime.

While we have formulated this problem in terms of the drag force exerted on the particle, by a fluid with imposed far-field velocity, we could equally have analysed the particle velocity, in response to an external force applied to it, say $\tilde {\boldsymbol{F}}_g$ (dimensional). For example, one may be interested in the settling velocity due to the submerged weight of the particle. For this statement of the problem, it is natural to non-dimensionalise stresses by $\tilde {F}_g/a^2$ and velocities by $\tilde {F}_g/\mu a$ . The dimensionless yield stress is then the Oldroyd number ${\textit{Od}} = a^2\tau _Y/\tilde {{F}}_g$ , and the unknown is the dimensionless particle speed $\boldsymbol{U}_s$ , in response to a unit-magnitude dimensionless force $\hat {\boldsymbol{F}}_g$ . The two dimensionless problems are clearly related. Suppose from the former that we obtain a dimensionless force $\boldsymbol{F}$ acting on a particle due to a uniform stream of unit velocity. Then the dimensional force is $\mu a U_0 \boldsymbol{F}$ , and the far-field velocity is $\boldsymbol{U}_0$ . In the alternative formulation, this corresponds to a driving force of magnitude $\tilde {F}_g=\mu a U_0\,|\boldsymbol{F}|$ , and a dimensionless particle velocity of magnitude $U_s=\mu aU_0/\tilde {F}_g=1/|\boldsymbol{F}|$ . The direction of the external force and particle velocity are opposite to those of the force and far-field velocity in the former problem. The Oldroyd number can be related to the Bingham number as

(2.6) \begin{equation} {\textit{Od}} = \frac {a^2\tau _Y}{\tilde {{F}}_g}\equiv \frac {\textit{Bi}}{\left |\boldsymbol{F}\right |}. \end{equation}

Our focus is on the low yield stress regime, ${\textit{Od}}\ll 1$ . However, we note that there is a minimum force below which the fluid remains unyielded and the particle remains static, denoted by a non-vanishing critical Oldroyd number ${\textit{Od}}_c$ . While we do not concern ourselves with this limit in this study, we comment that in some cases, ${\textit{Od}}_c$ may itself be quite small (e.g. ${\textit{Od}}_c\sim O(10^{-2})$ ; see Beris et al. Reference Beris, Tsamopoulos, Armstrong and Brown1985; Hewitt & Balmforth Reference Hewitt and Balmforth2018). Thus the study of flows with low yield stresses corresponds to Oldroyd numbers significantly smaller than the critical value, ${\textit{Od}}\ll {\textit{Od}}_c$ .

3. Flow past a sphere

We analyse uniform flow past a sphere, and employ axisymmetric spherical polar coordinates $(r,\theta )$ , with the far-field flow aligned with $\theta =0$ . We use the radius of the sphere as the length scale with which the problem is rendered dimensionless, and introduce the Stokes streamfunction $\psi (r,\theta )$ such that the velocities in the radial and polar direction, $(u,v)$ respectively, are given by

(3.1) \begin{gather} u = \frac {1}{r^2\sin \theta }\frac {\partial \psi }{\partial \theta }, \quad v=-\frac {1}{r\sin \theta }\frac {\partial \psi }{\partial r}. \end{gather}

No-slip boundary conditions on the surface of the sphere are then given by

(3.2) \begin{equation} \psi = \frac {\partial \psi }{\partial r}=0 \quad \text{on }\ r=1, \end{equation}

and a unit uniform flow in the far field corresponds to

(3.3) \begin{equation} \psi \sim -\frac {r^2}{2}\sin ^2 \theta \quad \text{as }\ r\to \infty . \end{equation}

In the regime ${\textit{Bi}}\ll 1$ , the problem close to the sphere reduces to the flow of a Newtonian fluid, since the ‘plastic’ contribution to the deviatoric shear stress is negligible (2.2). Then up to $O({\textit{Bi}})$ , the solution satisfying (3.2) is given by

(3.4) \begin{equation} \psi = -\frac {A}{4}\left (2r^2-3r+\frac {1}{r}\right )\sin ^2\theta +O({\textit{Bi}}). \end{equation}

In this expression, the coefficient $A$ depends on the Bingham number, and cannot be determined straight away since, as shown by Beris et al. (Reference Beris, Tsamopoulos, Armstrong and Brown1985) and below, the expression becomes non-asymptotic at large distances from the sphere. This means that we cannot impose the far-field boundary condition (3.3). Instead, the determination of $A$ comes through matching asymptotic expansions in the near and far fields. One approach here is to expand $A$ in an asymptotic series in ${\textit{Bi}}$ , and solve order by order. For the problem in this section with spherical geometry, this approach gives two terms below $O({\textit{Bi}})$ . In the later problems, however, we will show that that there is an infinite asymptotic series with terms larger than $O({\textit{Bi}})$ . Moreover, since every term has the same functional form, and the leading-order problem is linear, we can instead treat $A({\textit{Bi}})$ as a single ${\textit{Bi}}$ -dependent function, and perform the matching procedure with the far field to determine an algebraic relation for $A({\textit{Bi}})$ . This algebraic expression can subsequently be expanded for ${\textit{Bi}}\ll 1$ if desired.

The force on the sphere due to the solution (3.4) is $6\pi A + O({\textit{Bi}})$ in the flow direction, and the strain rate scales like $A/r^2$ when $r\gg 1$ . Thus when $r=O(\sqrt {A/Bi})$ , the strain rate is $O({\textit{Bi}})$ , and the neglected, plastic, term in (2.2) becomes of the same order as the viscous term. This is the viscoplastic equivalent to the Whitehead paradox for flows with small but non-vanishing Reynolds numbers. We make the rescaling to an outer region in which viscous and plastic stresses are comparable via

(3.5) \begin{equation} r = \sqrt {\frac {6\pi A}{\textit{Bi}}}\,R, \quad \psi = \sqrt {\frac {(6\pi A)^3}{\textit{Bi}}}\,\varPsi . \end{equation}

Note the use of capital letters to represent quantities in the outer region. We further subtract the free stream velocity (i.e. change to the frame of reference of the unyielded fluid at infinity), writing

(3.6) \begin{equation} \varPsi = -\frac {R^2}{2\sqrt {6\pi A\, {\textit{Bi}}}}\sin ^2\theta + \varPsi _0. \end{equation}

We rescale the strain rate tensor, the deviatoric stress tensor and the pressure as

(3.7) \begin{equation} \dot {\boldsymbol{\gamma }}=Bi\,\dot {\boldsymbol{\varGamma }}, \quad \boldsymbol{\tau }=Bi\, \unicode{x1D64F}, \quad p = Bi\, P. \end{equation}

The scaled velocity fields in the radial and polar directions associated with $\varPsi _0$ are given by

(3.8) \begin{equation} U=\frac {1}{R^2\sin \theta }\frac {\partial \varPsi _0}{\partial \theta }, \quad V=-\frac {1}{R\sin \theta }\frac {\partial \varPsi _0}{\partial R}, \end{equation}

and the far-field boundary condition requires that $U,\,V\to 0$ at $R\to \infty$ . The near-field boundary condition matches the streamfunction to the leading contribution from the inner field (3.4), thus

(3.9) \begin{equation} \varPsi _0\to \frac {R}{8\pi }\sin ^2\theta \quad \hbox{as}\ \ R\to 0. \end{equation}

This condition may be compactly included as a point forcing in the balance of momentum in the outer region. In this way we find that

(3.10) \begin{gather} -\boldsymbol{\nabla }\boldsymbol{\cdot }\unicode{x1D64F}+\boldsymbol{\nabla }P = \delta ^3({\boldsymbol{R}})\,\boldsymbol{e}_z \quad \text{and} \quad \unicode{x1D64F} = \left (1+\frac {1}{\dot {\varGamma }}\right )\dot {\boldsymbol{\varGamma }}, \end{gather}

where $\delta ^3$ is the Dirac delta function in three dimensions, $\boldsymbol{R}$ is the rescaled position vector, and $\boldsymbol{e}_z$ is the unit vector in the $\theta =0$ direction.

Thus, the outer problem (3.10) corresponds to a parameter-free, fully nonlinear viscoplastic problem, with vanishing velocity in the far field, driven by a point force of unit magnitude at the origin, which we term the ‘viscoplastic Stokeslet in three dimensions’. This problem requires numerical computation, but its solution leads to an asymptotic prediction of the flow fields and drag forces applicable for ${\textit{Bi}}\ll 1$ , thus providing considerable advantages over the requirement for individual simulations at each value of ${\textit{Bi}}$ . We note that rescaling to a parameter-free problem is a strategy employed in other studies of viscoplastic flows in the regime of small Bingham number. Examples include the flow far from a swimming sheet (Hewitt & Balmforth Reference Hewitt and Balmforth2017) and the flow in the neighbourhood of a stagnation point on a no-slip boundary (Taylor-West & Hogg Reference Taylor-West and Hogg2025); both of these examples correspond to two-dimensional flows (as in § 5 of the current paper), but are driven by boundary conditions different from those considered here.

In order to integrate (3.10) numerically, we must regularise the point force, for which we write

(3.11) \begin{equation} \delta ^3(\boldsymbol{R}) \approx \frac {1}{\varepsilon ^3}\exp \left (-\frac {\pi R^2}{\varepsilon ^2}\right )\!, \end{equation}

where $\varepsilon \ll 1$ is a regularisation parameter. This expression means that we have introduced a small length scale $\varepsilon$ at which the regularised approximation of the point force is no longer accurate. As will be explained below, asymptotic matching between the inner and outer regions to determine the effects of viscoplasticity on the purely viscous solution requires the evaluation of the numerical solution for $\varPsi _0$ in the regime $R\ll 1$ . Evidently some care is required in selecting an appropriate value of $\varepsilon$ that allows the near-field behaviour to be consistently evaluated, and this is reported below. Our numerical strategy to compute the solution of (3.10) does not regularise the rheological model. Instead, we adopt the FISTA algorithm proposed by Treskatis, Moyers-González & Price (Reference Treskatis, Moyers-González and Price2016), implemented in FEniCS (Alnæs et al. Reference Alnæs2015). Further details of the numerical method are provided in Appendix A.

We seek the limiting behaviour of this outer problem (3.10), when $R\ll 1$ , for which we consider the series in $R$

(3.12) \begin{equation} \varPsi _0 =\frac {1}{8\pi }\left (R\,f_0(\theta ) + R^2\,f_1(\theta )+\cdots \right )\!, \end{equation}

where the scale of the first term is fixed by the point force at the origin. Up to this order in $R$ , the plastic term does not enter the constitutive equation, thus both terms satisfy the Newtonian Stokes flow problem. This produces four possible independent solutions for $f_0$ and $f_1$ , but only $f_i=c_i\sin ^2\theta$ satisfies the regularity and symmetry conditions at $\theta =0$ . We have $c_0=1$ as a consequence of force balance, while $c_1$ requires numerical determination (see figure 1).

Figure 1.

Numerical evidence for the limiting behaviour (3.13) of the viscoplastic Stokeslet in three dimensions. In each panel, the dashed line is the prediction of (3.13), and the colours are from numerical simulations of the viscoplastic Stokeslet problem. (a) Scaled velocities $R U$ and $R V$ as functions of $\theta$ at fixed $R=0.001$ . The numerical solutions are for regularisation parameter $\varepsilon =10^{-5}$ . (b) Difference between $R U$ and the leading-order Stokeslet contribution to (3.13) as a function of $R$ along the axis of symmetry $\theta =0$ . The values of the regularisation parameter are shown in the legend. Note that the numerical results deviate from the asymptotic prediction at small $R$ due to the regularisation of the point force, but the agreement persists to lower $R$ as the regularisation parameter is reduced.

After determining $c_1=-7.6$ (to two significant figures) from the numerical solutions, we thus deduce that the limiting behaviour of this outer problem is

(3.13) \begin{equation} \varPsi _0\sim \frac {1}{8\pi }\sin ^2\theta (R - 7.6R^2)+\cdots \quad \text{as }\ R\to 0. \end{equation}

We reiterate that the first term in (3.13) corresponds to the Newtonian Stokeslet, and the second term is a uniform velocity of $7.6/4\pi$ in the far-field flow direction, arising due to the viscoplasticity. Figure 1(a) plots the variation of $RU$ and $RV$ at $R=10^{-3}$ , evidencing that (3.13) provides an accurate representation of the flow field at fixed $R\ll 1$ . In figure 1(b), we show the effects of the regularisation of the point force given by (3.11) on the computed velocity fields. In particular, in this panel we plot the difference between $RU(R,\theta =0)$ and $1/(4\pi )$ , the latter arising as the first term of (3.13). The power series expansion when $R\ll 1$ anticipates that this difference must vary linearly with $R$ , and this is what we find for $\varepsilon \ll R\ll 1$ . When $R=O(\epsilon )$ , we observe that the effects of the regularisation impact the numerically measured $1/(4\pi )-RU(R,\theta =0)$ . We therefore use the result with $\varepsilon =10^{-5}$ , and apply linear regression over the range $10^{-3}\leqslant R\leqslant 10^{-2}$ to determine $c_1=-7.6$ .

We can now carry out a matching between the outer and inner regions by examining the limits $R\to 0$ and $r\to \infty$ , respectively. In the regime $R\ll 1$ and using (3.5), we have found that

(3.14) \begin{align} \psi &= \sqrt {\frac {(6\pi A)^3}{\textit{Bi}}}\left ( - \frac {\textit{Bi}}{6\pi A}\frac {r^2}{2\sqrt {6\pi A\, {\textit{Bi}}}} + \frac {1}{8\pi }\sqrt {\frac {\textit{Bi}}{6\pi A}}\,r-\frac {7.6}{8\pi }\frac {\textit{Bi}}{6 \pi A}r^2 \right )\sin ^2\theta +\cdots \nonumber\\ &=-\frac {1}{4}\left (\left (2 +3.8\sqrt {\frac {6 A\, {\textit{Bi}}}{\pi }}\right )r^2- 3Ar\right )\sin ^2\theta +\cdots \end{align}

from the outer expansion, and

(3.15) \begin{equation} \psi = -\frac {1}{4}\left (2Ar^2 - 3Ar\right )\sin ^2\theta + \cdots \end{equation}

from the inner expansion. Thus, matching the two, we obtain the relation

(3.16) \begin{equation} 1 = A - 1.9\sqrt {\frac {6A {\textit{Bi}}}{\pi }}. \end{equation}

Now rearranging and expanding up to $O({\textit{Bi}})$ , we find

(3.17) \begin{equation} A = 1 + 1.9\sqrt {\frac {6 {\textit{Bi}}}{\pi }} + O({\textit{Bi}}), \end{equation}

and the force on the sphere is given by

(3.18) \begin{equation} F=6\pi A = 6\pi + 11.4\sqrt {6 \pi {\textit{Bi}}} + O({\textit{Bi}}). \end{equation}

We note that since this arises from a scalar multiple of the Newtonian solution, the partition of the leading-order correction to the force, $11.4\sqrt {6\pi {\textit{Bi}}}$ , into contributions from normal and shear stresses, is precisely the same as in the Newtonian case (i.e. the shear and normal stress contributions are in the ratio $2:1$ ).

As discussed in § 2, we may alternatively reinterpret this result as the velocity of the sphere in response to an external force. In this way, we write $U_s=1/F$ and ${\textit{Od}}=Bi/F$ to obtain

(3.19) \begin{equation} U_s = \frac {1}{6\pi }\left (1 - 11.4\sqrt {\textit{Od}}\right )+O({\textit{Od}}). \end{equation}

This shows how weak yield stress effects retard the particle, compared to the same force acting on a sphere in Newtonian fluid.

We verify our asymptotic result for the drag force (3.18) through further numerical computations. To this end, we numerically compute the solution to the full problem of viscoplastic flow around the sphere, driven by a unit far-field velocity at various values of dimensionless yield stress ( $10^{-5}\lt Bi\lt 10^{2}$ ). We again utilise the FISTA algorithm (see Appendix A). These computations supplement and extend the results due to Beris et al. (Reference Beris, Tsamopoulos, Armstrong and Brown1985), who had employed a regularised rheology. In figure 2, we plot the excess dimensionless force $F-6\pi$ as a function of ${\textit{Bi}}$ . We note that there is excellent agreement between the numerical results and the asymptotic solution (3.18). We have overlain the numerical results of Beris et al. (Reference Beris, Tsamopoulos, Armstrong and Brown1985), which also agree well but do not extend to as small values of ${\textit{Bi}}$ (figure 2).

Figure 2.

Viscoplastic correction to the drag force, $F-6\pi$ , on a sphere as a function of the Bingham number ${\textit{Bi}}$ . Included are the asymptotic prediction $11.4\sqrt {6\pi {\textit{Bi}}}$ (dashed line), results from numerical solutions of the full problem (red circles), corresponding numerical results of Beris et al. (Reference Beris, Tsamopoulos, Armstrong and Brown1985) (solid line), and experimental data from Ansley & Smith (Reference Ansley and Smith1967) (blue circles) and Arabi & Sanders (Reference Arabi and Sanders2016) (green circles). The latter have been grouped into those with Bingham numbers smaller and larger than the square of the Reynolds number, ${\textit{Re}}^2$ .

Also shown in figure 2 are experimental results from Ansley & Smith (Reference Ansley and Smith1967) and Arabi & Sanders (Reference Arabi and Sanders2016). These experimental results provide reasonable agreement with the numerical solutions at larger Bingham numbers, but they diverge from the theory somewhat when ${\textit{Bi}}\leqslant 1$ . We note that the experimental measurements may have been influenced by weak inertial effects. Proudman & Pearson (Reference Proudman and Pearson1957) showed that when $Re\ll 1$ , the drag exerted by a uniform flow aorund a stationary sphere is modified by $O(Re)$ , and that inertial effects on the velocity field become comparable with viscous processes when $r\sim 1/Re$ . In this subsection, however, we have shown that the viscoplastic effects become non-negligible when $r=O(Bi^{-1/2})$ , and the correction to the drag force is $O(Bi^{1/2})$ . Thus the neglect of inertia requires that

(3.20) \begin{equation} {\textit{Re}}^2\ll Bi\ll 1. \end{equation}

The experimental data are not all collected in this regime, so are not expected to be accurately predicted by the theory developed here.

Figure 3(a) shows an example of a numerical solution for the full problem of flow around a sphere at ${\textit{Bi}}=10^{-2}$ . The colour plot of strain rate shows that the fluid is unyielded (coloured black) in the far field, and yielded in an envelope around the sphere. Streamsurfaces are also shown in the frame of reference of the unyielded fluid. Figure 3(b) shows the boundary of the yielded envelope around the sphere from three such numerical simulations of the full problem, at different Bingham numbers, after rescaling to the outer coordinates via (3.5). Also shown is the yielded envelope from the viscoplastic Stokeslet simulation in three dimensions, with $\varepsilon =10^{-5}$ , evidencing that the outer problem effectively captures the behaviour far from the sphere. This is further supported by the streamsurfaces for the Stokeslet problem shown in figure 3(b), which closely mimic those in figure 3(a) for the flow around the sphere. In the rescaled coordinates (3.5), we find that the yielded envelope extends a distance approximately 0.31 upstream (and downstream) from the sphere, and a distance approximately 0.44 to the sides of the sphere. In unscaled dimensional units and correct to $O(Bi^{1/2})$ , after substituting (3.17) for $A$ , this corresponds to distances

(3.21) \begin{equation} 0.31\left (\sqrt {\frac {6\pi }{\textit{Bi}}}+5.7\right )a \quad \text{and} \quad 0.44\left (\sqrt {\frac {6\pi }{\textit{Bi}}}+5.7\right )a, \end{equation}

where $a$ is the radius of the sphere.

Figure 3.

(a) Example numerical simulation for the full problem of flow around a sphere at ${\textit{Bi}}=1/100$ . The colour plot shows the strain rate on a logarithmic scale (with black representing unyielded regions), and the black dashed lines indicate a selection of streamsurfaces in the frame of reference of the unyielded, far-field fluid. A quarter of the unit-radius sphere is visible in white at the bottom left corner of the plot.(b) Outer yield surface (plotted at $\tau =1.001\,Bi$ ) from numerical simulations for the three-dimensional Stokeslet problem (solid black) and the full problem of viscoplastic flow around a sphere at Bingham numbers shown in the legend. The latter have been rescaled to the outer coordinates, via (3.5). Also shown are the streamsurfaces for the Stokeslet solution (black dashed lines). In both panels, the flow is in the anticlockwise direction.

4. Arbitrarily shaped particles in three dimensions

We now outline how the matched asymptotic expansion for the drag force exerted by a viscoplastic flow around a sphere may be extended to an arbitrary particle in three dimensions. In this case, the length scale $a$ featured in the Bingham number (2.3) is set by an appropriate choice for the particle, e.g. the radius of the smallest sphere containing the particle. The argument follows as above. Below $O({\textit{Bi}})$ , the governing equation is the linear Stokes equation, thus the velocity field takes the form

(4.1) \begin{equation} \boldsymbol{u}_i = A_1({\textit{Bi}})\,\boldsymbol{u}_1+A_2({\textit{Bi}})\,\boldsymbol{u}_2+A_3({\textit{Bi}})\,\boldsymbol{u}_3, \end{equation}

where the $\boldsymbol{u}_i$ are three linearly independent solutions to Stokes flow around the particle. We can choose these to be solutions with unit far-field velocity in each of three orthogonal coordinate directions, $\{\boldsymbol{e}_1,\boldsymbol{e}_2,\boldsymbol{e}_3\}$ . The force on the particle, $\boldsymbol{F}$ , due to this flow is linear in the far-field velocity, related by a resistance matrix (or inverse mobility matrix), $ \unicode{x1D648}^{{\kern1pt}-1}$ :

(4.2) \begin{equation} F_i = M^{-1}_{\textit{ij}}A_{\!j}. \end{equation}

We write the drag force as $\boldsymbol{F}=F\hat {\boldsymbol{e}}_{\!f}$ , where $|\hat {\boldsymbol{e}}_{\!f}|=1$ , and substitute the magnitude $F$ into the rescaling of the outer region:

(4.3) \begin{equation} r=\sqrt {\frac {F}{\textit{Bi}}}\,R \quad \hbox{and}\quad \boldsymbol{u}=\sqrt {F\,{\textit{Bi}}}\,\boldsymbol{U}. \end{equation}

Crucially, this outer problem is the same as that for the sphere, namely the viscoplastic Stokeslet in three dimensions. Thus we can directly take the near-field limit of the outer solution, from (3.13) (after reorienting along the direction of the force, $\hat {\boldsymbol{e}}_{\!f}$ ), as

(4.4) \begin{align} \mathbf{u}&=\hat {\boldsymbol{U}}_0-\sqrt {F\,{\textit{Bi}}}\left (\frac {1}{8\pi }\left (\frac {\hat {\boldsymbol{e}}_{\!f}}{\sqrt {Bi/F}r}+\frac {(\hat {\boldsymbol{e}}_{\!f}\boldsymbol{\cdot }\boldsymbol{r)}\boldsymbol{r}}{\sqrt {Bi/F}r^3}\right )-\frac {7.6}{4\pi }\hat {\boldsymbol{e}}_{\!f}\right )+\cdots \nonumber \\ &=\hat {\boldsymbol{U}}_0-\frac {1}{8\pi }\left (\frac {\boldsymbol{F}}{r}+\frac {(\boldsymbol{F}\boldsymbol{\cdot }\boldsymbol{r})\boldsymbol{r}}{r^3}\right )+\frac {7.6\sqrt {F\,{\textit{Bi}}}}{4\pi }\hat {\boldsymbol{e}}_{\!f}+\cdots , \end{align}

where $\hat {\boldsymbol{U}}_0$ is the direction of the imposed far-field velocity. The far-field limit of the inner solution is the uniform stream $\boldsymbol{u}=\boldsymbol{A}({\textit{Bi}})\equiv (A_1({\textit{Bi}}),A_2({\textit{Bi}}),A_3({\textit{Bi}}))$ , plus the Stokeslet of strength $-\boldsymbol{F}$ :

(4.5) \begin{equation} \boldsymbol{u}=\boldsymbol{A}({\textit{Bi}})-\frac {1}{8\pi }\left (\frac {\boldsymbol{F}}{r}+\frac {(\boldsymbol{F}\boldsymbol{\cdot }\boldsymbol{r})\boldsymbol{r}}{r^3}\right )+\cdots . \end{equation}

Thus the matching gives

(4.6) \begin{equation} \hat {\boldsymbol{U}}_0 = \boldsymbol{A}-\frac {7.6}{4\pi }\sqrt {\textit{Bi}}\, \frac { \unicode{x1D648}^{{\kern1pt}-1}\boldsymbol{A}}{\sqrt {| \unicode{x1D648}^{{\kern1pt}-1}\boldsymbol{A}|}}+O({\textit{Bi}}). \end{equation}

This implies

(4.7) \begin{equation} \boldsymbol{A} = \hat {\boldsymbol{U}}_0 + \frac {7.6}{4\pi }\sqrt {\textit{Bi}}\,\frac { \unicode{x1D648}^{{\kern1pt}-1}\hat {\boldsymbol{U}}_0}{\sqrt {| \unicode{x1D648}^{{\kern1pt}-1}\hat {\boldsymbol{U}}_0|}} + O({\textit{Bi}}), \end{equation}

and the force on the particle is

(4.8) \begin{align} \boldsymbol{F}= \unicode{x1D648}^{{\kern1pt}-1}\boldsymbol{A}+O({\textit{Bi}}) = \unicode{x1D648}^{{\kern1pt}-1}\hat {\boldsymbol{U}}_0 + \frac {7.6}{4\pi }\sqrt {\textit{Bi}}\,\frac { \unicode{x1D648}^{{\kern1pt}-2}\hat {\boldsymbol{U}}_0}{\sqrt {| \unicode{x1D648}^{{\kern1pt}-1}\hat {\boldsymbol{U}}_0|}} + O({\textit{Bi}}). \end{align}

Note that this reduces to the result for the spherical particle, for which the resistance matrix is $ \unicode{x1D648}^{{\kern1pt}-1}=6\pi \unicode{x1D644}$ , so for a unit magnitude far-field velocity, the force is parallel to the imposed velocity and of magnitude given by (3.18). However, for a more general particle, (4.8) establishes that $\boldsymbol{F}$ is not necessarily aligned with $\hat {\boldsymbol{U}}_0$ .

While (4.8) applies for any particle in the ${\textit{Bi}}\ll 1$ regime, we can make no conclusions about the rate at which the asymptotic solution is approached. Particles with complex shapes may have Stokes solutions exhibiting regions of low strain rate, which could result in significant plugs at small but finite Bingham numbers. In general, however, the strain rate can only completely vanish at isolated points in the flow, as a consequence of the regularity of solutions to Stokes flow (e.g. Ladyzhenskaya Reference Ladyzhenskaya1969, p. 42). Hence any plugs near the particle will become infinitesimally small as ${\textit{Bi}}\to 0$ , and have a negligible effect on the leading-order solution. Nonetheless, for some particles (e.g. those with significantly non-convex shapes), these local plugged regions may only become negligible for very small Bingham numbers, resulting in a slow approach to the asymptotic result (4.8).

Now we again switch perspective to considering the velocity of the particle $\boldsymbol{U}_s=-\hat {\boldsymbol{U}}_0/|\boldsymbol{F}|$ in response to a unit-magnitude imposed force $\hat {\boldsymbol{F}}_g=-\boldsymbol{F}/|\boldsymbol{F}|$ . This gives

(4.9) \begin{equation} \boldsymbol{U}_s = \left ( \unicode{x1D648}-\frac {7.6}{4\pi }\sqrt {\textit{Od}}\, \unicode{x1D644}\right )\hat {\boldsymbol{F}}_g+\cdots . \end{equation}

Thus the velocity is reduced slightly by the weak yield stress, in an isotropic manner. A result of this is that any preferred direction is enhanced. For example, for a slender prolate ellipsoid, of unit length and aspect ratio $\epsilon \ll 1$ , the mobility matrix $ \unicode{x1D648}$ is approximately (Oberbeck Reference Oberbeck1876; Happel & Brenner 1983)

(4.10) \begin{align} \unicode{x1D648}=\frac {1}{8\pi }\begin{pmatrix} 2\left (\ln (2/\epsilon )-\dfrac {1}{2}\right ) & 0 & 0\\ 0 & \ln (2/\epsilon )+\dfrac {1}{2} & 0 \\ 0 & 0 & \ln (2/\epsilon )+\dfrac {1}{2} \end{pmatrix}, \end{align}

in a coordinate system aligned with the major axis and two perpendicular directions. Thus when $\epsilon \ll 1$ , there is approximately a factor of two between the settling velocities of particles with their major axes aligned with or perpendicular to the direction of gravity. In the presence of a weak yield stress ( ${\textit{Od}}\ll 1$ ), the matrix relating the settling velocity to the imposed force becomes

(4.11) \begin{align} \frac {1}{8\pi }\!\begin{pmatrix} 2\left (\!\ln (2/\epsilon )-\dfrac {1}{2}\!\right )\!-\!15.2\sqrt {\textit{Od}}& 0 & 0\\ 0 & \ln ({2}/\epsilon )+\dfrac {1}{2}\!-\!15.2\sqrt {\textit{Od}} & 0 \\ 0 & 0 & \ln (2/\epsilon )+\dfrac {1}{2}\!-\!15.2\sqrt {\textit{Od}} \end{pmatrix}\!, \end{align}

so the ratio of settling velocities becomes

(4.12) \begin{equation} 2\frac {\ln (2/\epsilon )-1/2}{\ln (2/\epsilon )+1/2}+\frac {15.2\left (\ln (2/\epsilon )-3/2\right )}{(\ln (2/\epsilon )+1/2)^2}\sqrt {\textit{Od}}+O({\textit{Od}}). \end{equation}

This represents a small increase in the tendency of the particle to settle in a direction parallel to its major axis. This is consistent with the results of Hewitt & Balmforth (Reference Hewitt and Balmforth2018), who demonstrated that in the regime of large Oldroyd number, long cylindrical bodies have a strong tendency to slide through the fluid along their axis, due to significant drag anisotropy. The result here demonstrates how this effect emerges at small Oldroyd numbers, albeit for slender ellipsoids, rather than cylindrical objects.

5. Flow past a circular cylinder

We now analyse two-dimensional uniform flow around a circular cylinder in the regime of low Bingham number (where the Bingham number is still given by (2.3), with $a$ now representing the dimensional radius of the cylinder). In this section, we employ plane polar coordinates $(r,\theta )$ , with $\theta =0$ along the direction of the flow, and the origin at the centre of the cylinder, noting that these independent variables now differ from the definitions in § 3. We also introduce the streamfunction for the flow, $\psi (r,\theta )$ , and in what follows will match between an inner, viscous region close to the cylinder and a distant outer region in which viscoplasticity affects the dynamics. In this way, the methodology is similar to that in § 3, but the changes of geometry and coordinate systems imply that the variables represent slightly different physical quantities – however, the commonality of the approach means that there are advantages in adopting this overlapping notation.

One significant difference in the analysis for circular cylinders is that no solution for a purely viscous problem can be constructed with no-slip boundary conditions on the surface of the cylinder and a uniform flow in the far field. Instead, the Stokes solution features logarithmic terms, leading to the Stokes paradox that is classically resolved by the introduction of weak inertial effects (Proudman & Pearson Reference Proudman and Pearson1957; Hinch Reference Hinch1991). Our analysis shows how this difficulty is resolved by weak yield stress effects for viscoplastic fluids.

We introduce the streamfunction $\psi$ in two dimensions such that the velocities in the radial and polar directions, $u$ and $v$ , respectively, are given by

(5.1) \begin{equation} u = \frac {1}{r}\frac {\partial \psi }{\partial \theta }, \quad v = -\frac {\partial \psi }{\partial r}. \end{equation}

The uniform far-field velocity is given by

(5.2) \begin{equation} \psi \sim -r\sin \theta \quad \text{as }\ r \to \infty , \end{equation}

while no-slip on the surface of the cylinder requires $\psi =\partial \psi /\partial r=0$ at $r=1$ . Up to $O({\textit{Bi}})$ , the problem is the Newtonian Stokes flow problem, which in two dimensions is given by ${\nabla} ^4\psi =0$ . The solution satisfying a no-slip boundary condition with the appropriate symmetry about $\theta =0$ is given by

(5.3) \begin{equation} \psi = A({\textit{Bi}})\left (r\ln \left (\frac {1}{r}\right ) +\dfrac {1}{2}r - \frac {1}{2r}\right )\sin \theta . \end{equation}

As in § 3, we take $A$ to depend on the Bingham number. The solution (5.3) has a strain rate that decays like $A/r$ , thus the viscous and plastic stresses become of the same order when $r=O(A/{\textit{Bi}})$ . We will show that the leading order of $A$ is not independent of ${\textit{Bi}}$ when ${\textit{Bi}}\ll 1$ , so we note that this scaling differs from the distance at which inertia becomes significant in the classical resolution of the Stokes paradox, and therefore also differs slightly from the suggestion of Hewitt & Balmforth (Reference Hewitt and Balmforth2018) (namely that the outer region is for $r=O(1/{\textit{Bi}})$ ). The force on the cylinder in this case is $\boldsymbol{F}=-4\pi A\, \boldsymbol{e}_x$ , where $\boldsymbol{e}_x$ is the unit vector along the $\theta =0$ axis. Similarly to the spherical case (§ 3), we now make a rescaling to an outer region,

(5.4) \begin{equation} R=\frac {\textit{Bi}}{4\pi A} r, \quad \varPsi _0=\frac {\textit{Bi}}{(4\pi A)^2}\left (\psi + \frac {4\pi A}{\textit{Bi}}R\sin \theta \right )\!, \end{equation}

in which $\varPsi _0(R,\theta )$ is the streamfunction for the perturbation to the uniform flow in the far field. In this way, our rescaling of the outer region diverges from the suggestion of Hewitt & Balmforth (Reference Hewitt and Balmforth2018), and also diverges from the analogous resolution of the Stokes paradox by the introduction of weak inertial effects (Hinch Reference Hinch1991). The problem for $\varPsi (R,\theta )$ reduces to that of an otherwise static Bingham fluid, driven by a point force per unit length of unit strength at the origin. We refer to this problem as the ‘viscoplastic Stokeslet in two dimensions’. In this case, neglecting inertial terms in the outer problem requires that $Re\ll Bi/A$ , which we will show corresponds to $Re\ll Bi\ln (1/{\textit{Bi}})$ .

The outer flow field must match the inner flow field when $r\gg 1$ , thus $\varPsi _0=(R\ln (1/R))/(4\pi )$ when $R\ll 1$ . This matching may be built in by including a two-dimensional point forcing in the balance of momentum. Thus the outer problem is given by

(5.5) \begin{gather} -\boldsymbol{\nabla }\boldsymbol{\cdot }\unicode{x1D64F}+\boldsymbol{\nabla }P = \delta ^2({\boldsymbol{R}})\,\boldsymbol{e}_x, \end{gather}

where, as in § 3, $ \unicode{x1D64F}$ , $P$ and $\boldsymbol{R}$ are respectively the stress tensor, pressure and position vector in the outer variables, and $\delta ^2$ is the Dirac delta function in two dimensions. To facilitate numerical integration, we must regularise the point force, which we do via

(5.6) \begin{equation} \delta ^2(\boldsymbol{R})\approx \frac {1}{\varepsilon ^2}\exp \left (-\frac {\pi R^2}{\varepsilon ^2}\right )\!. \end{equation}

Also as in § 3, we have introduced a small length scale $\varepsilon \ll 1$ over which this approximation is a poor representation of the point force. We integrate (5.5) numerically using the FISTA algorithm implemented in FEniCS (Alnæs et al. Reference Alnæs2015) as described in Appendix A. In this case, the near-field behaviour as $R\ll 1$ is given by the series solution

(5.7) \begin{equation} \varPsi _0 = R\ln (1/R)\,f_0(\theta )+R\,f_1(\theta )+\cdots . \end{equation}

The leading-order term is required to balance the point force at the origin, as described above. Plasticity does not enter until higher order in $R$ , thus we can find $f_0$ and $f_1$ by solving the Newtonian Stokes equations, and find that only $f_i=c_i\sin \theta$ satisfies the regularity and symmetry conditions at $\theta =0$ . Again, $c_0=1/4\pi$ is directly set by force balance, while $c_1$ requires numerical determination. Ultimately, we find that the leading-order behaviour is given by a Newtonian Stokeslet term plus an adjustment to the uniform stream:

(5.8) \begin{equation} \varPsi _0 = \frac {1}{4\pi }R\left (\ln \left (\frac {1}{R}\right )-3.4\right )\sin \theta +\cdots \quad \text{as } R\to 0. \end{equation}

The numerical evidence for this is given in figure 4. The coefficient 3.4 (to two significant figures) was determined from the intercept of a linear fit to the numerical results (with $\varepsilon =10^{-5}$ ) for the radial velocity along $\theta =0$ , over the range $10^{-4}\leqslant R\leqslant 10^{-3}$ , after subtracting the leading-order term proportional to $\ln (1/ R)$ (figure 4 b). We plot the outer velocity fields $U$ and $V$ at $R=10^{-4}$ as functions of $\theta$ , and show very close comparison with the fitted form (figure 4 a).

Figure 4.

Numerical evidence for the limiting behaviour, (5.8), of the viscoplastic Stokeslet in two dimensions. In each panel, the dashed line is the prediction of (5.8), and the colours are from numerical simulations of the viscoplastic Stokeslet problem. (a) Velocities $U$ and $V$ as functions of $\theta$ at fixed $R=0.0001$ . The numerical solutions are for regularisation parameter $\varepsilon =10^{-5}$ . (b) The radial velocity $U$ as a function of $R$ along the axis of symmetry $\theta =0$ . The values of the regularisation parameter are shown in the legend. Note that the numerical results deviate from the asymptotic prediction at small $R$ due to the regularisation of the point force, but that the agreement persists to lower $R$ as the regularisation parameter is reduced.

Now we use (5.8) to match between the inner and outer expansions, determining $A({\textit{Bi}})$ . From the outer, we have

(5.9) \begin{align} \psi = -r\sin \theta + \frac {(4\pi A)^2}{\textit{Bi}}\varPsi _0 +\cdots = \left (-1 + A \ln \left (\frac {4\pi A}{Bi\, r}\right )-3.4A\right )r\sin \theta +\cdots , \end{align}

which we match to the leading order from the inner (5.3):

(5.10) \begin{equation} \psi = \left (A\,\ln \left (\frac {1}{r}\right )+ \frac {A}{2}\right )r\sin \theta +\cdots . \end{equation}

Thus $A$ satisfies the algebraic relation

(5.11) \begin{equation} \ln \left (\frac {1}{\textit{Bi}}\right )+\ln (4\pi )-3.4-\dfrac {1}{2} =\frac {1}{A} +\ln \left (\frac {1}{A}\right ) \implies \frac {1}{A}=\mathrm{W}\left (\frac {4\pi \exp (-3.9)}{\textit{Bi}}\right )\!, \end{equation}

where $\mathrm{W}$ is the Lambert $\mathrm{W}$ function. This gives the expression for the dimensionless force on the cylinder:

(5.12) \begin{equation} F=4\pi A + O({\textit{Bi}}) = \frac {4\pi }{\mathrm{W}\!\left (4\pi \exp (-3.9)/Bi\right )} + O({\textit{Bi}}). \end{equation}

Alternatively, expanding the Lambert $\mathrm{W}$ function for large argument (or directly analysing the first expression of (5.11)), one obtains

(5.13) \begin{equation} A^{-1}=\ln (1/{\textit{Bi}}) -\ln (\ln (1/{\textit{Bi}})) + \ln (4\pi ) - 3.9 +\cdots , \end{equation}

thus

(5.14) \begin{equation} F = \frac {4\pi }{\ln (1/{\textit{Bi}})}+\frac {4\pi \ln (\ln (1/{\textit{Bi}}))}{[\ln (1/{\textit{Bi}})]^2}-\frac {4\pi (\ln (4\pi )-3.9)}{[\ln (1/{\textit{Bi}})]^2}+\cdots . \end{equation}

The leading term of (5.14) is as given by Hewitt & Balmforth (Reference Hewitt and Balmforth2018). However, the convergence of (5.14) is very slow for ${\textit{Bi}}\ll 1$ , and as shown in figure 5, the first term significantly underpredicts the force. Moreover, there are an infinite number of terms in the expansion before reaching the $O({\textit{Bi}})$ correction. All these terms are encapsulated in the algebraic expression (5.12).

Figure 5.

Force $F$ on a cylinder in a uniform flow of Bingham fluid, against the Bingham number ${\textit{Bi}}$ . Points show values obtained from numerical simulations of the full problem (circles, current; triangles from Hewitt & Balmforth Reference Hewitt and Balmforth2018). The solid curve is the asymptotic prediction (5.12), while the dotted curve retains only the leading-order term (as given by Hewitt & Balmforth Reference Hewitt and Balmforth2018), and the dashed curve indicates the three-term expansion (5.14).

If we are interested, instead, in (say) the settling velocity $U_s$ , in response to an imposed (unit) driving force, then we write $U_s=1/F$ , ${\textit{Bi}}=F\, {\textit{Od}}$ and $A=F/4\pi +O({\textit{Bi}})$ in (5.11) to obtain

(5.15) \begin{equation} U_s = \frac {1}{4\pi }\left (\ln \!\left (\frac {1}{\textit{Od}}\right )-3.9\right ) + O({\textit{Od}}). \end{equation}

Figure 5 compares the force on a circular cylinder in a Bingham fluid, evaluated from full numerical simulations with the predictions of this asymptotic theory. Here, we plot the results from new simulations that we have performed, along with those from Hewitt & Balmforth (Reference Hewitt and Balmforth2018). This figure shows that the asymptotic solution (5.12) agrees very well with the numerically determined results. We have also plotted the asymptotic series (5.14), retaining the first term and the first three terms. The one-term approximation is somewhat different from the numerical results over the range that has been investigated ( ${\textit{Bi}}\gt 10^{-5}$ ). This is improved by the inclusion of the second and third terms in the series, but both fail as ${\textit{Bi}}$ approaches unity since the series are divergent; this difficulty is not inherited by (5.12). Coincidentally, taking ${\textit{Bi}}\gg 1$ in (5.12) gives $F\sim \exp (3.9)\,Bi$ , following the correct scaling (but different prefactor) to the large Bingham result predicted by plasticity theory, $F\sim 4(\pi +2\sqrt {2})\,Bi$ (Randolph & Houlsby Reference Randolph and Houlsby1984; Hewitt & Balmforth Reference Hewitt and Balmforth2018). This is, of course, just a coincidence, but it can be observed in figure 5(a) as the asymptotic solution lying approximately parallel to the numerical results at the larger Bingham numbers.

Figure 6.

(a) Example numerical simulation of the full problem for flow around a cylinder at ${\textit{Bi}}=1/32$ . The colour plot shows the strain rate on a logarithmic scale (with black representing unyielded regions), and the black dashed lines indicate a selection of streamlines in the frame of reference of the unyielded, far-field fluid. A quarter of the unit-radius cylinder is just visible in white at the bottom left corner of the plot. (b) Outer yield surface (plotted at $\tau =1.001\,Bi$ ) from numerical simulations for the two-dimensional viscoplastic Stokeslet problem (solid black) and simulations for the full problem at Bingham numbers shown in the legend. The latter have been rescaled to the outer coordinates, via (5.4). Also shown are the streamlines for the Stokeslet solution (black dashed lines). In both panels, the flow is in the anticlockwise direction.

Figure 6(a) shows an example of a numerical solution for the full problem of flow around a circular cylinder at ${\textit{Bi}}=2^{-5}$ . The colour plot of strain rate shows that the fluid is unyielded (coloured black) in the far field, and yielded in an envelope around the cylinder. Streamlines are also shown in the frame of reference of the unyielded fluid. Figure 6(b) shows the boundary of the yielded envelope around the cylinder for three such simulations, at different Bingham numbers, after rescaling to the outer coordinates via (5.4). Also shown is the yielded envelope from the two-dimensional viscoplastic Stokeslet calculation, with $\varepsilon =10^{-5}$ , evidencing how the outer problem captures the behaviour far from the cylinder. Again, the streamlines for the Stokeslet flow (figure 6 b) closely match those for the flow around the cylinder (figure 6 a). In this case, we find that the yielded envelope around the cylinder extends to a distance approximately 0.10 upstream (and downstream) of the cylinder, and approximately 0.17 to either side. Translating to unscaled dimensional coordinates, and substituting (5.11) for $A$ , this corresponds to distances

(5.16) \begin{equation} 0.10\times \frac {4\pi a}{Bi\,\mathrm{W} ({4\pi \exp (-3.9)}/{\textit{Bi}})}+O(1) \!\!\quad \text{and} \!\!\quad 0.17\times \frac {4\pi a}{Bi\,\mathrm{W}({4\pi \exp (-3.9)}/{\textit{Bi}})}+O(1), \end{equation}

where $a$ is the radius of the cylinder.

6. Flow past an elliptic cylinder

Finally, we analyse two-dimensional viscoplastic flow past an elliptic cylinder, with (after non-dimensionalising) unit semi-major axis aligned with the $x$ -axis ( $\theta =0$ ), and semi-minor axis of length $b\lt 1$ , aligned with the $y$ -axis ( $\theta =\pi /2$ ). Once again, the Bingham number is defined by (2.3), with $a$ now representing the dimensional semi-major axis. We assume that the far-field velocity is at an angle $\alpha$ to the $x$ -axis, giving a far-field streamfunction, in cylindrical coordinates,

(6.1) \begin{equation} \psi \to -r\sin (\theta -\alpha ) \quad \text{for }\ r\to \infty . \end{equation}

This flow is therefore characterised by three dimensionless parameters: ${\textit{Bi}}$ , $\alpha$ and $b$ , representing the Bingham number, the angle of inclination of the uniform flow, and the aspect ratio of the ellipse, respectively. Following Shintani et al. (Reference Shintani, Umemura and Takano1983) (see also Berry et al. Reference Berry, Swain and Richmond1923), Stokes flow past an ellipse is compactly derived using elliptic coordinates

(6.2) \begin{equation} x=\epsilon \cosh \xi \cos \eta , \quad y=\epsilon \sinh \xi \sin \eta , \end{equation}

where $\epsilon \equiv \sqrt {1-b^2}$ is the eccentricity of the ellipse, and the boundary of the ellipse is given by $\xi =\xi _0=\tanh ^{-1}b=\ln \sqrt {(1+b)/(1-b)}$ . The solution local to the ellipse is then

(6.3) \begin{align} \psi &= \epsilon \Bigl [{\sin \tilde {\alpha }}\cos \eta \cosh \xi \left \{\xi -\xi _0+\cosh ^2\xi _0\left (\tanh \xi _0 -\tanh \xi \right )\right \} \nonumber\\ &\quad -\,\cos \tilde {\alpha }\sin \eta \sinh \xi \left \{\xi -\xi _0 -\sinh ^2\xi _0\left (\coth \xi _0-\coth \xi \right )\right \}\Bigr ], \end{align}

where $A$ and $\tilde {\alpha }$ are now two constants that depend on ${\textit{Bi}}$ , $\alpha$ and $b$ . The resulting force on the ellipse (see Berry et al. Reference Berry, Swain and Richmond1923; Shintani et al. Reference Shintani, Umemura and Takano1983) is

(6.4) \begin{equation} \boldsymbol{F}=- 4 \pi A \left (\cos (\tilde {\alpha })\, \boldsymbol{e}_x + \sin (\tilde {\alpha })\,\boldsymbol{e}_y\right )\!. \end{equation}

In this way, $|\boldsymbol{F}|=4\pi A$ and the force is orientated at angle $\tilde {\alpha }$ to the $\theta =0$ axis. There is no torque on the ellipse (see Berry et al. Reference Berry, Swain and Richmond1923; Shintani et al. Reference Shintani, Umemura and Takano1983). In the far field, $\xi \gg \xi _0$ , we find that to leading order, $\xi =\ln (2r/\epsilon )+\cdots$ and $\eta = \theta +\cdots$ , which are accurate up to $O(1/r)$ . The far-field behaviour of the streamfunction is therefore

(6.5) \begin{align} \frac {\psi }{Ar} = -\sin (\theta -\tilde {\alpha })\ln \left (\frac {2r}{1+b}\right ) +\frac {1}{1+b}\left (b\cos \tilde {\alpha } \sin \theta - \sin \tilde {\alpha } \cos \theta \right )+\cdots . \end{align}

As for flow around a circular cylinder, the viscous and plastic stresses become comparable when $r=O(A/{\textit{Bi}})$ . Thus we scale the radial coordinate and the streamfunction exactly as for the circular cylinder (§ 5), but further make the rotation $\varTheta =\theta -\tilde {\alpha }$ , so that in the outer $(R,\varTheta )$ coordinates, the force due to the cylinder is aligned with the $x$ -axis. The rescaled problem for the correction to the uniform stream $\varPsi _0$ is then precisely the viscoplastic Stokeslet in two dimensions, for which we already have the behaviour at $R\ll 1$ , (5.8), giving (after the rotation back to $\theta =\varTheta +\tilde {\alpha }$ )

(6.6) \begin{equation} \frac {\psi }{Ar} \sim -\frac {1}{A}\sin (\theta -\alpha ) + \left (\ln \left (\frac {4\pi A}{Bi \,r}\right )-3.4\right )\sin (\theta -\tilde {\alpha })+\cdots . \end{equation}

Asymptotic matching then entails equating (6.5) and (6.6). We expand the trigonometric terms, and equate coefficients of $\sin \theta$ and $\cos \theta$ in (6.5) and (6.6) (noting that the $\ln (1/r)$ terms automatically cancel), giving two algebraic relations for $A(Bi,\alpha ,b)$ and $\tilde {\alpha }(Bi,\alpha ,b)$ , which are accurate up to, but not including, $O({\textit{Bi}})$ :

(6.7) \begin{align} \ln \left (\frac {1}{\textit{Bi}}\right ) + \ln \left (\frac {8\pi }{1+b}\right )-3.4- \frac {1}{1+b} &= \frac {\sin \alpha }{\sin \tilde {\alpha }}\frac {1}{A}+\ln \left (\frac {1}{A}\right )\!, \end{align}

(6.8) \begin{align} \ln \left (\frac {1}{\textit{Bi}}\right ) + \ln \left (\frac {8\pi }{1+b}\right )-3.4- \frac {b}{1+b} &= \frac {\cos \alpha }{\cos \tilde {\alpha }}\frac {1}{A}+\ln \left (\frac {1}{A}\right )\!. \end{align}

The predictions of (6.7)–(6.8) for the angle $\tilde {\alpha }$ and the magnitude of the force $F=4\pi A$ are shown in figure 7 for ${\textit{Bi}}=1/64$ , and a range of values of $b$ . The results of numerical simulations are also shown for $b=1$ and $b=0.2$ , showing close agreement.

Figure 7.

(a) Difference between the angle of the force $\tilde {\alpha }$ and the uniform flow angle $\alpha$ as a function of $\alpha$ for ${\textit{Bi}}=1/64$ at various values of the dimensionless semi-minor axis $b$ . The lines show the results of the asymptotic prediction (6.7)–(6.8) for $b$ between $0$ (top) and $1$ (bottom) in increments of $0.2$ . Numerical simulations of the full problem were computed for $b=0.2$ , and are shown as blue circles. The red circles show the numerical simulation for $b=1$ (this is a single simulation plotted at a number of different $\alpha$ values, since the circular cylinder has no dependence on $\alpha$ ). (b) As in (a), but showing the magnitude of the force on the elliptic cylinder. Again, the lines correspond to $b$ between $1$ (now top) and $0$ (now bottom) in increments of $0.2$ , and the blue and red circles are from numerical solution of the full problem with $b=0.2$ and $b=0$ , respectively.

Seeking asymptotic series for $A$ and $\tilde {\alpha }$ , and expanding up to the first term at which the series is affected by the aspect ratio, $b$ , gives

(6.9) \begin{align} A&= \frac {1}{\ln (1/{\textit{Bi}})}+\frac {\ln (\ln (1/{\textit{Bi}}))}{[\ln (1/{\textit{Bi}})]^2}\nonumber\\&\quad +\left (-\ln \!\left (\frac {8\pi }{1+b}\right )+3.4 + \frac {\sin ^2\alpha + b\cos ^2\alpha }{1+b}\right )\frac {1}{[\ln (1/{\textit{Bi}})]^2}+\cdots , \end{align}

(6.10) \begin{align} \tilde {\alpha }&=\alpha + \frac {1}{\ln (1/{\textit{Bi}})}\left (\frac {1-b}{1+b}\right )\sin \alpha \cos \alpha +\cdots , \end{align}

which reduce to the circular cylinder result, (5.14) (and $\tilde {\alpha }=\alpha$ ), when $b=1$ , as they must do. The limit $b\to 0$ is also well defined and corresponds to the flow around a flat plate. The difference between the drag on the cylinder and the drag on the plate, $\Delta F$ , is

(6.11) \begin{equation} \Delta F=|\boldsymbol{F}(b=1)|-|\boldsymbol{F}(b=0)|=\frac {2\pi }{[\ln (1/{\textit{Bi}})]^2}\left (\ln 4 + 1 - 2\sin ^2\alpha \right )\!, \end{equation}

so that the drag is always reduced, and is naturally minimised when the flow is aligned with the plate ( $\alpha =0$ ). This case, of a flat plate cutting through a viscoplastic fluid in the same direction as its length, has been studied in the context of the large Bingham number regime, and the resulting boundary layer structure that arises, by Oldroyd (Reference Oldroyd1947) (in the case of a semi-infinite plate) and Balmforth et al. (Reference Balmforth, Craster, Hewitt, Hormozi and Maleki2017) (for a finite plate). The results above provide an extension of this classical problem to the low Bingham number regime, as well as arbitrary flow direction.

Again, it is also of interest to consider the settling velocity in response to an imposed force. This is achieved, as before, by writing $U_s=1/F$ , and substituting for $A=F/4\pi +O({\textit{Bi}})=1/4\pi U_s + O({\textit{Bi}})$ and ${\textit{Bi}}=F\, {\textit{Od}} = {\textit{Od}}/U_s$ at the matching stage. We then view ${\textit{Od}}$ and $\tilde {\alpha }$ (the angle the force makes with the major axis) as the independent variables, and $U_s$ and $\alpha$ (the angle the velocity makes to the major axis) as dependent variables. The relations (6.7)–(6.8) become (neglecting $O({\textit{Od}})$ terms)

(6.12) \begin{align} \ln \left (\frac {U_s}{\textit{Od}}\right ) + \ln \left (\frac {8\pi }{1+b}\right )-3.4- \frac {1}{1+b} &= \frac {\sin \alpha }{\sin \tilde {\alpha }}4\pi U_s+\ln \left (4\pi U_s\right )\!, \end{align}

(6.13) \begin{align} \ln \left (\frac {U_s}{\textit{Od}}\right ) + \ln \left (\frac {8\pi }{1+b}\right )-3.4- \frac {b}{1+b} &= \frac {\cos \alpha }{\cos \tilde {\alpha }}{4\pi U_s}+\ln \left (4\pi U_s\right )\!, \end{align}

hence

(6.14) \begin{align} U_s\sin \alpha &= \frac {1}{4\pi }\left (\ln \left (\!\frac {2}{(1+b)\,{\textit{Od}}}\right ) -3.4- \frac {1}{1+b}\right )\sin \tilde {\alpha }, \end{align}

(6.15) \begin{align} U_s\cos \alpha &= \frac {1}{4\pi }\left (\ln \left (\!\frac {2}{(1+b)\,{\textit{Od}}}\right ) -3.4- \frac {b}{1+b}\right )\cos \tilde {\alpha }, \end{align}

representing a weak preference to settle in a direction parallel to the major axis of the ellipse. For the flat plate, $b=0$ , the ratio of the settling velocity when aligned with the driving force, to that when perpendicular to it, is thus approximately given by

(6.16) \begin{align} \frac {U_s(b=0;\alpha =0)}{U_s(b=0;\alpha =\pi /2)} &=\frac {\ln (2/{\textit{Od}})-3.4}{\ln (2/{\textit{Od}})-4.4}+\cdots \end{align}

(6.17) \begin{align} &= 1+\frac {1}{\ln (2/{\textit{Od}})}+\frac {4.4}{[\ln (2/{\textit{Od}})]^2}+\cdots . \end{align}