2 Mobility and friction matrices

We consider $N$ spherical particles in an incompressible fluid of viscosity $\unicode[STIX]{x1D702}_{0}$ at a low Reynolds number, immersed in an external fluid flow $\boldsymbol{v}_{\infty }(\boldsymbol{r})$ , with stick boundary conditions on the particle surfaces. Linearity of the Stokes equations implies a linear relation (Brenner & O’Neill Reference Brenner and O’Neill1972; Cichocki, Felderhof & Schmitz Reference Cichocki, Felderhof and Schmitz1988; Kim & Karrila Reference Kim and Karrila1991)

(2.1) $$\begin{eqnarray}\left(\begin{array}{@{}c@{}}\tilde{\boldsymbol{F}}\\ \tilde{\boldsymbol{T}}\\ \tilde{\boldsymbol{S}}\end{array}\right)=-\left(\begin{array}{@{}ccc@{}}\unicode[STIX]{x1D73B}^{tt} & \unicode[STIX]{x1D73B}^{tr} & \unicode[STIX]{x1D73B}^{td}\\ \unicode[STIX]{x1D73B}^{rt} & \unicode[STIX]{x1D73B}^{rr} & \unicode[STIX]{x1D73B}^{rd}\\ \unicode[STIX]{x1D73B}^{dt} & \unicode[STIX]{x1D73B}^{dr} & \unicode[STIX]{x1D73B}^{dd}\end{array}\right)\boldsymbol{\cdot }\left(\begin{array}{@{}c@{}}\tilde{\boldsymbol{v}}_{\infty }-\tilde{\boldsymbol{U}}\\ \tilde{\unicode[STIX]{x1D74E}}_{\infty }-\tilde{\unicode[STIX]{x1D734}}\\ \tilde{\boldsymbol{E}}_{\infty }\end{array}\right),\end{eqnarray}$$

which defines the grand friction matrix $\unicode[STIX]{x1D73B}$ . The elements $\unicode[STIX]{x1D73B}^{pq}$ (with $p,q=t,r,d$ ) are the Cartesian tensors and the superscripts $t,r,d$ denote translational, rotational and dipolar components, respectively. Here, $\tilde{\boldsymbol{F}}=(\boldsymbol{F}_{1},\boldsymbol{F}_{2},\ldots ,\boldsymbol{F}_{N})$ are the $3N$ -dimensional vectors of forces with which particles act on the fluid and analogously for the torques $\tilde{\boldsymbol{T}}$ , as well as translational and rotational velocities of the particles ( $\tilde{\boldsymbol{U}}$ and $\tilde{\unicode[STIX]{x1D734}}$ ). Next, $\tilde{\boldsymbol{v}}_{\infty }=(\boldsymbol{v}_{\infty }(\boldsymbol{R}_{1}),\ldots ,\boldsymbol{v}_{\infty }(\boldsymbol{R}_{N}))$ are the values of external flow velocity calculated at the centres of the particles, $\unicode[STIX]{x1D619}_{i}$ . Similarly, $\tilde{\unicode[STIX]{x1D74E}}_{\infty }$ gives the vector of vorticities at the centres of the particles, with $\unicode[STIX]{x1D74E}_{\infty }=(1/2)\unicode[STIX]{x1D735}\times \boldsymbol{v}_{\infty }$ . Finally, $\tilde{\boldsymbol{E}}_{\infty }=(\boldsymbol{E}_{\infty }(\boldsymbol{R}_{1}),\ldots ,\boldsymbol{E}_{\infty }(\boldsymbol{R}_{N}))$ is the vector of strain rates, with $\boldsymbol{E}_{\infty }=(1/2)(\unicode[STIX]{x1D735}\boldsymbol{v}_{\infty }+(\unicode[STIX]{x1D735}\boldsymbol{v}_{\infty })^{\text{T}})$ whereas $\tilde{\boldsymbol{S}}=(\boldsymbol{S}_{1},\ldots ,\boldsymbol{S}_{N})$ are the particle stresses. Both are $5N$ -dimensional, due to the symmetric and traceless character of strain and stress tensors.

The grand mobility matrix $\unicode[STIX]{x1D741}$ is defined by the relation

(2.2) $$\begin{eqnarray}\left(\begin{array}{@{}c@{}}\tilde{\boldsymbol{U}}-\tilde{\boldsymbol{v}}_{\infty }\\ \tilde{\unicode[STIX]{x1D734}}-\tilde{\unicode[STIX]{x1D74E}}_{\infty }\\ -\tilde{\boldsymbol{S}}\end{array}\right)=\left(\begin{array}{@{}ccc@{}}\unicode[STIX]{x1D741}^{tt} & \unicode[STIX]{x1D741}^{tr} & \unicode[STIX]{x1D741}^{td}\\ \unicode[STIX]{x1D741}^{rt} & \unicode[STIX]{x1D741}^{rr} & \unicode[STIX]{x1D741}^{td}\\ \unicode[STIX]{x1D741}^{dt} & \unicode[STIX]{x1D741}^{dr} & \unicode[STIX]{x1D741}^{dd}\end{array}\right)\boldsymbol{\cdot }\left(\begin{array}{@{}c@{}}\tilde{\boldsymbol{F}}\\ \tilde{\boldsymbol{T}}\\ \tilde{\boldsymbol{E}}_{\infty }\end{array}\right),\end{eqnarray}$$

