2.1. Shear-thinning rheology: the Carreau–Yasuda model

Shear-thinning fluids experience a loss in apparent viscosity with applied strain rates, a property that results from changes in the fluid microstructure. As the rate of strain exceeds the rate of structural relaxation, one observes microstructural ordering in the fluid (Brown & Jaeger Reference Brown and Jaeger2011). Here, we capture the change in apparent viscosity due to this ordering using the Carreau–Yasuda model for generalized Newtonian fluids (Bird, Armstrong & Hassager Reference Bird, Armstrong and Hassager1987). The variation of viscosity with applied strain rate is given by

(2.1) $$\begin{eqnarray}{\it\eta}={\it\eta}_{\infty }+({\it\eta}_{0}-{\it\eta}_{\infty })[1+{\it\lambda}_{t}^{2}|\dot{{\it\gamma}}|^{2}]^{(n-1)/2},\end{eqnarray}$$

where ${\it\eta}_{0}$ and ${\it\eta}_{\infty }$ are the zero- and infinite-shear-rate viscosities respectively. The power-law index $n$ characterizes the degree of shear thinning ( $n<1$ ) and the relaxation time ${\it\lambda}_{t}$ sets the crossover strain rate at which non-Newtonian behaviour starts to become significant. The magnitude of the strain rate tensor is given by $|\dot{{\it\gamma}}|=({\it\Pi}/2)^{1/2}$ , where ${\it\Pi}=\dot{{\it\gamma}}_{ij}\dot{{\it\gamma}}_{ij}$ is the second invariant of the tensor. As an example, measured values for human cervical mucus can be fitted by the Carreau–Yasuda model with values of ${\it\eta}_{0}=145.7~\text{Pa}~\text{s}$ , ${\it\eta}_{\infty }=0~\text{Pa}~\text{s}$ , ${\it\lambda}_{t}=631.04~\text{s}$ and $n=0.27$ (Hwang, Litt & Forsman Reference Hwang, Litt and Forsman1969; Vélez-Cordero & Lauga Reference Vélez-Cordero and Lauga2013).

We non-dimensionalize the flow quantities taking the first mode, $B_{1}$ , of the surface actuation as the scale for velocity and the radius, $a$ , of the squirmer as the characteristic length scale. The strain rates are scaled with ${\it\omega}=B_{1}/a$ and the stresses by ${\it\eta}_{0}{\it\omega}$ , such that the constitutive equation takes the dimensionless form

(2.2) $$\begin{eqnarray}{\bf\tau}^{\ast }=\left\{{\it\beta}+(1-{\it\beta})[1+Cu^{2}|\dot{{\it\gamma}}^{\ast }|^{2}]^{(n-1)/2}\right\}\dot{{\bf\gamma}}^{\ast },\end{eqnarray}$$

where ${\bf\tau}$ is the deviatoric stress tensor, and dimensionless quantities are denoted by stars ( $^{\ast }$ ). The Carreau number $Cu={\it\omega}{\it\lambda}_{t}$ is the ratio of the characteristic strain rate, defined by the surface actuation ${\it\omega}$ , to the crossover strain rate, defined by the fluid relaxation $1/{\it\lambda}_{t}$ . The viscosity ratio is given by ${\it\beta}={\it\eta}_{\infty }/{\it\eta}_{0}\in [0,1]$ .

It is evident from (2.2) that when the actuation rate ${\it\omega}$ is much smaller or much larger than the fluid relaxation rate $1/{\it\lambda}_{t}$ , i.e. when $Cu\rightarrow 0$ or $Cu\rightarrow \infty$ , the shear-thinning fluid reduces to a Newtonian fluid of constant viscosity ${\it\eta}_{0}$ (dimensionless viscosity 1) or ${\it\eta}_{\infty }$ (dimensionless viscosity ${\it\beta}$ ) respectively. Recalling that for a given surface actuation the swimming speed of a squirmer in the Newtonian regime is independent of the fluid viscosity, we therefore expect the swimming speed of a squirmer in a shear-thinning fluid to converge to its Newtonian value in the limits $Cu\rightarrow 0$ or $Cu\rightarrow \infty$ . Non-monotonic variation of the swimming speed with $Cu$ is expected for any swimmer with prescribed kinematics and has been observed by Montenegro-Johnson et al. (Reference Montenegro-Johnson, Smith, Smith, Loghin and Blake2012) for some two-dimensional model swimmers. In this study, we employ both asymptotic analysis and numerical simulations to investigate these effects of shear-thinning rheology on swimming at low Reynolds numbers.

2.2. Asymptotic analysis

