A Method for Determining Stokes Flow Around Particles Near a Wall or in a Thin Film Bounded by a Wall and a Gas-Liquid Interface

A method for determining Stokes ﬂow around particles near a wall or in a thin ﬁlm bounded by a wall on one side and a nondeformable gas-liquid interface on the other side is developed. The no-slip boundary conditions at the wall are satisﬁed by constructing an image system based on Lamb’s multipoles. Earlier results for the image systems for the ﬂow due to a point force or a force dipole are extended to image systems for force or source multipoles of arbitrary orders. For the case of a ﬁlm, the image system consists of an inﬁnite series of multipoles on both sides of the ﬁlm. Accurate evaluation of the ﬂow due to these images is discussed, including the use of Shanks transforms. The method is applied to several problems including chains of particles, radially expanding particles, drops, and porous particles. © 2008 American Institute Physics


I. INTRODUCTION
The problem of determining velocity field around particles in close proximity of a wall or a gas-liquid interface arises in the analysis of many physical and biological phenomena including cell adhesion, slurry transport, and particulate flows in microfluidic devices. Understanding these phenomena require numerical simulations of particle-wall and particle-particle interactions involving hundreds or even thousands of particles. For example, there is considerable interest in simulating biofilms 1-3 formed by bacteria such as Escherichia coli. These are microcolonies in which bacteria adhere to each other and to a substrate to protect them from harsh environment while still providing adequate nutrients for them to grow. Transport phenomena-flow and mass transfer-play a role equally important as the cell biology in determining the morphology of such microcolonies, and this complex interplay of transport and biological factors can be studied by numerical simulations of a flow involving many particles near a wall.
Mo and Sangani 4 and Sangani and Mo 5 described an efficient method for computing low Reynolds number hydrodynamic interactions in a system consisting of many particles. The velocity induced by a spherical particle is expressed in terms of Lamb's multipoles, i.e., the coefficients appearing in Lamb's spherical harmonic solution of the Stokes equations of motion. 6,7 Repeated applications of translation formulas for Lamb's multipoles are then used to combine the velocity induced by groups of particles resulting in an algorithm for which the total computational effort increases only linearly with the number of particles. The method was devised for particles confined to a periodic unit cell or particles in an unbounded space.
The present study is intended to give expressions that are necessary so that the method of Sangani and Mo 5 can be extended to account for the presence of a rigid plane wall, a nondeformable gas-liquid interface, or both. The method is based on determining an image system for Lamb's multipoles. Expressions for the velocity induced by a point force or an array of point forces in a bounded or an unbounded fluid medium have been given by several investigators in the past. [8][9][10][11][12][13][14][15] Blake 10 gave an image system for a point force near a plane. This image system consists of a point force, a force doublet, and a source dipole. Blake and Chwang 16 derived image systems for higher-order singularities due to a point torque, a source, or a source dipole. The present study gives the image system for a general singularity of arbitrarily high order. Results for the image system due to force multipoles of arbitrary order were also given by Cichoki et al. 17,18 The main difference is that the present study gives the image system in terms of Lamb's multipoles that can be conveniently incorporated into the fast multipole algorithm described by Sangani and Mo. 5 The method is used to solve a number of simple flow problems involving either a single particle, a chain of particles, or particles in a thin liquid film bounded on one side by a wall and by a gas-liquid interface on the other side.
The organization of the article is as follows. In Sec. II, we derive the image system for a particle near a wall. In Sec. III, we consider several problems of single or multiple particle interactions with a wall. The results for a single particle are shown to be in excellent agreement with those published in the literature. In Sec. IV, we describe the image system for a film bounded by a wall and a gas-liquid interface. In Sec. V, we present results for spherical objects-rigid particles, porous particles, and drops and bubbles-and in Sec. VI, we summarize the important results.

II. LAMB'S MULTIPOLES AND THEIR IMAGE SYSTEMS
The velocity induced by a particle in an incompressible, Newtonian fluid at small Reynolds number is given by Lamb's solution: where r = x − x p , x p denotes the center of the particle, r = ͉r͉ is the distance of the point x from the center of the particle, c n s = 2 − n 2n͑2n − 1͒ , b n s = n + 1 n͑2n − 1͒ , ͑2͒ and p n s , n s , and n s are decaying spherical harmonics of order −n − 1 with their singularity at x p . These spherical harmonics are expressed in terms of spherical coordinates centered at where Y nm 0 = P n m ͑ p ͒cos m p and Y nm 1 = P n m ͑ p ͒sin m p , P n m being the associated Legendre function. The spherical polar components are related to the Cartesian components r i by the usual relations: r 1 = r p , r 2 = r ͱ 1 − ͑ p ͒ 2 cos p , ͑4͒ r 3 = r ͱ 1 − ͑ p ͒ 2 sin p .
The constants P nm k , T nm k , and ⌽ nm k will be referred to as Lamb's multipoles. Note that the components of the hydrodynamic force F i , torque L i , and the particle-induced stresslet S ij are related to the first few Lamb's multipoles as given by Mo and Sangani: 4  where is the viscosity of the fluid. For the case when the particle is a source of fluid, the mass flow rate ejected by the particle is given by A. An image system for P nm k Let us now determine the image system that accounts for the presence of a wall. We choose the Cartesian coordinate system with its center on the wall and with its x 1 axis passing through the center of the particle. Thus, we take h being the distance between the particle center and the wall, and ê 1 is the unit normal along the x 1 axis. Let u im represent the velocity induced by an image system so that the combined velocity u = u p + u im satisfies the no-slip boundary condition at the wall, i.e., at x 1 = 0. The velocity induced by the image system will also be expressed in Lamb's solution form ͓Eqs. ͑1͒-͑3͔͒ with Lamb's multipoles denoted by P nm k,im , T nm k,im , etc., and ͑r , p , p ͒ replaced by ͑r im , im , im ͒. Here, r im = x − x im , and x im =−hê 1 is the location of the image singularities.
Since the governing equations of motion are linear, we can determine the image system by considering one multipole at a time. As explained in the Appendix, the image system for the P nm k multipole requires, in general, a set of seven image multipoles given by Let us now consider a few special cases. For n = 1 and m = k = 0, corresponding to P 10 0 =−F 1 / ͑4͒, or a point force normal to the wall, the image system consists of only three nonzero multipoles consisting of a point force, a force dipole, and a source dipole ͑P 10 0,im =−P 10 0 , P 20 0,im =−4hP 10 0 , and ⌽ 10 0,im =−h 2 P 10 0 ͒. The resulting velocity due to the point force and its image system is given by

