FLUID MECHANICS IN CHEMICAL ENGINEERING
R. E. GEE and J. 6.
LYON
Polychemicals Department, E. I. du Pont d e Nemours & Co.,Inc., Experimental Station, Wilmington, Del.
Nonisothermal Flow of Viscous Non-Newtonian Fluids A mathematical model accurately describes nonisothermal laminar flow of Lucite acrylic resin through circular channels b In many chemical engineering operations, viscous non-Newtonian fluids are forced through channels or tubes which are at a temperature different from that of the fluid. Application of accepted engineering design procedures to nonisothermal laminar flow of this type of fluid i s very difficult. Effect of various factors on steady-state flow can be determined by solving the equations presented here.
SOME
of the reasons for the difficulties in applying conventional fluid flow theory to the problem of nonisothermal flow of viscous non-Newtonian fluids are : The viscosity of almost all viscous fluids changes markedly with temperature and, therefore, heat transfer affects the velocity distribution. Fluid viscosity changes with shear stress (non-Newtonian). Frictional heat generated during flow and cooling due to expansion must be included in many cases. Variations of thermal diffusivity and heat capacity with temperature must be considered. One of the earliest applications of the unsteady-state heat transfer equation involving heating or cooling of a flowing fluid was given by Graetz (7). His treatment was restricted to fluids having constant viscosity and constant thermal diffusivity flowing in undisturbed laminar motion. Yamagata (74) extended the treatment to include fluids whose viscosity varied in a simple manner with temperature. Brinkman (3) and Bird ( 2 ) further extended solutions of the flow equations to include heat generated by friction, and Toor (72) included cooling resulting from expansion during flow. In the present work, computer methods have been used to solve the energy and flow equations for non-Newtonian fluids which may be heated or cooled by heat transfer from the walls during flow. T h e variations of viscosity and thermal
956
diffusivity with temperature, the heat generated by friction, and cooling resulting from expansion have all been included. These equations have been solved for several cases of practical interest. Physical System
The physical system for bvhich the mathematical model is derived is as follows: The fluid is forced through a tube of circular cross section, in which the wall temperature is constant but may be different from that of the initial fluid temperature. The fluid enters the tube a t all times a t a uniform constant temperature, TO. The solution of the problem is started with the tube full of fluid at temperature TO. At time zero, a pressure, PO,is applied a t the inlet of the tube, and simultaneously the wall temperature, T,, is changed to the desired value. For times greater than zero, the applied pressure and wall temperature are held constant. A sketch of the system is shown in Figure 1. Mathematical Model
-dL/dx
a
=
f(~)
u = velocity
x = distance from boundary surface 7 = shear stress Newton formulated the basic concept of viscosity when he proposed that the shear rate is directly proportional to shear stress. The constant of proportionality is called viscosity, 7, which was thought to be a true constant for any liquid at constant temperature. Thus Equation 1 becomes: 1 D
= -
(7)
which is Newton’s law.
., c
Figure 1 .
AT T I M E Physical system
(2)
In later work
(RESERVO I R )
AT T I M E 0 = 0
(1 1
where
- d o /d x
Although the mathematical model developed is generally applicable to
INDUSTRIAL AND ENGINEERING CHEMISTRY
laminar nonisothermal flow of nonNewtonian fluids, the specific fluid considered in this paper is a high pollmer in the temperature range of practical interest. Viscosity-Shear-Temperature Dependency. When any fluid flows between fixed boundaries, the rate of shear that it undergoes is a function of the shear stress applied-that is, the change of velocity with distance from the boundary is a function of the shearing force:
0>0
it was found that many substances did not obey this law, and these substances were termed “non-Newtonian.” Plastic melts fall in this class. I n the case of viscoelastic fluids-Le., plastic melts-hydrodynamic theories by Bueche (4) and by Pao (9, 70) shed light on the reasons for deviations from Newton’s law. However, without regard for the fundamental mechanism causing the nonlinear dependence of shear rate and shear stress, a semiempirical function relating the two has been used, which fits the data for plastic melts very well. This function is similar to the function suggested by Dexter ( 6 ) , but has only one temperature dependent term. Substituting this function in Equation 1, instead of Newton’s law we have :
-
-du/dx
= CT(1
+k~”)
(3)
where C is a temperature-dependent term and k and n are characteristic (temperature-independent) constants for the particular fluid. This relationship, which does not require the viscosity concept, is perhaps a more technically correct way of looking a t the general flow. However, the same result can be accomplished by using Newton’s law and saying that the viscosity is a function of shear stress : At low shear stress, where k ~ isn negligible compared to 1, 7 = Q O and this becomes Newton’s law. At low stresses (where k ~ ”is negligible in Equation 3), the value of C is equal to 1 / q o in Equation 4. T h e constant QO is called the Newtonian viscosity and is a function of tempera-
Table 1. 1. Temperature
ture according to the following relationship : qo = A e E / R T
Equations 4 and 5 describe the effect of shear stress and temperature on the viscosity of plastic melts over a wide range of conditions and are thought to be applicable to most high-viscosity fluids. Velocity Distribution. The relationship for velocity distribution is given in Equation 1. However, in the case of non-Newtonian fluids, viscosity is a function of shear stress (given by Equation 4) and temperature (Equation 5 ) . This, of course, means that Equation 1 cannot be integrated by routine procedures. I n order to make the problem more specific, a tube of circular cross section has been selected as the flow path. For this case, the velocity, u, of the fluid a t any reduced radius, y, in the channel is given by the following equation (74):
where i, is the average velocity in the tube. Radial flow is neglected with respect to axial flow in the development of Equation 6. Although this equation is rigorous only for steady-state flow of an incompressible fluid, it is assumed that it approximates the unsteady-state velocity distribution if a correct material .balance is maintained a t each interval of time (quasi steady-state a t each time increment).
Equations Describing Flow of Non-Newtonian Fluids
a = ao/[l
+ a ( T - 273)] + 6(T - 273)l
cP = ( c , ) ~ [l 2.
(5)
Velocity
Temperature Distribution. When the walls of the tube are a t a different temperature than the fluid, heat is transferred radially by conduction through the walls. If radial flow of fluid is assumed to be negligible with respect to axial flow, heat is transferred along the length of the tube only by movement of the fluid. With viscous fluids, very large shear stresses are often encountered. In these cases, the temperature of the fluid rises at points of high shearing action. Accompanying the high shear stress is necessarily high pressure drop (sometimes higher than 5000 pounds per square inch). In this pressure range, most fluids are somewhat compressible, which means that cooling resulting from expansion must also be included. Variations of thermal diffusivity, a, and heat capacity, cp, with temperature have also been included in this treatment. All of these physical phenomena have been described mathematically (2, 8, 72, 74) and are given in Equation 7 :
a = a o / [ lf a ( T
cP = ( c p l 0 [I
4.
(7a)
rg)
equal to zero in Equation
7 or by continuing the incremental time calculations for a sufficiently long time.
Pressure Drop
Viscosity-Shear-Temperature Dependency
where
273)l
(7b) where the first term on the left represents the change in temperature along the tube length as fluid of changing temperature makes its way down the tube; the second term on the left expresses the possible transient conditibn of temperature with time. The first term on the right expresses the temperature change due to conduction through the walls. T h e second term on the right represents the temperature change due to localized areas of high shear, and the last term on the right expresses the cooling effects because of expansion. Linear relations (Equations 7a and 7b) have been used to express the thermal diffusivity and heat capacity dependence on temperature. Reference is made below to steadystate conditions. Here steady state is defined as the condition where the velocity and temperature distributions are independent of time (not length of tube). This is accomplished either by making
3.
-
+ b ( T - 273)l
1/11 = (1/110)(1 v o = AGEIRT
+ k+)
T h e two procedures have been checked and found to agree. Pressure Drop. The pressure drop at a given position in the tube is given by (74):
Although this equation is rigorously VOL. 49, NO. 6
0
JUNE 1957
957
derived on the basis of steady-state flow of an incompressible fluid, it is assumed that it approximates the unsteady-state pressure gradients if a correct material balance is maintained a t each interval of time in the velocity equation (which is used to obtain urnax).
265
Summary of Equations The equations that describe the flow characteristics of viscous non-Newtonian fluids (Equations 4 through 8) are summarized in Table I. Assumptions made to arrive at these equations are:
260
The flow path has a circular cross section. Conductive heat transfer in the axial direction is negligible compared to that in the radial direction. For momentum and continuity considerations, the fluid is incompressible. Radial flow is small compared to axial flow and can be accounted for in the velocity and pressure drop equations by maintaining a correct material balance at each increment along the tube. A quasi stead>--state is applicable for transient calculations. Constants in the equations have been obtained for Lucite 140 acrylic resin (Table 11). Viscosity constants ( A , E, k , and n of Equations 4 and 5) were determined from the extensive data of Tordella ( 7 3 ) . The coefficient of expansion \vas obtained from the Lvork of Spencer and Gilmore ( 7 7). Thermal diffusivity a t room temperature was taken from the results of Chung and Jackson (5). Thermal diffusivity data a t higher temperatures were determined by the present authors.
Computer Solution Because of the complexity of the above set of nonlinear partial differential equations, solutions by an analytical approach would be almost impossible. However, bvith the advent of high speed computers, particular solutions may be obtained by numerical integration methods. One of the difficulties encountered in preparing the equations for machine solution was the choice of a suitable numerical method. After many hand calculations and trials using different methods, a Gauss-Seidel iterative-type method was finally used.
Table 11. Constants Used in Computer Solution of Flow Equations A = 6.3 x 10-13p0i~ e o = 12.6 X lO-'sq. cm. per second a = 0.925 X IO-* OK.-' ( c , ) ~ = 1.02 X 10-7 ergs per gram K. b = 3.4 X 10-3 OK.-' E / R = 19,500 O K . ir = 2.4 X 10-12 (sq. cm./dyne)2 n = 2.0 z = 0.200 p = 1.09 grams per cc.
riiow
958
255
.-CENTER 2 50
0.2
0.4 y=
,
Figure 2.
(E)
a t that
position was then determined. The procedure was repeated with this same flow rate at successive intervals along the tube until the end of the tube was reached. At this point the pressure gradients were summed and compared with the total pressure drop specified in the problem. If the two agreed, the solution was satisfied; if they did not agree, a new mass flow rate was assumed and calculations were repeated until the calculated pressure drop did agree with the specified value. This entire procedure was then repeated a t the next increment of time and continued in this manner until steady-state 11 as reached, or until flow essentially stopped.
Calculated Results Shown in Figures 2 to 8 are examples of steady-state temperatures, velocities, and pressure gradients predicted for the flow of Lueite acrylic resin. The tube is l/g inch in diameter and 4 inches long, the applied pressure is 3000 pounds per square inch, and the
INDUSTRIAL AND ENGINEERING CHEMISTRY
0.8
1.0
Steady-state temperature distribution
The general procedure followed on the computer to obtain a solution for a specified applied pressure was to assume a mass flow rate and then calculate the temperature, viscosity and velocity distribution across the tube a t the entrance. The pressure gradient
0.6 r/R,
li
initial plastic temperature is 252" C. At this temperature, Lucite has a Newtonian viscosity of about IO4 poises, Lvhich increases by a factor of 2 for a temperature decrease of about 3" C. For the case where the wall temperature is approximately equal to the initial fluid temperature (Figure 2), the position of the maximum temperature, resulting from frictional heat generation, shifts slightly toward the center of the tube for increasing distances down the tube. I t can also be seen in Figure 2 that rhe radial position corresponding to the average flow temperature-Le., temperature averaged according to the flow rate--lies between 0.62 and 0.66. This value is in agreement with the measured value of 0.6 reported by Beyer and Dah1 (7). In Figure 3 the effect of wall temperature may be seen. Of interest here is the prediction that the wall temperature must be reduced about 60" C. before the center temperature is reduced below its initial value. This large temperature difference is necessary before it is "felt" in the center because of the combination of frictional heat generation and the increase in velocity near the center line with cooling at the wall. Figure 4 shows average flow temperatures and Figure 5 shows velocity profiles. In Figure 5, it can be seen for
FLUID M E C H A N I C S 260 v
a: 5 W
t
9 LL
I
I
T, =250°C.
1I 1
3
v ,
I
250
d W
2 240 0
1.0
I
0.4
0.6
I
0.8
I
1.0
z = 1 /Lo
I .o
0
I
I
0.2
Figure 4.
Steady-state average flow temperatures
y = r/Ro Figure 3.
I n the cold wall cases (Figures 5 and 6) the center velocity increases significantlv with distance down the tube. F o r this to occur, the pressure drop per unit length must be increasing along the tube, since the plastic is cooling as to flows (Figure 7). Similar results have been calculated for several tube diameters and lengths, a t various applied pressures, initial temperatures, and wall temperatures.
Steady-state temperature distribution
the "isothermal wall case ( T , = 250' C.) that the velocity distribution does not change appreciably with distance down the tu.be, even though there is a significant change in the temperature distribution and pressure gradient. This
occurs because these two effects tend to balance one another-the viscous heating decreases the viscosity near the wall and the decrease in the pressure gradient (and, therefore, shear stress) increases the viscosity near the wall.
1001
Comparison of Calculated and Experimental Results
Tw.i '
250OC
250°C.
20-
r
I
To check the validity of the mathematical model, experiments were carried I
0 0)
Tw =
rw= 200%
200°C
IO
B
k
0
sw
> 5
V- cm./sec. Tw = 160°C.
Tw = 180°C /
0.5
1.0 Figure
CENTER y - r/Ro
0.5
5. Steady-state velocity distribution
I .o
0
Z '0f/Lo Figure.6.
Steady-state velocity distribution VOL. 49, NO. 6
JUNE 1957
959
z
I= t
W J
R, = 1/16 IN.
5000
a z
0 t
4000
M
a \ L .-
a u)
3000
I I I
5 w0
2000
a
a a w
a
3
I000
v) 0)
W
U
n
b
I 0.12
0.14
016
018
i -
110
1.0
2 =.$/Lo
Figure 7.
EXPER. FLOW RATE
Acknowledgment
cy
The authors wish to thank J. K. Beasley and J. A. Cipolla for preparing and programming the equations for computer solution, and P. H. Squires for assistance in deriving the mathematical model.
The mathematical model accurately describes nonisothermal laminar flow of viscous non-Newtonian fluids through circular channels. The effect of pressure drop, fluid temperature, wall temperature, channel geometry, and fluid properties on the steady-state flow can be determined from solution of these equations. Solutions of the equations for transient conditions, although not proved, are felt to be valid. INDUSTRIAL A N D ENGINEERING CHEMISTRY
= = = =
=
= = =
= = = =
= = =
temperature coefficient of thermal diffusivity, ' K.-I Xewtonian viscosiiy extrapolated to 0' K., poises temperature Coefficient of heat capacity, ' K.-I heat capacity of fluid at constant pressure at 0 ' C., ergs/gram ' K. heat capacity of fluid at constant pressure at temperature T , ergsjgrams K. activation energy of flow, cal./ gram mole length of tube, cm. distance along tube to point under consideration, crn. weight of fluid leaving tube, grams pressure, dynes per sq. cm. gas constant, cal./gram mole O K. radius of tube, cm. radial distance to point under consideration, cm. temperature of fluid, K. average flow temperature = 2 J1 VT Y dY __U
= =
wall temperature, C. velocity, cm. per second
dmldb' average tubevelocity = __ .rrR:P = distance from boundary surface = reduced radius = r/Ro' = reduced length = Z/Lo =
cy0
i
vo 11
Nomenclature
Conclusions
960
Figure 8. Comparison of calculated and experimental flow rates
Steady-state pressure gradient
out in which molten Lucite acrylic resin was caused to flow through tubes cooled on the outside with circulating oil. Tube sizes ranged from ' / 8 inch in diameter by 4 inches long to ' / 4 inch in diameter by 16 inches long. The pressure range invesrigated was SO0 to 3000 pounds per square inch, the wall temperature range was 150' to 250' C., and the initial melt temperature was 250' C. Computer solutions were then obtained for the same conditions and comparisons of experimental and calculated flow rates were made (Figure 8). From the excellent agreement of all points. it has been concluded that the assumptions made are justified, and the numerical method of solution is correct. As local velocities for the material studied are very sensitive to point temperatures, and the calculated integrated velocitvLe., flow rate-agreed well with experimental results, it is felt that the predicted velocity and iemperature profiles are also essentially correct. Agreement of the predicted position corresponding to the average flow temperature with the reported value of Beyer and Dahl (7) tends to substantiate this belief
IO
- - g/sec.
0 p T
thermal diffusivity of fluid at temperature T , sq. cm. per second = thermal diffusivity of fluid at 0 ' C.: sq. cm. per second = coefficient of expansion at Tflo.,,, K.-l = Newtonian viscosity of fluid, poises = non-Xewtonian viscosity of fluid, poises = time, seconds = density of fluid, grams per cc. = shcar stress, dynes per sq. cm. =
Literature Cited (1) Beyer, C. E., Dahl, R. B., Modern PlQStiCS 30, 124 (1952). (2) Bird, R. B., SPE Journu! 11, 35 (1 955). ( 3 ) Brinkman,, H. C., Afipl. Sci. Research A-2, 120 (1951). (4) Bueche, F., J. Chem. Piiys. 22, 1570 (1954). (5) Chung, P. I
