2.1. The ambient Poiseuille flow

Consider (figure 1) a channel between two parallel walls $W_1,W_2$ separated by a distance $H$. We use a Cartesian coordinate system $(x,y,z)$ with axes $(x,y)$ along the walls and $z$ perpendicular to the walls, such that the walls $W_1$ and $W_2$ are in the planes $z=0$ and $z=H$, respectively.

Figure 1. Sketch of a single sphere with centre $O'$ and radius $a$ freely suspended in a Poiseuille flow between two parallel solid slip walls $W_1$ and $W_2.$

The Poiseuille fluid flow is driven by an imposed constant pressure gradient ${\rm \Delta} p_0/L$ along the $x$ direction, where $L$ is a length along $x$ that is large compared with $H$ (but small compared with a characteristic channel dimension in the $y$ direction) and ${\rm \Delta} p_0=p_0(L/2)-p_0(-L/2)$ is the negative pressure drop between the far away upstream and downstream sections.

In this article, we consider slip walls with equal slip lengths $\lambda$ on both walls. The expression for the Poiseuille flow velocity along $x$ is easily found to be (Feuillebois, Bazant & Vinogradova Reference Feuillebois, Bazant and Vinogradova2009; Ng Reference Ng2011)

(2.1)\begin{equation} u_0(z)=k_1(z+\lambda)+k_2 z^{2} \end{equation}

and the corresponding pressure may be written as

(2.2)\begin{equation} p_0(x)= 2\mu_0 k_2 x \end{equation}

with the constant coefficients

(2.3a,b)\begin{equation} k_1=-\frac{H}{2\mu_0} \frac{{\rm \Delta} p_0}{L} ,\quad k_2=-\frac{k_1}{H}, \end{equation}

which are independent of the slip length. The local shear rate at a position $z$ between the planes

(2.4)\begin{equation} \dot{\gamma}(z)=k_1 + 2 k_2 z \end{equation}

is thus also independent of $\lambda$.

The volume flow rate $\bar {u}_0 H$ (per unit channel length in the transverse direction $y$) is expressed in terms of the fluid velocity averaged across the channel width:

(2.5)\begin{equation} \bar{u}_0 = \frac 1H \int_{z=0}^{H} u_0(z) \, \textrm{d}z. \end{equation}

By inserting (2.1) and (2.3a,b) into (2.5) we obtain the expression

(2.6)\begin{equation} \frac{{\rm \Delta} p_0}{L} = -\frac{12}{H(6\lambda+H)} \, \mu_0 \bar{u}_0 \end{equation}

for the pressure drop per unit length of the channel.

2.2. Definition of the effective viscosity of the suspension

Consider now that the channel contains a dilute suspension of solid no-slip particles of typical length scale $a$. The suspension is entrained along $x$ by a constant pressure gradient that is, by comparison with the pure fluid, affected by the presence of particles.

Similarly to (2.9) in (I), the effective viscosity of the suspension is defined from the linear relationship between the volume flow rate $\bar {u}_0H$ and the average pressure gradient that is necessary to drive the flow field on a distance $L$ such that $L/H\to \infty$ with $H/a$ fixed. Like in (I), there are no effects of the far away container walls in the $y$ direction at distances large compared with $L$. Then, taking an average (denoted by $\left\langle \cdot \right\rangle$) over all possible positions and orientations of the particles

(2.7)\begin{equation} \left\langle \lim_{L\to\infty}\frac{{\rm \Delta} P}{L}\right\rangle = -\frac{12}{H(6 \lambda+H)} \, \mu_{\textit{eff}} \,\bar{u}_0 \end{equation}

defines the effective viscosity $\mu _{\textit{eff}}$, the search for which is the goal of this article.

2.3. Application of Lorentz reciprocal theorem

The flow field perturbed by the presence of particle $i$ is sought as the sum of the unperturbed flow (indicated by subscript $0$) and a perturbation (indicated by a prime). Accordingly, the perturbed flow velocity, stress tensor and pressure are written as

(2.8a–c)\begin{gather} {\boldsymbol{u}} = {\boldsymbol{u}}_0 + {\boldsymbol{u}}',\quad \boldsymbol{\sigma} = \boldsymbol{\sigma}_0 + \boldsymbol{\sigma}',\quad p = p_0 + p' . \end{gather}

It will not be necessary to calculate the perturbed flow field in full detail. Indeed, the method of solution uses the Lorentz reciprocal theorem, as in (I) (3.1a,b):

(2.9)\begin{equation} \int_{(S_i\cup S_L\cup S_{\textit{width}}\cup W)} {\boldsymbol{u}}_0\boldsymbol{\cdot} \boldsymbol{\sigma}'\boldsymbol{\cdot} {\boldsymbol{n}} \,\textrm{{d}}S = \int_{(S_i\cup S_L\cup S_{\textit{width}}\cup W)} {\boldsymbol{u}}'\boldsymbol{\cdot} \boldsymbol{\sigma}_0\boldsymbol{\cdot} {\boldsymbol{n}} \,\textrm{{d}}S , \end{equation}

where (see the control volume in figure 2) $S_i$ is the surface of a given particle labelled $i$, $S_L$ is the set of surfaces tending to infinity upstream and downstream, $S_{width}$ is the set of surfaces tending to infinity in the $y$ direction (not shown in figure 2) and $W=W_1\cup W_2$.

Figure 2. Control volume for the application of the Lorentz reciprocal theorem.

The derivation is as in (I):

1. (i) Integrals over $S_i$ take into account the no-slip boundary condition on the particle and their calculation is the same as in (I). They involve the force, torque, stresslet and quadrupole components $S_{i,xz}$ and $Q_{i,xzz}$. These quantities $S_{i,xz}$ and $Q_{i,xzz}$ are, respectively, the $(x,z)$ component of the second-order stresslet tensor ${\boldsymbol {S}}_i$ and the $(x,z,z)$ component of the third-order quadrupole tensor ${\boldsymbol {Q}}_i$ on particle $i$. These tensors are defined as

   (2.10)\begin{gather} {\boldsymbol{S}}_i = \int_{S_i} \left[\tfrac{1}{2}({\boldsymbol{r}}_i\, {\boldsymbol{f}}+ {\boldsymbol{f}} {\boldsymbol{r}}_i)-\tfrac{1}{3} ({\boldsymbol{r}}_i\boldsymbol{\cdot} {\boldsymbol{f}}) {\boldsymbol{I}} \right]\textrm{d}S, \end{gather}
   (2.11)\begin{gather} {\boldsymbol{Q}}_i = \int_{S_i} {\boldsymbol{f}} \; {\boldsymbol{r}}_i {\boldsymbol{r}}_i \,\textrm{{d}}S. \end{gather}
   Here, ${\boldsymbol {f}}$ is the stress of the perturbed flow field on particle $i$ and ${\boldsymbol {r}}_i$ is a vector pointing from the centre of particle $i$ to a point on its surface $S_i$.

2. (ii) Integrals on $S_L$ give the supplementary force difference ${\rm \Delta}\,f'$ (between the upstream and downstream far away sections) that is necessary to drive the particle $i$.

3. (iii) Integrals on $S_{width}$ (in the $y$ direction) vanish.

The main difference from (I) is in the calculation of the integrals on the slip walls $W_1$ and $W_2$ in (2.9), which cancel out, as we will see now.

For this purpose, we rewrite the integrands in terms of their components in the $(x,y,z)$ directions. The unperturbed flow velocity ${\boldsymbol {u}}_0$ on the wall is in the $x$ direction and so is the stress of the unperturbed flow field on the wall:

(2.12a,b)\begin{equation} {\boldsymbol{u}}_0 = u_0 \, {\boldsymbol{e}}_x,\quad \boldsymbol{\sigma}_0\boldsymbol{\cdot} {\boldsymbol{n}} = f_0 \, {\boldsymbol{e}}_x, \end{equation}

where $u_0$ and $f_0$ are, respectively, the velocity (2.1) and stress components in the $x$ direction. The perturbed flow velocity ${\boldsymbol {u}}'$ on the wall has components in the $x$ and $y$ directions since the wall is impermeable; now, the stress of the perturbed flow may have all its components $(\,f'_x,f'_y,f'_z)$:

(2.13a,b)\begin{equation} {\boldsymbol{u}}'=u'_x {\boldsymbol{e}}_x+u'_y {\boldsymbol{e}}_y,\quad \boldsymbol{\sigma}'\boldsymbol{\cdot} {\boldsymbol{n}} = f'_x {\boldsymbol{e}}_x+f'_y {\boldsymbol{e}}_y+f'_z {\boldsymbol{e}}_z. \end{equation}

The integral on $W_1$ on the left-hand side of (2.9) is then written

(2.14)\begin{equation} \int_{W_1} u_0\, f'_x \,\textrm{{d}}S \end{equation}

