Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

An analytical investigation of forced convective heat transfer to supercritical carbon dioxide flowing… Malhotra, Ashok 1977

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.

Item Metadata


831-UBC_1977_A1 M34.pdf [ 8.67MB ]
JSON: 831-1.0080805.json
JSON-LD: 831-1.0080805-ld.json
RDF/XML (Pretty): 831-1.0080805-rdf.xml
RDF/JSON: 831-1.0080805-rdf.json
Turtle: 831-1.0080805-turtle.txt
N-Triples: 831-1.0080805-rdf-ntriples.txt
Original Record: 831-1.0080805-source.json
Full Text

Full Text

AN ANALYTICAL INVESTIGATION OF FORCED CONVECTIVE HEAT TRANSFER TO SUPERCRITICAL CARBON DIOXIDE FLOWING IN A CIRCULAR DUCT Ashok MALHOTRA B.Tech., M.Tech., Indian I n s t i t u t e of Technology, New D e l h i , 1973 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY THE FACULTY OF GRADUATE STUDIES Department of Mechanical Engineering We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA October, 1977 ( 0 \ Ashok Malhotra, 1977 by in In presenting th i s thes i s in pa r t i a l fu l f i lment of the requirements f o r an advanced degree at the Univers i ty of B r i t i s h Columbia, I agree that the L ibrary sha l l make it f ree ly ava i lab le for reference and study. I fur ther agree that permission for extensive copying of th i s thes is for scho lar ly purposes may be granted by the Head of my Department o r by his representat ives. It is understood that copying o r pub l i ca t ion o f th i s thes i s for f i nanc i a l gain shal l not be allowed without my written permission. ASHOK MALHOTRA , MECHANICAL ENGINEERING artment of The Univers i ty of B r i t i s h Columbia 2075 Wesbrook Place Vancouver, Canada V6T 1W5 NOVEMBER 197 7 ABSTRACT A physical model and a numerical solution procedure has been developed to predict heat transfer behaviour in supercritical fluids. A major area of concentration was the modelling of the turbulent components of shear stress and heat flux. Traditionally, the turbulent fluxes are modelled by algebraic expressions such as the familiar mixing length methods. However, the use of this technique has not been entirely satisfactory. Newer methods for constant-property flows which model turbulent fluxes by considering the transport of quantities such as turbulent kinetic energy and the dissipation rate of turbulence have been extended to supercritical fluids. This involves the solution of two additional partial differential equations that are solved simultaneously with the equations of continuity, energy, and momentum. The numerical scheme has been developed on a com-pletely two-dimensional basis by extending the Pletcher-DuFort-Frankel f i n i t e difference method. Computed results for velocity and temperature profiles as well as wall temperature distributions exhibited reasonable agreement with previous experimental data and therefore indicate the v i a b i l i t y of the present method. Computations were carried out for supercritical carbon dioxide flowing through a circular duct in the reduced pressure range 1.0037 to 1.098. A consideration of the influence of buoyancy on the mean momentum balance permitted the calculation of unusual velocity profiles in this investigation. The existance of such velocity profiles had been accepted previously but the nature of their growth along a pipe has probably not been suggested previous to this work. No attempt was made to include buoyancy generated turbulence or additional fluctuating property correlations - i i -in this work, but suggestions are made regarding possible avenues of approach. Some of the incidental outcomes of this work were a new con-tinuous universal velocity profile implicit in cross stream distance an a new mixing length distribution for turbulent pipe flows. - i i i -TABLE OF CONTENTS Page ABSTRACT i i TABLE OF CONTENTS iv LIST OF TABLES v i i LIST OF FIGURES ...... v i i i LIST OF SYMBOLS x i I. INTRODUCTION . 1 1. .1 General 1 1.2 Review of Literature 8 1. i Need for the Present Work 16 II. Tllli CONSERVATION EQUATIONS 18 2.1 General 18 2.2 Statement of Problem 19 2.3 The Governing Equations 21 2.4 Algebraic Formulations for Turbulent Viscosity ............. 24 2.5 Boundary Conditions and Starting Profiles 26 2.6 Concluding Remarks 28 III. TURBULENCE MODEL 30 3.1 General 30 3.2 Reynolds Stress Equation 32 3.3 Turbulence Dissipation Equation 37 3.4 Algebraic Stress Modelling and Thin Shear Flow Simplifications 40 - iv -Page 3.5 Boundary Conditions for the Kinetic Energy-dissipation Model of Turbulence 45 3.6 Summary of Mathematical Formulation 53 3.7 Concluding Remarks 54 IV. THE NUMERICAL PROCEDURE .. 55 4.1 General 55 4.2 Required Capabilities of the Numerical Procedure 55 4.3 Finite Difference Formulation 57 4.4 Method of Solution 62 4.5 Preliminary Results 65 4.5.1 Laminar Flow Along a Flat PLate 65 4.5.2 Turbulent Flow Along a Flat Plate 67 4.5.3 Turbulent Flow in a Pipe 71 V. RESULTS AND DISCUSSION 79 5.2 Velocity and Temperature Profiles 80 5.3 Comparison of Wall Temperature Distributions 90 5.4 Scope for Further Results 101 5.5 Concluding Remarks 104 VI. FURTHER DISCUSSION 105 6.1 General 105 6.2 Turbulent Prandtl Number 105 6.3 Enhancement Factors 110 6.4 Free Convective Influences 115 VII. CONCLUSIONS AND RECOMMENDATIONS 121 - v -Page REFERENCES . . 123 APPENDIX A. Properties of Supercritical Fluids 130 APPENDIX B. Equations for the Turbulence Model 135 APPENDIX C. Computer Program 144 APPENDIX D. Listing of Computer Program 157 - v i -LIST OF TABLES Page TABLE 2.1 Value of F • 22 TABLE 3.1 Mixing Length Distribution in Turbulent Pipe Flow 49 TABLE 3.2 Summary of Mathematical Formulation 53 TABLE A.l Properties of Carbon Dioxide (Pressure = 1100.00 PSIA) ... 134 TABLE C l Description of Arrays Used in the Main Segment 150 TABLE C.2 Description of Input Variables ..... 152 TABLE C. 3 Description of Sub routines — 155 - v i i -LIST OF FIGURES Page FIGURE 1. Locus of Maximum C as a Function of Temperature for P Supercritical Carbon Dioxide . 3 FIGURE 2. Specific Heat (C p) of Carbon Dioxide as a Function of Temperature 4 FIGURE 3. Thermal Conductivity of Carbon Dioxide as a Function of FIGURE 4. Density of Carbon Dioxide as a Function of Temperature .. 6 FIGURE 5. Viscosity of Carbon Dioxide as a Function of Temperature . . 7 FIGURE 6. Schematic Diagram of the Problem .......... .... 20 FIGURE 7. Properties of Carbon Dioxide Near the Pseudocritical Temperature 41 FIGURE 8. Mixing Length Distribution i n Pipe Flow 50 FIGURE 9. Constant C ^ for the Dissipation Equation 52 FIGURE 10. The Finite Difference Grid 60 FIGURE 11. Blasius Velocity Distribution in the Boundary Layer Along a Flat Plate 68 FIGURE 12. Temperature Distribution on a Heated Flat Plate at Small Eckert Number for Pr = 2.57 ,69 - v i i i -Page FIGURE 13. Velocity Distribution i n a Turbulent Boundary Layer Along a Flat Plate . . 72 FIGURE 14. Growth of a Turbulent Boundary Layer and Decline of Shear Stress Along a Plate 73 FIGURE 15. Nusselt Number i n Pipe Flow 76 FIGURE 16. Turbulence Kinetic Energy i n Pipe Flow 76 FIGURE 17. Velocity Distribution i n Turbulent Pipe Flow ....... . 77 FIGURE 18. Temperature Profiles for Supercritical Carbon Dioxide — 82 FIGURE 19. Velocity and Temperature Profiles for Supercritical Carbon Dioxide 83 FIGURE 20. 'Unusual' Velocity Profile in Supercritical Carbon Dioxide 86 FIGURE 21. Development of 'Unusual' Velocity Profiles i n Super-c r i t i c a l Carbon Dioxide 87 FIGURE 22. Temperature and Density Profile in Supercritical Carbon Dioxide Corresponding to Data i n Figure 20...... 89 FIGURE 23. Wall Temperatures i n Supercritical Carbon Dioxide 93 I FIGURE 24. Wall Temperatures in Supercritical Carbon Dioxide ....... 93 FIGURE 25. Wall Temperatures i n Supercritical Carbon Dioxide 93 FIGURE 26. Wall Temperatures in Supercritical Carbon Dioxide ....... 95 - ix -Page FIGURE 27. Wall Temperatures in Supercritical Carbon Dioxide 97 FIGURE 28. Wall Temperatures in Supercritical Carbon Dioxide 99 FIGURE 29. Wall Temperatures in Supercritical Carbon Dioxide „. 100 FIGURE 30. Development of Unusual Velocity Profiles i n Supercritical Carbon Dioxide 102 FIGURE 31. Turbulence Kinetic Energy and Dissipation Profile Corresponding to Computations Illustrated i n Figure 19 ». 103 FIGURE 32. Turbulent Prandtl Number Variation Across a Pipe ........ 107 FIGURE 33. Wall Temperature Variations with Different Turbulent Prandtl Numbers '. 108 FIGURE 34. Enhancement Factors in Supercritical Carbon Dioxide ..... 113 FIGURE 35. Wall Shear Stress Variation 114 FIGURE 36. Wall Temperatures for Upflow and Downflow 117 FIGURE 37. Computed Radial Shear Stress Profiles for Carbon Dioxide at P/P = 1.027 at a Location x/D = 30 119 c FIGURE 38. Computed Values of Turbulence Kinetic Energy K for Carbon j Dioxide at P/P = 1.027 at a Location x/D = 15. (Data c Corresponding to Figure 36.) 119 FIGURE 39. Flow Chart of Computer Program 146 FIGURE 40. Nomenclature for Variables in Subroutines 156 - x -LIST 0F:SYMBOLS Constant in Van der Waal's equation Area Constant in Van Driest*s expression Constants i n starting profile for turbulence kinetic energy Constants in relationships determining f i n i t e difference grid spacing Constant in Van der Waal's equation Constants in pressure strain correlation Constants in the dissipation equation Constant in Reynolds' stress equation Constant in k - £ model Pipe diameter Eckert number for flow over that plate Representing function of quantity i n brackets Enhanced values of eddy d i f f u s i v i t i e s due to Hsu and Smith defined by Equation (1.7) Acceleration due to gravity Enthalpy A constant and a function respectively occurring in radial shear stress expression, defined by equations (1.11) and (1.12) Thermal conductivity Turbulence conductivity Turbulence kinetic energy Mixing Length Mass flow rate Nu , Nu Nusselt number as a function of streamwise location and for fully-developed flow respectively p Pressure p' Fluctuating component of pressure P. Production of turbulence kinetic energy Pr Prandtl number Pr Turbulence Prandtl number q Heat flux r Radial distance from axis of symmetry R Radius Re Reynolds number t Temperature t + Non-dimensionalised temperature t. Bulk temperature b t Outer wall temperature t Wall temperature in pipe w T Time u Velocity u' Fluctuating component of velocity u + Non-dimensionalised velocity = u/ / T IQ J W W U. Instantaneous velocity l i . . . . U Fluctuating Velocity component in i direction / \ Mean velocity component in i direction v Velocity in direction perpendicular to predominant flow direction v' Fluctuating component of velocity in direction perpendicular to predominant flow direction - x i i -x Strearawise distance measured from location where heating starts y Cross-stream distance measured from the wall y + Non-dimensionalised cross-stream distance = y A /p / u / p 3 w w/ w w T Shear stress T Turbulence shear stress y Viscosity y Turbulent viscosity p Density p ' Fluctuating component of density V Kinematic viscosity e Dissipation rate of turbulence kinetic energy e' Instantaneous value of £ Eddy d i f f u s i v i t y of momentum e, Eddy d i f f u s i v i t y of heat h 6 Boundary layer thickness 0' a' Constants occurring in k - £ model A Used in f i n i t e difference formulation as described in Figure 1Q T Quantity denoting diffusion of Reynolds stresses 3 Ratio of inner to outer diameter in a pipe Subscripts c Indicating c r i t i c a l or pseudocritical state w Indicating value of quantity at wall eff Effective value £,j,k Cartesian coordinates direction j , and k. £ corresponding to x, j to y. 1,2 Horizontal and vertical direction in two-dimensional flows t Turbulent - x i i i -ACKNOWLEDGEMENTS The author considers himself fortunate to have had an advisor in Dr. E.G. Hauptmann. His guidance and assistance came at every stage the work, without which this work would not have been possible. The author also wishes to express his gratitude to Dr. I.S. Gartshore for the interest and advice he rendered towards the present work. The University of British Columbia provided this wonderful opportunity to study at their institute as well as a much needed fellowship. Linda Blaine, who i s probably one of the best professionals available, did a l l the typing and much editing. She did this work, not as a professional but as an understanding sister would. Finally, the author wishes to thank his parents who have continuously provided encouragement from half-way across the world and others who are not recorded here but were nevertheless crucial towards the completion of this work. - xiv -1 I. INTRODUCTION 1.1 General Forced convection heat transfer to supercritical fluids has attracted the attention of many investigators during the past two decades. This interest has partly stemmed from the challenging nature of the problem and partly as a consequence of several engineering trends. There has been, a steady development of power plants towards supercritical conditions i n response to increasing demands for cheaper energy. Rocket motors are f r e -quently cooled by fuel at supercritical pressures. Supercritical fluids have been used as coolants for e l e c t r i c a l machines and have also been con-sidered as coolants for nuclear reactors. Very large changes of thermodynamic and physical properties near the c r i t i c a l state are responsible for the unique nature of supercritical  heat transfer. Effects which are normally approximated through constant  property idealizations become dominant near the c r i t i c a l state. Conse-quently, several commonly made assumptions and empirical correlations y based on constant property idealizations, become inapplicable in this region. At the present time, there i s insufficient knowledge of super-c r i t i c a l heat transfer mechanisms to produce completely satisfactory designs near the c r i t i c a l state. One of the goals of the present work is to improve this situation by an analytical examination of the phenomenon. The c r i t i c a l point i s often defined as the pressure and temperature at which no distinction can be made between the liquid and vapour phase of a fl u i d . At pressures higher than c r i t i c a l , the fl u i d i s termed supercritical. The term subcritical i s used for fluids below their c r i t i c a l pressure. At 2 supercritical pressures, heating of a f l u i d brings about a continuous change from a dense liquid - l i k e state to a lighter state resembling a gas (vapour- like) . No distinct interface occurs between a liquid and vapour at super-c r i t i c a l pressures. However, the pseudocritical temperature (the tempera-ture at which specific heat attains a local maxima) may be used as an approximate demarcating temperature between the liquid and vapour-like states. An i l l u s t r a t i o n of the preceeding terminology i s provided i n Figure 1. Property variations become less severe away from c r i t i c a l pressures and correspondingly the heat transfer mechanisms become simpler. The term supercritical heat transfer i s normally used only up to pressures where property variations are severe enough to make constant property heat transfer treatment inapplicable. The boundaries of this region are some- , what arbitrary. Below the c r i t i c a l pressure, boiling does not occur i f the f l u i d state does not encompass the saturation temperature. Therefore, single-phase heat transfer mechanisms apply. In a somewhat similar manner i f the' state of a supercritical f l u i d i s far removed from the pseudocritical tem-perature (so that property variations are not significant), conventional single-phase mechanisms apply. Abnormal heat transfer occurs only at or in the v i c i n i t y of the pseudocritical temperature. A brief discussion of the thermodynamic and transport properties of supercritical fluids i s provided in Appendix A. Figures 2 to 5 i l l u s t r a t e the variations of viscosity, density, specific heat '(C) and thermal con-ductivity of supercritical carbon dioxide at a pressure of 75.84 bars. 3 9.0 6.5 « a s a 304 305 306 307 TEMPERATURE °K 308 309 FIGURE 1. Locus of Maximum as a Function of Temperature for .'Supercritical Carbon Dioxide. 4 267.0 302.5 338.0 373.5 409.0 444.5 TEMPERATURE °K FIGURE 2. Specific Heat (C ) of Carbon Dioxide as a Function of P Temperature. I FIGURE 3. Thermal Conductivity of Carbon Dioxide as a Function of Temperature. 6 FIGURE 4. Density of Carbon Dioxide as a Function of Temperature. 7 FIGURE 5. Viscosity of Carbon Dioxide as a Function of Temperature. 8 1.2 Review of Literature Comprehensive surveys of near c r i t i c a l heat transfer with exten-sive bibliographies have been made by Hall"'", Hendricks, Simoneau* and. 2 3 Smith , and Hsu and Graham . The discussion in this section i s mainly restricted to forced convection aspects of supercritical heat transfer. Probably the earliest investigation of heat transfer to super-4 c r i t i c a l fluids was that of Schmidt, Eckert, and Grigull . They conducted their experiments in a natural convection loop and presented their results in terms of an apparent thermal conductivity, defined as the thermal con-ductivity required to transfer an equivalent amount of heat at the same temperature difference through a solid material of the same size. They found that apparent thermal conductivity becomes very large near the critical, point. A later study by Schmidt^ in a closed tube showed that apparent thermal conductivity at supercritical conditions can be as large as 10,000 times that of copper over a narrow range of temperatures. Con-trary to this and other studies indicating an improved behaviour at the c r i t i c a l point, several studies^ observed a sharp decrease i n the heat transfer coefficient near the c r i t i c a l point. This apparent contradiction was resolved when i t was demonstrated that maxima and minima in heat transfer can occur depending on the level of heat flux (Hauptmann^ and Styrikovich et a l . ). Differences between wall temperature and bulk temperature are 9 small in experiments where minima are observed. According to Hsu t the two results can be explained in boiling terms with instances of increased heat transfer coinciding with a mechanism similar to nucleate boiling, whereas a minimum heat transfer coefficient i s akin to film boiling. This 9 postulate was supported by the fact that maximum coefficients generally occurred at low temperature differences as in nucleate boiling. Some of the f i r s t direct observations in a supercritical f l u i d by G r i f f i t h and Sabersky^ indicated the existence of a bubble-like flow around a heating wire. Supercritical heat transfer i s sometimes accompanied by boiling-like 11 12 9 noise . Goldmann , therefore, was also led to conclude l i k e Hsu that under certain conditions of non-equilibrium, a boiling-like mechanism w i l l sometimes occur. Contrary to the preceeding evidence, flow visualisations carried out by Hauptmann^ and Green'" ^ have been in favour of single-phase forced convection mechanisms. The strongest argument against the boiling theories i s that two distinct phases are not possible under equilibrium at supercritical pressures and sufficient departure from con-15 dition:; of equilibrium has not yet been proved. Following Deissler , i t has been shown that conventional forced convection, theories that account for large property variations can account for most of the peculiarities including maxima and minima in supercritical fluids. This i s the approach adopted in most theoretical investigations and w i l l be pursued further i n the present work. Experimental investigations have mostly made use of smooth pipes of circular cross-section with uniform heat flux boundary conditions. The f i r s t detailed study was that of Wood"^, who made very careful measurements of velocity and temperature profiles. He also measured some unusual velocity profiles (which exhibited a maxima away from the center). Similar profiles have also been measured and studied carefully on a f l a t plate by Simoneau''"^ . A great deal of effort has been placed on developing 10 data correlations using either the Dittus Boelter type of equation or i t s modifications, d , t Nu x = C Re* Pr° (-f-) . (1.1) b Here, 'X' implies that properties are evaluated at a reference 18 temperature t . Suggestions for evaluating t v have been made by Eckert X X 19 and Bringer and Smith amongst others. Some other notable attempts to 20 21 develop correlations have been by Hess and Kunz , Krasnoschchekov et a l » 22 23 24 Gunson and Kellog , Swenson et a l . , and Mirapolsky and Shitsman Although these correlations are useful and one of the few means of predicting supercritical heat transfer, they have not been very success-3 f u l ; to quote Hsu and Graham : "It would be conceded by most investigators in this f i e l d that any of the correlations using the Dittus-Boelter type of equation are not f u l l y acceptable prediction schemes". Hall also concluded that existing correlations are inadequate as a means of predicting heat transfer i n the c r i t i c a l region. As the f l u i d state moves away from the c r i t i c a l point, property variations become less severe and then these correlations are more nearly applicable. A rather different approach to correlating and explaining forced convection heat transfer to supercritical fluids has been suggested by 25 Hendricks, Graham, Hsu, and Mediros , treating the problem as a pseudo-two-phase phenomenon. Another model i s the penetration model due to 26 Graham . The chief value in these latter approaches probably l i e s i n their heuristic description of mechanisms and not in their a b i l i t y to correlate data. 11 An alternative approach to predicting near-critical heat transfer i s application of turbulent forced convection theories with allowance for property variations. Accordingly, turbulent fluxes are defined as: < T(y) = (u + p O |H , (1.2) and q(7) = (k + P C p ^ ) |£ . (1.3) The principal unknowns in equation (1.2) and (1.3) are the flux variations T(y)» q(y)» and d i f f u s i v i t y e . The turbulent Prandtl number i s generally 15 27 regarded as equal to unity. Deissler ' pioneered attempts to analyse supercritical fluids and proposed that close to a wall d i f f u s i v i t y i s given by: £m 2 + + r, . 2 + + , . , ,„ — = n u y [1 - exp(-v n u y /v)] . (1.4) v w w 28 Goldmann redefined the non-dimensional parameter occurring i n the universal velocity profile i n order to make them more sensitive to property variations. 29 According to Tanaka , Goldman's proposal implies a redefinition of the same parameters occurring in the d i f f u s i v i t y expression which are, + ( f ) ' 1 / 2 dy + , (1.5) 0 w and 12 u J ^ p + (-£-) du . (1.6) Pw 0 30 Hsu and Smith modified Deissler's method by considering effect of density derivatives i n shear and heat transfer expressions as: e' - (1 + F)e ; e' = (1 + F )e : (1.7) m m h h m where, and, dy dy J n d£nC t + d£n£ 7 » (1.9) n , + , + dy dy The factors inside the parenthesis have been referred to as enhancement fac-31 tors. Another type of enhancement factor has been suggested by Yoshida 29 32 v Tanaka used Kato's expression for eddy d i f f u s i v i t y and Goldmann^s model for calculations. They were able to predict wall temperature peaks typical of deterioration regimes in supercritical heat transfer. Hess and 20 Kunz , on the other hand, adopted Van Driest's expression for eddy d i f f u s i -vity for calculations i n supercritical hydrogen and went on to develop a widely used correlation that involves an additional parameter v^/v^. Shear stress and heat flux variation have been assumed to be either constant or varying linearly in most investigations. This follows directly from the assumption that streamwise variations are not significant, 33 or equivalently that the flow i s one-dimensional. Shiralkar and G r i f f i t h found wall shear stress values by an additional mass flow constraint and were able to predict wall temperature trends in the deterioration regimes 13 34 for supercritical steam and carbon dioxide. Sastry and Schnurr , however, did not make the assumption of one-dimensionality except in the near wall region. They adopted the well known Patankar and Spalding scheme and used Van Driest's expression for d i f f u s i v i t y . Their computed results exhibited good agreement for supercritical steam i n the enhancement regime. In another group of studies i n which investigators having realized the importance of buoyancy even i n predominantly forced convection situations and have allowed the shear stress variations to be influenced by buoyancy. At least three attempts to make such calculations have been made and are 30 35 36 due to Hsu and Smith , Hall , and Tanaka Before undertaking the main part of investigation, a preliminary * one-dimensional study was carried out in order to gain familiarity with forced convection concepts in supercritical fluids. Calculations were made for pipe flow using a linear shear stress distribution as well as a shear 37 stress distribution that accounts for buoyancy in vertical flows , " »» -¥<^- i>+ ( 1- 1 0 ) where, I = PR pg(R-y)dy , (1.11) 0 The terms one-dimensional and two-dimensional in this investigation make reference to the minimum number of coordinates required to analyse the problem. 14 and, I y - pg(R-y)dy . (1.12) 0 Equations (1,10) to (1.12) follow directly by integrating the one-dimensional momentum equation. The remaining formulation i s essentially similar to that 33 of Shiralkar and G r i f f i t h . An important innovation i n the present inves-tigation was the use of temperature as an independent cross-stream variable in the numerical solution. This was made possible by rearranging equation (1.3) for 4^ rather than 4^ ". Since properties were tabulated as the function dt dy of temperature, this eliminated numerous interpolations for the properties which otherwise consume the bulk of computing time. An additional advantage was that inaccuracies associated with interpolaration were also reduced. It , also becomes possible to divide the grid in equal increments of enthalpy which yields an efficient grid distribution with suitably close points near the wall and in the v i c i n i t y of the c r i t i c a l point where properties change rapidly. Although i t may appear from the one-dimensional equations that computation time is not a severe limitation, this i s not the case due to time-consuming iterations for x and I. As a consequence of this study, w a need for an investigation at the two-dimensional level was reinforced. So" far, a l l analyses discussed have been quantitatively success-ful only in selected regions. Their main contribution, however, l i e s in the fact that they have predicted a principal peculiarity of supercritical fluids which i s 'enhancement' and 'deterioration' in heat transfer at low and high heat fluxes respectively. A one-dimensional analysis also succeeds in predicting some of the unusual velocity profiles. A more complete range 15 of velocity profiles has been calculated in this investigation (as discussed in Chapter V). In conclusion, calculation methods in supercritical fluids f a l l within two broad categories. The f i r s t of these i s through the use of empirical correlations and the second through the use of forced convection theories that make use of the mixing length type of algebraic formulations for d i f f u s i v i t y . The l a t t e r approach was adopted in this work in order to explore the possiblities that a more rigorous effort i n this direction can yield. However, an important change in this investigation was the replace-ment of the algebraic formulations for d i f f u s i v i t y by some of the more recent methods used for studying constant property turbulence flows. These methods formally referred to as turbulent models, d i f f e r according to the level of approximation and may involve one, two, or more additional d i f f e r e n t i a l equations for turbulence properties such as the Reynolds stress that are solved simultaneously with the three primary conservation equations (of mass, momentum, and energy). Therefore, this l a t t e r form of turbulence * modelling i s often referred to as multi-equation models of turbulence . Though s t i l l in their early stages of development, these models are regarded as an improvement over algebraic formulations (of the mixing length types) because of their success in application to many engineering situations. By accounting for the transport of turbulence properties, they place the art of turbulence modelling on a physically r e a l i s t i c l e v e l . See Launder and Spalding for multi-equation models of turbulence. 16 1.3 Need for the Present Work The designer of heat transfer equipment requires information on wall temperature distribution (or heat transfer coefficients) and wall shear stresses (or f r i c t i o n factors). It i s desirable that the prediction procedure available to him should be versatile enough to cover a wide variety of operating conditions, within which he may explore and optimise a design. Various attempts for supercritical fluids have been outlined in the l i t e r a -ture survey but as indicated, they are i n need of further development. However, they have pointed to a possible route which needs to be inve s t i -gated. The purpose of this dissertation i s to develop a numerical pro-cedure which can handle forced convection heat transfer calculation at super-critical conditions. The major areas of improvements i n the present work were the inclusion of a multi-equation model of turbulence and the extension of the analysis to a f u l l y two-dimensional lev e l . The procedmre i s general enough to be used (in most cases with only minor modifications) for any f l u i d flowing i n plane or axisymmetric configurations, but calcula-tions were restricted to flow of supercritical carbon dioxide through ci r c u l a r ducts. ** Carbon dioxide has a convenient c r i t i c a l range . Because of this and also because of i t s easy handling and av a i l a b i l i t y i n a f a i r l y pure form, carbon dioxide has been selected as a working f l u i d i n many experimental _ Heat transfer coefficients in supercritical fluids are heat flux dependent and therefore they do not have the same significance as i n constant property fluids. 'p = 73.8 bars, t = 304°K. c c 17 investigations. Its properties in the c r i t i c a l region are f a i r l y well known, * unlike most other fluids . Because of this and also because a variety of experimental data is available, carbon dioxide was selected as the f l u i d for testing the numerical procedure in this investigation. The problem investigated i n this study i s also interesting from an academic point of view in that the combination of d i f f i c u l t phenomena such as turbulent flow and large property variations are generalisations of phenomena that occur in most convective heat transfer problems. There i s also a need to inquire into the applicability of recent methods i n tur-bulence modelling and numerical analysis under the influence of such condi-tions ., The contributions of this investigation are twofold. F i r s t , the formulation of a physical model and second, the development of a numerical solution procedure. The three primary equations for conservation of mass momentum and energy are discussed i n Chapter II. These are f u l l y two-dimensional and include terms for buoyancy and viscous dissipation. The associated turbulence model discussed i n Chapter III i s a two-equation model of turbulence. Chapter IV describes the main features of the numerical solution procedure along with some preliminary results. The associated computer program i s described i n Appendix D. The remaining chapters dis-cuss the achievements and shortcoming of the prediction procedure. The properties of supercritical steam are also well known. 18 II. THE CONSERVATION EQUATIONS 2.1 General The main thrust of this investigation was concerned with flow i n circular ducts, but an attempt was made to keep the mathematical formulation suff i c i e n t l y general so that i t would apply to both plane and axisymmetric flow situations. This investigation was primarily concerned with time-averaged effects. The process of time averaging causes s t a t i s t i c a l corre-lations involving fluctuating velocities and temperature to appear i n the conservation equations. These unknown correlations are approximated from additional equations which generally take the forms of algebraic formulae (as in the mixing length models); they can also be determined from the solution of d i f f e r e n t i a l equations for time-averaged properties of turbu-lence. The latter approach i s relatively recent and has mostly been con-cerned with constant property fluids. Therefore, i t required a review before i t was introduced in this investigation. The complete mathematical model formulated in this work involves di f f e r e n t i a l equations for the fluctuating velocity correlations as well as the three primary conservation equations of mass, momentum, and energy. Only the latte r set of equations, along with algebraic formulae for turbu-lent viscosity, are discussed i n this chapter. The remaining equations are taken up i n the next chapter. This division has been made for the sake of c l a r i t y . 19 2.2 Statement of Problem The problem under consideration was forced convective heat trans-fer to supercritical fluids flowing past r i g i d surfaces in two-dimensional plane or axisymmetric configurations. Heat addition or removal took place only at the r i g i d boundaries. The flow was normally turbulent and the mean, flow was assumed to be steady (hydrodynamically as well as thermally). The f l u i d was assumed to be at supercritical pressures and the f l u i d properties were assumed to be known. Figure 6 shows a schematic diagram of the problem as applicable to pipe flow. The conservation equations considered in this chapter include a. term accounting for buoyancy. The necessity of including this term may be questioned since the primary interest in the present work is forced con-vection heat transfer. However, i t has been indicated by Hauptmann'7, 30 1 Hsu and Smith , and Hall that free convection effects can influence several predominantly forced convection situations. It i s d i f f i c u l t to define precisely when these effects are totally absent. Moreover, i t i s f e l t that the role of buoyancy in the mean momentum balance for forced convection situations i s not f u l l y understood. A consideration of this influence i n forced convection heat transfer provides an opportunity to con-tribute towards a further understanding of this situation. Nevertheless the mathematical model formulated i n this work is not generally applicable for both forced and free convection heat transfer since the turbulence model used in this work i s primarily meant for forced convection heat transfer. It should be pointed out that in previous analytical investiga-tions (treated as one-dimensional) the free convective influence has been conceived only i n terms of a buoyant body force contribution in the mean HEAT FLUX q w 21 momentum balance. 2.3 The Governing Equations The equations of continuity, energy and momentum governing the flow can be written in such a way that they apply to plane as well as to rotationally symmetric flow. For the problem under consideration, these * equations may be written as : 9^ (Pur ) + ^  (pvr ) - 0 (2.1) 9u , 9u P u 9x" + p v 97 dx a 9y a , 9u iV^ - pu'v') + F x (2.2) and. p u 9^ + p v -97 dp , 1 9 dx a 9y ra flc — - ph'v») + (u - 9u/9y' ^ 9y; (2.3) where: a = 0 (for plane flows), a = 1 (for axisymmetric flow) The body force term in the x-momehtum equation, F , may take the values indicated i n Table 1.1. * 39 See Eckert and Drake 22 TABLE 2.1 Value of F x FLOW CONFIGURATION F x I Non-buoyant or Horizontal Duct Flow 0 | II Upward Duct Flow ~Pg I III Downward Duct Flow +Pg 1 IV Upward Facing Flat Plate + TM pgdy Jy V Downward Facing Flat Plate For a horizontal pipe, the indicated body force i s zero. This i s actually only true for non-buoyant flows. The case of horizontal-pipe buoyant flow cannot be treated here since in that case the axisymmetric condition i s violated. The body forces for the f l a t plate IV and V of Table 1.1 are actually not the body forces; these terms arise due to body forces + pg in the y-momentum equation. However, since the y-momentum equation i s not stated spearately, these values can be seen as F^ with dp/dx taking a value equal to zero. The contribution of fluctuating property correlations has been neglected in previous investigations since i t i s f e l t that their overall influence i s not significant. As a preliminary step, this practice was continued in the present investigation in the preceeding conservation equa-tions. Hall''" has estimated these correlations and carried out a useful order of magnitude analysis. His conclusion was that these correlations are probably negligible, but further discussion on this issue i s provided in Chapter VI. 23 Since the pressure variations i n the cross-stream direction were not significant i n this investigation, then: , 3t k 3h m i \ P Whenever indicated by the turbulence model under use that the turbulent fluxes could be expressed as a product of a function (which was either known or could be determined) and the velocity gradient, then this function was called the turbulent viscosity (or conductivity). This i s the sense i n which the effective viscosity concept has been used i n this investigation: r - K £ . (2.5) pu'v' -t 3y K t 3h ph'v' C 3y * P (2.6) and in addition, u f - y + y t , (2.7) k c c = k + k . (2.8) ef f t The ratios of the two turbulent coefficients i s then the turbulent Prandtl number (Pr f c = C p y t / k t ) . The pressure gradient in the momentum equation for non-confined flows w i l l be known and for the confined flow, i t can be found by the additional constraint of known mass flow rate (see section 4.3 for details) (2.9) 24 2.4 Algebraic Formulations for Turbulent Viscosity No precise way i s yet available for completing the preceeding equation set. The conventional means however is through the use of alge-braic formulations for the turbulent viscosity. Predominant among these is Prandtl's mixing-length hypothesis. This type of formulation has been used in the past for an analytical treatment of supercritical f l u i d s . However, at the present time i t i s possible to make use of apparently improved models of turbulence and the intention in this investigation was to explore an alternative model. Nevertheless, the algebraic formulations were s t i l l necessary i n the laminar and transitional region close to the wall. In addition, they were useful for developing a solution procedure and for deriving certain boundary conditions and starting profiles. Some of the algebraic expressions that have found use in analy-t i c a l work dealing with supercritical fluids are due to Deissler, Kato, and Van Driest. These expressions were originally formulated for con-stant property fluids but are used with the assumption that i f the f l u i d properties appearing in the expressions are allowed to take their local values, they are s t i l l applicable. In addition, a l l the available 29 expressions may be redefined by the use of Goldmann's hypothesis However, the general experience i n this and other investigations ( a l l one-dimensional analysis) has been that whichever expression is used, the results are not altered significantly. The most widely used formulation i s due to Van Driest. It was also selected in the present work, u t = pl If • <2'10> I = Ky [1 - exp(-y +/A +)] , (2.11) where, K — 0.4 and A + = 26 . For y + > 26, i t i s possible to use % = Ky with comparable accuracy. However, this was not done in order to ensure a smooth merger of the transitional region with the turbulent region. Beyond y = 120, Van Driest's expression reduces to £ = ky within 1%. At least three different proposals exist for expressing the argument of the exponent in equation (2.11). None of these proposals i s advantageous under a l l circumstances. However, the procedure of using wall properties to evaluate the damping function has been shown to be superior to using local properties for fluids of moderate Prandtl number . Accordingly the expression chosen was: (2.12) Another expression for the mixing length £, away from the wall, used by 40 41 Pletcher and Nelson ' , was used in this investigation for some preliminary testing and i s for I > 0.0896 replace by I = .0898 , (2.13) where 6 i s the thickness of a developing boundary layer. Even for a f u l l y developed boundary layer this mixing length distribution has given good agreement for the velocity profiles and reasonable agreement for the 41 f r i c t i o n factors (about 5% higher). Once the turbulent viscosity is known, the turbulent conductivity is normally found by the assumption of unity turbulent Prandtl number. This _ — _ See Pletcher and Nelson and McEligot et a l 26 value has been advocated in previous investigations dealing with super-c r i t i c a l fluids. This assumption is probably a good one when the Prandtl number is close to unity. A value of 0.9 instead of unity has become popular i n recent years i n numerical studies of convective heat transfer. Accordingly, solutions were made with this value. This quantity i s d i s -cussed in greater detail i n Chapter VI. 2.5 Boundary Conditions and Starting Profiles The equations of continuity, energy, and momentum require starting profiles and boundary conditions before a solution can be effected. Starting enthalpy profiles were constant i n i:he problems considered. The velocity profiles, however, were not normally constant and these had to be generated by other means. In turbulent flow problems however, i t is not necessary to know these precisely since downstream results are not very sensitive to the starting profiles. It was nevertheless necessary to ensure that these profiles were smooth and continuous i n order to avoid i n s t a b i l i t i e s i n the numerical procedure. In a widely used numerical procedure , the starting velocity profiles normally follow from the 1/7 power law. This i s appropriate in numerical procedures where a one-dimensional analysis i s used in the near wall region. In the present solution procedure, a starting pr o f i l e was also required i n the near wall region. In principle i t should be possible to use a laminar law profile in the laminar sublayer and part of the transitional region (up to y + = 11),, a log law in the rest of the transitional region and most of the turbulent region with perhaps a 1/7*"*1 power law prof i l e near the edge of the boundary layer. Such a procedure w i l l however inevitably introduce discontinuities. 44 A continuous prof i l e due to Spalding is also available. However, i t s use 27 was not convenient for this work because i t appears to be implicit i n velocity. Continuous profiles e x p l i c i t l y defining u were devised for the * present work but the method of generating the i n i t i a l profiles by standard solution procedures has been preferred because the latter procedure i s more consistent with equation ( 2 . 2 ) . In a f u l l y developed constant property pipe flow, the shear stress varies as x = T r/R (2.14) w Using these values of shear stresses and noting that the Van Driest hypo-thesis can be translated to yield / [ 2 K 2py 2 {1 - exp(-y +/A +) 2}] . (2.15) The starting velocity profiles for pipe flow problems where the flow i s f u l l y developed before entering the solution regime can then be easily generated. In order to obtain the starting profile for f l a t plate, problems, the shear stress T i s assumed to be constant at the starting location and equation (2.15) i s solved for u up to a point where i t s til solution matches with a 1/7 power law pr o f i l e . In addition to starting profiles, the enthalpy and momentum equations require one boundary condition at each of i t s boundaries. These w i l l either be of known enthalpy and velocity or a known gradient. Thus, at the wall where y = 0, du 3y 2 2 2 -u + (u + 4TK py {1 - exp(-y +/A +) 2}) 1/2] *e.g., U + = B[l - exp(-y+/B - (y +) 1' 5/100)]; B = 2.5 *n(y + + 1) + 5.5. 28 u = 0; v = 0; t = 3 t t or -5— W dy y=0 w At a free stream boundary, u = uro; t = t Q and at an axis of symmetry, 37 =0 ; 3y" =0 • 2.6 Concluding Remarks The present model differs from a previous two-dimensional analysis" of supercritical fluids, f i r s t by the fact that additional terms have been added to account for viscous dissipation and buoyancy and second, the governing equations were used in their f u l l form a l l the Way to the wall. A one-dimensional analysis at the wall which i s an important feature of the 43 34 45 numerical procedure adopted by Sastry * , i s made possible by neglecting convection in the laminar and transitional near wall region. This assump-tion has been widely tested and proven to be good for constant or moderately variable property fluids. However, for supercritical fluids when the pseudo-critical temperature l i e s in this region i t appears that rapid changes i n density w i l l make this assumption unsuitable. What effect this has on the overall flow computation has not yet been sufficiently evaluated. One d i f f i c u l t y with using a one-dimensional analysis for supercritical fluids i s that their wall shear stress behaviour i s not known adequately either through standard correlations or experimental data. In the absence of this information, wall shear stresses in one-dimensional 29 methods can be found only by use of iterative schemes that do not always converge easily. With a two-dimensional analysis, this problem does not arise since in that case the wall shear stress can be found in a straight-forward manner from the computed velocity profiles. A higher order tur-bulence model that was incorporated in the analysis i s discussed in the next chapter. 30 III. TURBULENCE MODEL 3.1 General A major limitation i n prediction of turbulent flows i s an impre-cise knowledge of velocity and enthalpy correlations p U^ U^ . and p hU^ respectively. Several different approaches to modelling these correlations may be found in the literature. So far, the various different forms of algebraic formulations have found the widest use. However, in recent years, a different class of turbulence models has gained in popularity. These involve a simultaneous solution of d i f f e r e n t i a l equations for quantities such as the kinetic energy of turbulence, dissipation, and Reynolds stresses . Although a significant improvement over the algebraic formulations, the alternative models of turbulence considered here have their limitations as well. It should be pointed out that turbulence modelling i s not a rigorous theoretical effort. Exact equations are derived as a starting point and guide to modelling, but only those terms that can be exactly defined are retained. The remaining information i s supplied through approximation and empiricism. However, i t is f e l t that further improve-ments in turbulence modelling w i l l come from this direction considering the fact that these models are presently undergoing a rapid development. As already mentioned, only the high Reynolds number form of modelling w i l l be considered. That i s , the turbulence modelling w i l l be applicable only i n the f u l l y turbulent region of the flow. For laminar and transi-See for example, Launder and Spalding , Launder, Reece, and Rodi , and 47 Hanjalic for multi-equation approach to turbulence modelling. 31 tional regions, the present study makes use of the Van Driest hypothesis. This was f e l t j u s t i f i a b l e since modelling in this region has as yet made l i t t l e progress even for constant property fluids. As further progress takes place i n low Reynolds number regions, i t can be added to the methods proposed in this investigation since the remaining equations are i n any case solved in the laminar and transitional regions. Of the several different turbulence models comprising the multi-equation approach, the kinetic energy-dissipation (K - e) model appears to be adequate for the work carried out i n this investigation. It permits the analysis to stay at a reasonably simple level and appears to compare in accuracy with some of the more involved turbulence models. The kinetic energy-dissipation, model has already been formulated and much used for constant property fluids. The formulation of this model was re-examined here for applica-tion to supercritical fl u i d s . This examination starts at the Reynolds stress equation. Although this equation i s not used directly, i t i s required for introducing thin shear flow simplifications in the turbulence kinetic energy and dissipation equations. It also provides a necessary link between the desired correlations and the turbulent properties (kinetic energy and dissipation of turbulence). Further, the turbulence kinetic energy equation then does not have to be derived separately since i t i s a direct consequence of the Reynolds stress equation. A l l exact forms of the relevant equations were re-derived in this investigation for variable property flows since so * 3U. 3U. * 1 1 € = V 3X^ 3X^ 5 K " 2 U i U i • ** - 9 U 1 3 U i Thin shear flow simplifications imply >> U^ ; "§x~ > >'§x — 32 far, the literature has been concerned primarily with constant property situations. The turbulence correlations that arise because of the fluctua-tions in physical properties have been neglected i n the derivations of this chapter in the same manner as they were neglected in Chapter II. However, i t was deduced from the computed results (as discussed in Chapter VI) that a fluctuating property correlation that arises due to buoyancy (that i s , p'U_.g^ ) in the Reynolds stress equation may not be negligible under the influence of free convective effects. This correlation was neglected because the primary interest in this investigation i s forced convection and also because the science of turbulent free convection s t i l l appears to be at the development state for constant property fluids. The Cartesian tensor notation i s used in this section but towards the end of this chapter, the f i n a l equations employed in the numerical solution are also expressed in the conventional notation. 3.2 Reynolds Stress Equation A comprehensive discussion of Reynolds stress equation and i t s 46 closure for constant property fluids can be found in Launder et a l . The same equation was derived in this investigation for variable property fluids. The details are provided in Appendix B. An equation for fluctuating motion U i s obtained by subtracting the Navier-Stokes equation from the Reynolds equation. The resulting equation i s : An equation for the transport of Reynolds stresses U^ U_. i s often called the Reynolds stress equation. 33 3U 3T au. 1 + D, 1 Jk 3: 1 3pJ_ 1 _3_ , 8 U i . p 3X± " ^ W ^ v } p 3 X k - ^ 3U. -- u, ~^ + -=- — (PU.U, - piru,) k 3Xk p 3^ v ^ i " k ^ i T k ' (3.1) Multiplying equation (3.1) by U_. and adding to i t the corresponding equa-tion for U. multiplied by U^, leads after rearrangement and averaging (as shown in Appendix B) to the Reynolds stress equation: D U i U i DT au. = - (u.u, j"k 33^ + U i U k SX^ 3U, 3U. - 2 V — i - J -8 \ 3 X k p.- au 3u + — (—- + — P ^ ax. 3x.; 3 i 1 3 p 9 x k (Production) (Dissipation) (Pressure-strain or Redistribution) 3U.U. p u.u.u, - y - 5 J - 1 - p' (6.,u. + 6.. u.) x j k ax^ ] k 1 l k j (Diffusion) (3.2) The preceeding equations account for variation in properties p, y, and V but neglect the correlations with fluctuation in properties. The t i t l e s provided in brackets on the right-hand side are the physical interpretations usually associated with these terms. Except for some differences in the diffusion term, this equation i s very much like i t s constant property counterpart. In the latter case, diffusion i s normally represented as: ) Xk 3U.U-. — U.U.U. - V -—^ + V ( 6M U - + 6 - u U - ) i j k 3X, P ] k 1 i k j 3 X k (3.3) From a purely physical argument, this difference i s reasonable since variations in density from point to point in the flow should influence 34 the turbulent diffusion process associated with the Reynolds stresses. It is fortunate that the remaining terms are alike indicating that similar modelling practices for these terms appear suitable. Accordingly as a preliminary step, these terms were modelled as in earlier works 48 Viscous dissipation was described using Rotta's result for isotropic turbulence with the assumption that fine scale motions are isotropic: The pressure-strain correlation has been modelled by Launder et 46 49 al and Naot et a l using different approaches but with the same result. Namely: 2v (3.4) , dU. 9U. i j s l (3.5) where, (3.6) i (8 C2-2) - ( ^ . - f P c . . ) , (3.7) 11 where, 3U. 3U. 35 au au D i j " " Ul Dk 3X7 " Yk M l f » < 3 * 9 > 3 1 and au. p = " Yk ax; ' ( 3 - 1 0 ) 46 As pointed out by Launder et a l , one or more terms in the pressure-strain 37 correlation may be removed without altering the essential feature that the term as a whole be redistributive. Therefore, they proposed a simpler version of the pressure-strain correlation: <f>< . , - - Y(P. . - | P 6. .) . (3.11) i j , 2 i j 3 i j This form i s convenient because of i t s algebraic simplicity and appears adequate for algebraic stress modelling to be discussed later. Another reason for retaining the latter proposal i s that there i s considerable controversy regarding the value of in equation (3.7). However, when (3.11) i s used, y 1 S found directly from Crow's"*^ result for suddenly strained isotropic turbulence . When the f u l l pressure-strain correlation 46 i s used, Crow's result i s satisfied irrespective of the value of Presence of a nearby wall results i n additional effects on the pressure-strain correlation. This effect was not considered in the present work because i t was f e l t to be a minor consideration i n the thin shear flow model to be used. Only one of the shear stresses i s of interest in this model and i t s transport effects are found from the turbulent kinetic energy Equation (3.11) satisfies Crow's result when y = 0.6. 36 equation. The pressure-strain correlation drops out entirely when the Reynolds stress equation i s contracted to yield an equation for the turbulent kinetic energy. This i s also expected from a physical point of view since the effect of the pressure-strain correlation i s to redistribute turbulent kinetic energy from one component to another without influencing the overall balance of turbulent kinetic energy. 46 The diffusion terms have been treated by assuming that pressure and viscous diffusion are negligible; a result that has found further confir-mation by Irwin's experiments"'"'". The task then remains to model the t r i p l e correlation appearing i n equation (3.2). Various proposals for approxi-52 mating U^U^U^ are reviewed by Wolfshtein, Naot, and Lin"". The task of selecting or formulating this correlation in terms of known parameters was simplified in this analysis because the area of interest in this inves-tigation i s thin shear flows. A relatively simple proposal for this corre-53 lation i s due to Daly and Harlow , „ 3U.U. = ~ C ' r D, U n - r i - 1 • (3.12) U.U.U, s e k I 3X0 i J k I The general v a l i d i t y of this equation was not questioned i n this investiga-tion because the only place i t found use was in the reduction of equation (3.2) to the turbulent kinetic energy equation . In that case, this equation simplifies to a gradient diffusion hypothesis. That i s , the rate of transport of turbulence energy i s taken as proportional to the spatial gradient of the same quantity. This hypothesis appears to be physically 54 r e a l i s t i c and has been used by several investigators starting from Prandtl * During the algebraic stress modelling simplifications, an explicit state-ment of diffusion i s not required. 37 Therefore, equation (3.12) i s retained here for the time being. Equation (3.12) has been used with just about equal success i n thin shear flow computations as compared to another considerably more complex expression 46 for diffusion by Launder et a l Putting a l l these proposals together, the Reynolds stress equa-tion can be written asi D U.U. i .1 DT U.U. 9U. l 3U. j " k 9^ + u i u k 91^ i V - c i i-VJ • 3 5 i i K ) + «id.i + *u.2>+T^;(p!v£ ax 9U.U. ^ • ) -I (3.13) 3.3 Turbulence Dissipation Equation An equation for the dissipative correlation £ was obtained, by the following procedure. Equation (3.1) i s f i r s t differentiated by the „ 9U 9 i operator . ( ) and then multiplied by ^ v .. After some rearrangement 3 Y and averaging, this results in the desired equation. Details of this deri-vation are shown in Appendix B. The dissipation equation i s thent V D(£/V) DT V 9 p 9Xk V k 9U. u3U. (i) Z 3X^ 9X^ ^p 3X ' (ID 9U. 9U. 3U, _ 2 v ( — — 9Xk 9X£ 9XS 9U„ 3U, + 9X i 9 X k (III) (3.14) 38 - 2v U 3U. 3 U. 1 I k 3X^ SX^X^ 2v 1 3 3U. „ 3U. p 3Xk 3X£ 3X£ 3X k ; (IV) (V) 3U + 2v i 3 1 3 3U. 3X£ 3X£ p 3 X k ^ 3X^ > (VI) - 2V 3U 3U 3Ufc (VII) (3.15) The most outstanding difference between the present and constant property case i s that i t was possible to derive an equation for the ratio of d i s s i -pation to kinematic viscosity, e/v, rather than e by a derivation procedure that i s essentially similar to the constant property case. Of the terms involving source and sink terms III and IV i n equation (3.15) are the same as in the corresponding constant property dissipation equation. These have been shown to be negligible by Tennekes and Lumley"'"'. For the case of constant property fluids terms V and VI of equation (3.15) become equal to 2 3 U. - 2 (V 3X k3X £ ) ' Term VII i s similar to i t s constant property form. The source and sink terms were modelled collectively by Hanjalic 56 and Launder as follows: 39 .2 2 3U.8U.8U 3 U. 2 The same expression was adopted in this investigation for terms V, VI, and VII of equation (3.15). Since the right-hand side of equation (3.16) does not contain any d i f f e r e n t i a l operators, i t i s d i f f i c u l t to conceive a different result for variable property fluids. A more rigorous treatment of the source-sink combination i s d i f f i c u l t here because the modelling of equation (3.16) i s a result of empirical information. Terms I and II of equation (3.15) contribute towards the d i f f u -sion process. In these terms, the viscous diffusion and pressure d i f f u -sion were assumed to be negligible by Hanjalic and Launder"*^ and the following formula was derived for U^e': K — - 3e V ' - ^ l V i f ' (3-17> For the case of thin shear flows, this equation reduces to a simple gradient diffusion hypothesis essentially similar to the gradient diffusion hypothesis for turbulent energy. However, equation (3.17) i s not compatible with equation (3.15) since the variable under consideration is e/v rather than e. Therefore, i f equation (3.16) i s modified to include the gradient of (e/v) rather than e, the modelled dissipation equation for variable property fluids could be written as: 2 D(e/v) v 3 ( K rr-i— 3e/v. . P 1 W  V ~ D T ~ " p 3 i ^ ( p C e l ¥ i ^ + ( C e l 'l-Qa*T'. ( 3 ' 1 8 ) This compares with the constant property dissipation equation, which is 40 Without further examination, the additional modifications embodied i n equation (3.18) must be tentative. A choice had to be made i n this i n v e s t i -gation between equations (3.18) and (3.19). The task of numerical computa-tion i s essentially similar in both cases and i f the variations of kinematic viscosity are neglected in equation (3.18), i t too becomes an equation for the transport of the dissipative correlation e. In supercritical fluids, the properties y and p vary rapidly around the c r i t i c a l point but fortunately the variations of v is comparatively moderate as illu s t r a t e d in Figure 7 . Therefore, as a f i r s t step variations i n V were neglected in equation (3.18). The variations i n density however are s t i l l maintained inside the dif f e r e n t i a l operation. Equation (3.19) has already been used successfully ft* in moderately varying property situations 3.4 Algebraic Stress Modelling and Thin Shear Flow Simplifications An attractive simplification of the Reynolds stress equation i s a 58 procedure termed algebraic stress modelling by Rodi . This simplification i s effected in the Reynolds stress equation by introducing the assumption that variations in U.U./K are small compared to variations in O.U. i t s e l f . i j 1 3 Therefore, since M _ r(K) •= P - e , (3.20) Hence, Dotted lines i n Figure 7 arise because of uncertainty i n precise values of viscosity (see Appendix A). ** 57 See for example Khalil and Whitelaw . FIGURE 7. Properties of Carbon Dioxide Near the Pseudocritical Temperature. 42 n U.U. ^(u.u.)- r(u.u.) =~Yj-(P - e) . (3.2i) Introducing equation (3.21) in equation (3.13) yields an algebraic expression. for the Reynolds stresses U.U., 1 J U.U. - 2/3 6.. K . K " I fry " f P 6i-j> » (3.22) where A = 1 - Y (P/e - 1 + C x) In obtaining equation (3.22), the simpler version.of the pressure-strain correlation equation (3.11) was used. Further simplifications for thin shear flows yields: K 2 8 U 1 U XU 2 = f(P/e) , (3.23) where „ (C - 1 + y P/e) f(P/e) = ± (1- Y) — 5-. (3.24) (C - 1 + P/e) z Equation (3.23) makes possible the definition of an effective viscosity type of hypothesis. In the present study, a similar equation was derived using the f u l l pressure-strain correlation. The f i n a l result bears the same features as equation (3.22) but in this case the function f(P/e) i s somewhat more complicated: 43 1 f 2 1 8 1 5 (p/e - i + c r where, f 2 = J33(30 C 2 - 2) (P/e - 1 + + 10{(12 - 15 C 2)P/e + 11^ - 11}(5 - 9C2) • (3.26) Details of the derivation are provided i n Appendix B. Equation (3.23) is the basis of the kinetic energy-dissipation 59 model described by Launder and Spalding . In this model (which was at f i r s t probably conceived through physical arguments rather than through algebraic stress modelling), the function f(P /e) i s replaced by a constant* At f i t 3 t sight, this appears to be a drastic simplification of the algebraic stress model proposals. However, in actual practice, i t does not make a significant difference whichever procedure i s used. Preliminary calculations were made in this investigation using a constant value of f(P /e) (equal to 0.09) as well as with equation (3.24). In the latter case, was assigned a value equal to 2.4863 so that f(P /e) w i l l be equal to 0.09 when the pro-duction to dissipation ratio i s equal to unity. Computed results showed that there was no significant difference between the calculated velocity and temperature profiles. The reason for this i s that i n near wall flows, pro-duction to dissipation ratio retains a value close to unity i n most of the calculation regions. However, towards the edge of the boundary layer i t rapidly drops to zero but here the velocity gradients are very small and therefore errors in effective viscosity do not reflect i n the calculated 44 results. Therefore, for the sake of simplicity, a l l calculations reported in this investigation were made with a constant value of f(P/e) (equal to 0.09). A knowledge of turbulence energy K and the dissipation e i s necessary before the required Reynolds stress can be computed through equation (3.23). These quantities are determined through d i f f e r e n t i a l equations which can be obtained i n the following manner. By a contraction of the tensor and introducing thin shear flow simplifications in equation (3.14), leads directly to the turbulence energy equation for such flows: DK ~ ~ 1 , s 8 / K 2 3K ^  m = ~ Y 2 9X^- £ + T 9 x 7 ( p e U2 ^ V ( 3 ' 2 7> 2 Now substituting M in accordance with algebraic stress model proposals and using a new constant a', equation (3.27) becomes: K DK 1 3 ,_p_ K 2 3K. r — - 3 U 1 , ... Df = P3X7 to? ~ U1 U2 3X7" £ * ( 3 ' 2 8 ) Similarly, the dissipation equation becomes DT " p 3X0 K£ M? + ( C£l £ _ C £ 2 ) K * ( 3 , 2 9 ) 2 £ 2 Or in line with the notation used i n Chapter II, equations (3.28) and (3.29) can be rewritten as: 3K . 3K 1 3 / a K2 3K. —rr 3u u — + v - r - = — - - r - (r p — r -5-) - u'v' - r - - £ , (3.30) 3x 3y a 3y £0^ 3y 3x 45 and 2 2 9£ . 9£ 1 9 /_0l K 9E-, , /•„ P o \ E /o o i \ u — + v — = — (r p — r — ) + (C , — - C „) — . (3.31) 9x 9y a 3y ^ ea' 3y el e e2 K pr 7 £ Equations (3.30) and (3.31) involve four empirical constants. These are found by reference to standard sets of experimental data. They must necessarily have the same values as for constant property fluids i f these 60 c equations are to have general v a l i d i t y . A recent study by Stephenson for turbulent pipe flow uses the following values for the empirical constants: c£ = 1.0/C^, = 1«0/C U > C £ = 1.9, and C y = 0.09. These values were also used in the present study. The function f(P/e) i s written as since i t i s regarded as a constant for the purpose of compu-tations. Constant C is then found by reference to near wall data i n £1 plane flows which yields the relationship, * 2 - ^ 2 - C el» C f V < 3 - 3 2 ) In the present work, equation (3.32) was replaced by another similar equation. This i s discussed in the next section. Equations (3.23), (3.30)j and (3.31) comprise the kinetic energy-dissipation model of turbu-lence used for computations in this investigation. 3.5 Boundary Conditions for the Kinetic Energy-dissipation Model of  Turbulence Wall boundary conditions for the turbulence kinetic energy and dissipation equations are derived by noting that close to a wall convection 38 and diffusion of turbulent kinetic energy can be neglected , and production 46 and dissipation of this quantity are nearly in balance. This leads to the following relation. e = C 1 / 2 K —• . (3.33) U 8y In the laminar and transitional region close to the wall, use i s made of the Van Driest hypothesis but at the point where the turbulence model takes over, the two regions are matched by equating the respective effective viscosities. This leads to the relation for this region. K 2 2 2 du • " . Cy P T = P K y dy ' (3.34) The point at which the two solutions are matched i s sufficiently removed into the f u l l y turbulent region so the exponential damping term was not necessary + + in equation (3.34). The region of matching was between y = 75 and y - 150 for most computations carried out in this investigation. Another relation that applies i n this region i s 2 x - P<2y2(|^) • (3*35) Equations (3.33), (3.34), and (3.35) can be combined to yi e l d T = C 1 / 2 pK . (3.36) Further, to ensure complete consistency between the near wall and f u l l y turbulent region, the practice in the past has been to use equation (3.32) for obtaining the numerical value of constant C£^. However, occasional i n s t a b i l i t i e s in the solution caused a re-examination of equation (3.32) in the present study. It is possible that the problem never arose before 47 because of a different type of numerical solution procedure. Equation (3.32) was obtained with reference to plane flows where the shear stress close to a wall i s nearly constant in the absence of a pressure gradient. The relationship between the two empirical constants and C ^ was reformulated i n the present study with reference to f u l l y developed pipe flow. In this case, the shear stress distribution can be expressed as T - T.r/R . (3.37) Using this shear stress distribution and the results expressed i n equations (3.33), (3.34), (3.35), and (3.36) in the relevant dissipation, equation for constant property, f u l l y developed pipe flow yields 2 K 1 + ^ + %L R JI o'(C 0 - C _) C3 / 2 . (3.38) e e2 el u The procedure of obtaining equation (3.38) i s given in Appendix B. There was no indication during the modelling of the dissipation equation that i t is not generally v a l i d for both plane and axisymmetric flows. Therefore, the differences between equation (3.32) and (3.38) appear contradictory. The most l i k e l y cause of this difference should then not be the dissipation equation i t s e l f but equations (3.34) and (3.35) which follow from the near wall algebraic formulation. Let us therefore suppose that the near wall mixing length type distributions are different in plane and pipe flow and that in pipe flow i t i s expressed through 2 - U^UJ = K2yZ 0 f ( y ) 2 , (3.39) 48 where f(y) i s some unknown function. Now i f the preceeding derivation i s repeated using equation (3.39), then the value of f(y) could be recovered. This would s t r i c t l y involve solving a cumbersome dif f e r e n t i a l equation. Fortunately, this does not have to be done. The value of f(y) can be deduced approximately by an examination of equation (3.38). It should be f (7) = 1 — . (3.40) Therefore, the near wall mixing length distribution for pipe flow can then, be expressed approximately as Jo/R = K y / R . (3.41) / l + 3y/R + 3y 2/R 2 Equation (3.41) i s compared with a well known empirical result due to Nikuradse for f u l l y developed pipe flow, Jo/R - 0.14 - 0.08(1 - y/R) 2 - 0.06(1 - y/R) 4 . (3.42) Table 3.1 compares some calculated values from three different proposals in the near wall region. The last column of Table 3.1 contains values according to the standard plane flow proposal. Equation (3.41) reduces to the plane flow equation when the radius R i s large. The close agreement between the values contained in the f i r s t two columns of Table 3.1 appears to validate the argument presented in this section. 49 TABLE 3.1 Mixing Length Distribution in Turbulent Pipe Flow y/R A/R Equation (3.41) (This Investigation) Jt/R Equation (3.42) Nikuradse 5,/R = 0.4 y/R 0.01 0.0039 0.0039 0.0040 0.015 0.0059 0.0059 0.0060 0.02 0.0078 0.0078 0.0080 0.05 0.0186 0.0189 0.020 0.07 0.0253 0.0259 0.0280 ...... . A further comparison for the f u l l pipe i s illu s t r a t e d in Figure 8. The close agreement almost throughout the pipe i s surprising since the discussion and derivations of this section were only intended for the ;near wall region. The procedure used for obtaining equation (3.41) may be of significance in other investigations where a mixing length distribution i s required i n flows of arbitrary but known shear stress variation. However, for this work, the main consideration was formulating consistent constants for the dissipation equation. In line with the previous argument, this can be done by using either of the following procedures: a) by using equation (3.42) as the near wall algebraic formulation for pipe flows and equation (3.32) for obtaining the value of c e l J b) by continuing the use of plane flow algebraic formulation but using equation (3.38) for obtaining the value of C£^-50 FIGURE 8. Mixing Length Distribution in Pipe Flow. 51 The f i r s t procedure has the advantage that i t appears to be more correct. However, this i s not a significant advantage since the value of y/R at the point where the two solutions are normally matched i s small (about 0.05) and therefore the value of a n d t n e subsequent computed results should not be significantly different. The second procedure which is an approximation i s computationally simpler. While developing the numeri-cal solution procedure, computations were i n i t i a l l y made with both equations (3.32) and (3.38) and the results were found to be of comparable accuracy. The latt e r procedure was f i n a l l y selected for numerical solution because of i t s computational simplicity. Figure 9 shows the various possible values of Starting profiles for the turbulent kinetic energy equation were found by f i t t i n g an arbitrary but smooth monotonic function between a known value close to the wall (equation (3.36)) and another at the axis of symmetry found by reference to Laufer's experimental data*'"'*. Starting dissipation profiles are then found through equation (3.32). There appears to be as yet no standard procedure for obtaining these profiles. However, the downstream computations are not normally very sensitive to inaccuracies in the starting profiles. This i s illustrated in the next chapter i n the section on preliminary results. It may be pointed out that i n formulating the boundary conditions, the primary consideration was to obtain a workable procedure and that this area bears further inquiry. -52 FIGURE 9. Constant C for the Dissipation Equation. 53 3.6 Summary of Mathematical Formulation TABLE 3.2 Summary of Mathematical Formulation PRINCIPAL UNKNOWNS v, cross-stream velocity u, stream-wise velocity h, enthalpy p u v', turbulent shear stress p h'v', turbulent heat flux p, pressure K, kinetic energy of turbulence £, dissipation of turbulence energy p, density C , specific heat P u, viscosity k, thermal conductivity RELEVANT EQUATION IF ANY Equation (2.1) Equation (2.2) Equation (2.3) Equation (2.10) i n the near wall region and equation (3.23) i n the fu l l y turbulent region Equations (2.7) and (2.8) Equation (2.9) Equation (3.30) Equation (3.31) Tabulated property data Tabulated property data Tabulated property data Tabulated property data Chapters II and III describe the mathematical formulation used in the present study. The primary unknowns have been collected together in Table 3.2 along with the relevant relationships for these variables. 54 3c7 Concluding Remarks It i s apparent in review of this chapter that the multi-equation approach to turbulence modelling relies to a great extent on empiricism and approximation. Nevertheless, i t was decided to make use of this approach in preference to the algebraic formulations for the following reasons: a) this approach i s an improvement over the algebraic formulations, * particularly so when abnormal phenomenon are encountered ; b) i t i s f e l t that this approach holds good prospects of future improvement; c) since this approach involves empiricism and approximation, an important test of i t s usefulness i s i t s a b i l i t y to predict known experimental results. So far in constant property situa-tions, the agreement between predicted and measured results has been surprisingly good. It appeared worthwhile to test the same methods under the influence of severe property variations. A number of different proposals are available for adoption i n this area of turbulence modelling. Here the selection was governed by practical considerations. It was f e l t that the selected approach should be reasonably simple, at least as a f i r s t attempt of introducing this approach to an already complex phenomenon. An abnormal phenomenon encountered in supercritical heat transfer i s the existence of maximas in velocity away from the axis of symmetry. In such situations, the mixing length type algebraic formulations imply zero turbulent fluxes at the location of the maximas. A heat flux equal to zero does not appear to be r e a l i s t i c . 55 IV. THE NUMERICAL PROCEDURE 4.1 General Partial d i f f e r e n t i a l equations governing turbulent supercritical flows are non-linear and intri c a t e l y coupled. The velocity u occurs i n a l l the equations. The density p (which i s a function of enthalpy) i s also featured i n a l l the equations. The dissipation e appears i n the turbulence energy equation K and vice versa, while the turbulent viscosity depends on both. Besides this simultaneity, there is no obvious equation from which pressure can be extracted. In such a situation, approximate analytical methods cannot be profitably employed and only numerical procedures appear to have any hope of success. In the past decade a certain number of numerical schemes have proven to be successful in handling two-dimensional turbulent flow problems. A review of the various methods available for the analysis of turbulent flow heat transfer has been made by Patankar and 43 Spalding . The approach in this investigation was to select an existing numerical scheme and adapt i t s formulation for the problem being i n v e s t i -gated. 4.2 Required Capabilities of the Numerical Procedure Any numerical scheme for the problem must meet the following requirements: a) i t should be applicable to two-dimensional turbulent flows and preferably allow for the treatment of both axisymmetric and plane flows; 56 b) i t should permit calculation of confined flows; c) i t should be applicable to variable property flows; d) i t should be efficient in terms of computer time usage so as to be a practically useful tool for both research and engineering design. • The procedure selected to satisfy these requirements i s based on 62 the DuFort-Frankel explici t f i n i t e difference formulation . It has been 63 developed and used for two-dimensional turbulent flows by Pletcher „ and 40 41 Pletcher and Nelson ' . It was selected i n preference to the better known 43 Patankar-Spalding procedure . The Patankar-Spalding procedure normally makes use of a one-dimensional analysis close to the wall which does not appear to be suitable for supercritical fluids. The Pletcher scheme has the advantage that i t i s very flexible and permits easy modification for testing a variety of different proposals. Being based on an explicit formulation, i t permits an easier development of a computer code. An important point to note while selecting a numerical scheme is that existing computer codes cannot easily be used or modified and the only reasonable course i s to develop a new code. In this regard, a simpler formulation i s useful since already the problem i s a complex one because of the fact that five p a r t i a l d i f f e r e n t i a l equations and several auxiliary relationships constitute the physical model. A point in favor of the Patankar-Spalding procedure was that i t has already been tested on the kind of turbulence model being employed. However, this was not considered to be a major criterion for selection since i t was f e l t that i f a scheme is successful with the mixing length 5 7 type of algebraic formulations (such as the ones employed in Pletcher's scheme), i t would also be successful with higher order turbulence models. The numerical scheme selected has very good s t a b i l i t y performance. Further, i t does not make use of iterative or matrix-type solutions and i s therefore reasonably efficient in computer time usage. A more complete treatment of the f i n i t e difference scheme can be found i n the literature already quoted; only the important points are covered in this chapter. 4 . 3 Finite Difference Formulation The grid selected was orthogonal. Spacing between consecutive grid lines i n either x or y direction was varied but only gradually. If a central node i s denoted as (I, J ) , remaining nodes around i t are denoted as shown in Figure 1 0 . Numbers next to each node are used to denote these points and were used in the computer program, as consistently as possible. In order to solve for information at node 5 , a l l information about the flow has to be known at nodes 1, 2, 3 , and 4 and some at nodes 6 and 7 (density and velocities). The f i n i t e difference formulation for the energy equation i s as follows: ( h i + i , j " V i , ^ ( h i , j+ i - h i , j - i } P I , J U I , J Ax + + Ax_ P I , J V I , J Ay + + Ay_ ( p I + l ~ P I - 1 } , 1 = U T , 7— \—T~ + I,J Ax^ + Ax R T Ay, + Ay *r — J T ( RJ + i + V h I , J + l - °- 5 ( hI +l,J + hI-l,J> Ay + 58 ( R J + R J - 1 } J k e f f I,J , k e f f I.J-1 4 1 C _ T C 1 p I,J p I,J-1 Ay + u " i . J + l ^ I . J - l ) • eff I,J ^  Ay + + Ay_ 1 (4.1) The convective and diffusion terms in the momentum and turbulence equations are described i n a manner similar to that in the energy equation (4.1). The remaining terms in those equations are either algebraic or only contain first-order p a r t i a l derivatives. Therefore, the f i n i t e difference formulations for the momentum and turbulence equations are constructed similarly. The momentum equation i s : ( UI+1.J - U I - 1 , J } , n v ( UI.J+1 " " l . J - l *  P I , J U I , J Ax + Ax_ P I , J I,J Ay + + Ay_ ( PI+1 " P I - 1 } + 1 Ax + + Ax_ Rj Ay + + Ay_ ( RJ +1 + V 4Ay (y eff I,J+1 + ^eff l,J> { ( UI,J+1 - °- 5 ( uI+l,J + U I - 1 , J ) } " + ^eff I . J - l * { ° - 5 ( u I + l , J + V l . J * " U I S J - 1 } (Rj + R ^ ) 4Ay (y + gp I,J eff I,J (4.2) The f i n i t e difference formulation for the turbulence kinetic energy equation i s X ' J I , J AXj_ + Ax + P I , J V I , J Ay + Ay + — T — Rj Ay + + Ay_ ( R J + 1 + R*> 4Ay + ( yK I.J+l + y K I,J } { KI,J+1 " °- 5< Kl +l,J + K I - 1 , J ) } * ^ ' ^ - " ^ I,J + ^ I.J-1> { ° - 5 ( K l + l , J 59 + K I - 1 , J ) - K I , J - 1 ) } I + U and the dissipation equation ( UI,J+1 " U I , J - i r K I,J Ay + + Ay_ P I , J £ I , J ' (4.3) ( e i + l , J " h - l . J * ( £I,J+1 " £ I , J - 1 ) P I , J U I , J Ax + + Ax_ P I , J V I , J Ay + + Ay_ (R + R ) I , J + l + I,J> { ( £ I , J + 1 " 0 : 5 ( e l + l , J Rj Ay + + Ay_ 4Ay + £ I - 1 , J ) } " 1 - ' " ^ K 1,3 + \ I,J-1> { 0 ' 5 ( £ I + 1 , J + e i - l , J > " £ I , J - 1 } 4Ay ( u i , J + l - U I , J - i r - C e l - K I , J (Ay + + A y j p ^ j e ^ e2j 2 4 j p i j i i i ! — k i ± t (4.4) K where the parameter u = 0.09 p — . The f i n i t e difference formulation for the continuity equation i s : RJ + i + R J n 4Ax + + Ax_ L ( P U >I+1,J+1 ~ ( P U )I-1,J+1 + ( p u ) I + l , J ' ( P U ) l - l . j ] + ( p v R ) I + l , J + l " ( p v R ) I + l , J Q m Ay. (4.5) To obtain a pressure equation to be used in conjunction with the 41 partial d i f f e r e n t i a l equations, the method devised by Pletcher and Nelson i was used. In this method, the momentum equation i s rearranged to yield u e x p l i c i t l y and then multiplied by p . This equation i s then 1+1,J I,J integrated over the cross-section of the pipe or channel and equated to the known mass flow rate. Since the pressure i s only a function of the axial location, i t can be factored out of the integration to yield a 60 FIGURE 10. The Finite Difference Grid. 61 pressure equation, G -/ Pi, J U I + I , J D A / 1 ! d A + ( p i + i " pi-i> / 1 d A > < 4 - 6 > Jk JA JA where: PA = 4pj j u - j - j R j ( A y + + Ay ) Ay + Ay + (Ax + + Ax„) Ay_(R J + 1 + R j ) ( y e f ^ J + 1 + y g f f ^ ) + (Ax + + A x J Ay^CRj + R ^ ) ( y r f f ^ + p e f £ ^ , (4.7) P B = - 4 P i , j R j A y + A y - ( A x + + Ax-> v i , j ( u i , J + i - U I , J -P + 2 P I , J ( A X + + A x - > A y > e f f i , J + i + u e f f I , J > + < u I . J + l - 0 - 5 u I - l , J ) < R J * - l + V + 2 p I } J ( A x + + Ax_) A y + ( y e f f ^ + y £ f f I > J . 1 . ) ( u I > J _ 1 " °-5 " i - l , ^ ( R J + R J - 1 } + 4 p l , J u I - l , J u I , J R J ( A y + + A y-> A y+ A y-+ 4 p 2 ,RT(Ax, + Ax )(Ay, + Ay )Ay,Ay g (4.8) 1 , J J T — T — T — and PC = -4p I jRj(Ay + + Ay_)Ay+Ay_ (4.9) In equation (4.6), u , T i s multiplied with p instead of P T - T because the latter quantity remains unknown u n t i l PT.-, i s determined-Once p T i s determined, an iterative correction for the mass flow rate I I I j J can be carried out. However, this i s not necessary because the f i n i t e difference formulation requires a fine grid for convergence so the errors incurred by the use of p in equation (4.6) are negligible. 62 4.4 Method of Solution It should be pointed out that although the preceeding f i n i t e difference formulations were written for axisymmetric configurations, they are also valid for plane flows when a l l R's are assigned a value equal to unity. Briefly, the method of solution was as follows. F i r s t , a l l the equations were rearranged to exhibit the relevant unknown variable e x p l i c i t l y . Thus for example, equation (4.1) was rearranged so that only the unknown enthalpy j w a s o n the left-hand side. In this form, the equations are ready for a solution. Once a l l the variables are known at two previous grid lines (i.e., at I and 1-1 locations), then the equations can be solved one at a time to yield the required information for a l l J's at an '1+1' location. The order in which the equations are solved is pressure equation, momentum equation, energy equation, continuity equation, turbulence kinetic energy equation, and f i n a l l y dissipation equation. The turbulent fluxes are formu-1 lated through algebraic relationships prior to the solution of these equa-tions. Property values are computed at each '1+1* location after the solu-tion of the energy equation but prior to the solution of the continuity equation since the continuity equation requires the value of density at the *I+1' grid location. The solution then progresses to the next grid l i n e . More description about the solution procedure is provided i n Appendix C and D. It may be noted that information about the flow i s required at two previous grid lines. Normally such information is available, but i f i t is not (as sometimes may be the case at a starting location), then a stan-dard explicit formula that requires information at only one previous loca-63 tion was used at the f i r s t step. Standard explici t formulations are also embodied in the computer program but their use was minimised because they have poor s t a b i l i t y performance as compared to the DuFort-Frankel formula-tion. Properties are stored in tabular form at each pressure. The property tables were divided in equal increments of enthalpy rather than at equal temperature intervals. The former procedure i s superior because i t permits for greater number of divisions around the pseudo-critical tem-perature where the property variation i s steep. Property values at known enthalpies or temperature are found by linear interpolation. Once a l l variables are computed at internal nodes along ah 1+1 grid line, wall values are found by imposing the appropriate boundary conditions. These are mostly of the gradient type for the energy equation. Wall values in that case are found by, fcI,l = fcI,2 " k~7~ * <4-10> This simple relationship was f e l t to be sufficient since whenever any of the par t i a l d i f f e r e n t i a l equations contain a gradient type boundary condition, the y grid has to be very fine close to the wall (with y + equal to about unity) in order to secure convergence. Starting velocity and enthalpy profiles are found through relation-ships described in Chapter II in the section on boundary conditions. Equation (2.15) i s solved through a forward difference formula i n order to obtain the starting velocity pr o f i l e i n pipe flow problems. For the starting turbulence kinetic energy pr o f i l e , certain empirical relationships 64 * were tried , but f i n a l l y a simpler practice was adopted: K = AXU + A 2 . (4.11) The constants A^ and are found by a knowledge of the kinetic energy at two locations. Near the wall, kinetic energy is known through equation (3.35) and at the center of a pipe with f u l l y developed flow conditions, i t was assigned a value equal to 0.84 T w/P* found by reference to Laufer's experi-61 mental data . Equation (4.11) automatically satisfies the symmetry boundary condition and follows the correct trend. This simpler practice was f e l t to be adequate since downstream results are not very sensitive to starting pro-f i l e s i n turbulent flow problems. In the next section, the results of a computation are illustrated i n which the starting K p r o f i l e was intentionally kept somewhat inaccurate. The starting dissipation profile i s found through equation (3.32), a practice that has often been used i n other studies (e.g., 46 Launder et a l ). It may be recalled that the f l u i d i s unheated before entering the solution region and as a result behaves l i k e a constant property f l u i d at the starting location. The present solution procedure i s applicable to both laminar and turbulent flows. For laminar flows, the only change that has to be made is to equate the turbulent component of effective viscosity to zero. The main difference in solving various types of problems i s the way i n which grid spaces are arranged. If an extremely fine grid were to be used, various different kinds of grid spacings would yield convergence. This would be a An approximate empirical relationship designed in this investigation to describe Laufer's experimental d a t a 6 1 for f u l l y developed pipe flows i s : KA = 3.9 + 3.6(y/R) - 7.9(y/R)2 + 4.72(y/R)3 - 3.5(y/R)°- 5 6. 65 costly procedure. Therefore, grids lines are arranged i n a manner so as to secure convergence with the coarsest possible grid. The spacing between the grid lines i s varied in such a way that the flow regions with the steepest variations (such as near wall regions i n turbulent flows) contain the most nodes. The particular manner in which this i s done for different types of flows i s indicated alongside the preliminary results discussed next. 4.5 Preliminary Results The computer program was developed in stages and ver i f i e d thoroughly at each step. The program development can be approximately categorized in three steps: a) development of subroutines for solution of momentum, con-tinuity, and energy equations; b) addition of a turbulence viscosity subroutine based on algebraic formulations; c) addition of subroutines for solution of pressure, kinetic energy, and dissipation equations. Several tests were conducted during the development of the pro-gram. Some of the more important ones are reported here. During most of the preliminary testing, variation of properties was suppressed by using constant values through the property subroutine. 4.5.1 Laminar Flow_Along_a_Flat_Plate In order to test the formulation of the continuity and momentum equations, boundary conditions were simulated to correspond with the 66 Blasius flat-plate problem. The leading edge of the plate was at x=0, the plate being parallel to the x-axis. A solution was attempted starting from the leading edge and uniform streamwise velocities were imposed at the starting location. Before grid spacings can be decided, some idea of the length. scales involved i s necessary. In most boundary layer problems, this scale i s taken to be of the order of the boundary layer thickness at the starting location and then the grid lines in both x and y directions are spaced in some proportion of the boundary layer thickness. However, the f l a t plate problem i s unique because the starting location i s the leading edge where boundary layer thickness i s zero. Therefore,, an alternative length scale was devised by reference to laminar boundary layer theory. This length scale was taken to be equal to the thickness of the boundary layer at a downstream distance equal to the boundary layer thickness. Such a situation occurs only once very close to the leading edge during the development of the boundary layer and occurs at a distance of approximately 25-.V/U . Having obtained a length scale, the next step was to devise the grid spacing. Equal d i v i -41 sions in the cross-stream direction were adopted as suggested by Pletcher . The grid lines in the cross-stream direction were spaced at distances equal to 6 v/U . For the streamwise divisions, most of the previous numerical studies have suggested a step size proportional to the boundary layer thickness (Ax = A^6). The constant of proportionality in this example was assigned a value equal to 0.5, although this rule does not apply at the leading edge. At the leading edge, a step size equal to 9 v/U was *Since 6 ^ 5 /^p , the length scale i s found from x 'V 5 /—^ or x ^  25 •—-CO CO CO 67 selected. The main criterion for obtaining these divisions was to ensure that at least one nodal point l i e s within the boundary layer at the f i r s t grid line downstream from the leading edge. The manner of obtaining the f i n a l grid spacings was this: what was believed to be a conservative grid division along the X-axis was decided f i r s t , then the spacing along the Y-axis was narrowed successively u n t i l convergence was secured some dis-tance downstream from the leading edge. The f i n a l grid spacing thus found worked well for widely different values of free-stream velocity and f l u i d properties. The solution found in this way was in error by about 10% at the f i r s t step away from the leading edge and then progressively improved in the downstream direction u n t i l the Blasius velocity p r o f i l e was achieved. Figure 11 compares the computed values with the results due to * L. Howarth . The agreement here can be seen to be satisfactory. Hence, i t was concluded that the continuity and momentum equations are properly formulated and that the program was functioning well up to this stage. A subroutine for solving the energy equation was then introduced in the computer program. Thermal boundary conditions i n this case were a uniform i n i t i a l temperature profile and an isothermal wall. The previous grid spacing was retained in this case. The results are shown in Figure 12. The agreement here too was surprisingly good. The next step was to try the program for turbulent flows. 4.5.2 Turbulent Flow Along a FlatJPlate It w i l l be noted in Chapter II that the general form of the _ _ See Schlichting FIGURE 11. Blasius Velocity Distribution in the Boundary Layer Along a Flat Plate. 69 FIGURE 12. Temperature Distribution on a Heated Flat Plate at Small Eckert Number for Pr = 2.57. 70 governing equations i s the same in both laminar and turbulent flows. The only change i s i n the values of the effective viscosity. Accordingly, a subroutine based on the algebraic formulations of Chapter II was added to the program for computing the turbulent viscosity. The problem attempted was flow over a smooth f l a t plate. In this case, however, the solution was started some distance downstream from the leading edge. Starting profiles of velocity, u were found as discussed in Chapter II. The normal velocity v at the starting location was computed through,the use of the standard explicit formulation used i n conjunction with the starting u profiles. Since the variation of velocity near the wall i s very steep i n turbulent flows, the grid spacings i n the y direction are normally kept i n geometric progression rather than equal divisions as in the laminar case. 63 The nodal points i n the y direction are found from the relation : • * I . J * A y Y I , J - l * ( 4 « 1 2 ) A was prescribed a value of 1.1 in this case. It i s also necessary to y , ensure that the f i r s t node nearest the wall l i e s within the laminar sub-layer (y < 5). Grid spacing in the streamwise direction were found i n the same way as in the laminar flow problem, Ax = A <5 . (4.13) x A value of .65 for A^ was found to be adequate in this case. At the starting location about thirty nodal points lay within the boundary layer. These increased progressively with the growth of the boundary layer i n a manner implied by equation (4.12). The universal velocity profile given by 71 Kays^ was chosen as a reference for comparison. The universal velocity profile i s given by: u = y for y < 5, u + = -3.05 + 5.0 £ny + for 5 < y + < 30, (4.14) u + = 5.5 + 2.5 iny+ for y + > 30. The prof i l e obtained from equation (4.14) i s compared with the computed profile i n Figure 13. It can be seen that the result i s satisfactory. Wall shear stress values and the growth of the boundary layer were also checked to see that they satisfied the following relationships approximately (Figure 14): 6 x °' 8 —- = (—) , (4.15) 0 x o o T 0.2 — = (—) • (4.16.-X X o The next stage was to add a subroutine for the pressure equation and test the program on confined flows. 4.5.3 Turbulent_Flow_in a_Pipe At this stage, a subroutine for the solution of the pressure equation (4.6) was introduced in the computer program. After some i n i t i a l runs with the algebraic formulations for effective viscosity subroutines for the solution of kinetic energy and dissipation were added to the com-puter program. This completed the development of the solution procedure as FIGURE 13. Velocity Distribution in a Turbulent Boundary Layer Along a Flat Plate. 73 1 2 3 4 5 6 7 DISTANCE ALONG PLATE x/x FIGURE 14. Growth of a Turbulent Boundary Layer and Decline of Shear Stress Along a Plate. 74 required i n this investigation. Some examples from constant property pipe flow computations are described next. Starting profiles i n this case corresponded to hydrodynamically fully-developed flow conditions and constant starting enthalpy profiles. A constant wall heat flux boundary condition was imposed. Grid lines i n the cross-stream direction were divided according to equation (4.12). The only difference in this case i s that the constant A^ i s determined in such a way that the prescribed number of nodes cover the entire pipe. In most problems covered in this investigation, about 100 nodes along a diameter of the pipe were found to be adequate. The numerical solution has to be carried out at only half of these nodal points since the rest follow from symmetry. Divisions in the streamwise direction were found in accordance with equation (4.13). The value of A^ in pipe flow problems was kept between 0.3 and 0.65 depending on the particular problem. High heat flux cases require lower values of A^. However, in the i n i t i a l portion of the thermal entry length, A^ was kept much smaller and increased gradually in a geometric progression. This is necessary because of a rapid change in wall temperature close to the starting loca-tion. For the particular example illustrated in this section, A was x equal to 0.5. The f i r s t streamwise step was kept equal to y^ ^ a n <^ increased in ratios of 1.2 u n t i l the streamwise step became equal to a value indicated by equation (4.13). An expression for Nusselt numbers due to Bankston and McEligot 45 as reported by Sastry was used as a reference for comparison. They recommended the following equation for a i r : 7 5 Nu = 0.259 Re 0' 7 2 5 Pr°' 4 . (4.17} oo This comparison i s ill u s t r a t e d in Figure 15. The computed values of Nu^ are lower than that indicated by equation (4.17). However, i t appears that in this range of Reynolds numbers, equation (4.17) over-predicts the Nusselt number by about 6%. Experimental data of Kays*'"' as well as M i l l s * ^ suggest a value of Nusselt number approximately equal to 5 177 at Re = 1.026 x 10 . The computed value of Nu^ i s about 173. Similarly, at a Reynolds number of 50,000, equation (4.17) appears to s l i g h t l y over-65 predict the data (by about 3.5%) • Therefore, i t was concluded that the computed results are as satisfactory as can be expected. Next, a test run was made with a somewhat inaccurate starting profile of turbulence kinetic energy. The starting profile was determined in this manner: in order to find the values of constants and A^ of equation (4.11), the kinetic energy at the axis of symmetry was assigned a value equal to 0.1 T /p instead of the usual 0.84 T / p . S t r i c t l y i f the w w program i s functioning well, the correct values of kinetic energy K should be achieved at a downstream location starting from any arbitrary values. However, the starting p r o f i l e should not be such that i t causes large i n s t a b i l i t i e s i n the solution procedure. Figure 16 shows that the inaccuracies i n the starting pr o f i l e appear to be corrected at thirteen diameters downstream from the starting location. The f i n a l result included in this section i s the velocity d i s t r i -bution i n pipe flow. Equation (4.11) was again used as a reference for comparison. The results are shown in Figure 17. A correct prediction of the velocity distribution i s essential because in a way, i t plays a controlling 76 300 200 3 . COMPUTED • * "111 i / S Re = 1.026 x 10 Re = 5 x 10 • ! 0 j 4 8 12 15 STREAMWISE DISTANCE x/D FIGURE 15. Nusselt Number in Pipe Flow. RADIAL DISTANCE y/R FIGURE 16. Turbulence Kinetic Energy in Pipe Flow. 77 FIGURE 17. Velocity Distribution in Turbulent Pipe Flow. 78 role in a l l of the part i a l d i f f e r e n t i a l equations comprising the physical model. From these results, i t was ascertained that the program i s func-tioning well. In the next chapter, this solution procedure i s applied, to supercritical carbon dioxide. 79 V. RESULTS AND DISCUSSION 5.1 General Heat transfer predictions for supercritical carbon dioxide made with the procedure developed in this investigation are compared with previous experimental data i n this chapter. Computations were carried out in the reduced pressure range 1.0037 to 1.098. Although numerical predictions i n this range of pressures are comparatively d i f f i c u l t because of the proximity of the c r i t i c a l point, heat transfer within this range of pressures i s the least understood. A completely general prediction scheme should be able to make r e a l i s t i c predictions at these pressures. A comparison was made between measured and calculated temperature and velocity profiles i n pipe flow, followed by a comparison of wall tempera-ture distributions. The discussion of results i n this chapter mainly focuses attention on specific features as related to the flow under consideration. Approximations and assumptions of a general nature have already been stated. The most in f l u e n t i a l of these are discussed in greater detail i n the next chapter. Calculations were made for pipes because that i s where most of the experimental data and practical applications exist. There is one in v e s t i -gation by Hauptmann'' for supercritical carbon dioxide flowing on a f l a t plate . Part of that study was an optical examination of the flow For brevity, radial velocity and radial temperature profiles are referred simply as velocity and temperature profiles. 17 A f l a t plate study for supercritxcal nitrogen i s due to Simoneau . 80 which served to establish flow mechanisms. However, i t i s not very con-venient for comparison because temperatures were reported at only one location on the plate and the flow was partly laminar with transition occurring on the plate. However, i f need arises in a future investiga-tion to make such computations, the procedure presented in this i n v e s t i -gation can be adapted for such problems since both the numerical procedure and physical model i s applicable to plane flows as well. 5.2 Velocity and Temperature Profiles Compared to wall temperature measurements, experimental data for 16 velocity and temperature profiles i s scant. However, Wood has devoted an entire thesis to the measurement of such profiles and use i s made of his data here in order to test the theoretical model. Wood conducted experiments in an Inconel pipe with an inside diameter of 0.0229 m and total length of 1.435 m. It was mounted v e r t i -cally and the flow was upward. The heated section of the pipe was .7645 m long. The tube upstream of the heated section served to establish a f u l l y developed velocity p r o f i l e . He conducted his tests at two general conditions: 2 a pressure of 74.12 bars with a constant wall heat flux of 117350 W/m ; a 2 pressure of 75.84 bars with a heat flux of 92744.45 W/m . A series of runs were conducted by varying mass flow rates and inlet bulk temperatures. Temperature profiles were measured throughout, but velocity profiles were measured only during runs conducted at the lower pressure. Profile measure-ments were made at one location (30.7 diameters downstream from the beginning of the heated section) only. ' • Conditions corresponding to four of his experimental runs were 81 selected for computations. Fully developed velocity profiles and constant bulk enthalpies were used as starting conditions and computations were carried out i n step-by-step fashion up to a point where profiles were measured. Calculated and experimental temperature profiles are compared in Figure 18 at a pressure of 75.84 bars (p/P c =1.027, t = 305.6°K). Inlet temperature i n case (a) i s 293.15°K and for (b), i t i s 298.7°K. Both results compare within 0.5% with experimental data (within 8% based on wall-centerline temperature differences). The centerline temperatures are pre-dicted well with approximately correct amounts of temperature drops to various radial locations. Most of the temperature drop, however, takes place in the near wall region which cannot be shown on such a graph. It was not f e l t necessary to draw another graph on an expanded or logarithmic scale for the near wall region since there were no measurements i n that region for com-parison. However, Wood measured the temperature of the wall and these can be used as an indication of the accuracy of prediction i n the near wall region. The measured and computed wall temperatures for case (a) are 305.04°K and 304.4°K respectively. The same quantities for case (b) are 306.1°K and 305.55°K. The close agreement here suggests that the predictions are good in the entire flow region from wall to centerline. Figure 19 compares results at 74.12 bars (p/p £ = 1.0037, t = 304.3°K). The inlet temperature in this case i s 299.81°K. Property variations in this run were more severe due to the proximity to the c r i t i c a l point as well as higher heat flux, but agreement i s s t i l l good. Inaccuracies in the property data probably increase with the proximity of the c r i t i c a l state, but i f this has influenced the results, i t i s in the near wall region. Deviations (by about 10%) in velocity profile can be seen close to the wall. FIGURE 18. Temperature Profiles for Supercritical Carbon Dioxide. 83 FIGURE 19. Velocity and Temperature Profiles for Supercritical Carbon Dioxide. 84 Wall velocity in both cases i s zero, therefore, the deviation in the velocity profile should start reducing i n the near wall region. No such helpful constraint i s available for the prediction of temperature pr o f i l e s . There was some difference between measured (306.372) and computed (305.01) wall temperature i n this case. However, the agreement here i s s t i l l close enough to indicate that the prediction procedure has worked well under the conditions considered. A detailed comparison of the entire near wall region i s not possible because of lack of measurements for this part of the flow. In his measurements of temperature p r o f i l e , Wood noted that they become less "round1 as i n l e t bulk temperatures approached the pseudocritical temperature. The total temperature span from centerline to wall i s i n d i - , cative of the overall 'roundness' of temperature profiles. The same feature has been predicted computationally as well. This effect occurs because specific heat variations play an important role in determining the precise shape of temperature profiles. In constant property fluids, a similarity can often be noticed between temperature and velocity profiles when Prandtl number i s close to unity. This happens because equations for conservation of energy (with temperature as independent variable) and momentum have a similar form. The analogous situation for variable property fluids i s approximately a similarity between enthalpy and velocity profiles. In the present case, peaks at the pseudocritical temperature. Consequently, tem-perature varies gradually with enthalpy in the v i c i n i t y of the pseudocritical temperature. Therefore.the temperature profiles look f l a t t e r . While developing heat transfer correlations for supercritical fluids, i t may be worthwhile investigating the po s s i b i l i t y of a number based on enthalpy difference 85 rather than the usual practice of employing temperature difference. An attempt to develop such correlations was not made in this investigation because i t was f e l t that the prediction procedure is not yet sophisticated enough to be sufficiently reliable and accurate for this purpose. The f i n a l set of data due to Wood that were compared showed several interesting features. Computations corresponding to this data are shown in Figure 20. The trend appears to be approximately correct, but actual magnitudes are in error up to 15% at some locations. The unusual shape of the velocity profile here is a consequence of buoyancy force with, varying density. It appears that a variety of different shapes of velocity profiles are possible i n supercritical fluids. Figure 21 illustrates some of the formations that the velocity profile undergoes as the f l u i d progresses along the pipe. This profile could have been modified further at conditions of higher heat fluxes or i f the pipe had been longer as is illustrated later in Figure 30. At the time when Wood made these measurements, he found this to be a surprising discovery and had to argue in several different ways before he accepted them as a fact, even though their 30 ' ' existence had been postulated early by Hsu . The term unusual profiles was coined by him and i t is used here although these are now widely accepted. The entire range of shapes and their development in a pipe was however probably not known un t i l this investigation. Unusual profiles have 30,35,36 been predicted before by several investigators using a one-dimensional analysis but in that case, although a profile can be predicted in which velocity f i r s t rises and then f a l l s towards the axis of symmetry leading to a global maxima away from center. A profile such as in Figure 20 in which velocity rises, f a l l s , and rises again cannot be predicted by previous methods. In practice, a profile must transform gradually from i t s 86 ,90 75 ,60 .45 .30 « CO , 74.12 BARS \ .15 ...„..<. WOOD EXPERIMENTAL q = 117350 W/m2 D = 0.0029 m G = .789 Kg/s PREDICTION 3.0 3.4 3.8 VELOCITY m/sec 4.2 FIGURE 20. 'Unusual' Velocity Profile in Supercritical Carbon Dioxide. RADIAL DISTANCES y/R FIGURE 21, Development of 'Unusual' Velocity Profiles in Supercritical Carbon Dioxide, 88 usual constant property form to an approximately one-dimensional form. This gradual transformation (Figure 21) can be predicted by a two-dimensional analysis accounting for buoyancy. The actual profile i n Figure 20 has not yet attained i t s one-dimensional shape which i t could have had the pipe been long enough. The encouraging feature i n this computation i s the correct pre-diction of trend in a flow where a variety of different types of velocity profiles are possible (Figures 22 and 30). The discouraging feature i s that actual magnitudes of velocity were not well predicted. Although an error of 15% can sometimes be considered commonplace when predicting wall temperature distributions, better agreement was expected between measured and computed velocity profiles. The main reason why a better agreement could not be achieved in this flow i s the extreme proximity of the pseudocritical tempera-ture which led to uncertainties in both calculated and experimental results. In order to attain an unusual p r o f i l e , Wood had to take i n l e t temperatures within 0.005% of the pseudocritical temperature. This i s not a general requirement for unusual profiles (which can occur with lower i n l e t tempera-tures as well i f Q/G ratios are higher). However, due to experimental limitations, the value of Q/G was not raised further in Wood's experiments. At his particular operating conditions, the result was that the f l u i d travelled down the pipe without a significant change in temperature. At this temperature, density variations are very steep which means that the slightest errors in temperature w i l l result in. significant errors in density. This effect is shown in Figure 22. In most supercritical flows, similar differences normally occur over smaller regions and therefore effect the computation to a smaller degree. The slight difference in 89 FIGURE 2 2. Temperature and Density Profile in Supercritical Carbon Dioxide Corresponding to Data in Figure 20. temperature in Figure 22 could easily have been r e c t i f i e d by redesignating * certain empirical constants .that go into the mathematical model, but this was not done since the entire set of computations in this investigation should be made with fixed values of such parameters. Another important source of error in the present computation i s the possible inaccuracies i n property data caused as a result of the proximity of c r i t i c a l pressure and pseudocritical temperature - The same d i f f i c u l t y was faced by Wood in reducing his experimental data. He found that when the profile in Figure 20 was numerically inte-grated, i t did not yield correct mass flow rate. The computed pr o f i l e has this i n i t s favour. Given the density, i t ; Tields the correct mass flow rate since mass flow constraint i s an integral part of the mathematical procedure. 5.3 Comparison of Wall Temperature Distributions As mentioned in Chapter I, supercritical fluids can produce a variety of different wall temperature distributions varying from the enhancement regimes, where wall temperature rises less slowly than bulk temperature to deterioration regimes characterised by l o c a l peaks i n wall temperature. A r e a l i s t i c theoretical model must be able to predict both these regimes. Wall temperature measurements have been carried out i n a greater * These are empirical constants such as Von Kantian's constant for which values l i k e 0.4 and 0.42 have been suggested. ** Computations made after raising the pressure from 78.84 bars to 76.75 bars showed that velocity profiles were again normal. 91 number of investigations because i t i s relatively easier to measure. From a computational point of view, i t i s somewhat more d i f f i c u l t to predict wall temperatures accurately compared to velocity and temperature profiles because they require an accurate near-wall analysis. In the case of temperature and velocity profiles, inaccuracies of the near-wall region do not always influence predictions i n the rest of the flow. Moreover, temperature pro-f i l e s must always have the same general shape, that i s , decreasing towards the axis of symmetry. There i s no such helpful characteristic when predicting wall temperatures for supercritical fluids. A. Wood's Data Wood's experimental procedure has already been discussed. He also measured wall temperatures at various axial locations by soldering copper constantan thermocouples to the outer wall of his test section. Assuming uniform heat generation and no heat loss to the surroundings, he used the following equation relating inside and outside wall temperatures: g 2 - fa<*2? - -1] (5.1) 1 - 0 . • where 8 i s the ratio of inside diameter to outside diameter. Since much of Wood's data consists of outside wall temperatures, the same rela-tion was used for calculating inside wall temperatures. The f i r s t assump-tion involved in equation (5.1) was satisfied by the low thermal r e s i s t i v i t y of inconel and the second was checked by a heat balance (which indicated a loss of less than 0.3%). His measurements indicated that some thermocouples always read high and others always read low. This apparently happens because i t is not possible to attach a l l the thermocouples i n the same t - t = - r j — e w 2k 92 fashion. Therefore, he drew a smooth curve through the data for finding wall temperatures. His smoothed values of inside wall temperatures were used in the previous section but here the raw data reduced by equation (5.1) is used in Figures 23 to 25. Figures 23 to 25 show computed and experimental values of wall temperatures corresponding to graphs (a) and (b) of Figure 18; Figure 25 corresponds to the operating conditions in Figure 19. The agreement can be seen to be satisfactory at 75.84 bars, but at 74.12 bars, i t i s com-paratively poor. A possible reason for this was mentioned earlier that i s uncertainty of property data. Figures 24 and 25 display enhancement regimes found in supercritical fluids at moderate heat fluxes. It can be seen i n these figures that after an i n i t i a l length (about 9 diameters i n Figure 24 and 12 diameters in Figure 25), the wall temperatures are almost constant whereas the bulk temperatures rise linearly. The remaining computations of wall temperature in this investiga-tion are a l l in which deterioration (local peaks i n wall temperatures) was observed. B. Data of Jackson and Lutterodt Jackson and Lutterodt 6 7 conducted their experiments at a pressure of 75.8 bars (p/p c = 1.026) in a stainless steel tube (D = .01897 m). The tube was mounted ve r t i c a l l y and consisted of 64 diameters of entry length followed by a heated length of 129 diameters. Besides thermocouples at various axial locations, thermocouples were also attached around the tube to check for circumferential variations. 93 H C0 2 75.84 BARS 308 304 300 296 292 0 6 12 18 DISTANCE x/D FIGURE 23. Wall Temperatures in Supercritical Carbon Dioxide. WOOD,. EXPERIMENTAL Q = 927445 W/m2, G = 1.11 Kg/s —"PREDICTION H. H 306.5 304.5 302.5 300.5 298.5 C0 2 75.84 BARS M WOOD, EXPERIMENTAL Q = 927445 W/m2 G = 1.034 PREDICTION . j . . I I , I mi li immsstt 24 t^ 30 12 1.8 DISTANCE x/D FIGURE 24. Wall Temperatures i n Supercritical Carbon Dioxide 307 36 ai S3 H w H 12 18 DISTANCE x/D FIGURE 25. Wall Temperatures in Supercritical Carbon Dioxide. 94 Figure 26 compares computed and experimental data from their work. The shaded area shows the extent of circumferential variations of temperature. The shading was done in this investigation by joining the original data points. Computed results are approximately within this zone. Correct behaviour i s predicted in the thermal entrance length where wall temperature rises quickly. Buoyancy force was included corresponding to v e r t i c a l upflow, although i n this particular flow, buoyancy plays a minor role. (Jackson and Lutterodt confirmed this by changing the direction of flow and noting that there was no significant change in wall temperature distributions.) S t r i c t l y because of axial symmetry, there should not be any circumferential variation of temperature in this flow. In practice, this seems d i f f i c u l t to achieve due to unsymmetry in tube wall thickness and variation of i n l e t conditions. These may cause significant circumferential variations because heat transfer behaviour close to the c r i t i c a l region is very sensitive to inlet and boundary conditions. It may be pointed out that an even tube wall thickness is necessary i n order to obtain a uniform heat flux boundary condition, because heating in this type of experimental set-up i s provided by means of resistance heating. The measurements of circumferential varia-tions has not been given enough attention i n most of the other investigations concerned with wall temperature, measurements and could at least partly be responsible for observed differences between computed and experimental results i n those cases. C. Data of Schnurr 68 Schnurr conducted his experiments in a 0.0034 m outer diameter stainless steel tube with a 0.00038 m thick wall. An entrance length of 127 diameters was followed by an 82 diameter test length. Local tube wall 95 33 u o H 31 29 27 I 2 5 l 23 21 19' ZONE OF CIRCUMFERENTIAL VARIATION .^i^S^ q D G C0 2, 75.80 BARS EXPERIMENTAL JACKSON ANp LUTTERODT PREDICTION 17820 W/m2 .01897 m .159 Kg/s 10 20 30 40 STREAMWISE DISTANCE x/D 50 60 FIGURE 26. Wall Temperatures in Supercritical Carbon Dioxide. 96 temperatures were measured by chromel-alumel thermocouples. Schnurr did not make an explici t statement about the orientation of the tube but from his references to the top and bottom surfaces, i t i s inferred that the tube was horizontal, accordingly buoyancy force was dropped from the governing equations. Computations corresponding to his operating conditions at 75.29 bars (p/p c = 1.012) are shown in Figure 27. The agreement here i s not good in several respects. F i r s t , the temperatures are too low (by about 7%) at the beginning of the thermal entrance length. Second, a slight peak in wall temperature observed both experimentally and computationally has shifted considerably to the right. This i s probably a consequence of the former error since f l u i d temperatures must reach a certain level before deterioration i s affected. Last, at the termination of deterioration, temperatures are too high. It appears that a major part of the error i s due to neglected buoyancy effects. Heat flux i s f a i r l y high in this case and density gradients i n the v e r t i c a l direction are large. Therefore, as pointed out by Schnurr, buoyancy forces must play an important role i n this flow. The computation of this three-dimensional effect has not been taken up i n the present investigation. Nevertheless, this computation shows that this could be a worthwhile point to investigate in future work. D. Data of Ikraynnikov and Petukhov Ikraynnikov and Petukhov conducted their experiments i n a stain-less steel tube of 29 mm inside diameter with a wall thickness of 1.5 mm, 78 diameters of the tube served as the test length. This was preceeded by 25 diameters (calming section according to their terminology) for allowing 97 DISTANCE ALONG TEST SECTION, m FIGURE 27. Wall Temperatures in Supercritical Carbon Dioxide. 98 the flow to get f u l l y developed. Wall temperatures were measured by chromel-alumel thermocouples and the tube was oriented v e r t i c a l l y with the flow upward. They used relatively low mass flow rates and high heat fluxes. Therefore, deterioration as well as buoyancy effects are significant. Com-putations corresponding to their data shows the development of unusual velo-city profiles. Results for measured and computed wall temperature are plotted in Figures 28 and 29. Inlet temperatures for both runs are 15°C. Figure 28 compares results at the lower heat flux. Although computed temperatures are somewhat lower (by up to about 12%), the trend i s predicted well. This flow i s quite sensitive to inlet conditions. Computa-tions corresponding to the same conditions but with an inlet temperature of about 5°C higher (not plotted here) resulted in a double peak instead of the single peak found i n this case. Figure 29 shows a plot of another result from their investigation at a heat flux approximately nine times as high. Here the correct level of temperature peak is predicted as well as the location of deterioration (within 10%). However, computed wall temperatures drop far below predicted temperatures at the conclusion of deterioration. A reason for this appears to be a free convective effect that was excluded in this analysis and is discussed in the next chapter. Its effect i s to lower predicted tempera-tures for v e r t i c a l upflow whenever i t becomes dominant. A second cause which could be p a r t i a l l y responsible i s circumferential variations of * temperature which can be large just at the conclusion of such a deterioration. * , 67 See Jackson and Lutterodt 99 q D C0 2 81.06 BARS EXPERIMENTAL IKKRAYNIKOV AND PETURKOV 8955 W/m2 0.0029 m G/A = 190 Kg/nTs PREDICTION 20 40 STREAMWISE DISTANCE x/D 60 FIGURE 28. Wall Temperatures in Supercritical Carbon Dioxide. 100 AND PETUKHOV i q = 8432 W/m2 D = 0.0029 m j G/A = 190 Kg/m2s . PREDICTION t J , I 0 10 20 30 STREAMWISE DISTANCE x/D FIGURE 29. Wall Temperatures In Supercritical Carbon Dioxide. 101 Computed velocity profiles for this flow shown in Figure 30 demonstrate how the buoyancy forces modify the velocity profiles i n this flow. These modifications are much severer as compared to those i n Figure 21 because of strong heating at the walls. The near wall regions are not shown in Figure 30 for the sake of c l a r i t y . Since the velocity at the wall i s zero and since there i s no back flow or recirculation i n any of the problems considered here, the velocity always increases with radial distance in the near wall region i n this investigation. Figure 30(f) indicates the f i r s t signs of a return to a normal velocity p r o f i l e . If the pipe i s long enough, the flow eventually heats i t s e l f out of the c r i t i c a l region and behaves l i k e a constant property f l u i d . 5.4 Scope for Further Results This chapter was mostly restricted to a discussion of wall tem-perature distributions and radial velocity and temperature profiles. The calculation of other flow properties that comprises the mathematical formu-lation i s also a necessary part of each computation (such as those in Figure 31). The necessity of recording such data was not f e l t primarily because experimental data on these quantities i s not available for compari-son. Further, i t is f e l t that an adequate r e l i a b i l i t y cannot be placed on a l l the data at the present stage. It should be pointed out that the fact that some of the data compares well with experimental data i s not a sufficient criterion that a l l the associated data i s equally accurate. Another important point that deserves restating concerns the type of computations made. These were restricted to conform to the operating conditions in previous experimental studies. However, for the purpose of practical design work, several other operating conditions (such 103 nu, I *. ^ kmrnmrn u . ^ > r > ^ ' ' n^^ W r t^Wh^M,4 * i laiipiiim...*.. .25 .3 .35 .4 .45 .5 2 2 TURBULENCE KINETIC ENERGY, K x 10 m/s 1.1 1.9 2.7 3.5 4.3 5.1 2 3 DISSIPATION £ m /s Turbulence Kinetic Energy and Dissipation Profile Corresponding to Computations Illustrated i n Figure 19. 0.3 FIGURE 31. 104 as variable wall heat fluxes) w i l l probably be of interest. It i s hoped that the present procedure w i l l prove ammenable to further extensions. It i s f e l t that in the present form, the f u l l potential of the present procedure can be exploited only when used in conjunction with experimental data. Whenever an adequate agreement can be reached between experimental and predicted results, the present procedure can be used for generating similar results at other related operating conditions. Alternatively, the present procedure can be used for carrying out preliminary f e a s i b i l i t y studies before undertaking an actual experimental study. 5.5 Concluding Remarks To conclude, calculations made in this chapter indicate that there are several ways in which the theoretical analysis i s not completely satisfactory. However, overall agreement i s considered reasonable for the calculations illustrated. Predictions made in this chapter, though broadly covering the range of forced convection heat transfer phenomenon occurring in supercritical fluids are not exhaustive, more over computations in some cases were not carried out for the f u l l range of the heated section. For example, measurements corresponding to Figure 26 extend to 130 diameters, but calculations were terminated after 60 diameters. The only reason for doing so was to e f f i c i e n t l y u t i l i z e the available computer time. Computer time usage i s directly proportional to the number of diameters covered and therefore computations were terminated as soon as what was f e l t to be the important range was covered. 105 VI. FURTHER DISCUSSION 6.1 General Variations in turbulent Prandtl number, contribution of fluctu-ating property correlations, and the influence of free convective effects in the turbulence model are three major areas leading to uncertainty i n the formulation of the theoretical model. It was f e l t that these should be discussed more f u l l y since the further development of the theoretical procedure presented in this investigation w i l l probably depend on the cl a r i f i c a t i o n of these issues. 6.2 Turbulent Prandtl Number The turbulent Prandtl number was assigned a constant value equal to 0.9 in the present work. In actual practice, turbulent Prandtl number varies with molecular Prandtl number Pr and also with geometric location in the flow. A better practice was not adopted because there i s as yet no definite procedure for doing so, even in constant property flui d s . A large number of proposals, however, can be found. Reynolds7**' has examined more than thirty different ways of predicting turbulent Prandtl number. Calculations were made using some of these in order to study their effect on calculated heat transfer behaviour. A. Graber's Proposal 71 Graber took no account of the position within the flow but allowed for the dependence of molecular Prandtl number through the formula 106 P r t _ 1 = 0.91 + 0.13 P r 0 ' 5 4 5 . (6.1) The i l l u s t r a t i o n of Figure 32 shows the calculated Prandtl number using equation (6.1) for a supercritical flow. This flow i s the same as that in Figure 26 of Chapter V. The result shows the low sensitivity of turbu-lent Prandtl number to i t s molecular counterpart. This i s also apparent directly from equation (6.1). Even though values calculated from equation (6.1) are slig h t l y lower than 0.9, calculated temperatures do not d i f f e r significantly (See Figure 33). Equation (6.1), however, i n d i -cates that the practice of using a value less than unity for supercritical carbon dioxide appears to be better founded than that of other investiga-tions which employ a value of unity. B. Quaramby and Quirk's_Proposal 72 In view of the scatter in their data, Quaramby and Quirk f e l t that i t was impossible to isolate the dependence of turbulent Prandtl number on molecular properties. However, they proposed a relation for variation across the pipe and represented i t by Pr t = (1 + 400~ y / R) . (6.2) This gives Pr =1/2 near the wall and Pr f c = 1 in the core (Figure 32). Scatter around this result i s ±0.1 in the core and more in the wall layer. A value of 0.5 at the wall i s much lower than that suggested by other investigations and leads to lower calculated wall temperatures as seen in Figure 32. 107 FIGURE 32. Turbulent Prandtl Number Variation Across a Pipe. 108 FIGURE 33. Wall Temperature Variations with Different Turbulent Prandtl Numbers. 109 C. Cebici^s_Proposal 73 Cebici extended Van Driest's mixing length representation near a wall to turbulent heat fluxes. His arguments lead to the result - KD- - exp (-y"l"/A+)l Pr = c ^ v / ^ . (6.3) - = = • « =3S K h [ l - exp(-y +/B +), At the wall, Pr i s given by 1 ':h A" As y becomes larger, the exponential terms in equation (6.3) rapidly approach zero, leading to Pr = - . (6.5) h 74 Cebici's proposal and i t s later modification by Na and Habib allow for dependence of constants i n equation (6.3) on molecular Prandtl number. Therefore, i t combines both molecular and geometric effects. Molecular effects, however, are small in this flow as indicated already. Therefore, , calculations were made with equation (6.3) assigning constant values of + = 0.44 and B =35 (suggested by Cebici). Calculations made with equation (6.3) are shown in Figures 32 and 33. Calculated wall tempera-tures are higher because of higher Pr near the wall. Trends of variation in the last two porposals are completely opposite. A decreasing Pr with radial distance in the wall layer i s in agreement with Rotta's7""* proposal, whereas a trend of increasing Prandtl number in the high Reynolds number region i s in agreement with turbulence 110 modelling theory 7 6. A rigorous evaluation of the two proposals cannot be made because of a paucity of experimental data for this quantity. Further, a l l proposals discussed here were made for constant property f l u i d s . Additional effects due to variable properties cannot be ruled out. However, un t i l more i s known about the turbulent Prandtl number, the use of a con-stant value seems to be a reasonable practice since i t i s f a i r l y successful in predicting experimental behaviour. 6.3 Enhancement Factors As stated earlier, the fluctuating property correlations were neglected during the formulation of the mathematical model. Nevertheless, i t i s f e l t that in future investigations, the need may arise to question this assumption further. Probably many of the questions connected with these effects can be answered precisely only when experimental information on these quantities becomes available. At the present time, i t i s possible to make rough estimates of some of these additional terms i n the following manner. A s t r i c t averaging of the Navier-Stokes equations yi e l d for the turbulent shear stress? T = p u'v' + u p'v' + v p'u' + u'p'v' . (6.6) I II III Terms I, II, and III in equation (6.6) are additional fluctuating property correlations that arise in the momentum equation. Term I i s the most amenable to an estimate. It i s also probably the largest of three because term II contains v as a product and term III is a t r i p l e correlation. I l l Therefore, approximating equation (6.6) as: T = p u'v' + u h'v' , t (6.7) where, The additional assumptions embodied in equation (6.7) are: density... fluctuations due to pressure fluctuations are negligible and that higher 2 2 2 order effects are small (These involve 3 p/3h and h' v' ) .Equation (6.7) can now be rearranged as: , .' • T = p u'v' 1 + <(> u PPrx 3h/3y 3u/3y or, (6.8) T = p u'v' t 1 + u ,9p, p Pr ^3u;x (6.9) The factor inside the bracket may be regarded as an enhancement factor. Equation (6.9) i s equal i n magnitude to an enhancement factor proposed 30 * by Hsu and Smith when Pr i s equal to unity . The similarity here i s surprising because Hsu and Smith obtained their expression through physi-cal arguments based on an extension of the Prandtl's mixing length concept in order to account for density derivatives. Similar considerations for the energy equation w i l l lead to an analogous enhancement factor for the turbulent heat flux as was done by Hsu and Smith. Here the simpler procedure of regarding both these enhance-See Section 1.2. 112 ment factors as equal i s adopted in order to effect a calculation. The two enhancement factors should have the same order of magnitude on appeal to an analogy between velocity and enthalpy fields. A further simplification, used here i s to regard Pr as unity i n (6.9) for similarity with Hsu and Smith's proposal. The test flow chosen for computations i s the same as in the section on turbulent Prandtl numbers. Figure 34 illustrates computed values of enhancement factor. Calculations were effected by multiplying turbulent fluxes throughout by the enhancement factor. This factor has a value of unity at the wall since the velocity there i s zero and unity near the core because density gradients there are zero. Velocity gradients are also zero at the centre, but this fact did not cause any complication i n this flow because density f a l l s to zero faster than velocity. Calculated results showed less than 1% alteration in computed wall temperatures. This i s because enhancement i s nowhere i greater than 6% and influence, on calculated results must be much smaller. Of a l l calculated results, wall shear stresses were most effected as shown in Figure 35. Wall shear stresses are slightly lower than the unenhanced case at most points because normally i t tends to compensate for increases in shear stress away from the wall. Although in this particular flow, enhancement i s not large, i t may be larger whenever heat fluxes or flow rates are larger. This factor was not used in this investigation because i t i s f e l t that i t involves uncertainties and also because in the present work, i t s effect i s probably outweighed by other approximations. However, as prediction methods improve further i n accuracy, the need for an investi-gation in this area may arise. 113 114 115 6.4 Free Convective Influences The influence of free convective (or buoyant) effects i n forced convection situations of superciritical fluids can clearly be determined when experiments are carried out in pipe flow by changing the orientation. of the flow with respect to gravity, a l l other conditions remaining the same. Such experiments have been carried out by Jackson and Lutterodt*' 7. It was found that, as buoyant effects become significant, observed upflow wall temperatures are higher than downflow wall temperatures. Further experimental confirmation of these results has been provided by Bourke, 77 Pulling, G i l l , and Denton The prediction of the correct differences between upflow and downflow wall temperatures should s t r i c t l y be outside the capability of this investigation since the sole cause of this difference i s free convec-tive effects. However, the p o s s i b i l i t y of predicting this effect was 36,78 investigated i n this analysis because of a hypothesis due to Hall According to this hypothesis, buoyancy forces in upflow modify the shear stress distribution i n such a way as to result in fla t t e r velocity pro-f i l e s in upflow. Thereby the shear generated turbulence i s inhibited. This leads to higher wall temperatures i n upflow. Therefore, the implica-tion here i s that i f calculations are made accounting for such a modification in shear stress behaviour, then the correct differences between upflow and. downflow temperatures can be predicted. This hypothesis was supported by approximate one-dimensional calculations. Hall's contribution i s con-sidered invaluable in this investigation i n so far as emphasis and stimula-tion into the study of free convective effects i s concerned. However,., recently a closer look at the problem has become possible due to the 116 a v a i l a b i l i t y of two-dimensional analytical methods as i s explained next and i t appears that an important free convective effect may have been either overlooked or de-emphasized. It may be pointed out that the difference between upflow and downflow temperatures i s not a phenomenon unique to supercritical f l u i d s . It has also been noted i n constant property fluid s . Accordingly, Bates et a l in response to the preceeding hypothesis, accounted for buoyancy i n the mean momentum balance and computed wall temperatures in, constant property fluids for both upflow and downflow through a two-dimensional numerical procedure. Their results showed a complete reversal of observed behaviour. Calculated wall temperatures i n downflow were higher than i n upflow. Such a calculation for supercritical fluids was also made in this investigation. Figure 36 shows calculated results for carbon dioxide. Predicted wall temperature differences are again reversed from observed behaviour. That i s , wall temperatures are lower for upflow. In this flow, unusual velocity profiles associated with negative shear stresses were noted for upflow. Figure 37 compares the radial shear stress distribution i n this flow at a location (x/D =30) downstream from the start of heating. Although the shear stress i s lower at some points in upflow as expected, the wall shear stresses are higher. The wall shear stresses dominate heat transfer here and result in lower wall temperatures for upflow. Therefore, the conclusion drawn here i s that modified radial shear stress distribution i s not a sufficient criterion for obtaining the correct temperature differences and another explanation must be found. It should further be pointed out that due to a constantly changing velocity p r o f i l e and the contribution of the thermal entrance length, this flow i s clearly not 322 318 314 310 302 C0 2, 75.8 BARS q: = 49000 W/m2 D = .019 m G = .14 Kg/s DOWNFLOW \ / ** * UPFLOW •••o««r 0 10 20 30 DISTANCE x/D FIGURE 36. Wall Temperatures for Upflow and Downflow. 118 one-dimensional and therefore a one-dimensional calculation procedure that neglects flow acceleration cannot be expected to be correct. Figure 38 il l u s t r a t e s the computed turbulence kinetic energy profiles for this flow (at a location x/D = 15). However, such data was not used to deduce any conclusions here because the phenomenon under consideration appears to be dominated by the near wall region. Buoyancy forces may be regarded as influencing the flow in two ways. F i r s t , by a contribution to the mean-momentum balance, and second, * by influencing turbulence generation . Unless this second effect i s included in an analysis, i t i s believed that calculation procedures may continue to underpredict upflow wall temperatures and overpredict downflow temperatures solely because of this effect. Upward heated flows are stably s t r a t i f i e d and therefore, turbulence is suppressed in such flows. In downward flow, s t r a t i f i c a t i o n i s unstable and turbulent production enhanced. Stratification in ve r t i c a l flows i s in a direction perpendicular to the predominant heat flow direction and i s therefore not very obvious. It i s probably because of this reason that i t has not received sufficient atten-80 81 tion before. Recent Russian literature ' , however, indicates that the role of stratification-induced turbulence in ve r t i c a l confined, flows (both constant property and supercritical) i s now being studied. An adequate analytical treatment of free convective effect w i l l In the momentum equation, buoyancy forces appear as +pg^ when g^ and U^ are in the same direction. If this term is retained during the deriva-tion of the Reynolds stress equation, the corresponding contribution i n that equation i s [U.p'g. + U.p'g.] or approximately B[U.h'g. + U.h'g.]. This would further result in an additional production term in the kinetic energy equation equal to + 3 U.h'g.. 119 .4 .55 .7 RADIAL DISTANCE y/R FIGURE 37. Computed Radial Shear Stress Profiles for Carbon Dioxide at P/P = 1.027 at a Location x/D = 30. cn o x § Cd 2! Cd Cd U z Cd H (Data Corresponding to Figure 36.) \S* UP " DOWNFLOV, FLOW 1 BOTH LINES CC 1 INCIDE N ^ .1 .25 .4 .55 .7 RADIAL DISTANCE y/R .85 1.0 FIGURE 38. Computed Values of Turbulence Kinetic Energy K for Carbon Dioxide at P/P = 1.027 at a Locatl c x/D = 15. (Data Corresponding to Figure 36.) 120 probably depend on the progress in handling the same problem in constant property fluids*. The analytical treatment of such effects may be the principle interest of a future investigation on supercritical f l u i d s . *See Launder 8 2 and Gibson and Launder for recent methods i n constant property free shear flows. 121 VII. CONCLUSIONS AND RECOMMENDATIONS Conclusions and recommendations of a specific nature w i l l be found interspersed throughout this dissertation. Therefore, this section has been used for touching upon the broader points concerned with this work. It i s concluded from the computed results that both the physical modelling as well as the numerical solution procedure of the present work are viable procedures but in need of further development. A change was made from previous investigations by introducing the use of a com-pletely two-dimensional solution procedure and a two-equation model of turbulence in the f u l l y turbulent region of the flow. The mixing length type of algebraic formulations have been used in the past, mainly on the basis of how well they can predict experimentally determined behaviour. Although this s t i l l remains the criterion for testing the v i a b i l i t y of the present method, an additional advantage of physical realism and broader generality can be achieved with the multi-equation models of turbulence. The methods and approximations of the present work, though physically r e a l i s t i c , can only be regarded as a beginning to apparently promising extensions. A logical extension of the physical model presented in this investigation, appears to be the use of low Reynolds number modelling (in the near wall region). In fact, many of the improvements such as extension to combined free and forced convection (refer to Chapter VI) appear to be dependent upon such physical modelling. However, i t may be realized that in view of i t s complexity, supercritical heat transfer i s not the ideal phenomenon to select for any new radical changes in forced convection theories. Rather, i t i s f e l t that new theories 122 should f i r s t be thoroughly verified at the constant property level before being considered for supercritical heat transfer. The kind of turbulence model employed in this investigation appears to be receiving the atten-tion of several investigations at precisely such a level (that i s , con-stant property le v e l ) , and i t may be worthwhile to make use of this pro-gress with perhaps only minor modifications for application to super-c r i t i c a l fluids. It i s further concluded that the consideration of buoyancy i n the mean momentum balance i s a useful method of predicting unusual velocity profiles but not entirely sufficient to account for a l l of the free convective effects such as observed differences between upflow and downflow temperatures. The use of f u l l y two-dimensional numerical procedure was f e l t necessary in this investigation in order to account for the thermal entrance length and the near wall behaviour in supercritical f l u i d s . The numerical procedure used in this work w i l l probably prove amenable to further extensions. 123 REFERENCES 1. Hall, W.B.; HEAT TRANSFER NEAR THE CRITICAL POINT, Advances i n Heat  Transfer, VOL. 7, Academic Press, pg. 1 (1971). 2. Hendricks, R.C, Simoneau, R. J. , and Smith, R.V. ; "Survey of Heat Transfer to Near C r i t i c a l Fluids", Advances in Cryogenic Engineering, VOL. 15, pg. 197, Plenum Press, New York (1970). 3. Hsu, Y.Y. and Graham, R.W.; TRANSPORT PROCESSES IN BOILING AND TWO-PHASE SYSTEMS INCLUDING NEAR-CRITICAL FLUIDS", Hemisphere Publishing Corporation (1976). 4. Schmidt, E., Eckert, E., and Grigull, U.; "Warmetransport durch Flussigkeiten in der Nahe ihres Kritschen Zustandes", Jb. Dtsch.  Luft, VOL. 2, pg. 53 (1939). 5. Schmidt, E.; "Warmetransport durch Naturliche Konrektion i n Stoffen bei Kritischem Zustrand", Int. J. Heat Mass Transf., VOL. 1, pg. 92 (I960). 6. Powell, W.B.; "Heat Transfer to Fluids in the Region of the C r i t i c a l Temperature", Jet Propulsion, VOL. 27, pg. 776 (1957). 7. Hauptmann, E.G.; "An Experimental Investigation of Forced Convection Heat Transfer to a Fluid in the Region of i t s C r i t i c a l Point", Ph.D. Thesis, California Institute of Technology, Pasadena, California (1966). 8. Styrikovich, M.A., Shitsman, M.E., and Miroploski, Z.L.; "Certain Data on the Temperature Region of Vertical Boiler Pipes at Near C r i t i c a l Pressures", Teploenergetka, VOL. 3, pg. 32 (1956). 9. Hsu, Y.Y.; "Discussion of an Investigation of Heat Transfer of Fluids Flowing in Pipes Under Super-critical Conditions", International  Developments in Heat Transfer, ASME 188-189D (1961). 10. G r i f f i t h , J.D. and Sabersky, R.H.; "Convection i n a Fluid at Super-c r i t i c a l Pressures, Journ. A.R.S., VOL. 30, pg. 289 (1960). 11. Hines, W.S. and Wolf, H.; "Heat Transfer Vibrations with Hydrocarbon Fluids at Supercritical Pressures and Temperatures", presented at A.R.S.'Propellants, Combustion, and Liquid Rockets Conference (1961). 12. Goldmann, K.; "Heat Transfer to Supercritical Water at 5000 PSI Flowing at High Mass Flow Rates Through Round Tubes", International  Developments in Heat Transfer, pg. 565, A.S.M.E. (1963). 13. Sabersky, R.H. and Hauptmann, E.G.; "Forced Convection Heat Transfer to Supercritical Pressure Carbon Dioxide", Int. J. Heat Mass Transf., VOL. 10, pg. 1499 (1967). 124 14. Green, R. ; "Forced Convection Heat Transfer from a Cylinder i n Super-c r i t i c a l Carbon Dioxide", M.A.Sc. Thesis, University of B r i t i s h Columbia, (1970). 15. Deissler, R.G.; "Heat Transfer and Fluid Friction for Fully-developed Turbulent Flow of Air and Supercritical Water with Variable Fluid Properties", Trans. A.S.M.E., VOL. 76, pg. 73 (1954). 16. Wood, R.D.; "Heat Transfer in the C r i t i c a l Region: Experimental Investigation of Radial Temperature and Velocity Profiles", Ph.D. Thesis, Northwestern University, Evanston, I l l i n o i s (1963). 17. Simoneau, R.J.; "Velocity and Temperature Profiles i n Near-critical Nitrogen Flowing Over a Heated Flat Plate", Ph.D. Thesis, North Carolina State University, Raleigh, North Carolina (1970). 18. Eckert, E.R.G.; "Discussion of Heat Transfer and Fluid F r i c t i o n for Fully-developed Turbulent Flow of Air and Supercritical Water with Variable Fluid Properties", Trans. A.S.M.E., VOL. 76, p. 83 (1954). 19. Bringer, R.P. and Smith, J.M.; "Heat Transfer i n the C r i t i c a l Region", A.I.Ch.E. Journal, VOL. 3, pg. 49 (1957). 20. Hess, H.L. and Kunz, K.R.; "A Study of Forced Convection Heat Transfer to Supercritical Hydrogen", A.S.M.E. Paper Number 63-WA-205 (1963). 21. Krasnoshchekov, E.A. and Protopopov, V.S.; "Heat Exchange i n the. Supercritical During the Flow of Piped Carbonic Acid and Water'-', Teploenergetika, Vol. 6 (12), pg. 26 (1959). 22. Gunsen, W.E. and Kellog, H.B.; "Engineering Application Technique for Supercritical-Pressure Heat Transfer Correlations", A.S.M.E. Paper Number  66-WA/HT-ll (1966). 23. Swenson, H.S., Carver, J.R., and Kakarala, C.R.; "Heat Transfer to Supercritical Water in Smooth Base Tubes", J. Heat Transfer, VOL. 87, pg. 477 (1965). 24. Miropolski, Z.L. and Shitsman, M.E.; "Heat Transfer to Water and Stream with Varying Specific Heat in the C r i t i c a l Region", Journal of  Technical Physics, Vol. 27, pg. 2359 (1951). 25. Hendricks, R.C.,, Graham, R.W., Hsu, Y.Y., and Medeirbs, A.A.; "Correlation of Hydrogen Heat Transfer i n Boiling and Supercritical Pressure States", ARS Journ., VOL. 32, pg. 244 (1962). 26. Graham, R.W.; "Penetration Model Explanation for Turbulent Forced Convection Heat Transfer Observed in Near C r i t i c a l Fluids", NASA TND  5522 (1969). 27. Deissler, R.G. and Presler, H.F.; "Computed Reference Temperatures for Turbulent Variable Property Heat Transfer in a Tube for Several Common Gases", International Development in Heat Transfer, A.S.M.E., . pg. 579 (1963). 125 28. Goldmann, K.; "Heat Transfer to Supercritical Water and Other Fluids with Temperature Dependent Properties", Chem. Eng. Prog. Symp. Ser., VOL. 50, No. 11, pg. 105 (1954). 29..Tanaka, H., Nishiwaki, N., Hirata, M., and Tsuge, A.; "Forced Convection Heat Transfer to Fluid Near C r i t i c a l Point Flowing in Circular Tube", Int. J. Heat Mass Transfer, VOL. 14, pg. 739 (1971). 30. Hsu, Y.Y. and Smith, J.M. ; "The Effect of Density Variations on Heat Transfer i n the C r i t i c a l Region", J. Heat Transfer, VOL. 83, pg. 176 (1961). 31. Yoshlda, S., F u j i i , T., and Nisikawa, K.; "Analysis of Turbulent Forced Convection Heat Transfer to a Supercritical Fluid Flowing i n a Tube", J.S.M.E., VOL. 35, pg. 3185 (1972). 32. Kato, H., Nishiwaki, N., and Hirata, M.; "Turbulent Heat Transfer i n Smooth Pipes at High Prandtl Numbers", Proc. 11th Int. Cong. Appl.  Mech., pg. 1101 (1964). 33. Shiralkar, B.S. and G r i f f i t h , P.; "Deterioration in Heat Transfer to Fluids at Supercritical Pressure and High Heat Fluxes", A.S.M.E.  Paper Number 68-HT-39 (1968). 34. Sastry, V.S. and Schnurr, N.M.; "An Analytical Investigation of Forced Convection Heat Transfer to Fluids Near the Thermodynamic C r i t i c a l Point", J. Heat Transfer, VOL. 97C, pg. 226 (1975). 35. Hall, W.B.; "The Effects of Buoyancy Forces on Forced Convection Heat Transfer in a Vertical Pipe", Simon Engineering Lab. Report NE1, University of Manchester (1968). 36. Tanaka, H., Tsuge, A., Hirata, M. , and Nishiwaki, N.; "Effects of Buoyancy and Acceleration Owing to Thermal Expansion on Forced Turbulent Convection in Vertical Circular Tubes", Int. J. Heat Mass  Transfer, VOL. 16, pg. 1267 (1973). 37. Malhotra, A. and Hauptmann, E.G.; "Heat Transfer- to a Supercritical Fluid During Turbulent, Vertical Flow in a Circular Duct", International Seminar on Turbulent Buoyant Convection, I.C.H.M.T., Yugoslavia (1976). 38. Launder, B.E. and Spalding, D.B.; "Lectures in Mathematical Models of Turbulence", Academic Press, London and New York (1972). 39. Eckert, E.R.G. and Drake, Jr., Robert M.; ANALYSIS OF HEAT AND MASS TRANSFER, McGraw H i l l Book Company (1972). 40. Nelson, M.N.; "An Ex p l i c i t Finite Difference Analysis for Developing Internal Flows with Heat Transfer and Property Variations", M.Sc. Thesis, Iowa State University, Ames, Iowa (1972). . . 41. Pletcher, R.H. and Nelson, M.N.; "Heat Transfer to Laminar and Turbulent Flow in Tubes with Variable Fluid Properties", Proceedings of the F i f t h International Heat Transfer Conference, Tokyo, VOL. II, pg. 146 (1974). 126 42. McEligot, D.M., Smith, S.B., and Bankston, C.A.; "Quasi-developed Turbulent Pipe Flow with Heat Transfer", A.S.M.E. Paper Number 70-HT-8 (1970). 43. Patankar, S.V. and Spalding, D.B.; HEAT AND MASS TRANSFER IN BOUNDARY LAYERS, 2nd Ed., Intertext Books (1970). 44. Spalding, D.B.; "A Single Formula for the Law of the Wall", Journal of  Applied Mechanics, VOL. 28, Series E, pg. 455 (1961). 45. Sastry, V.S.; "An Analytical Investigation of Forced Convection Heat Transfer in a Circular Duct at Supercritical and Near-critical Pressures", Ph.D. Thesis, Vanderbilt University (1974). 46. Launder, B.E., Reece, G.J., and Rodi, W.; "Progress in Development of Reynolds Stress Closure", J. Fluid Mech., VOL. 68, pg. 537 (1975). 47. Hanjalic, K.; "Two-dimensional Asymmetric Turbulent Flow i n Ducts", Ph.D. Thesis, University of London (1970). 48. Rotta, J.; "Statische Theorie Nichthomogener", Turbulenz Z. Phys., VOL. 129, pg. 547 (1951). 49. Naot, D., Shavit, A., and Wolfshtein, M.; "Two-point Correlation Model and the Redistribution of Reynolds Stresses", Phys. Fluids, VOL. 16 (1973). 50. Crow, S.C.; "The Viscoelastic Properties of Fine-grained Incompressible Turbulence", J. Fluid Mech., VOL. 33 (1968). 51. Irwin, H.P.A.H.; "Measurements in Blown Boundary Layers and Their Prediction by Reynolds Stress Modelling", Ph.D. Thesis, McGill University (1974). 52. Wolfshtein, M. , Naot, D. , and Lin, C.C; "Models of Turbulence", Ben-Gurion University of the Negev, Department of Mechanical Engineering Report Number ME-746(N). 53. Daly, B.J. and Harlow, F.H.; "Transport Equations i n Turbulence", Physics of Fluids, VOL. 13, pg. 2634 (1970). 54. Prandtl, L.; "Uber ein neues Formelsystem fur die Ausgebildete Turbulenz", Nachrichlen von der Akad. der Wassenschaft in Gottingen (1945). 55. Tennekes, H. and Lumley, J.L.; A FIRST COURSE IN TURBULENCE, M.I.T. Press, Cambridge, Massachusetts (1972). 56. Hanjalic K. and Launder, B.E.; "A Reynolds Stress Model of Turbulence and i t s Applications to Symmetric Shear Flows", J. Fluid Mech., VOL. 52, pg. 609 (1972) 127 57. Khalil, E.E. and Whitelaw, J.H.; "The Calculation of Local Flow Properties in Two-dimensional Furnaces", Imperial College, Mech. Engng. Department Report HTS/74/38 (1974). 58. Rodi, W.; "The Prediction of Free Turbulent Boundary Layers by Use of a 2-equation Model of Turbulence", Ph.D. Thesis, University of London (1972). 59. Launder, B.E. and Spalding, D.B.; "The Numerical Computation of Turbulent Flow", Comp. Meth. in Appl. Mech. and Eng., VOL. 3, pg. 269 (1974). » . 60. Stephenson, P.L.; "A Theoretical Study of Heat Transfer in Two-dimensional Turbulent Flow in a Circular Pipe and Between Parallel and Diverging Plates" Int. J. Heat Mass Transfer, VOL. 19, pg. 413 (1976). 61. Laufer, J.; "The Structure of Turbulence of Fully-developed Pipe Flow", N.A.C.A. Report 1174 (1954). 62. DuFort, E.C. and Frankel, S.P.; "St a b i l i t y Conditions in the Numerical Treatment of Parabolic Differential Equations", Mathematical Tables  Aid Computation, VOL. 7, pg. 135 (1953). 63. Pletcher, R.H.; "On a Finite Difference Solution for the Constant Property Turbulent Boundary Layer", AIAA Journal, VOL. 7, pg. 305 (1969). 64. Schlichting, H.; BOUNDARY LAYER THEORY, McGraw H i l l Book Company (1960). 65. Kays, W.M.; CONVECTIVE HEAT AND MASS TRANSFER, McGraw H i l l Book Company (1966). / 66. M i l l s , A.F.; "Experimental Investigation of Turbulent Heat Transfer in the Entrance Region of Circular Conduits", J. Mech. Engr. Sci., VOL. 4, pg. 63 (1962). 67. Jackson, J.D. and Evans-Lutterodt, K. ; "Impairment of Turbulent Forced Convection Heat Transfer to Supercritical Pressure CO2 Caused by Buoyancy Forces", Res. Report N.E. 2, University of Manchester, England (1968). 68. Schnurr, N.M.; "Heat Transfer to Carbon-dioxide in the Immediate Vicinity of the C r i t i c a l Point, A.S.M.E. Paper Number 68-HT-32 (1968). 69. Ikryannikov, N.P., Petukhov, B.S. , and Protopopov, V.S.; "An Experimental Investigation of Heat Transfer in the Single-phase Near-critical Region with Combined Forced and Free Convection", Teplofizika Vysokikh  Temperatur, VOL. 10, pg. 96 (1972). 70. Reynolds, A.J.; "The Prediction of Turbulent Prandtl and Schmidt Numbers", Int. J. Heat Mass Transfer, VOL. 18, pg. 1055 (1975). 128 71. Graber, H.; "Der Warmeubergang in Glatten Rohren Zwischen Parallelen Platten in Ringspallen and Langs Rohrbundeln bei Exponentueller Warmeflussverteilung in Erzwungener Laminarer oder Turbulenter Stromung", Int. J. Heat Mass Transfer, VOL. 13, pg. 1645 (1970). 72. Quaramby, A. and Quirk, R.; "Axisymmetric and Non-axisymmetric Turbulent Diffusion in a Plain Circular Tube at High Schmidt Number", Int. J. Heat Mass Transfer, VOL. 7, pg. 45 (1973). 73. Cebici, T.; "A Model of Eddy Conductivity and Turbulent Prandtl Number", A.S.M.E. Paper Number 72-WA/HT-13 (1971). 74. Na, T.Y. and Habib, I.S.; "Heat Transfer in Turbulent Pipe Flow Based on a New Mixing Length Model", Applied Scientif. Res., VOL. 28, pg. 302 (1973). 75. Rotta, J.C; "Temperaturverteilungen i n der Turbulenlen Grenzschicht an der Ebenen Platte", Int. J. Heat Mass Transfer, VOL. 7, pg. 215 (1964). 76. Gibson, M.M. and Launder, B.E.; "On a Calculation of Horizontal Turbulent Free Shear Flows Under Gravitational Influence", A.S.M.E.  J. Heat Transfer, VOL. 98C (1976). 77. Bourke, P.J., Pulling, D.J., G i l l , L.E., and Denton, W.H.; "Forced Convective Heat Transfer to Turbulent CO2 in the Supercritical Region", Int. J. Heat Mass Transfer, VOL. 13, pg. 1339 (1970). 78. Hall, W.B. and Jackson, J.D.; "Laminarisation of Turbulent Pipe Flow by Buoyancy Forces", A.S.M.E. Paper Number 69-HT-55 (1969). 79. Bates, J.A., Schmall, R.A., Hasen, G.A., and McEligot, D.M.; "Effects of Buoyant Body Forces on Forced Convection i n Heated Laminarizing Flows", Proceedings of the F i f t h International Heat Transfer Conference Tokyo, VOL. II, Pg. 141 (1974). 80. Polyakov, A.F.; "Mechanism and Limites on the Formation of Conditions for Impaired Heat Transfer at a Supercritical Coolant Pressure", Teplofizika Vysokikh Temperatur, VOL. 13, pg. 1210 (1975). 81. Petukhov, B.S.; "Turbulent Flow and Heat Transfer i n Pipes Under Considerable Influence of Thermo-gravitational Forces", International Seminar on Turbulent Buoyant Convection, I.C.H.M.T., Yugoslavia (1976). 82. Launder, B.E.; "On the Effects of a Gravitational Field on the Turbulent Transport of Heat and Momentum", J. Fluid Mech., VOL. 68 (1975). 83. Hirschfelder, J.O., Curtiss, CF., and Bird, R.B.; MOLECULAR THEORY OF GASES AND LIQUIDS, Wiley (1954). 84. Michels, A., Sengers, J.V., and Van der Gulik, P.S.; "The Thermal Conductivity of Carbon Dioxide i n the C r i t i c a l Region", Parts I and II, Physica, VOL. 28, pg. 1201 (1962). 129 85. Sastry, V.S.; "Private Communication", Department of Mechanical Engineering, Tennessee State University, Nashville, Tennessee (1975). 86. Michels, A. and Michels, C.; "Isotherms of C0 2 Between 0° and 150° and Pressures From 16 to 250 Atm (Amagat Densities 18-206)", Proceedings of the Royal Society of London, VOL. A153, pg. 201 (1935). 87. Michels, A, Blaisse, B., and Michels, C.; "The Isotherms of C0£ i n the Neighbourhood of the C r i t i c a l Point and Round the Co-existence Line", Proc. Roy. Soc. of London (see 86 here), Vol. AJ.60, pg. 358 (1937). 88. Koppel, L.B. and Smith, J.M.; "Thermal Properties of Carbon Dioxide in the C r i t i c a l Region", Journal of Chemical and Engineering Data, VOL. 5, pg. 437 (1960). 89. Kestin, J., Whitelaw, J.H., and Zien, T.F.; "The Viscosity of Carbon Dioxide i n the Neighbourhood of the C r i t i c a l Points", Physica, VOL. 30, pg. 161 (1964). 90. Sengers, J.V. and Michels, A.; "The Thermal Conductivity of Carbon Dioxide in the C r i t i c a l Region", Second Symposium on Thermophysical Properties, A.S.M.E., New York, Academic Press, pg. 434 (1961). 130 APPENDIX A Properties of Supercritical Fluids 131 APPENDIX A. Properties of Supercritical Fluids Hsu and Graham3 have reveiwed the present state of knowledge concerning properties of near c r i t i c a l fluids. A useful treatment of the thermodynamics of the c r i t i c a l point can be found i n Hirschfelder, 83 Curtiss, and Bird . Classical modelling of a f l u i d near the c r i t i c a l state has been based on Van der Waal's proposal. His model i n which an allowance i s made for the attractive and repulsive forces between molecules leads to an equation of state i n the following form: (p + -%)(V - b) = Rt . ( A ' 1 ) V The Van der Waal's equation predicts a value of i n f i n i t y for C^ at the c r i t i c a l point. The Van der Waal's equation i s a useful tool for under-standing the thermodynamic behaviour near the c r i t i c a l state but i t does not succeed i n mapping the (P, p, T) surface with sufficient accuracy. Therefore, V i r i a l type equations requiring several empirical coefficients are usually employed, P - A(t )p + B(t ) p 2 + C(t ) p 3 + D(t ) p 4 . (A.2) A mathematically more analytic equation, based on the recently popular Ising model for ferromagnets has been proposed by Vicentini Missoni et a l . For numerical work, however, i t i s probably not an eff i c i e n t practice to use these equations directly in the computer program since these equations require an interative solution. Therefore, use was made of tabulated 132 values of properties. In most supercritical heat transfer systems, the pressure i s maintained above the c r i t i c a l pressure. Therefore, singularities such as an i n f i n i t e value for specific heat are avoided. However, the property variation i s s t i l l severe in the v i c i n i t y of the c r i t i c a l point. For example, the peak in specific heat i s s t i l l large even at pressures con-siderably removed from the c r i t i c a l pressure. Like the specific heat, the thermal conductivity also exhibits a peak around the c r i t i c a l point but because of the d i f f i c u l t i e s i n measuring this quantity near the c r i t i c a l state, there was considerable controversy regarding the behaviour of thermal conductivity near the c r i t i c a l state u n t i l very recent times. Some investigators believed that there i s no spike i n thermal conductivity. However, experiments due to Michels, 84 Sengers, and Van der Gulik have demonstrated that there i s a large spike, i n thermal conductivity right around the c r i t i c a l point. Subsequently, similar spikes have been noticed i n other fluids as well. Property values used in this investigation include this spike. Compared to thermal conductivity, variations i n viscosity around the c r i t i c a l point are less drastic. It i s l i k e l y that there i s a slight detectable spike i n viscosity as well. However, because of uncertainty i n the precise magnitude of this peak, i t i s normally neglected i n engineering applications. The viscosity data obtained in this investigation contained slight irregularities but these were smoothed out prior to feeding the property tables i n the computer program. The property data for C0 ? used in this investigation are the same 133 as that used with two previous numerical studies due to Schnurr"" and 9 85 8 86 87 Sastry ' . The density data i s from Michels and Michels ' * . Enthalpy 88 and specific heat data are from Koppel and Smith . Viscosity data are 89 from Kestin, Whitelaw, and Zien and thermal conductivity data are from 90 Sengers and Michels „ Table A.1 gives an example of how the property tables were set up at p/p c - 1.027 prior to feeding in the computer program. A good agreement between the experimental and predicted heat transfer data i n a number of test runs indicates that the property data used in this inves-tigation are probably accurate enough for the work carried out i n this investigation. The cases where such an agreement was not achieved appear largely to be a result of approximations i n the theoretical model rather than inaccuracies i n the property data. However, as theoretical prediction schemes improve and more accurate property data becomes available in future, i t should be employed in numerical prediction schemes. 134 TABLE A . l * Properties of Carbon Dioxide (Pressure = 1100.00 PSIA). ** N TEMPERATURE (RANKINE) VISCOSITY (LBM/FT-SEC) SP. HEAT (BTU/LBM-R) CONDUCTIVITY (BTU/SEC-FT-R) DENSITY (LBM/CU FT) 1 483.3 .7334E-04 .496 .1904E-04 61.505 ! 2 492.6 .6939E-04 .537 .1835E-04 60.017 3 502.0 .6530E-04 .572 .1759E-04 58.321 1 4 510.3 .6169E-04 .611 .1674E-04 56.455 1 5 517.5 .5866E-04 .651 .1584E-04 54.431 6 525.0 .5512E-04 .735 .1490E-04 51.850. 7 531.9 .5132E-04 .848 * .1418E-04 49.544 8 537.3 .4559E-04 .937 .1338E-04 47.603 9 542.0 .3937E-04 1.200 .1261E-04 45.320 10 545.5 .3405E-04 1.635 .1139E-04 42.410 11 547.8 .3037E-04 2.218 •1057E-04 39.564 12 549.0 .2828E-04 2.878 .1529E-04 36.927 13 549.4 .2744E-04 3.294 .1595E-04 34.250 fe 14 549.7 .2661E-04 3.709 .1662E-04 31.575 & 15 550.1 .2564E-04 3.861 .1726E-04 28.689 16 550.6 .2437E-04 3.397 .1785E-04 25.316 17 551.1 .2298E-04 2.930 .1804E-04 22.100 18 552.2 .2038E-04 2.427 .1433E-04 20.291 19 554.4 .1630E-04 1.892 .8252E-05 18.433 20 557.7 .1461E-04 1.408 .6588E-05 16.877 21 562.2 .1407E-04 1.068 .5256E-05 15.545 22 567.6 .1375E-04 .789 .4760E-05 14.446 23 574.9 .1352E-04 .618 .4543E-05 13.350 24 .584.1 .1333E-04 .514 .4458E-05 12.246 25 595.7 .1337E-04 .446 .4404E-05 11.255 26 607.0 .1347E-04 .409 .4389E-05 10.527 27 618.4 .1359E-04 .384 .4389E-05 9.951 28 631.8 .1373E-04 .362 .4422E-05 9.428 29 645.6 .1386E-04 .344 .4460E-05 8.945 30 659.7 .1398E-04 .330 .4499E-05 8.509 31 677.8 .1416E-04 .321 .4574E-05 8.117 32 694.3 .1433E-04 .309 .4653E-05 7.798 33 710.6 .1450E-04 .297 .4732E-05 7.489 34 726.8 .1468E-04 .287 .4811E-05 7.199 35 743.1 .1486E-04 .278 •4890E-05 6.935 * The data here is described in Br i t i s h Units because this i s the manner i n which, i t was fed to the computer program. The relevant conversion factors to S.I. units are contained in the computer program and are from Eckert and Drake . **ENTHALPY = 30.0 + 5(N - 1), (BTU/LBM). 135 APPENDIX B Equations f o r the Turbulence Model 136 Reynolds Stress Equation The starting point for deriving equations for turbulent quantities are the Navier-Stokes equations. These may be written as: 3U, „ 3U, . ^ • a 3lY. -3T + \ 83^ " ' p 3Xi + p 33^ ^ y 33^ * ^•" L J The instantaneous motions can be broken into mean and fluctuating components. The resulting equation i s : 3(U + U.) ^ 3U ^ 3U. 3U 3U. „ 3T 1 * U k 3 X ^ P k 3 X ^ + U k 3^ + U k 3 X ; = - P3x7 ( p + p t> (B.2) 4 - i 3 * 3 X ^ 1 + V p 9 X k The time averaged form of these equations are the Reynolds equations * 3U. * 3U. . , ^ - r , ~ 3U 1 + U. 1 i l ^ + r ^ - (y-s^-) ( P ¥ k ) • CB.3) W " °k 3X^ " p 3X. T p 3 ^ v^ 3^' p 3 ^ v" "1 Subtracting equation (B.3) from (B.2) results in the equation for the fluctuating component U (equation (3.1)). This i s repeated here for convenience: 3U. ^ 3U. . a , , a 3U 3U -I F + u k ^ - - ^ + ^ ( w ^ - \ < + i 4 ( p U A ~ p U i V (B.4) Multiplying equation (B.4) by U and adding to i t the corresponding equa-tion for U_. multiplied by U\, leads to: 137 3U.U. ~ 3U.U. 3U. 3U. U. a , U. ~ , 1 .1 + U 1 J = - uu - i - o — 1 1-^ 2 i J L . 3T ^ uk SX .^ k j 3 ^ k i 3 ^ p 3X± p 3X U. 3 3U. U. g 3U. 3 X ^ ( p U i U k - p W - T ^ ( p V k " p W - ( B*5 ) During the simplifications, use is made of the continuity equations, 3pU. 3pU. 3pU. 1x7-°- "3x7-°- a n d l x f = 0 ' ( B* 6 ) x i i The continuity equation (B.6) allows transformations such as -^|- <F) = ( PU kF) , - (B.7) 'k 9 X k 3 ^ where F i s any general flow property. Further rearrangements of the left-hand side of equation (B.5) lead tot /v 3U.U. - 3U U. 3U. 3U. , 3U , 3U -k2- + \ ' ^ " V j 3 ^ 7 V i ^ + V3xT + V < 1 3 p , u i i n + i x ( l i D !%• p i x f - p i r + P 3 \ ( u uj V (u 3U • 3U. 3U. • , ft ^ 1 3 • /„ „ _ 4 A _ 2 ^  — — 1 - ~ i r - (P U.U.U. ) ? K i d \ p 3 x k 3 x k p 3 x k 1 3 k 138 u. a u a  + -r|- (p u.u,) + -jdr (P U 4 U k } » p 82L l k P BX^ 3 k (B.8) or, 3U.U. ~ 3U.U. 1 3 + u — i - 1 3T + Uk 3X^ 3U. 3D. U.U. + U.U, j u k 3^ "i"k 3^ - 2v 3U. 3U. , 3D. 3D 3 X k 3 X k P * Y 3U.U. p SX^ . 1 P i j k ^ 3 \ p Jk x ik j U. „ U a  + —1 (p U.U, ) + - f a f - (P U.U ) p 3X^ p x k p 3X^ 3 k (B.9) Averaging equation (B.9) leads directly to che Reynolds stress equation which i s equation (3.2). Turbulence Dissipation Equation An equation for the transport of the dissipative correlation was obtained in the following manner. Equation (B.4) i s differentiated 3 with the operator -r^— leading to: a 3 U 4 -- d i 4f V + U k ^ V + 3X7 3X~ = - 3X7 P^ 3X±' ~ °k 3X, 3^ 3^ 3X£ 3U. ' , a  3U. 3D. 3U. k x „ . 3 U, 3U^ 3D. 3 (1 321) . u = t - A ^ c i v ^  u r 3D. 3 / 1 3 (B.10) 3D. Multiplying equation (B.10) with ^ results in 139 1 JL 2 3T 3D. U. k 3 9 X k 3U. 4 3U. 3U, 3U. l k l 9 U i 3 .1 3p' 3XA 3X& 3 ^ 3X£ 3X £ 3X± - U. 3U. 3D. l i 3D, 3D. 3D. i k l k 3X£ 3X £ 3Xk 3]^ 3X £ 3X £ 3D± 3D± 3Dk 3D. „ - D. k 3Xa 33^ ^ X ^ 3X £ 3X£ ^4 ( p U i V 3D i 3 3X£ 3X^  A 3 / I s p 3 ^ U i 3X,/ (B.ll) Rearranging the preceeding equation and multiplying throughout by 2v , V 3T 3D. ' v (3X, } + VD. k 9 X k  9 U i 3 3D. ' V (3x £> 3D. 3X£ 3 X j l (y 3D. -2v ) Xk ( 9 X k 3D. 3D v 3 • p 3 ^ 3D 3D. ' ^ VU, ( ^ i ) V k dX„ - 2V 3X£ 3XA p^ W±J 3X„ 3X£ 3X i3X k^ - 2v D. 2" 3D± 3 D ± k 3X r3X & ZX^ 1 3U. 3D. 3D -2v ^ ^ ^ + 2v 3D dX^ 3X£ 3 ^ " 3X„ 3X i 3 t P 3 X k ^ ^ . 2 v 3 p 3 x k 3D. 3D. 3X£ 3X, (y 9 x k +2V 3D. ~ 1 o 3X£ 3X£ (B.12) 140 3u 2 Replacing V(~-^) by d and averaging leads directly to the desired equation which i s equation (3.14). Algebraic Stress Modelling Equation (3.22) was derived by Rodi using only one of the compo-nents of the pressure-strain correlation. A similar expression using the f u l l pressure-strain correlation (equation (3.7)) i s derived here. Using D U.U. U U i JL _ r(U.U.) = —r1 (P - e) , (B.13) PT i V " i w j ' K the Reynolds stress equation can be written as Z £L ( p . e ) , 1 eS±i . C i | ( U . D . . | S i j K ) + P t j - UP,. - § « y P ) 3U. 3U. 9 3 1 where, C„ + 8 30 C. - 2 8 C - 2 L = A i ~ ' M 55 - a n d N I i - ' • ' ( B* 1 5 ) Equation (B.14) can be rearranged as: U.U. - 2/3 6..K A. 0 A, 9 K e 13 3 13 e 13 3 13 A- 3U. 3U. + T K ( 3 r + ^ ' • ( B - 1 6 ) J 1 where, 141 A l --(P / e-'l + V • A2 = (P/e - 1 + V ' a n d A3 = (P/e - 1 + y " (B.17) Equation (B.16) can then be expanded for thin shear flows to yield expressions for the stress components U ^ ' and U2- These are: A - 2 / 3 K 2 ^ 1 p 4 . 4 P ^ 2 (B.19) K 3 e ? + 3 e ' U 2 - 2/3 K 2 A l p + 4P ^2 (B.20> Substituting equations (B.19) and (B.20) in equation (B.18) yields 2 3U ^ = V A 3 ^ - ^ ( A , + A , ) ( 1 - A . | + 2 i L | ) • ( B " 2 I ) 1 2 e 3X2 3 e 1 2 / v 1 e 2 e dx £ or As 2 [3/2 M(P/£ - 1 + C,) + {(L - 2N)P/e + ^ - l } ( l - L - NJ] 3 ^ £ L — • ~ 2 K L fi - 3 X (P/e - 1 + C x) TT TT = - 9 OA U1 U2 3 e ZTj/c- _ i 4. r. ^ 2 (B.22) Further substitution for L, M, and N leads to equation (3.24). A Relation Connecting the Constants of the Dissipation Equation In order to derive this relation, use w i l l be made of equations 142 (3.32), (3.33), (3.34), (3.35), (3.36), and (3.37) of Chapter III. Replacing Tr- in equation (3.33) through equation (3.34) and K dy through (3.36) gives T 3 / 2 1 E = (1) -A- . (B.23) P <y Equation (3.30) for f u l l y developed constant property pipe flow can be written as e 2 1 Using equation (3.37) and noting that y = R - r, equation (B.23) can be written as .. •' T 3 / 2 , 3/2 Further, _3e 3r T 3/2 vpRy K 3 r 1/2 3/2 2 R - r (R - r)' (B.26) Using the results expressed in equations (B.26), (3.34), and (3.35), K 2 3e a'e 3r V a;c y b 2 + _ J 1 _ 2 R - r (B.27) and further, 143 Now putting R - r = y, and noting that r - R and further substituting the value of K from equation (3.35), e from (B.23) in equation (B.24) results i n T + 3y_ + 3y__ R R2 C )C e 2 u 3/2 (B.29) It should be pointed out that the entire treatment i n this section i s only valid for the f u l l y turbulent near wall region. APPENDIX C Computer Program 145 The computer program was compiled and executed on an IBM system 370, Model 168 computer under control of the MTS operating system. The Watfiv fortran compiler was used since the nature of the task undertaken required frequent debugging. However, the coding used in the program con-sists of commonly used fortran IV statements (with the exception of PRINT) and therefore the program should be compatible with other fortran IV com-pilers as well. The execution time for a run typical to the present work is about 40 seconds. The Main Segment The overall structure of the program is illustrated i n the adjoining flow chart (Figure 39). Control of operations is maintained through the main segment of the program. This section of the program con-tains a l l of the READ statements as well as most of the WRITE statements used in the program. Except for the subroutine VANDR, a l l other sub-routines are brought into operation only through this section of the program. For convenience of description and understanding, the main segment of the program has been divided into seven chapters. These divisions are marked on the computer l i s t i n g contained in Appendix D. However, i t was not possible to make s t r i c t demarcations of the main segment because of several small auxilliary calculations that are carried out in this section of the program. The approximate contents of the various chapters are described next. Chapter I. Dimensioning of arrays as well as statements containin'g input data. FIGURE 39 . Flow Chart of Computer Program. 148 CHAPTER II. READ statements for property data and conversion to appropriate units. CHAPTER III. I n i t i a l i s a t i o n of variables. CHAPTER IV. CALL statements for subroutines concerned with the starting profiles and Y-grid divisions. CHAPTER V. Call statements for subroutines concerned with solution of the d i f f e r e n t i a l equations and the pressure equation. Sub-routines used in assembling effective viscosities are also called in this section of the program. CHAPTER VI. Contains statements for imposing boundary conditions and the WRITE statements for dependent variables. CHAPTER_VII. Contains statements for rearrangement of array so that the computed variables do not exceed the allocated storage. Chapter V of the main segment u t i l i z e s a major portion of the total execution time (about 90%). Therefore, the programming approach for most of the other chapters gives preference to convenience rather than com-putational efficiency. For example, subroutine GRID i s called i n Chapter IV of the program. This subroutine i s actually required only for the Blasius f l a t plate problem. For other problems, i t i s superfluous. It was l e f t in because this subroutine was frequently used for testing the program. A l l computations are performed in SI units. Therefore, a l l inputs with the exception of physical property data are given in SI units. The appropriate conversion factors for the property data are included i n 149 the program. A description of the important aspects of the computer program is provided next. Storage Allocation A description of the arrays i s given in Table C l . The five dependent variables, U, V, T, DE, and ZK, as well as the properties RHP VISC, XK, and CP are stored i n arrays of size 10 x 100 (10 columns and 100 rows). The program w i l l also work without any change i n the logic of computation with arrays of a size as small as 3 x N, where N i s the number of grid divisions used i n a particular problem. Rather than changing N for every problem, 100 rows i n the arrays are provided in the DIMENSION statement since the Y-grid divisions rarely exceed this quantity. 100 rows imply 200 nodal points along each diameter of a pipe. A minimum of 3 columns are necessary i n this program for each of these arrays because values of variables from two previous steps are stored in the f i r s t and second column and the new computed values are stored in the third column. After the three columns are used up, variables from the second and third columns can be transferred into the f i r s t and second columns respectively. However, repeated replacements u t i l i z e more computational time. Therefore, a larger array size of 10 x 100 was used i n the program. In this way, computation proceeds up to the 10th streamwise step. Subsequently, the information from the 9th and 10th columns i s reassigned to 1st and 2nd columns of the arrays. Larger array sizes were not used i n order to keep the core usage within reasonable limits. 150 TABLE C .1. Description of Arrays Used in the Main Segment ARRAY DESCRIPTION CP(I,J) Specific heat C at node (I,J). '<• P DE(I,J) Dissipation e at node (I,J). I R(J) ' til Radial distance from axis of symmetry r, to J grid l i n e . RHO(I,J) Density p at node (I,J). RP(I,J) Production P at node (I,J). T(I,J) Enthalpy h (or temperature t) at node (I,J). TURV(I,J) Turbulent viscosity y t at node (I,J). U(I,J) Velocity U at node (I,J). UPL(J) ; + ,„ . , Tth Non-dimensional velocity U (location varies on the J grid l i n e ) . V(I,J) Velocity V at node (I,J). VISC(I,J) Viscosity u or V f f (varies) at node (I,J). XI (N) Temperature t on the n t h line of the property table. X2(N) til Viscosity ]X on the n line of the property table. X3(N) Specific heat C on the n line of the property table. P 1 X4(N) th •Thermal conductivity k on the n line of the property 1 tables. 1 X5(N) Density p on the n t h line of the property table. X7(N)* Enthalpies corresponding to XI(N). X(I) th Streamwise distance X to I grid. XK(I,J) Thermal conductivity k at node (I,J). XKE(I,J) 2 pK /e at node (I,J). .... continued *For arrays XI(N), X2(N), X3(N), X4(N), X5(N), and X7(N), the meaning given in this table should be regarded only as a preliminary meaning. During the course of computations, these w i l l normally be changed. 151 TABLE C l . Description of Arrays Used in the Main Segment (Continued) ARRAY XX(I) Y(J) YPL(J) ZK(I,J) DESCRIPTION Auxilliary variable containing predecided streamwise step divisions. th Distance from wall y to J grid l i n e . + th Non-dimensional distance y (location varies on J grid l i n e . Turbulence kinetic energy K at node (I,J). The remaining arrays are mostly one-dimensional and constitute a relatively minor proportion of the storage. Statements Defining Data and Operating Conditions Table C.2 describes the major input variables for the program. It w i l l be noted in Table C.2 that some input variables apply spe c i f i c a l l y to plane flow problems and others such as RADS are s p e c i f i -cally for pipe flow problems. If a pipe flow problem i s being solved, then the plane flow variables are assigned arbitrary values and vice versa. It w i l l be noted in the program l i s t i n g that a l l of the variables of Table C.2 (in Chapter I) are defined by ARITHMETIC statements. The same purpose can also be served through READ statements and this change may easily be introduced i f f e l t necessary by a future user. 152 TABLE Co2. Description of Input Variables VARIABLE CF GQ JEU JPRINT LNTP NC QW RADS SSTF TO TW UO UD DESCRIPTION If CF = 1, then computation scheme skips calculations for kinetic energy K and dissipation e. If CF = 0, then these calculations are performed. Mass flow rate. Number of Y-grid divisions. Number of streamwise steps that are skipped before output i s printed. If LNTP = 1, calculations are performed for turbulent case. If LNTP = 0 and CF = 1, then laminar flow computations are I performed. "lumber which determines the location Y(NC) where the multi-equation turbulence model takes over from the near wall algebraic formulations. Wall heat flux. Radius of pipe. Constant F^ for calculating streamwise step increments. I n i t i a l temperature. Constant wall temperature used for isothermal boundary condition problems. Free stream velocity. Use i n plane flow problems. Unheated distance on a f l a t plate for plane flow problems. The f i r s t executable statement of subroutines UVEL and PRES both assign a value to variable GRAV. If the flow direction i s v e r t i -cally upward, then GRAV i s assigned a value equal to -g. If buoyancy force i s to be neglected, then GRAV i s assigned a value equal to zero and i f the flow direction i s downward, then GRAV i s assigned a value equal to +g-153 The program as described solves for enthalpy as the dependent variable of the energy equation. If i t i s required that the temperature be the dependent variable, then the following statements contained i n Chapter II of the program are removed. TO = SLOP DO 110 N = 1,35 X7(N) = X1(N0 X1(N) = (25.0 + 5.0*N)*C7 XA(N) = X4(N)/X3(N) X3(N) =1.0 110 CONTINUED. Chapter IV of the computer program contains the following statement, R(J) = RADS - Y(J) . " This i s only applicable to axisymmetric flows. For plane flow problems, this i s changed to read, R(J) = 1.0. Another important change that w i l l be required for each problem i s i n the calculation of starting profiles. The program l i s t i n g s described i n Appendix D makes use of subroutine CHAMP for obtaining the Y-grid division and starting velocities. This subroutine i s applicable only for turbulent fully-developed flow starting conditions. For other problems, this sub-routine i s replaced by other problem-dependent subroutines. Another sub-routine CHANGE included in the program l i s t i n g applies to turbulent f l a t 154 plate problems. The program as described i n Appendix D i s set up for pipe flow problems where starting conditions require uniform enthalpy profiles and f u l l y developed velocity profiles. After understanding the working of the program, a prospective user should be able to make the necessary problem-dependent changes. Description of Subroutines A description of the subroutines Is provided i n Table C.3. A l l variables are handled as non-subscripted variables i n the subroutines U, V, E, and ZKIN. In order to simplify the notation used for different variables i n these subroutines, the following nomenclature i s used. A number which stands for the nodal location (as described i n Chapter II) i s attached to the variable name. For example, the velocity u w i l l be named around and at a central node as described in Figure 40. 155 TABLE C.3. Description of Subroutines pujjjMu^ Hinv'r 11111 ifrr TIT SUBROUTINE BCT2 CHAMP CHANGE GRID INITIA PRES PROP TURB UVEL VANDR WEL ZKIN ZZK DESCRIPTION Generates starting profiles of turbulence energy K and dissipation e. Generates the Y-grid and fully-developed starting velo-cit y profiles for pipe flow problems. Generates the Y-grid and starting velocity profiles for turbulent flow over a f l a t plate. Computes enthalpies or temperatures through the f i n i t e difference formulations for the energy equation. Generates the Y-grid for laminar flow over a f l a t plate, starting at the leading edge. I n i t i a l i z e s the arrays U(I,J), V(I,J), and T(I,J). Computes pressures through the mass flow constraint. Computes properties by linear interpolation of values i n the property tables. Computes and assembles effective viscosities and effec-tive conductivities through algebraic formulations. Computes velocities U through the f i n i t e difference for-mulation for momentum equation. Computes velocity gradients in one-dimensional flows through the use of Van Driest hypothesis. Computes velocities V through the f i n i t e difference formulations for the continuity equation. Computes turbulence energy k and dissipation e through the f i n i t e difference formulations of their respective equations. Computes and assembles effective viscosities and effective conductivities according to the K - £ turbulence model. 1 5 6 1-1 I 1+1 FIGURE 40. Nomenclature for Variables in Subroutines. 157 APPENDIX D Listing of Computer Program LISTING OF ASHN+... COMMON INDEX DIMENSION XPC100),X{100),Y{ 100) ,U<10, 100) ,VU0,100) , 1RH0U0,100)tVISCCI Of 100)tR(100) 1*XX1300) It XKtl 0,100) , CPU 0,100) ,T(10,100 ) 1,TURVI 10,100) DIMENSION X l O 5 ) , X2(35),X3{35),X4(35) , X5 ( 35 ) , X7( 35) DIMENSION YPL MOO) , UPLl 100 ) DIMENSION XKEMO, 100) ,ZK{ 10, 100 ) , RP U 0, 1 00) ,DE( 10,100) JPRINT=10 JPR=JPRINT T0=292.55 PI=3.1416 XXXP=7$10.0**6.0 RADS=.01 G0=.2 QW=5000 TW=TO UO =.45 UD=.i NC=35 CF=0.0 JEU=50 DELF=1 SSTF=0.5 DELTA=RADS LNTP=1 CHAPTER 2 INDEX=2 XP(1)=XXXP LISTING OF ASHN+ XPC 2)=XXXP Cl=5,0/9.0 C2=1.4882 C3=1000.0/0.23884 C4=3600.0/0.5778 C5=l.60184*10.0 C 7=1.89/1356.0/6.0*10.0**7.0 PEAK=-3.0 PEAK=3.0 00 102 N=i,35 READ' 5f 101) Xl(N),X2{Ni,X3tN),X4lN),X5(N) 101 FORMAT(5G9.4) X1(N)=C1*XKNI X2(NJ=C2*X2(N) X3(N)=C3*X3(N) X4IN)=C4*X4tNI X5(N}=C5*X5<N) 102 CONTINUE IF(PEAK.GE «0»0 I GO TO 104 DO 103 N=ll,22 READt5,101) X4IN) X4(N) = C4*X4(N) 103 CONTINUE 104 CONTINUE DO 109 N=i,35 X7(N)=(25.0+5.0*N)*C7 109 CONTINUE CALL PROP{T0,XXXP,SRHO,SVISC,SXK,SCP«SLOP,X1,X2,X3.X4,X5,X7) TO=SLOP DO 110 N=l,35 Ln LISTING OF ASHN+... X7(N)=X1(N) XI(N)={25.0+5.0*N)*C7 X4{N)=X4(N)/X3{N) X3IN)=1.0 110 CONTINUE CHAPTER 3 PRT=SCP*SVISC/SXK PRINT,PRT CALL INITIA(U,V,T,U0,T0) 00 12 1=1,10 DO 12 J=l, 100 XKEII,J)=0.0 ZKU,J>=0.0 RP(I,J> = 0.0 DEU,J)=0.0 TURViI,Ji=0.0 RHOd, J)=SRHO VISCC I,J) = SVISC CPU,J)=SCP XK(I,J)=SXK X K{I,J)= SXK/SC P C P U , J » = 1.0 12 CONTINUE W=25.0*VISC«lt ll/UO/RHOUtl ) xxm=o.o FACTOR=0.37 XX(2)=FACT0R*W DO 152 1=2, 299 XX{ 1 + 1 )=XXU ) + 5.0*FACT0R/f UO*RHOU , 1)/V I SC t 1, 1) /XX (I ) )**0.5 152 CONTINUE M ON o L I S T I N G OF ASHN+ >. DO 151 1=1,10 X i i ) = X X l I ) 151 CONTINUE N=JEU CHAPTER 4 CALL GRIDCWtNtYl C A L L CHAMP(RAOS,NtYfGQ.SRHO-SVISC»U,SHMAl 1 = 1 CALL TURB(U,VISC,RHO,Y,TURV,I,N,CP,XK,LNTP,DEITA,SHWA) CALL BCT24U ,Y » DE ,ZK , NC » N, SHWA.SRHO) DO 1500 J = l , 9 9 XK{ I, J }=XK ( I, J J-TURVU , J ) * C P t I, J) 1500 V I S C U , J ) = V I S C U , J } - T U R V U , J J X ( U = X X < i l X(2 )=XX<2) STF=Y( 2)/RADS DEL TA=RADS DO 13 J=1,N R{ J ) = RADS-Y(J> 13 CONTINUE C 2=1.9 - . 16/.3*(1.+3.*Y(NC)*Y{NC)/R< NC)/RINC) 1*3.*Y{NC >/R'NC)> CHAPTER 5 DO 154 L=lt70 DO 15 1=2,9 I F C C F . G T . 0 . 5 ) GO TO 77 7 CALL ZZKCXKE,ZK , DE , RP ,U , VIS C, RHQ, Y , TU RV 1, I »N, N C R ADS, CP, XK, SHWA.CF ) GO TO 775 777 CALL TURB(U,VISC,RHO,Y,TURV,I,N,CP,XK,LNTP,DELTA,SHWA) LISTING OF ASHN+.. 775 PRINT,SHWA XII+1) = STF*DELTA+XU) CALL PRESJ U,V,VISC,RHO,R,Y,GQ,N,I,X(I-1),X{I),X<I+U, lXPCI -n ,XP( 1+1 ) ,SHWA) NN=N-1 DO 14 J=2,NN CALL U V E L ( U < I , J + 1 ) , U ( I - 1 , J ) » U ( I , J ) , U ( I +1, J ) , U ( I , J-1 ) • iRHoUt j» tvisci i, j+i i, vi sea , J I , vi sen, j-n ,vi I .J I , 1R1J+1),RIJ)»R{J-l)» Y(J + l)» Y{J),Y(J-1) 1,XP( I - H.XPC I) .XPU+11 , X11 -1J ,X(I ) ,X{ 1 + 1) } CALL E(U{ I , J-1 I ,U< I , J > ,U(I , J - r l l ,TII , J - i ) ,TU , J) ,TI I , J + l ) , 1TC1-1,J),T(1+1,J),XP(I),XP(1-1),XP(1 + 1 ) , Y { J - i ) , Y ( J ) , Y { J + l ) , R { J - l ) 1RU)»R{J+1) , R H O ( I , J I , C P ( I , J ) , V I S C ( I , J ) , X K { I , J - l ) , X K U , J ) , X K ( I , J ^ l 1, X ( I ), X ( I-1),X< I + 1),V( N J ) ) CALL PROP (T{ 1 + 1 , J) ,XPM+1 ) , RHO {I +1 , J ),VISCU +1, J ) ,XK( 1 + 1, J ) , 1CP (1+1,J),SL0P,X1,X2,X3,X4,X5,X7) CALL VVELIUC1-1, J) ,U ( I , J ), U U +1, J ) »U ( 1-1, J - 1) ,U( I , J - l ) , U U + 1 , J - l ) 1, RHOl I -1 , J) » RHO ( I, J ) , RHOC I+1, J I ,RHOH-l, J - l ), RHOt I » J-1)»R HO{I+ 1, 1J-1),R( J+l),Ri J),R( J - l ) ,VM + 1,JI ,V( 1+1, J - l ) ,X(I-1 I ,X( I) ,XU + 1) , 1 Y ( J + l l , Y I J ) T Y ( J - l ) ] IF(J.LE.NC) GO TO 105 IF(CF.GE.0,5) GO TO 105 CALL ZKIN { ZK( I, J - i ) , ZK{ I, J ) , ZK{ I , J + l ), ZK 11-1 , J ) , Z KU+i , J ) ,DE( IIrJ-l),DE{ I,J),DE(ItJ^l)»DE(I-l,J)» DE(I+1» J),XKE(I,J-l),XKE*I , J ) , 1XKE{ I, J+l ) ,Y( J-1),Y {J) ,Y{ J + l ) ,R(J-1 ), R( J) ,R C J + l ) , X{ 1-1) , IX ( I I , XI 1+1) ,RHOJI, J) ,V(I , J J ,RPU, J ) ,U( I , J ),C2 ) 105 CONTINUE 14 CONTINUE CHAPTER 6 IF(STF.LE.SSTF)STF=STF#1.2 L I S T I N G O F A S HM + • • 1 0 J C O N T I N U E J C = 3 I P 1 = I + I DO 1 0 0 0 J = l , 9 9 X M I , J I = X K ( I # J i - C P { I , J I * T U R V C I * J ) 1 0 0 0 V I S C ( I t J ) = V I S C ( I , J ) - T U R V { I , J ) T U + l » n = T { I + l t 2 ) + 0 W * Y ( 2 } / X K { I , l l C A L L P R O P ' T t 1 + 1,1 ) , X P < 1 + 1 ) , R H Q { I + 1, I ) , V I S C U + 1, 1) , X K { I + 1, 1) , C P 1( 1 + 1 , 1 ) , S L O P , X I , X 2 , X 3 , X 4 , X 5 , X 7 ) D I A M = X ( I - l l / R A D S / 2 . 0 U { I + i , N ) = U U + l , N N ) V ( 1 + 1 , N ) = 0 . 0 R H O ( I * i , N >=RHO< I * i , N N > V I S C < I + l , N ) = V I S C U + i , N N ) X K { 1 + 1 , N ) = X K ( I + 1 , M N ) C P ( I * 1 , N I = CP< t + I . N N ) T ( I + l , N ) = T ( I + 1 , N N ) E N T H = 0 . 0 D O 2 0 0 0 K = 1 , N N E N T H = E N T H + { T ( I + l , K ) * U ( I + i , K ) * R H Q < I + 1 » K ) + T ' I + 1 , K + 1 ) * U f I + 1•K+1)* 1RH0 < I - 1 , K * l » J l M P ( K ) + R ( K + l ) ) * ( Y < K + ] * - Y ( K ) ) * P I / 2 . 0 2 0 0 J C O N T I N U E E N T H = E N T H / G Q W R I T E ( 6 , 7 5 ) D I A M , S L O P , E N T H , S H W A 1J F 0 R M A T { 4 X , ' 0 I A M S = « , F 6 . 2 , » W A L L T E M P = ' , P 7 . 3 , • B U L K E N T H . = • , G 11 . 4 , i « S H E A R S T . = ' , F 6 . 2 > I N 0 € X = 2 J P R = J P R + 1 I F ( J P R . L T . J P R I N T " G O T O 7 7 9 ^ O N L 1 S T i N G O F A S H N * . . . W R I T E ( 6 , 7 0 ) X I I ) W R I T E ( 6 , 7 1 ) W R I T E ( 5 , 112 I f U ( I f K l , K = 1 , N ) W R I T F ( 6 , 7 2 ) W R I T E ( 5 , l i 2 ) ( T ( I , K ) , K = 1 , N ) W R I T E ( 6 , 7 3 ) WRITEC5 ,112» I Z K ( I , K I , K = N C , N I W R I T F ( 6 , 7 4 ) W R I T E ( 5 , 1 1 2 ) ( D E ( I » K ) , K= N C » N ) W R I T F ( 5 , 1 L2I( T U R V ( 1 , K ) , K = l , N l 7 u F Q R M A T J 4 X , * X = ' » F 1 1 . 4 ) 7 i F 0 R M A T ( 4 X , ' V E L O C I T I E S ' ) 7<i F O R M A T ( 4 X , « E N T H A L P Y " ) 7J> F O R M A T ( 4 X K I N . E N E R ) 7-* F O R M A T ( 4 X , « 0 I S P . ' ) J P R = 0 7 7 * C O N T I N U E lo C O N T I N U E C H A P T E R 7 D O 1 8 1=1,2 X P ( I } = X P ( 1 + 8 ) DO 1 8 J = 1 , N Z K H , J ) = Z K ( I + 3 , J ) D E ( I ,Jl=OE< 1 + 8 , J l X K E I I , J ) = X K E ( 1 + 8 , J ) R P ( I , J ) = R P ( 1 + 8 , J ) UU • J » = U ( 1 * 8 , J I V ( I , J ) = V ( I + 8 , J ) R H O I I , J ) = R H 0 ( I + 8 , J ) V I S C ( I , J ) = V I S C ( I - « , J » L I S T i N G O F A S H N . . . X K ( i , J ) = X K U + 8 , J ) X ( I >=X{ 1 + 8 > C P ( I , J ) = C P ( T +8 , J ) T ( I » J ) = T ( I • 8 » J ) l b C O N T I N U E P R I N T , X ( 1 + 1 ) 1 1 ^ F O R M A T { I X , L O G l 1 . 4 ) 1 5 - t C O N T I N U E 77<i C O N T I N U E S T O P E N D S U B R O U T I N E Z K l N ( Z K l , Z K 2 , Z K 3 , Z K 4 , Z K 5 , D E i , D E 2 , D E 3 , 0 E 4 , D E 5 , 1 X K E 1 , X K E 2 , X K F 3 , Y 1 , Y 2 , Y 3 , R 1 , R 2 , R 3 , X 4 , X 2 , X 5 , R H 0 2 , V 2 , K P 2 , U 2 , 1 C 2 ) SIGK=1.0/0.09 Z Z 1 = 0 . 5 / R 2 / { Y 3 - Y 1 ) * { R 3 - R 2 » / ( Y 3 - Y 2 1 Z Z 2 = 0 . 5 / R 2 / { Y 3 - Y I ) * ( R 1 + R 2 ) / ( Y 2 - Y 1 ) Z l = Z Z i * ( X K E 2 + X K E 3 ) / S I G K Z 2 = Z Z 2 * ( XKE1 + X K E 2 I / S I G K C0EF=RH02*U2/(X5-X4 ) + 0 . 5 * { Z 1 + Z 2 ) T E R M = R P 2 - D E 2 * R H 0 2 * Z 1 * { Z K 3 - 0 . 5 * Z K 4 ) - Z 2 * { 0 . 5 * Z K 4 - I K I ) 1 - R H 0 2 * V 2 * ( Z K 3 - Z K 1 ) / { Y 3 - Y 1 ) +RH*)2*U2*ZK 4 / ( X 5 - X 4 ) ZK 5 = T F R M / C . O E F C A S S U M I N G Z , S A N D C Q E F A R E S A M E T E R M i = R P 2 * C 2 * Q E 2 / Z K 2 - 1 . 9 * R H Q ^ 0 E 2 * D E 2 11 K 2 1 - Z 1 * ( D E 3 - 0 . 5 * D F 4 ) - Z 2 * ( G . 5 * 0 E 4 - Q E 1 ) - R H 0 2 * V 2 * 11) E 3 - 0 E l I / t Y 3 - Y I ) l + R H C 2 * U 2 * D E 4 / { X 5 - X 4 ) D E 5 = T E R M 1 / C 0 E F R E T U R N E N D ON Ul L I S T i N G O F A S H N + . . S U B R O U T I N E B C T 2 ( U , Y , D E , Z K , N C , N , S H WA» S R H O I D I M E N S I O N U ( 1 0 , L O O ) , Y ( 1 0 0 ) , Z K ( 1 0 , 1 0 0 ) , D E ( 1 0 , L O O ) N O N C - 1 0 U S T = S H W A / S R H O DU=<U( l . N C + 1 ) - U { 1 , N C - 1 ) ) / ( Y { N C + l ) - Y ( N C - l ) ) U K l = . l 6 / . 3 * Y < N C ) * Y C N C ) * C U * D U U K 2 = 0 . 8 4 0 * U S T U 2 = U ( 1 , N ) U l = U ( l . N C > C 1 = I U K 2 * * L . O - U K 1 * * 1 . 0 I / ( U 2 - U 1 ) C 2 * U K 1 * * 1 . 0 - U l * C l N N = N - 1 D O 1 0 J = N C , ! M N 0 U = ( U J 1 » J + 1 ) - U ( 1 . J - l I I / ( Y ( J * l ) - Y { J - l ) ) Z K ( 1 , J ) = ( C l * U ( 1 , J l * C 2 1 * * 1 . 0 Z K ( 2 , J ) = Z K ( 1 , J ) D E ( 1 , J J = 0 . 3 * Z K ( 1 , J ) * n u 0 E ( 2 , J ) = D E I 1 , J ) U C O N T I N U E N C = N C + 1 0 R E T U R N E N D S U B R O U T I N E Z Z K ( X K E , Z K , O E , R P , U , V I S C , R H O , Y , T U R V 1 , I , N , N C , R A D S , C P , X K . S H W A , X X ) D I M E N S I CN X K E ( 1 0 , 1 0 0 ) , Z K { 1 0 , 1 0 0 ) , D E ( 1 0 , 1 0 0 ) , P P ( 1 0 , 1 0 0 ) , 1 U ( 1 0 , 1 0 0 ) , V I S C ( 1 0 , 1 0 0 ) , R H O ( 1 0 , 1 0 0 I , T U 4 V ( 1 0 , 1 3 0 » , C P ( 1 0 , I 0 0 » , I X K { 1 0 , 1 0 0 ) , Y ( 1 0 0 ) , Y P L ( I C O ) , U P L I 1 0 0 ) D E ( I , N ) = 0 E ( I , N - l ) Z K ( I , N ) = Z K ( 1 , N - 1 ) S H W A = V!SC( t , I ) *U( I » 2 ) / Y ( 2 ) L I S T I N G D F AS H N + * . N C F = 1 0 U S T = ( S H W A / R H O C 1 , 1 ) ) * * 0 . 5 O O 1 J = 1 , N C Y P L ( J ) = Y ( J ) * L ! S T / V I S C U , 1 ) * R H O ( I , 1 ) U P L ( J ) = U ( I , J ) / U S T 1. C O N T I N U E DO 2 J = 2 , N C R P L = 0 . 4 * Y ( J I I F I Y P L ? J ) . L E . 1 2 5 . 0 ) R P L = R P L * (1 . O - E X P ( - Y P L ( J ) / 2 8 - 0 ) ) T U R V U , J ) = P H O < I , J I * R P L * * 2 . 0 * ( U ( I , J « - l • - U i I . J - l ) l / i Y < J - L l - Y U - l > ) I F ( T U R V ( I , J ) . L E . 0 . 0 ) T U R V ( I , J ) = - T U R V { I , J ) 2. C O N T I N U E Z K W A = 3 , 3 3 3 * U S T * * 2 . 0 o N C F = N C F + l I F C Y P L ( N C F ) . L E . 7 5 . 0 1 G O T O 6 ! ) U = « U ( I , N C F + 1 ) - I J ( I , N C F - l ) ) / { Y < N C F + 1 ) - Y < N C F - 1 ) ) Z K.( I , N C F ) = , 1 6 / . 3 * Y ( N C F ) * Y ( N C F ) * D U * D U 0 E { I , N C F | = . 3 * Z K ( I , N C F ) * C U D E I I + 1 , N C > = D E ( I , N C ) Z K C I + l , N C ) = Z K ( I , N C ) N C C = N C + 1 D O 1 0 J = K C C , N X K F ( I , J » = R H O < I , J ) * Z K ( I , J ) * Z K ( I , J ) / D E C I , J I 1 U C O N T I N U E X K E J I , M C ) = P H 0 ( I , N C F ) * Z K ( I , N C F ) « 7.K 1 I , N C F ) / 9 F ( I , N C F ) F A C T = ( Y C N C C ) - Y { N C > I / ( Y < N C C l - Y ( N C F ) ) Z K { I , N C ) = Z K ( I , N C C ) - F A C T * ( Z K { I , N C C ) - Z K { I , N C F ) ) D E I I , N C ) = D E ( I , N C C ) - F A C T * ( D E C I , N C C ) - D E ( I , N C F ) ) N N = N - 1 D O 3 5 J=NCC,NN L I S T i N G O F A S H N + . . . R P t I , J t = X K E ( I, J » * ( U ( I , J * l » - U ( I , J - l l l * ( U ( I . J * l ) -U { [ , J - l ) i l / ( Y C J + l ) - Y ( J - l ) ) / ( Y ( J * l ) - Y ( J - l ) 1 * 0 . 0 9 T U R V { I , J ) = Q . 0 9 * X K 5 C I , J ) 3o C O N T I N U E T U R V l I , N ) =T .URV « I »t\iN ) N N N = N - 2 D O 3 J = 2 . N N N j> T U R V ( I » J ) =( T U R V U . J - i ) + T U R V ( I , J )+ T U R V ( I , J +1 H / 3 .0 P R T = 0 . 9 D O 4 J = 1,!M V I S C ( I t J ) = V I S C C I . J ) + T U R V ( I , J » X K { I , J > = X K { I , J » - T U R V U , J » * C P C I , J l / P R T * C O N T I N U E R E T U R N E N D 00 169 i r o * a -> ro U o o > eg > + a cv I C J cc OO I 2" * > f* II •« • ««. sj" — i d « • a o i » r H i — t o > r - l eg a: r - l X o j E f* o co > Z fl •—i < t u 2" »• *—< •* fl » X • r H ' w O i n UJ » r - l • r - r-t • II < Z> t o "~* • 0 0 Z o II • a t 1 u. II Z — i - t UJ H-<T IX. a UJ <••> z II r - l — t - t ~ o <r X II u *—* Z r - Ot rs »—i w O w o at: O o OO Q > LL c > O 3 LL og o oo m u oo — >-. i n > x CV eg O X X » a: ^ • x r-1 •» D un » a L O x ™o, * * eg O J CL ZZ X m> *> •4- ~t z> a • X f O * —< — >. _ J • ILI eg > > a ; > x t u o IT (Ni > CD eg > # > — I co c\j a: > + •— eg cc — * t _ og • | * CM ro —« > > I »r* > — L O i-« (NI ZZ + I ro f\j z> x w I II i n « ZZ X eg a •<-» <x LU N. + 0 0 og O oo (.5 H- or 3 * a t ce x> m • Co I II > < a: o pg a C L ' O X I . I f, O LTi eg • a > r - X I U-r g og — > IT. f • CC ~* U UJ ng h - ZD r g + (N! c X a -\ eg 3 fNl X I LO X *• pg r~ > r~ a : U J o t— i — i l un o •—» og i ~ t —< > > I I > 1 1 t r o > >- > *—» *m** » N . r o og ft o e t o oo- 1 >~« 1—1 r o > > r> + + rvt fi - *• o o M og oO t o N J > ^ » f fi { ¥. r - t > > > ft ."NI n> <J eg c i a — X <r C + ac Z3 *• CNJ rvi WO | •« r v oC cc • — ' LO CJ + + O • X r o ft + X O CL sC ac 1 w + > t LO -tt X X Og 1 — r»g X ft ft L O " V t 1 X — — t o 1 1 <t •4- X m r o \ a ' > >-r\j og ID a t of 3 >«» ^ 2^  Lf> LO r -« * t o o z II II O <-* eg O N J N J ~> LO og x ' * # on eg a C x X -f a- I il II u s; U J a: O LU O K o I * N . LO * rsj ro- rs eg * O X NJ CC 0 0 - J L I S T I N G O F . . . + A S H 2 + . . . U 5 = T E R M / C O E F 7 7 C O N T I N U E R E T U R N E N D S U B R O U T I N E V V E L i U 4 , U 2 » U 5 . U 7 » U 1 » U 6 , R H 0 4 , R H 0 2 , R H 0 5 , R H 0 7 » R H U 1 , 1 R H 0 6 . R 3 » R 2 » R l , V 5 , V 6 , X 4 , X 2 1 X 5 , Y 3 , Y 2 , Y 1 ) C O M M O N I N D E X I F ( I N D E X . G T . 1 ) GO T O 50 T E R M = - { R 2 + R 1 ) * ' R H 0 5 * U 5 - R H 0 2 * U 2 + R H 0 6 * U 6 - R H O I * U 1 ) / ( X 5 - X 2 ) T E R M = T E R M / 4 . 0 V 5 = ( T E R W * { Y 2 - Y 1 ) + R H 0 6 * V 6 * R i ) / R H Q 5 / R 2 G O TO 7 8 5w C O N T I N U F T E R M = { R 1 + R 2 ) / 4 . 0 / ( X 5 - X 4 ) * { R H 0 5 * U 5 - R H Q 4 * U 4 * R H 0 6 * U & - R H 9 7 * U 7 ) T E R M = - T E R M V 5 = ( T E R M * < Y 2 - Y 1 l - R H 0 6 * V 6 * R l J / R H 0 5 / R 2 7 o C O N T I N U E R E T U R N E N D S U B R O U T I N E E { U 1 , U 2 , U 3 , T I , T 2 • T 3 , T 4 , T 5 , X P 2 , X ? 4 t X P 5 , Y l , Y 2 , Y 3 , l R l f R 2 t R 3 t R H C 2 t C P 2 , V I S C 2 . X K l , X K 2 t X K 3 , X 2 , X * t X 5 , V 2 ) C O M M O N I N D E X C J = l . 0 I F ( I N D E X . G T . 1 ) G O TP 5 0 T E R M = U 2 * ( X P 5 - X P 2 l / ( X 5 - X 2 ) * C J - 0 . 5 / R 2 / ( Y 3 - Y U M ( R 2 + R3 ) * < X K 2 +X K 1 * I T 3 - T 2 ) / < Y 3 - Y 2 ) ~ ( R 1 + R 2 ) * { X K 1 + X K 2 ) * ( T 2 - T 1 )/ ( Y 2 - Y 1 I ) i - V l S C 2 * ( { U 2 - U 1 I / ( Y 2 - Y L » > * * 2 . 0 - R H O 2 * V 2 * C P 2 * < T 2 - T 1 » / < Y 2 - Y 1 ) T 5 = T E R M * { X 5 - X 2 ) / R H 0 2 / U 2 / C P 2 + T 2 G O T O 7 9 5 u C O N T I N U F 171 eg X CO < z og — i > > I I m eg > > r\j — i id X x + + rO PO x x a: cn + • C O O J * * a : ft » > > I I r o r o eg eg of a; CO LO o o II I! ft og og eg -fr —' - J > o g i + C O > r~« — r > - v I — ro —4 > > -v CO — > I — fO i r-l t— x — I # co eg X a — O * «• — r * ^ C \ J r—I eg > X og -» | + eg m —* G ZD *J X — —> a t *• # I eg LO ~~) C_) • O CO O * » -+ — > — v f + ~$ X —> X I — I I LO t~ LO X I •v v~ — Kr eg -j- to a a • u x o * i — eg to ¥ O Q . N x eg * eg O X a x> II Ii U X LU CC a L U eg eg u_ a. UJ c_> O * c_ eg ZD X -J? a eg LU O r -X II a LOI + Ui Z 2 cc r — X LO X v f X C O X »> eg X r - i m X r o a r~ a x - j • co —• • t o a , co o ~-CO to - X V •> X — I/O LO » r o u —• CO Nt — ' X > * 00 — -- L O C co 1 *-PL ro CO X —- »• CC ~ X IO CO £ eg C X X » a t o a r o ct. — C L r i X u i z z i — o X I CO a 2 Ct UJ o o o eg o L O C O o Ui ro a. a r - C X + w I — og w II u. II C L a I— O 0 t l X l u 0 Q » - , » - , » - i » - , O r—. r-i r~-4 •WW yf DOT LO eg CO r~ X X X — X c + f + + — i f* >*-X i — f + o f* H C X 1 r-» - c C — i rr« rr —» —« X -r-l * o + f* eg t—4 M ft hr- r l It-w. tO- X W w > > r-l X 1 CO gj- X X 1 f X X l X r*. fl 1 1 f — o r<4 + — f. fi <T O + f—1 —I fi + i—< fi — rr f + + H - • f* w eg r-. fi ~r M O "-rr co X r~ Z — i r - t X or ro X r - i c X — * X X r~l 1 » — *rrf w UJ r 5. r- z z r-l C" a < I - <_> — c~ II X < X u o <r )— r-i tr-4 LL II <r < X Z X to, —1 || II c X X II or. C •£ o e,i | h-o C O II II a X) QC LU t- r - l t-i O X • r - l a o t— II II < aC > L J X _ J UJ Z X »-< IS o C O C O C O C O C O c t UJ co O |r. o o~ L I S T I N G O F . . . + A S H 2 * . D O 12 J = l , 1 0 0 U I I t J ) = U O V ( I , J 1 = 0 . 0 T l I , J ) = T C i^ C O N T I N U F D O 1 4 1 = 1 , 1 0 T { I ,1 ) =TW U I I , 1 1 = 0 . 0 1 4 C O N T I N U E R E T U R N E N D S U B R O U T I N E T UR B { U , V I S C , R H O , Y , T ' J R V , I » N »C P , X K , L N T P , D E L T A , SHWA I D I M E N S I O N U i n . l O J I , V I S C < 1 0 , 1 0 3 ) , R H Q U O . 1 0 0 ) t T U R V < 1 0 , 100) , I C P U O , 1 0 0 » , X K I 1 0 , 1 0 0 1 , Y ( 1 0 0 » , Y P H 1 0 0 « , U P L ( 1 0 0 I S H W A = V I S C { I , 1 ) * U ( I , 2 ) / Y < 2 ) U S T = ( S H W A / R H O i I , 1 ) ) * * 0 . 5 D O 1 J = 1 , N Y P L ( J ) = Y { J ) * U S T / V I S C ( I , L ) * R H O ( I , i ) U P L ( J » = U< I t J l / U S T 1 C O N T I N U E N N = N - 1 DO 2 J = 2 , N N R P L = 0 . 4 * Y ( J ) I F ( Y P L ! J ) . L E . 1 0 0 . 0 ) R P L = P P L * ( I . 0 - F X P t - Y P L ( J ) / 2 8 . 3 ) ) I F C R P L . G E . ( 0 . 0 8 9 * O E L T A ) I R P L = 0 . 0 8 9 * 0 E L T A T U R V ( I , J ) = R H O { I , J ) * R P L * * 2 . 0 * ( U { I » J + I ) - U { I , J - l ) ) / ( Y ( J + l ) - Y ( J - l ) ) I F ( T U R V ( I , J > . L E . 0 . 0 i T U R V ( I , J l = - T U R V ( I , J l 2 C O N T I N U E 3 * N N N = N - 2 DO 3 J=2,NNN 173 L I S T i N G O f . . . - A S H 2 * . . . Y { 2 ) = S V I S C / S R H O / I SH W A / S R H O ) * * 0 . 5 D Y M = { 0 . 9 9 * R A D S / Y ( 2 ) I * * < 1 . 0 / ( N - 3 I ) P R I N T , D Y M U T { 1 ) = 0 . 0 N N = N - 1 D O 1 J = 3 , N N 1 Y ( J ) = DY M * Y I J - l ) DO 2 J = 1 , N U ( 1 » J ) =U M A X * { Y ( J ) / R A D S ) * * ( I . O / S N ) U ( 2 , J ) = U C 1 , J ) Z C O N T I N U F D O 3 J = 2 , N N Y D = ( Y ( J ) + Y { J - l ) ) / 2 . 0 S H= SHW A * { R A D S - Y ( J ) J / R A D S C A L L V A N D R . S H . Y B . S V I S C , S R H - 3 , D U > U T ( J I = U T { J - I I - 1 Y ( J ) - Y ( J - l ) ) * D U a C O N T I N U E D O 4 J = 2 , N N I F I U T C J J . L E . U f l . J I I U { 1 , J J = U T { J ) U I 2 , J ) = U ( 1 , J ) * C C N T I N U E P R I N T , S H W A , R E Y N . U M A X t U B A R W R I T E ( 5 , 1 1 2 ) ( U ( 2 , K ) , K = 1 , N ) W R I T E ( 5 , 1 1 2 » ( Y ( K I , K = l , N ) 1 1 * . F O R M A T ( I X , 1 0 G 1 1 . 4 ) D O 6 J = l , N o R U ) = R A O S - Y ( J J U ( 1 , N ) = U ( 1 , N N ) U I 2 » N ) = U ( 2 , N N ) S U M = 0 . 0 M 4> lib o x Q CO > I > X oo •> LO a x -* a x •> LT> IO X 4»-C \ J CM X * -* X o o o o o X — s: or _J + — — •— UJ — —> —* X tf + x • I ai — ~3 J " T C ? -ULI » + x a o •-! a u-' — jj i / l w x Q » a i— o a. » ~ a, + i—. — x ~3 O ~ > # m » •V —• 4 -s- O — i — + — _ j w I D. ~J u to U J Ci + a. _j u i * — + a x U J 1 o CD — - 3 4> 4—4 f—4 U J 1—4 44 X ) r- r-H — a 4—4 — » 1 —. a Z o + o _J L O ~3 l —4, Q. CO CO L U • o 4—4 o _J I—1 r~> D n o CO 1 cc —> o o L U > > 4f X 1 4—4 - 3 « w 4* to c (X > r ^ > 4* # r-4 ~ r - CC ~3 4> _ > — — + * + •> # X a. 4* • - 3 a — r—4 , — ^ •—4 4r — j - 3 " — ' -» c -» o —• 1 • 1—4 L U — 1 O o r-4 X o o lY. "7 -5 — n Cf. X X 4 ex. o 4—4 + «t<r» #• X5 at cr ~ > 4* I—1 — * or r»H r-4 KJ »r4" o \ w I * X 4—» (XI CO » — c -5 X — » - 3 - 3 •—»- 4—4 a r-< •• U J r •» ID > •> x> o w cc i—• (~l » 4—4 4—4 4 o D or * —<• + — — " > o 4—4 M — 1—4 a 1 o o o o _ — » •K 1 — X r-' X X X J .-4 o CM a U J - — ' cc DC • <\l — 4 L L •> -J _) r-4 r~4 t O o X ii • to cc < r-4 U J U J — t 4—4 — #. o o > 4 — ^ U J a o — — 1 O o o > * IT. • • <r (M QC 4—4 o > -5 c - 4f X + • I # a. ~ o <t (VI 1 I <•*. QC. 1—4 ~3 O + CO z m o > <r <\l —. >" 5. W r-4 » 1 — <r z • U J o — . X X 1 * X X o X 1 + 4—1 — . _ j - 3 4 •> o o z 4—4 1 1 + o U J U J • - 3 W 4—1 w L U «-r' _J • ~ » t a o i— . . c I L • 4—4 CM L P • o •> C 1 1 (£- O CD U.' • U £ 4r r - 4-4 <* f\l L O |i •— ~- + 4 1 — 4—4 X - 3 a a c • "J 2 * r> a to — 1 r-4 X X " 7 > II a. Q. II QC •r II n * — \ <X a z o II • II II I II II X X » o —1 a. a tn II 00 z Z3 cc U J o > X Q L z OL s: U J U J - > i—• CO W - 3 III _J o II 4—4 t - C cc •3L 4—4 < II X X II —J _J c a — 4—4 • — ' o «-» U J a => o a U J z r> U J 4—4 4—4 L U U J z o L U U J —' C D > > CD QC c <J5 O co O a U J CO CC o > o a. Ci o z c D C a + 4 a * + * a •44» C L # 2 —1 •—4 •—4 f-4 4—4 r-4 4—1 r-4 L I S T I N G O F . ..+ASH2+... P C ( J ) = - 4 . 0 * R H O ( I f J ) * P ( J ) * t D E L P + O E L M ) * D E L P * D E L M 1 C O N T I N U E P A { 1 1=0.0 P B ( 1 J-0.0 PC(11=0.0 P A ( N . = P A ( N N I P B < N ) = P B { N N ) P C < N ) = P C ( N N ) F H 1) =0.0 F2(11= 0.0 DO 2 J=2 , N P l ! J i = P P ( J ) / P A ( J ) F 2 { J ) = P C ( J ) / P A { J ) I C O N T I N U E SUM 1=0 .0 SUM2=0.0 D O 3 J=1 ,N\J S U M 1 = S U M 1 * ( F 1 I J J * F 1 ( J * 1 M * ( R ( J » * R ( J - 1 H * ( Y 0 1 ) - Y ( J M * P I / 2 . 0 SUM2 = SUM2*(F2( J ) + F 2 ( J + l » ) * < P U ) * R ( J + L ) ) * < Y ( J * l ) - Y C J ) ) * P I / 2 . 0 3 C O N T I N U E DELPR=(GC-SUM1 I/SUM2 XP5=DELPR+XP4 R E T U R N E N D ON L 1 S T i N G O F . , . * * S O U R C E * S U B R O U T I N E C H A N G E ( S R H O » S V I S C , U D , U O . Y , X X . U , V , N , S H W A , D E L T A ) 0 I M E N S I ON X X I 3 0 0 ) , Y ( 1 0 0 ) , U { 1 0 , l O O ) . V ( 1 0 . 1 0 0 ) . Y P L C i 1 0 ) , Y P L C l I 1 0 ) D I M E N S I O N U T < 2 , 1 0 0 ) R X = S R H O * U O * U D / S V I S C S H W A = S R H P * U G * * 2 . 0 * 0 . 1 8 5 / f A L 0 G 1 0 { R X ) ) * * 2 , 5 8 D E L T A L = 5 . 3 * U 0 * R X * * « - 0 . 5 ) D E L T A = 0 . 3 7 * U D * R X * * ( - 0 . 2 ) F A C L = D E L T A / D E L T A L D E L T A = O E L T A / F A C L V ( i ) = 0 . 0 Y ( 2 ) = 7 . 0 / < S H W A / S R H 0 ) * * 0 , 5 * S V I S C / S R H C Y < 2 ) = Y ( 2 I / F A C L F A C T 0 H = 0 . 6 5 X X ( 2 ) = U H X X < 3 1 = X X ( 2 ) + 0 . 3 7 *F A C TO R * X X { 2 ) *{ U O * X X { 2 I * S R H 0 / S V I S C ) * * ( - 0 . 2 ) X X ( 1 ) = U 0 - X X ( 3 ) + X X { 2 ) R X l = S R H n * U O * X X l D / S V I S C D E L T A 1= 0 . 3 7 * " X X < 1 ) * R X I * * ( - 0 . 2 ) D E L T A I = D E L T A 1 / F A C L S H W A 1 = S R H 0 * U C 1 * * 2 . 0 * 0 . 1 8 5 / 1 A L O G I C { R X l l 1 * * 2 . 5 8 D O 1 0 0 J = 3 T N 1 0 J Y l J ) = l . l * Y ( J - l ) Y P L M = ( S H W A / S P H 0 , * * 0 . 5 / S V t S C * S R H Q D O 1 0 2 J = 2 , N I F ( Y ( J ) . L E . D E L T A > U< 2 , J l = U O * { Y ( J ) / D E L T A | * * l 1 , 0 / 7 . 0 ) I F { Y< J ) . L E . D E L T A l ) IJ { 1 , J ) = U 0 * ( Y ( J ) / D E L T A 1 ) * * { 1 . 0 / 7 . 0 ) 1 0 _ . C O N T I N U E U T ( 1 , 1 I = 0 . 0 U T ( 2 , 1 ) = 0 . 0 D O 1 0 4 J = 2 . 3 0 X X ! rs! X X I X ~ X — c x c x N c * X ~ - ~ 1 I -3 C" I or 00 #> t ; 00 i-4 X - H I (X I — ) 0 0 w f> >• t_> >» I 00 I *—' > > w OO O OO > - > • #> ««* o —* NO t > • N > - • — "sn*. 4». t—4 r4-"t «—4 — <I | <S I - 4 S. -) S -3 ! X - X • - j O0 CM O0 •—t > rv i— or. + Q X O " 7 « Z —)<£—<[ — >-)>-) 0 \ i *nr-> I -3 — I r-4 ,i£ I  II • x; v: —4 • <• n — — XT ^ y « — Ovi ( V X. • » — C — — ~ O J II 1 1 f-* a: in a; • r— Z tn — • _ J 0\1 J r-4 t~ r-4 D <£ I— < > C J X O z o o o o t o II - J 4J> mm* > o O > > - 4 — — — O oj rv oj r-i in m m — »»r rrrr |— Z U J LU U J < or. h h hi J r — i * — i r—• 0C I — O cc cc cc O U J Z x r< *2 u a L U o a. < z t—« s: UJ z U - O U . »-4 o *-X "3 o o r-4 U _ 00 X UJ 


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items