2.1. The squirmer model

We consider an ellipsoidal microswimmer propelling through a viscoplastic fluid, as illustrated in figure 1(a). The swimmer has the aspect ratio $a_r = b_x/b_z$ , with semi-major axis $b_x$ and semi-minor axis $b_z$ , as shown in figure 1(b). We denote half of the focal length by $c = \sqrt {b_x^2 - b_z^2}$ , which yields the eccentricity $e=c/b_x$ ( $0 \leqslant e \lt 1$ , with $e\,{=}\,0$ describing a sphere). The prolate spheroidal coordinate system ( $\zeta, \tau, \varphi$ ) is related to the body-fixed Cartesian coordinates $(x^{\prime},y^{\prime},z^{\prime})$ via (Dassios, Hadjinicolaou & Payatakes Reference Dassios, Hadjinicolaou and Payatakes1994)

(2.1) \begin{equation} \zeta = \frac {1}{{2c}}\left ( {\sqrt {{{(x^{\prime} + c)}^2}+{y^{\prime 2}} + {z^{\prime 2}}} - \sqrt {{{(x^{\prime} - c)}^2}+{y^{\prime 2}} + {z^{\prime 2}}} } \right), \end{equation}

(2.2) \begin{equation} \tau = \frac {1}{{2c}}\left ( {\sqrt {{{(x^{\prime} + c)}^2}+{y^{\prime 2}} + {z^{\prime 2}}} + \sqrt {{{(x^{\prime} - c)}^2}+{y^{\prime 2}} + {z^{\prime 2}}} } \right), \end{equation}

(2.3) \begin{equation} \varphi = \arctan \left ( {\frac {y^{\prime}}{z^{\prime}}} \right). \end{equation}

Figure 1.

(a) Computational set-up of an ellipsoidal microswimmer in a Bingham fluid environment. (b) Sketch of an ellipsoidal microswimmer, where $b_x$ and $b_z$ are the semi-major and semi-minor axes of the ellipsoid. Here, $r_{eq}=(b_xb_z^2)^{1/3}$ refers to the equivalent radius of the ellipsoidal microswimmer.

In the equations, $-1 \leqslant \zeta \leqslant 1$ , $1 \leqslant \tau$ and $0 \leqslant \varphi \lt 2\pi$ . All points with $\tau = \tau _0 \equiv e^{-1}$ (i.e. ${x^{\prime 2}}/b_x^2 + ( {{{y^{\prime}}^2} + {{z^{\prime}}^2}})/b_z^2 = 1$ ) lie on the ellipsoid’s surface.

The microswimmer is modelled as a spheroidal squirmer in this research. In the squirmer model, the propulsion is generated by imposing surface squirming velocities. The surface velocity is assumed to be tangential, steady and axisymmetric. The general form of the surface velocity with an azimuthal component is expressed as (Pak & Lauga Reference Pak and Lauga2014; Pedley, Brumley & Goldstein Reference Pedley, Brumley and Goldstein2016)

(2.4) \begin{equation} \boldsymbol{u}_{sq}=-\sum _{n=1}^{\infty } \frac {2 B_n}{n(n+1)} P_n^{\prime }(\boldsymbol{p} \cdot \boldsymbol{n}) \boldsymbol{p} \cdot (\unicode{x1D644}-\boldsymbol{n n})-\sum _{n=1}^{\infty } C_n P_n^{\prime }(\boldsymbol{p} \cdot \boldsymbol{n}) \left | \boldsymbol{p} \cdot (\unicode{x1D644}-\boldsymbol{n n})\right | \boldsymbol{t}, \end{equation}

where ${\boldsymbol{p}}=\boldsymbol{e}_p$ is the unit orientation vector, $\boldsymbol{n}=\boldsymbol{e}_\tau$ is the unit normal vector and $\boldsymbol{t}=\boldsymbol{e}_\varphi$ is the unit azimuthal tangent vector; $\unicode{x1D644}$ is the unit tensor and $P^{\prime}_n$ is the derivative of the $n$ th-order Legendre polynomial. The coefficients $B_n$ and $C_n$ are often called polar modes and azimuthal (or swirling) modes, respectively. Typically, only the first two polar squirming modes ( $B_1$ and $B_2$ ) are considered. The first polar mode $B_1$ determines the swimming speed $U_N$ in an infinite Newtonian fluid in the Stokes flow regime as (Keller & Wu Reference Keller and Wu1977; Theers et al. Reference Theers, Westphal, Gompper and Winkler2016)

(2.5) \begin{equation} U_N=B_1 \tau _0\left(\tau _0-\left(\tau _0^2-1\right) \operatorname {coth}^{-1} \tau _0\right)\kern-2.5pt. \end{equation}

The swimming speed of a spherical swimmer is recovered for the spherical limit ( ${\tau _0} \to \infty$ ) as $U_N=2B_1/3$ . Therefore, $B_1$ is chosen as the characteristic velocity scale. The second polar mode $B_2$ regulates the strength of a force dipole, i.e. the stresslet, in an infinite Newtonian fluid (Ishikawa, Simmonds & Pedley Reference Ishikawa, Simmonds and Pedley2006). Recently, researchers have started to consider the effect of higher-order polar modes (Datt et al. Reference Datt, Zhu, Elfring and Pak2015; De Corato & D’Avino Reference De Corato and D’Avino2017; Pietrzyk et al. Reference Pietrzyk, Nganguia, Datt, Zhu, Elfring and Pak2019) or the azimuthal swirling modes (Pak & Lauga Reference Pak and Lauga2014; Pedley et al. Reference Pedley, Brumley and Goldstein2016; Binagia et al. Reference Binagia, Phoa, Housiadas and Shaqfeh2020; Fortune et al. Reference Fortune, Worley, Sendova-Franks, Franks, Leptos, Lauga and Goldstein2021) on the swimming of spherical squirmers. Nevertheless, the effect of the azimuthal modes on the swimming dynamics of non-spherical squirmers in a complex fluid has not been reported.

We consider the prescribed surface velocity (2.4) with the first two polar modes ( $B_1$ and $B_2$ ) and the second azimuthal mode ( $C_2$ ) as being non-zero, then it is simplified as

(2.6) \begin{align} {{\boldsymbol{u}}_{{{sq}}}} & = - {B_1} ( {{{\boldsymbol{e}}_\zeta }\cdot {{\boldsymbol{e}}_p}} ){{\boldsymbol{e}}_\zeta } - {B_2}\zeta ( {{{\boldsymbol{e}}_\zeta }\cdot {{\boldsymbol{e}}_p}} ){{\boldsymbol{e}}_\zeta } - 3{C_2}\zeta ( {{{\boldsymbol{e}}_\zeta } \cdot {{\boldsymbol{e}}_p}} ){{\boldsymbol{e}}_\varphi } \nonumber\\& =-B_1 \tau _0 (1-\zeta ^2)^{\frac {1}{2}} (\tau _0^2-\zeta ^2)^{-\frac {1}{2}} [(1+\beta \zeta) \boldsymbol{e}_\zeta +3 \chi \zeta \boldsymbol{e}_{\varphi }]. \end{align}

