2.1. Soft microswimmer

In this study, swimming microorganisms such as ciliates are modelled as soft microswimmers driven by propulsive torques. The original model was developed by Ishikawa et al. (Reference Ishikawa, Tanaka, Imai, Omori and Matsunaga2016) and expanded to two-body interactions by Matsui et al. (Reference Matsui, Omori and Ishikawa2020b ). The model consists of a capsule with a two-dimensional hyperelastic thin membrane encasing an incompressible Newtonian fluid with the same viscosity $\mu$ and density $\rho$ as the ambient fluid (cf. figure 1 a). The reference shape of this model is a prolate with aspect ratio $ \textit{AR}$ , varied from a sphere of radius $a$ under constant volume (i.e. $4\pi a^3/3$ ). The thrust force at the cilium tip and the reaction force at the cilium root during the ciliary beating are modelled as the propulsive torque. Based on the assumption that the microbial ciliary beat occurs along the meridional direction, the propulsive torque is distributed as a vector oriented in the azimuthal direction over the torque surface $A_t$ offset by a distance of $\varepsilon$ from the membrane surface $A$ . The propulsive torque per unit area in the reference state $\boldsymbol{L}_r$ is given by

(2.1) \begin{align} \boldsymbol{L}_r = \left [L_0 + \kappa \left ( \theta - \frac {\pi }{2} \right )\right ]\boldsymbol{n}\times \boldsymbol{t}, \end{align}

where $\kappa$ determines the torque strength of the asymmetry, $\theta$ is the polar angle with respect to the orientation vector $\boldsymbol{e}$ defined by the posterior and anterior material points, $\boldsymbol{n}$ is the unit outward normal vector and $\boldsymbol{t}$ is the unit tangential vector in the $\theta$ direction. In our model, $|\boldsymbol{L}_r|$ is assumed to be time invariant and is determined by adjusting the value of $L_0$ so that the swimmer in a steady state has a swimming speed $U_0$ under the free swimming condition. The swimmer in the reference state satisfies the torque-free condition on the overall torque surface due to the axisymmetric torque distribution. In contrast, in the deformed and steady states, this torque-free condition can be approximately maintained by locally adjusting the torque strength in response to the deformation, as follows:

(2.2) \begin{align} \boldsymbol{L} = \frac {\textrm{d}A_{t,r}}{\textrm{d}A_t}\boldsymbol{L}_r , \end{align}

where $\textrm{d}A_{t,r}$ and $\textrm{d}A_t$ are the local torque surface area in the reference and deformed states, respectively. This correction means that the number of cilia does not change with deformation.

The stresslet due to the presence of the capsule swimmer can be expressed by applying the reciprocal theorem proposed by Pozrikidis (Reference Pozrikidis1992) to the membrane interface

(2.3) \begin{align} \unicode{x1D64E} = \int { \left [ \frac {1}{2}(\boldsymbol{q}\boldsymbol{x} + \boldsymbol{x}\boldsymbol{q}) - \frac {1}{3}\boldsymbol{q}\boldsymbol{\cdot }\boldsymbol{x}\unicode{x1D644} \right ]}\textrm{d}A, \end{align}

where $\boldsymbol{q}$ is a stress jump in the interfacial membrane with $\boldsymbol{n}$ at position $\boldsymbol{x}$ , given by using viscous stress $\boldsymbol{\mathsf{\sigma }}$ as

(2.4) \begin{align} \boldsymbol{q}(\boldsymbol{x}) = \left [ \boldsymbol{\mathsf{\sigma }}_{\textit{out}}(\boldsymbol{x}) - \boldsymbol{\mathsf{\sigma }}_{\textit{in}}(\boldsymbol{x}) \right ] \boldsymbol{\cdot }\boldsymbol{n}(\boldsymbol{x}). \end{align}

Our model employs the swimming mode $\beta$ of the classical squirmer model (Blake Reference Blake1971), such as pusher ( $\beta \lt 0$ ) and puller ( $\beta \gt 0$ ) types, to represent the differences in fluid motion of various swimming microorganisms. Here, $\beta$ is calculated from the stresslet component in the swimming direction $\unicode{x1D64E}_{\textit{ee}}$ (Matsui et al. Reference Matsui, Omori and Ishikawa2020b ; Kubo et al. Reference Kubo, Omori and Ishikawa2025)

(2.5) \begin{align} \beta =\frac {\unicode{x1D64E}_{\textit{ee}}}{4\pi \mu U_{0} a^2}. \end{align}

The value of $\beta$ in our model is indirectly determined by adjusting $\kappa$ , which creates the swimmer’s anterior–posterior asymmetry. The correspondences between $\kappa$ , $\beta$ and $\unicode{x1D64E}_{\textit{ee}}$ for the ellipsoidal swimmers with $ \textit{AR} = 2$ are shown in figure 2. In this study, $\beta$ was set to −1 for pushers and 1 for pullers.

Figure 2.

The relation between the asymmetry of torque strength $\kappa$ , the stresslet component in the swimming direction $\unicode{x1D64E}_{\textit{ee}}$ and $\beta$ of the soft swimmer with $ \textit{AR} = 2$ and each $ \textit{Ca}$ . Here, $ \textit{Ca}$ is the ratio between the viscous force due to swimming and the elastic force of the membrane defined by (2.22).

2.2. Membrane mechanics

In order to describe local deformation on the membrane, our model employs thin-shell theory of capsules developed by Barthès-Biesel, Pozrikidis, and their co-workers, based on continuum mechanics. Here, we briefly explain the outline of thin-shell theory. For further details, please refer to Barthès-Biesel et al. (Reference Barthès-Biesel, Diaz and Dhenin2002), Lac et al. (Reference Lac, Barthes-Biesel, Pelekasis and Tsamopoulos2004), Walter et al. (Reference Walter, Salsac, Barthes-Biesel and Tallec2010), Pozrikidis (Reference Pozrikidis2010).

To compute membrane tension, we first define the local curvilinear coordinate ( $\xi ^1, \xi ^2$ ) on the two-dimensional membrane. The covariant basis vectors $\boldsymbol{A}_{\alpha }$ ( $\alpha = 1, 2$ ) in the reference state are defined as

(2.6) \begin{align} \boldsymbol{A}_1 = \frac {\partial \boldsymbol{X}}{\partial \xi ^1}, \hspace {10.0mm} \boldsymbol{A}_2 = \frac {\partial \boldsymbol{X}}{\partial \xi ^2}, \end{align}