͑17͒
For n = m = 1 and k = 0, corresponding to a force parallel to the wall, the image system consists of four multipoles corresponding to a point force, a force dipole, a torque, and a source dipole ͑P 11 0,im =−P 11 0 , P 21 0,im =2hP 11 0 , T 11 1,im =−hP 11 0 , and ⌽ 11 0,im = h 2 P 11 0 ͒. The resulting velocities due to the point force and its images are given by

͑18͒
Both of the above expressions agree with those given by Blake and Chwang. 16 B. An image system for T nm k Using the same method as outlined in the Appendix, it can be shown that the image multipoles for T nm k are given by The last equation applies to the case when n −1ജ m. ⌽ nm k,im must be set to zero when n Ͻ m or when k = 1 and n = 0. The case n = 1 corresponds to the image system for a point torque. It can be shown that the image system given above also agrees with that given by Blake and Chwang. 16 C. An image system for ⌽ nm k Finally, the image system for ⌽ nm k is given by The special case n = m = k = 0, corresponding to a mass source of fluid, requires an image system consisting of a point force dipole, a mass source, and a source dipole ͑P 20 0,im =−8⌽ 00 0 , ⌽ 00 0,im = ⌽ 00 0 , and ⌽ 10 0,im =−2h⌽ 00 0 ͒. The resulting velocity is, once again, in agreement with that given by Blake and Chwang. 16 D. An image system for nondeformable gas-liquid interface The image system for a particle near a surfactant-free gas-liquid interface in the limit of large surface tension and vanishingly small gas viscosity is relatively straightforward.
The boundary conditions at the gas-liquid interface are as follows: ͑i͒ the velocity of the liquid normal to the interface is zero and ͑ii͒ the tangential stress components are zero. For a particle placed at a distance h below an interface, the image multipoles at a distance h above the interface are simply given by The motion of a sphere near a planar interface of two immiscible, viscous fluids of arbitrary viscosities has been examined by Lee et al. 19 by using the method of images and reflections. A general solution of the Stokes flow near a planar interface was also given by Palaniappan 20 who gave explicit expressions for the first few image singularities. The expressions given above are in agreement with those of these investigators for the special case when the viscosity ratio of the two fluids vanishes.

E. Determination of Lamb's multipoles
With the image system described above for satisfying the boundary conditions at a wall or at a gas-liquid interface, we are now in a position to determine the hydrodynamic force on the particles or their velocities. For this purpose, we first need to determine Lamb's multipoles for each particle by satisfying the boundary conditions at the surface of all the particles. The velocity at a point in the fluid can be expressed as where u ϱ ͑x͒ is the undisturbed flow velocity, i.e., the velocity field in the absence of the particles, N is the total number of particles, and u p,im is the velocity due to the image system for particle p ͑denoted earlier for the sake of brevity by simply u im ͒.
To satisfy the boundary conditions on the surface of particle p, the flow induced by the images, undisturbed flow, and all the other particles are expanded near the center of the particle, i.e., near x p . Since the velocity due to these flows is regular at x p , it can be expressed in terms of the regular part of Lamb's solution: ⌽ nm k,r ΅ r n Y nm k ͑ p , p ͒.

͑32͒
The coefficients P nm k,r , T nm k,r , and ⌽ nm k,r are related to the image multipoles, the other particles' Lamb's multipoles and their images, and the imposed flow parameters. The transla- .͔ Finally, an application of the boundary conditions on the particle surface leads to the following set of equations ͑taken from Mo and Sangani 4 with slight modification͒: g n s P nm k + na 2n+3 ͑n + 1͒͑2n + 1͒͑2n + 3͒ g n r P nm k,r = 0, For rigid particles, the coefficients f n s , f n r , g n r , g n s , and h n s are all unity and V nm k and ⍀ nm k are related to the translational velocity V i and angular velocity ⍀ i through the relations V 10 0 = V 1 , V 11 0 =−V 2 , ⍀ 10 0 = ⍀ 1 , ⍀ 11 0 =−⍀ 2 , and ⍀ 11 1 =−⍀ 3 . Mo and Sangani 4 gave expressions for the coefficients f n s , f n r , etc., for other spherical objects ͑drops, bubbles, or porous par-ticles͒ and these will be given in Sec. V where we consider some problems concerning nonrigid particles. Finally, we also note that the above set of equations ͓Eqs. ͑33͒-͑35͔͒, with their right-hand side modified, may also be used for computing flow around charged particles with thin double layers as was done by Kang and Sangani. 21 Since the coefficients ⌽ nm k,r , T nm k,r , and P nm k,r for a given particle depend on the Lamb's multipoles of all the other particles and images, Eqs. ͑33͒-͑35͒ represent an infinite set of equations. This set is truncated as suggested by Mo and Sangani 4 to result, in general, a total of ͑3N s 2 −1͒N equations, where N s represents the maximum value of n used in P nm k multipoles. Upon solving these equations, one determines the particle multipoles and hence the force and torque acting on the particles. In suspension problems, the force and torque acting on the particles are specified and one determines the translational and rotational velocities of the particles.

III. PARTICLES NEAR A WALL
To test the validity of the image system derived here and to assess how the results depend on the truncation order N s , we have solved a number of problems involving one or more particles near a wall or in an unbounded fluid and compared the results to those available in the literature. [22][23][24][25][26][27][28][29][30][31][32][33] We first consider the case of a single spherical particle of radius a held fixed in a parabolic flow given by

