Drag and Torque on Clusters of N Arbitrary Spheres at Low Reynolds Number

Journal of Colloid and Interface Science 229, 184–195 (2000) doi:10.1006/jcis.2000.6981, available online at http://www.idealibrary.com on

Drag and Torque on Clusters of N Arbitrary Spheres at Low Reynolds Number A. V. Filippov1 High Temperature Chemical Reactions Engineering Laboratory, Department of Chemical Engineering, Yale University, New Haven, Connecticut 06520-8286 E-mail: [email protected] Received February 9, 2000; accepted May 16, 2000

Hydrodynamics of particle clusters suspended in viscous fluids is a subject of considerable theoretical and practical importance. Using a multipole expansion of the flow velocity in a series of spherical harmonics, Lamb’s fundamental solution of the Stokes flow outside a single sphere is generalized in this work to the case of N nonoverlapping spheres of arbitrary size with slip boundary conditions. The expansion coefficients are found by transforming the boundary conditions to the Lamb form and by transforming the spherical coordinates and solid spherical harmonics centered at different spheres. The problem is reduced to the solution of the linear system of equations for the expansion coefficients, which is carried out numerically. Based on the developed theory, the relation between the hydrodynamic and gyration radius of fractal-like aggregates with different structure is established. In another application, an asymptotic slip-regime dependence of the aggregate hydrodynamic radius on the Knudsen number and the number of particles is found by performing calculations of drag forces acting on the gas-borne fractal-like and straight chain aggregates. A good agreement is shown in comparing predictions of the described theory with available experimental and theoretical results on motion of various small sphere clusters in viscous fluid. °C 2000 Academic Press Key Words: Stokes flow; multipole expansion; aggregated particles; slip regime; mobility radius.

1. INTRODUCTION

The theory of hydrodynamic interactions between suspended particles and the surrounding fluid phase plays a fundamental role in understanding phenomena in a wide variety of systems such as colloids, polymers, and aerosols. For example, a proper account of hydrodynamic forces is required for the calculation of transport properties, such as the viscosity coefficient of suspensions and sedimentation velocity of aggregated particles. Incorporation of the long-range, many-body hydrodynamic forces between particles poses a major difficulty in theoretical analysis of these systems. In many cases, the particle shape is nearly 1 Present address: Mason Lab, P.O. Box 208286, Yale University, New Haven, CT 06520-8286.

0021-9797/00 $35.00

C 2000 by Academic Press Copyright ° All rights of reproduction in any form reserved.

spherical and the flow Reynolds number is small enough to assume the creeping flow model, which substantially simplifies the problem. Hydrodynamic forces in a collection of N spheres in Stokes flow have been studied by many authors starting from Smoluchowski (1), who used the method of reflections. Because of its complexity and poor convergence, the application of this method was restricted to cases involving few spheres (2). A great practical demand for effective and simple evaluation of hydrodynamic forces acting on multisphere configurations for N À 1 motivated creation of simplified models such as the “porous sphere” model of Debye (3) and the “modified porous sphere” model of Rosner and Tandon (4). However, the validity of these simple methods could be fully proven only by comparison with mathematically rigorous models for creeping flow around arbitrary conglomerates of spheres, which were not available until recently. The theories of Kirkwood and Riseman (5) and Rotne and Prager (6) described the hydrodynamic interactions between widely separated spheres, but were less accurate for close distances between spheres. The boundary element method in the general case requires solution of a two-dimensional integral equation and is effective only if the conditions of either aggregate symmetry or slenderness are met (7). In 1987, Durlofsky et al. (8) proposed an effective scheme for rigorous calculation of hydrodynamic forces on a set of N arbitrary spheres. It was based on the expansion of the Oseen equation relating the velocity field in the fluid to forces induced on the particle surfaces (9) and the Fax´en theorem for a sphere immersed in a flow field. Exact solution of a two-sphere problem (10) was used to account for the lubrication effect in the nearfield. As a result, a linear system of equations for hydrodynamic forces and their moments (torques and stresslets) was obtained. This scheme was later improved by Ladd (11) and by Cichocki et al. (12), who took into account the higher order multipole moments of hydrodynamic forces. A different approach to the problem has been developed by Hassonjee et al. (13) who used a superposition of N fundamental solutions of Lamb (14), expanding the velocity in series of

184

185

DRAG AND TORQUE ON CLUSTERS OF N SPHERES

spherical harmonics centered around individual spheres. The resulting linear system of equations for the coefficients of these expansions was found by transformation of spherical coordinates from one center to another, numerical integration of the equations with respect to the azimuthal angle, and collocation in longitude to force the no-slip boundary condition on each of the N spheres. In comparison to the approach of Durlofsky et al. (8) and Cichocki et al. (12) the velocity and pressure fields were readily found, along with the hydrodynamic forces and torques, with any desired accuracy, controlled by the number of terms retained in the series and the number of collocation points. However, the necessity of preliminary numerical integration for determining the coefficients of the system appeared to be a major drawback of the method of Hassonjee et al. (13). The approach developed in this work also uses the multipole expansion of the velocity in series of solid spherical harmonics centered at individual spheres. The expansion coefficients of the boundary conditions are also expanded in series of spherical harmonics, similar to the one-sphere case (2), truncated at some level L. However, in contrast to the scheme of Hassonjee et al. (13) the direct origin-to-origin transformation of spherical harmonics is used along with the transformation of spherical coordinates. This enables us to avoid the integration and collocation steps of Hassonjee et al. (13) and yields a system of 3 × N × L × (L + 2) linear equations for the expansion coefficients. The accuracy of the multipole expansion calculations is easily controlled by choosing the truncation level L, and the results converge very rapidly even for touching spheres. A similar technique was developed by Jeffrey and Onishi (10) and Mackowski (15) in their analyses of hydrodynamics of twosphere aggregates. More recently, Mo and Sangani (16) used multipole expansions of the fluid velocity field around particles to model infinite concentrated suspensions. In contrast to present work, Mo and Sangani (16) dealt with particle motion in cells with periodical boundary conditions and did not account for velocity slip effects. In the present paper, to facilitate application of the method to the cases of small, but nonzero Knudsen number, general slip boundary conditions are imposed on the sphere surfaces. This is practically important, for example, in calculations of hydrodynamic interactions between submicrometer aggregated aerosol particles and air at the sea level. The developed technique was applied to find the drag force and slip correction factor for the gas-borne fractal-like and chain aggregates under the slipregime conditions. The method is tested by comparison with available theoretical and experimental data on motion of sphere aggregates in viscous fluids.