and the integral on $W_1$ on the right-hand side of (2.9) is written

(2.15)\begin{equation} \int_{W_1} u'_x\, f_0 \,\textrm{{d}}S. \end{equation}

The slip conditions on the wall $W_1$ are written for the unperturbed and perturbed flow fields in the $x$ direction as

(2.16a,b)\begin{equation} u_0 = \frac{\lambda}{\mu} f_0,\quad u'_x=\frac{\lambda}{\mu}\, f'_x, \end{equation}

where $\lambda$ is the slip length on $W_1$. Substituting these conditions into (2.14) and (2.15) shows that these integrals are equal. Thus they compensate in (2.9). Likewise, it can be shown that the integrals on $W_2$, where the same slip length $\lambda$ applies, compensate in (2.9). Therefore, even though the integrals on the walls do not vanish as in the no-slip case in (I), the same formal approach may be carried out here for the slip walls.

Finally, from (2.9) the supplementary necessary force difference ${\rm \Delta}\,f'$ between the upstream and downstream far away sections that is necessary to drive the particle $i$ is (as in (3.20a,b) in (I))

(2.17)\begin{equation} u_0(z_i)F_{i,x}+ \dot{\gamma}(z_i)C_{i,y}+ \dot{\gamma}(z_i) S_{i,xz} + k_2 Q_{i,xzz} + {\rm \Delta} f' \,\bar{u}_0 =0 . \end{equation}

Here, $F_{i,x}$ and $C_{i,y}$ are, respectively, the force in direction $x$ and the torque in direction $y$ on particle $i$.

As in (I), (2.17) is generalized to a suspension of $N$ particles ($i=1\ldots N$) contained in a control volume $V$. We then obtain the supplementary force difference ${\rm \Delta} F' = \sum _{i=1}^{N} f'$ that is necessary to drive the suspension. Dividing ${\rm \Delta} F'$ by the volume $V$ containing the $N$ particles and taking the thermodynamic limit for $L\to \infty$ (with $N\to \infty$ and $n=N/V$ kept finite) yields the following result (like (3.22) in (I)):