where $\boldsymbol{X}$ is the material position vector in the reference state. Similarly, the covariant basis vectors $\boldsymbol{a}_{\alpha }$ in the deformed state are given by

(2.7) \begin{align} \boldsymbol{a}_1 = \frac {\partial \boldsymbol{x}}{\partial \xi ^1}, \hspace {10.0mm} \boldsymbol{a}_2 = \frac {\partial \boldsymbol{x}}{\partial \xi ^2}. \end{align}

The associated contravariant base vectors $\boldsymbol{A}^{\alpha }$ are defined by $\boldsymbol{A}^\alpha \boldsymbol{\cdot }\boldsymbol{A}_\beta = \boldsymbol{\delta }^\alpha _\beta$ (likewise for $\boldsymbol{a}^{\alpha }$ ). Here, $\boldsymbol{\mathsf{\delta }}^\alpha _\beta$ is the Kronecker tensor. The covariant and contravariant metric tensors are then defined as

(2.8) \begin{align} \unicode{x1D63C}_{\alpha \beta } = \boldsymbol{A}_\alpha \boldsymbol{\cdot }\boldsymbol{A}_\beta , \hspace {10.0mm} \unicode{x1D63C}^{\alpha \beta } = \boldsymbol{A}^\alpha \boldsymbol{\cdot }\boldsymbol{A}^\beta . \end{align}

Applying the same rule to the deformed state, the covariant and contravariant metric tensors $\unicode{x1D656}_{\alpha \beta }$ and $\unicode{x1D656}^{\alpha \beta }$ in the deformed state are defined. The relation between the membrane material points in the reference and deformed states is described as follows:

(2.9) \begin{align} \boldsymbol{\textrm{d}x} = {\unicode{x1D641}}_s \boldsymbol{\cdot }\boldsymbol{\textrm{d}X}, \end{align}

where $\unicode{x1D641}_s$ is the surface deformation gradient tensor: $\unicode{x1D641}_s = \partial \boldsymbol{x}/\partial \boldsymbol{X}=\boldsymbol{a}_\alpha \otimes \boldsymbol{A}^\alpha$ . Local deformation of the membrane can be described by the Green–Lagrange strain tensor as

(2.10) \begin{align} \unicode{x1D640}_s &= \frac {1}{2} \left [\unicode{x1D641}_s^T \unicode{x1D641}_s - \left (\unicode {x1D644} - \boldsymbol{N}\boldsymbol{N}\right )\right ]\!, \end{align}

and covariant form as

(2.11) \begin{align} \unicode{x1D640}_{s,\alpha \beta } & = \frac {1}{2} \big(\unicode{x1D656}_{\alpha \beta } - \unicode{x1D63C}_{\alpha \beta } \big) \unicode{x1D63C}^\alpha \otimes \unicode{x1D63C}^\beta . \end{align}

The strain invariants are given by the principal extension ratios $\lambda _1$ and $\lambda _2$ as

(2.12) \begin{align} I_1 &=2tr\unicode{x1D640}_s \notag \\[3pt]&= \unicode{x1D656}_{\alpha \beta } \unicode{x1D63C}^{\alpha \beta } - 2\notag \\[3pt]&= \lambda _1^2 + \lambda _2^2 - 2, \end{align}

and

(2.13) \begin{align} I_2 &=det\unicode{x1D63E} - 1 \notag \\[3pt]&= |\unicode{x1D656}_{\alpha \beta }| \big|\unicode{x1D63C}^{\alpha \beta } \big| - 1 \notag \\[3pt]&= \lambda _1^2 \lambda _2^2 -1 \notag \\[3pt]&= J^2_s - 1, \end{align}

where $\unicode{x1D63E} = \unicode{x1D656}_{\alpha \beta }\unicode{x1D63C}^\alpha \otimes \unicode{x1D63C}^\beta$ is the right Cauchy–Green tensor and $J_s$ = $det$ ( $\unicode{x1D641}_s)=\lambda _1 \lambda _2$ is the Jacobian expressing the areal dilatation. Assuming the two-dimensional hyperelastic membrane, the membrane tension $\unicode{x1D64F}_s$ can be described by using a strain energy function $w_s(I_1, I_2)$ in the reference state, as

(2.14) \begin{align} \unicode{x1D64F}_s &= \frac {1}{J_s} \unicode{x1D641}_s \boldsymbol{\cdot }\frac {\partial w_s}{\partial \unicode{x1D640}_s} \boldsymbol{\cdot }\unicode{x1D641}^T_s. \end{align}

Using the chain rule, the expression for the contravariant representation of $\unicode{x1D64F}_s$ is obtained as

(2.15) \begin{align} \unicode{x1D64F}_s^{\alpha \beta }&= \left [\frac {2}{J_s} \frac {\partial w_s}{\partial I_1} \unicode{x1D642}^{\alpha \beta } + 2 J_s \frac {\partial w_s}{\partial I_2} \unicode{x1D656}^{\alpha \beta }\right ]\boldsymbol{a}_\alpha \otimes \boldsymbol{a}_\beta . \end{align}

This study has employed the constitutive law derived by Skalak et al. (Reference Skalak, Tozeren, Zarada and Chien1973) to express the behaviour of a biological membrane

(2.16) \begin{align} w_s = \frac {G_s}{4} \left ( I^2_1 + 2I_1 - 2I_2 + CI^2_2 \right ), \end{align}

where $G_s$ is surface shear elastic modulus and $C$ is a non-dimensional parameter determining the areal incompressibility. In this work, the value of $C$ was set to ten with reference to the erythrocyte model (Omori et al. Reference Omori, Ishikawa, Imai and Yamaguchi2013), which characterises strain hardening behaviour and the incompressible areal dilatation of the membrane. The Poisson ratio of the Skalak law is $\nu _s = C/(1+C)$ .

The bending energy of the membrane $w_b$ is calculated based on a model proposed by Helfrich (Reference Helfrich1973)

(2.17) \begin{align} w_b = \frac {E_b}{2}\int \left ( 2H-c_0 \right )\textrm{d}A, \end{align}

where $H$ is a mean curvature of the membrane, $c_0$ is the spontaneous curvature and $E_b$ is the bending modulus. By assuming a symmetric membrane property, (i.e. $c_0 = 0$ ), the bending energy can be converted to the load due to a bending $\boldsymbol{q}_b$ as (Matsui et al. Reference Matsui, Omori and Ishikawa2020a )

