## Abstract

Despite their many advantages over other numerical methods, boundary integral formulations still fail to provide accurate predictions of mesoscale motion in dense suspensions of rigid particles because the nearly singular flow between surfaces in close proximity cannot be resolved accurately. A procedure for incorporating analytical solutions for the lubrication flow within a large-scale boundary integral equation method is shown. Although the method is applied to the case of spherical particles, in conjunction with the completed double layer boundary integral equation, it can be developed further to treat more complex geometries and can be adapted to other numerical techniques. In contrast to other apparently similar approaches, the present method does not resort to effective medium approximations, and in principle retains all the advantages typical of boundary integral approaches. The framework also allows for forces other than those due to hydrodynamic lubrication between particles, provided that they are a linear function of the relative velocity or at least can be linearized; for example, forces due to sub-continuum fluid behaviour or forces resulting from surface chemistry. It is shown using several benchmarks that the relative motion between two particles in various flows is captured accurately, both statically and dynamically, in situations where uncorrected simulations fail. Moreover, the computational effort is reduced substantially by the application of the method.

## 1. Introduction

Concentrated suspensions of rigid particles are an important and interesting class of materials, whose properties are not completely understood. Dense suspensions are especially interesting, from both the technological and scientific points of view. The ability to model their flow would allow the full exploitation of the potential of these materials. For this reason, tools for investigating the constitutive behaviour of concentrated suspensions have been sought by many, through theoretical and experimental means.

Historically, monomodal suspensions of spherical particles in a Newtonian fluid have received much attention, due to the simplicity with which they can be characterized. Nevertheless, even such simple ‘model’ systems exhibit strikingly complex behaviour (see Jánosi *et al*. 1997). Although, governed by Stokes's equations, and thus reversible, strong irreversibility is observed almost invariably in experiments. Examples are clustering, particle migration and anisotropy. The numerical simulation of these systems, therefore, is of high interest and could lead to new understanding.

The mesoscale simulation of flows of suspensions to obtain macroscopic properties is not a new endeavour. There have been three principal approaches: Stokesian dynamics (SD), finite element analysis (FEA) and boundary element analysis (BEA), each method having strengths and weaknesses relative to the others.

To date, SD (Bossis & Brady 1984; Brady & Bossis 1985, 1988) is perhaps the most successful method. In this, the near-field particle interactions are treated as pairwise additive forces which are a linear function of particle velocities. The collective far–field interactions are also modelled simultaneously, generating a moderately sized linear system that can be solved quickly enough to allow large-scale dynamic simulations. An early two-dimensional version of the SD was used to determine constitutive model parameters by Nott & Brady (1994) and Brady & Morris (1997). More recently, three dimensional large-scale simulations were used to investigate the evolution of particle structure and its relationship to the rheology of the suspension by Sierou & Brady (2002). The main drawbacks of SD are that they are only suited to spherical particles, it is difficult to model container boundaries and they do not fully account for coupled intermediate-field interactions between particles.

Another important numerical tool, FEA, has been used in two- and three-dimensional simulations, particularly in connection with inertial and non-Newtonian behaviour. However, the most physically realistic simulations are three-dimensional. Because of the complexity inherent in discretization, most of these were performed with fine regular Eulerian meshes of hexahedral cells, where spherical particles are approximated by the elements contained in its geometrical definition. While these models appear to provide a good qualitative representation of the dynamic interactions between particles and surfaces (see Glowinski *et al*. 1994, 1999; Patankar *et al*. 2001; Hu *et al*. 2001), it is not clear whether close-range interactions are modelled accurately. This is especially important at high concentrations, where such interactions become increasingly dominant. Less common are studies where the geometry is tracked by Lagrangian meshes, updated at every time step. Johnson & Tezduyar (1999) used such a method to simulate the sedimentation of 1000 spheres at relatively low volume fraction. Although the systems generated in FEA simulations are much larger than with SD, they are sparse. In the case of regular meshes, solutions could be obtained efficiently with iterative solvers. In addition, the particle shapes and container wall geometry are arbitrary.

Finally, many have sought to exploit the ease with which complex geometries could be described inherent in BEA. In some cases, direct formulations, where the primary variables are tractions and velocities on the surface of particles, were used (see Ingber & Mondy 1993; Li & Ingber 1994; Ingber *et al*. 1999). In these simulations, Gaussian elimination-type linear solvers were used, principally because the more efficient iterative solvers do not converge readily due to the high condition number. For this reason, indirect formulations which result in Fredholm integral equations of the second kind have gained popularity for large systems. The completed double-layer boundary integral equation method, described by Karrila *et al*. (1989), Phan-Thien & Kim (1994), Phan-Thien & Tullock (1994) and Power & Wrobel (1995), where the solution is sought in terms of a double-layer hydrodynamic potential in combination with a series of forces and torques located at particle centres, has gained widespread acceptance. Solutions with iterative solvers such as GMRES (Saad & Schultz 1986) can be obtained with small numbers of iterations, independently of the problem size. The fact that Krylov subspace iterative solvers do not require explicit knowledge of individual terms in the coefficient matrix also allows the use of Barnes–Hut (Barnes & Hut 1986) multipole methods, or the more effective fast multipole methods of Greengard & Rokhlin (1988), capable of reducing the operation count from (*N*^{2}) to (*N* log *N*) and (*N*), respectively. Thus, the principal drawback of BEA with respect to FEA, namely the difficulty in solving the fully populated linear systems inherent in BEA, has largely been overcome.

