Dynamics of Interacting Colloidal Particles Using the IIR Recursive Digital Filter Method

: This paper focuses on the numerical study of spherical particle sedimentation, taking into account hydrodynamic interactions. Infinite impulse response (IIR) digital filters, specially tailored to solve the sedimentation dynamics, were used in the present study to numerically solve the coupled ordinary differential equations with the time-dependent coefficients of the problem. Hydrodynamic interactions are modeled using the Rotne–Prager–Yamakawa (RPY) approximation, to which a correction is made to better account for short-range interactions. In order to validate both the proposed numerical resolution method and the RPY correction, this paper begins with the study of two interacting spherical particle sedimentation methods. Comparisons with previously published analytical or numerical results confirm the relevance of the present approach.


Introduction
The sedimentation of micro-particles (in the form of spheres, filaments, sheets, etc.) can represent a major obstacle to the effective use of colloidal suspensions in industrial environments.Colloidal particles experience gravitational, thermal, and hydrodynamic forces while moving in a fluid.The influence of the latter is a crucial aspect in various problems related to the mechanics and transport of suspensions.Important properties, such as the terminal sedimentation velocity of a suspension and the agglomeration rates of aerosols in the atmosphere, all depend on the relative motion of suspended particles [1].The same applies to the effectiveness of spray purification devices (aimed at removing suspended particles from gas streams).In addition, knowledge of the velocity distribution of microparticles within a suspension is essential for example in granulometry, where it can be used to determine the size distribution of the considered particles.
From a particle point of view, colloidal suspensions are usually studied using the Stokesian dynamics (SD) method introduced by Brady and Bossis [2].The SD method is based on solving the coupled Langevin equations describing the dynamics of N colloidal particles, which are written as follows when only translation is considered [2,3]: where m i and v i are the mass and the velocity vector of the ith particle, respectively; F h i is the hydrodynamic force exerted by the fluid on the ith particle; F nh i are the non-hydrodynamic forces (interparticle forces F int and/or external forces F e ) acting on the ith particle; and F b i is the Brownian force acting on the ith particle due to thermal agitation.The SD approach is implemented in well-known numerical codes such as Hoomd-Blues, Lammps, PyStokes, and many others.Langevin's equations are generally solved in these codes using the over-damped assumption, which consists of neglecting the inertial terms m i dv i dt in front of the viscous contribution to the hydrodynamic force F h i , leading to the following simplified equation of motion: In most unconstrained situations, the over-damped hypothesis yields remarkable results, provided that suspension behavior is only considered over long time scales, which is generally the case.However, when it comes to characterizing the local behavior of suspensions, i.e., at short times, the over-damped hypothesis may prove to be too imprecise.This is, for example, the case in optical two-point microrheology measurements [4,5].It is therefore within the framework of a global dynamics study that the present modeling is carried out, by solving Equation (1) over the whole time-scale, using first-order infinite impulse response digital filters (IIR).
Infinite impulse response (IIR) digital filters are usually limited to the field of linear time-invariant systems (LTISs), of which their dynamics are described by linear ordinary differential equations with constant coefficients, or equivalently by Laplace transfer functions expressed as polynomial fractions.In a previous paper, the use of first-order IIR filters was extended to the numerical resolution of Basset-type integro-differential equations [6], of which their Laplace transfer functions cannot be expressed as rational fractions.A similar approach was adopted here in order to extend the IIR-Filter method to the numerical resolution of system (1), which is composed of N coupled ordinary differential equations (ODEs), satisfied based on the particle velocities: v 1 , v 2 , • • • , and v N .These EDOs are linear as functions of velocities but non-linear as functions of particle positions r 1 , r 2 , • • • , and r N .It is therefore an original extension of the IIR-Filter method that is proposed in the present paper.
This paper first examines the sedimentation dynamics of two identical spherical microparticles to validate the method of solving the equations of the motion of colloidal particles in interactions using the IIR-Filter method.The results are compared with those obtained analytically or numerically, in other research works [7][8][9][10].Stimson and Jeffery [7] delivered, for the first time, an exact solution of the axisymmetric translational motion for two spherical particles moving at the same speed along their centerline.The authors obtained the solution in the form of an infinite series, which converges rapidly using bipolar coordinates.This is not the case when the particles are very close to each other.The calculated values of the forces are for identical spheres.O'Neill [11] then Goldman et al. [8] and finally Wakiya [12] also used this technique.The authors aimed to determine the resulting motion of two identical, homogeneous-free spheres, moving in an unlimited fluid, with a low Reynolds number.They determined the linear and angular velocities of the two spheres as a function of their relative separation and the orientation of the line of centers concerning the direction of gravity g.
Ganatos et al. [13] numerically solved the problem using the collocation technique.This technique was previously developed by other authors to process and study unbounded [14] and bounded [15] axisymmetric Stokes flows.

