Stokes flow in spheroidal particle-in-cell models with rappel and kuwabara boundary conditions

Int. J. Engng Sci. Vol. 33, No. 10, pp. 1465-1490, 1995

Pergamon 0020-7225(95)00010--0

Copyright ~) 1995 Elsevier Science Ltd Printed in Great Britain. All rights reserved 0020-7225/95 $9.50+ 0.00

STOKES FLOW IN SPHEROIDAL PARTICLE-IN-CELL MODELS WITH HAPPEL AND KUWABARA BOUNDARY CONDITIONS G. DASSIOS, M. H A D J I N I C O L A O U , F. A. COUTELIERIS and A. C. PAYATAKES Institute of Chemical Engineering and High Temperature Chemical Processes, and Department of Chemical Engineering, University of Patras, GR 265 00 Patras, Greece Abstract--Particle-in-cell models are useful in the development of simple but reliable analytical expressions for heat and mass transfer in swarms of particles. Most such models consider spherical particles. Here the creeping flow through a swarm of spheroidal particles, that move with constant uniform velocity in the axial direction through an otherwise quiescent Newtonian fluid, is analyzed with a spheroid-in-cell model. The solid internal spheroid represents a particle of the swarm. The external spheroid contains the spheroidal particle and the amount of fluid required to match the fluid volume fraction of the swarm. The boundary conditions on the (conceptual) external spheroidal surface are similar to those of the sphere-in-cell Happel model [1], namely, nil normal velocity component and shear stress. The stream function is obtained in series form using the recently developed method of semiseparation of variables. It turns out that the first term of the series is sufficient for most engineering applications, so long as the aspect ratio of the spheroids remains within moderate bounds, say - 1 / 5 < a 3 < - 5 . Analytical expressions for the streamfunction, the velocity components, the vorticity, the drag force acting on each particle, and the permeability of the swarm are obtained. Representative results are presented in graph form and they are compared with those obtained using Kuwabara-type boundary conditions. The Happel formulation is slightly superior because it leads to a particle-in-cell that is self sufficient in mechanical energy.

1. I N T R O D U C T I O N

Flow through a swarm of particles arises in many processes of practical importance, such as fluidization, sedimentation, flow in packed beds, etc. Determination of the flow field in each case is important for two reasons. First, the flow field in itself can be used to determine important engineering quantities, such as the drag force exerted on each particle, the macroscopic pressure gradient, the permeability of the swarm, etc. Second, the flow field is a necessary basis for the analysis of a host of transport processes, such as mass transfer with adsorption or reaction, heat transfer, fine suspended particle motion and deposition, etc. The flow in many of these processes is creeping, as the Reynolds number is smaller than unity. In engineering type analyses it is not usually necessary to have a detailed solution of the flow field over the entire swarm of particles, taking into account the exact positions of the particles (such solutions are cumbersome to use and unnecessarily detailed for many engineering applications). Rather, it is sufficient to obtain a relatively simple analytical expression that takes into account the effects of the neighboring particles on the flow field around a single particle of the swarm. This, in turn, can be used to develop relatively simple, yet reliable, models for heat and mass transfer. This rationale has led to the development of particle-in-cell models. Uchida [2] proposed a cell model for a sedimenting swarm of particles. In his model the spherical particle is surrounded by a fluid envelope with cubic outer boundary. The cubic shape of this boundary offers the advantage that it is space filling, but the difference in geometry between the inner and outer boundaries leads to a three-dimensional flow problem. Uchida proposed a simplified flow solution which, however, is not in good agreement with data [3, p. 376]. Brenner [4] developed an accurate solution of creeping flow for this geometry using a collocation method [3, p. 376]. Happel [1] and Kuwabara [5] proposed cell models in which both particle and outer envelope are spherical. "['his formulation has the significant advantage that it leads to an axially symmetric flow that has a simple analytical solution in closed form, and thus can be used 1465

1466

G. DASSIOS et al.

readily for heat and mass transfer calculations. On the other hand, it has the disadvantage that the outer envelope is not space filling, a difficulty which must be dealt with, when one tries to pass from the single unit cell to the assemblage of particles. (This point will be discussed further, in the section dealing with the determination of the permeability of the swarm.) The difference between the Happel and Kuwabara formulations lies in the boundary conditions. The Happel model assumes that the inner sphere--while at the center--moves with a constant velocity t7 along the axis. Assuming pseudosteady state, the following boundary conditions are imposed: non-slip flow on the inner sphere, nil radial velocity and nil shear stress on the outer envelope. The Kuwabara model assumes that the inner sphere is stationary and that fluid passes through the unit cell. The following boundary conditions are imposed: nil radial and tangential velocity on the inner sphere, velocity with axial component equal to a constant approach velocity t~ on the outer envelope, and nil vorticity on the outer envelope. Both formulations give essentially the same velocity fields (with the appropriate change of frame of reference) and approximately equal drag forces, Tien [6]. However, the Happel formulation has a significant advantage in that it does not require an exchange of mechanical energy between the cell and the environment. On the contrary, the Kuwabara model requires a small but discernible exchange of mechanical energy with the environment. The mechanical power given by the sphere to the fluid is not all consumed by viscous dissipation in the fluid layer. Rather, a small part is given to the environment, because on the external boundary the product of the shear stress and the tangential velocity is everywhere negative, except at the two polar points where it is nil [3, p. 390]. Neale and Nader [7] improved the formulation of Happel and Kuwabara by considering that the unit cell under consideration is embedded in an unbounded, continuous, homogeneous and isotropic permeable medium, which has the same permeability with that of the swarm of spheres. The flow in the exterior permeable medium is governed by Brinkman's equation. Flow continuity, and pressure and shear stress continuity are imposed at the interface between the fluid layer and the surrounding permeable continuum. An analytical solution in closed form is obtained for this model, too. Permeabilities predicted by this model are in good agreement with experimental data over a wide range of porosity values [6]. The Happel and Kuwabara models also give good agreement, but somewhat less so than the Neale and Nadar model. Pfeffer and Happel [8] and Pfeffer [9] used Happel's unit cell to model mass transfer to an adsorbing sphere in a swarm, for the case of low Reynolds and high Peclet number. Using the Levich [10] method for Pe >> 1, they obtained an expression of the form Sho = 0.997 f(3")Pe 1/3, where Sho is the overall Sherwood number and f(3') is a simple analytical function of the solid volume fraction 3'. Tardos et al. [11] obtained similar expressions for the low-Re/high-Pe case, based on the Kuwabara, and the Neale and Nader models. They concluded that all three unit cell models give comparable predictions, but that the Neale and Nader version gives somewhat better agreement with data. Epstein and Masliyah [12] proposed a useful generalization of the sphere-in-cell model by considering a spheroid-in-cell model for swarms of spheroidal particles. However, they had to solve the creeping flow problem numerically (a fact that has impeded the use of this flow model for heat and mass transfer calculations). The difficulty encountered in this case (as well as in other cases involving Stokes flow in spheroidal coordinates) is that the governing equation in terms of the stream function, E4~b = 0, is n o t separable (whereas it is separable in Cartesian, cylindrical and spherical coordinates, and R-separable in bispherical coordinates). This difficulty was resolved recently in Dassios et al. [13] with the introduction of the method of semiseparation of variables, which was then used to obtain an analytical solution of Stokes flow in a spheroid-in-cell with Kuwabara-type (see above) boundary conditions. The leading term of this series solution turns out to be satisfactory for engineering applications so long as the solid volume fraction 3' is not too large (say, 3' < - 0 . 3 ) , and the aspect ratio in the range, say, 0.20 < a 3 < 5. Couterlieris et al. [14] have already used this flow solution to model mass transfer

Stokes flow in spheroidal particle-in-ceU models

1467

to an adsorbing spheroid in a swarm of similar spheroids, in the case of low-Re/high-Pe. They obtained for the overall Sherwood number an expression of the form Sho = 0.997 g(y, a3)Pe u3, where g(Y, a3) is a given function of the solid volume fraction 7 and the aspect ratio a3, and give an extensive: table of g values. In the present work the solution to the Stokes flow problem in a spheroid-in-cell with Happel-type boundary conditions is obtained. The incentive for this is that the Happel-type boundary conditions are more compatible with the physics of flow in a swarm, since they ensure that each unit cell is self-sufficient in mechanical energy. In contrast, each unit cell with Kuwabara-type boundary conditions is required to give work (to an unspecified sink in the environment) at a small, but not entirely negligible rate, especially for large values of the solid volume fraction. Analytical expressions for the drag force and the permeability are also obtained, and the leading term of the solution is used to give results in graph form for wide ranges of values of the solid volume fraction and the aspect ratio. Both types of spheroids, prolate and oblate are considered. The corresponding results based on the Kuwabara-type boundary conditions are also given for comparison. An error in Kuwabara's expression for the drag force on a sphere-in-cell is detected and the corrected expression (which turns out to be much simpler than the erroneous one) is given in Appendix C.

2. F O R M U L A T I O N

OF THE PROBLEM:

HAPPEL

T Y P E BCs

Spheroidal particles can be either prolate or oblate. Since by using a simple transformation, as we will see later on, one obtains the results for the oblate spheroid from those for the prolate, we begin by treating the case of prolate spheroids only. Let the long semiaxis of the solid spheroid be denoted by a3 and the short semiaxis by al. The surface of the solid spheroid Ss is described by a-----~l + a'-~= 1,

(:~x, g2, $3)

on

Ss

(1)

where ($,, $2, x3) are the dimensional Cartesian coordinates, Fig. 1. Using the short semiaxis t71 as characteristic length, the following dimensionless quantities are defined: a3 --=~3 ~,

al = 1,

£i xi = ~-

c - e -- +X/=aa~-a~ - +Vaa 2 - 1,

(2)

where ~ is the semifocal length. Note that a 3 is the aspect ratio of the spheroid. Because of the geometry of the particles, it is more appropriate to introduce prolate spheroidal coordinates (77, 0, ~b), Fig. 2 [15]. The dimensional Cartesian coordinates (x~, xz, x3) are related to the prolate spheroidal ones (~/, 0, ~b) through the equations x~ = c sinh 77 sin 0 cos ~b] / x2 = c sinh 77 sin 0 sin ~b ~ x3 = c cosh ~7cos 0 J