which is a partial inverse of the relation (2.1). The two upper rows of the above relation describe the dynamics of the system. The third row (stresslet) is needed to compute the stress tensor in the suspension; however, it has no bearing on the movement of particles.

Another important hydrodynamic tensor is the full inverse of the grand friction matrix, introduced by Durlofsky et al. (Reference Durlofsky, Brady and Bossis1987), $\unicode[STIX]{x1D662}=\unicode[STIX]{x1D73B}^{-1}$ which links force and velocity multipoles

(2.3) $$\begin{eqnarray}\left(\begin{array}{@{}c@{}}\tilde{\boldsymbol{v}}_{\infty }-\tilde{\boldsymbol{U}}\\ \tilde{\unicode[STIX]{x1D74E}}_{\infty }-\tilde{\unicode[STIX]{x1D734}}\\ \tilde{\boldsymbol{E}}_{\infty }\end{array}\right)=-\left(\begin{array}{@{}ccc@{}}\unicode[STIX]{x1D662}^{tt} & \unicode[STIX]{x1D662}^{tr} & \unicode[STIX]{x1D662}^{td}\\ \unicode[STIX]{x1D662}^{rt} & \unicode[STIX]{x1D662}^{rr} & \unicode[STIX]{x1D662}^{td}\\ \unicode[STIX]{x1D662}^{dt} & \unicode[STIX]{x1D662}^{dr} & \unicode[STIX]{x1D662}^{dd}\end{array}\right)\boldsymbol{\cdot }\left(\begin{array}{@{}c@{}}\tilde{\boldsymbol{F}}\\ \tilde{\boldsymbol{T}}\\ \tilde{\boldsymbol{S}}\end{array}\right).\end{eqnarray}$$

The inverse friction matrix will play a central role in the construction of the generalized Rotne–Prager–Yamakawa approximation (GRPY), as detailed in the next section.

3 Generalized Rotne–Prager–Yamakawa approximation

The relation between the velocities of particles moving in a Stokes flow and the induced force density localized on particle surfaces ${\mathcal{S}}_{i},i=1,\ldots ,n$ , can be written as

(3.1) $$\begin{eqnarray}[\boldsymbol{U}_{i}+\unicode[STIX]{x1D734}_{i}\times \unicode[STIX]{x1D746}_{i}-\boldsymbol{v}_{\infty }(\boldsymbol{r})]_{\boldsymbol{r}\in {\mathcal{S}}_{i}}=\mathop{\sum }_{j}\int \unicode[STIX]{x1D64F}_{0}(\boldsymbol{r}-\boldsymbol{r}^{\prime })\boldsymbol{\cdot }\boldsymbol{f}_{j}(\boldsymbol{r}^{\prime })\,\text{d}^{3}\boldsymbol{r}^{\prime }.\end{eqnarray}$$

Here $\unicode[STIX]{x1D746}_{i}=\boldsymbol{r}-\boldsymbol{R}_{i}$ with $\boldsymbol{R}_{i}$ denoting the position of particle $i$ . On the other hand, $\boldsymbol{f}_{j}(\boldsymbol{r}^{\prime })$ is the density of the forces with which the particle $j$ is acting on the fluid. Finally, $\unicode[STIX]{x1D64F}_{0}$ is the Oseen tensor

(3.2) $$\begin{eqnarray}\unicode[STIX]{x1D64F}_{0}(\boldsymbol{r})=\frac{1}{8\unicode[STIX]{x03C0}\unicode[STIX]{x1D702}_{0}r}(\mathbf{1}+\hat{\boldsymbol{r}}\hat{\boldsymbol{r}}),\end{eqnarray}$$

which is the Green function for the Stokes problem for the unbounded space. In the above, $r$ is the length of the vector $\boldsymbol{r}$ and $\hat{\boldsymbol{r}}=\boldsymbol{r}/r$ . The key idea behind the GRPY approximation is to describe the force density in terms of its three multipoles only:

(3.3) $$\begin{eqnarray}\boldsymbol{f}_{j}=\boldsymbol{w}_{j}^{t}\boldsymbol{\cdot }\boldsymbol{F}_{j}+\boldsymbol{w}_{j}^{r}\boldsymbol{\cdot }\unicode[STIX]{x1D64F}_{j}+\boldsymbol{w}_{j}^{d}\boldsymbol{ : }\boldsymbol{S}_{j},\end{eqnarray}$$

where $\boldsymbol{w}_{j}^{p}$ are operators associated with different multipoles ( $p=t,r,d$ )