Resolution Method
Consider a finite number N of identical particles suspended in a base fluid with density ρ f and dynamic viscosity µ f (see Figure 1).In this study, the fluid analyzed was assumed to be incompressible, Newtonian, unbounded, and at rest at significant distances from the colloidal particles.The particles themselves are assumed to be rigid, spherical, and have a radius of a.They are supposed to be small enough to make inertial forces negligible at long-time dynamics in comparison to viscosity forces, but large enough to make the effects of thermal agitation (Brownian motion) negligible.Fluid properties (density and viscosity) play a crucial role in the sedimentation processes, affecting the behavior of particles as they settle through a fluid medium.These properties are related to each other by the Reynolds number Re.In the whole study, viscous flow regimes with Re less than 1 were considered and, in the case of water, at room conditions (density ρ f = 10 3 kg/m 3 and shear viscosity µ f = 1 mPa • s).Suspended particles were assumed to be silica, with density ρ PC = 2 × 10 3 kg/m 3 and radii of the order of a few micrometers, depending on the Reynolds numbers considered.

The Non-Hydrodynamic Forces
The dynamics of N > 1 colloidal particles, including hydrodynamic interactions, were examined.The effect of thermal agitation was neglected, and only the translational motion was considered.The initial positions were chosen in such a way that the particles did not come into contact (no contact forces) and were not subjected to any surface tension forces (unbounded fluid).This was, therefore, the problem of the fall of N solid spheres into a viscous fluid at rest.The external forces F e i (other than hydrodynamic) considered here were the weight P i and the Archimedes thrust Π i , which were exerted on particle (i): where ρ PC is the colloidal particle (PC) density.In addition, the non-hydrodynamic forces F nh i may also contain contact forces F c ij (of the Hertz-Voigt type for example), which was not taken into account in the present study.

The Hydrodynamic Interactions
Hydrodynamic interactions (HIs) are transmitted by the fluid between the particles that move through it.These interactions explain various observed effects.It is, for example, possible to understand the tendency of small objects to fall vertically or at different speeds if close enough during sedimentation.The hydrodynamic interactions F h i originate from the movement of the solvent (host liquid), which is described by the Stokes equation and is influenced locally and at a distance by the motion of the PCs.These disturbances are then transmitted by the solvent to the other PCs that compose the suspension, modifying its trajectories.These interactions are particularly significant in the case of concentrated suspensions, where the ϕ solid volume fraction is high.The disturbances caused by the fluid flow are assumed to be rapid enough (at the scale of colloidal particles) in comparison to the movement of the PCs.This means that flow induced disturbances propagation can, to a good approximation, be considered almost instantaneous [16].
Under these conditions, the flow of the fluid at a given time t is only a function of the positions and velocities of the N PCs that compose the suspension.Thus, the hydrodynamic force F h i exerted by the solvent on the PC (i) is a function of the N position vectors r j and the N velocity vectors v j of the PCs: The linearity of the Stokes equation allows us to establish the relationships between the translational velocities of particles and the forces acting on them.The hydrodynamic interactions between colloidal particles are expressed using the mobility tensor and the associated mobility matrix (M d,d ).The components M d,d ij can be computed using a variety of methods [17,18].The following relation defines the coefficients of the mobility matrix in the case of pure translation and a uniform imposed bulk flow v ∞ : The mobilitymatrix (M d,d ) is a symmetrical square matrix d × d, with d being the dimensionality of the suspension (i.e., d = 3 here).There is no general analytic expression of the mobility matrix (M d,d ) for N > 2. However, the expression of the matrix (M d,d ) can be simplified when only binary interactions are considered.Equation ( 5) is then written in the following simplified form, which will be systematically used in this work: where the flow velocity v ∞ of the fluid in the absence of particles is assumed to be zero throughout this study (fluid initially at rest).In assuming that the suspension has spherical symmetry, the M d,d ij r i , r j coefficients can be written as M d,d ij r ij with r ij = ∥r i − r j ∥.
The hydrodynamic forces F h j can then be determined by reversing the relation (6), i.e., by calculating the matrix (M d,d ) −1 .