The deviatoric stress tensor ${\bf\tau}^{\ast }$ in (2.2) is a nonlinear function of the strain rate tensor $\dot{{\bf\gamma}}^{\ast }$ . Assuming only a weak nonlinearity, we may uncouple the Newtonian and non-Newtonian contributions, writing

(2.3) $$\begin{eqnarray}{\bf\tau}^{\ast }=\dot{{\bf\gamma}}^{\ast }+{\it\varepsilon}\unicode[STIX]{x1D63C}^{\ast },\end{eqnarray}$$

with ${\it\varepsilon}\ll 1$ as a dimensionless measure of the deviation from the Newtonian case ( ${\it\varepsilon}=0$ ).

We observe that in the limits $Cu=0$ or ${\it\beta}=1$ , (2.2) reduces to a Newtonian constitutive equation. Thus, one may expect weakly nonlinear behaviour when the fluid relaxation rate is much faster than the surface actuation rate ( ${\it\varepsilon}=Cu^{2}\ll 1$ ), or when the zero-shear-rate viscosity is very close to the infinite-shear-rate viscosity ( ${\it\varepsilon}=1-{\it\beta}\ll 1$ ).

Henceforth we shall work in dimensionless quantities and therefore drop the stars ( $^{\ast }$ ) for convenience.

2.2.1. Expansion in Carreau number

Expanding all fields in regular perturbation series in ${\it\varepsilon}=Cu^{2}$ , we obtain, order by order, the constitutive equations

(2.4) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}{\bf\tau}_{0}=\dot{{\bf\gamma}}_{0},\\ {\bf\tau}_{1}=\dot{{\bf\gamma}}_{1}+\displaystyle \frac{(n-1)}{2}(1-{\it\beta})|\dot{{\it\gamma}_{0}}|^{2}\dot{{\bf\gamma}}_{0},\end{array}\right\}\end{eqnarray}$$

hence $\unicode[STIX]{x1D63C}=(n-1)(1-{\it\beta})|\dot{{\it\gamma}_{0}}|^{2}\dot{{\bf\gamma}}_{0}/2$ to leading order in (2.3). It should be noted that the first correction to the Newtonian behaviour is linear in $n$ , which points to a linear dependence of the swimming speed on $n$ upon using (2.8), elucidating the trend suggested by the two-dimensional numerical findings in Montenegro-Johnson et al. (Reference Montenegro-Johnson, Smith and Loghin2013). We also remark that this expansion is valid only when $Cu^{2}|\dot{{\it\gamma}}|^{2}$ is $o(1)$ and is therefore not uniformly valid across all values of strain rates.

2.2.2. Expansion in viscosity ratio

Expanding in perturbation series with ${\it\varepsilon}=1-{\it\beta}$ gives us, order by order, the constitutive equations

(2.5) $$\begin{eqnarray}\displaystyle & {\bf\tau}_{0}=\dot{{\bf\gamma}}_{0}, & \displaystyle\end{eqnarray}$$

(2.6) $$\begin{eqnarray}\displaystyle & {\bf\tau}_{1}=\dot{{\bf\gamma}}_{1}+\left\{-1+(1+Cu^{2}|\dot{{\it\gamma}}_{0}|^{2})^{(n-1)/2}\right\}\dot{{\bf\gamma}}_{0}, & \displaystyle\end{eqnarray}$$

where in this limit $\unicode[STIX]{x1D63C}=\{-1+(1+Cu^{2}|\dot{{\it\gamma}}_{0}|^{2})^{(n-1)/2}\}\dot{{\bf\gamma}}_{0}$ to leading order in (2.3). It should be noted that this asymptotic expansion is uniformly valid for all strain rates or Carreau numbers, which permits a full-range study of the non-monotonic swimming behaviour.

2.3. The reciprocal theorem

Stone & Samuel (Reference Stone and Samuel1996) demonstrated the use of the Lorenz reciprocal theorem in low-Reynolds-number hydrodynamics (Happel & Brenner Reference Happel and Brenner1983) to obtain the swimming velocity of a squirmer for a given prescribed surface actuation $\boldsymbol{u}^{S}$ without calculation of the unknown flow field, provided one can solve the resistance/mobility problem for the swimmer shape (with surface $S$ ). Lauga (Reference Lauga2009, Reference Lauga2014) then developed integral theorems extending this method for use with complex fluids. We use these methods in the subsequent calculations to obtain the swimming velocity of a squirmer in a shear-thinning fluid; the methodology adopted below closely follows the formulation in Elfring & Lauga (Reference Elfring, Lauga and Spagnolie2015).