(3.4) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\boldsymbol{w}_{j}^{t}(\boldsymbol{r})={\displaystyle \frac{1}{4\unicode[STIX]{x03C0}a_{j}^{2}}}\mathbf{1}\unicode[STIX]{x1D6FF}(\unicode[STIX]{x1D70C}_{j}-a_{j}),\quad \boldsymbol{w}_{j}^{r}(\boldsymbol{r})={\displaystyle \frac{3}{8\unicode[STIX]{x03C0}a_{j}^{3}}}\unicode[STIX]{x1D750}\boldsymbol{\cdot }\hat{\unicode[STIX]{x1D746}}_{j}\unicode[STIX]{x1D6FF}(\unicode[STIX]{x1D70C}_{j}-a_{j})\\ \boldsymbol{w}_{j}^{d}(\boldsymbol{r})={\displaystyle \frac{3}{4\unicode[STIX]{x03C0}a_{j}^{3}}}\unicode[STIX]{x1D6FF}(\unicode[STIX]{x1D70C}_{j}-a_{j})\hat{\unicode[STIX]{x1D746}}_{j}\boldsymbol{\cdot }\boldsymbol{{\mathcal{I}}}.\end{array}\right\}\end{eqnarray}$$

In the above $[\unicode[STIX]{x1D750}\boldsymbol{\cdot }\hat{\unicode[STIX]{x1D746}}]_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}=\unicode[STIX]{x1D716}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FE}}[\hat{\unicode[STIX]{x1D746}}_{j}]_{\unicode[STIX]{x1D6FE}}$ and $\boldsymbol{{\mathcal{I}}}$ , is the fourth-rank isotropic tensor, traceless and symmetric in its first and last index pairs:

(3.5) $$\begin{eqnarray}{\mathcal{I}}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D6FE}}={\textstyle \frac{1}{2}}(\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FF}}\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FE}}+\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FE}}\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FF}}-{\textstyle \frac{2}{3}}\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D6FE}}).\end{eqnarray}$$

Next, the velocity field in the left-hand side of (3.1) is approximated by a linear flow. The multipoles characterizing this flow: velocity, $\boldsymbol{U}_{i}-\boldsymbol{v}_{\infty }(\boldsymbol{R}_{i})$ , vorticity, $\unicode[STIX]{x1D734}_{i}-\unicode[STIX]{x1D74E}_{\infty }(\boldsymbol{R}_{i})$ , and strain rate, $-\boldsymbol{E}_{\infty }(\boldsymbol{R}_{i})$ , can be obtained from $[\boldsymbol{U}_{i}+\unicode[STIX]{x1D734}_{i}\times \unicode[STIX]{x1D746}_{i}-\boldsymbol{v}_{\infty }(\boldsymbol{r})]_{\boldsymbol{r}\in {\mathcal{S}}_{i}}$ by integration of this expression multiplied by the transposed $\boldsymbol{w}_{i}^{p}$ operators. This leads to the following relation for the elements of the inverse friction matrix

(3.6) $$\begin{eqnarray}\unicode[STIX]{x1D662}_{ij}^{\,pq}=\langle \boldsymbol{w}_{i}^{p}|\unicode[STIX]{x1D64F}_{0}|\boldsymbol{w}_{j}^{q}\rangle =\iint [\boldsymbol{w}_{i}^{p}(\boldsymbol{r}^{\prime })]^{\text{T}}\boldsymbol{\cdot }\unicode[STIX]{x1D64F}_{0}(\boldsymbol{r}^{\prime }-\boldsymbol{r}^{\prime \prime })\boldsymbol{\cdot }\boldsymbol{w}_{j}^{q}(\boldsymbol{r}^{\prime \prime })\,\text{d}\boldsymbol{r}^{\prime }\,\text{d}\boldsymbol{r}^{\prime \prime },\end{eqnarray}$$

where we have used the relation (2.3) together with an approximation (3.3). In the above, T stands for the transpose. Equation (3.6) can be applied to other geometries (periodic boundary conditions, presence of a wall etc.) by substituting $\unicode[STIX]{x1D64F}_{0}$ with a relevant Green’s function for a particular geometry, as outlined in Wajnryb et al. (Reference Wajnryb, Mizerski, Zuk and Szymczak2013).

A similar construction of the GRPY approximation on the level of the inverse friction matrix has been carried out by Durlofsky et al. (Reference Durlofsky, Brady and Bossis1987) and Brady & Bossis (Reference Brady and Bossis1988); however, their focus was on the suspension of freely moving, non-overlapping particles. As a result, they combined the RPY approach with the lubrication corrections which come into play as the particle surfaces approach one another – an effect not present in the rigidly connected beads in the macromolecular models.

Note that the translational and rotational components of the inverse friction matrix $\unicode[STIX]{x1D662}_{ij}^{\,pq}$ (with $p,q=t,r$ ) are of the same form as the respective elements of the grand mobility matrix $\unicode[STIX]{x1D741}^{ij}$ in the Rotne–Prager approximation. The explicit formulae for these components are not given here, since they can be found in Zuk et al. (Reference Zuk, Wajnryb, Mizerski and Szymczak2014).