Expressions of the Mobility Matrix Coefficients
The main purpose of this section is to define the approximate expressions of the mobility matrix coefficients.Two common approximations are presented: the quasi-point approximation (AQP) and the Rotne-Prager-Yamakawa approximation (RPY).

The Quasi-Point Approximation (QPA)
In the quasi-point approximation, which is the most rudimentary approach, each sphere is assimilated to a source point, also called a Stokeslet.The so-called reflection method makes it possible to iteratively take into account hydrodynamic interactions.To analyze the hydrodynamic interactions between two spheres of radii a 1 and a 2 , we can limit ourselves to the first reflection between two PCs, (1) and ( 2), and apply Faxen's laws [19].This helps us write the mobility matrix in the form of (7), which takes into account the Oseen-Burgers tensor S. For two spheres in three dimensions, this equation can be expressed as follows: where I 3 is the identity matrix (for d = 3) and r 12 = ∥x 1 − x 2 ∥.We recall here that r 12 • r 12 is not a scalar but a tensor, with components that are given in Cartesian coordinates by As an illustration, the translational mobility matrix is expressed here for two identical spheres and d = 3: where r 12 • r 12 is given by (8).Even if an analytic inversion of this 2d × 2d matrix is possible, it leads to humanly unusable expressions.Numerical inversion is therefore preferable for determining hydrodynamic forces F h i .

Rotne-Prager-Yamakawa (RPY) Approximation
This approximation is more accurate than the quasi-point approximation developed earlier.In the present study, the expression of the Rotne-Prager-Yamakawa tensor (RPY) proposed by Wajnryb et al. [20,21] is used for two different spheres of radii a 1 and a 2 in the absence of penetration (r 12 > a 1 + a 2 ).The matrix representation of mobility for two different spheres is given in this case by Let us recall the main properties of this approximation: • At the limit where r 12 ≫ sup(a 1 , a 2 ), the Rotne-Prager-Yamakawa mobility matrix tends to be the mobility matrix of the point approximation.This turns out to be legitimate and corresponds to the so-called far-field hypothesis.

•
When both spheres have the same radius, the following expression for the generalized mobility matrix is obtained: The literatureevaluates the coefficients of the mobility matrix for different configurations.Kim et al. [22] performed the computation of coefficients for two identical spheres, while Jeffrey and Onishi [23] considered two different spheres, using the Lamb expression [24].Perkins and Jones [25] were interested in the case of a single particle close to a flat wall, while Beenakker et al. [26] examined the case of several particles in the presence of a wall.Finally, Van Saarloos and Mazur [27] studied the case of several particles suspended in an infinite fluid.

