1. Introduction
The study of drag forces on objects moving through fluid media is one of the canonical problems in fluid mechanics (Batchelor Reference Batchelor2000). Today, the study of low- Reynolds-number (
$\textit{Re}$
) drag coefficients is critical in fields such as biomedical engineering, nanotechnology, drilling muds and slurry transport (Elgaddafi et al. Reference Elgaddafi, Ahmed, George and Growcock2012). The focus of this study is on drag measurements on objects moving in yield-stress fluids, a class of non-Newtonian fluids characterised by a finite stress threshold that must be exceeded before they flow. Yield-stress fluids are commonly found in both natural and industrial settings, including in the food industry, biomedical applications and drilling operations (Balmforth, Frigaard & Ovarlez Reference Balmforth, Frigaard and Ovarlez2014).
The settling motion of a single sphere or cylinder in yield stress fluids has been studied fairly well both in experiments (Ansley & Smith Reference Ansley and Smith1967; Jossic & Magnin Reference Jossic and Magnin2001; Merkak, Jossic & Magnin Reference Merkak, Jossic and Magnin2006; Tabuteau et al. Reference Tabuteau, Coussot, de and John2007; Jossic & Magnin Reference Jossic and Magnin2009; Tokpavi et al. Reference Tokpavi, Jay, Magnin and Jossic2009; Elgaddafi et al. Reference Elgaddafi, Ahmed, George and Growcock2012; Holenberg et al. Reference Holenberg, Lavrenteva, Shavit and Nir2012; Ahonguio, Jossic & Magnin Reference Ahonguio, Jossic and Magnin2014; Murch et al. Reference Murch, Krishnan, Shaqfeh and Iaccarino2017) and numerical simulations (Blackery & Mitsoulis Reference Blackery and Mitsoulis1997; Tokpavi et al. Reference Tokpavi, Magnin and Jay2008, Reference Tokpavi, Jay, Magnin and Jossic2009). It has long been observed that when an object with a density different from the surrounding yield-stress fluid is introduced, a critical force threshold (yielding criterion) must be exceeded for the object to settle at a constant terminal velocity (Ansley & Smith Reference Ansley and Smith1967; Jossic & Magnin Reference Jossic and Magnin2001; Merkak et al. Reference Merkak, Jossic and Magnin2006; Tabuteau et al. Reference Tabuteau, Coussot, de and John2007; Jossic & Magnin Reference Jossic and Magnin2009; Ahonguio et al. Reference Ahonguio, Jossic and Magnin2014). This yielding criterion was suggested in early studies of Ansley & Smith (Reference Ansley and Smith1967), and first determined analytically and numerically by Beris et al. (Reference Beris, Tsamopoulos, Armstrong and Brown1985) for a sphere settling in a Bingham fluid model. Later studies examined the yielding limit with numerical simulations (Blackery & Mitsoulis Reference Blackery and Mitsoulis1997) and in experiments (Jossic & Magnin Reference Jossic and Magnin2001; Tabuteau et al. Reference Tabuteau, Coussot, de and John2007; Jossic & Magnin Reference Jossic and Magnin2009).
Beyond yielding limit, studies have investigated the settling dynamics of spherical and cylindrical objects in yield-stress fluids, and have assessed the drag coefficient (Jossic & Magnin Reference Jossic and Magnin2001; Deglo de Besses et al. Reference Deglo de Besses, Magnin and Jay2004; Mitsoulis Reference Mitsoulis2004; Chafe & de Bruyn Reference Chafe, de and John2005; Tokpavi et al. Reference Tokpavi, Jay, Magnin and Jossic2009; Jossic & Magnin Reference Jossic and Magnin2009; Ouattara et al. Reference Ouattara, Jay, Blésès and Magnin2018) and details of the flow field around the settling object both in experiments (Tokpavi et al. Reference Tokpavi, Jay, Magnin and Jossic2009; Holenberg et al. Reference Holenberg, Lavrenteva, Shavit and Nir2012; Ahonguio et al. Reference Ahonguio, Jossic and Magnin2014; Sgreva et al. Reference Sgreva, Davaille, Kumagai and Kurita2020) and numerical simulations (Blackery & Mitsoulis Reference Blackery and Mitsoulis1997; Mitsoulis & Galazoulas Reference Mitsoulis and Galazoulas2009; Chaparian & Frigaard Reference Chaparian and Frigaard2017; Koblitz, Lovett & Nikiforakis Reference Koblitz, Lovett and Nikiforakis2018). Experiments have mainly used polymeric gels, particularly those based on aqueous Carbopol solutions as model yield stress fluids (Jossic & Magnin Reference Jossic and Magnin2001; Chafe & de Bruyn Reference Chafe, de and John2005; Holenberg et al. Reference Holenberg, Lavrenteva, Shavit and Nir2012; Fraggedakis, Dimakopoulos & Tsamopoulos Reference Fraggedakis, Dimakopoulos and Tsamopoulos2016; Fraggedakis et al. Reference Fraggedakis, Dimakopoulos and Tsamopoulos2016).
While the flow of yield-stress fluids past stationary or translating particles is relatively well understood, scenarios involving the combined translational and rotational motion of a single object (e.g. a sphere or cylinder) in yield-stress fluids remain less explored (Hewitt & Balmforth Reference Hewitt and Balmforth2018; Balmforth & Hewitt Reference Balmforth and Hewitt2025). In a small portion of their study, Hewitt & Balmforth (Reference Hewitt and Balmforth2018) investigated the axial motion with rotation of a cylinder in a Bingham fluid model, and showed in the limit of high and low Bingham numbers, drag force could be analytically estimated. The combined rotation and sedimentation of objects in yield-stress fluids have practical significance, particularly in biological contexts. For example, Helicobacter pylori (H. pylori), a gram-negative bacterium and leading cause of gastric ulcers, penetrates gastric mucus by flagellar propulsion (Hooi et al. Reference Hooi2017; Li et al. Reference Li, Choi, Leung, Jiang, Graham and Leung2023). Experiments showed that under certain force conditions, H. pylori can become immobilised despite continuous flagellar and head rotation (Celli et al. Reference Celli2009). This result also hints at some critical yield limit for locomotion of H. pylori that is likely due to a balance between thrust generated by the tail and drag on the rotating head (Nazari, Shoele & Mohammadigoushki Reference Nazari, Shoele and Mohammadigoushki2023). Therefore, the drag on the rotating head is particularly critical in shaping the bacterium’s locomotion, and understanding the yielding criterion, drag coefficient and flow fields around rotating and translating objects is crucial for a potential therapeutic strategies involving H. pylori’s infection. These fluid dynamics features remain unmeasured for bodies undergoing simultaneous translation and rotation in yield-stress fluids.
This study aims to systematically determine the yielding criterion, drag coefficient and flow field characteristics around objects undergoing combined translational and rotational motion in a simple, non-thixotropic yield-stress fluid formulated with Carbopol-980. Two geometries, a sphere and a cylinder, with both smooth and roughened surfaces, are examined. Additionally, numerical simulations based on the Herschel–Bulkley fluid model are performed and compared with experimental results. This paper is organised as follows. Section 2 describes the experimental methods. Section 3 presents the numerical model and simulation details. In § 4, the relevant dimensionless numbers governing both the experiments and the numerical simulations are introduced and discussed. Section 5 is devoted to the validation of the numerical model. In § 6, the experimental results are presented and compared with the numerical simulations. Finally, the main conclusions of the study and perspectives for future work are provided.
2. Experiments
2.1. Material
The yield-stress fluid used in this study is prepared by mixing Carbopol-980 (sourced from Lubrizol) with triethanolamine. The concentration of Carbopol is maintained at 0.1 wt
$\,\%$
for all experiments reported in this paper, with a fixed triethanolamine-to-Carbopol ratio of 1.5. Additionally, a Newtonian fluid based on corn syrup is prepared to provide a basis for comparison with the results obtained in the yield-stress fluid. The settling objects used in this experiment include a cylinder and a sphere, fabricated using a three-dimensional (3-D) printer (FormLab 3B). The sphere used had a diameter of
$d=$
8.07 mm, while the cylinder measured
$d=$
5 mm in diameter and
$L=$
14 mm in length. The dimensions of these objects are selected to have the same effective volume. For the roughened objects, the surface roughness was designed via 3-D printing at an average of 250 μm (see supplementary figure S1 available at https://doi.org/10.1017/jfm.2026.11782). These objects are hollow, allowing for the insertion of additional weight. By adjusting the mass of the objects with weight, we could precisely control the gravitational force acting on them under specific experimental conditions.
2.2. Experimental set-up
Fluid rheology is characterised using an HR-10 TA rheometer coupled with a 40 mm cone and plate configuration, on which sandpaper was attached to the surface of the plate to avoid wall-slip. To measure yield criterion, the drag as well as to analyse the flow field around settling and rotating objects, a custom-built Helmholtz coil set-up, as shown in figure 1, is used. The Helmholtz coil system consists of two large circular coils connected to a DC power supply. When a DC current passes through the coils, it generates a uniform magnetic field in the central region. A motor applies a controlled angular velocity
$\varOmega$
, causing the Helmholtz coils to rotate. Each object is embedded with a permanent magnet and when placed in a chamber containing the fluid, it is rotated at an imposed rotational velocity by the rotating magnetic field of the Helmholtz coil. The settling speed is measured using particle tracking velocimetry (PTV). Additionally, the details of the flow field around the object are analysed using a two-dimensional particle image velocimetry (PIV) method. For experiments related to PIV, the fluid was seeded with 50 ppm particles with an average size of 60-micron particles (from Potters Industries). The flow is illuminated in two planes (r–
$\theta$
), named as the orthogonal plane from now on, and the r–z plane, which is noted as the flow plane from now on. Further details on the experimental set-up can be found in our previous publications (Wu et al. Reference Wu, Solano, Shoele and Mohammadigoushki2022; Nazari et al. Reference Nazari, Shoele and Mohammadigoushki2023).
A schematic of the rotational Helmholtz coil set-up used to investigate the sedimentation and rotational dynamics of cylindrical and spherical objects under a uniform magnetic field. Here,
$F_d$
and
$F_g$
denote the drag and the gravitational forces, respectively.
Figure 1. Long description
Panel A: A photograph of the experimental setup showing a rotational Helmholtz coil apparatus. The setup includes a camera positioned to capture the motion of objects within the coil. A mirror is placed to reflect the view of the objects. The coil is rotating, as indicated by the omega symbol and arrow. Panel B: A diagram illustrating the forces acting on a cylindrical object within the coil. The object experiences drag force (Fd) upwards and gravitational force (Fg) downwards. The diagram includes a coordinate system with axes labeled z, r, and theta. The object is shown rotating, as indicated by the omega symbol and arrow.
3. Numerical simulations
Here, we consider an exterior flow problem of a steady settling axisymmetric object with rotation in an infinite viscoplastic fluid given by the regularised Herschel–Bulkley model (Saramito Reference Saramito2016; Saramito & Wachs Reference Saramito and Wachs2017), written in the dimensional form as
\begin{equation} -\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma } + \boldsymbol{\nabla }\!p = \boldsymbol{0}, \quad \boldsymbol{\sigma } = 2 \eta (\dot {\gamma })\, \boldsymbol{D}(\boldsymbol{u}), \quad \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u} = 0, \end{equation}
where
\begin{align} \eta (\dot {\gamma }) = k \, \dot {\gamma }^{\,n-1} + \sigma _y \, (\dot {\gamma }^2 + \epsilon ^2)^{-1/2} \end{align}
is the effective viscosity which depends on the shear rate
$\dot {\gamma }=|2\mathbb{D}|$
. The parameter
$n \gt 0$
is the power-law index and
$k$
is the associated consistency index. In addition,
$p$
is the pressure, and
$\boldsymbol{\sigma }$
and
$\boldsymbol{D}$
are the deviatoric stress and rate-of-strain tensors, respectively. Moreover,
$\boldsymbol{u}$
is the fluid velocity and
$\epsilon$
is a small regularisation parameter. The tensor norm is defined as
$|\boldsymbol{A}|^2 = (\boldsymbol{A} : \boldsymbol{A})/2$
for any tensor
$\delta \in \mathbb{R}^{3\times 3}$
.
It is worth noting that the use of regularisation methods for viscoplastic flows is known to involve inherent mathematical limitations, particularly regarding the convergence of the stress field and the accurate localisation of rigid (unyielded) regions. However, the model is shown by Glowinski, Lions & Trémolières (Reference Glowinski, Lions and Trémolières1981) to exhibit monotonic convergence, with the velocity field approaching the exact solution in both the
$H^1(\varOmega )$
and
$L^\infty (\varOmega )$
norms, even when it results in a smooth transition between yielded and unyielded regions over a narrow band as opposed to sharply defined rigid zones. For problems with imposed kinematics, such as the swirl rotation considered here, viscoplastic fluid flow is dominated by strong boundary-induced shear near the particle surface. In such cases, regions where
$D(\boldsymbol{u}_\epsilon ) \to 0$
(
$\boldsymbol{u}_\epsilon$
is the estimated flow velocity for an assumed regularised parameter
$\epsilon$
) do not control the global dynamics, and the large quiescent zones typically observed in internal flows, such as Poiseuille or Couette configurations, are absent. Consequently, this model is found to be suited for the current study to examine velocity fields and flow topology, hydrodynamic drag, and the far-field flow structure as a function of the Bingham number.
A Galerkin finite element method (FEM) is used to numerically solve (3.1) (Eastham & Shoele Reference Eastham and Shoele2020). Its axisymmetric weak formulation at the
$n{\text{th}}$
iteration is
\begin{align} \int _D2\eta \big (\dot {\gamma }^{n-1}\big )\left [\mathbb{D}(\boldsymbol{u}^n):\mathbb{D}(\boldsymbol{v}) + \frac {u_r^n v_r}{r^2}\right ]&r\,\text{d}S - \int _D{p^n(\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{v}) r}\,\text{d}S= 0, \end{align}
\begin{align} \int _D q(\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}^n) r \,\text{d}S &= 0, \end{align}
where
$D$
is the axisymmetric domain,
$r$
is the radial distance,
$(\boldsymbol{u}^n, p^n)$
refers to unknown velocity and pressure fields to be solved at the
$n{\text{th}}$
iteration, while
$\dot {\gamma }^{n-1}$
is calculated from the previous iteration’s solution. Here,
$u_r$
denotes the radial component of the velocity vector. The test functions
$(\boldsymbol{v}, q)$
are chosen to be Taylor–Hood Q
$_2$
–Q
$_1$
shape functions. More details on the axisymmetric weak form for the Stokes equation and the validation of the model can be found from Eastham & Shoele (Reference Eastham and Shoele2020) and Eastham, Mohammadigoushki & Shoele (Reference Eastham, Mohammadigoushki and Shoele2022).
All viscoplastic simulations reported in this study were performed in a circular computational domain with an outer radius equal to 20 times the equivalent object radius, except in the model validation section, where a larger domain radius of 40
$R$
was employed. The smallest Bingham number considered in the validation study was
$\textit{Bi}$
= 0.01, whereas all other simulations were conducted for
$\textit{Bi}\geqslant$
0.1. Based on the theoretical predictions of Taylor-West & Hogg (Reference Taylor-West and Hogg2025), these domain sizes are sufficiently large to ensure that the yield surface remains fully contained within the computational domain. The computational set-up and boundary conditions, along with the mesh employed in this study, are shown in figure 2. Because the yield surface causes velocity field deformations to decay rapidly to negligible values, more elements can be concentrated near the object to better capture the yield surface. Additionally, the mesh is constructed to preserve discrete back–forward symmetry. Furthermore, to maintain consistency across all cases, simulations are performed for normalised equations by the particle radius and the maximum surface velocity, accounting for the combined effects of swirling rotation and sedimentation.
Computational domain and imposed boundary conditions for the axisymmetric solver. The domain is partitioned into three concentric regions to optimise mesh density: an inner domain (thickness
$0.15R$
) for high-resolution capture of the yield surface, a middle region (thickness
$5R$
) and an outer region extending to
$20R$
to minimise far-field effects. Here,
$R$
is the radius of the sphere, while a similar set-up is also used for the cylindrical case. Boundary conditions are defined as follows: at the object surface, a no-slip condition is applied, where
$u_z = U$
(sedimentation velocity),
$u_r = 0$
and
$u_{\theta } = r(z)\varOmega$
(rotational velocity). Along the axis of rotation, symmetry conditions are enforced (
$\partial u_z/\partial r = 0$
). Here,
$z$
and
$r$
denote axial and radial coordinates, while
$\theta$
represents the azimuthal angle.
We use a Picard iteration scheme to solve the nonlinear system (3.3). The procedure starts with assuming a yield stress
$\sigma _y$
and a relatively large regularisation, and the values of both parameters are continuously updated in the solution procedure towards their final values. Each step with a given
$\textit{Bi}$
and
$\epsilon$
also involves subiterations in which, at each iteration, the fluid equation (3.1) is interpreted as a variable-viscosity Stokes equation with the current estimate of the effective constitutive relation. The system is iteratively solved until convergence with an absolute difference tolerance of
$10^{-6}$
. The procedure is repeated to determine the next step with a higher
$\sigma _y{}$
and a smaller regularisation parameter. Overall, 10 logarithmic increments of
$\textit{Bi}$
from the starting point of
$\textit{Bi}=0.1$
and 10 logarithmic increments of
$\epsilon$
from the starting value of
$\epsilon =0.1$
are considered. The overall solution procedure continues until the final convergence criteria of
$10^{-6}$
is achieved.
In classical perfect viscoplasticity models, the yield surface is defined by the condition
$|\boldsymbol{\sigma }| = \sigma _y$
, with the rigid region corresponding to
$|\boldsymbol{\sigma }| \lt \sigma _y$
, in accordance with the perfect viscoplasticity assumption. In the present work, however, we adopt a definition of the yield surface that is more consistent with a regularised formulation, identifying the rigid region based on yield term contribution in the stress as
$\dot {\gamma }/\sqrt {\dot {\gamma }^2+\epsilon ^2}\,=\,0.999$
.
The drag force and swirl torque are calculated from the pressure and deviatoric stress defined in (3.3) as
\begin{align} F_d &=\int _S \quad \,\,\,\,\left [\left (-p\boldsymbol{I}+\tau (\dot {\gamma }) \right )\boldsymbol{\cdot }\boldsymbol{n}\right ] \boldsymbol{\cdot }\boldsymbol{e}_z \, \text{d}S, \end{align}
\begin{align} T_d &= \int _{S} \boldsymbol{r} \times [(-p\boldsymbol{I} + \boldsymbol{\tau }(\dot {\gamma })) \boldsymbol{\cdot }\boldsymbol{n}] \boldsymbol{\cdot }\boldsymbol{e}_z \, \text{d}S ,\end{align}
where
$S$
is the surface settling object and
$\boldsymbol{r}$
represents the position vector from the axis of rotation to a point on the surface of the object. A systematic validation of the model is presented in § 5.
4. Dimensionless numbers
When a particle translates at a constant terminal velocity in a yield-stress fluid and inertial effects are negligible, the plastic drag coefficient,
$C^*_d$
, is used to quantify the resistance encountered by objects, and is defined as (Jossic & Magnin Reference Jossic and Magnin2001, Reference Jossic and Magnin2009; Mitsoulis & Galazoulas Reference Mitsoulis and Galazoulas2009)
\begin{equation} C^*_d = \frac {F_d}{A \, \sigma _y}, \end{equation}
where
$A$
is the projected area of the falling object and is given as
$A = \pi d^2/4$
. In addition, in numerical simulations, it is often convenient to normalise the drag force using the Stokes drag for a shear-thinning Herschel–Bulkley fluid. We therefore introduce Stokes’ drag as
\begin{equation} C_s = \frac {F_d}{6 \pi \eta \textit{UR}}, \end{equation}
where
$\eta$
is the viscosity, defined as
$k(U/R)^{n-1}$
with
$U$
being the characteristic sedimentation velocity, and
$R$
the radius of cross-section of the object. Here,
$k$
and
$n$
are the consistency factor and shear-thinning index of the Herschel–Bulkley model, respectively. In addition, we define
$C_s^V$
and
$C_s^P$
for the contributions of pressure and shear stress in
$F_d$
according to (3.4). The dimensionless torque on a rotating cylinder is characterised using
\begin{equation} T_s = \frac {T_d}{2\pi \eta \textit{UR}^2}. \end{equation}
Another dimensionless group that appears in the literature is the gravitational yield number, defined as the ratio of the yield stress to buoyancy forces (Jossic & Magnin Reference Jossic and Magnin2001, Reference Jossic and Magnin2009; Mitsoulis & Galazoulas Reference Mitsoulis and Galazoulas2009):
\begin{equation} Y_G = \frac {\sigma _y}{g d_e\Delta \rho }, \end{equation}
where
$g$
is the gravitational acceleration, and
$\Delta \rho$
is the density difference between the object and the fluid. The density of the particle is calculated as
$\rho = m/V$
, where
$m$
and
$V$
are the particle mass and volume, respectively. Additionally,
$d_e = (6V/\pi )^{1/3}$
, and for the sphere,
$d_e = 2R$
; while for the cylinder,
$d_e$
is the equivalent diameter of a sphere of the same volume as the cylinder. The dimensionless gravitational yield number has been widely used (albeit in the absence of rotation) to characterise the stability (or yield limit) of falling objects in yield-stress materials. Specifically, the yield limit
$Y_{G,max }$
defines the threshold beyond which an object will not fall in a yield stress material (cessation of motion).
For an object that undergoes both translation and rotation, two distinct viscous stresses arise: (i) a translational viscous stress associated with sedimentation (
$U$
); and (ii) a rotational viscous stress generated by the imposed angular velocity (
$\varOmega$
). A meaningful characterisation of the flow therefore requires two independent dimensionless groups to differentiate and quantify the relative contributions of translation and rotation. The appropriate choice of dimensionless numbers depends on which quantities are externally imposed and which are measured. In the experiments, the effective weight of the object
$m$
(and, hence, the driving force for sedimentation) and the angular velocity
$\varOmega$
are prescribed, while the sedimentation velocity
$U$
is measured. In contrast, in the numerical simulations, both
$U$
and
$\varOmega$
are directly imposed. Since our objective is to compare experiments and simulations, the dimensionless groups must be defined in a consistent manner across both approaches. To achieve this consistency, the experiments were designed such that
$ U$
is effectively imposed in an indirect but controlled manner. In particular, for a given imposed rotational rate
$\varOmega$
, the object’s effective weight was adjusted iteratively until the target sedimentation velocity was attained (see supplementary figure S2). The experimental tests are repeated for at least three iterations to converge to the desired sedimentation velocity, wherein the results are obtained and reported. This procedure allows the experimental conditions to be mapped onto the same (
$U$
,
$\varOmega$
) parameter space as the simulations. To characterise the relative importance of translation, we therefore define a Bingham number that compares the material yield stress to the viscous stresses generated by translation as
\begin{equation} \textit{Bi} = \frac {\sigma _y}{k (U/L_{c})^n}, \end{equation}
where the characteristic shear rate is
$\dot {\gamma }= U/L_{c}$
and
$L_{c}$
is the characteristic length (sphere radius or cylinder length). To quantify the relative importance of rotation, we non-dimensionalise the rotational velocities using
\begin{equation} \hat {\varOmega } = \frac { R\varOmega }{U}. \end{equation}
Large values of
$\hat {\varOmega }$
correspond to rotation-dominated flows, while small values indicate translation-dominated behaviour.
Finally, a Reynolds number can be defined for this study as
$\textit{Re} = \rho U^{2-n}d^{n}/{k},$
where
$\rho$
is the particle density. The Reynolds number assessed for this study is very small (
$\textit{Re}_{{max}} \approx 5.77\times 10^{-3}$
) and, therefore, inertia does not affect the results.
5. Model validation
5.1. Regularisation parameter and mesh sensitivity
The mesh resolution and the regularisation parameter are systematically examined to ensure that the numerical model accurately captures the flow dynamics for the present configuration and across the intended parametric regimes. Three representative values of the Bingham number
$\textit{Bi}$
are selected for two distinct rotational velocities
$\hat {\varOmega }$
to assess the sensitivity of the solution (see figure 3). In these tests, the settling particle is spherical and the surrounding fluid follows a Herschel–Bulkley rheology with
$n=0.36$
(matched with experimental values), consistent with the parameters adopted throughout the remainder of this study.
The regularisation parameter
$\epsilon$
is progressively reduced from
$1$
to
$10^{-5}$
and the sensitivity of the Stokes’ drag to the choice of
$\epsilon$
is evaluated. For the base finite-element mesh (red curves in figure 3), the results are found to converge for
$\epsilon \lt 10^{-3}$
. The value of
$\epsilon =5\times 10^{-4}$
is selected for all the numerical tests conducted in this work. A similar convergence behaviour is observed when the mesh is refined such that the average element size is reduced to
$66\,\%$
of that of the base mesh, yielding nearly identical predictions. Based on this sensitivity analysis, the base mesh shown in figure 2 is adopted for all subsequent simulations.
Stokes’ drag
$C_s$
of a settling rotational sphere in a Herschel–Bulkely fluid model as a function of the regularisation parameter
$\epsilon$
for
$\textit{Bi}=2$
(
$\triangle$
),
$50$
(
$\square$
) and
$200$
(
$\circ$
) assuming different scaled rotational velocities of
$\hat {\varOmega } = 10, 10^{3}$
. Two meshes are compared: the base mesh shown in figure 2 (mesh 1 with red symbols in this figure) and a more refined mesh with isotropic refinement and an average mesh size of
$66\,\%$
of the base mesh, shown with the blue colour symbols. Here, only every fourth data point is shown with symbols.
5.2. A rotating sphere in a Newtonian fluid
The first validation case considers a rotating sphere in a Newtonian fluid with
$\hat {\varOmega }=1$
. In this regime, the swirl and translational components of the flow are decoupled, and the problem admits an analytical solution expressed in spherical coordinates
$(r,\theta ,\phi )$
, representing the radial, polar and azimuthal directions, respectively, centred on a sphere of radius
$R$
(Stokes Reference Stokes1851),
\begin{equation} \frac {u_r}{U} = -\cos \theta\! \left ( 1 - \frac {3R}{2r} + \frac {R^3}{2r^3} \right )\!, \quad \frac {u_\theta }{U} = \sin \theta\! \left ( 1 - \frac {3R}{4r} - \frac {R^3}{4r^3} \right )\!, \quad \frac {u_\phi }{U} = \frac {1}{\hat {\varOmega }}\!\left (\frac {R}{r}\right )^3 \! \sin \theta . \end{equation}
The top panel of figure 4(a) presents the contour plot of the swirl velocity component, and the bottom panel of figure 4(a) shows the norm of the difference between the analytical and computational velocity vectors. An excellent agreement is achieved between numerical results and those of analytical solution, with the numerical error remaining below
$5 \times 10^{-4}$
.
(a) Comparison of flow for rotating sphere with
$\hat {\varOmega } =1$
in Newtonian fluid against analytical solution. The colour shows the norm of the difference between the flow velocity vectors. (b) Comparison of shear stress on a long cylinder with tangential axial motion and swirl rotation against analytical results for a wide range of
$\textit{Bi}$
and
$\hat {\varOmega }$
values.
Figure 4. Long description
Panel A: A contour plot compares the flow for a rotating sphere in a Newtonian fluid against an analytical solution. The horizontal axis is labeled z/R and the vertical axis is labeled r/R. The color scale indicates the norm of the difference between the flow velocity vectors, ranging from 0 to 4 times 10 to the power of −4. Panel B: Two line graphs compare the shear stress on a long cylinder with tangential axial motion and swirl rotation against analytical results. The horizontal axis of the top graph is labeled 1/Ω̂ and the vertical axis is labeled τ_D/(μU/R). The horizontal axis of the bottom graph is labeled Bi and the vertical axis is labeled τ_D/(μU/R). The graphs include data points for analytical (blue circles) and numerical (red triangles) results, with different line styles representing different values of Bi and Ω̂.
5.3. A rotating infinite cylinder in a Bingham fluid model
The second validation is performed for a rotating and axially translating infinitely long cylinder in a Bingham viscoplastic fluid (with
$n=1$
). In this case, the rotation axis is perpendicular to direction of translation. This problem can be solved analytically as shown by Hewitt & Balmforth (Reference Hewitt and Balmforth2018), resulting in analytical flow velocity that can be written in the cylindrical coordinate (
$z, r, \phi$
) as
\begin{align} \begin{aligned} \frac {u_z}{U} &= \frac {r_{\!p} {\textit{Bi}}}{C}\left [C^2 \log \left (\frac {r_{\!p}}{r}\right )-1+G(r) \right ]\!,\quad \frac {u_r}{U}= 0, \\[3pt] \frac {u_\phi }{U} &= \frac {r\,\textit{Bi}}{2} \left \{ S\left (\frac {r_{\!p}^2}{r^2}-1\right ) +\ln \!\left [ \frac {(1+S)\left (G(r)-S\right )} {(1-S)\left (G(r)+S\right )} \right ] \right \}\!, \end{aligned} \end{align}
where
$G(r) \equiv \sqrt {C^2 ({r/r_{\!p}} )^2+S^2}$
. Here,
$r_{\!p}$
denotes the radius of the yielded surface and
$(C,S) = (\cos \varGamma , \sin \varGamma )$
. The parameter
$\varGamma$
, together with
$r_{\!p}$
, is determined by satisfying the boundary conditions, which yields the tangential frictional shear on the cylinder as
$(({\tau _D R})/({\mu U})) = r_{\!p}\, C\, {\textit{Bi}}.$
This solution is based on the analytical results presented by Hewitt & Balmforth (Reference Hewitt and Balmforth2018), with a necessary correction in the expression for
$u_\phi$
((3.6) of Hewitt & Balmforth (Reference Hewitt and Balmforth2018)). For the computational model, the same mesh resolution used throughout the paper is employed here. The infinite geometry is approximated using a cylinder with an aspect ratio of
$L/R=40$
and the reported frictional drag forces are calculated from the middle section of the cylinder to avoid end effects. Figure 4(b) shows the comparison between computational and analytical results of surface shear forces for a different range of
$\textit{Bi}$
and
$\hat {\varOmega }$
, where a good match is found for a wide range of conditions.
5.4. A non-rotating sphere in a Bingham fluid model
Finally, the computational technique is validated for a settling sphere in a Bingham fluid model with yield stress of
$\sigma _y$
and viscosity of
$\mu$
. Here,
$\epsilon$
is selected based on a convergence study to be
$0.0005$
. Figure 5(a) compares the predicted Stokes drag on a sphere from the current method with previously reported results across different
$\textit{Bi}$
ranges, while figure 5(b) shows a similar comparison across varying gravitational yield number defined in (4.4). The simulations are carried out in a semicircular domain of
$40R$
, using a gradually coarsening mesh with element sizes ranging from
$0.001R$
near the surface to
$0.2R$
in the far field. The close agreement confirms the validity of the present computational approach.
Dimensionless Stokes drag as a function of (a)
$\textit{Bi}$
number and (b) gravitational yield number. The simulations are compared with results of Beris et al. (Reference Beris, Tsamopoulos, Armstrong and Brown1985) and experiments of Ansley & Smith (Reference Ansley and Smith1967).
6. Results and discussion
6.1. Fluid characterisation
Figure 6 displays the steady-state flow curve of the yield-stress fluid measured over a broad range of applied shear rates. These flow curves were obtained by ramping shear rates up and down to assess potential thixotropic behaviour in the fluid. As shown in figure 6, the ramp-up and ramp-down curves overlap, indicating the absence of hysteresis and confirming that the fluid exhibits no measurable thixotropy. Furthermore, our measurements do not detect any first normal stress differences
$(N_1)$
across the entire range of shear rates, suggesting that the fluid behaves as a simple viscoplastic, non-thixotropic material within the range of applied shear rates. The steady-state flow curve data are well described by the Herschel–Bulkley fluid model (represented by the fitted red curves in figure 6, with the corresponding fitting parameters provided in the figure caption).
Steady shear stress
$\sigma$
as a function of shear rate
$\dot {\gamma }$
for the yield stress fluid. The experimental data are fitted to the Hershel–Bulkley fluid model
$\sigma = \sigma _y + k\dot {\gamma }^n$
. Here,
$\sigma _y$
= 5.2 [Pa] is the yield stress,
$k$
= 8.5 [Pa.s
$^{n}$
] is the consistency index and
$n$
= 0.36 is the shear-thinning index.
6.2. Drag coefficient
In the following, we present the drag coefficient measurements and calculations as functions of rotation rate and sedimentation velocity.
6.2.1. Effect of rotation
Figure 7 presents the measured plastic drag coefficient,
$C_d^*$
, for cylinders and spheres with varying surface roughness and
$\hat {\varOmega }$
numbers at a fixed Bingham number
$\textit{Bi} \approx 2.6$
for the cylinder and
$\textit{Bi} \approx 2.1$
for the sphere. Note that in the experiments reported in figure 7, across all imposed angular velocities,
$U$
is controlled and fixed (see the corresponding raw data for object’s trajectory in supplementary figure S2). Three clear trends emerge from these experiments. First,
$C_d^*$
decreases monotonically with
$\hat {\varOmega }$
for all surface conditions and shapes. Second, roughened objects consistently exhibit higher plastic drag coefficients than smooth ones at a fixed rotation rate. Finally, at high
$\hat {\varOmega }$
values (or high angular velocities), the difference in
$C_d^*$
between smooth and roughened objects becomes negligible.
Plastic drag coefficient as a function of dimensionless rotation rate
$\hat {\varOmega }$
for (a) cylinder and (b) sphere. The measured drag coefficients with a roughened surface are denoted as filled square and the smooth surface as filled circles. Smooth, roughened cylinders, smooth and roughed spheres are denoted as SC, RC, SS and RS, respectively. In addition, the hollow symbols are the effective
$\hat {\varOmega }_{\textit{eff}}$
estimated based on the PIV data when the effect of wall-slip is removed. The dashed curves show the results of numerical simulations.
To further elucidate the nature of the above-mentioned trends, we performed PIV in a plane orthogonal to the cylinder’s longest axis (or direction of sedimentation). Figure 8 presents the two-dimensional (2-D), time-averaged, normalised fluid velocity fields (
$\langle U_\theta \rangle /(R\varOmega )$
) around the cylinder at a fixed
$\textit{Bi} \approx 2.6$
and varying rotation rates, for both smooth and roughened cylinders. At low
$\hat {\varOmega }$
numbers, despite a weak cylinder rotation, the surrounding medium shows negligible plastic deformation for both surface conditions. Two plausible scenarios may account for this observation. First, the yield-stress material adjacent to the cylinder may remain entirely unyielded, behaving as a rigid solid that permits slip of the cylinder. Alternatively, a yielded zone may exist, but remain too small to be resolved within the spatial resolution of our measurements. Due to experimental constraints, we cannot conclusively distinguish between these two possibilities.
Two-dimensional averaged normalised velocity profile around the cylinder in an orthogonal plane for (a) smooth and (b) roughened cylinders at a fixed
$\textit{Bi} \approx 2.6$
and various
$\hat {\varOmega }$
values (0.04, 1 and 4 from left to the right). (c) One-dimensional normalised averaged rotational velocity measured for the yield stress fluid in the radial direction. Filled and open symbols denote roughened and smooth surfaces, respectively.
Figure 8. Long description
The figure contains three panels labeled (a), (b), and (c). Panel A and Panel B: Two sets of contour plots show the two-dimensional averaged normalized velocity profile around a cylinder in an orthogonal plane for smooth and roughened cylinders, respectively. The x-axis represents the radial distance (r) in millimeters, and the y-axis represents the vertical distance (r) in millimeters. Each set includes three subplots for different values (0.04, 1, and 4 from left to right). The color scale on the right indicates the normalized velocity (U theta / (R Omega)) ranging from 0 to 1. Panel C: A scatter plot shows the one-dimensional normalized averaged rotational velocity measured for the yield stress fluid in the radial direction. The x-axis represents the radial distance (r) in millimeters, and the y-axis represents the normalized rotational velocity ((U theta) / (R Omega)). Filled symbols denote roughened surfaces, while open symbols denote smooth surfaces. Different symbols represent different values (0.04, 1, and 4).
At higher
$\hat {\varOmega }$
numbers (
$\hat {\varOmega } = 1$
and
$4$
), the plastic deformation zone shifts outward from the cylinder surface and becomes clearly resolvable in our experiments. Moreover, as
$\hat {\varOmega }$
increases, this zone extends progressively farther into the surrounding material. These PIV results suggest that the variation in plastic drag coefficient with
$\hat {\varOmega }$
is directly related to the spatial extent of the plastic deformation zone surrounding the rotating body. As the rotation rate increases (or
$\hat {\varOmega }$
increases), a larger volume of material undergoes plastic deformation, which reduces the effective resistance of the surrounding medium. Consequently, a smaller drag force (or object mass) is required to maintain a fixed value of
$\textit{Bi}\approx 2.6$
. It is also evident from the PIV data that wall slip occurs in these experiments. We will come back to the effect of wall-slip later in the text.
Figure 7 also shows that the plastic drag coefficient for roughened surfaces is consistently higher than the measured values for the smooth surface. This result is consistent with the experiments of Jossic & Magnin (Reference Jossic and Magnin2001, Reference Jossic and Magnin2009), who also observed higher plastic drag coefficients for roughened cylinders relative to smooth ones at
$\hat {\varOmega } = 0$
. PIV measurements in figure 8(c) further indicate that, at a fixed
$\hat {\varOmega }$
, the extent of the yielded zone does not differ markedly between roughened and smooth surfaces. Instead, the key distinction lies in the degree of wall slip: it is reduced for roughened surfaces compared with smooth ones. These findings suggest that wall slip effectively reduces plastic drag in yield-stress fluids. This result is consistent with predictions of Supekar, Hewitt & Balmforth (Reference Supekar, Hewitt and Balmforth2020), who showed that for a Bingham plastic fluid, the wall-slip reduces the plastic drag coefficient of an infinitely long cylinder that sediments laterally. The plastic drag coefficient for the sphere exhibits trends similar to those observed for the cylinder, although the differences between roughened and smooth surfaces are somewhat less pronounced. In addition, the plastic drag coefficient on spheres is, in general, smaller than that of cylinders. Corresponding PIV results for the sphere in the orthogonal plane are provided in supplementary figure S3.
To further elucidate the above-mentioned experimental trends, we performed numerical simulations using the Herschel–Bulkley viscoplastic fluid model with parameters matching the experimental conditions, including object dimensions, the rheological properties of the surrounding fluid and imposed translational velocity. The results, shown as dashed curves in figure 7, generally capture the experimental trend of increasing plastic drag coefficient with
$\hat {\varOmega }$
. However, the simulations consistently underestimate the plastic drag coefficients for both the cylinder and the sphere. Previous numerical simulations with the Herschel–Bulkley viscoplastic fluid model also showed that the numerically predicted plastic drag coefficient
$C^*_d$
underpredicts the experimental measurements for non-rotating cylinders in line with our results reported here (Mitsoulis & Galazoulas Reference Mitsoulis and Galazoulas2009; Tokpavi et al. Reference Tokpavi, Jay, Magnin and Jossic2009). Then, the natural question that arises is why simulations underpredict plastic drag coefficients? Several factors may contribute to this discrepancy. One likely cause is the presence of wall slip in the experiments, which is absent in the simulations of this study. To evaluate the effect of wall-slip, we calculated an effective rotation rate
$\hat {\varOmega }_{\textit{eff}}$
for the experiments based on the rotation rate of the first fluid element adjacent to the object surface. For example, at an imposed
$\hat {\varOmega } = 1$
and roughened cylinder, the effective rotational velocity of the fluid is
$\varOmega _{\textit{eff}} \approx 0.2 \varOmega$
, and therefore,
$\hat {\varOmega }_{\textit{eff}}\approx 0.2$
. Here,
$\varOmega$
denotes the imposed rotation of the solid body. This shifts the experimental data to the left in figure 7 (depicted as empty symbols), thereby removing the possible influence of wall slip on the drag coefficient comparison. Interestingly, when the effect of wall slip is removed, the agreement between experiments and simulations improves substantially. At lower
$\hat {\varOmega }$
values, however, discrepancies remain, with simulations consistently under-predicting the experimentally measured drag coefficients for both spheres and the roughened cylinder. Another subtle consequence of the wall-slip in experiments is that in experiments where significant slip occurs, the yielded (plastic) region is observed to be thinner (or closer to the object surface) than predicted numerically (see predictions of numerical simulations in supplementary figure S4). As a result, the stress gradients in the vicinity of the rotating object can be significantly altered compared with simulations that enforce a no-slip boundary condition. This modification of the stress field can change the effective (force–velocity) relationship, which may lead to differences in the plastic drag coefficients.
Beyond wall-slip, other factors may affect the drag in experiments in ways that are not captured by the viscoplastic simulations. Figures 9(a) and 9(b) show the averaged 2-D velocity vectors and velocity magnitude around the falling cylinder. In the absence of any rotation (
$\hat {\varOmega } = 0$
), the velocity field is asymmetric around the cylinder with the upstream flow clearly showing a stagnation point flow accompanied by a negative wake (see also the averaged velocity at the centreline of the cylinder in figure 9
c). As the rotational velocity increases, the stagnation point flow becomes less significant and eventually disappears at the highest
$\hat {\varOmega }$
number. Changing the surface roughness does not significantly affect the details of the flow field around the cylinder at low
$\hat \varOmega$
numbers. A similar trend is observed for flow past the smooth and roughened spheres, as shown in supplementary figure S5. In contrast, the counterpart results in the Newtonian fluid (corn syrup) show a symmetric flow field, highlighting how the non-Newtonian behaviour leads to asymmetries in the flow profiles (see supplementary figures S6−S7).
Two-dimensional and 1-D averaged normalised velocity profile in flow plane (
$v$
is fluid velocity component in z direction and
$U$
is the sedimentation velocity) of the (a) smooth and (b) roughened cylinder at a fixed
$\textit{Bi}\approx 2.6$
and various
$\hat {\varOmega }$
values: (i) 4; (ii) 1; (iii) 0.04 and (iv) 0). The averaged normalised velocity along the axis of the cylinder (
$r=0$
) for (c) smooth and (d) roughened cylinder.
It is worth noting that numerical simulations using the Herschel–Bulkley viscoplastic fluid model do not capture the flow asymmetry observed experimentally around the falling object. Interestingly, as the negative wake feature is dampened in the experiments at higher
$\hat \varOmega$
numbers, the discrepancy between drag coefficients of simulations and experiments diminishes (see empty symbols corrected for the effect of wall-slip in figure 7). This trend suggests a strong correlation between the formation of flow asymmetry (particularly the negative wake) around the falling object and the resulting drag coefficient. A similar asymmetry in the fore-and-aft direction has been reported for a sphere settling in yield stress fluids based on carbopol (Putz et al. Reference Putz, Burghelea, Frigaard and Martinez2008; Tokpavi et al. Reference Tokpavi, Jay, Magnin and Jossic2009; Holenberg et al. Reference Holenberg, Lavrenteva, Shavit and Nir2012; Sgreva et al. Reference Sgreva, Davaille, Kumagai and Kurita2020). This flow phenomenon has been suggested to be associated with either the viscoelasticity nature of the yield stress fluid or presence of ageing or thixotropy (Putz et al. Reference Putz, Burghelea, Frigaard and Martinez2008; Tokpavi et al. Reference Tokpavi, Jay, Magnin and Jossic2009; Sgreva et al. Reference Sgreva, Davaille, Kumagai and Kurita2020; Holenberg et al. Reference Holenberg, Lavrenteva, Shavit and Nir2012). Recent finite element simulations using elasto-viscoplastic (EVP) fluid models predicted the presence of fore–aft asymmetry around a falling sphere and showed that yield stress significantly enhances the negative wake behind the objects in an EVP fluid model (Fraggedakis et al. Reference Fraggedakis, Dimakopoulos and Tsamopoulos2016; Yazdi & D’Avino Reference Yazdi and D’Avino2023). Following the preparation steps noted by Dinkgreve et al. (Reference Dinkgreve, Fazilati, Denn and Bonn2018), the flow curves of the yield stress fluid were characterised both via shear rate ramp up followed by a shear rate ramp down, and showed no sign of thixotropy or hysteresis. In addition, the nonlinear rheological measurements of the yield stress fluid used in this study do not show any signs of normal stress differences or viscoelasticity, which suggests that the plastically deformed fluid around the object is unlikely to form normal stresses associated with viscoelasticity. These rheological characterisations further support the expectation that no fore–aft asymmetry should arise in the flow field for the fluid used in this study. Nevertheless, the origin of the reported negative wake around the object in this work remains unclear. Despite this, prior numerical simulations based on the EVP fluid model have suggested that, at fixed terminal velocity or Stokes’ drag, the plastic contribution to the drag decreases with increasing viscoelasticity of the EVP fluid model (Fraggedakis et al. Reference Fraggedakis, Dimakopoulos and Tsamopoulos2016). In this study, an EVP fluid model based on Oldroyd-type formulations was used and it is possible that drag reduction was a consequence of that model choice. It is worth noting that experimental evidence isolating the role of elasticity on drag in yield-stress fluids, particularly in the presence of imposed rotation, is currently lacking. This is the subject of our future work.
Simulation results for (a) normalised drag forces
$C_s$
, (b) ratio of the pressure drag to shear-induced drag, (c)
$C^*_d$
and (d) normalised torque as a function of
$\textit{Bi}$
numbers for different
$\hat {\varOmega }$
. The continuous and dashed curves correspond to simulations on the sphere and the cylinder, respectively. The arrows show the results for non-rotating objects at different
$\textit{Bi}$
values.
Figure 10. Long description
Panel A: A line graph shows the normalized drag forces as a function of dimensionless numbers for different Bingham numbers. The x-axis represents the dimensionless numbers, and the y-axis represents the normalized drag forces. Continuous and dashed curves correspond to simulations on the sphere and the cylinder, respectively. The arrows indicate the results for non-rotating objects at different values. Panel B: A line graph shows the ratio of the pressure drag to shear-induced drag as a function of dimensionless numbers for different Bingham numbers. The x-axis represents the dimensionless numbers, and the y-axis represents the ratio of the pressure drag to shear-induced drag. Continuous and dashed curves correspond to simulations on the sphere and the cylinder, respectively. Panel C: A line graph shows the normalized torque as a function of dimensionless numbers for different Bingham numbers. The x-axis represents the dimensionless numbers, and the y-axis represents the normalized torque. Continuous and dashed curves correspond to simulations on the sphere and the cylinder, respectively. The arrows indicate the results for non-rotating objects at different values. Panel D: A line graph shows the normalized torque multiplied by dimensionless numbers as a function of dimensionless numbers for different Bingham numbers. The x-axis represents the dimensionless numbers, and the y-axis represents the normalized torque multiplied by dimensionless numbers. Continuous and dashed curves correspond to simulations on the sphere and the cylinder, respectively.
To gain deeper insights into the dependence of the drag coefficients on the rotation, we carried out a systematic series of numerical simulations. Figure 10 presents the variation of drag force across a wide range of rotational conditions represented by
$\hat {\varOmega }$
. The line colours indicate different values of
$\textit{Bi}$
, with solid lines corresponding to the sphere and dashed lines to the cylinder. Two drag coefficients,
$C_s$
and
$C^*_d$
, are applied to illustrate the changes in drag force in figures 10(a) and 10(c). Additionally, figure 10(b) isolates the contribution of pressure versus viscous frictional forces, while figure 10(d) displays the induced torque on the rotating object. For each curve, the corresponding results for non-rotating objects are indicated with arrows on the rightmost side of the figures.
The most pronounced variations in drag force occur at large
$\hat {\varOmega }$
numbers, where azimuthal rotation controls the formation and shape of the yielded flow region around the rotating object. In this regime, both
$C^*_d$
and
$C_s$
exhibit a clear power-law dependence on
$\hat {\varOmega }$
, with a slope that strongly depends on
$\textit{Bi}$
. At smaller
$\hat {\varOmega }$
, the scaling approaches an approximately constant value of
$C_s \sim \hat {\varOmega }^{-0.1}$
, except for cases with thick yield regions associated with
$\textit{Bi}\lt 1$
. Therefore, for rotational rates small relative to the sedimentation velocity, the drag coefficient is a weak function of
$\hat {\varOmega }$
, consistent with the limited experimental observations reported in figure 7. A notable feature is the transition to a stronger power-law dependency for all cases, as shown in figure 10(a). The transition between the two regimes shifts to larger
$\hat {\varOmega }$
as
$\textit{Bi}$
increases. The intersection points of the transition at small and large
$\hat {\varOmega }$
are well approximated by a curve with slope
$1.64$
and is shown by the black line in figure 10(a).
Figure 10(b) presents the ratio of the pressure contribution (
$C^P_s$
) to the viscous contribution (
$C^V_s$
) in the the drag force, i.e.
$\int _S -p\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{e}_z\,\mathrm{d}S$
relative to
$\int _S (\boldsymbol{\tau }\boldsymbol{\cdot }\boldsymbol{n})\boldsymbol{\cdot }\boldsymbol{e}_z\,\mathrm{d}S$
, as a function of
$\hat {\varOmega }$
. The pressure-to-viscous drag ratio reaches a maximum near the intersection of the two
$C_s$
power-law regimes discussed earlier, attaining values of approximately
$28$
for
$\textit{Bi}=1000$
. The variation of this ratio with
$\hat {\varOmega }$
is closely linked to the shape of the yielded region surrounding the settling body and will be discussed later. This observation is also aligned with the scaling argument developed in Appendix A, where we discuss a physical interpretation of this trend in the high-swirl limit. When the surface traction becomes primarily azimuthal, only its meridional projection contributes to axial drag, giving
\begin{equation} C^{*}_{d,\tau } \sim C_{\tau }\big (1+\hat {\varOmega }^{2}\big )^{-1/2}. \end{equation}
A confined yielded layer can also generate a pressure contribution,
\begin{equation} C^{*}_{d,p} \sim C_{\!p} \textit{Bi}^{2/(n+1)} \hat {\varOmega }^{-(3n+1)/(n+1)}, \end{equation}
which reduces to
$C_{\!p} \textit{Bi}/\hat {\varOmega }^{2}$
for a Bingham material. Thus, the pressure-to-viscous drag ratio in figure 10(b) is controlled not only by the magnitude of rotation, but also by how rotation reshapes and confines the yielded region.
Our simulations suggest that scalings of plastic drag coefficient
$C^*_d$
(as shown in figure 10
c) with respect to
$\hat {\varOmega }$
follow similar trends to those reported for
$C_s$
. In the small
$\hat {\varOmega }$
regime, the plastic drag coefficient is only marginally affected by the object rotating condition and results effectively approach those of non-rotating objects below
$\hat {\varOmega } \lt 0.1$
. The experimental results of figure 7 that are performed at a fixed Bingham number of
$\textit{Bi}\approx 2.6$
and
$0\lt \hat {\varOmega } \lt 5$
suggest that
$C_d^* \sim \hat {\varOmega }^{-0.24}$
for the roughened cylinder and
$C_d^* \sim \hat {\varOmega }^{-0.22}$
for the smooth cylinder over a range of
$0.5\leqslant \hat {\varOmega } \leqslant 5$
. In addition,
$C_d^* \sim \hat {\varOmega }^{0.18}$
for roughened sphere and
$C_d^* \sim \hat {\varOmega }^{0.14}$
for smooth sphere in the range of
$ 0.2\leqslant \hat {\varOmega } \leqslant 8$
. These scalings are within the range predicted by the numerical simulations.
Finally, figure 10(d) shows the variation of the imposed torque with
$\hat {\varOmega }$
through the quantity
$T_s/\hat {\varOmega }$
. Similar trends are observed for both the sphere and the cylinder. In the rotation-dominated regime, the torque decreases with
$\hat {\varOmega }$
following a power-law behaviour for all tested
$\textit{Bi}$
, where the curves become steeper with increasing
$\textit{Bi}$
. This approximately linear power-law trend starts from
$\hat {\varOmega }\ge 1$
, after which the torque/rotation coefficient continues to decrease. The deviation from power law occurs at low rotational rates where the translational motion is the primary factor that controls the yielded region near the object.
To elucidate the origin of the transition between these two distinct regimes reported in figure 10, we plot the shear rate
$\dot {\gamma }$
and pressure distributions at a fixed
$\textit{Bi} = 1.0$
and over a broad range of
$\hat {\varOmega }$
in figure 11. In these visualisations, the sedimentation direction is from right to left (indicated by the arrow). At low
$\hat {\varOmega }$
, a large yielded region emerges at the front and rear of the object. However, as
$\hat {\varOmega }$
increases, rotation induces lateral expansion of the yielded zone, consistent with our PIV measurements. This modification of the yielded region shape results in a sudden modification of the local pressure field.
This can be seen clearly in
$\hat {\varOmega } \approx 10^{1/3}$
, wherein the expanded midplane yield zone becomes sufficiently large to allow flow circulation from the front to the back of the body. Therefore, the pressure difference between the front and back of the object reduces (figure 11
b). The modification of the pressure field at the front and rear of the object explains the rapid decline in total drag observed at higher rotation rates. Overall, the pressure field here is strongly linked to the shape of the yielded region, similar to previous observation by Eastham et al. (Reference Eastham, Mohammadigoushki and Shoele2022). Whether the yield region is controlled by the translational motion or the swirling rotation, the pressure contribution to the drag changes nature and that results in a different power law shown in figure 10.
(a) Shear rate
$\dot {\gamma }$
around the sphere and cylinder at various
$\hat {\varOmega }$
numbers and at a fixed
$\textit{Bi} = 1$
. (b) Pressure field around spherical body for different
$\hat {\varOmega }$
.
6.2.2. Effect of sedimentation speed
Figure 12 presents the variation of the measured drag coefficient
$C^*_d$
, as a function of
$\textit{Bi}$
for different imposed rotational velocities with panels (a) and (b) representing smooth and roughened cylinders and panels (d) and (e) denoting smooth and roughened spheres, respectively. First, at a fixed rotation rate, as
$\textit{Bi}$
increases, the measured drag coefficient decreases with the drag coefficient being higher for roughened surfaces than the smooth ones. These experimental trends are consistent with existing data reported in the literature when objects do not undergo rotation (or
$\hat {\varOmega }=0$
). Jossic & Magnin (Reference Jossic and Magnin2001, Reference Jossic and Magnin2009) conducted an experimental study on the settling of smooth and roughened spheres and cylinders in yield stress fluids based on Carbopol 940 and reported that
$C^*_d$
follows the same trend as a function of the
$\textit{Bi}$
. Moreover, they observed that
$C^*_d$
is consistently higher for roughened cylinders compared with smooth ones. A similar observation, with a decreasing trend of
$C^*_d$
at higher
$\textit{Bi}$
, was reported by Merkak et al. (Reference Merkak, Jossic and Magnin2006) for non-rotating spheres and cylinders by Tokpavi et al. (Reference Tokpavi, Jay, Magnin and Jossic2009).
Second, at a fixed
$\textit{Bi}$
, our results indicate that rotation produces a lower drag coefficient for all surface types and geometries. Note, to keep
$\textit{Bi}$
fixed and vary
$\varOmega$
, the mass of the object is systematically and carefully adjusted to maintain a constant sedimentation velocity. The observed reduction in drag with increasing rotation rate is due to the development of a plastic deformation zone that expands around the object as rotation increases. The enlargement of this plastically deformed zone should reduce the overall resistance to motion.
Drag coefficient as a function of
$\textit{Bi}$
for different imposed
$\varOmega$
: (a) smooth cylinder; (b) roughened cylinder; (c) smooth sphere and (d) roughened sphere. Symbols denote experimentally measured values and curves denote the results of numerical simulations.
We also performed numerical simulations to evaluate the drag coefficient as a function of
$\textit{Bi}$
for all cases reported in figure 12 as solid curves. In the absence of rotation, the simulations overpredict the drag coefficient for the smooth cylinder, whereas the agreement for the roughened cylinder is noticeably better. For the sphere, measurements without rotation could not be obtained because the required mass exceeded the practical limits imposed by the sphere’s size. Once rotation is introduced, the simulations systematically underpredict the experimentally measured drag coefficients as discussed in earlier results as well. To further assess the influence of the Bingham number on the negative wake and the surrounding flow field, we conducted PIV measurements in the flow plane at
$\textit{Bi} \approx 6$
and 5 (constant
$U\approx 0.025$
[mm s−1]) for both cylinders and spheres. Figure 13 presents the 2-D and 1-D averaged velocity profiles measured around the cylinders. Interestingly, although the flow field remains asymmetric, the negative wake vanishes even in the absence of rotation. Moreover, increasing the rotation rate does not significantly alter the flow field, regardless of whether the cylinders are smooth or roughened. These findings suggest that the emergence of a negative wake is strongly dependent on the intensity of the flow around the object. In contrast, the flow around spheres at the same
$\textit{Bi}$
still exhibits a negative wake, albeit weaker than that observed at lower Bingham numbers (cf. supplementary figures S5 and S8).
Two-dimensional and one-dimensional averaged normalised velocity profiles in the flow plane for cylinders. (a) Smooth and (b) roughened cylinders at a constant Bingham number (
$\textit{Bi}\approx 6$
) and
$\hat \varOmega = (i)\,40$
; (ii) 10; (iii) 0.4 and (iv) 0. (c) and (d) Corresponding averaged normalised velocity profiles along the
$z$
axis for smooth and roughened cylinders, respectively.
Figure 13. Long description
The image contains four panels showing velocity profiles for smooth and roughened cylinders in yield-stress fluids. Panel A: Four vertical bar graphs display two-dimensional normalized velocity profiles in the flow plane for smooth cylinders at different Bingham numbers. The x-axis represents the radial distance (r) in millimeters, and the y-axis represents the vertical distance (z) in millimeters. The color scale indicates the normalized velocity (v/U). Panel B: Four vertical bar graphs show two-dimensional normalized velocity profiles in the flow plane for roughened cylinders at different Bingham numbers. The axes and color scale are the same as in Panel A. Panel C: One line graph presents the averaged normalized velocity profiles along the axis for smooth cylinders. The x-axis represents the vertical distance (z) in millimeters, and the y-axis represents the normalized velocity (v/U). Different line styles and symbols represent different Bingham numbers. Panel D: One line graph shows the averaged normalized velocity profiles along the axis for roughened cylinders. The axes and line styles are the same as in Panel C.
To evaluate the sensitivity of the drag coefficient to
$\textit{Bi}$
, a systematic numerical simulation is performed. Figure 14 shows how
$C_s$
,
$C_s^P/C_s^V$
,
$C^*_d$
and
$T_s$
change as a function of
$\textit{Bi}$
for a wide range of
$\hat {\varOmega }$
numbers. A similar trend is observed in both the sphere and cylinder cases.
We begin with the behaviour of Stokes’ drag. For cases with large
$\textit{Bi}$
, it is noted that for all
$\hat {\varOmega } \lt 0.1$
, the Stokes drag coefficient (
$C_s$
) is best scaled by sedimentation velocity and viscosity, showing almost no dependence on the rotational state (see figure 14
a). Under these conditions,
$C_s$
exhibits a slope of approximately
$C_s\sim \textit{Bi}$
. These cases correspond to translation-dominated flow fields where the shear rate is strongly localised and the pressure field exhibits a higher value at the frontal and back stagnation regions, resulting in larger plastically deformed (yielded) zones at the front and back of the object, and only a narrow yielded layer along its flanks (see figure 15
a). For larger values of
$\hat {\varOmega }$
, we observe a
$C_s\sim \textit{Bi}^{\alpha }$
, where
$\alpha \gt 1$
, and
$\alpha$
increases slightly with
$\hat {\varOmega }$
. In this regime, the azimuthal velocity induced by rotation generates strong shear near the body surface and expands the surrounding yielded region, allowing the in-plane velocity to circulate more freely around the object (figures 15
b and 15
c). In contrast, at smaller
$\textit{Bi}$
, the variation of
$C_s$
is relatively weak, since the rotation largely controls the effective yielded regions surrounding the body.
Simulation results for (a) normalised drag forces
$C_s$
, (b) ratio of the pressure drag to shear-induced drag, (c)
$C^*_d$
and (d) normalised torque as a function of Rossby number for various
$\textit{Bi}$
. The continuous and dashed curves correspond to simulations on the sphere and the cylinder, respectively.
The ratio of pressure to frictional viscous drag is shown in figure 14(b). The pressure contribution to the drag coefficient increases as
$\textit{Bi}$
increases. In addition, this ratio exhibits a non-monotonic variation with
$\hat {\varOmega }$
across the entire range of
$\textit{Bi}$
. For instance,
$\hat {\varOmega }=10$
produces the largest
$C_s^P/C_s^V$
ratio for
$\textit{Bi}\le 100$
, whereas for larger
$\textit{Bi}$
, the maximum shifts to the higher rotational rate
$\hat {\varOmega }=100$
. This behaviour results from the evolution of the yielded region generated by the combined translational and swirl motions around the object (figure 15). As the rotation increases, the azimuthal flow expands the yielded envelope surrounding the body, which modifies the pressure distribution between the front and rear stagnation regions and alters the relative contributions of viscous and pressure drag.
Flow field calculated via numerical simulations for
$\hat {\varOmega } =(a)\, 0.1$
, (b)
$1$
and (c)
$10$
at three representative Bingham numbers:
$\textit{Bi} = (i)\,1$
, (ii)
$10$
and (iii)
$100$
. The object sediments from right to left. For each case, the panels represent: (1) normalised azimuthal velocity
$u_{\phi }/U$
; (2) normalised shear rate
$\dot {\gamma }/U$
; (3) normalised in-plane velocity magnitude
$(\boldsymbol{|u|}/U) \text{sgn}(u_x)$
, where the sign indicates direction relative to sedimentation; and (4) normalised pressure
$P/KU^n$
.
Figure 15. Long description
Panel A: Three heat maps showing different characteristics of the flow field for a Bingham number. Each heat map is divided into four quadrants labeled 1 to 4. The color scale ranges from blue to red, indicating varying values. Panel i shows normalized azimuthal velocity, Panel ii shows normalized shear rate, Panel iii shows normalized in-plane velocity magnitude, and Panel iv shows normalized pressure. Panel B: Three heat maps showing different characteristics of the flow field for a different Bingham number. Each heat map is divided into four quadrants labeled 1 to 4. The color scale ranges from blue to red, indicating varying values. Panel i shows normalized azimuthal velocity, Panel ii shows normalized shear rate, Panel iii shows normalized in-plane velocity magnitude, and Panel iv shows normalized pressure. Panel C: Three heat maps showing different characteristics of the flow field for another Bingham number. Each heat map is divided into four quadrants labeled 1 to 4. The color scale ranges from blue to red, indicating varying values. Panel i shows normalized azimuthal velocity, Panel ii shows normalized shear rate, Panel iii shows normalized in-plane velocity magnitude, and Panel iv shows normalized pressure.
Figure 14(c) illustrates the variation of the dimensionless drag coefficient
$C_d^*$
as a function of the Bingham number
$\textit{Bi}$
. For both spherical and cylindrical geometries, two distinct regimes emerge depending on the intensity of the imposed rotation. At low rotation rates,
$C_d^*$
decreases as
$\textit{Bi}$
increases and levels off at
$\textit{Bi}\gt 1$
. In this regime, the yield surface is governed primarily by sedimentation-driven stresses rather than rotational effects, as illustrated in the flow fields of figure 15(a), and the overall shape of the yield region remains similar. In contrast, for high rotation rates (
$\hat {\varOmega } \geqslant 10$
, corresponding to rotationally dominated flows),
$C_d^*$
exhibits a non-monotonic trend: it initially decreases to a local minimum before rising again at higher
$\textit{Bi}$
. This behaviour reflects a fundamental transition in the flow topology, as can be seen by comparing figures 15(b) and 15(c), where rapid rotation initially generates an expansive yielded region that envelops the body and reduces the overall drag.
However, as
$\textit{Bi}$
increases, the rotational stress becomes insufficient to overcome the fluid’s yield stress. Near the local minimum of
$C_d^*$
, the yield surface thickness reduces at the front and rear of the body (for example,
$\textit{Bi} \approx 3.16$
for
$\hat {\varOmega }=10$
) and eventually becomes pinned, progressively contracting towards the solid boundary. Interestingly, in the low-
$\textit{Bi}$
range, a nearly uniform power-law dependence of
$C^*_d$
on
$\textit{Bi}$
is observed across all tested
$\hat {\varOmega }$
curves, where
$C^*_d\sim \textit{Bi}^{-1}$
. This is consistent with the observation that as
$\textit{Bi}$
approaches zero, the drag force should be inversely proportional to
$\textit{Bi}$
. At this limit, different rotation rates primarily modulate the spatial extent of the already enlarged yield domain. By altering this effective confinement, the rotation modifies the resulting hydrodynamic resistance, while its relation to
$\textit{Bi}$
remains almost unchanged. The physical reason behind the observed local minimum in
$C_d^*$
is discussed in Appendix A.
Heat map showing the prediction of (a) non-dimensional drag
$C_s$
, (b) torque
$T_s$
and (c) the ratio of pressure drag to viscous drag
$C_s^P/C_s^V$
as a function of
$\hat {\varOmega }$
number and shear thinning index
$n$
for the cylinder with
$L/D=2.8$
.
6.3. Shear-thinning effects
The drag on a sedimenting and rotating object is expected to depend not only on the yield stress but also on the fluid’s shear-thinning behaviour. Although the effect of the shear-thinning power-law index on the plastic drag coefficient was not investigated experimentally, we examined it systematically through numerical simulations. For
$n=1$
, the fluid behaves as a Bingham plastic, whereas for
$n\lt 1$
, the apparent viscosity decreases with increasing shear rate according to a power-law relationship. Figure 16 presents contour maps of the dimensionless Stokes drag
$C_s$
, torque
$T_s$
, and the pressure-to-viscous drag ratio
$C_s^P/C_s^V$
for the cylinder as functions of the dimensionless rotation rate
$\hat {\varOmega }$
and the shear-thinning index
$n$
, for four representative values of
$\textit{Bi}$
. The corresponding plots for the sphere are shown in supplementary figure S9.
In principle, simulations suggest two distinct regimes based on the topology of the yielded zone: (i) rotationally induced yielding (high
$\hat {\varOmega }$
) and (ii) translationally induced yielding (small
$\hat {\varOmega }$
). In the rotation-dominated regime, the rapid rotational motion generates intense local shear rates in the immediate vicinity of the object’s surface. Because the shear rate is elevated, the fluid’s apparent viscosity in this rotationally induced yield region becomes exceptionally sensitive to
$n$
. Enhancing the shear-thinning effect (decreasing
$n$
) sharply reduces the local effective viscosity, effectively creating a highly mobilised, lubricated layer that envelops the body. Consequently, as shown in figures 16(a) and 16(b), both the drag
$C_s$
and torque
$T_s$
exhibit their highest sensitivity to
$n$
in this limit. Interestingly, for very strong rotation rates (e.g.
$\hat {\varOmega }\gt 100$
), the complex coupling between the yielded envelope size and the lubricated boundary layer leads to non-monotonic variations in the Stokes’ drag coefficients with respect to
$n$
.
Conversely, when translation dominates (
$\hat \varOmega \leq 1$
), the yield surface and the local shear-rate distributions are primarily dictated by the sedimentation of the object. Here, shear is concentrated primarily at the equator, driven by the front-to-back displacement of the fluid. Because the intense shear is less uniformly distributed around the body than in the rotational case, the overall hydrodynamic resistance shows a weaker, more monotonic dependence on the shear-thinning index. Still, as
$n$
increases towards Newtonian behaviour, the drag increases steadily, but with a lower rate than rotationally dominated flows, reflecting the loss of the shear-thinning effect. Additionally, the contour maps of
$ C_s^P/C_s^V$
(shown in figure 16
c) illustrate that shear-thinning significantly alters the balance of pressure and viscous drag forces acting on the object. The localised decrease in viscosity caused by shear-thinning (i.e. lower
$n$
values) leads to an increased relative contribution of the pressure form drag, especially in areas where rotation is significant.
6.4. Yield limit
An important aspect of our experiments is that, below a certain critical mass (or buoyant force), the object ceases to sediment and effectively becomes trapped within the yield-stress medium. In the following, we explore how introducing rotational motion can alter this threshold and trigger sedimentation. Central to this analysis is establishing an experimental criterion for identifying the onset of sedimentation. In the experiments reported here, the yield limit is determined by measuring the settling velocity as a function of added weight to the object. The object’s weight is incrementally adjusted until no detectable sedimentation occurs within a 20 min observation period (see supplementary figure S10 for object trajectories over time at different added weights). Note that experimental duration cannot be extended beyond 20 min due to overheating of the Helmholtz coil. The transition occurs from no motion (within 20 min) to slight motion with
$U \approx 0.1$
[mm min−1]. Moreover, the mass difference between these two conditions is minimal (
$\leqslant$
0.01 g), which approaches the resolution limit of our analytical balance (1 mg), making precise adjustments below 0.01 g challenging and less reliable. Therefore, the onset of sedimentation is defined as the maximum object mass for which no measurable translation is observed over the 20-minute duration of the experiment (i.e.
$U = 0$
).
Figure 17 presents the measured yield limit as a function of imposed
$\varOmega$
(in this case,
$\hat {\varOmega }/\textit{Bi}^{1/n}$
) for both spherical and cylindrical objects with smooth and roughened surfaces. In defining the yield limit, we use conditions under which
$U = 0$
. To eliminate the dependence of the dimensionless groups
$\textit{Bi}$
and
$\hat {\varOmega }$
on the translational velocity
$U$
, we therefore choose to present the x-axis using the combined dimensionless group
$\hat {\varOmega } / \textit{Bi}^{1/n}$
. This rescaling removes the explicit dependence on
$U$
at the yield limit. Several key observations emerge from these experiments. First, the measured yield limit decreases with decreasing rotation for both the cylinder (figure 17
a) and the sphere (figure 17
b), eventually levelling off for
$ \hat {\varOmega }/\textit{Bi}^{1/n} \lt 10^{-1}$
. Second, across nearly all tested rotational velocities, the yield limit is consistently higher for the cylinder than for the sphere. Third, surface roughness plays a significant role: at low
$\hat \varOmega /\textit{Bi}^{1/n}$
numbers, roughened surfaces exhibit a lower yield limit, whereas at higher
$\hat \varOmega /\textit{Bi}^{1/n}$
values (
$\hat \varOmega /\textit{Bi}^{1/n}\gt 0.1$
), smooth surfaces show lower yield limits.
Dimensionless yield limit as a function of
$\hat {\varOmega }/\textit{Bi}^{1/n}$
for (a) smooth and roughened cylinders and (b) spheres in the yield stress fluid. Here, SC and RC refer to smooth and roughened cylinders. In addition, SS and RS are smooth and roughened spheres. The continuous curve in panel (a) shows the high-
$\textit{Bi}$
infinite-cylinder scaling of Hewitt & Balmforth (Reference Hewitt and Balmforth2018). This curve is included as a qualitative asymptotic reference, not as a finite-cylinder fit.
Figure 17. Long description
Panel A: A line graph shows the dimensionless yield limit as a function of a variable for smooth and roughened cylinders. The x-axis is labeled with a variable and the y-axis is labeled with the dimensionless yield limit. The graph includes data points for smooth cylinders (SC) and roughened cylinders (RC), with a continuous curve representing the high-infinite-cylinder scaling. Panel B: Another line graph depicts the dimensionless yield limit as a function of the same variable for smooth and roughened spheres. The x-axis is labeled with the same variable and the y-axis is labeled with the dimensionless yield limit. The graph includes data points for smooth spheres (SS) and roughened spheres (RS).
Although there is no experimentally obtained data in the published literature on the effect of rotation on the yield limit, there are several reports on yield limit values in the absence of rotation. For the non-rotating limit (
$\hat {\varOmega }/\textit{Bi}^{1/n} =0$
), our measured yield limit values (
$Y_{G,\textit{max}} = 0.07$
for RC,
$Y_{G,\textit{max}} = 0.11$
for SC,
$Y_{G,\textit{max}} = 0.03$
for RS and
$Y_{G,\textit{max}} = 0.035$
for SS) are slightly different from the published results. For example, Jossic & Magnin (Reference Jossic and Magnin2001) reported
$Y_{G,\textit{max}} \approx 0.062$
for a roughened sphere and
$Y_{G,\textit{max}} \approx 0.088$
for a smooth sphere. Consistent with our experiments, their results also indicate that the yield limit is higher for smooth surfaces. Other experimental studies have reported yield limits of approximately 0.048 and 0.061 for spheres in yield stress fluids (Atapattu, Chhabra & Uhlherr Reference Atapattu, Chhabra and Uhlherr1995; Tabuteau et al. Reference Tabuteau, Coussot, de and John2007). Likewise, the asymptotic analysis of Beris et al. (Reference Beris, Tsamopoulos, Armstrong and Brown1985), using a Bingham fluid model, estimated
$Y_{G,\textit{max}} \approx 0.048$
for a sphere. Similar consistency is observed for falling cylinders. For cylinders oriented along their longest axis, the yield limit depends on both the aspect ratio and surface roughness. Jossic & Magnin (Reference Jossic and Magnin2001) measured
$Y_{G,\textit{max}} \approx 0.15$
for a smooth cylinder and
$Y_{G,\textit{max}} \approx 0.075$
for a roughened cylinder with an aspect ratio of 5. Numerical work by Iglesias et al. (Reference Iglesias, Mercier, Chaparian and Frigaard2020) likewise reported comparable yield limit values for cylinders with varying aspect ratios. The observed discrepancy between the yield-limit values obtained in our experiments and those reported in the literature may be due to variations in cylinder aspect ratio, surface roughness and/or the rheological characteristics of the surrounding fluid.
The above-mentioned experiments indicate that rotation has little effect on the yield limit at low imposed rotations. At higher rotation rates, however,
$Y_{G,max }$
increases markedly, indicating that rotation facilitates the onset of sedimentation in yield-stress fluids. This trend can be interpreted through the combined shear and pressure contributions to the drag. Rotation generates strong azimuthal shear near the particle surface and produces a yielded envelope around the body. In the high-swirl limit, the surface traction becomes increasingly aligned with the azimuthal direction, so that only its meridional projection contributes to the axial resistance. The shear contribution to the plastic drag therefore decreases as
\begin{equation} C^{*}_{d,\tau } \sim C_{\tau }\big (1+\hat {\varOmega }^{2}\big )^{-1/2}. \end{equation}
If the drag were purely shear-dominated, this would imply
\begin{equation} Y_G \sim \big (1+\hat {\varOmega }^{2}\big )^{1/2}, \end{equation}
consistent with the high-
$\textit{Bi}$
infinite-cylinder result of Hewitt & Balmforth (Reference Hewitt and Balmforth2018). For finite objects, however, Appendix A shows that confinement of the yielded layer can also produce a pressure contribution,
\begin{equation} C^{*}_{d,p} \sim C_{\!p} \textit{Bi}^{2/(n+1)} \hat {\varOmega }^{-(3n+1)/(n+1)}. \end{equation}
Thus, the relevant resistance is
\begin{equation} C_d^* = C^{*}_{d,\tau } + C^{*}_{d,p}, \quad Y_G \sim \frac {1}{C_d^*}. \end{equation}
The observed increase in
$Y_{G,max }$
therefore reflects a reduction in the axial shear resistance by rotation, while deviations from the shear-only scaling can arise from pressure drag associated with finite geometry and yielded-layer confinement.
To examine this interpretation experimentally, we performed detailed PIV measurements around the cylinder and sphere near the onset of sedimentation. Figure 18 shows the flow field around smooth and roughened cylinders at different values of
$\hat {\varOmega }/\textit{Bi}^{1/n}$
at the onset of sedimentation. At small values of
$\hat {\varOmega }/\textit{Bi}^{1/n}$
, the surrounding material barely deforms and no clear plastic deformation zone is captured within the spatial resolution of the PIV measurements. At
$\hat {\varOmega }/\textit{Bi}^{1/n}=0.07$
, a thin plastic deformation zone becomes visible near the rotating object. As
$\hat {\varOmega }/\textit{Bi}^{1/n}$
increases further, this yielded zone grows and moves farther away from the cylinder surface. Similar behaviour is observed for the sphere, as shown in supplementary figure S11. These observations support the idea that the increase in
$Y_{G,max }$
at high rotation is associated with the formation and expansion of a rotation-induced yielded region around the object.
Two-dimensional averaged velocity profile around the cylinders. (a) Smooth and (b) roughened with various
$\hat {\varOmega }/\textit{Bi}^{1/n}$
values from left to right (0.003, 0.07 and 0.28). (c) Radially averaged normalised velocity magnitude for cases shown in panels (a) and (b). In panel (c), the filled and open symbols denote the results for roughened and smooth surfaces, respectively.
Figure 18. Long description
Panel A: Three contour plots show the two-dimensional averaged velocity profile around cylinders. The left plot (i) represents a smooth cylinder with a value of 0.003, the middle plot (ii) represents a smooth cylinder with a value of 0.07, and the right plot (iii) represents a smooth cylinder with a value of 0.28. The color scale on the right indicates the velocity magnitude. Panel B: Three contour plots show the two-dimensional averaged velocity profile around roughened cylinders. The left plot (i) represents a roughened cylinder with a value of 0.003, the middle plot (ii) represents a roughened cylinder with a value of 0.07, and the right plot (iii) represents a roughened cylinder with a value of 0.28. The color scale on the right indicates the velocity magnitude. Panel C: A scatter plot shows the radially averaged normalized velocity magnitude for the cases shown in panels A and B. The x-axis represents the radial distance (r) in millimeters, and the y-axis represents the normalized velocity magnitude. Filled symbols denote the results for roughened surfaces, while open symbols denote the results for smooth surfaces. Different symbols represent different values (0.003, 0.07, and 0.28).
Another feature of experimental results reported in figure 17 is that the roughened objects have higher yield limits compared with the smooth surfaces at high rotation rates (or high
$\hat \varOmega /\textit{Bi}^{1/n}$
values). This suggests that onset of sedimentation requires less mass for roughened surfaces. These results are presumably correlated with the plastic deformation zone and/or wall-slip behaviour around the rotating objects. Our PIV results of figure 18(b) demonstrate that at high
$\hat \varOmega / \textit{Bi}^{1/n}$
numbers, the wall slip is less pronounced around the roughened cylinder. The larger plastically deformed zone around roughened cylinder, which may be caused by less wall-slip, may create a corridor of smaller resistance and lead to a higher yield limit (or equivalently less mass needed to achieve onset of motion).
We note that the numerical simulations should be interpreted with caution in the context of yield limits. The experiments define the yield limit operationally as complete arrest within experimental resolution,
$U=0$
. In contrast, the simulations impose a finite sedimentation velocity and therefore approach the yield limit only asymptotically as
$\textit{Bi}$
becomes large. This distinction already exists for non-rotating particles, where the computed resistance approaches the critical condition discussed by Beris et al. (Reference Beris, Tsamopoulos, Armstrong and Brown1985), but does not directly represent a static arrest calculation. For rotating bodies, the high-
$\textit{Bi}$
infinite-cylinder solution of Hewitt & Balmforth (Reference Hewitt and Balmforth2018) predicts a limiting force
${F_d}/{\textit{Bi}} \sim {2\pi }{\sqrt {1+\hat {\varOmega }^{2}}}$
, implying
$ Y_G \sim \sqrt {1+\hat {\varOmega }^{2}}$
in the limit of very high
$\textit{Bi}$
. The comparison shown in figure 17(a) indicates that the experimental cylinder data follow the same increasing trend with rotation, although quantitative agreement is not expected because the present cylinders are finite, the experiments may involve wall slip, and pressure drag becomes important at high rotation and high
$\textit{Bi}$
. Appendix A provides an analogous local scaling argument for a rotating sphere, showing that the shear contribution gives the same
$\sqrt {1+\hat {\varOmega }^{2}}$
increase in
$Y_G$
, while a confined yielded layer introduces an additional pressure contribution. The latter explains why the apparent
$Y_G$
inferred from velocity-controlled simulations can vary non-monotonically with
$\textit{Bi}$
. Therefore, the simulations and scaling arguments are used here to interpret trends in the onset of motion, not to claim an exact numerical prediction of the arrested yield limit.
It is worth noting that the comparison between the simulation-derived
$Y_G$
and the scaling in Appendix A should be restricted to the portions of the curves that satisfy the asymptotic validity window. For the present Herschel–Bulkley index,
$n=0.36$
, this window is
\begin{align} \hat {\varOmega }^{0.36}\ll \textit{Bi} \ll \hat {\varOmega }^{1.04}. \end{align}
The highlighted intervals in figure 19 identify where this condition is approximately satisfied for each imposed rotation. Within these intervals, the high-rotation numerical curves follow the trend predicted by the local scaling: after the rotation-assisted maximum,
$Y_G$
decreases with increasing
$\textit{Bi}$
as the yielded layer becomes increasingly confined and the pressure contribution becomes more important. In the pressure-dominated limit, the predicted dependence for
$n=0.36$
is
\begin{align} Y_G\propto \textit{Bi}^{-1.47}\hat {\varOmega }^{1.53}, \end{align}
which is consistent with the downward branches observed for the
$\hat {\varOmega }=10^2$
and
$\hat {\varOmega }=10^3$
curves. The lower-rotation cases and the portions of the curves outside the highlighted windows are not expected to follow this asymptotic trend quantitatively, because the assumptions of rotation-dominated yielding and a thin near-surface yielded layer are no longer fully satisfied. Moreover, because the simulations impose a finite settling velocity, they provide an apparent
$Y_G$
based on the computed drag rather than the unique arrested yield limit
$Y_{G,max }$
measured experimentally.
Gravitational yield number as a function of Bingham number for various imposed rotations (
$\hat {\varOmega }$
). Continuous curves correspond to the sphere and dashed curves correspond to the cylindrical cases. The highlighted segments mark the high-swirl thin-layer regime. At fixed
$\hat {\varOmega }$
, increasing
$\textit{Bi}$
thins the rotationally yielded layer and strengthens the pressure contribution to drag. Because
$Y_G$
is inversely proportional to drag, the numerical curves decrease after the peak in the highlighted ranges. Larger
$\hat {\varOmega }$
shifts the valid window to larger
$\textit{Bi}$
showing clear consistency with the scaling.
7. Conclusion and outlook
We investigated the sedimentation of rotating objects in a yield-stress fluid using a combination of experiments and numerical simulations. The results of this study can be summarised as follows.
The drag coefficient was found to be highly sensitive to rotation rate. Our results indicate that increasing rotation rate or
$\hat {\varOmega }$
lead to lower drag for both cylinders and spheres, with roughened surfaces consistently exhibiting greater drag than smooth ones. Numerical simulations reproduced the general trends in drag coefficient with respect to
$\hat {\varOmega }$
, but systematically underpredicted the measured drag, likely due to wall slip and nonlinear flow features such as the formation of negative wakes behind the falling object. Furtheremore, increasing the Bingham number reduced the drag coefficient, reflecting the growing dominance of viscous forces over yield stress. The scaling laws are identified based on simulations over a wide range of
$\textit{Bi}$
and
$\hat {\varOmega }$
.
We further examined the yielding criterion and its dependence on rotation. The results reveal that the yield limit decreases significantly with decreasing rotation rate
$\hat {\varOmega }$
, but levels off for
$\hat \varOmega /\textit{Bi}^{1/n}\lt 0.1$
. Surface roughness also plays a central role: at low
$\hat {\varOmega }$
, smooth surfaces exhibit a higher yield limit, whereas at high
$\hat {\varOmega }$
, roughened surfaces show a higher yield limit. The increase in the measured yield limit with rotation is consistent with high-
$\textit{Bi}$
theory for a rotating/translating cylinder and with the lubrication scaling in Appendix A: rotation reduces the axial projection of the shear resistance, while finite-body pressure effects and wall slip can modify the quantitative threshold. Finally, from a computational perspective, future improvements to the predictive model could include adopting an elastoviscoplastic framework and incorporating surface slip conditions.
Supplementary material
Supplementary material is available at https://doi.org/10.1017/jfm.2026.11782.
Funding
We gratefully acknowledge partial support for this work through National Science Foundation award CBET 2512810.
Declaration of interests
The authors report that there is no conflict of interest.
Appendix A. Lubrication scaling for high-swirl drag and gravitational yield number
The simulations show that, at large
$\hat {\varOmega }$
, the relative contributions of pressure and shear stresses to the drag are strongly affected by the structure of the yielded region surrounding the particle. We next develop a compact lubrication estimate to rationalise this behaviour and to connect the high-swirl drag scaling to the gravitational yield number. The analysis is intended as a local scaling argument rather than a closed solution of the full free-boundary viscoplastic problem.
Consider a sphere of radius
$R$
translating with velocity
$U\boldsymbol{e}_z$
and rotating about the same axis with angular velocity
$\varOmega \boldsymbol{e}_z$
. In spherical coordinates
$(r,\theta ,\phi )$
, where
$\theta$
is measured from the settling direction, the velocity imposed on the particle surface is
\begin{equation} u_r^s=U\cos \theta , \quad u_\theta ^s=-U\sin \theta , \quad u_\phi ^s=\varOmega R\sin \theta . \end{equation}
The tangential surface velocity is therefore
\begin{equation} \boldsymbol{u}_t^s = -U\sin \theta \,\boldsymbol{e}_\theta + \varOmega R\sin \theta \,\boldsymbol{e}_\phi , \quad |\boldsymbol{u}_t^s| = U\sin \theta (1+\hat {\varOmega }^2)^{1/2}, \end{equation}
where
$\hat {\varOmega }=R\varOmega /U$
. In the high-swirl regime,
$\hat {\varOmega }\gg 1$
, the local shear is dominated by the azimuthal motion, while the meridional component associated with sedimentation gives a smaller projection onto the settling direction.
At large
$\textit{Bi}$
, the rotationally yielded region may be localised near the particle surface. We denote its outer boundary by
\begin{equation} r_y(\theta )=R+h(\theta ), \quad h(\theta )\ll R, \end{equation}
and introduce the local normal coordinate
$y=r-R$
, with
$0\lt y\lt h(\theta )$
. In this thin layer, the dominant rotational shear rate is estimated as
\begin{equation} \dot {\gamma }_\varOmega \sim \frac {\varOmega R\sin \theta }{h}. \end{equation}
For a Herschel–Bulkley material, the corresponding viscous stress scale is
\begin{equation} \tau _\varOmega ^v \sim k \left ( \frac {\varOmega R\sin \theta }{h} \right )^n \!. \end{equation}
The yielded layer terminates where the stress magnitude approaches the yield stress. For a rotating sphere, the azimuthal shear stress decreases away from the surface because of curvature. Thus, the excess wall stress needed to sustain a yielded layer of thickness
$h$
scales as
\begin{equation} \tau _\varOmega ^v \sim C_\kappa \sigma _y\frac {h}{R}, \end{equation}
where
$C_\kappa =O(1)$
is a geometric curvature constant (e.g. for purely azimuthal Stokes flow around a sphere,
$C_\kappa =3$
). Combining (A5) and (A6) gives
\begin{equation} \frac {h(\theta )}{R} \sim \left [ \frac {\sin ^n\theta }{C_\kappa \textit{Bi}_\varOmega } \right ]^{1/(n+1)}, \end{equation}
where the rotation-based Bingham number is
$\textit{Bi}_\varOmega =\textit{Bi}/\hat {\varOmega }^n$
. This thin-layer approximation requires
$\textit{Bi}_\varOmega \gg 1$
, or equivalently,
$\hat {\varOmega }^n\ll \textit{Bi}$
.
A.1. Shear contribution
Let
$\boldsymbol{t}_\tau =\boldsymbol{\tau }\boldsymbol{\cdot }\boldsymbol{n}$
denote the deviatoric or shear traction acting on the particle surface. In the large-
$\textit{Bi}$
thin-layer limit, the magnitude of the plastic part of this traction is approximately
$\sigma _y$
, with a smaller Herschel–Bulkley correction. Using (A6),
\begin{equation} |\boldsymbol{t}_\tau | \sim \sigma _y + k \left ( \frac {\varOmega R\sin \theta }{h} \right )^n = \sigma _y \left [ 1+O\left (\frac {h}{R}\right ) \right ]\!. \end{equation}
The tangential traction opposes the local surface motion, so that
\begin{equation} \boldsymbol{t}_\tau \sim - |\boldsymbol{t}_\tau | \frac {\boldsymbol{u}_t^s}{|\boldsymbol{u}_t^s|}. \end{equation}
Only the meridional component of
$\boldsymbol{t}_\tau$
contributes to the axial drag, because
$\boldsymbol{e}_\phi \boldsymbol{\cdot }\boldsymbol{e}_z=0$
and
$\boldsymbol{e}_\theta \boldsymbol{\cdot }\boldsymbol{e}_z=-\sin \theta$
. Hence,
\begin{equation} -\boldsymbol{t}_\tau \boldsymbol{\cdot }\boldsymbol{e}_z \sim |\boldsymbol{t}_\tau | \frac {\sin \theta }{(1+\hat {\varOmega }^2)^{1/2}}. \end{equation}
The shear contribution to the drag is therefore
\begin{equation} F_{d,\tau } = -\int _S \boldsymbol{t}_\tau \boldsymbol{\cdot }\boldsymbol{e}_z\,\text{d}S \sim 2\pi R^2\sigma _y(1+\hat {\varOmega }^2)^{-1/2} \int _0^\pi \sin ^2\theta \,\text{d}\theta . \end{equation}
Using
$\int _0^\pi \sin ^2\theta \,\text{d}\theta =\pi /2$
, we obtain
\begin{equation} F_{d,\tau } \sim \pi ^2R^2\sigma _y(1+\hat {\varOmega }^2)^{-1/2}. \end{equation}
For a sphere,
$A=\pi R^2$
, so the shear contribution to the plastic drag coefficient is
\begin{equation} C_{d,\tau }^* = \frac {F_{d,\tau }}{\pi R^2\sigma _y} \sim \pi (1+\hat {\varOmega }^2)^{-1/2} \sim \pi \hat {\varOmega }^{-1}, \quad \hat {\varOmega }\gg 1 . \end{equation}
The results show that the leading shear drag decreases inversely with
$\hat {\varOmega }$
. This reduction is a geometric projection effect: at high rotation rates, the surface traction is directed mainly azimuthally and only a fraction proportional to
$U/(\varOmega R)=\hat {\varOmega }^{-1}$
contributes to the settling drag.
A.2. Pressure contribution
We next estimate the pressure contribution under the assumption that the unyielded material behaves locally as a stationary plug or
$ u_\theta ^y\simeq 0,\,\,\, u_r^y\simeq 0 \,\,\, \text{at} \,\,\, y=h(\theta ),$
where the superscript
$y$
denotes the velocity at the yield surface. Let
$u_\theta (y,\theta )$
be the meridional velocity in the thin layer and let
$p=p(\theta )$
be the leading-order pressure, which is approximately uniform across the film thickness. The meridional momentum balance is
\begin{equation} \frac {\partial \tau _{\theta y}}{\partial y} = \frac {1}{R}\frac {\text{d}p}{\text{d}\theta }. \end{equation}
Because the azimuthal shear dominates the local deformation, the meridional perturbation may be represented using an effective transverse viscosity,
\begin{equation} \eta _\perp \sim k\dot {\gamma }_\varOmega ^{\,n-1} + \frac {\sigma _y}{\dot {\gamma }_\varOmega }. \end{equation}
In the large-
$\textit{Bi}_\varOmega$
thin-layer limit, the plastic contribution dominates, giving
\begin{equation} \eta _\perp \sim \frac {\sigma _y h}{\varOmega R\sin \theta }. \end{equation}
With
$u_\theta (0,\theta )=u_\theta ^s=-U\sin \theta$
and
$u_\theta (h,\theta )=0$
, the depth-integrated meridional flux is
\begin{equation} q_\theta = \int _0^h u_\theta \,\text{d}y = \frac {h}{2}u_\theta ^s - \frac {h^3}{12\eta _\perp R} \frac {\text{d}p}{\text{d}\theta }. \end{equation}
Mass conservation in the thin layer gives
\begin{equation} \frac {1}{R\sin \theta } \frac{\text{d}}{\text{d}\theta } \left ( \sin \theta \,q_\theta \right ) = u_r^s-u_r^y. \end{equation}
Using
$u_r^s=U\cos \theta$
and
$u_r^y=0$
, substitution of (A17) yields the Reynolds-type equation
\begin{equation} \frac {1}{R^2\sin \theta } \frac{\text{d}}{\text{d}\theta } \left [ \sin \theta \frac {h^3}{12\eta _\perp } \frac {\text{d}p}{\text{d}\theta } \right ] = -\frac {U}{2R\sin \theta } \frac{\text{d}}{\text{d}\theta } \left [ h\sin ^2\theta \right ] - U\cos \theta . \end{equation}
The first term on the right-hand side represents tangential Couette pumping, whereas the second term represents normal squeezing due to sedimentation relative to the stationary plug. In the thin-layer limit, the squeeze term provides the dominant pressure scale.
The pressure contribution to the axial drag is
\begin{equation} F_{d,p} = 2\pi R^2 \int _0^\pi p(\theta )\cos \theta \sin \theta \,\text{d}\theta , \end{equation}
where an arbitrary constant added to
$p$
does not contribute to the drag. A scaling estimate follows from (A19):
\begin{equation} p \sim \eta _\perp U\frac {R^2}{h^3}, \quad F_{d,p} \sim \eta _\perp U\frac {R^4}{h^3}. \end{equation}
Using (A16) gives
\begin{equation} F_{d,p} \sim \sigma _y R^2 \left ( \frac {U}{\varOmega R} \right ) \left ( \frac {R}{h} \right )^2 . \end{equation}
Explicitly written for a sphere,
\begin{equation} C_{d,p}^* = \frac {F_{d,p}}{\pi R^2\sigma _y} \sim C_{\!p}\hat {\varOmega }^{-1} \left ( \frac {R}{h} \right )^2, \end{equation}
where
$C_{\!p}=O(1)$
depends on the detailed pressure distribution. Substituting the film-thickness estimate in (A7) gives, away from the polar regions,
\begin{equation} C_{d,p}^* \sim C_{\!p} \textit{Bi}^{2/(n+1)} \hat {\varOmega }^{-(3n+1)/(n+1)} . \end{equation}
In the Bingham limit,
$n=1$
, this reduces to
\begin{equation} C_{d,p}^* \sim C_{\!p}\frac {\textit{Bi}}{\hat {\varOmega }^2}. \end{equation}
A.3. Total drag and gravitational yield number
Combining the shear and pressure estimates gives the lubrication approximation of
\begin{equation} C_d^* = C_{d,\tau }^* + C_{d,p}^* \sim C_\tau (1+\hat {\varOmega }^2)^{-1/2} + C_{\!p}\hat {\varOmega }^{-1} \left ( C_\kappa \textit{Bi}_\varOmega \right )^{2/(n+1)} , \end{equation}
where
$C_\tau$
and
$C_{\!p}$
are order-one constants. The first term is the direct shear contribution and decreases with
$\hat {\varOmega }$
because of the projection of the surface traction onto the settling direction. The second term is the pressure contribution generated by lubrication flow through a thin yielded layer bounded by a stationary outer plug.
The estimate in (A26) should be interpreted as the large-
$\textit{Bi}$
confined-layer limit. At fixed large
$\hat {\varOmega }$
, the pressure term increases with
$\textit{Bi}_\varOmega$
because the yielded layer becomes thinner and the fluid displaced by sedimentation must pass through a more confined region. Thus, in the asymptotic stationary-plug limit, increasing
$\textit{Bi}$
is expected to increase the pressure contribution to
$C_d^*$
. However, the fully resolved simulations of figure 14(c) show that, for high-
$\hat {\varOmega }$
cases,
$C_d^*$
may first decrease and reach a local minimum at intermediate
$\textit{Bi}$
before increasing again at larger
$\textit{Bi}$
. This behaviour indicates a transition in the shape of the yielded region as discussed in the paper. At intermediate
$\textit{Bi}$
, rotation can expand the yielded region around the particle and partially relieve the pressure build-up, thereby reducing the total drag. This effect may be represented phenomenologically by replacing
$C_{\!p}$
with
$C_{p0}\, \mathcal{M}(\textit{Bi},\hat {\varOmega })$
, where
$\mathcal{M}\lt 1$
represents the increased mobility of the connected yielded region.
At larger
$\textit{Bi}$
, the yielded region contracts towards the particle surface, the locally confined-layer assumption becomes more appropriate and the pressure contribution increases. The observed local minimum in
$C_d^*$
therefore marks the transition from a rotation-assisted, pressure-relieved regime to a high-
$\textit{Bi}$
confined lubrication regime.
In either of the cases, the relative importance of the two terms is obtained from (A13) and (A24):
\begin{equation} \frac {C_{d,p}^*}{C_{d,\tau }^*} \sim C_{p\tau } \left ( C_\kappa \textit{Bi}_\varOmega \right )^{2/(n+1)} . \end{equation}
This estimate explains why pressure drag can become large when the yield surface is close to the particle. We can now estimate the gravitational yield number for a sedimenting rotating sphere,
\begin{equation} Y_G(\textit{Bi},\hat {\varOmega }) \sim \frac {2/3} { C_\tau (1+\hat {\varOmega }^2)^{-1/2} + C_{\!p}\hat {\varOmega }^{-1} \left ( C_\kappa \textit{Bi}_\varOmega \right )^{2/(n+1)} }. \end{equation}
This expression provides a compact scaling connection between rotation, the Bingham number and the apparent gravitational yield number. Because
$Y_G=(2/3)/C_d^*$
for a sphere, any local minimum in
$C_d^*$
corresponds to a local maximum in
$Y_G$
. Therefore, the local minimum in
$C_d^*$
observed in the fully resolved high-
$\hat {\varOmega }$
simulations implies that the apparent yield number can increase at intermediate
$\textit{Bi}$
and then decrease again as
$\textit{Bi}$
becomes large. In the shear-dominated high-swirl limit, the pressure term is negligible and
\begin{equation} Y_G \sim \frac {2}{3 C_\tau } (1+\hat {\varOmega }^2)^{1/2} \sim \frac {2}{3 C_\tau }\hat {\varOmega }, \quad \hat {\varOmega }\gg 1 . \end{equation}
In contrast, if the pressure term dominates, then for a Bingham material,
\begin{equation} Y_G \sim \frac {2}{3\,C_{\!p}} \frac {\hat {\varOmega }^2}{\textit{Bi}}. \end{equation}
Thus, the variation of
$Y_G$
with
$\textit{Bi}$
at fixed high
$\hat {\varOmega }$
need not be monotonic. At intermediate
$\textit{Bi}$
, rotation-assisted yielding can reduce
$C_d^*$
and increase
$Y_G$
, whereas at sufficiently large
$\textit{Bi}$
, the yielded layer becomes thin and confined, pressure drag increases, and
$Y_G$
decreases. This interpretation is consistent with the local minimum in
$C_d^*$
observed in the fully resolved simulations (figure 14
c). Thus, the apparent gravitational yield number increases with rotation because rotation reduces the drag required to maintain a given settling velocity. However, the rate of this increase depends on whether the resistance is governed primarily by shear projection or by pressure build-up in the confined yielded layer.
A.4. Range of validity
These scaling relations are derived based on four assumptions.
-
(i) The flow must be rotation-dominated:
(A31)
\begin{equation} \hat {\varOmega }=\frac {R\varOmega }{U}\gg 1 . \end{equation}
-
(ii) The yielded layer must be thin and close to the interface:
(A32)which requires
\begin{equation} \frac {h}{R}\ll 1, \quad \frac {h}{R} \sim \left [ \frac {\sin ^n\theta } {C_\kappa \textit{Bi}_\varOmega } \right ]^{1/(n+1)} , \end{equation}
$\textit{Bi}_\varOmega \gg 1$
away from the polar regions or, equivalently,
$\hat {\varOmega }^n\ll \textit{Bi}$
.
-
(iii) The rotation-controlled estimate of
$h$
assumes that the pressure stress does not reorganise the yield surface. A useful consistency requirement is(A33)Combining this with
\begin{equation} \hat {\varOmega } \gg \left ( C_\kappa \textit{Bi}_\varOmega \right )^{2/(n+1)} . \end{equation}
$\textit{Bi}_\varOmega \gg 1$
gives(A34)which for Bingham flow translates to
\begin{equation} {\textit{Bi}}^{\frac {2}{3n+1}}\ll \hat {\varOmega } \ll \textit{Bi}^{\frac {1}{n}} , \end{equation}
(A35)Within this window, the yielded layer is thin, its thickness is primarily set by rotational shear, and the stationary-plug pressure estimate is asymptotically consistent.
\begin{equation} \sqrt {\textit{Bi}}\ll \hat {\varOmega } \ll \textit{Bi} . \end{equation}
-
(iv) The analysis is not valid near
$\theta =0$
and
$\theta =\pi$
, where
$\sin \theta \rightarrow 0$
and the rotational surface velocity vanishes. Those polar regions require a separate local description.
For scaling analysis, inside these small polar regions, the local flow can be approximated as squeeze flow between the particle surface and the locally stationary outer plug. Introducing the local polar coordinate
$\rho =R\theta$
near
$\theta =0$
and
$\rho =R(\pi -\theta )$
near
$\theta =\pi$
, the normal velocity remains
$O(U)$
, while the rotational velocity scales as
$\varOmega \rho$
. If the polar cap has radius
$\ell _{\!p}$
and representative gap thickness
$h_{\!p}$
, the local Reynolds-type balance gives
\begin{equation} F_{d,p}^{\textit{polar}} \sim \eta _{\textit{eff},p} U\frac {\ell _{\!p}^4}{h_{\!p}^3}. \end{equation}
The effective viscosity is evaluated using the squeeze-induced shear rate,
\begin{equation} \dot {\gamma }_{\!p} \sim \frac {U\ell _{\!p}}{h_{\!p}^2}, \quad \eta _{\textit{eff},p} \sim k\dot {\gamma }_{\!p}^{\,n-1} + \frac {\sigma _y}{\dot {\gamma }_{\!p}}. \end{equation}
When the yield-stress contribution dominates, this reduces to
\begin{equation} F_{d,p}^{\textit{polar}} \sim \sigma _y\frac {\ell _{\!p}^3}{h_{\!p}}, \quad C_{d,p}^{*,\mathrm{polar}} \sim C_{p0} \left (\frac {\ell _{\!p}}{R}\right )^2 \left (\frac {\ell _{\!p}}{h_{\!p}}\right )\!, \end{equation}
where
$C_{p0}=O(1)$
. The polar cap size follows by matching the squeeze-induced shear to the outer rotational shear,
\begin{equation} \frac {\varOmega \ell _{\!p}}{h_{\!p}} \sim \frac {U\ell _{\!p}}{h_{\!p}^2} \quad \Rightarrow \quad \frac {h_{\!p}}{R}\sim \hat {\varOmega }^{-1}. \end{equation}
Using the outer estimate
$h/R\sim [\theta ^n/(C_\kappa \textit{Bi}_\varOmega )]^{1/(n+1)}$
gives
\begin{equation} \theta _c=\frac {\ell _{\!p}}{R} \sim \left (C_\kappa \textit{Bi}\right )^{1/n} \hat {\varOmega }^{-(2+1/n)}, \end{equation}
and, for a Bingham material, results in
$ \theta _c\sim C_\kappa \textit{Bi}\,\hat {\varOmega }^{-3}$
. This means that when
$\theta _c\ll 1$
, the polar caps mainly regularise the outer solution and modify the order-one pressure coefficient. If
$\theta _c$
is not small, the pressure field and yielded-layer topology must be obtained from the full viscoplastic solution.
Fd
Fg
0.15R
5R
20R
R
uz=U
ur=0
uθ=r(z)Ω
∂uz/∂r=0
z
r
θ
Cs
ϵ
Bi=2
△
50
◻
200
∘
Ω^=10,103
66%
Ω^=1
Bi
Ω^
Bi
σ
γ˙
σ=σy+kγ˙n
σy
k
n
n
Ω^
Ω^eff
Bi≈2.6
Ω^
v
U
Bi≈2.6
Ω^
r=0
Cs
Cd∗
Bi
Ω^
Bi
γ˙
Ω^
Bi=1
Ω^
Bi
Ω
Bi≈6
Ω^=(i)40
z
Cs
Cd∗
Bi
Ω^=(a)0.1
1
10
Bi=(i)1
10
100
uϕ/U
γ˙/U
(|u|/U)sgn(ux)
P/KUn
Cs
Ts
CsP/CsV
Ω^
n
L/D=2.8
Ω^/Bi1/n
Bi
Ω^/Bi1/n
Ω^
Ω^
Bi
YG
Ω^
Bi
Open access