We represent the velocity field and the associated total stress tensor for a force- and torque-free swimmer with $\boldsymbol{u}$ and ${\bf\sigma}$ respectively. We consider the corresponding resistance problem in a Newtonian fluid to simplify the calculation of the swimming velocity. The resistance problem (denoted with a hat) involves the rigid-body motion with translational velocity $\hat{\boldsymbol{U}}$ and rotational velocity $\hat{{\it\bf\Omega}}$ , and the corresponding velocity field and associated stress tensor are represented by $\hat{\boldsymbol{u}}$ and $\hat{{\bf\sigma}}$ respectively. Due to the linearity of the Stokes equation, we may write $\hat{\boldsymbol{u}}=\hat{\unicode[STIX]{x1D647}}\boldsymbol{\cdot }\hat{\unicode[STIX]{x1D650}}$ , $\hat{{\bf\sigma}}=\hat{\unicode[STIX]{x1D64F}}\boldsymbol{\cdot }\hat{\unicode[STIX]{x1D650}}$ and $\hat{\unicode[STIX]{x1D641}}=-\hat{\unicode[STIX]{x1D64D}}\boldsymbol{\cdot }\hat{\unicode[STIX]{x1D650}}$ . Here, for compactness, both the translational and the rotational components of velocity are contained in $\hat{\unicode[STIX]{x1D650}}$ and, similarly, the corresponding matrices contain both the translational and rotational terms. In weakly nonlinear complex fluids, the six-dimensional translational and rotational velocity of the swimmer $\unicode[STIX]{x1D650}$ is then given by

(2.7) $$\begin{eqnarray}\unicode[STIX]{x1D650}=\hat{\unicode[STIX]{x1D64D}}^{-1}\boldsymbol{\cdot }\left[\int _{S}\boldsymbol{u}^{S}\boldsymbol{\cdot }(\boldsymbol{n}\boldsymbol{\cdot }\hat{\unicode[STIX]{x1D64F}})\,\text{d}S-{\it\varepsilon}\int _{V}\unicode[STIX]{x1D63C}\boldsymbol{ : }\boldsymbol{{\rm\nabla}}\hat{\unicode[STIX]{x1D647}}\,\text{d}V\right].\end{eqnarray}$$

The integral over the volume of fluid $V$ external to $S$ in the equation measures the change in swimming dynamics due to the non-Newtonian behaviour of the fluid. For a spherical squirmer with axisymmetrical tangential surface distortions, there is no rotational motion and the translational velocity is given simply by

(2.8) $$\begin{eqnarray}\boldsymbol{U}=-\frac{1}{4{\rm\pi}}\int _{S}\boldsymbol{u}^{S}\,\text{d}S-\frac{{\it\varepsilon}}{8{\rm\pi}}\int _{V}\unicode[STIX]{x1D63C}\boldsymbol{ : }\left(1+\frac{1}{6}{\rm\nabla}^{2}\right)\boldsymbol{{\rm\nabla}}\unicode[STIX]{x1D642}\,\text{d}V,\end{eqnarray}$$

where $\unicode[STIX]{x1D642}=(1/r)(\unicode[STIX]{x1D644}+\boldsymbol{r}\boldsymbol{r}/r^{2})$ is the Oseen tensor (or Stokeslet). The first term on the right-hand side is the result of swimming in a Newtonian fluid (Stone & Samuel Reference Stone and Samuel1996), and the last term in the equation contains the weakly nonlinear effect, which can be evaluated analytically in some special cases and can be computed in general by numerical quadrature readily.

2.4. Numerical solution

The numerical simulations of the momentum equations at zero Reynolds number with the Carreau–Yasuda constitutive relation (2.1) are implemented in the finite element method software COMSOL. We use a square computational domain of size $500a\times 500a$ , discretized by approximately 30 000–50 000 Taylor–Hood ( $P2$ – $P1$ ) triangular elements. The mesh is refined near the squirmer in order to properly capture the spatial variation of the viscosity. Since slowly decaying flow fields are expected at low Reynolds numbers, a large domain size is important to guarantee accuracy. The simulations are performed in a reference frame moving with the swimmer, and the far-field (inlet) velocity is varied to obtain a computed zero force on the squirmer. In addition to comparing with the asymptotic analysis in this paper, we have validated our implementation against the analytical results for a three-dimensional squirmer in a Newtonian fluid (Lighthill Reference Lighthill1952; Blake Reference Blake1971) and a two-dimensional counterpart in a shear-thinning fluid (Montenegro-Johnson et al. Reference Montenegro-Johnson, Smith and Loghin2013).

3.1. Drag and thrust