Equations of Motion
The centers of inertia C j of spherical particles are assumed to be initially located in the z = 0 plane.The initial position of C j is supposed to be OC i (0) = (x 0 i , y 0 i , 0).In applying Newton's second law to each particle, the resulting system of coupled ordinary differential equations is as follows: where (M d,d ) −1 is the inverse of the mobility matrix defined above.
The present study sought to solve system (12) with zero initial conditions on the velocities: v i (0) = 0, ∀i = 1, 2, . . ., N. System ( 12) is rewritten as follows: where is a d × d matrix, with coefficients given by ( 7) or (10), depending on the considered approximation.To our knowledge, there is no systematic exact analytic approach (such as a Laplace transformation or variation of the constant) for solving the system analytically (13).Thus, the resolution exposed in the present study is based on a numerical approach that implements recursive IIR digital filters (the IIR-Filter method).This type of numerical approach has been successfully used to solve linear parabolic partial differential equations from heat conduction problems [28,29].It has also been used to solve linear convolutional integro-differential equations in fluid mechanics [6].In the former case, the time complexity of the IIR-Filter method is O(n), which is lower than numerical methods such as a piecewise linear approximation, which features a time complexity of O(n 2 ).
The present study extends the scope of the IIR-Filter method to the case of systems of linearly coupled ordinary differential equations with time-dependent coefficients.

Let us first write NI
as the non-time-invariant term in the system of Equation ( 13) for particle (i).The mobility matrix depends on the positions of the particles that are in motion.Therefore, these positions change continuously over time, and the components of the mobility tensor also vary over time.This represents a considerable difficulty for the numerical solution of the equations of motion.The application of the Laplace transform L to system (13) leads to the following system of algebraic equations: where v i (s), F e i (s), and NI i (s) are the Laplace transforms of v i , F e i , and NI i , respectively.Since there is no known analytic expression of velocities v i that verifies (13), the Laplace transform NI i (s) of the non-invariant term NI i cannot be computed analytically.Any direct analytic or numerical inversion of v i (s) (i.e., using an analytic expression of v i (s)), thus becomes impossible.In this case, the numerical resolution must use methods that do not require knowledge of NI i (s).The IIR-Filter method allows for this.This is the main interest of the digital approach developed in this work.We recall that, for simplicity, this study is based on the assumption that v i (0) = 0.In the framework of the IIR-Filter method, the excitatory force F e i (t) must be sampled into a regular set of discrete values in time F e i (nT e ) , where T e = ∆t is the sampling period, in order to calculate the velocity v i (nT e ).
The application to Equation ( 14) of the bilinear transformation of the type s → 2 T e z−1 z+1 [30], allows the transition from Laplace space to discrete time space, with z being the complex number introduced in the definition of the z-transformation (Z) of a discrete time signal x(nT e ): 1+z −1 in Equation ( 14) yields resulting in the following algebraic system, with the intervention of the z-transforms of the different functions of the problem: In noting that the non-invariant term system ( 15) is rewritten in the following form: From these algebraic z-equations, N coupled recurrence relations are deduced to numerically solve the problem of the dynamics of N spherical colloidal particles at time t n = nT e , taking into account hydrodynamic interactions: where Key result (18) is applied, in the following paragraphs, to the calculation of the dynamics of systems composed of two, and then three, colloidal particles, when hydrodynamic interactions are not negligible.The flowchart shown in Figure A1 summarizes the process of numerically solving the coupled recursion Equations (18) of the Stokesian dynamics of a colloidal suspension using the IIR-Filter method.

Full Dynamics versus Long-Time Dynamics
To test the relevance of numerically solving Equation (1) rather than (2), we consider here an elementary system in which two identical particles drive through hydrodynamic interactions, a third particle initially at rest (see Figure 2).The equations of motion are, in this case, neglecting the influences of gravity and thermal agitation: where F e 2 and F e 3 are the external forces that must be exerted on C 2 and C 3 , respectively, in order for them to move with uniform rectilinear translational motion.The dynamics, Equation (19), are solved using the IIR-Filter method (full dynamics) and also using the usual over-damped assumption (long-time dynamics).
We considered the case of three identical spherical particles (radius a = 10 µm) made of silica (density ρ PC = 2 × 10 3 kg/m 3 ), immersed in water at rest, at room conditions (dynamic viscosity µ f = 1 mPa • s and density ρ f = 10 3 kg/m 3 ).C 2 and C 3 are assumed to move at constant velocity U = Ue x , with U = 10 cm/s, and it is assumed that, initially, ∆y = 100 µm and ∆x = 10a.Equation ( 19) has been solved with a time step T e = 2 µs.
Figure 3a compares the short-time dynamics of C 1 , obtained without over-damping (solid line) and with over-damping (dotted line).It is clear from this figure that short-time dynamics are not well described in the present situation by the over-damped approximation.It is also clear from Figure A2 that the two approaches converge at the long-time scale of observation.
In the same way, Figure 3b shows that the external force F e 2 required to ensure uniform rectilinear motion of C 2 (or C 3 ) is not precisely calculated for short time scales using the over-damped method.
It can be concluded from this elementary study that measurements based on longtime results can easily be obtained by the over-damped method and possibly also by the IIR-Filter method (which is undoubtedly more complex to implement).If, on the other hand, it is necessary to evaluate the dynamics of colloidal particles at short time scales, it is imperative in this case to resort to solving the equations of motion without making the over-damped approximation.In this case, the IIR-Filter method is an interesting approach, as illustrated in the following paragraphs.