We now turn to the calculation of the dipolar elements of the $\unicode[STIX]{x1D662}$ matrix. The diagonal or ‘self’ ( $i=j$ ) elements of the $\unicode[STIX]{x1D662}^{dd}$ component are obtained from (3.6) by performing both integrals over the same surface:

(3.7) $$\begin{eqnarray}\unicode[STIX]{x1D662}_{jj}^{dd}=\frac{3}{20\unicode[STIX]{x03C0}\unicode[STIX]{x1D702}_{0}a_{j}^{3}}\boldsymbol{{\mathcal{I}}}.\end{eqnarray}$$

The off-diagonal ( $i\neq j$ ) components have the following form, stemming from the symmetry of the problem:

(3.8) $$\begin{eqnarray}\unicode[STIX]{x1D662}_{ij}^{dd}(\boldsymbol{R}_{ij})=\frac{3}{40\unicode[STIX]{x03C0}\unicode[STIX]{x1D702}_{0}a_{i}^{3}a_{j}^{3}}[f^{0}(R_{ij})\unicode[STIX]{x1D659}^{0}(\hat{\boldsymbol{R}}_{ij})+f^{1}(R_{ij})\unicode[STIX]{x1D659}^{1}(\hat{\boldsymbol{R}}_{ij})+f^{2}(R_{ij})\unicode[STIX]{x1D659}^{2}(\hat{\boldsymbol{R}}_{ij})],\end{eqnarray}$$

where $\boldsymbol{R}_{ij}=\boldsymbol{R}_{i}-\boldsymbol{R}_{j}$ and $\unicode[STIX]{x1D659}^{i}$ are the following tensors (Kim & Karrila Reference Kim and Karrila1991)

(3.9) $$\begin{eqnarray}\displaystyle & \displaystyle [\unicode[STIX]{x1D659}^{0}(\hat{\boldsymbol{r}})]_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D6FF}}={\textstyle \frac{3}{2}}\left(\hat{r}_{\unicode[STIX]{x1D6FC}}\hat{r}_{\unicode[STIX]{x1D6FD}}-{\textstyle \frac{1}{3}}\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}\right)\left(\hat{r}_{\unicode[STIX]{x1D6FE}}\hat{r}_{\unicode[STIX]{x1D6FF}}-{\textstyle \frac{1}{3}}\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D6FF}}\right), & \displaystyle\end{eqnarray}$$

(3.10) $$\begin{eqnarray}\displaystyle & \displaystyle [\unicode[STIX]{x1D659}^{1}(\hat{\boldsymbol{r}})]_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D6FF}}={\textstyle \frac{1}{2}}\left(\hat{r}_{\unicode[STIX]{x1D6FC}}\hat{r}_{\unicode[STIX]{x1D6FE}}\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FF}}+\hat{r}_{\unicode[STIX]{x1D6FD}}\hat{r}_{\unicode[STIX]{x1D6FE}}\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FF}}+\hat{r}_{\unicode[STIX]{x1D6FC}}\hat{r}_{\unicode[STIX]{x1D6FF}}\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FE}}+\hat{r}_{\unicode[STIX]{x1D6FD}}\hat{r}_{\unicode[STIX]{x1D6FF}}\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FE}}-4\hat{r}_{\unicode[STIX]{x1D6FC}}\hat{r}_{\unicode[STIX]{x1D6FD}}\hat{r}_{\unicode[STIX]{x1D6FE}}\hat{r}_{\unicode[STIX]{x1D6FF}}\right), & \displaystyle\end{eqnarray}$$

(3.11) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D659}^{2}=\boldsymbol{{\mathcal{I}}}-\unicode[STIX]{x1D659}^{0}-\unicode[STIX]{x1D659}^{1}. & \displaystyle\end{eqnarray}$$

For non-overlapping particles, the integral formula (3.6) can be shown to be equivalent to the following differential formula

(3.12) $$\begin{eqnarray}\unicode[STIX]{x1D662}_{ij}^{dd}=-\left[1+\frac{a_{i}^{2}+a_{j}^{2}}{10}\unicode[STIX]{x1D735}^{2}\right]\boldsymbol{{\mathcal{I}}}\boldsymbol{ : }\overrightarrow{\unicode[STIX]{x1D735}}\unicode[STIX]{x1D64F}_{0}(\boldsymbol{R}_{ij})\overleftarrow{\unicode[STIX]{x1D735}}\boldsymbol{ : }\boldsymbol{{\mathcal{I}}},\end{eqnarray}$$

where the arrow over the operator stands for the direction of its action, e.g. $[\unicode[STIX]{x1D64F}(\boldsymbol{r})\overleftarrow{\unicode[STIX]{x1D735}}]_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FE}}=\unicode[STIX]{x2202}_{\unicode[STIX]{x1D6FE}}T_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}(\boldsymbol{r})$ . This yields for $f^{i}$

