u-“ . y ; _ . urn-mi. W4!” >. ‘ [w ‘¢|(.“¢ s “I I . ”’4‘,‘ A . , u r :n ; .4! .- ,_ -Z‘E fff“ 9i. , r ‘ '3?! o . v 4 ‘3‘ , _ '5‘ «1‘0 ‘ “8'3", .2. I31 \y u . . old" -. u; 4‘ ‘ 'ml ‘ $13, ’ v)". '2' 5).. ‘1"! 1"; r ‘1. . In. a a, , if)“ u r a . .r (: l' v‘ ":33- *' .. ., 1‘ J J. I I .0 .‘a’zW/*“' , I u u: 5 S A.) - I :r 4 n . i .- al.‘n‘.f‘ zirsw .. .. “first-11. A %( V A 7A ‘- ":2! .‘ . i \- 'rl ‘ ‘ . . . ~r " ' ,'.’.-.\,; v n '3': . o -, 1132:," - s. v r ‘ m 4... u w h} u: ,» |'|- i‘rzf r J” ""1...- ”(”3 ' ' “kafiwé ‘ z “ "W, "55)". an , v ((er ’éw ' I I - . wt -'r r «1).. . Q' 3:4,; _. , l' v .57 . war'ne {w uh. . Jim .fi - 2‘5" 0 .fifiwzu '~ Date MICHIGAN TATE UNIVERSITY LIBRARY IUIIIHIWIIIW "H"H""1H""NIHIIIIUIIWIIIHII ‘1 3 1293 10599 1933 113MB? , mcmuafitate University I“ fl This is to certify that the thesis entitled OPTIMUM ECONOMIC TUBE DIAMETER FOR PUMPING HERSCHEL-BULKLEY FL UI DS presented by Edgardo Jose Garcia Caes has been accepted towards fulfillment of the requirements for M.S. degree in Food Science 7 Major professor May 3, 1985 James F. Steffe 0-7639 MS U is an Affirmative Action/Equal Opportunity Institution MSU RETURNING MATERIALS: Place in book drop to LIBRARIES remove this checkout from —3—_ your record. FINES will be charged if book is returned after the date stamped below. sap/36 296A .1 z 7 I '1 OPTIMUM ECONOMIC TUBE DIAMETER FOR PUMPING HERSCHEL-BULKLEY FLUIDS By Edgardo Jose Garcia Caes A THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of MASTER OF SCIENCE Department of Food Science and Human Nutrition 1985 IN MEMORY OF Willie and Jaime‘ DEDICATED TO My Family and Makito ACKNOWLEDGMENTS I wish to express my deepest appreciation to Dr. James F. Steffe, my major professor, for his guidance, suggestion, and continuous support throughout the course of this study. Sincere appreciation is extended to the guidance committee members: Dr. Ajit Srivastava, Dr. Kris Berglund, and Dr. Eric Grulke. A very special thanks to to my parents, Edgardo and Nancy; to Nancy Jeannette , Maria Elena, Mario, Patricia, Patricia Lorena, and Jaime Jose, for their support and for being a wonderful family; to Makito for her patience, understanding, and words of encourage- ment. A most sincere gratitude goes to Mr. Philip Lehner for his kindness and assistance during my studies in the United States. Thanks is also given to Nancy Heath for the typing of this thesis. iv TABLE OF CONTENTS Page LIST OF TABLES . . . . . . . . . . . . . . . ‘. vii LIST OF FIGURES . . . . . . . . . . . . . . . ix LIST OF APPENDICES . . . . . . . . . . . . . . X LIST OF SYMBOLS . . . . . . . . . . . . . . . xi Chapters 1. INTRODUCTION . . . . . . . . . . . . . . 1 1.1 General Remarks 1 1.2 Objectives 2 2. LITERATURE REVIEW 3 2.1 Herschel-Bulkley (H—B) Model 3 2.2 Optimization . . . 6 2.3 Power Requirements . 7 2.4 Economic Considerations . . 14 2.5 Optimum Economic Pipe Diameter 18 3. THEORETICAL DEVELOPMENT . . . . . . . . . . 20 3.1 Flow Behavior of Herschel— —Bulkley Fluids . . . 20 3.1.1 Laminar Flow . . . . . . 20 3.1.2 Laminar- Turbulent Transition . . . . . 26 3.1.3 Transitional and Turbulent Flow . . . '. 27 3.2 Total Annual Cost of a Pumping System . . . . 29 3.3 Optimum Economic Tube Diameter . . . . . . 35 3.4 Limitation of Design Method . . . . . . . 4O 4. RESULT AND DISCUSSION . . . . . . . . . . . 44 4.1 Model Verification . . 44 4.2 Cost Parameters and Other Pertinent Data for a Pumping System . . . 47 4.3 Optimum Diameter for a System Handling Tomato Ketchup . . . . . . . . . . . . . . 51 Chapter Page 4.4 Sensitivity Analysis . . . . 57 4.5 Optimum Diameter Using Apparent Viscosities . . 62 4.6 Optimum Diameter for a System Handling a Herschel-Bulkley Fluid in Turbulent Flow . . . 66 5. SUMMARY AND CONCLUSIONS . . . . . . . . . . . 71 APPENDICES . . . . . . . . . . . . . . . . . 73 REFERENCES . . . . . . . . . . . . . . . . . 122 vi Table 10. LIST OF TABLES Fluid properties and other pertinent data for the optimum diameter example problem given by Skelland 1967) . . . . . . . . . . Fluid properties and other pertinent data for the opti— mum diameter example problem given by Darby and Melson (1982 . . . . . . . . . . Variation of the installed cost at a tube system per meter length of tube with tube diameter Cost parameters for the tube system presented in, Table 3 . . . . . . . . Variation of the installed pump station cost with power requirements . . . . . . . Cost parameters for the pump station presented in Table 5 . . . Electrical energy cost, hours of operations per year, combined pump and motor efficiency, summation of the fittings resistance coefficients, tube length, pres— sure and elevation change, and mass flow rate used to estimate the costs, and optimum diameter for the pumping system presented in Table 3 and 5 . Rheological properties (Higgs and Norrington, 1971) and density (Lopez, 1981) for tomato ketchup at 25°C Optimum economic tube diameters, pumping system costs, and work and power requirements, at optimum for a system (Tables 4, 6, and 7) transporting tomato ketchup with properties given in Table 8. Flow condition of Do opt = 0.06907 m for the pumping system handling 4.0 kg/s of tomato ketchup with proper- ties given in Table 8 vii Page 45 46 48 50 52 54 54 55 55 Table Page 11. Percent change of D with :10% change of imput opt variable for the example presented in Section 4.3 . . 59 12. Cost constants of Equation (3.37) for the pump station presented in Table 5 based on Figure 10 . . . . . 62 13. Optimum economic tube diameter, pumping system costs, and work and power requirements at optimum estimated assuming Newtonian flow behavior and an apparent viscosity of 4.723 Pa 5 . . . . . . . . . . 64 14. Pumping system costs, and work and power requirements estimated using the H-8 model for D = 0.1138 m . . . 64 15. Optimum economic tube, diameter, pumping system costs, and work and power requirements at optimum estimated assuming Newtonian flow behavior and an apparent vis- cosity of 1.715 Pa 5 . . . . . . . . . . . 65 16. Pumping system costs, and work and power requirements estimated using the H-8 model for D = 0.09271 m . . 67 17. Rheological properties and density for a hypothetical Herschel-Bulkley fluid. . . . 67 18. Optimum economic tube diameter pumping system costs and work and power requirements, power at optimum for a system (Table 4, 6, and 7) transporting a H-B fluid with properties given in Table 17 . . . . . 68 19. Flow condition at 00 Mp = 0.04065 m for the pumping system handling 4.0 kg/s of a H- B fluid with proper- ties given in Table 17 . . . . 68 viii LIST OF FIGURES Figure 10. 11. 12. 13. Shear stress - shear rate relationship for time- independnet non-Newtonian and Newtonian fluids Optimum economic pipe diameter for minimum total cost at a fixed mass flow rate Velocity profile for Herchel-Bulkley fluid in a tube Calculation scheme to estimate the friction factOr . Calculation scheme to estimate the annual costs of a pumping system . Calculation scheme to estimate the optimum diameter Variation of the installed cost of the tube system with tube inside diameter Variation of the installed cost of the pump station with poer requirements Variation of tube system cost, pump station cost, operating cost, and total cost with tube inside diameter for a system transporting tomato ketchup with properties given in Table 8 . . . . . Variation of the installed cost of the pump station with poer requirements . . . . . . Variation of tube system cost, pump station cost, operating cost, and total cost with tube inside diameter for a system transporting a H-B fluid with properties given in Table 17 . . Friction factor versus tube inside diameter for the fluid properties given in Table 8 and a mass flow rate of 4.0 kg/s . . . . . . . . Friction factor versus tube inside diameter for the fluid properties given in Table 17 and a mass flow rate of 4.0 kg/s . . . . . . . . ix Page 21 30 36 41 49 53 56 61 69 81 82 LIST OF APPENDICES Appendix A. Alternative Equations for f and df/dD for H- B Fluids in Turbulent Flow . 8. Verification of df/dD Equation for Laminar Flow C. Approximation of df/dD for Turbulent Flow 0. Listing of Computer Program Page 74 77 83 86 bl LIST OF SYMBOLS annual fixed cost of the tube system expressed as a fraction of the initial installed cost of the tube system, l/yr, Equation (2.16) annual fixed cost of the pump station expressed as a fraction of the initial installed cost of the pump station, 1/yr, Equation (2.16) empirical wall effect parameter in mixing length theory, Equation (3.28) annual maintenance cost of the tube system expressed as a fraction of the initial installed cost of the tube system, 1/yr annual maintenance cost of the pump station expressed as a fraction of the initial installed cost of the pump station, 1/yr total installed cost of the tube system including the cost of fittings, valves, installation, etc., $/m, Equation (2.17) unknown cost of equipment of size 02, Equation (2.18) or (2.19) known cost of equipment of size 02 purchase cost of a new pipe, $/m, Equation (2.14) empirical constant for the pump station cost, S/WSI total direct cost of equipment of size 02 cost of electrical energy, S/N hr empirical constant for the pump station cost, $ total indirect cost of equipment of size 02 total annual operating cost per unit length of tube, $/yr m, Equation (2.20) xi C . p1 Cps Cpu est empirical constant for the tube system cost, $/m1+S total annual cost of installed tube system per unit length of tube, $/yr m, Equation (2.15) or (3.36) total cost of installed pump station, $, Equation (3.37) total annual cost of installed pump station per unit length of tube, $/yr m, Equation (3.38) total annual cost of a pumping system per unit length of tube $/yr m, Equations (3.7),(3.39), or (3.40) minimum total cost of a pumping system per unit length of pipe. $/yr m, Equation (3.39) or (3.40) tube/pipe inside diameter, m best estimate of D0 to start numerical search, m pt value of D at CT min combined fractional efficiency of pump and motor energy loss due to friction, J/kg, Equation (2.10), (2.11), and (2.12) ratio the total cost for fittings and installation of pipe and fittings to purchase cost of new pipe Fanning friction factor, Equation (3.15) laminar-turbulent transition valve of f best estimate of f to start numerical search acceleration due to gravity (9.8 m/sz) hours of operation per year generalized Hedstrom number, Equation (3.19) interest rate (fraction) consistency coefficient, Pa sn dimensionless fittings resistance coefficient Prandtl's universal mixing length constant = 0.36 tube/pipe length, m xii equivalent length of pipe for fittings frictional loss modified Prandtl's mixing length mass flow rate, kg/s life-time of equipment, yr flow behavior index, dimensionless power, Watts, Equation (2.13) constant for purchase cost of pipe dependent on the pipe material, dimensionless pressure at point 1, Pa pressure at point 2, Pa size of equipment with unknown cost size of equipment with known cost cost capacity factor turbulance parameter, Equation (3.27) radial coordinate laminar-turbulent transition valve of R generalized Reynolds number, Equation (3.18) laminar-turbulent transition valve of Re, Equation (3.21) tube/pipe radius, m exponent of tube system cost equation, dimensionless exponent of pump station cost equation, dimensionless dimensionless velocity, v/v dimensionless plug velocity, vo/v axial velocity component, m/s plug velocity, m/s mass average velocity, m/s xiii value of v at point 1, m/s value of 9 at point 2, m/s. work per unit mass, J/kg, Equations (2.9) (3.33) or (3.35) purchase cost of one.inch diameter pipe per unit meter of pipe length, $/m inp a small positive number elevation at point 1, m elevation at point 2, m axial coordinate Greek Letters kinetic energy correction factor at point 1 kinetic energy correction factor at point 2 -du/dg value of F at the pipe wall (r r 1 W1 rate of shear (~dv/dr),s' value of i at the pipe wall pressure change between points 1 and 2 (p2 - p1), Pa pressure drop due to friction, Pa elevatiOn change between point 1 and 2, (22 - 21), m dimensionless rate of shear,r/rw value of g calculated with A. and g < g. < 1 for j-1,2,3.. 3 °"J—' pastic viscoisty, Pa 5 laminar stability parameter, Equation (3.20) maximum value of K at g = g , 404 dimensionless mixing length, St/rw xiv value of A calculate for go 5_ gj.: 1 for j = 1, 2, 3 . . . Newtonian viscosity, Pa 5 apparent viscosity, Pa 5, Equation (2.6) dimensionless radial coordinate, r/rw value of g where K=1E value of g for 50-: Ej-i 1 j = 1, 2, 3 . . ., use in to evaluate the integral in Equation (3.31) dimensionless unsheared plug radius, IOI’TW laminar—turbulent transition valve of go Pi(3.1415) fluid density, kg/m3 parameter in the df/dD Equation for laminar flow, Equation (3.45) yield stress, Pa shear stress, Pa time average momentum flux, Equation (3.23) molecular flux, Equation (2.1) turbulent flux or Reynolds stress, Equation (3.24) value of Trz at r = rw, Equation (3.4) parameter in mixing length, Equation (3.26) laminar flow function, Equation (3.14) laminar-turbulent transition value of w XV ABSTRACT OPTIMUM ECONOMIC TUBE DIAMETER FOR PUMPING HERSCHEL-BULKLEY FLUIDS By Edgardo Jose Garcia Caes The optimum tube diameter, for which the total cost of a pumping system is a minimum, has been derived for the case of Herschel-Bulkley fluids in both laminar and turbulent flow. The method accounts for the tube system cost as a function of diameter, as well as the pump station and operating costs as a function of the power requirements. The optimum diameter can be estimated given the rheological properties, density of the fluid, mass flow rate, and economic parameters. The elevation and pressure difference in the system are irrelevant when a linear relationship is used for the pump station cost. The friction loss in fittings can be ignored when the tube length is much greater than the tube diameter. The pump station cost has less influence than the operating cost in determining the optimum diameter. The use of apparent viscosity and Newtonian flow behavior for non-Newtonian fluids may lead to severe errors in pipe sizing. 1. INTRODUCTION 1.1 General Remarks A problem associated with the design of fluid handling systems is the selection of tube or pipe size. The installed cost of a process piping system can vary between 7% and 60% of the total fixed investment (Wright, 1950). It is, therefore, important to choose the tube size that would result in the greatest economy while maintaining the designated operating conditions and performance requirements. Three criteria often control the selection of tube size; the pressure drop available, velocity allowable, and least annual cost. The first criterion is usually used when a given pressure drop must be absorbed by the tube. Limits in velocity may be encountered in the handling of slurries in which a minimum velocity must be maintained to keep the particles in suspension. Conversely, quality degradation of the product may restrict high velocities. The least annual cost applies when a given amount of fluid is to be pumped through a tube system. It is based on an economic balance of the capital and operating cost to give a tube size that will result in the least annual charge (Nolte, 1978; Kent, 1978). In this study, techniques to estimate the optimum tube diameter based on the least annual cost are developed for tube systems trans- porting non-Newtonian fluids. The Herschel-Bulkley model was selected due to its generality and wide application in fluid foods, as well as other fluid materials (Holdsworth, 1971; Higgs and Norrington, 1971; Steffe et al., 1983; Boger and Tiu, 1974). Non-Newtonian characteristics must be considered in the design of pumping systems when handling fluids of this type (Cheng, 1975; Johnson, 1982). Failure to do so may lead to under or over sizing, resulting in a system inefficienttx>operate or more costly to erect as suggested by Steffe (1983) and Nolte (1978). 1.2 Objectives The specific objectives of this study are as follows. Objective 1: Develop an equation to predict the total annual cost of a pumping system as a function of the tube diameter. Objective 2: Develop an equation to estimate the optimum economic tube diameter for pumping systems handling Herschel-Bulkley fluids. Objective 3: Demonstrate the design errors caused by using apparent viscosity and Newtonian flow behavior to design pumping systems handling non-Newtonian fluids. 2. LITERATURE REVIEW 2.1 Herschel-Bulkley (H-B) Model The flow behavior of many fluid foods and other industrially important fluids may be described by the H-8 model which can be written as (Herschel and Bulkley, 1926). .. '0 Trz — To + K y (2.1) where T = shear stress, Pa rz To = yield stress, Pa K = consistency coefficient, Pa 5n n = flow behavior index, dimensionless i = rate of shear (-dv/dr), 3’1 This model simplifies to other well-known models. The power law or Ostwald-de Naele model is written as _ °n Trz - K y (2.2) where A power law fluid is called pseudoplastic when 0 < n < 1 and dilatant when n > 1. Equation (2.1) reduces to the Bingham plastic model when n = 1 and n = K as O Trz = TO + n y (2.3) where n = plastic viscosity, Pa 5 Newtonian fluids are described by Equation (2.1) when To = O, n = 1, and p = K as where p = Newtonian viscosity, Pa 5 The shear stress—shear rate relationships for the above models are shown graphically in Figure 1. It is common practice to use an apparent viscosity (pa) and assume Newtonian fluid behavior to estimate the frictional pressure losses for the flow of non-Newtonian fluids in tubes. Apparent vis- cosity is defined as t = u i (2.5) Since equation (2.5) is used to describe H-B fluids, pa may be written in terms of the H-8 parameters using Equation (2.1) and (2.5) as _ --1 ~n-1 “a — To y + K y (2.6) From Equation (2.6), it is evident that pa is defined at a particular rate of shear. Therefore, the use of an apparent viscosity may lead to over- or under-estimation of the pressure losses and power SHEAR STRESS Figure l. (l) (2) (3) (S) (4) SHEAR RATE Shear stress - shear rate relationship for time-independent non-Newtonian and Newtonian fluids: (1) Herschel-Bulkley model, (2) Bingham plastic model, (3) Pseudoplastic model, (4) Dilatant model, (5) Newtonian model. requirements depending the rate of shear at which the apparent vis- cosity is measured. This, in turn, may lead to improper sizing of pipe, pump, and motor (Steffe, 1983). 2.2 Optimization The selection of a value for a given design variable to mini- mize the total cost of a project is possible whenever a change in this variable causes some costs to increase while other costs decreases (Skelland, 1967). For a pumping system, the total cost can be divided into three components: the tube system cost, the pump station cost, and the operating cost (Darby and Melson, 1982). The tube system cost primarily consist of the installed cost of tube, fittings, and values. This increases with increasing tube diameter (Skelland, 1967; Darby and Melson, 1982; Jelen, 1970). The pump station cost mainly consist of the installed cost of pump and motor while the operating cost parimarily consists of the cost of electrical power required to pump the fluid through the system. Both of these costs are directly pro- portional to the power requirements which decrease with increasing tube diameter since the pressure drop due to friction decreases with increasing tube diameter. Consequently, the pump station cost and operating cost decrease with increasing tube diameter (Skelland, 1967; Darby and Melson, 1982; Downs and Tait, 1953). This is shown graphi- cally in Figure 2. Clearly, the optimum value for the diameter can be obtained when the sum of these costs is at a minimum. Mathematically, the total cost C can be expressed as func- T tion at the tube diameter (D) with the following algebraic equation (Skelland, 1967). Total cost, CT CTmin --- ——----------- Pipe system cost, Cpi COST ($/(yr)(unit length of pipe)) ~“ Pump station cost, Cpu Operating cost, C I I I I 0P I Popt PIPE DIAMETER Figure 2. Optimum economic pipe diameter for minimum total cost at a fixed mass flow rate. CT(D) = C 1.(D) + C (D) (2.7) p um) + cO P P where CT = total annual cost of pumping system per unit length of tube, $/yr m i total annual cost of installed tube system per p unit length of tube, $/yr m C u total annual cost of installed pump station per unit p length of tube, $/yr m COp = gotal annual operating cost per unit length of tube, /yr m The analytical method for optimization of a function of a single variable involves differentiating with respect to the variable and equating the result to zero. So the result, for D, in the total cost equation is .9. d0 c =-d—c T dD (”HEP (ow-(15c (D)=O (2.8) pi Solving Equation (2.8) for 0 gives the optimum diameter for which the total cost is at minimum (Skelland, 1967; Jelen, 1970; Reklaitis et al., 1983). 2.3 Power Requirements The work per unit mass required to pump an incompressible isothermal fluid through a tube system from point 1 to point 2 under steady state conditions is given by the mechanical energy balance equation (Heldman and Singh, 1981) written as I: Ii p2 ' p1 w - Ef + a2 - a1 +--——E;——-+ 9(22-21) (2.9) where w = work per unit mass, J/kg Ef = energy loss due to friction, J/kg p = pressue, Pa p = fluid density, kg/m3 9 = acceleration due to gravity (9.8 m/sz) z = elevation, m v = mass average velocity, m/s a = kinetic energy correction factor 1,2 = subscripts referring to points 1 and 2, respectively Osorio and Steffe (1984) developed an equation for the kinetic energy correction factor (a) for Hershel-Buikley fluids in laminar flow. a is equal natwo for turbulent flow. For the purpose of tube/ pipe selection, however, the change in kinetic energy can be assumed to be zero SlflCG'UEttUbe has a constant diameter (171 = 92) and point one and two have been located far enough from any entrance,bend, or fitting to have the same velocity profile (dl = dz) (Skelland, 1967). The energy loss due to friction in a straight pipe can be written in terms of the Fanning equation as (Govier and Azis, 1972) 10 E1. = .;ngi. (2.10) where f = fanning friction factor L = tube/pipe length, m D = tube/pipe inside diameter, m The pressure drop due to friction depends on the flow char- acteristic, as well as the fluid properties. At slow flow, the fluid velocity is parallel to the tube axis and the pattern is smooth. This condition is known as laminar or streamline flow. As the veloc— ity of flow increases, there is a point where the flow becomes unstable, eddies develop, and cause the fluid to swirl in all directions to the line of flow. The flow is then turbulent. The region from the end of laminar to fully turbulent flow is known as transitional region. The theoretical relationship between pressure drop due to friction and flow rate for a H-B fluid in laminar flow can be optained by integrating Equation (2-1) as shown by Cheng (1970), Charm (1978), Skelland (1967), and Govier and Azis (1972). This relationship can be rewritten in term of the friction factor and generalized Reynolds number (Hands, 1978; Heywood and Cheng, 1982) and will be outlined later in this study. The transitional flow of non-Newtonian fluids has been subject to research for many years. Various criteria of transition has been developed based on the end of the laminar flow regime (Metzner and Reed, 1957; Ryan and Johnson, 1959; and Mishra and Tripathi, 1973). Le Fur and Martin (1967) applied the Ryan and Johnson criterion for ll Bingham and power law fluids. This criterion was also used by Hanks and Christiansen (1962) for nonisothermal flow of pseudoplastic fluids and by Cheng (1970) for H-B fluids. Hanks (1963) developed a more general stability criterion and applied itto the transitional flow of Bingham plastic fluids (Hanks, 1963). More recently Hanks and Ricks (1974) presented the transition flow behavior of H-B fluids based on his theory of laminar flow stability (Hanks, 1969). Numerous equations have been developed to calculate the fric— tion factor of power law (Dodge and Metzner, 1959; Shaver and Merrill, 1959; Kemblowski and Kolodziejski, 1973; Tomita, 1959; Szilas et al., 1981; Clapp, 1961; Hanks and Ricks, 1975), Bingham plastic (Tomita, 1959; Thomas, 1962; Hanks and Dadia, 1971; Darby and Nelson, 1981) and H-B (Torrance, 1963; Hanks, 1978) fluids in turbulent flow. Good reviews of these equations are found in articles by Heywood and Cheng (1982), Cheng (1975), Kenshington (1974), Govier and Azis (1972), and Skelland (1967). Unlike Newtonian flow, the friction factor prediction for non-Newtonian fluids varies greatly, depending on the equation used. This deviation increases with decreasing flow behavior index, but is not very sensitive to the yield stress, up to a Hedstrom number of 104. This is the motivation for using equations based on the power law model to predict friction factor for H-B fluids (Heywood and Cheng, 1982). However, this may lead to over estimation of the friction factor (Cheng, 1970). For all the methods developed for transitional and turbulent flow, the work of Hanks and Ricks (1974) and Hanks (1978) are the most comprehensive in describing the flow 12 behavior of H-B fluids in laminar, transition, and turbulent flow. This work will be presented later in this study. So far only the friction loss in a straight tube has been considered. However, to determine the total pressure drop in a tube system, one must add the friction loss arrising from any fittings, valves, and any other devices in the line. The total energy loss due to friction then can be written in terms of Equation (2.8) and the summation of the energy loss in fitting and other devices (Steffe et al., 1984) as (2.11) m -h I N .4, C31/"+1] (3 9a) (1 - £0) for g>g u - F” < “ > (1 - a )1/" + 1 (3-9b) o (1 _ £0)1/n n + 1 o for g $.50 In terms of the defined dimensionless variable, the expres- sion for the flow rate is given by l E 2 J 0 Euo dg + 2 gudg = 1 (3.10) 0 E 24 Substituting Equation (3.8) into Equation (3.10) yields, since u0 is constant for 0.: g_: :0 N r a c (E',EO) da'1-2 a c (5'.£0) dg' d5 = 1 a a a (3.11) Integrating the double integral in Equation (3.11) by parts, using Leibnitz' rule (Hanks and Ricks, 1974) gives 1 P :2 c (E, to) dz = 1 (3.12) g0 Substituting Equations (3.7) into (3.12) and integrating results in (1- a) 1‘3: 3 n (3.13) w (1':*§fii where (1-.;)2 2, (1-1) 2:2 “ w=(1+3mH1-€dhn[TTrfih'i I1+an°*'1$"EI (3.14) By combining Equation (3.6) and the definition of the Fanning friction factor as 25 f= W = ° (3.15) Equation (3.13) can be written in terms of the friction factor as f = if? (3.16) where w is given by Equation (3.14) and Re, the generalized Reynolds number, is (by definition) n D Re = 8 (Ti-7371')” pl (3.17) If one eliminates 9 using equation (3.6), and the definitions of f and Re, Equations (3.13) may be rearranged (Hanks, 1978) to give 2-n )2 (51)”— (3.18) 0 where w is given by Equation (3.14) and 2 T He = 14(3)? (3.19) Equation (3.19) is a generalization of the Hedstrom number. Equation (3.18) defines go as an implicit function of Re and He for He > 0. go = 0 when He = 0, i.e., To = 0. 26 3.1.2 Laminar-Turbulent Transition Laminar instability starts when the ratio K, the rate of change of angular momentum of a deforming fluid element to its rate of loss of frictional drag momentum, exceeds a critical valve E (Hanks, 1969). For rectilinear pipe flow, the stability parameter can be written as Re pr u; = 9L Ji_ 2 = w K dr (V I 16 (3.20) ZAPf where PW, ;, u, and u are given by equations (3.6), (3.7), (3.9a), and (3.14), respectively. K is a function of the radial position a having the value of zero at g = 1 and g = g0, and a maximum value at some point in the field (5 = E , K = E) where maximum instability occurs. The transitional critical Reynolds number (Rec) is obtained from Equation (3.20) when one sets a = E and E = 404 (Hanks and Ricks, 1974). This valve will give Rec = 2100 for Newtonian pipe flow. The radial position of maximum instability E is found by setting dK/dg = 0. For H-B fluids, the critical Reynalds number is then given by the following expression (Hanks and Ricks, 1978) 2+n 6464 n pz/“‘1 (2 + n) 1+n Rec = 2 1 + 2/n (3°21) (1 + 3n) (1 - 50C) where we is given by Equation (3.14) with :0 = 50c Equation (3.18) is also valid at Re = Rec. Now, by eliminating ReC with Equations (3-18) and (3.21), the relationship between He and goc can be obtained as 27 2 + n 3232 (2 + ”)1 + n EOCZ/n-l He = E/IH'I (3.22) n (1 - 50c) which defines Soc as in implicit function of He (Equation (3.19)) and n. The critical friction factor, fc, can be estimated from Equations (3.16) with w = we and Re = Rec. 3.1.3 Transitional and Turbulent Flow For transitional and turbulent flow, the time average momen- tum flux can be expressed as (Hanks, 1968) -_L-T Trz - Trz + Trz (3.23) where 1:2 is the molecular flux, given by Equation (2.1), and 1:2 is the turbulent flux (or Reynolds stress). This latter flux is given by Hanks and Dadia (1971), Hanks and Ricks (1975), Hanks (1978), as -T _ .2 th - pty (3.24) where t is a modified Prandtl's mixing length (Hanks, 1968) and is given in terms of the dimensionless variable I = t/rw as (Hanks and Ricks, 1975; Hanks, 1978) A = k (1 - F.){1- exp [-¢(1- €)]} (3-25) where 28 R - RC 0 = (3.26) «GE's k = Prandtl's universal mixing length constant = 0.36 and _§;fl l/n R = ( 1 + 3n) Re (i:) 2 (3°27) 16 R is a working parameter and reduces to R = Re/f'for Newtonain fluids (Hanks, 1968). The parameter RC is estimated from Equation (3.27) with Re = ReC and f = fc. The parameter 8 is given by the following empirical relation- ship for the H-8 model (Hanks, 1978). B =.%% 1 + 0.00352 He 2 (3.28) (1 + 0.000504 He) Substituting equations (2.1) and (3.24) into (3.23), Equation (3.23) can be rewritten in dimensionless form as -n n -2 2 2 _ Kv Fw n pv 1 PW 2 g - g + —————-— g + —-—-——- C (3.29) T F Tw where Fw is given by Equation (3.6) since g = 1 and A = 0 when a = 1. Using Equation (3.6), (3.15), (3.17), and (3.27), Equation (3.29) can be written as 29 2 + (1- to) 1;" + 38—0 - a )2/n 122:2 (3.30) 5=6 o 0 Combining equations (3.6),(3.12), (3.15),(3.17), and (3.27), Equa- tions (3.12) can be written in equivalent form 2-n 1 2‘" T n nR2 2 (1 - 60) (j—gfgfii fig' 5 C (EIdE € 1 (3.31) o where t(g) is given implicity by Equation (3.30). Finally, from the definitions of f, Re, R and He, it can be shown that R = 2 (3.32) The methodology to estimate the friction factor is outlined in Figure 4. A computer program (Appendix 0) written in FORTRAN 77 as developed to accomplish these calculations. 3.2 Total Annual Cost of a Pumpingggystem Assuming negligible kinetic energy change and substituting Equation (2.11), the work per unit mass, Equation (2.9), can be written as w =-3£§—L + Kf-7T + 7? + gAz (3.33) Input Variables K,11,to, p M 0, Des or DO Fluid Properties, Mass Flor Rate Pipe/TUbe Inside Diameter t pt 1. 9 Calculate 9 from Equation (3.34) 2. Re Calculate Re from Equation (3.17) 3. He Calculate He from Equation (3.19) 4. 50c Calculate 50c from Equation (3.22) through 1terat10n O_: 50c < 1 5. wC Calculate wC from Equation (3.14) with 50 = 50C ReC Calculate ReC from Equation (3.21) fc Calculate fc from Equation (3.16) with If Re < ReC then laminar flow. 8.-a g0 9.—a w 10.-a f Alternative for laminar flow 8.-b E C’o 9.-b w IO.-b f wc=wC and Re = ReC If transition or turbulent go to Step 11. Calculate go from Equation (3.18) through interation. (£0 = 0 if r0 = 0 (He = 0)). 60C: to < 1 Calculate w from Equation (3.14) Calculate f from Equation (3.16) Calculate go from Equation (3.15) guessing f > 2 IO/'(p v2). Eocg £0 < 1 Calculate w from Equation (3.14) Calculate f from Equation (3.16) and compare with the guess value in Step 8-b If Re > Rec, then transitional/turbulent flow 11. fest Guess a value for f > 210/(pV2) Figure 4. Calculation scheme to estimate the friction factor. 31 12. R 13. R If R < Rc’ then go to step 11. 14. go 15. B 16. 17. NJ 18. Cj 19. 20. 21. f Calculate RC from Equation (3.27) with Re = Re and f = f c c Calculate R from Equation (3.27) Guess a h1gher valve for fest Calculate go from Equation (3.32) or (3.15). 0 §.€0 < goc Calculate B from Equation (3.28) Calculate w from Equation (3.26) Generate values of xi from Equation (3.25) with go §_g §_1 (j = 1, 2, 3 ...) Generate values of Cj from Equation (3.30) with values of Ni and RiiyE1U=133-~) Evaluate the integral of Equation (3.31) by numerical methods with Cj and 50:53: _<_1 (j = 1,2,3 ...) Calculate Equation (3.31). If result # 1, then go to Step 11. If Equation (3.31) is equal to one, then f = fest Figure 4. Continued. where the friction factor f is obtained using Hanks' method described in the previous section and the friction loss coefficients, Kf, for fittings can be approximated with the Newtonian data for turbulent flow (Crane, 1982; Perry and Chilton, 1973; Govier and Azis, 1972) and the relationship given by Iwanami and Suu (1970), Steffe et al. (1984), and Soto and Shah (1976) for laminar flow. The mass average velocity may be written as g = 2 (3.34) Substituting Equation (3.34) into (3.33) gives the result . 02 02 _ 32fM L 8M E: 02 n p D n p D The annual cost of a pipe system can be estimated using Equation (2.17) as cpi = (a + b)cp0S (3.36) where C total annual cost of installed tube system per unit P‘ length of tube, $/yr m a = annual fixed cost of the tube system expressed as a fraction of the initial installed cost of the tube system, 1/yr b = annual maintenance cost of the tube system expressed as a fraction of the initial installed cost of the tube system, 1/yr C = empirical constant for the tube system cost, $/m1+S s = exponent of tube system cost equation, dimensionless 33 As stated before, Cp and s can be estimated from a log-log plot of the installed csot of the tube system (tube, fittings, valves, etc.) versus the tube inside diameter. Notice that Equation (3.36) can also be interpreted as Equation (2.15) if one lets s = p' and cp = (F + 1) x (39.37)P'. annual cost of the tube system as a function of the diameter from the This permits one to obtain the installed knowledge of the costs of one-inch tube and fittings. However, some error may be introduced by assuming F to be independent of D and extrapolating from the cost of one tube size. Therefore, this method should only be used in preliminary tube sizing when data and knowledge of the system are limited. More accurate results can be obtained if the variables Cp and s are calculated from the installed cost of the tube system for various tube diameters. Even in this case, extra- polating beyond the diameter range used should be avoided. That is, Equation (3.36) should be estimate using a range of diameters where the optimum diameter is expected. Notice that the annual fixed cost and the annual maintenance cost ratios are assumed to be independent of tube diameter. However, these costs, as well as other costs associated with the tube system which may depend on the tube diameter, may be included in the estimation of the installed cost of the tube system. Then the fixed and maintenance cost will be included in the varialbes Cp and S. If this is done, the term (a + b) in Equation (3.36) can be set equal to one. The cost of a pump station can be written, in a manner similar to Equation (2.19) as 34 cps = Co PS + CI (3.37) where CpS = total cost of installed pump station, $ CD = empirical constant for the pump station cost, $/WS CI = empirical constant for the pump station cost, $ exponent of the pump station cost equation, dimensionless (A II The value of CD, CI, and s' can be obtained from a plot of the installed cost of the pump versus the power require- ments. The installed pump station cost includes the purchase cost of pump, motor, and other costs dependent on the size of the pump. Equation(3.37) permits the use of a linear (s' = 1) or power (CI = 0) relationship for Cps versus P. Notice also that this equation can be interpreted as Equation (2.19) if one lets CI = Ci, P = Q1’ 51 = q, and CD = Céz’qu. The annual pump station cost per unit length of tube can then be expressed as - 1 1 5' Cpu -(a + b )(CDP + CI)/L (3.38) where a' = annual fixed cost of the pump station expressed as a fraction of the initial installed cost of the pump station, l/yr b' = annual maintenance cost of the pump station expressed as a fraction of the initial installed cost of the pump station, l/yr 35 Again, the fixed cost and the maintenance costs may be included in the estimation of the installed cost of the pump station for the different pump sizes accounting for these costs in the variables CI, CD, and S'. Then, the term (a' + b') in Equation (3.38) could be set equal to one. The total annual cost of a pumping system per unit length of tube can be obtained by adding Equations (2.20), (3.36), and (3.38) which gives, after rearrangement, S CehP (a' + b') (cDPS' + CI) CT = (a + 0) CPD + L Ceh P + 1 (3.39) Substituting Equation (2.13) into (3.39) yields c hnw E(a' + b')(C nS'wS 5'5 + c ) e D I + 1 c = (a + b)C US + _ T p LE CJIMW (3.40) where W = work per unit mass, J/kg, Equation (3.35) The procedure to estimate the pumping system costs is outlined in Figure 5. The computer program developed to accomplish these calculations is given in Appendix D. 3.3 Optimum Economic Tube Diameter As stated before, the optimum tube diameter, Dopt’ can be obtained by setting dCT/dD = 0 assuming that CT is only a function of D, i.e., 36 Imput Variables h, L, Ap,Az, 2Kf nsKsTsp o Cp,s, a, b CI, CD, 5', a', b' Ce’ h, E D or Dopt 1 l 2 f 3. W 4. P 5. Cpi 6. Cpu 7 Cop CT Pumping system parameters Fluid properties Tube system cost parameters Pump station cost parameters Operating Tube/pipe Calculate Calculate Calculate Calculate Calculate Calculate Calculate Calculate cost parameters inside diameter 9 from Equation (3.34) f from scheme in Figure 4 W from Equation (3.35) P from Equation (2.13) Cpi from Equation (3.36) Cpu from Equation (3.38) C0p from Equation (2.20) CT from Equation (3.39) or (3.40) Figure 5. Calculation scheme to estimate the annual costs of a pumping system. 37 d ( ) s—l 32n3ceh (a' + b')s'cDM5"1 —— C = a + b SC 0 — , + 1 ' dD T p opt TT2p2D6E Ceh E5 -1 2 d Dopt Dopt d 5f - DoptEfiIf + L sz - 4L ED'ZKf = 0 (3.41) For laminar flow, df/dD can be obtained from the derivative of Equation (3.16) with respect to the diameter which gives d _ 16 dRe 16 £1.11». — f — - — — (3.24) dB wRe2 dD Rewz dD Replacing V with Equation (3.34) in Equation (3.17) and taking the derivative of Re with respect to 0 gives 91%: = _l(3" [3 4 Re (3.43) Similarly, substituting go with Equation (3 15) in Equation (3.14) and taking the derivative of w with respect to 0 gives. dw psio d 440 £0 dD'=-.T——‘ ED f — "‘ff—-' (3.44) where C = [(1 + 3n)(1 + n)(1 - 5°)z + 2(1 + 2n)(1 + 3n)£o(1 - £0) + (1 + 3n)(1 + 2n)(1 + n)£§ I 1 (1 + 2n)(1 + n)(l - 5°)3 + 2(1 + 3n)(l + n)g° (1 - 5012+(1 + 3n)(1 + 2n):§ <1 - an) (3.45) 38 Substituting Equations (3.43) and (3.44) into Equation (3.42) and solving for df/dD gives d 64goo + 16 (4 - 3n) an I = RewD (1+ OED (3'46) 01" d 4f§oo + (4 - 3n)f d—D' f = D (1 + 05°) (3'47) where f, go and o are given by equations (3.16), (3.18), and (3.45), respectively. Equation (3.46) was confirmed for the special cases of the power law, Bingham plastic and Newtonian fluids by comparing it to independent analytical solutions for these fluids. It was also confirmed numerically for two examples of a H-B fluid (Appendix B). A numerical integration was required to estimate the friction factor for turbulent flow as seen in Equation (3.31); hence, the derivative of the friction factor with respect to the diameter must be approximated numerically for this flow condition by = f(D + x) - f(D—x) E UP) 2X (3.48) where f(0) = the friction factor expressed as a function of D x = a small positive number The backward difference method is used to evaluate the derivative for diameters just below the critical diameter where turbulent flow starts (Appendix C). Examples of f(D) versus 0 are shown in Appendix B. An alternative equation for df/dD for turbulent flow is given in Appen- dix A when the friction factor is estimated with the relationship developed by Torrance (1967). Since no general equation exists for the fitting resistance coefficient (Kf), it must be assumed to be independent of the diameter. That is, dKf/dD = o. In addition, when L > > o, L/D and 02/4L will be small numbers having a small influence on 00 Hence, constant pt' Kf values for Newtonian fluids for turbulent flow can be used as an approx1mat1on to evaluate Dopt' Then, the DOpt 1s g1ven 1mpl1c1ty by eliminating the dKf/dD term and solving for DOpt from Equation (3.41) as Ds+5 32M3ceh (a' + b')s'-cD PS"1 = . + 1 Opt (a + b)stn2p2E Ceh 5f - o -94 +-EQEE—zm< (3 49) opt dD L f ' where P = power, watts, Equation (2.13) By letting CD = 0, ZKf = 0, go = 0 and cp = (1 + F)(39.37)SX, and substituting df/dD for Equation (3.46) or (3.47), Equation (3.49) reduces to 40 1 (331.71) 110 4(1 + “IcthK (8:4(1 + 3n)) ’1 pt . s(a + b)(F + 1)(39.37)S x QE "8" (3.50) which is an equivalent form of Skelland's equation (Skelland, 1967; p. 245) for power law fluids in laminar flow, but with the variables expressed in SI units. The procedure to estimate the optimum diameter is shown in Figure 6. The computer program to do these calculations is given in Appendix D. 3.4 Limitation of Design Method Some assumptions are inherent in the design method presented. Even though some of these assumptions were stated previously, they will be summaried here: 1. Use of Newtonian values for the fittings resistance coefficients 2. Constant fluid density (incompressible fluid) Homogeneous or pseudohomogeneous fluids #0) No slip or apparent slip at the wall 5. No elastic or time-dependent behavior 6. Smooth wall (for turbulent flow only) 7. Isothermal flow 8. Steady state flow 9. Negligible kinetic energy change 41 boom Imput variables M, L, Ap, Az, 2Kf II, K9 T0: D C , s, , b p a C13C35959b L Pumping system parameters Fluid properties Tube system cost parameters Pump station cost parameters Operating cost parameters Guess DOpt Calculate 9 from Equation (3.34) Calculate f from scheme in Figure 4 Calculate df/dD from Equation (3.46) or (3.47) if the flow is laminar (Re < Rec) or from Equation (3.48) if the flow is turbulent (Re > Rec) Calculate W from Equation (3.35) Calculate P from Equation (2.13) If Equation (3.49) is true then Dopt = Dest’ otherw1se go to Step 1 Figure 6. Calculation scheme to estimate the optimum diameter. 42 The first assumption may be violated when handling non- Newtonian fluids under laminar flow. Under this flow condition, the fittings resistance coefficient increases with decreasing Reynolds number (Steffe et al., 1984). The frictional loss in fittings may become significant in a complex tube system with a great number of fittings. This may cause errors in estimating the optimum diameter, particularly with short tube systems. The next three conditions may be violated in the handling of heterogeneous or multiphase fluids. In these systems, a particle- free layer may form at the pipe wall creating a variation of solids concentration. The lubricating action of this liquid layer is known as effective slip. These systems cannot be accurately described by the H-B fluid model. More complex models are also required to describe the flow behavior of viscoelastic and time-dependent fluids. Viscoelastic fluids show partial elastic recovery on removal of deforming shear stresses. Such materials exhibit both viscous and elastic properties. Time-dependent fluids exhibit reversible decrease (thixotropic), irreversible decrease (rheomalaxis) or reversible increase (rheopectic) in shear stress with time at constant rate of shear (Skelland, 1967). These and various time-independent rheologi- cal models not described by Equation (2.1), such as the Ellis and Casson models, are not considered in this study. For turbulent flow, wall roughness leads to increased pressure drop (Cheng, 1975). Therefore, the power requirements will be underestimated for this condition since the pressure-drop/flow- rate relation used in this study for turbulent flow (Sections 3.1.3) is applicable only for smooth walls. This, in turn, will lead to under estimation of D Nonisothermal conditions will cause errors 0 t' in the design method since the consistency coefficient (K) depends on temperature. It decreases with increasing temperature according to the Arrhenius relationship (Cheng, 1975). Nonisothermal conditions may be caused by changes in environmental temperature or by mixing of various streams at different temperature. In addition, unsteady flow conditions may be encountered in start-up operations. Also, pressure surge waves may develop in long pipeline due to fluid inertia and compressibility (Cheng, 1975). Finally, appreciable kinetic energy changes may be found in complex tube systems with variation of tube diameter, entrances, fittings, etc. The current design method is not applicable for such systems. 4. RESULT AND DISCUSSION 4.1 Model Verification To validate the model, the optimum diameter (00 ) was first t estimated for the example given by Skelland (1967) (IllUstration 7.1 (c) pp. 253) for a power law fluid. His data, in terms of the vari- ables and units of the model developed in this study, are given in Table 1. The Dopt using this model was found to be 0.1653 (0.5425 ft) which is the same as the one obtained from Skelland's optimum diameter equation for power law fluids (as defined by Equation (2.2)). Notice that the answer in his illustration is different and is due to round- off error in his numerical constants. Direct calculation using his original equation (pp. 245) gave the same result. This is to be expected since Equation (3.49) is a general form of Equation (3.50) which is an equivalent form of Skelland's equation. The model was also tested using the example from Darby and Melson (1982) for a Bingham plastic fluid in turbulent flow. Their data are given in Table 2. for this case was found to be 0.865m which is 0.82% higher pt than their value. This small deviation is probably due to the fact The D0 that Darby and Melson assumed df/dD = 0 and used an approximation (Darby and Melson, 1981) of the friction factor relationship of Hanks;and Dadia (1971) in the derivation of their model. Even though the result of Darby and Melson was close to the one obtained using the model of this study, their assumption of df/dD = 0 may introduce 44 45 Table 1. Fluid properties and other pertinent data for the optimum diameter example problem given by Skelland (1967) Fluid Properties 0.5; K = 3.02 Pa 5”; p = 977.29 kg/m3 3 II Pipe Cost Parameters 354.58 s/m2'5; s = 1.5; a = 0.14; b = 0.06 C') II Power Cost Parameters Ce = 2.0 x 105 $/w hr; n = 6570 hrs/yr Other Pertinent Data n = 13.83 kg/s; L = 1523.93 m; E = 0.3 Optimum Pipe Diameter D = 0.1653 m opt 46 Table 2. Fluid properties and other pertinent data for the optimum diameter example problem given by Darby and Melson (1982) Fluid Properties n = 1.0; n = 0.03 Pa 5; To = 4.2 Pa; 0 = 1400 kg/m3 Pipe Cost Parameters cp = 409.58 $/m2'2; s = 1.2; a - 0.05 Pump Station Cost Parameters CI = 173800 $; CD - 0.6 $/W; s' = 1.0; a' = 0.05 Power Cost Parameters ce = 4.0 x 10"5 $/w hr; n = 8640 hrs/yr Other Pertinent Data M==729.16 kg/s; L = 1m; E = 0.6 Optimum Pipe Diameter D = 0.858 m opt 47 significant error for other fluid properties and flow conditions. This assumption is based on df/dD being much smaller than 5f/Dopt for the term 5f/Dopt-df/d0 which appears in Equation (3.49) if the equation is divided by D However, for power law fluids in laminar flow, 0 t' df/dD can vary ftom 20% to 5f/D for n = 1 to 74% of 5f/D for n = 0.1. For H-B fluids, df/dD was found to be as much as 68% of Sf/D. Therefore, the assumption of df/dD << 5f/D is questionable. In addition, the model was verified by comparing the Dopt obtained analytically (Equation (3.49)) and graphically (Figure 2) as will be shown later. 4.2 Cost Parameters and Other Pertinent Data for a Pumping System Consider a pumping system consisting of 100m of 304 stainless steel tubing with both ends at the same pressure and elevation. The tube system includes three tees (used as elbow), three 90° elbows, twenty-one union couplings, and two plug valves giving an overall fittings resistance coefficient of 10. A close coupled sanitary centrifugal pump is to be used with pump and motor (combined) effi- ciency of 70%. The variation of the installed tube system costs per meter length of tube with tube inside diameter are shown in Table 3. These are plotted on log-log coordinates in Figure 7. As seen, a straight line described by Equation (2.17) gives the constant Cp and 5 shown in the figure and a regression coefficient of 0.95. The values of CD and s are also given in Table 4. In addition, the fixed (a) and maintenance (b) annual cost ratios for the tube system are presented. The values of a was estimated from Equation (2.16) assuming 48 Table 3. Variation of the installed cost of a tube system per meter length of tube with tube diameter. Estimated from the purchase cost (January 1985) of 100 m of Tri-Clover 304 Stainless Steel Tubes (1-3 in tubes are gauge 16, 4 in tube in gauge 14), polished ID + 00; 3 tees (7MP); 3 90° Elbows (2CMP); 3 Caps (16AMP); 2 plug valves (DIOMP); 30 Gaskets (40MP-U); 30 Clamps (13MHHM); 36 Furrales (14RMP). Ladish Company, Tri-Clover Div., Kenosha, Wis- consin. Installation costs approximated with 1.5 man-hr/ joint-diameter (in) relation (Jeler, 1970), and labor cost of $35/man hr. Diameter Installed Cost 00 (in) ID (m) $/M 1 0.0221 36.53 1% 0.0348 44.24 2 0.0475 56.57 2% 0.0602 76.35 3 0.0729 94.65 4 0.0974 144.71 49 ~93 I C = 1097 0° for 0.0221 < D < 0.0974m Correlation Coefficient = 0.95 g 100'- E b \ D (D E" b 0: <3 0 . Q m a q I < 9 in E b 10 L l 1 L l 1 1 1 l 0.01 0.1 TUBE INSIDE DIAMETER (m) Figure 7. Variation of the installed cost of the tube system (Table 3) with tube inside diameter. 50 Table 4. Cost parameters for the tube system presented in Table 3. Based on Figure 7 and Equations (2.18) and (3.36) cp ($/m$+1) 1097 s 0.93 a (l/yr) 0.18 b (1/yr) 0.10 51 an interest rate of 12% and lifetime of 10 years. The value of b was taken as 10% of the installed cost of the tube system. The variation of the installed costs of the pump station with pump size are shown in Table 5. These are plotted in Figure 8. As seen, a straight line (5' = 1), described by Equation (3.37) gives the constants CI and CD shown in the figure with a regression coeffi- cient of 0.96. The valves of CI, CD, and s' are also shown in Table 6 along with fixed (a') and maintenance (b') annual cost ratios for the pump station. An interest rate of 12% and lifetime of 5 years were used to estimate a'. As for the tube system, b' was taken as 10% of the installed pump station. In addition, the system is to be operated 75% of the year (6,570 hrs/yr) and the electrical energy cost if 0.06 $/kW hr. These and other pertinent data are tabulated in Table 7. 4.3 Optimum Diameter for a System Handling Tomato Ketchup It is desired to determine the most economical diameter (D ) for transport of tomato ketchup at a mass flow rate of 4.0 kg/s. opt This fluid can be considered to be a homogeneous non-Newtonian fluid described by the H-B model (Higgs and Norrington, 1971). The fluid properties at 25°C are given in Table 8. Using given variables (Tables 4, 6, 7, and 8), the 00 minimum pt cost (C ) power (P), work (W), and pumping system costs (Cpi’ C 9 T . .m1n Cop) at optimum were estimated using the procedure outlined in Figures 5 and 6. The results are summarized in Table 9. Figure 9 [N 52 Table 5. Variation of the installed pump station cost with power requirements. Estimated from the purchase1c09t(January, 1985) of Tri-Flo close-coupled sanitrary centrifugal pumps, C216 (with water cooled rotary seal) and electric motor (60 cycle 230/460 volt-3 phase), 1750 rpm for 1-2 Hp pumps and 3500 rpm for 1-15 Hp pumps. ("Easy- Clean" totally-enclosed motor), Ladish Company, Tri- Clover Division, Denosha, Wisconsin. Installation cost taken as 25% of the total purchase cost (Peters and Timmerhaus, 1968). Power Installed Cost Hp watts ($) 1/2 372.9 1313 3/4 559.3 1348 1 745.7 1366 1 l/2 1118.6 1384 2 1491.4 1484 2 1491.4 1414 3 2237.1 1791 5 3728.5 1876 7 1/2 5592.8 2175 10 7457.0 2225 15 11185.5 2686 INSTALLED COST (5) 53 C = 1308 + 0.13 P ps Coorelation Coefficient = 0.96 I 3000 2500 2000 1500 1000 L 4 A I 1 2000 4000 6000 8000 10000 POWER (watts) Figure 8. Variation of the installed cost of the pump station (Table 5) with power requirements (linear relationship). 54 Table 6. Cost parameters for the pump station presented in Table 5. Based on Figure 8 and Equations (3.37) and (3.38) CI ($) 1308 CD $/W 0.13 s' 1.0 a' (l/yr) 0.28 b' (l/yr) 0.10 Table 7. Electrical energy cost, hours of operations per year, combined pump and motor efficiency, summation of the fittings resistance coefficients, tube length, pressure and elevation change, and mass flow rate used to estimate the costs, and optimum diameter for the pumping system presented in Table 3 and 5 ce ($/w hr) 6.0 x 10'5 h (hrs/yr) 6570 E 0.70 2k, 10.0 L (m) 100.0 49 (Pa) 0 02 (m) 0 M (kg/s) 4.0 55 Table 8. Rheological properties (Higgs and Norrington, 1971) and density (Lopez, 1981) for tomato ketchup at 25°C n 0.27 K (Pa s“) 18.7 To (Pa) 32.0 p (kg/m3) 1110.0 Table 9. Optimum economic tube diameters, pumping system costs, and work and power requirements, at optimum for a system (Tables 4, 6, and 7) transporting tomato ketchup with properties given in Table 8 Dopt (m) 0.06907 CTmin ($/yr m) 45.70 Cpi ($/yr m) 25.58 Cpu ($/yr m) 6.66 COp ($/yr m) 13.46 W (J/kg) 597.81 P (k Watts) 3.42 COST ($/ yr m) 56 80 70 Total cost, CT 60- ..( - ---- --------- C . +2 E:CT . ‘17 Tmin % min | 40" l Operating cost, COp I I 30- I ' Tube system cost, Cpi I 20- I I I 10_ Pump station cost, Cpu D 1 1 . 09th 1 1 . 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.10 TUBE INSIDE DIAMETER (m) Figure 9. Variation of tube system cost, pump station cost, operating cost, and total cost with tube inside diameter for a system (Table 4, 6, and 7) transporting tomato ketchup with properties given in Table 8. 57 shows the variation of the costs with tube inside diameter. As seen, Dopt obtained graphically and analytically was 0.06907m for a CT of 45.70 $/yr m or 5¢/ton of tomato ketchup pumped annually. m1” P was found to be 3.42 kW (4.58Hp). It can also be seen that tube diameter between 0.0585m and 0.0835m results in total costs which does not exceed the minimum value by more than 2% and that the deviation from minimum increases more rapidly as the diameter decreased. The Reynolds number was found to be less than the criti- cal Reynolds number for Dopt’ hence the flow was laminar as seen in Table 10. These values (Table 10) were estimated following the scheme shown in Figure 4 with the program given in Appendix 7.4. 4.4 Sensitivity Analysis This section is devoted to study the sensitivity of Dopt on the various input variables shown in Figure 6. This was done by estima— tion the percent change of Dopt obtained using a 110% value of each variable. Even though the analysis is mostly based on the example of Section 4.3, some general insight can be obtained on the relative importance of the cost components of the pumping system other variables in determining D The percent changes of 00 for each.of the opt' pt variables are shown in Table 11. The change in the Cpu variables, with the excpetion of 5', resulted in small changes of Dopt’ This is due to the small variation of CpS with P obtained in Figure 8. Changing s' just changed the nature of this relationship. The varia- tion of the variables of Cpi and COp resulted in considerable changes on the 00 can also be of Dopt° The greater 1nfluence of Cpl and C0 pt P 58 Table 10. Flow condition of Dopt = 0.06907 m for the pumping system handling 4.0 kg/s of tomato ketchup with properties given in Table 8 Re 107 He 8.85 f 0.2214 go 0.2815 Re 2754.0 c f 0.007506 Table 11. Percent change of 00 pt for the example presented in Section 4.3 with +10% change of imput variable Percent change of variable -10% +10% Variables Percent Change of Dopt Fluid Properties 0 -5.11 +5.17 K -3.04 +2.88 To -0.80 +0.78 0 +5.07 -4.39 Pumping System M -4.91 +4.65 L +0.07 -0.07 ZKf -0.09 +0.07 AP -- -- AZ -- -- Tube system cost C +4.11 -3.58 (Cpi) sp -5.52 +5.82 a +2.56 -2.36 b +1.39 -1.33 Pump Station Cost CI -- -- (Cpu) CD -0.43 +0.42 5' -2.62 +5.94 a' -0.32 +0.30 0' -0.12 +0.10 Operating Cost E +4.11 -3.58 Cop) Ce -3.49 +3.30 -3.49 +3.30 60 not1ced 1n Figure 9. Cpu has less 1nfluence on Dopt compared to CO An economic balance on Cpi and C0 . P° p__along gave a DOpt of 0.06602 m wh1ch is 4.42% lower than the value found when Cpu was considered. For the example to Table 2, DOpt excluding Cpu was found to be 0.855 m, 1.17% lower than the value found when Cpu was included. The changes of a and 0 did not effect the results as much as Cp and s. A change of +10% in a represents an appropriate change of +20% in the life-time (N) or 125% change in interest rate (1). These variations on N and i result in less than :3% change on Dopt' From Equation (3.49), it can be observed that where using a linear relationship (5' = 1) for Cps vs. P, the DOpt is independent of the pressure energy change (Ap) and the elevation change (82). Even if s' is not equal to one, D t can 0P be assumed to be independent of Ap and 82, since C generally varies pu little with P. To show this, the data in Table 5 were fitted to the curve shown in Figure 10. The value of CD, CI’ and s' for this curve are given in Table 12. The Dopt using these constants and A2 = Ap = 0 was found to be 0.06902 m which is only 0.07% lower than the value found using the linear relationship of Figure 8. For 82 = 20 m and pt was found to be 0.46% lower. This variation was also obtained for A2 = 0 and Ap = 217.78 kPa (2.15 atm) which shown Ap = 0, the Do the small influence of Ap and 82 have on Dopt’ The changes on L and ZKf also produced small changes on Dopt; hence, using the constant Kf value of Newtonian fluids in turbulent flow for approximating non—Newtonian fluid behavior will introduce negligible error. As seen from Equation (3.49), Dopt can also be INSTALLED COST (5) 3000 2500 2000 1500 1000 61 0-6 cps: 1075 + 5.53 P Correlation Coefficient = 0.95 l I k l L 2000 4000 6000 8000 10000 POWER (watts) Figure 10. Variation of the installed cost of the pump station (Table 5) with power requirements (non-linear relationship). 62 Table 12. Cost constants of Equation (3.37) for the pump station presented in Table 5 based on Figure 10 CI ($) 1075 CD ($/w5') 5.53 s' 0.60 assumed to be independent of Kf if the tube length is much greater than the DO expected. Otherwise, the summation of the resistance pt coefficient per unit length of pipe (ZKf/L) can be used as an appro- ximation without greatly effecting the results. Notice also that if s = 1 and L>>DOpt or an approx1mat10n of ZKf/L 15 used, 00 t 15 also P independent of L. In other words, Dopt can be estimated from the costs of a unit length of tube/pipe (e.g., one meter). The small effect of the error of Ap, Az, ZKf, and L on Do is of great value pt since these variables are usually not well known in preliminary sizing of a pipe system. Finally, as seen in Table 11, :10% change in the fluid properties and mass flow rate, except for the yield stress, resulted in considerable change on Dopt‘ 1 5 0 I' D' I H . g I M' '!i s The optimum diameter (Dopt)’ pumping system costs, power requirements and work requirements to transport tomato ketchup were estimated assuming Newtonian flow behavior and apparent Viscosities (pa) of 4.723 Pa 5 and 1.715 Pa sto show the problems that may arise 63 from this practice. The value of pa were calculated from Equation 1 and 50 s.1 (2.6) using a rates of shear of 15s- , respectively, to simulate point measurement (such as those which might be made with a Brookfield viscometer) at these shear rates. The results, using pa = 4.723 Pa 5, are shown in Table 13 and compared to the results of Table 9 (Section 4.3). As seen, as 00 was over estimated signifi- pt cantly given a CT 20.63% higher and a P value 37.56% lower. However, this Dopt does no$13ive the actual minimum as seen in Table 14 which illustrates the actual pumping system costs, and work and power requirements estimated using the H-8 model at D = 0.1138 m (the Dopt obtained using pa = 4.723 Pa 5). As seen, 0 = 0.1138 m gives a total cost which deviates from the minimum by 15.4%. In addition, the actual power requirement is 33.78% lower than the one estimated using the apparent viscosity. If the pumping system was designed using this apparent viscosity (4.723 Pa 5), the tube system cost (Cpi) and pump station cost (Cpu) would be estimated to be 40.7 $/yr m and 6.02 $/yr m for a tube size and pump size of 0.1138 m and 2.13 kWatts, respectively (Table 13). However, the operating cost would be 6.28 $/yr m since the actual power requirement at D = 0.1138m is 1.59 kWatts (Table 14). So the total cost for this system would be 53.0 $/yr m which is 15.97% higher than the one in Table 9. In addition, the system would have a oversized (hence, less efficient) pump. Table 15 shows the results obtained using pa = 1.715 Pa 3 and a comparison with the values of Table 9. The Dopt for this case was found to be 0.09271 m, 34.23% higher than the one obtained in Table 9. 64 Table 13. Optimum economic tube diameter, pumping system costs, and work and power requirements at optimum estimated assuming Newtonian flow behavior and an apparent viscosity of 4.723 Pa 5 Results for % Difference with _ the results of Dopt (m) 0.1138 +64.76 CT . ($/yr m) 55.13 +20.63 m1n Cpi (S/yr m) 40.7 +59.11 Cpu ($/yr m) 6.02 - 9.61 Cop ($/yr m) 8.41 -37.56 W (J/kg) 373.25 -37.56 p (kWatts) 2.13 -37.56 Table 14. Pumping system costs, and work and power requirements estimated using the H-8 model for D = 0.1138 m CT ($/yr m) 52.74 Cpi ($/yr m) 40.7 Cpu ($/yr m) 5.76 COp ($/yr m) 6.28 w (J/kg) 279.0 P (kWatts) 1.59 65 Table 15. Optimum economic tube, diameter, pumping system costs, and work and power requirements at optimum estimated assuming Newtonian flow behavior and an apparent vis- cosity of 1.715 Pa—s Results for % Difference u = 1.715 Pa-s with Results 3 of Table 9 Dopt (m) 0.09271 +34.23 CT ($/yr m) 46.42 + 1.58 min Cpi ($/yr m) 33.63 +31.47 Cpu ($/yr m) 5.84 -12.31 Cop ($/yr m) 6.95 -48.40 W (J/kg) 308.49 -48.40 P (kWatts) 1.6 -48.40 66 Again, this 00 t does not give the actual minimum as seen in Table 16. In this case, the total cost deviated from minimum by 5.51%. From Tables 15 and 16, it can also be seen that the actual power requirement at D ==0.09271 m is 18.7% higher than the one estimated using the apparent viscosity. If the tube size and pump size were to be selected, base or Table 15, the pumping system would have a undersized pump uncapable of meeting the actual operating conditions. Therefore, the pump would have to be replaced or the operating time would have to be increased resulting in a more expensive system. As these two examples show, the use of pa and Newtonian flow behavior to design non-Newtonian handling systems may lead to errors depending on the rate of shear at which pa was measured. 4.6 Optimum Diameter for a System Handling a Herschel-Bulkley Fluid in Turbulent Flow A problem was selected to test the optimum diameter model for a H-B fluid that resulted in a DOpt for which the flow was turbulent. For this purpose, the Dopt was estimated for a hypthetical H—B fluid with properties given in Table 17 and the cost data shown in Tables 4, 6. and 7. The results at C are shown in Table 18. The flow T . m1n condition for Dopt was found to be turbulent (Table 19) and the variation of the costs with D is shown in Figure 11. The Dopt obtained graphically and analytically was found to be 0.04065 m for a C of 24.16 $/yr m confirming Equation (3.49) for turbulent T . m1n flow. The diameter range, which CT did not exceed the CT by more min 67 Table 16. Pumping system costs, and work and power requirements estimated using the H-8 model for D = 0.09271 m CT ($/yr m) 48.22 Cpi ($/yr m) 33.63 Cpu ($/yr m) 6.04 C0p ($/yr m) 8.55 W (J/kg) 379.46 P (kWatts) 2.17 Table 17. Rheological properties and density for a hypothetical Herschel-Bulkley fluid n 0.70 K (Pa s”) 0.03 To (Pa) 2.0 9 (kg/m3) 1400.0 68 Table 18. Optimum economic tube diameter pumping system costs and work and power requirements, power at optimum for a system (Table 4, 6, and 7) transporting a H-B fluid with properties given in Table 17 DOpy (m) 0.04065 CT ($/yr m) 24.16 min Cpi ($/yr m) 15.62 Cpu ($l/yr m) 5.37 C0p ($/yr m) 3.17 W (J/kg) 140.68 P kWatts 0.80 Table 19. Flow condition at Dopt = 0.04065 m for the pumping system handling 4.0 kg/s of a H-B fluid with properties given in Table 17 Re 2.40 x 104 He 1.88 x 105 f 0.004883 50 0.1207 ReC 6.34 x 104 f 0.08200 90 80 7O 60 50 40 cosr ($/yr m) 30 20 69 Total cost, CT Tube system cost, Cpi I I I I I I I l I I I I I L 10 Pump station cost, Cpu DOpt Operating cost, COp 1+ 1 J 3 ==:_. 0.02 0.03 0.04 0.05 0.06 0.07 TUBE INSIDE DIAMETER (m) Figure 11. Variation of tube system cost, pump station cost, operating cost, and total cost with tube inside diameter for a system (Tables 4, 6, and 7) transporting a H-B fluid with properties given in Table 17. 70 than 2%, was 0.0365 m to 0.046 m, which is smaller than the example of Section 4.3. As observed in the previous example, the total cost deviated from minimum most slowly as the diameter increased (Figure 11). This rate of increase is practically given by Cpi as seen from the similarity of the sl0pes of the CT and Cpi curves for D > Dopt' 5. SUMMARY AND CONCLUSIONS 1. An equation to determine the total annual cost of a pumping system as a function of tube diameter (based on the costs of the tube system, pump station, and operation) has been developed for system handling Herschel-Bulkley fluids under laminar, transitional, or turbulent flow condition. 2. An equation to determine the optimum economic tube diameter has been developed for pumping systems handling Herschel- Bulkley fluids under laminar, transitional, or turbulent flow condi- tions. 3. The pump station cost had less influence than the operating cost in determining the optimum economic tube diameter. 4. The optimum economic tube diameter is independent of any elevation difference (A2) and pressure energy difference (Ap) in the system if a linear relationship (5‘ = 1) is used to correlate the pump station cost to power requirements. In addition, 82 and Ap do not have to be known accurately if the variation of the pump station cost with power is small. 5. The optimum ecnomic tube diameter can be obtained from the pumping system costs of a unit length of tube if a linear rela- tionship is used to correlate the pump station cost to the power reQuirements (s' = 1) and the length of the tube system is much 71 greater than the tube diameter (L >> Dopt) or the frictional loss in fittings and values is approximated as the summation of the fittings resistance coefficient per unit length of tube. 6. The use of apparent viscosity and Newtonian flow behavior for non-Newtonian fluids caused significant errors in the estimation of the optimum economic tube diameter, total annual cost, and power requirements of the pumping system. APPENDICES 73 APPENDIX A ALTERNATIVE EQUATIONS FOR f AND df/dD FOR H-B FLUIDS IN TURBULENT FLOW 74 APPENDIX A ALTERNATIVE EQUATIONS FOR f AND df/dD FOR H-B FLUIDS IN TURUBLENT FLOW Torrance (1967) developed a friction factor relationship for H-B fluids in turbulent flow as _l_ = _ 2.275 1.97 _ e J1? 0.45 -—n +-——n tn (1 go) n n n + 1.97 in [Re 1+3n fI-I'I/Z] (A.I) where Re and 50 are given by Equations (3.17) and (3.15), respectively Combining the definition of f, Re, and He, Equation (3.15) can be rewritten as 81. 2 2-n 2-n < n ) £0 = 16 (2He) 2 1+3n . (A.2) Rez:fi- f Equations (A+1) and (A+2) gives the friction factor as a function of Re, He, and n. Replacing V with Equation (3.34) in Equation (3.17) and go with Equation (3.15) in Equation (7.1), the derivative of f with respect to D from Equation (7.1) gives 75 76 3.94 [4 - 3n (1 - g0)]f3/2 [394+i + (1.-1.97ii)n(1-g0):| 0 ale. :3 -h This equation can be used instead of Equation (3.48) if the Torrance relationship (Equation (A.1)) is used to estimate the fanning fric— tion factor in turbulent flow. However, it is not clear what laminar- turbulent criterium should be used with the Torrance equation. Equation (A.3) was confirmed in the same manner as Equation (3.46). APPENDIX B VERIFICATION OF df/dD EQUATION FOR LAMINAR FLOW 77 APPENDIX B VERIFICATION OF df/dD EQUATION FOR LIMINAR FLOW The friction factor relation for power law fluids in laminar flow is given by Equation (3.16) with wIE = 1 as o=1 where Re is defined by Equation (3.17) When Equation (8.1) is differentiated with respect to D, it yields —0 ReD (3'2) When a value of zero for go (To = 0) is substituted into Equation (3.46), it reduces to Equation (B.2), indicating that Equation (3.46) is correct for the special case of the power law fluid. Equation (8.2) was also obtained by Darby and Melson (1982), and indirectly by Skelland (1967). For the Bingham plastic fluid, the friction factor is given by Equation (3.16) with K = n, and Re and w evaluated at n = 1 as 78 79 = Re|n=1 w'n=1 (8.3) where ReIn=1 = BEEF (3'4) and wIn=1 = 1 "3 £0 +'3 6: (8'5) When Equation B.6) is differentiated with respect to D, it gives 4 (1 - 52) df - o . ' 0 ReIn=1 LlJIn=1 D 1+ [3 wl _ :lgo n—l If one evaluates Equation (3.45) at n=1, it can be shown that 4 <1- e3) 0 n=1 _ 3 wln=1 (3-7) Then if n=1 is substituted hiEquation (3.46), it reduces to Equa- tion (B.6) which shows that Equation (3.46) is correct for the special case of the Bingham plastic. Similar results are found when considering the solution for a Newtonian fluid. In addition to the method just outlined, Equation (3.46) was confirmed numerically using Equation (3.48). The properties of tomato ketchup (Table 8) and the H-B fluid given in Table 17 along with mass flow rate of 4.0 Kg/s 80 were considered. These fluid properties and flow condition were also used for the examples given in Section 4.3 and 4.6, respectively. The variations of the friction factor with diameter for tomato ketchup and the H-B fluid are given in Figures 12 and 13, respectively. These curves can be obtained using the scheme shown in Figure 4. The analytical value obtained for df/dD at D = 0.06907m (the Dopt found in Section 4.3) was found to be 11.0585. Using x = 0.0001 in Equation (3.48), the numerical value was found to be 11.0586 which is only 0.001% higher than the analytical one. For the H-8 fluid (TablelJO, the analytical value of df/dD at D = 0.1m was found to be 1.045, 0.02% lower than the numerical value (1.0452) obtained using x = 0.001. It is clear that the analytical results are very close to the numercial results and small differences can be attributed to the limitations associated with the numeriCal solutiOn technique.- 81 .m\0x 0.3 no name 30am mmme m 0:6 0 magma cw co>mm mofiuumaoua omsam mnu no“ youwsmwo mummcm onau mamum> Ecuumm :omuumpm .NH whammm AEV mmbmzmm mowuumaouu omaau onu no“ nouoEmfio mowmcw onau msmum> EOuUMH cofluowum .mH musmwm A50 mmemz/TDL ER COMMON/RDOTNO/NOROOT CALCULATION OF U FROM 50. *MFR/(PI-OE‘DI--2) CALCULATION OF RE FROM EQ.(3-17) . 3.*N+1.))-~N = DI 2. )EFN =UF'(2. N) 8. *DERRET-RE2 RES/K i ) 1 F TOE N/Y YS.DE/FLWIDX/N 2/T OLC (3-34) U=4. 000 wwwwwwwwwwnnnmwnnn 000 0 Otom4mmwa-somm4mmbww unuuuvuuunuuuuuuuunuunnun 4444444444444444444 h m 444 1:515 (ON-t fl mggflwM-Aomm4mmbww4 uInuiunullnnIIunutnuuvnu 000 000 00040 44444 444mm l\)..bO(.O(X’\l nunnu" O 0 O 97 CALCULATION OF HE FROM EO (3-19) IF(Y$. E0. 0. ) THEN EOCOO. RE=(DE/YS)'DI--2'(YS/K)-*(2./N) END IF IF(HE.E0.0.) GO TO 10 CALCULATION OF EOC THROUGH ITERATION FROM EO.(3'22) HXsHE A1=O. A2=.999999999 CALL BISNEwT(A1 A2.40.TOLC.FUN1.DFUN1.EOC) IF(NOROOT.EQ. 0) GO TO 10 PRINT- PRINT~ gg%m1-.' THE DIMENSIONLESS UNSHEARED PLUG RADIUS.EOC.wAS NOT' ggIfiT-.' FOUND IN THE RANGE ’.A1.’ <= EOC <= ’,A2 PRINT'.’ ENTER A NEw RANGE FOR EOC: 0. <= EOC < 1.0 ...... ' READF.A1,A2 NOROOT=o GO TO 5 CONTINUE CALCULATION OF REC FROM Eo.(3-21) REC1=33600. -SORT(1. /27. )- N/(1.+3.*N)'*2 REC2=E2. +N)-'((2. +N)/ 1 +N)) REC3= -E0c)-'(1 +2. N) CALCULATION OF PSI-CRITICAL FROM EO.(3-14) wITH EO=EOC PSIC=Y(EOC) CONTINUE CALCULATION OF REC ES};REC2 PSICF'(2. /N)/REC3 R EC IT= 1. -REC/RE RECP = REC= R DIMC CR RETURN END 77 98 SUBROUTINE BISECT1(XA,XB.MAX.ERROR.NEWX) THIS SUBROUTINE COMPUTES THE OPTIMUM DIAMETER FROM EO.(3-49) (COMPUTES THE ROOT OF THE FUNCTION OPTDIAM). IT IS A COMBINATION OF THE BISECTION AND SECANT ITERATION METHOD. THE BISECTION INTERVAL IS USED TO START THE SECANT ITERATION. THE PROGRAM CONTINUES WITH THIS METHOD UNTIL THE SOLUTION IS FOUND OR THE FOLLOWING SITUATIONS OCCUR: 1- X FALLS OUTSIDE THE INTERVAL KNOWN TO CONTAIN THE SOLUTION; 2' X IS OUT OF RANGE OR INDEFINITE; 3- TOO FAR AWAY FROM THE SOLUTION; 4- THE NUMBER OF ITERATIONS EXCEEDS MAX. IF THESE SITUATIONS OCCUR. THE PROGRAM SWITCH TO THE BISECTION METHOD TO OBTAIN A SMALLER INTERVAL. REFERENCE: MOORE.E. 1982. ”INTRODUCTION TO FORTRAN AND ITS APPLICATION“. ALLYN AND BACON. INC.. BOSTON.MASS. LIST OF VARIABLES: DIFF= DIFFERENCE BETWEEN TWO ITERATION POINTS DIX= DIAMETER ARRAY GENERATED DURING THE OPTIMUM DIAMETER ITERATION ERROR= TOLERANCE ERROR FA= VALUE OF OPTDIAM AT XA FB= VALUE OF OPTDIAM AT XB FFX= FANNING FRICTION FACTOR FFY= FRICTION FACTOR ARRAY GENERATED DURING THE OPTIMUM DIAMETER ITERATION FM= VALUE OF OPTDIAM AT XM F03 VALUE OF OPTDIAM AT XO F1= VALUE OF OPTDIAM AT X1 LM= -1 IF X IS INDEFINITE; +1 1‘ OUT OF RANGE; O OTHERWISE MAX: MAXIMUM NUMBER OF SECANT ITERATION NEWX= ROOT OF OPTDIAM NO= NUMBER OF DATA POINTS GENERATED DURING THE OPTIMUM DIAMETER ITERATION NOROOT= ROOT INDICATOR: O-YES: 1-NO: 2-TOO MANY ITERATION X= POINT FROM THE SECANT ITERATION EQUATION XA= LOWER BOUND POINT USED IN THE BISECTION METHOD X8g UPPER BOUND POINT USED IN THE BISECTION METHOD XM= MIDPOINT BETWEEN XA 8 X8 X0: SECOND POINT REQUIRED FOR THE SECANT PROCESS X1= FOCUS POINT FOR THE SECANT ITERATION LIST OF SUBROUTINES: SWAP= INTERCHANGE THE VALUE OF TWO VARIABLES LIST OF FUNCTIONS: _ FRICFAC= FRICTION FACTOR FUNCTION OPTDIAM= OPTIMUM DIAMETER FUNCTION.EQ.(3-49) REWRITTEN AS OPTDIAM(X)=O. REAL NEWX COMMON/RD COMMON/FR COMMON/FF NO=2 IF(XA.GT.XB CALL SWAP(XA.X END IF FFx-FRICFAC(XA FFY§1 sFFx ,FFY; DTNO ICFC/ VSDI/ ) T 2 sXB :TDIAM(XB) XA+XB)/2. F0~F1.GT.0.) GO TO 90 ND.GT.99) GO TO 95 ABS(F1/FO).GT.5..OR.ABS(FOéF .GT. ALOG(ABS(F1)) GT 0. DR.ALD F 1) 5.) GO TO 40 ABS( O)).GT.O.) GO TO 40 000 50000:) O mmmwmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmoocomw mwmmmmmmmmmmmmmmmmm4444444444mmmmmmmmmmmmmmmmm 0040101RundomeImmbwn40mm4mmwa—Aomm4mmhww—Aomm4mmboo nunnunuunuuunuuuuuuunuuunnuunnnunnuuuuuununuuu 1.0 (D 000 _AOLO ll 0| CI 4 O) O O O U! 0 99 SECANT ITERATION DO 30 d= 1 MA IF(AB S(F15. GE. ABS(FO)) THEN CALL SWAP(FO.F CALL sw AP(XO x1 END IF X=X1-F1-(X1- xo)/(F1— -FO) LM=LEGVAR(X IF (LM NE. 0 GO TO 40 IF(X LT. XA. DR.x .XB) GO TO 40 DIFF=ABS(x-x1) IF(DIFF LE.ABS(x-ERROR)) GO TO 80 XO=X1 Fo=F1 x1=x FFx=FRICFAC(x) ND=ND+1 FFYéNO§=FFX DIx NO =x F1=DPTDIAM(x) CONTINUE BISECTIDN ITERATION FFX=FRICFAC(XM) NO=ND+1 FFYENO)=FFX DIx NO)=XM FM=OPTDIAM(XM) IF EFM.EQ.O.) GO TO 70 IF FAxFM.LE.o ) GO TO so XA=XM FA=FM Fo=FA XO=XA F1-FB X1=XB GO TO 60 XB=XM FB=FM X1=XB F1=FB XO=XA Fo=FA XM= EXA+XB)/2. iFXMABS(XA- X8). GT. ABS(XM*ERROR)) GO TO 10 Nwa=x FFX=FRICFAC(NEWX) N0=NO+1 FFYENO}=FFX DIx NO =NEWX RETURN NEWX=(XA+XB)/2. NOROOT=1 RET NEWX= (XA+XB)/2. NOR OOT =2 RETURN END IDLOKHDKMD ”NJ—544 400040) 1910mm MMMM 0!wa untnnnulnuuuu 000000000 9268C 944: 945: 100 FUNCTION OPTDIAM(DI) THIS FUNCTION EXECUTES THE OPTIMUM DIAMETER EQUATION. .(3- ~49). REWRITTEN AS OPTDIAM(DI)=O LIST OF VARIABLES: API= ANNUAL FIXED COST OF THE PIPE SYSTEM EXPRESS 0 AS A FRACTION OF THE INITIAL INSTALLED COST OF TH IPE SYSTEM.1/YR.EQ.(2-16) ANNUAL FIXED COST OF THE PUMP STATION EXPRES E FRACTION OF THE INITIAL COST OF THE PUMP STA I N. ANNUAL MAINTAINANCE COST OF THE PIPE SYSTEM EX RE AS A FRACTION OF THE INITIAL INSTALLED COST OF TH SYSTEM.1/YR ANNUAL MAINTAINANCE COST OF THE PUMP STATION EXPRES AS A FRACTION OF THE INITIAL INSTALLED COST OF THE P STATION. 1/YR CD= EMPIRICAL CONSTANT FOR THE PUMP STATION COST. $/W**PPP CEP= COST OF ELECTRICALM ENERGY. S/W HR CHEL= ELEVATION CHANGE. APP= BPI= E E P S D AS T O 1/YR.EQ .(2-16) P SSED E PIPE BPPB SED UMP CHPS= PRESSURE CHANGE. PA CI= EMPIRICAL CONSTANT FOR THE PUMP STATION COST. CP= EMPIRICAL CONSTANT FOR THE PIPE SYSTEM COST. $/M**(1+PPI) D = FLUID DENSITY. KG/Mfi'3 DIs TUBE/PIPE INSIDE DIAMETER.M EFF= COMBINED FRACTIONAL EFFICIENCY OF PUMP AND MOTOR FFx= FANNING FRICTION FACTOR HR= HOURS OF OPERATION PER YEAR K= CONSISTENCY COEFFICIENT PA s--N LEGT: TUBE/PIPE LENGTH MFR: MASS FLOW RATE KG/S OP1.DP2.OP3.DP4= WORKING CALCULATIONAL PARAMETERS PI: 3.141593 POWER= POWER REQUIREMENTS.EQ (2-13) PPI= EXPONENT IN THE PIPE SYSTEM COST EQUATION PPP= EXPONENT IN THE PUMP STATION COST EQUATION SUFFC= SUMMATION OF THE FITTINGS RESISTANCE COEFFICIENT WORK: WORK PER UNIT MASS.U/KG.EQ.(3-35 WORK1= WORK REQUIRED DUE TO PIPE FRICTION WORK2= WORK REQUIRED DUE TO FITTINGS FRICTION WORK3= WORK DUE TO PRESSURE AND ELEVATION DIFFERENCE IN THE SYSTEM YS= YIELD STRESS.PA LIST OF FUNCTIONS: DFFWD= DERIVATIVE OF FFX WITH RESPECT TO DI REAL MFR. PARAMETER COMMON/MF COMMON/PI COMMON/PU COMMON/FRIC CALCULATION OF WORK FROM EQ. (3-35) WORK1=32. I"LEGTRMFRI‘IMFR'"FFX/(PIWPII"DE"IDE*DI"HIE:) WORK2=8. *MFR*MFR*SUFFC/(PI*PI*OE*DE*OI**4) WORK3=CHPS/DE+9. B'CHEL WORKtWORK1+WORK2+WORK3 CALCULATION OF POWER FROM EQ.(2-13) POWEReMFRPWORK/EFF OPTIMUM DIAMETER EQUATION OP1=(API+BPI)IPPI*CPtPI-PI*OE*OE*EFF*DI**(PPI+5.) OP2=32. *CEPFHRPMFR**3 OP3=(APP+BPP)‘PPPPCOFPOWER**(PPP- 1. )/(CEP-HR) OP4-5. -FFx- DIFDFFWOEDI)+DI*SUFFC/LW OPTDIAM=1. -OP2*OP4- OP3+1 )/OP1 K I RC PE MP OHHmHm L p C C C F 994: l! 0 00000000000 J-AJ-b-L-b-L-h-L—L—fl-fi—A—A—L-Ld_‘J—A-b—L-L-b—‘J 000 OOOOOOOOOOOOOOOOOO00000000 wwwwwwwnwnnnmnnnna.‘444...“...4 4610'bww40004mmwa—POQOQ4QUIwa—n O 000-4» 0 101 FUNCTION OFFWD(DI) THIS FUNCTION COMPUTES THE DERIVA FACTOR WITH RESPECT TO THE TUBE/P IPE FOR LAMINAR OR TURBULENT FLOW LIST OF VARIABLES: DI= TUBE/PIPE INSIDE DIAMETER. DICg LAMINAR- TURBULENT TRANSITION VALUE OF DIM. F8 DERIVATIVE OF THE FRICTION FACTOR WITH RESPECT TO DIAMETER(DI) 1.DF2= WORKING VARIABLES TO CALCULATE OFF-LAMINAR = DIMENSIONLESS UNSHEARED PLUG RADIUS X: FANNING FRICTION FACTOR FH1‘ VALUE OF THE FRICTION FACTOR AT DI+H FH2= VALUE OF THE FRICTION FACTOR AT OI-H = SMALL POSITIVE NUMBER 1= -1 IF OFF 15 INDEFINITE; +1 IF OUT OF RANGE; o OTHERWISE N= FLOW BEHAVIOR INDEX SIGMA: PARAMETER IN THE DFFx/OO EQUATION FOR LAMINAR FLow.EQ.(3—47) T OF FUNCTIONS: 1= EQ.(3-45} CFAC= FRICTION FACTOR FUNCTION IV E OF THE FRICTION INSIDE DIAMETER LIs FFL FRI REAL N COMMON/FLWIDX COMMON/UNSPLG IF(DI.LT.DIC) DERIVATIVE FOR LAMINAR FLDw.EQ.(a-47) =FFL1(EO) FFXFSIGMA- EO- FFx- (3 *N 4. ) 1. .ESIGMA- E0) \\ NUMERICAL APPROXIMATION FOR TURBULENT FLOW H=o.1 IF(DI. LT. O. 01) THEN HsO QQO1 ELSE IF(DI. LT O 1) THEN H=O 001 ) THEN ELSEOIF(DI. LT. 1. EL8E1IF(DI. LT. 10. ) THEN END IF FH1=FRICFA ACIDI+HI FH2=FRICFA C DI-H UIUIUIUIU‘ wa—t 0000 01010101 mm4m n n n n 000 60= 618 .L-A—L-A-A—A‘A—L—b-b-A (DOCHDOCHDOCHDO 000000 IIIIIINIIIIIINNUIIIOIII 555555500000 0 000000004444 4mmbmnaomm4m N O 102 BACKWARD OR FORWARD DIFFERENCE ARE USED NEAR DIC FOR DI < DIC IFEFH1.GT.FFX.AND.FH2.GT.FFX) THEN IF ABS(FFX-FH1).LE.ABS(FFX-FH2)) TH FORWARD DIFFERENCE DFF=(FH1-FFX)/H ELSE BACKWARD DIFFERENCE DFF=(FFx-FH2)/H END IF ELSE QUADRATIIC APPROXIMATION.EQ.(3-48) DFF=1FH1~FH2)/(2.-H) END IF IzLEGVAR(DFF) EN IF THE DERIVATIVE IS INDEFINITE OR OUT OF RANGE. THEN IT IS NEGLECTED. THIS MAY ONLY HAPPEN NEAR DIC FOR DI < DIC WHERE THERE MAY BE AN ABRUPT INCREASE OF FFX OCHDOCJ 50000000000MMMMFJMMMMM-A-AAAA-AA-A-‘40000000000000000000000 01.004 0000000000000000000 ‘dd‘J‘JJ‘JJJA-bd—LJJ—LJ—AJ-A—A—A—b-A—L-AJ‘J-A—L—A—L—AJ—l—AA—A-‘JJ—A—L—A—L—AA—A.‘AJAAJJJJ‘JAAA-‘JJJJJJAJJ‘ "III!!!"llIlllllnullllnllllllIlllnnllllununllnllnnllllll"llllllllilfl 000000000000000b55555b~bb 000000000000000000000000000000000000000000000000000000000 #wnaomm4mmwaAOmm4mmbwn4 103 FUNCTION FRICFAC(DI) THIS FUNCTION CALCULATES THE FRICTION FACTOR ACCORDING TO THE SCHEME OF FIGURE 4. TORRANCE RELATIONSHIP. EQ.(A-1 . IS USED TO OBTAIN THE INITIAL GUESSES OF THE FRICTION FACTOR FOR TURBULENT FLOW. LIST OF VARIABLES: A1= LOWER BOUND GUESS OF EO OR EOC A2= UPPER BOUND GUESS OF ED OR EOC CONDIT= FLOW CONDITION I. E. LAMINAR.TURBULENT OR CRITICAL DE= FLUID DENSITY.KG/M'-3 DI= TUBE/PIPE INSIDE DIAMETER.M ED= DIMENSIONLESS UNSHEARED PLUG RADIUS EOC= LAMINAR-TURBULENT TRANSITION VALUE OF EO FC= LAMINAR-TURBULENT TRANSITION VALUE OF FF FF: FANNING FRICTION FACTOR FTORR= TORRANCE’ S FRICTION FACTOR FOR TURBULENT FLOW. EQ. (A-1) FFO= LOWER BOUND FOR FF OR FTORR. .(3- 15) WITH EO=1. O FF1= FFWENDBETDRRGUESS FOR THE CALCULATION OF THE TURBULENT FF2= FINAL LOWER BOUND GUESS FOR THE CALCULATION OF FTORR FF3= FINAL UPPER BOUND GUESS FOR THE CALCULATION OF FTORR FF4= LOWER BOUND GUESS FOR THE CALCULATION OF THE TURBULENT FF. EQ..(3-27) WITH R=RC FF4P= WORKING VARIABLE TO CALCULATE FF4 FF5.FF7= UPPER BOUND GUESS FOR THE CALCULATION OF TURBULENT FF FF6= LOWER BOUND GUESS FOR THE CALCULATION OF TURBULENT FF FF8= FINAL LOWER BOUND GUESS FOR THE CALCULATION OF THE TURBULENT FF FF9= FINAL UPPER BOUND GUESS FOR THE CALCULATION OF THE TURBULENT FF n§=HgENERALIZED HEDSTROM NUMBER.EQ.(3-19) K= CONSISTENCY COEFFICIENT.PA S**N LIMLow= COUNTER LIMUP= COUNTER MFR= MASS FLOW RATE KG/S N= FLOW BEHAVIOR INOE NOROOT= ROOT INDICATOR: D-YES; 1-NO; 2-TDO MANY ITERATION NOTIME= COUNTER PI= 3.141593 EQ.(3-21) e HERSCHEL-BULKLEY GENER YNOLDS NUMBER= RECFPSIC REC2,REC3= WORKING VARI REC E2,RE3= WORKING VARIABL TOLERANCE ERROR FOR EQ. TOLERANCE ERROR FOR THE OLERANCE ERROR FOR E8. 0. S I" D 3 H 2 D D _g C 1) m C I" m Z .4 .4 I) D Z M mprm-un 410m mrno TOLT= TOLERANCE ERROR FOR E TOLV= TOLERANCE ERROR FOR E U= MASS AVERAGE VELOCITY.M/ YS= YIELD STRESS.PA LIST OF SUBROUTINES: BISECT2= ROOT FINDING SUBROUTINE: BISEC BISNEWT= ROOT FINDING SUBROUTINE: BISEC LIST OF FUNCTIONS: ‘SECANT METHODS N-NEWTON METHODS DFUN1= DERIVATIVE OF FUN1 WITH RESPECT TO EOC DTORREN= DERIVATIVE OF TORREN WITH RESPECT TO FTORR DUNPGRA= DERIVATIVE OF USPGRA WITH RESPECT TO EO FFTM= EQ. (3- 31; REWRITTEN As FFTMEFF)=O. FUN1: EQ. (3—22 REWRITTEN AS FUN1 EOC)-Q. R: TURBULENT PARAMETER. EQ (3-27) TORREN= TORRANCE EQUATION. EQ.(A-1) REWRITTEN AS TORREN(FTORR)=O UNPGRA: EQ.(3-18) REWITTEN As UNPGRA(E0;=O Y: LAMINAR FLOW FUNCTION (PSI) EQ.(3-14 0(H1 (H30 "fllllllllllllfllllllflflIUMIOII'IIIOIIIIINNIII (H30 0 can»qounuunOJOunnqounnonm40unngonnaunmaounnNOHMAOMOACMom~unm nnIIIquIInunuunun11nuu11nlnIrnunlnunutuuu1lunuutuunununuun (WOCL‘O m (H30 0(H1 (H10 0 5OHDQOHJUOHDQQHJMKNONKHONKHONdh*44wh‘d~h*JCHDOCHDOCHDOCNDmdNDmahomunomaHDmflflbmflflnm~h44~b4Q~FJQ~HDQOHDQ I” O() AJJJ‘JJJJJ-AJA-A—AJ-A-A-A-A—b—L—LJ-L—A—DA—L—A—A—b—l..A—t—b-A—L-A—L—S-A-‘u-L-A-B-‘J—t—A—L—bJJ‘J—b—L—A—A-‘d-‘J-‘J—h—AJJJJAJAJA ANDHKHOMKHOMmHOMmNOMKMOMKHO”BHOMKMOMKHOMmHOMRHOMKHOMKHOMahAddhAJ-MAA-MAA—MAJAHMAJ.haaanLAAmMAA—h¢A—hA4 -AOUHDQOHflbQHOACMDQ\thQWMUACMomxflflmthOA O 104 REAL MFR.K.N.P PARAMETER(PI=3.141593) CHARACTER CONDIT-1O EXTERNAL FFTM EXTERNAL FUN1.DFUN1 EXTERNAL USPGRA.DUSPGRA EXTERNAL FFT2.DFFT2 EXTERNAL TORREN.OTORREN COMMON/MFRCCF/MFR.K/YSTDEN COMMON/CRICON/REC.FC.EOC/F COMMON/UNSPLG/EO COMMON/TOLER1/TOLV.TOLI/TO COMMON/ROOTNO/NOROOT/BLOCK COMMON/CODFLW/CONDIT COMMON/HEDSTR/HX CALCULATION OF U FROM EQ.(3-34) U=4.~MFR/(PI:DE-DIa-2) CALCULATION OF RE FROM EO.(3-17) RETF}N/(3.*N+1.))-*N R52: DI/2.)"N REasu--(2.-N) RE=8.*DE‘RE1*RE2*RE3/K CALCULATION OF HE FROM EO.(3-19) IF(YS.EQ.O.) THEN HE=O. I DE/FLWIDX/N/AVEVEL/U N/RE.HE /YS. LWCO LER2/TOLC/TOLER3/TOLL.TOLT 1/RC E SE HE=(DE/YS)‘DI'fi2'(YS/K)“(2./N) END IF IF(HE.EQ.O.) GO TO 10 CALCULATION OF EOC THROUGH ITERATION FROM EQ.(3-22) HXIHE A1=O. A2=.999999999 CALL BISNEWT1A1.A2.40.TOLC.FUN1,DFUN1,EOC) IF(NOROOT.EQ.O) GO TO 10 PRINT- PRINT' PRINT." THE DIMENSIONLESS UNSHEARED PLUG RADIU$.EOC.WAS NOT’ PRINT‘.I FOUND IN THE RANGE ’.A1.' <= EOC <= ’.A2 PRINT'.’ ENTER A NEW RANGE FOR EOC: 0. <= EOC < 1.0 ...... ’ READ‘.A1.A2 ‘ NOROOT=O GO TO 5 CONTINUE CALCULATION OF REC FROM EO.(3-21) REC1-33600.*SORT(1./27.)tN/(1.+3.*N)*-2 REC2'E2.+N)**((2.+N)/§1.+N)) REC3- 1.-EOC)**(1.+2. N) CALCULATION OF PSI-CRITICAL FROM EQ.(3-14) WITH EO=EOC PSIc-Y(EOC) CONTINUE CALCULATION OF REC RECPIREC1'REC2'PSIC*'(2./N)/RECS RECIRECP/PSIC CALCULATION OF FC FROM EQ.(3-16) WITH PSIIPSI-C & RE=REC FC=16./RECP ...AA bubb1>b oomqmmbww—AommxlmmuwnAowmqmmuww40mmqmmAwM—Aowmumwhwn llI'll"llllll"ll"llllllllll"II"llllllll"IIllll""ll"llllnllllllllI'll"llll'lnlllillllllllllllllllllllllllu MNMMMNMMNNMMMMMMMM 555 40 000 0 U" 000 000 M O U" 000M 0 (ADMFOMMMMMNMMMMMMNMNMMM 0000000 nonun mmmmmmmmwwmmmmmmmmm\Iqxlquq\Iqqqmmmmmmmmmmmmmmwwmmmm Ow waowmqmmuuua (00000 b O 105 IF(RE.NE.REC) GO TO 15 THE FLOW Is CRITICAL CONDIT=’CRITICAL’ FF=FC ED=EOC GO TO 60 CON IF(RE. GT. REC) GO TO 30 THE FLOW IS LAMINAR CONDIT=’LAMINAR’ ESTIMATION OF EO THROUGH ITERATION FROM EQ.(2-18) IF(HE.E0.0.) GO TO 25 A1=EOC A2=O. 999999 CALL BISNEWT(A1 A2. 50. TOLL USPGRA DUSPGRA ED) IF(NOROO EC. 0) GO PRINT- PRINTP S:§fi;-.’ THE DIMENSIONLESS UNSHEARED DLUG RADIUS.EO.WAS NDT’ ERINT'II FOUND IN THE RANGE ’,A1.’ <= EO <= ’.A2 PRINT-.’ ENTER A NEW RANGE FDR ED: 0. <= EO < 1.0 ........ ’ READ'.A1.A2 NOROOT=O GO TO 20 CONTINUE ESTIMATION OF FF FROM EQ.(3-14) 8 (3-16) FF=16/(RE-Y(EO)) GO TO 60 CONTINUE THE FLOW IS TURBULENT CONDIT=’TURBULENT’ ESTIMATION OF LOWER UPPER BOUND FOR FTORR. FFO CALCULATED FROM EQ. (3-15) WITH ED 0. HEN TH USED TO ALCULATE THE LOWER BOUND GUESS FOR vTHED FACTOR FOR WHICH THE CONDITION ED < .0 IS LI IS ALSO USED TO ESTIMATE THE FINAL TURBULENT FF FFO=2.-vs/(DE-Unu) FF1=FFO+1.E-5 FF2=FF1 FF3=1.O CALCULATION OF FTORR THROUGH ITERATION FROM EO. (A-1) THIS VALUE IS USED TO GET A RANGE FOR THE FF CALL BISNEwT(FF2 FF3é SO. TOLT TORREN DTORREN FTORR) IF(NOROO O0) GO 4O PRINT- PRINT: PRINTI'I THE INITIAL GUESS FOR THE FRICTION FACTOR (FTORR) ' PRINT " WAS NOT FOUND IN THE RANGE ’.FF2.’ <= FTORR <= ’.FF3 . PRINT-i’ ENTER A NEW RANGE FOR FTORR > ’.FFO.’ ...... , REA D-. FF2 F3 NOROOT 0 TO 35 CONTINUE JAJJA wwwwwwwwwwwwwww 0mmmmw0‘mwmmhbhhhbhhbbwwwmwwuwaMMMMMMNNNM-a-I-—-L II II II nu II II II II II 000000 000 IIIllllllIIIIIIIIIIIIIIIIIIIIII'II'II'IIII'II wwwwwuwwww ommqmmAwMAmeamwbwwa0mmqmmhwro—Aommqmwbwnaommqm HA} qqqqqqqqmmmmmmmmm \Immbmw—Aoomqmmbww—A 0" 000b- IIIIIIIIIIII"llIlllllllllIIIIVIIlll"IIIIIIIIMIIIIIIIIIII 106 CALCULATION OF RC FROM EQ.(3'27) 'ITH RE=REC 8 FF=FC RC=R(REC.FC) ESTIMATION OF A LOWER 8 UPPER BOUND FOR FF IN TURBULENT FLOW. F4 IS ESTIMATED FROM EQ.(3- 27) WITH R=RC. THIS WIL GIVES A FRICTION FACTOR GUESS FOR WHICH R > RC DR EO < . A CONDITION FOR TURBULENT FLOW N/(1-N)) -(FF4P;R C**N/RE)--(2 /(2 -N)) - 4..." w} _A N.GLO.7)THEN OR _A. N.GE.O.6) THEN T R _A E.O.5) THEN _A F N T F N.GLO.4)THEN F N _L- GE.O.3) THEN w. . . . - H- mH- mH- \IH' mH- un—u 50- 63A (Db N—ATIIITILDITIQImml‘HUTIT‘M (N.GE.O.2) THEN nmrmHmrflHmr“Hrrflflzmmflmmmmmmmmmmnflmmmmmmmmmmmmm OZHHmZHmMZHmmummnOZfimrmflrmnrmmrflnrfimrmmrmmmmfln ZDzwAUZmAUZmAz:wm—uoummqmmqmm\ImmqmmqqummqumbIS II umn llmll umu um" urnn "mu III’TIII "le II'U M4 1 2. IF IME=1 =FF1 =FF5 LOW=1 UP=1 FF4.GT.FF8) THEN =FF4 LOW=2 IF FF6.GT.FF8) THEN =FFG Low=3 IF FF7 LT.FF9) THEN =FF7 UP=2 IF TINUE ESTIMATION OF FF THROUGH ITERATION FROM EQ.(3-31) CALL BISECT2(FFSo FF9 20 TOLT FFTM FF) IF (NOROO .NDROOT E0 2) GO TO 55 IF(LIMLO TEOO1OAND LIMUP EQ.1) GO TO 50 95'= 98‘ 99: 02= 03= 107 IF THE ROOT IS NOT FOUND. IT MAKES SEVERAL ATTEMPS WITH A WIDER FF RANGE BEFORE IT ASKS THE USER TO ENTER A NEW RANGE N h- 1" II M an .5 4r- ".5 _A ~50! (NOTIME.EO.3.AND.LIMUP.EO.2) THEN "TH-‘71"?! “£11 (flbHflt _LUI‘T'I iNT*.’ THE FRICTION FACTOR (F.F.) WAS NOT FOUND IN THE RANGE’ INTF.' ’.FF8.’ <= F.F. k: ’.FF9 INT-.’ ENTER A NEW RANGE FOR F.F. > I.FFO.’ ...... ' 3 F9 ZDVVTJ'OTJ 0171111311102) EO=2.*YS/(FF‘DE‘U'U) CONTINUE FRICFAC'FF RETURN END nucuu 00000000 Abba-wwwwww (AIM—AOIDIDQOIUIIS u 0 his 0101 an 00 Abbh5hhb:hb:bbbbhbbb§b \I 15» II II 0 JJJ‘JJJ‘JJ‘J—L—L—A—L-fi—L—A-fi-A—L UIUIUIUIUIUIbb 0'!wa 401.0 000000000()0000 4..444444444444444444444444 #bbbbbbbbbbhbbbhbbbbb§bbbb (Dmxququl‘luqqqmmmmmmmmmmmmmm II II II II II II II II II II II II II II II II II II II II II II II II II II 0 108 FUNCTION FUN1(EC) THIS FUNCTION IS EO.(3-22) REWRITTEN AS FUN1(EC)=O. LIST OF VARIABLES: EC= DIMENSIONLESS UNSHEARED PLUG RADIUS AT THE LAMINAR-TURBULENT TRANSITION ( 9) 3-1 HX= GENERALIZED HEDSTROM NUMBER.EQ. N= FLOW BEHAVIOR INDEX P1.P2.P3= WORKING CALCULATIONAL PARAMETERS REAL N COMMON/FLWIDX/N/HEDSTR/Hx P1=E16800.*SORT(1 /27. g/N)-(2. +N )-~((2. +N)/(1. +N)) P2= EC/(1 -EC)*'(1. +N) --((2 -N )/N) P3=1./(1. -EC)‘*N FUN1=HX-P1*P2-P3 RETURN END FUNCTION DFUN1(EC) THIS FUNCTION IS THE DERIVA ATIVE OF FUN1EO(3-22). WITH RESPECT TO EC. I. E. OF UN1(EC) = D FUN1(EC)/ D EC. IT IS USED IN THE NEwTON S ITERATION METHOD IN SUBROUTINE BISNEWT. LIST OF VARIABLES: EC= DIMENSIONLESS UNSHEARED PLUG RADIUS AT THE LAMINAR-TURBULENT TRANSITION N= FLOW BEHAVIOR INDEX P1.P2.P3.P4,P5= WORKING CALCULATIONAL PARAMETERS REAL N COMMON/FLWIDX/N . P1='(16800. *SORT(1 /27. g/N;'(2.+N)-'((2.+N)/(1.+N)) P2= 2.-N)-ECr-((2. -2. *N /N /N P3= 1.-EC)**((2. +N)/N) P4=_2.+N)' EC~F((2. -N)/N)/N P5=(1. -Ec )-*((2. +2 2.!Ng/N) DFUN1=P1u(P2/P3+P4/P5 RETURN END mmqmmbwwaowmqmmbww 000000000000()000 03: O4: AA...L.A.A.;.A.A.A.A.A.L.A.L_A.L.A‘44.;4‘...‘ mmmwmwmmbbbbbbhbhhhbbhbbb# (D qmmbwwéommummbwwaomm 0 .54.;444444444444—54444.‘ UIUIUIUIU‘UIUIUIUIUIUIUIUIUIUIUIUIUIUIUIUI 0()00000000 THIS FUNCTION IS EO.(3-18) REWITTEN AS USPGRA(EX)=O. FUNCTION USPGRA(EX) LIST OF VARABLES: EX= DIMENSIONLESS UNSHEARED PLUG RADIUS HE= N= FLOW BEHAVIOR INDEX P1 .P2= 109 GENERALIZED HEDSTROM NUMBER.EO.(3-19) WORKING CALCULATIONAL PARAMETERS RE= GENERALIZED REYNOLDS NUMBER.EQ.(3-17) LIST OF FUNCTIONS: Y: LAMINAR FLOW FUNCTION THIS FUNCTION EXECUTES EQ. mDC‘D‘O‘UOD Znunka¢0n1 U4UNNIZ> r N MON/FLWCON/ REnsx-- (2. FUNCTION Y(EX) LIST OF VARIABLES: EX= N: P1,P2.P3.P4= DIMENSIONLESS UNSHEARED PLUG RADIUS FLOW BEHAVIOR INDEX REAL N COMMON/FLWIDX/N P1:( P2: 2. 1.-EX)-*2 -Ex)- (1. +3. aN)/(1. +2. -N) 'EX~(1. P3= EX-~2- ~(1. +3. ~N P4=(1 -Ex)-*(1. +N g/(1+N v=P4-(P1+P2+P3)-uN RETURN END (PSI).EO.(3-14) ‘FLWIDX/N -N)/N) (3-14). I .E. PSI=Y(EX) WORKING CALCULATIONAL PARAMETERS 000000000000000000000 4444444444.:AAA—54444444444444.5444 UIUIUIUIUIUImmmmmmmmmmmmmmmmmmmwmmmmmm OIUIUIUIUIUIUIUIUImmbbhb##«hbbbwwwwwwwwwwn 0 (3000000000000 munmmmmmmmmmmmmmmmmmwmmmmwm cooncnoooocnon\quqqqqqqmmmmmmmmm mmbwn—Aoomqmmth-somm 401015004 uTuuullnuuuuunounuunun11uu11nuluuu AJJJ-L—L-DA—LJA-‘JJJ—L-L-L—L-A—A—A—A-‘J—L 0 110 FUNCTION DUSPGRA(EX) THIS FUNCTION IS THE DERIVATIVE OF FUNCTION USPGRA. E0.(3*18). WITH RESPECT TO EX. I.E. DUSPGRA(EX ) D USPGRA(EX)/D EX. IT IS USED IN THE NEWTON’ S ITERATION METHOD IN SUBROUTINE BISNEWT. LIST OF VARABLES: EX= DIMENSIONLESS UNSHEARED PLUG RADIUS FFL1= SIGMA E0 (3 45) HE= GENERALIZED HEDSTROM NUMBER. EO. (3-19) N= FLOW BEHAVIOR INDEX P1.P2.P3.P4= WORKING CALCULATIONAL PARAMETERS RE= GENERALIZED REYNOLDS NUMBER.EO.(3-17) LIST OF FUNCTIONS: FFL1= EQ.(3-45) v= LAMINAR FLOW FUNCTION (PSI).EO.(3-14) L N MgN/N )WCON/RE HE/FLWIDX/N (N/(1. +3. *N))- gi‘HE: P2 P1 FFL1(EX; YSEX)‘*P1 p UR 1-RE Ex-*((2 -2 -N /N GRA=P3+P4 FUNCTION FFL1(EX) THIS FUNCTION EXECUTES EO. (3-45). I.E. SIGMA=FFL1(EX) LIST OF VARIABLES: EX= DIMENSIONLESS UNSHEARED PLUG RADIUS IEEA‘IR‘ZI’TMIESMSISESRP 45) IN LOW BEHAVIOR INDEX 2. P3.P4.P5.P6= WORKING CALCULATIONAL PARAMETERS 'U‘I‘I REA SOMMON/FLWIDX/N P2= 1. +2 P3=1.+N P4=1.-EX P5=PTFP3-P4-*2+2.*PZIP1FEX-P4+P1*P2‘P3'EXT'2 P6=P2~P3~P4-:3+2.-P1-P3-Ex-P4u-2+P1-P2-P4-Ex-*2 FFL1=P5/PB RETURN END ..n...n,..s ..n,..s ..g...A...A...n...n...n..as ._A.-.A...A>..lb ..A .41 O)01010)OiuncntfltntncntfltfltflcnLfltfltfl (If)()(3(3()()()()()(7(3(7(7(3(7 ()()(3 ()(3(3 .4A..A._A._A .A..45 .A .a..‘..as _A..Jl._AI_A.—‘ ounmannmounmounmannw (3 _A-b—‘J—b—Ad—b mannmounmow t:Is(0600000000000(0(JCOBJBJFJAJRJBJBJ k)-tCDUDOD~JOTU1t>u)hJ—bC)“300~J(DL”1>(O unlununluuttnuulunntuuucunlluu ()(3()(3(7()(3(7()(3(3(3(3(3(3(3(1 _L .5 .A .A .Ar—A .A .2... .5 .5 monnmounmaunma» -A O) 0‘ (3(3() 4 O) (0(D~d(b(nl>(ohJ-AC) (3()(5 II '1 II {I II 'l 'l Cl I. II II (flLflLfltflLflCflLflLfltfl (3 .A .2 _A..; .5 .2 .L .A 030101010301010) 111 FUNCTION TORREN(FX) THIS FUNCTION EXECUTES THE FRICTION FACTOR EQUATION OF TORRANCE FOR TURBULENTrFLOW.EO.(A-1). REWRITTEN AS TORREN(FX)=O. LIST OF VARIABLES: EO= DIMENSIONLESS UNSHEARED PLUG RADIUS WORKING VARIABLES T0 CALCULATE E0 ED1.E02.E03= FX= FANNING FRICTION FACTOR HE= GENERALIZED HEDSTROM NUMBER.EO.(3-19) N= FLOW BEHAVIOR INDEX P1.P2.P3= WORKING CALCULATIONAL PARAMETERS RE= GENERALIZED REYNOLDS NUMBER.EO.(3-17) REAL N COMMON/FLWIDX/N/FLWCON/RE.HE CALCULATION OF EO FROM EQ.(A-2I =(N/(1. +3. *N))*'(2. *N/(2.’N)I =16. I"(2. ‘HE)‘*(N/(2 NI) :Fx X*RE"(2 /(2 N) E01*E02/EOB TORRANCE EQUATION P1=O.45-2.75/N+1.97FALOG(1.-EO)/N P2=((1.+3.-N)/(4.:N)):-N P3=1.97nALOG(RE:P2-Fx-*(1.-N/2.))/N TORREN=P1+P3-1./SQRT(FX) RETURN END E01 E02 E03 EO= FUNCTION DTORREN(FXI THIS FUNCTION IS THE DERIVATIVE 0F FUNCTION TORREN. EQ. (A-1). WITH RESPECT TO FX. I. E. DTORREN(F X): DTORREN(FX)/D FX. IT IS USED IN THE NEWTON’ S ITERATION METHOD IN SUBROUTINE BISNEWT. LIST OF VARIABLES: E0= DIMENSIONLESS UNSHEARED PLUG RADIUS E01.E02.E03= WORKING VARIABLES TO CALCULATE E0 EX= FANNING FRICTION FACTOR HE= GENERALIZED HEDSTROM NUMBER.EO.(3-19) N= FLOW BEHAVIOR INDEX P1,P2.P3= WORKING CALCULATIONAL PARAMETERS RE= GENERALIZED REYNOLDS NUMBER.EO.(3-17I REAL N COMMON/FLWIDX/N/FLWCON/RE.HE CALCULATION OF E0 FROM EO.(A-2I EO1=(N,/(1. +3. I"f\'))"'”*(2. Ii'N/(2. -N)) EO2=16.'(2.‘HE)*’(N/(2. I) E03=FX*RE'*(2./(2.-NII EO=EO1 E02/EO THE DERIVATIVE OF TORREN WITH RESPECT TO FX P1=3.94*EO~ SORT(FX)+N-(1. O) P2=3.94-1.-N/2.)' -EO)-SQRT(Fx) P3=2. *N -EO)*FX£11. 5 DTORREN= P1+P2)/P3 RETURN END 0000000000000000000 000 000 000 4...4444.5444‘444444444444‘4‘.‘444444.544444 mmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmm mmmmwmmmwmmmmmmmmmmqqqqq\I‘lflQNlmmmmmmmmmm 0 \Jxlxlxldxldxlflxlxldxl‘dxlxlfl‘lflm J—t-‘AJ-PAAJOOOOOOOOOOCO mummuwwaommqmmuww—som nuununnunuuunuuunnuu 00000f300000 JJAJJAJJJ-‘J—A—‘J—b—‘A—L-‘J 0 112 FUNCTION FFTM(FT) THIS FUNCTION EXECUTES EO.(3-31) REwRITTEN AS FFTM(FT)=O. LIST OF VARIABLES: EO= DIMENSIONLESS UNSHEARED PLUG RADIUS FT: FANNING FRICTION FACTOR HE= GENERALIZED HEDSTROM NUMBER.EO.(3-19) g; FEOW BEHAVIOR INOEx P2.P3= wORKING CALCULATIONAL PARAMETERS RE= GENERALIZED REYNOLDS NUMBER.EO.(3-17) LIST OF FUNCTIONS: R: TURBULENT PARAMETER. EO.(3-27I ROMBIN= INTEGRAL IN EO.(3-31I REAL N COMMON/FLwIDX/N/FLwCON/RE.HE CALCULATION OF RP FROM EO.(3-27) PR=R(RE.FT) CALCULATION OF E0 FROM EO.(3-32) EO= 2. *HE/PRFt2)*-(N/(2. -N)) P2= -EO)*'((2. N)/N ) P3= N/(1. +3. :N)) CALCULATION OF EQ.(3—31) FFTM=11.-PR*'2*P2 P3 ROMBIN(EO PR)"(2. -N)/RE EEBURN N FUNCTION R(RX,FX) THIS FUNCTION EXECUTES E0 .(3- 27). I.E. ’ ’=R(RX FXI LIST OF VARIABLES: FX= FANNING FRICTION FACTOR N= FLOW BEHAVIOR INDEX RX= GENERALIZED REYNOLDS NUMBER P1.P2= WORKING CALCULATIONAL PARAMETERS REAL N COMMON/FLwIDx/N P1=ERX*(FX/16.)'*((2.-N)/2.))"(1./NI P2= 1.+3.*N)/N R=P1-P2 RETURN END uwuaowmummbwwaowmqmwwaAOmmqmmbww40m 00000000000000000000000000000000000 U" 190040101500»40mmqmmhwwaowmqmm A O co‘IqqqqqqqquammmmmmmmmmmmmmmmmmmhbISAIn:InsAhwwwwwwwwwwnnMMMMMMMM-s flflflflflddflflfl44444444444~IQ~1444QQQQQQQQQQQQQ\lxlxl\l\IQQQQQQQQQ\IQQQQQQQQQQQQQQQ 40 0| ummuwwao m an l0 0 mmmmmcocoooco l0 .5 ll 4....44444....‘4‘4‘44444444J44.54444444444‘444444.1....AA‘AJAAAAJAAJJ44‘A‘A4444A4444J4—AJ mmcmmmo £0 qmmbwn Om N 0| '1 II N '0 0| flxlflxlfl‘lxlxlfl 113 FUNCTION ROMBIN(EX.RYI THIS FUNCTION COMPUTES THE INTEGRAL IN E0 BY THE ROMBERG METHOD. THE INTEGRATION IS IN SEVERAL SECTIONS DEPENDING THE FLOW BEH INDEX AND THE INTEGRAL IS INDEPENDENTLY CO FOR EACH SECTION. THE VALUES OF THE DIMENS SHEAR RATE ARE SAVED IN ARRAY FS1 TO BE US LOWER & AND UPPER BOUND GUESSES IN CONSECU ROOT ITERATION OF FUNCTION FFT2(EO.(3-30). MILLER. A.R. 1982. "FORTRAN PROGRAMS". SY INC..BERKELEY.CA. LIST OF VARIABLES: DELTA = INTERVAL VALUE EX= DIMENSIONLESS UNSHEARED PLUG RADIUS FS1.FS2.FS3.FS4= DIMENSIONLESS SHEAR RATE ARRAYS FU1= VALUE OF FFT1 AT EX OR 0 FU2= VALUE OF FFT1 AT UPPER FU3= VALUE OF FFT1 AT X LOWER= LOWER LIMIT OF INTEGRATION N= FLOW BEHAVIOR INDEX PIECES= NUMBER OF INTERVALS RY= TURBULENCE PARAMETER.EO.(3-27) T= VALUES 0F ROMBERG TABLEAU TOLI= TLERANCE ERROR FOR THE INTEGRATION TOSUM= FINAL VALUE OF THE INTEGRATION IN EO.(3-31I UPPER= UPPER LIMIT OF INTEGRATION X= TRAPEZOIDAL POINTS LIST OF FUNCTIONS: FFT1= FUNCTION INSIDE THE INTEGRAL IN EO(3-31) ”-H'NHKDUA INTEGER PIECES.Nx(13) REAL N REAL LowER.T(92) REAL FS1(4097).F53(4097).F52(2049) COMMON/FLWIDX/N COMMON/TOLER1/TOLV.TOLI DO 5 KV=1,4097 FSTEKV§=O. F53 KV =0. CONTINUE DO 10 KV=1.2049 FS2(KV)=O. CONTINUE TOSUM=O. FU1=O. Fs1(1)=o. LOWER=EX 6F(§.LE.O.1) THEN 5L5; IF(N.LE.O.2) THEN 5L5; IF(N.LE 0.3) THEN ELSE V=.5 END IF UPPERev-(1.-LOVER)+LOVER CONTINUE FS1(2)=1 FU2=FFT1(UPPER.FS1(1),FS1(2).Ex.RY) FST(2)=FU2/UPPER~*2 FSA-FS1(2) PIECES=1 NX(1)=1 BEL;A=(UPPER-LOWERI/PIECES c=(FU1+FU2)/2 SUM=C T(1)=DELTA-C NM=1 NN=2 L1=1 1798825 1799= 3O 35 000 wwaommqmmhwaOmmqmmwaAOmm "ll (H15 0 (fits ”ll 0 oomcooooocooncomcnmmmmcooooooocooooooooooooooooooooocooocooooomonco h 01 Just:bwwmwwwwwwwnnmwwwwnnnaa..-4.A-4-..oo wnaommqm b a n1 0! O U'IUIUIUIUIUIU'IUIUIUIb 64: 66c JJJ-A—b—L-b-A-L-‘JAJ‘J—A—b-‘J-‘J—‘J—AJ—b-‘J-fi.L—‘J-‘J-‘a-‘d-‘A-‘J_fi-fidd-‘JJ-‘d‘dJ-fiJ-D-AJ-‘Jd mmmmmmmmmmmmmmmmmmmmmmmm 114 NM=NM+1 FOTOM=4. NX(NMI=NN PIECES=PIECEs-2 LL=PIECES-1 L2=(LL+1)/2 DELTA=(UPPER-LOWERI/PIECES COMPUTE TRAPEZOIDAL SUM FOR 2"(NM-1I+1 POINTS DO 30 II=L1.L2 I: II- 2- x=LowER+DELTA-I FU3=FFT1(x. FS1(II) FS1(II+1). Ex RY) SUM: SUM+F U3 FS2§III=FU3/X¥*2 F53 II 2)=FS2(II) CONTINUE NZ=NZ+L2 DO 35 KV=1L2+1 FS3(Kv-2 -1I=FST(KV) CONTINUE T(NN)=SUM-OELTA NTRA=NX(NM-1I KK=NM-1 COMPUTE T ARRAY DO 40 MM=1.KK U=NN+MM NT=NXENM-1)+MM-1 = FOTOM-T(U-1)-T(NT))/(FOTOM—1.) FOTOM=FOTOM-4. CONTINUE NEW ORDERED VALUES OF THE DIMENSIONLESS SHEAR RATE .L + .O. I GO TO 5 TENTRA+1I'T(NN:1II .L OLIII GO TO 60 T NN-1I'T(UII E. ABS 0 60 T. 12) GO TO 55 CONTINUE TOSUM=TOSUM+T(U) LOWEngPPER GO TO 75 9999999999) GO T0 70 I GO TO 65 .I 9 'fl‘flHHHTI‘TI 11 C M C) .4. 01.0 (010 E UPPER‘(%.'LOWERI/2.+LOWER CONTINUE TOSUM=TOSUM+(1. -UPPERI*(1.+FU2I/2. CONTINUE ROMBIN‘TOSUM RETURN END O000000000000OOOOOOOOOOOOOOOOOOOOOOOO ooooooooootmommwmmmmmmmmmmmmmmmqqqqqqquqqmm ommqmmwaJOmmqmmth.Aommqmmbwn-nomooqmmawn-sown: 0000 000 000 000 JAJJJJJJJAJ-‘J—A-fi-‘JJ—‘J-‘J—t—Acub-AA_b_A_A.AJul-hJub-A-AA-L-b-l-h—L—b-A-A—‘T—A-A-A‘J—A—‘J—b-A..A-L-b—L-L-A—A-AJJJJA 000 mmmmmwmwmmmmmmmomtomatomwmmmwmmIOmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmcomma) wwwwwwwwwnnnnmnMMMn-A‘a‘aa...»4.. mummwa-Aomooqmmbww40mmqmmhww4 O 115 FUNCTION FFT1(EY.C1.C2.EX,RY) THIS FUNCTION EXECUTES THE FUNCTION INSIDE INTEGRAL IN EQ.(3-31). I.E. FFT1=(2ETA)-( LIST OF VARIABLES: B: EMPIRICAL WALL EFFECT PARAMETER IN MIXING LENGTH THEORY.Eo.(3-2e) CV= DIMENSIONLESS RATE OF SHEAR (ZETA) C1: LOWER BOUND GUESS FOR CV C2= UPPER BOUND GUESS FOR CV X THE II‘F2 E= DIMENSIONLESS RADIAL COORDINATE (XI) EX= EIMENSIONLESS UNSHEARED PLUG RADIUS 8 Ez= EX HE= GENERALIZED HEDSTROM NUMBER EO.(3-19) L= DIMENSIONLESS MIXING LENGTH (LAMBDA). Q.(3-25I N= FLOW BEHAVIOR INDEX NOROOT= ROOT INDICATOR: O-YES; 1-N0; 2-TOO MANY ITERATION Q= PARAMETER IN MIXING LENGTH (PHI). EQ.(3-26) RC= LAMINAR-TURBULENT TRANSITION VALUE OF RY RE= GENERALIZED REYNOLDS NUMBER.EO.(3-17) 3;: EURBULANCE PARAMETER. EQ.(3-27) = Y TOLI= TOLERANCE ERROR FOR THE INTEGRAL OF EO.(3-31I TOLV= TOLERANCE ERROR FOR EQ.(3-30I LIST OF SUBROUTINES: BISNEWT= ROOT FINDING SUBROUTINE: BISECTION-NEWTON METHODS LIST OF FUNCTIONS: DFFT2= DERIVATIVE OF EFFT2 WITH RESPECT TO CV FFT2= EQ.(3'3OI REWRITTEN AS FFT2(CVI=O. HE EZ.R .L OTNO NOROOT \N CALCULATION OF B FROM EQ.(3-28) 8:22.-(1.+.00352rHE/(1.+.000504-HE)*-2)/N CALCULATION OF 0 (PHI) FROM EQ.(3-26) O=(RZ-RCI/(B*SORT(8.II CALCULATION OF L (LAMBDA) FROM EO.(3-25I L=.36-(1.-E)*(1.-EXP('Q‘(1.-EIII CALCULATION OF CV (zETA) THROUGH ITERATION FROM EQ.(3-3OI CALL BISNEWT(C1.C2.100.TOLV,FFT2.DFFT2,CV) NOROOTFO ' FUNCTION INSIDE THE INTEGRAL OF EO.(3-31I FT1=CV¢E~F2 ' F RETURN END 00000000000000 AAJJJJJAJJAJJJAAJJJAJAJ wwwmmmmmmmmmwmwwmmmwmmm mmmmmwmmwmmmaubbabhahbw aommqmmAwM—Aowmqmmawnaom égunuuuucunMunuuuuluuulluuulnn 00000000000000000 _LJ.A—A-A—A-A—L—A—A-A—L-AJJAAAJ—A‘J—L—L-LA mmmmmmmmmmmmmmmmmmmmmmmmww mmmmmmmmqqqqqqqqqqmmmmmmmm qmmAwMJOmmumwhwwaommqmmhww IINIIIIIIIIllllllfl"IlllllflllflllflfllllllIllIIfl 0 116 FUNCTION FFT2(CX) THIS FUNCTION IS EO.(3-3OI REWRITTEN AS FFT2(CXI=O. LIST OF VARIABLES: cx= DIMENSIONLESS RATE OF SHEAR (ZETA) E: DIMENSIONLESS RADIAL COORDINATE (XI) Ez= DIMENSIONLESS UNSHEARED PLUG RADIU L: DIMENSIONLESS MIXING LENGTH (LAMBDA). E0. (3 25) N: FLOW BEHAVIOR INDEX P1. P2= WORKING CALCULATIONAL PARAMETERS Rz= TURBULANCE PARAMETER. EO. (3-2 7) REA AL N CDMMON/FLWIDX/N/BLOCK2/E. EZ. RZ. L 333% 55% '55 93” ) I! 1: th2n 1 .E2 a:- 2 N . FFT2=P1+P2 ( ( / )/8 RETURN END FUNCTION DFFT2(CXI THIS FUNCTION IS THE DERIVATIVE OF FUNCTION FFT2(CXI. EO.(3T3OI WITH RESPECT TO CX. I.E. DFFT21CXI= D FFT2(CXI/D CX. IT IS USED IN THE NEWTON’S ITERATION METHOD IN SUBROUTINE BISNEWT. LIST OF VARIABLES: CX= DIMENSIONLESS RATE OF SHEAR (ZE TA IMENSIONLESS RADIAL COORDINATE (X IMENSIONLESS UNSHEARED PLUG RADI MENSIONLESS MIXING LENGTH (LAMBD ) I) D US I A). Bow BEHAVIOR INDEX T D g EO.(3-25I P WORKING CALCULATIONAL PARAMETERS URBULANCE PARAMETER. E0. (3' 27) DOZr'mm N-A II II N ll REAL N. L COMMON/FLWIDX/N/BLOCK2/E. EZ. RZ. L P1=N (1. -EZI-CXP*(NT P2= RZ**2*L-‘2 CX (1. -EZI"(2. /NI/4. DFFT2=P1+P2 RETURN END AAA—AJA—A (mmmmmom wmmommm awn-Foam: uunlluulu 000000 2040 117 SUBROUTINE BISNEWT(XA.XB.MAX.ERROR.FUN.DFUN.NEWX)I THIS SUBROUTINE COMPUTES THE ROOT OF THE FUNCTION FUN. IT IS A COMBINATION OF THE BISECTION AND NEWTON’ S ITERATION METHOD. THE MIDPOINT OF THE INTERVAL (PART OF THE BISECTION PROCESS) IS USED TO START THE NEWTON’ S ITERATION. THE PROGRAM CONTINUES WITH THIS METHOD UNTIL THE SOLUTION IS FOUND OR THE FOLLOWING SITUATIONS OCCUR; 1- DFUN‘O. 2- X FALLS OUTSIDE THE INTERVAL KNOWN TO CONTAIN THE SOLUTION: 3' THE DIFFERENCE IN SUCCESIVE APPROXIMATION DOES NOT DECRASED; 4- THE NUMBER OF ITERATIONS EXCEEDS MAX. IF THESE SITUATIONS OCCUR. THE PROGRAM SWITCH TO THE BISECTION METHOD TO OBTAIN A SMALLER INTERVAL. REFERENCE: MOORE. E. 1982. "INTRODUCTION TO FORTRAN AND ITS APPLICATION". ALLYN AND BACON. INC..BOSTON.MASS. LIST OF VARIABLES: DIFF= DIFFERENCE BETWEEN TWO ITERATION POINTS ERROR= TOLERANCE ERROR FA= VALUE OF FUN AT XA FB= VALUE OF FUN AT XB FMg VALUE OF FUN AT XM FUNC= VALUE OF FUN AT X FPRIME= VALUE OF DFUN AT X = -1 IF X IS INDEFINITE: +1 IF OUT OF RANGE: MAXIMUM NUMBER OF NEWTON’S ITERATION NEWX= ROOT OF FUN NDROOTF ROOT INDICATOR: 0“ YES: 1-NO OLDDIFF= DIFFERENCE OF PREVIOUS ITERATION OLDX= PREVIOUS VALUE OF X X¢ POINT FROM NEWTON’S ITERATION EQUATION XA= LOWER BOUND POINT USED IN THE BISECTION METHOD XB= UPPER BOUND POINT USED IN THE BISECTION METHOD XM= MIDPOINT BETWEEN XA 6 X8 LIST OF SUBROUTINES: SWAPs INTERCHANGE THE VALUE OF TWO VARIABLES LIST OF FUNCTIONS: FUN: FUNCTION WHOSE ROOT IS COMPUTED O OTHERWISE DFUN: DERIVATIVE OF FUN wITH RESPECT TO THE ROOT VARIABLE REAL NE COMMON/RDOTNO/NDROOT IF(XA. GT. XB)T EN CALL SVAP(XA XBI END IF FA-FUN(XA) IF(FA.EO.O .) GO TO 90 FBsFUN(XB) IFEFB.EO. o. ) OGO T0100 IF FA*FB.GT. .)GO TO 80 XM=(XA+XB)/2O OLDBIFxABS(xA-XB)/2. 2074: 2075=40 2076=50 20778 118 NEWTON’S ITERATION DD 20 U=1, MAX « OLDX=X . FPRIME= DFUN(X) M=LEGVAR(FPRIM IF(M.NE.O) GO TO 30 IF FPRIME EO. O.)GO TO 30 FUNC=FUN(X) x=X-FUNC/FPRIME M=LEGVAR(X) IF(M.NE.O) GO TO 30 IF (ABS(X-XM). GT. ABS(XA- XB)/2 ) GO TO 30 DIFF=ABS(X-OLD Dx) IF(DIFF. LE. ABS(X*ERROR)) GO TO 70 IF DIFF GE. OLDDIF) GO TO 30 OLDDIF=DIFF CONTINUE BISECTION ITERATION FM=FUN(XM) IFEFM.E0.0.IGO TO 60 IF FAtFM.LE.O.IGO TO 40 XA=XM FA=FM GO TO 50 XB=XM XM= (XA+XB)/2 §E&A85(XA- XB). GT. ABS(XM- ERROR)) GO TO 10 RETURN NEWX= X RETUR NEWX= (XA+XB)/2. NOR OOT RETURN NEWX=XA RETURN NEWX=XB RETURN END 4.444 . .AJJAJ 0155551:I:uhbbwwwwwuwwwwwwnnnwwnwmAAA44.-AA»; 0000000000000000000000000 _A II II II II II II II II II II II II II II I! II II II II II II II II II II II II II II II II II II II II II OI II II II II II OQDODQOIU'IbUM-POQOQQO’IUIAQMdomaidmwwa—tommflwlfiwa-AO 119 SUBROUTINE BISECT2(XA,XB.MAX,ERROR.FUN,NEWXI THIS SUBROUTINE COMPUTES THE ROOT OF THE FUNCTION FUN. IT IS COM MBINA TION OF THE BISECTION AND SECANT ITERATION METHOD. THE BISECTION INTERVAL IS USED TO STA AR T E SECANT ITERATION. THE PROGRSM CONTINUES WITH THIS METHOD 0 OCCUR: 1- X FQLLS OUTSIDE THE INTERVAL KNOWN TC CONTAIN THE SOLUTION 2- T OF RANGE OR INDEFINI E FAR AWAY M THE SO TION: - THE U ER OF ITERATIONS EXCEEDS X F E SITUATIONS O CUR. TH SWI CH TO THE BISECTION MET DO TO OBTAIN A SMA I ERVAL R ER NCE: MO R E "INTRODUCTION TO . E E O E 198 FORTRAN AND ITS APPLICATION“. ALLYN AND BACON. BOSTON MASS. LIST OF VARIABLES: DIFF= DIFFERENCE ERROSEN TWO ITERATION POINTS ERROR= TOLERANCE FA VALUE OF FUN 2T XA F LUE OF FUN AT XO F1 VALUE OF U T X L I INDEFINITE IF OUT OF RANGE; O OTHERWISE ITERATION M: -1 X 15 : + MAX: MAXIMUM NUMBER OF SECAN NE R F v mo... Z HMO m0 41> U TION O ECT ION METH DD XB= UPPER BOUND POINT XUSED IN THE BISECT ION METHOD WEEN 8 X8 XO= SECOND POINT REQUIRED FOR THE SECANT PROCESS X1= FOCUS POINT FOR THE SECANT ITERATION LIST OF SUBROUTINES: SWAP= INTERCHANGE THE VALUE OF TWO VARIABLES LIST OF FUNCTIONS: FUN= FUNCTION WHOSE ROOT IS COMPUTED REAL COMMON/ROOTNO/NOROOT (XA. GT. XBI THE CALL SWAP(XA. XBIN F1: 5 XM= (XA+XB)/2 Fo-F1 GT. 0. ) GO TO 90 /F OI .GT. 5T .OR. ABS(F O/F1I. GT I? TO 40 BS S(F 1)) -2. LOG(A BS S(FOI .-2.I GO TO 40 000 owmqmmwaAmeqmmwaJ bOOOQ O O m 0 m 0 m4 00 m MMMMMMMMMMMMMMMMMMMMNMMMMMMMMMMMMMMMMMMMMMMNMMMM O ‘4-AJ—L-AJ—A-b—b—fi—lA—A-A—b—A.‘_L—A—A_AJAAA—A...A—AJ‘JJJ—L—A—A—A—AJ—A—L-‘A—LAA mmmwmmmmwmmmmmmmmmmflu“flQflflQQflmmmmmmmmmmmmmmmmmmm mummbwwéowmqmmhwnaOmmqmmwa» flflflNII"llllllflflflll'lfllll'"“llllflflflflflflllII""ll"IIIIIIIINIINHIIIIIII! O 120 SECANT ITERATION DO so u=1 MAX IF(ABS(F1I.GE.ABS(FOII THEN CALL SWAP(FO.F1 CALL SWAP(XO,X1 END IF x=x1-F1n(X1-x0)/(F1-FO) LM=LEGVAR(X; LM.NE.O GO TO 40 LT.XA.0R.X.GT.XB) GO TO 40 ABS(X-X1I FF.LE.ABS(X*ERROR)I GO TO so OflXflXHUHH ZununAflA AflXMXOfixA BISECTION ITERATION XM) E0.0.I GO TO 70 F LE.O.I GO TO 50 GO TO 60 XM=2XA+XBI/2. iFXNABS(XA-XB).GT.A8$(XM*ERROR)) GO TO 10 NEWX=X RETURN NEWX=(XA+XBI/2. NOROOT=1 RETURN END 121 gégg=c SUBROUTINE SORTING(X.Y.NU) 2201=C SHELL-METZNER SORT FOR ARRAYS X 8 Y OF 2202=C INCREASING ORDER OF X. MAX. SIZE = 100. 2203=C MILLER. A.R. 1982. "FORTRAN PROGRAMS". 2204=C BERKELEY. CA. 2205=C 2206=C 220"= REAL X(100). Y(1OOI 2208: JUMP 2209=1O dUMPszMP/2 221o= IF(JUMP E0. 0) GO To 99 2211= d2=NU- dUM 2212: DO 30 d= 1. puz 2213= = 2214=2o U3=I+UUMP 2215= IF(X(I). LE. X(d3II so To 30 2216= CALL swApEx II ;.Xéd33g 2217: CALL SWAP .v 3) 2218= 1:1- 2219= IF(I. GT. P0) GO To 20 2220=3o CONTiNUE 2221: 60 T0 10 2222:99 RETURN 2223= END 2224=C §§§g=c SUBROUTINE SWAP(AA.BBI 2227=C THIS SUBROUTINE INTERCHANGE THE VALUE 2228=C 2229=C 2230= H0LD=AA 2231= AA=BB 2232= BB=HOLD 2233= RETURN 2234: END 2235=C SIZE NU. IN REFEREN E: C SYBEX N OF TWO VARIABLES REFERENCES 122 REFERENCES Barrett, 0.H. 1981. Installed cost of corrosion-resistance piping. Chem. Eng. 88(22):97-102. Boger, 0.0. and Tiu, C. 1974. Rheological properties of food products and their use in the design of flow system. Food Technol. in Australia 26:325-335. Boger, D.V. and Tiu, C. 1974. Rheology and its importance in food processing. Proceedings of 2nd National Chemical Engineering Conference. July 10-12. p. 449-460. Sunters Paradise, Queensland, Australia. Braca, R.M. and Happel, J. 1953. New cost data bring economic pipe sizing up to date. Chem. Eng. 60(1):180-187. Charm, S.E. 1978. "The fundamentals of food engineering." Third Edition. AVI Publishing Co., Westport, CT. Cheng, D.C.-H. 1970. A design procedure for pipe line flow of non- Newtonian dispersed systems. In "Proceedings of the First International Conference