The ratios $\beta = B_2/B_1$ and $\chi =C_2/B_1$ represent the squirmer’s dipolarity and chirality, respectively. Dipolarity $\beta$ helps in categorising the types of microswimmer as pusher ( $\beta\,{\lt}\,0$ , e.g. E. coli), puller ( $\beta \gt 0$ , e.g. Chlamydomonas) and neutral ( $\beta = 0$ , e.g. Volvox) types, while chirality $\chi$ represents the intensity of a microswimmer with rotating flagella and a counter-rotating body (Pedley et al. Reference Pedley, Brumley and Goldstein2016; Binagia et al. Reference Binagia, Phoa, Housiadas and Shaqfeh2020; Fadda, Molina & Yamamoto Reference Fadda, Molina and Yamamoto2020). We refer to the squirmer with $\chi \ne 0$ as a swirling squirmer. Note that the surface velocity for spherical squirmer is recovered in the spherical limit ( $\zeta \rightarrow \cos \theta$ ) (Housiadas et al. Reference Housiadas, Binagia and Shaqfeh2021; Kobayashi, Molina & Yamamoto Reference Kobayashi, Molina and Yamamoto2024) as

(2.7) \begin{equation} {{\boldsymbol{u}}_{sq}} = \left ( {{B_1}\sin \theta + \frac {{{B_2}}}{2}\sin 2\theta } \right){{\boldsymbol{e}}_\theta } + \frac {3}{2}{C_2}\sin 2\theta {{\boldsymbol{e}}_\varphi }, \end{equation}

where $\theta =\arccos (\boldsymbol{p} \cdot \boldsymbol{r} / r)$ represents the polar angle between the position vector $\boldsymbol{r}$ and the swimming direction $\boldsymbol{p}$ , while ${\boldsymbol{e}}_\theta$ is the unit tangent vector along the $\theta$ direction. In the present study, we investigate three types of swimmers: neutral swimmers ( $\beta =\chi =0$ ), pullers or pushers ( $\beta \ne 0,\chi =0$ ) and neutral swirling swimmers ( $\beta =0,\chi \ne 0$ ).

2.2. Basic equations

The locomotion of a squirmer through a viscoplastic fluid is simulated by the direct forcing/fictitious domain method (Yu & Shao Reference Yu and Shao2007; Yu & Wachs Reference Yu and Wachs2007). Let $\Omega$ and $P$ represent the entire domain and the domain occupied by the squirmer, respectively ( $P\,{\subset}\,\Omega$ ); $P$ is bounded by the squirmer surface $S$ . The governing equations are written as

(2.8) \begin{equation} \rho _{\!f} \frac {\mathrm{\ d} \boldsymbol{u}}{\mathrm{\ d} t}= \boldsymbol{\nabla } \cdot \left (-p \unicode{x1D644}+\boldsymbol{\tau }\right) + \boldsymbol{\lambda }\;\; {\textrm {in}} \;\;{\Omega }, \end{equation}

(2.9) \begin{equation} \boldsymbol{u}=(\boldsymbol{U}+\boldsymbol{\omega } \times \boldsymbol{r})+\boldsymbol{u}_{sq}\;\; {\textrm {in}} \;\;{S}, \end{equation}

(2.10) \begin{equation} \rho _s V_p \frac {\mathrm{\ d} \boldsymbol{U}}{\mathrm{\ d} t}=\int _S \boldsymbol{n} \cdot \left (-p \unicode{x1D644}+\boldsymbol{\tau }\right) \mathrm{d} s, \end{equation}

(2.11) \begin{equation} \frac {\mathrm{\ d}(\unicode{x1D645} \cdot {\boldsymbol{\omega }})}{\mathrm{d} t}=\int _S \boldsymbol{r} \times [\boldsymbol{n} \cdot (-p \unicode{x1D644}+\boldsymbol{\tau })] \mathrm{d} s. \end{equation}

In these equations, $\boldsymbol{u}$ , $p$ and $\rho _{\!f}$ represent the fluid velocity, pressure and density, respectively, $\boldsymbol{\tau }$ is the stress tensor and $\boldsymbol{\lambda }$ is a pseudo-body force (Lagrange multiplier) defined in the particle domain. The density, volume, moment of inertia tensor, translational velocity and angular velocity of the squirmer are denoted by $\rho _s$ , $V_p$ , $\unicode{x1D645}$ , $\boldsymbol{U}$ and $\boldsymbol{\omega }$ , respectively, and $\boldsymbol{r}$ is the position vector with respect to the particle centre of mass. A solenoidal velocity field $\boldsymbol{u}_a$ inside the squirmer is imposed (see Appendix A for the details of derivation). Therefore, the velocity field inside the squirmer is given by $ \boldsymbol{u}=(\boldsymbol{U}+\boldsymbol{\omega } \times \boldsymbol{r})+\boldsymbol{u}_a$ .

As a viscoplastic fluid, we assume a Bingham fluid model as one of the simplest viscoplastic fluids, for which the constitutive equation reads

(2.12) \begin{align}\boldsymbol{\tau} & =\left (2 \mu +\frac {\tau _y}{|\unicode{x1D640}|}\right) \unicode{x1D640}\quad \text{ for }|\boldsymbol{\tau }|\gt \tau _y, \nonumber\\ \unicode{x1D640} & =\boldsymbol{0}\quad \text{ for } |\boldsymbol{\tau }| \leqslant \tau _y. \end{align}

Here, $\tau _y$ is the yield stress, $\mu$ is the viscosity constant and $\unicode{x1D640}$ is the rate of strain tensor. The second invariant of a tensor is defined as $|\boldsymbol{\delta }|= \sqrt {(\boldsymbol{\delta }: \boldsymbol{\delta })/2}$ for any tensor $\boldsymbol{\delta } \in \mathbb{R}^{3 \times 3}$ . The Bingham model was employed to investigate the viscoplastic effect on the self-propulsion of a spherical squirmer (Eastham et al. Reference Eastham, Mohammadigoushki and Shoele2022), but it did not consider non-spherical body shapes.

We employ the augmented Lagrangian method proposed by Dean & Glowinski (Reference Dean and Glowinski2002) for the solution of the constitutive equation (2.12), which are rewritten as

(2.13) \begin{equation} \boldsymbol{\tau }=2 \mu \unicode{x1D640}+\tau _y \unicode{x1D63D}, \end{equation}

where $\unicode{x1D63D}$ is a non-dimensional tensor-valued function which satisfies $|{\unicode{x1D63D}}| \leqslant 1$ , ${\unicode{x1D63D}}:{\unicode{x1D640}} = |{\unicode{x1D640}}|$ , as well as the following relation (Dean & Glowinski Reference Dean and Glowinski2002)

(2.14) \begin{equation} \unicode{x1D63D}=P_{\Lambda }\left (\unicode{x1D63D}+\sigma \tau _y \unicode{x1D640}\right) \quad \forall \ \sigma \gt 0, \end{equation}

