2.1 Governing equations and boundary conditions

The flow field around the moving porous particles depends on their permeability. When the two particles are impermeable (zero permeability), the fluid has to travel around the particles. If the permeability is non-zero, the fluid can pass through the two particles (region marked in figure 1 with a dashed rectangle). In such a case, at low Reynolds number, the governing equations for the model are the Stokes equation for the external steady Newtonian fluid flow outside the particles ( $r_{a}>a$ , $r_{b}>b$ ) and the Brinkman equation for the internal porous region ( $0<r_{a}<a$ , $0<r_{b}<b$ ):

(2.1) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FB}^{2}\boldsymbol{v}^{e}=\frac{1}{\unicode[STIX]{x1D707}}\unicode[STIX]{x1D735}p^{e},\quad \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{v}^{e}=0, & \displaystyle\end{eqnarray}$$

(2.2) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FB}^{2}\boldsymbol{v}^{i}=\frac{1}{\unicode[STIX]{x1D707}}\unicode[STIX]{x1D735}p^{i}+\frac{1}{k}\boldsymbol{v}^{i},\quad \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{v}^{i}=0, & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D707}$ is the viscosity of the fluid and $k$ is the permeability of the porous material. The superscripts ‘ $e$ ’ and ‘ $i$ ’ represent the external and internal flows of the porous particles, respectively. The flow field $\boldsymbol{v}$ and induced stress tensors components $\unicode[STIX]{x1D635}_{rr}$ and $\unicode[STIX]{x1D635}_{r\unicode[STIX]{x1D703}}$ are continuous at the fluid-particle interface. The boundary conditions are

(2.3) $$\begin{eqnarray}\displaystyle & (\boldsymbol{v}^{e}-\boldsymbol{U})\boldsymbol{\cdot }\boldsymbol{n}=\boldsymbol{v}^{i}\boldsymbol{\cdot }\boldsymbol{n}, & \displaystyle\end{eqnarray}$$

(2.4) $$\begin{eqnarray}\displaystyle & (\boldsymbol{v}^{e}-\boldsymbol{U})\boldsymbol{\cdot }\boldsymbol{s}=\boldsymbol{v}^{i}\boldsymbol{\cdot }\boldsymbol{s}, & \displaystyle\end{eqnarray}$$

(2.5) $$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D635}_{rr}^{e}=\unicode[STIX]{x1D635}_{rr}^{i}, & \displaystyle\end{eqnarray}$$

(2.6) $$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D635}_{r\unicode[STIX]{x1D703}}^{e}=\unicode[STIX]{x1D635}_{r\unicode[STIX]{x1D703}}^{i}, & \displaystyle\end{eqnarray}$$

where $\boldsymbol{n}$ is the unit vector normal to the porous-fluid interface and pointing to the external fluid and $\boldsymbol{s}$ is the tangential vector at the interface. It should be noted that $\boldsymbol{v}^{i}$ represents the relative velocity of the fluid inside the porous particle with respect to the origin of the particle and $\boldsymbol{v}^{e}$ refers to the velocity in an inertial reference frame. Therefore, the particle velocity must be subtracted from $\boldsymbol{v}^{e}$ to construct the boundary conditions (2.3)–(2.4). On the wall of the cylindrical tube, the no-slip boundary condition applies:

(2.7) $$\begin{eqnarray}\boldsymbol{v}^{e}=0.\end{eqnarray}$$

2.2 Semi-analytical implementation

2.2.1 Calculation procedure based on the method of reflections

Since the governing equations and boundary conditions are linear, the velocity and pressure fields can be expressed as a summation of fields:

(2.8) $$\begin{eqnarray}\displaystyle & \boldsymbol{v}=\boldsymbol{v}_{c}^{(0)}+\boldsymbol{v}_{a}^{(1)}+\boldsymbol{v}_{b}^{(2)}+\boldsymbol{v}_{c}^{(3)}+\boldsymbol{v}_{a}^{(4)}+\boldsymbol{v}_{b}^{(5)}+\cdots , & \displaystyle\end{eqnarray}$$

(2.9) $$\begin{eqnarray}\displaystyle & p=p_{c}^{(0)}+p_{a}^{(1)}+p_{b}^{(2)}+p_{c}^{(3)}+p_{a}^{(4)}+p_{b}^{(5)}+\cdots , & \displaystyle\end{eqnarray}$$

where the superscript ( $j$ ) refers to the $j$ th reflection ( $j=1,2,3\ldots$ ) and subscripts $a$ , $b$ and $c$ refer to the boundary conditions satisfying particle  $A$ , particle $B$ and the cylinder, respectively. Each term of $\boldsymbol{v}$ and  $p$ separately satisfies the governing equations (2.1)–(2.2) and the boundary conditions resulting from the previous terms. For example, if the subscript of the last term is  $b$ , the summation of the terms satisfies the boundary condition of particle  $B$ . It should be noted that the reflections representing the flow fields due to the porous particles (reflections with subscripts ‘ $a$ ’ and ‘ $b$ ’) consist of two components: the flows inside and outside the porous particles with superscripts ‘ $i$ ’ and ‘ $e$ ’, respectively. Therefore, equation (2.8) becomes

(2.10) $$\begin{eqnarray}\displaystyle & \boldsymbol{v}^{e}=\boldsymbol{v}_{c}^{(0)}+\boldsymbol{v}_{a}^{(1,e)}+\boldsymbol{v}_{b}^{(2,e)}+\boldsymbol{v}_{c}^{(3)}+\boldsymbol{v}_{a}^{(4,e)}+\boldsymbol{v}_{b}^{(5,e)}+\cdots , & \displaystyle\end{eqnarray}$$

(2.11) $$\begin{eqnarray}\displaystyle & \boldsymbol{v}_{a}^{i}=\boldsymbol{v}_{a}^{(1,i)}+\boldsymbol{v}_{a}^{(4,i)}+\boldsymbol{v}_{a}^{(7,i)}+\cdots , & \displaystyle\end{eqnarray}$$

(2.12) $$\begin{eqnarray}\displaystyle & \boldsymbol{v}_{b}^{i}=\boldsymbol{v}_{b}^{(2,i)}+\boldsymbol{v}_{b}^{(5,i)}+\boldsymbol{v}_{b}^{(8,i)}+\cdots . & \displaystyle\end{eqnarray}$$

Based on the boundary conditions (2.3)–(2.6), the following relations can be obtained:

(2.13) $$\begin{eqnarray}\displaystyle & \boldsymbol{v}_{c}^{(0)}=[U_{0}(1-R^{2}/R_{0}^{2})]\boldsymbol{Z}, & \displaystyle\end{eqnarray}$$