(3.13a-c ) $$\begin{eqnarray}f^{1}(R_{ij})=a_{i}^{3}a_{j}^{3}\left(\frac{5}{R_{ij}^{3}}-\frac{8(a_{i}^{2}+a_{j}^{2})}{R_{ij}^{5}}\right),\quad f^{2}(R_{ij})=2a_{i}^{3}a_{j}^{3}\frac{a_{i}^{2}+a_{j}^{2}}{R_{ij}^{5}},\quad f^{0}=-2f^{1}-2f^{2}.\end{eqnarray}$$

The above formulae, for equal-sized particles, have been derived previously by Durlofsky et al. (Reference Durlofsky, Brady and Bossis1987), whereas the equivalent results for different-sized particles can be found in Kim & Karrila (Reference Kim and Karrila1991). However, the case of overlapping particles has not been considered in the literature. Here, the differential formula (3.12) is no longer valid, whilst the integral approach based on (3.6) can still be applied. The algebraic forms of the $f^{i}$ functions for overlapping particles are given in the Appendix. As the smaller sphere gets completely immersed in the larger one, the dipolar component of the inverse friction matrix approaches a constant value:

(3.14) $$\begin{eqnarray}\unicode[STIX]{x1D662}_{ij}^{dd}=\frac{3}{20\unicode[STIX]{x03C0}\unicode[STIX]{x1D702}_{0}a_{{>}}^{3}}\boldsymbol{{\mathcal{I}}},\end{eqnarray}$$

where $a_{{>}}$ is the radius of the larger (external) sphere. Note that the above is equal to the ‘self’ element of the larger particle, as given by (3.7). This is analogous to the earlier results of Wajnryb et al. (Reference Wajnryb, Mizerski, Zuk and Szymczak2013) for $t$ and $r$ components.

Next, we turn to the analysis of the translational-dipolar components. Symmetry considerations lead then to the following form of $\boldsymbol{m}^{td}$ :

(3.15) $$\begin{eqnarray}\displaystyle [\unicode[STIX]{x1D662}_{ij}^{td}]_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FE}}=[\unicode[STIX]{x1D662}_{ji}^{dt}]_{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D6FC}} & = & \displaystyle \frac{3}{20\unicode[STIX]{x03C0}\unicode[STIX]{x1D702}_{0}a_{j}^{2}}\left[h_{0}(R_{ij})\left(\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}\hat{R}_{ij,\unicode[STIX]{x1D6FE}}+\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FE}}\hat{R}_{ij,\unicode[STIX]{x1D6FD}}-\frac{2}{3}\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FE}}\hat{R}_{ij,\unicode[STIX]{x1D6FC}}\right)\right.\nonumber\\ \displaystyle & & \displaystyle \left.+\,h_{1}(R_{ij})\hat{R}_{ij,\unicode[STIX]{x1D6FC}}\left(\hat{R}_{ij,\unicode[STIX]{x1D6FD}}\hat{R}_{ij,\unicode[STIX]{x1D6FE}}-\frac{1}{3}\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FE}}\right)\right].\end{eqnarray}$$

For non-overlapping particles, the functions $h^{i}$ are given by

(3.16a,b ) $$\begin{eqnarray}h_{0}(R_{ij})=\frac{a_{j}^{2}(5a_{i}^{2}+3a_{j}^{2})}{6R_{ij}^{4}},\quad h_{1}(R_{ij})=\frac{5}{6}a_{j}^{2}\left(\frac{3}{R_{ij}^{2}}-\frac{5a_{i}^{2}+3a_{j}^{2}}{R_{ij}^{4}}\right).\end{eqnarray}$$

Finally, the $\boldsymbol{m}^{rd}$ and $\boldsymbol{m}^{dr}$ components are of the form

(3.17) $$\begin{eqnarray}[\unicode[STIX]{x1D662}_{ij}^{rd}]_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FE}}=[\unicode[STIX]{x1D662}_{ji}^{dr}]_{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D6FC}}=\frac{3}{16\unicode[STIX]{x03C0}\unicode[STIX]{x1D702}_{0}a_{j}^{3}}g(R_{ij})\hat{R}_{ij,\unicode[STIX]{x1D6FF}}(\unicode[STIX]{x1D716}_{\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}\hat{R}_{ij,\unicode[STIX]{x1D6FE}}+\unicode[STIX]{x1D716}_{\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FE}}\hat{R}_{ij,\unicode[STIX]{x1D6FD}}),\end{eqnarray}$$

where, for non-overlapping particles, $g(R_{ij})=a_{j}^{3}/R_{ij}^{3}$ . The results for overlapping particles for both $td$ and $rd$ components are given in the Appendix.

3.1 Hydrodynamic matrices of a rigid body

Here we apply the formalism introduced in the preceding section to the description of a macromolecule, immersed in an external linear flow

