Roles of Particle-Wall and Particle-Particle Interactions in Highly Confined Suspensions of Spherical Particles Being Sheared at Low Reynolds Numbers

The roles of particle-wall and particle-particle interactions are examined for suspensions of spherical particles in a viscous ﬂuid being conﬁned and sheared at low Reynolds numbers by two parallel walls moving with equal but opposite velocities. Both particle-wall and particle-particle interactions are shown to decrease the rotational velocity of the spheres, so that in the limit of vanishingly small gaps between the spheres and the walls, the spheres acquire a rotational slip relative to the walls. The presence of the walls also increases the particle stresslet and, therefore, the total viscous dissipation. In the limit of vanishingly small gaps, the increased viscous dissipation in the gaps between pairs of spheres aligned in the ﬂow direction is largely compensated by the reduction in the dissipation in the gaps between the spheres and the walls due to a reduction in the rotational velocity of the spheres. As a result, the effect of short-range particle interactions on the stresslet is generally insigniﬁcant. On the other hand, the channel-scale particle interactions in the shear ﬂow induced by the moving walls decrease the particle stresslet, primarily because the fraction of pairs of spheres that are aligned parallel to the ﬂow (the presence of which in a shear ﬂow reduces the stresslet) is relatively higher than in unbounded suspensions. Expressions are also derived for the total stress in dilute random suspensions that account for both the particle-wall and the channel-scale particle-particle interactions in determining the rotational velocities and stresses. The latter are shown to be consistent with recent numerical [Y. Davit and P. Europhys . Lett. 83 , 64001 (2008)] and Peyla and C. . Lett. , (2011)] ﬁndings according to which, for a range of sphere radius to gap width ratios, the effect of particle-particle interactions is to decrease the total dissipation. V C 2011 American Institute of Physics .


I. INTRODUCTION
Flows of suspensions of neutrally buoyant particles through channels of width comparable to the particle dimension are of considerable interest because of their occurrence in many experimental, biological, and technological systems including blood flow in capillaries, porous media, and microfluidic devices. Previous studies have generally focused on the hydrodynamic mobility of either isolated or pairs of spheres and on the pressure drop driven suspension flow through a channel under conditions of vanishingly small Reynolds numbers, 1-7 with relatively little attention being given to the case of a sphere or of pairs of spheres in another type of flow, namely, the shear flow induced by two parallel plane walls moving with equal and opposite velocities. Such a flow occurs, for example, in a parallel plate viscometer widely used for rheological measurements.
Ganatos et al. 2 determined the force on a sphere of finite size held fixed in a shear flow between two plane walls, but did not address the effect of the walls on the rheology of suspensions. Recently, Davit and Peyla 8 and Swan and Brady 9,10 determined the relative change in the viscous dissi-pation in such systems as a function of the particle volume fraction / and a, with a being the radius of the sphere divided by the half-width of the channel. Davit and Peyla used their numerical results for various / and a to estimate the Oð/Þ and Oð/ 2 Þ terms in the expansion of the relative dissipation in powers of / and found that although, as expected, the Oð/Þ coefficient increased monotonically as a increased from 0 to 1, the Oð/ 2 Þ coefficient decreased. For example, at a ¼ 0:5, their estimated Oð/Þ and Oð/ 2 Þ coefficients were found to equal, respectively, about 5 and À 5 compared with the well-known values for the bulk suspensions ða ¼ 0Þ of 5=2 and about 5.0. 11 More recently, Pelya and Verdier 12 have determined experimentally the viscosity of suspensions at various / and a and confirmed the finding by Davit and Peyla that the Oð/ 2 Þ coefficient decreases and becomes negative as a is increased. The existence of a local maximum relative dissipation at / % 0:45 À 0:50 has also been found recently by Yeo and Maxey, 13 which they attributed to the particle ordering near the walls.
In order to gain more insight into how the extra dissipation of a suspension is affected by a and by the particle-particle interactions, we examine, using semi-analytical tools, a few relatively simple problems. Although, it is of course possible today with modern day computers and the powerful algorithms that have been recently developed (e.g., Refs. 4,5,14,and 15) to perform direct numerical simulations that fully account for multiparticle interactions in the presence of walls, we believe that the analysis of the simple problems chosen here provide valuable insight into the role of particlewall and particle-particle interactions in determining the rheology of highly confined sheared suspensions in addition to providing estimates of the Oð/Þ and Oð/ 2 Þ coefficients of the relative dissipation.
The first problem concerns the motion of a single sphere freely suspended in a simple shear flow bounded by two parallel plane walls a distance apart equal to 2, i.e., twice the half-width of the channel, which has been set equal to unity with no loss of generality and moving with equal and opposite velocities of unit magnitude. This case was examined previously by Ganatos et al., 2 who determined the torques on a sphere held fixed in a simple shear flow and on a rotating sphere with stationary walls but did not evaluate the extra dissipation as given by the stresslet 13 induced by the presence of the sphere. Even in this simple case, the analysis to be presented leads to a rather interesting counter-intuitive result concerning the way in which the rotational velocity of the sphere changes with increasing confinement, in addition to yielding an expression for the Oð/Þ coefficient of the extra dissipation over the full range of a. For example, consider the special case when the sphere is placed at the center of the channel where its translational velocity is zero owing to symmetry. One expects that, in the limit of vanishingly small gaps between the sphere and the walls, the velocities of the points on the surface of the sphere closest to the walls will be equal to that of the walls and that the rotational velocity X will, therefore, approach unity, i.e., twice its value at a ¼ 0. In other words, one expects the walls to enhance the value of X. But, as will be seen shortly, exactly the opposite is the case in that the rotational motion is hindered by the walls with its value approaching 1=4 in the limit of vanishingly small gap. This result arises from a known, but not well appreciated, consequence of the lubrication analysis of the forces in small gaps, according to which, the torque on a sphere moving with translational velocity V near a stationary plane wall is roughly one-fourth of the torque on the same sphere rotating with angular velocity V=a. We also examine the case when the sphere is located off-center and determine, via matched asymptotic expansions, the average rotational velocity and stresslet in dilute suspensions of randomly distributed spheres. The predictions of this analysis are shown to be in good agreement with the results obtained numerically. The result for the average stresslet also corrects one given recently by Swan and Brady, 9 as explained later (cf., see the discussion following Eq. (11)).
We next examine the role of particle-particle interactions by analyzing separately the effect of short-and longrange interactions. We first consider a pair or a row of an infinite number of equi-spaced spheres arranged parallel to the flow direction. We show that the effect of particle interactions is to further reduce the rotational mobility X of each sphere and, hence, to enhance the rotational slip between each sphere and the walls. This additional reduction in X, however, decreases the viscous dissipation in the gap between the walls and the sphere and offsets the viscous dissipation arising from the hydrodynamic interaction between the spheres, thereby yielding an overall stress per sphere which is approximately the same as that for a single sphere between the two walls. Thus, short-range, lubrication-dominated particle interactions only modestly affect the stress induced in the suspension. We then determine the effect of long-range particle interactions on the relative viscous dissipation in dilute random suspensions and show that their contribution to the Oð/ 2 ) coefficient in the expansion of the relative viscous dissipation decreases as a is increased from 0 to about 0.6 and, in fact, becomes negative for a > 0:3 in agreement with the findings by Davit and Peyla. 8 This negative influence of the pair interactions on the overall stress in random suspensions is shown to arise primarily from the relative increase in the number of pairs of spheres that are aligned parallel to the flow compared to those in unbounded suspensions.