(2.14) $$\begin{eqnarray}\displaystyle & \boldsymbol{v}_{a}^{(1)}=\left\{\begin{array}{@{}ll@{}}(\boldsymbol{v}_{a}^{(1,e)}+\boldsymbol{v}_{c}^{(0)}-\boldsymbol{U}_{a})\boldsymbol{\cdot }\boldsymbol{n}=\boldsymbol{v}_{a}^{(1,i)}\boldsymbol{\cdot }\boldsymbol{n}, & \text{at}~r_{a}=a,\\ (\boldsymbol{v}_{a}^{(1,e)}+\boldsymbol{v}_{c}^{(0)}-\boldsymbol{U}_{a})\boldsymbol{\cdot }\boldsymbol{s}=\boldsymbol{v}_{a}^{(1,i)}\boldsymbol{\cdot }\boldsymbol{s}, & \text{at}~r_{a}=a,\\ \unicode[STIX]{x1D635}_{a,rr}^{(1,e)}=\unicode[STIX]{x1D635}_{a,rr}^{(1,i)}, & \text{at}~r_{a}=a,\\ \unicode[STIX]{x1D635}_{a,r\unicode[STIX]{x1D703}}^{(1,e)}=\unicode[STIX]{x1D635}_{a,r\unicode[STIX]{x1D703}}^{(1,i)}, & \text{at}~r_{a}=a,\\ 0, & \text{at}~Z=\pm \infty ,\end{array}\right. & \displaystyle\end{eqnarray}$$

(2.15) $$\begin{eqnarray}\displaystyle & \boldsymbol{v}_{b}^{(2)}=\left\{\begin{array}{@{}ll@{}}(\boldsymbol{v}_{b}^{(2,e)}+\boldsymbol{v}_{a}^{(1,e)}+\boldsymbol{v}_{c}^{(0)}-\boldsymbol{U}_{b})\boldsymbol{\cdot }\boldsymbol{n}=\boldsymbol{v}_{b}^{(2,i)}\boldsymbol{\cdot }\boldsymbol{n}, & \text{at}~r_{b}=b,\\ (\boldsymbol{v}_{b}^{(2,e)}+\boldsymbol{v}_{a}^{(1,e)}+\boldsymbol{v}_{c}^{(0)}-\boldsymbol{U}_{b})\boldsymbol{\cdot }\boldsymbol{s}=\boldsymbol{v}_{b}^{(2,i)}\boldsymbol{\cdot }\boldsymbol{s}, & \text{at}~r_{b}=b,\\ \unicode[STIX]{x1D635}_{b,rr}^{(2,e)}=\unicode[STIX]{x1D635}_{b,rr}^{(2,i)}, & \text{at}~r_{b}=b,\\ \unicode[STIX]{x1D635}_{b,r\unicode[STIX]{x1D703}}^{(2,e)}=\unicode[STIX]{x1D635}_{b,r\unicode[STIX]{x1D703}}^{(2,i)}, & \text{at}~r_{b}=b,\\ 0, & \text{at}~Z=\pm \infty ,\end{array}\right. & \displaystyle\end{eqnarray}$$

(2.16) $$\begin{eqnarray}\displaystyle & \boldsymbol{v}_{c}^{(3)}=\left\{\begin{array}{@{}ll@{}}-\boldsymbol{v}_{a}^{(1,e)}-\boldsymbol{v}_{b}^{(2,e)}, & \text{at}~R=R_{0},\\ 0, & \text{at}~Z=\pm \infty ,\end{array}\right. & \displaystyle\end{eqnarray}$$

(2.17) $$\begin{eqnarray}\displaystyle & \boldsymbol{v}_{a}^{(4)}=\left\{\begin{array}{@{}ll@{}}(\boldsymbol{v}_{a}^{(4,e)}+\boldsymbol{v}_{b}^{(2,e)}+\boldsymbol{v}_{c}^{(3)})\boldsymbol{\cdot }\boldsymbol{n}=\boldsymbol{v}_{a}^{(4,i)}\boldsymbol{\cdot }\boldsymbol{n}, & \text{at}~r_{a}=a,\\ (\boldsymbol{v}_{a}^{(4,e)}+\boldsymbol{v}_{b}^{(2,e)}+\boldsymbol{v}_{c}^{(3)})\boldsymbol{\cdot }\boldsymbol{s}=\boldsymbol{v}_{a}^{(4,i)}\boldsymbol{\cdot }\boldsymbol{s}, & \text{at}~r_{a}=a,\\ \unicode[STIX]{x1D635}_{a,rr}^{(4,e)}=\unicode[STIX]{x1D635}_{a,rr}^{(4,i)}, & \text{at}~r_{a}=a,\\ \unicode[STIX]{x1D635}_{a,r\unicode[STIX]{x1D703}}^{(4,e)}=\unicode[STIX]{x1D635}_{a,r\unicode[STIX]{x1D703}}^{(4,i)}, & \text{at}~r_{a}=a,\\ 0, & \text{at}~Z=\pm \infty ,\end{array}\right. & \displaystyle\end{eqnarray}$$

(2.18) $$\begin{eqnarray}\displaystyle & \boldsymbol{v}_{b}^{(5)}=\left\{\begin{array}{@{}ll@{}}(\boldsymbol{v}_{b}^{(5,e)}+\boldsymbol{v}_{a}^{(4,e)}+\boldsymbol{v}_{c}^{(3)})\boldsymbol{\cdot }\boldsymbol{n}=\boldsymbol{v}_{a}^{(5,i)}\boldsymbol{\cdot }\boldsymbol{n}, & \text{at}~r_{b}=b,\\ (\boldsymbol{v}_{b}^{(5,e)}+\boldsymbol{v}_{a}^{(4,e)}+\boldsymbol{v}_{c}^{(3)})\boldsymbol{\cdot }\boldsymbol{s}=\boldsymbol{v}_{a}^{(5,i)}\boldsymbol{\cdot }\boldsymbol{s}, & \text{at}~r_{b}=b,\\ \unicode[STIX]{x1D635}_{b,rr}^{(5,e)}=\unicode[STIX]{x1D635}_{b,rr}^{(5,i)}, & \text{at}~r_{b}=b,\\ \unicode[STIX]{x1D635}_{b,r\unicode[STIX]{x1D703}}^{(5,e)}=\unicode[STIX]{x1D635}_{b,r\unicode[STIX]{x1D703}}^{(5,i)}, & \text{at}~r_{b}=b,\\ 0, & \text{at}~Z=\pm \infty ,\end{array}\right. & \displaystyle\end{eqnarray}$$

$$\begin{eqnarray}\vdots\end{eqnarray}$$

For illustration, the calculation is conducted from the background flow having a parabolic profile, $\boldsymbol{v}_{c}^{(0)}$ . The following reflection, $\boldsymbol{v}_{a}^{(1)}$ , due to the presence of particle  $A$ , is required to make the summation $\boldsymbol{v}_{c}^{(0)}+\boldsymbol{v}_{a}^{(1)}$ satisfy the boundary condition of particle  $A$ (as shown in (2.3)–(2.6)). To acquire the boundary condition of $\boldsymbol{v}_{a}^{(1)}$ , we transform the expression of $\boldsymbol{v}_{c}^{(0)}$ in the cylindrical coordinates into the spherical coordinates (details of the coordinate transformation are provided in appendix A). The expression of $\boldsymbol{v}_{a}^{(1)}$ can be calculated by means of a stream function, which will be described in § 2.2.2. Similarly, the following reflection $\boldsymbol{v}_{b}^{(2)}$ is included to make the summation $\boldsymbol{v}_{c}^{(0)}+\boldsymbol{v}_{a}^{(1)}+\boldsymbol{v}_{b}^{(2)}$ satisfy the boundary condition of particle  $B$ . The expression can be solved using a similar procedure in the spherical coordinates of particle  $B$ . The added reflection $\boldsymbol{v}_{c}^{(3)}$ is needed to make the summation satisfy the boundary condition of the tube. The expressions derived in spherical coordinates of particle  $A$ , $\boldsymbol{v}_{a}^{(1)}$ , and particle  $B$ , $\boldsymbol{v}_{b}^{(2)}$ , should be transformed into cylindrical coordinates, such that $\boldsymbol{v}_{c}^{(3)}$ can be solved by utilizing cylindrical harmonics shown in § 2.2.3.

In general, the reflection $\boldsymbol{v}_{a}^{(1)}$ arising from $\boldsymbol{v}_{c}^{(0)}$ brings in the representative geometry dimensions $a$ (the radius of particle  $A$ ) and $R_{0}$ (the radius of the tube) included in the form $a/R_{0}$ by using a stream function. This indicates that the error at least scales as $a/R_{0}$ (Kim & Karrila Reference Kim and Karrila2005). The next reflection $\boldsymbol{v}_{b}^{(2)}$ is calculated from $\boldsymbol{v}_{a}^{(1)}$ and $\boldsymbol{v}_{c}^{(0)}$ . The velocity field introduces the term $b/R_{0}$ and $b/L$ , where $L$ is the separation distance between the two particles. This indicates that the error at least scales as the maximum value of $b/R_{0}$ and $b/L$ . The following reflection $\boldsymbol{v}_{c}^{(3)}$ is required to satisfy the boundary condition: $\boldsymbol{v}_{c}^{(3)}=-\boldsymbol{v}_{a}^{(1)}-\boldsymbol{v}_{b}^{(2)}$ at the cylinder wall. The magnitude of $\boldsymbol{v}_{c}^{(3)}$ scales as the same order of $\boldsymbol{v}_{a}^{(1)}$ and $\boldsymbol{v}_{b}^{(2)}$ at the wall, which is smaller than that near $\boldsymbol{v}_{a}^{(1)}$ and $\boldsymbol{v}_{b}^{(2)}$ . This indicates the error of $\boldsymbol{v}_{c}^{(3)}$ is smaller than the error of $\boldsymbol{v}_{a}^{(1)}$ and $\boldsymbol{v}_{b}^{(2)}$ . Therefore, the subsequent reflection can correct the error of the flow field to the higher order. By repeating the procedures and adding more reflections, the results are expected to converge. This convergence is verified in all our calculations (see appendix B).

2.2.2 Solution based on stream function

The solution of the external and internal flows of the porous particle in spherical coordinates can be expressed using the stream functions (Saad & Faltas Reference Saad and Faltas2014)

(2.19) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6F7}^{e}=\mathop{\sum }_{n=2}^{\infty }(A_{n}r^{-n+1}+B_{n}r^{-n+3})C_{n}(\cos \unicode[STIX]{x1D703}), & \displaystyle\end{eqnarray}$$

(2.20) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6F7}^{i}=\mathop{\sum }_{n=2}^{\infty }(\tilde{A}_{n}r^{n}+\sqrt{r}\tilde{B}_{n}I_{n-1/2}(\unicode[STIX]{x1D6FC}r))C_{n}(\cos \unicode[STIX]{x1D703}), & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D6FC}=\sqrt{1/k}$ ; $I_{n-1/2}$ is the modified Bessel function of the first kind; $C_{n}$ is the Gegenbauer polynomial of the first kind of order $n$ and degree $-1/2$ ; $A_{n}$ , $B_{n}$ , $\tilde{A}_{n}$ and $\tilde{B}_{n}$ are unknown coefficients which can be determined from the boundary conditions as shown in (2.3)–(2.6).

The velocity components in spherical coordinates can be derived from the stream functions

(2.21) $$\begin{eqnarray}\displaystyle & \displaystyle v_{r}^{e,i}=-\frac{1}{r^{2}\sin \unicode[STIX]{x1D703}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6F7}^{e,i}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}, & \displaystyle\end{eqnarray}$$

(2.22) $$\begin{eqnarray}\displaystyle & \displaystyle v_{\unicode[STIX]{x1D703}}^{e,i}=\frac{1}{r\sin \unicode[STIX]{x1D703}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6F7}^{e,i}}{\unicode[STIX]{x2202}r}, & \displaystyle\end{eqnarray}$$

By substituting (2.19)–(2.20) into (2.21)–(2.22), the velocity components can be expressed as

(2.23) $$\begin{eqnarray}\displaystyle & \displaystyle v_{r}^{e}=\mathop{\sum }_{n=2}^{\infty }-r^{-n-1}(A_{n}+B_{n}r^{2})P_{n-1}(\cos \unicode[STIX]{x1D703}), & \displaystyle\end{eqnarray}$$

(2.24) $$\begin{eqnarray}\displaystyle & \displaystyle v_{r}^{i}=\mathop{\sum }_{n=2}^{\infty }-\frac{1}{r^{2}}[r^{n}\tilde{A}_{n}+\sqrt{r}\tilde{B}_{n}I_{n-1/2}(\unicode[STIX]{x1D6FC}r)]P_{n-1}(\cos \unicode[STIX]{x1D703}), & \displaystyle\end{eqnarray}$$

(2.25) $$\begin{eqnarray}\displaystyle & \displaystyle v_{\unicode[STIX]{x1D703}}^{e}=\mathop{\sum }_{n=2}^{\infty }-r^{-n-1}[(n-1)A_{n}+(n-3)B_{n}r^{2}]\csc \unicode[STIX]{x1D703}C_{n}(\cos \unicode[STIX]{x1D703}), & \displaystyle\end{eqnarray}$$

(2.26) $$\begin{eqnarray}\displaystyle & \displaystyle v_{\unicode[STIX]{x1D703}}^{i}=\mathop{\sum }_{n=2}^{\infty }-\frac{1}{r^{2}}\{nr^{n}\tilde{A}_{n}+\sqrt{r}[nI_{n-1/2}(\unicode[STIX]{x1D6FC}r)+\unicode[STIX]{x1D6FC}rI_{n+1/2}(\unicode[STIX]{x1D6FC}r)]\tilde{B}_{n}\}\csc \unicode[STIX]{x1D703}C_{n}(\cos \unicode[STIX]{x1D703}), & \displaystyle\end{eqnarray}$$

where $P_{n-1}$ is the Legendre polynomial of order $n-1$ . The components of the stress tensor can be derived from the velocity components:

(2.27) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D635}_{rr}^{e,i}=-p^{e,i}+2\unicode[STIX]{x1D707}\frac{\unicode[STIX]{x2202}v_{r}^{e,i}}{\unicode[STIX]{x2202}r}, & \displaystyle\end{eqnarray}$$

(2.28) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D635}_{r\unicode[STIX]{x1D703}}^{e,i}=\unicode[STIX]{x1D707}\left[r\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}\left(\frac{v_{\unicode[STIX]{x1D703}}^{e,i}}{r}\right)+\frac{1}{r}\frac{\unicode[STIX]{x2202}v_{r}^{e,i}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}\right]. & \displaystyle\end{eqnarray}$$

The pressure can be obtained from the governing equations:

(2.29) $$\begin{eqnarray}\displaystyle & \displaystyle p^{e}=-2\unicode[STIX]{x1D707}\mathop{\sum }_{n=2}^{\infty }\frac{2n-3}{n}B_{n}r^{-n}P_{n-1}(\cos \unicode[STIX]{x1D703}), & \displaystyle\end{eqnarray}$$

(2.30) $$\begin{eqnarray}\displaystyle & \displaystyle p^{i}=\unicode[STIX]{x1D707}\mathop{\sum }_{n=2}^{\infty }\frac{\unicode[STIX]{x1D6FC}^{2}}{n-1}\tilde{A}_{n}r^{n-1}P_{n-1}(\cos \unicode[STIX]{x1D703}). & \displaystyle\end{eqnarray}$$

By substituting these expansions into the boundary conditions (2.3)–(2.6), the coefficients in the stream functions can be obtained. It should be remembered that our implementation in this work is based on the method of reflections. The stream function solutions developed here are applicable for the reflections involving particles  $A$ and  $B$ . For example, equation (2.14) is used to satisfy the boundary condition of particle  $A$ . To construct the algebraic equations for solving the unknown coefficients $A_{n}$ , $B_{n}$ , $\tilde{A}_{n}$ and $\tilde{B}_{n}$ for this reflection $\boldsymbol{v}_{a}^{(1)}$ , we transform the expression of $\boldsymbol{v}_{c}^{(0)}$ from the cylindrical coordinates into the spherical coordinates (details of the coordinate transformation are provided in appendix A). It should be noted that this procedure is not the re-expansion of the solution. It is a direct transformation of the unit vector and variables from cylindrical coordinates to spherical coordinates. Since (2.14) is a boundary condition at the surface of particle $A$ ( $r_{a}=a$ ), the transformed $\boldsymbol{v}_{c}^{(0)}$ is only a function of  $\unicode[STIX]{x1D703}_{a}$ . An arbitrary function of $\unicode[STIX]{x1D703}$ can be expanded into a series using either the Gegenbauer polynomial or Legendre polynomial (Happel & Brenner Reference Happel and Brenner1983):

(2.31) $$\begin{eqnarray}f(\cos \unicode[STIX]{x1D703})=\mathop{\sum }_{n=2}^{\infty }\tilde{a}_{n}C_{n}(\cos \unicode[STIX]{x1D703}),\end{eqnarray}$$

where

(2.32) $$\begin{eqnarray}\displaystyle & \displaystyle \tilde{a}_{n}=\frac{1}{2}n(n-1)(2n-1)\int _{-1}^{1}\frac{f(\cos \unicode[STIX]{x1D703})C_{n}(\cos \unicode[STIX]{x1D703})}{\sin ^{2}\unicode[STIX]{x1D703}}\,\text{d}(\cos \unicode[STIX]{x1D703}), & \displaystyle\end{eqnarray}$$

