Power Law Fluid Flow Through Beds of Spheres at Intermediate Reynolds Numbers

0263–8762/04/$30.00+0.00 # 2004 Institution of Chemical Engineers Trans IChemE, Part A, May 2004 Chemical Engineering Research and Design, 82(A5): 642–652

www.ingentaselect.com=titles=02638762.htm

POWER LAW FLUID FLOW THROUGH BEDS OF SPHERES AT INTERMEDIATE REYNOLDS NUMBERS Pressure in Fixed and Distended Beds S. D. DHOLE1, R. P. CHHABRA1,* and V. ESWARAN2 1 Department of Chemical Engineering, Indian Institute of Technology, Kanpur, India Department of Mechanical Engineering, Indian Institute of Technology, Kanpur, India

2

E

xtensive numerical predictions of friction, pressure and total drag coefficients for the steady incompressible flow of power law liquids through assemblages of spherical particles are presented in this paper. The governing equations (continuity and momentum) have been solved numerically for the unknown pressure and velocity components which, in turn, have been used to infer the values of the derived variables like vorticity, stream function, pressure and frictional components of the overall drag coefficient. The inter-particle hydrodynamic interactions have been mimicked by the two commonly used concentric spherein-sphere cell models. The numerical results presented herein span the following ranges of kinematic and physical conditions: power law index, 0.5, 0.6, 0.8, 1, 1.2, 1.4, 1.6, 1.8, thereby including both shear-thinning and shear-thickening flow behaviour; the particle Reynolds number, 1–500; and bed voidages of 0.4, 0.5 and 0.6, thereby covering packed and distended bed conditions. The accuracy of the numerical solution procedure has been validated by carrying out extensive comparisons with the previously available analytical and scant numerical results, mainly limited to the creeping flow conditions. The utility of such a simple approach is demonstrated by presenting detailed comparisons with the experimental results available in the literature. Keywords: power law; drag; spheres; pressure drop; packed beds.

INTRODUCTION

features, the simplest and possibly also the commonest type of non-Newtonian characteristics encountered in chemical and process engineering applications is the shear-thinning behaviour, although in recent years there has been a growing recognition of the importance of shear-thickening behaviour encountered in the handling and processing of concentrated suspensions and thick paste-like systems (Chhabra and Richardson, 1999). In spite of such an overwhelming pragmatic significance, the flow of even the simplest type, namely, power-law liquids, in packed and fluidized beds has received only very scant attention (Chhabra et al., 2001). Notwithstanding the importance of the detailed kinematics of flow, the key parameter of central interest in all these applications is the frictional pressure drop as a function of the system (voidage, size of packing), kinematic conditions (flow rate=Reynolds number) and non-Newtonian flow behaviour parameters. This work is thus concerned with the steady flow of power-law liquids (encompassing both shear-thinning and shear-thickening behaviours) through beds composed of spherical particles over a wide range of porosity values. It is, however, instructive and useful first to present a concise account of the previous pertinent studies available in the literature, especially those dealing with the flow of powerlaw fluids in such systems.

The steady incompressible flow of non-Newtonian fluids through beds of spherical particles represents an idealization of many industrially important processes. Typical examples include the use of polymer solutions in enhanced oil recovery, the use of sand pack beds for the filtration of polymer melts prior to processing, in the filtration of sewage sludges, and in fixed and fluidized beds employed extensively to carry out a range of biological and polymerization reactions (Mauret and Renaud, 1997; Chhabra et al., 2001). Consequently, considerable research efforts have been expended in understanding the fluid mechanical aspects of the non-Newtonian fluid flow in such multi-particle assemblages of spherical particles. Admittedly most liquids of industrial significance including polymer solutions, melts, suspensions and emulsions display a range of nonNewtonian flow characteristics including shear-thinning, shear-thickening, yield stress, time-dependence and viscoelasticity, etc. Amongst the various types of non-Newtonian *Correspondence to: Dr R.P. Chhabra, Department of Chemical Engineering, Indian Institute of Technology, Kanpur 208 016, India. E-mail: [email protected]

642

POWER LAW FLUID FLOW PREVIOUS WORK From a theoretical standpoint, it is readily conceded that, in addition to the governing equations, a mathematical description of inter-particle hydrodynamic interactions is also needed to define the problem completely. There are essentially three approaches available to model these hydrodynamic interactions in such multi-particle systems, namely, the capillary bundle, the submerged objects and the volume averaging techniques. Since detailed descriptions of these approaches are available in previous review papers (Agarwal et al., 1988; Chhabra et al., 2001), only the salient features are recapitulated here. While the various versions of the so-called capillary bundle model are available in the literature (Kemblowski et al., 1989; Chhabra, 1993a, b; Chhabra et al., 2001), perhaps the one developed by Comiti and Renaud (1989), which has been subsequently extended to the flow of power law fluids (Sabiri and Comiti, 1995, 1997; Ciceron et al., 2002a), is the most generalized one, as it accommodates wall effects and also applies to beds of particles of all shapes and=or of mixed sizes. The evaluation of the two structural parameters of the bed (tortuosity and specific surface area exposed to the fluid) which are determined using experiments with a Newtonian medium are germane to the success of this model. These parameters are assumed to be the intrinsic characteristics of the bed and are thus believed to be independent of the type of the fluid, at least for inelastic power law liquids. This approach has been successfully extended to the flow of power law liquids in the beds of mixed size spheres also (Ciceron et al., 2002b). Some of these and=or similar ideas have also been extended to the flow of viscoplastic liquids (Al-Fariss and Pinder, 1987; Dolejs et al., 1998) and of viscoelastic liquids in beds of mixed size spheres (Tiu et al., 1997). Finally, while this model works well over the full range of the Reynolds number (i.e. viscous and transitional regimes), it is limited to only moderate values of bed porosity, approximately
0.5–0.60 (Chhabra, 1993a; Mauret and Renaud, 1997; Ciceron et al., 2002a). In the second approach to submerged objects, the central problem is similar to that of the estimation of drag force on a sphere in the presence of other particles. Within the framework of the submerged objects model, two distinct approaches are available in the literature. In the first approach, the field equations are solved for a preconceived periodic arrangement of spheres and thus the results obtained relate to the specific values of the lateral and axial spacing between the particles, and the resulting drag force on a particle in the array is related to that on an isolated sphere through a function of the mean voidage of the system, the maximum packing density and=or of appropriate non-dimensional non-Newtonian model parameters. Since the maximum possible packing density is strongly dependent on the type of the arrangement, it is thus not possible to extrapolate these predictions to any other arrangements with different geometrical spacings. This approach has not been extended to the flow of generalized Newtonian fluids, but numerical solutions are needed which tend to be computation-intensive even for the creeping flow of Newtonian liquids. In the second approach, although somewhat less rigorous, the interparticle interactions are approximated by postulating each sphere to be surrounded by a hypothetical concentric