where $P_{\Lambda }(\unicode{x1D666})$ is a symmetry-preserving orthogonal projection operator defined by

(2.15) \begin{equation} P_{\Lambda }(\unicode{x1D666})= \begin{cases}\unicode{x1D666} & \text{ if }|\unicode{x1D666}| \leqslant 1, \\ \unicode{x1D666}{/}|\unicode{x1D666}| & \text{ if }|\unicode{x1D666}|\gt 1.\end{cases} \end{equation}

Here, $\unicode{x1D666}$ represents a non-dimensional tensor. Equation (2.14) implies that $\unicode{x1D640} = \boldsymbol{0}$ when $| \unicode{x1D63D}+\sigma \tau _y \unicode{x1D640} | \leqslant 1$ and $\unicode{x1D63D} = \unicode{x1D640}/|\unicode{x1D640}|$ when $| \unicode{x1D63D}+\sigma \tau _y \unicode{x1D640} | \gt 1$ , and thus (2.13) is equivalent to (2.12).

The following characteristic scales are used for the non-dimensionalisation scheme: $L_c = r_{eq}$ ( $r_{eq}= (b_xb_z^2)^{1/3}$ being the radius of a sphere with the same volume as the ellipsoid, i.e. the equivalent radius) for length, $U_c = B_1$ for velocity, $t_c=L_{c}/U_c$ for time, $\mu U_c/L_c$ for stress and $\mu U_c/L^2_c$ for body force. For convenience, we write the dimensionless quantities in the same form as their dimensional counterparts, unless otherwise specified. By employing the variations $\boldsymbol{v}$ , $\boldsymbol{V}$ and $\boldsymbol{\xi }$ for the fluid velocity and the squirmer translational and angular velocities, respectively, the complete set of dimensionless governing equations in the weak form comprise the following three parts.

(i) Equations of motion

(2.16) \begin{equation} Re \int _{\Omega }\left (\frac {\partial \boldsymbol{u}}{\partial t}+\boldsymbol{u} \cdot \nabla \boldsymbol{u}\right) \cdot \boldsymbol{v} \mathrm{d} \boldsymbol{x}=\int _{\Omega }(-\nabla p+ \nabla ^2 \boldsymbol{u}+Bi \nabla \cdot \unicode{x1D63D}) \cdot \boldsymbol{v}\,\mathrm{d} \boldsymbol{x}+\int _P \boldsymbol{\lambda } \cdot \boldsymbol{v}\ \mathrm{d} \boldsymbol{x}, \end{equation}

(2.17) \begin{equation} Re \left (\rho _r-1\right)\left (V_p\frac {\mathrm{d} \boldsymbol{U}}{\mathrm{d} t} \cdot \boldsymbol{V}+\unicode{x1D645} \cdot \frac {\mathrm{d} \boldsymbol{\omega }}{\mathrm{d} t} \cdot \boldsymbol{\xi }\right)=\int _P \!\!\left (Re {\frac {{\textrm {d}{\boldsymbol{u}_{a}}}}{{\textrm {d} t}}\textrm {d}{\boldsymbol{x}}}-\boldsymbol{\lambda }\right) \cdot (\boldsymbol{V}+\boldsymbol{\xi } \times \boldsymbol{r}) \mathrm{d} \boldsymbol{x}, \end{equation}

(2.18) \begin{equation} \int _P \left [\boldsymbol{u}-(\boldsymbol{U}+\boldsymbol{\omega } \times \boldsymbol{r})-\boldsymbol{u}_a \right] \cdot \boldsymbol{\gamma}\ \mathrm{d} \boldsymbol{x}=\mathbf{0}. \end{equation}

(ii) Continuity equation

(2.19) \begin{equation} \int _{\Omega } q \nabla \cdot \boldsymbol{u}\ \mathrm{d} \boldsymbol{x}=\mathbf{0}. \end{equation}

(iii) Constitutive equation

(2.20) \begin{equation} \unicode{x1D63D}=P_\Lambda \left (\unicode{x1D63D}+\sigma Bi \unicode{x1D640}\right) \quad \forall \ \sigma \gt 0. \end{equation}

In (2.18) and (2.19), $\boldsymbol{\gamma }$ and $q$ represent the corresponding variations. In this study, $\rho _r=\rho _s/\rho _{\!f}$ is the swimmer–fluid density ratio, which is set to be unity. The following dimensionless numbers are introduced:

(2.21) \begin{equation} {\textrm {Reynolds number}}:\ Re = \frac {{{\rho _{\!f}}{B_1}{r_{eq}}}}{{\mu }}, \end{equation}

(2.22) \begin{equation} {\textrm {Bingham number}}:\ Bi = \frac {{{\tau _y}{r_{eq}}}}{{\mu {B_1}}}. \end{equation}

The Bingham number here represents the non-dimensional yield stress.

2.3. Numerical method

Throughout this study, we set $\sigma =Re/(2Bi^2)$ in (2.20) as suggested by Dean & Glowinski (Reference Dean and Glowinski2002). The problem (2.16)–(2.20) is solved in an iterative way with the fractional step time scheme, as detailed in Yu & Wachs (Reference Yu and Wachs2007). The fictitious domain method for dealing with the swimming of a spherical or ellipsoidal squirmer in a Newtonian fluid has been validated in our previous studies (Lin & Gao Reference Lin and Gao2019; Xia et al. Reference Xia, Yu, Lin, Lin and Hu2025a ,Reference Xia, Yu, Zhang, Lin and Ouyang b ). The accuracy for particle sedimentation in Bingham fluids was verified in Yu & Wachs (Reference Yu and Wachs2007). In this study, we validate the swimming speed, power dissipation and swimming efficiency in a Newtonian fluid against the theoretical results derived by Demir et al. (Reference Demir, van Gogh, Palaniappan and Nganguia2024). A convergence test is performed in Appendix B.

As shown in figure 1(a), a cubic periodic computational domain is adopted with $L_x=L_y=L_z=32r_{eq}$ . This configuration allows the squirmer to swim periodically through the domain. The computational domain is sufficiently large so that the effect of domain size completely disappears for a viscoplastic fluid, and it can reproduce an unbounded domain. According to the convergence test in Appendix B, 16 grid points are used across the squirmer’s equivalent radius (i.e. $r_{eq}/\Delta x=16$ ) and the time step is set to be $\Delta t = 10^{-4}r_{eq}/B_1$ . The definitions and values of the non-dimensional parameters are summarised in table 1. In this paper, we consider a low Reynolds number of $Re=0.01$ , which is small enough to neglect inertia. The aspect ratios range up to 8, corresponding to eccentricities up to 0.99. For each simulation, the squirmer was started from rest and allowed to accelerate until it reached a steady swimming speed. Simulation results were obtained after reaching a steady time-averaged swimming speed.

Table 1.

List of the non-dimensional simulation parameters.

3.1. Neutral squirmer