(2.33) $$\begin{eqnarray}\displaystyle & \displaystyle g(\cos \unicode[STIX]{x1D703})=\mathop{\sum }_{n=2}^{\infty }\tilde{b}_{n-1}P_{n-1}(\cos \unicode[STIX]{x1D703}), & \displaystyle\end{eqnarray}$$

where

(2.34) $$\begin{eqnarray}\displaystyle \tilde{b}_{n-1}=\frac{1}{2}(2n-1)\int _{-1}^{1}g(\cos \unicode[STIX]{x1D703})P_{n-1}(\cos \unicode[STIX]{x1D703})\,\text{d}(\cos \unicode[STIX]{x1D703}).\end{eqnarray}$$

The $\boldsymbol{v}_{c}^{(0)}$ , which is now in spherical coordinates, can be expressed into either the Gegenbauer or Legendre polynomial, and substituted into (2.14). The first boundary condition of (2.14) can then be transformed into

(2.35) $$\begin{eqnarray}\displaystyle & & \displaystyle \mathop{\sum }_{n=2}^{\infty }-a^{-n-1}(A_{n}+B_{n}a^{2})P_{n-1}(\cos \unicode[STIX]{x1D703}_{a})+\mathop{\sum }_{n=2}^{\infty }\tilde{b}_{n-1}P_{n-1}(\cos \unicode[STIX]{x1D703}_{a})\nonumber\\ \displaystyle & & \displaystyle \quad =\displaystyle \mathop{\sum }_{n=2}^{\infty }-\frac{1}{a^{2}}[a^{n}\tilde{A}_{n}+\sqrt{a}\tilde{B}_{n}I_{n-1/2}(\unicode[STIX]{x1D6FC}a)]P_{n-1}(\cos \unicode[STIX]{x1D703}_{a}).\end{eqnarray}$$

The first summation term on the left-hand side of (2.35) is the general solution of the external flow of particle  $A$ , $\boldsymbol{v}_{a}^{(1,e)}$ , given by (2.23). The second summation term on the left-hand side is $\boldsymbol{v}_{c}^{(0)}$ expressed in a Legendre polynomial. The right-hand side represents the general solution of the internal flow of particle  $A$ , $\boldsymbol{v}_{a}^{(1,i)}$ , given by (2.24). Equating term by term under the summation sign, an algebraic equation can be constructed. Similarly, the other three algebraic equations can be constructed from (2.14). The coefficients  $A_{n}$ , $B_{n}$ , $\tilde{A}_{n}$ and $\tilde{B}_{n}$ in the stream functions for this reflection $\boldsymbol{v}_{a}^{(1)}$ can then be obtained. Similarly, other reflections with subscripts ‘ $a$ ’ and ‘ $b$ ’ can be solved by using the corresponding boundary conditions satisfying the surface of the two particles.

2.2.3 Solution based on cylindrical harmonics

In cylindrical coordinates, the velocity components are shown in terms of $(R,\unicode[STIX]{x1D6F7},Z)$ . Without loss of generality, the general solution to solve the flow field inside the cylindrical tube is given by Happel & Brenner (Reference Happel and Brenner1983):

(2.36) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{v}=\mathop{\sum }_{k=-\infty }^{\infty }\left[\unicode[STIX]{x1D735}\times (\unicode[STIX]{x1D6FA}_{k}\boldsymbol{Z})+\unicode[STIX]{x1D735}\unicode[STIX]{x1D6F9}_{k}+R\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}R}\unicode[STIX]{x1D735}\unicode[STIX]{x1D6F1}_{k}+\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6F1}_{k}}{\unicode[STIX]{x2202}Z}\boldsymbol{Z}\right], & \displaystyle\end{eqnarray}$$

(2.37) $$\begin{eqnarray}\displaystyle & \displaystyle p=-2\unicode[STIX]{x1D707}\mathop{\sum }_{k=-\infty }^{\infty }\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D6F1}_{k}}{\unicode[STIX]{x2202}Z^{2}}, & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D6F1}_{k}$ , $\unicode[STIX]{x1D6F9}_{k}$ and $\unicode[STIX]{x1D6FA}_{k}$ are cylindrical harmonic functions. In our previous work (Yao, Marcos & Wong Reference Yao and Wong2017), we have developed a framework to obtain such cylindrical harmonic functions with the following standard forms:

(2.38) $$\begin{eqnarray}M_{k}=\frac{1}{\unicode[STIX]{x03C0}}\left(\begin{array}{@{}l@{}}\displaystyle \cos k\unicode[STIX]{x1D6F7}\int _{0}^{\infty }\frac{1}{\unicode[STIX]{x1D706}}m_{ak}I_{k}(\unicode[STIX]{x1D706}R)\sin \unicode[STIX]{x1D706}Z\,\text{d}\unicode[STIX]{x1D706}\\ \quad +\displaystyle \sin k\unicode[STIX]{x1D6F7}\int _{0}^{\infty }\frac{1}{\unicode[STIX]{x1D706}}m_{bk}I_{k}(\unicode[STIX]{x1D706}R)\sin \unicode[STIX]{x1D706}Z\,\text{d}\unicode[STIX]{x1D706}\\ \quad +\displaystyle \cos k\unicode[STIX]{x1D6F7}\int _{0}^{\infty }\frac{1}{\unicode[STIX]{x1D706}}m_{ck}I_{k}(\unicode[STIX]{x1D706}R)\cos \unicode[STIX]{x1D706}Z\,\text{d}\unicode[STIX]{x1D706}\\ \quad +\displaystyle \sin k\unicode[STIX]{x1D6F7}\int _{0}^{\infty }\frac{1}{\unicode[STIX]{x1D706}}m_{dk}I_{k}(\unicode[STIX]{x1D706}R)\cos \unicode[STIX]{x1D706}Z\,\text{d}\unicode[STIX]{x1D706}\end{array}\right),\end{eqnarray}$$

where $M$ refers to $\unicode[STIX]{x1D6F1}$ , $\unicode[STIX]{x1D6F9}$ or $\unicode[STIX]{x1D6FA}$ and $m$ refers to  $\unicode[STIX]{x03C0}$ , $\unicode[STIX]{x1D713}$ or  $\unicode[STIX]{x1D714}$ . $\unicode[STIX]{x03C0}_{ik}$ , $\unicode[STIX]{x1D713}_{ik}$ and $\unicode[STIX]{x1D714}_{ik}$ are unknown functions of $\unicode[STIX]{x1D706}$ , where $i=a,b,c,d$ . Since the flow is axisymmetric and the boundary conditions are homogeneous in our study, the velocity expression will not be a function of $\unicode[STIX]{x1D6F7}$ , which requires $k=0$ and only  $\unicode[STIX]{x03C0}_{a0}$ , $\unicode[STIX]{x03C0}_{c0}$ , $\unicode[STIX]{x1D713}_{a0}$ and $\unicode[STIX]{x1D713}_{c0}$ contribute to the velocity expression. Hence, we have the following simplified forms:

(2.39) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6F1}_{0}=\frac{1}{\unicode[STIX]{x03C0}}\left(\begin{array}{@{}l@{}}\displaystyle \int _{0}^{\infty }\frac{1}{\unicode[STIX]{x1D706}}\unicode[STIX]{x03C0}_{a0}I_{0}(\unicode[STIX]{x1D706}R)\sin \unicode[STIX]{x1D706}Z\,\text{d}\unicode[STIX]{x1D706}\\ \quad +\,\displaystyle \int _{0}^{\infty }\frac{1}{\unicode[STIX]{x1D706}}\unicode[STIX]{x03C0}_{c0}I_{0}(\unicode[STIX]{x1D706}R)\cos \unicode[STIX]{x1D706}Z\,\text{d}\unicode[STIX]{x1D706}\end{array}\right), & \displaystyle\end{eqnarray}$$

(2.40) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6F9}_{0}=\frac{1}{\unicode[STIX]{x03C0}}\left(\begin{array}{@{}l@{}}\displaystyle \int _{0}^{\infty }\frac{1}{\unicode[STIX]{x1D706}}\unicode[STIX]{x1D713}_{a0}I_{0}(\unicode[STIX]{x1D706}R)\sin \unicode[STIX]{x1D706}Z\,\text{d}\unicode[STIX]{x1D706}\\ \quad +\,\displaystyle \int _{0}^{\infty }\frac{1}{\unicode[STIX]{x1D706}}\unicode[STIX]{x1D713}_{c0}I_{0}(\unicode[STIX]{x1D706}R)\cos \unicode[STIX]{x1D706}Z\,\text{d}\unicode[STIX]{x1D706}\end{array}\right). & \displaystyle\end{eqnarray}$$

By expanding (2.36), the flow field can be described by

(2.41) $$\begin{eqnarray}\boldsymbol{v}=\displaystyle \left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6F9}_{0}}{\unicode[STIX]{x2202}R}+R\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D6F1}_{0}}{\unicode[STIX]{x2202}R^{2}}\right)\boldsymbol{R}+\displaystyle \left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6F9}_{0}}{\unicode[STIX]{x2202}Z}+R\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D6F1}_{0}}{\unicode[STIX]{x2202}R\unicode[STIX]{x2202}Z}+\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6F1}_{0}}{\unicode[STIX]{x2202}Z}\right)\boldsymbol{Z}.\end{eqnarray}$$

We note the boundary condition at the wall of the tube $(R=R_{0})$ has the following form:

(2.42) $$\begin{eqnarray}\boldsymbol{V}(R_{0},Z)=V_{R}\boldsymbol{R}+V_{Z}\boldsymbol{Z},\end{eqnarray}$$

where

(2.43) $$\begin{eqnarray}\displaystyle & V_{R}=\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6F9}_{0}}{\unicode[STIX]{x2202}R}+R_{0}\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D6F1}_{0}}{\unicode[STIX]{x2202}R^{2}}, & \displaystyle\end{eqnarray}$$

(2.44) $$\begin{eqnarray}\displaystyle & V_{Z}=\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6F9}_{0}}{\unicode[STIX]{x2202}Z}+R_{0}\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D6F1}_{0}}{\unicode[STIX]{x2202}R\unicode[STIX]{x2202}Z}+\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6F1}_{0}}{\unicode[STIX]{x2202}Z}. & \displaystyle\end{eqnarray}$$

We now proceed to construct equations using the boundary conditions to solve the unknown coefficients in (2.39)–(2.40). By substituting (2.39)–(2.40), (2.43) can be written as

(2.45) $$\begin{eqnarray}\displaystyle V_{R}=\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6F9}_{0}}{\unicode[STIX]{x2202}R}+R_{0}\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D6F1}_{0}}{\unicode[STIX]{x2202}R^{2}}=\frac{1}{\unicode[STIX]{x03C0}}\left\{\begin{array}{@{}l@{}}\displaystyle \int _{0}^{\infty }[\unicode[STIX]{x1D713}_{c0}{I_{0}}^{\prime }(\unicode[STIX]{x1D706}R_{0})+\unicode[STIX]{x03C0}_{c0}\unicode[STIX]{x1D706}R_{0}{I_{0}}^{\prime \prime }(\unicode[STIX]{x1D706}R_{0})]\cos \unicode[STIX]{x1D706}Z\,\text{d}\unicode[STIX]{x1D706}\\ \quad +\,\displaystyle \int _{0}^{\infty }[\unicode[STIX]{x1D713}_{a0}{I_{0}}^{\prime }(\unicode[STIX]{x1D706}R_{0})+\unicode[STIX]{x03C0}_{a0}\unicode[STIX]{x1D706}R_{0}{I_{0}}^{\prime \prime }(\unicode[STIX]{x1D706}R_{0})]\sin \unicode[STIX]{x1D706}Z\,\text{d}\unicode[STIX]{x1D706}\end{array}\right\}.\end{eqnarray}$$

By utilizing a Fourier integral, we can eliminate variable $Z$ and (2.45) can be written as