͑36͒
␥ being the shear rate of the undisturbed flow at the wall, the viscosity of the fluid, and G the pressure gradient. In this case, all force and torque components vanish except for F 2 and L 3 . Since these quantities must be linear in ␥ and G, we express them according to ͑37͒ Table I shows the convergence of the results for the nondimensional force and torque coefficients f ␥ , f G , t ␥ , and t G with N s for the special case of a touching particle, i.e., for h = a. For this case, only Lamb's multipoles ͑P n1 0 , T n1 1 , ⌽ n1 0 ͒ are nonzero, and the total number of unknowns reduces to 3͑N s −1͒. Chaoui and Feuillebois 28 used a bipolar coordinate expansion technique to obtain accurate estimates ͑to 16 significant digits͒ for the shear flow coefficients, i.e., for f ␥ and t ␥ . The coordinate system chosen by these investigators did not allow direct evaluation for the touching case and therefore the results were presented for h = a͑1+͒. Their results for =2ϫ 10 −6 and 5 ϫ 10 −6 extrapolated to = 0 are also shown in Table I. The results for the parabolic flow coefficients, i.e., for f G and t G , were obtained previously by Pasol et al. 29 As can be seen in Table I, our technique yields accurate estimates of all these four coefficients. It may be noted that prior to Chaoui and Feuillebois, 28 O'Neill 25 also obtained estimates for f ␥ and t ␥ . His estimates ͑1.7009 and 0.943 993͒ are also in agreement with the results presented here-a slight discrepancy is observed only in the fifth digit for the value of f ␥ .
Next, we consider the translational motion of a particle in a fluid at rest at infinity. The hydrodynamic force and torque acting on the particle are given by

͑38͒
where V is the translational velocity of the particle and n is the unit vector perpendicular to the wall and pointing into the liquid. Table II shows the results for the coefficients f v w , f h w , and t h w for selected values of h / a. All the three coefficients diverge as h → a because of the lubrication forces ͑see Kim and Karrila 7 or Jeffrey and Onishi 34 ͒ in the narrow gap region between the particle and the wall. The table shows the results for N s = 2 and the smallest N s for which the results of numerical computations have converged to four significant digits ͑the fifth digit being rounded͒. We have included the results for N s = 2 in the table for comparison sake since the large-scale simulations of interacting particles are often carried out with N s = 2. For h Ͻ 1.5a, it may be necessary to add lubrication forces explicitly and the force dipoles to account for the velocity they induce as was done by Sangani and Mo. 35 Since the force dipoles due to lubrication effects by Sangani and Mo 35 were expressed in terms of Lamb's multipoles, the image system derived here could be readily used to incorporate the lubrication effects while keeping N s small.
Accurate estimates of f v w for the indicated values of h / a in Table II were previously determined by Brenner 22 and for f h w and t h w by Chaoui and Feuillebois. 28 Results presented here are in perfect agreement with their investigations.
The corresponding results for the force and torque on a particle translating near a gas-liquid interface are also presented in Table II, where the superscript G is used to indicate the presence of a gas-liquid interface. For this case, n is defined as the unit normal pointing into the gas phase. Once again, our results are in perfect agreement with those reported in the literature-the case of a particle moving towards the gas interface was treated by Brenner 22 and the case of a particle moving parallel to the gas-liquid interface, being the same as two particles moving parallel to each other in an unbounded fluid, was treated by Stimson and Jeffery. 30

A. Chain of particles
Next, we consider two problems that are motivated by their potential application to biofilms. As mentioned in the Introduction, biofilms are formed by cells that adhere to each other or to a substrate. The biofilms are generally heterogeneous with morphology that is governed by several variables including the "stickiness" of the cells and the fluid shear stress at the wall.
We shall explore the growth of biofilms through dynamic simulations of many particles in a future study. Here, we consider a simple system of touching particles, which form a straight rod as shown in Fig. 1, placed in a simple shear flow ͓cf. Eq. ͑36͒ with G =0͔. The particles are spherical and it is assumed that the bonds between a particle and its neighbors can withstand specified shear or tensile forces as will be explained below. Figure 2 shows the hydrodynamic force acting on the particles for two different cases. The forces are nondimensionalized by 6au 2 ϱ ͑x 1 j ͒, x 1 j being the distance between the center of particle j and the wall. In the first case, a chain of N particles is aligned perpendicular to the wall ͑ =0͒. The particle closest to the wall is labeled 1, and the farthest is labeled N. In the second case, a chain of 2N particles is aligned along the x 1 axis with no wall. The particles at the end of the chain are labeled as −N and N, and the symmetry consideration requires that the force on particle −j is simply the opposite of the force on particle j. In Fig. 2, we see that for both cases, the particles at the end of the chains experience essentially the same force. This suggests that the effect of the wall is unimportant for the particles far away from the wall. On the other hand, the force on the particle closest to the wall ͑j =1͒ is significantly greater than the force on a particle in the middle of a chain twice the length in an unbounded flow. The image system, which consists of a nega-TABLE II. Results for the force and torque coefficients for a particle moving parallel or perpendicular to the wall or the free surface. The numbers in the brackets indicate N s required to obtain the result accurate for four significant digits. Also shown are the results for N s =2.