We separate the swimming problem into a drag problem (a sphere undergoing rigid-body translation $\boldsymbol{U}$ inducing hydrodynamic drag) and a thrust problem (a sphere held fixed undergoing only tangential surface distortions thereby generating thrust). The superposition of these two subproblems gives the entire swimming problem in a Newtonian fluid in the Stokes regime; this is obviously not the case in a shear-thinning fluid due to its nonlinear constitutive equation. However, by looking at the thrust and the drag problems, one may gain insight into the more complex non-Newtonian swimming problem (Montenegro-Johnson et al. Reference Montenegro-Johnson, Smith and Loghin2013; Qiu et al. Reference Qiu, Lee, Mark, Morozov, Münster, Mierka, Turek, Leshansky and Fischer2014).

We derive the expressions for drag and thrust in a shear-thinning fluid again via the reciprocal theorem approach (§ 2.3) by utilizing the solution to the resistance problem in a Newtonian fluid. The drag force on a sphere moving with a velocity $\boldsymbol{U}$ in weakly shear-thinning fluid is given by $\boldsymbol{F}_{D}=-6{\rm\pi}\boldsymbol{U}-(3/4){\it\varepsilon}\int _{V}\unicode[STIX]{x1D63C}_{D}\boldsymbol{ : }(1+{\rm\nabla}^{2}/6)\boldsymbol{{\rm\nabla}}\unicode[STIX]{x1D642}$ , where $\unicode[STIX]{x1D63C}_{D}$ is formed by the solution to the Newtonian drag problem. Similarly, the thrust force generated by a sphere held stationary with surface actuation $\boldsymbol{u}^{S}$ in a weakly shear-thinning fluid is given by $\boldsymbol{F}_{T}=-(3/2)\int _{S}\boldsymbol{u}^{S}\,\text{d}S-(3/4){\it\varepsilon}\int _{V}\unicode[STIX]{x1D63C}_{T}\boldsymbol{ : }(1+{\rm\nabla}^{2}/6)\boldsymbol{{\rm\nabla}}\unicode[STIX]{x1D642}$ , where $\unicode[STIX]{x1D63C}_{T}$ is formed by the solution to the Newtonian thrust problem.

One could expect a drag reduction when a rigid sphere is pulled with a constant velocity through a shear-thinning fluid, since the fluid viscosity is reduced by the fluid straining motion. However, it is interesting to see in figure 2(a) that the thrust reduction caused by the shear-thinning rheology is larger than the drag reduction for a large range of $Cu$ . This more severe reduction in thrust than in drag then suggests slower swimming speeds compared with the Newtonian case, which correctly predicts the trend found by detailed calculations (figure 1). In addition, for very small or large values of $Cu$ , the difference between the magnitudes of drag and thrust ( $F_{D}-F_{T}$ ) vanishes, as shown in figure 2(b), respecting the limits where the swimming speed should recover the Newtonian value (figure 1).

Figure 2. In (a), the symbol ${\it\Delta}$ denotes the difference between drag or thrust in a shear-thinning fluid and a Newtonian fluid: the dash–dot curve represents the drag reduction in a shear-thinning fluid compared with a Newtonian fluid, while the solid and dashed curves represent the loss in thrust for squirmers with ${\it\alpha}=0$ and ${\it\alpha}=\pm 5$ respectively ( $n=0.25$ , ${\it\varepsilon}=0.1$ ). In (b), we show that the difference between drag and thrust is positive both for a neutral squirmer (solid) and for a pusher/puller (dashed) in a shear-thinning fluid. All quantities are dimensionless.

Although conceptually intuitive, the drag and thrust decomposition is not complete as it neglects the contribution of nonlinear products in the non-Newtonian stress to the full swimming problem, namely $\unicode[STIX]{x1D63C}\neq \unicode[STIX]{x1D63C}_{T}+\unicode[STIX]{x1D63C}_{D}$ , due to the nonlinearity in the constitutive equation. We will give a counterexample in § 3.2 below to illustrate a scenario when these intuitive arguments fail.

3.2. Addition of other squirming modes

The results from the detailed asymptotic and numerical analysis as well as the intuitive model for a two-mode squirmer seem to suggest that the shear-thinning rheology acts to hinder the locomotion. This raises the simple question of whether this conclusion still holds if other modes of surface actuation are present. The picture is clear for a Newtonian fluid: only the $B_{1}$ mode contributes to swimming and the addition of other modes does not alter the swimming speed. However, can the shear-thinning rheology render other modes, typically not considered in the Newtonian analysis, effective for propulsion? We address these questions using the asymptotic and numerical tools developed in the previous sections.