(2.46) $$\begin{eqnarray}\displaystyle & \displaystyle \int _{-\infty }^{\infty }V_{R}\cos \unicode[STIX]{x1D706}Z\,\text{d}Z=\unicode[STIX]{x1D713}_{c0}{I_{0}}^{\prime }(\unicode[STIX]{x1D706}R_{0})+\unicode[STIX]{x03C0}_{c0}\unicode[STIX]{x1D706}R_{0}{I_{0}}^{\prime \prime }(\unicode[STIX]{x1D706}R_{0}), & \displaystyle\end{eqnarray}$$

(2.47) $$\begin{eqnarray}\displaystyle & \displaystyle \int _{-\infty }^{\infty }V_{R}\sin \unicode[STIX]{x1D706}Z\,\text{d}Z=\unicode[STIX]{x1D713}_{a0}{I_{0}}^{\prime }(\unicode[STIX]{x1D706}R_{0})+\unicode[STIX]{x03C0}_{a0}\unicode[STIX]{x1D706}R_{0}{I_{0}}^{\prime \prime }(\unicode[STIX]{x1D706}R_{0}). & \displaystyle\end{eqnarray}$$

Similarly, the following equations can be obtained from the boundary condition in the $Z$ -direction as given by (2.44):

(2.48) $$\begin{eqnarray}\displaystyle & \displaystyle \int _{-\infty }^{\infty }V_{Z}\cos \unicode[STIX]{x1D706}Z\,\text{d}Z=\unicode[STIX]{x1D713}_{a0}I_{0}(\unicode[STIX]{x1D706}R_{0})+\unicode[STIX]{x03C0}_{a0}\unicode[STIX]{x1D706}R_{0}{I_{0}}^{\prime }(\unicode[STIX]{x1D706}R_{0})+\unicode[STIX]{x03C0}_{a0}I_{0}(\unicode[STIX]{x1D706}R_{0}), & \displaystyle\end{eqnarray}$$

(2.49) $$\begin{eqnarray}\displaystyle & \displaystyle \int _{-\infty }^{\infty }V_{Z}\sin \unicode[STIX]{x1D706}Z\,\text{d}Z=-\unicode[STIX]{x1D713}_{c0}I_{0}(\unicode[STIX]{x1D706}R_{0})-\unicode[STIX]{x03C0}_{c0}\unicode[STIX]{x1D706}R_{0}{I_{0}}^{\prime }(\unicode[STIX]{x1D706}R_{0})-\unicode[STIX]{x03C0}_{c0}I_{0}(\unicode[STIX]{x1D706}R_{0}). & \displaystyle\end{eqnarray}$$

By solving (2.46)–(2.49) simultaneously, $\unicode[STIX]{x03C0}_{a0}$ , $\unicode[STIX]{x03C0}_{c0}$ , $\unicode[STIX]{x1D713}_{a0}$ and $\unicode[STIX]{x1D713}_{c0}$ can be obtained with the following forms and finally the flow field in cylindrical harmonics can be obtained using (2.41):

(2.50) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D713}_{c0}=\frac{[I_{0}(\unicode[STIX]{x1D706}R_{0})+\unicode[STIX]{x1D706}R_{0}{I_{0}}^{\prime }(\unicode[STIX]{x1D706}R_{0})]\displaystyle \int _{-\infty }^{\infty }V_{R}\cos \unicode[STIX]{x1D706}Z\,\text{d}Z+\unicode[STIX]{x1D706}R_{0}{I_{0}}^{\prime \prime }(\unicode[STIX]{x1D706}R_{0})\displaystyle \int _{-\infty }^{\infty }V_{Z}\sin \unicode[STIX]{x1D706}Z\,\text{d}Z}{\unicode[STIX]{x1D706}R_{0}{I_{0}}^{\prime \prime 2}(\unicode[STIX]{x1D706}R_{0})+I_{0}(\unicode[STIX]{x1D706}R_{0})[{I_{0}}^{\prime }(\unicode[STIX]{x1D706}R_{0})-\unicode[STIX]{x1D706}R_{0}{I_{0}}^{\prime \prime }(\unicode[STIX]{x1D706}R_{0})]} & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

(2.51) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D713}_{a0}=\frac{[I_{0}(\unicode[STIX]{x1D706}R_{0})+\unicode[STIX]{x1D706}R_{0}{I_{0}}^{\prime }(\unicode[STIX]{x1D706}R_{0})]\displaystyle \int _{-\infty }^{\infty }V_{R}\sin \unicode[STIX]{x1D706}Z\,\text{d}Z-\unicode[STIX]{x1D706}R_{0}{I_{0}}^{\prime \prime }(\unicode[STIX]{x1D706}R_{0})\displaystyle \int _{-\infty }^{\infty }V_{Z}\cos \unicode[STIX]{x1D706}Z\,\text{d}Z}{\unicode[STIX]{x1D706}R_{0}{I_{0}}^{\prime \prime 2}(\unicode[STIX]{x1D706}R_{0})+I_{0}(\unicode[STIX]{x1D706}R_{0})[{I_{0}}^{\prime }(\unicode[STIX]{x1D706}R_{0})-\unicode[STIX]{x1D706}R_{0}{I_{0}}^{\prime \prime }(\unicode[STIX]{x1D706}R_{0})]}, & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

(2.52) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x03C0}_{c0}=\frac{I_{0}(\unicode[STIX]{x1D706}R_{0})\displaystyle \int _{-\infty }^{\infty }V_{R}\cos \unicode[STIX]{x1D706}Z\,\text{d}Z+{I_{0}}^{\prime \prime }(\unicode[STIX]{x1D706}R_{0})\displaystyle \int _{-\infty }^{\infty }V_{Z}\sin \unicode[STIX]{x1D706}Z\,\text{d}Z}{-\unicode[STIX]{x1D706}R_{0}{I_{0}}^{\prime \prime 2}(\unicode[STIX]{x1D706}R_{0})+I_{0}(\unicode[STIX]{x1D706}R_{0})[-{I_{0}}^{\prime }(\unicode[STIX]{x1D706}R_{0})+\unicode[STIX]{x1D706}R_{0}{I_{0}}^{\prime \prime }(\unicode[STIX]{x1D706}R_{0})]}, & \displaystyle\end{eqnarray}$$

(2.53) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x03C0}_{a0}=\frac{I_{0}(\unicode[STIX]{x1D706}R_{0})\displaystyle \int _{-\infty }^{\infty }V_{R}\sin \unicode[STIX]{x1D706}Z\,\text{d}Z-{I_{0}}^{\prime \prime }(\unicode[STIX]{x1D706}R_{0})\displaystyle \int _{-\infty }^{\infty }V_{Z}\cos \unicode[STIX]{x1D706}Z\,\text{d}Z}{-\unicode[STIX]{x1D706}R_{0}{I_{0}}^{\prime \prime 2}(\unicode[STIX]{x1D706}R_{0})+I_{0}(\unicode[STIX]{x1D706}R_{0})[-{I_{0}}^{\prime }(\unicode[STIX]{x1D706}R_{0})+\unicode[STIX]{x1D706}R_{0}{I_{0}}^{\prime \prime }(\unicode[STIX]{x1D706}R_{0})]}. & \displaystyle\end{eqnarray}$$

Therefore, the unknown coefficients can be solved based on the boundary conditions. For example, $\boldsymbol{v}_{c}^{(3)}$ is required to make the summation $\boldsymbol{v}_{c}^{(0)}+\boldsymbol{v}_{a}^{(1)}+\boldsymbol{v}_{b}^{(2)}+\boldsymbol{v}_{c}^{(3)}$ satisfy the boundary condition of the tube (as shown in (2.16)). We first transform the expression of $\boldsymbol{v}_{a}^{(1)}$ and $\boldsymbol{v}_{b}^{(2)}$ from the spherical coordinates into the cylindrical coordinates (details on the coordinate transformation are provided in appendix A) and set $R=R_{0}$ to acquire the boundary condition of $\boldsymbol{v}_{c}^{(3)}$ . This boundary condition (2.16) is substituted into (2.50)–(2.53) and the unknown coefficients $\unicode[STIX]{x03C0}_{a0}$ , $\unicode[STIX]{x03C0}_{c0}$ , $\unicode[STIX]{x1D713}_{a0}$ and $\unicode[STIX]{x1D713}_{c0}$ can be solved. Similarly, other reflections with subscript ‘ $c$ ’ can be solved by using the corresponding boundary conditions satisfying the wall of the tube.

2.2.4 Drag force and velocities of particles

After obtaining the velocity fields using both stream functions and cylindrical harmonics for the particles and the cylinder, respectively, the total velocity flow field can then be obtained using (2.10)–(2.12). However, only the flow fields derived from the stream functions (the reflection terms with subscripts  $a$ and  $b$ ) contribute to the drag forces acting on the two particles (Happel & Brenner Reference Happel and Brenner1983). Therefore, the drag force on the two particles follows:

(2.54) $$\begin{eqnarray}\displaystyle & \displaystyle F_{a}=F_{a}^{(1)}+F_{a}^{(4)}+F_{a}^{(7)}+F_{a}^{(10)}+\cdots , & \displaystyle\end{eqnarray}$$

(2.55) $$\begin{eqnarray}\displaystyle & F_{b}=F_{b}^{(2)}+F_{b}^{(5)}+F_{b}^{(8)}+F_{b}^{(11)}+\cdots . & \displaystyle\end{eqnarray}$$

On the other hand, the drag force acting on the surface of a porous particle by the external fluid can be obtained from the stream function (Happel & Brenner Reference Happel and Brenner1983):

(2.56) $$\begin{eqnarray}F=\unicode[STIX]{x1D707}\unicode[STIX]{x03C0}\int _{0}^{\unicode[STIX]{x03C0}}r^{3}\sin ^{3}\unicode[STIX]{x1D703}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}\left(\frac{E^{2}\unicode[STIX]{x1D6F7}^{e}}{r^{2}\sin ^{2}\unicode[STIX]{x1D703}}\right)\text{d}\unicode[STIX]{x1D703},\end{eqnarray}$$

where $E^{2}$ is the axisymmetric Stokesian operator

(2.57) $$\begin{eqnarray}E^{2}=\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}r^{2}}+\frac{\sin \unicode[STIX]{x1D703}}{r^{2}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}\left(\frac{1}{\sin \unicode[STIX]{x1D703}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}\right).\end{eqnarray}$$

By substituting the stream function given by (2.19) into (2.56), we can obtain the simple relation (Saad & Faltas Reference Saad and Faltas2014)

(2.58) $$\begin{eqnarray}F=-4\unicode[STIX]{x1D707}\unicode[STIX]{x03C0}B_{2},\end{eqnarray}$$

where $B_{2}$ is the only remaining coefficient in the stream function contributing to the drag force due to the orthogonality of the Gegenbauer polynomial.

We can now apply (2.56) to each term in (2.54)–(2.55) and solve the corresponding coefficient $B_{2}$ using the boundary conditions as presented in (2.13)–(2.18). Since the boundary conditions given by (2.13)–(2.18) are functions of particle velocities  $U_{a}$ , $U_{b}$ and background flow  $U_{0}$ , the drag force will also be functions of  $U_{a}$ , $U_{b}$ and  $U_{0}$ . Therefore, the drag force (2.54)–(2.55) can also be expressed as

(2.59a,b ) $$\begin{eqnarray}\boldsymbol{F}_{a}=6\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}a(D_{a}U_{0}+K_{a}U_{a}+K_{ab}U_{b})\boldsymbol{Z},\quad \boldsymbol{F}_{b}=6\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}a(D_{b}U_{0}+K_{ba}U_{a}+K_{b}U_{b})\boldsymbol{Z},\end{eqnarray}$$

where $\boldsymbol{F}_{a}$ and $\boldsymbol{F}_{b}$ are the drag forces on particles  $A$ and  $B$ ; $D_{a}$ , $D_{b}$ and  $K_{a}$ , $K_{b}$ are the coefficients related to the background flow and particle velocities, respectively; $K_{ab}$ and $K_{ba}$ are the coefficients referring to hydrodynamic interaction between the two particles.