(2.18) \begin{align} \boldsymbol{q}_b = -2E_b\big [\varDelta _sH + 2H(H^2 - K)\big ]\boldsymbol{n}, \end{align}

where $\varDelta _s$ is the Laplace–Beltrami operator and $K$ is the Gaussian curvature on the membrane surface.

2.3. Fluid mechanics

Considering the typical small size of swimming microorganisms, the surrounding flow field can be regarded as a Stokes flow. Thus, inertial effects on the internal and external flow fields can be neglected. The fluid velocity $\boldsymbol{v}$ at an arbitrary position $\boldsymbol{y}$ produced by $N$ swimmers in an ambient flow $\boldsymbol{v}_{\infty }$ can be written as a boundary integral on the deformed membrane $A$ and torque surface $A_t$ (Pozrikidis Reference Pozrikidis1992)

(2.19) \begin{align} \boldsymbol{v}(\boldsymbol{y}) = \boldsymbol{v}_{\infty }(\boldsymbol{y}) - \frac {1}{8\pi \mu }\sum _i^N \left [ \int _{A} \unicode{x1D64D}(\boldsymbol{x},\boldsymbol{y}) \boldsymbol{\cdot }\boldsymbol{q}(\boldsymbol{x}) \textrm{d}A^i + \int _{A_{t}} \unicode{x1D64D}(\boldsymbol{x},\boldsymbol{y}) \boldsymbol{\cdot }\boldsymbol{L}(\boldsymbol{x}) \textrm{d}A^i_t \right ]\!, \end{align}

where $\mu$ is viscosity inside and outside the membrane. The ambient flow is imposed as a shear flow in the $x$ – $z$ plane with a shear rate $\dot {\gamma }$ , given by $\boldsymbol{v}_{\infty } = (\dot {\gamma }z, 0, 0)$ , in this work. The Green’s functions $\unicode{x1D64D}$ and $\unicode{x1D64D}$ correspond to the stokeslet and rotlet in a triply periodic lattice, respectively (Beenakker Reference Beenakker1986).

2.4. Numerical method

This section provides briefly discretisation of the soft swimmer and the numerical method combining membrane mechanics and fluid mechanics in the Stokes flow. The swimmer’s membrane and torque surface are respectively discretised by 1280 triangular elements with 642 node points, as shown in figure 1(a). Local membrane deformations are tracked over time in the Lagrangian frame using these node points over time.

In order to solve the fluid–structure interaction between the membrane and surrounding fluid in the Stokes regime, this study has employed a numerical method similar to the previous work (Walter et al. Reference Walter, Salsac, Barthes-Biesel and Tallec2010). For the membrane mechanics, there is an equilibrium relation between the load due to a membrane in-plane stretch $\boldsymbol{q}_s$ and membrane tension, which can be written in weak form based on the principle of virtual work as

(2.20) \begin{align} \int _A \hat {\boldsymbol{u}}\boldsymbol{\cdot }\boldsymbol{q}_s\textrm{d}A = \int _A \hat {\boldsymbol{\mathsf{\varepsilon }}} : \unicode{x1D64F} \textrm{d}A, \end{align}

where $\hat {\boldsymbol{u}}$ and $\hat {\boldsymbol{\mathsf{\varepsilon }}}=({1}/{2})(\boldsymbol{\nabla} _s \hat {\boldsymbol{u}} + \boldsymbol{\nabla} _s \hat {\boldsymbol{u}}^T)$ are a virtual displacement and strain, respectively. This equilibrium equation is solved numerically by using a finite-element method to obtain $\boldsymbol{q}_s$ from membrane tension defined geometrically according to the Skalak constitutive law. The viscous load $\boldsymbol{q}$ on the membrane is given by a superposition of loads due to the stretch and bending (i.e. $\boldsymbol{q} = \boldsymbol{q}_s + \boldsymbol{q}_b$ ). For the fluid mechanics of the surrounding flow, substituting $\boldsymbol{q}$ into the boundary integral equation (2.19) results in a new velocity field $\boldsymbol{v}$ , which was calculated by a boundary-element method. This study has used a Gaussian quadrature method with six numerical integration points in each element to compute $\boldsymbol{v}$ . A material point on the torque surface is given by $\boldsymbol{x}_t = \boldsymbol{x}+\varepsilon \boldsymbol{n}(\boldsymbol{x})$ , and the flow due to the torque $\boldsymbol{L}$ is calculated in the same integral manner. The material point on the torque surface is defined as the point located a distance $\varepsilon$ away from the membrane in the normal direction and the flow induced by the torque $L$ is calculated using the same integration method. Here, $\varepsilon$ was set to $0.05a$ , which corresponds to the half-length of the cilium. Assuming a no-slip boundary at the membrane, $\textrm{d}\boldsymbol{x}/\textrm{d}t = \boldsymbol{v}$ , the material point $\boldsymbol{x}$ is updated by a second-order Runge–Kutta method. The above cycles are iteratively calculated to simulate hydrodynamic interactions in the suspension.

In order to accelerate the computational process, a multipole expansion technique (Pozrikidis Reference Pozrikidis1992) has been employed, classifying interactions into short-, middle- and long-range regimes depending on the distance between swimmers’ centres (Matsunaga et al. Reference Matsunaga, Imai, Yamaguchi and Ishikawa2016). Further details and settings are described in the Supplementary material. Mesh overlaps between swimmer surfaces can lead to numerical instabilities when the swimmers come into close proximity. Thus, pairwise repulsive force is introduced as an exponential function (Brady & Bossis Reference Brady and Bossis1985; Ishikawa & Pedley Reference Ishikawa and Pedley2007), defined as

(2.21) \begin{align} \boldsymbol{F}_{rep}(\boldsymbol{r}_t) = \begin{cases} \alpha _1 \frac {\alpha _2 \exp (-\alpha _2r_t)}{1 - \exp (-\alpha _2r_t)}\frac {\boldsymbol{r}_t}{|\boldsymbol{r}_t|} & (|\boldsymbol{r}_t| \leq a) \\[3pt] \boldsymbol{0} &(|\boldsymbol{r}_t|\gt a) \end{cases}, \end{align}

