2.1.1. Matrix-vector formulation of hydrodynamic interactions

In the Stokes flow regime, the equations of motion for particles are described by the Langevin equation, which represents the balance of forces affecting particle motion. The Stokesian dynamics is one approach that accounts for multibody hydrodynamic interactions through the effect of matrix inversion pioneered by Brady & Bossis (Reference Brady and Bossis1988). The changes in force density within a suspension inevitably lead to fluctuations at all points, and this, in turn, affects the changes in force density. This infinite sum of reflections is analogous to the inversion of the mobility matrix, which contains the information on force density (Durlofsky, Brady & Bossis Reference Durlofsky, Brady and Bossis1987). Therefore, in this section, we start from the Langevin equation and briefly describe the conventional methods for formulating multibody hydrodynamic interactions in particle suspensions using matrix-based approaches.

The number of N colloid particles suspended in the fluid that have viscosity $\mu $ and density $\rho$ satisfies the following force and torque balance (Brady Reference Brady1993):

(2.1)\begin{gather}{\boldsymbol{m}}\boldsymbol{\cdot }\frac{{\textrm{d}{\boldsymbol{U}}}}{{\textrm{d}t}} = {{\boldsymbol{F}}^H} + {{\boldsymbol{F}}^P} + {{\boldsymbol{F}}^B},\end{gather}

(2.2)\begin{gather}{{\boldsymbol{I}}_m}\boldsymbol{\cdot }\frac{{\textrm{d}{\boldsymbol{\varOmega} }}}{{\textrm{d}t}} = {{\boldsymbol{T}}^H} + {{\boldsymbol{T}}^P} + {{\boldsymbol{T}}^B}.\end{gather}

In the above equations, the left-hand sides can be ignored in the overdamped condition where the mass and relaxation time of the particles are sufficiently small. Here, ${\boldsymbol{m}}$ and ${{\boldsymbol{I}}_m}$ denote mass and moment of inertia, respectively, while ${\boldsymbol{U}}$ and ${\boldsymbol{\varOmega} }$ denote the translational and rotational velocity of the particle, respectively. Therefore, the forces and torques that make balances are hydrodynamic contributions denoted by superscript H, deterministic non-hydrodynamic force or torque with superscript $P$ and stochastic force or torque with B. Hydrodynamic force acts upon the particle due to the relative motion of the fluid with respect to the particle, and deterministic non-hydrodynamic force acts between a pair of particles based on the inter-particle potential. However, stochastic force and torque describe the thermal motion of particles through Brownian random walk and are uncorrelated with time.

To obtain the hydrodynamic contributions acting on a particle, a correlation between the velocity moments and force moments needs to be established. In the Stokes flow regime, characterized by the negligible fluid inertia, a linear relationship exists between an infinitude of higher-order velocity derivatives and hydrodynamic traction moments. Mediating this relationship is the Green's function ${\boldsymbol{J}}$, a fundamental solution of the Stokes equation. The point force ${\boldsymbol{f}}(\kern 1.5pt {\boldsymbol{y}})$, located on the surface element $\textrm{d}{S_{\boldsymbol{y}}}$ of the particle suspended in the fluid where ${\boldsymbol{y}}$ is a position vector, is integrated to determine the fluid velocity ${\boldsymbol{u}}({\boldsymbol{x}})$ at point ${\boldsymbol{x}}$:

(2.3)\begin{equation}{\boldsymbol{u}}({\boldsymbol{x}}) = {{\boldsymbol{u}}^\infty }({\boldsymbol{x}}) - \frac{1}{{8{\rm \pi} \mu }}\sum\limits_{\alpha = 1}^N {\int_{{S^\alpha }} {[{\boldsymbol{J}}({\boldsymbol{x}} - {\boldsymbol{y}})\boldsymbol{\cdot }{\boldsymbol{f}}({\,\boldsymbol{y}})]\,\textrm{d}{S_{\boldsymbol{y}}}} } .\end{equation}

The multipole expansion of the Green's function about the centre ${{\boldsymbol{x}}^\alpha }$ of each particle yields

(2.4)\begin{equation}{J_{ij}}({\boldsymbol{x}} - {\boldsymbol{y}}) = {J_{ij}}({\boldsymbol{x}} - {{\boldsymbol{x}}^\alpha }) + (x_k^\alpha - {y_k}){\nabla _k}{J_{ij}}({\boldsymbol{x}} - {{\boldsymbol{x}}^\alpha }) + \cdots .\end{equation}

The above expansion is substituted into (2.3), and after simplification through some replacements, it yields

(2.5)\begin{align}{u_i}({\boldsymbol{x}}) = u_i^\infty ({\boldsymbol{x}}) - \frac{1}{{8{\rm \pi} \mu }}\sum\limits_{\alpha = 1}^N {[{J_{ij}}({\boldsymbol{x}} - {{\boldsymbol{x}}^\alpha })F_j^H + {R_{ij}}({\boldsymbol{x}} - {{\boldsymbol{x}}^\alpha })T_j^H - {K_{ijk}}({\boldsymbol{x}} - {{\boldsymbol{x}}^\alpha })S_{jk}^H + \cdots ].}\end{align}

The replaced terms are equal to (2.6)–(2.10), while $\varepsilon$ is the Levi-Civita symbol in three dimensions. Here, ${{\boldsymbol{F}}^H}$ and ${{\boldsymbol{T}}^H}$ are the total force and torque on the particle exerted by the fluid. The propagator for the torque ${\boldsymbol{R}}$ is called the rotlet and ${{\boldsymbol{S}}^H}$ is the stresslet on the fluid exerted by the particle. In (2.5), each moment of the force density induces a velocity disturbance of the fluid through propagators. For these moments, ${{\boldsymbol{F}}^H}$ corresponds to the monopole or zeroth moment, the asymmetric part of the dipole or first moment becomes ${{\boldsymbol{T}}^H}$, while the symmetric part becomes ${{\boldsymbol{S}}^H}$. In this study, terms beyond the first moment are truncated in subsequent expansions, except for two higher multipole contributions arising from the finite size of particles. There are two main reasons for this truncation. The first reason is for computational efficiency. The second reason is that while the force moments and velocity moments are linearly related by the proportionality tensor with the same degree of correlation, in this study, background velocity moments are obtained up to the first moment through a fluid dynamics solver, thus coupling up to this degree of accuracy is considered. Therefore, higher force moments lose their significance as they would not introduce additional coupling effects beyond the first moment.

(2.6)\begin{gather}F_i^H = \int_{{S^\alpha }} {{f_i}(\kern 1.5pt {\boldsymbol{y}})\,\textrm{d}{S_y}} ,\end{gather}

(2.7)\begin{gather}{R_{ij}} ={-} \frac{1}{2}{\varepsilon _{ijk}}{\nabla _k}{J_{il}},\end{gather}

(2.8)\begin{gather}T_i^H = \int_{{S^\alpha }} {{\varepsilon _{ijk}}(\kern 1.5pt {y_j} - x_j^\alpha ){f_k}(\kern 1.5pt {\boldsymbol{y}})\,\textrm{d}{S_y}} ,\end{gather}

(2.9)\begin{gather}{K_{ijk}} = \frac{1}{2}({\nabla _k}{J_{ij}} + {\nabla _j}{J_{ik}}),\end{gather}

(2.10)\begin{gather}S_{ij}^H = \int_{{S^\alpha }} {\left[ {(\kern 1.5pt {y_i} - x_i^\alpha ){f_j} + (\kern 1.5pt {y_j} - x_j^\alpha ){f_i} - \frac{2}{3}{\delta_{ij}}(\kern 1.5pt {y_k} - x_k^\alpha ){f_k}} \right]\textrm{d}{S_y}} .\end{gather}

Assuming the spherical particles with radius a, each particle generates fluid disturbances due to its finite size in the fluid. The higher terms of the force moments that are added correspond to the quadrupole contribution $(1/6){a^2}{\nabla ^2}{J_{ij}}$ and octuplet contribution $(1/10){a^2}{\nabla ^2}{K_{ijk}}$ from straining motion (Durlofsky et al. Reference Durlofsky, Brady and Bossis1987).

(2.11)\begin{align}{u_i}({\boldsymbol{x}}) = u_i^\infty ({\boldsymbol{x}}) - \frac{1}{{8{\rm \pi} \mu }}\sum\limits_{\alpha = 1}^N {\left[ \begin{array}{@{}l@{}} \left( {1 + \dfrac{1}{6}{a^2}{\nabla^2}} \right){J_{ij}}({\boldsymbol{x}} - {{\boldsymbol{x}}^\alpha })F_j^H + {R_{ij}}({\boldsymbol{x}} - {{\boldsymbol{x}}^\alpha })T_j^H\\ \quad - \left( {1 + \dfrac{1}{{10}}{a^2}{\nabla^2}} \right){K_{ijk}}({\boldsymbol{x}} - {{\boldsymbol{x}}^\alpha })S_{jk}^H \end{array} \right]} .\end{align}