063301-5
A method for determining Stokes flow around particles Phys. Fluids 20, 063301 ͑2008͒ tive force and higher-order multipoles, induces a net flow that is along the positive x 2 axis, causing the drag on the particle to increase in the case of a particle attached to a wall.
In the absence of a wall, the other particles induce a flow along the negative x 2 axis, causing a decrease in the drag. Finally, it may be noted that the nondimensional drag on a single particle in a shear flow was about 1.7 ͑cf. Table I͒ compared to about half that value for the first particle in a chain of 15 particles. Thus, the drag coefficient decreases with the number of particles in both wall-bounded and unbounded flows. The force balance on particle N requires that the bond between particles N − 1 and N must exert a −F N force on particle N and F N on particle N − 1. The force balance on particle N − 1 then requires that the bond between particles N − 1 and N − 2 must withstand a force equal to F N + F N−1 . Thus, in general, the bond between particles j and j + 1 must be able to withstand force B j given by where ê s and ê t are unit vectors perpendicular and parallel to the chain as given by ê s = − sin ê 1 + cos ê 2 , ê t = cos ê 1 + sin ê 2 , ͑40͒ and b j s and b j t are the nondimensional shear and tensile forces on the bond between particles j and j + 1. The required strength for the bond connecting particle 1 to the wall will be denoted by B 0 . Since this bond is always aligned along the x 1 axis, the tensile and shear forces for this bond are the components along, respectively, the x 1 and x 2 axes.
When N is large, the slender body theory 36 may be used to estimate the force on each particle and hence the bond forces. For the = 0 case, the shear component of the bond forces is expected to scale as N 2 / log N. Figure 3 shows b 1 s as a function of N obtained by first determining the force on each particle and then using Eq. ͑39͒. Results for b 0 s are similar. For large N, these results are well fitted by the relations If we treat the chain of spheres as equivalent to a cylinder of radius a and length 2aN and neglect the effect of the wall, then the use of the slender body theory would have given us the coefficient of log N in the above expressions to equal 3 4 , close to the fitted values. Let us now suppose that all the bonds between the particles are identical and can withstand a maximum shear force S p , and that between the particle and the substrate can withstand S s . Similarly, let the maximum tensile strengths be given by, respectively, T p and T s . If the shear rate ␥ is progressively increased, then when it exceeds a critical value ␥ c , to be determined, one of the bonds will break. For = 0, the bond will break by exceeding the shear strength. If the bond with the substrate is relatively weak, i.e., if B 0 s / S s Ͼ B 1 s / S p , then the bond with the substrate will break first, and the entire chain of particles will eventually be carried away by the fluid ͑assuming that the particles are rough so that the lubrication force remains finite for nearly touching particles͒. Otherwise, the bond connecting particles 1 and 2 will break. The critical shear for this bond to break is given by If the shear rate is progressively increased beyond ␥ c given by the above expression, then the chain of particles 2 through N will rotate around an axis passing through the center of particle 1 until the chain inclination angle is such that the shear force exerted by the bond equals S p . This assumes that the strong lubrication forces will keep the center to center particle distance close to 2a and that the time scale for forming new bonds is relatively short compared to the time required for the chain to rotate by an angle . The critical shear rate for the chain to remain at an angle is therefore given by To determine the inclination of chain as a function of the shear rate, we therefore need to compute b 1 s ͑͒. The results are shown in Fig. 4. We see that as the shear rate is progressively increased, the chain inclination angle also increases.
In order that the chain simply rotates around particle 1 when all bonds among particles can withstand equal strength, the tensile force on bond 1 must not exceed its strength. This requires that 6␥ c ͑͒a 2 b 1 t ͑͒ Ͻ T p . The results for b 1 t are shown in Fig. 5 for three different N. As one would expect, b 1 t goes through a maximum at = 45°, the principal extension axis for a simple shear flow. The condition for the chain to rotate is given by According to the slender body theory, one expects b 1 t ͑͒Ϸ͓b 1 s ͑0͒ / 2͔sin cos and b 1 s ͑͒Ϸb 1 s ͑0͒cos 2 , where b 1 s ͑0͒ is the nondimensional shear force on the bond between the first two particles when the chain is aligned perpendicular to the wall that can be estimated using Eq. ͑41͒. Figure 6 shows the computed values of b 1 t / b 1 s for a few selected values of N and the predicted asymptote for N → ϱ. Clearly, the chain of spheres will snap off if the shear rate is sufficiently high so that the chain rotates to an angle s given by s Ϸ tan −1 ͑2T p / S p ͒ ͑in radians͒. Figure 4 may then be used to determine the shear rate beyond which the chain snap off will occur. This maximum shear rate is given by Substituting for ␥ p , s , and b 1 s , the last being estimated according to the slender body theory, we find that the maximum shear rate for the model chain considered here to remain attached to the wall is given by

B. Radially expanding particle
The second problem we consider is the problem of determining the translational velocity of a growing spherical particle near a wall. The boundary condition on the particle surface is given by where n is the unit outward normal on the particle surface. We want to determine V such that the net force on the particle vanishes. The cells in biofilms consume nutrients ͑referred to sometimes as substrates͒ and grow. In a quiescent fluid, the morphology of the biofilm will be controlled by how the cells grow, divide, and rearrange. This provides a motivation for first calculating in detail the flow due to a single growing particle.
For this problem, ⌽ 00 0 =−a 2 U 0 ͓cf. Eq. ͑8͔͒, and since the force on the particle vanishes, P 10 0 = 0. The results of numerical computations for V / U 0 for selected values of h / a are given in Table III which shows the minimum N s required to determine V / U 0 correct to four significant digits. The lubrication effects due to an expanding particle are similar to that of a translating particle to leading order as the gap between the particle and the wall approaches zero. As a result, V / U 0 → 1 as = h / a −1→ 0. The results of numerical computations are well described by where the coefficient of the log term is derived by using the lubrication analysis, while the coefficients of and 2 log are obtained by fitting the numerical results. The negative coefficient of the log͑1 / ͒ term in the above expression suggests that the distance between the wall and the closest point on the surface of the particle will decrease with time indicating that the particle surface will come arbitrarily close to the wall as the time progresses.
When the particle is at great distances from the wall, V / U 0 decays as ͑a / h͒ 2 . The large distance asymptote is related to the velocity induced by the image system for ⌽ 00 0 and can be shown to be given by

063301-7
A method for determining Stokes flow around particles Phys. Fluids 20, 063301 ͑2008͒ Figure 7 shows the results of numerical computations and the above two limiting expressions.

