Computational Mechanics 26 (2000) 388±397 Ó Springer-Verlag 2000

Lubrication approximation in completed double layer boundary element method S. Nasseri, N. Phan-Thien, X.-J. Fan

388

Abstract This paper reports on the results of the numerical simulation of the motion of solid spherical particles in shear Stokes ¯ows. Using the completed double layer boundary element method (CDLBEM) via distributed computing under Parallel Virtual Machine (PVM), the effective viscosity of suspension has been calculated for a ®nite number of spheres in a cubic array, or in a random con®guration. In the simulation presented here, the short range interactions via lubrication forces are also taken into account, via the range completer in the formulation, whenever the gap between two neighbouring particles is closer than a critical gap. The results for particles in a simple cubic array agree with the results of Nunan and Keller (1984) and Stoksian Dynamics of Brady et al. (1988). To evaluate the lubrication forces between particles in a random con®guration, a critical gap of 0.2 of particle's radius is suggested and the results are tested against the experimental data of Thomas (1965) and empirical equation of Krieger-Dougherty (Krieger, 1972). Finally, the quasi-steady trajectories are obtained for time-varying con®guration of 125 particles.

1 Introduction Suspensions provide an economical way of transporting large quantities of solid particulate material in industry as well as in medical ®eld (e.g., pulp handling in paper manufacture, petroleum processing in ¯uidized beds and self diffusion of cells in blood). However, suspension mechanics are not well understood, mainly because of the multi-particle interactions. Exact solutions for two spheres in a shear ¯ow are available in the literature. Lin et al. (1970), obtained an exact solution of the Stokes equations for the motion of two spheres of arbitrary size and arbitrary orientation with respect to a shear ®eld by using spherical bipolar coordinates. Arp and Mason (1977), presented a general method of calculating forces, torques and translational and rotational velocity components of a pair of equally-sized rigid spheres in a viscous ¯uid undergoing uniform shear ¯ow. The method is based on the matrix formulation of the hydrodynamic resistances originally documented by Brenner and O'Neill (1972).

S. Nasseri (&), N. Phan-Thien, X.-J. Fan Department of Mechanical & Mechatronic Engineering, Sydney University, NSW 2006 Australia

Far-®eld interactions for three or more spherical particles have been considered by Mazur and Van Saarloos (1982). These authors developed a general scheme for evaluating the mobility tensor for any three-dimensional con®guration of spheres and derived explicit expressions up to order R 7 , where R is the interparticle spacing. Their expressions are identical with those previously obtained by Kynch (1959) using the method of re¯ections. Applications of the theory have been reported by Kamel and Tory (1989) for special planar con®gurations of sedimenting spheres and by Ladd (1988) for periodic con®gurations. Durlofsky et al. (1987), developed a very ef®cient simulation technique, termed the ``multipole moment method,'' which permits the rapid evaluation of the force and torque on ®nite clusters of spherical particles. The method accounts for the many body interactions at large spacing (in a manner equivalent to the method of re¯ections) and lubrication forces at close spacings by considering pair interactions. The sphere mobility matrix is ®rst formed by expanding the integral formulation for Stokes ¯ow for a J sphere system in conjunction with FaxeÂn's law for the particle velocities in the moments of the force distribution on the surface of each particle. The mobility matrix is inverted to yield a far-®eld approximation to the resistance matrix. Then lubrication is introduced in a pairwise additive manner using the exact two-body resistance calculated by Arp and Mason (1977). However, the method does not readily permit evaluation of the higher order terms in the series expansion for the velocity ®eld, nor does it readily permit evaluation of the ¯uid velocity ®eld. The method has been applied to study the resuspension of spheres in shear ¯ow in a channel ¯ow by Durlofsky and Brady (1989). Hassonjee et al. (1988, 1992), developed a method of solution, called the ``multipole truncation technique,'' for treating three-dimensional clusters of arbitrary sized spherical particles and accounting for the multiparticle interactions at any given spacing. Their study addressed the mobility problem for a multiparticle con®guration settling freely under gravity and the resistance problem for clusters of particles ®xed in space in a uniform ¯ow (1988) as well as the mobility and resistance problems in shear or Poiseuille ¯ow (1992). The multipole truncation techniques uses a linear superposition of Lamb's spherical harmonic solution of the creeping ¯ow equations capable of describing an arbitrary disturbance on the surface of a sphere. The solution is expressed in spherical coordinates in the form of a Fourier series and using the orthogonality property of the eigenfunctions in the azimuthal direction