II. A SINGLE SPHERE AT THE CHANNEL CENTER
Let us first consider the simple case of a neutrally buoyant spherical particle freely suspended in a viscous fluid, the viscosity l of which is set equal to unity without loss of generality ðl ¼ 1Þ. The fluid is being sheared by two plane walls separated by a distance equal to 2 and moving with velocities 61, thereby generating a simple shear flow with shear rate c ¼ 1 away from the sphere. Furthermore, the sphere is placed midway between the two walls so that, because of symmetry, its translational velocity is zero. We shall determine the rotational velocity and the stresslet induced by the sphere as functions of a, its radius divided by the half-width of the channel, by modifying the method recently developed by Ozarkar and Sangani, 17 who determined the mobility and resistivity of a sphere placed in a thin film of liquid bounded on one by a rigid wall and a stress-free planar interface on the other, by expressing the flow induced by the sphere in terms of Lamb's multipoles located at the center of the sphere plus an image system that ensures that the boundary conditions at the wall and at the gas-liquid interface are satisfied. The accuracy of the method in the present case depends on two parameters for a given value of a; the order N s of Lamb's multipoles used to represent the flow induced by the particle, and the number of reflections N r of these multipoles on the two sides of the film to account for the no-slip boundary conditions on the plane walls. 17 We used up to N s ¼ 20 and N r ¼ 32 which gave accurate results, say within 5%, for a up to about 0.9.

A. Results for the rotational velocity
The computed results for the angular velocity X of a freely suspended particle are shown in Fig. 1, where, as expected, in the limit a ! 0, the rotational velocity is seen to approach 1 = 2 , i.e., the vorticity of the imposed shear flow. On the other hand, the results for larger a are surprising since one expects intuitively that the velocity of the points on the surface of the sphere closest to the walls should approach that of the walls due to the no-slip condition, i.e., one anticipates that X ! 1 as a ! 1. Instead, the effect of the walls is to hinder the rotation of the particle, and in fact, as shown below, X ! 1=4 as a ! 1. This hindrance effect of the walls is better understood by examining separately the two limiting behaviors of X for small and for large a, respectively.
We start by recalling that, as e 1 À a approaches zero, the flow is dominated by the lubrication stresses in the narrow gap between the walls and the sphere. But since, due to symmetry, the net force on the particle is zero, we need to consider only the lubrication contributions to the torque on the sphere. The torque on a sphere rotating with an angular velocity X in a quiescent liquid film bounded by two walls is given by  17 Similarly, the torque on the sphere held fixed in a shear flow with c ¼ 1 is given by e ln e À1 þ 0:26e; where the coefficients multiplying the O(ln e À1 ) and O(e ln e À1 ) terms were obtained from the analysis of torque on a fixed sphere near a moving wall by O'Neill and Stewartson. 20 Figure 2 shows a comparison between the numerical results for R LR and R LS and the expressions given by Eqs. (1) and (2). The coefficients R LR and R LS determined here are also in very good agreement with those obtained previously by Ganatos et al. 2 The angular velocity of a freely rotating sphere obtained by requiring that the net torque on the sphere must be equal to zero is, therefore, given by e ln e À1 þ 0:26e 4 5 ln e À1 À 0:47 þ 132 125 e ln e À1 þ 1:18e from which it is evident that, since the coefficient of the leading Oðln eÞ term in R LR for a rotating sphere is four times that of a sphere fixed in the shear flow, we have that X ! 1=4 as e ! 0. In other words, the hindrance due to walls arises because the magnitude of the torque on a sphere rotating with angular velocity X is, to leading order, four times greater than that due to a sphere translating with velocity Xa near a stationary wall. It is interesting to note that the opposite is the case for the force on a sphere, in that the magnitude of the force on a sphere rotating with angular velocity V=a near a stationary wall is, to leading order, four times smaller than that on a sphere translating with velocity V. Now, let us consider the opposite limit, i.e., a ! 0, for which the rotational velocity can be determined from the velocity gradient near the center of the channel induced by the images of a point stresslet which, in this limit, is determined from Eq. (33) in Ozarkar and Sangani 17 with n ¼ 2; m ¼ 1; k ¼ 0. One must, however, first determine numerically the coefficients / 0;r 21 (related to the rate of strain), P 0;r 21 , and T 1;r 11 (related to the vorticity) in terms of the stresslet coefficient P 0 21 , use Eq.

083302-3
Roles of particle-wall and particle-particle interactions Phys. Fluids 23, 083302 (2011) then substitute the result into Eq. (35) in Ref. 17. This analysis yields with c 1 ¼ 0:213, c 2 ¼ 0:221, c 3 ¼ 1:78, and c 4 ¼ 1:86. The leading order correction in the above equation arises from the coefficient c 1 which is related to the vorticity induced by the wall images of a point stresslet placed at the center of the channel. The coefficient c 3 is related to the rate of strain induced by these images and, as we shall see later, this coefficient is important in determining the leading order correction to the stresslet (cf. Eq. (6)). The coefficients c 2 and c 4 arise from two separate contributions of equal magnitude. The first is related to the Laplacians of the vorticity and the rate of strain induced by the wall images, while the second is due to the leading order contributions to the vorticity and strain rate from the source quadrupole (/ 0 21 ¼ ða 2 =10ÞP 0 21 in Ref. 17). As seen in Fig. 1, the computed results are in excellent agreement with the two limiting cases described above within their respected ranges of applicability.
Several studies have investigated the behavior of a particle near a single wall or between two walls. Goldman et al. 21 and Chaoui and Feuillebois 22 considered a freely suspended sphere placed in a simple shear flow near a plane wall and found that, in the limit of vanishingly small gap between the particle and the wall, the ratio Xa=V of the rotational to the translational velocity V of the sphere approaches 0.4218=0.7431 ¼ 0.5676. On the other hand, according to Ozarkar and Sangani, 17 for a freely suspended sphere in a pressure gradient driven flow in a film bounded by a plane wall on one side and a stress-free plane interface on the other, the ratio Xa=V approaches 0.25=0.775 ¼ 0.290. In both cases, the determination of this ratio required that the flow field be first obtained in the outer region, i.e., away from the gap between the wall and the sphere, which can only be achieved from a detailed numerical analysis. This is in contrast to the present case where the value of the leading coefficient, 1=4, follows solely from the lubrication analysis. Finally, Staben et al. 6 examined the case of a sphere freely suspended at the center of two plane walls in a pressuredriven flow, where the rotational velocity of the particle is zero because of symmetry but the translational velocity is not. In all three cases, however, the velocity of the sphere becomes vanishingly small as the gap between the sphere and the wall decreases to zero, so that there is no slip between the sphere and the wall in contrast to the problem examined here where the corresponding slip velocity between the particle surface and the walls approaches a nonzero value.
To assess the accuracy of the Fluid Particle Dynamics (FPD) method proposed by Tanaka and Araki 23 and further developed by Peyla 24 and Davit and Peyla, 8 calculations were also carried out for different values of a=d, d being the grid size. In the study cited earlier (Ref. 8), Davit and Peyla used a=d ¼ 3 to simulate confined suspensions at various particle volume fractions and sphere radii to wall gap ratios. It was found that much higher values of a=d, equal to about 12 or greater, were required to obtain accurate results in agreement with our method, especially for large a.
Figures 3(a) and 3(c) show the velocity profiles obtained using the FPD method for two different a. The disturbance flows generated by the presence of the sphere obtained by subtracting the imposed shear flow are shown in Figs. 3(b) and 3(d). We see clearly the formation of vortices ahead and behind the sphere at the higher value of a ( Fig. 3(b)). (The upper wall moves from right to left so that the vorticity of the imposed shear is counter-clockwise.) These vortices act counter to the imposed vorticity, which explains why the presence of the walls leads to the reduction in the angular velocity of the sphere.
The velocity fields shown here are similar to those shown by Bikard et al., 25 who used a finite element method for the purpose of computing the velocity and stress distributions. These investigators also determined the rotational speed of the sphere for a < 0:9. Their results, presented in terms of the time required for a sphere to complete one rotation, given by curve 2 in their Figure 10 are also in reasonable agreement with ours, although their expression X ¼ 1=½2ð1 À 2:212a þ 0:64a 2 Þ (cf. their Eq. (9)) is inconsistent with their own results in Figure 10 and predicts an increase in the rotational velocity with an increase in a for small particle radii. The rotational velocity of the sphere was also computed by D'Avino et al. 26 for a 0:83 for both Newtonian as well as viscoelastic fluids using a method of finite elements. Once again, their results for Newtonian fluids are in agreement with ours.