∇ ·v=0

[2]

ri → ∞ : v → v∞

[3]

αv λ Πr t = Ui , i = 1, . . . , N , [4] µ where v is the local fluid velocity, p is the pressure, µ is the fluid dynamic viscosity, ri is the radius vector with origin at the center of the ith sphere, v∞ is the flow velocity at infinity, Ui is the local velocity on the surface of the ith sphere, αv is the velocity slip coefficient, and λ is the mean free path. Πr t is the tangential component of hydrodynamic stress Πr , acting on the surface element of a sphere (2): µ ¶ ri ∂v v m − [5] + ∇(ri · v). Πr = − p + m ri ∂ri ri ri |ri | = ai : v −

The velocity slip boundary condition [4] accounts for fluid velocity discontinuity near the particle surface, which becomes negligible in the continuum limit l/a j → 0. Following Happel and Brenner (2) and Mackowski (15), an equivalent set of three scalar boundary conditions can be constructed by applying the operators ri ·, −ri ∇· and ri · ∇× to both sides of Eq. [4],

∂ ∂vr − αv λ ∂ri ∂ri · ∂ × sin θi ∂θi µ 1 − αv λ

·

vr = Uri

¸

1 ∂ ¡ 2 ¢ αv λ ri vr + 2 2 2 ∂r ri i ri sin θi µ ¶ 2 ¸ ∂ vr ∂vr sin θi + = −∇ · Ui ∂θi ∂ϕi2 ¶ ∂ (∇ × v)r = (∇ × Ui )r , ∂ri

[6]

[7] [8]

where ri , θi , and ϕi are the spherical coordinates centered about the ith sphere. No restrictions will be imposed on the particle motion, so that the following results are applicable to N contacting particles in an aggregate, and to N independently moving separated spheres. 3. THE SOLUTION OF THE N-SPHERE PROBLEM

3.1. General Form of the Solution The general solution of the external boundary problem [1]– [3], [6]–[8] can be given in terms of solid spherical harmonics: v = v∞ +

N X

Vj

[9]

j=1

2. FORMULATION OF THE PROBLEM

Consider a low Reynolds number steady flow of viscous fluid past a set of N arbitrary nonoverlapping spheres with radii a j , described by equations µ1v = ∇ p

[1]

Vj =

∞ · X

¢ ¡ j j ∇ × r j χ−(n+1) + ∇8−(n+1)

¸ (n − 2) (n + 1) j j 2 − r ∇ p−(n+1) + r j p−(n+1) . µ2n(2n − 1) j µn(2n − 1) n=1

[10]

186

A. V. FILIPPOV

Equation [10] has the form of Lamb’s general solution (14) of the creeping flow Eqs. [1]–[2] (here in notation of Happel and j j Brenner (2)) outside a single sphere j with χ−(n+1) , 8−(n+1) , j and p−(n+1) being the following linear combinations of solid harmonics u i− mn with origin at the center of the jth sphere, j

8−(n+1) =

n X

j j− amn u mn

[11]

n X 1 j j j− p−(n+1) = bmn u mn µ m=−n

[12]

m=−n

j

χ−(n+1) =

n X

j j− cmn u mn

[13]

Pnm (cos θ j ) exp(imϕ j ),

[14]

m=−n j− = u mn

1 r n+1 j

where Pnm is the associated Legendre function. The velocity field [9] satisfies Eqs. [1]–[2], as well as the boundary condition [3] at infinity. Each of the functions V j can be considered as a perturbation of the flow due to the presence of the jth particle.

Cmnkl = (−1)k+l

The boundary conditions [6]–[8] can be used to find unknown i i i , bmn , and cmn , determining the solution. An elcoefficients amn egant way to do this is to express contributions V j of all spheres in the total velocity field [9] in terms of the harmonic functions centered at any given sphere i. For this purpose, Eq. [10] can be modified, ∞ · X

¡ ¢ ¡ ¢ j j ∇ × ri χ−(n+1) + ∇ × Ri j χ−(n+1)

n=1

−

¢ (n − 2) ¡ 2 j r + 2ri · Ri j + Ri2j ∇ p−(n+1) µ2n(2n − 1) i

(n + 1) j ri p µn(2n − 1) −(n+1) ¸ (n + 1) j Ri j p−(n+1) + µn(2n − 1)

[15]

where Ri j is the vector, connecting the centers of spheres j and i as shown in Fig. 1. The next step is to use the transformation of solid harmonics from origin j to origin i j− u mn (r j , θ j , ϕ j ) =

l ∞ X X

ij

Cmnkl u i+ kl (ri , θi , ϕi )

[17]

l=0 k=−l l k u i+ kl = ri Pl (cos θi ) exp(imϕi ),

where Ri j , θi j , ϕi j are the coordinates of vector Ri j in the spherical system centered at particle j (here, the misprint of Ref. (17) in the exponent of (−1) is corrected). The equations for the coefficients in expansions [11]–[13] can be now obtained by substituting Eqs. [9]–[10] in the boundary conditions [6]–[8] and taking into account Eq. [17]. The nontrivial part of these lengthy calculations is connected with expressing the products of the components of vector Ri j , solid harmonics u i+ kl , and their . After some elementary, but gradients in terms of functions u i+ kl lengthy algebra, the recurrence relations for Legendre polynomials (18) yield the equations Ri j ·

[16]

[18]

where u i+ kl are the positive spherical harmonics of degree l. If

(n + l − m + k)! j− u (Ri j , θi j , ϕi j ), (n − m)!(l + k)! (m−k)(n+l) [19]

µ

j

+ ∇8−(n+1) +

Ri j = r j − ri ,

the length of vector Ri j is greater than r j , the transformation ij coefficients Cmnkl can be expressed in terms of negative spherical harmonics ij

3.2. Finding the Expansion Coefficients from the Boundary Conditions

Vj =

FIG. 1. Definition of vectors in transformation of spherical coordinates and spherical harmonics from origin i to origin j.

=

ri ri

¶ u i+ mn