to satisfy the no-slip boundary condition on the surface of each particle. The method can be used to obtain solution of multiparticle con®gurations for up to 64 spheres ®xed in space in a uniform ¯ow and can also be used to obtain the ¯uid velocity ®eld around multiparticle con®guration. Recent studies have led to a partial understanding of the shear-induced migration of particles across streamlines in monodisperse suspension of non-colloidal particles at low Reynolds number. Geometries which have been studied include tube (Chapman, 1990), channel (Leighton and Acrivos, 1987b; Koh et al., 1994), Poiseuille ¯ow, wide-gap Couette ¯ow (Abbott et al., 1991), Couette reservoir ¯ow (Leighton and Acrivos, 1987b; Chow et al., 1994), resuspension ¯ow (Leighton and Acrivos, 1986; Chapman and Leighton, 1991; Scha¯inger et al., 1990) and parallel plate ¯ow (Chapman, 1990; Chow et al., 1994). In general, particles are found to preferentially concentrate in regions of lower shear rates and/or lower curvature. Sangani and Mo (1994) and also Mo and Sangani (1994) have developed a technique based on the periodic Green's function for the Stokes equations, which can be used to calculate transport properties in random suspensions. Stokesian dynamics (Brady and Bossis, 1988) includes exact two body interactions and results in an O…N 3 † problem (with N being the number of particles involved). This method requires a FaxeÂn relation and lubrication results for its formulation, however, as FaxeÂn laws are not available for shapes other that spheres or ellipsoids, the method seems limited in application to particles of other geometries. Therefore, when the particles shape is complex, or when the number of particles is moderately large, the boundary element method seems to be very popular in recent years, because of a reduction in the dimensionality of the problem, and the ease of mesh generation (only surface meshing is required without the need of remeshing at subsequent time steps). There are two boundary element formulations: the direct formulation where the variables of interest are the surface velocity and traction ®elds (Youngren and Acrivos, 1975; Rallison, 1981, 1984; TranCong and Phan-Thien, 1989; Phan-Thien et al., 1991; Ingber, 1990; Ingber and Mondy, 1993; Stone, 1994; Nasseri and Phan-Thien, 1996, 1997a, b, c), and the indirect formulation in which one works with a ®ctitious surface potential (Brebbia et al., 1984; Kim and Karrila, 1991; Pozrikidis, 1992). The direct formulation is not generally suitable for dense systems as the large number of boundary elements necessary will cause the system too illconditioned. Although Ingber et al. (1995) have reported a massively parallel implementations of a ®rst kind method which solve a dense system of 6240 equations using an iterative technique, however, the volume fraction that was reached is still very low (of the order of a few percent). The completed double layer boundary element method (CDLBEM) (Pakdel and Kim, 1991; Phan-Thien et al., 1992; Phan-Thien and Tullock, 1993; Phan-Thien and Kim, 1994a; Phan-Thien and Fan, 1995; Fan and Phan-Thien, 1997; Rider and Phan-Thien, 1997), allows solution of many body problems using an iterative strategy, which removes the problem of solving a large and possible illconditioned linear system. Furthermore, if domain de-

composition is used, the memory requirement will be reduced. This also makes the use of parallel virtual machine (PVM) message passing system feasible (Phan-Thien and Tullock, 1994; Seeling and Phan-Thien, 1994; Rider and Phan-Thien, 1997). Therefore, the velocity and stresslet of each particle are obtained which in turn can be used to calculate transport properties such as effective suspension viscosity. Fan and Phan-Thien (1997), formulated a traction based CDLBEM for periodic suspensions via distributed computing (PVM). The hydrodynamic interaction of particles at in®nity was handled by using O'Brien method (1979), using a mean ®eld values of the traction and the background ¯ow. After a de¯ation of the extreme eigenvalue 1 of the adjoint double layer operator, an iterative solution strategy was implemented, which solves for the traction ®eld on the surfaces of a group of near-by particles sequentially. The implementation was tested on a periodic suspension of spheres in a cubic array and proved to be accurate. Rider and Phan-Thien (1997), used CDLBEM, the PVM message passing system and an asynchronous iterative technique to calculate the viscosity of concentrated system of large numbers of particles con®ned by a rigid spherical container. Although their research was successful, in order to revolve the traction ®eld in the narrow gap between two neighboring particles, a large number of boundary elements were needed. In this paper, we develop a CDLBEM that incorporates the short range interaction by using lubrication forces whenever the gap between the two generic particles is less than a critical value. The formulation is benchmarked with Nunan and Keller's results (1984) for a periodic array of spheres in a simple cubic lattice. The use of PVM and domain decomposition has made this study successful (the most computationally intensive system studied was one involving 729  96 element spheres in a system of about 50% volume fraction). Using the average stresslet on the surfaces of all particles, we de®ned the effective intrinsic shear viscosity, and the results obtained for different con®gurations have been examined against the results of Nunan and Keller (1984) and Brady et al. (1988), for simple cubic array of particles, and against the results of Thomas (1965) and Krieger (1972), for randomly distributed particles.

2 Short-range interactions Knowledge of interactions between particles in close proximity is essential in understanding the suspension phenomena. To accurately determine particle interactions in a suspension, particularly a concentrated suspension, both the far-®eld and near-®eld physics must be modelled correctly. Both far and near-®eld interactions have been incorporated in the CDLBEM. However, the near-®eld interaction requires much more number of elements to capture. If a denotes the characteristic length of the particles (here radius of the spherical particle which is considered as unity), then the separation between the closest points on two particle surfaces, or gap, will be written as a, with   1.

389

390

Ty For two rigid particles in relative motion, the dominant 1 1 …4† contributions to the force and torque may be determined 8pla2 U ˆ 8 ln  : by the methods of lubrication theory. In particular, by using lubrication equations it is shown that the leading order terms in the force and torque are singular as  1 and 2.2 Squeezing motion ln  1 for squeezing and shearing motions respectively. Now consider a moving sphere (again denoted by Sa ) of radius a moving towards a second stationary sphere (of 2.1 radius b ˆ ba, and denoted by Sb ). This motion generates Shearing motion Consider a sphere of radius a (denoted by Sa ) moving in a ``squeezing'' ¯ow (see Fig. 1b). Similarly, the lubrication viscous ¯uid velocity U (along x) past a second stationary result for the force on sphere A is sphere of radius b ˆ ba (denoted by Sb ), in a direction Fz b2 …1 ‡ 7b ‡ b2 † 1 1 transverse to the sphere-sphere axis (see Fig. 1a). This  ln ‡ K…b† ˆ 6plaU  …1 ‡ b†2 5…1 ‡ b†3 motion generates a ``shearing ¯ow'' and the lubrication 2 3 4 results for the force and torque on sphere A are thus (for …1 ‡ 18b 29b ‡ 18b ‡ b † 1 detailed calculations see Kim and Karrila, 1991)  ln ‡ O…† ; 4 …5†