B. Induced stresslet
Highly confined suspensions also lead to large stresses. The stress resulting from the presence of a rigid spherical particle can be determined from the induced stresslet S ij which is related to the traceless, symmetric part of the first moment of the traction at the particle surface 13 where r i is the position vector of a point on the surface of the sphere relative to its center, f j is the traction (force per unit area) exerted by the fluid, and @D is the surface of the particle. For the special case of an unbounded simple shear flow given by u 1 i ¼ x 2 d i1 , the only nonzero components of the stresslet for a force-free sphere are S 12 and S 21 , both of which are equal to 10pa 3 =3. 13 Now, in the presence of the walls, we let these non-zero components of the stresslet be denoted by S Ã times this infinite dilution value, so that S Ã depends only on a with S Ã ð0Þ ¼ 1 when particle-particle interactions are negligible. Consequently, the dissipation rate in the channel divided by that in the absence of the spherical particles (sometimes referred to as the relative effective viscosity of the suspension) equals 1 þ ð5/=2ÞS Ã , with / being the volume fraction of the spheres, which reduces to Einstein's well known expression for the relative viscosity when a ¼ 0. Note that here, and henceforth, the volume fraction is defined  Sangani, Acrivos, and Peyla Phys. Fluids 23, 083302 (2011) as the total volume occupied by the spheres divided by the channel volume. Figure 4 shows the numerical results for S Ã , as obtained using the method of Ozarkar and Sangani, 17 as well as their comparison with Eqs. (6) and (9) given below for a small and a ! 1, respectively.
The small a analysis that led to Eq. (4) also yields where c 3 and c 4 are as in Eq. (4). In the opposite limit, a ! 1, the stresslet is determined, once again, by examining separately two problems: (i) the sphere held fixed in the shear flow generated by the two moving walls and (ii) in which the sphere rotates in the presence of stationary walls. For the first case, the stresslet is given by 10p a 3 S s =3, with 10 log e À1 À 0:9 þ 198 125 e log e À1 þ 0:24e !
and for the second case by 10pXa 3 S r =3, where, as can be shown by applying the reciprocal theorem, S r is given by with R LR and R LS being the coefficients appearing in the expressions for the torques introduced earlier (cf. Eqs. (1) and (2)). The coefficients of the log e À1 and e log e À1 terms in Eq. (7) were taken from the expression given by Jeffrey 27 for the stresslet induced on a fixed sphere in the presence of a moving wall. Combining the results for the above two problems, we see that the stresslet induced by a freely suspended sphere between the two moving walls is given by 10pa 3 S Ã =3 with  6) and (9).

083302-5
Roles of particle-wall and particle-particle interactions Phys. ð9Þ with X being given by Eq. (3). As seen in Fig. 4, the numerical results for S Ã compare very well with the limiting expressions given by Eqs. (6) and (9).
Note that, since X ! 1=4 in the limit of vanishingly small e, S Ã diverges as ð9=5Þ log e À1 . Note also that, in this limit, the contribution from the wall motion is 14 times that due to the particle rotation. Finally, it is also interesting to point out that, although the contributions to the total torque from the wall shear and the particle rotation are of opposite signs and cancel each other for a freely suspended sphere, both contribute positively to the total stress which becomes smaller as X decreases due to the presence of the wall. In particular, the total stress would have been greater if X had been equal to unity-a case one might have thought would contribute the least to the lubrication stresses given the absence of slip between the velocities of the walls and the points on the sphere closest to the walls.

