Rheology of a concentrated suspension of spherical squirmers: monolayer in simple shear flow

In this study, we consider a monolayer suspension of squirmers in a thin fluid film. The film is assumed as infinitely periodic, and the two surfaces are assumed to be flat and stress free. The film thickness $L_{z}$ is set as $L_{z}=2.1a$ throughout the study. The reason why we consider a monolayer suspension, instead of a 3-dimensional suspension, is because it allows us to handle a larger system size with a limited number of particles. Moreover, a film suspension can in principle be generated experimentally by using a soap film. As shown in Fig. 1, in-plane shear flow is induced in the monolayer suspension of squirmers; the $x$-axis is taken in the flow direction, the $y$-axis is taken in the velocity gradient direction, and the $z$-axis in taken perpendicular to the film plane. The background shear flow can be expressed as $(u_{x},u_{y},u_{z})=({\gamma}y,0,0)$. Gravity also acts in the $x-y$ plane, and $\alpha$ is the angle the gravitational acceleration ${\bm{g}}$ makes with the $-y$ axis, as given by Eq. (5) and shown in the figure. We assume that body centres and orientation vectors of squirmers are placed in the centre-plane of the fluid film. Thus, due to the symmetry of the problem, the squirmers remain in the same centre plane all the time.

The hydrodynamic interactions among squirmers in an infinite periodic monolayer suspension are calculated in the same manner as in our former studies (Ishikawa & Pedley, 2008; Ishikawa et al., 2008), so we explain the methodology only briefly here. The Stokesian dynamics method (Brady & Bossis, 1988; Brady et al., 1988) is employed. The hydrodynamic interactions among an infinite suspension of particles are computed by the Ewald summation technique. By exploiting the Stokesian dynamics method, the force ${\bm{F}}$, torque ${\bm{T}}$ and stresslet S of squirmers are given by:

where R is the resistance matrix, ${\bm{U}}$ and ${\bm{\Omega}}$ are the translational and rotational velocities of a squirmer, $\langle{\bm{u}}\rangle$ and $\langle{\bm{\omega}}\rangle$ are the translational and rotational velocity of the bulk suspension, and $\langle\mbox{\bf{E}}\rangle$ is the rate of strain tensor of the bulk suspension. ${\bf Q}_{sq}$ is the irreducible quadrupole providing additional accuracy, which is approximated by its mean-field value (cf. (Brady et al., 1988)). Index far or near indicates far- or near-field interaction, and 2B or Sq indicates interaction between two inert spheres or two squirmers, respectively.

The computational region is a square of side $L$, and a suspension of infinite extent is modelled with periodic boundary conditions in the $x$ and $y$ directions. In order to express the stress free surfaces of the fluid film, the monolayer is periodically replicated also in the $z$-direction with $2.1a$ intervals. For the simple shear flow in the $x$, $y$-plane, the periodic conditions in the $x$ and $z$-directions are straightforward. In the $y$-direction, the periodicity requires a translation in the $x$-direction of the spheres at $y=L$, relative to those at $y=0$, by an amount $L\gamma t$ in order to preserve the bulk linear shear flow, where $t$ is time. The number of squirmers $N$ in a unit domain is set as $N=128$. The areal fraction of squirmers is $\phi$ (not the volume fraction $c$), and it is varied in the range $0.1\leq\phi\leq 0.75$. Parameters in Eq. (6) for the repulsive force are set as $F_{0}=10$ and $\tau=100$. The initial configuration of the suspended squirmers is generated by first arranging them in a regular array and then applying small random displacements and random orientations. The dynamic motions are calculated afterwards by the fourth-order Adams-Bashforth time marching scheme. The computational time step is set as $\Delta t=5\times 10^{-4}/\gamma$.

The apparent shear viscosity $\eta$ of the film suspension can be calculated from the suspension averaged stresslet $\langle\mbox{\bf{S}}\rangle$. The excess apparent viscosity is given by