45b ‡ 58b2 ‡ 45b3 ‡ 16b4 † 1  ln 4  375…1 ‡ b† ‡ O…† ; …1† 4…16

Ty b…4 ‡ b† 1 ˆ ln ‡ B…b† 8pla2 U 10…1 ‡ b†2  ‡

…32

33b ‡ 83b2 ‡ 43b3 † 1  ln ‡ O…† ; 3  250…1 ‡ b† …2†

where 6plaU and 8pla2 U are the drag force and torque exerted on sphere A in the absence of the second sphere. Also, the O(1) terms A…b† and B…b† should be determined by matching with the outer solution (see Kim and Karrila, 1991). For two equal spheres …b ˆ 1†, considering the leading terms only, we will reach at

Fx ˆ 6plaU

1 1 ln ; 6 



21…1 ‡ b†

4b…2 ‡ b ‡ 2b2 † 1 ln ‡ A…b†  15…1 ‡ b†3

Fx ˆ 6plaU

…3†

where K…b† should be determined by matching with outer solution (Kim and Karrila, 1991). For equal spheres considering the leading terms only, we obtain

Fz ˆ 6plaU

1  4

1

:

…6†

Now consider two equal size spheres i and j, moving with the velocities Ui and Uj in viscous ¯uid respectively. The force exerted from particle j on i is then expressed as sq

Fij ˆ Fij ‡ Fsh ij ;

…7†

sq Fij

where and Fsh ij are the squeezing and shearing force respectively. We de®ne a set of coordinates as (lij ; pij ; qij ), with lij being the unit vector along the line of the particles' centres and the two others de®ned as y

lxij lzij dx ‡ lij lzij dy pij ˆ  1=2 y …lxij †2 ‡ …lij †2



y

…lxij †2 ‡ …lij †2

1=2

dz ;

…8†

y

lij dx ‡ lxij dy qij ˆ  1=2 ; y …lxij †2 ‡ …lij †2

…9†

where dx , dy and dz are the unit vectors of ®xed frame of coordinates. The squeezing and shearing forces between two particles is now written as

3 sq Fij ˆ pla 1 Uij
lij lij ; 2 sh Fij ˆ pla ln  1 Uij
…pij pij ‡ qij qij † ;

…10† …11†

where Uij ˆ Uj Ui . Since particle i is moving in viscous ¯uid with which an ambient ¯ow is associated, the drag force exerted on this particle is

Fi ˆ

Fig. 1. Flow produced by shearing motion of a sphere past another (a), and by squeezing motion of a sphere approaching another ®xed sphere (b)

6pla…Ui

rUT
Ri † ;

…12†

where Ri is the position vector of particle i in ®xed frame and rUT is the gradient of ¯uid velocity in shearing motion. Equilibrium requires that the total forces on particle i must be zero, hence

N…i† X

Fij ‡ Fi ˆ 0 ;

…13†

jˆ1…i6ˆj†

with N…i† being the total number of particles in close proximity of particle i. Considering all the above mentioned equations, we reach at

Fij

ln  6

1

2 4

N…j† X

Fjk

kˆ1…j6ˆk† 1

N…i† X

3

Fim 5Q

mˆ1…i6ˆm† T

ˆ pla ln  QrU
Rij ;

…14†

where Rij ˆ Rj Ri and Q ˆ …3  1 =2 ln  1 †…lij lij ‡ pij pij ‡ qij qij †. It is proposed that, whenever the gap between particles are less than a certain critical value, these particles are excluded from the CDLBEM formulation, and their effects are modelled by lubrication forces, which are obtained by solving Eq. (14). This linear system of equations can be solved by a direct Gauss elimination method. However, this quickly becomes too memory and CPU intensive to be of practical interest. In the present work, we have used an iterative strategy based on the generalized minimum residual (GMRES) procedure originally proposed by Sadd and Schultz (1983) and later developed by Shakib et al. (1989). This method which is particularly useful for solving linear nonsymmetric systems, is based on the generation of Krylov spaces. Within the GMRES iteration, phase and element-by-element preconditioning (Shakib et al., 1989) is accommodated. The algorithm allows for a domain decomposition into explicit element groups and implicit element groups with either direct or iterative solvers. For more information, the reader is referred to Shakib et al. (1989).

3 CDLBEM A group of particles are subjected to shear ¯ow at low Reynolds number. The equations of motion are Stokes equations, lr2 u ˆ rp;

r
u ˆ 0;

x2D ;

…15†

where u is the velocity vector, p is the hydrostatic pressure, l is the viscosity of suspending ¯uid and D is the ¯ow domain in which the position vector x is de®ned. These equations are typical solved by converting them into a set of boundary integral equations, through the Green functions (Pozrikidis, 1992). For the suspension problem considered here, the resulting integral equations are of the ®rst kind, and tend to be ill-conditioned. Moreover, iterative solution strategies do not work. The CDLBEM is designed so that one always has a set of second integral equations. Essentially, the methods represents the disturbance velocity ®eld by a double layer distribution (PhanThien and Kim, 1994),