III. DILUTE RANDOM SUSPENSIONS
When the center of the sphere is no longer at the midplane of the channel, its stresslet, S Ã , as well as its translational and rotational velocities will be functions of its position in addition to a. In order to evaluate the extra dissipation of a very dilute random (hard sphere) suspension, we consider, therefore, the problem of determining the average values of S Ã for the special case when the probability density for the sphere position is uniform throughout the channel except, of course, for distances from the walls that are less than one radius. We shall examine the limiting cases of small and large spheres separately and compare the analytical results with those obtained numerically for selected values of a.
A. The small sphere limit: a ( 1 The correction to the stresslet of a sphere, relative to its value when a ! 0, is either Oð1Þ or Oða 3 Þ depending on whether its center is OðaÞ or Oð1Þ from one of the walls. Therefore, for a suspension having a uniform probability of finding a sphere anywhere in the channel, the coefficients of the leading OðaÞ corrections to the average stresslet and rotational velocity can be determined from the detailed calculations involving the interaction of the sphere with a single wall, and the method of matched asymptotic expansions may then be used to evaluate the next few corrections as shown in Appendix A to yield hXi ¼ 1=2 À 0:078a À 0:078a 2 þ 0:109a 3 À 0:312a 4 þ Oða 5 Þ: Swan and Brady 9 have also obtained the coefficient of the OðaÞ term for the average stresslet, but unfortunately their estimate 0.12 for the value of this coefficient is incorrect owing to two errors in their calculations. First, in determining the stresslet from their formula for a linear shear flow, these authors substituted an incorrect strain rate tensor-they for the simple shear flow. Second, to determine the average stresslet, these authors integrated the Oða 3 Þ outer region approximation (cf. Appendix A) for the stresslet all the way from the center of the channel to a distance y ¼ a from the wall. This outer region approximation, however, is not valid for distances from the wall that are comparable to a. The first error can be corrected by replacing g 3 in their Eq. (64) by h ES 3 , whereas the second one requires that an inner region approximation be constructed as we have done in Appendix A. Note that the Oða n Þ outer region approximations are singular and behave as (a=y) n as y ! 0 for n ! 3, so that, in principle, all higher-order outer region approximations must be determined to obtain a correct estimate of the coefficient of the OðaÞ correction to hXi from the outer region analysis alone. B. The large sphere limit: 1 À a ( 1 In this limit, the expressions for hS Ã i and hXi are similar to Eqs. (3) and (9), respectively, provided they are modified to account for the fact that the center of the sphere is, generally, no longer on the mid-plane. Let then the gap between the sphere and one of the walls be eð1 À aÞ with À1 < a < 1 and e ¼ 1 À a ( 1. Now the translational velocity of the sphere is nonzero and we must use both the force and torque balances on the sphere to determine its translational and rotational velocities. For this purpose, we need to first determine the force F À6paf t U on a sphere translating with velocity U along the center-plane of a channel with stationary walls, a problem already examined by Ganatos et al., 2 Staben et al., 6 and Bhattacharya et al. 4 Results obtained with N r ¼ 20 using the method of Ozarkar and Sangani 17 are in very good agreement with those in the literature and can be fitted quite well by means of the following expression for a ranging from 0.3 to 0.9: e ln e À1 À 0:9e: The translational and rotational velocities of the off-center sphere are then given by with Here, X, R LS , and R LR correspond to the results for the sphere placed at the mid-plane, while those denoted with the subscript o refer to the corresponding 083302-6 Sangani, Acrivos, and Peyla Phys. Fluids 23, 083302 (2011) quantities when the sphere is placed off-center. The nondimensional stresslet may be subsequently estimated using where S Ã is given by Eq. (9). To leading order, f t ¼ ð16=15Þ ln e À1 and 4 À Xa ¼ Oð1= ln eÞ, and therefore, the off-center translational velocity given by Eq. (13) is Oðln eÞ À2 . Since the center of the sphere is shifted from the channel center-plane by ae, the coefficients of the OðeÞ terms in R LS , R LR , and f t etc., which were determined for the sphere placed at the center of the channel, will now be the functions of a. Nevertheless, neglecting this dependence on a will lead to relative errors of only Oðe= ln eÞ in the estimates of the angular velocity and of the stresslet. For the special case of randomly placed spheres with uniform probability distribution for a, we obtained the average values by numerically integrating Eqs. (13) and (14) with a varying from zero to unity. For a equal to 0.85, 0.9, 0.95, and 0.99, we found that the ratio of hXi to its center line value equals, respectively, 0.92, 0.93, 0.95, and 0.97, indicating that the average rotational velocity of large spheres randomly placed in off-center positions is slightly smaller than if the spheres were placed at the center of the channel. The values for the ratios of the corresponding stresslets are slightly higher and are given by, respectively, 1.1, 1.09, 1.08, and 1.06. Note that this differs from what was found earlier in the case of the sphere centered along the mid-plane, where a decrease in X was followed by a decrease in S Ã .

C. Comparison with numerical results
For selected values of a, the stresslet and rotational velocities of the sphere were also evaluated using the method of Ozarkar and Sangani 17 at ten different positions of the sphere and the corresponding average values were thereby computed. As shown in Figs. 5 and 6, the numerical results are in good agreement with those obtained analytically for small and for large spheres. Also shown are the results for a sphere placed at the channel center. The channel-centered spheres have smaller stresslets and greater rotational velocities, although the differences are minor for a greater than say 0.7.

IV. PARTICLE-PARTICLE INTERACTIONS
We shall now investigate the effect of particle-particle interactions by first analyzing, in Sec. IV A, the lubricationdominated, short-range interactions of relatively large spheres for the special case when the spheres are placed on the mid-channel plane, to be followed, in Sec. IV B, by an analysis for small spheres, separated by a distance comparable to the channel width but, otherwise, arbitrarily placed within the channel. The results of such an analysis will then be compared in Sec. IV C with those obtained by direct numerical computations for a few selected cases of mediumsized spheres. Finally, in Sec. IV D, we shall examine the more general case of random dilute suspensions and thereby arrive at an explanation of why the Oð/ 2 Þ coefficient in the expansion of the relative viscous dissipation in powers of / may be expected to decrease as a is increased.

A. Lubrication-dominated, short-range interactions
Let us first consider two freely suspended equi-sized spheres, placed at the mid-plane with their center to center distance equal to R, the channel width being equal to 2 as before. To assess the importance of the lubrication-dominated, short-range interactions, we examine the limit a ! 1 and R ! 2. Since, when the pair of spheres is aligned in the vorticity direction, there are no lubrication effects in the gap between the spheres, we focus our attention mainly to the case when the pair of spheres is aligned in the flow direction. Now, the lubrication effects in the gap between the two

083302-7
Roles of particle-wall and particle-particle interactions Phys. Fluids 23, 083302 (2011) spheres cause an additional torque on each sphere both of which are rotating with the same angular velocity X. This extra torque equals 8pa 3 XL p with L p ! ð1=4Þ ln e À1 p as e p R=2 À a ! 0. 19 Therefore, to the leading order, the torque balance on a sphere gives (cf. Eq. (3)) that X ! ð1=5Þ ln e À1 ð1=4Þ ln e À1 p þ ð4=5Þ ln e À1 with e 1 À a representing the gap between a sphere and the walls as before. For the special case, when the gap between the spheres equals that between the sphere and the walls, i.e., e p ¼ e, the coefficient of the torque contribution due to particle-particle interaction is 1=4 compared to 4=5 from the moving walls and X ! 4=21. Comparing this result with X ! 1=4 for the single particle case, we see that the short-range, lubrication-dominated, particle-particle interaction reduces the rotational velocity by about 24%. On the other hand, the corresponding effect of the particle-particle interactions to the total stresslet is negligible because, by adding the relevant contribution from the rotating particles to the expression given by Eq. (7), one finds that with X given by Eq. (15). For the special case e p ¼ e, this yields S Ã ! c ln e À1 with c ¼ 64=35 compared to 9=5 for the single sphere case, a net increase of only 1.6%. Thus, we see that the increase in the stresslet due to the particle-particle interaction is offset by the corresponding reduction in the contribution from the particle-wall gaps as a result of the reduced rotational velocity.
In the extreme case of nearly touching spheres, i.e., e p ( e, the rotational velocity of each sphere vanishes as X ! ð4=5Þ ln e À1 = ln e À1 p and S Ã ! ð48=25Þ ln e À1 , a modest increase of about 12% from the single sphere result. The lubrication forces arising from the rotation of the spheres also induce a translational motion of each sphere towards the opposite walls which, by means of a force balance in the direction normal to the walls, is found to be Oðe ln eÞ, too small to influence the leading order analysis for the rotational velocity and the stresslet given above by Eqs. (15) and (16).
Next, let us consider the case of a periodic row of spheres aligned in the flow direction with the centers of the spheres on the mid-plane and with the center-to-center distance between two adjacent spheres fixed at R ¼ 2. The lubrication analysis for this case yields X ! 2=13 and S Ã ! ð24=13Þ ln e À1 as a ! 1, and hence, here, the particle interactions reduce the rotational velocity very significantly (38%) while increasing the stresslet, although by less than 5%. The increase in the viscous dissipation in the gap between the particles is again mostly compensated by the decrease in the dissipation in the gaps between the particles and the walls.

B. Long-range interactions
Now, let us consider the effect of long-range interactions by examining the opposite limit a ( R ¼ Oð1Þ, R being, once again, the center-to-center distance between the two spheres, although now we do not restrict their centers to lie on the mid-plane. As in the case of the single sphere analysis which yielded Eqs. (4) and (6), the disturbance flow created by each sphere can be approximated in this limit by that due to a point stresslet of magnitude 10pa 3 =3. The rotational velocity and the normalized stresslet of each sphere can then be shown to be given by with c 1 and c 3 representing, as before, the contributions from the walls, and X p and S p the contributions from the other sphere. For the special case of spheres located at the midplane, c 1 ¼ 0:218 and c 3 ¼ 1:78. For more general case, they are functions of the x 2 -coordinate of the center of the sphere and their values can be determined from the outer region analysis presented in Appendix A. X P and S p , on the other hand, are functions of the positions of the centers of the two spheres, which will be denoted, henceforth, by x and X (cf. Fig. 7). Let us first consider the limiting case a ( R ( 1 for which the spheres are separated by a distance that is large compared to their radius but small compared to the channel width. Furthermore, let the centers of the spheres be also at large distances from the either wall. This limit corresponds to the well-examined case of two widely separated freely suspended spheres in an unbounded shear flow for which 11 whereR i are the components of the unit vector along the line joining the centers of the two spheres ðrecall that R i X i À x i Þ. For the special case of the pairs of spheres aligned in the flow direction, we have thatR i ¼ d i1 and, therefore, S un p ¼ À5=ð2R 3 Þ and X un p ¼ 5=ð4R 3 Þ. Thus, the  stresslet of a sphere is decreased due to the presence of a distant second sphere in the flow direction. As in the case of the short-ranged, lubrication dominated interactions, the presence of another sphere in the vorticity direction ðR i ¼ Rd i3 Þ does not affect the stresslet or the rotational velocity of a sphere. Finally, the maximum contribution to S p (for fixed R) results from the presence of a second sphere along one of the principal axes of the extensional flow, i.e., at 645 to the flow direction in the plane of shear jRj 1 ¼ jRj 2 ¼ R= ffiffi ffi 2 p Á À , where it is positive and equals 15=ð4R 3 Þ.
Next, let us consider the other limiting case corresponding to R ) 1 where the presence of the walls plays a significant role as noted by Liron and Mochon, 28 who obtained an expression for the velocity induced by a point force in a fluid bounded by two plane walls. These investigators showed that the velocity due to a point force is screened by the walls, so that the channel-width averaged velocity decays as R À2 for large R compared to the slower decay, proportional to R À1 , in an unbounded Stokes flow. To obtain an expression for S p for a sphere centered at X due to a stresslet at the center x of another sphere, one must determine the gradients of the velocity induced by a point force at x with respect to both X and x. This results in an expression for S p which, as shown in Appendix B, decays as R À2 rather than as R À3 (cf. Eq. (18)) as is the case of an unbounded Stokes flow interaction of two freely suspended spheres. In other words, the interactions of freely suspended particles in wall-bounded sheared flows are longer-ranged than their counterparts in unbounded flows. This is a consequence of the fact that the velocity field varies more rapidly in the direction normal to the walls owing to the no-slip condition at the walls than in the direction parallel to them. Since the walls are separated by an O(1) distance, both the velocity and its gradient normal to the walls are of the same order of magnitude. As shown in Appendix B, the detailed analysis gives where q 2 ¼ R 2 1 þ R 2 3 , ; h ¼ cos À1 ðR 1 =qÞ, and R 2 ¼ q 2 þðX 2 À x 2 Þ 2 . The exponential decay constant b equals p=2 when neither of the spheres is located on the mid-plane, and equals 2.1062 otherwise. It is interesting to note that S p and X p for both spheres are equal in the two extreme limits R ( 1 and R ) 1. For intermediate separations, this is not the case. Figures 8 and 9 show S p and X p as functions of q for two special cases of pairs of spheres aligned in the direction of the flow (h ¼ 0; x 2 ¼ X 2 ¼ 0, 1=2). The symbols represent the results obtained by the method described by Ozarakar and Sangani 17 with N r ¼ 24 (the number of reflections of the point stresslet at x on the either side of the channel), while the dashed curves represent the small q approximations as given by Eqs. (18) and (19). (Note that q ¼ R for the special case corresponding to X 2 ¼ x 2 .) The computed results for X p and S p are seen to agree well with the above asymptotes. We also see that, in accordance with Eq. (20), both X p and S p decrease more rapidly with increasing q for the case in which the spheres are centered at the mid-plane, compared to that in which they are off-centered. Also, when the spheres are centered at the mid-plane, X p changes its sign at q ffi 1:95; from positive values at smaller q to small negative values for larger with the minimum value of about À 0.017 attained at q ffi 2:3: (Only the positive values of X p for q < 1:95 are shown in Fig. 9 because of the log scale used for plotting.) A similar behavior was also observed for S p for the mid-planecentered spheres with the sign change observed at q ffi 3:8.
Equation (20) suggests that both S p and X p must change sign at sufficiently large q if the sign of either X 2 or x 2 is changed. This prediction was verified numerically for q ¼ 4:1 for which S p equals À0.0454 for x 2 ¼ X 2 ¼ 1=2 and 0.0453 for x 2 ¼ ÀX 2 ¼ 1=2. Note that this is not the case for FIG. 8. Reduction in the normalized stresslet of a sphere (divided by a 3 ) due to the presence of a second sphere at a distance q from it in the flow direction. The circles represent the numerical results when the spheres are centered at the mid-plane (x 2 ¼ X 2 ¼ 0) and the plus signs for the case when the spheres are at x 2 ¼ X 2 ¼ 1 2 . The dashed curve represents the small q approximation given by Eq. (18), and the solid curve represents the large q approximation given by Eq. (20) for FIG. 9. Reduction in the rotational velocity of a sphere (divided by a 3 ) due to the presence of a second sphere at a distance q from it in the flow direction. Refer to Fig. 8 for the captions.

083302-9
Roles of particle-wall and particle-particle interactions Phys. Fluids 23, 083302 (2011) q ¼ 2 for which the corresponding values were, respectively, À 0.216 and 0.134. Finally, the prediction that S p ¼ À2X p for large q was also verified from the results for q > 3:8.

C. Comparison with the numerical results for pairs of spheres
Detailed numerical calculations for selected values of a for a pair of freely suspended spheres with R ¼ 2 and aligned in the flow direction x 2;3 ¼ X 2;3 ¼ 0 À Á were carried out using the method developed by Ozarkar and Sangani 17 using multipoles up to N s ¼ 9, and the results were compared with those obtained for the single sphere with the same a and N s . It was found that the presence of the second sphere changed the rotational velocity of the sphere by less than 1% for a 0:8. This is consistent with our analysis for small a which gave X p ffi À0:005 for R ¼ 2 and c 1 ¼ 0:213, suggesting thereby that the effect of the presence of a second sphere at this distance should have a negligible influence on the rotational velocity of a sphere. The effect on the stresslet was somewhat greater: the stresslets of pairs of spheres at a ¼ 0:6 and 0.8 were found to be lower than those for the single sphere by, respectively, 6% and 8%. Note that for this case, c 3 ¼ 1:78 and S p ¼ À0.4 (for R ¼ 2Þ; hence, the small a analysis predicts a modest reduction in the stresslet due to the presence of a second sphere. The large a analysis based on the lubrication effects, on the other hand, predicted a modest increase in the stresslet. As a consequence, we expect the difference between the stresslet of a pair of spheres aligned in the flow direction and that for the single sphere with same a to remain relatively modest for the entire range of values of a for R ¼ 2.
Similar computations were also carried out for pairs of spheres aligned in the vorticity direction. Both the shortrange lubrication analysis and the long-range asymptotic expressions show that the presence of a second sphere should has no effect on the stresslet or on the rotational velocity of a sphere in accord with the computed values of both X p and S p which were found to be generally small. For example, for R ¼ 2 we obtained X p ffi À0:008 and S p ¼ 0:06. Therefore, we expect the presence of a sphere in the vorticity direction to have a negligible effect on both the rotational velocity and the stresslet of a sphere. This was confirmed from the results obtained by direct computations. For example, for a ¼ 0:8, the stresslet was found to be only 2% greater in the presence of a second sphere placed along the vorticity vector while the angular velocity was 2% lower.

V. AN APPROXIMATE CALCULATION OF THE STRESSLET TO O (/)
Let us now consider the problem of determining the Oð/Þ correction to the average stresslet or, equivalently, the Oð/ 2 Þ correction to the relative viscous dissipation in dilute random suspensions. Detailed calculations involving arbitrarily positioned pairs of spheres between two plane walls are extremely cumbersome, and therefore, we shall be content with only an approximate calculation using Eq. (17) which ignores the effect of higher-order multipoles induced in each sphere. Even this approximate calculation is nontrivial, but, as we shall see, provides an interesting qualitative insight into the role of particle interactions in highly confined suspensions.
To determine the average stresslet of a sphere in a random suspension, we must evaluate the following integral where Pðq; h; x 2 X 2 Þ is the probability density for finding another sphere near ðq cos h; x 2 ; q sin hÞ given that the test sphere is present at 0; X 2; 0 À Á . In trying to evaluate integrals such as Eq. (21), one encounters the difficulty, well-known in the suspensions literature, that such integrals are not absolutely convergent. For example, for the case of unconfined, infinitely extended suspensions examined by Batchelor and Green, 11 S p is given by Eq. (18), hence the integrand in Eq. (21) decreases as R À3 while the volume of the suspension over which the integral must be evaluated increases as R 3 , so that although the integral of S p over a large but finite volume approaches a constant, it depends, in general, on the shape of this domain. Batchelor and Green 11 observed that even though the average stresslet is dependent on the shape of the outer boundaries, the rheology of bulk suspensions is not since the average rate of strain is also similarly shape dependent. This observation was then exploited by them to determine the relation between the average stress and the rate of the strain that involved an integral that is absolutely convergent and independent of the shape of the outer boundary ("infinity") of the suspension or the flow imposed at infinity.
In the present problem, the integrand decreases as q À2 (cf. Eq. (20)) while the volume of the suspension grows as q 2 , so that, once again, the integral is, in general, non-absolutely convergent and the integral as expressed in Eq. (21) is ill-defined. However, for the special case in which the pair probability density is symmetric around the mid-plane for q ) 1, the integration over x 2 of the leading term in Eq. (20) will vanish irrespective of how the integration over a finite domain involving the other two dimensions is carried out. Thus, the remaining integral-which depends on the remainder, i.e., the exponentially decaying terms-is absolutely convergent and the expression (21) which extends the integration to infinity is justified. We must stress, however, that for the more general case, in which the pair probability distribution is not symmetric around the mid-plane, e.g., when all the spheres are confined in one, off-center plane, one must examine carefully the implication of the long-range interactions before proceeding further.
To simplify matters further, we shall consider here only the case in which the probability density is uniform and equal to 3/=½4pa 3 1 À a ð Þ for all possible positions of the second particle which do not lead to an overlap with either of the walls or with the first particle.
The integral in Eq. (21) was evaluated numerically as follows. First, we note that the dependence of the integrand on h is of the form A þ Bcos 2 h, with A and B being the functions of q; x 2 ; and X 2 . Therefore, the integration with respect to h simply equals 2pðA þ B=2Þ), which is the same, incidentally, as the integrand evaluated at h ¼ p=4 multiplied by 2p. Next, since the integrand in Eq. (21), upon integrating over x 2 , decays exponentially with q for h ¼ p=4 (cf. Eq. (20)), the integration with respect to q can be truncated to a finite q c ; in the results presented below, we used q c ¼ 3:5 given that no noticeable change was observed for a few selected calculations with q ¼ 3 or q c ¼ 4:5. Furthermore, we note that, since the pair probability is zero when the spheres overlap, the lower limit on the integration with respect to q can be replaced by q l ðx 2 ; X 2 Þ which is given by q 2 l ¼ 4a 2 Àðx 2 À X 2 Þ 2 for x 2 j À X 2 j < 2a and zero otherwise. We then decomposed S p into two parts: S un p given by Eq. (18) which can become very large for very small a and q and the remainder. The former was integrated analytically with respect to q from q ¼ q l to q ¼ q c and subsequently integrated numerically with respect to x 2 and X 2 using Simpson's rule with Dx 2 ¼ DX 2 ¼ 0:02. The triple integration of the remainder part which varies over the channel width scale was performed numerically with a coarser grid involving about 10 divisions with respect to each integration variable. Note that this part is computationally intensive. We, therefore, used the method of Ozarkar and Sangani with somewhat smaller N r (equal to 12) to evaluate the integrand.
The results of the numerical integration of Eq. (21) are expressed as equal to k p /, and the computed values of k p as a function of a are shown in Fig. 10. We see that these results extrapolated to a ¼ 0 are in agreement with the result k p ¼ 1 obtained by Batchelor and Green 11 (cf. their Eq. (7.2) but without the integral which accounts for the effects or the higher-order pair interaction reflections) for unbounded suspensions using their renormalization method. Our computations, therefore, are consistent with this very significant benchmark result without having to apply a renormalization, albeit for a very specific geometry of confined suspensions with the flow at "infinity" generated by the moving walls. We should note that the modification of the particle interactions due to the presence of the walls (for a ( R ¼ 0ð1ÞÞ played a very important role in obtaining the correct result since we would have obtained k p ¼ À3=2 as a ! 0, if we had simply used S un p , as given by Eq. (18), instead of S p to evaluate Eq. (21)-which is still an absolutely convergent integral since S un p decays as q À3 . Since S p is related to the gradient of the velocity induced at X due to a sphere located at x, the volume integration over x (i.e., over q; :h, and x 2 ) can be can be converted to surface integrals involving the surfaces of the excluded volumes using the divergence theorem as shown in Appendix C. The alternate expression derived there is shown to lead more readily to the result k p ! 1 as a ! 0 and may have offered a slight computational advantage in determining k p for other values of a as well.
As seen in Fig. 10, the contribution from the long-range pair-interactions becomes negative at a % 0:3. This is because, the number of pairs aligned parallel to the walls increases relative to the pairs that are obliquely oriented as a increased and, as noted earlier, the contribution from such pairs is negative, causing the overall contribution to the stresslet from the pair-interactions to decrease. This trend reverses at a % 0:6 since, for larger a; the volume over which pair interactions are significant decreases due to an increase in the excluded volume.
The computed results for k p are well fitted by the relation k p ¼ 1 À 3:3a 1 À 6a 2 þ 9a 3 (22) which is denoted by the solid line in Fig. 10. The coefficient 3.3 in the above equation was chosen such that k p vanishes at a % 0:3, and the other two coefficients were subsequently determined by trial and error. Finally, we note that our analysis assumed that the normalized stresslet of the sphere at x was simply equal to unity. The normalized stresslet of the individual sphere in a random suspension actually varies, of course, with x 2 and this fact may be approximately accounted for by multiplying k p by the average stresslet in the absence of the other spheres. Thus, the relative viscous dissipation may be approximated by Our calculation of the Oð/ 2 Þ coefficient is only approximate since it does not account for the detailed pair interactions among the spheres that are separated by distances R that are comparable to the radius of the spheres. Batchelor and Green 11 have carried out such detailed calculations for the special case of unconfined suspensions ða ! 0Þ and obtained k p ffi 5:0=ð5=2Þ ¼ 2:0, about twice the value we obtained. We expect that the error for moderately large values of a will be smaller since for such suspensions the inter-particle distances will be comparable to the channel width and we have shown that on such length scales: (i) the particle-particle interactions do not affect the lubrication effects significantly and (ii) the average stresslet induced by other spheres decays rapidly. For larger a, for which k p is negative, a

083302-11
Roles of particle-wall and particle-particle interactions Phys. Fluids 23, 083302 (2011) possibility exists for the relative viscous dissipation to have a local maximum at / = À 1=(2kp), which for a ¼ 0.5 occurs at / % 0:5. Of course, this estimate is very approximate considering that our calculation of the Oð/ 2 Þ coefficient only accounts for the leading order term in the particle-particle interactions and that the suspension at / % 0:5 is anything but dilute. Finally, it is also interesting to compare the above result with that given by Davit and Peyla 8 and quoted in the introduction of the present study. The numerical results obtained by these investigators for a ¼ 0:5 correspond to hS Ã i ¼ 2 and k p ¼ À1 which are in good agreement with the values hS Ã i % 1.5 and k p ¼ À1 (cf. Figs. 5 and 10) obtained here.

VI. CONCLUSIONS
Particulate suspensions sheared in the presence of walls show rheological and mobility behavior that is strikingly different from those in the absence of walls. Particle-wall interactions reduce the rotational velocity of the individual particles and increase the stress induced by the particles, with the short-range, lubrication-dominated, particle-particle interactions appearing to play a relatively insignificant role in highly confined suspensions. The particle-particle interactions in the presence of walls are long-ranged in that, the presence of a second particle changes the induced stresslet of a test particle by an amount that is proportional to the inverse of the distance squared compared to the inverse cubed relation in unbounded suspensions. For the special case considered here of random suspensions with uniform probability density, these very long-range interactions, however, do not contribute to the average stress in the suspension. In contrast, when the separation distance between pairs is comparable to the channel width, pair interactions play the most significant role in determining the pair contribution to the total stress and may contribute negatively to the total stress.

ACKNOWLEDGMENTS
We are grateful to Mr. Levan Jibuti for providing the results presented in Figure 3. Ashok Sangani acknowledges the support from the IR=D program while serving at the National Science Foundation.

APPENDIX A: MATCHED ASYMPTOTIC ANALYSIS FOR SMALL SPHERES
Let the center of the sphere lie in the lower half of the channel at a distance y from the lower wall, with a y 1. Note that y is related to x 2 introduced in the main text by y ¼ x 2 þ 1. We shall develop expansions for the stresslet and the rotational velocity in: (i) the outer region, where y ¼ Oð1Þ, and (ii) the inner region, where y ¼ OðaÞ, a distance from the wall comparable to the particle radius. Let us denote the terms of Oða k Þ in the inner and outer expansions of the normalized stresslet by, respectively, S k;in and S k;out and those of the corresponding terms for X by X k;in and X k;out , respectively.
The leading behavior in the outer region corresponds to the trivial case where the walls exert no influence on the freely suspended sphere. This yields S 0;out ¼ 1 and X 0;out ¼ 1=2. Let us, therefore, consider the inner region in more detail where the stretched distance Y y=a between the center of the sphere and the lower wall is Oð1Þ and where one must account for the influence of the lower wall.
Calculations for a simple shear flow around a sphere freely suspended near a single wall have been carried out by several investigators (e.g., Refs. 1, 2,14,[17][18][19]. Particularly noteworthy are those by Chaoui and Feuillebois 22 who computed, with remarkably high precision to 16 significant digits, the translational and rotational velocities of the sphere but not the stresslet. We used the method of Ozarkar and Sangani 17 to determine both the rotational velocity and the stresslet as functions of Y. The results for the stresslet and the rotational velocity were found to be well approximated by the following expressions: For 0 < 1=Y < 0:85, For 0:85 < 1=Y ¼ 1=ð1 þ eÞ < 1, S 0;in ¼ 0:847 ln e À1 À 0:41 þ 1:44e ln e À1 À 0:3e ð1=5Þ ln e À1 þ 0:6376 ; (A3) X 0;in ¼ 0:4218 À 0:025e ð1=5Þ ln e À1 þ 0:6376 : Note that the stresslet approaches an Oð1Þ constant equal to 0:847 Â 5 % 4:23 in the limit of vanishingly small gap between the sphere and the wall while the rotational velocity approaches zero. The expression for the stresslet in the inner region expressed in terms of the outer region variable y is given by S in ! 1 þ 15a 3 =ð16y 3 Þ þ Oða 5 Þ, hence the matching requirement yields that S 1;out ¼ S 2;out ¼ 0 and S 3;out ! 15=ð16y 3 Þ as y ! 0. Similarly, X t 1;out ¼ X t 2;out ¼ 0 and X t 3;out ! À5=ð32y 3 Þ as y ! 0. The Oða 3 Þ corrections to the stresslet and rotational velocity are determined next by treating the particle as a point stresslet with magnitude corresponding to its value for an infinitely wide channel and computing at that point the rate of strain and vorticity induced by its images on the two sides of the channel. Such computations were already made for the special case of the sphere placed at the mid-plane in Sec. II where we gave the corrections in terms of c 1 and c 3 (cf. Eqs. (4) and (6)). These calculations were repeated now as functions of y, and the results were expressed as given by S 3;out ¼ 15 16y 3 þ C 3;reg ðyÞ; X 3;out ¼ À 5 32y 3 À C 1;reg ðyÞ: (A5) Figure 11 shows the results for C 3;reg and C 1;reg obtained using 32 reflections on each side of the channel and subtracting from them the contribution from the first reflection from the lower wall which has already been accounted for by the first terms on the right-hand side of the above expressions.
At y ¼ 1, we obtained C 3;reg ð1Þ ¼ 0:84 and C 1;reg ð1Þ ¼ 0:057, so that S 3;out ð1Þ ¼ ð15=16Þ þ 0:84 ¼ 1:78 ¼ c 3 and X 3;out ð1Þ ¼ Àð5=32Þ À 0:057 ¼ À0:213 ¼ Àc 1 in agreement with the result for the sphere placed at the mid-plane quoted in Sec. II. At the other end, we find that C 3;reg ! 0:729 À 0:96y; C 1;reg ! À0:364 þ 0:99y as y ! 0: (A6) The above results expressed in terms of the inner region variable Y ¼ y=a suggest that the matching requirement is S 1;in ¼ S 2;in ¼ 0 and S 3;in ¼ 0:729 as Y ! 1. Similarly, X 1;in ¼ X 2;in ¼ 0 and X 3;in ¼ 0:364 as Y ! 1. To determine the Oða 3 Þ term in the inner expansion, we must solve, therefore, for a sphere freely suspended at a distance Y from a wall in an undisturbed flow having a rate of strain 0.729 times that of the undistrubed shear flow and a vorticity, that is, 0:364=0:5 or 0.728 times that of the undisturbed shear flow. The two numbers are essentially the same and within the numerical accuracy used in the computations, so that the undisturbed flow is essentially the same as the simple shear flow. It can be shown, therefore, that the Oða 3 Þ terms of the inner expansions are given by Finally, we note that Eq. (A6) suggests that the Oða 4 Þ corrections in the inner region will be nonzero and require solving for the case in which the rate of strain and vorticity of the undisturbed flow increase linearly with Y. We shall not pursue this calculation here but note that it can be shown that its solution will force a correction of only Oða 6 Þ in the outer region. Therefore, the Oða 4 Þ corrections in the outer region are zero.
We now turn to the evaluation of the average stresslet and rotational velocity. For this purpose, we need expressions for these quantities that are applicable for an arbitrary position of the particle. A good approximation to these quantities for the entire range of particle positions can be obtained by adding the inner and outer expansion results and subtracting from the sum the overlapping part. Keeping in mind that the inner region expression is in terms of the stretched coordinate, we write Let V j i ðx; XÞ be the x i -component of the velocity at X due to an applied force of unit magnitude at x along the x jaxis. A complete expression for this velocity may be found in Liron and Mochon. 28 In turn, the velocity induced by a freely suspended sphere with its center at x can be approximated by that due to a point stresslet S jk when the radius of the sphere is very small (compared to the channel width and the distance from any other particle) and is given by For the case of shear flow of interest to us in the present study, The stresslet induced in a freely suspended sphere at X due to this velocity induced by the sphere at x can be determined from the symmetric part of the velocity gradient. In particular, the 1,2-component of the induced stresslet is given by S p 12 ¼ ð20=3Þpa 3 e 12 with e 12 ¼ 1 2 083302-13 Roles of particle-wall and particle-particle interactions Phys. Fluids 23, 083302 (2011) Note that S p , the particle-particle interaction contribution to the normalized stresslet introduced in the main text (cf. Eq. (17)) is related to S p 12 by means of the relation S p a 3 ¼ 3S p 12 =ð10pa 3 Þ, and therefore, For jX À xj ( 1, V j i can be approximated by the velocity induced by a point force in an unbounded medium which is given by where R i ¼ X i À x i . On substituting Eq. (B4) into Eq. (B3), one obtains the expression for the induced stresslet in the unbounded shear flow as given by Eq. (18) in the main text. For jX À xj ) 1, the point-force velocity is given by (cf. Eq. (51) in Liron and Mochon 28 ) Oðq À1=2 e À2:1062q Þ þ ðd i2 d ja þ d j2 d ia Þ Oðr a q À1=2 e À2:1062q Þ þ d ia d jb Oðr arb q À1=2 e Àpq=2 Þc; (B5) where f ðx 2 Þ ¼ 1 À x 2 2 , r a ¼r a q ¼ X a À x a , q 2 ¼ R 2 1 þ R 2 3 , and a and b are indices that take values 1 and 3 (but not 2). Now, noting that V 1 1 ¼ 3=ð16pÞ bf ðx 2 Þ f ðX 2 Þ cos 2h q À2 þ Oðq À1=2 e Àpq=2 Þc; V 2 2 ¼ 3=ð16pÞ bOðq À1=2 e À2:1062q Þc V 2 1 ¼ V 1 2 ¼ 3=ð16pÞ bOðq À1=2 e À2:1062q Þc; where h ¼ cos À1 ðc 2 =qÞ, and substituting for the various velocity components in Eq. (B3), one obtains Eq. (20) in the main text. Finally, the corresponding expression for X p , the particle-interaction contribution to the rotational velocity of a particle (cf. Eq. (17)), can be derived in a similar manner to yield X p ¼ ÀS p =2, a relatively simple relationship since the leading terms, of oðq À2 Þ for large q, in both S p and X p are related to the partial derivative of ð@ 2 V 1 1 =@x 2 @X 2 Þ. The complete expression for V j i involves two infinite series, both of which decay exponentially with q-one series decays as q À 1 2 e Ànpp=2 n ¼ 1; 2; :::Þ while the other decays as q À 1 2 e ÀImðz m Þp=2 ðm ¼ 1; 2; :::Þ, where z m satisfies the equation z ¼ sinh z and Im stands for the positive imaginary part. For large q, therefore, the leading exponentially decaying term is q À 1 2 e Àpp=2 . For the special case, x 2 ¼ X 2 ¼ 0, however, this term is identically zero, and the leading exponentially decaying term is Oðq À1=2 e ÀImðz 1 Þq=2 ) with Imðz 1 Þ=2 ¼ 2.1062.

APPENDIX C: AN ALTERNATE EXPRESSION FOR k p
We begin with the expression for the average particle contribution to the normalized stresslet as given by /k p ¼ a 3 2ð1 À aÞ ð 1Àa aÀ1 hS p iðXÞdX 2 (C1) with hS p iðXÞ ¼ ð V ½s p ðX; xÞpðx; XÞdV: (C2) Here, pðx; XÞ is the probability of finding a sphere at x given that a sphere is present at X and a 3 S p ðX; xÞ is the induced normalized stresslet on the sphere centered at X due to the stresslet on the sphere centered at x. The volume integration over x must be carried over the domain V which is the entire volume of the suspension minus the excluded volumes V w (defined by 1 À jx 2 j < aÞ near each wall and V p near the sphere centered at X (cf. Fig. 7). For the special case of a uniform pair probability distribution considered in the main text, we have pðx; XÞ ¼ 3/=½4pa 3 ð1 À aÞ p o and, therefore, upon substituting for s p from Eq. (B3) into Eq. (C2), we obtain where B represents outer boundary of the suspension and S w and S p represent the surfaces of the excluded regions (cf. Fig. 7). Now, at the outer boundary B, n 2 ¼ 0 and u 2 is exponentially small and, therefore, the middle integral on the righthand side of Eq. (C4) can be neglected. Thus, one needs to evaluate only the first and third integrals which must be done numerically as in the method described in the main text.
Note that in the limit of small particles ða ! 0Þ, the first integral being OðaÞ (since u 1 ¼ OðaÞ near the walls) can be neglected, and if 1 À jX 2 j ) a, then the integral over S p can be evaluated using Eq. (B4) to determine u i , and hence, hS p iðXÞ in a straightforward manner to yield k p ¼ 1-the same result as one obtained by Batchelor and Green 11 in the absence of higher order interactions between the spheres, i.e., in the absence of the integral in their Eq. (7.2).