UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Prediction of combined heat and mass transfer for a fully developed flow of supercritical water Xu, Li 2000

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

Item Metadata

Download

Media
831-ubc_2000-0619.pdf [ 3.39MB ]
Metadata
JSON: 831-1.0080975.json
JSON-LD: 831-1.0080975-ld.json
RDF/XML (Pretty): 831-1.0080975-rdf.xml
RDF/JSON: 831-1.0080975-rdf.json
Turtle: 831-1.0080975-turtle.txt
N-Triples: 831-1.0080975-rdf-ntriples.txt
Original Record: 831-1.0080975-source.json
Full Text
831-1.0080975-fulltext.txt
Citation
831-1.0080975.ris

Full Text

Prediction of Combined Heat and Mass Transfer for a Fully Developed Flow of Supercritical Water By L i Xu B.Sc, Tsinghua University, Beijing, 1984 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE In THE F A C U L T Y OF GRADUATE STUDIES DEPARTMENT OF MECHANICAL ENGINEERING We accept this thesis as conforming To the required standard THE UNIVERSITY OF BRITISH COLUMBIA May, 2000 © L i Xu, 2000 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. The University of British Columbia Vancouver, Canada D E - 6 (2/88) ABSTRACT Supercritical water oxidation (SCWO) is a technology for destroying organic wastes. In this process, water, organic waste and oxidant are brought together, and organics are oxidized to carbon dioxide and water. Under supercritical conditions, water has different features compared to normal water. Inorganic salts become insoluble at supercritical condition. Salts dissolved in organic wastes and produced from base and acid reactions in the system precipitate out of solution in the reactor. Fouling from salt is a serious problem for SCWO development. The salt deposition rate is proportional to the mass transfer coefficient and concentration difference between the wall and the bulk solution. In this work, a mass transfer model using concepts introduced by Deissler is developed based on an analogy among momentum, heat and mass transfer. The study focuses on supercritical water flow in a heated horizontal pipe. We ignore entrance length effects and assume that it is fully developed turbulent flow. Buoyancy is neglected in present research, so the flow is axisymmetric. The mixing length model with Van Driest's expression is used to calculate an eddy viscosity. We call this method to calculate heat and mass transfer coefficient the "Deissler-Van Driest" model. The model is validated by comparing laminar parabolic velocity profiles and turbulent velocity profiles for constant properties. The agreement in both cases is very good. Different mesh sizes are tested to check the effect on the results. Further validation is performed for variable properties: results are compared to the velocity profiles and heat transfer coefficients available in the literature, and also heat transfer coefficients predicted by Swenson et al's correlation (Swenson, H.S., Carver, J.R., and Kakarale, C.R., J. of Heat Transfer, November 1965.) Mass transfer coefficients calculated using the Deissler-Van Driest model are very similar to those predicted by the Swenson et al. correlation (with Prandtl number replaced by Schmidt number). Mass transfer coefficient increases steadily with bulk temperature. There is no peak in the pseudocritical region (where heat transfer coefficient peaks). Under supercritical conditions, the diffusion coefficient increases with temperature, but density and viscosity decrease with temperature. Higher diffusion coefficient, lower density and viscosity accelerate the diffusion of i i molecules and the transport of bulk fluid. This increases the mass transfer coefficient, especially at pseudocritical temperature. iii Table of Contents A B S T R A C T ii T A B L E OF C O N T E N T S iv L IST OF F I G U R E S vi A C K N O W L E D G E M E N T S vii i 1. INTRODUCTION A N D L I T E R A T U R E R E V I E W 1 1.1. Supercritical water oxidation 1 1.2. Heat transfer in supercritical water 2 1.2.1. Experiments 3 1.2.2. Computation 4 1.3. Mass transfer in supercritical water 7 1.4. Scope and outline of thesis 9 1.5. References 10 2. M O D E L F O R M U L A T I O N 14 2.1. Statement of problem < 14 2.2. The governing equations 14 2.3. Turbulence modeling 15 2.4. The Deissler-Van Driest model 19 2.5. The numerical procedure 23 2.6. References 25 3. M O D E L V A L I D A T I O N 27 3.1. Constant properties 27 3.2. Velocity and temperature profiles with heat tranfer 28 3.3. Heat transfer coefficient predictions 30 3.4. References 35 4. M A S S T R A N S F E R PREDICTIONS A N D DISCUSSION 60 4.1 Concentration profile 60 4.2 Mass transfer coefficients 60 5. C O N C L U S I O N S AND R E C O M M E N D A T I O N S 72 5.1. Conclusions 72 5.2. Recommendations 74 A P P E N D I X A : C O M P U T E R C O D E 75 iv A P P E N D I X B: S A M P L E INPUT A N D OUTPUT F ILES Chapter 1 Figure 1.1. List of Figures Water properties at different pressure 13 Chapter 2 Figure 2.1. Supercritical water flowing in a horizontal pipe with boundary conditions 26 Chapter 3 Figure 3.1. Laminar velocity profile for different mesh size (constant property) 36 Figure 3.2. Turbulent velocity profile (constant property) 37 Figure 3.3. Turbulent velocity profile for different mesh size (constant property, mesh=737, dy + =0.1, dy,+ = LOW)/,../) 38 Figure 3.4. Turbulent velocity profile compared with Deissler calculation at 34.47MPa, wall temperature 482°C, R=3.15mm 39 Figure 3.5. Effect of heat flux and mass flow rate on velocity profile 40 Figure 3.6. Effect of heat flux and mass flow rate on temperature profile 41 Figure 3.7. Velocity profile for different varation of shear stress 42 Figure 3.8. Temperature profile for different variation of heat flux 43 Figure 3.9. Influence of mesh size on velocity profile (mesh=549, dy* =0.1, dy,+ =1.0kfy M + ) 44 Figure 3.10. Influence of mesh size on temperature profile (mesh=549, dy,+ =0.1, dy,+ =1 .0k iy M + ) 4 5 Figure 3.11. Influence of mesh size on heat transfer coefficient (mesh=554, dy+ =0.1, dy,+ = 1 . 0 1 ^ M + ) 46 Figure 3.12. Heat transfer coefficient compared with Yamagata's measurement (horizontal pipe) 47 Figure 3.13. Heat transfer coefficient compared with Yamagata's measurement (horizontal pipe) 48 Figure 3.14. Heat transfer coefficient compared with Yamagata's measurement (vertically upward) 49 Figure 3.15a. Heat transfer coefficient compared with UBC experiment (bulk temperature).... 50 vi Figure 3.15b. Heat transfer coefficient compared with UBC experiment (bulk enthalpy) 51 Figure 3.16a. Heat transfer coefficient compared with UBC experiment (bulk temperature).... 52 Figure 3.16b. Heat transfer coefficient compared with UBC experiment (bulk enthalpy) 53 Figure 3.17a. Heat transfer coefficient compared with Bazargan model 54 Figure 3.17b. Heat transfer coefficient compared with UBC experiment (bulk temperature).... 55 Figure 3.17c. Heat transfer coefficient compared with UBC experiment (bulk enthalpy) 56 Figure 3.18. Heat transfer coefficient for different turbulent Prandtl number 57 Figure 3.19. Effect of heat flux and mass flow rate on heat transfer coefficient 58 Figure 3.20. Heat transfer coefficient for different variation of heat flux 59 Chapter 4 Figure 4.1. Effect of heat flux and mass flow rate on concentration profiles (molecule diameter = 5.52e-10m) 64 Figure 4.2. Effect of molecule diameter on concentration profiles 65 Figure 4.3. Mass transfer coefficient for molecule diameter =5.52e-10m 66 Figure 4.4. Diffusion coefficient of molecule using Stokes-Einstein relation 67 Figure 4.5. Mass transfer coefficient for different molecule diameter 68 Figure 4.6. Effect of heat flux and mass flow rate on mass transfer coefficient (molecule diameter = 5.52e-10m) 69 Figure 4.7. Influence of turbulent Schmidt number on mass transfer coefficient (molecule diameter = 5.52e-10m) 70 Figure 4.8. Mass transfer coefficient for different variation of mass flux (molecule diameter = 5.52e-10m) 71 vii A C K N O W L E D G E M E N T First, I'd like to thank my supervisor Dr. S. Rogak. His valuable guidance and assistance came through my study and research at UBC. I am indebted to him for helping me in the direction of my research and the process of writing my thesis. He spent a lot of time to discuss every detail of my simulation results and improve my thesis writing. The most important thing I learned beyond the Degree is the method and skill how we should analyze and deal with a difficult problem. I am grateful to my supervisor for providing this great opportunity to study at UBC. I also wish to thank the other two members of my thesis committee. Dr. Bushe provided me with guidance and assistance during my simulation and thesis writing. I am grateful to Dr. Green for his interest and advice in the present work. I wish to thank Dr. Ollivier-Gooch for his guidance and help in CFD course project. Dr. Fraser and Majid Bazargan must be acknowledged. Majid Bazargan provided me with water properties, and UBC experimental data for comparison. I got a lot of information from Majid's collection and discuss some problem during my work. I also thank Sanja Boskovic, and other individuals in our SCWO group. I thank Dr. Liqiu Wang for providing me useful information. I acknowledge many other friends who are not recorded here. Finally, I thank my mother who provided encouragement and support. I am indebted to her in all of my life. V l l l 1. INTRODUCTION AND L I T E R A T U I R E R E V I E W 1.1 Supercritical water oxidation The critical point of water is 374° C and 22.1 MPa. At pressures and temperatures higher than this critical point, water is supercritical. At supercritical pressures, water changes continuously from a liquid-like state to a vapor-like state with increasing temperature. There is no distinct interface between a liquid and vapor at supercritical pressures. The temperature at which — ap is a maximum is known as the pseudocritical temperature . The properties of water change dramatically near pseudocritical temperature (Figure 1.1). We can see from Figure 1.1 that the density of water drops to a tenth of the normal water density. The other properties such as heat capacity, thermal conductivity, and viscosity also change significantly. The specific heat peaks at or very near the pseudocritical temperature. Supercritical water oxidation (SCWO) is a process for destroying organic wastes. At supercritical conditions, water, organic waste and oxidant are brought together, and organics are rapidly oxidized to carbon dioxide and water. Because reaction occurs in a single-phase, the destruction efficiency is very high. Gloyna et al. (1994) reviewed the engineering aspects of SCWO technology, including organic destruction, inorganic solubility, corrosion and heat transfer. The performance of the first commercial SCWO industrial waste facility (McBrayer and Griffith, 1995) demonstrated the large-scale commercial viability of the process. These reviews suggest that SCWO technology is an environmentally attractive, safe and innovative wastewater treatment option. Inorganic salts become insoluble at supercritical conditions. Salts dissolved in the waste and produced from acid-base reactions in the system precipitate out of solution in a SCWO reactor. In a tubular reactor (eg., the UBC pilot plant), salts deposit on the tube wall and inhibit heat transfer. It is possible that salt deposition could lead to plugging of the pipe. An approach for controlling salt build-up is to flush the reactor with cold water This is very close to where the maximum in Cp occurs. 1 to dissolve the deposited salts. This would result in substantial downtime for a commercial SCWO plant. Fouling may occur in various parts of a SCWO facility including preheaters where rising bulk fluid temperatures decreases salt solubility and at the same time results in large wall -bulk temperature differences. This fouling problem motivates the present research. Salt deposition rates on the pipe wall depend on the mass transfer coefficient and mass transfer driving force (concentration difference between wall and bulk solution). The objective of this research is to model mass transfer in the UBC/NORAM supercritical water oxidation facility in order to predict how fast a reactor can be plugged. The approach is to use existing knowledge about heat transfer to develop mass transfer models. 1.2 Heat transfer in supercritical water Over the past few decades, as Hall et al. (1967-68) and Jackson and Hall (1979) reviewed, heat transfer to supercritical fluids has been studied intensively. Existing correlations based on constant properties were found to be inadequate to predict heat transfer in the critical region. The upstream conditions affect local heat transfer coefficients more than for constant property flows. The specific heat influences the peak of the heat transfer coefficient near the pseudocritical temperature. Jackson and Hall wrote the basic conservation equations-continuity, momentum and energy-in dimensionless form by using the inflow properties as reference. The entry conditions were specified as uniform temperature and some prescribed velocity distribution. They did not solve the governing equations themselves, however they did obtain expressions for certain dimensionless groups. These indicate that the Nusselt number must be a function of Prandtl number, Reynolds number, and other non-dimensional variables in the following form: d ph Ph Cph kh Where Nu is the Nusselt number 2 Re is the Reynolds number Pr is Prandtl number L is the length of the pipe [m] d is the diameter of the pipe [m] p is the density of fluid [ kg I nr' ] /u is the viscosity [kg/msj Cp is the specific heat capacity of the fluid [kJ/kgKJ k is the thermal conductivity of the fluid [kW/mKJ with the subscripts w and b referring to properties evaluated at wall temperature and bulk temperature. Note that for a given pressure, all fluid properties depend on temperature, so the last 4 dimensionless ratios above are coupled. From an experimental point of view, the heat transfer coefficient is often correlated with mass velocity G (kg/m2s). 1.2.1 Experiments Swenson et al. (1965) obtained heat transfer coefficients experimentally over a wide range of data for heat transfer to supercritical water flowing vertically upward. Reynolds number varied from 27,000 to 680,000. Nusselt number was from 60 to 3600. Prandtl number varied from 0.8 to 9.4. The inside diameter was 9mm, and heated length was 1.8288m. For fluid temperatures above (but near) the pseudocritical temperature, the heat transfer coefficient decreased along the pipe length due to the inlet effect, and then started to increase to the value for fully developed turbulent flow. The inlet effect was no longer significant after approximately L/d=97. Large property variations might extend the thermal entrance length. In the pseudocritical region, heat flux strongly affected heat transfer coefficient. At low heat fluxes, the heat transfer coefficient had a sharp maximum near the pseudocritical temperature. At high heat flux, the heat transfer coefficient curve was relatively flat and was much lower. This decrease in heat transfer coefficient with increase in heat flux was significant as the pseudocritical temperature was approached. The heat transfer coefficient was also found to be directly proportional to the mass velocity. The correlation developed from these experiments is discussed in section 3.3. Yamagata (1972) also conducted experimental investigations on forced convective heat transfer to supercritical water in horizontal and vertical pipes. The heat transfer coefficient reached a maximum at a bulk temperature slightly lower than the pseudocritical temperature, and its maximum value decreased with increasing heat flux. Compared to the Swenson correlation, Yamagata's data predicted the heat transfer coefficient to be considerably larger. At higher heat flux and lower flow rate, the measured wall temperature was remarkably higher than what was calculated by the correlation equation developed for low heat fluxes, i.e. heat transfer deteriorated. At some higher heat fluxes to horizontal pipes, the heat transfer coefficient was not uniform around the pipe periphery. It was higher at the bottom of the horizontal pipe and lower at the top. Adebiyi and Hall (1976) carried out an experimental investigation of heat transfer to carbon dioxide at supercritical pressure. Sufficiently large diameter pipes were used in their experiment to study the role of buoyancy. Peripheral temperature varied for a uniformly heated horizontal pipe flow. Their experimental results showed a marked deterioration in heat transfer at the upper part of the pipe and an improvement at the lower part. The buoyancy effect was still developing after 100 diameters of heated length. They suggested that abnormally long thermal entry lengths were a characteristic of supercritical pressure fluids. Harrson and Watson (1976) observed multiple peaks in wall temperature in their experiment at test section. They found that wall temperature varied in a regular manner. These variations induced local axial conduction heat fluxes, according to the authors. Other authors had reported multiple peaks in wall temperature at high heat fluxes. 1.2.2 Computation 4 Deissler (1954) analyzed heat transfer and fluid friction for a fully developed turbulent flow of supercritical water. He assumed that the eddy diffusivities for momentum and heat transfer were equal (Pr, = 1); the shear stress and the heat flux were constant across the pipe; and the static pressure could be considered constant across the pipe. He derived one-dimensional differential equations for shear stress and heat transfer which included molecular and turbulent transport. The eddy diffusivities were evaluated in two flow regions. For the region away from the wall, the Von Karman expression was used to calculate the eddy diffusivity. For the region close to the wall, the eddy diffusivity was taken to be a function of velocity u and distance from the wall y . This analysis indicated that the effect of variable shear stress and heat transfer across the pipe on the velocity and temperature distribution was negligible. It was also concluded that the effect of molecular shear stress and heat transfer in the region far from the pipe wall was slight. Goldmann (1954) proposed a method for analysis of heat transfer to supercritical water in fully developed turbulent flow. This method was based on an assumption which allowed an adaptation of the universal velocity profile to flow with variable properties. He used the analogy concept between momentum and heat transfer in fully developed turbulent flow. Goldmann suggested an integrated form of u+ and y+ for nonisothermal flow, which will be discussed later in section 3.3. Swenson et al. found that these models (particularly that of Deissler) compared poorly with their measurements, but suggested that discrepancies in the water properties might account for the differences. Malhotra (1977) investigated forced convective heat transfer to supercritical carbon dioxide flowing in a circular vertical duct. A two-dimensional numerical model was developed to predict heat transfer. The equations of continuity, momentum (x-direction) and energy were solved. The pressure gradient in the momentum equation was found for a specified mass flow rate. The kinetic energy-dissipation (k- s) model was then used to get a turbulent viscosity. Since the (k-s) turbulence model is applicable only in the fully 5 turbulent region, a mixing length model was still necessary in the laminar and transitional regions close to the wall. In Malhotra's study, the turbulent Prandtl number was assumed to be unity. To solve the equations of continuity, momentum and energy, boundary conditions are needed. At the axis of symmetry, the velocity and temperature gradients were set to zero. At the wall, the velocities were zero, and either the heat flux or the wall temperature were fixed. The inflow enthalpy profiles were assumed constant in his investigation, but the inflow velocity profiles were not normally constant. In turbulent flow problems, it is not necessary to have accurate inflow velocity profiles since downstream results are not very sensitive to the inflow profiles. In the numerical procedure, the inflow velocity profiles would normally follow from the I/7thpower law. It is possible to use a laminar profile in the laminar sublayer and part of the transitional region (up to y+=ll), a log law in the rest of the transitional region and most of the turbulent region with perhaps a \/l'h power law profile near the center of the pipe. Such a procedure would inevitably introduce discontinuities. In Malhotra's study, the mixing length model with the Van Driest expression was used to produce inflow velocity profiles. His numerical simulation used staggered grid spacing. A finite difference method was used. The code was validated against various standard uniform property test cases. The calculated radial velocity and temperature profiles in supercritical carbon dioxide flowing through vertical pipe were compared with the experimental data. A comparison was also made for wall temperature distribution between the simulation results and measured values. The relative difference was about 15%. Overall agreement was considered reasonable for the calculation. The prediction was correct for the development of 'unusual' velocity profiles in supercritical carbon dioxide. These 'unusual' velocity profiles were due to buoyancy. It was concluded that both the physical assumptions used as well as the numerical solution procedure were reasonable but in need of further development. The use of a fully two-dimensional numerical procedure was necessary in 6 order to account for the thermal entrance length and the near wall behavior in supercritical fluids. Lee and Howell (1998) performed a numerical simulation to study turbulent convective heat transfer to supercritical fluids in the entrance region of a vertical tube. The mixing length model included the effect of density fluctuations on turbulent diffusivity and it predicted differences of up to 10% in turbulent diffusivity over the standard mixing length model which does not include density fluctuations. The finite difference method was used to solve governing equations iteratively with the SIMPLE solution algorithm. For velocities and enthalpy, relaxation factors were 0.5-0.75 for stability of convergence. Properties were renewed with an underrelaxation factor 0.25-0.5. These low relaxation factors were used to prevent divergence. The grid was very fine near the wall and became larger toward the center of the pipe. It was found that the significant property variations near the pseudocritical temperature caused a slow approach to fully developed conditions, but they didn't report the entrance length under supercritical conditions. They compared with the Swenson and Dittus Boelter correlations and the Yamagata data, and their simulation results were higher than the correlations and measurements (about 10%) higher than the Yamagata's measurement). Due to the strong coupling of the heat transfer and fluid velocity through the pressure and temperature dependent properties, they concluded that no simple correlation could be expected to give useful predictions except over very limited ranges, and design of equipment operating at near critical conditions might require two-dimensional numerical analysis for adequate prediction. 1.3 Mass transfer in supercritical water Compared with heat transfer investigation in supercritical water, mass transfer has only been researched in recent years. Oh et al. (1997) did numerical analysis for the MODAR vertical-vessel-type reactor. They developed a CFD model (using Fluent, RNG k-e turbulence model) to calculate the detailed flow field, temperature distribution and salt-particle trajectories in the reactor flow domain. The wall model was not specified in the paper. They investigated the effect of particle sizes on salt deposition on the wall. Particles smaller than 20 jum reached the top head and inlet nozzle structure. Medium-7 size particles impacted the middle sections of the reactor. Particles greater than 500 jjrn were transported and dissolved into the cold brine pool at the bottom of the reactor. Chan et al. (1994) experimentally investigated salt formation and deposition in a supercritical water reactor. They also developed a simple deposition rate model and compared plugging time predictions to the measured plugging times. It was assumed that all salt beyond the solubility immediately deposited on the pipe walls at the location in the flow where solubility was exceeded and supersaturation did not occur. This model can not predict the correct starting deposition location in the axial direction. Due to temperature differences between the wall and the bulk solution, salt will deposit at the wall before the bulk temperature reaches the solubility temperature. Chan et al.'s solubility data was inconsistent with that from other studies. They did not determine the deposition mechanism in their experiments, the deposition mechanism and deposition rate are clearly related. Rogak and Teshima (1999) examined experimentally the deposition of sodium sulfate in a heated flow of supercritical water. It appeared that salt deposition depended on turbulent diffusion of salt molecules to a hot surface. Salt thickness was inferred from outside tube temperature measurement. They also developed a simple heat and mass transfer model. No nucleation of salt particles in the bulk solution was considered in the model. They predicted peak salt layer thickness to within a factor of two of measured values. They concluded that particle nucleation was not important in their experiments and better mass transfer modeling was needed, because mass transfer was predicted simply by extending empirical heat transfer correlations. Hodes (1998) investigated experimentally and numerically salt deposition rates from near-supercritical salt solution to a heated cylinder. He developed an understanding of salt deposition kinetics and nucleation mechanism at supercritical conditions for laminar flow. It was assumed that salt nucleation occurred only at the salt layer-solution interface, and there was no homogeneous nucleation of salt in the bulk solution. This assumption was validated by visual observations. He developed a model that predicted homogeneous 8 nucleation and supersaturation are unlikely in the bulk solution at the investigated conditions. The deposition mechanism was molecular diffusion of salt to the wall surface and subsequent heterogeneous nucleation there. Hodes concluded that Salt deposition rate is proportional to the mass transfer coefficient and the concentration differences between the wall and the bulk solution. A reduction in mass transfer coefficients may increase, decrease or not affect the deposition rate at the wall surface. At the conditions investigated in the deposition rate experiments, decreased mass transfer coefficients led to increase in the deposition rate on the wall surface due to an increased driving force for mass transfer. Protopopov et al. (1994) presented a numerical solution of the set of differential laminar equations for convective heat and mass transfer to supercritical water flow in a vertical heated tube. They proposed a correction factor to modify an abnormal reduction of the coefficients of diffusion of impurities in water in the region of pseudophase transition at supercritical pressure. Coefficients of diffusion of water-suspended particles O.Oljum'm diameter at supercritical pressure showed minima in the pseudocritical region. 1.4 Scope and outline of thesis Chapter 1 reviewed previous experimental and theoretical work related to heat and mass transfer in supercritical water. The mechanisms of mass transfer and salt deposition model were discussed in this section, but further work is needed to understand mass transfer in supercritical water. The purpose of this study is to simulate momentum, heat and mass transfer in supercritical water in the U B C pilot plant. The mixing length model is used to calculate eddy viscosity. Due to the similarity among these transport phenomena, some well -known correlations of momentum and heat transfer can be compared with the numerical results. The mass transfer coefficients should be plausible i f the heat transfer data are in good agreement with experimental and correlation results. 9 The continuity, momentum, energy, and concentration equations are written in Chapter 2. The model uses Deissler approach with Van Driest's mixing length expression. Buoyancy effects are neglected. We ignore the thermal entrance length effect, and assume that it is fully developed turbulent pipe flow. Based on the analogy among momentum, heat and mass transfer, the law of the wall is used to model heat and mass transfer. The solution procedure is explained in this chapter. The program code is validated through comparing with a laminar parabolic velocity profile and turbulent velocity distribution for pipe flow with constant properties in Chapter 3. Further validation is implemented for variable properties by comparing with velocity profiles, heat transfer correlations and experimental data. The effect of different turbulent Prandtl number on heat transfer coefficient is discussed. Based on the Deissler-Van Driest model, mass transfer coefficients are calculated and compared with correlation in Chapter 4. The influence of heat flux, mass flowrate, diffusion coefficient, and turbulent Schmidt number on mass transfer coefficient is discussed in this chapter. Conclusions from numerical simulation are discussed in Chapter 5, and finally, the recommendations are stated for future improvements. 1.5 References Adebiyi, G.A., and Hall, W.B., "Experimental Investigation of Heat Transfer to Supercritical Pressure Carbon Dioxide in a Horizontal Pipe", Int. J. Heat Mass Transfer. Vol . 19. Pp. 715-720, 1976. Chan, J.P.C., LaJeunesse, C.A., and Rice, S.F., "Experimental Techniques to Determine Salt Formation and Deposition in Supercritical Water Oxidation Reactors", HTD-Vol. 296, Fire, Combustion, and Hazardous Waste Processing, Book No. Goo913, 1994. Deissler, R.G., "Heat Transfer and Fluid Friction for Fully Developed Turbulent Flow of Air and Supercritical Water With Variable Fluid Properties", Trans, the A S M E , January, 1954. 10 Gloyna, E.F., Li , L., and McBrayer, R.N., "Engineering Aspects of Supercritical Water Oxidation", Wat. Sci. Tech. Vol.30, No. 9, pp. 1-10, 1994. Goldmann, K., "Heat Transfer to Supercritical Water and Other Fluids with Temperature-Dependent Properties", Nuclear Engineering, C.E.P. Symposium Series, Vol. 50, No . l l , 1954. Hall, W.B., Jackson, J.D., and Watson, A., "A Review of Forced Convection Heat Transfer to Fluids at Supercritical Pressures", Proc. Instn. Mech. Engrs, Vol 182, Pt 31 1967-68. Harrison, G.S., and Watson, A., "An Experimental Investigation of Forced Convection to Supercritical Pressure Water in Heated Small Bore Tubes", Proc. Instn. Mech. Engrs, Vol 190, 40/76, ImechE 1976. Hodes, M.S., "Measurements and Modeling of Deposition Rates from Near-Supercritical, Aqueous, Sodium Sulfate and Potassium Sulfate Solutions to a Heated Cylinder", Ph.D. Thesis, MIT, 1998. Jackson, J. D., and Hall, W.B., "Forced Convection Heat Transfer to Fluids at Supercritical Pressures", Turbulent Forced Convection in Channels and Bundles, Vol. 2 1979. Lee, S.H., and Howell, J.R., "Turbulent Developing Convective Heat Transfer in a Tube for Fluids near the Critical Point", Int. J. Heat Mass Transfer. Vol. 41, No. 10, pp.1205-1218, 1998 Malhotra, A., "An Analytical Investigation of Forced Convective Heat Transfer to Supercritical Carbon Dioxide Flowing in a Circular Duct", PhD Thesis, UBC, 1977. McBrayer, R., and Griffith, J., "Operation of the First Supercritical Water Oxidation Industrial Waste Facility", The International Water Conference, 56th Annual Meeting, 1995, Pittsburgh. Oh, C.LI., Kochan, R.J., and Beller, J.M., "Numerical Analysis and Data comparison of a Supercritical Water Oxidation Reactor", AIChE Journal, Vol. 43, No.6 June 1997. Protopopov, V.S., Kuraeva, I.V., Gorbacheva, E.V., Zav'yalova, T.E., Makarov, M.V., and Fokina, N.I., "Convective Diffusion of Iron-Corrosion Products in an Aqueous Coolant", Thermal Engineering, Vol. 41, No.3, 1994, pp. 170-173. Rogak, S.N., and Teshima, P., "Deposition of Sodium Sulfate in a Heated Flow of Supercritical Water", AIChE Journal, Vol.45, No.2, February 1999. Swenson, H.S., Carver, J.R., and Kakarala, C.R., "Heat Transfer to Supercritical Water in Smooth-Bore Tubes", J. Heat Transfer, November 1965. 11 Yamagata, K., Nishikawa, K., Hasegawa, S., Fujii, T., and Yoshida, S., "Forced Convective Heat Transfer to Supercritical Water Flowing in Tubes", Int. J. Heat Mass Transfer. Vol 15. pp.2575-2593. 1972. 12 Figure 1.1. Water properties at different pressure 13 2 MODEL FORMULATION 2.1 Statement of the problem The problem in this study is forced convective heat and mass transfer to turbulent supercritical water flowing in axisymmetric horizontal pipe (Figure 2.1). Because the salt solution is very dilute, we use pure water properties during modeling. Since water properties change significantly near critical region, we can not use constant property during simulation. At given temperature and pressure, we can use a library of water properties of lapws957.f (Wagner and Vxu(3, 1999) to find the water properties. The influence of buoyancy in the momentum balance for forced convection in supercritical water is not fully understood. For a horizontal tube, if a buoyancy force is included in the momentum equation, the axisymmetric flow assumption will be invalid. Forced convection is the dominant heat transfer mode. For higher mass flowrate and lower heat flux, the buoyancy is not important. We ignore the buoyancy force and simplify this already complicated problem in the current study. 2.2 The governing equations The equations of continuity, momentum, energy and species conservation governing the two-dimensional symmetric flow can be written as follows (Rohsenow et al. 1998). Continuity: d(pur) | d(pvr) = Q dx dr x-momentum: d(pu2r) d(puvr) _ ^dp + d dx dr dx dr .du <9vN rju(— + —) dr dx dx ^4 du 2 dv 2v^ 3 dx 3 dr 3r (2.2) r-momentum: d(puvr) d(pv2r) _ ^dp + d dx dr dr dx ,dv 3ws rju(— + —) dx dr + • dr A dv 2 du 2v x rju( ) 3 dr 3 dx 3r (2.3) Energy: 14 d(puCTr) d(pvCTr) d dx dr dr M d(CrTX Pr dr + • dx M d{CpT) Pr dx dp dp , + ur h vr h /urq) dx dr (2.4) Where the viscous dissipation function is: 0 = 2 dr r dx ,dv du.7 2 + (—+ —y - -dx dr 3 1 5 . . du (rv) + — r dr dx (2.5) Species mass conservation: d(puClr) d(pvCjr) d dx dr dr 2.3 Turbulence modeling The governing equations (2.1 to 2.6) are still precisely correct for turbulent flow by using instantaneous pressure, velocity, temperature and concentration in the equations, but there is no need for us to consider the details of turbulence. We are only concerned with its time-averaged effects. We take a time average of each governing equation. We assume that it is statistically stationary flow. The instantaneous pressure, velocity, temperature and concentration equal to the time average plus a fluctuation term. P = p + p U - u + u V - v + v T = T + t' C = c + c Where p,u,v,T,c are time averages. The fluctuating terms are functions of time, and by definition: p = 0, u = 0, v = 0, t' = 0, c = 0 . The process of time averaging causes statistical correlation's involving fluctuating velocities, temperature and concentration to appear in the above governing equations. We have no direct way of knowing the magnitudes of these terms. We depend on turbulence models to know their effect in terms of known or calculable quantities. We solve the mean-flow equations with 15 rpD d^ dr dx rpD dx (2.6) turbulence model equations, and compute the relevant correlations. Thus we can simulate the behavior of turbulent flow in many important respects. In a standard mixing length model, the turbulent fluxes can be expressed as: - ^ ^ , 7 (2.7) or k, d(CpT) - p C / ' v ' = - f - ^ . (2.8) Cp dr -pc^^pD,8^- (2.9) or Unlike the molecular viscosity, ju, is the eddy viscosity and it is not a property of the fluid. Its value will vary from point to point in the flow. This is a major assumption that the motion is somehow isotropic and the turbulent dissipation is directly related only to local mean velocity fields. In the flows of interested here, the eddy viscosity is generally many orders of magnitudes greater than the molecular viscosity, except in the viscous sublayer adjacent of a solid wall. Similar considerations apply to the turbulent conductivity kt and diffusivity Dt. The turbulent Prandtl and Schmidt number are defined as: C nu, Pr, = ^ - ^ (2.10) Sc,=^- (2.11) pD, The turbulent Prandtl number Pr, must be determined for calculation of turbulent convective heat transfer. For supercritical water, most of the numerical (theoretical) studies use the same value of turbulent Prandtl number as for the constant property case because there is little information about property variation effects on turbulent Prandtl number. The turbulent Prandtl number is assumed to be around 1.0. The turbulent Schmidt number Sc, is often approximated as unity, but the effect of varying these two numbers will be considered later. The governing equations for turbulent flow (time-average) can be written as following (all overbars are removed because all quantities are averaged): 16 d(pur) d(pvr) = dx dr (2.12) d(pu r) ^  d(puvr) dp d dx dr dx dr r(ju + jut) du dr + • dx r(ju + ju,) du dx + • dr dx + • dx 1 du 2 dv 2v^ 3 dx 3 dr 3r (2.13) d(puvr) d(pv2r) _ dp d dx dr dr dr dr dx r(ju + ju,) dv dx dx du rM(—) or + • dr A dv 2 du 2v. rju{ ) 3 dr 3 dx 3r (2.14) d(puCpTr) | d(pvCpTr) _ d dx dr dr dp dp + ur h vr h jury) dx dr Pr Pr, dr + • dx Pr Pr, dx (2.15) d(puCir) d(pvCjr) d dx + dr dr pbc, dr + • dx rp(D + ^ ) d C < pSct dx (2.16) The prediction of the eddy viscosity is a major problem in turbulent flow. The mixing length model is used in this study to calculate the turbulent viscosity, and then the above governing equations can be solved. This turbulence model approach is a practical way of treating turbulent flows. We require a turbulence model to be accurate, economical and simple. What constitutes the 'best' turbulence model will differ according to the problem under consideration. The Prandtl mixing length model is simplest and earliest turbulence model. It provides an adequate basis for many engineering applications (Launder and Spalding, 1972). Prandtl adopted Newton's viscosity law and applied it to turbulent flow. The turbulent viscosity is related to the Prandtl mixing length / and mean velocity gradient by 17 (2.17) du The mixing length / is not a quantity that can be measured directly. It is reasonable to assume that the mixing length may scale with the distance from the wall / = ky, (empirically k = 0.4), but this is not valid for regions very close to the wall-the viscous sublayer. The mixing length has to be modified to account for laminar effects. Van Driest (1956) introduced a damping function into the mixing length equation. Van Driest suggested . l = ky[l~exp{-y+ I A+)\ (2.18) When the argument of the exponential is close to zero, the mixing length / will be lowered. The constant A+ =26 is a non-dimensional empirically determined, effective sublayer thickness. Lee and Howell (1998) included the effect of density fluctuations on turbulent transport characteristics, expressed as: ju, = pl2 du dy R I d(CpT) R i d(CpT) "C Pr, dy Cp Pr, by (2.19) Where /? is compressibility, and y (2.20) Lee and Howell (1998) predicted differences of up to 10% in turbulent diffusivity over the standard mixing length model (equations 2.17 and 2.18). By applying the appropriate boundary conditions and using the turbulent viscosity equations (2.17) and (2.18), the governing equations can be discretized and solved in a computational domain. The velocity, temperature, and concentration are calculated in every control volume. Furthermore, the heat transfer coefficient, mass transfer rate and mass transfer coefficient can be computed along the pipe. The wall temperatures can then be compared with the measurements. The inflow temperature profiles are constant, and then heat is added. This is a thermally developing flow starting from the inlet. The hydrodynamics and thermal entrance lengths in 18 turbulent pipe flow are much shorter than their corresponding lengths in laminar pipe flow, so we neglect the entrance length effect and assume that it is fully developed turbulent flow (for 3u 3D 3D constant properties, —- = 0, v = 0, — = const, -^- = 0). We can largely simplify the above ox dx dr process and use a one-dimensional model. 2.4 The Deissler-Van Driest model The method of analysis for one-dimensional fully developed turbulent flow here is adapted from Deissler (1954). He assumed that the turbulent Prandtl number was 1.0, the variations of shear stress and heat flux across the tube had a negligible effect on the velocity and temperature distribution, and the static pressure was constant across the tube. From the x-momentum (2.13), energy (2.15), and concentration (2.16) equations, the shear stress, heat transfer and mass transfer equations are written in the following form r = (ju + /ul)~ (2.21) dy q = -(k + -L—)— (2.22) Pr, dy j = (D + ^-)^- (2.23) pSc, dy dC. Note that the term T—- , presented in equation 2.15, is neglected here, as in earlier works by dy Deissler (1954), Goldmann (1954) and others. We integrate these equations from the wall to the tube center. Due to a fine mesh size for every step, we use constant properties in each cell. Momentum, heat and mass transfer include contributions both molecular diffusion and turbulent mixing. In the laminar sublayer adjacent to the wall, the molecular viscosity is dominant. In fully turbulent region, the eddy viscosity is dominant and may be 10 times as large as the molecular values. 19 For equation 2.23, the Stokes-Einstein relation is used to calculate the diffusion coefficient D of a molecule: b T D = ^ ±~r (2-24) 37TjLldm Where bk = 1.38e — 23 [ J/K] is Boltzman's constant. dm is the diameter of a molecule. This diameter is difficult to estimate. Rogak and Teshima (1999) tested the diameter range from 2 to 60 Angstroms in their model. In this model, we use 2, 5, and 60 Angstroms in Chapter 4. The similarity of the mathematical form among the momentum, heat and mass transfer equations (2.21-2.23) defines this group of phenomena as analogous. There are no detailed experimental measurements on concentration fluctuations in a turbulent flow field, so we must depend on some assumed analogy among the mass, momentum, and heat transfer in order to relate mass transfer data to momentum or heat transfer data. We need to solve the equations (2.21) to (2.23) above in order to calculate the velocity, temperature and concentration profiles in fully turbulent pipe flow. It is important to have appropriate expressions for shear stress, heat flux and mass flux on the left-hand side of the above three equations. For constant properties, the laminar velocity profile is parabolic. The shear stress is proportional to the derivative of velocity, so the shear stress is a linear function of radius for laminar flow. For constant properties and fully developed statistically stationary turbulent flow, the pressure forces are balanced by the shear stress so that: nr2 dp = r27irdx (2.25) We get from equation (2.25) that = T- (2.26) 2dx r dp For fully developed turbulent flow, — is constant, therefore the shear stress varies linearly with dx radius as well. The mean turbulent shear stress in the center of the tube must be zero because the mean velocity gradient is zero. For turbulent flow under supercritical conditions, we assume a linear relationship exists for shear stress: 20 — = l - f (2.27) For similarity, we assume linear relationship exists for heat flux and mass flux as well: — = 1 - - (2.28) ci.., R ./, R (2.29) These assumptions will be discussed in more detail later, but it should be noted that Deissler and Eian (1952) demonstrated that the velocity, temperature, and concentration profiles are relatively insensitive to shear stress distribution for fully developed turbulent flow in a pipe. Hatton (1964) and Webb (1971) reached the same conclusions for a turbulent boundary layer and pipe flow. This is consistent with the fact that in calculating heat transfer and mass transfer rates, the primary resistance is concentrated near the wall, and this region near the wall is of the most essentia] consideration. Equations (2.21) to (2.23) are non-dimensionalized by dividing both sides by the shear stress, heat flux and mass flux at the wall. T (u+ a,) du *"„ rw dy qw P r , qw dy • i - = ( 0 + ^ - ) - L * ( 2 . 3 J , Jw pSc, j„ dy The eddy viscosity is computed from the Prandtl mixing length model together with the Van Driest expression. The characteristic properties are evaluated at wall conditions: 21 u (2.33) Goldmann (1954) proposed an integrated form to calculate u+ and y+: du = I- dy •o (2.34) Lee and Howell (1998) used both equations (2.33) and (2.34) in their calculation of the mixing length. There was less than 3% difference in the heat transfer coefficient between these two equations. There were very small differences between their predicted velocities and temperatures. The non-dimensional mixing length (from equation 2.18) is expressed as: / + =^ + [ l -exp(-^- ) ] Using the above characteristic properties and substituting the shear stress, heat flux and mass flux linear relationship (2.27) to (2.29) into equations (2.30) to (2.32), we get the non-dimensional velocity, temperature and concentration gradients as following: 22 du+ _-M + ^ M +4pl+2(l-y/R) dy+ ~~ 2pl+2 dT+ l-y/R k/Pv+pcl / P r w i n i + I dy dC+ _ \-ylR (2.35) d y + D/Sc+r2d^/Sc, dy+ where P = —> P = —, k = lT> c P = > D = r V Pw Pw w C pw Uv, r p + Pw^pW yj^W ' Pw ^rp rjr^ aw £~t+ *\j^"w ^ Pw £ The velocity and concentration at the wall are zero. The heat flux, mass flux and wall temperature are given. The preceding method is adapted from Deissler, and the Prandtl mixing length with the Van Driest expression is used to calculate the eddy viscosity. We call this method the "Deissler-Van Driest" model. Majid Bazargan (Bazargan, 1999) is currently working with professor Fraser in UBC on a fully developed heat transfer model which is based on the Goldmann's integrated definition of u+ and y + . This model is referred to as "Bazargan model" in the figures of the present thesis 2.5 The numerical procedure Appendix A provides the FORTRAN code for this model, with sample input and output files in Appeddix B. It is required that the numerical procedure for this problem should be applicable to variable property flows and should be efficient in terms of computing time. The program code is 23 Scwater.f and variable water properties under certain temperature and pressure are calculated from lapws957.f. The input data for this program are turbulent Prandtl and Schmidt numbers, tube radius R (mm), pressure P (MPa), mass velocities G (kglslm2), heat flux Qw(kW I m2), wall temperature TW(°C), and maximum wall temperature. The velocity at the wall is zero. The wall temperature is given. The concentration of solution at the wall surface is assumed to be zero. Because we don't know the velocity profile initially, we must make an initial guess for the wall shear stress. Since there are large velocity, temperature and concentration gradients near the wall, we use a very fine mesh at the viscous sublayer. For the first step, dy+ =0.1 is used; this is then increased by a constant factor each step to the tube center. The effect of mesh size on the results will be discussed and tested later. Starting from the wall, for a step dy*, we can get the velocity, temperature and concentration by integrating equation (2.35). At this temperature and pressure (input), using a library of water properties, we can find the water viscosity, density, conductivity, heat capacity, enthalpy and diffusivity. Then we can calculate the mass flow rate, product of mass flow rate and enthalpy, product of mass flow rate and concentration in this annular element of area. This process is continued step by step to the tube center to predict velocity, temperature and concentration profiles. The mass flow rate, bulk enthalpy and bulk concentration are calculated from following equations: R m= ^2n(R-y)u(y)p(T(y))dy o R mH{Th)= \2rc{R-y)u(y)p{T{y))H(T(y))dy o R mC{Th)=\2n{R- y)u(y)p(T(y))C(y)dy o From the bulk enthalpy and heat capacity, the bulk temperature is calculated. This integrated mass flow rate is compared with the mass flow rate from the mass velocity times the tube area, and the wall shear stress is adjusted. This process is repeated until the relative difference between these two mass flow rates less than 0.01. This means that the reported mass 24 velocities here could be in error by 1%, resulting in up to about 1% error in the heat transfer coefficients. From the heat transfer rate and mass transfer rate equations, we define: /,.= _ £ _ T -T w (2.36) where h is the heat transfer coefficient, similarly hm is the mass transfer coefficient. For a given heat transfer rate q and mass transfer rate/, we can calculate the heat transfer coefficients h and mass transfer coefficients hm from the above two equations. 2.6 References Bazargan, M , and Fraser, D., "Heat Transfer to Supercritical Water", Proceeding at CANCAM 99, 17th Canadian Con gress of Applied Mechanics, Held at McMaster University, Hamilton, Ontario, Canada, May 30-June 3, 1999. Deissler, R.G., "Heat Transfer and Fluid Friction for Fully Developed Turbulent Flow of Air and Supercritical Water With Variable Fluid Properties", Transactions of the ASME, January, 1954. Goldmann, K., "Heat Transfer to Supercritical water and Other Fluids with Temperature Dependent Properties", Nuclear Engineering, C.E.P. Symposium Series, Vol. 50, No.l 1, 1954. Launder, B.E. and Spalding, D.B., "Mathematical Models of Turbulence", 1972, Academic Press, London and New York. Lee, S.H., and Howell, J.R., "Turbulent Developing Convective Heat Transfer in a Tube for Fluids near the Critical Point", Int. J. Heat Mass Transfer. Vol. 41, No. 10, pp. 1205-1218, 1998. Rogak, S.N., and Teshima, P., "Deposition of Sodium Sulfate in a Heated Flow of Supercritical Water", AIChE Journal. Vol.45, No.2, February 1999. Rohsenow, W.M., Hartnett, J.P., and Cho, Y.I., "Handbook of Heat Transfer", Third Edition, 1998, McGraw-Hill. Wagner, W., and Prw/J, A., "The IAPWS Formulation 1995 for the Thermodynamic Properties of Ordinary Water Substance for General and Scientific Use", To be submitted to J. Phys. Chem. Ref. Data, 1999. 25 u(y) r x y u (0 )=0 C (0 )=0 T (0 ) = Tw r = R Figure 2 . 1 . Supercritical water flowing in a horizontal pipe with boundary conditions 26 f 3. M O D E L V A L I D A T I O N 3.1 Constant properties Validation will tell us whether the solutions we get for a series of simple test cases are correct. The test case we chose should have a known solution. Constant properties are used first. For fully developed laminar pipe flow, the velocity profile u(y) is known exactly: u = 2U[\-{^)2} (3.1) K where U is average flow velocity. From Chapter 2 equation (2.27), the shear stress is: ^ = l - i (3.2) R The friction factor / and wall shear stress r for laminar flow are (White, 1986): / = £ (33) Re rw = l~fpU2 (3.4) To compare these solutions to the numerical model, the eddy viscosity is set to zero. A Reynolds number of 753.7 was used arbitrarily; a constant step size (uniform mesh) was used. The numerical model is lower than the parabolic velocity near centerline at mesh =39 in Figure 3.1. The model velocity is increased when the mesh increases to 78. Our results and parabolic velocity agree very well at mesh=389. A further increase mesh to 1942 points has no significant effect. Turbulent velocity distributions are shown in Figure 3.2. The Van Driest universal velocity profile is compared with our model. The Van Driest expression for constant properties is as follows: 27 dy+ l + {l + 4/<V 2[l-exp(-v + / A+)]2}U1 where k=0.4, and A+ = 26. For both this "universal velocity profile" and the Deissler-Van Driest model, variable step size was used as discussed in Chapter 2. Both were calculated with 78 mesh points. Although this is too few, our model's prediction of velocity agrees very well with universal velocity profile because the same number of points was used for each. For shear stress varying linearly along the radius, we can see that in Figure 3.2 near the center of pipe, our model is a little bit lower than the universal velocity profile. For constant shear stress assumption, the agreement is very good, since the universal velocity is derived based on the assumption of constant shear stress. We can see from this figure, there is little difference between velocity profiles with different shear stress distributions. This means that for fully developed turbulent flow, the major resistance for momentum transfer is near the wall, as stated before. We compare the velocity profile with meshes 78, 508, 737 and 1652 points in Figure 3.3. Velocity increases as the number of nodes is increased from 78 to 737, but not from 737 to 1652. From the above laminar and turbulent velocity profile, we can see that the mesh size does affect the accuracy of results. During simulation, we always need compromise between computing time and accuracy. Since this is one-dimensional model, and it is not time consuming, we can choose a very fine mesh to calculate accurate results. The choice of mesh size is discussed further in section 3.3. 3.2 Velocity and temperature profiles with heat transfer There is no experimental data for velocity and temperature profiles across a flow of supercritical water flowing in a horizontal pipe. Therefore validation of our model is difficult. Deissler (1954) analyzed heat transfer to air and supercritical water in fully developed turbulent flow, and made measurements for air. Deissler used the Von Karman expression for eddy diffusivity in a different way for the region close to the wall and the region further from the wall. He calculated velocity and temperature distributions for both air and supercritical water. His experimental and calculated velocity and temperature distributions were in good agreement for air. 28 A comparison between our model and Deissler's predictions is shown in Figure 3.4 for a pressure P=34.47 MPa, and a wall temperature of 482°C. The integration step size is increased from dy+ =0.1 by a factor 1.01 for each step (about 500 points). The heat transfer parameter is given in the following equation: 1= 9 w, (3-6) At j5 = 0, the agreement is good between Deissler's calculation and our model results. For other beta values, our results are higher than Deissler's calculation. The reason is that we use a different turbulence model, and different properties. Deissler stated that the difference between his result and several other analyses of heat transfer to supercritical water is due principally to the difference in the physical properties used. We chose the values of P, q, G, R etc. for comparison with available measurements. In principle, this calculation method is applicable for any value of these parameters for turbulent flow. The following velocity and temperature profiles are for the UBC system. Figure 3.5 shows the velocity profile at different heat flux and mass flow rate for condition P = 24.2 MPa, R=3.15mm and Tw = 386°C. It is obvious that the higher the mass flow rate, the larger the average velocity, and the bigger velocity difference from the bulk solution to the wall. At lower heat fluxes, the temperature difference between wall and bulk solution becomes smaller, and the density difference becomes less. Smaller density difference means that the bulk solution density is less at lower heat flux. The average velocity is higher at smaller density for fixed mass flow rate, so the velocity difference is larger at lower heat flux. Figure 3.6 shows the effect of heat flux and mass flow rate on the temperature distribution at P = 24.2MPa, R = 3.\5mm, Tw =386°C. We can see from this graph that the temperature decreases significantly at higher heat flux. At zero heat flux, the temperature is uniform across the pipe radius. At low heat flux, the temperature difference between the wall and bulk solution is smaller, so the temperature profile is 2 9 flatter than that at high heat flux. Temperature decreases more at lower mass flow rate. The heat transfer coefficient increases with velocity or mass flow rate for forced convective heat transfer, so the temperature difference between the bulk solution and the wall will decrease for fixed heat flux and wall temperature. The temperature profile becomes flatter at higher mass flow rate. We compare the effect of shear stress distributions on the velocity profile in Figure 3.7. At the centerline, the velocity from linear variation of shear stress is 5% lower than that from uniform shear stress. The temperature profile can be seen in Figure 3.8 for different heat flux distribution. For fixed wall temperature, the center temperature from a linear variation of heat flux is about 1% higher than that from uniform heat flux. The higher bulk temperature will result in a slightly large heat transfer coefficient for linear variation of heat flux. From Figure 3.7 and 3.8, we find that the effect of variation of shear stress and heat flux across the pipe on the velocity and temperature profiles is slight for fully developed turbulent flow with variable fluid properties. The effects of mesh size on the velocity and temperature profiles are shown in Figure 3.9 and 3.10 at R = 3A5mm,P = 2A2MPa,G = S16kglm2s,q = \07kW/m2, 7w = 386°C. Velocity increases with mesh number increase in Figure 3.9. The bulk temperature decreases with mesh number increase in Figure 3.10. The velocity and temperature profiles do not change significantly when mesh points increase from 549 to 1277. During the simulation of heat and mass transfer coefficient, we use the fine mesh size ( about 550 points) in order to calculate more accurate results. 3.3 Heat transfer coefficient predictions Heat transfer coefficients can be calculated using the Deissler-Van Driest model and compared with empirical correlations and existing experimental data. Although there are dozens of data sets and correlations, we chose to compare the model with the Swenson et 30 al's correlation and Yamagata's data because these works are discussed fully by many authors. Normally, the Nusselt number is correlated with the Reynolds and Prandtl numbers, but the normal heat transfer correlation must be modified for variable properties at supercritical conditions. For example, the Swenson correlation (Swenson, et al. 1965) for supercritical water is: hd 'Gd~ 0.923 r-Nu = — = 0.00459 K Pw •T„ K 0.613 Pw Ph 0.231 (3.7) For constant properties, one may use the much simpler Dittus-Boelter equation (Lee and Howell, 1998) : Nu = — = 0 . 0 2 3 ( ^ ) 0 ' 8 ( ^ ^ ) 0 4 (3.8) K Ph kh where Nu is the Nusselt number [-] h is the heat transfer coefficient [ kW I m2 K'\ d is the diameter of the tube [m] k is the thermal conductivity of the fluid [ kW/mK ] G is the mass velocity [ kg I m2s] /LI is the viscosity [ kg/ms] Cp is the specific heat capacity of the fluid [ kJ/kgK ] / is the enthalpy of the fluid [ kJ/kg ] T is the temperature [ K ] p is the density of fluid [ kg I nr' ] with the subscripts w and b referring to properties evaluated at wall temperature and bulk temperature. For constant property flow, the heat transfer coefficient is constant for fixed Reynolds number. Near critical conditions, heat transfer will be enhanced. This is mainly because heat capacity is very high in this region. 31 We have shown that the mesh size affects the velocity and temperature profiles. We can see this influence on heat transfer. As shown in Figure 3.11, heat transfer coefficients decrease with mesh number increase. When mesh number increases from 60 to 554, the heat transfer coefficient decreases about 20%. Increasing the mesh size from 554 to 1281 changes the heat transfer coefficient by only 2%. The computing time is acceptable for this one-dimensional modeling, so we use mesh number of about 500 during simulations, calculated by dy* = 0.1; dy*=\.0]dyi_*. In Figure 3.12, we compare the Deissler-Van Driest model with Yamagata's (Yamagata, et al. 1972) measurements, and the Swenson correlation. The agreement is good at P = 24.5MPa,G = \260kg/m2s,q = 233kW/m2,R = 3.75mm. The heat transfer coefficient obtained with the Swenson et al. correlation is higher than that obtained with the Deissler-Van Driest model and measurement in the critical region, the Swenson correlation reaches its peak value before the pseudocritical temperature. In the subcritical region, our model and Dittus-Boelter's both overpredict the heat transfer coefficient, whilst the Swenson et al. correlation fits the experimental data better. In Figure 3.13, the mass velocity is increased to G = 1830% /m2s . Both the Deissler-Van Driest model and the Swenson correlation overestimate the heat transfer coefficient in the critical region. The agreement among the Deisser-Van Driest model, the Swenson correlation and measurements is good above the pseudocritical temperature. Our results are close to Dittus-Boelter's, both overestimate heat transfer coefficient in the subcritical region. Because our model does not include the buoyancy force in the momentum equation, our results can also be compared with Yamagata's measurement in a vertical tube as well. As shown in Figure 3.14 with R = 5mm, P = 24.5MPa, G = \200kg/m2s, q-465kW/m2, the Deissler-Van Driest model underestimates the heat transfer coefficient in the pseudocritical region and after this region, but the agreement is good in the subcritical region. We can see from these three graphs that the Dittus-Boelter correlation performs poorly in the critical region. It is used mainly for constant properties. Because water properties do 32 not change significantly in the subcritical region, we can compare our model results with the Dittus-Boelter correlation in the subcritical region. Our model results are also compared with UBC experimental data. Of the many UBC experiments, we selected only those believed to be no influence by buoyancy. At P = 24.2 MPa, G = 516kgIm2s, q = 107kW/m2, R = 3.15mm in Figure 3.15a and 3.15b. Near the pseudocritical temperature, the measurements fluctuate largely. The Deissler-Van Driest model appears to fit the measurement very well for this condition. In this experimental facility, there are two test sections in series, each with Z / D = 500 (ie. Entrance length effects should be minimal). At P = 25.2MPa, G = 955kg/m2s, q = 307kW/m2, R = 3A5mm, the heat transfer coefficients are plotted against bulk temperature and bulk enthalpy, shown in Figure 3.16a and 3.16b. The Deissler-Van Driest model's prediction is higher than the experiment. The Swenson et al. correlation fits the data better in the subcritical region. We have also compared predictions with a "Deissler- Goldmann" model developed by Majid Bazargan working with Dr. Fraser at UBC. Figures 3.17a,b, c, show heat transfer coefficients at P = 24.34MPa, G = 640kg/m2s, q = 300kW/m2, R = 3.15mm . Our model predicts lower heat transfer coefficients than Bazargan's model in Figure 3.17a, especially in the critical region. Our peak heat transfer coefficient is about 20% less than Bazargan's predictions. Comparison with experiment is given in Figure 3.17b and 3.17c. The agreement is good between the Deissler-Van Driest model and the measurement. Our model results are a little bit lower than the experiment data in subcritical region. The Swenson correlation performs poorly in this situation. Bazargan used different algebraic expressions to calculate the eddy viscosity. He also used a simple Prandtl mixing length model. This explains the discrepancy between the results. It should be noted that Bazargan's preliminary modeling shown in Figure3.17 uses a much coarser mesh than we have used. The influence of the turbulent Prandtl number is shown in Figure 3.18. The heat transfer coefficients decrease with turbulent Prandtl number especially the peak value. In numerical simulation, we normally chose the turbulent Prandtl number to be 0.9 or 1.0. In his review paper, Kays (1994) stated that direct numerical simulations indicate that turbulent Prandtl number is near 1.00 as the wall is approached. According to Kays, the experimental data confirmed this. For turbulent flow, the near wall region is most important for momentum, heat and mass transfer. We assume turbulent Prandtl number to be 1.0 in this study. The heat transfer coefficient is not a unique function of bulk temperature because the change in physical properties with temperature is not only large and nonlinear but also nonmonotonic. Heat flux and mass flow rate will affect the heat transfer coefficient as shown in Figure 3.19 at P = 24.2MPa, R = 3A5mm. The heat flux is increased from q = \07kW/m2 to 200kW/m2, and the maximum heat transfer coefficient near the pseudocritical temperature reduces dramatically for higher heat flux. But in the subcritical region, heat transfer coefficients are the same at different heat fluxes, and after the pseudocritical region, the heat flux affects the heat transfer coefficient slightly. We can see from this graph that the heat flux only has a significnantly influence on the heat transfer coefficient in the pseudocritical region and on the peak value. The mass velocity is increased from G = 576kg fm2s to lOOOkg/m2^. The higher the mass flowrate, the larger the heat transfer coefficient. The peak heat transfer coefficient increases significantly with mass flowrate. The variation of heat flux across the pipe has slight influence on the heat transfer coefficient as shown in Figure 3.20 at P = 24.34MPa,G = 640kg/m2s,q = 300kW/m2, R = 3A5mm . The heat transfer coefficient from linear variation of heat flux is higher by about 6.7% than that from the uniform heat flux in the subcritical region, and the peak value of the linear variation of heat flux is about 5.7% higher than the maximum heat transfer coefficient by using uniform heat flux. We use a linear variation of shear stress, heat flux and mass flux across the pipe during the simulation. 34 The heat transfer coefficients are compared with the correlations and experiments in the above figures. The Deissler-Van Driest model is as good as the Swenson correlation without adjusting any parameters. There are no measurements of mass transfer coefficients, but based on the analogy among the momentum, heat and mass transfer, if the calculated heat transfer coefficient is in good agreement with experimental results, we assume that the calculated mass transfer coefficients are acceptable. 3.4 References Bazargan, M. , and Fraser, D., "Heat Transfer to Supercritical Water", Proceeding at CANCAM 99, 17th Canadian Cong ress of Applied Mechanics, Held at McMaster University, Hamilton, Ontario, Canada, May 30-June 3, 1999. Deissler, R.G., "Heat transfer and Fluid Friction for Fully Developed Turbulent Flow of Air and Supercritical Water With Variable Fluid Properties", Transactions of the ASME, January, 1954. Kays, W.M., "Turbulent Prandtl Number-Where Are We?", Journal of Heat Transfer, May 1994, Vol. 116, pp284-295. Lee, S.H., and Howell, J.R., "Turbulent Developing Convective Heat Transfer in a Tube for Fluids near the Critical Point", Int. J. Heat Mass Transfer. Vol. 41.No. 10. pp. 1205-1218,1998. Swenson, H.S., Carver, J.R., and Kagarala, C.R., "Heat Transfer to Supercritical Water in Smooth-Bore Tubes", Journal of Heat Transfer, November 1965. White, F.M., "Fluid Mechanics", pp298-306, Second Edition, 1986. Yamagata, K., Nishikawa, K., Hasegawa, S., Fujii, T., and Yoshida, S., "Forced Convective Heat Transfer to Supercritical Water Flowing in Tubes", Int. J. Heat Mass Transfer. Vol. 15. Pp. 2575-2593, 1972. 35 0.25 R=3.15mm, Re=753.7 (constant property) i i i i i i 0.2 0.15 E g o parabolic velocity profile Z5 .. Deissler-Van Driest, mesh=39 0.1 p -. Deissler-Van Driest, mesh=78 0 — Deissler-Van Driest, mesh=389 r - Deissler-Van Driest, mesh=1942 0.05 n / i i I I i i QY. \ I i : 1 1 1 ' 0 0.5 1 1.5 2 2.5 3 3.5 y m (from wall to center) xio~ 3 Figure 3.1. Laminar velocity profile for different mesh size (constant property) 36 R=3.15mm ) Re=753739 (constant property) y+=y/y* Figure 3.2. Turbulent velocity profile (constant property) 37 Figure 3.3. Turbulent velocity profile for different mesh size (constant property, mesh=737, dyx+ =0 .1 , dy^ =1 .0h / y / _ 1 + ) 38 R=3.15mm, P=34.47MPa, Tw=482C Figure 3.4. Turbulent velocity profile compared with Deissler calculation at 34.47MPa, wall temperature 482°C,R=3.15mm 39 R=3.15mm, P=24.2MPa, Tw=386C — i — i — I I I I r— —i r—i—i—i—r— — q=107kW/m2, G=1200kg/m2s - q=107kW/m2, G=576kg/m2s - . q=400kW/m2, G=576kg/m2s 10* y+ Figure 3.5. Effect of heat flux and mass flow rate on velocity profile 40 R=3.15mm, P=24.2MPa, Tw=386C 390 i 1—<—1—• i i > • i 1—i—i—• * * *11 1 1—i—' i i i ' i 1 1—i—i i i 385 380 375 ^ 3 7 0 o 365 360 355 350 345 \ \ — q=107kW/m2, G=1200kg/m2s Nx -q=107kW/m2, G=576kg/m2s \ - . q=400kW/m2, G=576kg/m2s \ \ \ \ \ \ S N N I ' l I I 1 l l l l l 1 1 1 1 1—L I 1 I 10~1 10° 10 1 10 2 10 3 10 4 y+ Figure 3.6. Effect of heat flux and mass flow rate on temperature profile 41 R=3.15mm, P=24.2MPa, G=576kg/m2s, q=107kW/m2, Tw=386C y+ Figure 3.7. Velocity profile for different variation of shear stress R=3.15mm, P=24.2MPa, G=576kg/m2s, q=107kW/m2, Tw=386C f 1 | 1 1 1 1 1 1 1 1 | 1 i i—r i i i r | i i r—i—i—i t i | r r i i—r r T I - uniform heat flux -- . linear variation of heat flux V -V V v. — • ' 1 1 ' ,1,1 , I , , , , • l l u • • r i . • • 1 . . • • i • i • 1 • • i t i I 1 I 10" 10u 10 10 y+ 10J 10' Figure 3.8. Temperature profile for different variation of heat flux 43 R = 3 . 1 5 m m , P = 2 4 . 2 M P a , G = 5 7 6 k g / m 2 s , q = 1 0 7 k W / m 2 , T w = 3 8 6 C Figure 3.9. Influence of mesh size on velocity profile (mesh=549, dy* =0 .1 , d y j + = 1 . 0 1 ^ M + ) 44 Figure 3.10. Influence of mesh size on temperature profile (mesh=549, d y * =0.1, dy* =1.0kfr,._,+) 45 355 360 365 370 375 380 385 390 bulk temperature C Figure 3.11. Influence of mesh size on heat transfer coefficient (mesh=554, d y x + = 0.1, dy* = 1.0kfyM +) 46 Figure 3.12. Heat transfer coefficient compared with Yamagata's measurement (horizontal pipe) 47 1101 R=3.75mm, P=24.5MPa, G=1830kg/m2s, q=233kW/m2 TOO 0 90 OJ 1 80, • £ 701 CD - Deissler-Van Driest model - . Swenson + Dittus-Boelter o Yamagata experiment 280 300 320 340 360 380 400 420 bulk temperature C 440 460 480 Figure 3.13. Heat transfer coefficient compared with Yamagata's measurement (horizontal pipe) 48 R=5mm, P=24.5MPa, G=1200kg/m2s, q=465kW/m2 50 i—; 1 1 —i 1 1 :—r© 1 O 45 h 40 o CM E § 35 -»—> § 301 251 O O aj 201 U) c - Deissler-Van Driest model - . Swenson + Dittus-Boelter o Yamagata experiment 15 "+- , | + 4. 4- ^-zfc-^K o o CO CD 10 Vertically upward 260 280 300 320 340 360 380 400 420 bulk temperature C Figure 3.14. Heat transfer coefficient compared with Yamagata's measurement (vertically upward) 49 Figure 3.15a. Heat transfer coefficient compared with UBC experiment (bulk temperature) 50 70, R=3.15mm5 P=24.2MPa, G=576kg/m2s, q=107kW/m2 60 o CNJ E 50 CD O O , v_ 301 c 03 20 03 CD sz 10 - Deiss ler-Van Driest model - . Swenson o U B C experiment 1500 2000 2500 bulk enthalpy kJ/kg 3000 Figure 3.15b. Heat transfer coefficient compared with UBC experiment (bulk enthalpy) 51 451 R=3.15mm, P=25.2MPa, G=955kg/m2s, q=307kW/m2 1 1 1 1 + 40 h o <^35| 301 •4—1 c CD ; ^ 251 5= CD O O 2 0 | cc CD - Deissler-Van Driest model - . Swenson + Dittus-Boelter o U B C experiment 10 c^Po o ° o — O " ' 150 200 250 300 bulk temperature C 350 400 Figure 3.16a. Heat transfer coefficient compared with UBC experiment (bulk temperature) 52 R=3.15mm, P=25.2MPa, G=955kg/m2s, q=307kW/m2 451 1 1 1 1 1 1 : 1 :r_ 1 1— + 51 i L_ 1 1 1 1 1 1 1 1 600 800 1000 1200 1400 1600 1800 2000 2200 2400 2600 bulk enthalpy kJ/kg Figure 3.16b. Heat transfer coefficient compared with UBC experiment (bulk enthalpy) 53 301 R=3.15mm, P=24.34MPa, G=640kg/m2s ) q=300kW/m2 O CNJ E 25 h 20 h •+—> C CD O 0 15| O o CD H — CO co CD Deissler-Van Driest model - . Swenson — Bazargan model(1999) 200 250 300 350 400 450 500 bulk temperature C Figure 3.17a. Heat transfer coefficient compared with Bazargan model 54 Figure 3.17b. Heat transfer coefficient compared with UBC experiment (bulk temperature) 55 R=3.15mm, P=24.34MPa, G=640kg/m2s, q=300kW/m2 500 1000 1500 2000 2500 3000 3500 bulk enthalpy kJ/kg Figure 3.17c. Heat transfer coefficient compared with UBC experiment (bulk enthalpy) 56 R=3.15mm, P=24.34MPa, G=640kg/m2s, q=300kW/m2 1 1—: 1 1 : 1 1 i : i n l i | ' ' 1 1 —= 1 1 — 320 330 340 350 360 370 380 390 400 bulk temperature C Figure 3.18. Heat transfer coefficient for different turbulent Prandtl number R=3.15mm, P=24.2MPa o E 50 c CD ' o 40 3= CD O O 30 "Vi c CO < 20 co CD JZ — q=107kW/m2, G=1000kg/m2s - q=107kW/m2, G=576kg/m2s -. q=200kW/m2, G=576kg/m2s 320 340 360 380 400 bulk temperature C 420 440 Figure 3.19. Effect of heat flux and mass flow rate on heat transfer coefficient 58 25, R=3.15mm, P=24.34MPa, G=640kg/m2s, q=300kW/m2 O E 20 © 1 5 | 'o CD o o JD 10I c CO CTj 0 - . linear variation of heat flux - uniform heat flux 200 250 300 350 400 450 500 bulk temperature C Figure 3.20. Heat transfer coefficient for different variation of heat f l u x 59 4. MASS T R A N S F E R PREDICTIONS AND DISCUSSION 4.1 Concentration profile There are no measurements in the literature for comparison of concentration profiles. We have validated velocity and temperature profile already. There is similarity among these transport phenomena, and the same method should be suitable for concentration calculations. The concentration profiles at different heat fluxes and mass flow rates are shown in Figure 4.1 at P = 24.2MPa, R = 3.15mm, Tw = 386°C. The non-dimensional concentration (equation 2.34) is the concentration times friction velocity and divided by deposition mass flux. Concentration increases with an increase in molecule diameter as shown in Figure 4.2. The mass flux distribution across the tube has little influence on concentration profile. This is similar to the velocity profile. 4.2 Mass transfer coefficients In convection mass transfer problem, by analogy, we have the non-dimensional Sherwood number: Sh=f(Re,Sc) (4.1) Similar to the Swenson correlation for heat transfer, by substituting Prandtl number with Schmidt number, we can extend the Swenson correlation for mass transfer: Sh = ^ - = 0.00459(—)° 9 2 3 ( - ^ - ) ° 6 1 3 ( ^ ) ° 2 3 1 (4.2) Dw Mw PWDW ph Note that wall properties are used to evaluate Sc while integrated properties are used to calculate Pr in the original heat transfer correlation. This is arbitrary and other methods of evaluating Sc might be tried in the future. 60 In equation 4.2 is divided by equation 3.7, we can calculate the mass transfer coefficient from the heat transfer coefficient in the Swenson et al. correlation: ) 0.613 (4.3) K PwCpD] w where Cp = is the integrated heat capacity. As shown for P = 2A.34MPa, G =.6A0kg/m2s, q = 300kW I m2, R = 3.15mm , dm = 5.52e-\0m, in Figure 4.3, the mass transfer coefficient increases monotonically with bulk temperature. There is no peak in the pseudocritical region as for the heat transfer coefficient, but mass transfer coefficient increases significantly near the pseudocritical temperature. This is mainly attributed to water properties. From the Swenson et al. correlation, we can see that the mass transfer coefficient increases with a diffusion coefficient and decreases with density and viscosity. The diffusion coefficient increases with temperature at supercritical pressure in Figure 4.4. Density and viscosity decrease with temperature at supercritical pressure in Figure 1.1. Figure 4.4 shows diffusion coefficients at three different molecule diameters 2e-10, 5.52e-10, 60e-10 m. We can see that the diffusion coefficients decrease with increasing molecule diameter. We use a molecule diameter 5.52e-10 in our mass transfer simulations. We use the Stokes-Einstein equation to calculate the diffusion coefficient in mass transfer modeling. This equation is applied under normal conditions. Under supercritical condition, it is not clear how we should modify this formula? This is one major concern about the mass transfer model. The Swenson correlation shows the same tendencies that the mass transfer coefficient increases with temperature in Figure 4.3. The agreement is good between the Deissler-Van Driest model and the Swenson correlation in the subcritical region. The mass transfer coefficient by the Swenson correlation is about 50% higher than that computed 61 by the Deissler-Van Driest model near critical region and 6% higher after pseudocritical region. At P = 24.2MPa,G = 576kg/m2s,q = \07kWlm2,R = 3.15mm in Figure 4.5, the Deissler-Van Driest model agrees well with the Swenson correlation at molecule diameter =5.52e-10m. The Deissler-Van Driest model predicts the mass transfer coefficient 11% higher than the Swenson correlation in subcritical region and about 2.7% higher after the pseudocritical region. The mass transfer coefficients are comparable near the pseudocritical region. The Deissler-Van Driest model is higher than the Swenson correlation at molecule diameter = 2e-10m, but it is lower than the Swenson correlation at molecule diameter = 60e-10m. The mass transfer coefficient increases with mass flow rate in Figure 4.6 at P = 24.2MPa,R = 3.15mm . The higher average fluid flow velocity in the tube promotes transport of bulk fluid, therefore it results in higher mass transfer coefficient at larger mass flow rate. The mass transfer coefficient does not vary significantly with heat flux. In Figure 4.7 at P = 24.2MPa, G = 576kg/m2s, R = 3.15mm, q = \07kWIm2, the mass transfer coefficient decreases with the turbulent Schmidt number. This is obvious, since a smaller turbulent Schmidt number means a higher diffusion coefficient. Similar to heat transfer, we use a turbulent Schmidt number of 1.0 in the mass transfer modeling. The variation of mass flux across the tube affects the mass transfer coefficient slightly in Figure 4.8. The linear variation of mass flux and uniform mass flux across pipe are compared. The agreement is very good in the pseudocritical region. The mass transfer coefficient calculated from a linear variation of mass flux is 2.7% higher than that from uniform mass flux. We use a linear variation of mass flux in this modeling. From the above analysis and comparison, we believe that the Deissler-Van Driest calculation is applicable to mass transfer. It predicts heat transfer coefficients slightly better than the Swenson correlation, so we conclude that the Deissler-Van Driest model is 62 better for mass transfer modeling than the Swenson correlation. Nevertheless, extension of the Swenson correlation to predict mass transfer results in predictions quite close to the 1-D numerical model. Further comparison and measurements are needed to validate the calculated mass transfer coefficient. 63 R=3.15mm, P=24.2MPa s Tw=386C o i i i t i i i | i I I I i n i r——.T . , i i r T i ] r—„, , I T . , , T | - . q=400kW/m2, G=576kg/m2s - q=107kW/m2, G=576kg/m2s ^ " — q=107kW/m2, G=1200kg/m2s 10° 101 102 103 y+ ;ure 4.1. Effect of heat flux and mass flow rate on concentration profiles (molecule diameter = 5.52e-10m) 64 R=3.15mm, P=24.2MPa, G=576kg/m2s, q=107kW/m2, Tw=386C 140 120 100 80 + O 60 40 20 - . molecule diameter=2e-10m - molecule diameter=5.52e-10m — molecule diameter=60e-10m • • I L - J — L - J . J _l 1 • • • 10" 10 10 J y+ 10* Figure 4.2. Effect of molecule diameter on concentration profiles 65 R=3.15mm, P=24.34MPa, G=640kg/m2s, q=300kW/m2 E Es c CD "o 3=4 CD O O ^ 3 C/) c CO C/) 2 03 03 E 1 Deissler-Van Driest model - . Swenson 260 280 300 320 340 360 380 400 420 440 bulk temperature C Figure 4.3. Mass transfer coefficient for molecule diameter = 5.52e-10m 66 Figure 4.4. Diffusion coefficient of molecule using Stokes-Einstein relation 67 R=3.15mm, P=24.2MPa, G=576kg/m2s, q=107kW/m2 101 • 1 1 1 1 1 9h E E •*—> c CO o 'CD o o to Vi c CO 1— C/) w CO E dm=2e-10m - . Deissler-Van Driest model o Swenson dm=5.52e-10m - Deissler-Van Driest model + Swenson dm=60e-10m — Deissler-Van Driest model < Swenson < < < < < <_<<<<<< 1^ ' 320 340 360 380 400 420 440 bulk temperature C Figure 4.5. Mass transfer coefficient for different molecule diameter 68 R=3.15mm, P=24.2MPa - G=1200kg/m2s, q=107kW/m2 G=576kg/m2s, q=107kW/m2 -. G=576kg/m2s, q=250kW/m2 O 1 300 350 400 bulk temperature C 450 Figure 4.6. Effect of heat flux and mass flow rate on mass transfer coefficient (molecule diameter = 5.52e-10m) 69 Figure 4.7. Influence of turbulent Schmidt number on mass transfer coefficient (molecule diameter =5.52e-10m) 70 R=3.15mm, P=24.2MPa, G=576kg/m2s, q=107kW/m2 in E Es c CD 'o 3=4 0 O O £ 3 tf) C 03 i _ CO 2 CO CO E 1 - linear variation of mass flux -. uniform mass flux across pipe 320 340 360 380 400 bulk temperature C 420 440 Figure 4.8. Mass transfer coefficient for different variation of mass flux (molecule diameter = 5.52e-10m) 71 5. CONCLUSIONS AND R E C O M M E N D A T I O N S 5.1 Conclusions A heat and mass transfer model was developed based on concepts introduced by Deissler and Van Driest. The variation of heat flux and mass flux across the tube has a slight influence on the heat and mass transfer coefficients for fully developed turbulent flow. The linear variation of heat flux predicts heat transfer coefficient about 6.7% higher than uniform heat flux across the pipe, and a 5.7% higher peak heat transfer coefficient. The agreement of the mass transfer coefficient with different variations of mass flux is very good, especially near the pseudocritical temperature. The mass transfer coefficient is estimated with a linear variation of mass flux to be 2.7% higher than that estimated with uniform mass flux across the pipe. The heat and mass transfer coefficients decrease with increasing turbulent Prandtl and Schmidt numbers. The peak heat transfer coefficient increases significantly with a decrease in heat flux. The mass transfer coefficients also change with heat flux. Both heat and mass transfer coefficients increase with mass flow rate. We set the turbulent Prandtl and Schmidt numbers to 1.0 in these simulations. The heat transfer coefficient increases with bulk temperature, and reaches a peak value near the pseudocritical temperature. Beyond this point, the heat transfer coefficient decreases with bulk temperature. This is mainly due to high heat capacity at the pseudocritical temperature. In comparison with Yamagata's experiment, the Deissler-Van Driest model is better than the Swenson et al. correlation except in the subcritical region. The maximum heat transfer coefficient from the Swenson et al. correlation is reached before the pseudocritical temperature. The Deissler-Van Driest model agrees well with the UBC experimental data. The difference between measurement and calculation might be due to a variety of causes. Different physical properties used in earlier works produce different results. Also, fully developed turbulent flow is assumed in this study; entrance length effects may have affected some of the experiments. 72 Buoyancy is neglected, and axisymmetric pipe flow is assumed in this study. According to Bazargan's measurement in the UBC SCWO pilot plant, in many cases, buoyancy does play a role in supercritical water flow, especially at low mass flow rates and high heat fluxes. If buoyancy is considered for supercritical water flow in horizontal pipe, the axisymmetric pipe flow assumption is not valid anymore. Although we attempted to select buoyancy-free experiments, there is no proof that the experiments from UBC or other groups are truly free of buoyancy effects. Turbulent flow was a big challenge in science and technology during the last century and continues to be in this century. Which turbulence model is the best one? There is no obvious answer. It depends on which problem we study. The simple Prandtl mixing length model is good for constant property pipe flow. Is it good for variable property flow? Since complicated variations of water properties occur near the critical region, we use algebraic turbulence model to make this problem easy to study. The mass transfer coefficient increases steadily and monotonically with bulk temperature. There is no peak in the pseudocritical region. The mass transfer coefficient changes dramatically near the pseudocritical temperature. This mass transfer coefficient variation is due to changing water properties. Under supercritical pressures, the diffusion coefficient increases with temperature, while density and viscosity decrease with temperature. These properties vary significantly in the pseudocritical region. The total effect of a higher diffusion coefficient , a lower density and viscosity speeds the diffusion of molecules and transport of bulk fluid. This increases the mass transfer coefficient. The Deissler-Van Driest model and Swenson et al. correlation show the same tendency for mass transfer coefficient. Near the pseudocritical region, the mass transfer coefficient from the Swenson correlation is higher by about 50% than the Deissler-Van Driest model. The difference decreases to 6% after the pseudocritical region. One uncertainty in mass transfer modeling is the diffusion coefficient. The Stokes-Einstein formula is used under normal conditions. Is it applicable under supercritical conditions? And what modification we should have? 73 Although we don't have direct measurements for mass transfer coefficients now, based on an analogy between heat and mass transfer, comparison and analysis, we believe that the Deissler-Van Driest model calculation is applicable to mass transfer. It predicts the heat transfer coefficients very well, so predicted mass transfer coefficients should be acceptable. The Deissler-Van Driest model may be better than the Swenson et al. correlation in heat transfer modeling. 5.2 Recommendations New measurements of mass transfer and further comparison are needed to improve the mass transfer model. In order to study the effect of buoyancy and thermal entrance length, a three-dimensional model should be used. The Navier-Stokes equations, energy and concentration equation can be solved using commercial software. We can compare the wall temperature between simulation results and measurement, predict flow field, heat transfer, mass transfer and provide guidelines for experiments and design. However, this is still a very difficult problem. Aside from the need for an enormous computational mesh to resolve the wall region, it is not clear how to model the turbulence in the wall region or the diffusion of molecules in near-critical conditions. 74 Appendix A Computer Code SCWATER.F To simulate heat and mass transfer in supercritical water *************************************************************** PROGRAM main Program f o r S i m u l a t i n g Mass Transfer to SCW i n Tubes using the Deissler-Van D r i e s t model with v a r i a b l e p r o p e r t i e s w r i t t e n by Dr Rogak i n FORTRAN 77, modified by L i Xu *************************************************************** CALL LAMINAR CALL TURBULENT CALL SCWATER STOP END *************************************************************** subroutine laminar compare laminar v e l o c i t y p r o f i l e s between p a r a b o l i c p r o f i l e and using "Deissler-Van D r i e s t " i n t e g r a t i o n f o r constant p r o p e r t i e s *************************************************************** IMPLICIT NONE INTEGER n,m,i PARAMETER (n=500) REAL* 8 uu(n),ybar(n),y(n),u(n),ul(n) REAL*8 RADIUS,G,AREA REAL*8 UFRIC,RP,QW,DW,CPW,TW,P,KW,VIScW,TBULK, mflow REAL*8 HW,PRW,MFL0W1,FLOWERR, CPAVG REAL*8 DBULK,HBULK,CPBULK,VISBULK,KBULK,PRBULK REAL*8 UBAR,REY,FRIC,SHEARW,YFRIC,PI 0PEN(3,FILE='ulaminar') OPEN(10,FILE='parabolic') PI = 4\*ATAN(1.0) G=120 RADIUS=3.15e-3 DW=998 Ubar=G/DW mflowl=G*PI*RADIUS**2 viscW=l.003e-3 Rey=ubar*2.*RADIUS*DW/viscW fric=64./Rey shearw=fric*DW*Ubar**2/8. ufric=(shearw/DW)**0.5 yfric=viscW/DW/ufric RP=RADIUS/yfric CALL LPROFILE(RP,m,ybar, uu) do i=l,m y ( i ) = y b a r ( i ) * y f r i c u ( i ) = u u ( i ) * u f r i c ul(i)=2.*ubar*(1-((RADIUS-y(i))/RADIUS)**2) w r i t e ( 3 , 6 0 ) i , y ( i ) , u ( i ) w r i t e ( 1 0 , 6 0 ) i , y ( i ) , u l { i ) 75 enddo 60 format(16,2e20.5) RETURN END SUBROUTINE LPROFILE(RP,m,ybar,uu) * This r o u t i n e i n t e g r a t e s away from the w a l l i n t o the bulk * c o n d i t i o n s f o r laminar flow and constant p r o p e r t i e s ************************************************************* IMPLICIT NONE INTEGER n,m PARAMETER (n=500) REAL* 8 uu(n),ybar(n) REAL*8 UP, YP, VISBAR,DYP,DUDY,RP REAL*8 T,H,CP,D,VIS,K,PR,PI,CORR PI=4.*ATAN(1.0) c we i n t e g r a t e out from w a l l , i n non-c dimensional coordinates (y+=YP; T+=TP i s nomenclature) UP=0 YP=0 DYP=0.1 VISBAR=1 m=0 DO WHILE(YP.LT.RP) YP=YP+DYP IF(YP.GT.RP) THEN DYP=DYP-(YP-RP) YP=RP END IF m=m+l CORR=l-YP/RP dudy=CORR/VISBAR UP=UP+DUDY*DYP uu(m)=up ybar(m)=yp END DO RETURN END ************************************************* subroutine t u r b u l e n t * compare u n i v e r s a l v e l o c i t y p r o f i l e with Deissler-Van D r i e s t model * f o r constant p r o p e r t i e s ******************************************************************* IMPLICIT NONE INTEGER n,m,i PARAMETER (n=5000) REAL*8 uu(n),ybar(n),uplus(n),uvan(n) REAL*8 RADIUS,G,AREA REAL*8 UFRIC,RP,QW,DW,CPW,TW,P,KW,VIScW,TBULK,mflow REAL*8 HW,PRW,MFLOW1,FLOWERR,CPAVG REAL*8 DBULK,HBULK,CPBULK,VISBULK, KBULK, PRBULK 76 REAL*8 UBAR,REY,FRIC,SHEARW,YFRIC,PI 0PEN(11,FILE='uturbulent') 0PEN(12,file='universal') PI=4.*ATAN(1.0) G=120000. RADIUS=3.15e-3 DW=998 Ubar=G/DW mf lowl=G*PI*RADIUS**2 viscW=l.003e-3 Rey=ubar*2.*RADIUS*DW/viscW write(*,*)'Rey', Rey fric=0.316/Rey**.25 shearw=fric*DW*Ubar**2/8. C loop to f i n d c o r r e c t mass flow 80 ufric=(shearw/DW)**0.5 yfric=viscW/DW/ufric RP=RADIUS/yfric CALL TPROFILE(UFRIC,RP,MFLOW,m,ybar,uu, uvan) c convert from non-dimensional forms mflow=mflow*yfric**2*DW*UFRIC FLOWERR=MFLOW/MFLOWl IF(ABS(FLOWERR-1.).GT.0.01) THEN SHEARW=SHEARW/FLOWERR* * 2. GOTO 8 0 END IF do i=l,m w r i t e ( 1 1 , 6 0 ) i , y b a r ( i ) , u u ( i ) w r i t e ( 1 2 , 6 0 ) i , y b a r ( i ) , u v a n ( i ) enddo 60 format(I6,2E20.5) r e t u r n END ************************************************************** SUBROUTINE TPROFILE(UFRIC,RP,MFLOW,m,ybar,uu,uvan) * This r o u t i n e i n t e g r a t e s away from the w a l l i n t o the bulk * c o n d i t i o n s . The v e l o c i t y p r o f i l e i s determined * by the f r i c t i o n v e l o c i t y . ************************************************************** IMPLICIT NONE INTEGER n,m PARAMETER (n=5000) REAL*8 uu(n),ybar(n),uvan(n) REAL*8 UFRIC,RP,Q W , T W , P REAL*8 DW,hw,CPW,KW,VISW,PRW,TBULK,mflow REAL* 8 UP,YP,TP,VISBAR,RHOBAR,CONDBAR,CPBAR REAL*8 MIXL, DUDY, DTDY,DYP,DMASS,HFLOW REAL*8 T,H,CP,D,VIS,K,PR,HAVG,TOFH REAL*8 PI,upv,dudyv,CORR PI=4.*ATAN(1.0) we i n t e g r a t e out from w a l l , i n non-77 dimensional coordinates (y+=YP; T+=TP i s nomenclature) UP=0 upv=0 YP=0 DYP=0.1 VISBAR=1 RH0BAR=1 mflow=0 m=0 DO WHILE(YP.LT.RP) YP=YP+DYP IF(YP.GT.RP) THEN DYP=DYP-(YP-RP) YP=RP END IF m=m+l MIXL=0.4*(l-EXP(-YP/2 6))*YP CORR=l-YP/RP CORR=l. DUDY=(SQRT(VISBAR**2.+4*MIXL**2.*RHOBAR*CORR)-VISBAR) C / (2*MIXL**2.*RHOBAR) dudyv=2./(l+(l+4.*0.4**2*yp**2*(1-exp(-yp/26.))**2)**0 . 5) upv=upv+dudyv*DYP UP=UP+DUDY*DYP uu(m)=up uvan(m)=upv ybar(m)=yp The mass flow (non-dimensional) i n the l a s t annular c r o s s - s e c t i o n a l element i s DMASS=2.*PI*rhobar*(RP-YP)*DYP*UP mflow=mflow+DMASS DYP=DYP*1.01 END DO RETURN END **************************************************************** subroutine SCwater s i m u l a t i n g heat and mass t r a n s f e r i n s u p e r c r i t i c a l water using Deissler-Van D r i e s t model ************************************************* IMPLICIT NONE INTEGER NIN,NOUT,IT,NSTEP,n,m,i PARAMETER (n=5000) REAL*8 uu(n) , t t ( n ) , cbar(n),ybar(n),uplus(n) REAL*8 RADIUS,G,TMAX,DELTAT,AREA REAL*8 UFRIC,RP,QW,DW,CPW,TW,P,KW,VIScW,TBULK,mflow REAL*8 HW,PRW,MFLOW1,FLOWERR,CPAVG REAL*8 DBULK,HBULK,CPBULK,VISBULK,KBULK,PRBULK REAL*8 UBAR,REY,FRIC,SHEARW,YFRIC,heatcoef v a r i a b l e s r e l a t e d to mass t r a n s f e r : REAL*8 MASSCOEF,MCOEFl,MCOEF2,MCOEF3,DIFFW,SCW,QMASS REAL*8 CBULK, DIFFUSE, PI, Pr,Nub,Prt,Set,Reys REAL*8 TK,HK,HS,Q,Hdittus,Re,Hd,HSwenson,MS,MD REAL*8 Nu,Tguass input/output u n i t numbers must be passed to the eqtest r o u t i n e s 78 COMMON/COUT/NIN,NOUT OPEN(9, FILE='read') OPEN(13,FILE='heat' ) OPEN(14,FILE='mass') OPEN(15,FILE='upyp') OPEN(16,FILE='typ') OPEN(17,FILE='cyp') PI=4.*ATAN(1.0) NIN = 2 NOUT = 6 READ(9,*) P r t , S c t READ(9,*) RADIUS RADIUS=RADIUS/1000 . AREA=PI*RADIUS**2 . • READ(9,*) P READ(9,*) G READ(9,*) QW QW=QW*1000. READ(9,*) TW read(9,*)Tmax TW=TW+273. Tmax=Tmax+273. w r i t e (13,*)' Tw Tb h hs $ hd Hb Re' write(14,20) 20 FORMAT( ' Tw(C) Tb(C)', C ' h M Ms Hb Re' ) write(*,20) QMASS=1. DO WHILE(TW.LE.TMAX) IF((TW-27 3).LE.3 60.)then DELTAT=5. ELSEIF((TW-273).LT.38 0)THEN DELTAT=2. ELSEIF((TW-273).LT.420)THEN DELTAT=1. ELSEIF((TW-27 3).LT.4 50.)THEN DELTAT=2. ELSE DELTAT=5. ENDIF TW=TW+DELTAT MFLOWl=G*AREA CALL PROP (TW, P, DW, HW, CPW, VISCW, KW, PRW) DIFFW=DIFFUSE(TW,VISCW) SCW=VISCW/DIFFW/DW Ubar=G/DW Rey=ubar*2.*RADIUS*DW/viscW fric=0.316/Rey**.25 shearw=fric*DW*Ubar**2/8. C INPUT / OUTPUT ON PC: C loop to f i n d c o r r e c t mass flow 80 ufric=(shearw/DW)**0.5 yfric=viscW/DW/ufric RP=RADIUS/yfric 79 CALL PROFILE(UFRIC,RP,QW,QMASS,TW,P,TBULK, CBULK, MFLOW, $ m,ybar,uu,tt,cbar,Prt,Set) c convert from non-dimensional forms mflow=mflow*yfric**2*DW*UFRIC FLOWERR=MFLOW/MFLOWl c WRITE(*,*) 'FLOWERR=',FLOWERR IF(ABS(FLOWERR-1.).GT.0.001) THEN SHEARW=SHEARW/FLOWERR* * 2. GOTO 8 0 END IF CALL PROP(TBULK,P,DBULK,HBULK,CPBULK,VISBULK, KBULK, PRBULK) Ubar=mflow/DBULK/AREA Rey=ubar*2.*RADIUS* DBULK/visBULK HEATCOEF=QW/(TW-TBULK) MASSCOEF=QMASS/CBULK CPAVG=(HW-HBULK)/(TW-TBULK) Pr=CPW*viscW/KW Nub=HEATCOEF*RADIUS*2/KBULK TK=TBULK HK=HBULK C i t e r a t i o n to c a l c u l a t e the heat t r a n s f e r c o e f f i c i e n t from Swenson c c o r r e l a t i o n Reys=G*2.*RADIUS/visCW 30 Nu=0.00459*Reys**0.923*((HW-HBULK)/(TW-TBULK)* + viscW/KW)**0.613*(DW/DBULK)**0.231 HSwenson=Nu*KW/2./RADIUS Tguass=TBULK TBULK=(TBULK+TW-QW/HSwenson)12 . C WRITE(*,*)'Swenson',TBULK IF (ABS(Tguass-TBULK).GT.0.1)THEN CALL PROP(TBULK,P,DBULK,HBULK,CPBULK,VISBULK,KBULK,PRBULK) GOTO 30 ENDIF CPAVG=(HW-HBULK)/(TW-TBULK) MS=HSwenson*DIFFW/kw*(KW/(DW*CPAVG*DIFFW))**0.613 c i t e r a t i o n t o c a l c u l a t e the heat t r a n s f e r c o e f f i c i e n t from D i t t u s -c B o e l t e r c o r r e l a t i o n 40 Re=G*2.*RADIUS/VISBULK Hdittus=0.023*Re**0.8*(CPBULK*VISBULK/KBULK)**0.4* $ KBULK/2./RADIUS Tguass=TBULK TBULK=(TBULK+TW-QW/Hdittus)12 . c w r i t e ( * , * ) ' D i t t u s ' , T b u l k IF (ABS(TBULK-Tguass).GT.0.1)THEN CALL PROP(TBULK,P,DBULK,HBULK,CPBULK, VISBULK, KBULK, PRBULK) GOTO 4 0 ENDIF wr i t e (13,100) TW-27 3.,TK-27 3.,heatcoef/1000., $ HSwenson/1000.,Hdittus/1000., HK/1000., REY write(14,100) TW-27 3.,TK-27 3.,heatcoef/1000., c MASSCOEF*1000.,MS*1000., HK/1000., REY write(*,100) TW-273.,TK-27 3.,heatcoef/1000., c MASSCOEF*1000.,MS*1000.,HK/1000., REY 100 FORMAT(2F8.1,3F10.3,2F10.0) 80 ENDDO do i=l,m write ( 1 5 , 6 0 ) i , y b a r ( i ) , uu ( i ) write(16, 6 0 ) i , y b a r ( i ) , t t (i) w r i t e ( 1 7 , 6 0 ) i , y b a r ( i ) , c b a r (i) enddo 60 format(16,2e20.5) w r i t e (*,*) 'meshsize',m r e t u r n END SUBROUTINE PROFILE(UFRIC,RP,QW,QMASS,TW,P,TBULK, C B U L K , MFLOW, $ m, ybar, uu,tt,cbar,Prt,Set) C This r o u t i n e uses Deissler-Van D r i e s t ideas to i n t e g r a t e away from the c w a l l i n t o the b u l k ' c o n d i t i o n s . The v e l o c i t y p r o f i l e i s determined c by the f r i c t i o n v e l o c i t y , which i s an input to t h i s r o u t i n e , c F l u i d p r o p e r t i e s change through the p r o f i l e , which means that we must c c a l l the thermodynamics r o u t i n e to evaluate p r o p e r t i e s . The ro u t i n e c t h e r e f o r must use dimensionless AND dimensional v a r i a b l e s . (2*********************************************** ************** IMPLICIT NONE INTEGER n,m PARAMETER (n=5000) REAL* 8 uu(n),tt(n),cbar(n) ,ybar(n) REAL*8 UFRIC,RP,QW,TW,P REAL*8 DW,hw,CPW,KW,VISW,PRW,TBULK,mflow REAL* 8 UP,YP,TP,VISBAR,RHOBAR,CONDBAR,CPBAR REAL*8 MIXL, DUDY, DTDY,DYP,DMASS,HFLOW REAL*8 T,H,CP,D,VIS,K,PR,HAVG,TOFH c v a r i a b l e s r e l a t e d to mass t r a n s f e r REAL*8 QMASS,CBULK,DIFF,CONC,CONCP REAL*8 CFLOW,DIFFBAR,SCW,DIFFW,DIFFUSE,DCDY REAL*8 PI,Prt,Set,CORR PI=4.*ATAN(1.0) c we i n t e g r a t e out from w a l l , i n non-c dimensional coordinates (y+=YP; T+=TP i s nomenclature) UP=0 YP=0 TP=0 CONCP=0 DYP=0.1 VISBAR=1 RHOBAR=l CONDBAR=l DIFFBAR=1. CPBAR=1 CALL PROP (TW, P, DW, HW, CPW, VISW, KW, PRW) c For mass t r a n s f e r , need Schmidt number DIFFW=DIFFUSE(TW,VISW) SCW=VISW/DW/DIFFW mflow=0 HFLOW=0 CFLOW=0 c write(*,50) ufric,RP,qw,dw,prw, tw, p 50 format(7G11.3) 81 m=0 DO WHILE(YP.LT.RP) YP=YP+DYP IF(YP.GT.RP) THEN DYP=DYP-(YP-RP) YP=RP END IF m=m+1 MIXL=0.4*(l-EXP(-YP/2 6))*YP CORR=l-YP/RP c CORR=l. DUDY=(SQRT(VISBAR**2. +4*MIXL**2.*RHOBAR*CORR)-VISBAR) C / (2*MIXL**2.*RHOBAR) UP=UP+DUDY*DYP uu(m)=up ybar(m)=yp DTDY=CORR/(CONDBAR/PRW + RHOBAR*CPBAR*MIXL* *2.* DUDY/Prt) DCDY=CORR/(DIFFBAR/SCW + MIXL**2.*DUDY/Sct) TP=TP+DTDY*DYP CONCP=CONCP+DCDY*DYP C re-evaluate p r o p e r t i e s at the new temperature T=TW-TP*qw/(DW*CPW*Ufric) CONC=CONCP*QMASS/Ufric tt(m)=T-273 cbar(m)=CONCP IF(T.gt.360) then C eqtest w i l l work c otherwise assume p r o p e r t i e s need not be r e c a l c u l a t e d CALL PROP(T,P,D,H,CP,VIS, K,PR) VISBAR=VIS/VISW RHOBAR=D/DW CONDBAR=K/KW CPBAR=CP/CPW DIFF=DIFFUSE(T, VIS) • DIFFBAR=DIFF/DIFFW END IF C The mass flow (non-dimensional) i n the l a s t annular c r o s s - s e c t i o n a l c element i s DMASS=2.*PI*rhobar*(RP-YP)*DYP*UP mflow=mflow+DMASS C enthalpy c a r r i e d i n t h i s element of area i s HFLOW=HFLOW+DMASS*H CFLOW=CFLOW+DMASS*CONC DYP=DYP*1.01 END DO HAVG=H FLOW/MFLOW CBULK=CFLOW/MFLOW TBULK=TOFH(P,HAVG, T) C We expect that the c e n t e r l i n e temperature w i l l be close to the bulk c temperature, and only a small c o r r e c t i o n i s needed to f i n d the true Cbulk c temperature. RETURN END Q* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * -k -k -k * * * * * * * * * * * * * * * * * * * * * * REAL*8 FUNCTION TOFH(P,H,T) Q* ************************************************************* * 82 C This f u n c t i o n returns T(H,P) given an i n i t i a l guess f o r T c P i n MPa, T i n K, H i n J/kg REAL*8 P,H,T REAL* 8 DS,D,HGUESS,BDENS,DENS,HB,CP,CPB,DT 10 DS = BDENS(T,P,0) D = DENS (P, T, DS, 1 . D-6) HGUESS=HB(T,D)*1000. CP=CPB(T,D)*1000. c w r i t e ( * , * ) H,Hguess,T,CP DT=(H-Hguess)/CP T=T+DT IF(ABS(DT).GT.0.2) GOTO 10 TOFH=T RETURN END Q*********************************************************************** SUBROUTINE PROP(T,P,D,H,CP,VIS, K, PR) Q*********************************************************************** C This r o u t i n e c a l l s the eqtest and other r o u t i n e s to get c the usual thermodynamic p r o p e r t i e s as a f u n c t i o n of T and P c A l l are i n SI REAL*8 T,P,D,H,CP,VIS,K,PR REAL* 8 DS,VISC,THK DOUBLEPRECISION BDENS,DENS,VISCOSITY,TCONDUCT,CPB, HB DS = BDENS(T,P,0) D = DENS(P,T,DS,1.D-6) c a l l TRANSPORT(T,P,VIS,K) H=HB(T,D)*1000. CP=CPB(T,D)*1000. PR=vis*cp/K RETURN END REAL*8 FUNCTION DIFFUSE(T,VISC) c returns d i f f u s i o n c o e f f i c i e n t of a molecule using c S t o k e s - E i n s t e i n , with a l l v a r i a b l e s i n SI. Q* ******************************************************************* * REAL*8 DM,DIFF,T,VISC,PI,BOLTZ PARAMETER (PI = 3.14159265359, BOLTZ = 1.38E-23) C diameter of the molecule i n meters DM = 5.52e-10 DIFF = BOLTZ*T/(3*PI*VISC*DM) DIFFUSE = DIFF RETURN END ******************************************************************** * The f o l l o w i n g subroutine i s from Majid Bazargan f o r c a l c u l a t i o n of * v i s c o s i t y and thermal c o n d u c t i v i t y . The v i s c o s i t y i s from the IAPS * Formulation 1985 f o r the v i s c o s i t y of ordinary water substance, and * the thermal c o n d u c t i v i t y i s from the IAPS Formulation 1985 f o r the * v i s c o s i t y of ordi n a r y water substance, the tenth i n t e r n a t i o n a l 83 * conference on the p r o p e r t i e s of steam, Moscow, USSR., September, 1 ************************************************* ****************************************************************** subroutine TRANSPORT(T, P, VISC, THK) ************************************* C C GIVEN T AND P, THIS ROUTINE CALCULATES VISCOSITY AND C THERMAL CONDUCTIVITY OF WATER, "eqsub.f" OR "iapws957.f" C MUST BE ATTACHED TO THE MAIN PROGRAM USED IN CONJUNCTION C WITH THIS ROUTINE. C IMPLICIT REAL*8 (A-H,0-Z) IMPLICIT INTEGER (I-N) REAL*8 ACOEF(4),BCOEF(6,7) ,ECOEF(4) ,FC0EF(5,6) C COMMON/CCPEQ/TCEQ, PCEQ, DCEQ COMMON/CSUB2/R,XMOL,TC,PC,DC COMMON/CNORM/TNORM, DNORM COMMON/CSUB3/TTR,PTR,DLTR,DVTR,TBOYL,PBOYL,DLB,DVB COMMON/COUT/NIN,NOUT C BRO=BDENS(T,P,0) RO=DENS(P,T,BRO,1.D-6) C TSTAR = 647.27 ROSTAR = 317.763 PSTAR = 22.115E06 VISTAR = 55.071E-6 THKSTAR = 0. 4 94 5 C TBAR = T/TSTAR ROBAR = RO / ROSTAR PBAR = P/PSTAR Ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c C NEWTON'S FORMULA IS USED IN THE FOLLOWING TO CALCULATE C del(RHO)/del(P) C DELP =0.1 PO = P PI = P + DELP P2 = P + 2*DELP P3 = P + 3*DELP P4 = P + 4*DELP RO0 = RO BR01=BDENS(T,PI,0) R01=DENS(PI,T,BROl,1.D-6) BR02=BDENS(T,P2,0) R02=DENS(P2,T,BR02,l.D-6) BR03=BDENS(T,P3,0) R03=DENS(P3,T,BR03,1.D-6) BR04=BDENS(T,P4,0) R04=DENS(P4,T,BR04,l.D-6) C DELRO0=ROl-RO0 DELR01=R02-R01 84 DELR02=R03-R02 DELR03=R04-R03 C DEL2RO0=DELROl-DELRO0 DEL2R01=DELR02-DELR01 DEL2R02=DELR03-DELR02 C DEL3RO0=DEL2ROl-DEL2RO0 DEL3R01=DEL2R02-DEL2R01 C DEL4 RO0=DEL3ROl-DEL3RO0 C DRODP = (DELRO0-(0.5* DEL2RO0) + (DEL3RO0/3)-(DEL4RO0/4) )/DELP XTBAR = ROBAR* PSTAR*DRODP/ROSTAR Ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c ACOEF(l) = 1.000000 C ACOEF(2) = 0. 978197 ACOEF(3) = 0 . 579829 ACOEF(4) =- 0 . 202354 BCOEF(1, 1) = 0 . 5132047 BCOEF(1, 2) = 0 . 2151778 BCOEF(1, 3) =- 0 . 2818107 BCOEF(1, 4) = 0. 1778064 BCOEF(1, 5) =- 0 . 0417661 BCOEF(1, 6) = 0 . 0000000 BCOEF(1, 7) = 0. 0000000 BCOEF(2, 1) = 0. 3205656 BCOEF(2, 2) = 0 . 7317883 BCOEF(2, 3) =- -1. 070786 BCOEF(2, 4) = 0. 4605040 BCOEF(2, 5) = 0 0000000 BCOEF(2, 6) =--0 01578386 BCOEF(2, 7) = 0 0000000 BCOEF(3, 1) = 0 0000000 BCOEF(3, 2) = 1 2410440 BCOEF(3, 3) =--1 263184 BCOEF(3, 4) = 0 2340379 BCOEF(3, 5) = 0 0000000 BCOEF(3, 6) = 0 0000000 BCOEF(3, 7) = 0 0000000 BCOEF(4, 1) = 0 0000000 BCOEF(4, 2) = 1 476783 BCOEF(4 3) = 0 0000000 BCOEF(4 4) = -0 4924179 BCOEF(4 5) = 0 .1600435 BCOEF(4 6) = 0 . 0000000 BCOEF(4 7) = -0 . 003629481 BCOEF(5 1) = -0 .7782567 BCOEF(5 2) = 0 . 0000000 BCOEF(5 3) = 0 .0000000 BCOEF(5 4) = 0 . 0000000 BCOEF(5 5) = 0 . 0000000 BCOEF(5 6) = 0 . 0000000 BCOEF(5 7) = 0 . 0000000 BCOEF(6 1) = 0 . 1885447 85 BC0EF(6, 2) 0. ooooooo BC0EF(6, 3) = 0. ooooooo BC0EF(6, 4) = 0. ooooooo BC0EF(6,5) = 0 . ooooooo BCOEF(6, 6) = 0. ooooooo BCOEF(6,7) = 0 . ooooooo ECOEF(l) = 1. 000000 ECOEF(2) = 6. 978267 ECOEF(3) = 2 . 599096 ECOEF(4) =--0 998254 FCOEF(1, 1) = 1. 3293046 FCOEF(1, 2) =--0 . 40452437 FCOEF(1, 3) = 0 24409490 FCOEF(1, 4) = 0 018660751 FCOEF(1, 5) =--0 12961068 FCOEF(1, 6) = 0 044809953 FCOEF(2, 1) = 1 7018363 FCOEF(2, 2) =--2 2156845 FCOEF(2, 3) = 1 6511057 FCOEF(2, 4) =--0 76736002 FCOEF(2, 5) • = 0 37283344 FCOEF(2, 6) =--0 11203160 FCOEF(3, 1) = 5 2246158 FCOEF(3, 2) = • -10.124111 FCOEF(3, 3) = 4 9874687 FCOEF(3, 4) = -0 27297694 FCOEF(3, 5) = -0 43083393 FCOEF(3, 6) = 0 13333849 FCOEF(4, 1) = 8 .7127675 FCOEF(4, 2) = -9 . 5000611 FCOEF(4, 3) = 4 . 3786606 FCOEF(4, 4) = -0 .91783782 FCOEF(4 5) = 0 .ooooooo FCOEF(4 6) = 0 .ooooooo FCOEF(5 1) = -1 .8525999 FCOEF(5 2) = 0 .93404690 FCOEF(5 3) = 0 .ooooooo FCOEF(5 4) = 0 . ooooooo FCOEF(5 5) = 0 . ooooooo FCOEF(5 6) = 0 . ooooooo C CALCULATIONS FOR VISCOSITY; VISC=VISCO *VISCI*VISC2 C VO = 0.0 DO 21 k=l,4 KK=K-1 VO = VO + ACOEF(K)/(TBAR**KK) 21 CONTINUE VISCO = (SQRT(TBAR))/VO C VI = 0.0 DO 23 1=1,6 11=1-1 DO 22 J=l,7 86 JJ=J-1 VI = VI + BCOEF(I, J ) * ( (1/TBAR -1)* * 11) *{ (ROBAR-1)* *JJ) 22 CONTINUE 23 CONTINUE VISC1 = EXP(ROBAR*VI) C VISC2=1 c XBART=ROBAR* DRODP c IF( (TBAR.GE.0.997.AND.TBAR.LE . 1.0082) .OR. c +(ROBAR:GE.0.775.AND.ROBAR.LE.1.290))THEN c IF(XBART.GE.21.93)THEN c VISC2=0.992*(XBART**0 . 0263) C ENDIF c ENDIF C VISC=VISTAR*VISC0*VISC1*VISC2 Ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c C CALCULATIONS FOR THEMAL CONDUCTIVITY; THK=THK0*THK1 + THK2 C THO =0.0 DO 31 L=l,4 LL=L-1 THO = THO + ECOEF(L)/(TBAR**LL) 31 CONTINUE THK0 = SQRT(TBAR)/THO C TH1 = 0.0 DO 33 M=l, 5 MM=M-1 DO 32 N=l,6 NN=N-1 TH1 = TH1' + FCOEF(M,N)*((1/TBAR -1)**MM)*((ROBAR-1)**NN) 32 CONTINUE 33 CONTINUE THK1 = EXP(ROBAR*THl) C C TO CALCULATE THK2, de l ( P ) / d e l ( T ) AT CONSTANT DENSITY IS REQUIRED C LIKEWISE del(RO)/del(P) , NEWTON'S FORMULA IS USED: C DELT =0.1 TO = T T l = T+DELT T2 = T+2*DELT ; T3 = T + 3* DELT T4 = T+4 *DELT P0 = P PI = PB(Tl,RO) P2 = PB(T2/RO) P3 = PB (T3, RO) P4 = PB(T4,RO) C DELP0=P1-P0 DELP1=P2-P1 DELP2=P3-P2 DELP3=P4-P3 C 87 DEL2P0=DELP1-DELP0 DEL2P1=DELP2-DELP1 DEL2P2=DELP3-DELP2 DEL3P0=DEL2P1-DEL2P0 DEL3P1=DEL2P2-DEL2P1 C DEL4P0=DEL3P1-DEL3P0 C DPDT = (DELPO-(0.5*DEL2P0)+(DEL3P0/3)-(DEL4 PO/4) )/DELT C c ONE=(0.0013848*SQRT(ROBAR))/(VISC0*VISC1) ONE=(0.00138 4 8*SQRT(ROBAR))/VISC/10e-6 TWO=((TBAR/ROBAR)* (DPDT*TSTAR/PSTAR))**2 THREE=XTBAR** 0.4 67 8 FOUR=EXP((-18.66*(TBAR-1)**2)-(ROBAR-1)**4) THK2=ONE*TWO*THREE* FOUR THK=THKSTAR*(THK0*THK1+THK2) RETURN END 88 Appendix B Sample input and output files Sample input file <read> Input Explanations 1. 3.15 24 . 2 576 107 380 410 1 / / P r t Set //R(mm) //P(MPa) //G(kg/m2s) //q(kW/m2) //Tw(C) //Tmax(C) Output file <ulaminar> //laminar velocity output for constant properties Output file <parabolic> // laminar analytic parabolic velocity for constant properties Output file <uturbulent> // turbulent velocity for constant properties Output file <universal> // Van-Driest universal velocity profile for constant properties Output file <upyp> // velocity profile with variable properties Output file <typ> // temperature profile with variable properties Output file <cyp> // concentration profile with variable properties Sample output file <heat> Tw: wall temperature (C) Tb: bulk temperature (C), h: heat transfer coefficient using Deissler-Van Driest model (kW/m2C) hs: heat transfer coefficient from Swenson correlation (kW/m2C) hd: heat transfer coefficient from Dittus-Boelter correlation (kW/m2C) Hb: bulk enthalpy (kJ/kg) Re: Reynolds number Tw Tb h hs hd Hb Re 381.0 375 5 19.433 20. 530 13.792 1873. 64679 382 . 0 377 5 23.795 29. 600. 14.466 1911. 67917 383.0 379 6 31.127 36. 070 15.557 1969. 73491 384.0 381 1 37.428 36. 444 18.538 2050. 82092 385 . 0 382 0 35 .518 32 . 162 37.927 2127 . 90871 386.0 382 5 30.980 27 . 636 33.741 2186. 97520 387 . 0 383 1 27.141 23. 994 30.123 2232 . 102299 388 . 0 383 6 24.203 21. 422 26.454 2269. 105981 389.0 384 1 21.931 19. 268 23.868 2300. 108924 89 390 0 384 7 20 130 . 17 576 21 496 2328 . 111349 391 0 385 3 18 688 16 214 19 868 2352 . 113401 392 0 385 9 17 476 15 097 18 526 2374 . 115140 393 0 386 5 16 464 14 163 17 398 2394 . 1167 55 394 0 387 1 15 587 13 369 16 436 2413. 118032 395 0 387 8 14 827 12 724 15 602 2430 . 119163 396 0 388 4 14 161 12 126 14 877 2446. 120171 397 0 389 1 13 573 11 600 14 237 2461. 121068 398 0 389 8 13 049 11 133 13 668 2475. 121876 399 0 390 5 12 579 10 715 13 158 2489. 122606 400 0 391 2 12 153 10 339 12 631 2502 . 123267 401 0 391 9 11 767 . 9 998 12 216 2514 . 123867 402 0 392 6 11 414 9 687 11 837 2526. 124409 403 0 393 4 11 090 9 402 11 489 2538 . 124903 404 0 394 1 10 791 9 141 11 170 2549. 125358 405 0 394 8 10 515 8 899 10 875 2560 . 125771 406 0 395 6 10 258 8 675 10 601 2570 . 126147 407 0 396 3 10 019 8 466 10 346 2580. 126491 408 0 397 1 9 796 8 272 10 109 2590. 126808 409 0 397 8 9 594 8 097 9 886 2599. 127232 410 0 398 6 9 397 7 926 9 678 2609. 127498 411 0 399 4 9 213 7 765 9 482 2618 . 127740 Sample output file <mass> Tw: wall temperature (C) Tb: bulk temperature ( C ) h: heat transfer coefficient using Deissler-Van Driest model (kW/m2C) M: mass transfer coefficient using Deissler-Van Driest model (mm/s) Ms: mass transfer coefficient using Swenson correlation for mass transfer (mm/s) Hb: bulk enthalpy (kJ/kg) Re: Reynolds number Tw(C) Tb (C) h M Ms Hb Re 381 0 375. 5 19. 433 1 780 1 780 1873. 64679 382 0 377 . 5 23. 795 1 936 2 071 1911. 67917 383 0 379 6 31. 127 2 164 2 434 1969. 73491 384 0 381 1 37 . 428 2 475 2 789 2050. 82092 385 0 382 0 35. 518 2 805 3 035 2127 . 90871 386 0 382 5 30. 980 3 074 3 223 2186. 97520 387 0 383 1 27 . 141 3 286 3 386 2232. 102299 388 0 383 6 24 . 203 3 460 3 481 2269. 105981 389 0 384 1 21. 931 3 609 3 . 605 2300. 108924 390 0 384 7 20 . 130 3 739 3 .716 2328 . 111349 391 0 385 3 18 . 688 3 854 3 .817 2352 . 113401 392 0 385 9 . 17 . 476 3 958 3 . 908 2374 . 115140 393 0 386 5 16. 464 4 056 3 . 993 2394 . 116755 394 0 387 1 15. 587 4 142 4 .071 2413. 118032 395 0 387 8 14 . 827 4 221 4 . 119 2430 . 119163 396 0 388 4 14 . 161 4 295 ' 4 . 186 2446. 120171 397 0 389 1 13. 573 4 365 4 .250 2461. 121068 398 0 389 8 13. 049 4 430 4 .310 2475. 121876 399 0 390 5 12 . 579 4 .492 4 .367 2489. 122606 400 0 391 2 12 . 153 4 . 550 4 .421 2502 . 123267 90 401 0 391 9 11 767 4 606 4 473 2514 . 123867 402 0 392 6 11 414 4 658 4 523 2526. 124409 403 0 393 4 11 090 4 709 4 571 2538 . 124903 404 0 394 1 10 791 4 757 4 616 2549. 125358 405 0 394 8 10 515 4 803 4 660 2560. 125771 406 0 395 6 10 258 4 848 4 703 2570 . 126147 407 0 396 3 10 019 4 891 4 744 2580 . 126491 408 '0 397 1 9 796 4 932 4 783 2590 . 126808 409 0 397 8 9 594 4 977 4 807 2599. 127232 410 0 398 6 9 397 5 015 4 844 2609. 127498 411 0 399 4 9 213 5 052 4 880 2618 . 127740 91 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

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"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0080975/manifest

Comment

Related Items