0 -< 77 < ~ 0 -< 0 -< :r . 0 -< ~b < 27r

(3)

The equation describing any spheroidal surface S,~ (fixed 77) is +

)1, c 2 sinh 2 71 c 2 cos h2 ~7

(Xl, X2, X3) on

S,.

(4)

The solid surface, S~, in the system of spheroidal coordinates is obtained by setting )7 = a -= sinh -1 1. C

(5)

G. DASSIOS et al.

1468

/

/

/

/

/

._Ix_ \

Jb3

xx \

\

I I

i

t---o ....

\\\

'

X2

//

\x',,,q..../'/

Fig. 1. Solid spheroid in a spheroidal fluid layer (spheroid-in-cell).

According to the spheroid-in-cell model, the solid surface S, is contained in a conceptual spheroidal surface Sa. Sa is confocal with S, and its size is chosen so that the solid volume fraction of the cell is equal to that of the swarm. To obtain S~ we set 7/=/3 =- sinh_ 1 __bl c

(6)

where bl is the short dimensionless semiaxis of the external spheroidal surface S~ (bl =/~i/al). The value of bl is obtained from b~ = + ~ -

c2

(7)

where b3 is the long dimensionless semiaxis (b3 = b3/t~l), and is determined as follows. Let the solid volume fraction of the cell be denoted by 3/ ( = 1 - e). Now, the volume of the cell is (4g/3)~3b2b3 and the volume of the solid is (4g/3)53a3 (recall that a~ = 1). Consequently, in order for the solid volume fraction of the cell to be equal to 3' we must have 3/b~b3 = a3.

(8)

Equations (7) and (8) give b3 -

c2b3 -

a-23 = 0

(9)

3/ which has only one real positive root and is obtained as follows. Let c g be the value c2 = 3\~_~/

1 3/1" = 3COS[iCOS-'

(10)

Stokes flow in spheroidal particle-in-ceU models

1469

I.

X2

.(p=const.

i

v - - ° 1

Fig. 2. Spheroidal coordinates.

If c 2 -- c 2, then b3

r a 3 + .,/(a3~ 2

a3

.~/(a3) 2

(11a)

If c 2 -->c 2, then 2V]c

[1 cos_l(3V~ a3~ ]

b~=-]--c°st]

\ 2-~~-] J

(11b)

Substituting b3 in equation (7) from equation (11), one obtains bl, while/3 is obtained from equation (6). For reasons of convenience the following independent variables are introduced [3]. r = cosh 7/,

~ = cos 0.

(12)

Then xt = c ~ v T ~ xl = c ~ X / ] - - ~ x3 = c ~ At this point we also set

cos th1 sin ~b ~ J

1 -< J~l -1<-~+1. 0 -< ~b < 21r

aa ~-1= ~ / l + c 2 • ~ = cosh a =--=c c2 r a = cosh/3

=--b3 C

(13)

(14) (15)

1470

G. DASSIOS

et al.

where ¢ is the eccentricity defined by C

= --. aa

(16)

Two confocal prolate spheroids are considered: the inner one, S~, at "t = "t~, is solid and is moving with velocity of magnitude a in the positive direction of the x3-axis, inside an otherwise quiescent fluid spheroidal layer which is confined by the outer spheroid S0, at "t = To. Following the formulation of Happel [1], the velocity component normal to So and the tangential stresses are assumed to vanish on S0. Assuming that the flow is creeping, (Re = 2a~tT//2 << 1), the equation of motion becomes E4~t = E2(E2d/) = 0

(17)

where ~b(=~,/t~ 2) is the dimensionless stream function, and the operator E 2 is given by E2 -