IV. PARTICLES IN A THIN FILM
We now consider a more difficult problem of determining force and torque on a particle in a thin film bounded on one side by a rigid wall and by a gas on the other side. The wall is at x 1 = 0, and the gas-liquid interface is nondeformable and at x 1 = H. The gas is assumed to have negligible viscosity. The particle center is at ͑h ,0,0͒. Let us first consider the image system for the case when a point force in the direction of the x 2 axis is applied at ͑h ,0,0͒. We then require image multipoles at x 1 =−h to satisfy the no-slip conditions at the wall and at x 1 =2H − h to satisfy the conditions at the gas-liquid interface at x 1 = H. Introduction of the latter, however, results in the violation of the no-slip boundary condition at the wall and requires additional image multipoles at x 1 =−2H + h. Likewise, additional multipoles are required at x 1 =2H + h to correct the boundary condition at the gas-liquid interface due to the multipoles at x 1 =−h. This leads to an infinite number of image multipoles on both sides of the liquid film as shown in Fig. 8. Since the image system for ͓P nm , T nm , ⌽ nm ͔ multipoles consists, in general, of multipoles of order n + 2, every second reflection on the wall side involves multipoles that are two orders higher than the previous ones. Thus, the complete image system for a point force does not only involve infinite number of images but also multipoles that are increasingly higher order as the distance of the image point from the wall or the gas-liquid interface increases.
We should note here that a similar case of a fluid bounded on either side by rigid walls has been investigated by a number of investigators. 11,15,[37][38][39][40] This case also requires an infinite number of images. To avoid the evaluation of these images, Liron and Mochon 11 extended the Stokes equations of motion to the region outside the walls and used a method of Fourier transform to obtain Green's function, i.e., the velocity field due to a point force. The resulting expression, which contains an infinite series of Fourier-Bessel integrals, is quite cumbersome and contains terms that converge very slowly when the point force is near one of the walls. Although Staben et al. 37 recently regularized the slowly converging integrals to devise an efficient method for computing Green's function, the expressions for Green's function remain cumbersome, and the method presented here may provide an alternate approach to solve such problems.
The velocity induced by a point force singularity at ͑h ,0,0͒ may be expressed as where v ij ϱ corresponds to the unbounded flow, i.e., v ij ϱ ͑x͒ = 1

͑49͒
with s = x − hê 1 , s = ͉s͉, and v ij im is the contribution from the images on the two sides of the film. We describe first a method that can be used for determining accurately the velocity field due to point force or higher-order singularities. An approximate method for evaluating the same using Shanks transforms 41,42 is described in Sec. IV A.
To determine v 22 im , we need to determine the regular coefficient ⌽ 11 0,r at ͑h ,0,0͒ ͓cf. Eq. ͑32͔͒ when the multipole P 11 at that point is set to unity. We first determine the image multipoles at x 1 =−h and determine contribution to ⌽ 11 0,r from those multipoles as described in Sec. II. Next, we determine the multipoles on the gas-liquid side and determine their contribution to ⌽ 11 0,r . This is followed by the second reflection on the wall side, and so on. This leads to a series whose first few terms for the special case, h = H / 2, are given below:

͑50͒
We have evaluated the first 240 terms ͑120 reflections on each side of the film͒ in the above series. The 120th reflection involved an image system with multipoles P n1 with n = 61. The series in Eq. ͑50͒ follows a regular sign pattern ͑ϪϩϪϪϩϪϩϩ͒ repeating every ninth term. The computational effort increases roughly linearly with the square of the number of reflections. To estimate the contribution from higher reflections, we found correlations for c n from the  computed 240 terms by plotting every ninth term in the above series. An example of this is shown in Fig. 9 where we have plotted the product of the c 8p+1 coefficient and distance d 8p+1 of the image contributing to that coefficient from ͑h ,0,0͒ as a function 1 / p. The dashed line in the Fig. 9 represents a fit given by with A 1 = −0.645 33, B 1 = 0.221 56, and C 1 = −1.252 39. The R 2 of the fit is better than 0.995. Since d 8p+1 increases linearly with p, the above correlation suggests that c 8p+1 decreases as 1 / p for large p.
We obtained similar fits for the other coefficients c 8p+j ͑j =2,3, ... ,8͒ with an R 2 of 0.995 or better. The sum of the coefficients of leading terms, i.e., A j ͑i =1,2, ... ,8͒, was found to be −0.0066, or nearly zero, indicating that the sum of every eight terms decreases roughly as 1 / p 2 , and that series ͑50͒ will converge. We corrected each A j by 0.0066/ 8 so that the corrected sum of these coefficients is exactly zero and refitted the terms c 8p+j to determine modified B j and C j . These modified fits for the coefficients c n were then used to determine the contribution for n Ͼ 240 to series ͑50͒. This procedure yielded for the velocity induced by the images at the point of force singularity. The same procedure was used to determine Green's function at other points in the fluid. Table IV gives Green's function at several points along the x 2 and x 3 axes ͑with x 1 = h͒. We shall compare next the results of numerical computations to an asymptotic expression for the velocity at large distances from the point force.
At large distances from the force singularity, we expect the velocity to become parallel to the wall. Thus, we write u i = ū i ͑͒f 1 ͑x 1 ͒, i = 2,3, ͑53͒ with = ͑x 2 2 + x 3 2 ͒ 1/2 and f 1 = 3 Note that u i given by the above expression satisfies the no-slip condition at the wall and the stress-free condition at the gas-liquid interface. The average of f 1 over the film thickness is unity so that ū i represents the film-averaged velocity. The pressure satisfies Laplace equation in the x 2 -x 3 plane, and for large distances, it is expected to be given by the dipole approximation The film-thickness-averaged velocity satisfies Darcy's equation according to which the velocity is proportional to the pressure gradient. Thus, Substituting for the velocity and pressure fields into the Stokes equation of motion then yields These expressions are similar to those given by Bhattacharya et al. 43 for the flow between two walls. In fact, Green's function in the present problem is given by the sum of Green's functions for the case of two walls separated by distance 2H and with point forces applied at ͑h ,0,0͒ and ͑2H − h ,0,0͒. By comparing with their expression ͓Eq. ͑31͒ by Bhattacharya et al. 43 ͔, we find that  Table IV and the dashed line represents the predictions based on the above expressions. We have similarly computed the pressure along the x 2 axis and verified that the computed values are in excellent agreement with Eq. ͑55͒ with A given by Eq. ͑58͒.
In summary, the method for computing the velocity due to a point force or any other singularity consists in first evaluating as many reflections as possible ͑typically 120 on each side of the film͒, correlating these terms, and using the correlations to evaluate the contributions from the higher reflections.

