GASP - A Computer Code for Calculating Thermodynamic and Transport Properties for Ten Fluids, Lecture notes of Thermodynamics

Download GASP - A Computer Code for Calculating Thermodynamic and Transport Properties for Ten Fluids and more Lecture notes Thermodynamics in PDF only on Docsity! 1. Report No. 2. Government Accession No. NASA TN D-7808 4 Title and Subtitle GASP -A COMPUTER CODE FOR CALCULATING THE THERMODYNAMIC AND TRANSPORT PROPERTIES FOR TEN FLUIDS- PARAHYDROGEN, HELIUM, NEON, METHANE, NITROGEN, CARBON MONOXIDE, OXYGEN, FLUORINE, ARGON, AND CARBON DIOXIDE 7. Author(s) Robert C. Hendricks, Anne K. Baron, and Ildiko C. Peller 9. Performing Organization Name and Address Lewis Research Center National Aeronautics and Space Administration Cleveland, Ohio 44135 12. Sponsoring Agency Name and Address National Aeronautics and Space Administration Washington, D.C. 20546 3. Recipient's Catalog No. 5. Report Date Feb "'US, z":,' 1 : )r,"[2 6. Performing Organization Code 8. Performing Organization Report No. E-6501 10. Work Unit No. 502-24 11. Contract or Grant No. 13. Type of Report and Period Covered Technical Note 14. Sponsoring Agency Code 15. Supplementary Notes 16. Abstract A FORTRAN IV subprogram called GASP calculates the thermodynamic and transport prop- erties for 10 pure fluids: parahydrogen, helium, neon, methane, nitrogen, carbon monoxide, oxygen, fluorine, argon, and carbon dioxide. The pressure range is generally from 0. 1 to 400 atmospheres (to 100 arm for helium and to 1000 atm for hydrogen). The temperature ranges are from the triple point to 300 K for neon; to 500 K for carbon monoxide, oxygen, and fluorine; to 600 K for methane and nitrogen; to 1000 K for argon and carbon dioxide; to 2000 K for hydrogen; and from 6 to 500 K for helium. GASP accepts any two of pressure, temperature, and density as input conditions. In addition, pressure and either entropy or enthalpy are also accepted as input variables, a desirable flexibility for cycle analysis. The properties available in any combination as output include temperature, density, pressure, entropy, enthalpy, specific heats (Cp, Cv) , sonic velocity, (0P/_p)T, (0P/_T)p, viscosity, thermal conductivity, and surface tension. The subprogram design is modular so that the user can choose only those subroutines necessary to the calculations. GASP was written for the engineer user who requires both accuracy and speed in calculating thermodynamic and transport properties. 17. Key Words (Suggested by Author(s}) Fluids; Thermodynamic; Transport prop- erties; Computer; Hydrogen; Helium; Neon; Methane; Nitrogen; Carbon monoxide; Oxygen; Fluorine; Argon; Carbon dioxide 19, Security Classif. (of this report) Unc lassifi ed 18. Distribution Statement Unclassified - unlimited STAR category 34 (rev.) 20 Security Classif. {of this page) 21. No, of Pages Unclassified 210 ' For sale by lhe National Technical Information Service. Springfield, Virginia 22151 22. Price* $7...5 GASP - A COMPUTER CODEFOR CALCULATING THE THERMODYNAMIC AND TRANSPORT PROPERTIESFOR TEN FLUIDS: PARAHYDROGEN, HELIUM, NEON, METHANE, NITROGEN, CARBON MONOXIDE, OXYGEN, FLUORINE, ARGON, AND CARBON DIOXIDE by Robert C. Hendricks, Anne K. Baron, and Ildiko C. Pellet Lewis Research Center SUMMARY A FORTRAN IV subprogram, GASP, has been developed to calculate the thermo- dynamic and transport properties for 10 pure fluids: parahydrogen, helium, neon, methane, nitrogen, carbon monoxide, oxygen, fluorine, argon, and carbon dioxide. The pressure range is generally from 0.1 to 400 atmospheres (to 100 arm for helium and to 1000 arm for hydrogen). The temperature ranges are from the triple point to 300 K for neon; to 500 K for carbon monoxide, oxygen, and fluorine; to 600 K for methane and ni- trogen; to 1000 K for argon and carbon dioxide; to 2000 K for hydrogen; and from 6 to 500 K for helium. Two-phase liquid-vapor properties are available; however, fluid mix- tures are not at this time part of GASP. GASP accepts any two of pressure, temperature, and density as input conditions. In addition, pressure and either entropy or enthalpy are also accepted as input variables, a desirable flexibility for cycle analysis. The properties available in any combination as output include temperature, density, pressure, entropy, enthalpy, specific heats (Cp, Cv) , sonic velocity, (_P/_P)T' (_P/_T)p, viscosity, thermal conductivity, and surface tension. The subprogram design is modular so that the user can choose only those subrou- tines necessary to his calculations. Existing equations are used for all fluids except fluorine and hydrogen. The equations for these fluids were developed by the authors from National Bureau of Standards data and the PVT surface curve-fitting technique of Bender. GASP was written for the engineer user who requires both accuracy and speed in calculating thermodynamic and transport properties. INT RODUC T ION Cryogens are important fluids in many areas of modern technology. Among these are food preparation, power transmission, antipollution power-conditioning systems, neurosurgery, the H2-O 2 propulsion systems for land-based and space vehicles, and proposed high-energy fuels for aeronautics. Because their properties are known and gradients are readily established, cryogens are often the primary test fluids in heat- transfer and fluid dynamics research. A Lewis Research Center literature and computer library search revealed no single computer code that could provide a flexible medium for calculating internally consistent thermodynamic and transport properties of the desired cryogenic fluids. A subprogram (subroutine) exists for nitrogen, called NTWO (ref. 1), which has the desired input-output flexibility. However, its results for nitrogen are not the best currently available. We decided to construct a multifluid, multiproperty computer subprogram based on the NTWO input-output structure using, for eight fluids, the best available equations that were also compatible with the desired program logic structure. The program is called GASP, an acronym for gas properties (ref. 2). We then decided that the only feasible way to maintain a viable properties program is by modifying and updating a basic properties program at periodic intervals. For this report, the 10 fluids chosen are based on the current and projected requirements of cryogenics research and development for aerospace and environmental applications at the Lewis Research Center and other NASA facilities and in industry. The equation-of-state coefficients are programmed so that fluids may easily be added, replaced, or updated when new equations or additional data become available. As such, we suggest that this report represent an ongoing program of maintaining GASP, where future revisions can be made by issuing looseleaf inserts to replace pages of this document. The main report consists of two sections which are directed to the research-oriented GASP user who must know both what is included in GASP and how to use it. The first section, ANALYSIS OF GASP, explains the contents of the program. It includes the Bender method for pressure-volume-temperature (PVT) surface curve-fitting (refs. 3 to 5); the equations of state selected for eight of the fluids; the equations of state devel- oped by the authors for parahydrogen and fluorine; a discussion of the transport property calculations, including the prediction of the thermal conductivity anomaly in the critical region; comparisons to measured PVT data for each fluid; the question of extrapolation to higher pressures for parahydrogen in particular and other fluids in general; and com- parison of PVT data to other equations of state. The second section, USER'S GUIDE TO GASP, presents detailed instructions for in- put and output summarized in table form for handy reference. Three accuracy ranges for thermodynamic calculations using PVT measured data are discussed in relation to GASP. The use of a modular design in the programming of GASP simplifies dividing it 2 into smaller subprograms for users with limited interest or limited core storage. Cer- tain options are also provided for users wishing to make calculations in troublesome areas of the PVT surface. The mathematical and FORTRANsymbols used are defined in appendixA. The exact equationsusedare presented in appendixesB to E. Detailed information about the com- puter program is given in appendixesF to J. And the properties of parahydrogen are discussed in appendixK. ANALYSIS OF GASP THEORY (EMPIRICAL EQUATIONS) Selection of Equation of State We relied heavily on the published and private work of the German researcher Dr. Eberhard Bender (refs. 3 to 5). The equations of state for all but carbon monoxide, neon, and helium are those published by Bender or developed by us using Bender's com- puter program. The history of equation-of-state development is well delineated in the literature and only the important aspects of the authors' work will be mentioned. The modified Benedict-Webb-Rubin (BWR) virial equation (ref. 6) developed by Strobridge (ref. 7) was an important advance. The form of this equation is 4 2 P : P(T,p) : _ Ai(T)pi + _ Bj(T)p2j+l e-CP 2 i=l j=l (1) where P is the pressure and Ai(T) and Bj(T) are polynomials in T and T -1. (All symbols are defined in appendix A and details are found in appendix B. ) This equation enabled PVT calculations to be made in the liquid and vapor phases with a single set of coefficients. The derived properties of enthalpy, entropy, (_P/aT)p, and (_P/_P)T could also be obtained. The equation, however, had several serious short- comings: PVT predictions at the higher densities were greatly in error, calculations could not be made at the liquid-vapor boundary without an additional saturation equation, calculations could not be made in the two-phase or metastable region, and the _2p/aT2 did not give good values so that specific heats could not be calculated directly. Bender (ref. 4), whose major interest was fluid mixtures, used a modified Strobridge-BWR equation with additional terms in the "i" summation to improve PVT predictions at the higher densities. This equation is Equations for fluorine and parahydrogen. - Prydz and Straty (ref. 18) have recently measured PVT data for fluorine. Their report on the properties of fluorine and on the entire F 2 PVT data set used by Prydz (including his own and other researchers} was chosen as the best and most up-to-date source for fluorine. The coefficients which we obtained for fluorine, 2 table II, by using Bender's curve-fitting technique compare well with the data as reported in the next section. We found no single coefficient set for the fluid parahydrogen, although its properties are and have been of much interest. Many limited-range equations, especially the early work of Roder and Goodwin (refs. 8 and 19) have been available in computerized forms. Tables of smoothed hydrogen data and their extrapolations have been compiled and com- puterized, the most popular being the NBS TABCODE computer program. 3 A most re- cent publication by McCarty and Weber (ref. 20) uses several equations of state, depend- ing on the region of interest, and includes many useful parahydrogen properties over a wide PVT range. We wanted a coefficient set for a Bender parahydrogen equation to use in GASP that could also compare favorably to reference 19. Goodwin's parahydrogen data set and Bender's curve-fit program were used to develop a set of coefficients (table II) that pre- dict parahydrogen PVT data quite well and can also reproduce the NBS high-pressure PVT extrapolations of reference 20. The simplicity of this Bender-type equation (eq. (2)) compared to the computerization of reference 19 should be preferable for all but the most critical user. Calculation of Transport Properties Transport properties are related to the dynamics of the system (e. g., viscosity to fluid flow, and thermal conductivity to heat transfer}. Precise measurements in such systems are extremely difficult, and errors of over 5 percent in the gaseous phase and 15 to 20 percent in the dense-gas or liquid phase are not uncommon. Consequently, at this point in the development of GASP, we chose to abandon "precision" in favor of con- sistency. In order to do this, we chose generalized forms for fitting the transport prop- erties which, for the most part, are (1) functions of only temperature for the dilute gases and (2) functions of only density (excess data, i.e., _ - p* and k -k*) for dense fluid regions. Such representations can obviously lead to erroneous results. For example, thermal conductivity data near the critical density but at 1.5 times the critical temper- ature are not representative of all data along the critical isochore, as would be implied; and viscosity coefficients are temperature sensitive at reduced densities greater than 2. 2The derived properties of fluorine need further work. The authors recommend that the user check the results above 1 MPa (10 atm). The extensive PVT data were made available from NBS. 3Available from H. M. Roder, National Bureau of Standards, Boulder, Colorado. Viscosity. - The dilute-gas viscosity for all fluids in GASP, with the exception of hydrogen, is based on the simplified Chapman-Enskog model iref. 21) t_ = 0. 2669x10 -4 2gt*(2, 2) (4) The dilute-gas viscosity for hydrogen is based on the work of Diller, Hanley, and Roder (ref. 22). Extensions of this work to higher temperatures (1000 K) proved satisfactory. However, above 1000 K, the Chapman-Enskog model was again used. The Enskog dense-gas theory suggests that an excess in viscosity can be treated as a function of density alone. From the work of Jossi, Stiel, and Thodos (ref. 23), this seems justified for the fluids in GASP and was modified and adopted for use herein. 4 [_(_ _*)+10-4] 1/4 E a i_ = iPR i=0 (5) where T1/6 4- c MI/2 2/3 Pc The formulation appears valid except for dissociation and some regions beyond PR = 2 where some temperature dependence can be noted. For fluorine, a small correction factor was added, based on the work of Hanley 4. The regions covered and the specific formulation and coefficientsused are given in table X and appendix B. Thermal conductivity. - E ucken proposed that the dilute-gas thermal conductivity be the sum of translational and internal molecular energies times some factor # The Mason and Monchick model (ref. 24) gives the sum as _-'_fiCi=1.32Cv+ 1.77R-0.91qZ_ 4private communication with H. J. M. Hanley of National Bureau of Standards, Boulder, Colorado. (7) 7 where Z represents the collision number. Brokaw (ref. 25) determined a simple em- pirical form for the sum _-_ fiCi = 1.24 Cv + 2.54 While the preceding equations are the basis of dilute-gas theory and used to determine the thermal conductivity of fluorine, empirical formulations given by Stiel and Thodos (ref. 26) and by Roder and Diller (ref. 27) for hydrogen are used in GASP. gives a reduced form (8) Reference 26 (9) while reference 27 gives a slightly more elaborate form for the sum. See appendix B and table X for specific forms and coefficients and table III for the sources of the con- stants and coefficients. For the dense-gas region, excluding near critical, the excess thermal conductivity (like viscosity) is a function of density alone. The exceptions here are the thermal con- ductivities of helium and hydrogen near the saturation boundaries, which are also func- tions of temperature. These fluids are handled separately and use the fits determined by McCarty and Roder (refs. 16 and 28) and by Roder and Diller (ref. 27) for helium and hydrogen, respectively. For other fluids covered by GASP, the formulation of Stiel and Thodos (ref. 26) is used: lOgl0 [(k - k*)Zc5 _1 = Z ai(l°gPR )i (I0) i=O where P CZ - c RPcT c and 1/6 /-MTc X- p2/3 C other comparisons, those in the near-critical range are the most revealing since this portion of the PVT surface is the most difficult to represent by an equationof state valid for the whole surface. The comparisons for the fluids fit by the authors, fluorine and parahydrogen, and for a fluid of interest to the Lewis Research Center's Aerospace SafetyResearch and Data Institute, oxygen, are madewith most of the available data for the fluid as furnished to the authors by the NBSCryogenic Data Center. Comparisons of data to values calculated from a state equationfor carbon monoxideand neoncanbe found in references 9 and 10, respectively. The results of these comparisons are presentedas figures 1to 18and summarized in table IV. The measuredPVT data are the standard for comparison, and three plots are shownfor each set with the percent deviations plotted as the ordinate andthe reduced comparedparameter, T, P, or p in that order, as the abscissa for each fluid. Methane. - Figure 1 compares GASP-calculated PVT data for methane to the data of Jansoone, Gielen, DeBoelpaep, and Verbeeke (ref. 34). The temperatures are within +0.1 percent, and the pressures generally agree within -0.5 to 1.0 percent. The den- sities show some irregularities near PR = 0.98; but since the temperature range is within +2 K of the critical value (190.77 K), this is expected. In general, the GASP pro- gram will show significant deviations slightly below the critical point, but it is good above the critical point. This may be caused by the mismatch, near the critical point, between Bender's saturation data and the vapor pressure equation used in GASP (eq. (B18)). While reference 34 presents an equation that predicts these data very well, the equation is valid only for the near-critical range. The PVT data of Vennix (ref. 35) for methane are compared in figure 2. The GASP- calculated temperatures are within +0.25 percent o_ta, except for a few points which are within +1.0 percent; the pressures are generally within +1.0 percent, with more than one-half of the pressure values agreeing within +0.5 percent. A few pressures near the maximum density range (~2.5 pc } show deviations greater than 2 percent. The GASP densities are within _2.5 percent, with most of the values within +1.0 percent. One density value at 190.80 K showed a 6.5 percent deviation from the Vennix data. Vennix (ref. 35) presented an equation of state for methane with 24 co- efficients which predicts his PVT data with about the same deviation in pressure as the Bender equation in GASP (fig. 3}. The claimed density predictions of reference 35 are better than the GASP values. We tested the Vennix equation and found it unacceptable for PVT and derived property calculations at temperatures less than ~155 K and for some higher ranges in derived properties. The PVT data of Goodwin (ref. 13} for methane are compared in figure 4. The GASP-calculated temperatures are within -0.8 to 0.3 percent, with most data within _0.2 percent and major departures confined to 0.7 _ T R _ 1. The calculated pressures are most satisfactory for PR > 1, with deviations of -0.5 to 3.0 percent, and least 11 satisfactory at very low pressures, where there are five points with deviations above 20percent. The cal_'_[ateddensities are within +0.5 percent, except in the critical region (PR _ 1), where departures of -3 to 2 percent are noted. Nitrogen. - The Strobridge PVT values for nitrogen were calculated by using equa- tion (1) and necessary constants from reference 7. The equation was extrapolated to 1000 K and 50 megapascals (or MN/m2). As might be expected, these extrapolated PVT data show the largest deviations from Bender's results. The results in figures 6 to 8 show temperatures within 22 percent for densities less than 2.5 Pc" The temperature deviation in the dense liquid ranges to -4 percent, with the largest errors at the highest densities. The pressures are within +5 percent, except for the dense-liquid low- temperature region, where deviations to -10 percent and greater are found. The Strobridge fit is inadequate in this portion of the PVT surface, and we consider Bender's values more accurate than the Strobridge-BWR fit of the standard PVT surface. The densities, figure 6(c), are within 22.0 percent, except for deviations near the critical isochore. Figure 9 compares calculated PVT values with the liquid-saturation and near- critical data of Weber (ref. 36). The density comparison is most crucial in the near- critical area of the PVT surface. The values calculated by GASP are within +0.5 per- cent, except for data on two isochores at 0.265 and 0. 467 g/cm 3, where deviations are from -3.0 to 1.0 percent. The liquid values (high density) predicted by Bender are lower than those of Weber by approximately 0.003 g/cm 3, while the Bender density along the 0.265-g/cm 3 isochore is slightly higher. Figure 10 is discussed later. Coleman and Stewart (ref. 37), after an exhaustive compilation of the literature, derived an equation of state similar to that of Strobridge (ref. 7). A comparison of this wide-range PVT equation presented in reference 38 with the Bender fit indicated little difference in PVT results. However, the equation was judged unacceptable because of anomalous calculated derived properties and also because it was the more complicated form. Oxygen. - From the PVT data of Weber (ref. 39), 1352 data points were used by Bender (ref. 4) in his oxygen fit. We compare 1478 PVT data points from reference 39 with those computed from GASP in figure 11. Most temperatures are within 20.15 per- cent, including the PVT saturation data. The majority of the pressures are within 1 percent, except in the high-density region (p > 2.5 pc ), where 1 to 2 percent devia- tions are common and a few deviations to 20 percent are found at very low temperatures (_70 K). Most densities are within +0.2 percent, except in the critical region, where the maximum deviation is 22 percent. Figure 12 is discussed later. Argon. - Subcritical PVT data for argon reported by Van Itterbeek, Verbeeke, and Staes (ref. 40) are compared with GASP-computed values in figure 13. The temper- atures agree within _:0.3 percent, with most being within +0.1 percent. The pressures are within +5 percent, with a few errors to 10 percent for the 90 to 110 K range. The 12 densities are in excellent agreement, with more than one-haLf the Bender values being within +0. I percent of the data. We noted two discrepancies in the plotted and tabulated data of reference 40, so the temperature value selected for two of the isotherms may be in error by 0.05 K. Carbon dioxide. - The PVT data of Michels, Blaisse, and Michels (ref. 41) are com- pared to the GASP-calculated carbon dioxide values in figure 14. A few values listed in reference 41 were in the GASP two-phase region and are not considered. The temper- atures are within +0. i percent. The pressures are within +0.5 percent, except for three deviations to 5 percent at the highest measured density. The density deviations are well scattered within 14 percent. Helium. - Helium values for five isobars and a temperature range of 2.5 to 1500 K were calculated by using McCarty's three-region fit (ref. 17) as the best available source for comparison. The comparison to the Mann equation used in GASP is shown in fig- ure 15 (eqs. (BI) and (B2)). The temperatures are within +i percent, except for near- critical deviations to 3 percent and dense-liquid (P > 2.8 pc ) deviations to +i0 percent. The pressures are within +0. I percent for T _ T c and within +10 percent at subcritical temperatures. Large density deviations occur in the near-critical region (3 to 6 K), and elsewhere density agrees within 1 percent. The PVT range for helium based on Mann's equation from reference 11 (table i) is quite limited compared to reference 17. There- fore, part of the comparison is really an extrapolation and hence subject to uncertainty. We do not recommend the use of Mann's equations in GASP for T < 6 K. Fluorine. - All extensive compilation of the thermodynamic properties of fluorine has been made by Prydz and Straty (ref. 18), including 850 points of measured PVT data. These data and saturation data from reference 42 were used by the authors to obtain the fluorine coefficients in table If. The comparisons are shown in figure 16. The temper- atures are within +0.15 percent, except for very low temperatures, where a few devia- tions are near 1 percent. The pressure errors are within 1 percent, except for the re- gion where densities are greater than 2.5 Pc (melting locus) and temperature is less than 90 K, where larger errors occur. The densities are within +0.5 percent, except near the critical isochore or isotherm, where deviations to 3 percent occur. Parahydrogen. - The parahydrogen data compiled by Roder, Weber, and Goodwin (ref. 19) and by McCarty and Weber (ref. 20) are compared to GASP-calculated values in figures 17 and 18. The temperatures agree within +0.25 percent, except for a few points which deviate to +0.5 percent. The pressures are within +1 percent, except for several points in the near-critical area or near the melting locus. The densities are within +0.5 percent, ex- cept near the critical isochore, where deviations to +3 percent occur. Densities greater than 0.06 g/cm 3 are all within +0.1 percent. The critical density Pc of reference 19 is 0. 03142 g/cm 3, whereas the Pc we determined in the Bender fit is 0.0311 g/cm 3. 13 density and calculation of pressure gave meaningful results. Nitrogen - generated "data. " - While a consistent set of nitrogen data is not avail- able, several comparisons can be made in figure 7. The PVT values were compared by inputting either (P,T), (P,p), or (p, T) as generated by GASP into the equation of Cole- man and Stewart (ref. 37). The deviations are quite small, but a change in form of equation (BI) does not appear warranted. In figure 8 a similar comparison was made to the 33-term modified BWR equation of Stewart, Jacobsen, and Meyers (ref. 44). Generated P,T data gave a comparison of deviations in density, as shown in figure 8(b). The deviations are well within 1 percent, with increased deviations at higher densities. However, the generated density, when in- putted into the 33-term equation, did not always properly converge, and significant devia- tions for T < T c were noted. We then generated p,T and compared the deviations in calculated pressure. Again the deviations for T < T c or T -< Tsat in the liquid region were large, and the conver- gence not always assured. The technique of solution and the cause of these irregularities require further investigation. Nitrogen - measured data. - Figure I0 compares the calculated values of refer- ence 44 with the near-critical nitrogen data of Weber (ref. 36). The deviations in pres- sure (-i to 3 percent, fig. 10(a)) appear slightly better than those of figure 9(a). The deviations in density range from I. 5 percent at low densities to +0. 1 percent at higher densities (fig. 10(b)) and also appear slightly better than those of figure 9(b). Methane. - The recent methane data of Goodwin (ref. 13) represent an accurate con- sistent set of PVT measurements for fluid methane. Using these data and the techniques of Hust and McCarty (ref. 45), McCarty (ref. 46) developed a Stewart-type equation, and the relative errors in pressure are given as figure 5(a). For the most part, the relative error falls between +I percent, with the exception of a few points at lower pressure (with rather large deviations, _ i0 percent) and some near the critical pressure. The density deviations of figure 5(b) are for the most part within -0.25 to 0.5 percent, except for a few points near the critical density. Comparing figure 4(b) with 5(a) and figure 4(c) with 5(b) gives Parameter Pressure Density Region Low Near critical High Near critical Other High only Equation of state (eq. (B1)) N=6 M=2 N=9, M =6 (Bender-type equation, ref. 4) Stewart-type equation, ref. 44) A few points (5) to 40 perceut A few points to 10 percent; for the most part, lto5percent -1 to 3 Percent -2 to 1.5 Percent -0.3 to 0.1 Percent -0.4 to 0.05 Percent A few points (4) to 10 percent ±2 Percent ±1 Percent -0.7 to 1 Percent ±0.2 Percent -0.25 to 0.6 Percent -0.2 to 0.05 Percent 16 The McCarty fit (ref. 46)of the Goodwindata does represent the data more accurately than the Bender-type fit. However, a refit of methaneusing the Goodwindata may also prove fruitful. Oxygen. - The data of Weber (ref. 39) represent an accurate and consistent set of PVT data for fluid oxygen. These data and the Bender-type equation are compared in figure 11. When the 33-term equation and coefficients as determined by Stewart, Jacobsen, and Meyers (ref. 44) are used, the relative error in pressure appears as in figure 12(a). The relative error at low pressure appears systematically scattered about +2 percent, with a few points outside this range. At higher pressure the scatter appears about +1 percent, again a few points outside this range. The relative error in density is illustrated in figure 12(b). For the most part the scatter is within +0. 1 percent, except near the critical density and some points at rather low densities. Comparing figure ll(b) with 12(a) and figure ll(c) with 12(b) reveals Parameter Region Equationof state(eq. (BI)) N=6, M=2 (Bender-type equation, ref. 4) N=9, M=6 (Stewart-type equation, ref. 44) Pressure Low +2.5 Percent +2 Percent !High +1.5 Percent :el Percent Density Near critical +2 Percent +0.8 Percent Dther +0.3 Percent +0.1 Percent From these comparisons, either equation can be used with good results, except in pre- dicting density near the critical point, where the Stewart-type equation is better. Parahydrogen. - The parahydrogen data of Goodwin (refs. 8 and 19) represent prob- ably one of the most exhaustive PVT studies in the literature. We used these data to de- termine a Bender-type equation (eq. (B1)) for parahydrogen. The PVT comparisons are made in figure 17. When the 33-term equation and coefficients as detel mined by Stewart (ref. 44) are used, the relative error in pressure is given as in figure 18(a). The rela- tive errors lie for the most part between +1 percent at low pressures and +0. 5 percent at high pressures. The relative errors in density (fig. 18(b)) are for the most part within +0.2 percent, except near the critical density, where errors to s0.7 percent are noted. Cross-comparing the relative errors of figure 17(b) with 18(a) and figure 17(c) with 18(b) reveals 17 Parameter Region Equation of state (eq. (BI)) N=6, M=2 N=9, M=6 (Bender-type equation, ref. 4) (Stewart-type equation, ref. 44) Pressure Low -0.5 to 1.5 Percent -1.0 to 0.8 Percent Moderate :e0.2 Percent _-0.4 Percent High +0.3 Percent ±0.5 Percent Density Near critical +2 Percent +0.6 Percent Other -0.3 to 0.2 Percent ±0.2 Percent It would appear that the Bender-type equation is superior except in predicting density in the near-critical region, where clearly the number of constraints of reference 44 (13 additional terms) does appear helpful. The greatest advantage of having a single accurate and consistent set of data to de- rive a curve fit is the accuracy with which it reproduces the data. This is readily ob- vious from the parahydrogen results, where both McCarty (ref. 46) and the authors used Goodwin's data to develop a PVT surface. The results are self-evident. There is virtually no improvement in predicting the PVT surface, except for density near the critical point. Here the extra terms do seem to merit use. While similar remarks can be made for oxygen, methane serves as an example of not having a self-consistent set of data available (e. g., the comparison of the Bender-fit to the McCarty fit in the section Methane). Clearly, the McCarty fit is better. The question as to whether it would be better if all the methane data were used in the fit still remains to be investigated. Some current drawbacks in changing to a 33-term equation (Stewart type) are as follows: (1) The coefficients are available for only four fluids. Thus, several refits would be required. (2) Such an equation is limited to specific computing machines (a) Because it requires double precision, 10 +88 capability, additional storage, and data set modifications 00) Because exponentials can become sufficiently large as to cause overflow and underflow (3) Thirteen additional constants (i. e., 33 as opposed to 20 for the Bender-type equation) are required. (4) The first virial coefficient is not a polynomial. Some current advantages in changing to a 33-term equation are (1) Better PVT representation of the near-critical region and along the freezing locus (2) Better behavior of _2p/_T2 at low temperatures (One may not have to use num- erical determination of C v in this range, as is now done in GASP. ) 18 Since many coefficients are put into the program COMMONSby SETUP, it is recom- mendedthat all calculations for onefluid be madebefore changingfluids. (4) The controls KR, KS, andKP which tell GASPwhatvariables are to be used as input andwhat properties are requested for output must be correctly initialized in the call statement for subroutine GASP. The corresponding input variables in the call state- ment and COMMON/PROPTY/... must also be correctly initialized. The controls KS andKR determine which of the variables T, P, D, H, or Sor com- binations thereof are neededas thermodynamic input. The input control KP specifies which properties are soughtas output. KR is also an output variable since it gives the correct region number for the variables in a specific call to GASP, as shownby the sketch in table V. Dependingon the input for KS andKP, the other possible output var- iables are T, P, D, H, and all of COMMON/PROPTY/ except the control KU. As mentionedabove, KR is both input and output andmust be reset before each call to GASP. The input options are (1) KR=0 when the user wishes GASP to determine a value for KR (2) KR=I when the user wishes saturation conditions 7 The output for KR will be (1) KR=I for saturation (2) KR=2 for liquid (3) KR=3 for vapor The control KS specifies which variables are to be used as input for a call to sub- routine GASP. (In the remaining discussion on GASP input and output, the input var- iables are assumed to be in user's units specified by KU. Output is always returned in the KU system of units. ) The following table shows the input and output for all KS, KR combinations: 7Saturation or coexistence conditions exist on the PVT surface when pressure is a function of only temperature and the liquid and vapor states both exist at that pressure. Thus, when KR=I, two outputs for each property are available in COMMON/PROPTY/and only one independent vari- able is required for some input options, as shown in the KS--KR input/output chart (table V). To find saturation properties given temperature TI: prior to the call to GASP, set P-=-0, T=T1, and KR=I. Select the proper KP and CALL GASP(KS, KP, T, P, D, H, KR). The saturation pressure is returned as P. All other properties are in COMMON/PROPTY/etc. To find saturation properties given pressure Pl: follow the above procedure except set T=0 and P=Pl. The saturation temperature is returned as T. If you give GASP both T and P, it uses T and does not alter P. 21 KR KS Input 0 TandP TandD PandD PandH PandS 1 T or pa T P P P Outpat 1 T or pa, p T T,DL, and DV T,DL, and DV DL and DV 2 D P T Dand r DandT 3 D P T D and T D and T aSee footnote 7, p. 21. The input control KP specifies which derived and transport properties are requested by the user. It is the sum of the individual KP options and is described in table V. This binary sum allows GASP to uniquely identify any combination of requests. The following table shows the output locations for the specific KR, KP combinations: Value added to KP input 8 16 32 Output for KR=2 or KR=3 H S CP CV GAMMA C MU K EXCESK SIGMA Output for KR= 1 Liquid Vapor HL HV SL SV CPL CPV CVL CVV GAMMAL GAMMAV CL CVP MUL MUV KL KV EXCL EXCV SIGMA ........ Name of calculated property None requested Enthalpy Entropy Specific heat at constant pressure Specific heat at constant volume Specific-heat ratio Sonic velocity Viscosity Thermal conductivity Anomalous thermal conductivity Surface tension of fluid as func- tion of temperature TROUBLESHOOTING FOR USER ERRORS After experience with GASP, we have found that several common errors are easily detected and corrected: (1) Failure to set 1 -< KU -< 5 will cause a "division by 0.0" and/or no valid an- swers. Set KU to its proper value. (2) Failure to set 1 -< KS -< 5 will most likely cause a halt to the program because 22 of an execution error. The branching on KS in subroutine GASP is a computed "GO TO. " Simply set KS to its proper value. (3)Failure to set KP will return enthalpy ifKP is odd and no derived properties if KP is even. (4)Ifa wrong value is entered for KR, itis treated as KR=0. Ifa user enters KR=I when he does not want saturation properties, he will get them anyway for T < T c and otherwise will get a wrong answer. (5)Ifany T, P, D, H, or S is entered incorrectly, that value will be used and the answer willbe wrong. (6)Ifthe COMMON/PROPTY/is duplicated incorrectly, there are a variety of pos- sible errors, almost all serious. (7)Ifthe call to SETUP has an illegalfluidrequest (notone of i0 fluids,keypunch error, etc.), the execution stops after an error message. Other small problems may be encountered ifGASP is modified for different com- pilers or computers. The FORTRAN IV coding in GASP is machine independent except for a few Hollerith format statements which can be easily changed. The reader who needs more detailed information should read the appendixes. ADDITIONAL INFORMATION The approximate core storage for the complete GASP program is (23306)8= (9926)10 locations. The following problems have previously been encountered when converting to non- IBM machines or to different FORTRAN IV - FORTRAN V compilers: (1) Users of the IBM 360-67 should run in double precision by inserting IMPLICIT REAL*8 (A-H, O-Z) and REAL*8 MU, MUL, MUV, K, KL, KV in subprogram GASP and IMPLICIT REAL*8 (A-H, O-Z) in all other subroutines. COMMON/PROPTY/KU, KZ, DL, DV, etc., should be changed for proper alinement. (2) Data statements are found in subroutines BLOCK DATA, THERM, VISC, and SURF. Many compilers differ in formatting data statements. (3) The multiple-entry routine CHECK (TCHECK, PCHECK, DCHECK) has an entry point, DCHECK, whose input vector (KU, D) does not correspond in kind and number with the other entry points (KU, KR, T) or (KU, KR, P). To our knowledge, this has caused a problem on only one compiler, a FORTRAN IV for a CDC 3800. It was easily remedied by an equivalence statement. We adapted the code to fit all the compilers and machines 8 tried so far with the pre- ceding exceptions. 8The machines tried so far are the UNIVAC 1108, CDC 3600, CDC 3800, IBM 360/67, and IBM 7094-7044 DC. 23 Mathematical symbol k k* M ml_5 n1-24 P Pc Pmax Pmin P sat (ap/ )p (_P/ap)T FORTRAN sy mb ol K KCP KL KP KR KS KU KV NAMGAS P PC PS PSS PTV PTVL PTW PDT Description thermal conductivity, W/(cm)(K) thermal conductivity of dilute gas, W/(cm)(K) region delimiter used in numerical calculation of Cv thermal conductivity of saturated liquid, W/(cm)(K) thermodynamic and transport properties specification thermodynamic region specification state relation specification units specification thermal conductivity of saturated vapor, W/(cm)(K) molecular weight, g/g-mole set of constants for equation of specific heat at "zero" pressure Hollerith code used to specify fluid set of constants for equation of state pressure, MPa pressure at thermodynamic critical point upper pressure limit, MPa lower pressure limit, MPa pressure used internal to program, MPa saturation pressure used internal to program, MPa saturation pressure (eq. (3)), MPa partial derivative of pressure at constant volume, MPa/K partial derivative of pressure at constant volume of saturated liquid, MPa/K partial derivative of pressure at constant volume of saturated vapor, MPa/K partial derivative of pressure at constant temperature, J/g 26 Mathematical symbol S So T Tc Tmax Tmin Tt TO U V V' V" Z Z C E/k A FORTRAN sy mb ol PDTL PDTV R S SL SV T TC TS TSS GAMMA GAMMAL GAMMAV EPSK De scr iption partial derivative of pressure at constant temperature of saturated liquid, J/g partial derivative of pressure at constant temperature of saturated vapor, J/g gas constant, J/(g)(K) entropy, J/(g)(K) reference entropy, J/(g)(K) entropy of saturated liquid, J/(g)(K) entropy of saturated vapor, J/(g)(K) temperature, K temperature at thermodynamic critical point, K upper temperature limit, K lower temperature limit, K temperature used internal to program, K saturation temperature computed by function TSS, K triple-point temperature, K reference temperature, K internal energy, H - (P/p), J/g specific volume, cm3/g volume of saturated liquid, eq. (3) volume of saturated vapor, eq. (3) compressibility factor, P/pRT compressibility factor at critical point, eq. (10) ratio of specific heats, Cp/C v ratio of specific heats of saturated liquid ratio of specific heats of saturated vapor potential parameter, used to reduce temperature, K thermal conductivity parameter, _fM Tlc/6///p2/3=c _M, cal/(cm)(sec)(g) 27 Mathematical symbol r P P Pc Pswitch Pt _d 1)* _2(2, 2)* Subscripts: C exp R t total FORTRAN symbol EXCESK MU MUL MUV D DL DS DSL RHOC DV SIGMA Description reacting conductivity, W/(cm)(K) dynamic viscosity, g/(cm)(sec) dynamic viscosity of dilute gas, g/(cm)(sec) dynamic viscosity of saturated liquid, gJ(cm)(sec) dynamic viscosity of saturated vapor, g/(cm)(sec) TIJ66_ _/_-_p2/3 viscosity parameter, c /_vl c , where the units are T c in K, Pc in atm, and viscosity in cP density, g/era3 density of saturated liquid,g/cm 3 density used internal to program, g/cm 3 density of saturated liquid, used internal to program, g/cm 3 density at thermodynamic criticalpoint, g/cm 3 density where calculationof C v changes from numeric to analytic triple-pointdensity, g/cm 3 density of saturated vapor, g/cm 3 surface tension, dyne/cm hard-sphere collisiondiameter, A (or i0-I0 m) Lennard-Jones collisionintegral for diffusion Lennard-Jones collisionintegral for viscosity critical point experimental reduced triple point total 28 The derivative (_P/_P)T (eq. (B3)) is straightforward. From equation (B1), Z .2cp p2j e-Cp2/T2 i=l j=l (B9) and (_2p/_T2)p, it is more convenient to re-To calculate the derivatives arrange the virial equation: 7 3 _cp2/T 2 p = _ Ci(P)W3"i + _ Dj(p)e T -(j+l) i=i j=i (B 10) where 3 C 1 = n22P nlp2 3 4 5C 2 =Rp+ +n6P +ngp +nllP n2p2 3 4 5 6C 3 = +n7P +nl0P +n12P + nl3P n3p2 3 5C 4 = +n8P +n24P C 5 =n4 p2 +n23 p3 2 C 6 = n5P 2 C 7 = n21P D l= p3(n14+ nl7P 2) D2=p3(n15+n18 p2) D3= p3(n16+ n19 p2) (Bil) 31 (B12) For fluids other than helium, a = 0 and T 2 = 1. In order to determine the specific heats, the /\_82p/_T2)o is required, although it is never used in that form. After differentiating equation (B12) to give (B5), the integral (eq. (B6)) becomes T 2 -i (3 - i)(2 - i)Ci(p) 2 P + 2 T2 T_LT2 T P +_ P J j=l P T 2 e-cp2/T2 T -j :f_[Ei(T)o_-l+Fi(T)oei-le-COe/T_'] d°i=l i(T)__ _ I e Fi(T) -- i 2 "= j=1 (i - I)' 0 - I): 1 (B13) 32 E I(T) = 6n4 12n5 20n21 2n 3 + __ + -- + T T 2 T 3 T 2 2n 8 6n23 E2(T) = 2n22T + __ + T 2 T 3 E3(T) = 0 2n24 E4(T ) - T 2 F I(T) : 12n15 20n16 6n14 + --+ T T 2 T 3 -- F2(T ) = -2n20 ot k k TT 2 U _14 2) nlS/_ 2 3) + n16 (_ 2 t) 1 /6 12n18 20n19 _+_ +T + _ + + "'_+T+T-) T 3 F3(T ) - -2n20 n617 -- 2 +2)+1118(1 3)n14(_ 2 4)1 +(n20)2_(n n15 n16__ TT 4 17 + n18 n19_ T+7) (B14) In order to facilitate evaluating equation (B7), rearrange the virial equation to the form IR i i-1 i p2j -cp2/T2]P = pT + Ai(W ) k + Bj(T) -- e = pW[_] i=2 T j =1 T (B15) After differentiation of equation (B15) with respect to T and substitution into equa- tion (B7), the integral becomes 33 Specific heat at constantpressure: Cp = Cv + T \_T/p (B23) An alternate method for Cv used at high density is C V - AT AT (B24) AH (B25) Cp - AT Sonic velocity: c= C_/_va(_p_)T (B26) With (aP/aT)p, (0P/ap)T, P, p, T, and the Bridgeman Tables (appendix E), the user can calculate any other thermodynamic relation of his choice. TRANSPORT PROPERTIES In general, the thermal conductivity, viscosity, and surface tension are calculated by the techniques outlined in the main text under Calculation of Transport Properties. Here we specify these calculations in greater detail. Viscosity For fluids other than hydrogen, the dilute-gas viscosity is given by * (B27)= 0.26693x10 -4 2_2"(2 , 2) 36 and the excess viscosity function is specified by . =[0.0093324 p4- 0.040758 p3+ 0.058533 p2 + 0. 023364 PR + 0. 10230] 4 - lx10 -4 I00 If the fluid is fluorine, a small correction factor is used (P-_*)F2 =(_- _*)[ 1-1sin(2"38081npR)]4 (B28) (B29) The viscosity then becomes p = (# - #*) + p (B30) If the fluid is hydrogen, a two-part fit is used which for the most part is from H. M. Roder of NBS, Boulder, Colorado: Let and /_H2 + 19.55/\T + 1175.9/ * = 1.799 T O. 6835 T > 100 K (B32) PH 2 where R -> -80 (B33) A : exp [5. 7694 + In p + 65p 3/2- 6x10 °4 e127" 2p] (B34) B = 10 + 7.2 \0.07/ - 17.63e (B35) The viscosity for hydrogen becomes = + Ae (B/T) × 10 -6 /_H 2 H 2 (B36) 37 At elevated temperatures, the value of * predicted from equation (B32)dominates, PH 2 and it becomes unnecessary to correct for density effects (i. e., one can assume that A_0). While it is known that the excess viscosity of cryogens is temperature dependent (refs. 47 to 55), it is not yet known how to represent this effect for several fluids by using corresponding-states principles. Gibbon and Kuebler (ref. 56) present a pseudo- corresponding-states analogy to predict the viscosities of argon, nitrogen (limited), neon, and helium. Haynes (ref. 55) presents a "virial" equation for argon viscosity. However, attempts to extend these equations to other fluids were met with limited suc- cess. Further, as will be shown later, the data of various investigators differ substan- tially for the same fluid region. These data appear to be dependent on the apparatus, the measurement technique, and the fluid region where the measurements are made (refs. 57 to 59). (Recall that thermal conductivity measurements suffer the same prob- lems; so perhaps do all transport properties. ) It is difficult to make generalizations, however, the authors have found some trends. In figure 22 we plot the viscosity data of Diller (ref. 50) as a function of reduced density. Similar plots were made for nitrogen and argon by using the data of Zhdanova (refs. 60 and 61). See also reference 47. (1) The i§otherms are nearly linear and intersect the saturation locus for TR < 1; the extended isotherms appear to "focus" at PR = 2. For T R > 1 the isotherms tend to "focus" near PR = 1. (2) The saturation locus appears to be continuous; however, at T R _ 1, some dis- continuities in the surface are noticeable. (3) Viscosity in the region PR -->2 can be described by p = 0.63(PR - 2)tanlm(TR) ] T R 1.3 --- TR-< 1 2 0.08(9+TR3 ) TR>I A/_ = (1. 616x10 p-5) (B37) However, the formulation is limited to this region; subsequent comparisons found equa- tion (B37) better, but not significantly better, than assuming that the temperature de- 38 Fluorine. - First, compute the dilute-gas thermal conductivity 1.77R+ 6 ' _ 2 _2(2,2)* 2 _(1, 1)* Cv'O _- _2(1, D*J v, 2 * 0.02 T O. 83333 (B43) Then add the excess function (k -k*) +k* F 2 (1344) Parahydrogen. - From an interpolation routine provided by Roder, In k* = [(Interpolation of tabulated values) - (Dependent on T alone)]H 2 (B45) Let B =39.6 - 2 65 _-_] J (B46) Then =0 lk" {C [C ( C3 )] 2)kH2 . H2 exp 1Bp + 2 + _ P T <250K (B47) where C 1 = 0. 9885311 C 2 = 32.0887 C 3 = -910. 141 For temperatures between 150 and 250 K, the thermal conductivity is given by k= (k -k*)Eq. (B41) x 4.184+ 3.383xi0 -5x T O.72872 for T > 250 K (B48) 41 and k=x(k)Eq" (B48)+ (1 -X)kEq" (B47) for 150-<T-< 250K (B49) where x = 0.01(T - 150 K) (B50) Above 2000 K, departures due to dissociation begin to become noticeable. Helium. - From an interpolation routine provided by Roder, compute function of CONZ(T): k* as a • CONZ(T) kHe - 1000 (BSl) Let lOgl0 S = -621. 369 p3 + 224.256 p2 _ 29.485 p + 2. 094 196 Then the excess function becomes (B52) - × 10-3(k. k*)H e l°gl0 (T) 10S (B53) and finally, the thermal conductivity is given by k =(k-k*)H e * (B54)+ kHe Thermal Conductivity Anomaly The prediction of the anomalous thermal conductivity of simple substances has been investigated by Sengers (refs. 29, 30, and 32), Brokaw (ref. 25), and Hendricks and Baron (ref. 31). While each technique reproduces the CO 2 data of references 29, 30, and 32 with good agreement, the simplified technique of reference 32 is used herein. 42 Define the parametric groupings: (B55) (k - _F ) _ (B56)y= 3.o_1o-5I1- TRI-°"6 where and _F _F = f(p) + k* is the frozen or nonanomalous part. (B57) Within the region 0.4 -<PR -< i.6, the anomalous thermal conductivity can be com- puted by using the following formulas: For x2 < 0.4: For 0.4-<x fl-< 3: For x/_ > 3: _r = k - XF = 3.05'<10 -5 p_ (B58) (1- pR) 1" 71 y[109]_0.6(x_)t/o.35 (B59) A r 3.05x10 -5 _ (B60) =k -_F- 11-TRI°'6 43 In addition, the critical constraints usedare Pc = Pc(Pc'Tc ) (C6) C O2p = 0 = ac as determined by the saturation data (C7) (C8) (C9) Two additional constraints = 0 (C 10) = 0 (Cll) are sometimes needed to obtain a good fit in the critical region. These constraints are referred to as the fifth-order attachment, while equations (C7) and (C8) give the usual third -order attachment. The coefficient n20 in equation (B2) is not determined in the fit but must be varied systematically until an acceptable set of coefficients is determined. The data needed for this method then are (1) PVT data over the pressure and temperature range of interest (2) Smoothed PVT data along the saturation boundary (3) Estimates of the critical parameters and ac (Initial estimates can be found from the saturation loci. ) 2_1. (4) Value of n20 in equation (B2) (Initial estimate satisfies n20 Pc The items (3) and (4) along with the weight factors can be varied to obtain the best set of coefficients for equation (B2). Techniques to actually solve the system of equations are in the literature; see, for example, Hust and McCarty (ref. 45) and Morsy (ref. 65). 46 APPEND_ D STATE EQUATION OF STEWART, JACOBSEN, AND MEYERS Stewart, along with Jacobsen and Meyers (ref. 44), extended the work of Coleman and Stewart (ref. 37) to nitrogen and oxygen. As Bender's analysis and program was made available to Stewart and his coworkers, they utilized Bender's technique to develop another state equation. The state equation is essentially equation (B1) where N = 9 and M = 6, with a %/T term. The explicit forms for Ai(T) and Bi(T) become AI(T) = RT n 4 A2(T) =nlT + n 3 + T n 5 + -- + E I(T) T 2 n 8 n9 A3(T) =n6T+ n 7 +-- +-- T T 2 nl 2 A4(T) = nl0T + nll+ __ T A5(T) = n 13 A6(T) = A7(T) = A8(T) = A9(T) = n14 n15 T T 2 n16 T nl 7 nl 8 + T T 2 n19 T 2 (D1) 47 B I(T) n20 + n21 T 2 T 3 n22 B2(T ) - T 2 B3(T ) n24 + n25 T 2 T 3 n23 +-- T 4 n26 n27 B4(T ) - + T 2 T 4 B5(T) = n2___88+ n29 T 2 T 3 B6(T) = no + n___1 + n3____2 T 2 T 3 T 4 (D2) where El(T) = n 2 The exponential term of equation (B1) becomes 2 -n33P e (D3) (D4) The explicit forms for Ci(p) and Dj(p) become 48 APPEND_ E THERMODYNAMIC RELATIONS AND DERIVATIVES The symbols Cp, Cv, H, P, R, S, T, and p have the same meaning as defined elsewhere in this report. The other symbols used exclusively in this appendix are de- fined as follows: A = E - TS Helmholtz free energy or work content E internal energy F = H - TS Gibbs free energy or free energy K equilibrium constant V specific volume To illustrate the facility of the partial derivatives, Roder and Weber (ref. 66) give five which are useful to engineers: Specific heat input: (a_P avJp pCp Energy derivative: _-_/V pC v , __T_'p Isothermal bulk modulus: \aV/T \bo/T 51 Volume expansivity: • aP The background material necessary to derive these and other parameters as the Joule-Thomson coefficient aT)H = 1 can be found in most thermodynamic texts. GASP provides the partial derivatives (aP/_p) T and (_P/_T)p. With the aid of the following thermodynamic derivatives and the Bridgeman Tables, any thermodynamic parameter can be found. The following thermodynamic tables were taken from refer- ence 67. Differential energy formulas: dE = TdS- P dV dH = TdS +VdP dA=- S dT- P dV dF =- S dT+VdP Maxwell relations: S ,as Iv 52 T V T P Energy-function derivatives: V P S T P V Heat-capacity relations: p \_T]p _ ___v_ _(._ % Cv [_JP[_JT Cp-Csat = T(_V_ (_P_ \aT/p\aT/sa t 53 (aS)v = - (aV)s = T[ P_-_]T \aT/pJ = =c /av_ +T(_] _ (aE)v -(aV)E P t'_-P)T \aT/p : :C I_V_ T(_]__ vpV_ (aH)v -(aV)H P_aP/T + \aT/p t_-T/p ] (aF)V -- - (aV)F L \-_/p taP/T] = =_spV] (aA)V - (aV)A _-_'# T P Pc {av'_ (aE)s:-(_S)E:7 [ PbP]m VCp (aH)S = -(aS)H = T _[v -sT(_'_I(aF)S = -(aS)F = - -_ Cp \aTipJ (aA)S = -(aS)A = T i"\aP]T \aT/pJ [-_/pJ <_.)..:-<_..),=-v[c_-(m,_,_jl-_ro[i>t77)T/.av_+m(aV]"]k7711>J <_>..:-<_.,_:_v[c_,(_v_1+_h(_v_+_(_v_1t> t77/t>J L _77/p t-D7/_j (aA)E=-(aE) A P Cp t-_l T taTIpJ _'TIp 56 _- __ ms(_V_ (a;)H -(aH); -v(cp+s)+ [_TIp JTipl xaT/pl Pt_)_ (aA)F -(_F) A -S[ taP/TJ _'T/p 57 APPENDIX F DESCRIPTION OF SUBROUTINES IN GASP This appendix includes a discussion of the input and output and important features of the major subroutines in GASP. The equations are presented in appendix B. The FORTRAN IV variables mentioned correspond to the program listing in appendix H (see also symbol list, appendix A). GENERAL MATHEMATICAL ROUTINES Function SOLVE(X1, F, DF) Function SOLVE performs a Newton-Raphson iteration given the initial estimate X1, the function F, and the derivative function DF. The convergence is determined when The value of TOL is 1. E-5 for iterations 1 to 40, 1. E-4 for 41 to 60, 1. E-3 for 61 to 80, and 1. E-2 for 81 to 100. In all cases studied, the convergence was usually obtained in fewer than 40 iterations. For the exceptions, usually in the near-critical region of the PVT surface, the values returned with the increased tolerance are the best obtainable us- ing equation (B3). The maximum number of iterations is 100, and an appropriate mes- sage is written if this number is reached. Subroutines ROOT(X0, X2, FOFX, FUNC, Xl) and ROOTX(X0, X2, FOFX, FUNC, Xl) Subroutines ROOT and ROOTX are identical except for name. The duplication is necessary for the double iterations in the solutions for temperature and density given pressure and enthalpy (KS=4) or pressure and entropy (KS=5) as input (see also table V). The solution method is a modified, half-interval, search technique for a monotonic function, FUNC, with a root between X0 and X2 such that FUNC(X1) = FOFX, where Xl is the answer returned. The number of iterations does not exceed 100, and the tolerance is varied in the same manner as in function SOLVE. In addition, both the root and the 58 (52, element number) (53, element number) DIST(element number) EPSOK(elementnumber) HSCHI HSCH2 (2, element number) (3, element number) OM2OMl(element number) WTMOL TCCOF(15, element number) is stored as (I, element number) XLAMB_ f zc5 ) TCR (4, element number) RHOCR (5, element number) (14, element number) (15, element number) AKSTCO(18, KODE) CKMKST KODE TCSTAR ZETA FF DTRIPL SWT KSWT minimum enthalpy, J/g maximum enthalpy, J/g intermolecular distance potentialparameter used to reduce temperature, k/_ correction factor for r molecular weight used in k - k* reduced correlation for thermal conductivity; Z also c used in surface tension temperature at thermodynamic criti- cal point, Tc, K density at thermodynamic critical point, pc , g/era 3 (k - k*)),Z 5 coefficients for POLY C index for k /k T curve C k T , cal/sec-cm-K C k /k T curve coefficients for POLY c KODE= 1 for inert gas =2 for diatomie =3 for carbon dioxide =4 for methane (CH4) viscosity normalizing parameter, ( surface tension multiplier, C density at the triple point lower temperature switch to numeri- cal specific heats use numerical techniques on neon, carbon monoxide, and helium where temperature or density switch is set 61 DIFTT RHOSWT AMUC, BMUC temperature increment for evaluating numerical specific heats density delimiter for switching to numerical specific heat calculation viscosity correction parameters at high density (eq. (B38)) Function CHECK Entry: DCHECK(KU, D). - At this entry point the density D is converted to g/cm 3 and checked to see ifD is out of range. Entry: PCHECK(KU, KR, P). - This entry converts the pressure P to MPa and checks for out of range. IfKR=I, P is checked to see ifitis out of saturation range. If P is out of range, the program writes an out-of-range note and continues. Entry: TCHECK(KU, KR, T). - This entry converts the temperature T to Kelvin and checks for out of range. If KR=I, T is checked to see if it is out of saturation range. If T is out of range, the program writes an out-of-range note and continue. These entry points convert the variables from the user's units to those represented by KU=I and check for out-of-range variables. Appropriate messages are written for any out-of-range input, but the calculation is allowed to continue. SUBROUTINES TO OBTAIN STATE VARIABLES(KS OPTIONS) The following subroutines use the mathematical routines, function CHECK, and sub- routines listed with each (see the modular structure of GASP, appendix G and table XI). Their use is determined by the KS ahd KP options (table V), and they are called by sub- routine GASP. If a user wants to use only a few of these subroutines, he can disassem- ble the GASP program by following the instructions in appendix G and the discussion for the routine of interest. These subroutines are called twice by GASP for saturation prop- erties, once with DL and once with DV as input for D. Subroutine DENS(KU, T, P, D, DL, DV, KR) Subroutine DENS solves equation 03I)for the density, given T and P in units indi- cated by KU. The region number KR is returned, and the density is returned in D for KR=2 or KR=3. For KR=I, the saturation values are returned in DL and DV. IfKR=I 62 for input and either T=0. or P=0. for input, the saturated value is calculated and re- turned for the variable which was input as 0. The solution is obtained by ROOT for subcritical pressures and by SOLVE for sat- uration or supercritical pressures. Special initial estimates were found necessary for convergence with SOLVE in the near-subcritical temperature region when 0.9T c <T <0.9999T c. Subroutine PRESS(KU, T, D, P, KR) Subroutine PRESS solves equation (B 1) for the pressure as a function of T and D in regions KR=2 and KR=3 and as a function of only T in region KR=I (using subroutine PSSS). The result, P, is returned in user's units indicated by KU. The correct value of KR is also returned, and the calculation is direct. Equation (B1) is repeated in its expanded form for N= 6, M = 2: P = RpW + A(T)p 2 + B(T)p 3 + C(W)p 4 + D(T)p 5 + n13 p6 + p3 EE(T ) + p2F(T)_e-n20p2/T2 where if the fluid is helium, T 2 = T, otherwise T 2 = 1; and where n3 n 4 n 5 n21 A(T) =nl T +n 2+u+-+n+ T T 2 T 3 T 4 n 8 n23 B(T) =n22T2+n6T +n 7+_+ T T 2 C(T) = n9T + nl0 n24 D(T) =nll T +n12 +- T E(T) = 63 solve for temperature TS, given pressure PS, and density DS: DTSF - 8(TSF) 8T SubroutineTEMPPH(KU, P, H, T, D, DL, DV, KR) SubroutineTEMPPH solves equation (B20)by using equation (B10}for the temper- ature parameter T anddensity D, given P and H as input in user's units indicated by KU. The doubleiteration is performed by using ROOTand ROOTXwith function TSHF for re- gions KR=2and KR=3. In region KR=I the saturation values are determined for DL and DV by DENS, andT is foundby function TSS(using SOLVE). KR is also returned. SubroutineTEMPS(KU, P, S, T, D, DL, DV, KR) SubroutineTEMPS solves equations (B21)and (B16)for T and D in the samemanner as TEMPPH, by using P andS as input and function TPSF for the double iteration with ROOT and ROOTX. Function TSHF (TS) Function TSHF is a multiple-entry routine used in conjunction with TEMPPH. It ob- tains a trial value of DS by using the given PS and a trial TS. Then it obtains a trial H, which is compared with the input enthalpy within an iteration in TEMPPH. Entry TPSF(TS) Entry TPSF is a function used in conjunction with TEMPPS. It obtains a trial value of DS from the given PS and a trial TS. Then it finds a trial S, which is compared to the input entropy within an iteration in TEMPPS. SUBROUTINES TO OBTAIN DERIVED THERMODYNAMIC PROPERTIES The following routines assume that the variables T and D have been input or pre- viously calculated in the user's units. This condition is satisfied in subroutine GASP. 66 WhenKR=I is input or hasbeenso determined, the corresponding saturated variable (HL, SL, etc., or HV, SV, etc.) is output. Subroutine ENTH(KU, KR, T, P, D, H, HL, HV) Subroutine ENTH calculates enthalpy H in user's units (KU) by using equation (B20). SL and SV are returned for saturation conditions, KR=I. Subroutine ENT(KU, KR, T, P, D, S, SL, SV) Subroutine ENT calculates entropy S in user's units (KU) by using equation (B21). SL and SV are returned for saturation conditions, KR=I. Function HDINT(DS, DSL) Function HDINT is a multiple-entry routine which computes the integral used in the enthalpy computation from density DSL to density DS. DS Entry SDINT(DS, DSL) Entry SDINT computes the entropy integral from the density DSL to the density DS. DS 67 Function HDINTF(DS) Function HDINTF is a multiple-entry routine which evaluatesthe integrand of func- tion HDINT, where DSis the variable of integration. 2 2 P P Entry SDINTF (DS) Entry SDINTF evaluates the integrand of function SDINT, where DS is the variable of integration. Function HSS(PS, DS) Function HSS is a multiple-entry routine which computes the enthalpy in the region KR=3 (table V) or saturated vapor enthalpy (KR=I), given pressure PS, density DS, and temperature TS. H =H 0+ Cp,0dT + (Z - I)RT + P T 0P dp TO p2 T where Cp, 0 = _ mi Ti-1 j=l and the values m 1 to m 5 are stored by subroutine SETUP, except those for hydrogen. 68 While direct calculation of CV may be valid at densities greater than the switching den- sities for some fluids in GASP, the numeric method is preferred because of the some- times erratic behavior of aP2/ST 2. In addition, if T < Tc, the specific heats for neon, carbon monoxide, and helium are computed by SPCHV. The P-T diagram given in fig- ure 29 will give the user a better feel for the physical regions involved. The first partial derivatives from subroutine PTRHO (eqs. (B9) and (B12)) are made available through COMMON/DERIV/PDT, PDTL, PDTV, PTVL, and PTVV. These are with PDTV for saturated vapor and PDTL for saturated liquid and PTV = 8(_T) p with PTVL for saturated liquid and PTVV for saturated vapor. Function C PPRLF(DS) Function CPPRLF evaluates the integral used in computing specific heats: Subroutine PTRHO(D, T) Subroutine PTRHO evaluates two partial derivatives. COMMON/PARTLS/PTV, PDT. PTV = a(_)p The answers are returned in 71 Subroutine SPCHV(KU, KR, T, P, D, CV, CVL, CVV) Subroutine SPCHV computes the specific heat at constant volume CV given temper- ature T, pressure P, and density D. If KR=I, the saturated liquid or vapor specific heat, CVL or CVV, respectively, is computed as requested by GASP. CV- AT AT See subroutine CPPRL. Subroutine CVPS(KVP, KR, CVS) Given temperature, pressure, and density, this routine is used to determine five values of internal energy U for a spline fit used in SPCHV to compute the specific heat at constant volume. Subroutine SONIC(KU, KR, T, D, GAMMA, C) This routine calculates the sonic velocity C, in user's units KU, by using equa- tion (B26) and the ratio GAMMA=CP/CV, where GAMMA is calculated in subroutine GASP. Subroutine VISC(KU, KR, T, D, MU) Subroutine VISC uses T and D as input in user's units KU. Dynamic viscosity is calculated by using one or more of equations (B27) to (B36) depending on the fluid, as explained in appendix B. Calculations of dynamic viscosity are direct evaluations of curve fits but may be altered as by equation (B38). A special form is used for hydrogen and helium. 72 Function VISCD(DS, TS) Function VISCD, developed by McCarty of NBS, computes the dynamic viscosity of helium. Subroutine THERM(KU, KR, T, P, D, EXCESK, FK) Subroutine THERM uses T, P, and D in user's units KU to calculate the thermal conductivity FK in user's units KU. An optional coding section calculates the critical excess thermal conductivity EXCESK associated with the critical anomaly in the PVT region, O. 4 < P/Pc < 1.6. See also references 25 and 29 to 32, the subroutine listing in appendix H, and equations (B55) to (B60). EXCESK will be computed for the density regime 0.4 < P/Pc < 1.6. To obtain the thermal conductivity of a near-critical fluid, the user must add FK and EXCESK. The equations used for the thermal conductivity of the different fluids are (B39) to (B54). Subroutines CONZ and CONC are used for the special forms for parahydrogen and helium (eqs. (B45) to (B47) and (B51))o Subroutine SURF(KU, KR, T, SIGMA) Subroutine SURF uses T in user's units KU to calculate the surface tension SIGMA of the liquid by using equations (B61) and (B62). Subroutine DGUESS(TS, TCR, DST) Subroutine DGUESS provides an estimate of near-critical density that is used by subroutine DENS. Function CONZ(TEMP) Function CONZ computes the thermal conductivity of dilute gaseous helium as a function of temperature by the technique of H. M Roder of NBS, Boulder, Colorado (private communication). 73 P, p, T calculations (KS=l, 5) KP = 1 E nthalpy KP= 2 E ntropy Storage saved, 66318 (3) E nthalpy-entropy calculations P, p, T calculations (KS=I,3) = t I IStorage saved, 131378 (4) Specific user requirement - specific heats, sonic velocity, and derivatives (a) If p <--Pswitch (see table on p. 70 of appendix F), P, p, T calculations (KS=l,3) KP=4 Specific heats, sonic velocity, and derivatives ri co it \ Storage saved, 111758 76 (b) If p > Pswitch' P, p, T calculations (KS=l, 3) KP= 1 E nthalpy K-P=2 Entropy KP=4 Specific heats, sonic velocity, and derivatives (5) Transport property calculations Storage saved, 46508 P, p, T calculations (KS=I,3) KP=8 Viscosity KP= 16 Thermal conductivity KP=32 Surface tension Storage saved, 76128 The specific illustration of removing all but one fluid from subroutine SETUP and BLOCK DATA is given as follows: 77 TEST PROGRA_ PREPAFES TABLES FOR EACH FLUID OVER A RANGE OF PVT COMMDN/PROPTY/KU_DLvDV_HLwHV_S_SLtSVmCVmCVL!'CVVtCP_CPLtCPVmGAMMAt [GAqMAL i GAMMAV i C i CL _CVP t MU_ MULL MUV_ K _,KL t KV_ SI _ MA pEXCES K, EXCLt'EXCV COWMDN IDER IV/PDT_ PTV_ PDTL tPDTV, P TVL t PTVV REAL MUtMUL _MUVt K_KL_ KV CO_4M3N/CHECKS/DCHlpDCH2_PCH I_PCH2_PCH3, TCHI tTCHZtTCH3tDST ,TST m H ISCHI_HSCH2 DATA NAM/3HCH_t2HN2_ZHO2_2HARI3HCOZ,2HNE_ZHCDIZHHEt2HFZtZHH2/ DIMENSION NAq( IOIBP(3I_T(250}IYPL(2501[O} pTSIART(3tIO)tDTA(3elO) DATA Pll._ [.P 100.I DATA TSTAPTIg5.tZOS.t4OO. m65.,200, t400. _bO. p225.t_OO.,BS.t225.,400 _'.l, 2ZO.t, 420.,50O*l'25.,80._,400.,70-tZOO.'_O0° _3.tt3. tI30*_70.,,203._, *_O0._ 16._40,t 200o/ DATA CTAIiO.t, 20.t [DO.t, IO. _,ZO. _,1DO. _lO. m20. _,I)O. _I0. s,20. tl30.t * 2., [0._ I00.I KS=[ KP=3[ KU=I WPlTE( 6_ II I FORMAT( IHII DO 1000 t=lOt I0 CALL SETUP (NAN[LIt P( 2 I=PCMZ/. lO 1325 T( I )=TSTART( l_t | OT=DTA(ILL) T S=TSTAF, T( 2,L ) N=Z K PT =0 DO 20 J=IvSO K PT =J T(J÷II=T(JI+DT TJ=T(J+I) IF(TJ,GE°TCH3] E_TO 25 IF(L.EO. IO .AND. TJ._E.Igg9. S) GOTO 25 I F( TJ-TS+,05I ZO_ 210_, 210 ZlO T(J+I)=TS DT=DTA(N_[ ) N=N +]. TS=TCH3 IF(N.LI._) TS=TSTART( N,L l ZO CONTINUE Z5 KPT=KPT+I T(KPT ) =TCH3 DO I00 I=I_3 _= Pl I 1,_.101325 DO 50 J=ItKPT KR=O CALL C_SP (KS_KP_T(JJ_ZtD_H_KR) YPL(J_ I)=D YPL (J_, _)=H YPL(J_ 3)=S YPL (Jr 4 )=CP YPL(J_ 5)=CV YPtlJtE|=C YPL(J_7)=K YPL( Jr 8 )=MU YPL(J_ 9)=PPT YPL(J, IO )=PTV _.6 FORMAT(IH tIb_bE14.6) 50 CONTINUE C TABLE OF C CENSITYm ENTROPY_ ENTHALPY_ CP_ CV, S(_NIC VELOCITY C THERMAL C3NDLCTIVITY, VISCUSITY AND PARTIAL DERIVATIVES C FOR EACH T P WRITE( 6,800( NAM(LItKU 800 FORMAT( IHI _5X_34H THERMODYNAMIC PROPERTY TABLE -- ,_A_ _3X,3HKU= 1 13} WI_ITE(E,80ll P[I ) BOl FORMAT[LHO _5Xt FLU,Z ,IIH ATM ISOBAR) WRITE( 6,701 I WRITE( 6_ 702} 701 FORMAT( 3HO _ 4HTEMPt 5Xt 7HDENSI TY _6X _ BHE NTHALP Yt_tX t 7HENT ROPY_TX t 2 HC =ptloxt 21,_cvtgX_THSDN VELI_X_,tOHTHERM CONOISXI_HVISCtTXtSHDP/DDITX,_5 $HDP/DT_/] 702 FORMAT ( 4X _ IHK t 8X,_4HG/CC t 8X_, 3H JIG _ 8X tSHJ IG-K_I Xt 5H JIG-K vTX _ 5HJIG-K_ ,_=8X,_.HCMIStTX_8HJIE'4-S-K_SX_bHG/CM-S_4X_IiJHMN'CC/G'MZ_3X_IHMN/K-N2m *Ill DO [50 J=I-KPT WRITE (6,803} I IJI,,(YPL(Jt NI_,N=J.,IO! IF (4OC(J,501.NE°O| S{1 TO 15_, IF ( J.EO. KPT) GO TT_ lbu WRITE( 6t80Ol NAM(L } WRIT_(6_801) P(1) '78 U 3 0 p. N IE O, eq I I$ I I I _ZZZZZ u_O 00 O.U.U. U_ _ _ i,'- uJ a._ t_ W _ U,. I.L U_ U__ _;. un _n _ m u_ _mmmmm ,,,75' X: t-- I I < u It" 3: I,,- ooo __ooooooo ,,, ,?,,,?, ?,,?, ,, ,_ ,,,, ,0,,,,,,,? ??? ??? o .... o ?? dd;dggddgdgddgg;4ggggd;g_gdggdgdddggdddggdgg_g4 I II,llllllll?OOlllllllllllllllllll,??O0?llllllllll II N 00000000000000000000000000000000000000000000000 • oo0000000000 8_ Z'' ' ''_Z_±__' '_&_, , • i i !*1tlt1111_ii.I IIIiii I i_ I JgJgdoddggoJgdgodgdoooooooooooooooJooo_ooogoooo ooo7oo_oooooo ooo?ooo_oooooo_oo_oooo ooo 00000000000000000000000000000000000000000000000 co u_ __.=_ o• _Z _. _u_ _NN_ oooo o_o?_oo_ ÷++÷*÷+÷* ÷+$ III 00o0000 O000000_O0000000 000_0000000 O000 _O0 00o00_00000_00000_00000_0000000000000000000000 000000000000000000000_0000000000000000000000000 _ .....__R_oo_ooo__ _o~_°°°°_°°°°°°°_°°°_ooooooo o_o_oooo°°°°_ 81 _ ?? ??? ??? ????_???_???????,,, ,,?, O_O0000000000000000000000000000000000_O00000000 ?_?_?_?°°° ,,, ,,,,,,_°°°_°°° ,*=,_= , • •ooooooggo;ogo;gd;'''=oo ooooooooooo_ _ooo oo I I II oogoooooooooooooogoo ooo oooo*oo _ooooeeeeeeeeeee_e #eee =e=eeeeee=e=e==$*eee 00000000000000000000 =, ,_,,=e ,,= II|IIIIII O00000000000000000000000000=m_m,_*,,, ,,m, ,** m*m eeeeeee_eeee_eeeeee 000000000000000000000000000 00000000000000 _NN_NN_O00000000_O00 00000000000000000000 &_,_, O0000000000000000000 ooooooocooooooooooo_oooooo_ U l M _u UJ lu I _1 I-- I_ I- Ll u I1. Z Z Z ooooooooooo_oooooooo ; ; ooo ooooooooooooooooo o o 82 z / l,l.I _ _-.> z_ .l QI,.l "_ I/I ,¢I I-- 8 i??l???? *+1+.11++11.1+!1+.***.1,.**111+++_..++oo-oooooooooooooooooooooooooooooooo_oodddd_dddo I ÷ . _,,,_,,,_,,,_,,, ,,, ,,,,,,,,,,_,,,,,,,,_,,,_ ;_ooooooo_oo_oo_;;o;;_;;;;_;oo;;o;;oo............. __g_ 3_ooo ooog_ooooooo_g_ooo ooooooooooogoooo_oooooo_ liil _g___g____33_33_3_o _oooo ooooooo ooo?ooooooo oo oooo oo ooo oo _ ...... 000000000000 IIIIIIIIII g _ &_lll_lll_ .... IlI_'IIIA_=IIIAIII_ll_IIIlIIIIIII 00000000000000000000000000000000000000000000000 83 TO I01 2O 25 3O 35 SATUR ATIflt,l CURVE COLFFIC I_ NTS Cf:MM(]NICOSAT/ CPSI( I)tCPSZ,CPS3tCPS4,CPS5tCPSb,CPS7 CP,Q CURVE COEFFICIENTS CO_MON/COCP_ITO( IItCCOP_,CCIIP2_CCOP31CCOP4tCCF]P5 REFERENCE ENTROPY_ ENTHALPY t CP CORRECTION FACTOR COMMDN/REFND/SOT_ mHOTn tCPOCnR 7HTERM_ STERM PARAMETERS FOR CHECKING REGION ANO RANGE LIMITS ON DENSITY PRESSUPE_, TEMPERATURE, EHTHALPY. DENSITY AND TEMPERATURE ESTIMATES FOR _EWTON-RAPHSC)N ITERATIONS COMMON/CHECKSIDCHI( I)_,{)CH2_PCHI,PCH2,pPCHO,TCHImTCH2pTCH3t, DSTt, TSTt, H IS CHI, HSCH2 CONSTANTS F._R THERMAL COND(_TIVITY CALCULATION. COMMr]N/TCt_ND/CKMKST(QIt XLA_B_7. C51, RHOCR_ TCRt, TCSTARt. CKSTAR(I81 CJNSTANTS FOR VISCOSITY CALCULATION. CO_'40N/COFMU/EPSK vWMt,D IS pRHOCRT _ZETAA t,AMUX ,BMUX CuNSTANTS FIR SURFACE TENSION CALCULATIUN. COM_ONISURCOV./PCTCt T'..RIT,FI XI Tt ZET S_ITCF FOR CP-CV CALCULATION. VALUES IN CERTAIN REGIONS ARE CA. CULATEO BY _UME_ICAL )IFFERENTIATIUN ...WHENEVER THE ANALYTICAL DERIVATIVES ARE GO]Dr THEY ARE USED ......... CO_]N/SW ITSIK SWI T, TSWIT,D IF T,R SWI T CO_MON/SOL ID/DTRIP < GAS=O 00 I0 l=Ip I0 IF(NAMC.AS.E_.MATCH(II ) KGAS=I CO_T INUE IF(KGAS.FQ.O! GO TO TO I HE=O IF (,',GAS.EQ.B) IHE=L IFL=O IF (KGAS.EQ.9| IFL=I IF (IFL.EO.I! WRITE(6,72| I HY=O I FI KGAS.EQ.IO I IHY=I WRITE(6_,IOll IMESSAGI ItKGASItI=I,15) FORMAT ( IHO, 15A6) STORE CONSTANTS FOR BENDER-S EQ. OF STATE DO 20 I=1t25 R(I)= CO F( It, KGAS I STORE SATURATION CURVE FOR GAS DO 25 I=26,32 CPSI{ 1-25I =COF(I,KGAS) STORE CPO COEFF IC I_:_TS ANO REFERENCE TO DO 30 I=33t, 38 TO(1-32) = COF[In, KGAS) STORE REFEREN"- _ E NTROPYt, E NTHALPY, CPO CORRECt I CN FACTOR TX=TD ! II COMPUTE LOWER BOUND F]R THE ENTHALPY AND ENTROPY CALCULATIONS IN HSS A I=CCOP 5".25 A 2= CCOP to*.3333333 A3=CCOP3*.5 STERM =(((AI*TX+A2)tTX+A31*TX+CCOP2)*TX+CCOPItALOGiTX) A I=CCOPS*.2 A 2= CCOP 4't'.25 A3= CCOP__*. 3333333 A4=CCOP 2. ,5 HTER_ =(( ((AIW'TX+A2)W'TX+A3)*TXc-A_)tTX+CCOPI)*TX 10,312 = _.q68_,18_/2.01572 IF(IHY.EQ.I) STERM=I_.312,WALOG(TXI IF( IHY.EQ.I! HTERM=IO.312*TX SOTO= COF(39_,KGAS) HOTO= COF(_O_,KGAS) CPOCO R=CO F( _I_KGAS! STORE CRITICAL VALUES AND REGION BOUNDARY CONSTANTS DO 35 I=_2,53 DCHI{ l-_il = COFI I,KGAS) XLA_B= TCCOF( IoKGAS| ICS-- TCCOF(2, KGASI TCR= TCCOF( 3,KGAS l RHOCR = TCC3FI4_KGAS) LOAD K-K* CURVE COEFFICIENTS DO t¢5 I=5, 13 86 • 5 CKMKSTI I-6| z TCC_F| ItKGAS) TCSTAR- TCCDF( [SpKGAS| I m TCCOF(14_KGASI_.[ LOAD KSTAR/KSTARTC CURVE COEFFICIENTS DO 50 J"IB 18 50 CKSTAR(J l= AKSTCO|J_I I STORE CONSTANTS FOR U-Ui, AND VISCOSI TY CALCULATION WM=WTMOL|KGAS) Z ETAA-ZETA( KGAS ) RHO CRT-TCCD F( 4t KGAS | EPSK zEPSOK(K GAS ! O IS:,,O IST( KGAS ) STORE CONSTANTS FOR SURFACE TENSION ( RF IS RIEDEL FACTOR) T CR IT -TC H2 PCTC=| PCPZ/.LOI325)i,4,.6(_6666T*TCH2i_,I ,.3333333. FIXIT-FF(KGAS) Z ET,,Z C5"* .2 ENTER SWITCHING PARAMETER FOR CP,CV CALCULATION K SW IT=KSWT(KGAS) T SWIT'SWT(KGAS I IF (KSWIT.EQ.2) TSWIT=TCH2 D IFT_DIFTT(KGAS ) RSWIT-RHCISWT(KGAS) *RHOCRT DTR IP =CTR IPL( KGAS ) AMUX=AMUC( KGAS ) BMUX= BMUC( KGA S) RETURN 70 WRITE(AtTI} 71 FDRMAT(IHOo 6QH ERROR IN CODE FOR NAMGAS - NO CnNSTANTS STORED- 2 PRDGRAM STOP. ) ?2 FORMAT(IHOtl[SH Th1_ REGION 125-L65K FOR T AND P GREATER THAN 10 IATM YIELDS POOR RESULTS FOR THE DERIVED PROEPRTIES--BEWARE. STOP END C STORES T_E COEFFICIENTS FOR ALL FLUIDS FF)R TH_ EQUATION OF STATE C AND THE TRANSPORT EQUATIONS. C STORES CHNVER SION CONSTANTS NEEDED BY ALL FLU) f3S C BLOCK rATA COMMON IGASES/MATC H( ]._ ItME SSAG( 15, L O! COMMON/ CON123/ DCONV(S|e TCONV( 51tPCIJNV(51 COMMON/CONV4/SC()N V(5) C OMMON/CONVS/CCNN V(5) COMMON/CONVA/HCONV( 51 COMMON/ALLCF)F/COF(53t [O}_TCCOF( [5tlO! tAKSTCI)(18 _4) pDIST(I} It IWTHOL ( IO)_EPSOK( TO! t(ETA( IOItFF (IO) ,SWT(IO),KSWT(IO)tDIFTT (IO 1, 2RHOSwT( IOItDTRIPL( IO)tAMUC(IO),BMUC(IO) rtIMENS ION MCH4( 15!,MW2( L5)tMOZ( 15! tMAR([5) tM=O2 (15) tMNE(L5 ) tMCO(15 I!,MHE| 15)tMF2| 151tMH2(15) EOUIVALENCE(MESSAG( It I|,MCH4(I) ! _(MESSAG(It2| _MN2(I!| t ( MESSAG| ],3) It M]2( l))t (MESSAG( Lt4)tMAR( Ill t|MES_AG(I t5| ,M=02 |11 I t(MESSAG| I_O)t 2MNE|) I ItIMESSAG( lpTItMCO( I! )o(MESSAG(ItS! tMH--- (I)! t (MESSAG( [tg!t MF2 3([) ), | MESSAG( It 101,MH2 ( 11! DATA (MATCH( I|_ I=ItBII3HCHt, t2HNZt2H02 t2HARt3HCCZtZHNEo2HC(]tZHHE/ DATA MATCH( 9 )/2HF2/ OATA MATCI'( lO )/2HH2/ DATA wCH¢| ]) /(9OH THERMODYNAMIC ANL) TRANSPCRT PROPERTIES FOR ME ITHANE PC=45.66 ATM,T:=I9O.T7 K_ROC=.I02 C,/CC / DATA MN2( 1! /90"I TH-R_DOYNAHIC AND TRANSPORT PROPERTIES FOR N IITROGEN RC=33.72ATM,TC=126.3 KtROC=.3105G/CC / DATA M02(I) /90H THERMODYNAMIC AND TRANSPCRT PROPERTIES FOR OX IYGEN PC=SO.IbATMtTC=154.78 KtROC=.4325C,/CC / DATA MAR| I! /C)OH THERMODYNAMIC AND TRANSPORT PF_OPERTIES FOR AR tGON PC= z,B.O14AT_t T:=ISO.7 K,RPC=.531 G/CC / DATA MC02( [l /9DH THERMODYNAMIC AND TRANSPCQT PR(]PERTIES FOR CO 2 PC=72.B6qATMt TC=304.21. K, RUC=.464 G/CC / 'ATA MNE| II /qOH THERMOr_yNAMIC AND TRA_SPC_RT PROPERTIES FOR NE 87 ].ON PC= 26.19ATM ,Tc=64-6 K, RHOC=o483 G/CC / DATA MCOI 11 IgOH THERMODYNAMIC AND TRANSPORT L PC=36°52gATMt T.C=132.91 K, RHOC=.2997 G/CC I DATA MHE(ll /90H THERMODYNAMIC AND TRANSPORT I PC = 2.265ATM, TC = 5°2K, RHOC = .0693G/CC / DATA MF2(|| /90HTHERMDDYNAMIC AND TRANSPORT IORINE PC=BI°47ATM, T==1¢¢.31K,RHOC=-5738G/C C DATA MHZ( L)/ 90HTHERMODYNAMIC AND TRANSPORT LG E_IPC" 12.759& TM t TC =32 • g76Kt RHOC =. 03163G/CC DATA TCONV /1o11. ,1.8,2_'1./ DATA PCONV/t.tD-86923Z7,165-03776 ,24'1./ DAtA DCONVlZSl.t62.62796t2w_1-/ DATA SCONVI2WW1.,O.23900574,2*l./ OATA HCONVI2el.tO.6302|03,2*1./ DATA CCONV/2_I.tO-O3ZBO84,2_l./ MET PAN E- .................. BENDER 1971 DATA (COFI I, 1 ), I: 1, 5;11 .518251 1- .B6366_O2E3 , .25005236E 5 ,-.125338_8E8 2 .T5523689 ,-.12111233E3 , .20547188E6 3- .61948317E6 t- .25603803E 2 , .11556713E6 6 ._949_1630E8 t-.I1956135Ell , .88302981E12 5 .12302028E13 t-.64499295E14 , .37E2 ,4_'0. t 6 .29976933E3, .70773219 , -.T23_gLZZE-Z 7-.1213063gE-0, .I_8489295E-9 , .11167E3 8 .39225023E-2 ,-.56781803E-4 , .234_2607E-S 98.355173_6 , .803207E3 , .260815 I .57 , .01167 ,_,.62T 2 .gO66E2 t .19077E3 , .bE3 3 .25E3 t 220.28242 ,2007. / NITROGEN .................. BENDER IgTl PROPERTIES FOR CO PROPERTIES FOR HE PROPERTIES FOR FLU / PROPERTIES FCR HYDRO / t .ITI91020EL t .34169567E9 t , .32337543E2 , . 27_25297 E5 t-. 31713_86 EIO -. 32097996E2 . zt0851153E-4 , _1.9698625Z ,-. 223950[T E-'9 , . I E--*, , . 5067E2 , .45 DATA (COF(I,21,I=1,531/ .296797 , .48058175 , 1-.1506TOIOE3 , .26071365E4 ,-.12792742E7 , .29_36228E8 , Z .37500265 t-.SOT38465E2 , .1_699236E5 , .1_*0998E1 , 3-.2_136776E3 t .2895/_771 ,-.27613799E2 , .360/'826_E3 , 6- .200 83357E7 t .4326518q_: 9 t-.16513521Ell t-. I0141_39E8 , 5 .4721_396E10 ,-.12016418E12 , .LE2 , 4"0. , -.66869126E1 , 6-.20086668E3 , .37077481 ,-.62602129E-__ , .56187159E-4 , 7-.25986366E-6 t ._gO25967E-9 , .77364 E2 , .29093035E+2 8 .8993017OE-3, -.92391t038E-5• .26008592E-7, -.I4102926E-13, 9 .Z3351167 E1 , .231188578E3 , .35693888E-1 • o14e19 E-3 , 11.1208 , .010132 _3._17 o .5067 E2 , 2 .6_ E2 , .1263 E3 , .100 E_ , .92 , 3 .18 E3 ,2.1159267 , 1262. / OXYGEN .................... BENDER 1970 DATA(COF(1,3), 1=l,531 / .259832 , .36811077 , 1-.1_07C678E3 , .250617_4E _' ,-.10081345E7 , .IgOT_I6_E8 , 2-._.OI3_'g66E-2 , .65172112E2 , .10962206E5 , .6972158 3-.262_2649E3 , .191378 , .29_16771 EZ , .78932076E2 , _" o1923 158E7 , .46 1082'_tE 9 ,-. 39936263E11 t-. 5668995E7 , 5 .136_6286E10 , .919771gBEI1 ,5._ , _0. , -.5150_18E1 6-.25626822E3, .26280757 ,-.36222681E-2 , .26516635E-4 , 7-.10006386E-6 t .15423155E-9 , .c)OE2 , .29165189E2 8-.5691_86E-3 t .6803263_E-5 ,-.48_56604E-_ , . 13822647E-9 , 94.72426938 t .35642E3 , .03125117 , .15E-3 , 1 1.50 , .010 1325 t 5. 083 ,I.OO3OE2 , 2 .5635C7E2 , .15478E3 , .5E+3 ,1.3 , 3 .23E3 ,88.46_226 ,739.86 / ARGON ..................... BENDER 1970 DATA (COF( 1,41, I=l, 531 / .208128 , .19825921 t I-.BIT3__IlRE2 , .17777_7E4 ,.-.8240654_tE6 , .31666098E8 , 2-.4_202671E-1 , .62161_2E2 , .11_.32_8E_ , ._779752 t 3-.196_5227E3 ,-.21572754 , .165_@141E3 ,-. 281 42112E2 , @ .8 253205gE5 ,- .gI538377E 7 ,-. [ 83_0752E13 ,-. 33858130E7 , 5 .15532886E10 ,-.67_79568E11 ,3.5 , _'0. _ .1381_648E+2 , 6-.573g(:331E+3 t -.18072036_+0 t .18350292E-2 ,-.10957549E-4 , 7 .36123787E-T , -.50397228E-IO_ .8728E2 ,2.5 , 8 .0 , .0 , .O , .0 t 9 2.757177 , .237932E3 , .208152 , .1E-3 , 11._8 , .0101325 ,4.865 , .5066E2 , 2 .8378E2 , .1537E3 , .i E4 ,1.3 , 3 .25E3 _ 70.0510g2 ,725.240@3 l 88 139.98,2.C159/ DATA ZETA /.468 gO513E-l,. 40786245E_ i, .-$0i 15154r:_ i ' 1" 2762 8636E-1_ • 224073g4E-I t. _7495F-1,. 402544E-1 , .38 43, • 268 525E-1, 2 .40786245E- i/ DATA FF/ loO15t 5_l°Op 1.02, 1.1_ 1., 2.27, 1.0, 1.0/ DATA DTRIPL/°57tI.121,1.400_I.415,1.17,1.247,.83{_,.21,1.71,.[I/ DATA SgT/l13op2_68°t86.,216.t25. t70. tb.3e60.,68. / DATA KSWT 15"1,3"2,2.1/ DATA CI FTT/7_ l.t. 11 2" 1./ DATA RHOSWT/2.5t2.2t6*2.4,2.2,2.2/ DATA AMUC/- -06, O. i .04, 2_'0., 2.. 1,3.0. / DATA BMUC/.O7t°2, O.t.2tO.,.2,4*O. / END C3MPUTE P_ESSURE P GIVEN TEMPERATURE T AN) OENSITY O. U_IITS ARE SPECIFIED BY KU. IF KR IS RETURNED OR SPECIFIED AS It P IS =OMPUTED AT SATURATION AS A FUNCTION OF T ONLY,. SUBROUTINE PRESS|KUtT,DtP,KRI COMM3N/ CON123/ DCONV(5),TCONV(SItPCONV(51 COMMON/BEND/ R tCPItCP2,CP3tCP4tCP5tCP6tCP7 _CP8tCPg,CPI3, ICPIIICPI2tCPI3tCP14t:P15tCPI6tCPIT_CPI8sCPIgtCP20 2 ,CP21tCP22, CP23tCP26 C(IqMON/CHECK S/DCH liD= H2,PCH l t PCH21, PCH3 wTCHI t,TCH2 t TCH3, DST, TST p IHSCHI, t_SC H2 COqNON/BENDIT/TS CDMHON/HELFLU/ IHEt IFL9 IHY CO_IHON / TERRDR/IROUT I ROUT=2 DETERMINE REGION TS=TC HECK(KUt KR_T ) T2=I.0 IF (IHE.EO.I) T2=TS IF (KR-1I IO, 7OtlO 10 DSfOChECK(KUp D) IF (TS-TCH2I 20,20,50 20 CALL DENS(I, TS,ZE,ZE,DSLtDSVtl) IF {DS-DSLI 30,60t40 30 IF (DS-DSVl 50t60_60 &O KR=2 GO T3 80 50 KIQ=3 GO TO 80 REGION 1 60 KR=I 70 CALL PSSS(PSI GO T3 gO REGIONS 2 AND 3 80 AI= R*TS A2 = 1 ./TS A3 = A2*A2 A4=(( (CP21*A2*CPSI*A2÷CP41*A2+CP3},_AZ4.CPZ+CPL.TS AS=( CP23_A2+CP8 I* AZ+CP T+( CPZZ_TS+CP6 )*TS A6 = cPg*TS*CPIO A7= CP11'kTS+CPI2 +CPZ4_A2 A8 = ( (CP16=A2÷CP 15)=A2+CPI4)_A3 Ag = ( (CPIg'I=A2+Cp I8)*A2_-CPITI=wA3 .._i= EX P (-CP 20_'OS* D S IT 2 ) PS={( [ ((CP13_DS+A7+Ag_WB].)*DS+Ab)_,DS,,.ASeABW_BI) w,[-)S+A4I,OS+AI)_OS go P=PS*PCONV(KU } R ET JR N END 91 C CONPUTE TENPER&TU_E T GIVEN PRESSURE P AND DENSITY O. C UNITS ARE SPECIFIED BY KU° IF KR IS RETURNED OR C SUBROUTINE TENPIKUtPtDBTt KR l CC_INON/ CON123/ DCONV(5|BTCONV(5|tPCONV(5! COqMON/BENO/ R pep ItCP2 _C P3eCP6t CP5 tCP6 tCP/_CP8 _CP9 *C PlO I [CPI tt CPI2eCPI3tCP I6t_P 15tCP[6tCP[T_CPt8tCP[9tCP20 2 tCPZ[tCP22vCP23tCP26 CORNON/HELFLU/ IHEt IFkt |HY COqMON/Ct_ECKS/DCHItDCHZtPCH ItPCH2 tPCH3 t TCHI tTCH2 ITCH3 eDST eTST, 1HSCHI, HSCH2 COMHON/BEND3/DSoA IrA 3tA4tAStA68 tATB eA8B IAI2 COqNDN IIERROR/IROUT I ROUT_' 3 EXTERNAL TSFt DTSF PS-PCHECK (KUt KR tP | C C DET ERN [NE REGION C IF (KR-I) lOl70tIO 10 DS=OCHECKIKUvDI IF (PS-PCH2) 20e20. 50 20 TS=TSS(PS ! CALL DENS(IpTS_ZEtZE_OSLeDSVtll IF |DS-DSLI 30t60_40 30 IF (DS-DSVI 50t60e60 40 KR=2 T S:,TS- 10o IF (IHE.EQoll TS=3.O GO TO 80 50 KR=3 TS=TST GO TO 80 C C REG ION 1 C 60 KR=L GO T_3 110 70 TS=TSS(PS) GO TO II0 C C REGIONS Z AND 3 C 80 A1 = DS*CS A3 =( ( ( (CPI3_'DSeCP IZI_,DSeCPIOI'_DS+CPTI*DSo'CPZ I*AI -PS A4 = ( ( ( ( Cp I I_DS+CP q) _D S+CP 6) _3 S+CP/)*D S_RI *Z) S A5=(( CPE_*A leCP8 )*D S_CP 3 I*A 1 A6fi= (CP23tDS÷CP4|*A1 A7_=CP 5*A 1 A8B=CP 21_'AI A12 = 2**CP22_A I*DS TS=SOLVE( TSm TSF tDTSF | C C VER IFY PEGION C ]F (PS-PCH2) 1 lOt l 10190 90 IF (TS-TCH2I lOOt lOOtllO IO0 KR=2 110 T =T S,(,TCONV( KU I RETURN END 92 C C C COMPUTE DENSITY O GIVEN TEMPERATURE T ANO PRESSURE P. UNITS ARE SPECIFIED BY KU. IF KR IS RETURNED OR SPECIFIED AS It THE SATURATED LIQUID AND VAPOR DENSITIESt 01. AND OV RESPECTIVELYt ARE COMPUTED AS A FUNCTION OF T OR P. THE OTHER VALUE MUST BE INPUT AS 0o0 ° 7O 75 SUBROUTINE DENS(K Ut T_ P _D t DL tD Vl, KR) COMM3N/ CDNI23/ OCONVIS|tTCONVI5)tPCONV(51 COMMON /CHECK I/NI COqMON / BEND/ R tCPltCP2_C P3tCP6tCP5 tCP6 tCP7 tCP8 tCP9 ICPI3 o |CPllt CP12tCPI3_CPI4_PIS_CPI6_CP17tCP18tCPIgtCP20 2 _CPZIICP22t CP23_CP 26 CO_MON IBENO !.31A lwA6pA 5,A6tATtA8_Agt PS COMMON/BENDISIAIOeAII_AI2tAI3tAI6 COWN]N/BENDIT/TS COWMON/CHECKS/OCHItDCHZ_PCHLtPCHZt PCH3_ TCHI _TCH2 tTCFS,{)ST,TST, 1HSCH1, HSCHZ C0_ NON/CO FMU/EP SK pWMt D I S ! RHOC R T t ZE TAA CO_MON/FELFLU/ IHE. IFLt, IHY COMMON/SOL ID/DTR IP COWMON / IERROR / IROUT EXT_NAL DSFz, DDSF I ROUT=I IF (KR.EOol) GO TO 70 T S=TCI_ ECK (KUt KRt T | GO T_ 5 IF (T°GT°O°O} GO TO 75 PS=PCHECK(KUt KR tP | TS=TSS(PSI IF (T.LE.O°I T=TS*TCONVIKU) GO T] 5 TS=TChECK(KUtKRtT| CALL PSSS(PS) IF iP..LE.O..! P=PS*PCONVIKU) T2=I.O IF IIHEoEO,I| T2=TS COMPUTATIONS COMMON TO ALL REGIONS At= R*TS A2= I.ITS A3= A2*A2 A4=(( (CP21*AZ_CP51*AZeCP6I*A2+CP31*A2eCPZeCPI*TS A5= ( CP23*A2_CP8 I_A2*CP 7+[CP22*TS_CP61*TS A6= cPg*TS*CP 10 AT= CPII_TSt-CPI2 _CP24_A2 A8=( ( CP 16*AZ_-CP 15 _*A Z ¢'CP 16l _A3 A9=( ( CP xg_vA 2*CP 18 )_A Z _'CP £ 7l *A3 A I0=-2 .*CP 20/T2 AIL=6**CP 13 A12=5._A7 AI3= 4.*A6 A 16= 2 .*A4 DETERMINE REGION IF IKR-II 10t80,10 lO PS=PCHECKIKU_KRtP) IF (PS-PCH21 llO,110, L30 I00 IF (TS-TCH2|I20,12CtI30 120 KR=2 EST=R HOCPT*.9O CALL ROOT(DCH2tEST_O. tDSF_OS) GO TO 150 130 KR=3 EST=R HOCRT*3. CALL ROOT( ESTt DCHlt 0° _DSF _D S) GO T3 150 110 IF (TS-TCH21 20_20t50 20 CaLL PSSS(PSS| IF (ABS((PSS-PS)/PSSI-I.E-4| 60_30,30 30 IF {PS-PSS| 50t6.qt60 93