We compute the swimming speed, power dissipation and swimming efficiency relative to their corresponding Newtonian values as a function of Bingham number $Bi$ for both spherical and ellipsoidal neutral squirmers in figure 2. The results for the spherical squirmer ( $a_r=1$ ) agree well with results obtained from the finite element simulation of Eastham et al. (Reference Eastham, Mohammadigoushki and Shoele2022). From figure 2(a,b), we find that the squirmers for all aspect ratios swim slower and require more power than their Newtonian counterpart as the Bingham number increases. While the squirmer experiences a reduction in speed, more elongated squirmers are able to maintain speeds closer to that of their counterpart in Newtonian fluids. In the range of $10 \leqslant Bi \leqslant 100$ , the power dissipation increases approximately exponentially with $Bi$ , following $P/P_N \sim Bi^{0.83}$ from fitting, and is nearly independent of the aspect ratio. Although the power dissipation increases monotonically with the Bingham number, figure 2(c) shows that the variation in the swimming efficiency is non-monotonic as $Bi$ increases. The scaled swimming efficiency, $\eta /\eta _N$ , reaches a maximum value before decreasing to a low value. This behaviour has been reported by Eastham et al. (Reference Eastham, Mohammadigoushki and Shoele2022) for a spherical squirmer. The optimal Bingham number for maximum efficiency is smaller in the case of a long ellipsoid. At $Bi \leqslant 0.1$ , the elongated squirmer swims more efficiently, while the spherical squirmer has higher efficiency for moderate to high values of the Bingham number.

Figure 2.

(a) Swimming speed $U$ , (b) power dissipation $P$ and (c) swimming efficiency $\eta$ , of a neutral squirmer, as a function of the Bingham number $Bi$ . The dashed lines represent the numerical results obtained using the continuous Galerkin finite element method for spherical squirmers, as presented in Eastham & Shoele (Reference Eastham and Shoele2020). Here, the results are scaled by the corresponding values for a squirmer in a Newtonian fluid.

The observations in figure 2 can be made more clearly by plotting the results as a function of the aspect ratio for various $Bi$ . Figure 3(a) suggests that ellipsoidal squirmers swim much faster compared with spherical squirmers, and the influence of body shape on swimming speed becomes significantly more pronounced in environments with higher viscoplasticity. As illustrated in figure 3(b,c), ellipsoidal squirmers expend less energy and are more efficient swimmers compared with spherical squirmers for the same value of the Bingham number.

Figure 3.

(a) Normalised swimming speed, (b) normalised power dissipation and (c) swimming efficiency, of a neutral squirmer, as a function of aspect ratio. Curves are obtained from Demir et al. (Reference Demir, van Gogh, Palaniappan and Nganguia2024) for swimming in Newtonian fluids, while the symbols denote numerical results.

The results in figure 2 have been obtained for ellipsoidal microswimmers with constraints on volume, while some studies fixed the semi-major length $b_x$ when the aspect ratio changes (van Gogh et al. Reference van Gogh, Demir, Palaniappan and Pak2022; Demir et al. Reference Demir, van Gogh, Palaniappan and Nganguia2024). In figure 4, we replot figure 2 as a function of the modified Bingham number $Bi^*$ , which is defined as

(3.3) \begin{equation} Bi^* = \frac {{{\tau _y}{b_x}}}{{\mu {B_1}}}=a^{2/3}_r Bi. \end{equation}

Our results show that the constraint of constant semi-major length does not qualitatively alter the results obtained under the fixed volume assumption. One interesting observation is that the scaled power dissipation has a much weaker dependence on the aspect ratio under the constant semi-major length constraint, as observed in figure 4(b). A similar weak dependence of particle shape on power dissipation has also been reported by Demir et al. (Reference Demir, van Gogh, Palaniappan and Nganguia2024) for a squirmer swimming through a heterogeneous medium with fixed semi-major axis length. Figure 4(c) exhibits a similar non-monotonic variation of the scaled efficiency with respect to $Bi$ , as observed in figure 2(c).

Figure 4.

(a) Swimming speed $U$ , (b) power dissipation $P$ and (c) swimming efficiency $\eta$ , as a function of the modified Bingham number $Bi^*$ , defined by using the semi-major length $b_x$ . The results are normalised by their corresponding values for a neutral squirmer in a Newtonian fluid.

The locomotion behaviour of spherical squirmer closely correlated with the confinement effects induced by the yielded region surrounding the swimming body (Eastham et al. Reference Eastham, Mohammadigoushki and Shoele2022). Figure 5(a–c) shows the yield-surface profiles around a neutral swimmer for various aspect ratios. The yield surface is given by the iso-surface of $\left |\boldsymbol{\tau }\right |=Bi$ . Indeed, increasing the Bingham number for a fixed aspect ratio results in a reduction of the yielded region around the squirmer and the plasticity becomes more dominant, as previously demonstrated by Eastham et al. (Reference Eastham, Mohammadigoushki and Shoele2022). The decreased swimming speed of neutral swimmers with increased confinement is also consistent with results from other models involving confinement (Reigh & Lauga Reference Reigh and Lauga2017; Nganguia et al. Reference Nganguia, Zhu, Palaniappan and Pak2020). By comparing figure 5(a–c), increasing the squirmer’s aspect ratio at a fixed Bingham number results in an expansion of the yielded region around the squirmer. Consequently, the ellipsoidal squirmer moves faster than the spherical squirmer for a given Bingham number. To gain further insight, we display in figure 5(d) the volume of the yielded region $V_Y$ normalised by the squirmer volume $V_p$ versus the Bingham number. It can be observed that the volume of the yielded region $V_Y$ decreases with $Bi$ and increases as the squirmer’s aspect ratio becomes larger. At low yield stress, the volume of the yielded regime declines as $V_Y \sim Bi^{-3/4}$ for all aspect ratios. This scaling can be derived by assuming ${V_Y} \sim {r^3_Y}$ , with $r_Y$ being the radius of the yielded regime (i.e. ${\left | \boldsymbol{\tau } \right |_{r = {r_Y}}} = Bi$ ). As will be discussed in the next paragraph, the velocity disturbance induced by a neutral squirmer in the far field decays as $r^{-3}$ in the limit $Bi \to 0$ , where $r$ is the distance. The derivative of the velocity disturbance decays as $r^{-4}$ , so $\left | \boldsymbol{\tau } \right |$ also decays as $\left | \boldsymbol{\tau } \right | \sim {r^{ - 4}}$ , leading to $r_Y^{ - 4} \sim Bi$ . As a result, we obtain the scaling of ${V_Y} \sim B{i^{ - 3/4}}$ . However, as the Bingham number increases, the yield surface progressively shrinks toward the squirmer’s surface, and the adopted scaling assuming the far-field decay tendency fails to hold for $Bi \gg O(10^{-1})$ . Specifically, the more elongated squirmer has a greater impact on the shape of the yield surface, producing a much larger yielded region. Thus, in high $Bi$ environments, the yielded region around an ellipsoidal squirmer is significantly larger than that around a spherical squirmer, resulting in a much higher swimming speed for the ellipsoidal squirmer compared with the spherical one.

Figure 5.