(2.18)\begin{equation} \lim_{ L\to\infty} \frac{{\rm \Delta} P'}{ L} = -\frac{1}{\bar{u}_0} \lim_{ L\to\infty} \frac{1}{V} \sum_{i=1}^{N} \left[ u_0(z_i)F_{i,x}+ \dot{\gamma}(z_i)C_{i,y}+ \dot{\gamma}(z_i) S_{i,xz} + k_2 Q_{i,xzz} \right]. \end{equation}

From (2.7), the average $\left\langle {.} \right\rangle$ over all possible positions and orientations of the particles of this quantity provides the effective viscosity $\mu _{\textit{eff}}$. For force-free and torque-free particles, $F_{i,x}=C_{i,y}=0$, the result is

(2.19)\begin{equation} \mu_{\textit{eff}} = \mu_0 + \frac{H(H+6\lambda)}{12\bar{u}_0^{2}} \left\langle \lim_{ L\to\infty} \frac{1}{V} \sum_{i=1}^{N} \left[ \dot{\gamma}(z_i) S_{i,xz} + k_2 Q_{i,xzz} \right] \right\rangle . \end{equation}

This formula is quite general, in the sense that it applies to any particle shape in any polydisperse dilute suspension.

2.4. Dilute monodisperse suspension of spheres

The considered suspension is monodisperse. It is also assumed to be dilute in the sense that its volume fraction $\phi =Nv/V$ is small compared with unity, where $v$ is the volume of one particle. It is standard for this suspension to introduce the intrinsic viscosity defined as $[\mu ]=(\mu_{\textit{eff}}-\mu _0)/(\mu _0 \phi )$.

In a monodisperse suspension all particle contributions are statistically equivalent, so that (2.19) yields the intrinsic viscosity:

(2.20)\begin{equation} [\mu] = \frac{H(H+6\lambda)}{12 \bar{u}_0^{2} v \mu_0} \left\langle \dot{\gamma}(z) S_{xz} + k_2 Q_{xzz} \right\rangle . \end{equation}

Here the index $i$ has been dropped for a test particle that is identical to the other ones. The centre of this particle is denoted $z_c$. Equation (2.20) replaces (4.5) in (I); $H$ there is replaced here by $H(H+6\lambda )$ in the numerator of the prefactor; moreover, the stresslet $S_{xz}$ and quadrupole $Q_{xzz}$ components both depend on the slip length $\lambda$ on the walls. Recalling (2.3a,b) and writing all constants in terms of $k_1$, the final result is

(2.21)\begin{equation} [\mu] = \frac{3}{\mu_0 v k_1(1+6\lambda/H)} \left\langle \left( 1 -2 \frac{z_c}{H} \right) S_{xz} -\frac{Q_{xzz}}{H} \right\rangle . \end{equation}

This equation replaces the result (4.6) for no-slip walls in (I).

Consider now spherical particles of radius $a$ and volume $v=\frac 43 {\rm \pi} a^{3}$. We define the dimensionless channel width, particle centre coordinate and slip length as

(2.22a–c)\begin{equation} \tilde{H}=\frac{H}{a} ,\quad \tilde{z}_{c}=\frac{z_c}{a},\quad \tilde{\lambda}=\frac{\lambda}{a} \end{equation}

and the dimensionless stresslet and quadrupole as

(2.23a)\begin{gather} \tilde{S}_{xz} = \frac{S_{xz}}{\mu_0 a^{3} k_1},\quad \tilde{Q}_{xzz} = \frac{Q_{xzz}}{\mu_0 a^{4} k_1}. \end{gather}

Equation (2.21) can then be recast in dimensionless form as

(2.24)\begin{equation} [\mu] = \frac{9}{4{\rm \pi}(1+6 \tilde{\lambda}/ \tilde{H})} \left\langle \left( 1 -2 \frac{ \tilde{z}_{c}}{ \tilde{H}} \right) \tilde{S}_{xz} -\frac{1}{ \tilde{H}} \tilde{Q}_{xzz} \right\rangle . \end{equation}

As in (I) and in Navardi & Bhattacharya (Reference Navardi and Bhattacharya2010), a uniform distribution of spherical particles is considered in the channel, in the sense that particle centres occupy all possible positions (excluding overlap with walls) with an equal probability. In (I) that was concerned with elongated particles, it was assumed that a ‘high-frequency’ flow field is applied to the suspension, so that the orientation angle of each elongated particle relative to the flow field is practically kept constant during the period of the external flow oscillation. Now for spheres, the orientation problem is irrelevant and for a particle distribution to stay unmodified, it is permissible for particles to move a significant distance along the flow field provided their centre stays at the same position $z$. In that case, a uniform particle distribution set in a steady Poiseuille flow field stays as such later on. A possible migration of particles across streamlines (across $z$) would be attributed to either lift forces due to fluid inertia, concentration effects or Brownian motion. All these phenomena are out of the scope of the present article, provided: (i) the Reynolds number responsible for migration (Ho & Leal Reference Ho and Leal1974) is low enough and also the distance travelled along the flow field is not large (otherwise a low migration velocity across streamlines could give a significant migration distance); (ii) the volume fraction is low enough for interactions between particles to be negligible; and (iii) the particles are non-Brownian, that is, larger than a few micrometres. The average $\left\langle \cdot \right\rangle$ then only amounts to calculating an integral over all possible positions of the test particle across the channel, excluding the possible overlap of the sphere with walls:

(2.25)\begin{equation} [\mu] = \frac{9}{\displaystyle 4{\rm \pi}\left(1+6\frac{ \tilde{\lambda}}{ \tilde{H}} \right) \left( \tilde{H}-2\right)} \int_1^{ \tilde{H}-1} {\mathcal{J}}( \tilde{z}_{c})\, \textrm{d} \tilde{z}_{c}, \end{equation}

where the integrand ${\mathcal {J}}( \tilde{z}_{c})$ contains a stresslet part ${\mathcal {J}}_S( \tilde{z}_{c})$ and a quadrupole part ${\mathcal {J}}_Q( \tilde{z}_{c})$:

(2.26a)\begin{gather} {\mathcal{J}}( \tilde{z}_{c}) = {\mathcal{J}}_S( \tilde{z}_{c})+ {\mathcal{J}}_Q( \tilde{z}_{c}), \end{gather}

(2.26b,c)\begin{gather} {\mathcal{J}}_S( \tilde{z}_{c})= \left( 1 -2 \frac{ \tilde{z}_{c}}{ \tilde{H}} \right) \tilde{S}_{xz} , \quad {\mathcal{J}}_Q( \tilde{z}_{c})= -\frac{1}{ \tilde{H}} \tilde{Q}_{xzz} . \end{gather}

3.1. Interactions based on individual walls

In general, the hydrodynamic interactions of the sphere with both slipping walls should be taken into account. As stated in the introduction, this is a difficult task. Indeed, for various calculation methods such as the boundary element method (BEM) (see Pozrikidis Reference Pozrikidis1992; Sellier Reference Sellier, Aliabadi and Wen2010) and the method of multipoles (see Bhattacharya, Bławzdziewicz & Wajnryb Reference Bhattacharya, Bławzdziewicz and Wajnryb2005a,Reference Bhattacharya, Bławzdziewicz and Wajnrybb; Ekiel-Jeżewska & Wajnryb Reference Ekiel-Jeżewska, Wajnryb, Feuillebois and Sellier2009) it would first require the derivation of the Green tensor between two parallel slip walls. This very challenging issue is left for future works.

For a large enough gap $\tilde{H}$ between walls, in the case of no-slip walls, the earlier calculation by (I) has shown that the interactions of particles with the nearest wall give a good approximation to the effective viscosity, by comparison to the accurate multipole method which takes both walls simultaneously into account. For the case of a freely moving particle between no-slip parallel walls, the earlier study by Pasol et al. (Reference Pasol, Martin, Ekiel-Jeżewska, Wajnryb, Bawzdziewicz and Feuillebois2011) showed that the nearest wall approximation and the wall superposition approximation both gave excellent results for the particle velocity for a wide range of parameters. Based on these experiences, we consider here the following.

1. (i) The zero wall approximation, denoted as 0w, in which the stresslet and quadrupole in (2.26b,c) are calculated ignoring the interactions of the particle with the walls, that is (see (I))

   (3.1a,b)\begin{equation} \tilde{S}_{xz, \textit{0w}} = \frac{10{\rm \pi}}{3} \left(1-\frac{2 \tilde{z}_{c}}{ \tilde{H}}\right) , \quad \tilde{Q}_{xzz, \textit{0w}} = -\frac{8{\rm \pi}}{3 \tilde{H}} , \end{equation}
   and injected into the expression (2.26) for the integrand, to be denoted as ${\mathcal {J}}_ \textit {0w}$; the result for the intrinsic viscosity calculated with (2.25) is then denoted as $[\mu ]_ \textit {0w}$.

2. (ii) The nearest wall approximation, denoted as nw, in which the stresslet and quadrupole in (2.26b,c) are calculated by considering interactions with the nearest wall only; the integrand is then denoted as ${\mathcal {J}}_ \textit {nw}$ and the intrinsic viscosity $[\mu ]_ \textit {nw}$ is calculated by integration using (2.25). By symmetry of ${\mathcal {J}}_ \textit {nw}$ with respect to the channel mid-plane, a simpler formula is

   (3.2)\begin{equation} [\mu]_ \textit{nw} = \frac{9}{\displaystyle 2{\rm \pi}\left(1+6\dfrac{\tilde{\lambda}}{ \tilde{H}} \right) \left( \tilde{H}-2\right)} \int_1^{ \tilde{H}/2} {\mathcal{J}}_ \textit{nw}( \tilde{z}_{c})\, \textrm{d} \tilde{z}_{c} . \end{equation}

3. (iii) The wall superposition approximation, denoted as ws. In this case, consider first the influence of the wall at $\tilde{z}_{c}=0$ over the whole gap (note that ${\mathcal {J}}_ \textit {1w}= {\mathcal {J}}_ \textit {nw}$ only for $1\leq \tilde{z}_{c} \leq \tilde{H}/2$). The integrand, denoted as ${\mathcal {J}}_ \textit {1w}( \tilde{z}_{c})$, here takes into account the interactions with the wall at $\tilde{z}_{c}=0$. Integrating ${\mathcal {J}}_ \textit {1w}( \tilde{z}_{c})$ for $\tilde{z}_{c}$ from $1$ to $\tilde{H}-1$ provides the resulting contribution to the intrinsic viscosity:

   (3.3)\begin{equation} [\mu]_ \textit{1w} = \frac{9}{\displaystyle 4{\rm \pi}\left(1+6\dfrac{ \tilde{\lambda}}{ \tilde{H}} \right) \left( \tilde{H}-2\right)} \int_1^{ \tilde{H}-1} {\mathcal{J}}_ \textit{1w}( \tilde{z}_{c})\,\textrm{d} \tilde{z}_{c} . \end{equation}

   The calculation of the contribution of the second wall at $\tilde{z}_{c}= \tilde{H}-1$ to the intrinsic viscosity is carried out in the same way. Without entering into detail, it is clear that both walls will give equal contributions. These contributions are then added. In this process, the common part for large distance from the wall should not be counted twice, thus we subtract the 0w case. Finally, the formula for the wall superposition approximation is

   (3.4)\begin{equation} [\mu]_ \textit{ws} = 2\, [\mu]_ \textit{1w} - [\mu]_ \textit{0w} \end{equation}
   in which $[\mu ]_ \textit {0w}$ is calculated with (3.1a,b) and (2.25) as explained above. The integrand for this case is
   (3.5)\begin{equation} {\mathcal{J}}_ \textit{ws} = 2\, {\mathcal{J}}_ \textit{1w} - {\mathcal{J}}_ \textit{0w} . \end{equation}

In any case (either nw or ws), the quantities required for averaging are the stresslet and quadrupole on a freely suspended sphere in the ambient flow (2.1) in the presence of a single slip wall. The bipolar coordinates method (BCM) will be used in this article to obtain very accurate results for these quantities. This borrows results from previous papers to be quoted below and also provides novel derivations for the quadrupole, for four basic flow fields presented in § 3.2. Then § 3.3 presents the new calculation of the quadrupole for a sphere held fixed in ambient linear and quadratic shear flows.

3.2. Relevant normalized quantities

The vector ${\boldsymbol {r}}_i$ in the stresslet operator (2.10) and quadrupole operator (2.11) is now for a single particle (dropping the $i$ subscript) denoted as ${\boldsymbol {r}}= {\boldsymbol {x}}- {\boldsymbol {x}}_C$, where ${\boldsymbol {x}}$ is a current point on the sphere and ${\boldsymbol {x}}_C$ is the position of the sphere centre $C$ as measured from a point located on the lower wall $W_1$.

The translational and rotational velocities of a freely moving sphere in the flow field (2.1) are, by linearity of Stokes equations, respectively along $x$, say $U_x$, and along $y$, say $\varOmega _y$. By linearity of Stokes equations, they are obtained as the superposition of the corresponding freely moving velocities in a laminar shear flow $k_1z$ (see Loussaief et al. Reference Loussaief, Pasol and Feuillebois2015) and in a quadratic shear flow $k_2z$ (see Ghalia et al. Reference Ghalia, Feuillebois and Sellier2016):

(3.6a,b)\begin{gather} U_x = U^{s}_x + U^{q}_s,\quad\varOmega_y = \varOmega^{s}_y + \varOmega^{q}_y \end{gather}

(3.6c)\begin{gather} \frac{U^{s}_x}{k_1(z_c+\lambda)}= u^{s}_x=\frac{\displaystyle\frac{1}{2( \tilde{z}_{c}+ \tilde{\lambda})} f^{r}_{xy}c^{s}_{yx} + f^{s}_{xx} c^{r}_{yy} } { f^{t}_{xx}c^{r}_{yy} - f^{r}_{xy}c^{t}_{yx}} , \end{gather}

(3.6d)\begin{gather} \frac{ \varOmega^{s}_y}{k_1/2}= \omega^{s}_y=\frac{2( \tilde{z}_{c}+ \tilde{\lambda}) f^{s}_{xx}c^{t}_{yx} + f^{t}_{xx} c^{s}_{yx} } { f^{t}_{xx}c^{r}_{yy} - f^{r}_{xy}c^{t}_{yx}} , \end{gather}

(3.6e)\begin{gather} \frac{U^{q}_x}{k_2\left(z_c+\dfrac{a^{2}}{3}\right)}= u^{q}_x= {{ \tilde{z}_{c}^{2} c_{yy}^{r}\, f_{xx}^{q}+ \tilde{z}_{c}\, f_{xy}^{r} c_{yx}^{q}}\over{[ \tilde{z}_{c}^{2}+\frac 13] [c_{yy}^{r}\, f_{xx}^{t}-f_{xy}^{r} c_{yx}^{t}]}}, \end{gather}

(3.6f)\begin{gather} {{ \varOmega_q}\over{k_2 z_c}}=\omega^{q}_y={{c_{yx}^{q} f_{xx}^{t}+ \tilde{z}_{c}\, f_{xx}^{q} c_{yx}^{t}}\over{c_{yy}^{r}\, f_{xx}^{t}-f_{xy}^{r} c_{yx}^{t}}}. \end{gather}

We here also define dimensionless velocities. The ‘$f$'s and ‘$c$'s on the right-hand sides are dimensionless friction coefficients which either have the limit of unity (when subscripts are identical) or vanish (when subscripts are different) when the sphere is infinitely far from the wall. The results (3.6c) and (3.6d) were obtained (in Loussaief et al. Reference Loussaief, Pasol and Feuillebois2015) by writing that the sums of forces and torques for a sphere translating and rotating in a fluid at rest and a sphere held fixed in a linear shear flow $k_1z$ vanish. Likewise, (3.6e) and (3.6d) were obtained (in Ghalia et al. Reference Ghalia, Feuillebois and Sellier2016) from vanishing total force and torque for the quadratic shear flow $k_2z$. Adding the results in (3.6a) and (3.6b) amounts to expressing the condition that the total force and torque in the flow field (2.1) vanish. By total force and torque, we mean those due to translation and rotation in a fluid at rest and particle held fixed in the flow fields $k_1z$ and $k_2z$.

Expressions for the normalized stresslet and quadrupole defined in (2.23) are displayed in table 1, in terms of dimensionless stresslet coefficients $s^{\bullet }_{xz}$ and quadrupole coefficients $q^{\bullet }_{xzz}$, where the superscript bullet ($\bullet$) represents $t, r, s$ and $q$ for translation, rotation, ambient linear shear flow and ambient quadratic shear flow, respectively. The coefficients $s^{\bullet }_{xz}$ and $q^{\bullet }_{xzz}$ for $(t, s, q)$ tend to unity for $\tilde{z}_{c}\to \infty$ whereas $s^{r}_{xz}$ and $q^{r}_{xzz}$ vanish (the coefficients $6{\rm \pi}$ and $8{\rm \pi}$ in the definitions of these last quantities are chosen arbitrarily).

Table 1. Normalized stresslet and quadrupole (defined in (2.23)) for elementary flow fields: a translating sphere, a rotating sphere in a fluid at rest and a sphere held fixed in linear shear and quadratic shear flows.

In the last column of table 1, the normalized quadrupoles $q^{t}_{xzz}$ due to translation and $q^{r}_{xzz}$ due to rotation were obtained by applying the Lorentz reciprocity theorem (see appendix A):

(3.7)\begin{gather} q^{t}_{xzz} = 3 \tilde{z}_{c} \left( \tilde{z}_{c} + 2 \tilde{\lambda} \right) f^{t}_{xx} -6 \tilde{z}_{c} \left( \tilde{z}_{c} + \tilde{\lambda} \right) f^{s}_{xx}+3 \tilde{z}_{c}^{2}\, f^{q}_{xx}, \end{gather}

(3.8)\begin{gather} q^{r}_{xzz} = \frac 34 \tilde{z}_{c} \left( \tilde{z}_{c} + 2 \tilde{\lambda} \right) f^{r}_{xy}+ \tilde{z}_{c} c^{s}_{yx} - \tilde{z}_{c} c^{q}_{yx}. \end{gather}

All dimensionless coefficients, in these formulae and in table 1, except $q^{s}_{xzz}$ and $q^{q}_{xzz}$, have been calculated previously using bipolar coordinates for a sphere interacting with a slip plane wall. The results of Loussaief et al. (Reference Loussaief, Pasol and Feuillebois2015) and Ghalia et al. (Reference Ghalia, Feuillebois and Sellier2016) will be used here. For instance, the dimensionless stresslet coefficient for translation, $s^{t}_{xz}$, is obtained from earlier results as follows:

1. (i) From (2.23a) and the first line, first column in table 1, the stresslet component due to translation is $S_{xz}^{t}=6{\rm \pi} \mu _0 aU_x s^{t}_{xz}$, which is the original definition of $s^{t}_{xz}$ in Ghalia etal. (Reference Ghalia, Feuillebois and Sellier2016) (see (66) there).

2. (ii) In Ghalia et al. (Reference Ghalia, Feuillebois and Sellier2016), (B 4a–c) gives $s^{t}_{xz}$ in terms of coefficients $B_n^{t}, C_n^{t}, E_n^{r}$.

3. (iii) These coefficients appear in the expressions of harmonics (describing the flow field) as series in bipolar coordinates (see (3.3) in Loussaief et al. (Reference Loussaief, Pasol and Feuillebois2015)); $B_n^{t}, C_n^{t}$ are for the translation problem and $E_n^{r}$ is for the rotation problem. These coefficients are calculated accurately in bipolar coordinates in Loussaief et al. (Reference Loussaief, Pasol and Feuillebois2015).

The determination of the novel coefficients $q^{c}_{xzz}$ and $q^{q}_{xzz}$ is the topic of the next subsection.

The expressions of the stresslet and quadrupole for translation and rotation in table 1 are valid for any translation velocity $U_x$ and rotational velocity $\varOmega _y$, respectively. Consider now a freely moving sphere in the flow field (2.1) (with (2.3a,b)). Its translational and rotational velocities are obtained by adding the contributions of the linear and quadratic shear flows shown in (3.6):

(3.9a,b)\begin{gather} U_x=k_1(z_c+\lambda)u^{s}_x+ k_2\left(z_c^{2}+\frac{a}{3}\right)u^{q}_x, \quad \varOmega_y=\frac{k_1}{2} \omega^{s}_y + k_2 z_c \omega^{q}_y. \end{gather}

These expressions are then introduced in table 1. The stresslet and quadrupole for a freely moving sphere in the flow field (2.1) are obtained as superpositions of those for a sphere translating and rotating in a fluid at rest and those for a sphere held fixed in the shear flows $k_1z$ and $k_2z$:

(3.10)\begin{align} \tilde{S}_{xz} &= 6{\rm \pi} \left[ \left( \tilde{z}_{c}+ \tilde{\lambda}\right) u^{s}_x -\frac{1}{ \tilde{H}} \left( \tilde{z}_{c}^{2}+\frac 13 \right) u^{q}_x \right] s^{t}_{xz}, \nonumber\\ &\quad +6{\rm \pi} \left[ \frac{1}{2} \omega^{s}_y -\frac{ \tilde{z}_{c}}{ \tilde{H}} \omega^{q}_y \right] s^{r}_{xz} +\frac{10{\rm \pi}}{3} s^{s}_{xz} -\frac{20{\rm \pi}}{3} \frac{ \tilde{z}_{c}}{ \tilde{H}} s^{q}_{xz}, \end{align}

(3.11)\begin{align} \tilde{Q}_{xzz} &= -2{\rm \pi} \left[ \left( \tilde{z}_{c}+ \tilde{\lambda}\right) u^{s}_x -\frac{1}{ \tilde{H}} \left( \tilde{z}_{c}^{2}+\frac 13 \right) u^{q}_x \right] q^{t}_{xzz} \nonumber\\ &\quad +8{\rm \pi} \left[ \frac{1}{2} \omega^{s}_y -\frac{ \tilde{z}_{c}}{ \tilde{H}} \omega^{q}_y \right] q^{r}_{xzz} +2{\rm \pi} \left( \tilde{z}_{c}+ \tilde{\lambda}\right) q^{s}_{xzz} -2{\rm \pi} \frac{ \tilde{z}_{c}^{2}}{ \tilde{H}} q^{q}_{xzz} . \end{align}

3.3. Determination of the normalized quadrupoles $q^{s}_{xzz}$ and $q^{q}_{xzz}$

The reader only interested in the results for the effective viscosity may skip this subsection and go directly to § 3.4.

We derive here the normalized quadrupoles $q^{s}_{xzz}$ and $q^{q}_{xzz}$ (recall table 1) for a sphere, with radius $a$ and centre $O',$ held fixed near the $z=0$ plane slip wall in an ambient linear or quadratic shear flow $({\boldsymbol {u}}_a,p_a)$ with stress tensor $\boldsymbol {\sigma }_a$. The disturbance flow $({\boldsymbol {u}}',p')$ about the sphere has stress tensor $\boldsymbol {\sigma }'$ and it is recalled that the sphere distance to the wall is $z_c>a$ while the wall slip length is $\tilde{\lambda} a.$ The quadrupole $Q_{xzz}$, exerted by the disturbed flow on the sphere, is $Q_{xzz}=Q^{a}_{xzz}+Q^{d}_{xzz}$ with

(3.12a,b)\begin{equation} Q^{a}_{xzz}=\int[{\boldsymbol{e}}_x\boldsymbol{\cdot}{} \boldsymbol{\sigma}_a\boldsymbol{\cdot}{\boldsymbol{n}}](z-z_c)^{2}\,\textrm{d}S,\quad Q^{d}_{xzz}=\int[{\boldsymbol{e}}_x\boldsymbol{\cdot} \boldsymbol{\sigma}'\boldsymbol{\cdot}{\boldsymbol{n}}](z-z_c)^{2}\,\textrm{d}S. \end{equation}

From the definitions (2.23) and table 1, it follows that the dimensionless quadrupoles $q_{xzz}^{s}$ and $q_{xzz}^{q}$ are defined by

(3.13)\begin{gather} Q_{xzz}=2{\rm \pi}\mu_0 k_1 a^{4}( \tilde{z}_{c}+\tilde{\lambda})q^{s}_{xzz}\quad \textrm{for} \ {\boldsymbol{u}}_a={\boldsymbol{u}}_s:=k_1(z+a\tilde{\lambda}){\boldsymbol{e}}_x\ \textrm{and}\ p_a=0, \end{gather}

(3.14)\begin{gather} Q_{xzz}=2{\rm \pi}\mu_0 k_2 a^{3} z_c^{2}q^{q}_{xzz}\quad \textrm{for}\ {\boldsymbol{u}}_a={\boldsymbol{u}}_q:=k_2 z^{2}{\boldsymbol{e}}_x\ \textrm{and}\ p_a=2\mu_0 k_2 x. \end{gather}

Let us now derive the integrals in (3.12a,b). Following the original idea suggested by Dean & O'Neill (Reference Dean and O'Neill1963) for a sphere rotating parallel to a no-slip wall, it is convenient to use spherical polar coordinates $(r,\theta ,\varphi )$ with the sphere centre as the origin and $r=O'M, x=r\sin \theta \cos \varphi , y=r\sin \theta \sin \varphi , z=z_c+r\cos \theta .$ Here $\theta$ and $\varphi$ run in $[0,{\rm \pi} ]$ and $[0,2{\rm \pi} ],$ respectively. In addition, the disturbance velocity ${\boldsymbol {u}}'$ polar components are $(u'_r,u'_{\theta },u'_{\varphi }).$ On the $r=a$ sphere surface $S,$ with unit normal ${\boldsymbol{n}}=\boldsymbol{O}'\boldsymbol{M}/a$, we write for convenience of calculation ${\boldsymbol {e}}_x\boldsymbol {\cdot } \boldsymbol {\sigma }'\boldsymbol {\cdot }{\boldsymbol {n}}=\mu _0(T_1+T_2)$ with

(3.15)\begin{gather} aT_1=\cos\theta\cos\varphi\left[\left({\frac{\partial u'_r}{\partial\theta}}\right)-u'_{\theta}\right] +\sin\varphi[u'_{\varphi}]-{\frac{\sin\varphi}{\sin\theta}} \left({\frac{\partial u'_r}{\partial\varphi}}\right), \end{gather}

(3.16)\begin{gather} T_2=\left[2\left(\frac{\partial u'_r}{\partial r}\right)- \frac{p'}{\mu_0}\right]\sin\theta\cos\varphi +\cos\theta\cos\varphi\left({\frac{\partial u'_{\theta}}{\partial r}}\right) -\sin\varphi\left({\frac{\partial u'_{\varphi}}{\partial r}}\right). \end{gather}

The contribution of $T_e$ for $e=1,2$ to the quadrupole $Q^{d}_{xzz}$ is denoted $Q^{d,e}_{xzz}$. On the sphere boundary $(z-z_c)^{2}\,\textrm {d}S=a^{4}[\cos \theta ]^{2}\sin \theta \,\textrm {d}\theta \,\textrm {d}\varphi$ and ${\boldsymbol {u}}'=-{\boldsymbol {u}}_a$. Using (3.15) easily gives

(3.17)\begin{equation} Q^{a}_{xzz}=Q^{d,1}_{xzz}=0\quad \textrm{if}\ {\boldsymbol{u}}_a={\boldsymbol{u}}_s, \quad (Q^{a}_{xzz},Q^{d,1}_{xzz})={\frac{8{\rm \pi}\mu_0k_2 a^{5}}{5}}\left({\frac{2}{3}},{\frac{1}{7}}\right) \ \textrm{if}\ {\boldsymbol{u}}_a={\boldsymbol{u}}_q. \end{equation}

The calculation of $Q^{d,2}_{xzz}$ is tricky. It employs the cylindrical coordinates $(\rho ,z,\varphi )$ with $\rho =\{x^{2}+y^{2}\}^{1/2}$ being the distance to the $(O,{\boldsymbol {e}}_z)$ axis. The cylindrical components, $(v_{\rho },v_z,v_{\varphi })$, of the disturbance velocity ${\boldsymbol {u}}'$ satisfy

(3.18a–c)\begin{equation} u'_r=\sin\theta[v_{\rho}]+\cos\theta[v_z],\quad u'_{\theta}=\cos\theta[v_{\rho}]-\sin\theta[v_z],\quad u'_{\varphi}=v_{\varphi}. \end{equation}

The relations (3.18a–c) provide the following expressions for the radial derivatives of $(u'_r,u'_{\theta },u'_{\varphi }),$ required in (3.16) (note that $\theta$ is kept constant in these derivatives):

(3.19a)\begin{gather} \displaystyle{\frac{\partial u'_r}{\partial r}}= \sin\theta\left({\frac{\partial v_{\rho}}{\partial r}}\right) +\cos\theta\left({\frac{\partial v_z}{\partial r}}\right), \end{gather}

(3.19b,c)\begin{gather} \displaystyle{\frac{\partial u'_{\theta}}{\partial r}} = \cos\theta\left({\frac{\partial v_{\rho}}{\partial r}}\right) -\sin\theta\left({\frac{\partial v_z}{\partial r}}\right),\quad {\frac{\partial u'_{\varphi}}{\partial r}}={\frac{\partial v_{\varphi}}{\partial r}}. \end{gather}

For each considered ambient velocity field, either ${\boldsymbol {u}}_a^{s}$ (linear shear) or ${\boldsymbol {u}}_a^{q}$ (quadratic shear), the disturbance flow $({\boldsymbol {u}}',p')$ about the sphere is written (see e.g. Ghalia et al. Reference Ghalia, Feuillebois and Sellier2016)

(3.20a,b)\begin{gather} v_{\rho}=\frac{1}{2}\left\{\frac{{\rho} Q_1}{c}+U_0+U_2\right\} \cos\phi,\quad v_{\phi}={\frac{1}{2}} (U_2-U_0)\sin\phi, \end{gather}

(3.21a–c)\begin{gather} v_z=\frac{1}{2}\left\{\frac{zQ_1}{c}+2W_1\right\}\cos\phi, \quad p'=\mu_0 Q_1\cos\phi,\quad c=(z_c^{2}-a^{2})^{1/2}, \end{gather}

with functions $W_1, U_0, U_1$ and $U_2$ solely depending upon $(\rho ,z)$ (or equivalently upon $(r,\theta$)). Injecting (3.19a–c)–(3.21a–c) into (3.16) provides, by elementary manipulations not detailed here for conciseness,

(3.22)\begin{align} Q^{d,2}_{xzz}&=\mu_0\int_S a^{2} T_2[\cos\theta]^{2}\,\textrm{d}S={\rm \pi}\mu_0 a^{4}\int_0^{{\rm \pi}}\tilde{T}_2 [\cos\theta]^{2}\sin\theta \,\textrm{d}\theta,\end{align}

(3.23)\begin{align} \tilde{T}_2&={\frac{\sin\theta}{\sinh\alpha}} \left[1+{\frac{\cosh\alpha\cos\theta}{2}}\right] \left({\frac{\partial Q_1}{\partial r}}\right)+ \sin\theta\cos\theta\left({\frac{\partial W_1}{\partial r}}\right)\nonumber\\ &\quad +\left[1+{\frac{\sin^{2}\theta}{2}}\right] \left({\frac{\partial U_0}{\partial r}}\right) +{\frac{\sin^{2}\theta}{2}}\left({\frac{\partial U_2}{\partial r}}\right). \end{align}

In finding the disturbance flow, it is more convenient to determine the functions $W_1, U_0, U_1$ and $U_2$ in terms of the bipolar coordinates $(\xi ,\eta )$ which are related to $(\rho ,z)$ as follows (see e.g. Happel & Brenner Reference Happel and Brenner1973):

(3.24a–c)\begin{equation} \rho=\frac{c\sin\eta}{\cosh\xi-\cos\eta},\quad z=\frac{c\sinh\xi}{\cosh\xi-\cos\eta},\quad c=(z_c^{2}-a^{2})^{1/2}. \end{equation}

For each angle $\varphi$ the liquid domain above the slip wall is obtained taking $0\leq \eta \leq {\rm \pi}$ and $0\leq \xi \leq \alpha$ with $\alpha >0$ given by $\cosh \alpha = \tilde{z}_{c}.$ In the $(\xi ,\eta ,\varphi )$ coordinates the $\xi =0$ and $\xi =\alpha$ surfaces are the slip $z=0$ plane wall and the sphere boundary, respectively. It is then found (see Ghalia et al. Reference Ghalia, Feuillebois and Sellier2016) that

(3.25)\begin{gather} W_1(\xi,\eta)=k_lc^{l}(\cosh\xi-t)^{1/2}\sin\eta\sum_{n\geq 1}A_n\sinh( \gamma_n\xi) P'_n(t), \end{gather}

(3.26)\begin{gather} Q_1(\xi,\eta)=k_lc^{l}(\cosh\xi-t)^{1/2}\sin\eta\sum_{n\geq 1}[B_n\cosh(\gamma_n\xi) +C_n\sinh(\gamma_n\xi)]P'_n(t), \end{gather}

(3.27)\begin{gather} U_0(\xi,\eta)=k_lc^{l}(\cosh\xi-t)^{1/2}\sum_{n\geq 0}^{\infty}[D_n\cosh(\gamma_n\xi) +E_n\sinh(\gamma_n\xi)]P_n(t), \end{gather}

(3.28)\begin{gather} U_2(\xi,\eta)=k_lc^{l}(\cosh\xi-t)^{1/2}\sin^{2}\eta\sum_{n\geq 2}[F_n\cosh(\gamma_n\xi) +G_n\sinh(\gamma_n\xi)]P''_n(t), \end{gather}

with $c=a\sinh \alpha , \gamma _n=n+1/2, t=\cos \eta , P_n$ the usual Legendre polynomial of order $n$ and primes designating differentiations. In (3.25)–(3.28) the superscript $l$ that is either $l=1$ or $l=2$ refers to either the ambient linear shear flow or quadratic shear flow, respectively. The dimensionless coefficients $A_n, B_n, C_n, D_n, E_n, F_n$ and $G_n$ in (3.25)–(3.28) are required to vanish as $n$ becomes infinite. They are obtained by enforcing the condition ${{\nabla }}\boldsymbol {\cdot }{\boldsymbol {u}}'=0$ in the $0\leq \xi \leq \alpha$ liquid domain, the Navier slip and impermeability conditions on the $\xi =z=0$ slip wall and the no-slip condition ${\boldsymbol {u}}'=-{\boldsymbol {u}}_a$ on the $\xi =\alpha$ sphere surface (see details in Loussaief et al. (Reference Loussaief, Pasol and Feuillebois2015) and Ghalia et al. (Reference Ghalia, Feuillebois and Sellier2016)).

Recall that in (3.23) the partial derivatives with respect to $r$ are to be performed at constant angle $\theta .$ Then (see also Dean & O'Neill Reference Dean and O'Neill1963), on the $\xi =\alpha$ sphere boundary

(3.29a,b)\begin{align} \left({\frac{\partial f}{\partial r}}\right) & ={\frac{t-\cosh\alpha}{c}} \left({\frac{\partial f}{\partial \xi}}\right), \quad \textrm{d}\theta={\frac{\sin\theta \,\textrm{d}\eta}{\sin\eta}}, \end{align}

(3.29c,d)\begin{align} \sin\theta & ={\frac{\sinh\alpha\sin\eta}{\cosh\alpha-t}},\quad \cos\theta={\frac{t\cosh\alpha-1}{\cosh\alpha-t}}, \end{align}

with $t=\cos \eta$ and $\partial f/\partial \xi$ to be taken at constant value of $\eta$ for $f(r,\theta )=f(\xi ,\eta )$. Appealing to (3.29a–d) when using (3.22)–(3.23) then finally yields

(3.30)\begin{align} Q^{d,2}_{xzz}/\mu_0=\int_S a^{2} T_2[\cos\theta]^{2} \,\textrm{d}S=-{\frac{{\rm \pi} k_l a^{2} c^{l+1}} {4}}\{S_A+S_B+S_C+S_D+S_E+S_F+S_G\}, \end{align}

where $S_L$ designates a sum involving the coefficients $L$ introduced in (3.25)–(3.28). Each sum $S_L$ is expressed in terms of 10 integrals $I_m(n,\alpha )$ in appendix B. For example,

(3.31)\begin{equation} S_A=\sinh\alpha\sum_{n\geq 1}\left\{\sinh\alpha[1-\textrm{e}^{-2\gamma_n\alpha}]I_4(n,\alpha) +2\gamma_n[1+\textrm{e}^{-2\gamma_n\alpha}]I_3(n,\alpha)\right\}A_n \end{equation}

with the definitions

(3.32)\begin{equation} I_m(n,\alpha)=\textrm{e}^{\gamma_n\alpha}\int_{-1}^{1}\frac{(1-t^{2})[(\cosh\alpha)t-1]^{3} P'_{n}(t)\,\textrm{d}t}{(\cosh\alpha-t)^{m+3/2}}\quad \textrm{for}\ m=3,4. \end{equation}

Appendix B both defines and explains how to calculate analytically all 10 integrals $I_m$. Since they are too lengthy to be reproduced in this article, formulae are displayed only for $I_3$ and $I_4$ in appendix B while the analytical results for the other integrals $I_m$ are given in the supplementary material available at https://doi.org/10.1017/jfm.2020.429. Finally, combining (3.17) with (3.30) provides the quadrupole $Q_{xzz}=Q^{a}_{xzz}+Q^{d,1}_{xzz}+Q^{d,2}_{xzz}$ for the ambient linear and quadratic shear flows. Using (3.13) and (3.14) then gives the desired coefficients $q^{s}_{xzz}$ and $q^{q}_{xzz}$.

3.4. Accuracy issues and comparisons against a BEM approach

The BCM is known to be very accurate provided its numerical implementation is adequately handled. The encountered series, as in (3.25)–(3.28), are truncated beyond a given integer $N_t$. Of course, the value of $N_t$ required to ensure a given accuracy level depends on $(z_z/a,{\tilde{\lambda} })$. Moreover, while coding in Fortran double precision $(D)$ was sufficient for some domains of parameters (see Ghalia et al. Reference Ghalia, Feuillebois and Sellier2016), it has been found here necessary to code in Fortran quadruple precision $(Q)$ (as in Loussaief et al. Reference Loussaief, Pasol and Feuillebois2015) to accurately compute the new terms $q^{s}_{xzz}$ or $q^{q}_{xzz}$ when $\tilde{z}_{c}$ becomes large whatever ${\tilde{\lambda} }\geq 0.$ Moreover, the quadruple precision code uses another numerical method, i.e. the Thomas algorithm (as in Loussaief et al. Reference Loussaief, Pasol and Feuillebois2015), instead of the classical LU factorization (as in Ghalia et al. Reference Ghalia, Feuillebois and Sellier2016) to invert the tridiagonal matrix providing the series coefficients. The interest of the Thomas algorithm is that it allows the calculation of more terms in the series at a reduced cost in memory and CPU calculation time. Displaying results using these two different methods also provides here a supplementary check for the validity of the numerical results. The comparison is illustrated for $q^{q}_{xzz}$ in table 2 taking $z_z/a=1.005, 1.5, 100$ and ${\tilde{\lambda} }=0,1,10$. Clearly, for $z_z/a=100$ the calculations in quadruple precision are more accurate and efficient (smaller required value of $N_t$) whatever ${\tilde{\lambda} }.$ For instance, taking $N_t=2000$ in quadruple precision is sufficient to reach an eight-digit accuracy when $z_z/a=100$ even for the large slip case ${\tilde{\lambda} }=10$. When the sphere is very close to a wall with large slip length it is necessary to take $N_t$ large for computations in both double and quadruple precision.

Table 2. Normalized quadrupole $q^{q}_{xzz}$ computed with the BCM using either Fortran double precision $(D)$ or quadruple precision with the Thomas algorithm $(Q)$ and truncating the series in (3.25)–(3.28) beyond $N_t$. For example, the column $(D,2000)$ gathers the values of $q^{q}_{xzz}$ computed in Fortran double precision taking $N_t=2000$.

In view of the great amount of calculations (see the previous subsection, appendix B and the supplementary material), it is useful for validation to compare the BCM predictions, for both dimensionless quadrupoles $q^{s}_{xzz}$ and $q^{q}_{xzz},$ with the ones obtained by the quite different BEM implemented in Ghalya, Sellier & Feuillebois (Reference Ghalya, Sellier and Feuillebois2012) for a solid arbitrary-shaped particle interacting with the $z=0$ slip wall. This method uses the Green tensor fulfilling the boundary conditions on the wall and implements six-node (P2 quadratic approximation) isoparametric curved triangular meshes on the particle surface. More details can also be found in Sellier (Reference Sellier, Aliabadi and Wen2010). Some numerical comparisons for $q^{s}_{xzz}$ and $q^{q}_{xzz}$ are displayed in table 3 for a few settings $(z_z/a,{\tilde{\lambda} })$ and taking $N_t=20\,000$ for the BCM (with quadruple precision) and $1058$ nodal points on the sphere boundary for the BEM. Since the BEM requires much more refined meshes on the sphere surface as $\tilde{z}_{c}-1$ vanishes the results in table 3 are given for $\tilde{z}_{c}\geq 1.1.$ As observed in table 3, a very good agreement (up to a four-digit accuracy) is found between the BEM and the quadruple precision BCM (with the Thomas algorithm). Finally, comparisons between the BCM and BEM techniques are given in table 4 for the key quantities $\tilde{S}_{xz}$ and $\tilde{Q}_{xzz}$ taking $H=5a.$ Again, a good agreement is observed.

Table 3. Computed normalized quadrupoles $q^{s}_{xzz}$ and $q^{q}_{xzz}$ using either the BCM with $N_t=20\,000$ or the BEM spreading $1058$ collocation nodal points on the sphere surface.

Table 4. Computed normalized quantities $\tilde{S}_{xz}$ and $\tilde{Q}_{xzz}$ using the BEM ($\star$ symbol) or the BCM taking either $N_t=2000$ (a) or $N_t=20\,000$ (b). Here $H=5a$ and $1058$ collocation nodal points are distributed on the sphere surface for the BEM.

4.1. Comparison and validation for no-slip walls

The intrinsic viscosity for no-slip walls was calculated in (I) using the accurate multipole method which fully takes into account the interactions with both walls. That method was corrected for lubrication to speed up the convergence rate. Plots in (I) were in terms of the dimensionless ratio $H/d= \tilde{H}/2$. We use also $H/d$ in this article for comparison with (I) and also for all further plots. This choice is also motivated by possible comparisons with experiments in which the particle diameter is provided rather than the radius. A comparison was made in (I) with a ‘simplified model’ which took into account only the stresslet term and in this case only the nearest wall. That comparison clearly showed that both models match for $H/d\geq 10$ while the behaviour is quite different for $H/d$ down to unity. Indeed, the intrinsic viscosity predicted by the ‘simplified model’ drops down to zero for $H/d\to 1$.

The calculation of the quadrupole term in bipolar coordinates, presented in § 3.3, will be used here. First, the full nearest-wall model nw which takes into account both the stresslet and the quadrupole is employed. The comparison is shown in figure 3. It is observed that the curves for the two-walls exact model (say tw) and the nw model have the same shape for small $H/d,$ which is not the case when only stresslet is taken into account (compare with figure 2 in (I)). It is also observed that the value of $[\mu ]$ for the nw model is smaller than that for the two-walls model. This is expected, since some friction effects from both walls at the same time are discarded by the nw model, especially in the middle of the gap. For comparison, figure 3 also displays the results for the no-wall (0w) approximation already shown in (I).

Figure 3. Comparison and validation for no-slip walls. Results for two walls (tw) are from (I) and were calculated for two parallel walls with the multipole method. The wall superposition (ws, dash-dotted line, red) and nearest wall (nw, lines with larger dots, blue) results are obtained with the present bipolar coordinates calculations. The no wall (0w, lines with circles) curves are obtained from the analytical equations (3.1a,b). Dashed lines for the various models (see legend) are for stresslet only. The horizontal dotted line at $[\mu ]=2.5$ recalls the result of Einstein (Reference Einstein1906) for the intrinsic viscosity in unbounded fluid.

Results for the wall-superposition (ws) model, taking into account both the stresslet and quadrupole in the interaction with each wall separately, are also plotted in figure 3. Curves for the nw and ws models are remarkably close together in the range $H/d\geq 1.2$.

The ability of the nw and ws models to predict the effective viscosity is quantified by considering their relative ‘error’ $({\rm \Delta} [\mu ])/[\mu ]$ by comparison with the exact tw model. For instance for the nw model, $({\rm \Delta} [\mu ])/[\mu ]= |1-[\mu ]_{nw}/[\mu ]_{tw}|$. These errors are plotted in figure 4. They are small and of comparable magnitude, the ws approximation being a little better that the nw one in the range $1.5 \leq H/d \leq 6$. In the same figure, we plot the horizontal dashed line which indicates the 5 % error level. It is remarkable that for both nw and ws approximations the relative error remains below that level in the wide range $H/d\geq 2$. For $H/d\geq 3$, both relative errors are even less than 2 %.

Figure 4. Comparison and validation for no-slip walls: relative differences between the nw and ws models and the exact two-walls (tw) results of Feuillebois et al. (Reference Feuillebois, Ekiel-Jeżewska, Wajnryb, Sellier and Bławzdziewicz2016), i.e. $|1-[\mu ]_ \textit {nw}/[\mu ]_ \textit {tw}|$ and $|1-[\mu ]_ \textit {ws}/[\mu ]_ \textit {tw}|$.

In conclusion, both nw and ws models are validated in the range $H/d\geq 2$, this lower bound being represented as a vertical dash-dotted line in figure 4.

It is expected that the influences of walls decrease with increasing normalized slip length $\tilde{\lambda}$. Consequently, for $\tilde{\lambda}>0$ the accuracy of both nw and ws models for slip walls is likely to be even better than 5 % in the whole range of distances $H/d\geq 2$.

The counterpart to restricting the domain to $H/d\geq 2$ because of accuracy issues is that the interesting increase of $[\mu ]$ for small gaps due to lubrication (see figure 3) will not be reconsidered for slip walls. This is unavoidable in this article, since considering interactions with both walls at the same time would be compulsory for the range of small gaps $H/d< 2$; that is, using only lubrication with individual walls for $H/d$ close to unity would not be sufficient.

To make this idea more precise, consider first no-slip walls. For $H/d$ close to unity, lubrication interactions with separate walls are of order $\log (H/d-1)$ and the next order term for small gaps (numerically very close to the first one) would include interactions with both walls, as pointed out in Ekiel-Jeżewska et al. (Reference Ekiel-Jeżewska, Wajnryb, Bawzdziewicz and Feuillebois2008). Using the precise multipole method (see Bhattacharya et al. Reference Bhattacharya, Bławzdziewicz and Wajnryb2005a,Reference Bhattacharya, Bławzdziewicz and Wajnrybb), it was shown there that the second-order lubrication terms (constant and of the same order as the leading logarithmic terms) are essentially different from the superposition of the corresponding single-wall terms.

This problem carries over to slip walls. When calculating the force on a sphere near a single slip wall, Davis, Kezirian & Brenner (Reference Davis, Kezirian and Brenner1994) found a weak log-singularity in the lubrication regime at low slip and no singularity for other values of the slip. A similar behaviour is expected for the stresslet and quadrupole. These quantities are actually calculated in this lubrication regime from our accurate bipolar coordinates results. As for the next order lubrication terms, they are expected to be numerically very close to the first-order ones and to involve both slip walls. These terms cannot be calculated at the present time.

4.2. Integrands in the expressions of the effective viscosity for slip walls

This subsection is concerned with the variations of the integrands ${\mathcal {J}}_ \textit {0w}, {\mathcal {J}}_ \textit {nw}$ and ${\mathcal {J}}_ \textit {ws}$ (defined in § 3.1) in the expressions for the intrinsic viscosity and is motivated by twoissues:

1. (i) It is expected on physical grounds that the regions close to the slip walls would have a larger contribution than the region in the bulk of the channel and that these contributions would be more sensitive to variations in the slip length $\tilde{\lambda}$.

2. (ii) It was noticed for the no-slip case in figure 3 that the quadrupole term significantly contributes to the shape of the curve for the effective viscosity in the range $H/d \leq 5$. This is an incentive for studying the relative contribution of the quadrupole to the integrand for the case of slip walls.

The integrands ${\mathcal {J}}_ \textit {0w}$, ${\mathcal {J}}_ \textit {nw}$ and ${\mathcal {J}}_ \textit {ws}$ depend on $H/d, \tilde{\lambda}$ and $z_c/H$. We consider first the very confined case $H/d=2$. The integrands are plotted together versus $z_c/H$ in figure 5(a) for $\tilde{\lambda}=0, 0.3, 1$. It is observed that they are positive, symmetric with respect to the channel mid-plane and monotonic over each half-channel. Such a smooth behaviour allows calculation of the integral for the intrinsic viscosity by a standard Gauss–Legendre quadrature with 16 points on a half-channel. The smallest Gauss–Legendre integration point is at $z_c/H=0.2513$, not too close to the contact with wall at $z_c/H=0.2500$, i.e. higher than the lowest point $\tilde{z}_{c}=1.005$ considered in table 2, so that using $N_t\leq 20\,000$ in the bipolar coordinates calculations is sufficient for accuracy.

Figure 5. (a) Integrands ${\mathcal {J}}_ \textit {0w}$, ${\mathcal {J}}_ \textit {nw}$ and ${\mathcal {J}}_ \textit {ws}$ (represented by circles, dots and solid lines, respectively) for the very confined case $H/d=2$. The sphere is in contact with a wall at $z_c/H=0.25$ and $z_c/H=0.75$. The integrands ${\mathcal {J}}_ \textit {nw}$ and ${\mathcal {J}}_ \textit {ws}$ look practically superimposed. Values of the normalized slip length $\tilde{\lambda}=0, 0.3$ and $1$ correspond to the upper, middle and lower solid lines (black, green and red), respectively. (b) Relative contribution of the quadrupole ${\mathcal {J}}_Q/ {\mathcal {J}}$ for the cases $\textit {0w}$, $\textit {nw}$ and $\textit {ws}$, with the same notation as in (a). (c) Zoom of the preceding plots near the wall, with the same notation as in (a).

Note that this integration procedure was already applied for the no-slip $\tilde{\lambda}=0$ case in § 4.1.

The zero-wall integrand ${\mathcal {J}}_ \textit {0w}$ (see figure 5a) is independent of $\tilde{\lambda}$ and always smaller than the other integrands. As expected, the contribution of the wall regions is larger than that of the bulk (say, of the domain $| \tilde{z}_{c}/H -0.5| \leq 0.15$). In the bulk region, each integrand is practically insensitive to $\tilde{\lambda}$. By contrast, ${\mathcal {J}}_ \textit {nw}$ and ${\mathcal {J}}_ \textit {ws}$ significantly decrease with increasing $\tilde{\lambda}$ in the wall regions. Moreover, ${\mathcal {J}}_ \textit {nw}$ and ${\mathcal {J}}_ \textit {ws}$ are very close together in the whole domain for all values of $\tilde{\lambda}$.

The relative contribution of the quadrupole term ${\mathcal {J}}_Q/ {\mathcal {J}}$ (recall (2.26b,c)) is displayed in figure 5(b), still for $H/d=2$ and $\tilde{\lambda}=0, 0.3, 1$. It is larger than 15% and nearly insensitive to $\tilde{\lambda}$ in the bulk region. It varies slightly with $\tilde{\lambda}$ in the wall regions, as shown by the zoom in figure 5(c).

The influence of increasing the normalized channel width $H/d$ is investigated in figure 6 for $H/d=5$ (moderately confined case) and in figure 7 for $H/d=50$ (widely separated walls), for comparison with figure 5 for $H/d=2$. The bulk region referred to above increases with $H/d$. Near the walls, the sensitivity of ${\mathcal {J}}$ to $\tilde{\lambda}$ is about the same for all values of $H/d$. The quadrupole part ${\mathcal {J}}_Q/ {\mathcal {J}}$ decays significantly as $H/d$ increases, except in the very vicinity of the channel mid-plane. Note that ${\mathcal {J}}_Q/ {\mathcal {J}}=1$ at mid-channel, i.e. $z_c/H= \tilde{z}_{c}/ \tilde{H}=1/2$, because the local shear rate $1-2 \tilde{z}_{c}/ \tilde{H}$ and thus also ${\mathcal {J}}_S$, see (2.26b,c), vanishes there. On the other end, as $H/d$ increases, the sensitivity of ${\mathcal {J}}_Q/ {\mathcal {J}}$ to $\tilde{\lambda}$ is enhanced near walls.

Figure 6. (a) Variation of the integrand ${\mathcal {J}}_ \textit {0w}, {\mathcal {J}}_ \textit {nw}, {\mathcal {J}}_ \textit {ws}$. (b) Zoom of (a) near the wall at $z_c=0$. Curves are for the moderately confined case $H/d=5$. The sphere is in contact with a wall at $z_c/H=0.1$ and $z_c/H=0.9$. Values of the normalized slip length $\tilde{\lambda}=0$, 0.3 and 1 correspond to the upper, middle and lower solid lines (black, green and red), respectively. (c)Ratios ${\mathcal {J}}_{ \textit {0w},Q}/ {\mathcal {J}}_{ \textit {0w}}, {\mathcal {J}}_{ \textit {nw},Q}/ {\mathcal {J}}_{ \textit {nw}}, {\mathcal {J}}_{ \textit {ws},Q}/ {\mathcal {J}}_{ \textit {ws}}$. (d) Zoom of (c) near the wall at $z_c=0$.

Figure 7. (a) Variation of the integrand ${\mathcal {J}}_ \textit {0w}, {\mathcal {J}}_ \textit {nw}, {\mathcal {J}}_ \textit {ws}$. (b) Zoom of (a) near the wall at $z_c=0$. Curves are for widely separated walls, $H/d=50$. The sphere is in contact with a wall at $z_c/H=0.01$ and $z_c/H=0.99$. Values of the normalized slip length $\tilde{\lambda}=0$, 0.3 and 1 correspond to the upper, middle and lower solid lines (black, green and red), respectively. (c)Ratios ${\mathcal {J}}_{ \textit {0w},Q}/ {\mathcal {J}}_{ \textit {0w}}, {\mathcal {J}}_{ \textit {nw},Q}/ {\mathcal {J}}_{ \textit {nw}}, {\mathcal {J}}_{ \textit {ws},Q}/ {\mathcal {J}}_{ \textit {ws}}$. (d) Zoom of (c) near the wall at $z_c=0$. (e) Zoom of the peak near the channel mid-plane.

4.4. An interpolation formula for the effective viscosity

The preceding numerical results were obtained from a lengthy quadruple precision Fortran code enforcing the bipolar coordinates calculations. Building this code required substantial efforts. For a practical prediction of the effective viscosity, we therefore propose a handy and accurate formula interpolating our results obtained with the ws model.

Considering that $[\mu ]$ has the Einstein limit of $5/2$ at $H/d\to \infty$ for any $\tilde{\lambda}$, it is convenient to fit, in the sense of least squares, the numerical results for $[\mu ]-5/2$ with a polynomial in $d/H$ and $\tilde{\lambda}$ as follows:

(4.2)\begin{equation} [\mu]_{\textit{fit}} = {\frac{5}{2} }-\sum_{i=0}^{5} \sum_{j=0}^{8} c_{i,j} \left( \frac dH \right)^{5-i} \tilde{\lambda}^{8-j} . \end{equation}

The selected upper bounds $(5$ for $i$ and $8$ for $j$) were chosen as a compromise between a reasonable number of coefficients $c_{i,j}$ and a sufficient relative accuracy of at most $0.5\,\%$ (see later in figure 15). This value is consistent with expected accuracy of the $\textit {nw}$ and $\textit {ws}$ approximations.

We used (4.2) for fitting the numerical results for $[\mu ],$ given by the ws model, up to $H/d=100$ and added, from physical considerations, the value $[\mu ]=5/2$ at $d/H=0$ for all values of $\tilde{\lambda}$. Moreover, even though the range of slip lengths $\tilde{\lambda}\leq 1$ was kept for physical reasons (see above), calculated values of $[\mu ]$ in the range $1 < \tilde{\lambda}\leq 2$ were used to improve the fitting precision of the polynomial in the range $0\leq \tilde{\lambda}\leq 1$. The coefficients $c_{i,j}$ of the fitting polynomial are displayed in table 6.

Table 6. Coefficients $c_{i,j}$ of the polynomial in formula (4.2) fitting the numerical results obtained with the ws model.

Results for $[\mu ]$ using the fitting formula (4.2) are displayed as solid lines and compared to the exact results displayed as crosses in figure 13. Values for $[\mu ]$ are also shown as a contour plot in figure 14(a) and a zoom in the range $2\leq H/d \leq 10$ is provided in figure 14(b). Values of the relative error between the fitted and exact results are displayed in figure 15. Only values $H/d\leq 5$ are represented in this latter figure in order to show details of the relative error contour lines (for $H/d$ from $5$ to $100$ the contour lines become nearly parallel to the $\tilde{\lambda}=0$ vertical axis). It is remarked that, as discussed after (4.2), the error is mostly smaller than 0.4 %, except for small $\tilde{\lambda}$ and $H/d$ close to 2 where it is less than 0.5 % (this contour line is very close to $H/d=2$ in figure 15). For $H/d>3.5$, the error is smaller than 0.2 % for all values of $\tilde{\lambda}$.

Figure 13. Variations of the intrinsic viscosity $[\mu ]$ (as predicted with the ws model) with $H/d$ from 2 to 100, for (from top to bottom): $\tilde{\lambda}=0, 0.01, 0.05$ and $0.1$ to $1$ in steps of $0.1$. Crosses: exact results; solid lines: fitting formula (4.2).

Figure 14. (a) Contour plot of the variations of the intrinsic viscosity $[\mu ]$ with $\tilde{\lambda}$ and $H/d$ extracted from the exact results displayed in figure 13 (crosses). (b) Zoom plot for $2\leq H/d\leq 10$.

Figure 15. Contour plot of the relative error between the exact results and the fitted formula for $[\mu ]$.

Finally, note that for the particular case $\tilde{\lambda}=0$, the fitting formula (4.2) complements the results of (I) in the case of a suspension of spheres.