where $\boldsymbol{r}_t$ is the distance between node points on the torque surface of each swimmer. The parameters related to the strength and decay rate of the repulsive force have been tuned as $\alpha _1 = 20$ and $\alpha _2 = 25$ .

In this study, two non-dimensional parameters were introduced to characterise the system dynamics. The first is the capillary number, defined as

(2.22) \begin{align} \textit{Ca} = \frac {\mu U_0}{G_s} . \end{align}

which quantifies membrane deformability as the ratio of viscous forces due to swimming with $U_0$ to the elastic restoring force of the membrane. The second parameter, $\dot {\gamma }a/U_0$ , characterises the relative strength of the imposed shear flow compared with the swimmer’s self-propulsion. All computations were normalised by using the swimmer’s reference radius $a$ , swimming speed $U_0$ and surface shear elastic modulus $G_s$ .

2.5. Rheological properties

The particle bulk stress tensor can be obtained as the average of the stresslet over $N$ particles in the unit domain

(2.23) \begin{align} \boldsymbol{\mathsf{\varSigma }}^{(p)} = \frac {1}{V} \sum _{i=1}^{N} \unicode{x1D64E}_i = \frac {3c}{4\pi a^3} \left \langle \unicode{x1D64E} \right \rangle , \end{align}

where $V$ is the volume of the unit fluid domain and $\left \langle \right \rangle$ means the ensemble average. This study focuses on the contributions to viscosity and normal stress differences as fundamental rheological properties of the suspension. The particle contribution to viscosity is given by

(2.24) \begin{align} \frac {\mu _a - \mu }{\mu } = \frac { \mathsf{\varSigma }^{(p)}_{xz}}{\dot {\gamma }\mu } = \frac {3c}{4\pi a^3\dot {\gamma }\mu } \left \langle \unicode{x1D64E}_{xz} \right \rangle . \end{align}

Here, $\mu _a$ is the apparent viscosity and $(\mu _a - \mu ) / \mu$ is the specific viscosity. The dimensionless first and second normal stress differences are defined, in a similar manner, as

(2.25) \begin{align} \frac {\mathsf{\varSigma }^{(p)}_{xx} - \mathsf{\varSigma }^{(p)}_{zz}}{\dot {\gamma }\mu } = \frac {3c}{4\pi a^3\dot {\gamma }\mu } \left \langle \unicode{x1D64E}_{xx} - \unicode{x1D64E}_{zz} \right \rangle \end{align}

and

(2.26) \begin{align} \frac {\mathsf{\varSigma }^{(p)}_{zz} - \mathsf{\varSigma }^{(p)}_{yy}}{\dot {\gamma }\mu } = \frac {3c}{4\pi a^3\dot {\gamma }\mu } \left \langle \unicode{x1D64E}_{zz} - \unicode{x1D64E}_{yy} \right \rangle . \end{align}

3.1. Problem setting

The ellipsoidal swimmer with $ \textit{AR} = 2$ was initially aligned with a polar angle of $\varTheta = \pi /4$ and azimuthal angle of $\varPhi = \pi /4$ , as shown in figure 3. The capillary number $ \textit{Ca}$ was fixed at 0.1 while the imposed shear rate was varied in each condition.

Figure 3.

Problem setting of a simulation using the solitary swimmer model. The ellipsoidal swimmer with $ \textit{AR} = 2$ is initially oriented towards $\varPhi = \varTheta = \pi /4$ . A shear flow is set in the $x$ – $z$ plane as $\boldsymbol{v}_\infty = (\dot {\gamma }z, 0, 0)$ . A red (blue) marker shows the anterior (lateral) node point on the membrane.

3.2. Second normal stress difference

A characteristic of the swimming behaviour of elliptical capsule swimmers in shear flows is the reorientation of their swimming direction, which differs from the Jeffery orbit. For example, swimming direction, which was initially oriented at $\varPhi = \pi /4$ , converges towards the vorticity direction $(\varPhi = 0)$ whilst undergoing rotation due to shear. Once convergence has occurred in the vorticity direction, the swimmer achieves steady motion even under shear flow. The reorientation of the elliptical capsule swimmer towards the vorticity direction produces significant differences between the pusher and puller with respect to the second normal stress difference. Thus, in this section, we mainly discuss the second normal stress difference of the suspension. For further details on the results regarding reorientation, shear viscosity and first normal stress difference are provided in the Supplemental Material.

We evaluated the second normal stress difference of the dilute suspension from the stresslet acting on the solitary swimmer during steady rotational motion. In the dilute limit, the suspension rheology is often discussed by evaluating the coefficient of the first-order term in the polynomial expansion with respect to $c$ : the value of (2.26) divided by $c$ . The result is shown in figure 4, where all data were averaged over the time for one rotation. In order to compare with the former results of ellipsoidal passive capsules, instead of microswimmers, we introduce another capillary number, $ \textit{Ca}^*$ , defined using the viscous force due to the shear flow (Zhang & Graham Reference Zhang and Graham2020)

(3.1) \begin{align} \textit{Ca}^* = \frac {\mu \dot {\gamma }a}{G_s}. \end{align}

Due to the deformability, the swimmer tends to reorient in its rotational motion. As a result, the second normal stress difference shows a different tendency depending on the swimming mode, especially in the low shear rate regime. At low shear rates, the swimming directions of both the pusher and the puller are aligned with the vorticity direction. Even when the swimming directions are the same, the stresslet in the orientational and lateral directions is opposite for pushers and pullers. Thus, the second normal stress difference yields markedly different results for the pusher and puller. Conversely, at the high shear rates ( $ \textit{Ca}^*\gt 0.5$ ), this stress difference asymptotically approaches that of the passive capsule regardless of the swimming mode because the swimmer deforms mainly due to the background flow.

Figure 4.

Second normal stress differences of the dilute suspensions estimated by the solitary swimmer with $ \textit{AR} = 2$ and $ \textit{Ca} = 0.1$ . The values are averaged over one rotational period. Each marker shows the difference in the rotational mode in the steady motion. For comparison, a blue dashed lines indicate the results for a passive prolate capsule with $ \textit{AR} = 2$ (Zhang & Graham Reference Zhang and Graham2020). Here, $ \textit{Ca}^*$ is the capillary number based on the shear flow defined as $ \textit{Ca}^* = \mu \dot {\gamma }/G_s$ .

4.1. Problem setting