Yield-surface profiles around a neutral swimmer ( $\beta =0$ ) for various $Bi$ . Panels show (a) $a_r=1$ , (b) $a_r=3$ , (c) $a_r=6$ . (d) Volume of the yielded region as a function of $Bi$ . The dashed line in (d) is added to represent the slope $Bi^{-3/4}$ . The yield surface is given by the iso-surface of $|\boldsymbol{\tau }|=Bi$ .

The flow decay around a microswimmer is important for discussing the interaction of microswimmers and tracer diffusion. Since the stresslet of a neutral squirmer is zero, the flow decays as $\sim{\kern-1pt}r^{-3}$ in a Newtonian fluid in the far field (Theers et al. Reference Theers, Westphal, Gompper and Winkler2016). figure 6(a,b) shows the component of fluid velocity $u$ in the swimming direction versus the position relative to the squirmer centre along the semi-major axis for different aspect ratio for fixed Bingham numbers of $Bi=0.01$ and 1.0, respectively. The decay tendency of $r^{-3}$ is marginally observed for a spherical squirmer with $Bi = 0.01$ . The flow decreases significantly with increasing distance, rapidly becoming negligible in the vicinity of the yield surface. For a mild viscoplastic environment of $Bi=1.0$ , the velocity reduces to a negligible value (i.e. $u/U \lt 10^{-4}$ ) at $x/b_x \lt 1.6$ , implying that near-field fluid mechanics dominates. When we fix $a_r=3$ instead, and vary the Bingham number (figure 6 c), we observe that the flow decreases significantly with increasing Bingham number. This trend is consistent with the results in figure 2(a) and figure 4(a) that show that the propulsion speed decreases monotonically with increasing Bingham number. The results in figure 6 differ from the observed behaviour in Newtonian fluids and heterogeneous media (Demir et al. Reference Demir, van Gogh, Palaniappan and Nganguia2024). In viscoplastic fluids, the flow decays much faster than in Newtonian fluids and heterogeneous media, where the flow decays as $\sim{\kern-1pt}r^{-3}$ in the far field, regardless of the shape of the microswimmer.

Figure 6.

Laboratory-frame velocity component in the swimming direction scaled by the swimming speed $U$ versus the position relative to the centre of squirmer along the semi-major axis. (a,b) Effect of aspect ratio for (a) $Bi=0.01$ and (b) $Bi=1.0$ . (c) Effect of the Bingham number on velocity for $a_r=3$ . Here, $x/b_x=1$ corresponds to the particle’s surface. The black dash line is added to represent the slope $\sim{\kern-1pt}r^{-3}$ . The vertical dash-dot lines indicate the corresponding positions of the yield surface, where the flow properties undergo rapid transition.

To better understand the flow convergence near the yield surface, we further examine the fluid velocity profiles in this region. As described figure 7, the flow increases from the yield surface following $r^3$ apart from the vicinity of the squirmer, where $r$ is the distance from the wall. This scaling is regardless of either the Bingham number or squirmer’s aspect ratio, offering a consistent description of the near-yield-surface flow behaviour. From figure 7(b,c), the fluid velocity can be reasonably approximated as $u=A(x_Y-x)^3$ in the whole yielded region for $Bi \geqslant 1$ . Here $A$ is a coefficient related to the squirmer’s aspect ratio and Bingham number. To understand the $r^3$ growth, consider how the velocity disturbance induced by a point force in a Newtonian fluid grows in the vicinity of an infinite flat plate. The Blakelet solution (Blake Reference Blake1975) provides $r^2$ growth of the velocity disturbance in a Newtonian fluid near a flat plate, in the direction perpendicular to the wall, which is slower than in the Bingham fluid. Since the apparent viscosity of the Bingham fluid increases as the velocity gradient decreases, the viscosity increases near the wall and the velocity decays rapidly. This can be one reason for the faster growth of the velocity in the Bingham fluid. While universal theory is not currently available, it will be an important focus of our future work.

Figure 7.

Laboratory-frame velocity component in the swimming direction scaled by the swimming speed $U$ versus the position relative to the position of the yield surface along the semi-major axis. (a,b) Effect of aspect ratio for (a) $Bi=0.01$ and (b) $Bi=1.0$ . (c) Effect of the Bingham number on velocity for $a_r=3$ . Here, $x_Y$ represents the relative coordinate of the yield surface in the x-direction, thus $(x_Y -x)/(x_Y-b_x) = 0$ corresponds to the position of the yield surface. The black dash line is added to represent the slope $\sim{\kern-1pt}r^{3}$ .

3.2. Pullers and pushers

In Stokes flow of a Newtonian fluid, the swimming speed of a squirmer is independent of the swimming dipolarity $\beta$ . In contrast, in a viscoelastic fluid, the swimming speed depends on the value of $\beta$ . For a spherical squirmer in an upper convected Maxwell fluid, the speed is given by $U/U_N=1-0.2\beta Wi$ , where $Wi$ is the Weissenberg number (Datt & Elfring Reference Datt and Elfring2019). This equation indicates that the swimming speed decreases for positive values of $\beta$ and increases for negative values. It is reasonable to discuss the swimmer type dependence on the locomotion of the squirmer in a viscoplastic fluid. To address this question, we study cases of a puller/pusher ( $\left | \beta \right |=3$ ) with higher polar modes $B_{n\gt 2}$ and $C_n$ set to zero. We examine the velocity fields around the pushers, pullers and neutral swimmers in a Bingham environment in figure 8, where the yield surface is shown with the dash lines. By comparing figure 8(a,b,c) or (d,e, f), one can observe that the flow field and the yield surface are symmetrical for $\beta =3$ and −3, but different for neutral swimmers. We compare the swimming speed, power dissipation and efficiency among pullers, pushers and neutral swimmers in table 2. We find that these values differ between neutral swimmers and pullers or pushers, but are independent of the sign of $\beta$ (i.e. puller versus pusher). In some cases, the puller/pusher may swim faster than the neutral swimmer, but in others, it may swim slower. However, its swimming efficiency is much lower than that of a neutral swimmer, as the second polar mode contributes to mixing of the surrounding fluid and increases power dissipation.

Figure 8.

The velocity fields (colours) and streamlines (black lines) for (a,b,c) $Bi=1$ and $a_r=3$ , (d,e, f) $Bi=10$ and $a_r=6$ . Panels show (a,d) $\beta =-3$ , (b,e) $\beta =3$ and (c, f) $\beta =0$ . The green or blue lines represent the yield surface, while the black lines are the streamlines.

Table 2.

List of the normalised swimming speed, power dissipation and swimming efficiency for specific cases of a pusher, puller and neutral swimmer in a viscoplastic fluid.

Figure 9 shows the normalised propulsion speed, power dissipation and swimming efficiency of a puller ( $\beta =3$ ) as a function of aspect ratio. The results show that the puller uses more energy and swims more slowly in viscoplastic fluids. In addition, ellipsoidal pullers swim faster with increasing aspect ratio for a given Bingham number. By comparing figures 9(b) and 3(b), we observe that ellipsoidal pullers expend more energy compared with spherical pullers, while ellipsoidal neutral swimmers expend less energy compared with spherical swimmers. As illustrated in figure 9(c), the hydrodynamic efficiency in viscoplastic fluids is non-monotonic: it first decreases as $a_r$ approaches 1.5–2, and then increases for $a_r\gt 2$ .