A. Approximate methods for evaluating sums
We have seen that the series produced by the images on the two sides of the film converges rather slowly. In largescale simulations, it will be necessary to seek methods that may be used to determine Green's function or its derivatives at lower computational cost. The convergence of a series may often be accelerated through an application of a suitable transform. A number of transforms for evaluating a series are described by Shanks. 41,42 Let be the sum of the first n terms of a series. Then one may transform the original sequence ͕A n ͖ to a new sequence ͕B k,n ͖ ͑n = k , k +1,k +2, ...͒ by the kth order transform defined by · · · · · · · · · · · · ⌬A n−1 · · · · · · ⌬A n+k−1 Έ Έ · · · · · · · · · · · · ⌬A n−1 · · · · · · ⌬A n+k−1 Έ , ͑60͒ where ⌬A n = A n+1 − A n . The resulting sequence generated through such a transform is denoted in the operator form by e k ͑A n ͒. The simplest and well-known transform corresponds to k = 1 for which It is often found that the transformed sequence ͕B k,n ͖ converges more rapidly than the original sequence ͕A n ͖ ͑provided, of course, that the series is convergent͒ as the transformation may reduce oscillations or "noise" from the original series. Further noise from the transformed series may be reduced by repeated application of the transform. Thus, for example, e k p ͑A n ͒ refers to a sequence generated by applying the e k transform p times on the sequence ͕A n ͖. Note that a finite sequence ͕A n ͖ with n Յ N reduces to a sequence of N −2kp numbers upon application of e k p . Given a finite number of terms in a series, one may experiment with different combinations of k and p to determine an efficient algorithm for estimating an infinite series. Table V shows the results of applying various transforms to the series for evaluating the regular part of Green's function v 22 im . For the purpose of comparison, we took the contributions from the first 33 reflections on the wall and gas sides. c n corresponds to the combined contribution from the first n reflections on the two sides of the film. Due to space limitations, only nine terms are shown in the table even though the results indicated used n up to 33. Also, only five significant digits are shown in the table for each result even though we used a double precision arithmetic. The table shows the sequence resulting from the e k ͑k =1,2,4͒ transform applied several times. With the total of 33 terms, it is possible to apply e k transform 16/ k times to yield a single estimate. The estimates thus obtained are shown in Table VI. Also shown in that table are the errors in the resulting estimates when compared to the accurate estimate obtained by the method described earlier. We see that the least error occurs when the e 4 transform is used. It is possible that this particular transform is better since the sign pattern for the coefficients c n has a frequency of 4, corresponding to a regular sign pattern repeating after every four reflections on each side, or a total of eight terms.
We have also compared the estimates of Green's function obtained by using the e k p transform with a limited number of reflections to the more accurate estimates obtained  earlier at other positions along the x 2 and x 3 axes ͑cf. Table  IV͒. It was found that the e 6 transform gives slightly better estimates than e 4 at most positions. With only 25 reflections on each side, e 6 2 gives estimates that are accurate to within 2% for the most cases examined.
Although the use of the e k transform reduces the computational time substantially ͑note that, as mentioned earlier, the number of multipoles required by the image system increases with the number of reflections, and hence the computational cost of c n roughly increases in proportion to n and the total computing time of Green's function increases as N 2 , N being the total number of reflections͒, we were initially surprised to find that the error did not decrease exponentially with p, as was the case in the number of series reviewed by Shanks. 42 We surmise that these nonlinear transforms work extremely well when the functional dependence of c n on n is relatively simple. In the present problem, the image system continues to become more complex with increasing n and these trends are apparently not captured very efficiently with very few terms. To support this hypothesis, we have tested the case when all the higher-order multipoles are ignored and only the point force images are kept in all the reflections. In this case, we found that even with 17 reflections, the repeated application of the transforms gave results that were accurate to six significant digits.

V. RESULTS FOR THE FORCE AND MOBILITY OF SPHERICAL OBJECTS
To determine the force or torque on particles in a film, we need to determine the regular coefficients ͑P kl r , T kl r , ⌽ kl r ͒ at the center of a particle in terms of Lamb's multipoles ͑P nm , T nm, ⌽ nm ͒ at the center of the same particle or any other particle in the liquid film. The procedure described in the previous section for determining ⌽ 11 k,r due to a P 11 k multipole at ͑h ,0,0͒ can be used to compute all regular coefficients in terms of Lamb's multipoles. Once these coefficients are determined, it is relatively straightforward to compute the force or torque on spherical objects. The contributions from the lower-order multipoles ͑n Յ 1͒ to lower-order regular coefficients were obtained by the accurate method described in Sec. IV in which up to 240 terms were evaluated and correlated to estimate the remaining terms in the series. For higher-order multipoles, we simply used the estimates from the first 240 reflections. All calculations to be presented below were carried out with N s = 9 which was sufficient for most cases except for a / h close to unity.

A. Rigid particles
We first consider the problem of a particle held fixed at x 1 = h = H / 2 in a film ͑cf. Fig. 8͒ with parabolic flow at infinity given by where U is the film-thickness-averaged velocity of the fluid. The force and torque acting on the particle are written as with 3 =3U / 4h being the vorticity of the undisturbed flow at the particle center. f 1 s and h 1 s are the coefficients that appear in Eqs. ͑33͒-͑35͒. For rigid particles, both these coefficients are unity but we include them in the expressions for the force and torque as they will be useful later when we give results for porous particles and drops. The results of computations for the force and torque coefficients for selected values of a / h are shown in Table VII. The force coefficient increases monotonically from 9 8 at a / h = 0 to about 1.9 at a / h = 1. For small values of a / h, the force coefficient can be shown to be given by the asymptotic relation The factor 9 8 in the numerator corresponds to the undisturbed flow evaluated at the center ͑x 1 = h = H / 2͒ of the particle while v 22 im is related to the velocity induced by the image system at the particle center ͑cf. Sec. IV͒. As seen in Fig. 11, the asymptotic relation gives accurate results up to a / h = 0.3.
To check the results of our analysis by an independent method, we also obtained the force coefficient by using a standard computational package due to FLUENT 6.3 which employs a finite volume method. The computational domain assumed periodicity along the x 3 axis with a period S. The TABLE VII. The force and torque coefficients f p and t p for a rigid particle placed in a parabolic flow in a film