(n + m)ζi j ri i+ (n − m + 1)ζi j i+ u m(n−1) + u m(n+1) (2n + 1) (2n + 1)ri à i+ ! u (m+1)(n+1) ηi j i+ + − ri u (m+1)(n−1) (2n + 1) ri +

(n + m)(n + m − 1)ξi j ri i+ u (m−1)(n−1) (2n + 1)

−

(n − m + 1)(n − m + 2)ξi j i+ u (m−1)(n+1) (2n − 1)ri

[20]

i+ Ri j · ∇u i+ mn = (n + m)(n + m − 1)ξi j u (m−1)(n−1) i+ + (n + m)ζi j u i+ m(n−1) − ηu (m+1)(n−1) [21]

187

DRAG AND TORQUE ON CLUSTERS OF N SPHERES

¡ ¢ i£ ri · ∇ × u i+ (n + m)(n − m + 1)ξi j u i+ mn Ri j = − (m−1)n ri ri ¤ i+ − mζi j u i+ [22] mn + ηi j u (m+1)n ξi j =

1 1 (Ri j x + i Ri j y ); ηi j = (Ri j x − i Ri j y ); 2 2

[23]

ζi j = Ri j z ,

where the Ri j x , Ri j y , and Ri j z are the projections of the vector Ri j to the coordinate axes of Cartesian frame. The resulting infinite i i i , bmn , and cmn has set of linear equations for the coefficients amn the form (n + 1) i (n + 1)ai i − amn + b ai 2(2n − 1) mn + ai(2n+1)

∞ X l N X X ¡

ij

j

ij

j

ij

j

Dklmn akl + E klmn bkl + Fklmn ckl

¢

j=1 l=1 k=−l j6=i i = ai2n+1 X mn

[24]

1 i (1 + 2nαv Kni )(n + 1)(n + 2)amn ai2 − [n + 2(n 2 − 1)αv Kni ] ∞ X l N X X ¡ (2n+1)

+ ai

(n + 1) i b 2(2n − 1) mn ij

j

ij

j

ij

j

G klmn akl + Hklmn bkl + L klmn ckl

[32]

¢

N X

Fj

[33]

j=1 N X (T j + R j × F j ),

[34]

j=1

[25]

[1 + (n + 2)αv Kni ]n(n +

where R j is the vector connecting a central point of the aggregate (e.g., the centre of mass) and the centre of the jth particle.

i 1)cmn

∞ X l N X X ¡

ij

j

ij

j

Mklmn bkl + Nklmn ckl

¢

3.3. Expansion Coefficients for N Rigid Spheres

j=1 l=1 k=−1 j6=i i = ai(n+1) Z mn

[26] Kni =

λ , ai

[27]

where Kni is the Knudsen number based on the radius of ij ij ij ij ij sphere i, the coefficients Dklmn , E klmn , Fklmn , G klmn , Hklmn , ij ij ij L klmn , Mklmn , and Nklmn are all linear combinations of the elij ements of the transformation matrix Cklmn [19], given in the i i i Appendix, and X mn , Ymn , and Z mn are the coefficients in the expansions n ∞ X X ri i · Ui = X mn Pnm (cos θi ) exp(imϕi ) [28] ri m=−n n=1 n ∞ X X

i Ymn Pnm (cos θi ) exp(imϕi )

[29]

n=1 m=−n

ri · ∇ × Ui =

[31]

where ex , e y , and ez are unit vectors in the direction of each of the three Cartesian coordinate axes. After finding the torques and forces acting on each of the spheres, the determination of the total hydrodynamic force F and torque T is straightforward,

T=

i = ain Ymn

−ri ∇ · Ui =

·µ ¶ 1 j j F j = −4π µ b11 − b−11 ex 2 µ ¶ ¸ 1 j j i + i b11 + b−11 e y + b01 ez 2 ·µ ¶ 1 j j T j = −8π µ c11 − c−11 ex 2 µ ¶ ¸ 1 j j i + i c11 + c−11 e y + c01 ez , 2

F=

j=1 l=1 k=−l j6=i

+ ai(2n+1)

In the case of two spheres (N = 2), Eqs. [24]–[26] can be reduced to the form given by Mackowski (15) for two-sphere aggregates by choosing the z axis of the Cartesian coordinate system parallel to vector R12 . Once the flow velocity field is found, the hydrodynamic frictional force F j on sphere j and hydrodynamic torque T j about its center can be expressed in terms of the expansion coefficients (2)

∞ X

n X

n=1 m=−n

i i i , Ymn and Z mn on the right-hand side of the The terms X mn system [24]–[26] can be easily found in a specific case of rigid spheres, when the particle surface velocity Ui is described by

Ui = Ui0 + ωi × ri ,

where Ui0 is the velocity of the center of the ith sphere and ω i is its angular velocity. For any vector A, the following equation holds due to the definition of Legendre polynomials Pnm (in notation of Ref. (18)): ¡ ¢ 1 ri · A = P1−1 (cos θi ) exp(iϕi ) Aix − i Aiy + P10 (cos θi )Aiz ri 2 ¡ ¢ + P11 (cos θi ) exp(−iϕi ) Aix + i Aiy [36] 1 P1−1 (cos θi ) = − sin θi ; P10 (cos θi ) = cos θi ; 2 P11 (cos θi ) = sin θi .

i Z mn Pnm (cos θi ) exp(imϕi ). [30]

[35]

[37]

Substituting Eq. [35] in the left-hand side of the system

188

A. V. FILIPPOV

[28]–[30] and equating the coefficients of spherical harmonics term-by-term results in the equations ¢ 1 1¡ i i δn U − iU0y 2 0x ¡ i ¢ i = − U0x + iU0y δn1

i = X 1n i X −1n

[38] [39]

i i 1 X 0n = U0z δn i =0 Ymn i Z 1n

=

¡

ai ωix

¡

[40] −

iωiy

[41]

¢

δn1

[42]

¢ i

i = −2ai ωix + iω y δn1 Z −1n

[43]

i Z 0n = 2ai ωzi δn1 ,

[44]

where δnk is the Kronecker delta. 4. RESULTS