Figure 9.

(a) Normalised swimming speed, (b) normalised power dissipation and (c) swimming efficiency, of a puller ( $\beta =3$ ), as a function of aspect ratio. Curves are obtained from Demir et al. (Reference Demir, van Gogh, Palaniappan and Nganguia2024) for swimming in Newtonian fluids, while the symbols denote numerical results.

In order to highlight the dependence of the squirmer’s aspect ratio and Bingham number on the swimming speed variation due to the second polar mode, we plot a phase map in figure 10(a), which shows the regimes where speed-enhanced ( $U_{\beta =3}/U_{\beta =0} \gt 1$ ) and hindered ( $U_{\beta =3}/U_{\beta =0} \lt 1$ ) swimming occur in the $a_r$ – $Bi$ space. By applying data fitting, the transition condition for $U_{\beta =3}/U_{\beta =0}=1$ is represented by a dashed curve in the figure. It is observed that the second polar mode only reduces the swimming speed of elongated squirmers in a strongly viscoplastic environment when $a_r$ exceeds a critical value, while in other scenarios, it enhances the speed. The ratio $U_{\beta =3}/U_{\beta =0}$ as a function of the aspect ratio $a_r$ is shown in figure 10(b). As a general trend, $U_{\beta =3}/U_{\beta =0}$ decreases with increasing $a_r$ for a fixed value of Bingham number. Specifically, the speed of a spherical puller is approximately 1.7 times greater than that of a spherical neutral swimmer at $Bi=100$ , whereas the speed of an elongated puller ( $a_r \gt 4$ ) is only half that of the neutral squirmer. To understand this behaviour, we compare the volume of the yielded regime for pullers and neutral swimmers in figure 10(c) as a function of $Bi$ . The results indicate that the volume of the yielded region depends on both the Bingham number and the body shape. In general, the second polar mode consistently increases the volume of the yielded region for a spherical squirmer, leading to a speedup of the squirmer. However, for ellipsoidal squirmers with aspect ratios exceeding a certain value, the yielded region created by the polar mode is much smaller than that of neutral squirmers at large Bingham numbers ( $Bi \geqslant 10$ ). As a result, the increased confinement slows down the ellipsoidal squirmer.

Figure 10.

Ratio of swimming speed between a puller ( $\beta =3$ ) and a neutral squirmer in a viscoplastic fluid. (a) Phase diagram over the aspect ratio $a_r$ and Bingham number $Bi$ . Symbols are coloured by the speed ratio $U_{\beta =3}/U_{\beta =0}$ , and the dashed line gives the boundary between these states by data fitting. (b) Value of $U_{\beta =3}/U_{\beta =0}$ as a function of the Bingham number $Bi$ as a function of the aspect ratio. (c) Normalised volume of yielded region $V_{Y,\beta =3}/V_{Y,\beta =0}$ as a function of $Bi$ .

3.3. Neutral swirling squirmer

The swirling mode plays a role in complex fluids such as viscoelastic fluids. Theoretical and numerical studies have indicated that the swimming speed of squirmers is enhanced due to the coupling between chirality and viscoelasticity (Binagia et al. Reference Binagia, Phoa, Housiadas and Shaqfeh2020; Housiadas et al. Reference Housiadas, Binagia and Shaqfeh2021). In particular, previous studies highlight the critical role of the first normal stress difference induced by the swirling flow in boosting swimming speeds (Binagia & Shaqfeh Reference Binagia and Shaqfeh2021; Kobayashi et al. Reference Kobayashi, Molina and Yamamoto2024). To date, however, no one has considered how the presence of the azimuthal modes in the squirmer model impacts the swimming kinematics in a viscoplastic fluid, and studies of swirling squirmers in complex fluids have been limited to spherical shapes.

In the $a_r$ – $Bi$ plane plotted in figure 11(a), we report the change in swimming speed when introducing the second swirling mode ( $\chi =3$ ) relative to its non-swirling counterpart (i.e. $U_{\chi =3}/U_{\chi =0}$ ) for different parameter points. Similar to our observation in figure 10(a), while the swirling motion generally contributes to an increase in the swimming speed of a squirmer for small $a_r$ , beyond a given threshold in aspect ratio and/or Bingham number, the swirling squirmer swims slower than the non-swirling swimmer with the same $a_r$ and $Bi$ . The transition condition for $U_{\chi =3}/U_{\chi =0}=1$ is fitted and represented by a dashed curve in the figure. In figure 11(b), we show that $U_{\chi =3}/U_{\chi =0}$ decreases with the aspect ratio, similar to the results observed in figure 10(b). Specifically, the speed enhancement of a spherical swirling squirmer relative to a non-swirling squirmer is approximately 1.7 times at $Bi=100$ , while the speed of ellipsoidal squirmers with large aspect ratios is significantly reduced by the swirling mode. We plot the variation in the normalised volume of the yielded region ( $V_{Y,\chi =3}/V_{Y,\chi =0}$ ) in figure 11(c). For all aspect ratios, $V_{Y,\chi =3}/V_{Y,\chi =0}$ first increases with increasing $Bi$ , reaching a local maximum at $Bi=O(10^{-1})$ , then becomes smaller as $Bi$ continues to increase. For a spherical squirmer, the swirling mode always increases the volume of the yielded region. However, for ellipsoidal squirmers, the yielded region generated by swirling squirmers can be much smaller than that of non-swirling squirmers at large Bingham numbers. As a result, the increased confinement reduces the speed of the squirmer with a large aspect ratio.

Figure 11.

Ratio of swimming speed between a swirling squirmer ( $\chi =3$ ) and a non-swirling neutral squirmer in a viscoplastic fluid. (a) Phase diagram over the aspect ratio $a_r$ and Bingham number $Bi$ . Symbols are coloured by the speed ratio $U_{\chi =3}/U_{\chi =0}$ , and the dashed line gives the boundary between these states by data fitting. (b) Value of $U_{\chi =3}/U_{\chi =0}$ as a function of the aspect ratio. (c) Normalised volume of yielded region $V_{Y,\chi =3}/V_{Y,\chi =0}$ as a function of the Bingham number.

The stress distributions in figures 12 and 13 for $a_r=1.5$ and 6 can help us to understand the observed dependence of the yielded region on the swirling mode. In the case of $a_r=1.5$ and $Bi=10$ , the magnitude of both normal and shear stresses around a swirling squirmer is much larger than that around a non-swirling squirmer, as shown in figure 12. Therefore, the yielded region around a swirling squirmer is larger than that around a non-swirling squirmer. Unlike the case of $a_r=1.5$ and $Bi=10$ , for the case of $a_r=6$ and $Bi=100$ , the magnitude of both normal and shear stresses around a swirling squirmer is smaller than that around a non-swirling squirmer, as shown in figure 13, resulting in a reduced yielded region for the swirling squirmer. The mechanism by which the swirling mode reduces stress around a squirmer with a large aspect ratio remains unclear and warrants further investigation.