063301-11
A method for determining Stokes flow around particles Phys. Fluids 20, 063301 ͑2008͒ upstream condition of uniform flow ͑u 2 = U͒ was applied at x 2 =−S / 2 while the downstream condition on the plane x 2 = S / 2 was specified to correspond to the "outflow condition." The Reynolds number was set to 0.005 and a / h was set to 0.5. We carried out computations for two cases: S =10h and 20h. For the latter case, the flow domain was divided into 400 000 cells and the computational time on Dell Precision 360 was about 7 h. The computed force coefficients for these two cases were, respectively, 1.79 and 1.70. Assuming that these forces vary linearly with h / S, we estimated f p to equal to about 1.61 for the case of a single particle in an unbounded ͑in the x 2 and x 3 directions͒ film. This compares well to the 1.55 obtained by the present method.
The results for the torque coefficient t p are given in Table VII. This coefficient decreases monotonically from its value of unity at a / h → 0 to about 0.3 at a / h = 1. The image system produces vorticity at the particle position that is opposite to that of the undisturbed flow, and this is responsible for the reduction in the torque as a / h increases. The asymptotic formula with f p given by Eq. ͑64͒ is shown to be in excellent agreement with the computed results as can be seen in Fig. 12.
It is interesting to compare the results for the force and torque on a particle in a thin film to those for a particle placed in an unbounded parabolic flow studied in Sec. III. Taking ␥ = 3 2 and G =− 3 4 in Eq. ͑36͒ yields the same undisturbed flow as that given by Eq. ͑62͒. If the flow for x 1 Ͼ H is entirely neglected, then f p computed here will equal 3 2 Using the values of f ␥ and f G from Table I yields   3 2 f ␥ − 3 8 f G = 1.82 for a particle resting on a plane wall. This is less than 5% lower than the computed value of f p = 1.9 for a particle in a film, suggesting that the bounded nature of the film flow is not very important as far as the force is concerned. The torque on the particle in a film on the other hand is significantly different. If we had used the results from the unbounded parabolic flow, we would have obtained t p =2t ␥ − t G = 0.89 compared to the 0.3 obtained for a particle in the thin film. As mentioned above, the image system for the thin film induces vorticity that is opposite to the vorticity of the undisturbed flow. This may be responsible for the significant difference in the results for the two cases.
Next, we consider the problem of determining the force and torque on a translating and rotating particle through a quiescent film. The force and torque acting on the particle with velocity V 2 ê 2 and angular velocity ⍀ 3 ê 3 are expressed as As mentioned earlier, f 1 s and h 1 s are unity for rigid particles. Note that the reciprocal theorem requires that f r = 4 3 t t . Thus, we need to determine only f t , t t , and t r . The results of computations for these coefficients are shown in Figs. 13 and 14. The lines in these figures correspond to the asymptotic results that will be discussed next.
In the limit = ͑h − a͒ / a → 0, the force and torque coefficients diverge due to lubrication effects in the gap between the wall and the surface of the particle. No lubrication effects arise in the gap between the particles and the gas-liquid interface. Chaoui and Feuillebois 28 give asymptotic expres- sions for the force and torque coefficients in the limit → 0 for the case of a particle translating or rotating parallel to a wall. Those expressions are applicable here but with different O͑1͒ constants. Combining their results with our numerical results, we obtain In Eqs. ͑67͒-͑69͒, the coefficients of log and log are taken from Chaoui and Feuillebois 28 while the O͑1͒ and O͑͒ are estimated from our numerical results for a / h equal to 0.85 and 0.9. In the limit a / h → 0, it can be shown that As seen in Figs. 13 and 14, the above asymptotic results compare well to the computed results for f t , t r , and t t .
The above results can be combined to determine the translational and rotational velocities of a freely moving ͑force and torque free͒ particle due to imposed parabolic flow ͓cf. Eq. ͑62͔͒. The results are given in Figs. 15 and 16. We see that the velocity of the particle is greater than the mean velocity of the fluid for a / h less than about 0.7. This is because the center line velocity is greater than the mean velocity of the fluid. The rotational velocity is roughly equal to half the vorticity of the undisturbed flow at the particle center ͑3U / 8h͒. We also see that V remains greater than ⍀a for the entire range of values of a / h indicating that the surface on the particle closest to the wall will have a positive velocity component along the x 2 axis. For larger values of a / h, both translational and rotational velocities begin to hinder due to the presence of the wall. The exact expressions and the asymptotic expressions in the limit → 0 ͑setting f 1 s = h 1 s = 1 for rigid particles͒ are given by The errors on the extreme right side of the above expressions are O͑1 / log ͒. The lines in Figs. 15 and 16 represent various asymptotic results. For small a / h, the asymptotic expressions given earlier for the coefficients f p , f t , t p , etc., were substituted in the middle expressions in Eqs. ͑73͒ and ͑74͒. These are indicated by the dotted lines. In the opposite limit, i.e., for → 0, we have shown two asymptotic results. The first is obtained again from the middle expressions in the above equations by substituting f p = 1.9 and t p = 0.3 corresponding to the results for these coefficients for a / h = 1 and the asymptotic expressions Eqs. ͑67͒-͑69͒ for the rest of the force and torque coefficients. These asymptotic results are given by the dashed line. The dot-dashed line represents the simplified form of the resulting asymptotic expressions given on the extreme right-hand side of Eqs. ͑73͒ and ͑74͒.

B. Drops and bubbles
Having determined the regular coefficients ͑P kl r , T kl r , ⌽ kl r ͒ in terms of multipoles ͑P nm , T nm , ⌽ nm ͒ through the image system, it is relatively easy to determine the force or mobility of other spherical objects or slightly nonspherical particles. We shall present here the results for porous particles and for drops although they could be equally applied with slight modification to moderately deforming drops or charged particles.
For the motion of drops or bubbles under conditions of small Capillary numbers, for which the surface tension is large enough to keep the fluid particles nearly spherical, the boundary conditions simplify to continuity of tangential stress and velocity components at the drop-fluid interface. The normal component of the velocity satisfies u i n i = V i n i at the interface, V i being the translational velocity of the drop.