4.1. Numerical Implementation Based on the theory described above, a FORTRAN computer code for determining the flow around N arbitrary rigid spheres has been developed. An IMSL subroutine DLSLCG using the LU-factorization was applied for solution of the linear system of Nex = 3 × N × L × (L + 2) equations for the expansion coefficients obtained by truncating the series in Eqs. [24]–[26] after l = L terms. The truncation level L was defined by observing the convergence of the results and depended on the required accuracy, the kind of calculated parameters, and the intersphere distances. For example, Table 1 shows the calculated total drag force acting on quiescent aggregates of 8 equal contacting spheres of radius a in multiples of the drag force 6πµav∞ on an isolated sphere in the same flow. The flow is directed normal to the plane of the “rectangular” aggregate (two touching parallel straight chains of 4 spheres) and to the major axis of the straight chain of 8 spheres. Although the number of expansion coefficients Nex grows rapidly with the truncation level L, so does the accuTABLE 1 Convergence of the Calculated Drag Force for Various Configurations of 8 Spheres Calculated total drag force Maximal order Number of expansion of harmonics coefficients “Cube” “Rectangle” Straight chain L Nex 2×2×2 2×4×1 8×1×1 1 2 3 4 5 6 7

72 192 360 576 840 1152 1512

2.2861 2.3342 2.3408 2.3414 2.3417 2.3418 2.3418

2.8623 2.8782 2.8824 2.8836 2.8837 2.8838 2.8838

3.4231 3.4803 3.4814 3.4821 3.4822 3.4822 3.4822

racy of the calculation, and 1% accuracy for total drag force can be achieved already with L = 2 ÷ 3. Similar convergence has been also found for distributed parameters such as drag force and torque on individual particles in aggregates of equal spheres, as well as for total drag and torque on arbitrary aggregates. Convergence of individual drag and torque on smaller spheres in aggregates of unequal spheres can be worse, but this does not change the rapid convergence of total parameters because of the smaller contribution of these particles. Due to the efficiency of this method, the solution of creeping flow problems for aggregates of up to some hundred spheres can be performed on average workstations and even PCs. In the case of colliding particles and aggregates, the development of extremely high flow velocity gradients may impede the convergence of the method. In these situations, well-developed methods involving the lubrication theory “patches” (19) can be used to modify the algorithm to restore convergence. In the present paper, we consider dynamics of solid aggregates without reciprocal spherule motion, where the lubrication effects are absent and the developed theory can be used in its simplest form. All numerical calculations described in this paper have been carried out on a Dell Dimension XPS R450 personal computer (Pentium II 450 MHz with 256 Mb RAM). The minimal truncation level L = 3 was determined by considerations of accuracy. The computation time strongly depended on the order Nex of the linear system and in the following examples varied from some seconds to some minutes. For instance, a calculation for 8 spheres with L = 5 involved a solution of Nex = 840 equations and took 18 s. Hydrodynamic forces and torques, as well as the flow field have been calculated for various aggregate configurations including numerically generated fractal-like aggregates. 4.2. Particle Agglomerates in Continuum Regime The developed theory will be applied to calculated hydrodynamic parameters of two important classes of gas-borne clusters which can be encountered in nature: chains of spheres and random fractal-like aggregates. Fractal-like aggregates of N spherules with equal radius a are characterized by the scaling dependence µ N = kf

RG a

¶D f ,

[45]

where D f and k f are the aggregate fractal dimension and prefactor, respectively, and RG is the gyration radius, which can be defined as RG2 =

N 1 X R2 + a2 N i=1 i

[46]

where Ri is the distance between the aggregate center of mass and the sphere center. The gyration radius [46] has a geometrical meaning of the mean square distance between the aggregate

DRAG AND TORQUE ON CLUSTERS OF N SPHERES

189

where K is the translation tensor for a given aggregate, the components of which can be easily found by performing calculations of the drag force for three independent (e.g., perpendicular) directions of the translation velocity. An average parameter, characterizing the aggregate mobility in the flow, is the aggregate hydrodynamic radius RH µ ¶ 1 1 1 1 = 2π + + , [48] RH K1 K2 K3

FIG. 2. Distribution of the magnitude of hydrodynamic drag acting on individual spheres in a fractal-like aggregate of 100 spheres with fractal dimension 1.8.

center of mass and all points on the particle surfaces. This definition has an advantage that the gyration radius of one sphere coincides with its geometric radius, while for larger N it quickly converges to the usual definition of the aggregate gyration radius, which neglects the term a 2 on the right-hand side of Eq. [46]. Intensity of gray in Fig. 2 shows the calculated distribution of the hydrodynamic force F j (in multiples of 6πµaU ) acting on spherules in a random aggregate of N = 100 particles with fractal dimension 1.8 and prefactor 2.3, generated numerically using a modified cluster–cluster aggregation algorithm of Ref. (20). Also shown is the flow velocity field (relative to the aggregate) in one of the planes passing through the aggregate center of mass parallel to the translation velocity U. In this case values of the drag force on spheres situated along the aggregate perimeter may exceed those for particles in the aggregate core by a factor of about 10. This significant shielding effect is characteristic also for other elliptic problems, such as electrostatic analysis and continuum regime heat and mass transfer to/from aggregated particles (21). Figure 3 shows the distribution of the hydrodynamic torque, calculated (in multiples of a · 6πµaU ) for the same aggregate under the same flow conditions. The shielding effect here is still more pronounced than for the case of drag force (Fig. 2), and the torque on perimeter-located spheres can be more than 20 times higher than that for core spheres. In general, the total drag and torque acting on aggregates are strongly dependent on the direction of the flow. For “stiff” aggregates translating with velocity U in respect to the fluid, the drag force is given by (2) F = −µK · U,

[47]

where K 1 , K 2 , and K 3 are the eigenvalues of the translation tensor K. The hydrodynamic radius of relatively dense aggregates can be close to the gyration radius RG . For example, if a is the sphere radius, the gyration radius of the cluster shown in Figs. 2 and 3 equals 8.325a, while the hydrodynamic radius calculated using the described theory is 8.186a. This agrees with assumptions made by Schmidt-Ott et al. (22), as well as with calculations of Ref. (23) using Kirkwood–Riseman theory for very large random aggregates for aggregates with a fractal dimension about 1.8. Figure 4 shows the dependencies of the average ratio RH /RG calculated for fractal-like aggregates with different spherule numbers N , fractal dimensions D f , and prefactor k f . Each symbol represents the result of averaging of RH /RG over 30 simulated fractal geometries. The aggregate fractal dimension is the major factor determining the slope of the curves RH /RG (N ), which are close to straight lines in the log–log scales, whereas the prefactor k f is mainly responsible for the shift of the lines in the vertical direction. For dense aggregates with relatively high D f and k f values, the ratio RH /RG weakly depends on the spherule number, and the hydrodynamic radius may exceed the