At low Reynolds number, the inertial forces are negligible, thus the total drag force on the particle is zero and using (2.58), the velocities of particles  $U_{a}$ , $U_{b}$ can be obtained. Subsequently, the total flow field can be obtained by substituting  $U_{a}$ , $U_{b}$ into (2.10)–(2.12). The practical implementation is numerically conducted for the given geometrical parameters, including particle radii, permeability, cylinder radius and separation distance. We can calculate any arbitrary number of reflections in our solution when needed for convergence. As described earlier, the calculation procedure for the method of reflections is conducted based on an iterative process by adding more terms. We terminate the calculation when the solution reaches a convergence, i.e. when the magnitude of the last term is three orders smaller than the sum of the previous terms (details of the convergence are provided in appendix B). We provide a detailed example on how to conduct the calculation for the first four reflections in appendix C.

3.1 Approximate formula for drag force

Since the drag force given by (2.58) only includes coefficient  $B_{2}$ , we approximate the stream function outside the porous particle (2.19) as

(3.1) $$\begin{eqnarray}\unicode[STIX]{x1D6F7}^{e}=\displaystyle rB_{2}C_{2}(\cos \unicode[STIX]{x1D703}).\end{eqnarray}$$

By substituting (3.1) into (2.56) and using (2.14), the first reflection of drag force acting on particle  $A$ (i.e. the first term in (2.54)) can be obtained:

(3.2) $$\begin{eqnarray}\boldsymbol{F}_{a}^{(1)}=6\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}a(D_{a}^{(1)}U_{0}+K_{a}^{(1)}U_{a})\boldsymbol{Z},\end{eqnarray}$$

where

(3.3) $$\begin{eqnarray}\displaystyle & \displaystyle D_{a}^{(1)}=\displaystyle \frac{2a^{2}[a\unicode[STIX]{x1D6FC}(3R_{0}^{2}\unicode[STIX]{x1D6FC}^{2}-2a^{2}\unicode[STIX]{x1D6FC}^{2}-12)\cosh a\unicode[STIX]{x1D6FC}+(6a^{2}\unicode[STIX]{x1D6FC}^{2}+12-3R_{0}^{2}\unicode[STIX]{x1D6FC}^{2})\sinh a\unicode[STIX]{x1D6FC}]}{3R_{0}^{2}[a\unicode[STIX]{x1D6FC}(3+2a^{2}\unicode[STIX]{x1D6FC}^{2})\cosh a\unicode[STIX]{x1D6FC}-3\sinh a\unicode[STIX]{x1D6FC}]},\qquad & \displaystyle\end{eqnarray}$$

(3.4) $$\begin{eqnarray}\displaystyle & \displaystyle K_{a}^{(1)}=\displaystyle \frac{2a^{2}\unicode[STIX]{x1D6FC}^{2}(\sinh a\unicode[STIX]{x1D6FC}-a\unicode[STIX]{x1D6FC}\cosh a\unicode[STIX]{x1D6FC})}{a\unicode[STIX]{x1D6FC}(3+2a^{2}\unicode[STIX]{x1D6FC}^{2})\cosh a\unicode[STIX]{x1D6FC}-3\sinh a\unicode[STIX]{x1D6FC}}.\qquad & \displaystyle\end{eqnarray}$$

Here, $K_{ab}^{(1)}$ is zero, since (2.14) has not included the boundary condition of particle  $B$ . Similarly, by using (2.15), the first reflection of drag force on particle  $B$ (i.e. the first term in (2.55)) is given by

(3.5) $$\begin{eqnarray}\boldsymbol{F}_{b}^{(2)}=6\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}a(D_{b}^{(2)}U_{0}+K_{ba}^{(2)}U_{a}+K_{b}^{(2)}U_{b})\boldsymbol{Z},\end{eqnarray}$$

where

(3.6) $$\begin{eqnarray}\displaystyle & & \displaystyle \!\!\!\!\!\displaystyle D_{b}^{(2)}=\nonumber\\ \displaystyle & & \displaystyle \!\!\!\!\!\quad \displaystyle \frac{\displaystyle \left\{\begin{array}{@{}l@{}}\displaystyle \left[\begin{array}{@{}l@{}}3a^{4}b^{4}\unicode[STIX]{x1D6FC}^{4}(6+b^{2}\unicode[STIX]{x1D6FC}^{2}-3L^{2}\unicode[STIX]{x1D6FC}^{2})\\ (12+8a^{2}\unicode[STIX]{x1D6FC}^{2}-3R_{0}^{2}\unicode[STIX]{x1D6FC}^{2}+2a^{4}\unicode[STIX]{x1D6FC}^{4}-a^{2}R_{0}^{2}\unicode[STIX]{x1D6FC}^{4})\\ \quad +\,3ab^{4}L^{3}\unicode[STIX]{x1D6FC}^{4}(3+2a^{2}\unicode[STIX]{x1D6FC}^{2})(12+2b^{2}\unicode[STIX]{x1D6FC}^{2}-3R_{0}^{2}\unicode[STIX]{x1D6FC}^{2})\end{array}\right]\!\cosh a\unicode[STIX]{x1D6FC}\cosh b\unicode[STIX]{x1D6FC}\sinh a\unicode[STIX]{x1D6FC}\\ \quad +\,\displaystyle \left[\begin{array}{@{}l@{}}-4.5a^{3}b^{4}\unicode[STIX]{x1D6FC}^{3}(6+b^{2}\unicode[STIX]{x1D6FC}^{2}-3L^{2}\unicode[STIX]{x1D6FC}^{2})(4+2a^{2}\unicode[STIX]{x1D6FC}^{2}-R_{0}^{2}\unicode[STIX]{x1D6FC}^{2})\\ \quad +\,L^{3}(-9b^{6}\unicode[STIX]{x1D6FC}^{5}-54b^{4}\unicode[STIX]{x1D6FC}^{3}+13.5b^{4}R_{0}^{2}\unicode[STIX]{x1D6FC}^{5})\end{array}\right]\!\cosh b\unicode[STIX]{x1D6FC}\sinh ^{2}a\unicode[STIX]{x1D6FC}\\ \quad +\,\displaystyle \left[\begin{array}{@{}l@{}}13.5a^{3}b^{3}\unicode[STIX]{x1D6FC}^{2}(2+b^{2}\unicode[STIX]{x1D6FC}^{2}-L^{2}\unicode[STIX]{x1D6FC}^{2})(4+2a^{2}\unicode[STIX]{x1D6FC}^{2}-R_{0}^{2}\unicode[STIX]{x1D6FC}^{2})\\ \quad +\,L^{3}(27b^{5}\unicode[STIX]{x1D6FC}^{4}+54b^{3}\unicode[STIX]{x1D6FC}^{2}-13.5b^{3}R_{0}^{2}\unicode[STIX]{x1D6FC}^{4})\end{array}\right]\!\sinh b\unicode[STIX]{x1D6FC}\sinh ^{2}a\unicode[STIX]{x1D6FC}\\ \quad -\,\displaystyle \left[\begin{array}{@{}l@{}}4.5a^{4}b^{3}\unicode[STIX]{x1D6FC}^{3}(2+b^{2}\unicode[STIX]{x1D6FC}^{2}-L^{2}\unicode[STIX]{x1D6FC}^{2})\\ (12-3R_{0}^{2}\unicode[STIX]{x1D6FC}^{2}+2a^{4}\unicode[STIX]{x1D6FC}^{4}+8a^{2}\unicode[STIX]{x1D6FC}^{2}-a^{2}R_{0}^{2}\unicode[STIX]{x1D6FC}^{4})\\ \quad +\,4.5ab^{3}L^{3}\unicode[STIX]{x1D6FC}^{3}(3+2a^{2}\unicode[STIX]{x1D6FC}^{2})(4+2b^{2}\unicode[STIX]{x1D6FC}^{2}-R_{0}^{2}\unicode[STIX]{x1D6FC}^{2})\end{array}\right]\!\sinh b\unicode[STIX]{x1D6FC}\sinh 2a\unicode[STIX]{x1D6FC}\\ \quad -\,\displaystyle \left[\begin{array}{@{}l@{}}0.5a^{5}b^{4}\unicode[STIX]{x1D6FC}^{5}(3+2a^{2}\unicode[STIX]{x1D6FC}^{2})(6+b^{2}\unicode[STIX]{x1D6FC}^{2}-3L^{2}\unicode[STIX]{x1D6FC}^{2})\\ (12+2a^{2}\unicode[STIX]{x1D6FC}^{2}-3R_{0}^{2}\unicode[STIX]{x1D6FC}^{2})\\ \quad +\,0.5a^{2}b^{4}L^{3}\unicode[STIX]{x1D6FC}^{5}(3+2a^{2}\unicode[STIX]{x1D6FC}^{2})^{2}(12+2b^{2}\unicode[STIX]{x1D6FC}^{2}-3R_{0}^{2}\unicode[STIX]{x1D6FC}^{2})\end{array}\right]\!\cosh b\unicode[STIX]{x1D6FC}\cosh ^{2}a\unicode[STIX]{x1D6FC}\\ \quad +\,\displaystyle \left[\begin{array}{@{}l@{}}1.5a^{5}b^{3}\unicode[STIX]{x1D6FC}^{4}(3+2a^{2}\unicode[STIX]{x1D6FC}^{2})(2+b^{2}\unicode[STIX]{x1D6FC}^{2}-L^{2}\unicode[STIX]{x1D6FC}^{2})\\ (12+2a^{2}\unicode[STIX]{x1D6FC}^{2}-3R_{0}^{2}\unicode[STIX]{x1D6FC}^{2})\\ \quad +\,1.5a^{2}b^{3}L^{3}\unicode[STIX]{x1D6FC}^{4}(3+2a^{2}\unicode[STIX]{x1D6FC}^{2})^{2}(4+2b^{2}\unicode[STIX]{x1D6FC}^{2}-R_{0}^{2}\unicode[STIX]{x1D6FC}^{2})\end{array}\right]\!\sinh b\unicode[STIX]{x1D6FC}\cosh ^{2}a\unicode[STIX]{x1D6FC}\end{array}\right\}}{\displaystyle 3a\unicode[STIX]{x1D6FC}^{2}L^{3}R_{0}^{2}[(3b\unicode[STIX]{x1D6FC}+2b^{3}\unicode[STIX]{x1D6FC}^{3})\cosh b\unicode[STIX]{x1D6FC}-3\sinh b\unicode[STIX]{x1D6FC}][(1.5a\unicode[STIX]{x1D6FC}+a^{3}\unicode[STIX]{x1D6FC}^{3})\cosh a\unicode[STIX]{x1D6FC}-1.5\sinh a\unicode[STIX]{x1D6FC}]^{2}},\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

(3.7) $$\begin{eqnarray}\displaystyle & & \displaystyle \!\!\!\!\!K_{ba}^{(2)}=\nonumber\\ \displaystyle & & \displaystyle \!\!\!\!\!\quad \displaystyle \displaystyle \frac{\displaystyle 2a^{3}\unicode[STIX]{x1D6FC}^{2}b^{3}\displaystyle \!\left(\begin{array}{@{}l@{}}3a^{2}\unicode[STIX]{x1D6FC}^{2}\cosh ^{2}a\unicode[STIX]{x1D6FC}+2a^{4}\unicode[STIX]{x1D6FC}^{4}\cosh ^{2}a\unicode[STIX]{x1D6FC}\\ \quad +\,3\sinh ^{2}a\unicode[STIX]{x1D6FC}-3a\unicode[STIX]{x1D6FC}\sinh 2a\unicode[STIX]{x1D6FC}\\ \quad -\,a^{3}\unicode[STIX]{x1D6FC}^{3}\sinh 2a\unicode[STIX]{x1D6FC}\end{array}\right)\!\displaystyle \left(\begin{array}{@{}l@{}}6b\unicode[STIX]{x1D6FC}\cosh b\unicode[STIX]{x1D6FC}-5b^{3}\unicode[STIX]{x1D6FC}^{3}\cosh b\unicode[STIX]{x1D6FC}\\ \quad +\,15bL^{2}\unicode[STIX]{x1D6FC}^{3}\cosh b\unicode[STIX]{x1D6FC}-6\sinh b\unicode[STIX]{x1D6FC}\\ \quad +\,3b^{2}\unicode[STIX]{x1D6FC}^{2}\sinh b\unicode[STIX]{x1D6FC}-15L^{2}\unicode[STIX]{x1D6FC}^{2}\sinh b\unicode[STIX]{x1D6FC}\end{array}\right)}{\displaystyle 5aL^{3}[(3a\unicode[STIX]{x1D6FC}+2a^{3}\unicode[STIX]{x1D6FC}^{3})\cosh a\unicode[STIX]{x1D6FC}-3\sinh a\unicode[STIX]{x1D6FC}]^{2}[b\unicode[STIX]{x1D6FC}(3+2b^{2}\unicode[STIX]{x1D6FC}^{2})\cosh b\unicode[STIX]{x1D6FC}-3\sinh b\unicode[STIX]{x1D6FC}]},\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