where $c$ is the volume fraction. The areal fraction $\phi$ and $c$ satisfy the relation $c=4a\phi/(3L_{z})$, where $L_{z}$ ($=2.1a$) is the film thickness. When squirmers are torque-free, the stresslet is symmetric and $\eta_{xy}=\eta_{yx}=\eta$. When squirmers are bottom-heavy, on the other hand, the stresslet becomes asymmetric and $\eta_{xy}\neq\eta_{yx}$. The dimensionless first and second normal stress differences are defined as

In order to obtain suspension-averaged properties, stresslet values are averaged over all particles during the time interval $15\leq t\gamma\leq 30$, given that the ensemble values reach the steady state.

We introduce a dimensionless number $Sq$, which is the swimming speed, scaled with a typical velocity in the shear flow, i.e.

In simulation cases with $Sq>1$, parameters are modified as $F_{0}=10Sq$, $\Delta t\gamma=5\times 10^{-4}/Sq$, to prevent overlapping of particles. Stresslet values are averaged during the longer time interval $75/Sq\leq t\gamma\leq 150/Sq$, to obtain the steady state values.

We first calculated the apparent viscosity $\eta$ of the suspension, for both inert and squirming spheres. The results are plotted as a function of areal fraction $\phi$ in Fig. 2(a). We see that $\eta$ of a suspension of squirmers with $Sq=1$ and $\beta=+1$ (i.e. pullers) increases rapidly with $\phi$. By comparing with the results for inert spheres in the present study, it is found that the motility of a squirmer considerably increases the viscosity. In the case of inert spheres, layers of particles are formed along the flow direction, and near-field interactions between particles are reduced. In the case of squirmers, on the other hand, such layers are destroyed by the motility of the squirmers, which move around irregularly, leading to frequent near-field interactions. Near-field interactions generate large lubrication forces, which result in large particle stresslets as well as the large apparent viscosity.

In Fig. 2, Einstein’s equation in the dilute limit as well as the former numerical results of Singh & Nott (2000), for inert spheres, are also plotted for comparison. We see that all results agree well with Einstein’s equation in the dilute regime. The main differences between the present study and Singh & Nott (2000) are the boundary condition and the repulsive interparticle force. In Singh & Nott (2000) a monolayer suspension is sheared between two walls, while in the present simulation the shear is generated without a wall boundary. The coefficients of repulsive force given by Eq. (6) are $\eta_{0}a^{2}\gamma F_{0}=10^{-4}[N],\tau=100[-]$ in Singh & Nott (2000), while $F_{0}=10[-],\tau=100[-]$ in the present simulation. We see that the apparent viscosity of a suspension of inert particles reported by Singh & Nott (2000) is larger than that in the present study. The discrepancy may come from the differences in the boundary condition and the repulsive interparticle force. By introducing the wall boundary, motions of spheres in the $y$-direction are restricted. The restricted motions may enhance the near-field interactions between particles, which results in the larger apparent viscosity. Moreover, the repulsive force may be larger in the present study than in Singh & Nott (2000) (it is hard to find the value of the dimensionless coefficient analogous to $F_{0}$ in that paper). This is because squirmers easily overlap if the repulsive forces are not strong enough, given that lubrication flow between two squirming surfaces cannot mathematically prevent the overlapping. The strong repulsive force may reduce the apparent viscosity, which may be another reason for the discrepancy.

The stresslet can be decomposed into the hydrodynamic contribution and the repulsive contribution, as shown in Fig. 2(b). The repulsive contribution to the particle bulk stress can be calculated as (Batchelor, 1977)

where $V$ is a unit volume, $\mbox{\boldmath$r$}^{ij}$ is the centre-centre separation of squirmers $i$ and $j$, and $\mbox{\boldmath$F$}^{ij}$ is their pairwise interparticle force given by Eq. (6). We see that the hydrodynamic contribution is dominant in the low $\phi$ regime. When $\phi\geq 0.6$, on the other hand, the repulsive contribution becomes larger than the hydrodynamic contribution, which is the reason why the apparent viscosity rapidly increases in the high $\phi$ regime.