643

envelope of fluid. The size of the hypothetical envelope is chosen such that the void fraction of each cell is equal to the overall mean void fraction of the system. This approach has been shown to yield satisfactory predictions of gross macroscopic parameters like drag on swarms of bubbles and droplets moving in power-law and other generalized Newtonian fluids (Jarzebski and Malinowski, 1987; Gummalam and Chhabra, 1987; Tripathi and Chhabra, 1994; Chhabra, 1998; Zhu, 2001). Similarly, scant analytical and numerical predictions are also available for the flow of power law liquids (both n > 1 and n < 1) through beds of spherical particles (Mohan and Raghuraman, 1976a, b; Kawase and Ulbrecht, 1981a, b; Chhabra and Raman, 1984; Satish and Zhu, 1992; Jaiswal et al., 1991, 1992, 1993, 1994). However, most of these studies are limited to the creeping flow region (Re < 1), except for the limited results of Jaiswal et al. (1992), which extend up to Re ¼ 20 and for shear-thinning flow behaviour, i.e. n < 1. Detailed comparisons between these predictions and the corresponding experimental results both for n < 1 and n > 1 show a fair to moderate degree of agreement (Jaiswal et al., 1993, 1994; Chhabra, 1993a, b). It is also appropriate to mention here that the predictions based on this approach are in excellent agreement with the experimental values of friction factor for the flow of shear-thinning liquids past a bundle of circular cylinders both for Newtonian and power law liquids (Tripathi and Chhabra, 1992, 1996; Prasad and Chhabra, 2001; Shibu et al., 2001; Malleswara Rao and Chhabra, 2003) up to about Re  100–200. Finally, a hybrid approach, based on a combination of the capillary bundle and the submerged objects model, has also been developed by Machac and Dolejs (1981) and Dolejs and Machac (1987). This approach also hinges on the determination of two parameters, namely, a characteristic linear dimension and a bed factor, both of which must be evaluated from experimental data with Newtonian fluids. Unlike in the approach of Comiti and Renaud (1989), in this case the so-called bed factor is a function of the shape and orientation of particles in the creeping flow regime and it shows additional dependence on the Reynolds number beyond the creeping flow conditions. Thus, these parameters are not intrinsic properties of the bed. Nor is it clear whether these show additional dependence on nonNewtonian characteristics of the liquid. In the volume averaging method, the central problem is to identify the effective characteristics of the fluid-porous medium such that the resistance to flow is equal to that in the real medium. While this approach has had a degree of success for the flow of Newtonian fluids in a porous medium, it has not proved to be very successful for the flow of power law fluids and also is not completely devoid of empiricism (Chhabra et al., 2001). From the foregoing description, it is thus safe to conclude that very little theoretical information is known about the drag force on an assemblage of spheres in power law fluids at intermediate Reynolds numbers, especially for shear-thickening fluids. This work sets out to fill this gap in the existing literature. In particular, the field equations have been solved numerically for the wide ranges of physical (0.4
e
0.6), rheological (0.5
n
1.8) and kinematic parameters (1
Re
500). The paper is concluded by presenting detailed comparisons for the frictional pressure drop in fixed and distended beds with the other numerical and experimental results available in the literature.

Trans IChemE, Part A, Chemical Engineering Research and Design, 2004, 82(A5): 642–652

644

DHOLE et al.

PROBLEM STATEMENT AND THE CELL MODEL IDEALIZATION Consider the steady and incompressible flow of an inelastic power law fluid past an ensemble of spheres (of radius R) as shown schematically in Figure 1(a). Owing to the f-symmetry, Vf ¼ 0 and no flow variable depends upon the f-coordinate. In spherical coordinates, the equation of continuity, r- and y-components of the momentum equations in dimensionless form are written as (Bird et al., 1987): Continuity 1 @ 2 1 @ (r Vr ) þ (V sin y) ¼ 0 2 r @r r sin y @y y r-component
 @Vr @V V @V V2 þ Vr r þ y r  y @t @r r @y r  n @p 2 1 @ 2 1 (r trr ) þ ¼ þ @r Re r2 @r r sin y  tyy þ tff @  (try sin y)  @y r

(1)

Furthermore, the diffusion terms appearing in the r- and y-components, equations (2) and (3), are expressed in a conservative form and, substituting for the extra stress components using equations (4) and (5), one gets: r-component V2 @Vr 1 @ 2 2 1 @ (Vr Vy sin y)  y þ 2 (r Vr ) þ r @r @t r  r sin y @y  @p 2n Z 2 2Vr 2 @Vy 2Vy cot y D Vr  2  2  ¼ þ @r Re r @y r2 r   nþ1 2 @Z e @Z þ err þ ry (6) @r Re r @y y-component

(2)

@Vy 1 @ 2 1 @ 2 VV (Vy sin y) þ r y þ 2 (r Vr Vy ) þ r @r @t  r sin y @y r 1 @p 2n Z 2 2 @Vr Vy þ D Vy þ 2 ¼  r @y Re r @y r2 sin2 y   2nþ1 @Z e @Z þ ery þ yy @r Re r @y where

