Flow of a viscoplastic fluid around a particle

Abstract

We study the force exerted by the uniform flow of a Bingham fluid around two- and three-dimensional particles in the regime of slow creeping flow and relatively weak yield stress. Matched asymptotic expansions are employed to couple a viscously dominated Stokes flow close to the particle with a far field in which the yield stress and viscous stresses are comparable. The far-field region is therefore modelled as a Bingham fluid driven by a point force at the origin (i.e. a viscoplastic Stokeslet). It features the full nonlinearity of the viscoplastic rheology, and its solution is computed through direct numerical simulation. Asymptotic matching then leads to a quasi-analytical expression for the drag force in terms of the dimensionless Bingham number ${\textit{Bi}}$, which measures the magnitude of the yield stress relatively to viscous effects at the particle scale. We deploy this methodology to determine the drag force on a sphere in three dimensions, and circular and elliptic cylinders in two dimensions, confirming our asymptotic predictions by comparison with full numerical simulations of the motion. We also generalise the three-dimensional result to arbitrary particles. The viscoplastic correction to the Newtonian drag in three dimensions scales as ${\textit{Bi}}^{1/2}$. In two dimensions, however, the effects of viscoplasticity are non-negligible at leading order. The drag varies with $[\ln (1/{\textit{Bi}})]^{-1}$, but this asymptotic result is only approached very slowly. Instead, an accurate representation of the drag is derived in terms of a single algebraic relation between the drag and the Bingham number.

1. Introduction

We analyse the uniform flow of a yield stress fluid around two- and three-dimensional particles, and compute the drag that is exerted on them. When the particles are sufficiently small, or the yield stress sufficiently weak – collectively here measured through a Bingham number, defined explicitly below, which is the ratio of the yield stress to the viscous stresses on the scale of the particle – we determine the drag force asymptotically. This entails matching between a far-field region in which the viscous and plastic stresses both contribute to the leading-order force balance, and a region around the particle in which significant viscous stresses are generated, and dominate the weak plastic stresses. The far-field region exhibits a yield surface, delineating unyielded fluid translating with the uniform far-field velocity, and a yielded envelope in which the oncoming material is deformed. The calculation leads to quasi-analytical expressions for the drag, akin to classical results for the weak inertial modification of viscous drag (Oseen 1927; Proudman & Pearson 1957), but in this case revealing the effects of a non-vanishing yield stress, expressed in terms of the Bingham number. The drag on particles in uniform flows, and conversely the sedimentation dynamics of particles in response to imposed forces, is a problem that has many applications in natural settings and industrial processes. For viscoplastic fluids, the yield stress can be exploited, for example in fractionation of pulp fibres in the paper-making industry (Madani et al. 2010), or to enhance the transport of coarse solids by pipeline (du Plessis & Ansley 1967).

The drag forces exerted by uniform flows of Newtonian fluids on spheres or cylinders have presented important problems in fluid dynamics (Oseen 1927; Batchelor 1967) and have spawned mathematical advances to enable their calculation. In the regime of low Reynolds number, here expressing the ratio of the inertia of the fluid motion to the viscous stresses on the particle, Proudman & Pearson (1957) have shown how to employ matched asymptotic expansions to couple an inner, ‘Stokes’ region close to the particle in which inertial effects are negligible to an outer, ‘Oseen’ region in which viscous and inertial effects are comparable. In particular, for the case of two-dimensional flow around a cylinder, this construction provides a resolution to the ‘Stokes paradox’, in which there is no solution to the uniform flow around a cylinder with vanishing Reynolds number (Hinch 1991). The treatment of an elliptic cylinder with axes not necessarily aligned with the oncoming uniform flow was provided by Shintani, Umemura & Takano (1983) (see also Berry, Swain & Richmond 1923; Hasimoto 1953; Imai & Lighthill 1954), who showed how the small Reynolds number results in corrections to the drag and lift forces on the ellipse. In particular, the force on the ellipse need not be aligned with the uniform flow direction. Similarly, in three dimensions, the treatment of ellipsoidal particles in a uniform flow at low Reynolds number has been considered by Oseen (1927) and Breach (1961). The methods of matched asymptotic expansions in the regime of small particle Reynolds numbers have also been applied to shear flows (Saffman 1965) and flows in channels (Hogg 1994), in each case leading to an inertial migration perpendicular to the flow due to the modified drag.