First and second normal stresses also appear in the suspension of inert spheres and squirmers, as shown in Fig. 3. $N_{1}$ is negative in sign, so particles are basically compressed in the flow direction. $|N_{1}|$ increases as $\phi$ is increased, similar to the apparent viscosity. The value of $|N_{1}|$ in a suspension of inert particles, as reported by Singh & Nott (2000), is larger than in the present study. The discrepancy may again come from the differences in the boundary condition and the repulsive interparticle force.

$N_{2}$ in the suspension of inert spheres is also negative in sign, so the particle stresses satisfy $\Sigma^{(p)}_{xx}<\Sigma^{(p)}_{yy}<\Sigma^{(p)}_{zz}$. In the case of squirmers, the sign of $N_{2}$ changes at around $\phi=0.55$. The hydrodynamic contribution to $N_{2}$ is positive and steadily increases with $\phi$, while the repulsive contribution to $N_{2}$ is negative and rapidly increases with $\phi$. Since the repulsive contribution overwhelms the hydrodynamic contribution when $\phi\geq 0.6$, $N_{2}$ becomes negative in the large $\phi$ regime.

The difference between pushers and pullers can be investigated by changing the swimming mode parameter $\beta$, which is positive for pullers and negative for pushers. Figure 4(a) shows the effective viscosity $\eta$ for $Sq=1$ and various values of $\beta$. We see that $\eta$ increases as the absolute value of $\beta$ is increased. This is probably because strong squirming velocity with large $|\beta|$ enhances near-field interactions between squirmers, which induces a strong stresslet. The effect of $\pm\beta$ is not symmetric, but pullers generate larger viscosity than pushers. In order to confirm these tendencies, we calculated the normalised probability density function of squirmers, defined as

where $P(r_{0}|r_{0}+r)dr$ is the conditional probability that, given that there is a squirmer centred at $r_{0}$, there is an additional squirmer centred between $r_{0}+r$ and $r_{0}+r+dr$. The results are plotted in Fig. 4(b). We see that $I(r)$ with $\beta=3$ has a considerably larger value than that with $\beta=0$ in the small $r$ regime, while $I(r)$ with $\beta=-3$ does not. On the other hand, $I(r)$ with $\beta=-3$ has a larger peak than that with $\beta=0$. Thus, near-field interactions between squirmers are increased as $|\beta|$ is increased. The difference between pushers and pullers is obvious in Fig. 4(b): pullers tend to come closer to each other than pushers. A former study reported that the face-to-face configuration is stable for puller squirmers while unstable for pusher squirmers, though the face-to-face configuration was expressed as a stress-free surface in the paper (Ishikawa, 2019). Thus, puller squirmers facing each other can come closer more easily than pushers.

We also investigated the effect of $Sq$, i.e. the ratio of swimming velocity to the shear velocity. The results are shown in Fig. 5. We see that the apparent viscosity $\eta$ increases as $Sq$ is increased. This is because large swimming velocity, relative to the shear velocity, destroys the formation of layers by the particles, and enhances their near-field interactions. This tendency can be confirmed by Fig. 5(b), in which $I(r)$ with $Sq=10$ has larger values than with $Sq=0.1$ and 0.5. Thus, increase in the apparent viscosity can be understood as coming from the increase of near-field interactions.

When squirmers are bottom-heavy, cells tend to align in the direction opposite to gravity. To discuss the effect of bottom-heaviness, we introduce a dimensionless number $G_{bh}$, defined as

$G_{bh}$ is proportional to the ratio of the time to swim a body length to the time for the axis orientation $\bm{p}$ of a non-swimmer to rotate from horizontal to vertical. When $G_{bh}$ is sufficiently large, the gravitational torque due to the bottom-heaviness balances the hydrodynamic torque due to the background shear flow, and squirmers tend to align in a particular direction.