FIG. 3. Distribution of the magnitude of hydrodynamic torque acting on individual spheres in a fractal-like aggregate of 100 spheres with fractal dimension 1.8.

190

A. V. FILIPPOV

FIG. 4. The ratio of the hydrodynamic and gyration radius as function of the aggregate size, calculated for different values of the aggregate fractal dimension D f and prefactor k f .

gyration radius. In the limiting case of spherically symmetric 3dimensional aggregate consisting of a large number of densely packed spherules the gyration radius is close to (3/5)1/2 of both the geometric radius and hydrodynamic radius (which are equal in this case), so that the ratio RH /RG is close to 1.29 for large N . For aggregates with fractal dimension below 1.7 ÷ 1.5 the ratio RH /RG is decreasing with N , and the slope of the lines in Fig. 4 is steeper for aggregates with smaller fractal dimension. In the limiting case of straight chains, which are onedimensional aggregates, the slender body theory predicts the following asymptotic behavior for ratio of hydrodynamic and gyration radius (24) at N → ∞: √ µ ¶ 1 RH 3 +O = . [49] RG ln(2N ) ln2 (2N )

or, equivalently, in the form αv Kn, where Kn is the Knudsen number based on the radius of the considered spherule. The value of the velocity slip coefficient depends on the momentum accommodation of gas molecules at the surface (26) with typical values of order 1 ÷ 1.5. In various aerosol applications, the boundary condition [4] is applicable at least for values of αv Kn below 0.1 and, in some cases, up to 0.25 (27). The limiting case αv Kn → ∞ corresponds to ideal slip conditions on the surface, and the developed theory can be used in modeling bubble dynamics in liquids. However, because we are considering dynamics of gas-borne rigid clusters of touching spheres, the presented results will be mainly restricted to the region of αv Kn not higher than 0.3. As a result of the velocity slip effect, the total drag force and torque acting on the cluster must be corrected. Figure 5 shows the size dependence of the hydrodynamic radius calculated for straight chains of spheres and fractal-like clusters of equal spheres in the continuum (solid lines) and slip regime with αv Kn = 0.3 (dashed lines). In this and other examples in this section, the data for fractal-like aggregates are obtained by numerical generating the random clusters with fractal dimension 1.8 and prefactor 1.2 and statistical averaging the calculation results performed for 30 different realizations for each cluster size. The resulting statistical spread of the calculated data was not high and did not exceed a few percent of the corresponding mean values. Qualitatively, the velocity slip effect is characterized by the slip correction factor, which is equal to the ratio of hydrodynamic radius at nonzero Knudsen number to its continuum limit value: C(N , Kn) =

RH (N , Kn) . RH (N , 0)

[50]

Qualitatively, the presented results agree with conclusions of Rogak and Flagan (25), who used the approximate Kirkwood– Riseman theory to analyze the relation between the hydrodynamic and gyration radius of various numerically generated fractal-like aggregates with different fractal dimensions, but did not specify and study the influence of the fractal prefactor. 4.3. Particle Clusters in the Slip Regime Calculations of the aggregate dynamics in the slip regime of small Knudsen numbers are important, for example, in describing the deposition of submicrometer fibrous and random organic and inorganic clusters in lungs and filters. The presented theoretical description of the near-continuum regime is based on the slip boundary condition [4], where the velocity slip coefficient and the mean-free path enter only in the form of product lαv

FIG. 5. Size dependence of the hydrodynamic radius for straight chains of spheres and fractal-like clusters with fractal dimension 1.8 and prefactor 1.2 in the continuum regime (Kn = 0, solid lines) and slip regime (αKn = 0.3, dashed lines).

DRAG AND TORQUE ON CLUSTERS OF N SPHERES

191

In the case of one sphere (N = 1), the slip correction factor is given by a classical solution of the Stokes problem with velocity slip boundary conditions [3] (14): C(1, Kn) =

1 + 2αv Kn . 1 + 3αv Kn

[51]

The calculated dependencies of the slip correction factor of fractal-like aggregates containing from 5 to 35 spheres on the product αv Kn are shown in Fig. 6. Qualitatively, the behavior of all curves is similar and looks as if the solution [51] for one sphere, shown by a dotted line, were stretched along the αv Kn axis. The hydrodynamic radius RH (N , Kn) is equal to the radius of the sphere having the same average mobility in the continuum limit as a given aggregate at Knudsen number Kn. This parameter has sense in the considered near-continuum regime. A “true” hydrodynamic or mobility radius RM is the radius of the sphere having the same average mobility as a given aggregate under the same conditions (28). Taking into account Eq. [50] for the slip correction factor, the following equation for the mobility radius RM can be obtained C(1, KnM )RM = C(N , Kn)RH (N , 0) KnM =

λ a = Kn , RM RM

[52] [53]

where KnM is the “adjusted” Knudsen number based on the mobility radius. In the continuum limit, the mobility radius coincides with the aggregate hydrodynamic radius RH (N , 0), but in general case, it is a function of the Knudsen number and the aggregate shape and size.

FIG. 7. Slip-regime independence of the mobility radius of straight chains of spheres, shown normalized with the chain gyration radius, on the product αKn.

Using the described theory, it is possible to calculate the slip correction factor and the continuum-regime hydrodynamic radius RH (N , 0) on the right-hand side of Eq. [52]. Then, using inversion of the simple analytical solution [51] for one sphere, the slip-regime dependence of the mobility radius on the parameters N and Kn for aggregates of a given class can be found. The results of such calculations for straight chains of up to 35 spheres of different size are shown in Fig. 7, where the mobility radius RM , normalized by the aggregate gyration radius RG , is plotted as function of the product αv Kn. Remarkably, the mobility radius is independent on the Knudsen number throughout the whole slip region. Similar calculations for random fractal-like aggregates with fractal dimension 1.8 yield still more interesting features shown in Fig. 8. Due to insensitivity of the ratio RM /RG in the continuum regime, which can be noticed in Fig. 4, all the curves nearly collapse into one already for aggregates of more than 5 spheres, although some weak dependence of the ratio Rad /RG on the Knudsen number is now observed. These properties can be used in developing and testing of new interpolation formulas for the mobility parameters in the intermediate regime of moderate Knudsen numbers. 4.4. Comparison with Published Theoretical and Experimental Results