Due to the flow disturbance in the summation of (2.11) without contribution from the particle $\alpha$, denoted by ${\boldsymbol{u}^{\prime}}({{\boldsymbol{x}}^\alpha })$, the particle velocities at the centre of the particle ${{\boldsymbol{x}}^\alpha }$ are given by Faxen's formulae (2.12)–(2.14). In these equations, ${{\boldsymbol{U}}^\alpha }$ and ${{\boldsymbol{\varOmega} }^\alpha }$ refer to the translational and rotational velocities of the particle $\alpha$, respectively.

(2.12)\begin{gather}U_i^\alpha - u_i^\infty ({{\boldsymbol{x}}^\alpha }) ={-} \frac{{F_i^{H,\alpha }}}{{6{\rm \pi} \mu {a_\alpha }}} + \left( {1 + \frac{{a_\alpha^2}}{6}{\nabla^2}} \right){u^{\prime}_i}({{\boldsymbol{x}}^\alpha }),\end{gather}

(2.13)\begin{gather}\varOmega _i^\alpha - \varOmega _i^\infty ={-} \frac{{T_i^{H,\alpha }}}{{8{\rm \pi} \mu a_\alpha ^3}} + \frac{1}{2}{\varepsilon _{ijk}}{\nabla _j}{u^{\prime}_k}({{\boldsymbol{x}}^\alpha }),\end{gather}

(2.14)\begin{gather}- E_{ij}^\infty = \frac{1}{2}\left[ { - \frac{{S_{ij}^{H,\alpha }}}{{(10/3){\rm \pi} \mu a_\alpha^3}} + \left( {1 + \frac{{a_\alpha^2}}{{10}}{\nabla^2}} \right)({\nabla_j}{u^{\prime}_i}({{\boldsymbol{x}}^\alpha }) + {\nabla_i}{u^{\prime}_j}({{\boldsymbol{x}}^\alpha }))} \right].\end{gather}

Using Faxen's equations, the velocity moments that include particle translational (${\boldsymbol{U}} - {{\boldsymbol{u}}^\infty }$) and angular velocity against the background flow (${\boldsymbol{\varOmega} } - {{\boldsymbol{\varOmega} }^\infty }$), as well as the rate of strain of the background flow $- {{\boldsymbol{E}}^\infty }$, are represented by the linear combination of force moments ${{\boldsymbol{F}}^H}$, ${{\boldsymbol{T}}^H}$ and ${{\boldsymbol{S}}^H}$. By representing the corresponding proportionalities in a matrix-tensor formulation, the far-field mobility matrix ${{\mathcal{M}}^\infty }$ is constructed. However, to account for lubrication interactions, all terms of the multipole expansion are necessary. Therefore, for particles within the lubrication cutoff radius, the exact resistivity is reflected in a pairwise additive manner. In this process, the far-field corrective term is excluded to prevent double-counting of far-field interactions, resulting in the contribution of the near-field resistance matrix ${{\mathcal{R}}_{nf}}$ (Jeffrey & Onishi Reference Jeffrey and Onishi1984; Jeffrey Reference Jeffrey1992). The near-field and far-field resistivity are combined to give the grand resistance matrix ${\mathcal{R}} = {({\mathcal{M}^\infty })^{ - 1}} + {{\mathcal{R}}_{nf}}$.

(2.15)\begin{equation}\left( {\begin{array}{@{}c@{}} {{{\boldsymbol{F}}^H}}\\ {{{\boldsymbol{T}}^H}}\\ {{{\boldsymbol{S}}^H}} \end{array}} \right) ={-} {\mathcal{R}}\boldsymbol{\cdot }\left( {\begin{array}{@{}c@{}} {{\boldsymbol{U}} - {{\boldsymbol{u}}^\infty }}\\ {{\boldsymbol{\varOmega} } - {{\boldsymbol{\varOmega} }^\infty }}\\ { - {{\boldsymbol{E}}^\infty }} \end{array}} \right).\end{equation}

The inverse problem is to obtain velocity moments from force moments, and the proportionality tensor ${\mathcal{M}}$ in this relation is called the grand mobility matrix, which has a relation ${\mathcal{M}}= {{\mathcal{R}}^{ - 1}}$ with the grand resistance matrix. Each component of the grand resistance matrix is important for subsequent developments and is given as shown in (2.16). Another characteristic of the proportionality tensor is their symmetry and positive definiteness, which arise from the dissipative nature of the system (Kim & Karrila Reference Kim and Karrila2005).

(2.16)\begin{equation}{\mathcal{R}} = \left( {\begin{array}{@{}ccc@{}} {{{\boldsymbol{R}}_{FU}}}&{{{\boldsymbol{R}}_{F\varOmega }}}&{{{\boldsymbol{R}}_{FE}}}\\ {{{\boldsymbol{R}}_{TU}}}&{{{\boldsymbol{R}}_{T\varOmega }}}&{{{\boldsymbol{R}}_{TE}}}\\ {{{\boldsymbol{R}}_{SU}}}&{{{\boldsymbol{R}}_{S\varOmega }}}&{{{\boldsymbol{R}}_{SE}}} \end{array}} \right).\end{equation}

Returning to (2.1) and (2.2), the non-hydrodynamic inter-particle force ${{\boldsymbol{F}}^P}$ is pairwise additive and is obtained by taking a derivative of the analytic form of the potential based on the inter-particle distance. Though various analytical potentials can be employed to study fouling dynamics, in our case study in § 3.2, we use a Morse potential due to its finite-sized derivative, including the potential well. In the definition in (2.17), ${r_{d,}}_{ij}$ is the inter-particle distance between particles i and j, equal to $|{{{\boldsymbol{x}}_i} - {{\boldsymbol{x}}_j}} |$, while ${a_i}$ and ${a_j}$ are the radius of particles i and j, respectively. However, the attraction range parameter $\kappa$ controls the magnitude of attraction and is suggested in accordance with the specific case study. The potential well depth ${V_0}$ is represented as $6{\rm \pi} \mu u_{avg}^\infty {a^2}/{V_0}$, signifying the dimensionless ratio of the characteristic hydrodynamic force $6{\rm \pi} \mu u_{avg}^\infty a$ to the attractive force ${V_0}/a$ in the specific simulation system. In this ratio, $u_{avg}^\infty$ is the surface average of the background fluid velocity at the inlet surface.

(2.17)\begin{gather}{\varphi _{Morse}}({r_{d,}}_{ij}) ={-} {V_0}\left[2\,{\textrm{e}^{ - \kappa \{ {r_{d,}}_{ij} - ({a_i} + {a_j})\} }} - {\textrm{e}^{ - 2\kappa \{ {r_{d,}}_{ij} - ({a_i} + {a_j})\} }}\right],\end{gather}

(2.18)\begin{gather}F_{ij}^P ={-} \frac{{\partial {\varphi _{Morse}}({r_{d,}}_{ij})}}{{\partial {r_{d,}}_{ij}}}.\end{gather}

Despite the singular lubrication forces, a small amount of particle overlap may occur due to the finite time step of the simulation (Ekiel-Jezewska, Tomasz & Szymczak Reference Ekiel-Jezewska, Tomasz and Szymczak2008). The same Morse potential is also used to prevent unphysical overlap at ${r_{d,ij}} < {a_i} + {a_j}$.

Finally, the Brownian random force and torque terms ${{\boldsymbol{F}}^B}$ and ${{\boldsymbol{T}}^B}$ originate from the thermal fluctuation of the fluid and are given by the following equations, with angle brackets representing an ensemble average over the entire suspension:

(2.19)\begin{gather}\langle {{\boldsymbol{F}}^B}(t)\rangle = \langle {{\boldsymbol{T}}^B}(t)\rangle = 0,\end{gather}

(2.20)\begin{gather}\langle {{\boldsymbol{F}}^B}(t){{\boldsymbol{F}}^B}(t^{\prime})\rangle = 2{k_B}T{{\boldsymbol{R}}_{FU}}\delta (t - t^{\prime}),\end{gather}

(2.21)\begin{gather}\langle {{\boldsymbol{T}}^B}(t){{\boldsymbol{T}}^B}(t^{\prime})\rangle = 2{k_B}T{{\boldsymbol{R}}_{T\varOmega }}\delta (t - t^{\prime}).\end{gather}

In this section, each term that constitutes the force balance in the Langevin equation has been introduced. However, there are additional considerations when describing particle dynamics inside a confined channel. For example, as wall particles are immobile, a constraint must be imposed to fix the hydrodynamic velocity moments to zero. Similarly, since there are no thermal fluctuations for the wall particles, (2.19)–(2.21) should also be appropriately modified to account for the Brownian motion of interior particles only while preserving multi-body hydrodynamic interactions. The derivation of the modified formulation, including these constraints, is carried out in the next section.

2.1.2. Stokesian dynamics with velocity constraint