However, Mammoli (2002) showed that, despite the accuracy with which arbitrary geometries could be represented, the macroscopic properties predicted by BEA were reliable only for low to moderate concentrations of the rigid phase (up to *ϕ*=0.30). The possible causes for this are two. First, the particle distributions are not random, but a function of the flow. Second, the forces between neighbouring particles could be underpredicted by the solution of the boundary integral equation. The latter problem is a consequence of the inability to perform the nearly singular boundary integral in the regions of surfaces in near-contact with sufficient accuracy. Various numerical schemes were proposed to overcome this problem (see Mustakis & Kim 1995, 1998), but none provides a satisfactory solution for arbitrarily small surface separations, especially in the case of limited mesh density. An alternative possible solution for this problem, which made use of known analytical solutions for pairwise interactions between two spherical particles, was outlined by Mammoli (2002).

The importance of considering lubrication forces when dealing with short-range interactions has already been recognized elsewhere. A seemingly similar correction procedure for BIE formulations was proposed by Qi *et al*. (2000). The present formulation, however, is more basic, and does not require the modification of the fluid viscosity using effective media relations, which inevitably biases the problem solution as a consequence of the viscosity–concentration relation chosen. Lubrication corrections were also considered recently by Dance & Maxey (2003), in the context of a force coupling method (FCM) formulation. In what follows, the methodology is formalized, and incorporated within the framework of the completed double layer boundary integral equation (CDL-BIE), giving rise to the corrected CDL-BIE (CCDL-BIE). The results from the new formulation are compared to several static and dynamic analytical benchmarks. The flows are analysed in detail, and some implications of the results are discussed.

## 2. Formulation

### (a) The CDL-BIE

For large problems, a boundary integral formulation that results in a second-kind Fredholm integral equation is desirable, because it facilitates the use of iterative linear solvers and acceleration methods to solve the fully populated linear system that results from the discretization of the BIE.

In the present work, the CDL-BIE is used as a basis for the development of the proposed method. The derivation of the CDL-BIE has been discussed at length in the literature (see Kim & Karrila 1991; Kim & Power 1993; Phan-Thien & Kim 1994; Power & Wrobel 1995). The resulting equations are summarized below.

Let *Ω* be the fluid volume between a container and *n* particles, and let *Γ*_{0} be the surface of the container. Also denote the surface of particle *p* by *Γ*_{p}, and the totality of the particle surfaces by *Γ*. The velocity at a point is(2.1)where the kernel function is(2.2)in which *r*=|* x*−

*| and*

**y***is the unit normal to the surface pointing into the fluid domain. Because the double layer representation can only describe flow fields generated by force- and torque-free particles, a series of Stokeslet/Rotlet pairs (the second term on the right-hand side of (2.1), also known as the range completer) each originating from point*

**n**

**y**^{p}inside each particle

*p*, is necessary to ‘complete’ the flow field. In these, the Oseen tensor is given by(2.3)where , while each

**F**^{p}and

**T**^{p}is the applied force and torque on particle

*p*, known in the mobility problem.

Using a limiting process to bring the point * x* to the boundary, and assuming no-slip boundary conditions on all surfaces, the boundary integral equation for is(2.4)where

*is the velocity on the boundary. This is a Fredholm equation of the second kind for the unknown double layer density*

**U***ϕ*. In the mobility problem, the velocity on particle surfaces is unknown, but it can be related to the translation

**u**^{p}and rotation

*ω*

^{p}about point

**y**^{p}inside particle

*p*to obtain(2.5)where , . Because

**u**^{p}and

*ω*

^{p}are unknown, an additional six degrees of freedom per particle are introduced. Furthermore, the integral equation admits a double layer potential field given by rigid body motion at each of the particle surfaces and a velocity field equal to the container surface normal as eigensolutions. One method to simultaneously constrain out the particle eigensolutions and remove the rank deficiency in the system is to link the rigid body motion of each particle with the double layer potential on the same particle, as suggested by Kim & Karrila (1991). Another possibility is to relate the double layer density to the force and torque exerted on each particle, with(2.6)and(2.7)This is the method usually chosen for resistance problems, where the forces and torques exerted on the particles are unknown. In the present formulation, it is applied to solve the mobility problem because it results in a matrix structure which can be manipulated more easily for preconditioning purposes, particularly when the lubrication correction is introduced, as will be shown later. The container eigensolution is removed by adding the term to (2.4) and setting(2.8)