We applied triply periodic boundary conditions proposed by Lees & Edwards (Reference Lees and Edwards1972) to the computational cubic lattice, to express hydrodynamic interactions in the infinite domain. Initially, twenty-seven swimmers were randomly assigned positions and orientations within the domain, as shown in figure 5(a). By adjusting the size of the domain, the desired volume fraction $c$ can be controlled. The shear flow was set in the $x$ – $z$ plane as $\boldsymbol{v}_\infty = (\dot {\gamma }z, 0, 0)$ (cf. figure 5 b).

Figure 5.

Problem setting of a simulation using multiple swimmer models. (a) Twenty-seven swimmers are suspended in the computational unit, and the triply periodic boundary condition is applied to the domain. The shear flow is set in the $x$ – $z$ plane as $\boldsymbol{v}_\infty = (\dot {\gamma }z, 0, 0)$ . (b) Instantaneous snapshot of 27 pusher swimmers in the sheared suspension ( $c = 0.1$ ). Red and blue square markers show trajectories of the swimmers.

Figure 6.

(a–b) Time evolution of the correlation function between the swimmer’s orientation and vorticity axis of the shear flow in the suspension. (c–d) Time evolution of the correlation function between the projected swimmer’s orientation in the shear plane and elongational axis of the shear flow. The swimmers with $ \textit{AR} = 2$ and $ \textit{Ca} = 0.1$ are exposed to the shear flow with $\dot {\gamma }a/U_0 = 0.05$ .

Figure 6. Long description

The image contains four line graphs labeled (a), (b), (c), and (d), each showing the time evolution of correlation functions in a suspension of swimmers under shear flow. Panel A: A line graph titled Pusher, c = 0.1 shows the time evolution of the correlation function between the swimmers’ orientation and the vorticity axis of the shear flow. The x-axis is labeled tU0/a and ranges from 0 to 1000. The y-axis is labeled Iv and ranges from -0.5 to 1.0. Two lines are plotted: one in red for initial polar orientation and one in blue for initial random orientation. The red line starts at 1.0 and decreases rapidly, while the blue line starts near 0 and fluctuates around it. Panel B: A line graph titled Puller, c = 0.05 shows the time evolution of the same correlation function as in Panel A. The x-axis and y-axis labels and ranges are identical to Panel A. The red line starts at 1.0 and decreases rapidly, while the blue line starts near 0 and fluctuates around it. Panel C: A line graph titled Pusher, c = 0.1 shows the time evolution of the correlation function between the projected swimmers’ orientation in the shear plane and the elongational axis of the shear flow. The x-axis is labeled tU0/a and ranges from 0 to 1000. The y-axis is labeled Ie and ranges from -0.5 to 1.0. Two lines are plotted: one in red for initial polar orientation and one in blue for initial random orientation. Both lines fluctuate around 0. Panel D: A line graph titled Puller, c = 0.05 shows the time evolution of the same correlation function as in Panel C. The x-axis and y-axis labels and ranges are identical to Panel C. Both lines fluctuate around 0.

4.2. Distribution of swimmers’ orientations

First, we investigate how the orientation of microswimmers changes in a more concentrated suspension, where interactions between the swimmers cannot be ignored. We define the orientational correlation function $I_v$ between the swimmer’s orientation and the vorticity axis as follows:

(4.1) \begin{align} I_v = \frac {3}{2}\left \langle (\boldsymbol{e}_i\boldsymbol{\cdot }\boldsymbol{n}_{v})^2 - \frac {1}{3} \right \rangle _N, \end{align}

where $\boldsymbol{n}_{v}$ is the vorticity axis: $\boldsymbol{n}_{v} = (0, 1, 0)$ . Specifically, $I_v = 1$ corresponds to the swimmer being parallel to $\boldsymbol{n}_{v}$ , $I_v = -0.5$ corresponds to a perpendicular orientation and $I_v = 0$ indicates an isotropic distribution with respect to $\boldsymbol{n}_{v}$ . Typical results of ellipsoidal pushers and pullers with $ \textit{AR} = 2$ are shown in figures 6(a) and 6(b), respectively. The shear strength is set as $\dot {\gamma }a/U_0 = 0.05$ and the capillary number is 0.1 for both cases. The correlation function reached a developed state when $tU_0/a \gt 250$ for both swimmer types. In the case of the dilute limit without the cell–cell interactions, the swimmer in low shear rate tends to align to the vorticity axis of the shear flow (cf. Supplementary Material). However, when the interparticle interaction is considered, the time-averaged value in the developed state converges to nearly zero. This result suggests that reorientation towards the vorticity direction is being eliminated by fluid interactions between swimmers. Furthermore, to verify initial condition dependence, the simulations were initiated from a polar-order state where the initial swimming direction was aligned with the vorticity axis. Even in this case, it was confirmed that the value of the correlation function finally approaches zero. These results suggest that hydrodynamic interactions prevent the reorientation of microswimmers along the vorticity axis.

The question arises as to whether the swimmer is swimming isotropically within three-dimensional space. We subsequently define an additional correlation function, $I_e$ , in a manner analogous to (4.1). By changing the reference axis from $\boldsymbol{n}_{v}$ to the extension direction vector $\boldsymbol{n}_{e} = (1/\sqrt {2}, 0, 1/\sqrt {2})$ , we obtain a quantity $I_e$ that represents the correlation between the projected swimmer’s orientation and elongational axis of the shear flow. As shown in figures 6(c) and 6(d), the value of $I_e$ is positive and no longer equal to zero. This suggests that the mean swimming directions manifest an asymmetric distribution, exhibiting a bias towards the direction of shear elongation.

In order to gain deeper insight into these results, we evaluate the averaged orientational distribution of the swimmers using the generalised orientation distribution function $P$ , defined as

(4.2) \begin{align} P = \frac {\left \langle N_l(\varPhi , \varTheta ) \right \rangle _{t}}{\Delta S(\Delta \varPhi , \Delta \varTheta ) n_{mean}}, \end{align}