In the framework of Stokesian-based methods, various attempts have been made to incorporate a channel's confined geometries in describing particle-laden flow. One of them is constructing a new Green's function in the presence of wall geometry using Blake's solution to obtain images of Stokeslets based on the distance H from the wall (Blake Reference Blake1971; Bossis, Meunier & Sherwood Reference Bossis, Meunier and Sherwood1991; Swan & Brady Reference Swan and Brady2007). This additional Green's function induced from wall geometry ${{\boldsymbol{J}}_{wall}}({\boldsymbol{x}},{\boldsymbol{y}};H)$ is superposed onto the fundamental function for particles in (2.22) to form the total Green's function ${\boldsymbol{J}}({\boldsymbol{x}},{\boldsymbol{y}};H) = {{\boldsymbol{J}}_{non\text{-}peri}}({\boldsymbol{x}},{\boldsymbol{y}}) + {{\boldsymbol{J}}_{wall}}({\boldsymbol{x}},{\boldsymbol{y}};H)$.

(2.22)\begin{equation}{{\boldsymbol{J}}_{non\text{-}peri}}({\boldsymbol{x}},{\boldsymbol{y}})\frac{1}{{8{\rm \pi} \mu }}\left( {\frac{{\boldsymbol{I}}}{r} + \frac{{{\boldsymbol{rr}}}}{{{r^3}}}} \right).\end{equation}

The above function is also called the Oseen tensor (Oseen Reference Oseen1927), while ${\boldsymbol{I}}$ denotes the identity tensor, ${\boldsymbol{r}} = {\boldsymbol{x}} - {\boldsymbol{y}}$ and ${r^2} = {\boldsymbol{r}}\boldsymbol{\cdot }{\boldsymbol{r}}$. For the fluid velocity at ${\boldsymbol{x}}$ induced by the periodic images at ${\boldsymbol{y}}$, the Green's function derived by Hasimoto (Reference Hasimoto1959) is used while ${\boldsymbol{k}}$ represents lattice vector on a reciprocal space:

(2.23)\begin{equation}{{\boldsymbol{J}}_{peri}}({\boldsymbol{x}},{\boldsymbol{y}}) = \frac{1}{{\mu V}}\sum\limits_{{\boldsymbol{k}} {\ne} 0} {{\textrm{e}^{\textrm{i}{\boldsymbol{k}}\,\boldsymbol{\cdot}\,({\boldsymbol{x}} - {\boldsymbol{y}})}}\frac{1}{{{k^2}}}({\boldsymbol{I}} - {\hat{\boldsymbol{k}}\hat{\boldsymbol{k}}})} .\end{equation}

When this total Green's function is integrated through a surface boundary integral that takes wall geometry into account, it allows for the evaluation of fluid disturbances at arbitrary points (Ladyzhenskaya Reference Ladyzhenskaya1963). For the extension in the case of multiple wall geometries, the same method of continuously superposing Green's functions has been used for studying the dynamics and rheology of particles placed between the parallel and perpendicular walls (Swan & Brady Reference Swan and Brady2003, Reference Swan and Brady2011; Staben, Zinchenko & Davis Reference Staben, Zinchenko and Davis2003). While it has the advantage of accurately reflecting the boundary effects where point forces are located, it also has certain limitations. The rapid increase in complexity with each superposed boundary makes it challenging to derive an explicit form of the Green function; an explicit form of the Green function does not already exist in the case of the flow between parallel slits according to Swan & Brady (Reference Swan and Brady2011). Expandability to accommodate changed boundary conditions is also limited, as the strategies to integrate force densities on each surface boundary have to be changed.

In this study, to circumvent these difficulties, the confined channel wall is simulated using a particle-wall scheme, as shown in figure 2. The primary advantage of this scheme is that it allows for the utilization of the same form of Green's function for both interior particles and wall particles, eliminating the need for calculating complex Green's functions. The placement of imaginary particles on the wall has been previously employed in a suspension balance model (Singh & Nott Reference Singh and Nott2000; Nott, Guazzelli & Pouliquen Reference Nott, Guazzelli and Pouliquen2011). In our study, we have developed an extension to arbitrary channel geometries by coupling a fluid dynamic solver to obtain background velocity moments while satisfying the microscopic force balances on the suspended particles and deriving macroscopic properties from the overall force balance on the channel. This approach develops into an efficient method for solving the interior dynamics of particles, leveraging the constraint properties on wall particles. Furthermore, by analysing the hydrodynamic force moments applied to the wall particles, it becomes possible to calculate the stress applied to a specific location on the channel wall, a capability challenging to achieve when using Oseen's superposition approach.

Figure 2. Schematic pictures of boundary treatment methods using wall particles. In panel (a), the excluded volume scheme is displayed, where the wall particles are entirely excluded from the fluid domain where fluid and particle velocity are evaluated. In panel (b), the included volume scheme is displayed, where the wall particles are arranged in a way that removes the excluded volume.

However, the accompanying limitations of the particle-wall scheme are examined. Since it models the boundary using wall particles, and point forces are located on the surfaces of spherical particles, it practically simulates a rough or bumpy wall. Moreover, as it considers particles as a patch on the wall approximating the wall resistivity, the computational cost increases rapidly as the system size becomes larger, in proportion to the cube of the total number of particles used. To mitigate this cubic dependence on the total number of particles, the effects from the wall particles are separated and treated as predetermined hydrodynamic force moments acting on the interior particles, originating from the properties of the prescribed velocity of wall particles. This treatment allows an efficient description of particle dynamics inside the channel, especially at the beginning of the simulation. However, as the volume fraction inside the channel increases in the late stage of the simulation, hydrodynamics within the channel becomes dominated by particle–particle interactions, and the influence of the wall on particle motion becomes less significant (Swan & Brady Reference Swan and Brady2011). The deviatoric effects of the particle-wall scheme on the transport properties and scaling performance of the constrained SD by the system size are discussed in § 3.1.

There are two ways of arranging particles to replace the wall boundary, which are the excluded volume scheme in figure 2(a), where the particles completely retreat from the target boundary, and the included volume scheme in figure 2(b), where the particles intrude over the target boundary. The excluded volume scheme has the advantage that macroscopic properties converge to those of pure fluid when there are no interior particles. However, the included scheme focuses on filling the place of the wall to the greatest extent. If the background flow is uniform, or in simulations where both force-rotational velocity and force-strain rate coupling do not exist (Nott & Brady Reference Nott and Brady1994), there is no practical difference between the two methods as they do not differ in the hydrodynamic forces experienced by interior and wall particles. However, in this study, the background velocity within the channel is calculated through the fluid dynamics solver, allowing for the application of any form of background flow field within the Stokes flow region. Due to its generality, force moments up to the dipole term, which corresponds to the same degree of accuracy as velocity moments, are considered. This increases the accuracy of the mobility relation governing HI. Thus, as depicted in figure 2(b), if particle centres are located on the target boundary, additional stresses are exerted on the channel wall even when there are no suspended particles due to the velocity gradient at the centre of wall particles. Therefore, in this study, we will proceed with the development according to the excluded volume scheme in figure 2(a) to avoid this situation.

If we expand the hydrodynamic forces on interior particles I and wall particles W through the resistance formulation in (2.15), each force can respectively be expressed as follows:

(2.24) \begin{align}{{\boldsymbol{F}}^I} &={-} {\boldsymbol{R}}_{FU}^{II}\boldsymbol{\cdot }({{\boldsymbol{u}}^I} - {{\boldsymbol{U}}^\infty }) + {\boldsymbol{R}}_{FU}^{IW}\boldsymbol{\cdot }{{\boldsymbol{U}}^\infty } - {\boldsymbol{R}}_{F\varOmega }^{II}\boldsymbol{\cdot }({{\boldsymbol{\omega} }^I} - {{\boldsymbol{\varOmega} }^\infty })\nonumber\\ &\quad + {\boldsymbol{R}}_{F\varOmega }^{IW}\boldsymbol{\cdot }{{\boldsymbol{\varOmega} }^\infty } + {\boldsymbol{R}}_{FE}^{II}\boldsymbol{\cdot }{{\boldsymbol{E}}^\infty } + {\boldsymbol{R}}_{FE}^{IW}\boldsymbol{\cdot }{{\boldsymbol{E}}^\infty },\end{align}

(2.25) \begin{align}{{\boldsymbol{F}}^W} &= {\boldsymbol{R}}_{FU}^{WW}\boldsymbol{\cdot }{{\boldsymbol{U}}^\infty } - {\boldsymbol{R}}_{FU}^{WI}\boldsymbol{\cdot }({{\boldsymbol{u}}^I} - {{\boldsymbol{U}}^\infty })\nonumber\\ &\quad + {\boldsymbol{R}}_{FU}^{WW}\boldsymbol{\cdot }{{\boldsymbol{\varOmega} }^\infty } - {\boldsymbol{R}}_{F\varOmega }^{WI}\boldsymbol{\cdot }({{\boldsymbol{\omega} }^I} - {{\boldsymbol{\varOmega} }^\infty }) + {\boldsymbol{R}}_{FE}^{WW}\boldsymbol{\cdot }{{\boldsymbol{E}}^\infty } + {\boldsymbol{R}}_{FE}^{WI}\boldsymbol{\cdot }{{\boldsymbol{E}}^\infty }.\end{align}