In the case of a collection of particles in an unbounded fluid, the velocity on the particle surfaces is given by(2.9)where is the velocity that would exist at * x* in the absence of particles. The system can be solved using a GMRES solver with preconditioning based on the known matrix structure, as discussed by Mammoli & Ingber (1998).

The particle velocities obtained from the solution of the BIE allow the tracking of the particle positions, and consequently the evolution of the particle structure as a function of the imposed flow is possible. Complementary to this is the determination of the macroscopic constitutive behaviour of the suspension. Together with the deformation rate that is applied to the RVE, determination of the resultant stress tensor, in a volume-averaged sense, yields the constitutive response. The volume-averaged stress tensor can be obtained from the double layer potential solution by matching coefficients of like terms in the multipole expansions of the velocity field from a direct formulation, as suggested by Kim & Karrila (1991) and Phan-Thien & Kim (1994).

### (b) Lubrication forces between two particles

It is evident from previous benchmarks (see Mammoli 2002) that the numerical solution of the BIE does not adequately capture the lubrication flow between neighbouring particles. To a good approximation, it can be assumed that lubrication flow between two nearly touching particles is due to the relative motion of the particles only, and is not affected either by the mean velocity of the particle pair or by the presence of neighbouring particles.

Any relative motion between two particles can be decomposed into four basic motions, illustrated in figure 1: rotation about an axis normal to the line of centres (mode I), translation along an axis normal to the line of centres (mode II), rotation about the line of centres (mode III) and translation along the line of centres (mode IV).