(3.8) $$\begin{eqnarray}\displaystyle & \displaystyle K_{b}^{(2)}=\displaystyle \frac{2b^{3}\unicode[STIX]{x1D6FC}^{2}(\sinh b\unicode[STIX]{x1D6FC}-b\unicode[STIX]{x1D6FC}\cosh b\unicode[STIX]{x1D6FC})}{(3ab\unicode[STIX]{x1D6FC}+2ab^{3}\unicode[STIX]{x1D6FC}^{3})\cosh b\unicode[STIX]{x1D6FC}-3a\sinh b\unicode[STIX]{x1D6FC}}. & \displaystyle\end{eqnarray}$$

Thus far, we have derived the first term of (2.55), given by (3.5), as the approximate drag force to the first leading order acting on particle  $B$ , in which the boundary conditions of particle  $A$ , particle  $B$ and the cylindrical tube have been taken into account. Following the same procedure as above, the approximate drag acting on particle  $A$ is simply (3.5) by exchanging the letters $b$ and  $a$ .

3.2 Model validation

Figure 2.

Drag force coefficients (a)  $D_{a}$ or $D_{b}$ and (b)  $K_{a}$ or $K_{b}$ when $a/R_{0}=0.1$ , $b/R_{0}=0.2$ . (c)  $|K_{a}+K_{ab}|$ when $a/R_{0}=b/R_{0}=0.1$ . (d)  $|K_{a}-K_{ab}|$ when $a/R_{0}=b/R_{0}=0.01$ . All the drag force coefficients are plotted as a function of the normalized distance $L/R_{0}$ .

In this section, we validate our semi-analytical solution with previous works, where the motion of two solid particles has been studied. As $\unicode[STIX]{x1D6FC}$ goes to infinity, the permeability approaches zero and the particle becomes a solid sphere. We first validate our full model solution with our previous work (Yao et al. Reference Yao and Wong2017), which is also based on the method of reflections. We consider the case $a/R_{0}=0.1$ , $b/R_{0}=0.2$ and $\unicode[STIX]{x1D6FC}R_{0}=50\,000$ by using the procedures introduced in § 2. At such a large  $\unicode[STIX]{x1D6FC}$ , the fluid going through the internal region is extremely limited and can be approximated as the solid particle case. Figure 2(a,b) shows that the drag force coefficients due to the background flow ( $D_{a}$ and  $D_{b}$ ) and particle velocities ( $K_{a}$ and  $K_{b}$ ) evaluated at $\unicode[STIX]{x1D6FC}R_{0}=50\,000$ are in good agreement with our previous work (Yao et al. Reference Yao and Wong2017). This has confirmed the consistency of the two frameworks. Moreover, Sonshine & Brenner (Reference Sonshine and Brenner1966) theoretically investigated the motion of two or more solid particles travelling along the axis of a cylindrical tube without the background flow by developing the point force formula based on the stream function. The particles are equally spaced along the axis and translate with a constant velocity. We compare our solution with their work for the case that ( $a/R_{0}=b/R_{0}=0.1$ , $\unicode[STIX]{x1D6FC}R_{0}=50\,000$ ) (figure 2 c). We can clearly see that the drag force coefficients calculated from our full semi-analytical results are in good agreement with the previous work except when the two particles are very close ( $L/R_{0}<0.3$ ). In addition, the drag force of two identical particles approaching each other with a constant velocity in an unbounded fluid based on the lubrication theory has been previously solved (Jeffrey Reference Jeffrey1982). We then consider the case $a/R_{0}=b/R_{0}=0.01$ and $\unicode[STIX]{x1D6FC}R_{0}=50\,000$ by using the procedures introduced in § 2 to compare with the lubrication theory solution. In such a case of a small particle size, the wall effect in our solution is negligible. It should be noted that when the two particles are very close, our solution cannot match with the lubrication theory solution (figure 2 d). This could be due to the fact that the method of reflections is not accurate to solve cases when the two solid particles ( $\unicode[STIX]{x1D6FC}R_{0}=50\,000$ ) are very close, consistent with what is reported in the literature (Bhattacharya et al. Reference Bhattacharya, Mishra and Bhattacharya2010).

Figure 3.

(a) Drag force coefficient $D_{b}$ as a function of normalized distance $L/R_{0}$ for $\unicode[STIX]{x1D6FC}R_{0}=5,10,20$ . (b) Relative error as a function of $\unicode[STIX]{x1D6FC}R_{0}$ .

Here, we present the validation of the approximate solution derived in § 3.1. In figure 3(a), we compare our full semi-analytical results with the approximate solution given by (3.6) when $a/R_{0}=0.1$ , $b/R_{0}=0.2$ and $\unicode[STIX]{x1D6FC}R_{0}=5,10,20$ . As $\unicode[STIX]{x1D6FC}R_{0}$ increases, $D_{b}$ increases, which indicates the porous particle experiences a larger drag force due to the background flow. Since (3.6) is the approximate expression of the drag force, it is not unexpected that the magnitude is smaller than that of the full semi-analytical solution. As $\unicode[STIX]{x1D6FC}R_{0}$ decreases, the difference between (3.6) and the semi-analytical result becomes smaller (figure 3 a), which implies that the particles’ contribution to the drag force is not as important. To better highlight this, we investigate the effect of $\unicode[STIX]{x1D6FC}R_{0}$ on the relative error between the full semi-analytical solution and the approximate solution for $a/R_{0}=0.1$ , $b/R_{0}=0.2$ and $L/R_{0}=0.301$ . Figure 3(b) demonstrates that the decrease in $\unicode[STIX]{x1D6FC}R_{0}$ leads to the decrease in the relative error. When $\unicode[STIX]{x1D6FC}R_{0}$ is smaller than 5, the relative error is less than 10 % which indicates that (3.6) is a good approximate solution of the drag force acting on the porous particle. We repeated our calculations for (3.7)–(3.8) and found again that these approximate solutions are valid when $\unicode[STIX]{x1D6FC}R_{0}$ is smaller than 5. We can also infer that the particle–particle interaction decreases as $\unicode[STIX]{x1D6FC}R_{0}$ decreases.

3.3 Two identical particles moving along the axis of the tube

Previous works have investigated the permeability of porous media theoretically and experimentally. Kapur et al. (Reference Kapur, Charkoudian, Kessler and Anderson1996) measured the permeability of polymer gels when they are stabilized within a microporous support. Mika & Childs (Reference Mika and Childs2001) developed a model in which gel networks are considered as an assemblage of spheres to calculate the hydrodynamic permeability. They have provided the equation to calculated the radius $a_{p}$ of such spheres. Theoretical works have provided the equations to calculate the hydrodynamic permeability based on $a_{p}$ and porosity  $\unicode[STIX]{x1D700}$ . One of the most widely used equations to calculate the permeability is from Brinkman (Reference Brinkman1949):

(3.9) $$\begin{eqnarray}k=\displaystyle \frac{a_{p}^{2}}{18}\left(3+\frac{4}{1-\unicode[STIX]{x1D700}}-3\sqrt{\frac{8}{1-\unicode[STIX]{x1D700}}-3}\right).\end{eqnarray}$$

Happel & Brenner (Reference Happel and Brenner1983) have also derived the relationship based on the cell model method (Mika & Childs Reference Mika and Childs2001):

(3.10) $$\begin{eqnarray}k=\displaystyle \frac{2a_{p}^{2}}{9(1-\unicode[STIX]{x1D700})}\left[\frac{3-4.5(1-\unicode[STIX]{x1D700})^{1/3}+4.5(1-\unicode[STIX]{x1D700})^{5/3}-3(1-\unicode[STIX]{x1D700})^{2}}{3+2(1-\unicode[STIX]{x1D700})^{5/3}}\right].\end{eqnarray}$$

Mika & Childs (Reference Mika and Childs2001) have shown the comparison between the theoretical work and experimental data conducted by Kapur et al. (Reference Kapur, Charkoudian, Kessler and Anderson1996) and both equations (3.9)–(3.10) give a good prediction of the permeability of the gel. In addition, the permeability can be predicted with the Carman–Kozeny formula (Krainov, Smirnov & Tereshonok Reference Krainov, Smirnov and Tereshonok2014):

(3.11) $$\begin{eqnarray}k=\displaystyle \frac{\unicode[STIX]{x1D700}^{3}a_{p}^{2}}{45(1-\unicode[STIX]{x1D700})^{2}}.\end{eqnarray}$$

In our study, we consider a material with a large porosity, such as a polymer gel with porosity $\unicode[STIX]{x1D700}=0.96$ . Based on the equation and parameters from Mika & Childs (Reference Mika and Childs2001), we can calculate $a_{p}=2.66\times 10^{-10}/\sqrt{0.372(1-\unicode[STIX]{x1D700})}$  m. Using (3.9)–(3.11), $\unicode[STIX]{x1D6FC}=\sqrt{1/k}$ can be obtained. In our main discussion, we choose $\unicode[STIX]{x1D6FC}R_{0}=40$ and $\unicode[STIX]{x1D6FC}R_{0}=100$ . Using these values, we summarized the calculated values of  $a_{p}$ , $\unicode[STIX]{x1D6FC}$ and $R_{0}$ in tables 1 and 2. The experiments of Kapur et al. (Reference Kapur, Charkoudian, Kessler and Anderson1996) reported a value of $\unicode[STIX]{x1D6FC}=2.26\times 10^{8}~\text{m}^{-1}$ . Using this value, we also calculated the corresponding $R_{0}$ in tables 1 and 2.

Table 1.

Variables corresponding to $\unicode[STIX]{x1D6FC}R_{0}=\sqrt{1/k}R_{0}=40$ .

Table 2.

Variables corresponding to $\unicode[STIX]{x1D6FC}R_{0}=\sqrt{1/k}R_{0}=100$ .

It can be seen that $R_{0}$ is at least one order larger than  $a_{p}$ . In our study, the minimum size of the porous particle is $a/R_{0}=0.1$ and $a$ is at least six times larger than  $a_{p}$ . Although the Brinkman equation is only fully accurate for $a_{p}/a\ll 1$ , it can be used as a first-order approximation in $a_{p}/a$ provided that $a_{p}/a<1$ (Auriault Reference Auriault2009). It should be noted that $R_{0}$ is in the nanometer scale based on our calculation result. This nanoscale of fluidic structures is of interest for a great variety of applications, ranging from drug delivery (Sinha et al. Reference Sinha, Valco, Sharma, Liu and Ferrari2004) and energy transfer (Daiguji et al. Reference Daiguji, Yang, Szeri and Majumdar2004; van der Heyden et al. Reference van der Heyden, Bonthuis, Stein, Meyer and Dekker2006) to biochemical analysis (Foquet et al. Reference Foquet, Korlach, Zipfel, Webb and Craighead2002). Previous works have shown the procedure to fabricate the nanochannel. Lee et al. (Reference Lee, Yang, Myung and George2003) have developed a nanochannel fabrication technique by using chemical–mechanical polishing and thermal oxidation. They have successfully fabricated narrow channels (i.e. 10 mm long, 25 nm wide and 100 nm deep) with smooth inner surfaces. Song et al. (Reference Song, Kim, Kwon, Kim and Kim2012) have fabricated nanochannels by utilizing the crystallinity-dependent anisotropic wet etching of silicon. The dimension of the nanochannel is down to 16 nm in depth and 175 nm in width and up to 8 cm in length. Even at such a scale, the Stokes equation and Brinkman equation are still valid (Urbakh, Daikhin & Klafter Reference Urbakh, Daikhin and Klafter1995; Bocquet & Charlaix Reference Bocquet and Charlaix2010).