where $N_l$ is defined by the number of swimmers oriented within the surface element $\Delta S$ on the unit sphere defined by angular intervals of $\Delta \varPhi$ and $\Delta \varTheta$ , and $n_{mean} = N/4\pi$ is the global surface number density. Figure 7(a–d) shows contour plots of $P$ mapped onto the spherical surface. An anisotropic distribution concentrated in the elongational direction is observed at low shear rates $\dot {\gamma }a/U_0 = 0.05$ . As the shear strength increases, the orientational distribution becomes biased towards the streamwise direction. At very high levels of shear, the distribution concentrates on the shear plane. Such shear plane concentrations can be attributed to the predominance of external flow, which, in conjunction with the vorticity inherent in the shear flow, results in an augmentation of swimmer’s rotation. These results indicate that the orientational distribution in the low shear regime becomes more asymmetric with the bias in the elongational quadrant. Changes in the swimming direction due to fluid interactions exert an effective randomness analogous to thermal fluctuations in Saintillan (Reference Saintillan2010). In the rotation of rod-like objects within a shear flow field, orientational fluctuations lead to symmetry breaking between the times spent in the extensional and compressional quadrants. This effect becomes dominant in the low shear rate regime because the time scale required for alignment due to the interactions is likely shorter than that for flow-induced reorientation towards the vorticity axis. Although hydrodynamic interaction also induces local swimmers’ aggregation and weak nematic order in the near field, these effects are not strong enough to form a global coherent polar order.

Figure 7.

(a–d) Orientational distribution for pusher suspensions with $ \textit{AR} = 2$ , $ \textit{Ca} = 0.1$ , $c = 0.1$ and various $\dot {\gamma }a/U_0$ . Here, $P(\varTheta ,\varPhi )$ is defined as the time-averaged probability density of the orientation over one rotational period, which has been normalised by the homogeneous one.

4.3. Shear viscosity and normal stress differences

Based on the results shown in figure 6, this study evaluated the rheological properties of the suspension of the soft microswimmers as time-averaged values over one rotational period after the microstructure had reached a sufficiently developed steady state.

Results of the pusher suspensions with $ \textit{Ca} = 0.1$ are shown in figure 8. Similar to the former experiments (Rafaï et al. Reference Rafaï, Jibuti and Peyla2010; López et al. Reference López, Gachelin, Douarche, Auradou and Clément2015), negative specific viscosity was measured when $\dot {\gamma }a/U_0 \sim 0.1$ and the value decreases with higher $c$ (see figure 8 a). In contrast, for high $\dot {\gamma }a/U_0$ , the results show the shear-thinning characteristic with positive values. Regarding the normal stress differences of the pusher suspension, the response to shear differs qualitatively between the first and second normal stress differences (cf. figures 8 b and 8 c). The first normal stress differences in the pusher suspension become negative in low shear rate regimes, because the pusher swimmer orienting closer to the streamwise direction generates a negative active stress $\mathsf{\varSigma }^{(p)}_{xx}$ . In the transverse direction, on the other hand, it tends to generate a positive active stress, meaning that the first normal stress difference tends to be negative. In the vorticity axis direction, the orientational distribution is isotropic. In the shear plane, however, the swimming direction is slightly offset from the streamwise direction towards the extensional direction. Consequently, a difference arises between $\mathsf{\varSigma }^{(p)}_{zz}$ and $\mathsf{\varSigma }^{(p)}_{yy}$ , generating a positive second normal stress difference. By contrast, at high shear rates, shear-induced deformation becomes more dominant than that of propulsive torque, as the swimmers are stretched primarily in the $x$ -direction and compressed in the $y$ - and $z$ -directions by the shear flow. Therefore, under high shear, the first normal stress difference becomes positive while the second normal stress difference asymptotically approaches zero.

Figure 8.

Rheological properties of suspensions of pusher swimmers with $ \textit{AR} = 2$ and $ \textit{Ca} = 0.1$ . (a) Specific viscosity. The contributions to (b) 1st and (c) 2nd normal stress differences. Error bars mean standard deviations for the three independent cases discussed throughout the paper.

Figure 9.

Rheological properties of suspensions of puller swimmers with $ \textit{AR} = 2$ and $ \textit{Ca} = 0.1$ . (a) Specific viscosity. The contributions to (b) 1st and (c) 2nd normal stress differences.

Figure 9. Long description

Three line graphs depict rheological properties of puller swimmers with varying parameters. Panel A: A line graph shows the specific viscosity of puller swimmers. The x-axis represents the product of gamma, a, and U0 on a logarithmic scale, while the y-axis represents the specific viscosity normalized by mu. Two data series are plotted: one with blue triangles for c = 0.05 and another with green diamonds for c = 0.01. The blue series shows a decreasing trend, while the green series remains relatively constant. Panel B: A line graph illustrates the contributions to the first normal stress difference. The x-axis is the same as in Panel A, and the y-axis represents the normalized first normal stress difference. The blue series shows a decreasing trend, while the green series remains relatively constant. Panel C: A line graph depicts the contributions to the second normal stress difference. The x-axis is the same as in the previous panels, and the y-axis represents the normalized second normal stress difference. Both the blue and green series show slight increases. A small diagram in the bottom right corner of Panel B shows a coordinate system with axes labeled x, y, and z, and an arrow indicating the direction of v_infinity.

Next, we investigate the rheology of the puller suspensions. The capillary number is again set as $ \textit{Ca} = 0.1$ and the result is shown in figure 9. Basically, the orientational distributions of pullers are similar to those of the pusher suspensions. However, the sign of the active stress generated by the puller is opposite to that of the pusher. Therefore, at low shear rates, the puller exhibits a high viscosity increase (cf. figure 9 a). This effect is observed at a volume fraction of around 0.05 and is barely noticeable at a volume fraction of 0.01, where the frequency of interactions also decreases. Similarly, the results for normal stress differences are the opposite of those for pushers. At low shear rates, the first normal stress difference is positive, while the second is negative.

4.4. Specific viscosity in dilute suspensions with or without hydrodynamic interactions

According to the scaling by Batchelor & Green (Reference Batchelor and Green1972), the frequency of $n$ -body interaction is proportional to $c^n$ . For $c = 0.01$ , the contribution of two-body interactions is expected to be approximately 1 % of that of a solitary swimmer. Thus, the rheological properties of the suspension with $c = 0.01$ can be expected to be well described by the leading-order term of $O(c)$ in this scaling. However, the comparison between the multi-body system and solitary swimmer reveals distinct differences arising from alignment of the swimmers at the low shear rate (cf. figure 10). When the hydrodynamic interactions between swimmers are considered, the viscosity contribution of the pusher decreases relative to the results of the solitary swimmer, and vice versa for the puller. These results demonstrate that interactions between swimmers play a decisive role in disrupting stable orbital motion in shear flow, even at a volume fraction of 1 % and with a low frequency of hydrodynamic interactions. In other words, the ideal approximation of the active suspension rheology based on the solitary motion might be untenable, even in low volume fractions. Indeed, the previous experiments (Rafaï et al. Reference Rafaï, Jibuti and Peyla2010; López et al. Reference López, Gachelin, Douarche, Auradou and Clément2015) reported inherent viscosity change at the low volume fraction with the range of $O(10^{-3})$ . It is expected that hydrodynamic interactions will be effective because the presence of flagella enables cells to interact over a wider range.