The effect of $G_{bh}$ and the gravitational angle $\alpha$ on the apparent viscosity is shown in Fig. 6. Since an external torque due to the bottom-heaviness is exerted on a squirmer, the partcle stress tensor becomes asymmetric, and the $xy$ and $yx$ components become different. We see that $G_{bh}$ considerably affects the value of $\eta$. For $\beta=0$ and $-3$, in horizontal flow ($\alpha=0$), $\eta$ is decreased by the bottom-heaviness, and $\eta_{yx}$ with $\beta=-3$ even becomes negative when $G_{bh}\geq 50$. The decrease in the apparent viscosity may be caused by two mechanisms: (a) squirmers with large $G_{bh}$ tend to swim in the same direction, and cell-cell collisions are suppressed; and (b) aligned squirmers induce a net squirming stresslet when $\beta\neq 0$, which directly contributes to $\eta_{xy}$ and $\eta_{yx}$. The effect of $\alpha$ is also significant: $\eta$ with $\beta=3$ takes its maximum values around $\alpha=\pi/8$, whereas that with $\beta=-3$ takes its minimum values around $\alpha=\pi/8$. So the tendencies are almost opposite between the pusher and the puller.

In order to clarify the mechanism of bottom-heavy effects, we here discuss the orientation of bottom-heavy squirmers in shear flow. Figure 7 shows normalised probability density distribution as a function of angle $\zeta$ that is defined from the $x$-axis in the counter clockwise direction as shown in Fig. 7(b). When $I(\zeta)=1$ for all $\zeta$, the orientation distribution is isotropic. We see that the peak of $I(\zeta)$ is higher, and at a larger value of $\zeta$, in the $G_{bh}=100$ case than in the $G_{bh}=30$ case due to the strong bottom-heaviness. Squirmers with $G_{bh}=100$ are oriented with approximately $\zeta=0.38\pi$ regardless of $\beta$. When $\beta=-3$, i.e. pushers, the stresslet with $\alpha=0$ is directed as in Fig. 7(b) bottom-left, and the apparent viscosity decreases. If the direction of gravity is rotated to $\alpha=\pi/2$, as in Fig. 7(b) bottom-right, the direction of the stresslet becomes opposite, and the apparent viscosity increases. These schematics can qualitatively explain the results in Fig. 6. When the squirmers are pullers, the sign of the stresslet is opposite to that shown in Fig. 7(b). Thus, the apparent viscosity increases with $\alpha=0$ whereas it decreases with $\alpha=\pi/2$, which is again consistent with Fig. 6. We note that the internal configuration, at $\alpha=0$, is more regular for $\beta=-3$ than for $\beta=+3$; this is indicated by the higher peak for $\beta=-3$ in Fig. 7(a). For $\beta=+3$, large numbers of the squirmers are aggregated and almost jam the system.

The first and second normal stress differences are affected by the angle of gravity $\alpha$, as shown in Fig. 8. The first normal stress difference with $\beta=3$ increases as $\alpha$ is increased, whereas that with $\beta=-3$ decreases. The second normal stress difference with $\beta=3$ takes its minimum values around $\alpha=3\pi/8$, whereas that with $\beta=-3$ takes its maximum values around $\alpha=3\pi/8$. So the tendencies are again almost opposite between the pusher and the puller. The opposite tendency can be explained by the sign and rotation of the stresslet, as schematically shown in Fig. 7(b).

Last, we investigate the rheology under constant $Sq\cdot G_{bh}$ conditions. This product represents the ratio of the gravitational torque to the shear torque, independently of the swimming speed $V_{s}$. Thus, a solitary squirmer under constant $Sq\cdot G_{bh}$ conditions is expected to have the same orientation angle relative to gravity. The results for excess apparent viscosity and normal stress differences under the condition of $Sq\cdot G_{bh}=30$ are shown in Fig. 9. We see that $\eta$ is considerably increased in the small $G_{bh}$ regime. The effect of swimming is larger than that of the background shear in the small $G_{bh}$ regime. The large swimming velocity enhances near-field interactions between squirmers, which results in the large apparent viscosity. The normal stress difference $N_{2}$ for $\beta=-3$ is also enhanced in the small $G_{bh}$ regime. This is because the active stresslet, caused by the squirming velocity, plays a dominant role compared to the passive stresslet, caused by inert spheres, in the small $G_{bh}$ regime. Hence, the rheology under constant $Sq\cdot G_{bh}$ conditions is strongly affected by $Sq$ especially when it is large.