y-component @Vy @V V @V VV þ Vr y þ y y þ r y @t @r r @y r  1 @p 2n 1 @ 2 1 þ (r try ) þ ¼ r @y Re r2 @r r sin y  try  tff cot y @  (tyy sin y) þ @y r

D2 

(3)

The power-law model in a dimensionless form is written as (Bird et al., 1987): tij ¼ 2Zeij

(4)

and the viscosity, in turn, is given by:
ðn1Þ=2 I Z¼ 2 2

(5)

where I2 is the second invariant of the rate of deformation tensor and its expression in terms of Vr, Vy and their derivatives are available in standard text books (Bird et al., 1987). In equations (1)–(5), the velocity terms have been scaled using the free stream velocity Vo, the radial coordinate r using the sphere radius R, pressure using, rVo2, the extra stress components using m(Vo=R)n, and finally the viscosity with m(Vo=R)n1.

Figure 1. Schematics of flow and the cell model idealization.



 1 @ 2 @ 1 @ @ r sin y þ r2 @r @r r2 sin y @y @y

(7)

(8)

Note that equations (6) and (7) are identical to those presented by Adachi et al. (1973) for the flow of a power law liquid past a sphere. The Reynolds number, Re, is based on the sphere diameter d ( ¼ 2R) and is defined as follows: rVo2n d n (9) m Furthermore, err, eyy, and ery are the components of the rate of deformation tensor and their expressions in terms of Vr and Vy are available, for instance, in Bird et al. (1987). The cell model approach postulates the each sphere of radius R to be surrounded by a hypothetical envelope of the fluid of radius R 1 , as shown in Figure 1(b). As detailed elsewhere (Chhabra, 1993a, b; Jaiswal et al., 1991; Shibu et al., 2001), the value of R 1 is chosen such that the void fraction of each cell is equal to the overall porosity, e, of the assemblage. Thus, in a dimensionless form: Re ¼

R1 ¼ r1 ¼ (1  e)1=3 (10) R Thus, by simply varying the value of r 1 , one can simulate the beds of various voidages including the limiting conditions of r 1 ! 1 , i.e. e ¼ 1 corresponding to the case of the flow past a single sphere. Conversely, one can readily calculate the value of r 1 for known values of bed porosity. The physically realistic boundary conditions for this flow are that of no-slip at r ¼ 1, i.e. on the surface of the solid sphere. However, much confusion exists regarding the identification of appropriate boundary conditions at the cell boundary, i.e. r ¼ r 1 (Happel and Brenner, 1965). In the free surface cell model, the cell boundary is postulated to be frictionless, i.e. try ¼ 0 thereby emphasizing the noninteracting nature of the cells whereas Kuwabara (1959) advocated the use of the zero-vorticity condition at r ¼ r 1 . Some other plausible boundary conditions together with their

Trans IChemE, Part A, Chemical Engineering Research and Design, 2004, 82(A5): 642–652

POWER LAW FLUID FLOW relative merits and de-merits have also been proposed in the literature (Slobodov and Chepura, 1982). Admittedly, all such boundary conditions represent approximations and it is not possible to offer a sound theoretical justification for any of these boundary conditions at the cell boundary. However, based on the previous extensive comparisons between predictions and experimental results (Chhabra, 1993a, b), the free surface cell model will be used here, although limited results were also obtained for the zero-vorticity cell model to validate the accuracy of the numerical solution procedure and to ascertain the effect of the boundary condition at the cell boundary on the values of drag coefficient. Thus, the boundary conditions used in this work can be expressed as: at r ¼ 1 Vr ¼ 0;

Vy ¼ 0

(11a)

at r ¼ r 1 Vr ¼  cos y

try ¼ 0

(11b)

It is worthwhile re-iterating here that the two cell models, namely the free surface and the zero vorticity models, are identical in all respects except that try ¼ 0 at r ¼ r 1 is prescribed in the former and ory ¼ 0 condition at r ¼ r 1 is implemented in the latter case. In addition to these, the flow is assumed to be symmetrical about y ¼ 0 and y ¼ p planes, i.e. the conditions of Vy ¼ 0 and @ Vr= @ y ¼ 0 are satisfied at these planes. Thus, equations (1), (6) and (7) subject to the boundary conditions outlined in equation (11) have been solved numerically for mapping the flow domain 1
r
r 1 and 0
y
p in terms of Vr (r, y), Vy(r, y) and p(r, y) for a range of values of e, n and Re. A detailed description of the numerical solution procedure is presented in a later section; it is sufficient to add here that, once the velocity and pressure fields are known, the individual and the total drag coefficients can be evaluated as described here. It is customary to introduce a dimensionless drag coefficient as: CD ¼

2FD rVo2 pR2

The pressure component CDP is evaluated as: ðp CDP ¼ 2 pjr¼1 sin 2ydy

(12)

(13a)

645

scheme while the viscous=diffusive terms were approximated using the standard centre difference method. The timestepping solution procedure involves two steps: a predictor step which involves the calculation of the unknown velocities using an assumed pressure field. This is iterative in nature as the viscous terms are implicitly discretized. The pressure and velocity values are then repeatedly corrected until the velocity field satisfies the equation of continuity within prescribed limits (
106). The implicit discretization helps avoid the numerical instability commonly encountered at the relatively small Reynolds numbers used in this work. The time-stepping is pseudo-transient; starting from arbitrary initial conditions, the time-stepping of the momentum equations is continued until steady-state solution is obtained (
103). Once the steady-state and the fully converged velocity profile is known, the true pressure field is obtained by solving the Poisson equation, together with appropriate boundary conditions. These in turn were used to evaluate the derived variables including vorticity, stream function, friction, pressure drag coefficients, etc. All results reported herein have been checked for grid independence and most of the computations are based on the use of a 52  52 mesh. Further refinements of the mesh altered the results by less than 0.5%, but at the expense of a significant increase in CPU time. Furthermore, the drag values stabilized at least up to four significant digits were finally accepted encompassing wide ranges of the power-law index, Reynolds number and the bed porosity. RESULTS AND DISCUSSION A simple scaling of the governing equations and the boundary conditions suggests the individual and total drag coefficients to be functions of the power-law index, Reynolds number and the bed voidage. In this work, theoretical estimates of the individual and total drag coefficients have been obtained as functions of the power-law index (0.5
n
1.8), the Reynolds number (1
Re
500) and for a range of values of e (0.4, 0.5 and 0.6) for both cell models. Thus, these results can be used to delineate the influence of n as well as that of the cell boundary conditions on the overall hydrodynamic parameters. First, it is useful to establish the accuracy of the numerical solution procedure used in this work.