The structure of the flow field for fluids that exhibit a yield stress – viscoplastic fluids – is rather different from their Newtonian counterparts. The stress falling below the yield stress at large distances from the particle results in unyielded material, thus it is only sufficiently close to the stationary particle that deformation occurs. Several studies have concerned the flow fields and forces resulting from uniform flow of a Bingham fluid past spheres and cylinders. A number of experimental studies have reported the drag force (or terminal settling velocity) of objects in viscoplastic fluids (Valentik & Whitmore 1965; Ansley & Smith 1967; Brookes & Whitmore 1969; Atapattu, Chhabra & Uhlherr 1995; Tabuteau et al. 2007; Arabi & Sanders 2016). The region of deformed fluid around a sphere was reported by Atapattu et al. (1995), using flow visualisation techniques to provide quantitative measurements of the dimensions of the envelope. Modelling strategies have included the use of variational extrumum principles to bound and estimate the drag (Yoshioka, Adachi & Ishimura 1971; Adachi & Yoshioka 1973), numerical simulation employing a regularised rheology to ease the computation of the yield surface (Beris et al. 1985; Deglo De Besses et al. 2003; Tokpavi, Magnin & Jay 2008; Nirmalkar, Chhabra & Poole 2012), and numerical simulations employing the augmented-Lagrangian algorithm to accurately capture the yield surface (Roquet & Saramito 2003; Yu & Wachs 2007; Putz & Frigaard 2010; Chaparian & Frigaard 2017; Hewitt & Balmforth 2018; Iglesias et al. 2020). The most relevant to the current study of flow of Bingham fluid around spheres is by Beris et al. (1985), who carried out a large number of finite element numerical simulations using a regularised rheology, and provided an outline of the asymptotic approach to determine the force on the sphere in the low Bingham number regime. In particular, they developed a scaling argument that the correction to the Newtonian force depends upon the square root of the Bingham number ${\textit{Bi}}$ . They evidenced this dependence numerically, but did not evaluate directly the magnitude of this correction, as we do here. Two-dimensional flow around a circular cylinder in the low Bingham number regime was considered by Hewitt & Balmforth (2018), who argued that the force on the cylinder varies with $\ln (1/{\textit{Bi}})^{-1}$ , analogously to the classical result for a Newtonian fluid at low Reynolds number. The agreement between their numerical results and this leading-order scaling was fairly poor; here, we show that the asymptotic series decays very slowly with decreasing Bingham number, and we derive an alternative formulation.

In this study, we show how to evaluate the drag forces exerted by flows of Bingham fluid past both spheres and circular cylinders in the regime of relatively weak yield stress, using quasi-analytical methods that are built upon matched asymptotic expansions. Moreover, we demonstrate that this analytical framework can be readily extended to flows past ellipses and three-dimensional particles, demonstrating how the yield stress modifies both the magnitude and alignment of the induced drag. Our asymptotic predictions are tested against numerical simulations of the flow, which are computed using the augmented-Lagrangian algorithm and thus avoid regularisation of the rheological model.

Our study is structured as follows. First, in § 2, we define the general problem. In § 3, we consider uniform flow past a stationary spherical particle, revisiting the scaling anticipated by Beris et al. (1985), and we evaluate asymptotically the modification to the drag exerted on the sphere. We then generalise this result to an arbitrary three-dimensional particle in § 4. Our study also tackles two-dimensional flows around circular and elliptical cylinders in §§ 5 and 6, respectively. In these cases, the effects of plasticity enter the leading-order expression for the force, even when the yield stress is relatively small. This is a consequence of the Stokes paradox for slow viscous flow around a cylinder (Hinch 1991), which is resolved by our analysis. In § 6, we also illustrate how the magnitude and orientation of the force are affected by the aspect ratio of the ellipse and its orientation relative to the oncoming flow. We provide conclusions in § 7, and there is one appendix, providing details of the numerical methods employed.

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. 1985; Hewitt & Balmforth 2018). 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. (1985) 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 2017) and the flow in the neighbourhood of a stagnation point on a no-slip boundary (Taylor-West & Hogg 2025); 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 (2016), implemented in FEniCS (Alnæs et al. 2015). 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. (1985), 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. (1985), 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. (1985) (solid line), and experimental data from Ansley & Smith (1967) (blue circles) and Arabi & Sanders (2016) (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 (1967) and Arabi & Sanders (2016). 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 (1957) 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 1969, 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 1876; 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 (2018), 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 1957; Hinch 1991). 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 (2018) (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 (2018), and also diverges from the analogous resolution of the Stokes paradox by the introduction of weak inertial effects (Hinch 1991). 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. 2015) 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 (2018). 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 2018). The solid curve is the asymptotic prediction (5.12), while the dotted curve retains only the leading-order term (as given by Hewitt & Balmforth 2018), 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 (2018). 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 1984; Hewitt & Balmforth 2018). 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. (1983) (see also Berry et al. 1923), 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. 1923; Shintani et al. 1983) 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. 1923; Shintani et al. 1983). 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 (1947) (in the case of a semi-infinite plate) and Balmforth et al. (2017) (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}

7. Discussion and conclusions