Here, ${{\boldsymbol{F}}^I}$ and ${{\boldsymbol{F}}^W}$ denote the hydrodynamic forces on the interior and wall particles exerted by the fluid, respectively. In subsequent developments, the additional superscript H will be omitted for the sake of simplicity since the grand resistance formulation itself is based on the common understanding that it represents relationships involving hydrodynamic force moments. Regarding the meaning of the proportionalities, ${\boldsymbol{R}}_{FU}^{II}$ corresponds to the resistance matrix, indicating the coupling between the force and translational velocity among the interior particles, while ${\boldsymbol{R}}_{FE}^{IW}$ represents the coupling between the force and strain rate among the interior and wall particles. The remaining components of the resistance matrix follow the same fashion of coupling notation. Here, ${{\boldsymbol{u}}^I}$ and ${{\boldsymbol{u}}^W}$ denote the translational velocities of interior and wall particles, respectively, while ${{\boldsymbol{\omega} }^I}$ and ${{\boldsymbol{\omega} }^W}$ represent the rotational velocities of interior and wall particles. Wall particles are assumed to be immobile, and velocity constraints are applied as ${{\boldsymbol{u}}^W} = {{\boldsymbol{\omega} }^W} = {\textbf{0}}$. Subsequently, as the resistance tensor is independent of the particles’ momentum and is pre-determined based on the configuration of the particles and the topology of the wall, the deterministic part of the hydrodynamic forcing term under the background velocity is substituted, as shown in (2.26) and (2.27), and the velocity moments of interior particles are reformulated, as in (2.28):

(2.26)\begin{gather}{\boldsymbol{A}} = ({\boldsymbol{R}}_{FU}^{II} + {\boldsymbol{R}}_{FU}^{IW})\boldsymbol{\cdot }{{\boldsymbol{U}}^\infty } + ({\boldsymbol{R}}_{F\varOmega }^{II} + {\boldsymbol{R}}_{F\varOmega }^{IW})\boldsymbol{\cdot }{{\boldsymbol{\varOmega} }^\infty } + ({\boldsymbol{R}}_{FE}^{II} + {\boldsymbol{R}}_{FE}^{IW})\boldsymbol{\cdot }{{\boldsymbol{E}}^\infty },\end{gather}

(2.27)\begin{gather}{\boldsymbol{B}} = ({\boldsymbol{R}}_{TU}^{II} + {\boldsymbol{R}}_{TU}^{IW})\boldsymbol{\cdot }{{\boldsymbol{U}}^\infty } + ({\boldsymbol{R}}_{T\varOmega }^{II} + {\boldsymbol{R}}_{T\varOmega }^{IW})\boldsymbol{\cdot }{{\boldsymbol{\varOmega} }^\infty } + ({\boldsymbol{R}}_{TE}^{II} + {\boldsymbol{R}}_{TE}^{IW})\boldsymbol{\cdot }{{\boldsymbol{E}}^\infty }{,_{}}\end{gather}

(2.28) \begin{gather}\left[ {\begin{array}{*{20}{@{}c@{}}} {{{\boldsymbol{F}}^P} + {{\boldsymbol{F}}^B} + {\boldsymbol{A}}}\\ {{{\boldsymbol{T}}^B} + {\boldsymbol{B}}} \end{array}} \right] = \left[ {\begin{array}{@{}cc@{}} {{\boldsymbol{R}}_{FU}^{II}}& {{\boldsymbol{R}}_{F\varOmega }^{II}}\\ {{\boldsymbol{R}}_{TU}^{II}}& {{\boldsymbol{R}}_{T\varOmega }^{II}} \end{array}} \right]\boldsymbol{\cdot }\left[ {\begin{array}{@{}c@{}} {{{\boldsymbol{u}}^I}}\\ {{{\boldsymbol{\omega} }^I}} \end{array}} \right],\end{gather}

(2.29) \begin{gather}\left[ {\begin{array}{*{20}{@{}c@{}}} {{{\boldsymbol{u}}^I}}\\ {{{\boldsymbol{\omega} }^I}} \end{array}} \right] = \left[ {\begin{array}{@{}cc@{}} {{{({\boldsymbol{R}}_{FU}^{II})}^{ - 1}}}&{{{({\boldsymbol{R}}_{F\varOmega }^{II})}^{ - 1}}}\\ {{{({\boldsymbol{R}}_{TU}^{II})}^{ - 1}}}&{{{({\boldsymbol{R}}_{T\varOmega }^{II})}^{ - 1}}} \end{array}} \right]\boldsymbol{\cdot }\left[ {\begin{array}{*{20}{@{}c@{}}} {{{\boldsymbol{F}}^P} + {{\boldsymbol{F}}^B} + {\boldsymbol{A}}}\\ {{{\boldsymbol{T}}^B} + {\boldsymbol{B}}} \end{array}} \right].\end{gather}

In writing (2.28), the governing equation in (2.1) substitutes the hydrodynamic force on the interior particles. Here, with regards to ${({\boldsymbol{R}}_{FU}^{II})^{ - 1}}$, it should be noted that it is not a result of direct inversion of the force-velocity coupling part of the resistance tensor, but rather the submatrix obtained after inverting the entire resistance matrix between interior particles in (2.28). The contribution from the inter-particle non-hydrodynamic interactions ${{\boldsymbol{F}}^P}$ is calculated in the pairwise additive fashion as in conventional SD, while the Brownian force and torque terms, represented as ${{\boldsymbol{F}}^B}$ and ${{\boldsymbol{T}}^B}$, can be simplified as follows:

(2.30)\begin{gather}\langle {{\boldsymbol{F}}^B}(t)\rangle = \langle {{\boldsymbol{T}}^B}(t)\rangle = 0,\end{gather}

(2.31)\begin{gather}\langle {{\boldsymbol{F}}^B}(t){{\boldsymbol{F}}^B}(t^{\prime})\rangle = 2{k_B}T{\boldsymbol{R}}_{FU}^{II}\delta (t - t^{\prime}),\end{gather}

(2.32)\begin{gather}\langle {{\boldsymbol{T}}^B}(t){{\boldsymbol{T}}^B}(t^{\prime})\rangle = 2{k_B}T{\boldsymbol{R}}_{T\varOmega }^{II}\delta (t - t^{\prime}).\end{gather}

While the Brownian random force has been explicitly added in (2.29), it alone does not fully account for the effects of thermal motion. In concentrated colloid systems, diffusivity varies depending on the configuration, and particle diffusion can be enhanced by random walks in structures with high diffusivity (Ermak & McCammon Reference Ermak and Mccammon1978). Therefore, a mean drift term with equal strength is added to counteract the erroneous particle depletion resulting from configurations that have high diffusivity (Brady & Bossis Reference Brady and Bossis1988).

(2.33) \begin{align}\left[ {\begin{array}{*{20}{@{}c@{}}} {{{\boldsymbol{x}}^I}({t_{n + 1}})}\\ {{{\boldsymbol{\theta} }^I}({t_{n + 1}})} \end{array}} \right] &= \left[ {\begin{array}{*{20}{@{}c@{}}} {{{\boldsymbol{x}}^I}({t_n})}\\ {{{\boldsymbol{\theta} }^I}({t_n})} \end{array}} \right] + ({t_{n + 1}} - {t_n})\nonumber\\ &\quad \times \left( \begin{array}{l} \left[ {\begin{array}{@{}cc@{}} {{{({\boldsymbol{R}}_{FU}^{II})}^{ - 1}}}&{{{({\boldsymbol{R}}_{F\varOmega }^{II})}^{ - 1}}}\\ {({\boldsymbol{R}}_{TU}^{II})^{ - 1}}&{{{({\boldsymbol{R}}_{T\varOmega }^{II})}^{ - 1}}} \end{array}} \right]\boldsymbol{\cdot }\left[ {\begin{array}{*{20}{@{}c@{}}} {{{\boldsymbol{F}}^P} + {{\boldsymbol{F}}^B} + {\boldsymbol{A}}}\\ {{{\boldsymbol{T}}^B} + {\boldsymbol{B}}} \end{array}} \right]\\ + \left[ {\begin{array}{*{20}{@{}c@{}}} {\boldsymbol{\nabla }\boldsymbol{\cdot }{{({\boldsymbol{R}}_{FU}^{II})}^{ - 1/2}}}\\ {\boldsymbol{\nabla }\boldsymbol{\cdot }{{({\boldsymbol{R}}_{T\varOmega }^{II})}^{ - 1/2}}} \end{array}} \right] \end{array} \right).\end{align}

Various algorithms have been proposed to account for the mean drift term simultaneously during time-stepping (Grassia, Hinch & Nitsche Reference Grassia, Hinch and Nitsche1995; Delong et al. Reference Delong, Usabiaga, Buscalioni, Griffith and Donev2014). Among these methods, the mid-point algorithm of Fixman (Reference Fixman1986) has been widely adopted and can be readily implemented within the current formalism. Employing this method, the evolution equations for the particle trajectories become