FIG. 6. Slip correction for a single sphere and fractal-like clusters of different sizes as function of the product of the velocity slip coefficient and the Knudsen number. The aggregate fractal dimension is 1.8; the prefactor equals 1.2.

Very few computational data on hydrodynamic forces on aggregates are available for slip-regime conditions of nonzero, but small Knudsen numbers Kn. Table 2 shows the dimensionless drag force on each of two equal spheres moving along their line of centers with equal velocities, calculated using the method of reflections (27) and Multipole Expansion Method (MEM) calculations described in this work. The data, calculated for two distances h between spheres and different values of the product

192

A. V. FILIPPOV

TABLE 3 Dimensionless Drag Fi /(6πµai v ∞ ) Acting on Each of Two Spheres in a Flow Perpendicular to the Symmetry Axis as Function of the Ratio of Their Sizes, Calculated by Davis (29) and by the MEM Method Described in This Work

FIG. 8. Dependence of the ratio of the aggregate mobility and gyration radius on the product αKn, calculated for fractal-like aggregates with fractal dimension 1.8 and prefactor 1.2.

F2 /(6π µa2 v∞ )

Dynamic shape factor K s

Ratio a2 /a1

Ref. (29)

MEM

Ref. (29)

MEM

Ref. (29)

MEM

1 2 5 10

0.7331 0.576 0.3516 0.2129

0.7331 0.576 0.3515 0.2102

0.7331 0.8523 0.95 0.9825

0.7331 0.8523 0.9501 0.9827

1.1637 1.0964 1.0176 1.0035

1.1637 1.0964 1.0177 1.0033

calculated for each of the spheres are normalized by the drag force experienced by an isolated sphere of the same radius in the same flow, given by Eq. [54] with C = 1. The aggregate dynamic shape factor K s is F 6π µaeq v∞ ! 13 Ã N X 3 = (a j ) ,

Ks =

αv Kn, are normalized by the drag force experienced by an isolated sphere in the same flow:

aeq

[55]

[56]

j=1

F1 = 6πµav∞ C(1, Kn).

[54]

The agreement between predictions of two different methods is very good. In contrast to the slip-regime situations, numerous data for the case Kn = 0 enable a detailed testing of the proposed method in the continuum limit. For example, Table 3 shows the drag force acting on two spheres of different radii a1 and a2 in their translation perpendicular to the symmetry axis of such two-particle aggregate, when the separation between spheres is 10% of the smaller radius. The results shown in Table 3 are based on the paper of Davis (29) who used the solution in bispherical coordinates and on the MEM predictions of this work. The data TABLE 2 Correction Factor for Two Equal Spheres Moving along Their Line of Centers with Equal Velocities Calculated for Two Distances h between Spheres and Different Values of the Product αv Kn h/a = 0.005

F1 /(6π µa1 v∞ )

h/a = 1.08616

α Kn

Meth. of reflections

MEM

Meth. of reflections

MEM

0 0.1 0.2 0.3 106

0.64572 0.65531 0.66185 0.66652 0.69388

0.64572 0.65531 0.66185 0.66652 0.69388

0.70245 0.71441 0.72447 0.72822 0.7618

0.70245 0.71442 0.72448 0.72822 0.7618

Note. Shown are the calculations using the method of reflections (27) and MEM predictions of this work.

where aeq is the radius of the equal volume sphere. For the ratio a2 /a1 = 10, an about 1.4% discrepancy for the drag on a smaller particle results, probably, from relatively slow convergence observed for this property in both methods. However, since the contribution of the smaller particle to the total drag in this case is small, the total drag and dynamic shape factor converge very well and agree to 0.01% in all considered cases. Obviously, most of the published experimental and theoretical data deals with aggregates of equal diameter spheres. Lasso and Weidman (30) performed experiments on sedimentation of closed packed aggregates of spheres in a silicone oil and determined the drag force exerted by the fluid. Some of these structures are shown in Fig. 9, while the measured dynamic shape factors are compared in Table 4 with the predictions of the MEM theory described in this work, as well as with the calculations

TABLE 4 Dynamic Shape Factor K s for Hexagonal-Closed-Packed Agglomerates of Spheres (Shown in Fig. 8), Measured by Lasso and Weidman (30), Calculated by Cihocki and Hinsen (31), and Predicted by the MEM Method of This Work N

Experiment (30)

Theory (31)

MEM

3 5 13 57

1.195 1.095 1.143 1.131

1.1959 1.0783 1.1356 1.144

1.1958 1.0778 1.1342 1.1372

193

DRAG AND TORQUE ON CLUSTERS OF N SPHERES

TABLE 6 Dynamic Shape Factor Ts for Straight Chains Rotating Parallel and Perpendicular to Major Axis, Calculated by the MEM Method of This Work, and by the BEM Method (7) Rotation perpendicular to major axis

FIG. 9. Hexagonal closed-packed sphere aggregates and direction of their translation in experiments of Lasso and Weidman (30).

of Cichocki and Hinsen (31). All three sets of data agree very well. Straight chains of spheres are, probably, the most theoretically well-investigated aggregate arrangements, allowing effective application of the Boundary Element Method (BEM). The MEM calculations of dynamic shape factor for chain motion perpendicular to the major axis are shown in Table 5 along with the data calculated by the BEM method (7) as well as with the calculations of Durlofsky et al. (8) and experiments of Kasper et al. (32). The agreement is very good, with MEM values situated between predictions of the BEM and Durlofsky et al. (8) theory. No experimental data are available for hydrodynamic torque acting on aggregates, and here accurate theoretical calculation is particularly important. A dimensionless shape factor for torque can be determined as Ts =

T , 8πµaωN L a2

[57]