063301-13
A method for determining Stokes flow around particles Phys. Fluids 20, 063301 ͑2008͒ The resulting equations from the application of the boundary conditions can be put in the form given by Eqs. ͑33͒-͑35͒ with where ␤ = / is the ratio of viscosities. The above expressions, taken from Mo and Sangani, 4 correct one important typographical error: expression ͑76͒ for f n r and g n s had the inverse sign missing. The case n =1 in h n s can be treated separately by requiring that T 1m = 0, which is equivalent to requiring that the net torque on the fluid particles must be zero. Table VIII gives the results for the translational velocity of a force and torque free fluid particle in a parabolic imposed flow. We see that the velocity variation with a / h reduces as the viscosity ratio is decreased. For bubbles ͑␤ Ϸ 0͒, the nondimensional velocity varies only by about 15% as a / h is varied from 0 to 0.9.

C. Porous particles
Porous particles are often used as models for proteins and macromolecules. We assume that the flow inside the porous particles satisfies the Brinkman equation where ij is the stress tensor inside the particles, and k is the Darcy permeability of the porous medium. The flow outside the particles, of course, satisfies the Stokes equations of motion. The boundary conditions at the particle surface are the continuity of velocity and traction. Using a solution in terms of spherical Bessel functions to describe the flow inside a particle and Lamb's solution for the flow outside the particle, and applying the boundary conditions at r = a yields Eqs. ͑33͒-͑35͒ with the following expressions for the coefficients f n s , f n r , etc., that appear in Eqs. ͑33͒-͑35͒ ͑Mo and Sangani 4 ͒: where ␤ = a / ͱ k and G n ͑␤͒ = i n ͑␤͒ / i n−1 ͑␤͒. Here,i n ͑␤͒ ϵ͑ / 2␤͒ 1/2 I n+1/2 ͑␤͒ is the modified spherical Bessel function. The limit ␤ → ϱ corresponds to the case of a rigid, nonporous particle. ͓Note that Eq. ͑81͒ corrects the expression given by Mo and Sangani 4 in which ␤ in the denominator in the second term on the extreme right side of Eq. ͑81͒ was missing.͔ The force and torque on a particle placed in an undisturbed parabolic flow given by Eq. ͑62͒ are given by Eq. ͑63͒. For a porous particle, f 1 s and h 1 s can be determined by substituting n = 1 in Eqs. ͑80͒ and ͑83͒, respectively. In some applications, it is of some interest to estimate the average velocity of the fluid through the porous particle. The average velocity is given by where V is the volume occupied by the particle. In deriving the last equality in the above expression, use was made of the fact that u i is related to the divergence of stress tensor ͓cf.
Eq. ͑79͔͒ and that the integral of the latter can be related to the force by means of the divergence theorem. Table IX gives the results for the force and torque coefficients f p and t p for a porous particle placed in a parabolic flow in a film for various values of a / h and ␤. For a translating and rotating porous particle through a quiescent film, the force and torque acting on a particle are given by Eq. ͑66͒ with f 1 s and h 1 s given by Eqs. ͑80͒ and ͑83͒, respectively.   Table X gives the results for f t and t t coefficients for a translating porous particle, and Table XI gives the results for f r and t r coefficients for a rotating porous particle in a quiescent film for various values of a / h and ␤. The asymptotic expressions for the force and torque coefficients in the limit a / h → 0 for the rigid particle given in Sec. V A are also applicable for a porous particle with appropriate expressions for f 1 s and h 1 s . In the limit of very high permeability, the force coefficient f p can be determined from the second equality in Eq. ͑84͒ with u i replaced by the undisturbed velocity u i ϱ . This yields the following for ␤ 1: Thus, for highly porous particles, the force coefficient decreases as a / h is increased. The opposite is the case for nearly nonporous particles for which the presence of the wall increases the drag. We therefore expect that for some intermediate values of permeability, the force coefficient as a function of a / h will go through a maximum. Indeed, this is what is found for ␤ =1.

VI. SUMMARY
We have described a method for accounting for the presence of a rigid wall or a nondeformable gas-liquid interface on the Stokes flow interactions of spherical particles. The flow induced by the particles and the image system that ac-counts for the boundary conditions at the wall or the interface are both expressed in terms of Lamb's multipoles. This will allow an easy integration into an O͑N͒ algorithm described earlier by Sangani and Mo. 5 The method is applied to a number of problems including chains of particles, drops, and porous particles.

ACKNOWLEDGMENTS
This work was partially funded by the Syracuse Center of Excellence under the EPA Grant No. X-83232501-0.

APPENDIX: DERIVATION OF THE IMAGE SYSTEM
The no-slip boundary condition at the wall is equivalent to the conditions where 1 is the x 1 component of the vorticity = ٌ ϫ u. Substituting Eq. ͑3͒ into Eq. ͑1͒ and using the properties of spherical harmonics, we find that where the summation over n, m, and k is the same as in Eqs. ͑1͒ and ͑2͒, and Y n,m k,p ϵ Y nm k ͑ p , p ͒. We consider separately the image system for each of Lamb's multipoles. By setting all multipoles except one, P nm k , to zero using  To satisfy the above conditions, the image system must involve, in general, seven multipoles: P nm k,im , P n+1,m k,im , T nm 1−k,im , T n−1,m 1−k,im , ⌽ nm k,im , ⌽ n−1,m k,im , and ⌽ n−2,m k,im , all of which must, of course, be proportional to P nm k . Furthermore, a dimensional analysis suggests that T nm 1−k,im and ⌽ n−1,m k,im must be proportional to h and ⌽ nm k,im must be proportional to h 2 . Substituting h = im r im and using the recursion relations among the Legendre functions and their orthogonality, we obtain eight linear equations among the seven multipoles. Using seven of these equations to determine the image multipoles and then verifying that the remaining equation is automatically satisfied, we obtain the image system for P nm k as given in the main text ͓cf. Eqs. ͑10͒-͑16͔͒.