(2.34) \begin{align}{\left. {\left[ {\begin{array}{*{20}{@{}c@{}}} {{{\boldsymbol{x}}^I}}\\ {{{\boldsymbol{\theta} }^I}} \end{array}} \right]} \right|_{{t_{n + 1/2}}}} &= {\left. {\left[ {\begin{array}{*{20}{@{}c@{}}} {{{\boldsymbol{x}}^I}}\\ {{{\boldsymbol{\theta} }^I}} \end{array}} \right]} \right|_{{t_n}}} + \frac{{({t_{n + 1}} - {t_n})}}{2}{\left. {\left[ {\begin{array}{@{}cc@{}} {{{({\boldsymbol{R}}_{FU}^{II})}^{ - 1}}}&\ {{{({\boldsymbol{R}}_{F\varOmega }^{II})}^{ - 1}}}\\ {{{({\boldsymbol{R}}_{TU}^{II})}^{ - 1}}}&\ {{{({\boldsymbol{R}}_{T\varOmega }^{II})}^{ - 1}}} \end{array}} \right]} \right|_{{t_n}}}\nonumber\\ &\quad \boldsymbol{\cdot }\left[ {\begin{array}{*{20}{@{}c@{}}} {{{ {({{\boldsymbol{F}}^P} + {\boldsymbol{A}})} |}_{{t_n}}} + {{ {{{\boldsymbol{F}}^B}} |}_{{t_n}}}}\\ {{{ {\boldsymbol{B}} |}_{{t_n}}} + {{ {{{\boldsymbol{T}}^B}} |}_{{t_n}}}} \end{array}} \right],\end{align}

(2.35) \begin{align}{\left. {\left[ {\begin{array}{*{20}{@{}c@{}}} {{{\boldsymbol{x}}^I}}\\ {{{\boldsymbol{\theta} }^I}} \end{array}} \right]} \right|_{{t_{n + 1}}}} &= {\left. {\left[ {\begin{array}{*{20}{@{}c@{}}} {{{\boldsymbol{x}}^I}}\\ {{{\boldsymbol{\theta} }^I}} \end{array}} \right]} \right|_{{t_n}}} + ({t_{n + 1}} - {t_n}){\left. {\left[ {\begin{array}{@{}cc@{}} {{{({\boldsymbol{R}}_{FU}^{II})}^{ - 1}}}&\ {{{({\boldsymbol{R}}_{F\varOmega }^{II})}^{ - 1}}}\\ {{{({\boldsymbol{R}}_{TU}^{II})}^{ - 1}}}&\ {{{({\boldsymbol{R}}_{T\varOmega }^{II})}^{ - 1}}} \end{array}} \right]} \right|_{{t_{n + 1/2}}}}\nonumber\\ &\quad \cdot \left[ {\begin{array}{*{20}{@{}c@{}}} {{{ {({{\boldsymbol{F}}^P} + {\boldsymbol{A}})} |}_{{t_{n + 1/2}}}} + {{ {{{\boldsymbol{F}}^B}} |}_{{t_n}}}}\\ {{{ {\boldsymbol{B}} |}_{{t_{n + 1/2}}}} + {{ {{{\boldsymbol{T}}^B}} |}_{{t_n}}}} \end{array}} \right].\end{align}

First, by using (2.34), particle displacement at ${t_{n + 1/2}}$ is computed, and subsequently, employing the updated resistance matrix, particle displacement at ${t_{n + 1}}$ is determined according to (2.35). This step corresponds to the final phase of the algorithm, as illustrated in figure 1, which entails determining the new particle positions through time marching.

When drawing an analogy between the present approach and conventional SD, the grand resistance matrix presented in (2.15) is replaced with the resistance matrix between interior particles in (2.28). Furthermore, the Brownian motions described by (2.19)–(2.21) are replaced by those designed for interior particles in (2.30)–(2.32). Consequently, in constrained SD, the procedure for calculating the square root or inverse of the grand resistance matrix in conventional SD is replaced with the process of determining the square root or inverse of the resistance matrix, which accounts for the interactions among the interior particles. In essence, as the system size increases, primarily due to the greater number of wall particles, the computational overhead associated with direct matrix inversion is confined to that from interior particles. This structural advantage ensures that the number of independent equations is directly proportional to the number of interior particles, effectively minimizing the computational cost stemming from the number of wall particles.

However, the influence of wall particles on computational cost is not entirely negligible throughout the entire algorithm. This is because when obtaining the grand resistance matrix in (2.16), it is necessary to invert the far-field mobility matrix containing information about all the particles to account for long-ranged hydrodynamics induced by wall particles. As a result of this process, the dynamics of particles in accordance with ambient flow and thermal fluctuation reflect the effects from the channel walls. However, since hydrodynamic effects from far-field interactions change slowly over time, this process is performed only once per single diffusional time of a particle $6{\rm \pi} \mu {a^3}/{k_B}T$ (Nott & Brady Reference Nott and Brady1994), where ${k_B}$ is the Boltzmann constant and T is the system temperature. In our study, the time step used for simulations is 10⁻⁴ times a single diffusional time, so the inversion of the grand mobility matrix is performed once every 10⁴ steps, significantly reducing its computational proportion in the simulation.

However, as (2.29) takes a fully generalized form independent of the number of wall particles, this formulation can also be applied in conventional SD where wall particles are not present. Interestingly, Wang & Swan (Reference Wang and Swan2019) employed a similar form of the governing equation to investigate the heterogeneity effect of patched beads on colloid gelation. In this study, each of the sphere particles comprising the beads is compelled into rigid body motion through the action of constraint force. Although developed for different systems, this constraint-containing Stokesian dynamics formulation shares the structural advantages presented in (2.29) and could find applications in various systems.

Not only essential aspects of particle motion, but also the rheological properties of interior particles and the hydrodynamic stress ${{\boldsymbol{S}}^H}$ exerted on the channel wall can be efficiently calculated. First, obtaining the hydrodynamic force ${{\boldsymbol{F}}^W}$ on the channel wall is straightforward, achieved by performing the matrix-vector multiplication in (2.25) using the velocity vector of interior particles obtained during the time marching. Similarly, the hydrodynamic stress experienced by the wall and interior particles is determined using the velocity vector of internal particles and (2.15). More specifically, the particle-phase stress is divided into non-hydrodynamic and hydrodynamic contributions, as illustrated in (2.36) (Batchelor Reference Batchelor1970, Reference Batchelor1977). The first two terms, ideal osmotic pressure and inter-particle elastic stress, constitute the non-hydrodynamic stress, while the third term, representing the hydrodynamic component, is elaborated upon in (2.37) for interior particles (${\langle {\boldsymbol{\Sigma} }\rangle ^H} = {\langle {\boldsymbol{S}}\rangle ^I}$), and (2.38) for the wall stress simulated by wall particles (${\langle {\boldsymbol{\Sigma} }\rangle ^H} = {\langle {\boldsymbol{S}}\rangle ^W}$).

(2.36)\begin{gather}\langle {\boldsymbol{\Sigma} }\rangle ={-} nkT{\boldsymbol{I}} - n\langle {\boldsymbol{r}}{{\boldsymbol{F}}^P}\rangle + {\langle {\boldsymbol{\Sigma} }\rangle ^H},\end{gather}

(2.37) \begin{gather}{{\boldsymbol{S}}^I} ={-} {\boldsymbol{R}}_{SU}^{II}\boldsymbol{\cdot }({{\boldsymbol{u}}^I} - {{\boldsymbol{U}}^\infty }) + {\boldsymbol{R}}_{SU}^{IW}\boldsymbol{\cdot }{{\boldsymbol{U}}^\infty } - {\boldsymbol{R}}_{S\varOmega }^{II}\boldsymbol{\cdot }({{\boldsymbol{\omega} }^I} - {{\boldsymbol{\varOmega} }^\infty })\nonumber\\ \quad +\,{\boldsymbol{R}}_{S\varOmega }^{IW}\boldsymbol{\cdot }{{\boldsymbol{\varOmega} }^\infty } + {\boldsymbol{R}}_{SE}^{II}\boldsymbol{\cdot }{{\boldsymbol{E}}^\infty } + {\boldsymbol{R}}_{SE}^{IW}\boldsymbol{\cdot }{{\boldsymbol{E}}^\infty },\end{gather}

(2.38) \begin{align}{{\boldsymbol{S}}^W} &= {\boldsymbol{R}}_{SU}^{WW}\boldsymbol{\cdot }{{\boldsymbol{U}}^\infty } - {\boldsymbol{R}}_{SU}^{WI}\boldsymbol{\cdot }({{\boldsymbol{u}}^I} - {{\boldsymbol{U}}^\infty }) + {\boldsymbol{R}}_{S\varOmega }^{WW}\boldsymbol{\cdot }{{\boldsymbol{\varOmega} }^\infty } - {\boldsymbol{R}}_{S\varOmega }^{WI}\boldsymbol{\cdot }({{\boldsymbol{\omega} }^I} - {{\boldsymbol{\varOmega} }^\infty })\nonumber\\ &\quad + {\boldsymbol{R}}_{SE}^{WI}\boldsymbol{\cdot }{{\boldsymbol{E}}^\infty } + {\boldsymbol{R}}_{SE}^{WW}\boldsymbol{\cdot }{{\boldsymbol{E}}^\infty }.\end{align}