where L is a characteristic dimension of the aggregate in respect TABLE 5 Dynamic Shape Factor K s for Straight Chains of Spheres Moving Perpendicular to Their Major Axes, Calculated by the MEM Method of This Work, by BEM Method (7), by Multipole Force Expansion (MFE) (8), and Measured in Experiments of Kasper et al. (32) N

MEM

BEM (7)

MFE (8)

Experiment (32)

2 3 4 5 6 10 15 20 25 30

1.15 1.276 1.386 1.485 1.576 1.889 2.209 2.484 2.729 2.953

1.147 1.272 1.381 1.48 1.567 1.875 2.194 2.468 2.713 2.935

n.a. n.a. n.a. 1.485 n.a. 1.894 2.226 2.505 2.762 2.993

1.14 1.27 1.35 1.48 1.52 1.86 2.19 2.4 2.61 2.8

Rotation parallel to major axis

N

This work (MEM)

Geller et al. (7)

This work (MEM)

Geller et al. (7)

2 3 4 5 6 7 9 10 15 20 30

0.468 0.338 0.278 0.243 0.219 0.202 0.179 0.17 0.143 0.128 0.111

0.468 0.338 0.278 0.243 0.221 0.2 0.177 0.169 0.142 0.127 0.11

0.902 0.865 0.845 0.833 0.825 0.819 0.811 0.808 0.799 0.795 0.79

0.906 0.868 0.853 0.841 0.832 0.777 0.769 0.766 0.756 0.751 0.747

to the rotation axis and ω is its angular velocity. Table 6 shows the values of the shape factor Ts for straight chains of spheres, calculated by the BEM method (7) and by the MEM method of this work. For rotation perpendicular to major axis L a = N a is assumed, while L a = a is taken for the case when the rotation and chain major axis coincide. The agreement within 1% is found in all cases beside the case of rotation parallel to the major axis for chains with N > 7, where the discrepancy is relatively high. However, an abrupt drop of the BEM-calculated shape factor at N = 7 indicates to a possibility of systematic error in BEM data for higher N . Smooth behavior of MEM predictions and their invariance to increasing the truncation level L observed in control calculations imply that the described theory provides adequate results in these situations too. In the case of more complicated, fractal-like aggregate morphology, most of the published theoretical and experimental data do not specify the value of the aggregate fractal prefactor k f , which makes a direct quantitive comparison difficult. In experiments, the ratio of hydrodynamic radius RH and gyration radius RG of fractal-like aggregates has been estimated experimentally, from static and dynamic light scattering measurements with suspensions of colloidal (33) and aerosol clusters (34). For example, diamonds in Fig. 10 show the ratio RH /RG based on the data obtained by Wiltzius (33) for colloidal aggregates of relatively small sizes of up to some hundred monomers. Calculations for numerically generated fractal-like aggregates with fractal dimension 1.8 in the same size range are shown by triangles (for prefactor k f = 1.3) and circles (k f = 2.3). These prefactor values have been reported by different authors for soot aggregates (35) with a morphology similar to that of colloidal clusters. The difference between the two theoretical predictions is not appreciably higher than the spread of the experimental data. Another uncertainty factor is the polydispersity of real aggregates, because of which the measurements can yield only

194

A. V. FILIPPOV

nonzero, but small Knudsen numbers, often realizable in practice, for example, in the deposition of submicrometer air-borne biological or inorganic clusters in human lungs. Another challenging problem is prediction of soot dynamic properties in a high-pressure environment, typical for diesel and aircraft engines, where the molecular mean-free path can become smaller than the primary particle size. The computational efficiency of the theory described can be a valuable property in view of possible future applications, including predictions of the thermophoretic, gravitational, and inertial deposition of flocculated spherical particles in chemical industry, combustion science, and environmental engineering.

APPENDIX FIG. 10. Ratio of hydrodynamic and gyration radius of colloidal clusters measured by Wiltzius (33) (diamonds) and calculated in this work for numerically generated fractal-like aggregates with fractal dimension 1.8 and prefactor 1.3 (triangles) and 2.3 (circles).

population-averaged values of both radii, RH and RG . In spite of these complications, the aggreement between the described theory and experiment is reasonably good.

The coefficients for the terms under summation signs in Eqs. [24]–[26] are ij

n ij C ai klmn n ij = 2 [n − 1 − 2(n 2 − 1)α Kni ]Cklmn ai

Dklmn = ij

G klmn ij

E klmn = 5. CONCLUSIONS

[A1] [A2]

[2(l + 1) − n(l − 2)]ai2 − n(l − 2)Ri2j 2l(2l − 1)ai

ij

Cklmn

ij

Hydrodynamic properties of arbitrary clusters of N spheres were studied using the multipole expansion technique developed in light scattering problems (15) and the general solution of the Stokes problem (14) in terms of solid spherical harmonics. The linear system of equations for the coefficients of the multipole expansions was found by direct transformation of spherical coordinates and spherical harmonics from one center to another, and expansion of the boundary conditions for each of the N spheres, as in Lamb’s theory. The method was applied to find dependence of the aggregate mobility radius on the number of particles for the gas-borne chain and fractal-like aggregates. In another application of the developed theory, the dependence of the cluster mobility radius on the Knudsen number was determined under the slipregime conditions. These results can be used in development of the interpolation formulas for the aggregate–gas dynamic interactions in the transition regime of moderate Knudsen numbers. Direct comparison with available experimental and theoretical data on hydrodynamic forces and torques acting on clusters of rigid spheres has shown a very good agreement and rapid convergence of results. This enables performing calculations for dense aggregates with N of up to some hundred using average personal computers and workstations. Because no additional approximations are made, the developed method can be considered as a straightforward extension of the Lamb’s general solution for Stokes flow (14) to the case of N arbitrary spheres. The general slip boundary conditions on the sphere surfaces allow application of the method to cases of

− (m + n + 1)3+ ai ζi j Cklm(n+1) − (n − m)3−

ζi j i j ηi j i j C − 3− Ckl(m−1)(n−1) ai klm(n−1) ai

ij

+ 3+ ηi j ai Ckl(m−1)(n+1) ij

− (n + m + 2)(n + m + 1)3+ ξi j ai Ckl(m+1)(n+1) + (n − m − 1)(n − m)3− ij

ξi j i j C ai kl(m+1)(n−1)

[A3]

