Ind. Eng. Chem. Res. 1991,30,517-524
517
GENERAL RESEARCH Viscoplastic Flow in Centered Annuli, Pipes, and Slots Edmund J. Fordham, Simon H. Bittleston,* and M. Ahmadi Tehrani Schlumberger Cambridge Research Ltd., P.O. Box 153, Cambridge CB3 OHG,U.K. A general method for the practical calculation of viscoplastic flows in centered annuli of arbitrary radius ratio is presented. The analysis is worked for the Casson, Herschel-Bulkley, and Robertson-Stiff rheological models, but the method is readily extended to other viscoplastic rheologies. A robust algorithm for the practical calculation of pressure drops and velocity profiles from given volumetric flow rates is presented, differing from other calculations where the pressure drop is the given parameter. Illustrative calculations are given for the Herschel-Bulkley model, together with a simple experimental test of pressure drop results for a particular case. Previously published results (also readily extended to other rheologies) for the simpler but related pipe and slot channel geometries are reviewed with some extensions for comparison with the annulus case. 1. Introduction Calculations of the fully developed,llaminar,axial flows of viscoplastic fluids through simple channels such as circular pipes, wide plane slots, and centered annuli of arbitrary aspect ratio are commonly needed in the design of equipment handling such materials, and they are important reference cases for studies of more complicated flows. Restricting our attention to fluids of the “generalized Newtonian” class, Le., non-Newtonian fluids possibly exhibiting a yield stress T~ but without thixotropy, we present in this paper a general method for the calculation of flow rates, pressure drops, and velocity profiles for such fluids flowing in a centered annulus of arbitrary aspect ratio. The calculations are illustrated by reference to the rheological models of Casson (1959)and Herschel and Bulkley (1926)and the less well-known model developed for cement slurries by Robertson and Stiff (1976), but are readily reworked for any generalized Newtonian rheological model where the shear rate i. may be expressed as a function of shear stress T in closed form. In particular the special cases of the pseudoplastic (“power law”-where T~ vanishes) and Bingham plastic rheologies (these cases were considered in the classical paper of Frederickson and Bird (1958))follow naturally. Such calculationshave been previously published (asgraphical results) for the Casson rheology by Shul’man (1970)and Jaisighani (1974)(reported by Shah (1980))and for the Herschel-Bulkley rheology by Hanks (1979)(the last missed by the review of Bird et al. (1983)).However, these results are specific to particular rheological models, and they take a pressure gradient P = -dp/dx (or head gradient if not in horizontal flow) as a given quantity and compute from P velocity profiles and volumetric flow rates Q. In experimental and design work it may well be Q that is the given parameter and a pressure drop and velocity profiles that are required to be calculated. We shall call these the “direct” and “inverse”problems, respectively. For example, volumetric flows may be easier to measure directly and accurately than pressure drops, particularly with pressurized systems,
* Author to whom correspondence should be addressed. 0888-5885/91/2630-0517$02.50/0
or Q may be specified in a design program or may be determined by a positive-displacement pump. The “inverse” problem is then the relevant calculation. If families of pressure drop/flow rate characteristics are calculated by the “direct” method, it is of course true that the “inverse” problem may then be solved by graphical interpolation (assumingthat a relevant region of parameter space has been covered by the “direct” calculations). However, it is much more efficient to have a turnkey method for practical solution of the “inverse”problem, and possibly essential if such calculations are embedded in iterative design work or the analysis of more complicated flows. An efficient computerized version of graphical interpolation of design charts generated by the “direct” calculation would in effect reduce to our “inverse”problem. A robust practical algorithm for the “inverse”problem and a pseudoprogram as a coding outline are included in our analysis. The principal computational difficulty is the iterative search for the positions of the two yield surfaces in the flow. The major practical utility is thus for flows where a significant plug flow region exists, corresponding to large values of the “plasticity number” as used, e.g., by Anshus (1974)in analyzing the Bingham plastic case. Although simpler methods are available where the plug width may be neglected, the method is nevertheless general for any nonzero plug width. The simpler slot and pipe geometries may be analyzed similarly and are summarized in Appendices 1 and 2. Several of these results are already in the literature (sources but not formulas are given in the review article of Bird et al. (1983))but are reviewed here for ease of reference and comparisons with extreme cases of the annulus geometry. The generalization to other rheological models in the pipe and slot geometries should also be obvious. We hope that this article will thus serve as a practical review and design guide for calculations involving such fluids in simple channels. We emphasize the generality of the analysis and its ease of extension. We do not give detailed design charts, since these can easily be generated by anyone with access to a personal computer equipped with some numerical software. Some graphical results for 0 1991 American Chemical Society
518 Ind. Eng. Chem. Res., Vol. 30, No. 3, 1991 10.0,
I
n 1.0
0.8
0.6 0.4
Yield
SU
0.2
o . o ~ " ' ' ~ ' " ' ~ " ' ' ~ " 0.8 "' " . J 0.2
0.0
0.4
Zero stress surface
1.o
0.6
r0
dr) -tTo
n
0.8
Figure 2. Annular (and pipe) flow notation. 0.6 0.4
0.2
0.0 0.0
0.2
0.4
0.6
2. Analysis: Centered Annulus Geometry Figure 2 shows a sketch of the velocity and shear stress radial profiles for axial laminar flow in an annulus (not necessarily narrow); essential notations are also given on the figure. We normalize most radii in the analysis by the outer radius r2 and write such normalized radii in upper case R = r / r 2 ,R1 = r 1 / r 2 ,etc. 2.1. Stress Profile and Yield Surfaces. We summarize here the analysis given, e.g., by Frederickson and Bird (1958)for a Bingham plastic fluid although in fact it is general for any generalized Newtonian fluid with a yield stress. For flow under a head gradient P (positive quantity), the radial shear stress profile is given by a force balance and an elementary integration as
1.o
0.8
WWr3
n
10.0
d 0.4
7.5
3 CT:
0.2
3
zI=
2.5
0.01, 0.0
"
I
0.2
'
'
'
'
I
0.4
"
"
0.6
'
'
'
"
0.8
'
'
'
the Herschel-Bulkley model are however included as illustrations (Figure l). These are in effect equivalent to the charts of Hanks (1979),although for different parameter ranges and with different nondimensionalizations, but similar charts can be generated for other rheologies and parameter ranges with essentially the same software.
I
1.o
WWr3 Figure 1. Nondimrneional prreaure gradient VI) nondimeneional flow rate for a Hernchel-Bulkley fluid. R1 = 0.2 (a), 0.4 (b),0.6 (c) and 0.8 (d). Valuee of n indicated.
where A = 2r0/Pr2and Ro is the zero stress radius (a constant of integration and thus far undetermined). We assume throughout that P is sufficient to overcome the yield stress at the walls so that fluid flows; where 1~(R)l > T ~ a, plug flow region exists bounded by two yield surfaces. The inner and outer yield radii are determined by the condition 17(R+)l = 70; solving for Ri, one finds that the zero stress radius is the geometric mean of the yield radii R+R- = Ro2 (2) and that the plug flow width is R+-R-= A (3) A further condition is required to solve for Ro and R,; this is obtained by matching the velocities a t the edge of the plug region. 2.2. Velocity Profiles and Volumetric Flows. We distinguish in our notation between the (signed) axial velocity gradient du/dr and the unsigned "shear rate" *(r) = Idu/drl which is related to the local shear stress magnitude IT(r)l by a rheological model. Clearly, du/dr = a g n ( d r ) )j+); simple integration gives u(r) once j+) is found.
Ind. Eng. Chem. Res., Vol. 30, No. 3, 1991 519 Applying no-slip boundary conditions at the walls and continuity of u at the yield surfaces gives (4)
which is the condition to be solved with ( 2 ) and (3) for Ro and R*. This can be accomplished by a combination of root-finding and a numerical quadrature if the pressure gradient P (and hence A) is given. Volumetric flow is then given by
(Ro2and A being related to Ri by (2) and (3)). The prefactor in (8) is written thus to ensure that the functions f ( R ) are well behaved at all values of the dimensionless plug width A, and r/Aq being independent of 70 and monotonic in P, to clarify the limiting cases of pseudoplastic or Newtonian fluids (rO 0) or of high pressure gradient flows where A 1 ro (Casson)
+
+ A+" = A'(+ + C)B
171 171
= r0
> r0 (Herschel-Bulkley) for 171 > r0 = A'CB for
171
(Robertson-Stiff) (6) which can be solved for + in terms of
7
as follows:
where m = l / n and b = 1/B, and I'cc) = rO/qo1 I'(HB)= (so/A)"', and I'(R9)= C are characteristic shear rates (inverse time scales) for the three rheologies. Note that for the Robertson-Stiff model 70 = A'CB. Substituting from (l), the radial profiles +(R) can be written generally as
c'
-f(R;R-,R+) *)I'A forR-5 R s R+
sd, R2f(R;R-,R+)dR -
OT
R,5R51
(8)
where r takes the form relevant to the model, q = 1 for the Casson model, q = m for the Herschel-Bulkley, q = b for the Robertson-Stiff; and the dimensionless shear rate profiles f ( R ) are given by
R2f(R$-,R+) dR (10)
R1
where Q* = QAq/rI'r$; note that Q* thus defined is independent of ?,, 3. Calculations 3.1. Calculating Flows Q from Given Pressure Drop P. Equations 10, 4,3,and 2 are the basic results of the flow analysis, with the functions f ( R ) chosen according to the rheological model. In the special case of Bingham plastic, the integrations in (10) and (4) are analytical. Calculations given for the Casson rheological model, e.g., by Shul'man (1970), use the fact that some terms of (10) and (4)can be integrated in closed form, but there is little point in writing these out if other terms remain for numerical computation. This would also be contrary to the objective of our analysis, which was to develop a code for practical calculation that did not depend on particular analytical features of a given rheology. The calculations of flows Q from given pressure drops P is straightforward; (10) gives the flow provided the yield radii R , are known, which are in general found from (3) and (4) by using standard quadratures and root-finding methods. 3.2. Calculation of Pressure Drops P from Given Flows &. This is the more awkward problem, as (4) is not sufficient to determine the yield radii R , when P is not known. The yield radii have to be determined in general by a numerical solution of (4)and (10) simultaneously (except for pseudoplastic flows where R* are equal). We present below a method that is robust and general for the centered annulus geometry. 3.3. A Robust Algorithm for the "Inverse"Problem. The computational problem reduces to finding simultaneously the zeroes of two functions of Ri:
B(R-,R+)= Q* - 1 1 R 2 f d R+ L R - R 2 f d R (12) R+
forR1SRSR-
JR-
1
(from (4)and (lo),respectively) where f means f(R;R-,R+) as given by (9), ( 2 ) , and (3). We use the robust root-finding method of Brent (1973) as described by Press et al. (1986) for finding the zeroes of 3 and 5'; two iterations were nested for R , simultaneously. Other interval-chopping algorithms could also be used. A two-dimensional version of the Newton-Raphson method was also tried, but experiment showed that convergence was not guaranteed. Brent's method, like similar interval-choppingmethods, requires outer bounds that are known to bracket the required root; however, in our nested procedure we must also make sure that all iterates generated in the outer iterations allow the inner iteration to work from bounds that are certain to bracket its root. This can in fact be done as follows. The outer iteration finds a zero of 8, iterating on R-, with R+ being determined by the inner iteration which delivers
520 Ind. Eng. Chem. Res., Vol. 30, No. 3,1991
a value for R , such that 3 = 0, given an iterate R-. The critical step is to supply correct outer bounds for both levels of iteration. These are for the outer iteration
R- =
(:
for which 9 1 0 for which 9 C 0
(13)
where R , is given by
R , is physically the limit to which Rt both tend as the flow rate becomes indefinitely large. In practice we take R , t (where t 0
< R,
(15)
The condition R- C R. is certainly required but is automatically satisfied by the bounds we have selected for the R- iterates. Technical demonstrations that the bounds in (13) and (15) correctly bracket the roots sought are given at Appendix 3. Once R* are found, A, Ro, and hence P and the profiles T ( R )and q ( R ) are known; the latter can be immediately integrated for the full velocity profile u(r). 4. Discussion 4.1. Pressure/Flow Curves for Herschel-Bulkley
Fluids. The method outlined above has been used to generate families of curves of nondimensional pressure (or head) gradients vs volumetric flow rates, for various Herschel-Bulkley rheologies and annulus radius ratios. These are shown in Figure 1. We assume throughout that the flows remain laminar and fully developed. Transition to turbulence is not addressed in this paper, and dimensional resulta shown should have sufficientlylow Reynolds numbers that the flows are laminar. Inspection of (10) and other results shows that, with the nondimensionalizations chosen for P and Q, the only free parameters are the flow behavior index n (6) and the radius ratio R1= rl/r2. The other rheological parameters appear in the dimensionless groups (70 only with P, all three rheological parameters (as with Q). Q* is not convenient for direct presentation of resulta that are eventually required in dimensional form, because of the Aq factor: we therefore plot dimensionless ~ ~nondimensional . pressure flow rates as Q / ~ r rThe gradient Pr& - R J / ~ =T (1 ~ - R J / A plotted on the figures includes a factor (1 - R , ) to ensure that its value at the starting flow is unity for all cases, which is presentationally convenient. (1- R l ) / Ais also, by (3),the reciprocal of the fractional plug flow width a c r w the channel. The annulus radius ratio is varied from frame to frame; a family of curves is given for various flow behavior indices n. For a known geometry and rheology, the pressure gradient P developed by a given flow Q (or vice versa) may be estimated from these curves by graphical interpolation; the fractional plug flow width is also immediate. A minimum fractional plug flow width of 10% has been taken in Figure 1since the method is of greatest utility where a significant plug exists. All curves show a power-law-like behavior for high flows, sublinear for n < 1. For highly shear thinning
0.0 0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
Q I L s-I Figure 3. Dimensional pressure gradient VB flow rate; experimental and calculated. 0.5% xanthan g u m in water; Herschel-Bulkley par?uneters: T~ = 1.59 Pa; n = 0.54; A = 0.143 Paas". Error bars indicate transducer accuracy. Bands to calculated curves indicate sensitivity of calculation to experimental errors in T~
fluids (small n),the dependence of P on flow is weak, and at low flows P may exceed that for a less shear thinning fluid, although the reverse is usually the case for P well above the starting flow value. 4.2. Experimental Verification: Pressure Drops. Although we shall see below that calculation of pressure gradient can be remarkably insensitive to the geometry, it is not insensitive to the rheological parameters, especially T@ We have made an experimental test of the calculations of pressure drop using a 0.5% aqueous solution of a commercial xanthum gum. The fluid exhibited an apparent yield stress, and ita rheology was found to be well represented by the Herschel-Bulkley model. Rheological parameters as determined from data taken with a Fann 35 viscometer at six different shear rates in the range 5.1-1022 s-l are given in the caption to Figure 3. As a check on the existence of a yield stress in this material, subsequent rheological measurements on similar gum at the Same concentration (0.5%) in aqueous solution have compared traditional viscometry with the vane rheometer technique as described by, e.g., Dzuy and Boger (1985). The yield stress wm determined by measuring the maximum torque developed over a yield surface within the fluid at very low rotation rates. A four-bladed vane (40.7 mm in diameter by 44.0 mm high) at a rotational speed (approximately 1.5 revolutions per hour) of 0.026 rad was used in a container of 100-mm diameter. For a 0.5% solution, the yield stress obtained was comparable to that obtained by fitting to the Herschel-Bulkley formula data obtained in a traditional concentric cylinder rheometer over a shear rate range of 5.1-1022 s-l (13 points). For a weaker solution (0.25%), the vane method gave a yield value about 50% lower than that found from the conventional rheometer and the Herschel-Bulkley model. These data indicate that, a t the concentrations used in the pressure drop test, the Herschel-Bulkley fit to the conventional rheometer data is likely to indicate a reliable value for the yield stress, although for weaker solutions discrepancies with the vane method may occur. The pressure drop measurements were made in an experimental rig designed for studying flow through annuli. The test section consisted of two vertical concentric cylinders 300 cm long, with inner and outer diameters of 4 and 5 cm, respectively, producing a 5-mm uniform annular gap. A differential pressure transducer (Honeywell, Model 41105) calibrated in the range 0-1400 mm of H20, with
Ind. Eng. Chem. Res., Vol. 30, No. 3, 1991 521 1.00 7 1.00.
a
El
0.20 1
hulus
,
1
0.75
0.15
-
v1
E
0.10
\
0
'S 3
* 31.0
32.0
33.0
'34.0
35.0
rlmm
b"--
0.05
0.00
-35
-25
-15
-5
5
15
25
35
rlmm
I
Figure 6. Velocity profdm for very wide annulus and a pipe. Model rheology and volumetric flow a for Figure 4. h u m gradienta 830 (annulus)and 778 Pa-m-' (pipe). Outer wall shear rata 49 (annulus) and 45 s-l (pipe). 3200 Annulur
-0.02 31.0
32.0
33.0
34.0
35.0
r/mm
Figure 4. (a) Velocity profdea for a very narrow annulus and a slot. Model rheology: so = 10 Pa; n = 0.66; A = 0.3 Pa@. Volumetric flow rate 2.0 m8 h-' in both cases. Calculated pressure gradients 21 721 and 21725 Paem-'; outer wall shear rates 1373 and 1411 s-l. (b) Relative velocity differences in profiles of (a).
pipe flow
Figure 7. Pressure gradient vs radius ratio of annulus, for fixed mean flow velocity of 0.25 ms-' in channel of r, = 35 mm. Equivalent slot and pipe calculationsshown for comparison. Model rheology a for Figure 4.
o.oo""""'''"'""'"''""''"'~ 0
5
10
15
20
25
30
35
rlmm Figure 6. Velocity profiiee for a wide annulua and a slot of the same gap size. Model rheology and flow rate as for Figure 4. Calculated pressure gradients 888 (annulus) and 921 Pe.m-' (slot); inner and outer wall shear rates (annulus)260 and 53 8-l; wall shear rate (slot) 75 6-1.
an accuracy of h0.590of its full-scale range, was used. The transducer and connecting tubing were flushed frequently with deionised water to expel any gas bubbles and inhibit ingress of polymer solution. The fluid was pumped through the annulus under computer control, and the pressure drop was measured at several flow rates. The results are given in Figure 3, with error bars indicating estimated transducer accuracy from the manufacturers' specification. The calculated curve was obtained by using the separately measured rheological parameters as input and is shown in the figure for comparison. The sensitivity of the calculation to uncertainty in the measured yield
stress is also shown. Agreement between experiment and calculation is generally within or close to the limits of the experimental uncertainties. 4.3. Velocity Profllee. Examples of calculated velocity profiles are shown in Figures 4-6 for centered annuli of various aspect ratios. In the case of Figure 4a, we show an extreme case of a very narrow annulus, which is intuitively reasonable to approximate by a slot. The calculated profile is only slightly skewed from the profile calculated for the same volumetric flow per unit width in an infinite plane slot (using the analysis outlined in Appendix l),and the pressure drop is for practical purposes identical. Parameters are given in the caption. Figure 4b shows the relative error in the slot approximation for the whole of the velocity profile. For "wide" annuli, however, the slot flow approximation is not acceptable for calculation of velocity profiles or derived parameters such as the wall shear rates required in some applications. Figure 5 shows a case where the full analysis is required for the velocity profde and the slot flow approximation clearly fails. Although the pressure gradients calculated by the two methods are still fairly close, the wall shear rates (and therefore stresses) are not. See the Figure 5 caption for data. It is in this type of intermediate radius ratio that our calculation achieves its full utility. The experimental verification of a calculated velocity profile for various flows, rheologies, etc., is a more involved project than pressure drop measurements and has yet to be done.
522 Ind. Eng. Chem. Res., Vol. 30, No. 3, 1991 z = fractional plug flow width y o / H
plug radius r*/rz (pipe flow)
(slot flow) or fractional
Greek Symbols
.y j
I=
-lo
0
+yo
I
Figure 8. Slot flow notation. Slot width (into paper) is W.
It is interesting to compare an extreme wide annulus case with the pipe flow calculation. The velocity profiles for such a case are shown in Figure 6 with the profile calculated (using the analysis outlined in Appendix 2) assuming a circular pipe but no center body. Obviously wrong within and close to the center body, the approximation is nevertheless satisfactory for most of the true flow region. Outer wall shear stresses agree to within 8%. The calculated pressure drops are a h in fair agreement. These are examined in more detail in Figure 7, which shows calculated pressure drops as a function of annulus radius ratio (for fixed mean flow velocity as the center-body radius is varied). The pressure drops on the slow flow approximation, and for pipe flow, are shown for comparison.
Nomenclature Dimensions are indicated by the standard S.I.unit in brackets; no entry indicates a dimensionless quantity. Numbers in parentheses are equation numbers. A = consistency index in Herschel-Bulkley model (6) [Pa@] A’= consistency index in Robertson-Stiff model (6) [Pa@] B = flow behavior index in Robertson-Stiff model (6) b = 1/B C = shear-rate offset in Robertson-Stiff model (6) [s-l] f = dimensionless shear rate (9) 9 = defined by (11) 9 = defined by (12) H = half-height of infinite plane slot [m] m = l/n n = flow behavior index in Herschel-Bulkley model (6) p = fluid pressure [Pa] P = pressure (or head) gradient -dp/dx [Pa m-’1 q = model-dependent index defined in (9) Q = volumetric flow rate [m3s-l] Q* = dimensionless flow rate QAq/rr 31‘ (annulus) QB* = dimensionless flow rate Q/Wdr (slot) Qp* = dimensionless flow rate Q/&r (pipe) r = radius coordinate [m] rl = inner radius of annulus [m] r2 = outer radius of annulus or pipe radius [m] r+ = yield radius in pipe flow [m] R = dimensionless radial coordinate r / r 2 R1 = annulus aspect ratio r1/r2 Ro = dimensionless zero-stress radius (annular flow) R, = dimensionless outer (+) and inner (-) yield radii (annular flow) R. = defined by (14) R’ = dimensionless radius, r/r* (pipe flow) u = flow velocity (in x direction) [m s-l] W = width of plane slot segment [m] x = axial (flow direction) coordinate [m] y = channel width coordinate (slot) [m] yo = distance of yield surface from center (slot flow) [m] Y = dimensionless coordinate y / y o
i. = (unsigned) shear rate Idu/drl (annulus, pipe) or Idu/dyl (slot) [s-ll = characteristic shear rate T O / V O (Casson model) [e-’] r(HB) = characteristic shear rate (ro/AIm(Herschel-Bulkley model) [s-l] = characteristic shear rate (4(Robertson-Stiff model) [s-ll r = characteristic shear rate for relevant model [s-l] A = dimensionless inverse pressure gradient 2s0/h2(annular flow) vo = high-shear viscosity in Casson model (6) [Pes] T = shear stress [Pa] T~ = yield stress of viscoplastic fluid [Pa]
Appendix 1. Analysis: Infinite Plane Slots Flow in a plane slot is a limiting case of the annulus but is best analyzed by noting the obvious symmetries. See Figure 8 for notation and compare Figure 2. The stress profile is linear across a slot and the yield surface distances &yo from the channel center (.y = 0) are immediate; the fractional plug flow width z in a slot of half-height H is z = yo/H = To/PH. Velocity profiles and volumetric flows are then calculated analogously to (10). Writing Q,* = Q/ W H T ,where Q/ W is the flow per unit width W of slot, one obtains (where Y = y / y o )
(xl”Y(Y
$={
il”U(Y
- 2Y1I2+ 1)dY - l)mdY
[ x 1 ” Y ( P - 1) dY
(Casson)
(Herschel-Bulkley) (16) (Robertson-Stiff)
which simplify to iQ8*z
= 1 - -z1/2 12 5
+ -z32
(m + l)(m + 2)Qa*zm= 2(1 - z)”’+l[z + m 1
-(b 2
+ 2)Q8*zb= 1 1 - & ( b+ 2)zb + ~ 2
1 10
- -23
(Casson)
+ 11 (Herschel-Bulkley)
b z ~ + ~ (Robertson-Stiff) (17)
where we write all left-hand sides to be independent of ro, analogously to (10). The Casson result was given by Kooijman and van Zanten (19721, and the HerschelBulkley result is implied in the paper of Grinchik and Kim (19721, although not given in this form. These formulas give practical solutions to the “direct” problem and these rheologies. The “inverse“ problem in general requires a suitable root-finding method to solve for z. However, in the well-known special case of the Bingham plastic: 3 3 1 ~ Q , * z 1 - -2Z + -z3 2 there is a little-known analytic solution: z = 2(~,*
+ 1)1/* sin
arcsin (Q,*+ 1)-3/2
which is fairly eaey to use, although we have not found it previously in the literature.
Ind. Eng. Chem. Res., Vol. 30, No. 3, 1991 523
Appendix 2. Analysis: Circular Pipes Similarly to the above, stress profiles are linear, the fractional plug radius z is given by 4 = r 4 r 2 = 2sO/Pr2, where r+ is the yield radius, and volumetric flows are given by ( ~ l ’ * R ’ z ( R ’ - 2R’1/2+ 1) dR‘ (Casson)
9 = {j1”Rt2(R’z3 (J1”R‘2(R’b
1)”’ dR’ (Herschel-Bulkley) (20)
- 1)dR’
(Robertson-Stiff)
where Qp* = Q / ~ r and ~ ~R’r = r/r*. These simplify to 16 1 4Qp*z = 1 - -zl/’ + -24 - -z4 7 3 21
1’
m+l
(Casson)
(Herschel-Bulkley)
Now because of the following inequalities
Re2
sR’f Ri
dR > S R ’ R 2 fdR (24) Rl
(where f = f(R;R.,R,) in the above integrands) (23) implies 9 < 0 as required. Note that A = 0 corresponds to infinite flow rate. Hence by lemmas 1 and 2 the bounds given in (13) correctly bracket a zero of 9. Lemma 3: For R+ = 1,3 > 0. Proof: Since f 1 0 for all R and f is not zero everywhere, 3 > 0. Lemma 4: For R+ = R-, and R- < R,, then 9 < 0. Proof: The expression for 3 becomes
and the integrands above reduce to
(b
3)Q,*Zb = 1 1 1 - - ( b + 3)zb + ~ b z * + ~ (Robertson-Stiff) (21) 3 The Casson result here was previously obtained by Kooijman and van Zanten (1972) and by Lih (1975);the Herschel-Bulkley result was obtained by Froishteter and Vinogradov (1980). Numerical root-finding is necessary in general for the “inverse” problem. The well-known Bingham plastic special case
where q is model-dependent as defined as in section 2.3 (but real and positive). Now since by assumption R- < R,, the following inequalities hold:
is a quartic equation that in principle has a closed-form solution for 2. However, unlike the slot flow case (a cubic equation for z ) , the solution involves several levels of nested fractional exponents and is thus unlikely to be more convenient in practical calculation than a numerical method.
Appendix 3. Details of the Numerical Scheme A3.1. Bounds for the Functions 3 and 9. This section shows that the bounds given for the zeroes of 9 and 3 in the numerical problem of section 3.3 are correct. This is shown by lemmas 1 and 2 for the outer iteration (13) and by lemmas 3 and 4 for the inner iteration (15). Lemma 1: For R- = R1, and 9 = 0 , 9 > 0. Proof: For R- = R1, 3(Rl,R+) = -&:f(R;Rl,R+) dR
3 = 0, so R+ = 1 and S(R1,I) = Q* 2 0 with equality when there is no flow ( Q = 0). Lemma 2: For R- = R, and 9 = 0 , 9 < 0. Proof Since 9 = 0, for R- = Re, 9
for which R+ = Re is a solution by the definition of R. (14). Thus A = 0 and
and so 9 < 0. Hence by lemmas 3 and 4 the bounds given in (15) correctly bracket a zero of 9 when R- < Re. A3.2. Pseudoprogram for the Numerical Solution. The simplest algorithm is set up with a “reverse communication” form of Brent’s method. On the first call to the routine “brent” the brackets for the root are provided with the corresponding function values. On each return from this routine a new estimate is provided for the root and the callin? routine is expected to provide the function value corresponding to this estimate. “Reverse communication” is a technique that reduces the complexity of iterative algorithms. With the more usual forward communication, a routine that needs to find the root of a transcendental equation will call an iterative root-findingscheme and will pass to that routine the name of the function that evaluates the equation. With reverse communication, the routine that needs to find the root also evaluates the transcendental equation while the iterative root finding scheme just returns new guesses for the root. To do this, the iterative scheme generally needs to hold some information from previous invocations. In particular,
Znd. Eng. Chem. Res. 1991,30, 524-532
524
it remembers either the initial bracketing interval or some updated smaller bracket for the root. The full algorithm follows, in a "pseudoprogram" format. The subroutine "brent" uses its arguments to bracket the root on the first call and thereafter updates its first argument to provide a new estimate of the root.
START: comment Find R, such that 3(R,,RR,)= 0 comment Root bracket R, E (Rl,l) Re := R1 lOOP call brent (R,,3(Re,R,),l,3(l,l)) until 3(R,,RR,)E 0 comment Find yield surfaces such that 3(R-,R+) = 0, S(R-8,) = 0 comment Outer root bracket is R- E (Rl,RR,-t)
R- := R1 loop comment Inner bracket is R+ E (RJ) R+ := Rloop call brent(R+,S(R-,R+),1,3(R-,l)) until 3(R-,R+)= 0 call brent(R-,O(R-,R+),R. - c,O(RI-e,RR.)) until $'(R-,R+) = 0
END. Literature Cited Anshus, B. E. Bingham plastic flow in annuli. Ind. Eng. Chem. Fundam. 1974,13,2,162-4. Bird, B. R.; Dai, G. C.; Yarusso, B. J. The rheology and flow of viscoplastic materials. Rev. Chem. Eng. 1983,1 (l),1-70. Brent, R. Algorithms for Minimisation without Derivatives; Prentice-Hall: Englewood Cliffs, NJ, 1973.
Casson, N. A flow equation for pigment oil suspensions of the printing ink type. In Rheology of Disperse Systems; Mill, C. C., Ed.; Pergamon: Oxford, UK, 1969. Dzuy, N. Q.; Boger, D. V. Yield stress measurements with the vane method. J. Rheol. 1985,29 (31,335-47. Frederickson, A. G.; Bird, R. B. Non-Newtonian flow in annuli. Ind. Eng. Chem. 1958,50(31,347-352. Froishteter, G. B.; Vinogadov, G. V. The laminar flow of plastic disperse systems in circular tubes. Rheol. Acta 1980,19,239-50. Grinchik, I. P.; Kim, A. Kh. Axial flow of a non-linear viscoplastic fluid through cylindrical pipes. J. Eng. Phys. 1974,23,1039-41 (published in Russian 1972). Hanks, R. W. The axial laminar flow of yield-pseudoplasticfluids in a concentric annulus. Ind. Eng. Chem. Process Des. Dev. 1979, 18 (3), 488-93. Herschel, W. H.; Bulkley, R. Konsistenzmessungen von GummiBenzollosungen. Kolloid 2. 1926,39,291-300. Jaisinghani, R. Annular flow of a Casson fluid. M.S. thesis, University of Wisconsin-Milwaukee, 1974. Kooijman, J. M.; van Zanten, D. C. The flow of a Cassonian fluid through parallel-plate channels and through cylindrical tubes. Chem. Eng. J. 1972,4,185-8. Lih, M. M A . Transport Phenomena in Medicine and Biology; Wiley: New York, 1975;Chapter 9. Press, W. H.; Flannery, B. P.; Teukolsky, S. A.; Vettering, W. T. Root-finding and non-linear sets of equations. In Numerical Recipes: the art of scientific computing; Cambridge University Press: Cambridge, UK, 1986,Chapter 5,see especially pp 251-4. Robertson, R. E.; Stiff, H. A. An improved mathematical model for relating shear stress to shear rate in drilling fluids and cement slurries. SOC.Pet. Eng. J. 1976,Feb, 31-36. Shah, V. L. Blood flow. Adu. Tramp. Processes 1980,1, 1-57. Shul'man, Z. P. Calculation of a laminar axial flow of a non-linear viscoplastic medium in an annular channel. J. Eng. Phys. 1973, 19, 1283-9 (published in Russian 1970). Received for review February 21, 1990 Revised manuscript receioed August 24, 1990 Accepted September 12, 1990
Applications of a Generalized Equation of State for Associating Mixtures S. Jayaraman Suresh and J. Richard Elliott, Jr.* Chemical Engineering Department, University of Akron, Akron, Ohio 44325-3906
A generalized equation of state that explicitly recognizes self-association and cross-association has been evaluated for systems with two associating species a t high pressures as well as low pressures. Forty-five systems were evaluated for components ranging from water to propanol. The Soave equation of state and the UNIFAC model were used for comparison. Percent average absolute deviations in bubble-point pressure were 3.54 for the Soave equation, 2.53 for UNIFAC, and 1.62 for the association model. Generalized correlations of the pure-component parameters and a single binary interaction parameter were used for all models. Some limitations are discussed, but the results show that the reaction equilibrium approach is a clear improvement over the Soave equation of state and competitive with UNIFAC. Introduction For many years, one of the most challenging problems in thermodynamics hae been prediction of phase equilibria in hydrogen-bonding systems. These systems are commonly encountered in the chemical industry, and the ability to make quantitative predictions would facilitate optimization of many commercial processes. Previous successful treatments of these mixtures have used activity coefficient methods. The Van Laar equation, the Wilson (1964) equation, and UNIFAC (Fredenslund et al., 1975) are examples of significant developments by this approach. The Van Laar and the Wilson equations are usually considered to be correlative models whereas
the UNIFAC model is considered to be a predictive model. Nevertheless, we have implemented UNIFAC to evaluate the activity coefficients in the liquid phase only and used the Soave equation (1972) to evaluate the fugacity coefficients in the vapor phase. The Soave equation for the vapor phase requires an adjustable binary interaction coefficient. Hence, our implementation of UNIFAC in conjunction with the Soave equation becomes more a correlation with one adjustable parameter. This treatment is typical of the way these equations are used in process simulators like ASPEN. Many sets of UNIFAC parameters have been compiled; therefore the specific set of UNIFAC parameters was taken from Reid et al. (1987) for
088&5885/ 91/2630-0624$02.50/ 0 0 1991 American Chemical Society