Figure 10.

Comparison of the rheological properties in the dilute regime obtained by multiple swimmer models with $c = 0.01$ , with those estimated by the solitary swimmer model. (a) Viscosity contribution. The (c) 1st and (d) 2nd normal stress differences. The aspect ratio of the swimmer is set to 2 and capillary number is set to 0.1.

4.5. Comparison with hydrodynamic interactions and Brownian rotational diffusion

From the discussion so far, we found that the interparticle hydrodynamic interaction induces effective orientational fluctuation and results in a biased distribution towards the elongational direction of the shear flow. Similar distributions have been reported by Saintillan (Reference Saintillan2010) using a time-dependent Fokker–Planck equation for orientation. In this model, orientational fluctuations are modelled as Brownian diffusion. Here, clarifying the differences between fluctuations arising from hydrodynamic interactions and those due to Brownian rotational diffusion may be worthwhile. In order to consider the effects of Brownian motion on the suspension rheology, we added Gaussian noise to the rotational motion of the soft swimmer. Assuming isotropic diffusivity, the rotational displacement due to the Brownian motion $\Delta \omega _b$ follows a Gaussian distribution with $\left \langle \Delta \omega _b(\Delta t)\right \rangle = 0$ and is defined by $D_R$ as

(4.3) \begin{align} \Delta \omega _b(\Delta t) = \sqrt {2 D_R \Delta t}. \end{align}

Here, the non-dimensional rotational diffusivity is set as $D_{R}^*=D_{R} a/U_0 = 0.004$ as estimated from the simulation results of a non-Brownian soft swimmer for a volume fraction of $c=0.1$ , without background flow ( $\dot {\gamma }a/U_0 = 0$ ), a capillary number of $ \textit{Ca}=0.1$ and an aspect ratio of $ \textit{AR}=2$ . The results of the Brownian soft swimmer are shown in figure 11. Similar to hydrodynamic interactions, a biased orientational distribution towards elongation was observed even in the solitary Brownian soft pusher at the low shear rate $\dot {\gamma }a/U_0 = 0.05$ . However, the higher peak value of $P$ on the line of $\varPhi = 0.5\pi$ for the Brownian swimmer is observed, compared with the result considering only the fluctuation induced by the interparticle hydrodynamic interaction. This quantitative discrepancy is likely attributable to the emergence of a weak alignment effect arising from the hydrodynamic interaction. As a result, the degree of specific viscosity becomes higher (cf. figure 11 c). These results imply that the effects of hydrodynamic interactions on suspension rheology cannot be fully replaced by those of Brownian noise.

Figure 11.

(a) Orientational distribution estimated by the solitary Brownian soft pusher. (b) Distribution of $P(\varTheta )$ in the shear plane ( $\varPhi =0.5\pi$ ). The red line indicates the results of non-Brownian soft pushers with hydrodynamic interactions, whilst the green line represents the results of a solitary Brownian soft pusher. (c) Specific viscosity estimated from a solitary Brownian soft pusher and non-Brownian soft pushers with hydrodynamic interactions. The results for Brownian swimmer assume a dilute limit ( $\mu _a \sim \mu - 7.16c$ ) and, for comparison, have been estimated using a volume fraction of $c = 0.1$ .

Figure 11. Long description

The image contains three panels with different types of graphs. Panel A shows a 3D surface plot representing the orientational distribution estimated by the solitary Brownian soft pusher. The axes are labeled x, y, and z, with a color gradient indicating the probability P, ranging from 0.5 to 2.0. Panel B is a line graph displaying the distribution in the shear plane. The x-axis is labeled with the variable Theta/pi, ranging from 0 to 1.0, and the y-axis is labeled P(Phi = 0.5pi), ranging from 0 to 4. The graph includes two lines: a red line for non-Brownian soft pushers with hydrodynamic interactions and a green line for a solitary Brownian soft pusher. The graph is divided into two quadrants: the Elongational quadrant (shaded in blue) and the Compressional quadrant (shaded in orange). Panel C is a scatter plot showing the specific viscosity estimated from a solitary Brownian soft pusher and non-Brownian soft pushers with hydrodynamic interactions. The x-axis labels two categories: Solitary soft Brownian swimmer and Suspension with hydrodynamic interactions. The y-axis is labeled (mu_a - mu)/mu, ranging from -1.0 to 0. The plot includes two data points with error bars representing the variability in the measurements.

4.6. Effects of the capillary number

In this section, we investigate the effects of $ \textit{Ca}$ on specific viscosity and normal stress differences. The shear rate is set as $\dot {\gamma }a/U_0 = 0.05$ for all cases. The change in $ \textit{Ca}$ under a fixed value of $a\dot \gamma /U_0$ is similar to the change in membrane stiffness.

The orientational distributions of the swimmers for various $ \textit{Ca}$ are shown in figure 12. Due to numerical instability arising from excessive membrane deformation, stable simulations in this study were limited to pusher suspensions with $ \textit{Ca} \geq 0.1$ and puller suspensions with $ \textit{Ca} \leq 0.1$ , respectively. In the pusher suspensions, the increase in $ \textit{Ca}$ resulted in the decrease in the peak value of $P$ . Consequently, as $ \textit{Ca}$ increases, the effect of active stress approaches zero (see figure 13). On the other hand, $ \textit{Ca}$ has little influence on the normal stress differences. For pullers, the peak value of $P$ does not change with $ \textit{Ca}$ , but $\varTheta$ at which $P$ is the maximum shifts slightly closer to the $x$ -axis ( $\varTheta =0.5\pi$ ) as $ \textit{Ca}$ increases. As the bias in the direction of elongation is reduced, the effect of active stress diminishes with increasing $ \textit{Ca}$ , much like with the pusher.