1 [ 02 ~,2 _0_2 ] . c2(z2 _ ¢2) ('t2 _ 1) ~ + (1 ) ~2

(18)

The boundary conditions are: BCI:

v,7 = (_x3" ¢/_) on

"t = "t,,

(19)

BC2:

Vo = (_~3" _0) on

t = "t,~

(20)

BC3:

v,~ = 0

on

"t = "to

(21)

BC4:

II~o = 0

on

"t = "to

(22)

where v, and vo are the 7/and 0 components of the dimensionless velocity p ( = f s / a ) , ¢1_ and 0 are the unit vectors of the spheroidal coordinates in the , / a n d 0 directions, respectively, ~3 is the unit vector of the Cartesian coordinates in the x3 direction, and II~o (=I'Ino~l//2a) is the dimensionless tangential stress. Note that for t ='t0 the unit normal vector r1 on the external spheroid S0, the Cartesian components of which are given by h- =

('t0 1~"~-~ ~ 2 COS (#, " t f l V T ~ sin ~b, X/-~ - 1¢) ~

'

(23)

coincides with the unit vector ~_. Equations (19) and (20) express the non-slip flow condition. Equation (21) implies that there is no flow across the boundary of the fluid envelope S0. Equation (22) expresses the assumption that the tangential stresses vanish on S0, in accordance with the argument advanced by Happel [1]. In order to express the BCs in terms of the stream function ~b, we use the relations 1

05

(24)

vn - - hoh~, O0 1

0~b

(25)

Ve - hnh6 0rl II,7 o = 2Ano

(26)

where A is the rate of deformation tensor defined by A = l[V_v + (V_v)'], and h . , ho, hd. art' the metric coefficients defined by ha =

tg~rr ,

K = 77, 0, ~b

(27)

Stokes flow in spheroidal particle-in-cell models

1471

where _r(=g/al) is the dimensionless position vector. The metric coefficients have the form:

h, = c ~

]

(28)

h o = c

h, = cV-d-zSv"i-

J

and the velocity and the shear stress are defined as follows: 1

0q, ~ _ 1 0~"

(29)

1 01p ve = v¢ = c2~X/--~_ ~21V]--Z_~,2 Or

(30)

v,~ = or = c 2 ~ V - -

IIno = I-I,¢

2

L(za - 1) ~ - - (1 - ~,2) O~"2J + z(1 - z 2) -~r - ~'(1 - ~.2) 0-~"

(31)

The BCs in terms of ~O(r, ~) are given below

BC4:

,2__ ,,2 [ -~

( ¢ - 1)

BCI:

OqJ c2(z z - 1)~" 0"-2=

on

z = ~,,

(32)

BC2:

0q,= _c2( 1 _ ~2)T 0T

on

z = r~

(33)

BC3:

~--~=0

on

z=zt3

(34)

~

.2\ Ol]/

.2, a2~b]

- (1 - g ) ~--~j = ( ¢ - 1) - ~ + ~r(l - g ) - ~

on

z=rt3.

(35')

Because of BC3, BC4 can be written as BC4:

( r 2 - ~r2) ~Ozqj = 2z rg~O Or on

z = rt3.

(35)

The corresponding formulation for Kuwabara type boundary conditions is given in Dassios et al. [13].

3. S O L U T I O N : H A P P E L - T Y P E BCs

The complete solution of equations (17), (18) was obtained in a previous publication [13], through the introduction of the concept of semiseparability. In terms of the standard spheroidal parameterization the solution assumes the semiseparable form: qt('t', ~r) = ~ [g,(,0G,,(~- ) + h,(z)H.(~')].

.=2

(36)

Here the functions g,(r) and hn(z) satisfy specified inhomogeneous O D E ' s of second order (see [13]). G,(ff) and Hn(~') denote Gegenbauer functions of the first and second kind, respectively, of degree ( - 1 / 2 ) and of order n. The expression (36) indicates that the equation of E4~b = 0 does not accept a complete separation of the z and ~ variables in the usual sense. Indeed, although ~0 is written as a sum of separable solutions, the z-dependence of the G,(~) eigenfunction does not coincide with the r-dependence of the H,(~') eigenfunction and, furthermore, the z-functions of order n are linear combinations of Gegenbauer functions of mixed order ( n - - 2 , n, n + 2). In addition, individual terms of the series are not solutions of E$ 33-10-G

G. DASSlOS et al.

1472

equation (17), whereas the complete expansion is [13]. This is exactly what is meant by the term "semiseparation". The Gegenbauer functions have the following characteristics: G.(~) are regular on the x3-axis; H.(~') are singular on the x3-axis; G.(~) are regular in the interior of the spheroid S.; H . ( r ) are regular everywhere (since T > 1), except on the focal segment (where Ir = 1). As we demand of our solution to be regular of the x3-axis and in the space between S= and Sa, we use for the representation of our solution only the terms g.(r)G.(~). Furthermore, taking into account the symmetry of the qJ-field on either side of the equatorial plane (ff = 0), we retain only the even terms of the solution. Thus we obtain

~b(~, ~)= ~

gn(r)G.(~).

(37)

n=2,4...

Now, using the expressions for g.(T) from [13], we obtain 2B

=6 +

+ ~

E

+

n=2,4...

[CnGn(T) + D.Hn(r)]Gn(~)

[~.+2G.(r) + ~n+EHn(T)]Gn+2(~) +

n=2,4-.,

(38)

n=2,4-..

where C2

~t.+2 -

2(2n + 1)

An+2Otn+2]

(39a)

[Bn/3n - B,,+2an+2]

(39b)

[A./3.

-

-

C2

9~.+2 =

2(2n + 1)

with (n - 3)(n - 2) a . = (2n - 3)(2n - 1)

/3. -

(n + 1)(n + 2) (2n - 1)(2n + 1)"

(40)

Here A., B., Cn and D. are constants to be determined from the BCs. Using the orthogonality of the Gegenbauer polynomials and the relation G" = -Pn-l(~')

n = 2, 3 , . . .

and also taking into account that = - G I ( ~ ) = ~ (4n + l)G2n+l(~) n=l

[GI(~) is not orthogonal to the other members of the family, G.(~), n -> 2], the BCs can be written as follows: BCI:

g2('r,,) = 2c2G2(~'~) g2.('c,~) = 0

BC2:

BC3: BC4:(4

- ~.2) ~ n=2,4..,

n = 2, 3 . . . .

(41)

g~('r~) = 2c2GI(T,~) g~.(r~) = 0

n = 2, 3 . . . .

(42)

g2.(ra) -- 0

n = 1, 2, 3

(43)

g~,(ra)G.(~ ) = Era ~ n =2,4...

g.(~:a)G.(~)"

(44)

Stokes flow in spheroidal particle-in-cellmodels Now, setting

c2B2/6= ~1,

1473

the gn(~') functions have the form:

g2(z) = ~tlGl(~') + C2G2(~) + D2H2(r) + ~4G4(r) + ~4H4(T)

(45a)

g,,('r)= ,~,,G,,-a('C)+ ~,,H,,_a('r)+ C,,G,,(v)+ D,,H,,('c)+ ~¢.+zG#+z('t')+ ~n+2nn+2('g), n = 4, 6 . . . . .

(45b)

Using the reccurence relations and the Gengebauer functions (Appendix A) to rewrite the BCs for all even values of n, one obtains a system of linear algebraic equations (Appendix B). Up to an arbitrary even value n = N this system has (2N + 1) unknowns and 2N equations. The existence of an extra unknown creates a problem, when one tries to solve the system by truncating it at n = N. An analogous problem was also encountered in the case of the spheroid-in-cell with Kuwabara type boundary conditions [13]. As in that case, we can provide an additional equation by demanding that the solution for the spheroid-in-cell should tend to the corresponding solution for the sphere-in-cell (that is, the Happel cell solution) as c ~ 0, in such a way that no r 3 sin 0 term appears for infinitesimal values of c. This condition fixes the value of the coefficient of .~/4('t')G2(~) in equation (44) to nil, ~4 = 0.

(46)

Thus, for n = N we now have a system of 2N equations with 2N unknowns, which can be solved readily. The only minor difficulty lies in the somewhat complex form of the coefficients of the unknowns, Appendix B. Once the stream function has been obtained, the velocity components, the vorticity, the flow rate, the drag force, and the permeability can be determined readily, as we will see below. The leading terra of the solution, is given by

~(2)( 't', ~) = [M1Gl('t') + C2G2('t') + M4G4(~) +

D2H2('~)]G2(~)

(47a)

where M1, C2, -'~4 and /)2 are constant coefficients the values of which depend on the geometrical parameters ,',, and ,'a. In any particular instance, these coefficients are obtained readily by solving the linear system

Gl(,ra )

G2('t',, )

G4( "t',~)

H2('t'a)

G~(T,,)

G~(T,~)

G~.('r,~)

H.~('r,,)

G,(,'/,)

G,(,-/,)

-- ~)GI(T/3) 2raG~(ra)

•

- ~) G2('t'/3) - 2TaG~,(ra)

I] C2

D2

- 2r/3G,~(r#)

- 2~aH.~(r#)

i)-] - 1) |

(47b)

where the primes and double-primes denote the first and second derivatives of the corresponding functions, respectively. The leading term seems to contain most of the important physics of the flow, so long as a 3 and 3' are not too large, while the higher order terms provide a slight correction. The same holds true for the first term of the series solution for the case of Kuwabara type BCs, which was developed in [13]. Comparison of the first term of the analytical solution with results obtained numerically show that the first term can be used with excellent to very good accuracy in the domain 3' < - 0 . 3 and -0.20 < a 3 < - 5 . In the remainder of the present work we confine our considerations to values of T and a 3 sufficiently small to use only the leading term ~b(2)(~", ~), equation (47).

G. DASSIOSet aL

1474 4. V E L O C I T Y

COMPONENTS, VORTICITY,

AND FLOWRATE

4.1 Happel-type B C s

Substituting ~0(2)('r, ~') from equation (47) into equations (29) and (30) we obtain the following expressions for the velocity components

v~

G,(~')

(48)

= c2~V'-~-~_ ~ ~X/-~L-i_ 1 g~C'r)

-

lVi-z-~_ ~g~('~" )

c~~

(49)

where g~(z) is the derivative of g2(v)The non-vanishing component of the vorticity, .(z) is obtained readily from the expression o~C2) _

I

[(z a - l)g::(z) - 2 g 2 ( z ) ] V l

_ c ~V-~-zX-: :xq:-~-:~ :~¢,c,~ _

- ~2

2c' ~x/-~-:- :(z' _ :2)

(50)

In order to determine the flowrate through the cell, 0c, (since the outer boundary is considered impermeable to the fluid), we change the system of reference by switching to a moving system of coordinates with its origin fixed on the center of the solid spheroid. Then, the flowrate through the cell can be obtained by integrating the velocity component ~ - a on the ring-shaped equatorial cross section,

qc =

-

f•l

2zr~[avc(z, O) - a] d 6 = - 2 n a S a

fb: ,o[vc(z,

O) - 1] doo

(51)

where ~o (= N / a : ) is the dimensionless radial coordinate in the cylindrical system. For ff = 0 it is w = c V ~ - 1, and equation (51) becomes

qo = - 2 ~ a c

z[vc(z, O) - 1] dz.

2

(52)

T=

Using equation (49) and integrating we obtain

Now, BC1 and BC3 give g2('G) = 1 and g2('r~) = 0. Consequently, we obtain O~2~= = ~ a .

(53)

4.2 Kuwabara-type B C s

Using the expression for qJ(2)(T, ~) developed in Dassios et al. [13], we obtain the following expression for the flowrate through a spheroid-in-cell with Kuwabara-type boundary conditions,

Stokes flow in spheroidal particle-in-cell models

1475

with

v(~)(r~, ~ ) = g~(~) - g2(~) =: -~ {A2[G~(r~) - G~(T.)] ÷ A3[SG4(I.Sfl)Ol(ra) ÷

1(~8)

G4('~a)-6G4(rfl)]÷

A,[H4(ra) - n4(rfl)]}

(55)

where the coefficients A2, A3, A4 are defined as , , _ _,t-/3 5 G4(~a)[H2(r~) - r,,n;(r~)] A 2 = a4(ra)n2(r,~) - G4(r,~)n2(ra)

A3 = Gi(r~)H2(T,,) - G2(r~)Hi(r,,) =

1

2

A4 = G2(r~)G;(%) - G;(r,,)G4(%) + 5 G4(re)[Gz(r,,) _ r~G;(r~)] re

(56)

and D is given by 1 D - - [AEG=(ro) + 6A3G,(ro) + &H2(ro)]. 2G2(r~)

(57)

5. D R A G F O R C E 5.1 H a p p e l - t y p e B C s

The drag force exerted on the spheroid in the x3 direction, Fo, is obtained conveniently from the line integral [3]

/~D= ~la~tL~ASB o93 -O0_n - , - (- £ ~ 7 E O1 2/ d s .

(58)

The integration follows a meridian from A to B, Fig. 1; n is measured normal to the meridian, outwardly, and s is measured along the meridian. In terms of (r, ~) we have to :: cV'-~ - 1V1 - ~2,

0

~-1

-~n - c V ~

0

- ~20z'

ds -

cVTrT ~

d.

~.

(59)

Substituting 0(2)(v, ~) from equation (47) into equation (58) and integrating we obtain

P~) = 6/ral/iaL(2)(,t-,~, r~)

(60)

where the factor L(2)(r,,, za) is given by

L(2)(r'~'r°)-r'~(~'31)3/z{-T2(')[-7~+(--~4

+ l~~/ l n ( Z\ ' t"' , ~+- l1~] /J

+\

~,,-1

~ 2% 11_2

~ - ~o)ln(~° ÷1) +1 } / \%-1/

where

T2(T ) = -2a¢1G1(~" ) ÷ 10~4G4(T)

(61)

1476

G. DASSIOS et al.

and T~(r) = -2~/, G'I(r) + 10~I4G~(~').

(62)

The drag coefficient is defined as CD--

2 1

(63)

2

and therefore C(D2,

24 ~2" = Ree L '(r~, rt3),

(Re_2al__Pt7 / /2 ,"

(64)

As c ~ 0 + (in which case the spheroid-in-cell tends to the sphere-in-cell with the same y value), the drag force given by equation (58) tends to the corresponding expression of the Happel [1] model. 5.2 Kuwabara-type BCs

The drag force on the spheroid, in this case, is obtained by inserting qjt2)(r, ~) from Dassios et al. [13] into equation (58) and integrating. We obtain

/~(D2) = 6~¢a1/2t2L(2)('t',~, "t't3)

(65)

where the factor L(2)(z~, zt~) is given by

5

{ 1

L~2'(v~,, rt3) = 6D ~V,-~Z-i_ 1.

q (~1)2

(5~-6~+

1)

[15"c3,~- T,~ - (15T4 - 6~2 - 1)coth-1 ~,,]}.

(66)

The corresponding drag coefficient is given, again, by equation (64). As c--*0+, the drag force given by equation (65) tends to the correct (see Appendix C) expression for the Kuwabara sphere-in-cell model.

6. PERMEABILITY 6.1 Happel-type BCs

An expression for the permeability of an homogeneous swarm, ~, is obtained as follows. We consider a control volume that has a r e a / ] normal to the x3-axis (which is the direction of the main flow), and length 7. The space-averaged velocity is denoted by g/and the space-averaged pressure drop along 7 is denoted by AP (P denotes the total pressure, i.e. pressure that accounts for the hydrostatic effect, or pressure at zero-gravity). The control volume contains Ns solid spheroids, with Ns given by 3y~l Ns - 4tra2~ a .

(67)

At steady state, a linear momentum balance on the fluid in the control volume gives .4aP - NsFD = 0

(68)

From equations (67) and (68) we obtain A[' [

33, 47r~ 2-----~3 FD.

(69)

Stokes flow in spheroidal particle-in-cell models

1477

Now, the permeability is defined as

~_

Oti (AP/7)"

(70)

Combining equations (58), (59), (68) and (69) and defining a dimensionless permeability k(=/~/ti 2) we obtain 2a3 0 k = 93~L t~

(71)

The effective (:ross sectional area .4~ that corresponds to a single spheroid-in-cell is needed to evaluate the ratio 0/tT. Since the spheroidal shape is not space-filling (a shortcoming that is also shared by the sphere-in-cell models; see discussion in Tien [6]), we proceed as follows. We assume that the control volume is prism-shaped with cross-sectional area ,4c and length ic in the direction of the macroscopic flow. The cross-section can be square or hexagonal, or can have any other convenient plane-filling shape. Furthermore, the volume of the space-filling unit prism must be equal to that of the spheroid-in-cell. At this point, we introduce a constant of proportionality Z, so that we have Ac = (2,~1) 2

~c = 2Zt~3.

(72)

This factor is to be determined from the requirement that --

4

2

3'Acl¢ = ~ ]['/~ 1/~ 3Equations (72) and (73) give

(73)

(717)1/3

Z = \~7/

•

(74)

The flowrate through a unit prism must be equal to that through a spheroid-in-cell, and thus ~c = .4c0 = 4

a20.

(75)

Combining equations (53) and (75) we obtain (76)

t7 Finally, equation (71) becomes

( 717)1/3 a3b2 k (2) =

\-~T/

3L (2)

(77)

where L (2) is given by equation (61).

6.2 Kuwabara-type BCs In a similar way we obtain (~)2/3/i;1/3 y(2 )

(78)

and k(2) = [ 7r ]l/S asy(2) where y(2) and L (2) are given by equations (55) and (66), respectively.

(79)

1478

G. DASSIOS et

al.

7. PROLATE SPHEROID IN AN UNBOUNDED FLUID (REVISITED) The solution to this problem was obtained by Payne and Pell [16] and, also, Happel and Brenner [3]. It is of interest to see how that solution is related to the present work. For the sake of continuity, however, this discussion is given in Appendix C.

8. THE CASE OF THE OBLATE SPHEROID-IN-CELL

8.1 Happel-type BCs The spheroid is oblate when a 3 ~ 1. In this case the semifocal distance is given by = V1 - a 2

(80)

whereas b 3 and b~ are obtained from expressions similar to (11) and (7),

[~2,311/3 12y

Vt~yy/ + t 3 )

J

+ [~y- Vt~y/

b, = V-b3z + ?2.

(81) (82)

As is shown in Payne and Pell [16] (also, Happel and Brenner [3]), the solution qJ(a, ~) for the oblate spheroid can be obtained by making the transformations

c~-ig

v---~iA and

(83)

which is the same as determining ? with use of equation (80) instead of (2). The values of A on the inner spheroid, A~, and the outer spheroid, At3, are given by a3

b3

;t~ = - - ~

a~ =--_.c

(84)

From equations (47), (48) and (49) we obtain O(2)(A, ~) = [M~G,(iA) + C~G2(iA) + D~H2(iA) + M4G4(iA)]G2(~) U(A2) =

v7 ) =

ax(~)

b_2A V ~

Az.#-A-r--~g2(iA)

G2(~')

c2 A2~-~ IVT~_~Eg~(iA)

(85) (86)

(87)

where the transformation (T, c)---~(iA, -iF) has been made in the expressions for A~, C~, D~, A~, B~ as well. The flowrate through the cell, ~2), is given again by equation (53). Finally, the drag f o r c e / ~ ) is given by F~ ) = 6zt2aal L(2)(A=, At3)

(88)

where L(2)(Aa, A B ) -

Aa(Aza + 1) 3/2

3

1

{.-A_~A~) [A_~+ (~_~_ A~)cot_ 1A~]

(89)

Stokes flow in spheroidal particle-in-cell models

1479

with A2(A.) = 2M'IGI(A.) + 10M~G4(A.) A~(A.) = 2M'1G~(A.) + 10M~G~(A.).

(90)

The permeability is given by equation (77) again, but here L (2) is given by equation (90).

8.2 Kuwabara-type BCs The conversion is done in a way entirely similar to the one used above, and it is omitted for the sake of brevity. 9. SAMPLE C A L C U L A T I O N S Sample streamlines and isovorticity lines are shown in Figs 3(a) and (b) for two typical cases: Fig. 3(a) shows the case of a prolate spheroid (a3 = 3) with moderate packing density (y = 0.05), whereas Fig. 3(b) shows the case of an oblate spheroid (a 3 = 1/2) with high packing density (y = 0 . 3 ) A comparison of streamlines obtained using Kuwabara-type BCs with streamlines obtained using Happel-type BCs is given in Fig. 4. A plot of Re C~ ) versus y for several prolate spheroid aspect ratio values (a3 = 1, 2, 3, 4, 5) is given in Fig. 5. As it can be seen, Re C~ ) increases slowly with increasing y, for small y values (say y < -0.005), and more rapidly at high y values. For any given value of y, the minimum value of Re C~ ~ iis obtained at a3 = 1 (sphere-in-cell) and increases as a3 increases, that is, as the deviation from sphericity increases. A similar plot for several oblate spheroid aspect ratio values (a3 = 1, 1/2, 1/3, 1/4, 1/5) is given in Fig. 6. Again, Re C~ ~ increases with increasing y, slowly for small y values, and rapidly for large ones. Here, however, the behavior of the Re C~ ~ versus y curves is more complicated than it is in the case of the prolate spheroid. Specifically, for small values of y, Re C(D2) decrease:s as the deviation from the spherical geometry increases, whereas, for relatively high values of y, Re C(D 2~ increases as the deviation from the spherical geometry increases. These observations can be explained as follows. For small y values the flow around a spheroid is not influenced very much by neighboring bodies. For fixed radius ~1, the surface of a prolate spheroid increases whereas the surface of an oblate spheroid decreases as the deviation from sphericity increases. Under creeping flow conditions this difference in surface area causes the observed difference in the drag coefficient for small y values (see also Happel and Brenner [3, pp. 145-149 and pp. 154-156]. For large 3/ values the effect of neighboring bodies becomes dominant, and an important difference between prolate and oblate spheroidal geometry expresses itself. For fixed y, the equatorial area available to the flow in a prolate spheroid-in-cell increases with increasing deviation from sphericity, whereas it decreases in an oblate spheroid-in-cell. It is this feature of the oblate cell that can cause an increase of Re Co, despite the reduction in the surface of the solid, provided that y is sufficiently large. This behavior of Re C~ ) is also observed in the case of the spheroid-in-cell model with Kuwabara-type boundary conditions. It is interesting to determine the magnitude of the drag force exerted on a spheroid-in-cell in relation to that exerted on a sphere-in-cell under comparable conditions, specifically, for equal volumes of the inner solids, the same solid volume fraction 7, the same fluid viscosity, and the same approach w~Aocityin both cells. To this end, we define the ratio

)2" = -/~D FD,s

(91)

where PD,s is the drag force on the sphere-in-cell with the same type of BCs on the outer boundary. The radius of the inner sphere, tL is given by a = (1 +

c2)1/6al .

(92)

1480

G. DASSIOS et aL

8 '

I

Axis Ratio "7 = 0 . 0 5

6

I

' ~

I

3

2

;<3

0 co¢= 2.0 --2

--

-4

--

2.7= 2.0

1.0 0.4

I

0.1 -6

--

-8 -8

I -6

I -4

0.2

0.01

0.05

I

I

I

I

2

4

6

-2

0

8

X1

(a)

2

I

I

I

I

Axis Ratio -/ =0.3

X3

I =

0.5

I

I

l 1

I

I

0 co~= 2 . 0 ~ 0.02

0.01 0.1

-1

-2 -2

I

I -1

I

I 0

2

X1 (b) Fig. 3. Sample streamlines and isovorticity lines around: (a) a prolate spheroid (a3 = 3), for 7 = 0.05, and (b) an oblate spheroid (a 3 = 1/2), for ~, = 0.3.

Stokes flow in spheroidal particle-in-cell models

I

1481

I

Axi$ R a t i o = 2 3v I 0 . 0 5

]

I

Kuwabara type BCs

Happel type BCs

2

X~

0

ii.

-2

,I, = 2 . 7 .

=

~,

~

-6 -6

~

I

I

-2

0

1.5

=

1.3

=

q/

-4

=

-oi

•

1.9 j

---- 0 . 9 ~ 0.4 0.1 /

-4

",1,

0.2

0.05

I

I

2

4

6

Xl

Fig. 4. Sample streamlines in two identical unit cells with a 3 = 3.0 and 3' = 0.1: (a) Kuwabara-type BCs; (b) Happel-type BCs.

H a p p e l - t y p e BCs, combining equations (60), (91), (92) and equation (10) f r o m [1] we tin 3 (2 - 3",/1/3 + 3T 5/3 - 272) X ~ ) = 2(1 + c2) 1'6 (3 + 275'3) L(2)(z~' ~a) (93) t L (2) given by equation (61). N o t e that 1 + c 2 =

z2/(~

-

1).

1000 B D

c~ (3

100 -

Q:

10 0.001

I

I

0.01

0.1

0.3

"7

Fig. 5. Plot of Re C ~ ) versus the solid volume fraction 3' for several prolate spheroid aspect ratio values.

1482

G. DASSIOS et al. 1000

~

100

•.f

10 0.001

I

i

0.01

0.1

0.3

"7 Fig. 6. Plot of Re C~ ) versus the solid volume fraction 7 for several oblate spheroid aspect ratio

values. For Kuwabara-type BCs, combining equations (65), (91), (92) and equation (3.18') in Appendix C, we obtain X~) -

(5 - 9y ~/3 + 57 - y 2) 5(1 + c2) 1/6 L(2)(z~, z~)

(94)

with L (2) given by equation (66). The corresponding expressions for the oblate spheroid-in-cell are obtained with the transformation (T, c) ~ (iA, -iF). In Fig. 7 calculated values of X (2) are plotted versus y for several prolate and oblate spheroids-in-cell. The reason for the difference in the behavior of the curves for the prolate and oblate geometry was discussed above. We observe that this difference is substantially more pronounced in the case of Kuwabara-type BCs. In Fig. 8 a plot of the dimensionless permability k (2) versus 3' for several aspect ratio values is given, for Happel-type BCs. As it can be seen, k (2) decreases rapidly with increasing y, and the slope increases as y increases. For both prolate and oblate geometries, and for any fixed value of y, k (2) increases monotonically with increasing aspect ratio a3. A comparison of Re C~ ~ and k ~2~ values obtained with Kuwabara-type BCs with the corresponding results obtained with Happel-type BCs is given in Figs 9 and 10, respectively. The drag coefficient calculated with Kuwabara type BCs is systematically higher than that with Happel type BCs (see also Appendix C). The difference is small (albeit discernible) for small values of the packing density, y, and increases for large values of y, because the difference caused by the different BCs on the outer envelope becomes more pronounced as the density of the swarm increases.

10. CONCLUSIONS We present a spheroid-in-cell model with Happel-type boundary conditions, which can be used as the basis for the analysis of heat and mass transport phenomena in swarms of spheroidal particles.

Stokes flowin spheroidal particle-in-cell models

1483

1.2 al P I

1.0 _

............ ..-:.'::"

0.9

0.8

0.7

0.6 0.001

I

I

O.Ol

o.1

0.3

(a) __

a s = 0.2,," ,, 0.25

1.4

i

i

1.2 1

1.0

:

,

:-:----::

/

I

I

"o 33s ~

I

"

"

~

"

I"

I

""

0.5

....

0.8

0.6

0.4 0.001

I 0.01

I O. 1

9'

(b) Fig. 7. Plot of the normalized drag force X (2) v e r s u s the solid volume fraction y for several prolate and oblate spheroid aspect ratio values. (a) Unit cell with Happei-type BCs; (b) unit cell with Kuwabara-type BCs. The equation of motion for creeping flow, E4~b = 0, is not separable in spheroidal coordinates in the usual sense, but the method of semiseparation of variables, developed recently by Dassios e t aL [13], is used to obtain the solution. The solution of the boundary value problem in hand is a series, but, fortunately, the leading term, ~(2)(~, ~), gives satisfactory approximation for engineering applications, provided that the solid volume fraction is not too high (say, 3/< 0.03) and the aspect ratio is in a moderate range (say, - 1 / 5 < a3 < -5). (This was established by comparing the first term of the expansion with the numerical solution.) We obtain expressions for the velocity components v,~ and vo, the vorticity coy, the drag force Pt,, the drag coefficient CD, the normalized drag force X and the

1484

G. DASSIOS et al. 1 0 "~

10

s

10 2

10'

10

10

o

-T

1 0 -=

0.001

0.01

0.1

0.3

9,

Fig. 8. Plot of the dimensionless permeability k (2) versus the solid volume fraction 3' for several aspect ratio values.

dimensionless permeability k. Sample calculations are made and the results are presented in graphs. All numerical calculations are based on the leading term of the solution ~t2)(r, ORe CtD2> increases with increasing 3' for prolate as well as oblate spheroids. For small values of 3' the rate of increase is relatively small, but it increases with increasing 7. For prolate spheroids, Re C(D 2~ increases as the deviation from sphericity increases for any value of 3'- The dependence of Re C~ ) on a 3 is more complex in the case of oblate spheroids. For oblate spheroids-in-cell, Re C~ > decreases as the deviation from sphericity increases, for 3' small, whereas it increases for 3' large. This difference in behavior between prolate and oblate

300

K u w a b a r a type B C s

....

/J/,

100 --

_-

i l

a3=

1

3

/

_~

1

Ira3= -~10

0.001

~

0.01

0.1

0.3

7 Fig. 9. Plot of Re C ~ ) versus solid volume fraction 3' for three typical aspect ratio values. Solid lines: unit cell with Happel-type BCs. Dotted lines: unit cell with Kuwabara-type BCs.

Stokes flow in spheroidal particle-in-cellmodels

~"~~ 101~

0.001

a3 =

1.0

0.01

1485

3.0

0.1

0.3

7 Fig. 10. Plot of the dimensionless permeability k (2) v e r s u s solid volume fraction y for three typical aspect ratio values. Solid lines: unit cell with Happel-type BCs. Dotted lines: unit cell with Kuwabara-type BCs. geometries is due to the interplay of two factors, namely, the surface area on the inner spheroid and the equatorial area available to the flow. This difference in behavior is also observed in the case of the normalized drag force, Fig. 7. The dimensionless permeability k decreases rapidly as y increases for prolate and oblate geometry as well. For all values of y, k increases with aj, for prolate as well as oblate geometries. The predictions of the spheroid-in-cell model with Kuwabara-type BCs, concerning Re CD and k as functions of y and aj, are in reasonable agreement with those based on the spheroid-in-cell model with Happel-type BCs, Figs 9 and 10. The latter version of the model, however, has the important feature that it does not imply exchange of mechanical energy with the environment. This is due to the Happel boundary condition, according to which the shear stress vanishes on the outer boundary, and therefore no work is exchanged with the environment; the rate of mechanical energy provided by the motion of the inner spheroid is exactly equal to the rate of viscous energy dissipation in the fluid envelope. In this sense, it is superior to the model with Kuwabara-type BCs, which implies a small, but discernible, rate of exchange of mechanical energy with the environment. The spheroid-in-cell models can be used for heat and mass transfer calculations in swarms of spheroids, in a way similar to that in which the Happel and Kuwabara sphere-in-cell models have been used for swarms of spheres. Acknowledgments--This work was supported by the EC, DG XII, Program BRITE, Project P-2289, Contract No.

RllB-0290-C(AM) and by the Institute of Chemical Engineering and High Temperature Chemical Processes (ICE/HT-FORTH). REFERENCES [1] J. HAPPEL, A.LCh.E. Jl 4, 197 (1958). [2] S. UCHIDA, Inst. Sci. Technol. Univ. Tokyo (in Japanese) 3, 97 (1949); Abstract, Ind. Engng Chem. 46, 1194 (1954) (translation by T. MOTAI). [3] J. HAPPEL and H. BRENNER, Low Reynolds Number Hydrodynamics. Prentice Hall, Englewood Cliffs, N.J. (1965) [also Martinus Nijholl, Dordrecht (1986)]. [4] H. BRENNER, Eng. Sc.D. thesis, New York University (1957). [5] S. KUWABARA, J. Phys. Soc. Japan 14, 527 (1959). [6] C. TIEN, Granular Filtration of Aerosols and Hydrosols. Butterworths, Boston (1989).

1486

G. DASSIOS et al.

[7] [8] [9] [10] [11] [12] [13] [14] [15] [16]

G. H. N E A L E and W. K. N A D E R , A.LCh.E. Jl 20, 530 (1974). R. P F E F F E R and J. H A P P E L , A.LCh.E. Jl 10, 605 (1964). R. PFEFFER, Ind. Engng Chem. Fundam. 3, 380 (1964). V. C. LEVICH, Physicochemical Hydrodynamics, p. 80. Prentice Hall, Englewood Cliffs, N.J. (1962). G. I. T A R D O S , C. G U T F I N G E R and A. A B U A F , A.LCh.E.JI 22, 1147 (1976). N. EPSTEIN and J. H. M A S L I Y A H , Chem. Engng J. 3, 169 (1972). G. DASSIOS, M. H A D J I N I C O L A O U and A. C. P A Y A T A K E S , Q. Appl. Math. 52, 157 (1994). F. A. C O U T E L I E R I S , V. N. B U R G A N O S and A. C. P A Y A T A K E S , J. Colloid Interface Sci. 161, 43 (1993). P. M O O N and D. E. SPENCER, Field Theory Handbook, pp. 28-29. Springer, New York (1971). L. E. P A Y N E and W. H. PELL, J. Fluid Mech. 7, 529 (1960).

(Received 7 September 1994; accepted 7 January 1995)

APPENDIX

A

Gegenbauer Functions A presentation of the Gegenbauer functions of the first and second kind and of their most useful properties is given in Dassios et al. [13]. For reasons of convenience and completeness we list, here, the expressions for G1, (32, G4,//2 and

H4, GI(x ) = - x

(A1)

G2(x ) = 1 (1 - x z)

(A2)

G4(x ) = 81-(5x 2 - 1)(1 - x 2)

(A3)

x+l

x

H2(x)=l (1-x2)ln ~ - 1 1 x+l H4(x) = ~ (5x 2 - 1)(1 - x2)in ~

+2 x 2 + ~ (15x - 13)

(An)

(A5)

where x = ~ or r.

APPENDIX

B

S y s t e m o f L i n e a r A l g e b r a i c E q u a t i o n s A r i s i n g f r o m the H a p p e l - t y p e B o u n d a r y Conditions For n = 2 [with ~4 = 0 see equation (46)]

M, G,(t~,) + C2G2(t~,) + ~4G4(":,~) + D2H2(t~,) = 2G2(t,~)/(r~ - 1)

(BI)

•dlG~(r,~ ) + C2G~('I'a) + ,.~¢4G,~('t',~)+ D2H~('t'o~) = 2Gl('r,~)/(t ~ - 1)

(B2)

M, GI(ro) + C2G2(t/3) + M4G4(rt0) + D2H2(ro) = 0

(B3)

1

2

1

2 fC4G4(T~I) + O4H~,(vt~) + M6G~,(vt~) + ~6H;(vt3) - M4t) = O. (B4) 35

Stokes flow in spheroidal particle-in-cell models

1487

For n = 4 [with ~4 = 0, see equation (46)]

M, G 2 ( ~ ) + M6G6(r~) + B6H6(~) + C4G4(~) + D4H4(~) = 0

(B5)

M4G~(~) + ~t6G~(r~) + ~sHg(~) + C4G~(~) + D4H~(~) = 0

(B6)

•~462(~'~) + ~666('t',~)--I-~6H6(~'0) q- C4G4(i'#) + D4H4(TIj ) = 0

(B7)

4 33 {"~6G~('t'#)+

~oH,~(z#) + C6G~('t',o) + D6Hg(z/~ ) + MsG~(T.) + ~8H[(z#)}

4

For n = 6, 8, 10 . . . .

MnGn 2('r,~)+ ~,,Hn-2('r,~)+ .~,,+2Gn+2('r,~)+ ~'.,+2H,,+2(r,~) + C.,Gn(.r~,)+ D.Hn(T,~) = 0

(B9)

•~,,G.'_z(v,~)+ ~,,H~_z(T,~) + .~,,+2G"+2(v,~) + ~.+2H'+2(v,~ ) + C,,G'(v,~)+ D.,H~('r~,)= 0

(BI0)

~,,G,,_ 2(~) + ~,,H,,_2('r#)+ ~.+2G,,+z(Tt0 + ~.,+zH,,.2(rp) + C.~G.('r#)+ D,,H,,('r#)= 0

(BII)

C.[('~# - y.)G;',(r#) - 2 . # G ' ( r # ) ] + D.[(r~ - yn)H;(r#) - 2 r a i l ' ( r # ) ] +

2

.

~t.+2[(z# - Y.)Gn+2(*a) - 2reG;+2(r~)] + ~.+2[(r~ - y.)n;;+2(r#) - 2 r a n ' + = ( r # ) ]

- a.+~{A.,+=G.(r a) + ~n+~H.('r#) + C.,+~G.+2(ra) + Dn+zH.+z(r#) + .~.+4Gn+2(r#) = -.~,,[(~

-- 1 , , , ) G " _ = ( ~ # )

- 2~#G'_2(~#)]

- ~.[(~

- 3,.)H'/,_2(r#)

+

~n+4Un+4(,r#)}

- 2r#H'_=(r#)]

+ ~8.+2[sf,,_=G._4(r#) + ~.-2H.-4(r/~) + C,, 2G._2(l'#) + Dn_2H._2(t#) + sf.G.(r#) + ~.H.(r#)]

(B12)

where 2n2- 2n-3 (2n + 1)(2n + 3)

3'~

(n - 1)n ~n-2

an+2

=

(2n - 5)(2n - 3)

(n + 1)n (2n + 1)(2n +3)"

(When the expansion is truncated at a value n = N -> 2, the corresponding expression in curly brackets, {-}, which is generated by higher order terms of the expansion, is set equal to nil to effect closure.) ES 33-10-H

1488

G. DASSIOS et al. APPENDIX

C

Corrected Equations of the Kuwabara (1959) Sphere-in-Cell Model (The Notation and Equation Numbers of the Original are Retained) Kuwabara derived an expression for the drag force X, [ref. his equation (3.16) and (3.18)] by calculating the total rate of viscous dissipation of energy, equation (3.15), and setting this equal to XU. In the framework of his formulation, this is not correct. The mechanical power given by the sphere to the fluid, XU, is not all consumed by viscous dissipation in the fluid layer; a small, but non-negligible, part is given to the external boundary. This is easy to see, since the product [-'r~o(b, O)ve(b, 0)], where ~,o is the (r, 0) component of the deviatoric stress tensor, is everywhere positive, except on the two polar points where it is nil (see also Happel and Brener 3, p. 390]). When this is taken into account properly, the correct expression for the drag force and the drag coefficient are obtained from equations (3.16) and (3.17) by setting 1 1 { , = a'~. (3.18') L(~:I=K'= 1 _ g9~ + ~ 3 _ ~ 61 The same expression is obtained by substituting the stream function from equations (3.8)-(3.9) into equation (56) (of the present work) and integrating. Equation (3.18') is also in complete agreement with the results obtained from the solution for the spheroid-in-cell, when the latter tends to become a sphere-in-cell [equation (63) and Fig. 4 of the present work, as 33---~1]. It is interesting to note that the correct expression, (3.18'), is much simpler than the erroneous one, (3.18). This error has remained unreported for a long time (even though Tardos et al. [11] use the correct expression, without commenting on it). One reason may be that the relative magnitute of the difference in the value of L resulting from the error is small (albeit discernible). The following corrected equations should also be used instead of those appearing in the original publication. C = - 1 U ( I + 21-/-3)/K

(3.9c)

D = ~ U/(13g)

(3.9d)

X = 67rgaUL(5).

(3.16)

APPENDIX

D

Prolate Spheroid in an Unbounded Fluid (Revisited) Let us consider the case of a stationary prolate spheroid in an unbounded fluid that flows in the direction of the axis with a uniform approach velocity ax3- This system is governed by the equations that govern the flow in a spheroid-in-cell, with Kuwabara type BCs, if we let b 1-~ ~, keeping c constant, in which case Sty, becomes a sphere of infinite radius, while S,, remains unchanged. The two boundary conditions on St3 are then to be replaced by the limit ~b(r, 5) =c2( r 2 - 1)G2(5)

as

r-~.

(Ol)

This implies that from the general solution, equation (85) in Dassios et al. [13], we should keep only the terms that behave as c2r2 + o ( ~ ) as t--~ 0% that is ~b(~', 5) = Fa.o(T)G2(5) + c2F2.2(T)G2(5)

(D2)

/~2 2A2 4 - 2 n2(r) F2.o(r) = ~ - cGI (r) + ~ - c G4(Z) + C2c G2('r) +/32 c

(D3)

with

r:,2(~) =

~ n4(r) 175 c3

i,, , 1--~c G,(r)

(D4)

where A2,/~2, C2, 1)2, A4, and/~4 are constants to be determined from the boundary conditions. Rearranging we obtain B2

2A24

-2

+ D2 ~ ] G 2 ( s r )

+

/~4 n4(I') 175

"~4 4 1--~c G4(r)] G2(5).

(05)

At this point we have six unknowns, but only three BCs (two on the spheroid and one at infinity), which means that the usual conditions do not suffice to formulate a well-posed problem. To proceed, we impose a geometrical limit condition, by requiring that the solution should tend to the Stokes solution, namely ~bs(r, O) = (Ar -1 + Br + Cr2)sin 2 0,

(D6)

Stokes flow in spheroidal particle-in-cell models

1489

as c ~ O + , that is, as the spheroid tends to become a sphere. Taking in account the asymptotic behavior of the Gegenbauer functions, in the limit as c--* 0+,

cGl(t)---~-r,

c2G2('/') --->- ~1r 2 ,

1H2(.r)_.1 '

54 c4G,(z)---~--~r

1

21

and renaming the constants of the non-vanishing terms as follows -

A 2 --- C2 c2,

B 2 = ~,

-- B 2 ¢

A1 - - - ~

we obtain ~/,(z, g-) = [AzG2(z) + B2H2(~ ) + A 1GI(~)]G2(g-).

(07)

Substituting into the BCs qs('r,~, g-) = 0 and 0~b(z,~, g-)/0z = 0, and into the limiting condition, equation (D1), we obtain (D8)

A 2 = --2C 2

_ ~2 G2(~) - r~c~(~)

B 2 - __ ~ A1 = c 2

-- ~

(D9)

1

(D10)

H2(t,~ ) - z~H~(t~)"

Equations (D7)-(D10) finally give i/J(z, 0 = lc2(z2 - 1)(1 - g-2)[1

[(z~ + 1 ) / ( ~ - 1)]coth -] z - [ z / ( z 2 - 1)] ] - [ ( ~ + 1 ) / ( ~ - 1)]coth -~ z,~ - [z,~/(r~ - 1)]J

(Dll)

which is the solution proposed by Payne and Pell [15], and by Happel and Brenner [3, pp. 154-156]. The solution for the ease when the spheroid translates with velocity (-a_~3), while the fluid at infinity rests, is obtained by subtracting ~ = 1/2c20"2- 1)(1 - g-2) from (Dll). Here, ~b~ denotes the stream function that corresponds to the uniform velocity field _v= _x3- Substituting ~b from equation ( D l l ) into equation (58) and integrating gives the following expression fi)r the drag force, FD = 6r.~.~aLG,(r,,)

(D12)

with 4 Lc,(z,,) = ~ [ ( ~

+ 1)coth -1 ~,~ - z~l-'

(D13)

which is the same with the expression obtained by Payne and Pell [16], and by Happel and Brenner [3], using the limit method developed by Payne and Pell. It is interesting to note that equation (D5) is not an exact solution of equation (17). It should be emphasized that exact solutions of (17) are obtained either in the form of a complete expansion, equation (36), or as closed-form solutions (generalized eigenfunctions) of which equation (D5) is not one, Dassios et aL [13]. However, elimination of the terms containing G4(z) and H4(z) based on the geometrical limit condition leads "fortuitously" to the exact solution. Here we should discuss the approach used by Happel and Brenner, because they, too, arrived at (D5) by assuming a solution in the form g(z)G2(g-). Direct substitution of this product into (17) led to a fourth order differential equation for g(r) that involved the variable g- in non-separable form. Two linearly independent solutions for g(r), namely Gz(r) and H2(ar) were readib' derived from E2~ = 0. Happel and Brenner managed to obtain a third one by independently setting two complementary parts of the expression for E4~ = 0 equal to nil and observing that there existed a common solution, namely G~(z). These three solutions were sufficient to satisfy the three boundary conditions (two on the boundary, one at infinity) of the single spheroid problem. In view of the special technique used to derive their solution, Happel and Brenner carefully stated that equation ( D l l ) is "a possible form for the stream function". The present analysis helps to clarify this point. The same technique could not be applied in the spheroid-in-cell problem, since the two finite boundaries require four linearly independent solutions and, as we have proved, (see also Dassios et al. [13]) there do not exist four linearly independent separable solutions with a g'-dependence expressed via the same Gegenbauer function.

NOMENCLATURE

Ac sO. A. a3, al

cross sectional area of the control volume cross sectional of the space-filling unit prism that corresponds to a spheroid-in-cell constants defined by equation (39a) constants, equation (39a) dimensionless major semiaxis and minor semiaxis (radius) of the

a3, al

~n

Bn b3, bl

solid spheroid, respectively (a 3 equals the aspect ratio of the solid spheroid) major semiaxis and minor semiaxis (radius) of the solid spheroid, respectively constants defined by equation (39b) constants, equation (39b) dimensionless major semiaxis and

1490

G. DASSIOS

CD, C(D 2)

C. c, t~

CO D

o.

minor semiaxis (radius) of the external spheroid, respectively major semiaxis and minor semiaxis (radius) of the external spheroid, respectively friction factor; value of CD based o n I//(2) constants, equation (38) dimensionless and dimensional semifocal length of a prolate spheroid, respectively critical value of c, equation (10) coefficient, given by equation (57) constants, equation (38) drag force exerted on the solid spheroid (in the direction of the x3-axis); value of ~'D based on

et al.

_v, v m vo

(2) ,v r(2) t , ,vO X, X (2) X ~ ), X ~ ) x l , x2, x 3

£1, £2, -r3

dimensionless velocity, and velocity, respectively "0 and 0 components of the dimensionless velocity p values of v, and v o based on q~(2) normalized drag force, equation (91); value of X based on ~O(2) value of X (2) for Happel-type and Kuwabara-type BCs, respectively dimensionless rectilinear coordinates, Fig. 1 unit vectors of the system of coordinates (xl, x2, x3) G r e e k letters

i//(2)

Gn(x)

g.(r) H.(x) h,,('r)

k=

~/a]

k (2)

L(2)(r., ~) 7

7o Us

Gegenbauer function of the first kind, of degree (-1/2) and of order n function appearing in the complete solution, equation (36) Gegenbauer function of the second kind, of degree (-1/2) and of order n function appearing in the complete solution, equation (38) permeability dimensionless permeability dimensionless permeability, based on I/t(2) function defined by equation (61), or equation (66) length of the control volume (in the macroscopic flow direction) length of the space-filling unit prism that corresponds to a spheroid-in-cell number of solid spheroids contained in the control volume

Mc

t~

P P.(x) Pe P.(X) Qn(x)

Re

s~

s~ S

0 /7

coordinate measured normal to the external spheroid, S o, outwardly unit vector normal on the external spheroid, So total pressure Legendre function of the first kind, of degree 0 and of order n Peclet number Legendre function of the first kind, of degree 0 and of order n Legendre function of the second kind, of degree 0 and of order n fiowrate through a spheroid-incell; value of tT¢ based on ~0(a) Reynolds number surface of the internal (solid) spheroid surface of the external spheroid arc length measured along a meridian of St3 superficial velocity through the swarm approach velocity, a is related to U through equation (76), or equation (78)

a an /3 /3n 7 AP e sr "q 0 A2, A3, A4

A2(Aa), A~(A~)

A /2 s¢

~,~ T2(~a), T~(%) r ~'~, "rt3 to, t~ to,p ~o

• = -qJ ~O(=~/tit~ 2) t~ ~b(2) to t%

parameter defined by equation (5) number given by equation (40) parameter defined by equation (6) number given by equation (40) solid volume fraction pressure drop along the control volume (P = total pressure) eccentricity of the solid spheroid spheroidal coordinate, defined by equation (12) spheroidal coordinate, Fig. 2 unit vector in the B-direction spheroidal coordinate, Fig. 2 coefficients, given by equation (56) combinations of Gegenbauer functions for oblate spheroid and of their derivatives respectively, equation (89) oblate spheroidal coordinate, equation (83) dynamic viscosity eccentricity of the solid spheroid, equation (16) dimensionless and dimensional stress tensor fluid density combinations of Gregenbauer functions and of their derivatives, respectively, equation (62) prolate spheroidal coordinate, defined by equation (12) parameters, defined by equations (15), (16) dimensionless and dimensional radial (cylindrical) coordinate, respectively ~b-component of the dimensionless vorticity vector spheroidal coordinate, Fig. 2 unit vector in the ~b-direction proportionality factor, defined by equation (72) opposite value of ~0, Figs 3 and 4 dimensionless stream function stream function leading term of q~ radial cylindrical coordinate, equation (59) ~-component of the vorticity, equation (50)