0

and the frictional component CDF is given by (Adachi et al., 1973): ð 2nþ2 p (Zo)jr¼1 sin2 ydy (13b) CDF ¼ Re 0 Thus, the total drag coefficient CD is simply the sum of CDP and CDF . NUMERICAL SOLUTION PROCEDURE The numerical solution methodology used here is a SMAC-type implicit scheme originally due to Harlow and Welch (1965), and it has been used in our previous work (Shibu et al., 2001; Mandhani et al., 2002; Shukla et al., 2004). This scheme was implemented on a staggered grid for the solution of the equations of continuity and momentum, i.e. equations (1), (6) and (7) in dimensionless form. The convective terms have been discretized using the upwind

Validation of Results In the limit of Re ! 0 (creeping flow region) and n ¼ 1, i.e. for Newtonian fluid behaviour, an analytical solution is available which can be used to estimate the values of CDP , CDF and CD (Happel and Brenner, 1965; Jaiswal et al., 1991). Table 1 presents a comparison of the analytical and numerical values of the individual and total drag coefficients where the two values of the individual and the total drag coefficients are seen to differ at most by 0.6%. This also suggests that, at Re ¼ 0.1, the creeping flow model is a good approximation. Furthermore, detailed comparisons between the point velocities based on the analytical solution and those obtained in the present work showed the two values to be within 0.75–1% of each other in the whole flow domain except for a small region near the front and rear stagnation points. Similarly, a comparison between the present and the literature values (Jaiswal et al., 1991) for n ¼ 1 and for a range of values of

Trans IChemE, Part A, Chemical Engineering Research and Design, 2004, 82(A5): 642–652

646

DHOLE et al. Table 1. Comparison between analytical and numerical results for n ¼ 1 and Re ¼ 0.1 (for free surface cell model). Analytical e 0.4 0.5 0.6 0.8

Numerical

CDP

CDF

CD

CDP

CDF

CD

14,400.00 5681.08 2476.38 550.93

6088.00 3448.10 2073.10 805.28

20,488.00 9131.18 4549.48 1356.21

14,356.5 5696.61 2479.68 552.18

6098.3 3492.5 2098.37 808.13

20,454.8 9189.11 4578.05 1360.31

e over the range 1
Re
100 revealed the two values to be in excellent agreement (within 0.8%). With reference to the zero vorticity cell model, Table 2 presents a typical comparison between the literature values and the present results. Again, the two values are seen to be generally within 5–6% of each other. This order of discrepancy is not at all uncommon in such numerical studies. Finally, the present results were contrasted with those of Jaiswal et al. (1992) for the flow of power law fluids in the range 1
Re
20 and 0.6 < n < 1.0 and the two values of drag coefficient were found to be generally within 4–5% of each other in most cases. From these representative comparisons with previous results coupled with our past experience (Shibu et al., 2001), it is believed that the present numerical predictions of drag coefficients are reliable and accurate to within 2–3%.