ui …x† ˆ

u1 i …x†

Z

‡

Kij …x; y†uj …y†dS…y†;

where u is the double layer density vector, u1 i …x† is the ambient velocity ®eld, and the double layer kernel, K(y, x), is

K…y; x† ˆ

n…y†


S

…16†

y†…x jx

y†…x yj5

y†

;

…17†

in which n is the normal unit vector is directed into the ¯uid. Now, it is well-known that the double layer suffers a jump on S:

u1 i …x† Z ˆ ui …x† ‡ Kij …x; y†uj …y†dS…y†; x 2 S :

391

ui …x†

…18†

S

In symbolic notation

u1 …x† ˆ …1 ‡ j†u…x† ; …19† R where j…
† ˆ s Kij …x; y†…
†dS…y†. The double layer can represent only force-free and torque-free solutions (i.e., the inclusions have to be forcefree and torque-free). While this is not a critical issue, it is desirable that the representation be completely general so that solutions with a ®nite force and torque over a close surface may also be represented. The process of adding a suitable combination of particular solutions to the righthand side of Eq. (19), is known as a range completer. Power and Miranda (1987), in their investigation of the second-kind formulation of Stokes' ¯ow, suggested that one could complete the range of Eq. (19) by point-force and point-torque solutions, suitably located at each of the inclusions (see Pozrikidis, 1992; Kim and Karrila, 1991; Phan-Thien and Kim, 1994). This method has become known as the CDLBEM. With this range completer, the boundary integral equations (19) takes the form u…x†

u…x†

u1 …x†

ˆ …1 ‡ j†u…x† ‡

X F…p† p

 1 …p† G…x; x…p† † …T  r† ; 2 8pl …20†

where F…p† and T…p† are the external force and torque acting on inclusion p, respectively, and G…x; x…p† † is the Stokeslet that corresponds to a point force located at x…p† ; this latter point can be anywhere within the inclusion, but a good choice would be the centre of mass of the surface of the inclusion.1 The Stokeslet is given by

G…x; y† ˆ

1 jx

yj

1‡

…x

! y†…x y† : jx yj2

On the surface of inclusion p, the displacement ®eld is a rigid body displacement: 1

x2D ;

3 …x 2p

The centre of mass of the inclusion surface is de®ned by Z 1 x…p† ˆ …p† x dS : S S…p†

u…x† ˆ U…p† ‡ x…p†  x ;

392

…21†

where U…p† and x…p† are the rigid translational and rotational displacements of the inclusion; they are both unknowns and must be found as parts of the solution procedure. Considering one single inclusion, there are six nontrivial eigenfunctions associated with the eigenvalue at 1 of the double layer, each is associated with a rigid body motion mode of the inclusion. These nontrivial eigenfunctions (e.g., u…p;f † ; p ˆ 1; . . . ; N; f ˆ 1; . . . ; 6) are the null functions to (1 ‡ j) and solution to Eq. (20) will not be unique, as any linear combination of the solution and eigenfunctions will be another solution. To render a unique solution, we impose the constraint that the projection of the double layer density, ui , onto the null space is equal to the components of rigid body motion:

U…p† ˆ

3 X iˆ1

x…p†  r…p† ˆ

…p†

u…p;i† Fi 3 X iˆ1

ˆ

3 X

3 X

ˆ

u…p;i† hu…p;i† ; ui ;2

where

8 1 < p ei on particle p (i = 1, 2, 3), …p;i† …24† u ˆ S…p† : 0 otherwise, 8 1 > < q ei  r…p† on particle p (i = 1, 2, 3), …p† …p;i‡3† u ˆ Ii > : 0 otherwise, …25† …p†

2  dS ;

…29† De®ning a new operator as:

H…
† ˆ j…
† ‡

X

u…p;f † hu…p;f † ;
i ;

…30†

f ;p

the boundary integral equations can be rewritten as

…1 ‡ H†u ˆ b ;

…31†

where H has the same eigenvalues as j, but the extreme eigenvalue 1 has been mapped to zero (for more detailed descriptions see Phan-Thien and Kim, 1994). We will be mainly concerned with shear ¯ow

_ vx ˆ cy …vy ˆ vz ˆ 0† ;

…32†

…23†

iˆ1

ri

p

 1 …p† G…x; x…p† † …T  r†
: 2 8pl

3.1 Domain decomposition and PVM In order to solve the Eq. (31), Tullock and Phan-Thien (1992) found out that a domain decomposition would accelerate the convergence rate. In this scheme, the boundary is divided into M subdomains (m ˆ 1; . . . ; M) and the unknown vector on a subdomain m is denoted by u…m† . Considering the Picard iterative method for solving Eq. (31), we have

…p†



X F…p†

where (vx ; vy ; vz ) are the components of u1 , and c_ is the shear rate.

u…p;i‡3† Ti

r… p† ˆ x x…p† ; Z  …p† r…p†
r…p† Ii ˆ

u1

…22†

iˆ1

u…p;i‡3† hu…p;i‡3† ; ui ;

bˆ

u…n‡1† ˆ b…n†

Hu…n† ;

…33†

in which n is the iteration number and the convergence is achieved if u ˆ ……ku…n‡1† k ku…n† k†=ku…n† k†  O…10 3 †. Hence, by employing the domain decomposition procedure, Eq. (33) can be written as

u…m† ‡ H…m;p† u…p† ˆ b…m† ;

…34†

…26†

with the summation implied. For each subdomain the set of equations

…27†