In this Section, we present a complementary method for calculating the rheological properties of a suspension of steady, spherical squirmers. The methodology for measuring the bulk viscosity is similar to that of a rheometer; a suspension of squirmers situated between flat parallel plates is sheared using Couette flow, and the drag force on the plates is measured. The overall dynamics of the active suspension are solved by summing pairwise lubrication interactions between closely-separated squirmers, along with hydrodynamic forces and torques associated with solitary squirmers in shear flow. This approach utilises elements of previous work (Brumley & Pedley, 2019), but with several key advances. Firstly, in the present formulation, an arbitrary areal fraction of squirmers is permitted, so that cells are no longer confined to a hexagonal crystal array. Secondly, the entire suspension will be subject to a background shear flow. The dependence of the rheology on a number of key suspension and external parameters are presented. Despite the key differences between this approach and that of Section 2, the results are in good quantitative agreement at sufficiently large areal fractions.

Two infinite no-slip plates situated at $y_{1}=+H/2$ and $y_{2}=-H/2$ move with velocities $\gamma y_{1}\bm{e}_{x}$ and $\gamma y_{2}\bm{e}_{x}$, respectively, so that the fluid between the plates is subject to a uniform shear with rate $\gamma$. That is, the fluid velocity is given by $\bm{u}=(u_{x},u_{y},u_{z})=(\gamma y,0,0)$. A suspension of $N$ identical squirmers is introduced into the fluid (see Fig. 10).

The position and orientation vectors of all squirmers are restricted to lie in the $x$-$y$ plane, so that each squirmer has $2+1$ degrees of freedom. Moreover, the domain is subject to periodic boundary conditions in the $x$-direction, with period $L$ such that the squirmer areal fraction $\phi=N\pi a^{2}/LH$, where $N=n_{x}n_{y}$. This can also be expressed in the following way:

where $\epsilon_{0}a$ is the spacing between adjacent squirmers in the case where they are arranged in a regular hexagonal array (as depicted in Fig. 10). The total dimensional force $\bm{F}_{i}$ and torque $\bm{T}_{i}$ on the $i^{\text{th}}$ squirmer are composed of several contributions, as outlined below. As in Brumley & Pedley (2019), these will be scaled by $\eta_{0}\pi a$ and $\eta_{0}\pi a^{2}$ respectively, so that they have units of velocity. As in Section 2, the parameter $\eta_{0}$ is the viscosity of the solvent around the squirmers. The resistance formulation for the spheres in Stokes flow is given by

where the resistance matrix is given by $\bm{R}=\bm{R}^{\text{sq-sq}}+\bm{R}^{\text{sq-wall}}+\bm{R}^{\text{drag}}$. The vectors $\bar{\bm{F}}=\bm{F}/(\eta_{0}\pi a)$ and $\bar{\bm{T}}=\bm{T}/(\eta_{0}\pi a^{2})$ are of length $2N$ and $N$ respectively, and are composed of the forces $[F_{x}^{(i)},F_{y}^{(i)}]$ and torques $[T_{z}^{(i)}]$ for all squirmers $i\in\{1,2,\dots,N\}$. The corresponding vectors $\bm{V}$ and $a\bm{\omega}$ contain the linear and angular velocities, respectively, of all squirmers in the suspension. At zero Reynolds number, all squirmers must be force- and torque-free. Equation (34) therefore reduces to

Explicit hydrodynamic coupling between squirmers is limited to lubrication interactions. For two squirmers $i$ and $j$, with minimum clearance $\epsilon_{ij}=(|\mathbf{x}_{i}-\mathbf{x}_{j}|-2a)/a$ (subject to periodic boundary conditions), lubrication interactions occur only if $\epsilon_{ij}<1$. This value is chosen because $\log\epsilon_{ij}$, the functional dependence of these forces and torques, is equal to zero at that point. These terms therefore increase continuously from zero as squirmers approach one another from afar. The matrix $\bm{R}^{\text{sq-sq}}$ captures the hydrodynamic forces and torques acting on every pair of spheres sufficiently close to one another, arising from their linear and angular velocities. These expressions, evaluated for no-slip spheres, can be found in Kim & Karrila (2005) and Brumley & Pedley (2019). In a similar fashion, $\bm{R}^{\text{sq-wall}}$ captures the forces and torques due to motion of the spheres in close proximity to the bounding walls. The matrices $\bm{R}^{\text{sq-sq}}$ and $\bm{R}^{\text{sq-wall}}$, composed of blocks for each pair of squirmers, are non-zero only for pairs that are sufficiently close to permit lubrication interactions. The sparsity pattern of these resistance matrices therefore depends on the physical configuration of the system, and must be computed at every time step. Conversely, the final component of the resistance matrix is always diagonal, and is given by