Figure 12.

Stress fields around the squirmer ( $a_r=1.5$ ) in a viscoplastic fluid with $Bi=10$ in the $x$ – $y$ plane: (a,b) swirling swimmer $(\chi =3)$ and (c,d) non-swirling swimmer. Panels (a,c) show the normal stress $\tau _{xx}$ , and panels (b,d) show the shear stress $\tau _{xy}$ . The dashed line represents the yield surface.

Figure 13.

Stress fields around the squirmer ( $a_r=6$ ) in a viscoplastic fluid with $Bi=100$ in the $x$ – $y$ plane: (a,b) swirling swimmer $(\chi =3)$ and (c,d) non-swirling swimmer. Panels (a,c) show the normal stress $\tau _{xx}$ , and panels (b,d) show the shear stress $\tau _{xy}$ . The dashed line represents the yield surface.

3.4. Force analysis

The dimensionless net force acting on a force-free squirmer is given by

(3.4) \begin{equation} \boldsymbol{F} = \underbrace {\int _S {\left ( { - p} \right){\unicode{x1D644}} \cdot {\boldsymbol{n}}\ {\textrm {d}}S} }_{{{\boldsymbol{F}}^{pres}}} + \underbrace {2\int _S {\unicode{x1D640}} \cdot {\boldsymbol{n}}\ {\textrm {d}}S}_{{{\boldsymbol{F}}^{visc}}} + \underbrace {{Bi}\int _S {{\unicode{x1D63D}} \cdot {\boldsymbol{n}}\ {\textrm {d}}S} }_{{{\boldsymbol{F}}^{plas}}} = {\boldsymbol{0}}, \end{equation}

where $\boldsymbol{F}^{pres}$ is the contribution due to pressure, $\boldsymbol{F}^{visc}$ results from the viscous stress and $\boldsymbol{F}^{plas}$ represents the plastic force.

In a Newtonian fluid, the force can also be decomposed into contributions of the thrust force and towing force as $\boldsymbol{F}={\boldsymbol{F}_{thrust}} + {\boldsymbol{F}_{tow}} = 0$ . Here, $\boldsymbol{F}_{thrust}$ is the thrust force generated by the squirmer when it is being held fixed, while ${\boldsymbol{F}_{tow}}=\alpha \boldsymbol{U}$ is the towing force required to tow the squirmer at the swimming speed $\boldsymbol{U}$ with $\alpha$ being the translational drag coefficient (Yang et al. Reference Yang, Lu, Zhao, Kawamura and Chalmers2017). However, this decomposition in complex fluids is not entirely rigorous due to the nonlinearity of the rheology, i.e. $\boldsymbol{F}\ne {\boldsymbol{F}_{thrust}} + {\boldsymbol{F}_{tow}}$ , as discussed in Datt et al. (Reference Datt, Zhu, Elfring and Pak2015). We have carefully compared the thrust and towing forces in our settings and found that the difference between them is approximately 10 % in magnitude at $Bi=10$ . Although such a decomposition is no longer rigorous in a viscoplastic fluid, we can still gain some insight into the propulsion behaviour by analysing the trust and drag forces separately. In particular, the comparison of the contributions of pressure, viscosity and plasticity and the comparison of the contributions of thrust and drag coefficients are helpful to facilitate our qualitative understanding.

In this subsection, we discuss changes in the thrust forces that act to generate propulsion of the squirmer and how they vary with aspect ratio. This is done by analysing the force contributions parallel to the swimming axis on a fixed squirmer. Figure 14(a–c) illustrates the variation of the decomposed thrust forces as a function of the aspect ratio for different swimmer types at $Bi=10$ . We find that the total thrust forces do not change much for various aspect ratios, suggesting that the variation in free-swimming speed, $U$ , is primarily attributed to changes in swimming drag associated with squirmer’s aspect ratio. Moreover, the viscous, plastic and pressure contributions reveal a strong dependence on the cell shape. For all swimmer types, the viscous and plastic forces significantly increase with the aspect ratio, while the pressure force shows a decrease with $a_r$ . By comparing figure 14(a–c), we observe that the neutral squirmer generates a larger thrust force than the puller/pusher or swirling squirmer. In figure 14(d), we provide the translational drag coefficient $\alpha ={F}_{tow}/{U}$ , which decreases with the aspect ratio. Meanwhile, in figure 14(e), we observe that the volume of the yielded region around the fixed squirmer increases with the aspect ratio. This is because the reduction in the confinement effect causes a decrease in the translational drag coefficient as the aspect ratio increases. This ultimately leads to a higher swimming speed as the aspect ratio increases.

Figure 14.

(a–c) Thrust force contributions as a function of aspect ratio at $Bi=10$ . Panel (a) neutral swimmer, $\beta =\chi =0$ , panel (b) puller, $\beta =3, \chi =0$ , panel (c) neutral swirling swimmer, $\beta =0,\chi =3$ . (d) Translational drag coefficient of different type of squirmers versus aspect ratio. The forces are normalised by $\mu B_1 r_{eq}$ . (e) Volume of yielded region of squirmers versus aspect ratio.

For a neutral squirmer with constant-volume constraints, the results in figure 14 are quite different from the observed behaviour in heterogeneous media (Demir et al. Reference Demir, van Gogh, Palaniappan and Nganguia2024). In heterogeneous fluids, the thrust force increases with the aspect ratio, and the drag coefficient shows a non-monotonic trend as a function of the aspect ratio: $\alpha$ generally increases with increasing the aspect ratio for $a_r \gt 2.3$ . In summary, the increase in swimming speed with aspect ratio in a heterogeneous fluid is primarily attributed to the enhanced thrust force, whereas in a viscoplastic fluid, the increase in swimming speed results from a decrease in the drag coefficient.

4.1. Motion of fluid particles

As depicted in figure 15(a), we consider an axially symmetric squirmer that moves in a straight line along the $x$ direction with a constant velocity $U$ in a viscoplastic fluid. Here, $L_x$ , $L_y$ and $L_z$ represent the domain sizes in the $x$ , $y$ and $z$ directions, respectively. We define $A=L_yL_z$ as the cross-sectional area of the region, which is equal to $ ( 32r_{eq})^2$ . The simulation is carried out to the extent that this cross-sectional area is larger than the yielded region. The area farther away from this region is solidified by plasticity, so the motion of the fluid particles, including Brownian motion, is assumed to be negligible.

Figure 15.

(a) Numerical settings for the calculation of squirmer-induced diffusion. (b–d) Typical trajectories of tracers as a neutral squirmer moves along the $x$ axis. Panel (b), $a_r=3,Bi=0.01$ ; (c) $a_r=3,Bi=1.0$ ; (d) $a_r=8,Bi=1.0$ . The starting and ending points of the tracer trajectories are marked by filled black circles and open red circles, respectively. Here, $y=0$ represents the swimming axis of the squirmer.