The hydrodynamic part of the stress is further divided into the stress from the background velocity field, the Brownian stress and the dissipative part of the inter-particle stress. If relative motion between particles is induced by inter-particle forces and they are hydrodynamically coupled, this motion induces stresslets (Chu & Zia Reference Chu and Zia2016). The individual contributions constituting the hydrodynamic part can be calculated by determining the slip velocity under the influence of the corresponding force using (2.29), followed by performing matrix-vector multiplications as described in (2.37) and (2.38).

One of the primary interests in fouling dynamics is the observation of how the deformation of clogs and resuspension kinetics is affected by the hydrodynamic stress resulting from the combined effects of particle–fluid interaction and particle–particle interaction. From a Lagrangian perspective, the inter-particle stress, which falls under the non-hydrodynamic contribution, has primarily been associated with the elastic energy stored within the given particle structure (Whittle & Dickinson Reference Whittle and Dickinson1997; Johnson, Landrum & Zia Reference Johnson, Landrum and Zia2018; Lee et al. Reference Lee, Jin, Kim, Myung and Ahn2021). However, a comprehensive understanding of the effects of hydrodynamics on agglomerate deformation has not yet been fully established. Therefore, the method developed in this section, which allows for an easy evaluation of the kinetics of individual particles influenced by the interaction of various opposing forces within the flow channel, can offer advantages in studying the impact of hydrodynamics on fouling dynamics.

While the constrained Stokesian dynamics formulation has indeed resulted in a significant enhancement of computational speed during the initial phases of the simulation, it fundamentally possesses a time complexity of $O({({N^I})^3})$ when the number of interior particles is ${N^I}$. Therefore, as fouling advances to later stages, the computational workload escalates rapidly. In the subsequent subsection, various acceleration schemes developed to facilitate efficient calculations for hydrodynamic coupling are presented, assuming constrained velocity on wall particles.

2.1.3. Acceleration scheme on the constraint Stokesian dynamics

The most time-consuming step in considering many-body hydrodynamics is matrix inversion, which becomes the rate-determining step in the SD paradigm. Therefore, to circumvent matrix inversion, a method of dividing hydrodynamic force into near-field lubrication and far-field contributions and obtaining the desired properties through an iterative scheme has been considered (Sierou & Brady Reference Sierou and Brady2001; Banchio & Brady Reference Banchio and Brady2003). According to this scheme, far-field force ${{\boldsymbol{F}}_{ff}}$, torque ${{\boldsymbol{T}}_{ff}}$ and stresslet ${{\boldsymbol{S}}_{ff}}$ are defined as follows:

(2.39) \begin{equation}\left[ {\begin{array}{*{20}{@{}c@{}}} {{{\boldsymbol{F}}_{ff}}}\\ {{{\boldsymbol{T}}_{ff}}}\\ {{{\boldsymbol{S}}_{ff}}} \end{array}} \right] ={-} {({{\mathcal{M}}^\infty })^{ - 1}}\boldsymbol{\cdot }\left[ {\begin{array}{*{20}{@{}c@{}}} {{\boldsymbol{u}} - {{\boldsymbol{U}}^\infty }}\\ {{\boldsymbol{\omega} } - {{\boldsymbol{\varOmega} }^\infty }}\\ { - {{\boldsymbol{E}}^\infty }} \end{array}} \right].\end{equation}

Therefore, the action of the grand resistance matrix is divided into contributions from near-field lubrication and far-field interactions:

(2.40) \begin{align}\left[ {\begin{array}{*{20}{@{}c@{}}} {\textbf{0}}\\ {\textbf{0}}\\ {{{\boldsymbol{S}}^H}} \end{array}} \right] &={-} \left[ {\left[ {\begin{array}{@{}ccc@{}} {{{\boldsymbol{R}}_{FU,nf}}}& \ {{{\boldsymbol{R}}_{F\varOmega ,nf}}}& \ {{{\boldsymbol{R}}_{FE,nf}}}\\ {{{\boldsymbol{R}}_{TU,nf}}}& \ {{{\boldsymbol{R}}_{T\varOmega ,nf}}}& \ {{{\boldsymbol{R}}_{TE,nf}}}\\ {{{\boldsymbol{R}}_{SU,nf}}}& \ {{{\boldsymbol{R}}_{S\varOmega ,nf}}}& \ {{{\boldsymbol{R}}_{SE,nf}}} \end{array}} \right] + {{({{\mathcal{M}}^\infty })}^{ - 1}}} \right]\boldsymbol{\cdot }\left[ {\begin{array}{*{20}{@{}c@{}}} {{\boldsymbol{u}} - {{\boldsymbol{U}}^\infty }}\\ {{\boldsymbol{\omega} } - {{\boldsymbol{\varOmega} }^\infty }}\\ { - {{\boldsymbol{E}}^\infty }} \end{array}} \right]\nonumber\\ &\quad + \left[ {\begin{array}{*{20}{@{}c@{}}} {{{\boldsymbol{F}}^P} + {{\boldsymbol{F}}^B}}\\ {{{\boldsymbol{T}}^B}}\\ { - {\boldsymbol{x}}{{\boldsymbol{F}}^P}} \end{array}} \right].\end{align}

In the above equation, the force balance for interior particles in the first and second rows is expanded as follows. The superscripts in the resistance matrix denote the submatrix of the near-field matrix corresponding to interior particle I or wall particle $W$.

(2.41) \begin{align} {\textbf{0}} & = {\boldsymbol{F}}_{ff}^I - {\boldsymbol{R}}_{FU,nf}^{II}\boldsymbol{\cdot }({{\boldsymbol{u}}^I} - {{\boldsymbol{U}}^\infty }) + {\boldsymbol{R}}_{FU,nf}^{IW}\boldsymbol{\cdot }{{\boldsymbol{U}}^\infty } - {\boldsymbol{R}}_{F\varOmega ,nf}^{II}\boldsymbol{\cdot }({{\boldsymbol{\omega} }^I} - {{\boldsymbol{\varOmega} }^\infty }) + {\boldsymbol{R}}_{F\varOmega ,nf}^{IW}\boldsymbol{\cdot }{{\boldsymbol{\varOmega} }^\infty }\nonumber\\ & \quad + {\boldsymbol{R}}_{FE,nf}^{II}\boldsymbol{\cdot }{{\boldsymbol{E}}^\infty } + {\boldsymbol{R}}_{FE,nf}^{IW}\boldsymbol{\cdot }{{\boldsymbol{E}}^\infty } + {{\boldsymbol{F}}^P} + {{\boldsymbol{F}}^B}, \end{align}

(2.42) \begin{align} {\textbf{0}} & = {\boldsymbol{T}}_{ff}^I - {\boldsymbol{R}}_{TU,nf}^{II}\boldsymbol{\cdot }({{\boldsymbol{u}}^I} - {{\boldsymbol{U}}^\infty }) + {\boldsymbol{R}}_{TU,nf}^{IW}\boldsymbol{\cdot }{{\boldsymbol{U}}^\infty } - {\boldsymbol{R}}_{T\varOmega ,nf}^{II}\boldsymbol{\cdot }({{\boldsymbol{\omega} }^I} - {{\boldsymbol{\varOmega} }^\infty }) + {\boldsymbol{R}}_{T\varOmega ,nf}^{IW}\boldsymbol{\cdot }{{\boldsymbol{\varOmega} }^\infty }\nonumber\\ & \quad + {\boldsymbol{R}}_{TE,nf}^{II}\boldsymbol{\cdot }{{\boldsymbol{E}}^\infty } + {\boldsymbol{R}}_{TE,nf}^{IW}\boldsymbol{\cdot }{{\boldsymbol{E}}^\infty } + {{\boldsymbol{T}}^B}. \end{align}

Similar to the previous subsection, since the resistance tensor is independent of the particle's momentum and determined solely by the particle's configuration and the topology of the wall, the predetermined hydrodynamic force from the wall can be substituted, as illustrated in (2.43) and (2.44). The remaining terms are then restructured into a matrix formulation, as presented in (2.45).

(2.43)\begin{gather}{\boldsymbol{A}^{\prime}} = {\boldsymbol{R}}_{FU,nf}^{IW}\boldsymbol{\cdot }{{\boldsymbol{U}}^\infty } + {\boldsymbol{R}}_{F\varOmega ,nf}^{IW}\boldsymbol{\cdot }{{\boldsymbol{\varOmega} }^\infty } + {\boldsymbol{R}}_{FE,nf}^{II}\boldsymbol{\cdot }{{\boldsymbol{E}}^\infty } + {\boldsymbol{R}}_{FE,nf}^{IW}\boldsymbol{\cdot }{{\boldsymbol{E}}^\infty },\end{gather}

(2.44)\begin{gather}{\boldsymbol{B}^{\prime}} = {\boldsymbol{R}}_{TU,nf}^{IW}\boldsymbol{\cdot }{{\boldsymbol{U}}^\infty } + {\boldsymbol{R}}_{T\varOmega ,nf}^{IW}\boldsymbol{\cdot }{{\boldsymbol{\varOmega} }^\infty } + {\boldsymbol{R}}_{TE,nf}^{II}\boldsymbol{\cdot }{{\boldsymbol{E}}^\infty } + {\boldsymbol{R}}_{TE,nf}^{IW}\boldsymbol{\cdot }{{\boldsymbol{E}}^\infty },\end{gather}