where $\bm{I}_{n}$ is the $n\times n$ identity matrix. This captures the drag on a solitary translating and rotating sphere in Stokes flow. Inclusion of this matrix ensures that even widely separated spheres that do not experience lubrication interactions, are subject to solitary Stokes drag.

There are a number of contributions to the forces and torques that are independent of the squirmer velocities. First, there are the contributions which depend on the arrangement of cells with respect to one another. The squirming motions generate contributions for all pairs that are sufficiently close to one another. We refer the reader to Brumley & Pedley (2019) for detailed expressions of these forces and torques, $\bar{\bm{F}}^{\text{sq}}$ and $\bar{\bm{T}}^{\text{sq}}$, respectively, noting that a change of reference frame must be made for each pair, in order for the expressions to apply. Closely separated pairs of squirmers and squirmers near the no-slip plate also experience a repulsive force $\bar{\bm{F}}^{\text{rep}}$ parallel to the vector joining their centres (or perpendicular to the wall). This follows the same functional form as in Eq. (6), but in principle we allow squirmer-squirmer and squirmer-wall interactions to have different strengths $F_{0}$ and interaction ranges $\tau^{-1}$.

There are several forces and torques which do not require pairwise geometries. Every squirmer is subject to a propulsive force parallel to its orientation vector, $\bm{p}_{i}$, according to

where $V_{s}=2B_{1}/3$ is the swimming speed of a solitary squirmer. Each squirmer also experiences a gravitational torque in the $z$-direction due to bottom-heaviness:

where $\bm{g}=-g(\sin\alpha,\cos\alpha,0)$ and $g$ is the strength of gravity. Finally, the effect of the background shear flow is to exert a hydrodynamic force and torque on each squirmer as follows:

where $y_{i}$ is the $y$ coordinate of the $i^{\text{th}}$ squirmer. Before proceeding, it is instructive to consider the dynamics of the system if interactions between squirmers and the bounding plates are completely neglected. Under these conditions, the matrix system Eq. (39) reduces to the following:

Since $\bm{R}^{\text{drag}}$ is diagonal, it is easily inverted, and the following results are obtained:

Equation (52a) reveals that a solitary squirmer will swim at speed $V_{s}$ in the direction of its orientation, and be advected by the background shear flow in the $x$-direction. Similarly, the orientation of the squirmer will evolve according to the gravitational torque and background vorticity (Jeffery, 1922). In addition to these solitary dynamics (Eq. (52)), the complete resistance formulation in Eq. (39) includes hydrodynamic and steric effects through pairwise interactions, but only for sufficiently close squirmers.

The ability for “solitary” squirmers – i.e. squirmers with no neighbours within lubrication range – to propel themselves is a critical advancement of the present model. This maintains realistic behaviour at low areal fractions of the suspension, when it is quite feasible that a squirmer $i$ will have $\epsilon_{ij}>1\ \forall\ j$, and ensures that the matrix system in Eq. (39) is well conditioned for all configurations.

The dynamics of the sheared squirmer suspension are calculated numerically by solving Eq. (39) with the Matlab solver ode15s. The squirmers, initially distributed on a hexagonal close-packed array with mean spacing $\epsilon_{0}$, are each subject to a random translational perturbation within a disk of radius $\epsilon_{0}/2$, ensuring that cells are non-overlapping. Initially, the squirmers’ orientations are taken to be random.