We first consider two identical porous particles moving along the axis of the tube ( $a/R_{0}=b/R_{0}=0.1$ and $a/R_{0}=b/R_{0}=0.3$ ) with $\unicode[STIX]{x1D6FC}R_{0}=40$ and $\unicode[STIX]{x1D6FC}R_{0}=100$ . The semi-analytical results such as drag force coefficients and particle velocities are plotted against the normalized separation distance between the two particles $L/R_{0}$ . Due to kinematic reversibility at low Reynolds number, the coefficients of the two identical particles are expected to be identical. Our semi-analytical results have also verified that this is the case. For convenience, the new variables $D(D=D_{a}=D_{b}),K(K=K_{a}=K_{b})$ and $K^{\ast }(K^{\ast }=K_{ab}=K_{ba})$ are utilized to represent the coefficients (figure 4 a–c).

Figure 4.

Drag force coefficients and normalized velocity of two identical porous particles. (a) Drag force coefficient  $D$ ; (b) drag force coefficient  $K$ ; (c) drag force coefficient  $K^{\ast }$ and (d) normalized particle velocities $U=U_{a}=U_{b}$ . Two particle sizes and two $\unicode[STIX]{x1D6FC}R_{0}$ are considered: $a/R_{0}=b/R_{0}=0.1,0.3$ ; $\unicode[STIX]{x1D6FC}R_{0}=40,100$ .

The coefficient $D$ (figure 4 a), describing the drag force due to the background flow on the particles when they are at rest, increases monotonically from a minimum to a constant value. At small separation distance $L/R_{0}$ , the effect of particle interaction is strong inducing a small velocity gradient between them and thus a small drag force. Details of the streamlines and velocity contour will be discussed below. As $L/R_{0}$ increases, the effect of particle interaction diminishes and each particle behaves as a single particle case, indicated by the constant  $D$ . In addition, bigger particles have larger drag force coefficients due to the stronger wall effect. When $\unicode[STIX]{x1D6FC}R_{0}$ is larger ( $\unicode[STIX]{x1D6FC}R_{0}=100$ ), the particle has a smaller permeability, hence it is harder for the fluid to go through the internal region of the particle. The porous particle experiences a larger drag force than that with a smaller $\unicode[STIX]{x1D6FC}R_{0}$ ( $\unicode[STIX]{x1D6FC}R_{0}=40$ ). Thus, as $\unicode[STIX]{x1D6FC}R_{0}$ increases, the coefficient  $D$ has a similar relationship with $L/R_{0}$ , but with a larger magnitude.

Unlike the coefficient $D$ , $K$ represents the drag force contribution due to the movement of one particle while the other is at rest in the absence of the background flow (figure 4 b). This drag force is against the direction of the particle motion, which explains the negative sign of  $K$ . When the particle is moving in the fluid, the presence of the neighbouring resting particle increases the velocity gradient in the vicinity of the moving particle and thus also increases the drag force acting on the moving particle. This effect is more pronounced as the two particles get closer due to stronger particle–particle interaction. As the separation distance $L/R_{0}$ increases, $K$ approaches the single particle case. In addition, smaller $\unicode[STIX]{x1D6FC}R_{0}$ allows the fluid to smoothly pass through the particle, represented by the smaller coefficient $K$ magnitude. Unsurprisingly, we also found that bigger particles have a larger drag force magnitude due to the stronger wall effect.

The drag coefficient $K^{\ast }$ represents the particle–particle interaction in the absence of background flow (figure 4 c). It describes the drag force acting on a resting particle due to the movement of the neighbouring particle. The coefficient $K^{\ast }$ is positive since the drag force is in the same direction as the moving particle. When the neighbouring particle is moving, it drags the surrounding fluid inducing drag force on the resting particle. As the two particles get very close, the coefficient $K^{\ast }$ reaches its maximum value, which arises from a large shear rate in the gap between the two particles. It is not unexpected that as $L/R_{0}$ increases, $K^{\ast }$ decays to zero since the hydrodynamic interaction becomes negligible.

As mentioned earlier, the velocities of particles can be calculated by applying zero total drag force due to the negligible inertia effect at low Reynolds number. The velocity of particle  $A$ is the same as particle  $B$ since the two particles are identical ( $U=U_{a}=U_{b}$ ), which implies the two particles will not attract nor repel. This is consistent with the kinematic reversibility argument. Figure 4(d) shows that the particle velocity is almost unaffected by the separation distance $L/R_{0}$ , with the smallest magnitude when the two particles are almost in contact. We also find that smaller $\unicode[STIX]{x1D6FC}R_{0}$ gives rise to larger particle velocity, as it is easier for the fluid to go through the particle. When the particle is bigger, it possesses a smaller velocity due to the non-uniformity of the background parabolic flow as well as the interaction with the cylinder wall.

To better highlight the hydrodynamic interaction between the two porous particles, we choose two representative separation distances ( $L/R_{0}=0.3,0.7$ ) with $a/R_{0}=b/R_{0}=0.1$ and $\unicode[STIX]{x1D6FC}R_{0}=40,100$ , to investigate the corresponding streamlines and velocity contour in the $Z$ -direction in the vicinity of the two particles (figure 5) when they are fixed. The pattern of the streamlines is related to the permeability of the particle. Large permeability (small $\unicode[STIX]{x1D6FC}R_{0}$ ) allows fluid to go through the internal region of the particle easily, and thus results in streamlines that are almost parallel to the axis of the tube (top half of figure 5 a,b). Due to higher fluid resistance to move through the porous particle, the $Z$ -direction velocity of the fluid inside the particle is smaller than that outside. When $\unicode[STIX]{x1D6FC}R_{0}=100$ , the particle is less permeable, inducing evident distortion of the streamlines in the vicinity of the two particles (bottom half of figure 5 a,b). The streamlines and velocity contour are squeezed towards the wall of the tube. The $Z$ -direction velocity of the fluid in the vicinity of the two particles is smaller than that of $\unicode[STIX]{x1D6FC}R_{0}=40$ because of the higher fluid resistance inside the particle. The velocity contour is denser when $\unicode[STIX]{x1D6FC}R_{0}$ is larger indicating larger shear rates, which explains why the particle with larger $\unicode[STIX]{x1D6FC}R_{0}$ experiences a larger drag force than that with smaller $\unicode[STIX]{x1D6FC}R_{0}$ . The effect of separation distance $L/R_{0}$ on the velocity field can also be seen from figure 5. At small $L/R_{0}$ , the boundary layers of the two particles interact and thus as the streamline exits the first particle, it does not recover to its upstream condition before entering the second particle (figure 5 a). Such distortion is not visible when the two particles are far apart (figure 5 b). The fluid in the gap between the two particles experiences smaller shear rate when they are close to each other.

Figure 5.

Streamlines (left) and normalized velocity contour in the $Z$ -direction $v_{Z}/U_{0}$ (right) around two identical particles $a/R_{0}=b/R_{0}=0.1$ . Since it is an axisymmetric model, we only plot half of the region of the flow field: $\unicode[STIX]{x1D6FC}R_{0}=40$ (top half) and $\unicode[STIX]{x1D6FC}R_{0}=100$ (bottom half). Two representative separation distances are considered: $L/R_{0}=0.3$ (a) and $L/R_{0}=0.7$ (b).

Figure 6 shows the streamlines and velocity contour in the $Z$ -direction around two identical particles when $a/R_{0}=b/R_{0}=0.3$ . We can see that the distortion of the streamlines is more significant since the two particles are bigger than those presented in figure 5. The velocity contour further reinforces this phenomenon. Here, the particles squeeze the streamlines towards the wall and the velocity contours in the vicinity around the boundary of the particle are denser than that in figure 5. As mentioned earlier, the denser velocity contour indicates a larger velocity gradient, which results in a larger drag force. The effect of $\unicode[STIX]{x1D6FC}R_{0}$ on the velocity field for the larger particles here is similar to that in figure 5. However, it is interesting to note that even though the background velocity is parabolic and the local fluid velocity moving at the axis of the tube is maximum, we find that the velocity of the fluid inside the particle is smaller than the velocity of the fluid outside the particle nearer to the wall, due to the fluid resistance inside the particle. This effect is more evident when the particle is bigger with sufficiently large $\unicode[STIX]{x1D6FC}R_{0}$ (i.e. sufficiently large fluid resistance). In the limit of $\unicode[STIX]{x1D6FC}R_{0}$ approaching zero, the particle is completely permeable (simply becomes a control volume), the velocity field will therefore be unaffected regardless of the size of the particle.

Figure 6.

Streamlines (left) and normalized velocity contour in the $Z$ -direction $v_{z}/U_{0}$ (right) around two identical particles $a/R_{0}=b/R_{0}=0.3$ . Since it is an axisymmetric model, we only plot half of the region of the flow field: $\unicode[STIX]{x1D6FC}R_{0}=40$ (top half) and $\unicode[STIX]{x1D6FC}R_{0}=100$ (bottom half). Two representative separation distances are considered: $L/R_{0}=0.7$ (a) and $L/R_{0}=1.5$ (b).

3.4 Two particles with unequal sizes moving along the axis of the tube

Figure 7.

Drag force coefficient $D_{a}$ of particle size ratio $b/a=2,5,7$ for $\unicode[STIX]{x1D6FC}R_{0}=40$ (a) and $\unicode[STIX]{x1D6FC}R_{0}=100$ (b). The size of particle  $A$ , $a/R_{0}=0.1$ . The horizontal line refers to the single particle case.

We now consider the case when the two particles are unequal. We keep the size of particle  $A$ as $a/R_{0}=0.1$ and vary the radius of particle $B$ with $b/a=2,5,7$ . Figure 7 shows the normalized drag force coefficient  $D_{a}$ . As mentioned previously, when the two particles are close to each other, a larger region where the velocity gradient of the fluid is small is induced between them that causes the smaller  $D_{a}$ . As the separation distance $L/R_{0}$ increases, $D_{a}$ approaches the single particle result. Such an effect is more pronounced with the bigger sizes of particle $B$ ( $b/a=2,5,7$ ) due to the stronger particle–particle interaction. Similar to figure 4(a), the permeability influences the drag force. When the particle’s permeability is small ( $\unicode[STIX]{x1D6FC}R_{0}=100$ , figure 7 b), it is harder for the fluid to go through the internal region of the particle resulting in a larger drag force on particle  $A$ . At larger permeability ( $\unicode[STIX]{x1D6FC}R_{0}=40$ , figure 7 a), the flow field is more parallel and the particle–particle interaction has relatively less effect on the particle’s motion, which causes smaller $D_{a}$ magnitude as well as smaller variation of $D_{a}$ against $L/R_{0}$ .

Figure 8.

Drag force coefficient $D_{b}$ of particle size ratio $b/a=2$ (a,b), $b/a=5$ (c,d), $b/a=7$ (e,f). The size of particle  $A$ , $a/R_{0}=0.1$ ; $\unicode[STIX]{x1D6FC}R_{0}=40$ (a,c,e) and $\unicode[STIX]{x1D6FC}R_{0}=100$ (b,d,f). The horizontal line refers to the single particle case.

Figure 9.

Streamlines (left) and normalized velocity contour in the $Z$ -direction $v_{z}/U_{0}$ (right) around two particles $a/R_{0}=0.1$ and $b/a=2$ . Since it is an axisymmetric model, we only plot half of the region of the flow field: $\unicode[STIX]{x1D6FC}R_{0}=40$ (top half) and $\unicode[STIX]{x1D6FC}R_{0}=100$ (bottom half). Three representative separation distances are considered: $L/R_{0}=0.301$ (a), $L/R_{0}=0.35$ (b), and $L/R_{0}=0.6$ (c). The three representative $L/R_{0}$ have been marked in figure 8(a,b).