Sedimentation of Two Identical Colloidal Particles
The configuration of two identical spherical particles is a simple and effective way to test the accuracy and convergence of the numerical resolution method presented in the present study (result (18)).The numerical results obtained using our approach are first compared with approximate or numerical analytical solutions available in the literature.The "ideal" configuration, where two identical spherical particles sedimenting at small Reynolds numbers in a viscous fluid at rest and unbounded, has numerous applications in colloid chemistry, rheology, meteorology, and the sedimentation of dilute suspensions.

Presentation
Let us consider the sedimentation problem of two identical spherical particles of radius a, density ρ PC , and suspended without initial velocity (v 1 (0) = 0 and v 2 (0) = 0) in a solvent initially at rest (v ∞ = 0), with viscosity µ f and density ρ f (see Figure 4).The centers of inertia C 1 and C 2 of the two spherical particles are placed in the z = 0 plane.The initial positions of C 1 and C 2 are OC 1 (0) = (x 0 1 , y 0 1 , 0) and OC 2 (0) = (x 0 2 , y 0 2 , 0).We need to solve the system of coupled differentials, Equation (12), which becomes, in the case of N = 2, identical to the following particles: with the initial conditions v 1 (0) = 0 and v 2 (0) = 0.By explaining the various vectors, we obtain here