For the configuration presented in Fig. 10, determining the effective suspension viscosity requires calculating the wall shear stress on each of the no-slip plates at $y=\pm H/2$. We emphasise that only squirmers which are sufficiently close to the walls to facilitate lubrication interactions will contribute to the wall shear stress. Of the full set of squirmers $S=\{1,2,\ldots,N\}$, the following subsets can be identified:

The sets $S_{\text{T}}$ and $S_{\text{B}}$ identify squirmers that have a clearance of less than $a$ with the top and bottom plates respectively (i.e. $\epsilon<1$), and therefore whose behaviour contributes to the lubrication forces and torques. The $x$-component of the force on the bottom plate is given by

where the three terms

represent the $x$-component of the force on the plate, due to the squirming motion of the $i^{\text{th}}$ sphere, the difference between the squirmer velocity, $\bm{V}_{i}$, and the plate velocity, $-(H\gamma/2)\bm{e}_{x}$, and the angular velocity of the squirmer, respectively. In Eq. (56a) we have made use of the fact that $B_{1}=\frac{3}{2}a\gamma Sq$ and $B_{2}=\beta B_{1}$. Similar expressions exist to calculate the tangential force, $\bar{F_{x}}^{\text{T}}$, on the top plate. Here, the value $\epsilon_{i}^{\text{B}}a$ represents the minimum clearance between the $i^{\text{th}}$ squirmer and the bottom plate. The repulsive force acts normal to the wall, and therefore does not contribute to the shear force. Although the background shear and gravitational torques influence the motion of the squirmers, their combined effects are encapsulated in Eqs. (56b) and (56c).

In dimensional form, the tangential shear stress on the top and bottom plates are given by $\sigma^{\text{T}}=\eta_{0}\pi a\bar{F_{x}}^{\text{T}}/(L\delta)$ and $\sigma^{\text{B}}=\eta_{0}\pi a\bar{F_{x}}^{\text{B}}/(L\delta)$, respectively, where $L=(2+\epsilon_{0})an_{x}$ and $\delta=2.1a$ is the thickness of the monolayer (see Section 2). The effective bulk viscosity is therefore equal to

The above expression utilises the mean value of the shear stress across both plates. The excess apparent viscosity is therefore equal to

Upon first glance, it appears that the viscosity depends on the system size through the denominator in Eq. (58). Although the number of squirmers interacting with the wall would depend on suspension micro-structure, to leading order the number of terms in Eq. (55) is proportional to $n_{x}$. Moreover, the terms in Eq. (55) have the dimensions of velocity, matching the dimensions of $a\gamma$ in Eq. (58). In what follows, we will compare the results of Section 2 using Eq. (21), with the present lubrication formulation Eq. (58).

We should note here that the lubrication-theory-based (LT) method does not give an obvious way to compute normal stresses, so the results that follow will concentrate on predicting the shear viscosity.

We first calculate the excess apparent viscosity in the absence of gravity. Suspensions of $N=132$ spheres with an areal fraction of $\phi=0.80$ were simulated over $t\in[0,150]$, for both active squirmers ($Sq=1$, $\beta=1$) and passive spheres ($Sq=0$). The suspension viscosity can be calculated at every time-step using Eq. (58), the results of which are plotted in Fig. 11(a). It is evident that the squirmer suspension (blue) has a higher mean viscosity than the suspension of passive spheres (black), and also exhibits greater excursions from the mean value.

The contribution to the mean squirmer suspension viscosity (blue curve in Fig. 11(b)) arising purely from linear and angular velocities is $\eta/\eta_{0}-1=20.6$, approximately 90% of the total mean viscosity ($\eta/\eta_{0}-1=22.8$). This is still considerably higher than the mean value for the passive sphere case $\eta/\eta_{0}-1=12.0$. Despite the propensity for active swimmers to redistribute themselves, the interior interactions still give rise to a higher suspension viscosity.