…1 ‡ H…m;m† †u…m† X  …p† ˆ b…m† H…m;p† u

S…p†

where r…p† is the position vector from the particle's surface …p† centre, Ii the surface moment and S…p† the surface area of

(no sum in m† ;

…35†

p6ˆm

 ˆ …u…m† u…m 1† †=2, is solved by a standard Gauss with u elimination and cycled through all subdomain for each iteration. The process is repeated until u is less than requested tolerance. In order to have a load balancing situation, the domain decomposition should be reasonably X …p;f † …p;f † …1 ‡ j†u…x† ‡ u hu ; ui ˆ b ; …28† uniform (Phan-Thien and Tullock, 1994), so the computational time of parallel sub-tasks will be uniform across f ;p the processors, otherwise, some processors will be idle where the sum is over p ˆ 1; . . . ; N and f ˆ 1; . . . ; 6, waiting for others to ®nish their sub-tasks. nothing that the null function u…p;f † is only nonzero on the PVM (Beguelin et al., 1995), a software library for particle p, and writing and executing parallel programs based on message passing, is chosen here as a means to illustrate the application of CDLBEM. In general, processes of a parallel ap2 The angle bracket denotes: plication distributed over a collection of processors must Z communicate problem parameters and results. In distribha; bi ˆ a
b dS : uted memory multiprocessors or workstations on a s particle p. This procedure is actually Weilandt's de¯ation as described by Kim Karrila (1991), which removes the eigenvalue 1. The resulting boundary integral equations read

network, the information is typically communicated with explicit message-passing subroutine calls. To send data to another process, a subroutine is usually provided that requires a destination address, message and message length. The receiving process usually provides a buffer, a maximum length and the sender address. PVM is made up of two parts: a daemon process that any user can initiate on a machine, and a user library that contains routines. Programs written in Fortran or C languages, are provided access to PVM through calling PVM library routines for functions such as process initiation, message transmission, reception and synchronization (for more information see Geist et al., 1994; Hockney, 1994). In the serial implementation of the domain decomposition described in previous section, the boundary is divided into subdomains, each is actually the boundary surface of a particle. Each subdomain is discreatised using superparametric elements (constant in u, up to quadratic for the geometry) leading to the set of linear algebraic equations [Eq. (35)]. These equations are then solved for the unknowns on the particle m; u…m† , assuming that the solution u…p† …p 6ˆ m† is known (from previous iteration) on the subdomain. This is recycled until we have a convergent solution. The algorithm can be converted easily to a master/slave model by farming out the calculations on each particle to each processor. The master task reads in the data ®le, calculates the mesh data and the null functions for each particle, and sends all the data to the slaves [O…N† calculations]. It then sends a particle ``id'' to each slave process, which calculates the system matrix as needed and solves for the solution on this particle, using the previous solution vector for the right-hand side [O…N 2 † calculations]. This solution is then sent back to the master task, which calculates the rigid body motions and the stresslets on the particles. A check for convergence is made and the current solution vector is sent to all slaves to start a new iteration. For a pseudocode for the master and the slave program see

Phan-Thien and Kim (1994) or Rider and Phan-Thien (1997). In this procedure, updating the double layer density [ u in Eq. (35), as adopted by Rider and Phan-Thien, 1997] as soon as a slave has ®nished its most recent computational task will improve the convergence rate.

3.1.1 CDLBEM with lubrication forces The CDLBEM has all the interaction terms incorporated. However, in order to capture the near-®eld solution accurately, mesh re®nement has to be carried out in the narrow gap between two neighboring particles. A good compromise between the CPU time and the accuracy of the solution is to use the near-®eld solution in lieu of the CDLBEM whenever the gap between two particles is less than a critical value. That is, we solve Eq. (14) for the lubrication forces (neglecting lubrication torques), for the particle concerned, delete the particle from consideration in the CDLBEM, and use these lubrication forces in its range completer. We expect that this simple, but effective treatment will provide an upper bound for the force, and the viscosity of the suspension. 4 Numerical results In this paper two different con®gurations have been considered. In the ®rst one, particles have been formed a simple cubic array, whereas in second the particles are distributed randomly (see Fig. 2). The discretization scheme used by Phan-Thien and Fan (1995), has been employed here. That is, the boundary is discretized into super-parametric 9-noded quadratic elements (see Phan-Thien and Kim, 1994; Tullock, 1993). On each boundary element Se , the density function is assumed constant, so that its local approximation can be written as u…x† ˆ N e …x†ue , with ue being the constant value of the density function on the element e, whereas the geometry is quadratic.

Fig. 2. Two different con®gurations for particles subjected to shear ¯ow

393

The iterative procedure is run on a distributed computing environment (PVM farm) consisting 28 DEC Alpha 500/266, 26 have a local memory of 128 Mb each, and 2 of them each with a local memory of 512 Mb, sharing the same disk space (through NFS automount). Communication is through a fast ethernet of 100 MBits per seconds. Message passing is supported through PVM 3.3.11, that enables a collection of heterogeous computers to be used as a coherent and ¯exible concurrent computational resource. 394

4.1 Particles in simple cubic array The effective viscosity of a periodic suspension of spheres in cubic lattices has well been investigated (Nunan and Keller, 1984; Brady et al. 1988; Claeys and Brady, 1993a). In addition to providing a constitutive framework for periodic suspensions, Nunan and Keller (1984) also calculated numerically and asymptotically the ``intrinsic effective viscosity'' in shear ¯ow (using the direct formulation that takes explicit account of the symmetry of the ¯ow ®eld). This intrinsic viscosity is de®ned as hSxy i ; …36† bˆ lc_ in which l is the ¯uid viscosity (here considered as unity) p and hSij i ˆ …U=Vp †Sij , is the stresslet of one generic particle. Here we de®ne the volume fraction of an array of spherical particles (which are neutrally buoyant) by the ratio of all particles' volume to the volume of surrounding cubic cell or NVp ; …37† Uˆ V where Vp is the volume of each particle, N is the number of particles and V is the cubic cell volume. Equation (37), can be expressed in terms of the gap3 between particles, with N ! 1 and   1, resulting in Uˆ

…4=3†p ; …2 ‡ †3

…38†

in which the radius of each particle or a is considered as unity. For particles in a simple cubic array, induced by shear ¯ow, the convergence in b is initially checked against the number of particles. For low volume fraction a fast convergence is obtained, whereas for high volume fraction this convergence is very slow. For instance, when particles are formed in a cubic array with volume fraction of 34.4% (see Table 1), the convergence in intrinsic shear viscosity is obtained when the number of particles reaches 729 (9  9  9 particles). Figure 3 illustrates the intrinsic effective shear viscosity of 729 particles in different volume fractions. The present numerical results are compared to those of Nunan and Keller (1984) and Brady et al. (1988). Compared to the CDLBEM, where no lubrication force is considered, the results obtained in the present work for effective shear 3

Notice that the de®nition for gap is different with what was de®ned by Nunan and Keller (1984).

Table 1. The averaged stresslet on particles (S) and the effective shear viscosity (b) in terms of number of particles (N), when F = 34.4% (e = 0.3) N

S

b

8 27 64 125 216 343 512 729

13.54 15.49 16.45 17.36 17.85 18.20 18.56 18.67

1.113 1.273 1.314 1.327 1.376 1.396 1.401 1.403

Fig. 3. The intrinsic effective shear viscosity, obtained by CDLBEM, in terms of volume fraction, for 729 particles in simple cubic array. The results are compared with the results calculated by Nunan and Keller (1984) and Brady et al. (1988)

viscosity are generally improved for high volume fraction suspensions. Also, there is good agreement between the results obtained here for high volume fraction suspension and the singular solution …b  p4 …ln  ln…2a ‡ ††, when U ! Umax † calculated by Brady et al. (1988). Note that with the de®nition explained here, the gap they considered is actually equal to =…2a ‡ †. In general, the total CPU time depends on the volume fraction. In fact, adaptive integration techniques are more employed when particles become closer together and this will increase the slave CPU time. Also, master's issuing and updating tasks, receiving the new iterates for u, and calculating various geometrical properties will be much more when the number of elements increases. When the gap between particles is about the size of each particle's radius ( ˆ 1; U ˆ 15:5%), we noticed that the total CPU time for one iteration, considering a group of 729 particles each discretised to 96 elements, is about 2 h, whereas this time will increase (about 3 h) when the gap is reduced to 0.2(U ˆ 39:3%) of particle's radius. The size of linear system for 729  96 elements spheres is 209 952  209 952, which would require 176 Gb of memory to store the linear system (using single precision) if domain decomposition were not being used.

between particles are much more important compared to cubic array of spheres. It is found that even in moderate to low volume fractions, there may be some near-touching particles for which the lubrication theory is valid. For a system with high volume fraction, considering a very small critical gap (i.e., lubrication forces are applicable for smaller gaps than this critical one) underestimates the effective viscosity. In addition, for a dispersed system the critical gap can not be chosen very big as this will overestimates the result (see Fig. 4, dotted line, hollow circle). From these results, a critical gap of 0.2a is recommended and has been tested for three different schemes of randomly distributed particles. When suf®cient number of particles are considered (729 particles compared to 125 particles, see the same ®gure), for all volume fractions the viscosity graph is located within Thomas empirical bounds (1965). The empirical equation of Krieger-Dougherty (Krieger, Fig. 4. The intrinsic effective shear viscosity in terms of volume 1972) suggests that the effective viscosity has the following fraction, for particles in random con®guration. Numerical results form: are compared with the results obtained by Thomas (1965) and Krieger (1972)

4.2 Particles in random configuration We next consider a suspension of particles, randomly placed in a spherical cell (Fig. 2). As previous section, each particle is discretised to 96 QUAD9 elements and all the neutrally buoyant particles are subjected to shear ¯ow. In this con®guration we notice that increasing the number of particles in random con®guration will result in an increase in effective viscosity. Solving the Eq. (35), for the double layer density of each particle, we then calculate the average stresslets on all particles and use them to evaluate the effective viscosity (b). The results are tested against the experimental results of Thomas (1965) and the empirical equation of Krieger (1972). In this con®guration, the role of lubrication forces

 gr ˆ 1

U Um



a

;

…39†

where Um  0:68 is the maximum volume concentration, and a ˆ 1:82. As shown in the inset of Fig. 4, our result is in agreement with Krieger-Dougherty equation.

4.3 Trajectories in shear flows Owing to the low computational cost of the method presented here, it can be used to simulate the time-dependent motion of large system of particles. Therefore, a numerical algorithm for tracking three-dimensional motion of particles in Stokes ¯ow has been employed (Tullock, 1993). Using the particles' velocities and an Euler kinematic scheme, the particles' trajectories and orientations will be determined at the end of each time step.

Fig. 5. Time-dependent motion of particles in shear ¯ow. The particles in simple cubic array have been traced for 30 iterations

395

396

Figure 5 shows the trajectories of 125 particles, initially forming a cubic array with a concentration of 20%, which are induced in a shear ¯ow and traced for 30 iterations. It is clear that a cubic lattice of spheres in a uniform shear ¯ow will be distorted by the ¯ow. As stated by Nunan and Keller (1984), the spheres will remain on a lattice, but the lattice will not remain cubic. The result obtained here is in agreement with the foregoing statement. The total CPU time, time spent by slaves solving Eq. (35) plus the time spent by master on the serial task, is approximately 12 h. This is obtained when two processors have been operating, with average slave CPU time of about 85% of total time.

5 Final remarks The CDLBEM is an ef®cient method for solving the stresslet or velocity on the surfaces of a large number of particles migrating in a viscous ¯uid associated with an ambient ¯ow. Indeed, the method is an ideal candidate to implement in a parallel computing environment, in a distributed system; and large-scale simulations ought to be performed in such environments. To accurately determine particle interactions in a concentrated suspension, the short-range lubrication forces have been evaluated in this work, and incorporated in the CDLBEM, resulting in a good agreement between what have obtained here and the exact solutions or experimental values for the intrinsic effective viscosity of suspension. References

Abbott JR, Tetlow N, Graham AL, Altobelli SA, Fukushima E, Mondy LA, Stephens TS (1991) Experimental observations of particle migration in concentrated suspensions: Couette ¯ow. J. Rheol. 5: 773±795 Arp PA, Mason SG (1977) The kinetics of ¯owing dispersions, VIII, Doublets of rigid spheres (theoretical). J. Colloid Interface Sci. 61: 21±43 Beguelin A, Dongarra J, Geist A, Manchek R, Sunderam V (1995) Recent enhancements to PVM. Int. J. Supercomputer Appl. 9(2): 108±127 Brady JF, Bossis G (1988) Stokesian dynamics. Ann. Rev. Fluid Mech. 20: 111±157 Brady JF, Phillips RJ, Lester JC, Bossis G (1988) Dynamic simulation of hydrodynamically interacting suspensions. J. Fluid Mech. 195: 257±280 Brebbia CA, Telles JCF, Wrobel LC (1984) Boundary Element Techniques: Theory and Application in Engineering. Springer-Verlag, Berlin, Germany Brenner H, O'Neill ME (1972) On the Stokes resistance of multiparticle systems in a linear shear ®eld. Chem. Eng. Sci. 27: 1421±1439 Chapman BK (1990) Shear Induced Migration Phenomena in Concentrated Suspensions. PhD thesis, University of Notre Dame Chapman BK, Leighton DT (1991) Dynamic viscous resuspension. Int. J. Multiphase Flow 17: 469±483 Chow AW, Sinton SW, Iwamiya JH, Stephens TS (1994) Shearinduced migration in Couette and parallel-plate viscometers: NMR imaging and stress measurements. Phys. Fluid A 6: 2561±2576 Claeys ILAM, Brady JF (1993a) Suspensions of prolate spheroids in Stokes ¯ow, Part I: Dynamics of a ®nite number of particles in an unbounded ¯uid. J. Fluid Mech. 251: 411±442

Durlofsky L, Brady JF (1989) Dynamic simulation of bounded suspensions of hydrodynamically interacting particles. J. Fluid Mech. 200: 39±67 Durlofsky L, Brady JF, Bossis G (1987) Dynamic simulation of hydrodynamically interacting particles. J. Fluid Mech. 180: 21±49 Fan XJ, Phan-Thien N (1997) Completed double layer boundary element method for periodic suspension. ZAMP 48: 1±12 Geist A, Beguelin A, Dongarra J, Jiang W, Manchek R, Sunderam V (1994) PVM3 user's Guide and Reference Manual. Oak Ridge National Laboratories, Oak Ridge, TN Hassonjee Q, Ganatos P, Pfeffer R (1988) A strong interaction theory for the motion of arbitrary three-dimensional clusters of spherical particles at low Reynolds number. J. Fluid Mech. 197: 1±37 Hassonjee Q, Pfeffer R, Ganatos P (1992) Behavior of multiple spheres in shear and Poiseuille ¯ow ®elds at low Reynolds number. J. Multiphase Flow. 18(3): 353±370 Hockney R (1994) The communication challenge for mpp. Parallel Computing 20: 389±398 Ingber MS (1990) Dynamic simulation of the hydrodynamic interaction among immersed particles in Stokes ¯ow. Int. J. Num. Method. Fluids 10: 791±809 Ingber MS, Mondy LA (1993) Direct second kind boundary integral formulation for Stokes ¯ow problems. Comp. Mech. 11: 11±27 Ingber MS, Womble DE, Mondy LA (1995) Simulations of hydrodynamics among immersed particles in Stokes ¯ow using a massively parallel computer. Boundary Elements 17: 515± 522 Kamel MT, Tory EM (1989) Sedimentation of clusters of identical spheres, I. Comparison of methods for computing velocities. Powder Technol. 59: 227±248 Kim S, Karrila SJ (1991) Microhydrodynamics: Principles and Selected Applications. Butterworth-Heinemann, Boston, USA Krieger IM (1972) Rheology of monodispersed lattices. Adv. Coll. Int. Sci. 3: 111±136 Koh CJ, Hookham P, Leal LG (1994) An experimental investigation of concentrated suspension ¯ows. J. Fluid Mech. 266: 1±32 Kynch GJ (1959) The slow motion of two or more spheres through a viscous ¯uid. J. Fluid Mech. 5: 193±208 Ladd AJC (1988) Hydrodynamic interactions in a suspension of spherical particles. J. Chem Phys. 88: 5051±5063 Leighton DT, Acrivos A (1986) Viscous resuspension. Chem. Eng. Sci. 41: 1377±1384 Leighton DT, Acrivos A (1987b) The shear induced migration of particles in concentrated suspensions. J. Fluid Mech. 181: 415±439 Lin CJ, Lee KJ, Sather NF (1970) Slow motion of two spheres in a shear ®eld. J. Fluid Mech. 43: 35±47 Mazur P, Van Saarloos W (1982) Many-sphere hydrodynamic interactions and mobilities in a suspension. Physica 115A: 21±57 Mo G, Sangani AS (1994) A method for computing Stokes ¯ow interactions among spherical objects and its application to suspensions of drops and porous particles. Phys. Fluids 6(5): 1637±1652 Nasseri S, Phan-Thien N (1996) On the path and ef®ciency of two micromachines with rigid tails. Comp. Mech. 18(3): 192±199 Nasseri S, Phan-Thien N (1997a) Modelling micromachines with elastic parts in a viscous environment. Comp. Mech. 20(3): 242±246 Nasseri S, Phan-Thien N (1997b) Geometric optimization of a micromachine with spiral tail. Comp. Mech. 20(3): 267±271 Nasseri S, Phan-Thien N (1997c) Hydrodynamic interaction between two nearby swimming micromachines. Comp. Mech. 20(6): 551±559 Nunan KC, Keller JB (1984) Effective viscosity of a periodic suspension. J. Fluid Mech. 142: 269±287

O'Brien J (1979) A method for the calculation of the effective transport properties of suspensions of interacting particles. J. Fluid Mech. 91: 17±39 Pakdel P, Kim S (1991) Mobility and stresslet functions of particles with rough surfaces: A numerical study. J. Rheology 35: 797±823 Phan-Thien N, Fan XJ (1995) Traction based completed adjoint double layer boundary element method in elasticity. Com. Mech. 16: 360±367 Phan-Thien N, Kim S (1994) Microstructures in Elastic Media: Principles and Computational Methods. Oxford University Press, New York Phan-Thien N, Tran-Cong T, Graham AL (1991) Shear ¯ow of periodic arrays of particles clusters. J. Fluid Mech. 228: 275±293 Phan-Thien N, Tullock DL (1993) Completed double layer boundary element method in elasticity. J. Mech. Phys. Solids 41: 1067±1086 Phan-Thien N, Tullock DL (1994) Completed double layer boundary element method in elasticity and Stokes ¯ow: Distributed computing through PVM. Comp. Mech. 14: 370±383 Phan-Thien N, Tullock DL, Kim S (1992) Completed double layer in half-space: a boundary element method. Comp. Mech. 9: 121±135 Power H, Miranda G (1987) Second kind integral formulation of Stokes ¯ow past a particle of arbitrary shape. SIAM J. Appl. Math. 47: 689±698 Pozrikidis C (1992) Boundary Integral and Singularity Methods for Linearized Viscous Flow, CUP, Cambridge Rallison JM (1981) A numerical study of the deformation and burst of a viscous drop in general shear ¯ows. J. Fluid Mech. 109: 465±482 Rallison JM (1984) The deformations of small viscous drops and bubbles in shear ¯ows. Ann. Rev. Fluid Mech. 16: 45±66 Rider PF, Phan-Thien N (1997) A numerical viscometer via distributed computing and the completed double layer boundary element method. Eng. Anal. Boun. Elem. 19: 57±65

Sadd Y, Schultz MH (1983) GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems. Research Report, YALEU/DCS/RR-254, Department of Computer Science, Yale University Sangani AS, Mo G (1994) Inclusion of lubrication forces in dynamic simulations. Phys. Fluids 6(5): 1653±1662 Scha¯inger U, Acrivos A, Zhang K (1990) Viscous resuspension of a sediment within a laminar and strati®ed ¯ow. Int. J. Multiphase Flow 16: 567±578 Seeling C, Phan-Thien N (1994) Completed double layer boundary element algorithm in many-body problems for a multi-processor: An implementation on CM5. Comp. Mech. 15: 31±44 Shakib F, Hughes TJR, Johan Z (1989) A multi-element group preconditioned GMRES algorithm for nonsymmetric systems arising in ®nite element analysis. Comp. Meth. App. Mech. Eng. 75: 415±456 Stone HA (1994) Dynamics of drop deformation and breakup in viscous ¯uids. Ann. Rev. Fluid Mech. 26: 65±102 Thomas DG (1965) Transport characteristics of suspensions: VIII, A note on the viscosity of Newtonian suspensions of uniform spherical particles. J. Colloid Sci. 20: 267±277 Tran-Cong T, Phan-Thien N (1989) Stokes problems of multiparticle systems: A numerical method for arbitrary ¯ows. Phys. Fluids A 1(3): 453±461 Tullock DL (1993) New Developments and Applications of the Boundary Element Method for some Problems in Elasticity and Viscous Flow. PhD thesis, Univ. of Sydney, Australia Tullock DL, Phan-Thien N (1992) A completed double layer boundary element method with domain decomposition. Proc. BETECH 92, Albuquerque, New Mexico, June 1992 Youngren GK, Acrivos A (1975) Stokes ¯ow past a particle of arbitrary shape: A numerical method of solution. J. Fluid Mech. 69: 377±403

397