(3.18) $$\begin{eqnarray}\boldsymbol{v}_{\infty }(\boldsymbol{r})=\unicode[STIX]{x1D646}_{c}\boldsymbol{\cdot }(\boldsymbol{r}-\boldsymbol{R}_{0}),\end{eqnarray}$$

where $\unicode[STIX]{x1D646}_{c}$ is a constant velocity gradient matrix. We will represent the macromolecule as a collection of rigidly connected beads of different radii, $a_{i}$ , which can potentially overlap. When the macromolecule moves as a rigid body, the translational and rotational velocities of individual constituents ( $\boldsymbol{U}_{i}$ , $\unicode[STIX]{x1D734}_{i}$ ) are linked with those of the body as a whole ( $\boldsymbol{U}_{c}$ , $\unicode[STIX]{x1D734}_{c}$ ) in the following way

(3.19) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\boldsymbol{U}_{i}=\boldsymbol{U}_{c}+\unicode[STIX]{x1D734}_{c}\times (\boldsymbol{R}_{i}-\boldsymbol{R}_{0})+\boldsymbol{E}_{c}\boldsymbol{\cdot }(\boldsymbol{R}_{i}-\boldsymbol{R}_{0}),\\ \unicode[STIX]{x1D734}_{i}=\unicode[STIX]{x1D734}_{c},\\ \boldsymbol{ E}_{\infty ,i}=\boldsymbol{E}_{c},\end{array}\right\} & & \displaystyle\end{eqnarray}$$

where $\boldsymbol{E}_{c}$ is the symmetric part of $\unicode[STIX]{x1D646}_{c}$ and $\boldsymbol{R}_{0}$ is the reference point with respect to which the translation and rotation are defined. The relations (3.19) can be written as

(3.20) $$\begin{eqnarray}\left(\begin{array}{@{}c@{}}\tilde{\boldsymbol{U}}\\ \tilde{\unicode[STIX]{x1D734}}\\ \tilde{\boldsymbol{E}}_{\infty }\end{array}\right)=\unicode[STIX]{x1D64F}_{N}(\boldsymbol{R}_{0})\boldsymbol{\cdot }\left(\begin{array}{@{}c@{}}\boldsymbol{U}_{c}\\ \unicode[STIX]{x1D734}_{c}\\ \boldsymbol{E}_{c}\end{array}\right),\end{eqnarray}$$

where $\unicode[STIX]{x1D64F}_{N}$ is $11N\times 11$ matrix. We can also define the rigid-body grand friction matrix $\unicode[STIX]{x1D73B}_{c}$ , by a relation analogous to (2.1) involving $\boldsymbol{U}_{c}$ , $\unicode[STIX]{x1D734}_{c}$ and $\boldsymbol{E}_{c}$ as well as the forces, torques and stresslets with which the macromolecule acts on the fluid. The relation between $\unicode[STIX]{x1D73B}_{c}$ and the $N$ -body grand friction matrix $\unicode[STIX]{x1D73B}_{N}$ is then given by

(3.21) $$\begin{eqnarray}\unicode[STIX]{x1D73B}_{c}=\unicode[STIX]{x1D64F}_{N}(\boldsymbol{R}_{0})^{\text{T}}\boldsymbol{\cdot }\unicode[STIX]{x1D73B}_{N}\boldsymbol{\cdot }\unicode[STIX]{x1D64F}_{N}(\boldsymbol{R}_{0}),\end{eqnarray}$$

with $\unicode[STIX]{x1D64F}_{N}$ defined in (3.20).

The construction of the hydrodynamic matrices in the GRPY approximation now proceeds as follows. First, we construct the $N$ -body inverse friction matrix $\unicode[STIX]{x1D662}_{N}$ , as outlined in § 3. Next, it is inverted to obtain the $N$ -body friction matrix $\unicode[STIX]{x1D73B}_{N}$ . Subsequently, the friction matrix of the macromolecule is obtained by projection (3.21). Finally, $\unicode[STIX]{x1D73B}_{c}$ is partially inverted to yield the rigid-body grand mobility matrix $\unicode[STIX]{x1D741}_{c}$ , which can then be used to calculate the intrinsic viscosity of the system:

(3.22) $$\begin{eqnarray}[\unicode[STIX]{x1D702}]_{\infty }=\mathop{\sum }_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}\frac{[\unicode[STIX]{x1D741}_{c}^{dd}]_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FC}}}{10\unicode[STIX]{x1D702}_{0}v}.\end{eqnarray}$$

In the above, intrinsic viscosity is normalized by the volume of the macromolecule, $v$ . To be precise, the above corresponds to the high-frequency limit of intrinsic viscosity, since we neglected the contribution from the Brownian relaxation of the particle orientation (Rallison Reference Rallison1978; Cichocki, Ekiel-Jezewska & Wajnryb Reference Cichocki, Ekiel-Jezewska and Wajnryb2012).