We have studied the drag force on objects in a uniform flow of Bingham fluid, in the regime of low Bingham number. By analogy with classical approaches to include inertia in such flows, we employed matched asymptotic methods to connect an inner solution, satisfying Stokes flow and no-slip on the particle boundary, with an outer solution, where the geometry of the object can be ignored, but the full nonlinear viscoplastic flow equations need to be solved. The advantage of this approach is that the outer solution is parameter-free and need be calculated only once (in each case of two or three dimensions). The limiting behaviour of this solution as the origin is approached then provides the matching condition through which to fix the undetermined constants from the inner solution.

We first extended the argument of Beris et al. (1985) that the correction to the force on the sphere scales with $\sqrt {\textit{Bi}}$ , and used our numerical solution of the outer problem to further determine the constant of proportionality. In dimensional quantities, we find that the force on the sphere is

(7.1) \begin{equation} F = 6\pi \mu a U_0\left (1 + 11.4 \sqrt {\frac {\mu a \tau _Y}{6\pi U_0}}+\cdots \right )\!. \end{equation}

This result is based upon the effects of the yield stress influencing the drag force at an asymptotic order before the effects of inertia, which we demonstrated requires that

(7.2) \begin{equation} \frac {{\textit{Re}}^2}{\textit{Bi}}= \frac {\rho ^2 U_0^3 a}{\mu \tau _Y}\ll 1, \quad \text{as well as} \quad {\textit{Bi}}=\frac {a \tau _Y}{\mu U_0}\ll 1. \end{equation}

In principle, both criteria can be met for a sufficiently small particle, given other material properties. However, we have not found experimental data from this dynamical regime in order to test the result. We have also extended the analysis to a general particle in three dimensions, deriving results in terms of the mobility matrix of the particle, $ \unicode{x1D648}$ , in Newtonian Stokes flow. In particular, we showed how the drag anisotropy on a slender ellipsoid is enhanced by the inclusion of a weak yield stress.

For the case of two-dimensional flows around cylindrical objects, the asymptotic analysis is similar but resolves the pathology of the Stokes paradox for Newtonian fluids with logarithmic divergence of the purely viscous solution. Our matched asymptotic formulation leads to much improved expressions for the drag force on a circular cylinder. In fact, we found that the first term neglected by Hewitt & Balmforth (2018) is only an order $\ln (\ln (1/{\textit{Bi}}))/\ln (1/{\textit{Bi}})$ smaller than the leading-order term, and thus only decays to a tenth of the leading-order term for extremely small values, ${\textit{Bi}}\lesssim 3\times 10^{-16}$ . For an elliptic cylinder, we found that the drag force deviates from that of the circular cylinder at $O([\ln (1/{\textit{Bi}})]^{-2})$ , but that the angle between the directions of the force and the velocity scales as $[\ln (1/{\textit{Bi}})]^{-1}$ . We further calculated the dependence of both of these quantities on the aspect ratio of the ellipse and the direction of the flow. A particular limit of the elliptic cylinder is that of a flat plate, which follows immediately from our results.

Beyond experimental validation, there are a number of interesting directions for further work. We have considered a uniform far-field velocity field in an unbounded domain. A natural extension would be to consider the role of boundaries in the flow. In general, the existence of unyielded fluid in viscoplastic flows can have the effect of cloaking boundaries (e.g. see Balmforth & Hewitt 2025); however, in the low Bingham number regime considered here, the yielded envelope is large compared to the lengthscale of the particle, so the inclusion of boundary effects could be particularly relevant. This effect has been demonstrated and studied through numerical solutions for the case of a sphere by Blackery & Mitsoulis (1997) and a cylinder by Zisis & Mitsoulis (2002) and Roquet & Saramito (2003). Second, we have considered particles under the action of an imposed force, for which the outer solution is driven to leading order by a point force (i.e. a Stokeslet). In the absence of a net force, as for force-free swimming microorganisms, the dominant behaviour at large distances from the particle corresponds to faster decaying singularities, such as the force dipole (the ‘stresslet’) and force quadrupole (Batchelor 1970; Lauga & Powers 2009). Several studies (Supekar, Hewitt & Balmforth 2020; Eastham, Mohammadigoushki & Shoele 2022; Hewitt 2024) have considered the effect of a yield stress on such biologically motivated problems, and the methodology employed in this paper provides a strategy by which to investigate the low Bingham number regime. Due to the faster decay of the flow disturbance, the outer region would occur closer to the particle than in this study, and additional numerical simulations would be required for the different parameter-free problems arising in the outer region (e.g. the viscoplastic stresslet). Another direction would be to consider other viscoplastic constitutive laws, such as a Herschel–Bulkley model, or viscoelastoplastic models, though the additional nonlinearity of such rheological models presents an additional challenge.