These results suggest that membrane deformability influences primarily viscosity through structural changes in the orientational distribution of swimmers. In particular, we infer that a softer membrane reduces the contribution of active stress to shear viscosity.

Figure 12.

Effects of $ \textit{Ca}$ on the orientational distribution of the swimmers. The left (right) graphs indicate results of pusher suspensions with $c = 0.1$ (puller suspensions with $c = 0.05$ ) and $\dot {\gamma }a/U_0 = 0.05$ . The aspect ratio of the swimmer is set to 2: (a) $ \textit{Ca} = 0.15$ and (c) $ \textit{Ca} = 0.2$ of pusher swimmers; (b) $ \textit{Ca} = 0.03$ and (d) $ \textit{Ca} = 0.1$ of puller swimmers. (e–f) Orientation distribution functions on the line with $\varPhi = 0.5\pi$ . The blue (red) area indicates that the swimmer spends in the extensional (compressional) quadrant.

Figure 13.

Effects of $ \textit{Ca}$ on the rheological properties. (a–b) Specific viscosity. (c–d) Contribution to 1st normal stress differences. (e–f) Contribution to 2nd normal stress differences. The left (right) graphs indicate results of pusher suspensions with $c = 0.1$ (puller suspensions with $c = 0.05$ ) and $\dot {\gamma }a/U_0 = 0.05$ . The aspect ratio of the swimmer is set to 2.

Figure 13. Long description

The image contains six line graphs, each illustrating different rheological properties for pusher and puller suspensions. Panel A and B show specific viscosity, Panel C and D show the contribution to the first normal stress differences, and Panel E and F show the contribution to the second normal stress differences. Each graph plots a variable against the capillary number (Ca). Panel A, C, and E represent pusher suspensions with a concentration of 0.1, while Panel B, D, and F represent puller suspensions with a concentration of 0.05. The aspect ratio of the swimmer is set to 2. Panel A shows an increasing trend in specific viscosity with increasing Ca. Panel B shows a decreasing trend in specific viscosity with increasing Ca. Panel C shows a fluctuating trend in the contribution to the first normal stress differences with increasing Ca. Panel D shows a relatively stable trend in the contribution to the first normal stress differences with increasing Ca. Panel E shows a fluctuating trend in the contribution to the second normal stress differences with increasing Ca. Panel F shows a decreasing trend in the contribution to the second normal stress differences with increasing Ca.

4.7. Effects of the aspect ratio

Finally, we investigate the effects of the aspect ratio $ \textit{AR}$ on the suspension rheology. The capillary number and shear rate are set to 0.1 and 0.05, respectively. For $ \textit{AR} = 1$ , the resultant orientation distribution is almost isotropic for both pushers and pullers, as shown in figure 14. Generally, the soft spherical swimmer rotates at an almost constant angular velocity during the one-periodic rotation (Matsui et al. Reference Matsui, Omori and Ishikawa2020b ). By contrast, ellipsoidal swimmers, like Jeffery orbits, have a non-constant rotational speed and tend to spend more time aligned with the streamwise direction. As the swimmer’s aspect ratio increases, a bias towards an orientation in the direction of shear elongation becomes observable (see figures 14 b and 14 c), due to uneven rotational velocity and fluid interactions between particles.

Changes in the aspect ratio also strongly influence rheology due to differences in orientation distribution (cf. figure 15). In pusher suspensions, specific viscosity changes from positive to negative values as the aspect ratio increases. In pullers, however, viscosity increases with the aspect ratio. For the first normal stress difference, the negative value increases with $ \textit{AR}$ for the pusher, but increases monotonically for the puller. The second normal stress difference is similarly affected by the aspect ratio, with values increasing in the pusher and decreasing in the puller. These results suggest that the swimmer’s body shape is one of the crucial factors causing the asymmetrical orientational distribution, resulting in the appearance of the active rheology of the microswimmers.

Figure 14.

Effects of the swimmer’s shape on the orientational distribution for the pusher suspensions with $c = 0.1$ and $\dot {\gamma }a/U_0 = 0.05$ . The capillary number is set to 0.1: (a) $ \textit{AR} = 1$ , (b) $ \textit{AR} = 1.5$ . (c) Distribution of $P(\varTheta )$ in the shear plane ( $\varPhi =0.5\pi$ ).

Figure 15.

Effects of swimmers’ shapes on the rheological properties. The capillary number is set to 0.1 in all swimmer and shear rate is set to 0.05. (a,b) Specific viscosity. (c,d) Contribution to 1st normal stress differences. (e,f) Contribution to 2nd normal stress differences. The left (right) graphs indicate results of pusher suspensions with $c = 0.1$ (puller suspensions with $c = 0.05$ ).

Figure 15. Long description

The figure contains six line graphs comparing the rheological properties of pusher and puller suspensions. Each graph plots data against the aspect ratio (AR) on the horizontal axis. The left column (a, c, e) represents pusher suspensions with a capillary number of 0.1, while the right column (b, d, f) represents puller suspensions with a shear rate of 0.05. Panel A: The line graph shows the specific viscosity ((μa - μ)/μ) for pusher suspensions. The vertical axis represents the specific viscosity, and the data points show a decreasing trend as the aspect ratio increases. Panel B: The line graph shows the specific viscosity ((μa - μ)/μ) for puller suspensions. The vertical axis represents the specific viscosity, and the data points show an increasing trend as the aspect ratio increases. Panel C: The line graph shows the contribution to the first normal stress differences (Σxx(p) - Σzz(p)/γμ) for pusher suspensions. The vertical axis represents the contribution to the first normal stress differences, and the data points show a decreasing trend as the aspect ratio increases. An inset diagram illustrates the coordinate system with velocity components. Panel D: The line graph shows the contribution to the first normal stress differences (Σxx(p) - Σzz(p)/γμ) for puller suspensions. The vertical axis represents the contribution to the first normal stress differences, and the data points show an increasing trend as the aspect ratio increases. Panel E: The line graph shows the contribution to the second normal stress differences (Σzz(p) - Σyy(p)/γμ) for pusher suspensions. The vertical axis represents the contribution to the second normal stress differences, and the data points show a fluctuating trend as the aspect ratio increases. Panel F: The line graph shows the contribution to the second normal stress differences (Σzz(p) - Σyy(p)/γμ) for puller suspensions. The vertical axis represents the contribution to the second normal stress differences, and the data points show a decreasing trend as the aspect ratio increases.