Drag Phenomena Figure 2 shows the dependence of the friction, pressure and total drag coefficients on the Reynolds number and porosity for the two boundary conditions, namely, the zero shear stress and zero vorticity condition for n ¼ 0.5, 0.8, 1.6 and 1.8, thereby encompassing highly shear-thinning (n ¼ 0.5) to extremely shear-thickening fluids (n ¼ 1.8). Based on a detailed analysis of these results, the key findings can be summarised as follows. (i) For a fixed value of porosity (e) and the Reynolds number, the friction, pressure and total drag coefficients increase as the value of the flow behaviour index is progressively increased from n ¼ 0.5 to 1.8. However, the rate of increase is strongly dependent on the values of e, Re and n. Thus, for instance, for e ¼ 0.4 and Re ¼ 1, the values of CDP increase from 117.06 for n ¼ 0.5 to 8717 for n ¼ 1.4 while the corresponding value for n ¼ 1 is 1438. The corresponding figures for e ¼ 0.4 and Re ¼ 100 are 2.37 for n ¼ 0.5, 21.34 for n ¼ 1 and 89 for n ¼ 1.4. The friction drag coefficient CDF shows qualitatively similar behaviour, although the relative change is not as steep. Thus, for instance, for the aforementioned sets of conditions (i.e. for Table 2. Typical comparison between the present and literature (LeClair and Hamielec, 1968) values of CD for n ¼ 1 for the zero vorticity cell model. e ¼ 0.5

e ¼ 0.6

Re

Present

Literature

Present

Literature

1 10 20 50 100 500

1085.17 108.30 55.87 24.74 14.46 6.79

1130.0 114.5 58.00 25.50 14.80 6.68

549.37 56.13 27.80 12.92 8.53 3.96

571.00 53.56 27.57 13.80 8.85 4.16

e ¼ 0.4, Re ¼ 1), the value of CDF ranges from 73 for n ¼ 0.5, 610 for n ¼ 1 to 3269 for n ¼ 1.4. Furthermore, the use of the zero vorticity condition at the cell boundary always yielded slightly larger values of the drag coefficient than that for the free surface cell model. The difference between these two values is also strongly dependent on the values of e, Re and n. It varies from 20% at low Reynolds numbers to about 10% at Re  100–500. Clearly, the higher value of drag is attributable to the additional dissipation of energy due to the non-zero value of shear stress at the cell boundary in this case. (ii) Irrespective of the values of n and Re, the values of drag coefficient decrease with the increasing value of the bed porosity. This is simply due to the reduction in the velocity gradients occurring in the bed. (iii) Some further insights can be gained by examining the relative contributions of the pressure and friction drag to the total drag coefficient. Representative results displaying the variation of CDP=CDF with n, e and Re are shown in Figure 3. Evidently, the ratio CDP=CDF is strongly influenced by all the three parameters, namely, e, n and Re. Firstly, both cell models lead to qualitatively similar trends. In low porosity systems, the pressure drag dominates over the friction drag and while this ratio changes very little up to about Re ¼ 20–50 for shear-thinning materials, it rises quite steeply, attaining maximum values which are larger than the corresponding values for the Newtonian fluid behaviour at high Reynolds numbers. This clearly shows that, as the value of the Reynolds number is progressively increased, the pressure forces increase much more rapidly than the viscous forces, for constant values of n and e. This behaviour is observed for all values of bed porosity. On the other hand, as the value of bed porosity increases, the ratio (CDP=CDF) decreases for constant values of n and Re. Likewise, for shear-thickening fluids, while the ratio (CDP=CDF) is always higher than the corresponding value for the Newtonian fluids, but the dependence on the Reynolds number weakens with the increasing value of the flow behaviour index, e.g. see the figure where for n ¼ 1.8, the value of the drag ratio changes only very little over the Reynolds number range 1
Re
500. This is presumably so due to the steeply rising viscous forces in shear-thickening systems. This suggests that the critical Reynolds number marking the cessation of the creeping flow conditions also tends to be higher under these conditions. In general, with lower bed porosity and=or higher the value of n, the creeping flow persists up to higher values of the Reynolds number. This is further supported by the relative constancy of the product (CD  Re) and by the detailed examination of the corresponding streamline plots. However, under very highly shear-thinning conditions (n
0.5) and especially at high Reynolds number, (Re 0 300) the pressure drag was seen to be negative and this was also consistent with the corresponding flow (vorti-

Trans IChemE, Part A, Chemical Engineering Research and Design, 2004, 82(A5): 642–652

POWER LAW FLUID FLOW

647

Figure 2. Variation of total (CD), pressure (CDP) and friction (CDF) drag coefficient with the particle Reynolds number (Re), power-law index (n) and voidage (e) for the two cell models (FS=ZV): (
) e ¼ 0.4, model FS; (s) e ¼ 0.4, model ZV; (m) e ¼ 0.5, model FS; (n) e ¼ 0.5, model ZV; (j) e ¼ 0.6, model FS; (u) e ¼ 0.6, model ZV. (The same legend has been used in all figures.)

city) patterns wherein the ‘reverse flow’ was seen to occur in the rear of the particle. This behaviour was seen only in one case of Re ¼ 500 and n ¼ 0.5, and it needs to be emphasized here that all possible remedial measures like mesh refinement, varying the relaxation factor, etc. were taken but this phenomenon was always present. Therefore, this is believed to be not a numerical artefact. While the exact reasons for

this behaviour are not immediately obvious, it is certainly a manifestation of the complex interplay between the inertial, pressure and viscous forces in such highly shear-thinning systems, and these interactions are believed to also be responsible for the ratio (CDP=CDF) going through a maximum value at about Re  50 in highly shear-thinning systems (n ¼ 0.5). Some further insights can also be

Trans IChemE, Part A, Chemical Engineering Research and Design, 2004, 82(A5): 642–652

648

DHOLE et al.

Figure 2. Continued

gained by an order of magnitude analysis (Jaiswal et al., 1991, 1993). Within the framework of the free surface cell model, and for moderately dense assemblages, an order of magnitude analysis suggests the following dependences: mean angular velocity  0(1=d) inertial forces  0(1=d2 ) shear forces  0(1=d2n ) pressure forces  0(1=dmo )

(14)

Figure 3. Functional dependence of the ratio (CDP=CDF) on the Reynolds number, power-law index and voidage.

where m0 ¼ 2 for inertia dominated flows (low n and=or high Re) and m0 ¼ 2n þ 1 (n < 1) for viscous flow, i.e. high n and low Re; d is the dimensionless radial gap, i.e. d ¼ (r 1 7 1). Clearly, as the porosity increases, d will increase and the inertial forces will increase as (1=d2), whereas the pressure forces will scale less steeply. This is

Trans IChemE, Part A, Chemical Engineering Research and Design, 2004, 82(A5): 642–652

POWER LAW FLUID FLOW

Figure 3. Continued

qualitatively in line with the results seen in Figure 3. Also, the shear forces will scale with e less steeply for shearthinning fluids but more steeply for shear-thickening fluids, though the above-noted expression for shear forces is strictly valid only at low Reynolds numbers. COMPARISON WITH EXPERIMENTS Many investigators have reported pressure drop results for the flow of inelastic power law liquids in packed and fluidized beds of spherical particles, and most of these have been critically reviewed by Kemblowski et al. (1989) and more recently by Chhabra et al. (2001). Based on their extensive studies, Kemblowski et al. (1989) recommend the following equation for estimating the pressure drop in a packed bed of spheres (e  0.37–0.4; 0.5
n
1.6 and Rek < 115): f ¼

150 1:75 mc K12 þ 2 Rek {m2c (K1  1)2 þ K12 }1=2

(15)

649

Thus, equations (17) and (18) allow a direct comparison between the present results and the predictions of equation (15), as shown in Figure 4 for e ¼ 0.4. Clearly, the correspondence between the two results is seen to be satisfactory over the range of conditions as: 0.5
n
1.6. Note that, since the two definitions of the Reynolds number, Re and Rek, are related through a function of n and e, this comparison is limited to the maximum value of Re ¼  10–15, which corresponds to the maximum value of Rek  110. Generally, the present predictions begin to deviate from that of equation (15) as the values of the power-law index decrease and=or the Reynolds number increases or both. Furthermore, while the average error associated with equation (15) was stated to be 15%, a detailed examination of their results reveals maximum errors of the order of 100% in some cases. Furthermore, some of the test fluids used in this study are now known to be viscoelastic, which is obviously not accounted for in any way in equation (15). In view of this, the correspondence seen in Figure 4 is about as good as can be expected in this kind of work. For the results shown in Figure 4, the mean error is 36% for the free surface cell model and 24% for the zero vorticity cell model, both of which are well within the accuracy of the experimental results. Similarly, the modified capillary model of Comiti and Renaud (1989) for power law fluids is given by the following equation: fpore ¼