−i £ ij (n + m + 1)(n − m) ξi j Ckl(m+1)n ai ¤ ij ij − mζi j Cklmn + ηi j Ckl(m−1)n [A4] · 2(l + 1) − n(l − 2) i j Cklmn = [n + 1 − 2n(n + 2)αKni ] 2l(2l − 1)

Fklmn =

ij

Hklmn

ij

ij

− (m + n + 1)3+ ζi j Cklm(n+1) + 3+ ηi j Ckl(m−1)(n+1) ¸ ij − (n + m + 2)(n + m + 1)3+ ξi j Ckl(m+1)(n+1) · n(l − 2)Ri2j i j 1 2 C − 2 [n − 1 − 2(n − 1)α Kni ] 2l(2l − 1)ai klmn ai ij

ij

+ (n − m)3− ζi j Cklm(n−1) + 3− ηi j Ckl(m−1)(n−1) ¸ ij − (n − m − 1)(n − m)3− ξi j Ckl(m+1)(n−1) [A5]

DRAG AND TORQUE ON CLUSTERS OF N SPHERES ij

L klmn = ij

1 ij [n − 1 − 2(n 2 − 1)α Kni ]Fklmn ai

1 ij [1 + α(n − 1)Kni ]Fklmn l £ 1 ij = [1 + α(n − 1)Kni ] n(n + 1)Cklmn ai

Mklmn = ij

Nklmn

[A6] [A7]

8. 9.

ij

+ ζi j (m + n + 1)nCklm(n+1) ij

+ ξi j n(n + m + 1)(n + m + 2)Ckl(m+1)(n+1) ¤ ij − ηi j nCkl(m−1)(n+1) , [A8] where 3+ and 3− are the combinations 3+ =

(n + 1)(l − 2) − (l + 1) ; (2n + 3)l(2l − 1)

(n − 1)(l − 2) − (l + 1) . 3− = (2n − 1)l(2l − 1)

3. 4. 5. 6. 7.

[A9]

According to definition of spherical harmonics [18] and Eq. [19], all transformation coefficients Cklmn are assumed to be zero, if |k| > l or if |m| > n.

10. 11. 12. 13. 14. 15. 16. 17. 18. 19. 20. 21. 22. 23. 24. 25.

ACKNOWLEDGMENTS 26. The author is thankful to D. E. Rosner, D. W. Mackowski, and A. S. Sangani for insightful discussions and comments. The financial support of NASA Grant NAG3-1951, AFOSR Grant 97-1-0266, and NSF Grant CTS-9980747 is gratefully acknowledged.

REFERENCES 1. Smoluchowski, M., Bull. Int. Acad. Sci. Cracovie. 1A, 28 (1911). 2. Happel, J., and Brenner, H. “Low Reynolds Number Hydrodynamics.” Kluwer, Dordrecht, 1983.

27. 28. 29. 30. 31. 32. 33. 34. 35.

195

Debye, P., and Bueche, A. M., J. Chem. Phys. 16, 573 (1948). Rosner, D. E., and Tandon, P., AIChE J. 40, 1167 (1994). Kirkwood, J. G., and Riseman, J., J. Chem. Phys. 16, 565 (1948). Rotne, J., and Prager, S., J. Chem. Phys. 50, 4831 (1969). Geller, A. S., Mondy, L. A., Rader, D. J., and Ingber, M. S., J. Aerosol Sci. 24, 597 (1993). Durlofsky, L., Brady, J. F., and Bossis, G., J. Fluid Mech. 180, 21 (1987). Ladyzhenskaya, O. A., “The Mathematical Theory of Viscous Incompressible Flow.” Gordon and Breach, New York, 1963. Jeffrey, D. J., and Onishi, Y., J. Fluid Mech. 139, 261 (1984). Ladd, A. J. C., J. Chem. Phys. 90, 1149 (1988). Cichocki, B., Felderhof, B. U., Hinsen, K., Wajnryb, E., and BlÃawzdziewicz, J., J. Chem. Phys. 100, 3780 (1994). Hassonjee, Q., Ganatos, P., and Pfeiffer, R., J. Fluid Mech. 197, 1 (1988). Lamb, H., “Hydrodynamics.” Dover, New York, 1945. Mackowski, D. W., J. Colloid Interface Sci. 140, 1139 (1990). Mo, G., and Sangani, A. S., Phys. Fluids 6, 1637 (1994). Mackowski, D. W., Appl. Opt. 34, 3535 (1995). Arfken, G., “Mathematical Methods for Physicists,” 3rd ed. Academic Press, Orlando, FL, 1985. Sangani, A. S., and Mo, G., Phys. Fluids 6, 1653 (1974). Thouy, R., and Jullien, R., J. Phys. 6I, 1365 (1996). Filippov, A. V., Zurita, M., and Rosner, D. E., J. Colloid Interface Sci. 219, 216 (2000). Schmidt-Ott, A., Baltensperger, U., G¨aggeler, H. W., and Jost, D. T., J. Aerosol Sci. 21, 711 (1990). Meakin, P., Chen, Z.-Y., and Deutch, J. M., J. Chem. Phys. 80, 2982 (1984). Batchelor, G. K., J. Fluid Mech. 44, 119 (1970). Rogak, S. N., and Flagan, R. C., J. Colloid Interface Sci. 134, 206 (1990). Rosner, D. E., and Papadopoulos, D. H., Ind. Eng. Chem. Res. 35, 3210 (1996). Reed, L. D., and Morrison, F. A., J. Aerosol Sci. 5, 175 (1974). Rogak, S. N., Flagan, R. C., and Nguyen, H. V., Aerosol Sci. Technol. 18, 25 (1993). Davis, M. H., Chem. Eng. Sci. 24, 1769 (1969). Lasso, I. A., and Weidman, P. D., Phys. Fluids 29, 3921 (1986). Cichocki, B., and Hinsen, K., Phys. Fluids 7, 285 (1994). Kasper, G., Niida, T., and Yang, M., J. Aerosol Sci. 16, 535 (1985). Wiltzius, P., Phys. Rev. Lett. 58, 710 (1987) Wang, G. M., and Sorensen, C. M., Phys. Rev. E 60, 3036 (1999). ¨ O., ¨ Xing, Y., and Rosner, D. E., Langmuir 11, 4848 (1995). K¨oyl¨u, U.