The projection operator $\unicode[STIX]{x1D64F}_{N}(\boldsymbol{R}_{0})$ defined in (3.19)–(3.20) involves an arbitrary reference point $\boldsymbol{R}_{0}$ . A natural question which can be raised in this context is whether the choice of $\boldsymbol{R}_{0}$ has an effect on the elements of the rigid-body matrices $\unicode[STIX]{x1D73B}_{c}$ and $\unicode[STIX]{x1D741}_{c}$ . As it turns out, the answer is different for different components of these matrices. Translational and rotational components of the grand mobility matrix have been analysed in Kim & Karrila (Reference Kim and Karrila1991) and it has been shown that $\unicode[STIX]{x1D741}_{c}^{tt}$ , $\unicode[STIX]{x1D741}_{c}^{tr}$ and $\unicode[STIX]{x1D741}_{c}^{rt}$ depend on the choice of the origin, whereas $\unicode[STIX]{x1D741}_{c}^{rr}$ does not. Generalizing this reasoning to dipolar components it can be shown that $\unicode[STIX]{x1D741}_{c}^{dd}$ , $\unicode[STIX]{x1D741}_{c}^{dr}$ and $\unicode[STIX]{x1D741}_{c}^{rd}$ are not origin sensitive. In some of the alternative approaches to the calculation of intrinsic viscosity one finds that the result is origin sensitive (de la Torre & Bloomfield Reference de la Torre and Bloomfield1978; Harding Reference Harding1997) and it is proposed $\boldsymbol{R}_{0}$ needs to be chosen in such a way as to guarantee the minimum energy dissipation (so called ‘viscosity centre’). A possible reason for this discrepancy is that the derivation by de la Torre & Bloomfield (Reference de la Torre and Bloomfield1978) involves approximations on the level of the friction matrix, the elements of which do depend on the choice of origin.

4 Examples

We demonstrate the performance of the GRPY approximation by calculating the intrinsic viscosity for two simple model systems. In the first example, we represent a spherical particle of radius $a_{c}$ as a collection of beads, either positioned on the surface of the sphere only or positioned throughout its entire volume. This corresponds to the shell model and filling model in the nomenclature of Carrasco & de la Torre (Reference Carrasco and de la Torre1999a ). In both cases, the starting point is a three-dimensional cubic lattice with lattice constant $\unicode[STIX]{x1D706}$ . For the filling model, we place the spheres in all the nodes that are at a distance $a_{c}$ or smaller from the centre. For shell model, we choose the outer layer of the respective volume filling (figure 1). Furthermore, to calculate the intrinsic viscosity we need to assign an effective radius to such a conglomerate in order to estimate its volume as needed in (3.22). As shown by Cichocki, Ekiel-Jezewska & Wajnryb (Reference Cichocki, Ekiel-Jezewska and Wajnryb2014), the correct procedure in such a case is to estimate $v$ based on the hydrodynamic (Stokes) radius, defined as $a_{hyd}=(2\unicode[STIX]{x03C0}\unicode[STIX]{x1D702}_{0}\text{Tr}\unicode[STIX]{x1D741}^{tt})^{-1}$ . For an almost spherical shape, this approach leads to the correct value of the intrinsic viscosity, up to quadratic terms in surface roughness. Note that $\unicode[STIX]{x1D741}^{tt}$ in the formula for the Stokes radius needs to be calculated using the same approach as the one used to compute $\unicode[STIX]{x1D741}^{dd}$ (i.e. bead model combined with the GRPY approximation).

Figure 1.

Sphere modelled as a conglomerate of beads with different radii: $a=0.5\unicode[STIX]{x1D706}$ (a) and $a=0.2\unicode[STIX]{x1D706}$ (b). The spheres are positioned on a cubic grid with a lattice constant $\unicode[STIX]{x1D706}$ (c).

Next, we analyse the accuracy of the GRPY approximation for different radii ( $a$ ) of the constituent spheres by comparing their prediction with the exact result for a single sphere ( $[\unicode[STIX]{x1D702}]_{\infty }=5/2$ ). As can be inferred from the data in table 1, both models perform well for a wide range of bead radii, with an accuracy not worse than approximately 4 %, even at $a/\unicode[STIX]{x1D706}=0.05$ , which corresponds to a very sparse filling. At the other extreme, for the shell model with touching spheres ( $a/\unicode[STIX]{x1D706}=0.5$ ) we get an accuracy of approximately 0.5 %. The accuracy of the results increases significantly with the decrease of the lattice constant $\unicode[STIX]{x1D706}$ , which results in an increase of the total number of spheres ( $N$ ) used to represent our molecule. As shown in table 1, an approximately twofold increase in $N$ results in an approximately twofold reduction in the deviation from the exact result. Importantly, the model continues to perform well for overlapping beads ( $a/\unicode[STIX]{x1D706}>0.5$ ), which demonstrates its applicability to macromolecular models, where overlapping beads are often used to represent complex biological structures. The shell model performs better than the filling model, since – for a given number of beads – it represents the surface of the sphere in a more precise way. The omission of the internal beads does not affect the accuracy while allowing for a reduction in the computational complexity of the model.