We first note that the $B_{3}$ mode alone leads to locomotion in a shear-thinning fluid, in stark contrast to a Newtonian fluid, as only the $B_{1}$ mode has a non-zero surface average (see (2.8)). Indeed, any odd mode alone may lead to locomotion in a shear-thinning fluid (even modes alone do not swim by symmetry). We also find quite distinctive behaviour when the $B_{3}$ mode is combined with other modes. In the $Cu\ll 1$ regime, we can derive an analytical expression allowing us to predict the values of ${\it\alpha}$ and ${\it\zeta}=B_{3}/B_{1}$ for faster or slower swimming. To quadratic order in $Cu$ , we find

(3.2) $$\begin{eqnarray}\frac{U}{U_{N}}=1+Cu^{2}(1-{\it\beta})\frac{(n-1)}{2}C_{1}[1+C_{2}(1+C_{3}{\it\zeta}){\it\alpha}^{2}+C_{4}(C_{5}{\it\zeta}^{2}+C_{6}{\it\zeta}-1){\it\zeta}],\end{eqnarray}$$

where the additional numerical constants are given by $C_{3}=0.51$ , $C_{4}=0.70$ , $C_{5}=0.18$ and $C_{6}=1.66$ . Again, the swimming speed is even in ${\it\alpha}$ and we recover (3.1) when ${\it\zeta}=0$ , as expected. From (3.2), we can predict that faster swimming ( $U/U_{N}>1$ ) will occur if

(3.3) $$\begin{eqnarray}\left.\frac{\partial ^{2}U}{\partial Cu^{2}}({\it\alpha},{\it\zeta})\right|_{Cu=0}>0,\end{eqnarray}$$

in other words when the term in the square brackets in (3.2) is negative. In figure 3(a), we plot the level set curve below which faster swimming occurs in the small- $Cu$ regime; we find that this can only occur when ${\it\zeta}$ is negative for any ${\it\alpha}$ . For example, when ${\it\alpha}=0$ we must have ${\it\zeta}<-10.11$ , while ${\it\alpha}=\pm 5$ , ${\it\zeta}<-2.22$ leads to faster swimming.

Figure 3. In (a), we show the level set curve below which faster swimming occurs for specific values of ${\it\alpha}$ and ${\it\zeta}=B_{3}/B_{1}$ when $Cu\ll 1$ . In (b), we show the swimming speed of a neutral squirmer (solid line) and a pusher/puller (dashed line) at values of ${\it\zeta}$ chosen below the curve in (a) so that the swimming speed is larger than the Newtonian value at small $Cu$ (upper solid line: ${\it\alpha}=0$ , ${\it\zeta}=-15$ ; upper dashed line: ${\it\alpha}=\pm 5$ , ${\it\zeta}=-4$ ); conversely, with values of ${\it\zeta}$ above the curve in (a), the swimmers always swim slower than in a Newtonian fluid, as shown (lower solid line: ${\it\alpha}=0$ , ${\it\zeta}=15$ ; lower dashed line: ${\it\alpha}=\pm 5$ , ${\it\zeta}=4$ ). Here, ${\it\varepsilon}=0.1$ , $n=0.25$ .

In figure 3(b), we show the variation of swimming speed for two swimmers with ${\it\alpha}$ and ${\it\zeta}$ chosen below the level set curve (the upper solid and dashed lines) and two swimmers with ${\it\alpha}$ and ${\it\zeta}$ chosen above the level set curve (the lower solid and dashed lines). We note that the swimming speeds of the faster swimmers in the small- $Cu$ regime experience a subsequent fall below the Newtonian value and then a rise above it as $Cu$ increases, before asymptoting to the Newtonian swimming speed at high $Cu$ . This indicates that micro-organisms (with a given swimming gait) can swim both faster and slower than in a Newtonian fluid depending on the actuation rate of that gait. In contrast, the two swimmers that swim slower in a Newtonian fluid in the small- $Cu$ regime remain slower for larger $Cu$ , with a non-monotonic variation similar to that observed previously (see figure 1). These results also hold qualitatively for large values of ${\it\epsilon}$ .

We emphasize that the thrust and drag reduction model is unable to explain the faster swimming speed with the addition of a $B_{3}$ mode, because in these cases the thrust still decreases more than the drag over a wide range of Carreau numbers, if they are considered separately. This serves as a counterexample demonstrating how analysing drag and thrust separately may not adequately describe swimming in complex media. This fact points specifically to the interaction between the thrust and drag fields, due to the nonlinearity in the constitutive equation, as the cause of faster than Newtonian swimming.