Figure 8 shows the variation of the coefficient $D_{b}$ against separation distance $L/R_{0}$ . Interestingly, $D_{b}$ does not vary monotonously when particle $B$ is bigger than particle $A$ except for $a/R_{0}=0.1$ , $b/R_{0}=0.2$ and $\unicode[STIX]{x1D6FC}R_{0}=40$ (figure 8 a). As $L/R_{0}$ increases, the coefficient $D_{b}$ decreases first to a minimum value and then increases to a constant representing the single particle case. Similar to  $D_{a}$ , smaller permeability ( $\unicode[STIX]{x1D6FC}R_{0}=100$ ) causes the porous particle to experience a larger $D_{b}$ (figure 8 b,d,f).

To investigate the detailed flow field around the two particles, we choose three representative separation distances ( $L/R_{0}=0.301$ , 0.35 and 0.6) with $b/a=2$ when $\unicode[STIX]{x1D6FC}R_{0}=40$ and 100 and plot the corresponding streamlines and velocity contour in the $Z$ -direction in the vicinity of the two particles when they are fixed along the axis (figure 9). The permeability of the particle and separation distance $L/R_{0}$ can significantly affect the flow field. At $\unicode[STIX]{x1D6FC}R_{0}=40$ (top half), it is easier for the fluid to pass through the two particles and the streamlines are almost parallel to the axis of the tube with smaller distortion. In addition, a bigger particle can induce larger distortion and the streamlines are squeezed towards the wall when the fluid passes through particle  $B$ . As discussed earlier, due to the resistance inside the porous particle, the fluid has a smaller $Z$ -direction velocity inside the particle than that outside. This phenomenon is more evident for larger particle  $B$ . At $\unicode[STIX]{x1D6FC}R_{0}=100$ (bottom half), the particle has smaller permeability and the fluid tends to move around the particle. Especially in the vicinity of particle  $B$ , the streamlines are extremely squeezed and more fluid flows outside the particle. Furthermore, the velocity contour at $\unicode[STIX]{x1D6FC}R_{0}=100$ is stretched outwards as compared with $\unicode[STIX]{x1D6FC}R_{0}=40$ due to the larger distortion. Such an effect causes the velocity contour to become denser around the two particles, which is more apparent around particle  $B$ . Since a dense velocity contour indicates a large velocity gradient, the two particles experience a large drag force, which agrees with that previously discussed (figures 5 and 6). Taking a closer look at the velocity contour of $\unicode[STIX]{x1D6FC}R_{0}=40$ , when the two particles are almost in contact (figure 9 a top half) and when the two particles are slightly apart (figure 9 b top half), the contour looks almost identical, except that the contour is slightly denser in figure 9(b). When the separation distance $L/R_{0}$ is even larger (figure 9 c top half), the velocity contour is even denser. This explains the monotonic increase in $D_{b}$ in figure 8(a). However, for small permeability ( $\unicode[STIX]{x1D6FC}R_{0}=100$ ), the velocity contour when the two particles are almost in contact is denser than that when the two particles are slightly apart (figure 9 a,b bottom half). When the separation distance $L/R_{0}$ is large (figure 9 c bottom half), the velocity contour is denser than that in figure 9(a,b) (bottom half). The corresponding drag coefficient $D_{b}$ is therefore observed to decrease at small $L/R_{0}$ and increase at large $L/R_{0}$ (figure 8 b).

Figure 10.

Drag force coefficient $K_{a}$ of various particle size ratios $b/a=2,5,7$ with $\unicode[STIX]{x1D6FC}R_{0}=40$ (a) and $\unicode[STIX]{x1D6FC}R_{0}=100$ (b). The size of particle  $A$ , $a/R_{0}=0.1$ .

Figure 11.

Drag force coefficient $K_{b}$ of various particle size ratios $b/a=2$ (a,b), $b/a=5$ (c,d), $b/a=7$ (e,f) with $\unicode[STIX]{x1D6FC}R_{0}=40$ (a,c,e) and $\unicode[STIX]{x1D6FC}R_{0}=100$ (b,d,f). The size of particle  $A$ , $a/R_{0}=0.1$ .

The drag coefficients $K_{a}$ and $K_{b}$ show the drag force due to the movement of one particle while the other is at rest in the absence of the background flow. Figure 10 shows the variation of $K_{a}$ , which is due to the movement of particle $A$ against $L/R_{0}$ . The magnitude of $K_{a}$ decreases monotonously to the value of the single particle case. Furthermore, the bigger the resting particle  $B$ is, the larger the magnitude of $K_{a}$ due to larger velocity gradient between the two particles. Also, when $\unicode[STIX]{x1D6FC}R_{0}$ is small, the flow field is more parallel and the particle–particle interaction has relatively less effect on the particle–particle motion. Therefore, both the magnitude and the variation of $K_{a}$ against $L/R_{0}$ are significantly smaller than that of larger $\unicode[STIX]{x1D6FC}R_{0}$ (figure 10). Similar behaviour is also observed in the coefficient $K_{b}$ (figure 11).

Figure 12.

Drag force coefficient $K_{ab}$ (lines) or $K_{ba}$ (symbols) of various particle size ratios $b/a=2,5,7$ with $\unicode[STIX]{x1D6FC}R_{0}=40$ (a) and $\unicode[STIX]{x1D6FC}R_{0}=100$ (b). The size of particle  $A$ , $a/R_{0}=0.1$ .

Figure 13.

(a,c) Normalized velocity of particle  $A$ . The horizontal line refers to the velocity of the single particle case. (b,d) Normalized velocity of particle  $B$ . $\unicode[STIX]{x1D6FC}R_{0}=40$ (a,b) and $\unicode[STIX]{x1D6FC}R_{0}=100$ (c,d).

The coefficients $K_{ab}$ and $K_{ba}$ represent the particle–particle interaction. The drag force has the same direction as the moving neighbouring particle, which makes the coefficients positive (figure 12). Our semi-analytical results have validated that the two coefficients are identical, consistent with what is reported in the literature (Purcel Reference Purcel1977; Reuland, Felderhof & Jones Reference Reuland, Felderhof and Jones1978). When two particles travel at low Reynolds number, the resistance matrix is symmetric based on Lorentz reciprocal theorem (Reuland et al. Reference Reuland, Felderhof and Jones1978; Higdon & Kojima Reference Higdon and Kojima1981; Happel & Brenner Reference Happel and Brenner1983). Thus, the two coefficients $K_{ab}$ and $K_{ba}$ (off diagonal components of the matrix) are identical. When the neighbouring particle is moving, it drags the surrounding fluid inducing drag force on the resting particle. In both cases, when the two particles are very close, the coefficients $K_{ab}$ and $K_{ba}$ have the maximum value, which arises from a large shear rate in the gap between the moving and the resting particles. As the separation distance increases, the coefficient decays to zero, which indicates that the hydrodynamic interaction can be neglected. When $\unicode[STIX]{x1D6FC}R_{0}$ is larger (figure 12 b), the particle is less permeable resulting in a larger magnitude of $K_{ab}$ and  $K_{ba}$ .

We investigated the effect of the separation distance $L/R_{0}$ and radii ratios $b/a$ on the velocities of particles $U_{a}$ and $U_{b}$ (figure 13). The two particles are moving freely along the axis of the tube. We can see that the particle interaction only has a significant effect on $U_{a}$ (figure 13 a,c) and not so much on $U_{b}$ (b,d). The effect is more significant when $b/a$ is larger. However, the effect of the particle interaction on particle $B$ is weak, which can be seen from the near constant $U_{b}$ (figure 13 b,d). In addition, the permeability can influence the particle motion. When $\unicode[STIX]{x1D6FC}R_{0}=40$ , the particle is more permeable and has larger velocity than that for the case when $\unicode[STIX]{x1D6FC}R_{0}=100$ . Furthermore, $U_{a}$ is always larger than $U_{b}$ for any given $\unicode[STIX]{x1D6FC}R_{0}$ . This implies the two particles are approaching when the smaller particle  $A$ is placed behind the bigger particle  $B$ . Using the kinematic reversibility, we can deduce that the two particles will repel each other if we reverse the direction of the Poiseuille flow, or when the smaller particle is in front of the bigger particle.

3.5 Effect of permeability on the particle–particle interaction

Figure 14.

(a) Effect of permeability parameter $\unicode[STIX]{x1D6FC}R_{0}$ on characteristic separation distance $L^{\ast }/R_{0}$ . The size of particle  $A$ , $a/R_{0}=0.1$ and particle size ratios $b/a=2,3,4$ . (b) Effect of particle size ratio on the characteristic permeability parameter $\unicode[STIX]{x1D6FC}^{\ast }R_{0}$ .

Our previous work on solid spheres (Yao et al. Reference Yao and Wong2017) has found that the characteristic distance $L^{\ast }$ to neglect particle–particle interaction is three times the sum of the particles’ radii. This characteristic distance is defined as the distance when the velocity of each particle is within 1 % difference of that of the single particle case. However, this characteristic distance is expected to change, since particle’s permeability changes the flow field and correspondingly the particles’ velocities. As previously discussed, the particle–particle interaction becomes stronger as $\unicode[STIX]{x1D6FC}R_{0}$ increases. This indicates the characteristic distance should also increase with $\unicode[STIX]{x1D6FC}R_{0}$ . Moreover, if $\unicode[STIX]{x1D6FC}R_{0}$ is small enough, the particle–particle interaction can be completely neglected even when the two particles are almost touching, i.e. in this case $L/R_{0}=L^{\ast }/R_{0}$ . We can then define this particular $\unicode[STIX]{x1D6FC}R_{0}$ as the characteristic $\unicode[STIX]{x1D6FC}^{\ast }R_{0}$ . Figure 14(a) shows the relationship between the normalized characteristic distance $L^{\ast }/R_{0}$ and $\unicode[STIX]{x1D6FC}R_{0}$ . When $\unicode[STIX]{x1D6FC}R_{0}$ is smaller than 100, $L^{\ast }/R_{0}$ increases significantly followed by a gentler linear increase at larger $\unicode[STIX]{x1D6FC}R_{0}$ . When $a/R_{0}=0.1$ and $b/a=2$ , $L^{\ast }/R_{0}$ has a linear relationship with $\unicode[STIX]{x1D6FC}R_{0}$ . Increase in $\unicode[STIX]{x1D6FC}R_{0}$ , i.e. a less permeable porous particle, leads to an increase in the overall internal resistance, which finally results in stronger particle–particle interaction together with larger characteristic distance. Since the particle–particle interaction becomes stronger when $b/a$ is larger, the characteristic $\unicode[STIX]{x1D6FC}^{\ast }R_{0}$ decreases as $b/a$ increases (figure 14 b).

3.6 Two particles with different permeabilities

In this section, we study the case when the two porous particles have different permeabilities and plot the streamlines and velocity contours in the $Z$ -direction in the vicinity of the two particles when they are fixed along the axis (figure 15). When the permeability of particle $A$ ( $\unicode[STIX]{x1D6FC}R_{0}=40$ ) is smaller than that of particle $B$ ( $\unicode[STIX]{x1D6FC}R_{0}=100$ ) (top half), it is easier for the fluid to pass through particle  $A$ and the streamlines are almost parallel to the axis of the tube. The streamlines around particle  $B$ have been pushed towards the wall. This phenomenon is more evident in the velocity contours, which have more uniform distribution for particle  $A$ as compared with particle  $B$ . These features are in accord with our previous discussion. Similarly, if we switch the permeability of the two particles (bottom half), the distribution of the streamlines and velocity contours will be changed accordingly.

Figure 15.

Streamlines (a) and normalized velocity contour in the $Z$ -direction $v_{z}/U_{0}$ (b) around two particles $a/R_{0}=0.1$ , $b/a=2$ and $L/R_{0}=0.35$ with different permeabilities. Since this is an axisymmetric model, we only plot half of the region of the flow field. Two combinations of different permeabilities are considered: $\unicode[STIX]{x1D6FC}R_{0}=40$ for particle  $A$ (smaller particle), $\unicode[STIX]{x1D6FC}R_{0}=100$ for particle  $B$ (bigger particle) (top half) and $\unicode[STIX]{x1D6FC}R_{0}=100$ for particle  $A$ (smaller particle), $\unicode[STIX]{x1D6FC}R_{0}=40$ for particle $B$ (bigger particle) (bottom half).