Ind. Eng. Chem. Fundam. 1083, 22, 183-186
183
Fully Developed Flow of Power-Law Fluids in Ducts TaJo Liu Research Laboratwles, Eastman K&k
Company, Rochester, New York 14650
The Galerkin flnite-element method has been applied to determine the pressure drop/flow rate relationship for fully developed flow of power-law fluids in a duct of arbRrary cross section. The duct shape can be characterlted by three geometric parameters, which can be evaluated accurately without experimentation. The pressure drop/fbw rate relationshipis completely determined if the three geometric parameters are known. Miller's empirical approach for estimating pressure losses is examined and found to be reasonably accurate provided that the power-law index is close to unity.
Introduction The accurate prediction of the pressure drop/flow rate relationship for non-Newtonian flow in ducts is of great importance to many industrial applications (Miller, 1972). Several researchers have studied the non-Newtonian flow behavior in ducts of arbitrary cross section. Schechter (1961)applied a variational method to solve for the velocity profiles of power-law fluids in a rectangular duct and then determined the relationship between the friction factor and the Reynolds number. Wheeler and Wissler (1965) solved the same problem by the finite-difference method and verified their predictions experimentally. Mizushima et al. (1965) and Mitsuishi et al. (1968) took a variational approach to study non-Newtonian flow in elliptical and triangular ducts. The methods used by the above authors are either too awkward or impossible to generalize for ducts of unusual cross section. Kozicki et al. (1966) and Kozicki and Tiu (1967) proposed a new approach to predict pressure losses. Their method requires the determination of two geometric parameters and a function of shear stress to characterize the behavior of the fluid model. The geometric parameters are available for limited duct shapes. It requires the Newtonian flow data to evaluate the geometric parameters for complex-shaped ducts. Miller (1972) suggested a much simpler empirical method. According to his analysis, there exists a universal curve if the average shear rate is plotted against the average shear stress for flow of a given non-Newtonian fluid in ducts of arbitrary cross section. Such a curve can be obtained by performing the flow experiment in a duct of simple geometry, such as a circular duct. To obtain the pressure drop/flow rate value for a particular duct from this curve, it is necessary only to know one shape factor for this duct. The shape factor is defined as the product of the friction factor and the Reynolds number for Newtonian flow in this duct. Hanks (1974) pointed out that if the shear rate is higher than a critical value or the fluid under consideration has significant yield stresses, Miller's method is no longer valid. In this paper, we describe the determination of the pressure drop/flow rate relationship for power-law fluids by the Galerkin finite-element method, which is applicable to duds of any geometry. Such a relationship depends on three parameters, which characterize the duct geometry; they can be computed if the finite-element solution of the flow model is available. We also rigorously examine the validity of Miller's method; his predictions for pressure losses are reasonably accurate for fluids with power-law index close to unity. 0196-4313/83/1022-0183$01.50/0
Flow Model We consider fully developed laminar flow of non-Newtonian fluids in a straight duct of arbitrary but constant cross section as shown in Figure 1. The cross-sectional area of the duct is D and the characteristic length of the duct h is defined as h2 E D. The flow is assumed to be isothermal and incompressible. It is generally agreed that the secondary flows in ducts are weak and have a very small effect on the pressure drop/flow rate relationship (Walters, 1979). Since we consider only the inelastic power-law fluids, we assume there is no secondary flow and the velocity is purely axial. The equations of motion reduce to a single equation
The rheological equation for power-law fluids is (Bird et al., 1960) = -ijA where A is the rate of the deformation tensor and
(2)
z
(g)]
=m[(g)'+ where m and n are material constants. ij
T,,
= -si
Tyz
= -ij(
(3) T,,
and
rYzare
$)
Substituting eq 4 into eq 1, we obtain
+
=
!dz E
The dimensionless variables are
Substituting eq 6 into eq 5, we obtain
where It =
[
( 5 )] z
+
0 1983 American Chemical Society
(8)
184
Ind. Eng. Chem. Fundam., Vol. 22, No. 2, 1983 i
t
Figure 1. Fully developed flow in a duct. D = constant cross-sectional area; z = axial direction; = axial velocity component.
and D is the dimensionless domain whose area is equal to 1.
The no-slip boundary condition requires
w = 0 (on dD)
(9)
where dD is the boundary of D. The volumetric flow rate Q in the duct is defined as
Introducing the dimensionless variables into eq 10, we obtain
where the shape factor X(n) is given as
X(n)
D
lw
dx dy
Figure 2. An isoparametric element and its correspondence in the local coordinates. Table 1. Geometric Parameters in Eq 19 for Various Duct Geometries geometry
a
1. circular 2J;; 3.7471 2. square (Wheeler and Wissler, 1965) 3. semicircular 4.1798 4. equilateral triangular 4.3821 4.5559 5. 0.25 of an annulus (radius ratio = 0.5) 6. L-shaped 4.5409 (which consists of 3 squares)
b
C
J;;
3J;; 1.7330 5.8606
1.3344 6.6023 1.0389 6.8785 1.0574 7.4281 1.5485 7.1521
interpolation function is a complete quadratic polynomial, and we have
6
(12)
Y = C4iYi i=l
which depends on the duct geometry and the power-law index n. Equation 11 expresses the pressure drop/flow rate relationship for flow of power-law fluids in ducts. We can predict pressure losses from this equation if X(n) is known. If we define the Reynolds number Re and the friction factor f as ,,Q2-nh3n-4
Re
m
2h6( -
(13)
$)
f =
pQ2S where S is the wetted perimeter, we obtain, after comparing the product o f f and Re with eq 11
Finite-Element Solution To determine X(n),we shall solve eq 7-9 for w and then perform the integration in eq 12. The Galerkin finiteelement method was adopted to solve eq 7-9 because of its flexibility in handling the irregular domain of integration. Detailed descriptions of this method are in various textbooks (Heubner, 1975; Zienkiewicz, 1977; Chung, 1978). The isoparametric element we use and its correspondence in the local coordinates (E,{) are shown in Figure 2. The
After the isoparametric element and the interpolation function are specified, it is possible to follow the procedure explained in the textbooks mentioned previously to solve eq 7-9. Since eq 7 is nonlinear, we use Newton's method to determine w iteratively. The efficient frontal technique (Hood, 1976) is applied to solve the overall system of equations. Once the convergent solution for w is obtained, we can determine the local shape factor X(")(n)for each element
where e is the area of the element. By substituting eq 16a into eq 18, we can perform the integration in the local coordinates. The value of X(n) is obtained by summing the contributions of all elements. Example Calculations We checked the accuracy of our finite-element solutions by solving for flow (a) in a circular duct and (b) in a square duct, so that we could compare our solutions with the exact solution available for (a) and with the numerical solution
Ind. Eng. Chem. Fundam., Vol. 22, No. 2, 1983
185
Table 11. Values of h ( n ) for Four Duct Geometries geometry 3
4
5
6
n
Aa
Bb
A
B
A
B
A
B
2.0 1.5 1.0 0.7 0.4
0.67252 0.51438 0.30143 0.15233 0.028060
0.67285 0.51440 0.30144 0.15233 0.028171
0.64491 0.49309 0.28824 0.14486 0.026276
0.64572 0.49323 0.28823 0.14487 0.026253
0.58759 0.44710 0.25869 0.12818 0.022491
0.58876 0.44737 0.25865 0.12820 0.022410
0.55021 0.41509 0.23679 0.11555 0.019569
0.55038 0.41514 0.23678 0.11555 0.019580
a A = h ( n ) X 10, h(n)computed by the finite-element method. constants in Table I.
a least-squares fit. We found that eq 19 is valid for the various duct geometries studied. The constants were determined by a least-squares procedure (Kuester and Mize, 1973, p 251) with eight values of X(n) for n between 1.3 and 0.6. The values of these geometric parameters are listed in Table I. Values of X(n) computed from eq 19 agree well with those obtained directly by the finite-element solutions over a wide range of n. Some representative results are shown in Table 11. Fredrickson and Bird (1958) studied fully developed flow of power-law fluids in an annulus. We found that h(n)for an annulus cannot be arranged exactly in the form in eq 19. However, the errors are small if we approximate the exact solution of X(n)by eq 19. The three parameters were evaluated by the same least-squares procedure with six values of X(n), ranging from n = 2 to n = 0.4, for various aspect ratios K . The three geometric parameters are tabulated in Table 111. Values of X(4) and X(0.2) computed by eq 19 were compared with those obtained by Fredrickson and Bird (1958) and with the approximate variational results by Bird et al. (1977, p 230). The results (Table IV) show that our approximations of X(n)are quite accurate.
Table 111. Geometric Parameters in Eq 1 9 for an Annulus K
a
b
c
0.01 0.1 0.3 0.5 0.7 0.9
3.6463 3.8375 4.7851 6.1184 8.4339 15.2898
1.3810 3.2173 4.5983 6.0412 8.3781 16.3777
7.3481 7.9237 9.6916 12.2920 16.8863 30.7895
B = h ( n ) x 10, h ( n )computed by eq 1 9 using the three
obtained by Wheeler and Wissler (1965) for (b). Owing to the nondimensionalization,the area of D is always equal to 1. By using a grid with 40-50 elements for a duct geometry, we can obtain the numerical solution accurate to the fifth decimal place. For a given duct geometry, we start the computation for w and X(n) by choosing n = 1and then decrease (or increase) the power-law index n and repeat the computation consecutively using the previous convergent solution as the initial guess. Therefore, we can generate a sequence of X(n) in one computer run. Wheeler and Wissler (1965) reported the difficulty of obtaining a convergent solution for power-law fluids of small values of n. Owing to the steep velocity gradients close to the duct boundary, we also fail to produce a convergent solution if n is less than 0.4. Higher order interpolation functions are required if such a solution is desired. Wheeler and Wissler (1965) devised an equation that expresses X(n) for a square duct as follows
Validity of Miller’s Method The average shear stress 7 and the average apparent wall shear rate 7 used by Miller have the following expressions
dP
h2
+ c)]-l
X(n) = [a’/.(!
y=-
where a, b, and c are constants that depend only on the shape of the duct. They obtained the three constants by
8
h2SX(1)
Table IV. Values of h ( n )for an Annulus
K
0.01 0.1 0.3 0.5 0.7 0.9
h(4 1 Fredrickson and Bird e t al. Bird (1958) (1977) 0.094920 0.082156 0.062428 0.046083 0.030916 0.014507
0.090237 0.080601 0.062056 0.045984 0.030902 0.014507
eq 19
Fredrickson and Bird (1958)
0.094063 0.081860 0.062367 0.046067 0.030912 0.014497
0.10826 X 0.49638 X 0.12156 X 0.27413 X 10-’ 0.39880 X 0.10507 X lo-’
h(0.2) Bird et al. (1977) 0.67798 0.39430 0.11239 0.26663 0.39563 0.10496
eq 1 9
X X X
X X X
lo-’
0.10885 X 0.50048 X 0.12197 X 0.27443 X 0.39841 X 0.10621 X
loT4 lo-’ lo-’
Table V. Values of R ( n )for Various Duct Geometries geometry
a
n
1
2
3
4
2.0 1.5 1.2 0.9 0.7 0.5 0.3
1.1428 1.0909 1.0435 0.97297 0.90323 0.80001 0.63160
1.0925 1.0590 1.0283 0.98237 0.93676 0.86914 0.75994
1.1020 1.0660 1.0320 0.97962 0.92534 0.84022 0.68758
1.0493 1.0320 1.0156 0.98997 0.96295 0.91967 0.83963
Geometries 7 and 8 correspond to annuli with
K
5 1.0885 1.0577 1.0282 0.98186 0.93277 0.85349 0.70408
= 0.9 and 0.1, respectively.
6 1.0816 1.0528 1.0256 0.98369 0.94022 0.87196 0.74934
7a 1.2037 1.1269 1.0596 0.96398 0.87441 0.75013 0.56619
Ba 1.1565 1.0988 1.0469 0.97118 0.89793 0.79244 0.62752
186
Ind. Eng.
Table VI.
Chem. Fundam., Vol. 22, No. 2, 1983
Ratio of for Various Duct Geometries (?i/?, ) geometry
n 2 .o 1.5 1.2 0.9 0.7 0.5 0.3
2
3
4
5
6
7
8
1.0942 1.0455 1.0178 0.99139 0.97481 0.95941 0.94602
1.0754 1.0353 1.0133 0.99389 0.98321 0.97578 0.97485
1.1864 1.0868 1.0330 0.98453 0.95618 0.93268 0.91814
1.1023 1.0475 1.0178 0.99185 0.97772 0.96817 0.96794
1.1165 1.0548 1.0209 0.99019 0.97230 0.95786 0.95001
0.90141 0.95244 0.98181 1.0084 1.0230 1.0327 1.0333
0.97650 0.98927 0.99609 1.0017 1.0041 1.0048 1.0019
Substituting eq 20-21 into eq 11, we obtain after some rearrangement
where R ( n ) is defined as
Therefore the validity of Miller's method depends on whether eq 22 is universal, which requires that the values of R(n) are close enough for ducts of different shapes. Values of R(n) were computed from eq 19 for eight geometries, and the results are given in Table V. Generally, the difference between R ( n ) values is more significant as n is remote from unity. It is also interesting to examine how large the variance will be if we try to predict pressure losses from a universal curve based on the flow data from a circular duct. Considering two duct geometries A and B, we obtain from eq 22 the exmession
Letting duct B be a circular duct, we can compute the ratio of 7 for the duct geometries in Table V. The results are listed in Table VI. In general, Miller's predictions are accurate over a range of n from 1.5 to 0.3; the errors are less than 9.0% even for the worst case, Le., geometry 4, which corresponds to the equilateral triangular duct. However, large deviations appear as n becomes greater than 1.5. Conclusion We have discussed the application of the Galerkin finite-element method for the determination of the pressure drop/flow rate relationship for fully developed flow of power-law fluids in a duct of arbitrary cross section. The method we propose here is superior to the previous methods in two aspects. (1)The Galerkin finite-element method is easily applicable to any duct geometry. The duct geometry can be characterized by three geometric parameters. It is necessary only to evaluate the three parameters by using the values of X(n) obtained by the finite-element method. (2) Pressure losses in ducts can be estimated from eq 11 with high accuracy. There is no need to evaluate any geometric parameters experimentally. Initial study of this method may be somewhat complicated. However, once the computer program for the flow model is available, only minor effort is needed to set up a grid for any given duct geometry and to perform the computation. Miller's observation was verified, and according to our calculations, his prediction is reasonably accurate, especially if n is near unity.
Nomenclature
a, b, c = geometric parameters, eq 19 D = cross-sectionalarea of the duct D = cross-sectionalarea of the duct, dimensionless dD = boundary of D e = area of a particular element f = friction factor, eq 14 h = characteristic length of the duct m,n = material constants of the power-law model, eq 3 P = fluid pressure Q = volumetric flow rate R(n) = defined by eq 23 Re = Reynolds number, eq 13 S = wetted perimeter of the duct @ = axial velocity component 10 = axial velocity component, dimensionless Z,7, z = Cartesian coordinates n, y = lateral Cartesian coordinates, dimensionless Greek Letters
9 = average shear rate, eq 2 1 A = rate of the deformation tensor [, F = local coordinates i j = apparent viscosity, eq 3 TI = apparent viscosity, dimensionless, eq 8 K = ratio of radius of inner cylinder to that of outer cylinder A(n) = shape factor, eq 12 X(e)(n)= local shape factor in each element, eq 18 p = fluid density T = stress tensor 7 = average shear stress, eq 20 T ~ T,,~ ~ =, components of T 4 = interpolation function Subscripts A, B = different duct geometries i = index for nodes
Literature Cited Bird, R. B.; Stewart, W. E.; Lightfoot, E. N. "Transport Phenomena"; Wiiey: New York, 1960. Bird, R . B.; Armstrong, R. C.; Hassager, 0. "Dynamics of Polymeric Liquids: Voi. 1 Fluid Mechanics"; Wiley: New York, 1977. Chung, T. J. "Finite-Element Analysis in Fluid Dynamics"; McGraw-Hill: New York, 1978. Fredrickson, A. G.; Bird, R. B. Ind. Eng. Chem. 1958, 5 0 , 347. Hanks, R. W. Ind. Eng. Chem. Fundam. 1974, 73, 62. Heubner, K. H. "The FlnitaEiement Method for Engineers"; Wiiey: New York, 1975. Hood, P. Int. J . Num. Meth. Eng. 1976, 70, 379. Kozicki, W.; Chou, C. H.; Tiu, C. Chem. Eng. Sci. 1966, 27,665. Kozicki, W.; Tiu, C. Can. J. Chem. Eng. 1967, 45, 127. Kuester, J. L.; Mize, J. H. "Optlmization Techniques with Fortran"; McGrawHill: New York, 1973. Miller, C. Ind. Eng. Chem. Fundam. 1872, 7 7 , 524. Mitsuishi, N.; Kitayama, Y.; Aoyagi, Y. I n t . Chem. Eng. 1966, 8 , 168. Mizushima, T.; Mitsuishi, N.; Nakamura, R. Kaguku Kogaku 1965, 28, 848. Schechter, R. S. AIChE J. 1961, 7 , 445. Waiters, K. J . Non-Newtonian Fluid Mech. 1979, 5 , 113. Wheeler, J. A.; Wissler, E. H. AIChE J. 1965, 7 7 . 207. Zlenkiewicz. 0. C. "The Finite-Element Method"; McGraw-Hill: London, 1977.
Received for review January 6, 1982 Revised manuscript received December 16, 1982 Accepted January 5, 1983