Analytical solutions which give the forces resulting from the motion of a spherical particle in the vicinity of a second particle exist for absolute motions corresponding to each of the four modes of motion, in the form of infinite series. These were developed by O'Neill and Majumdar (1970*a*) for modes I and II, by Jeffery (1915) for mode III, and in various forms by Stimson & Jeffery (1926), Brenner (1961) and by Cooley & O'Neill (1969) for mode IV. For small interparticle separations, it is advantageous to use asymptotic solutions, which were developed based on the infinite series solutions (see O'Neill and Majumdar 1970*b*; Jeffrey & Onishi 1984). The solutions can be cast as the forces and torques on each particle that result from their combined motions. In what follows, the notation is altered from the original source to render it common for all types of motion.

Consider first motion of type I, for two particles aligned in the *x*_{1}-direction. Particle A, of radius *r*_{A}, rotates about an axis through its centre parallel to the *x*_{3}-direction with angular velocity and particle B, of radius *r*_{B}, rotates about an axis through its centre parallel to the *x*_{3} direction with angular velocity . The flow field set-up by the particles' motion results in forces on each particle along the *x*_{2}, and couples in the *x*_{3} direction, given by(2.10)where *k*=*r*_{B}/*r*_{A}, *ϵ*=*d*/*r*_{A}, with *d* taken as the interparticle gap (the minimum distance between the spherical surfaces). The coefficients and are interpreted as follows: the subscript *i* is 1, if the target particle, for which the force or couple is being calculated, is the same as the source particle, whose motion generates the force or couple, 2 otherwise. The superscript refers to the relative motion mode (I in this case). For example, in the first equation in the above set, represents the force on particle A due to its own rotation combined with the presence of particle B, normalized by the magnitude of the Stokes force of the same isolated particle moving in an infinite fluid with speed .

Similarly, in mode II relative motion, with particle A moving in the *x*_{2} with velocity and particle B moving in the with velocity , the flow field applies forces to each particle in the *x*_{2}-direction, and couples in the *x*_{3}-direction. These are given by(2.11)As will become evident later, for the purposes of this work, only a portion of , , and given in (2.11) are of interest, namely the portion produced by the relative motion.

In mode III motion, the particles rotate about their line of centres, in the *x*_{1}-direction for the case illustrated in figure 1. The angular velocities of particles A and B are and , respectively. The resultant couples in the *x*_{1}-direction are given by(2.12)

Finally, in mode IV the spheres approach each other along their line of centres. In the case illustrated in figure 1, the spheres A and B are aligned in the *x*_{1}-direction, with velocities and . As before, only forces generated by relative motion are of interest. These are in the *x*_{1}-direction, and their magnitudes are given by(2.13)

The coefficients *f* and *t* in equations (2.10)–(2.13) are given in appendix A. It should be noted that the *O*(1) and *O*(*κ*) terms that appear in the expressions for the force and couple coefficients cancel out in the scheme that will be developed in the following sections, while the *O*(*e*) terms are assumed negligible.

### (c) Application of the lubrication correction to the BIE

Consider a suspension of monomodal spherical particles. The singular flow field in the volume of fluid between two surfaces in near-contact (interparticle gaps smaller than approximately 3×10^{−2} of the particle radius) is not resolved correctly by the numerical solution of the boundary integral equation due to the practical difficulty in using a sufficiently fine surface discretization. However, it is possible to utilise the analytical solutions between two isolated particles, discussed in the preceding section, to correct for the lack of resolution in the lubrication region.

The lubrication solution cannot simply be added to the numerical solution of the BIE because it is already captured, albeit inaccurately, by the solution of the BIE. The part of the singular flow implicit in the solution of the BIE must be identified, removed and replaced with the analytical solution. This is accomplished by adapting a technique suggested by Zinchenko (1998) for the simulation of heat transfer in a packed bed. With reference to figure 2, given a collocation point on particle A, the surface of particle B is modified so that a minimum critical gap between the two exists. The critical gap is set so that an acceptably accurate solution is achieved with a specified mesh refinement level. The accuracy of the numerical integration on each constant element is ensured by an adaptive Gaussian integration algorithm which uses a coordinate transformation where the Jacobian becomes small in regions of near-singularity (see Telles & Oliveira 1994). The lubrication force for the modified geometry, which is now correctly predicted by the BIE solution, is removed using the analytical lubrication force for the modified geometry and replaced with the analytical lubrication force for the original geometry. The discretized version of the ‘corrected’ BIE for an unbounded flow and a collocation node on particle *k* then becomes(2.14)where is a modified boundary, similar to the original boundary except that the neighbours of particle *k* are shrunk sufficiently to obtain particle separations equal to the critical gap *ϵ*_{crit}; *N*(*p*) is the number of neighbouring particles to *p*, and is the *i*th neighbour of particle *p*. The force and represent the analytical lubrication forces applied to particle *p* by its neighbour *q* for the original and modified geometry, and similarly for the torques.

In the calculation of the lubrication forces for a particular particle pair, it is convenient to operate in a local coordinate system defined by the relative position of the particles, as shown in figure 3. In this coordinate system, it is simple to show that the lubrication force and torque to which particle *p* is subjected, due to the presence of neighbouring particle *q*, can be written in terms of a linear mapping of the respective velocities: (2.15)(2.16)where the coefficient matrices and can be deduced easily from the expressions in §2*b*, and are presented in detail in appendix A. The local rates of translation and rotation can be expressed in terms of global coordinates by the coordinate rotation R^{pq} for the particle pair *pq*:(2.17)(2.18)Forces in the global coordinate system can be obtained by using the inverse rotation, so that finally one has(2.19)(2.20)where, for example, . This can now be substituted in (2.14), which in turn can be discretized and solved as normal. It must be noted here that the treatment of the lubrication forces represented in (2.19) and (2.20) is where the principal difference between this work and the work of Qi *et al*. (2000) is: in the present work, the lubrication forces on each particle are cast as variables, and obtained simultaneously with the entire problem, while in the work cited, a separate, decoupled force balance is solved to calculate the lubrication forces on each particle, which are then applied to the particles in the BEM problem. The force balance includes the force on an isolated particle resulting from its isolated motion in a fluid with effective properties obtained using an equation for the effective viscosity of a suspension. The advantage of the present method is twofold: lubrication forces are completely coupled with the entire problem, and no effective medium approximation is required.

It was postulated that the lubrication correction must be solely a result of the relative motion between two particles in a pair. Inspection of the expressions for coefficients relating the forces and torques on particles to their motion, given in (A 1)–(A 3), reveals that the expressions yield the total forces acting on a spherical particle in the presence of the other particle, and not just the forces due to the lubrication flow between particles. All effects other than those resulting from the lubrication flow must be stripped from the formulation, since they are already calculated accurately by the BIE solution. One possible method for doing so is by altering the expressions for the calculation of the coefficients at the source. A less direct but simpler method for obtaining the same result is to ensure that the forces yielded by the multiplication of the particle velocities with the lubrication matrices are independent of the translation or rotation of the reference frame.

Consider first a translating frame of reference. Its effect can be removed by subtracting the centroid velocity from the velocity of the each particle in the pair. For simplicity, the case of two identical particles is considered, however, similar considerations can be made for an arbitrary size ratio *k*. Clearly, only the matrices that multiply the velocity vectors are affected by this manipulation, since angular velocities are unaltered by a translating reference frame. Thus, the term in (2.19) is replaced with . Given that, for two equal particles, , it is easy to show that(2.21)The torque matrices are treated identically.

The second type of motion which does not generate any lubrication force, and which corresponds to frame invariance with respect to a rotating frame of reference, is the motion of a rigid dumbbell. Consider two particles, *p* and *q*, with relative velocity and separated by a distance * r*. The part of which is normal to

*is given by(2.22)where*

**r***r*is the magnitude of

*. The angular velocity*

**r**

*ω*_{d}which, in conjunction with , constitutes a rigid dumbbell motion, can be related to by(2.23)where the coefficients of the matrix D can be obtained with some tensor algebra.

The term in (2.19) is replaced by . By substituting the results of (2.23), the removal of the dumbbell angular velocity can be expressed in terms of quantities which appear explicitly in the matrix equations, namely(2.24)Similar treatment is given to the torque equation. Thus modified, the lubrication correction matrices produce lubrication forces that are frame invariant, while remaining in terms of quantities that are part of the original solution vector.

Finally, it is recognized that the lubrication correction forces presented here are only applicable for spherical particles, and that this is somewhat counter to the general spirit of BIE formulations, for which applicability to arbitrary geometry is cited as a strength. Analytical forms of the lubrication forces and torques for arbitrary smooth geometries exist in the literature (see Claeys & Brady 1988; Cox 1973), and can be applied in the present framework. For geometries with sharp edges or corners, analytical solutions do not exist to the best of the author's knowledge. However, empirical relationships could be established with the aid of numerical simulation, as suggested, for example, by Mustakis & Kim (1995, 1998).

### (d) Preconditioning

The linear system that results from the discretization of (2.5) or (2.14) is of the form(2.25)where * Φ* is the vector of nodal double layer densities,

*is the vector of centroidal velocities (translational and rotational),*

**u***is the vector of velocities due to the point forces and torques,*

**g***is the vector of applied forces and torques, and is the integral operator. One of the features of the correction to the CDL-BIE discussed here is that the large-scale matrix structure for the solution of the mobility problem, (2.25), remains unaltered. However, with the lubrication correction the sub-matrix U no longer consists of a diagonal array of 6×6 matrices, but is in fact fully populated. This change in sub-structure has implications in the formulation of a preconditioner. As noted previously, the principal reason for using CDL-BIE is the low condition number of the linear system that results from its discretization, which makes the system suitable for the application of very efficient iterative solution techniques, such as the combination of multipole acceleration with a GMRES solver. The presence of the U and F matrices disturbs the diagonal dominance of the system, and the number of iterations required to obtain a solution may increase substantially. Previous work (see Greenbaum*

**f***et al*. 1992, 1996; Mammoli & Ingber 1998) has shown that a left preconditioner of the form(2.26)is a good approximation to the inverse of the mobility matrix and is very effective at restoring the diagonal dominance of the CDL-BIE. The application of the preconditioner on a vector is described by(2.27)It is readily apparent that the formation and inversion of the matrix FU is the most computationally expensive operation in the preconditioning process. As such, its parallel implementation is desirable. To understand the process, it is instructive to inspect the structure of the F and U matrices. The F matrix results from the discretization of (2.6) and (2.7). In both the original and CCDL-BIE, it takes the form(2.28)where the submatrices F

^{i}have dimension 6×

*M*

^{i},

*M*

^{i}being the number of degrees of freedom on particle

*i*. The structure of the matrix U, on the other hand, is affected by the correction procedure. The form of the uncorrected and corrected matrices, U and U

_{C}, respectively, is(2.29)where each submatrix U

^{i}or U

^{ij}has dimension

*M*

^{i}×6. The product FU is a diagonal array of

*n*6×6 matrices. Its inverse is simply the array of inverses of the original matrices. It is ideally suited for parallel evaluation and inversion, with minimal operation count and communication overhead, because the preconditioner is calculated and applied locally.

The calculation becomes substantially more complex in the case of the preconditioner for the ‘corrected’ CDL-BIE. The matrix U_{C} is potentially fully populated because all particles (except for isolated ones) are subject to lubrication correction forces from their neighbours, which must be accounted for in the range completer, as can be inferred from (2.14). The product FU_{C} is also a fully populated matrix. The parallel evaluation of the matrix coefficients and the inversion of the matrix are now both more desirable and more difficult. The evaluation of the matrix coefficients is parallelized by splitting the number of particles as evenly as possible between processes. On each process, the submatrices F^{i} are evaluated for the particles allocated locally. Likewise, only the row blocks of the matrix U_{C} which correspond to the collocation nodes on the locally allocated particles are evaluated in each process. The row block of the product FU_{C} can be evaluated locally by multiplying the partial rows of the matrix F with the partial columns of U_{C}. The matrix FU_{C} has dimensions 6*n*×6*n*, substantially smaller than the overall dimensions of the CCDL-BIE system. However, if the number of particles *n* is large (thousands of particles) the evaluation of the inverse can still be computationally expensive, and should be done in parallel by using a suitable direct solver. For small to moderate values of *n*, say up to 100, it is simpler to gather the entire matrix to a single process, perform the inversion on this process, and scatter the results back to the individual processes. This latter option is used for the purposes of the work described here.

## 3. Benchmarks

### (a) Static

Analytical solutions are available for several two-particle interactions. Those used here are depicted in figure 4. Forces or torques are applied to the centres of two identical spherical particles, and the resulting motion is calculated. The motion predicted by the BIE solution, at various levels of discretization, with and without the asymptotic correction, is compared to the motion predicted using the analytical solution.

The first benchmark problem is the motion of two particles under the action of opposing forces along the line of centres, which in this case is parallel to the -axis. As a result of the forces, each particle approaches the other along the -direction with a velocity of magnitude *u*. The ratio of *u* with the speed *u*_{∞} of a single particle in an unbounded fluid under the action of the same force is shown in figure 5*a*, while the corresponding error is shown in figure 5*b*.

Several points emerge from this benchmark. Homogeneous meshes of 24 and 96 elements do not provide sufficient accuracy for *ϵ*<1.0. A homogeneous refinement of 384 elements is slightly less accurate than the adaptive mesh at *ϵ*<0.1, which uses 456 elements. Gaps much smaller than *ϵ*<0.1 are prevalent in dense suspensions (see Sierou & Brady 2002), and it is clear that obtaining accurate results for such simulations would require excessive mesh refinement. It is also apparent that if one is to maintain a reasonable computational cost, the standard BEM solution is accurate down to gaps larger than 0.04 of the particle radius, where the error in the relative velocity is approximately 6% with an element size to local gap ratio *s*/*d* of 0.5. On the other hand, the lubrication correction is very effective at small separations, where the error is very small. The largest error arises at moderately small gaps. For example, the error at *ϵ*=0.0316 is 22%. This is due in part due to the increasing inaccuracy of the BEM, and partly to the geometry modification. The problem is that the geometry modification affects the drag force in the intermediate field, and this is not corrected by the analytical solution. Thus, in order to gain orders of magnitude in accuracy for small gaps, one must accept a modest loss in accuracy for moderately small gaps. The consequences of this error on the predicted particle trajectory will be discussed in the dynamic benchmark section.

In the second benchmark problem, forces are applied to the centre of each particle in opposite directions along lines normal to the line of centres (the -direction in this case). The resulting motion is translation of the particles parallel to the -axis with velocity *u* and rotation about the -axis with angular velocity *ω*. The analytical solution, determined using O'Neill and Majumdar (1970*a*,*b*), is used to benchmark the numerical solution, again for various values of the interparticle gap and various discretization levels.

In the limit of a vanishing gap, the motion of the particles asymptotes to the motion of a rigid dumbbell. For two identical particles of unit radius separated by a gap *ϵ*, the following holds:(3.1a)(3.1b)where *A*(*ϵ*), *B*(*ϵ*), *C*(*ϵ*), *D*(*ϵ*) can be written in asymptotic terms for small *ϵ*. It can be shown with some algebraic manipulation that, although the individual terms *A*(*ϵ*), *B*(*ϵ*), *C*(*ϵ*), *D*(*ϵ*) are individually divergent as the gap *ϵ* vanishes, various combinations thereof are finite. In particular, it can be shown that(3.2)By manipulating (3.1*a*) and (3.1*b*), it can be shown that *δ*=*u*−*ω* behaves as for small *ϵ*. In other words, the angular velocity and the linear velocity converge slowly, as the reciprocal of ln(*ϵ*). By re-writing (3.1*a*) as(3.3)and noting that *B*(*ϵ*) behaves as −ln(*ϵ*)/3 for small *ϵ*, it is seen that the limiting linear and angular velocity for a vanishingly small gap is . This limiting velocity is subtracted from the analytical and numerical solutions, and the results are plotted in figure 6. Translation is normalized by the velocity *u*_{∞} of a single particle in an unbounded fluid, rotation by the angular velocity *ω*_{∞} due to a unit torque.

The error in the translation and rotation induced by the shearing motion is smaller than the error for squeezing motion, especially when using the adaptive mesh with element size to gap ratio of 0.5. However, the number of elements is still excessive and would result in inefficient simulations. In addition, because the geometry modification must be performed to correct the error in the squeezing motion, it is necessary to introduce the asymptotic correction for all other modes of motion to compensate for the geometry modification. The results for the corrected method are in good agreement with the analytical solution, and as in the case of squeezing flow, improve with diminishing gap. Again, the error is largest when the gap is moderate, approximately 0.03 of the particle radius.

In the third benchmark, the twisting motion resulting from the application of a torque along the line of centres is examined. It should be noted that this motion is not singular, and the conventional numerical solution of the BIE provides accurate results. The flow for two equal particles rotating in opposite directions at equal rates can be substituted with the flow of one particle in the half space delimited by a planar surface located halfway between the two particles. When the particle touches the surface, the torque required to sustain the twisting motion is 1.2 020 569 times the torque to sustain the same motion in an unbounded fluid (see Jeffrey & Onishi 1984). For a given torque, the angular velocity of a particle decreases with decreasing interparticle distance. The limiting angular velocity is subtracted from the particle velocity, and the difference is normalized by the angular velocity in an unbounded fluid. The results are shown in figure 7.

In fact, even with 24 elements, the error in the rotational velocity at the smallest gap is less than 10%. With adaptive mesh refinement, this error is reduced to approximately 1%. Using the lubrication correction, and reiterating the fact that the correction must be applied because of the squeezing flow, the error is less than 1% for the smallest gap.

The final static benchmark involves the settling motion of two particles under the action of gravity, with the line of centres oriented at an angle *π*/4 from the direction of gravity. The analytical solution for this situation was derived by Goldman *et al*. (1966).

In the case, where the particles centres are in the *x*_{1}*x*_{3}-plane and gravity acts in the *x*_{1}-direction, the motion will consist of a settling velocity *u*_{s} in the *x*_{1}-direction, a drift velocity *u*_{d} in the *x*_{3}-direction, and the two particles will rotate along the *x*_{2}-axis at equal rates *ω* and opposite directions. The particle centres have no relative velocity, thus the only singular flow is of the mode I type, and gives rise to the rotation of the particles. The settling and drift velocities are not singular and are captured adequately by the numerical solution of the BIE. In figure 8, the predicted drift velocity and rotation rate are plotted versus the analytical solution. The settling velocity is not plotted because the error there is negligible in all cases.

It is seen that the rotation is captured very well by the corrected method, while with the uncorrected method even extremely fine discretizations fail to predict it accurately. The drift velocity, on the other hand, is adequately predicted by all methods, and in some cases the uncorrected BEM provides better results. As stated previously, the lubrication correction is not expected to improve the drift velocity prediction, because this is not affected by the near-singularity in the flow. The modification of the geometry necessary for the lubrication correction is responsible for the small loss in accuracy observed. In the settling velocity, the error is very small compared to the magnitude of the velocity and can be neglected.

### (b) Dynamic

The first dynamic benchmark concerns the interaction of two particles in a linear shear field. The analytical solution for this case was first obtained by Batchelor & Green (1972). An approximation of this solution was obtained more recently by da Cunha & Hinch (1995). The latter solution is used here as a benchmark for the CCDL-BIE code. The shear field is described by *v*_{1}=*x*_{2}. The particles are located at initial positions (−5.0, 0.1, 0.0) and (5.0, −0.1, 0.0) in one case, and (−5.0, 0.5, 0.0) and (5.0, 0.3, 0.0) in the other. The relative motion of the two particles should be identical, since the flow field experienced by the particles in the second case is equivalent to the flow field experienced in the first case, as viewed from a reference frame moving with a velocity *v*_{1}=−0.4. The relative trajectories are shown in the central panel of figure 9, where the semi-analytical results obtained using da Cunha and Hinch's procedure are also plotted. Within numerical error, the two CCDL-BIE results are indeed identical to each other. In addition, they are very close to the semi-analytical trajectory, in contrast to the case for the uncorrected CDL-BIE, where the dynamic simulation results in the interference of the particle surfaces.

The second benchmark simulates two particles located in a fluid in rigid body rotation, which is equivalent to a stationary flow as viewed from a rotating frame of reference. The angular position of the particles, and the gap between them, are plotted against the angular position of the dumbbell formed by the two particles, in the rightmost panel of figure 9. Particle angular positions are equal to the overall rotation of the fluid, as expected, and the interparticle gap remains constant, confirming the independence of the correction method with respect to a rotating frame of reference.

## 4. Multi-body interactions

One of the often cited reasons for the irreversible behaviour observed in flowing concentrated suspensions of rigid spherical particles is the fact that in reality particles are not perfectly smooth. Their roughness, it is said, may be responsible for an asymmetry in the approach/separation phase of the interactions. It was shown by da Cunha & Hinch (1995) that such interactions may indeed result in particle self-diffusivity. From the results of the two-particle interaction in a shear flow presented in §3, it is clear that particles whose surface roughness is smaller than 10^{−5} should behave reversibly. This study does not account for the possible effect of other particles in the vicinity of the interacting pair, as da Cunha & Hinch (1995) clearly acknowledge. As a consequence, the validity of the results presented can only be proven in the dilute limit.

Using the CCDL-BIE, it becomes possible to study the effect of additional particles. In particular, it is possible to calculate the minimum gap that results from the interaction of multiple particles, which would then determine the surface smoothness required to maintain reversibility. In the case of three particles in a shear flow, the configuration which produces the smallest gap is that of three aligned equidistant particles with centres in the plane of shear.

Because of the symmetric configuration, the centre of the middle particle remains fixed in the direction of the velocity gradient, while the trajectory of the other two is symmetric. A plot of interparticle gap versus angular position is given in figure 10. The minimum gap in this case is extremely small, and in most cases of practical interest far smaller than the continuum limit.

It is interesting to see the effect of imposing a surface roughness model similar to the one imposed by da Cunha & Hinch (1995), where in this case the roughness is set to 10^{−6} times the particle radius. To prevent particles from approaching each other when their separation is below this threshold, the physical mechanism of roughness is mimicked: further approach results in a large repulsive force. This effect is obtained very simply by modifying the coefficients of (2.15). In the present case, the coefficients resulting from mode IV motion were increased by a factor of 100 in case the particles are approaching and the separation is less than 10^{−6} times the particle radius, and unaltered otherwise.

The effect is visible in the plot of gap versus angular position at the microscale, and the macroscale effect, in terms of deviation from the reversible trajectory, are visible in figure 10, in comparison with their smooth particle counterparts. While this simple roughness model may help in qualitatively explaining certain irreversible behaviour observed in experimental work, it is certainly inadequate to represent the true surface phenomena. One possible improvement would be to use solvation forces, as described by Douglas Frink & van Swol (1996, 1998), which account of effects below the continuum limit due to finite size of the fluid molecules. Similarly, other microscale or nanoscale surface effects, including those modelled by nano-scale simulations such as molecular dynamics (MD) or dissipative particle dynamics could easily be included within the same framework.

## 5. Conclusions

A method for correcting the underprediction of interparticle forces in the CDL-BIE is presented. Its incorporation into the matrix structure of the CDL-BIE is such that the solution procedure differs very little from the solution of the uncorrected system. The principal difference is that the direct inversion of a preconditioner matrix of non-trivial size by a direct method occurs prior to the solution of the entire system. This matrix, however, is very small compared to the size of the entire system, by a factor of at least 200, and should not be a significant factor in the overall solution time. The CCDL-BIE, in comparison with the original CDL-BIE, has fewer degrees of freedom, a substantially smaller condition number, and far greater accuracy in most conditions.

The CCDL-BIE was tested on several benchmark problems, static and dynamic. In all cases, the corrected method performed very well. The worst results are obtained with particles separated by moderate gaps, where the correction is small, but the solution of the BEM is slightly inaccurate. In dynamic simulations, two-particle trajectories are almost identical to those produced by the semi-analytical solution of da Cunha & Hinch (1995). In the case of three particles, for which analytical solutions are unavailable, the minimum separation is several orders of magnitude smaller than in the case of two-particle interactions, leading to the conclusion that multi-body interactions are much more sensitive to surface characteristics such as roughness, and that nanoscale surface–surface interactions should be taken into consideration. In fact, using continuum solutions for the simulation of multi-particle systems is almost certainly aphysical in practice, and serves to explain why irreversibility is always observed experimentally. The introduction of a tiny imperfection in a three-particle interaction results in a strong macroscopic irreversibility manifested by particle self-diffusion.

The present modifications to the CDL-BIE can be viewed as the combination of SD in the near-field with BEA in the intermediate and far-field. This combination weds the advantages of the two methods, namely the near-field accuracy of SD with the flexibility offered by BEA. The CCDL-BIE represents a substantial step towards the ability to perform dynamic simulations of the flow of concentrated suspensions. For a given accuracy, the improvement resulting from the lubrication correction is of orders of magnitude, as a consequence of the better condition number of the matrices, and the reduced number of degrees of freedom.

Finally, the lubrication correction algorithm could be easily adapted to numerical simulation methods other than BEA, for example FEA. Moderate Reynolds number flows could be investigated without modification of the lubrication solutions. Nonlinear fluid constitutive behaviour, on the other hand, would require the derivation of appropriate expressions. The derivation of analytical expressions could also be supplanted by MD simulations, which may be more appropriate than continuum analytical expressions for the length-scale typically encountered in reality.

## Acknowledgements

The concept of the analytical correction for close range interactions was inspired by a suggestion by Dr Alan Graham, of Los Alamos National Laboratory, and by further conversations with Prof. Alexander Zinchenko, of the University of Colorado at Boulder. These interactions were instrumental in the development of the theory and analysis code presented here. This work was supported in part by the United States Department of Energy (DOE) grant DE-FG03-97ER25332. This financial support does not constitute an endorsement by the DOE of the views expressed in this paper. Computations were performed on the LosLobos supercluster at the University of New Mexico Center for High Performance Computing.

## Footnotes

- Received March 16, 2005.
- Accepted October 31, 2005.

- © 2006 The Royal Society