Figure 2.

High-frequency component of the intrinsic viscosity $[\unicode[STIX]{x1D702}]_{\infty }$ for the rigid dumbbell as a function of the beads separation length $l$ measured in unit lengths of the bead radius, $a$ . Three different methods are used to obtain the results: the GRPY approximation with two beads (one per each sphere in the dumbbell), the multipole method and the GRPY approximation with 895 beads filling each sphere of the dumbbell. Red triangles mark the values of intrinsic viscosity obtained using the GRPY approximation for free beads, i.e. without imposing the constraint (3.20). The pictograms below the horizontal axis show schematically the geometry of the dumbbell for the respective $l/a$ values.

Table 1.

Intrinsic viscosity of a spherical particle represented by the shell model (with the beads covering the surface of the sphere) and the filling model (with beads filling the entire volume of the sphere) for different number of beads ( $N$ ).

In the second test of the GRPY approximation we consider a dumbbell with two equal-sized spheres of radii $a_{0}$ rigidly connected to each other. Figure 2 shows the intrinsic viscosity of such a system as a function of the distance between the centres of the spheres $(l)$ . Here, we compare the results of the simple GRPY approximation with two beads (one per each sphere in the dumbbell) against two more precise methods: the multipole series expansion method (Cichocki et al. Reference Cichocki, Felderhof, Hinsen, Wajnryb and Bławzdziewicz1994) as well as a representation of the dumbbell by a filling model with the lattice constant $\unicode[STIX]{x1D706}=a_{0}/6$ and $a=\unicode[STIX]{x1D706}/2$ (this corresponds to the representation of each sphere of the dumbbell by $895$ smaller, touching beads).

As expected, when the particles are widely separated, the results of the GRPY approximation coincide with those of the multipole method, since GRPY becomes exact at the far-field limit. However, the two-sphere GRPY model also works very well when the particles are overlapping, with a deviation from the result of the multi-bead model of not more than 5 %. Based on the data in table 1 we expect the multi-bead model constructed with this many beads to be within 0.5 % of the exact result. One of the reasons for such good performance of the GRPY approximation is its correct limiting behaviour at a complete overlap of the spheres (cf. (3.14)). The non-monotonic dependence of $[\unicode[STIX]{x1D702}]_{\infty }$ on $l/a$ is the result of normalization of the intrinsic viscosity by the volume of the dumbbell (3.22). Both the total dissipation and the volume increase with $l/a$ , but at a different rate, which leads to the overall non-monotonic behaviour of $[\unicode[STIX]{x1D702}]_{\infty }(l/a)$ .

An additional advantage of the present scheme is that it allows for the separation of two different contributions to intrinsic viscosity: that coming from the presence of internal constraints between the beads and that coming from the finite size of the beads themselves (Rallison Reference Rallison1979). Red triangles in Figure 2 mark the values of $\unicode[STIX]{x1D702}$ calculated without imposing the constraint (3.20). As observed, if the beads are close to each other or overlap, the viscosity of free beads is close to that of the constrained system, i.e. the contribution of constraints is relatively small. Then, as the distance between the beads is increased, constraints play an increasingly important role (for $l/a=4$ their contribution equals the one connected with the finite size of the beads). Importantly, the GRPY approach incorporates both contributions to the intrinsic viscosity, since it takes explicit account of the volume of the beads, in contrast to the Oseen-level approaches. The data in Figure 2 indicates that the GPRY approach will be particularly effective for modelling compact systems, such as proteins and protein complexes, whereas for long linear polymers the Oseen-based approaches might prove sufficient.

5 Summary

Using a reformulation of the RPY approximation based on the integral formalism (equations (3.6)) of Wajnryb et al. (Reference Wajnryb, Mizerski, Zuk and Szymczak2013) we have provided explicit formulae for the components of the inverse friction matrix, $\unicode[STIX]{x1D662}$ , both for non-overlapping and overlapping particles. These have then been used for the calculation of the intrinsic viscosity of complex-shaped macromolecules, using their representation as a collection of beads. By inverting the $N$ -particle $\unicode[STIX]{x1D662}$ matrix and then projecting the resulting friction matrix on the rigid-body motion of the molecule, we recover the friction and mobility matrices of the macromolecule, $\unicode[STIX]{x1D73B}_{c}$ and $\unicode[STIX]{x1D741}_{c}$ . The dipolar element of the latter is then used to calculate the high-frequency intrinsic viscosity, $[\unicode[STIX]{x1D702}]_{\infty }$ .

We have tested the RPY approximation for $[\unicode[STIX]{x1D702}]_{\infty }$ by comparing its predictions either to exact results (for a sphere) or to the predictions of the virtually exact multipole series expansion method (for a dumbbell). GRPY has been shown to perform well not only for separated beads but also for overlapping particles. Importantly, the GRPY tensors were shown to have the correct limiting behaviour at a complete overlap, with one sphere fully immersed in another.