The areal fraction of cells was systematically varied in the range $0.5<\phi<0.8$ for both squirmers and passive spheres. A sufficiently long averaging window ($10<t<150$) was taken when calculating the time-averaged suspension viscosity. Figure 11(b) shows the results of the lubrication simulations (circles) together with the findings using Stokesian dynamics (dashed, cf. Fig. 2). There is good quantitative agreement between the complementary simulation methods across all values of $\phi$ studied.

The effect of the squirmers’ swimming properties was also studied. The two dimensionless parameters are $Sq=V_{s}/{a\gamma}$, the swimming speed relative to background shear rate, and $\beta=B_{2}/B_{1}$, which effectively controls the stresslet sign and strength. Figure 12(a) for $\beta=1$ shows that the suspension viscosity increases for small $Sq$ before plateauing around $Sq=1$. Increasing $Sq$ further does not result in a noticeable shift in the viscosity, on average.

The dependence of the suspension viscosity on the swimming mode $\beta$ is shown in Fig. 12(b) for two different areal fractions. The viscosity is higher for $\phi=0.7$ than for $\phi=0.6$ for all values of $\beta$ studied. The viscosity for pushers ($\beta<0$) does not vary significantly with $\beta$, whereas for pullers ($\beta>0$) the viscosity increases dramatically with $\beta$. This increase is much more marked than the Stokesian dynamics findings (see Fig. 4(a)). Cross channel probability distributions for the squirmer positions reveal that pullers spend more time near the boundaries than pushers do, and this likely has a corresponding effect on the suspension viscosity.

The effect of bottom-heaviness is to provide an external torque on each squirmer, which tends to reorient the cell in the gravitational field (see Fig. 10). In the absence of hydrodynamic interactions or external shear, the orientation of each squirmer $\bm{p}_{i}$ will become anti-aligned with the gravity vector $\bm{g}$ over a timescale that is inversely proportional to the normalised strength of gravity, $G_{bh}$. A series of simulations were performed in which the strength of gravity and the angle with respect to the shear flow, were independently varied.

Figure 13(a) illustrates the excess apparent viscosity as a function of $G_{bh}$ for three different squirming modes, $\beta$ (all with $Sq=1$, $\phi=0.7$, $\alpha=0$). For small $G_{bh}$, the ordering matches the results of Fig. 12(b), since hydrodynamic effects alone result in the strongest accumulation of pullers near the wall. As $G_{bh}$ is increased, the viscosity differences become more pronounced, with both neutral squirmers and pushers exhibiting $\eta/\eta_{0}-1\approx 0$, i.e. no enhancement over the solvent viscosity. The internal squirmer motions are similar to those computed by the SD method: pullers with high $G_{bh}$ generate near-jamming aggregates. As the gravitational angle, $\alpha$, is varied (for fixed $G_{bh}=100$), the curves exhibit cross-over points (see Fig. 13(b)). For $-\pi/4\lesssim\alpha\lesssim\pi/4$, puller suspensions exhibit the greatest viscosity. However, for $\alpha\gtrsim\pi/4$ and $\alpha\lesssim-\pi/4$, pusher suspensions have a greater viscosity. In Fig. 13(b), the large peak in the $\beta=3$ curve around $\alpha=-\pi/4$, corresponds to a situation where the (clockwise) viscous torques due to background shear balance the (counterclockwise) gravitational torques for upward-swimming squirmers. This results in a system which is, to leading order, equivalent to the $G_{bh}=0$ case (1% difference between respective viscosities).

The dependence on $\alpha$ can be understood in terms of the fluid mechanics of propulsion in each case. For $\alpha\approx 0$, squirmers will tend to align vertically in the shear flow (i.e. perpendicular to the channel), so that their poles are closest to the walls. Moreover, since squirmers swim upwards, the lubrication interactions will be dominated by anterior poles interacting with the upper plate. Under these conditions, pullers will be drawn closer to the wall than pushers, and so will impart a greater lubrication drag on the top plate. For $\alpha\approx\pi/2$, gravity will tend to align squirmers parallel to the walls, so that equatorial regions of the squirmers are in the lubrication zones. The movement of fluid away from the poles in the case of pushers, means that under these circumstances, they will be drawn closer to the translating plates than their puller counterparts, and therefore exhibit a higher viscosity.