16 þ 0:194 Repore

where fpore and Repore are defined as follows:


fpore

 Dp de3 ¼ L 3rVo2 t3 (1  e)

(20)

and Repore ¼ Re  F(n,e)

where

(19)

(21a) 2(n1) 2n

K1 ¼ xRek Rek ¼

rVo2n d

c(1  e)
n m 3 9þ c¼ (150Se)(1n)=2 12 n d2 e3 S¼ 150 (1  e)2

(16a)

F(n,e) ¼

e t 22n3 {3(3n þ 1)=4n}n (1  e)n

(21b)

(16b) (16c) (16d)

and finally, mc and x were correlated with the flow behaviour index only involving polynomials of 5–6 degrees in power law index, n. The drag coefficient can be linked to the friction factor by equating the rate of energy dissipation per unit volume and this approach yields 3 f ¼ CD e3 4 where the friction factor f is defined as: f ¼

Dpde3 rLVo2 (1  e)

(17)

(18)

Figure 4. Typical comparison between the present results and the predictions of equation (15) for e ¼ 0.4 (open symbols for zero vorticity and solid symbols for free surface cell model).

Trans IChemE, Part A, Chemical Engineering Research and Design, 2004, 82(A5): 642–652

650

DHOLE et al.

Finally, for a bed of spheres, the tortuosity t is related to the bed porosity as: t ¼ 1  0:41 ln e

(22)

Equation (19) has been found to correlate the pressure drop data in beds of mono-size and mixed-size spheres (Chhabra et al., 2001; Ciceron et al., 2002a, b) over wide ranges of porosity conditions (0.31
e 9 0.6) and up to about Repore  600. Figure 5(a–c) shows comparisons between the present predictions (for both cell models) and that of equation (19) for three values of e, i.e. 0.4, 0.5 and 0.6, respectively, in the range of applicability of equation (19) in terms of fpore. Again, good correspondence exists between the experiments and predictions at low Reynolds numbers Re
50–100 depending upon the value of the power law index and porosity. For instance for e ¼ 0.4, the predictions and experiments begin to diverge at Re ¼ 25 for n ¼ 0.5, whereas for n ¼ 0.8, good match exists up to about Re  60. However, as the value of the porosity increases, good match is obtained under most conditions, as can be seen in Figure 5(b) and (c) for e ¼ 0.5 and e ¼ 0.6, respectively. For the results shown in Figure 5(a–c), the mean deviation is 38% for the free surface cell model and 39% for the zero vorticity cell model. In assessing the comparisons shown in Figures 4 and 5, the following factors must be borne in mind. Equation (15) is based on the flow of polymer melts and solutions which often show viscoelastic effects which are captured neither by equation (15) nor by equation (19), nor by the present theory. Indeed, the pressure drop results for viscoelastic fluids can vary by up to a factor of two or three depending upon the extensional viscosity behaviour of these fluids and due to other phenomena like slip and adsorption (Chhabra, 1993b). Furthermore, equation (15) has more than a dozen fitted parameters, especially the polynomial expressions for mc and x. This is probably difficult to justify as their data relates to only five or six values of the flow behaviour index and to rather narrow range of values of e. In contrast, equation (19) relies only on two structural parameters and uses an appropriate value of the flow behaviour index, depending upon the prevailing shear rate in the bed, whereas the present model predictions are based on the constant values of n. Similarly, while the mean errors of associated with the use of equation (19) is quite reasonable (10–12%), maximum errors of up to 30–40% are not uncommon in such experimental studies, which are limited to only shear-thinning (n < 1) fluids. Besides, except for the degree of arbitrariness inherent in the predictions due to the boundary condition at the cell boundary, no further adjustment is possible in the theoretical results. Bearing in mind these factors coupled with the highly idealized nature of the cell model, the match between experiments and the present predictions (especially for e ¼ 0.5 and e ¼ 0.6) is regarded as satisfactory and acceptable for shear-thinning and shearthickening fluids. These results also confirm the expectation that the cell model is perhaps more appropriate in the porosity range, e 0 0.5 or so and at Reynolds number not exceeding  100 or so (Chhabra et al., 2001; Ciceron et al., 2002a). Finally, it is appropriate to re-iterate here that equation (19) is based on experimental data for liquids with n < 1 and it is really not expected to apply to shearthickening fluids with n > 1. Indeed this explains in part the rather large discrepancies seen in Figure 5.

Figure 5. Typical comparisons between the present results and the predictions of Equation (19). (a) e ¼ 0.4; (b) e ¼ 0.5; (c) e ¼ 0.6. (Open symbols for zero vorticity and solid symbols for free surface cell model.)

Before leaving this section, it is desirable to develop a correlation based on the present numerical results. In order to retain the basic form of the well known Ergun equation, the following equation was obtained by using the regression procedure: f ¼

230:4 þ 0:33 Rek

(23)

Trans IChemE, Part A, Chemical Engineering Research and Design, 2004, 82(A5): 642–652