(2.45) \begin{gather}\left[ {\begin{array}{*{20}{@{}c@{}}} {{{\boldsymbol{u}}^I} - {{\boldsymbol{U}}^\infty }}\\ {{{\boldsymbol{\omega} }^I} - {{\boldsymbol{\varOmega} }^\infty }} \end{array}} \right] = \left[ {\begin{array}{@{}cc@{}} {{{({\boldsymbol{R}}_{FU,nf}^{II})}^{ - 1}}}& \ {{{({\boldsymbol{R}}_{F\varOmega ,nf}^{II})}^{ - 1}}}\\ {{{({\boldsymbol{R}}_{TU,nf}^{II})}^{ - 1}}}& \ {{{({\boldsymbol{R}}_{T\varOmega ,nf}^{II})}^{ - 1}}} \end{array}} \right]\left[ {\begin{array}{*{20}{@{}c@{}}} {{\boldsymbol{F}}_{ff}^I + {{\boldsymbol{F}}^P} + {{\boldsymbol{F}}^B} + {\boldsymbol{A}^{\prime}}}\\ {{\boldsymbol{T}}_{ff}^I + {{\boldsymbol{T}}^B} + {\boldsymbol{B}^{\prime}}} \end{array}} \right].\end{gather}

While (2.45) serves as an analogy to (2.29), the acceleration scheme includes the additional separation of far-field components into a force vector. This separation results in a sparse near-field matrix that, in most cases, becomes singular as contributions from the particles that have not entered the lubrication boundary layer become zero. To ensure positive definiteness and enhance convergence, the introduction of an isotropic tensor on the right-hand side leads to the transformation of the near-field resistance to ${\boldsymbol{R}}_{FU,nf}^{II} + \beta {\boldsymbol{I}} = \tilde{{\boldsymbol{R}}}_{FU,nf}^{II}$ and the far-field force on the interior particles to ${\boldsymbol{F}}_{ff}^I + \beta ({{\boldsymbol{u}}^I} - {{\boldsymbol{U}}^\infty }) = \tilde{{\boldsymbol{F}}}_{ff}^I$, where $\beta$ is a positive scalar. This treatment changes (2.45) into (2.46), making it a suitable form for practical implementation:

(2.46) \begin{equation}\left[ {\begin{array}{*{20}{@{}c@{}}} {\tilde{{\boldsymbol{F}}}_{ff}^I + {{\boldsymbol{F}}^P} + {{\boldsymbol{F}}^B} + {\boldsymbol{A}^{\prime}}}\\ {\tilde{{\boldsymbol{T}}}_{ff}^I + {{\boldsymbol{T}}^B} + {\boldsymbol{B}^{\prime}}} \end{array}} \right] = \left[ {\begin{array}{@{}cc@{}} {\tilde{{\boldsymbol{R}}}_{FU,nf}^{II}}& \ {\tilde{{\boldsymbol{R}}}_{F\varOmega ,nf}^{II}}\\ {\tilde{{\boldsymbol{R}}}_{TU,nf}^{II}}& \ {\tilde{{\boldsymbol{R}}}_{T\varOmega ,nf}^{II}} \end{array}} \right]\boldsymbol{\cdot }\left[ {\begin{array}{*{20}{@{}c@{}}} {{{\boldsymbol{u}}^I} - {{\boldsymbol{U}}^\infty }}\\ {{{\boldsymbol{\omega} }^I} - {{\boldsymbol{\varOmega} }^\infty }} \end{array}} \right].\end{equation}

As the far-field mobility remains undetermined, the translational and rotational velocities of the interior particles are expressed as functions of the far-field force and subsequently substituted into (2.40). Following several steps of rearrangement, it takes the form

(2.47) \begin{align} & \left( {{\mathcal{M}}^\infty }\boldsymbol{\cdot }\left( {{\boldsymbol{I}} - \beta \left[ {\begin{array}{@{}ccccc@{}} {{{(\tilde{{\boldsymbol{R}}}_{FU,nf}^{II})}^{ - 1}}}&{\textbf{0}}&{{{(\tilde{{\boldsymbol{R}}}_{F\varOmega ,nf}^{II})}^{ - 1}}}&{\textbf{0}}&{\textbf{0}}\\ {\textbf{0}}&{\textbf{0}}&{\textbf{0}}&{\textbf{0}}&{\textbf{0}}\\ {{{(\tilde{{\boldsymbol{R}}}_{TU,nf}^{II})}^{ - 1}}}&{\textbf{0}}&{{{({\tilde{{\boldsymbol{R}}}_{T\varOmega ,nf}^{II}} )}^{ - 1}}}&{\textbf{0}}&{\textbf{0}}\\ {\textbf{0}}&{\textbf{0}}&{\textbf{0}}&{\textbf{0}}&{\textbf{0}}\\ {\textbf{0}}&{\textbf{0}}&{\textbf{0}}&{\textbf{0}}&{\textbf{0}} \end{array}} \right]} \right)\right.\nonumber\\ &\left.\quad + \left[ {\begin{array}{@{}ccccc@{}} {{{(\tilde{{\boldsymbol{R}}}_{FU,nf}^{II})}^{ - 1}}}&{\textbf{0}}&{{{(\tilde{{\boldsymbol{R}}}_{F\varOmega ,nf}^{II})}^{ - 1}}}&{\textbf{0}}&{\textbf{0}}\\ {\textbf{0}}&{\textbf{0}}&{\textbf{0}}&{\textbf{0}}&{\textbf{0}}\\ {{{(\tilde{{\boldsymbol{R}}}_{TU,nf}^{II})}^{ - 1}}}&{\textbf{0}}&{{{(\tilde{{\boldsymbol{R}}}_{T\varOmega ,nf}^{II})}^{ - 1}}}&{\textbf{0}}&{\textbf{0}}\\ {\textbf{0}}&{\textbf{0}}&{\textbf{0}}&{\textbf{0}}&{\textbf{0}}\\ {\textbf{0}}&{\textbf{0}}&{\textbf{0}}&{\textbf{0}}&{\textbf{0}} \end{array}} \right] \right)\boldsymbol{\cdot }\left[ {\begin{array}{@{}c@{}} {\tilde{{\boldsymbol{F}}}_{ff}^I}\\ {\tilde{{\boldsymbol{F}}}_{ff}^W}\\ {\tilde{{\boldsymbol{T}}}_{ff}^I}\\ {\tilde{{\boldsymbol{T}}}_{ff}^W}\\ {{{\boldsymbol{S}}_{ff}}} \end{array}} \right]\nonumber\\ & \quad ={-} \left[ {\begin{array}{@{}c@{}} {{{(\tilde{{\boldsymbol{R}}}_{FU,nf}^{II})}^{ - 1}}\boldsymbol{\cdot }({{\boldsymbol{F}}^P} + {{\boldsymbol{F}}^B} + {\boldsymbol{A}^{\prime}}) + {{(\tilde{{\boldsymbol{R}}}_{F\varOmega ,nf}^{II})}^{ - 1}}\boldsymbol{\cdot }({{\boldsymbol{T}}^B} + {\boldsymbol{B}^{\prime}})}\\ { - {{\boldsymbol{U}}^\infty }}\\ {{{(\tilde{{\boldsymbol{R}}}_{TU,nf}^{II})}^{ - 1}}\boldsymbol{\cdot }({{\boldsymbol{F}}^P} + {{\boldsymbol{F}}^B} + {\boldsymbol{A}^{\prime}}) + {{(\tilde{{\boldsymbol{R}}}_{T\varOmega ,nf}^{II})}^{ - 1}}\boldsymbol{\cdot }({{\boldsymbol{T}}^B} + {\boldsymbol{B}^{\prime}})}\\ { - {{\boldsymbol{\varOmega} }^\infty }}\\ { - {{\boldsymbol{E}}^\infty }} \end{array}} \right]\nonumber\\ & \quad \quad + \beta {{\mathcal{M}}^\infty }\boldsymbol{\cdot }\left[ {\begin{array}{@{}c@{}} {{{(\tilde{{\boldsymbol{R}}}_{FU,nf}^{II})}^{ - 1}}\boldsymbol{\cdot }({{\boldsymbol{F}}^P} + {{\boldsymbol{F}}^B} + {\boldsymbol{A}^{\prime}}) + {{(\tilde{{\boldsymbol{R}}}_{F\varOmega ,nf}^{II})}^{ - 1}}\boldsymbol{\cdot }({{\boldsymbol{T}}^B} + {\boldsymbol{B}^{\prime}})}\\ {\textbf{0}}\\ {{{(\tilde{{\boldsymbol{R}}}_{TU,nf}^{II})}^{ - 1}}\boldsymbol{\cdot }({{\boldsymbol{F}}^P} + {{\boldsymbol{F}}^B} + {\boldsymbol{A}^{\prime}}) + {{(\tilde{{\boldsymbol{R}}}_{F\varOmega ,nf}^{II})}^{ - 1}}\boldsymbol{\cdot }({{\boldsymbol{T}}^B} + {\boldsymbol{B}^{\prime}})}\\ {\textbf{0}}\\ {\textbf{0}} \end{array}} \right]. \end{align}