The number of sample tracers used in the simulation is $N=1000^2$ , and their initial positions, ${\boldsymbol{r}}_i^{initial}$ , are set uniformly at random in the mid-plane, where the flow is initially unyielded. After the squirmer passes away, the sample tracers move to ${\boldsymbol{r}}_i^{final}$ in the unyielded region and stop moving. Particles are assumed to move with the velocity created by the squirmer and the effect of Brownian motion is ignored. The displacement of the $i$ th tracer is defined as $\Delta {{\boldsymbol{r}}_i} = {{\boldsymbol{r}}_i^{final}} - {\boldsymbol{r}}_i^{initial}$ , and the ensemble-averaged squared displacement tensor can be calculated by

(4.1) \begin{equation} {\unicode{x1D64D}}_2 = \frac {1}{N}\sum \limits _{i = 1}^N {\left ( {\Delta {{\boldsymbol{r}}_i} \otimes \Delta {{\boldsymbol{r}}_i}} \right)} = \left [ {\begin{array}{*{20}{c}} {r_\parallel ^2}&0&0\\ 0&r_ \bot ^2&0\\ 0&0&r_ \bot ^2 \end{array}} \right], \end{equation}

where $r_\parallel$ and $r_ \bot$ are the average displacements over the $N$ tracers parallel to (longitudinal), and in the plane perpendicular to (transverse), the swimming direction. The velocities of the tracers are interpolated from the velocity fields, and the time integration is performed using a second-order Runge–Kutta method. The accuracy of the calculated displacements is guaranteed by the number of particles and the database resolution used in this study. The representative trajectories of the tracer particles for specific cases ( $\beta =0$ ) are shown in figure 15(b–d). At a low Bingham number of $Bi=0.01$ , the tracers initially positioned near the swimming axis tend to move forward, while those farther from the axis are displaced backwards. In contrast, at $Bi=1$ , all tracers are displaced backwards. In general, the longitudinal average squared displacement $r^2_\parallel$ dominates over the transverse average squared displacement $r^2_ \bot$ , hence we focus on the longitudinal squared displacement in the following discussion.

Figure 16 presents the scaled longitudinal average displacement, ${r_\parallel ^2}/r^2_{eq}$ , as a function of $a_r$ and $Bi$ . Figure 16(a) reveals that, when the aspect ratio exceeds 2, $r^2_\parallel$ decreases as the aspect ratio increases. This can be understood by comparing the tracer trajectories in figure 15(c,d) for a neutral squirmer ( $Bi=1$ , $a_r=3$ and 8) in motion. As the aspect ratio increases, a decrease in the backward displacement is observed for tracers initially positioned near the swimming axis, leading to a reduction in the longitudinal average squared displacement. In figure 16(b), we find that $r^2_\parallel$ generally increases with the Bingham number. As presented in figure 15(b,c), an increase in the backward displacement of tracers is observed as $Bi$ increases, then the longitudinal average squared displacement increases.

Figure 16.

Dependence of the scaled longitudinal average displacement, ${r_\parallel ^2}/r^2_{eq}$ , on (a) the aspect ratio and (b) the Bingham number for a neutral squirmer.

4.2. Diffusion coefficient

The diffusivity is a measure of the spreading of the tracer particles. The definition of the diffusion tensor is given by

(4.2) \begin{equation} \unicode{x1D63F} = \lim _{t \rightarrow \infty } \frac {\langle \left [\boldsymbol{r}(t)-\boldsymbol{r}(0)\right] \otimes \left [\boldsymbol{r}(t)-\boldsymbol{r}(0)\right]\rangle }{2t}, \end{equation}

where $\boldsymbol{r}$ is the position of a fluid particle, and the angle bracket denotes an ensemble average over particles in the whole domain.

Here, we focus on fluid particle diffusion in a dilute suspension of squirmers, allowing us to neglect multi-body interactions of squirmers. Since the motion of tracers in a plastic fluid is restricted to the vicinity of the squirmer, the tracer particles only move when a swimmer comes close and otherwise remain stationary. Furthermore, we assume an isotropic suspension, where the orientation of the squirmers is isotropic. Under these assumptions, the tracer exhibits a three-dimensional random walk, taking steps only when the squirmer comes close. The diffusion tensor becomes isotropic, and given by

(4.3) \begin{equation} \unicode{x1D63F} = D \unicode{x1D644} = \left ( \left \langle \frac {r_\parallel ^2 + 2 r_ \bot ^2}{3} \right \rangle _{\! \! \! A} \frac {U \! A}{2 V_p} \phi \right) \unicode{x1D644}, \end{equation}

where $\langle \ \rangle _{\! A}$ indicates the average over the cross-sectional area $A$ , and $\phi$ is the volume fraction of squirmers in a dilute suspension. This definition indicates that the diffusivity is linearly related to the total squared displacement of tracers, the swimming speed $U$ and the volume fraction $\phi$ .

Figure 17 shows the numerically calculated squirmer-induced diffusivity, $D$ , as a function of $a_r$ and $Bi$ . The results reveal a non-monotonic behaviour of diffusivity with respect to both the aspect ratio and the Bingham number. As these parameters increase, the diffusivity first reaches a maximum and then decreases. The increase in $D$ with $a_r$ in figure 17(a) is attributed to the increase in swimming speed shown in figure 3(a), while the decrease in $D$ with further increased $a_r$ results from the reduction in the average squared displacement in figure 16(a). In contract, the increase in $D$ with $Bi$ for $Bi \leqslant 1$ is driven by the increase in the average squared displacement, while the decrease in $D$ with $Bi$ at higher Bingham numbers is due to the reduced swimming speed, as seen in figure 2(a).

We then investigated the effect of the swimming dipolarity $\beta$ ( $\left | \beta \right | \leqslant 1$ ) on the squirmer-induced diffusivity $D$ . As shown in figure 18, we see that both $r_\parallel ^2$ and $D$ reach minimum values when $\beta =0$ , i.e. for neutral squirmer, and increase as $\left | \beta \right |$ increased. This result is qualitatively consistent with previous studies for a Newtonian fluid (Pushkin et al. Reference Pushkin, Shum and Yeomans2013; Thiffeault Reference Thiffeault2015). Figure 19 shows the representative trajectories of tracer particles for the spherical cases ( $a_r=1$ and $Bi=0.1$ ). As $\left | \beta \right |$ increases, the tracers initially positioned near the swimming axis tend to move forward with larger displacement, thereby enhancing the average longitudinal squared displacement and diffusivity.

Figure 17.

Dependence of squirmer-induced diffusivity on (a) the aspect ratio and (b) the Bingham number for a neutral squirmer.

Figure 18.

Effect of $\beta$ on (a) the longitudinal average displacement and (b) the squirmer-induced diffusivity.

Figure 19.

Typical trajectories of tracers as a squirmer ( $a_r=1$ and $Bi=0.1$ ) moves along the $x$ axis. Panel (a) $\beta =0$ ; (b) $\beta =-1$ ; (c) $\beta =1$ . The starting and ending points of the tracer paths are marked by filled black circles and open red circles, respectively.