POWER LAW FLUID FLOW Equation (23) correlates the complete set of present results for n > 1 and n < 1 with an average error of about 31% for 250 numerical data points while 17 points deviate by more than 60%; the maximum deviation is 94% for the extreme case of n ¼ 1.8 and Re ¼ 500. There are no discernable trends with respect to n or e or Re. For shear thickening fluids, high values of the power law index are accompanied by extremely high values of the consistency coefficient, m, thereby resulting in small values of the Reynolds number. Similarly, the values of n 9 0.5 or so are usually associated with significant levels of viscoelasticity and thus are beyond the range of applicability of the present study. In summary, equation (23), although completely empirical, can be used with an average error of 31% under most conditions of practical interest. CONCLUSIONS In this study, the steady incompressible flow of power-law fluids through beds of spherical particles at moderate Reynolds numbers has been investigated numerically. The particle–particle hydrodynamic interactions have been approximated using the two commonly employed cell models, namely the free surface cell and the zero vorticity cell models. The field equations within the framework of the concentric spheres cell model have been solved for the unknown velocity and pressure fields using a finite difference method. Extensive numerical predictions of drag on spheres spanning wide ranges of physical and kinematic conditions as 1
Re
500, 1.8 n 0.5 and 0.4
e
0.6 have been obtained. In addition, the variation of the individual contributions of the pressure and friction to the total drag with the Reynolds number, power law index and the bed voidage has been studied. Irrespective of the values of Re, e and n, the resistance to flow in shear-thinning liquids is reduced below the corresponding value for the flow of a Newtonian liquid under otherwise identical conditions. Likewise, shearthickening behaviour increases the resistance to flow with reference to an equivalent Newtonian liquid. It is appropriate to highlight here that no prior results are available for shear thickening fluids. The accuracy and reliability of the numerical solution procedure and that of the new results have been ascertained by making detailed comparisons with the previous analytical and=or numerical results for similar systems. The utility of this simple approach has been established by showing fair to good match between the present predictions and the experimental results of pressure drop in beds of spheres available in the literature (0.4 < e < 0.6). This approach thus complements the traditional capillary bundle approach which has been found more suitable in dense systems (low values of e) and possibly at low Reynolds numbers. NOMENCLATURE CD CDF CDP d e FD f fpore K1

drag coefficient friction drag coefficient pressure drag coefficient sphere diameter, m bed porosity drag force, N friction factor pore friction factor, equation (20) modified Reynolds number, equation (16a)

I2 m m0 n p Dp=L R Re Rek Repore R1 r r1 S t Vo Vr, Vy

651

dimensionless second invariant of the rate of deformation tensor power-law consistency coefficient, Pa sn constant, equation (14) power-law index dimensionless pressure pressure gradient, Pa m1 sphere radius, m Reynolds number modified Reynolds number, equation (16b) pore Reynolds number, equation (21) cell radius, m dimensionless radial coordinate dimensionless cell radius constant, equation (16d), m2 time variable, s superficial velocity, m s1 radial and angular velocities, m s1

Greek symbols d (r 1 7 1) dimensionless components of the rate of deformation tensor eij F spherical coordinate Z apparent viscosity, Pa s y spherical coordinate function, equation (15) mC x function, equation (16a) r fluid density, kg m3 o vorticity, s1 t tortuosity factor components of extra stress tensor, Pa tij c function, equation (16a), Pa snm1n

REFERENCES Adachi, K., Yoshioka, N. and Yamamoto, K., 1973, On non-Newtonian flow past a sphere, Chem Eng Sci, 28: 2033–2043. Agarwal, P.K. and O’Neill, B.K., 1988, Transport phenomena in multiparticle systems-I. Pressure drop and friction factors: unifying the hydraulic radius and submerged object approaches, Chem Eng Sci, 43: 2487–2499. Al-Fariss, T. and Pinder, K.L., 1987, Flow through porous media of a shearthinning liquid with yield stress, Can J Chem Eng, 65: 391–405. Bird, R.B., Armstrong, R.C. and Hassager, O., 1987, Dynamics of Polymer Liquids, Vol 1. Fluid Dynamics, 2nd edn (Wiley, New York, USA). Chhabra, R.P., 1993a, Bubbles, Drops and Particles in non-Newtonian Fluids (CRC Press, Boca Raton, FL, USA). Chhabra, R.P., 1993b, Fluid flow, heat and mass transfer in non-Newtonian fluids: multiphase systems, Adv Heat Transfer, 23: 187–278. Chhabra, R.P., 1998, Rising velocity of a swarm of spherical bubbles in nonNewtonian power law fluids at high Reynolds numbers, Can J Chem Eng, 76: 137–140. Chhabra, R.P. and Raman, J.R., 1984, Slow non-Newtonian flow past an assemblage of rigid spheres, Chem Eng Commun, 27: 23–46. Chhabra, R.P. and Richardson, J.F., 1999, Non-Newtonian Flow in the Process Industries (Butterworth-Heinemann, Oxford, UK). Chhabra, R.P., Comiti, J. and Machac, I., 2001, Flow of non-Newtonian fluids in fixed and fluidised beds, Chem Eng Sci, 56: 1–27. Ciceron, D., Comiti, J., Chhabra, R.P. and Renaud, M., 2002a, NonNewtonian fluidisation of spherical particles, Chem Eng Sci, 57: 3225–3234. Ciceron, D., Comiti, J. and Chhabra, R.P., 2002b, Pressure drops for purely viscous non-Newtonian fluid flow through beds packed with mixed sized spheres, Chem Eng Commun, 189: 1403–1414. Comiti, J. and Renaud, M., 1989, A new model for determining mean structure parameters of fixed beds from pressure drop measurements: application to beds packed with parallelepipedal particles, Chem Eng Sci, 44: 1539–1545. Dolejs, V. and Machac, I., 1987, Pressure drop in the flow of a fluid through a fixed bed of particles, Int Chem Eng, 27: 730–736. Dolejs, V., Siska, B. and Dolecek, P., 1998, Modification of Kozeny-Carman concept for calculating pressure drop in flow of viscoplastic fluids through fixed beds, Chem Eng Sci, 53: 4155–4158. Gummalam, S. and Chhabra, R.P., 1987, Rising velocity of a swarm of spherical bubbles in a power law non-Newtonian liquid, Can J Chem Eng, 65: 1004–1008.

Trans IChemE, Part A, Chemical Engineering Research and Design, 2004, 82(A5): 642–652

652

DHOLE et al.