The Results
Equation ( 21), describing the fall dynamics of two identical spheres, is solved using recurrence relations (18).One of our objectives was to choose among the two approximations mentioned above the one that delivers the results closest to physical reality.The two spheres are initially placed at the same altitude x 0 1 = x 0 2 , with the centers being 2a 1 + dy apart.When the two particles settle, they remain at the same distance from each other by moving at the same speed, denoted v x .This is a consequence of the reversibility property of flows at low Reynolds numbers.This common velocity of the two particles is higher than the velocity that each of them would have taken separately.The smaller the initial distance between the particles, the greater the speed of the doublet.
Figure 5a shows the temporal variations in the v x component of the two particles.
Figure 5b compares the limit velocity of the fall of a single sphere x obtained numerically in the presence of a second identical sphere, for the initial conditions indicated above.The verticality of the fall is checked (v 1y = v 2y = 0 m/s) as well as its location in the z = 0 plane.x of the two spheres relative to the limit fall velocity V 0 s of one sphere, as a function of 2a/r 12 .
Figure 5b legitimately reveals that both the RPY and AQP approximations tend toward a common limit when the distance r 12 between the centers of the two spheres becomes large in front of their diameter d = 2a.The physical reason for this common limit lies in the fact that at a high distance from a small source (in this case, a sphere), the latter can be considered as point-like.The RPY approximation, on the other hand, yields higher-limit velocities than those of the quasi-point approximation when the two spheres are close.This discrepancy is due to the fact that the quasi-punctual approximation (AQP) is less valid the closer the spheres are.
It can also be seen that the fall velocity of each of the two spheres is higher than the fall velocity of a single sphere.The closer the spheres are initially, the higher the velocity.Hydrodynamic interactions thus accelerate the fall velocity at low Reynolds number in an incompressible, resting, unbounded Newtonian viscous fluid.This result is important for the modeling of the drying of complex suspensions by Stokesian dynamics.
The influence of a vertical shift between the two spheres on their relative motion was also analyzed.Positioning sphere (1) above sphere (2), during the movement of (1), results in the displacement of the fluid by "pushing" fluid in front of the sphere and the deflection to the right of sphere (2) (Figure 6a).The deflections are reversed when sphere (1) is positioned below sphere (2), as shown in Figure 6b.
The influence of a vertical shift between the two spheres on their relative motion was also analyzed.For example, if sphere ( 1) is initially positioned above sphere (2), then during its movement, (1) displaces the fluid by pushing the fluid in front of it, and consequently, sphere ( 2) is deflected to the right, as shown in Figure 6a.The deflections are reversed when sphere ( 1) is positioned below sphere (2), as shown in Figure 6b.The accuracy and convergence of the IIR-Filter method are now examined by comparing the numerical results obtained with those resulting from other analytical and/or numerical techniques available in the literature [7,9,13].In [7], the authors determined the motion created in a viscous fluid at rest at infinity using two solid spheres moving with small, constant, and equal velocities, parallel to their center line.Their solution is based on the determination of the Stokes current function for the movement of the fluid.The forces required to maintain the motion of the spheres were then calculated.In the case of identical spheres, the forces are equal, and their modulus F is written as follows: where r = a/ sinh α.The ratio r/2a of the distance between the centers of the two spheres to their diameter is written here as cosh α.According to [9], the coefficient λ ∥ can be written in the following form, which we denote as λ B ∥ : Thus, λ ∥ represents the ratio between the force required to maintain the motion of one sphere in the presence of the other and the force necessary to maintain its motion with the same velocity if the other sphere was at an infinite distance.Figure 7 shows the evolution of the drag correction factor λ ∥ of the two spheres when falling parallel to the line connecting their centers.
In the case of the motion of the two spheres being perpendicular to the line of their centers, Goldman et al. [8] found the following empirical formula: where ξ = a/(y 2 − y 1 ).
The results obtained agree with those of Brenner [9] for large separation distances in front of the diameter d = 2a of the spheres.In this case, the relative error is less than 10 −2 %.On the other hand, the relative error increases and becomes non-negligible as the separation distance between the two particles decreases.This behavior is explained by the influence of near-field interactions (lubrication effect), which are not properly taken into account in the mobility matrix.
Figure 8 shows the evolution of the drag correction factor λ ⊥ for two spheres falling perpendicular to the line connecting their centers.The results obtained using the chosen resolution method, with the use of the RPY approximation, are consistent with those of Goldman et al. [8] (refer to the insert in Figure 6 showing the relative error in %, defined as for the sedimentation of two identical spherical particles in a perpendicular direction at the line connecting their centers.
The relative error calculated in the case of the parallel fall is greater than that calculated in the perpendicular situation.The physical origin of this difference is due to the fact that in the first case, the dynamics of the particles downstream are strongly influenced by the disruption of the flow caused by the sphere upstream, just as when a ship moves in the wake of another rather than alongside it.
As shown both in Figures 7 and 8, the numerical results obtained by IIR filters are consistent with those obtained analytically for relatively large separation distances, as noted above.In the case of small separation distances, where deviations from analytical results are more pronounced, we found that the results obtained with the IIR-Filter method are very close to those obtained with the Euler method for the two approximations considered here (QPA and RPY).The observed discrepancies are therefore due more to the limitations of the QPA and RPY approximations than to the calculation methods.We therefore proposed in this work to correct the RPY mobility matrix in order to obtain better results at low separation distances.Error (in %)  Error (in %) Figure 8. Evolution of the drag correction factor λ ⊥ ; the two spheres fall perpendicular to the line connecting their centers.

Correction to the RPY Mobility Matrix
In order to correct the deviations observed at small separation distances, we propose to introduce, in the RPY mobility matrix, a corrective term of an order greater than 2 into a/r 12 .The corrected matrix RPY (M) c is thus written as  (25) where α(r 12 ) is a corrective function.The aim of the present correction is to minimize the discrepancy between our numerical results obtained using the RPY approximation and the analytical results of both H. Brenner [9], in the case of the motion of the two particles parallel to the axis connecting their centers, and of Goldman et al. [8], in the case of perpendicular motion.The expression of the function α(r 12 ) should be determined as a function of the separation distance r 12 .This function reduces the discrepancy between the results of the present study and the analytical solutions [8,9] for a given separation distance of r 12 .Figure 9 shows the α values leading to a minimum error for a given value of ζ = 2a/r 12 .
An interpolation of the results shown in Figure 9 led to the following expression of α as a function of the parameter ζ = 2a/r 12 , with a regression coefficient R 2 = 0.998: The numerical IIR-Filter solution of the two-particle sedimentation problem with hydrodynamic interactions was repeated using the corrected mobility matrix (25), with the α expression (26).As can be seen from Figure 10a,b, the introduced correction considerably improves the correspondence between numerical and analytical results, which confirms its relevance and allows the lubrication effect to be integrated very simply into the RPY mobility matrix.Error (in %)

Sedimentation of Three Particles
The study of the sedimentation of a triplet of particles (see Figure 11) is discussed here.The three identical spherical particles can be arranged evenly or irregularly.The literature examines these different configurations [13,31,32].
In the first configuration, the three spheres are positioned in such a way as to be equidistant along the horizontal axis, with a distance equal to 6a between each particle.The particle, initially at the center of this configuration, moves during sedimentation, at a higher speed than the two spheres positioned at the ends.It drags the two extreme spheres with it, causing them to come together.
Figure 12 illustrates the trajectories of the three particles obtained using the IIRbased resolution method with the corrected RPY mobility matrix.The velocities of the two particles (C 1 ) and (C 3 ) increase as they approach and exceed the velocity of the sphere (C 2 ).This pair of particles then approaches the sphere (C 2 ) until it reaches a critical distance, as described by Ganatos [13].The numerical results for this configuration are compared with those obtained by Ganatos et al. [13] who solved the steady-state creeping-motion governing equations of the flow (Stokes equation and mass conservation) for three sedimenting spheres.Figure 13 illustrates the variation in the horizontal distance y 21 between one of the extreme spheres (1 or 3) and the central sphere (2) as a function of dimensionless time t * = V 0 s t/a. Figure 14 shows the temporal variation in the vertical distance x 21 between one of the spheres at the ends and the central sphere.The results obtained from previous research conducted by Ganatos et al. [13] are included in both cases.The comparison between the numerical results of this study and the results of Ganatos et al. [13] reveals a global agreement in the analytical and numerical evolutions of the distances between the particles.For t * ⩽ 100, the curves overlap almost perfectly.However, discrepancies t * > 100, when the correction factor in the mobility matrix is not integrated.These differences are greatly reduced when the corrective function is integrated into the same matrix, thus confirming the effectiveness of the lubrication correction we introduced.
Figure 15 shows the temporal evolution of the vertical velocities of the three particles obtained using the IIR-Filter method, using the RPY c −corrected mobility matrix.The velocity of the central sphere remains higher than those of the two extreme spheres, up to t * = 106.From this point on, the particle doublet gradually approaches the central sphere until it reaches the critical distance.The physical explanation for this strange behavior lies in the hydrodynamic interactions between the three spheres.Initially, the C 2 sphere is pushed forward by C 1 and C 3 through the hydrodynamic interactions (incompressible fluid).As a result, C 2 falls faster than C 1 and C 3 .Before t * , spheres C 1 and C 3 are disturbed by the wake of sphere C 2 , and approach each other, forming a doublet, and its fall velocity gradually increases (see Figures 5 and 15) until it equals that of C 2 at t = t * .Finally, for t > t * , the doublet passes C 1 and continues to accelerate.According to the results of Ganatos et al. [13], the value of t * at which the vertical velocities of the three particles become equal is 107.8.This represents a deviation of 1.6% from the results obtained with the IIR-Filter method based on the corrected RPY mobility matrix, thus confirming the relevance of both the lubrication correction we of the IIR-Filter method.

Conclusions
The sedimentation dynamics ordinary differential equations (ODEs) of a number N > 1 of colloidal particles (PCs), coupled by hydrodynamic interactions and suspended in an incompressible, unbounded Newtonian fluid, were numerically solved in this study for the first time using infinite impulse response recursive filters (the IIR-Filter method).The N = 2 and N = 3 cases were analyzed in detail, and many analytical and numerical results already established with other approaches were accurately found numerically using the IIR-Filter method.
An effective correction to the mobility matrix was also proposed in this work, within the framework of the Rotne-Prager-Yamakawa (RPY) approximation.This correction takes better account of hydrodynamic interactions between colloidal particles in the near-field (lubrication phenomenon at low inter-particle separations).This correction significantly improves the accuracy of the numerical results provided by the RPY approximation in the case of two and three particles.The theoretical results, obtained by directly solving the partial differential equations (PDEs) of the Stokesian dynamics of the fluid flowing around the colloidal particles moving in it, can then be recovered with great precision by solving the PCs's ODEs with the mobility matrix approach.
The range of applications for the present approach is vast, and could extend, for example, to the study of microrheology or the drying of complex suspensions.

Figure 1 .
Figure 1.N spherical particles suspended in a viscous fluid initially at rest.

Figure 2 .
Figure 2. Two identical particles (C 2 and C 3 ), moving according to a uniform rectilinear translation at the same velocity U 2 = U 3 = U = Ue x , are driven by hydrodynamic interactions to a third particle (C 1 ), initially at rest.

Figure 3 .
Figure 3. Comparisons of C 1 short-time dynamics, obtained without making the over-damped assumption (solid line) and making the over-damped assumption (dashed line).(a) Changes in the U 1,x component of C 1 's velocity as a function of time.(b) Changes in the external force F e 2 exerted on C 2 , as a function of time.

Figure 4 .
Figure 4. Sedimentation of two identical spheres in a stationary fluid.

Figure 5 .
Figure 5. Fall velocities of two identical spheres taking into account hydrodynamic interactions.(a) Changes in the v x component of the two spheres as a function of reduced time t * .(b) Changes in the fall velocity v ℓ,HI

Figure 6 .
Figure 6.Deviations of the fall trajectories of two spheres in hydrodynamic interaction.(a) Trajectories of the two sphere when the blue one (C 1 ) is initially positioned above the red one (C 2 ).(b) Reverse situation of (a).

Figure 7 .
Figure 7. Evolution of the drag correction factor λ ∥ of the two spheres falling parallel to the line connecting their centers.

Figure 10 .
Figure10.Evolution of the drag correction factor λ: (a) the two spheres' sediment parallel to the line connecting their centers; (b) the two spheres' sediment perpendicular to the line connecting their centers.The generalized mobility matrix used is given by(25).

Figure 11 . 3 Figure 12 .
Figure 11.Three particles sedimenting in a viscous fluid at rest.

Figure 13 .
Figure13.Evolution of (y 2 − y 1 )/a as a function of dimensionless time t * = V 0 s t/a.The symbols (circles) show the results obtained by Ganatos et al.[13] (reproduced with permission from Peter Ganatos, Robert Pfeffer, Sheldon Weinbaum published by Cambridge University Press, 2024).

Figure 14 .
Figure14.Evolution of (x 2 − x 3 )/a as a function of dimensionless time t * = V 0 s t/a.The symbols (circles) show the results obtained using the results of Ganatos et al.[13] (reproduced with permission from Peter Ganatos, Robert Pfeffer, Sheldon Weinbaum published by Cambridge University Press, 2024)

v 3 v 2 Figure 15 .
Figure 15.Velocity of the three particles obtained using the IIR-Filter method with the corrected mobility matrix RPY c .

Figure A1 .Figure A2 .
Figure A1.Flowchart of IIR-Filter method applied to the resolution of system(1).Appendix B