The iterative form of the equation presented above includes velocity constraints for wall particles but remains in a fully generalized form. Consequently, it can be immediately reduced to the form used in the accelerated Stokesian dynamics scheme (ASD) for studying the systems such as dispersion between parallel plates or colloidal gels (Swan & Brady Reference Swan and Brady2011; Ouaknin, Su & Zia Reference Ouaknin, Su and Zia2021), even when the number of wall particles ${N^W}$ is zero. This adaptive nature does not require any changes to the governing equation. Alternatively, in systems where thermal motion is not significant, as indicated by ${{\boldsymbol{F}}^B} = {{\boldsymbol{T}}^B} = {\textbf{0}}$, it can be reduced into a form for obtaining deterministic force moments. In the presence of thermal motion, the relationships below replace (2.30)–(2.32), as thermal motion can also be divided into contributions from the near- and far-field parts. To avoid unnecessary repetition, only the force term is represented since torque is conventionally included in the force vector (Brady & Bossis Reference Brady and Bossis1988).

(2.48)\begin{gather}\langle {\boldsymbol{F}}_{nf}^B(t)\rangle = \langle {\boldsymbol{F}}_{ff}^B(t)\rangle = 0,\end{gather}

(2.49)\begin{gather}\langle {\boldsymbol{F}}_{nf}^B(t){\boldsymbol{F}}_{nf}^B(t^{\prime})\rangle = 2{k_B}T{\boldsymbol{R}}_{FU,nf}^{II}\delta (t - t^{\prime}),\end{gather}

(2.50)\begin{gather}\langle {\boldsymbol{F}}_{ff}^B(t){\boldsymbol{F}}_{ff}^B(t^{\prime})\rangle = 2{k_B}T({{\mathcal{M}}^\infty })_{FU}^{ - 1}\delta (t - t^{\prime}).\end{gather}

The near-field contribution, which encompasses the limited contributions from each particle pair $\gamma$ within the cutoff radius of $4a$, results in a sparse resistance matrix. Due to its pairwise additive nature, it can be factorized by Cholesky decomposition (Banchio & Brady Reference Banchio and Brady2003).

(2.51)\begin{equation}F_{nf,i}^B = \sum\limits_{(i,j) \in S(\gamma )} {(L_{ii}^\gamma \cdot {\varPsi _{nf,i}} + L_{ij}^\gamma \cdot {\varPsi _{nf,j}}).}\end{equation}

In this context, ${\boldsymbol{L}}$ represents the lower triangular matrix of ${\boldsymbol{R}}_{FU,nf}^{II}$, while ${{\boldsymbol{\varPsi} }_{nf}}$ is a stochastic vector with zero mean and unit variance. It is important to clarify that the superscript $II$ in the factorized resistance matrix does not imply that lubrication contributions are exclusively between the interior particles; rather, it denotes the submatrix corresponding to interior particles in the matrix. Indeed, while the near-field thermal motion of interior particles is influenced by contributions from the wall particles within the cutoff radius, the computational advantage lies in the fact that only the matrix components corresponding to interior particles undergo Cholesky decomposition.

In the next part, the far-field contribution unavoidably involves considering the far-field effects from all the particles, as mentioned in § 2.1.2, but this step needs to be executed only once per unit diffusion time. Additionally, an approximation of $({{\mathcal{M}}^\infty })_{FU}^{ - 1/2}{{\boldsymbol{\varPsi} }_{ff}}$ using the Lanczos process, a type of Krylov subspace method, can circumvent the direct inversion of the far-field mobility matrix. Since this process uses mobility components from all particles, this part aligns with conventional ASD, and the detailed procedure for the Lanczos process is omitted here but can be found in its pioneering works (Ando et al. Reference Ando, Chow, Saad and Skolnick2012; Chow & Saad Reference Chow and Saad2014; Ouaknin et al. Reference Ouaknin, Su and Zia2021).

Returning to the governing equation, (2.47) can be solved using an iterative method such as the generalized minimum residual method (Saad & Schultz Reference Saad and Schultz1986). By substituting the obtained far-field force back into (2.46), the particle velocity can be obtained. When expressed in a form that incorporates drift velocity, similar to those in (2.34) and (2.35), it takes the following form:

(2.52) \begin{align}{ {{{\boldsymbol{x}}^I}} |_{{t_{n + 1/2}}}} &= { {{{\boldsymbol{x}}^I}} |_{{t_n}}} + \frac{{({t_{n + 1}} - {t_n})}}{2}\left( \begin{array}{@{}l@{}} { {{{\boldsymbol{U}}^\infty }} |_{{t_n}}} + { {{{(\tilde{{\boldsymbol{R}}}_{FU,nf}^{II})}^{ - 1}}} |_{{t_n}}}\boldsymbol{\cdot }{ {(\tilde{{\boldsymbol{F}}}_{ff}^I + {{\boldsymbol{F}}^P} + {{\boldsymbol{F}}^B} + {\boldsymbol{A}^{\prime}})} |_{{t_n}}}\\ \quad +\, { {{{(\tilde{{\boldsymbol{R}}}_{F\varOmega ,nf}^{II})}^{ - 1}}} |_{{t_n}}}\boldsymbol{\cdot }{ {(\tilde{{\boldsymbol{T}}}_{ff}^I + {{\boldsymbol{T}}^B} + {\boldsymbol{B}^{\prime}})} |_{{t_n}}} \end{array} \right),\end{align}

(2.53) \begin{align}{ {{{\boldsymbol{x}}^I}} |_{{t_{n + 1}}}} &= { {{{\boldsymbol{x}}^I}} |_{{t_n}}} + ({t_{n + 1}} - {t_n})\nonumber\\ &\quad\times\, \left( \begin{array}{@{}l@{}} { {{{\boldsymbol{U}}^\infty }} |_{{t_{n + 1/2}}}} + { {{{(\tilde{{\boldsymbol{R}}}_{FU,nf}^{II})}^{ - 1}}} |_{{t_{n + 1/2}}}}\boldsymbol{\cdot }\{ { {(\tilde{{\boldsymbol{F}}}_{ff}^I + {{\boldsymbol{F}}^P} + {\boldsymbol{A}^{\prime}})} |_{{t_{n + 1/2}}}} + { {{{\boldsymbol{F}}^B}} |_{{t_n}}}\} \\ \quad +\, { {{{(\tilde{{\boldsymbol{R}}}_{F\varOmega ,nf}^{II})}^{ - 1}}} |_{{t_{n + 1/2}}}}\boldsymbol{\cdot }\{ { {(\tilde{{\boldsymbol{T}}}_{ff}^I + {\boldsymbol{B}^{\prime}})} |_{{t_{n + 1/2}}}} + { {{{\boldsymbol{T}}^B}} |_{{t_n}}}\} \end{array} \right).\end{align}

The theoretical time complexity of the algorithm outlined above corresponds to $O(N\,\textrm{log}\,N)$. Recently, Fiore & Swan (Reference Fiore and Swan2019) reported that the SD scheme can be accelerated up to $O(N)$ through a reformulation using the property of a saddle-point matrix. Appendix A demonstrates how to reformulate the constrained SD in a manner that exploits the properties of the saddle-point formulation. This is grounded in the common idea that both the saddle-point formulation and ASD share the concept of dividing hydrodynamic force moments into far-field and near-field interactions. This approach also enables the formation of the Schur complement required for the saddle-point formulation, even when there is a strong constraint on the velocity vector.

In summary, the algorithm developed in our study offers an advantage in terms of efficiency compared with Euler-based methods for systems with thousands of particles, as it does not require the process of remeshing after the application of new boundary conditions. While this study only implements the formulation in § 2.1.3 for the subsequent results, there are potential opportunities for future improvements for systems involving tens of thousands or more particles through GPU parallelization, as suggested by Fiore & Swan (Reference Fiore and Swan2019).

However, the process described in (2.52) and (2.53) corresponds to solving for the new particle positions through time marching, as depicted in the final step of figure 1. Before proceeding to the next section, in which we define the case system for fouling studies, we would like to provide a brief review of the coupling relationship between FEM and constraint SD. Obtaining the ambient velocity through FEM at the algorithm's outset and subsequently evaluating the velocity disturbance appears to be a one-way coupling, according to the form of the flowchart in figure 1. However, fluid and particles inherently interact through a two-way coupling. This is due to the continuous feedback loop that occurs between the fluid and particles through a series of processes: initially, the velocity disturbance is calculated at the centre of the particles, corresponding to the influence from the particles; subsequently, the changed positions of the particles are determined, corresponding to the influence from the fluid. Given the boundary conditions of the channel, FEM provides a steady-state solution of pure fluid as a reference in this feedback process. In the following section, the target geometry and boundary conditions are introduced, allowing for the acquisition of a unique FEM solution for the time-marching process.