Happel, J. and Brenner, H., 1965, Low Reynolds Number Hydrodynamics (Prentice Hall, Uppersaddle River, NJ, USA). Harlow, F.H. and Welch, J.E., 1965, Numerical calculation of time-dependent viscous incompressible flow of fluid with free surfaces, Phys Fluids, 8: 2182–2188. Jaiswal, A.K., Sundararajan, T. and Chhabra, R.P., 1991, Hydrodynamics of Newtonian fluid flow through assemblages of rigid spherical particles in intermediate Reynolds number regime, Int J Eng Sci, 29: 693–708. Jaiswal, A.K., Sundararajan, T. and Chhabra, R.P., 1992, Flow of power law liquids through particle assemblages at intermediate Reynolds numbers, Can J Chem Eng, 69: 1235–1241. Jaiswal, A.K., Sundararajan, T. and Chhabra, R.P., 1993, Hydrodynamics of creeping flow of power law liquids through particle assemblages, Int J Eng Sci, 31: 293–306. Jaiswal, A.K., Sundararajan, T. and Chhabra, R.P., 1994, Pressure drop for the flow of dilatant fluids through a fixed bed of spherical particles, Can J Chem Eng, 72: 352–353. Jarzebski, A.B. and Malinowski, J.J., 1987, Drag and mass transfer in slow non-Newtonian flows over an ensemble of Newtonian spherical drops or bubbles, Chem Eng Commun 49: 235–247. Kawase, Y. and Ulbrecht, J., 1981a, Drag and mass transfer in nonNewtonian flow through multi-particle systems at low Reynolds numbers, Chem Eng Sci, 36: 1193–1202. Kawase, Y. and Ulbrecht, J., 1981b, Motion of and mass transfer from an assemblage of solid spheres moving in a non-Newtonian fluid at high Reynold numbers, Chem Eng Commun, 8: 233–249. Kemblowski, Z., Dziubinski, M. and Sek, J., 1989, Flow of non-Newtonian fluids through granular media, Transport Phenom Polym Syst, 5: 117–175. Kuwabara, S., 1959, The forces experienced by randomly distributed parallel circular cylinders or spheres in a viscous flow at small Reynolds numbers, J Phys Soc Jpn, 14: 527–532. LeClair, B.P. and Hamielec, A.E., 1968, Viscous flow through particle assemblages at intermediate Reynolds numbers—steady state solutions for flow through assemblages of spheres, Ind Eng Chem Fundam, 7: 542–549. Malleswara Rao, T.V. and Chhabra, R.P., 2003, A note on pressure drop for the cross-flow of power-law liquids and air=power law liquid mixtures past a bundle of circular rods, Chem Eng Sci, 58: 1365–1372. Machac, I. and Dolejs, V., 1981, Flow of generalised Newtonian liquids through fixed beds of non-spherical particles, Chem Eng Sci, 36: 1679–1686. Mandhani, V.K., Chhabra, R.P. and Eswaran, V., 2002, Forced convection heat transfer in tube banks in cross flow, Chem Eng Sci, 57: 379–391. Mauret, E. and Renaud, M. 1997, Transport phenomena in multi-particle systems—I. Limits of applicability of capillary model in high voidage beds. Application to fixed beds of fibers and fluidised beds of spheres, Chem Eng Sci, 52: 1807–1817.

Mohan, V. and Raghuraman, J., 1976a, Bounds on the drag for creeping flow of an Ellis fluid past an assemblage of spheres, Int J Multiphase Flow, 2: 581–589. Mohan, V. and Raghuraman, J., 1976b, A theoretical study of pressure drop for non-Newtonian creeping flow past an assemblage of spheres, AIChE J, 22: 259–264. Prasad, D.V.N. and Chhabra, R.P., 2001, An experimental investigation of the cross flow of power law liquid past a bundle of cylinders and in a bed of stacked screens, Can J Chem Eng, 79: 28–35. Sabiri, N.E. and Comiti, J., 1995, Pressure drop in non-Newtonian purely viscous fluid flow through porous media, Chem Eng Sci, 50: 1193–1201. Sabiri, N.E. and Comiti, J., 1997, Experimental validation of a model allowing pressure gradient determination for non-Newtonian purely viscous fluid flow through packed beds, Chem Eng Sci, 52: 3589–3592. Satish, M.G. and Zhu, J., 1992, Flow resistance and mass transfer in slow non-Newtonian flow through multiparticle systems, J Appl Mech ASME, 59: 431–437. Shibu, S., Chhabra, R.P. and Eswaran, V., 2001, Power law fluid flow over a bundle of cylinders at intermediate Reynolds numbers, Chem Eng Sci, 56: 5545–5554. Shukla, R., Dhole, S.D., Chhabra, R.P. and Eswaran, V., 2004, Convective heat transfer for power law fluids in packed and fluidised beds of spheres, Chem Eng Sci, 59: 645–659. Slobodov, E.B. and Chepura, I.V., 1982, A cellular model of biphasal media, Theor Found Chem Eng, 16: 235–239. Tiu, C., Zhou, J.Z.Q., Nicolae, G., Tunnan, N.F. and Chhabra, R.P., 1997, Flow of viscoelastic polymer solutions in mixed bed of particles, Can J Chem Eng, 75: 843–850. Tripathi, A. and Chhabra, R.P., 1992, Slow power law fluid flow relative to an array of cylinders, Ind Eng Chem Res, 31: 2754–2759. Tripathi, A. and Chhabra, R.P., 1994, Hydrodynamics of creeping motion of an ensemble of power law liquid drops in an immiscible power law medium, Int J Eng Sci, 32: 791–803. Tripathi, A. and Chhabra, R.P., 1996, Transverse laminar flow of nonNewtonian fluids over a bank of cylinders, Chem Eng Commun, 147: 197–212. Zhu, J., 2001, A note on slow non-Newtonian flow over an ensemble of spherical bubbles, Chem Eng Sci, 56: 2237–2241.

The manuscript was received 14 July 2003 and accepted for publication after revision 18 February 2004.

Trans IChemE, Part A, Chemical Engineering Research and Design, 2004, 82(A5): 642–652