A FORTRAN PSLCHRCMLL Ibesisfor Liza? Be we; a? LL LL :TTTC’HT GAL! STAT E UTLH’ERSTTT LLOYD EUGELLE 2.ng ...- L972 "LIIIIILIHZLILILJIIIIILIIIJLLILI111MIILHISIIIILILI, LIB RARY " ['Mich‘gan Stat: L 1‘. umvchity I} 151251. 3' ABSTRACT A FORTRAN PSYCHROMETRIC MODEL BY LLOYD EUGENE LEREW Psychrometric properties are basic to many agricultural systems. A psychrometric model is therefore required for computer simulation and analysis of such systems. A collection of theoretical and empirical psychrometric equations have been programmed as a set of efficient and interconnected FORTRAN sub-programs. Dry-bulb, wet-bulb and dew-point temperature, absolute and relative humidity, vapor and saturated vapor pressure, latent heat, enthalpy, and specific volume are included in the model. Given any two of these prop- erties (if they are independent), it is possible to find the remaining properties with the model. Specific instructions for the use of the model are given for a num- ber of pairs of independent starting properties. Major Professor /ég;/;;;/zflg7fi,.. Approved Department Chairman A FORTRAN PSYCHROMETRIC MODEL BY LLOYD EUGENE-LEREW A THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of MASTER.0F SCIENCE Department of Agricultural Engineering 1972 464T 0‘ ACKNWLEDGMENTS The author gratefully acknowledges the guidance, trust and encour- agement provided by his major professor, Dr. F. W. Bakker-Arkema. The financial support provided by the Craps Research Division of the Agricultural Research Service of the United States Department of Agriculture is greatly appreciated. ii PREFACE The thesis, A Fortran Psychrometric Model, written by Lloyd E. Lerew represents only a small part of his research efforts over the past two years. Following is a complete list of the projects to which Lloyd E. Lerew has either contributed significantly or has been its sole re- searcher: 1. Grain Drying Bakker-Arkema, F. W., T. W. Evans, and L. E. Lerew (1970) ‘Michigan State Grain Drying Models. ASAE Paper No. 70-832. Bakker-Arkema, F. W. and L. E. Lerew (1970) Single particle cooling and thin-layer drying simulation. In Proceedings of the Institute for Simulation of Cooling and Drying Beds of Agricultural Products, Agricultural Engineering Department, Michigan State University, East Lansing, Michigan. Bakker-Arkema, F. W., R. J. Patterson, and L. E. Lerew (1971) Multiple-zone drying in stationary and moving bed dryers. ASAE Paper No. 71-302. Bakker-Arkema, F. W., L. E. Lerew, and R. J. Anderson (1972) Better quality grain through the use of concurrent flow grain drying with counterflow cooling. In Proceeding of Grain Damage Symposium, Ohio State University, Columbus, Ohio. Bakker-Arkema, F. W., L. E. Lerew, and S. F. DeBoer (1972) Grain dryer simulation. Research Report - Michigan Agricultural Experiment Station, Michigan State University, East Lansing, Michigan. Bakker-Arkema, F. W., M. G. Roth, and L. E. Lerew (1972). Gleichzeitig Warme-und Stoff-Austausch in einem getreidetrockner mit gegenstrom. Landtechnische Forschung (accepted for publication). Brooker, D. B., L. E. Lerew, and F. W. Bakker-Arkema (1971) Drying by alternating heated and unheated air: a mathematical analysis. ASAE Paper No. 71-307. iii Lerew, L. E. (1970) Comparison of thin-layer drying equations for corn. Unpublished Report. 2. Potatoes Lerew, L. E. (1972) Energy and time requirements for heating potatoes after storage. In Proceedings of the Midwest Potato Storage Workshop, Cooperative Extension Service, Michigan State University, East Lansing, Michigan. Rosenau, J. R., F. W. Bakker-Arkema, and L. E. Lerew (1970) Cooling of a bed of potatoes (or apples, tomatoes, sugarbeets, onions): a realistic model. In Proceedings of the Institute for Simulation of Cooling and Drying Beds of Agricultural Products, Agricultural Engineering Department, Michigan State University, East Lansing, Michigan. 3. Other Bakker-Arkema, F. W. and L. E. Lerew (1972) Heat exchangers, computer analysis. IFT Convention, Minneapolis, Minnesota. Lerew, L. E. (1972) Toward more effective use of the CDC 6500 (how to get the 6500 before it gets you). To be published, Department of Agricultural Engineering, Michigan State University, East Lansing, Michigan. Lerew, L. E., F. W. Bakker-Arkema, and G. H. Bulgrien (1972) User's guide to SGRAF6 and SPLOTB: printer and Calcomp plot packages. Unpublished report. It was the intention of the guidance committee to include in this thesis all the papers of which Lloyd E. Lerew had been the author or a co-author. Unfortunately current administration procedures within the Colleges of Engineering and of Agriculture and Natural Resources at Michigan State University did not allow this. F. W. Bakker-Arkema Chairman Guidance Committee of L. E. Lerew and Major Professor. iv TABLE OF CONTENTS INTRODUCTION . . . . . . . . . . . . . . . . BASIC PSYCHROMETRIC EQUATIONS . . . . FORMULATION OF THE MODEL . . . . . . . . USE OF THE MODEL . . . . . . . . . . . . . . ZEROIN . . . USER CONSIDERATIONS . . . . . . . . . . . MODEL VALIDATION . . . . . . . . . . . . . RECOMMENDATIONS FOR FURTHER STUDY . . SUMMARY . BIBLIOGRAPHY . . . . . . . . . . . . . . . . APPENDIX A: APPENDIX B: APPENDIX C: APPENDIX D: APPENDIX E: APPENDIX F: Original References . . . . . . Model Listing . . . . . . . . . Suggested Calling Sequence With 0 O O 0 DB Known Suggested Calling Sequence To Find DB . Validation Program . . . . . . Results Of Validation Program . Page 15 17 20 21 22 23 24 27 28 34 37 38 41 LIST OF TABLES Page Mnemonic Codes Used in the FORTRAN Psychrometric Model . . . . . 10 Psychrometric Properties by Independent Groups . . . . . . . . . 10 Relationships Given by the Equations . . . . . . . . . . . . . . 11 Pairs of Properties Which Require TWo-Dimensional Search . . . . 11 Pairs of Beginning Properties for the Psychrometric Model . . . 12 Function Subprograms in the Model . . . . . . . . . . . . . . . 13 vi LIST OF SYMBOLS Dry-bulb temperature, R Wet-bulb temperature, R Dew-point temperature, R Vapor pressure, psi Saturation vapor pressure at T, psi Saturation vapor pressure at wa, psi Humidity ratio, 1b water per 1b dry air Humidity ratio at saturation at Tw lb water per lb dry air 13’ Latent heat of vaporization of water at saturation, Btu per 1b Latent heat of vaporization of water at Btu per 1b wa, Latent heat of vaporization of water at po, Btu per 1b Heat of sublimation of ice, Btu per lb Heat of sublimation of ice at T Btu per 1b wb’ Heat of sublimation of ice at po, Btu per 1b Enthalpy of air-vapor mixture, Btu per lb dry air Relative humidity, decimal 3 Air specific volume, ft per 1b dry air vii INTRODUCTION Psychrometrics is basic to research and design of many agricultural processing, storage and environmental systems applications. A number of professional organizations have publiched psychrometric charts and tables to be used as standards for their industries (ASHRAE, 1963; ASAE, 1971). Charts and tables are not well suited to the application of computer aided design and simulation. The recent increase in the use of computer techniques has created the need for a standard computer model of the psychrometric chart. The psychrometric model presented here is based on a combination of theoretical and empirical equations. Ten thermodynamic prOperties are considered. In all but five cases, it is possible to find the eight remaining prOperties when given two independent prOperties which uniquely describe a state point of a dry air--water vapor mixture. The model has been designed to be easily implemented as well as easily modified. Several novice programmers have already successfully applied portions of the model (Farmer, 1971; Schluckbier, 1972). Equations for additional pr0perties or more accurate equations can be easily added. Portions of the model have been used in the simulation of grain dryers (Bakker-Arkema et a1., 1972). BASIC PSYCHROMETRIC EQUATIONS Brooker (1967, 1970) has compiled and tested a number of equations for psychrometric prOperties. Appendix A contains the original references given by Brooker for the equations used here. All error limits pre- sested here are taken from Brooker (1970). The following assumptions were made by Brooker (1970): (1) The properties of steam are prOperly described by Kennan and Keyes (1936). (2) Air-water vapor mixtures can be treated as a mixture of ideal gases. The saturation line can be described by 2 . 11 8: 6489 - 0.46057 Ln T (1) LnP8 = 23.3924 - 459.69 S T 5 491.69 2 3 4 and by Ln“, /R) _ A + BT + CT +2131: + ET (2) 5 FT - GT 491.69 S T 5 959.69 where R - 0.3206182232000000 D 04 D - 0.2153211916363544 D -04 A - -0.27405525836l4256 D 05 E - -0.4620266568199822 D -08 B = 0.5418960763289505 D 02 F a 0.2416127209874000 D 01 C - -0.4513703841126545 D -01 G - 0.1215465167060546 D -02 The form of equation (1) results from the integration of the Clausius-Clapeyron equation. Coefficients were obtained by fitting the equation to the data of Keenan and Keyes for a temperature range of 0 to 32 F. J range of -20 to 32 F is as follows: The deviation from Keenan and Keyes for a temperature Alll'ilT.lT|ill‘l.lT.l[ 0.00002 S P S 0.00069 error 0.2993 S % error S 0.8501 Equation (2) is empirical. The deviation from Keenan and Keyes data is shown below for a few data points where the error was large when compared to the error at other temperatures. At 32 F Absolute error = 0.0002 psia Percent error 8 0.228 At 40 F Absolute error = 0.0001 psia Percent error 8 0.081 At 470 F Absolute error - 0.5027 psia Percent error - 0.097 Latent heats are found by the following empirical equations: h18 - 1220.844 - 0.05077 (T - 459.69) (3) 459.69 S T S 491.69 hf8 I 1075.8965 - 0.56983 (T - 459.69) (4) 491.69 S T S 609.69 2 1/2 hfg = (1354673.214 - 0.9125275587 T ) (5) 609.69 S T S 959.69 Equation (3) can be used for temperatures as low as -20 F with an absolute error in his of 0.66 Btu per 1b. For a temperature range of 0 to 32 F, the error range is as follows: 0.0209 S h 5 0.1932 error 0.0017 S % error S 0.0158 Error limits for equations (4) and (5) from 32 to 480 F are as follows: 0.00 S h S 1.47 error 0.00 S % error S 0.20 From 480 to 500 F, the absolute error increases to 3.2 Btu/lb. At 500 F, the percent error is 0.45. The wet-bulb line is described by Pswb - Pv - B (wa - T) (6) where B a 0.2405(Pswb - Patm) (1 + 0.15577 Pv/Patm) 0.62194 h' fs 459.69 S T S 959.69 Equation (6) was derived on the assumption that(H' - H) 1b of water vapor is combined with one 1b of dry air and H 1b of water vapor is combined with one 1b of dry air and H lb or water vapor by an adiabatic process at wet bulb temperature. The derivation is shown by Brooker (1967). For temperatures below 32 F values of h' are used in equation is (6) rather then h'fg' There is no way to express the accuracy of equation (6) except to use the equation and then check the answer against some known standard. Examples are given by Brooker (1967). The equation for absolute humidity is 0.6219 P v P - P atm v H _ (7) 459.69 S T S 959.69 This equation is based directly on the definition of H and assumes that dry air and water vapor are ideal gases. Specific volume is found by the equation a 53.35 T (8) as 144 (P - P ) a V V tm 459.69 S T S 959.69 Equation (8) is the ideal gas law and needs no explanation. Enthalpy is best described as Enthalpy - Enthalpy of air + Enthalpy of water (or ice) at dew point temperature + Enthalpy of evaporation (or sublimation) at dew-point temperature + Enthalpy added to the water vapor (superheat) after evaporization. The results in the equations h - 0.2405 (Tdb - 459.69) -H [143.35 + 0.485(491.69 po)] + h 18 H + 0.448 H ('1‘db po) 459.69 s po s 491.69 (9) and h - 0.2405 (Tdb - 459.69) + H (po - 491.69) (10) I I - +h fg H + 0.448H (Tdb po) 491.69 s T s 671.69 dp The value of specific heat for water vapor is considered to be 0.448 Btu per lb per F in equations (9) and (10). This value is apprOp- riate for dew-point temperatures of 70 F or less. 'If for some reason the dew-point temperature is higher, the specific heat of the vapor should be considered a function of temperature and increased. Relative humidity is defined as rh - 53. (11) Equations (1) through (11) constitute the complete set on which the psychrometric model is based. FORMULATION OF THE MODEL In order to simplify the notation, the mnemonic codes listed in Table 1 were develOped. A number of the ten prOperties shown in Table l are not independent. To uniquely describe a point on the psychrometric chart, two independent prOperties must be given. Table 2 lists dependent prOperties together by groups; two pr0perties not in the same group in Table 2 do describe a point uniquely. A complete model of the psychrometric chart should be capable of finding the remaining eight prOperties regardless of which pair of independent prOperties are given. If relationships are known such that given any prOperty within a group all other prOperties within that group can be found, all that remains to completely model the chart is to specify the relationships among groups. Equations (1) through (5) represent relationships between DB and the other prOperties of group I. Equation (6) is a relationshipimmong grOUps I, II, and III (DB, PV, and WB respectively). A relationship between members of group II is given by equation (7). Equation (8) deals with groups I, II, and V, equations (9) and (10) with I, II, and III, and equation (11) with groups I, II, and IV. This is summarized in Table 3. The form of the equations makes explicit solution for some of the prOperties impossible or quite difficult. Equation (2), for example, is difficult to solve for DB, the dy bulb temperature,and equation (6) can not be solved explicity for WB. The use of a search technique helps to overcome these difficulties. .l’llli'll‘flll '1' ll 1‘ l Tu T l [1‘ .TI Careful examination of equations(1) through (11) reveals that relationships between groups III and IV, groups III and V, or groups IV and V with any other group are non-existent. This means that, given pr0perties from any two of these groups, the relationships to explicitly solve for the remaining prOperties are missing. The equations which express relationships among groups are all keyed to groups I and II (see Table 3). With the aid of a two-dimensional search, it is possible to find group I and group II prOperties when only group III, IV, or V prOperties are given. Table 4 lists the possible pairs of prOperties which would require a two-dimensional search. ! I [III [.6]- TABLE 1. Mnemonic Codes Used in the FORTRAN Psychrometric Model fl . Code DB DP EN PS PV RH VS WB Property dry bulb temperature dew point temperature enthalpy absolute humidity latent heat of vaporization (sublimation) saturation vapor pressure vapor pressure relative humidity search or special specific volume of dry air wet bulb temperature 921.29. R R Btu/lb dry air lb water/lb dry air Btu/1b water psia psia decimal ft3/1b R TABLE 2. Psychrometric PrOperties by Independent Groups Group I II III IV PrOperties DB, HL, PS DP, HA, PV EN, ws RH VS TABLE 3. Relationships Given by the Equations Equation Number 1 2 10 11 Relationship I I I I I I, II, III II I, II, v I, II, III I, II, III I, 11, IV TABLE 4. Pairs of Pr0perties which Require a Two-Dimensional Search EN E RH and RH and RH and VS and V8 and VS TABLE 5. Pairs of Beginning PrOperties for the Psychrometric Model DB DP EN E PS PV RH VS WB P - Possible to find remaining pr0perties starting with these two . N - Two dimensional search necessary to find remaining properties starting with these two. I - Infeasible; this pair does not uniquely describe a point on the psychrometric chart DB DP EN HA HL PS PV RH VS WB P P P I I P P P P P I P P I P P P P P P P N N I P P I P P P I P P P P P P P P P P P N N N 10 TABLE 6. Name DBPSS# DPDBENS# DPHAS# ENDBDP HADBRH HADP HAPV HLDB Function Subprograms in the Model Statement 0 - PS - PSDB (DB) .2405*(DB-459.69)-EN+HA*HLDB(DP)+ .448*HA*(DB-DP)-HA*(143.35+.485* (491.69-DP)) for DP < 491.69 .2405*(DB-459.69)-EN+HA*HLDB(DP)+ .448*HA*(DB-DP)+HA*(DP-49l.69) for DP 2 491.69 0 - PVHA(HA)-PSDB(DP) - [.2405*(DB-459.69)+HLDB(DP)*HA+.448* HA*(DB-DP)]-[HA*(143.35+.485*(49l.69 -DP))] for DP < 491.69 - E.2405*(DB-459.69)+HLDB(DP)*HA +.448*HA*(DB-DP)]+HA*(DP-49l.69) for DP 2 491.69 - HAPV(RH*PSDB(DB)) - HAPV(PSDB(DP)) - .6219*PV/(PATM-PV) = 1220.844 -..05077*(DB-459.69 for DB < 491.69 = 1075.8965-.56983*(DB-459.69) for 491.69 S DB < 609.69 11 Equation Number 1 10 10 11 11 12 - SQRT(1354673.214-0.912575587*DBZ) for 03 2 609.69 5 PSDB - EXP[23.3924-11286.6489/DB-.46057* ALOG(DB)] for DB<49l.69 1 a R*Exp[(A+B*DB+C*DB2+D*DB3+E*DBT)/ (F*DB-G*DBZ)] where R,A,B,C,D,E,F,G are constants (Brooker, 1970) 2 PVDBVS I PATM-S3.35*DB/VS/144 8 onsws - TPSDB(WB)*[.62194*HLDB(WB)*PATM] 6 -[.2405*(PSDB(WB)~PATM)*(WB-DB)] *PATMT/TE.62194*HLDB(WB)*PATMJ +{.15577*[.2405*(PSDB(WB)-PATM*(WB-DB])i PVHA a HP*PATM/(.6219+HA) 7 DBHA - PVHA(HA)/PSDB(DB) 11 RHPSPV - PV/PS 11 VSDBHA - 53.35*DB*(.6219+HA)/144./.6219/PATM 8 WBDBHAS# o - TWB-DB-(PSDB(WB)-PVHA(HA))/[.2405* (PSDB(TWB)-PATM)*(1.+.15577*PVHA(HA)/ PATM)]*(.62194*HLDB(WB)) 6 # These equations are solved via SUBROUTINE ZEROIN. l(1.. v'llll'll‘ll‘. 13 Since the occurrence of a situation which would provide only one of the pairs of preperties listed in Table 4 is considered to be very unlikely, a two-dimensional search is not included with this model. A number of such routines are available (Rosenbrock, 1966; Powell, 1964) and can be implemented easily if the need arises. Table 5 depicts fourty-five combinations of the ten properties taken two at a time, and shows which are infeasible due to dependence, which would require the use of a two dimensional search, and which are possible with this model. Note that of the thirty-eight physically possible starting pairs, thirty-three can be used without the aid of a two-dimensional search. The model contains nineteen function subprograms and one subroutine (a zero finding technique discussed later) and required less than 1300 octsl memory positions. Several of the function subprograms are used internally only and need not concern the user. Those with which the user is concerned have been named using the codes of Table l and a mnemonic scheme. The first two letters are the code of the prOperty that is sought; the second two and third two (if present) are the codes of the prOperties of which the first is a function. If the first is a function of two values, they are listed in alphabetical order. Examples: RHDBHA - relative humidity as a function of dry bulb temp and absolute humidity (humidity ratio). PVHA - vapor pressure as a function of humidity ratio. If the prOperty which is being sought is a temperature (dry-bulb, wet-bulb, or dew-point) the equations require the search routine for a solution;and the user will have to supply additional parameters l4 (discussed later). An "S" is added to the names of these routines to help remind the user to supply the additional values. Table 6 lists the sixteen function names which are accessible to the user and the principle statement along with the number of the equation_ from which it was develOped. Most of the develOpments are straight- forward. The zero finding search technique (subroutine ZEROIN) is used to solve equations of difficult form or, as in the case in DPHAS, to solve for a prOperty common to two equations. The computer programs are listed in Appendix B. USE OF THE MODEL The sixteen functions in Table 6 plus the search technique enable the user to find the eight remaining prOperties given DB and any other independent prOperty. In addition, it is possible to find DB given one prOperty from group II and one property from groups III, IV or V. Appendices C and D contain lists of suggested function calls and the order of calls and the order of calls to accomplish this. Use of the non-search functions (those without an 8 attached) is simple and best illustrated by examples: 1) If TABS - absolute dry bulb temperature and R:- relative humidity (decima1)* and it is desired to find the humidity ratio H, the following FORTRAN statement would be used a - HADBRH(TABS,R)) ii) If the vapor pressure VAP is desired, the following might be used VAP - PVHA(HADBRH(TABS,R)) This is simply a nested call similar to x - EXP(SIN(Y)) for esin(y) Use of the search functions (those with an 8 appended) X - s is only slightly more difficult. The general form is: FNAME(PAR1,PAR2,GUESSl,GUESS2,EPSLON) where FNAME - the function name (DBPS, WBDBHAS, etc.) PAR1,PAR2 - first and second (if applicable) parameters specified by the two letter codes in the function name. GUESSl - a low estimate of the answer GUESSZ - a high estimate of the answer EPSLON I allowable tolerance in the result The values of GUESSI, GUESSZ, and EPSLON should be chosen with care as they greatly affect the time required for the zero-finding search * Note that all temperatures are in degrees Rankine and that the relative humidity is in decimal form. 15 16 technique to arrive at a solution. In general, it should be remembered that input psychrometric quantities are seldom accurate to more than 2 or 3 decimal places and hence an EPSLON of .01 or .001 should be sufficient. The search routine guarantees convergence (usually better than linear, Dekker, 1967), but small EPSLONs require more time. If' the values of GUESSl and GUESSZ are chosen such that both are above or both are below the correct answer, an apprOpriate error message is printed but computation continues. Example: If the following are known TDB I dry bulb temperature, R TDP I dew point temperature, R and the wet bulb temperature is needed, Appendix C says that the absolute humidity must first be found from HADP and then the wet bulb can be found from WBDBHAS. Therefore the calls could take the form ' HAVS=HADP(TDP) . ' TWB-WBDBHAS(TDB,HABS,TDP,TDB,.01) or TWB-WBDBHAS(TDB,HADP(TDP,HADP(TDP),TDP,TDB,.01) if the absolute humidity is not desired separately. Note the use of TDP and TDB as GUESSl and GUESS2. ZEROIN ZEROIN is the root finding search technique included in the psychrometric mode. Although used by some of the functions in the model, this subroutine remains completely independent of them and can thus be freely used for other purposes. ZEROIN is based on an algorithm presented by Dekker (1967) which is a mixture of linear interpolation, extrapolation and bisection. Convergence to a small interval in which a zero exists is guaranteed. No derivates are used. ZEROIN must be given a function for which a zero exists on an interval J, end points A and B of this interval and EPS the largest acceptable final interval. ZEROIN then works with three points A, B, and C on the interval and the function values at these points FA, EB, and FC. Point B is the most recent ilerate, point A the previous one and point C the "contrapoint" or the last one satisfying SIGN (FB)# SIGN (FC). Initially, point B is taken as the end point of J which results in the smallest absolute value of the function; points A and C are taken as the Opposite end point. Interpolation or extrapolation is performed between A and B yielding a provisional value I for the next iterate according to the formula: I-J—l—‘B'A FB +3 FB-FA If this division results in a divide fault or overflow, bisection is used to determine I by the equation: B+C I M 2 l7 18 Computation continues using the newly calculated value. If no problem exists with the division, I and'M are checked to determine which is closer to B. I is then assigned this value. Having determined a suitable value for I, the difference between B and I is compared to EPS. If EPS is larger, I is re-evaluated as: I I SIGN (C-B) x EPS+B If EPS is smaller, the following reassignments are made in preparation for the next step: A-B B=I FAHFB FB-the function evaluated at I. One further check is necessary before beginning the following step. If SIGN FB I SIGN FC, C is reassigned to A and EC is reassigned the value of FA. Whenl C - BI S 2 x EPS the termination procedure is initiated. ZEROIN takes the mid-point value between B and C, assigns it to A and evaluates the function at this point, FA. The final values re- turned are A and B or C (whichever one has a function evaluation different in sign from FA). A call to ZEROIN takes the form: CALL ZEROIN (A,B,EPS,FUNC) where A and B a the end points of the interval containing a zero EPS - the largest acceptable final interval containing a zero FUNC - the name of the user supplied function subprogram containing the function for which a zero is to be found. Note that FUNC must be typed EXTERNAL in the calling routine. ZEROIN returns in A and B the endpoints of the final interval of length less than EPS. Either A or B or their average can be used as the zero. 19 ZEROIN checks the function values of the initial interval end points. If the function values have the same sign, ZEROIN sets the value of ICHK in the labeled common area ICHECK to one. Calculations Calculations then continue as this indicates that multiple roots or no roots occur within the given interval. In the case of multiple roots one of them will be found. In the case of no roots ZEROIN will converge to the point where the absolute value of the function is a minimum, whether this occurs at an end point or arbitrary value within the interval. Subroutine ZEROIN can easily be used to form additional functions from the sixteen given. For example, suppose that dew-point temperatures DP and wet-bulb temperatures WE are known in some cases and relative humidity RH and DP are known in other cases and that in both cases the dry-bulb temperature DB is to be found. Appendix D shows that in the first case absolute humidity and then vapor pressure should be found so that PVDBWB can be used with ZEROIN to solve for DB. This might be programmed as FUNCTION DBl(T) COMMON/BLOCK/PV,WB DBl I PV-PVDBWB(T,WB) RETURN END where T is a dummy variable representing the various values for DB tried by ZEROIN. Appendix D also shows that in the second case the absolute humidity should be found and that RHDBHA should be solved for DB using ZEROIN. This might be as follows: FUNCTION DBZ(T) COMMON/BLOCKZ/RH,HA D82 I RH-RHDBHA(T,HA) RETURN END 20 The calling routine might be PROGRAM SAMPLE(INPUT,etc.) COMMON/PRESS/PATM COMMON/BLOCKl/PV,WB COMMON/BLOCKZ/RH,HA EXTERNAL DB1 EXTERNAL 032 statements to assign values to PATM, EPS ESTl, EST2 and determine whether DP and WB or DP and RH are given when DP and WE are given CALL ZEROIN(EST1,EST2,EPS,DB1) DB I (EST1+EST2)/2. when DP and RH are given CALL ZEROIN(EST1,EST2,EPS,DB2) DBI(ESTl+EST2)/2. END Note the labeled common area PRESS. This must be specified in all programs which use the model and a value for atmospheric pressure in psia must be assigned to PATM. USER CONSIDERATIONS The following should help the user avoid possible errors and locate problems that do arise: 1. These routines are said to give accurate results from 440 to 960 R (-20 to 500 F). Accuracy within this region is generally better than can be read from a psychrometric chart. The values returned from these functions for a temperature of 491.69 R (32 F) are those for water vapor mixtures at that temperature. Lower temperatures will cause the heat of fusion to be included where apprOpriate. A11 temperatures must be in degrees Rankine. A11 relative humidities must be decimal. A MODE 0 error (on CDC 6000 series computers) indicates that the EXTERNAL card has been omitted. The following names are used by this package (in addition to the functions listed in Table 6) and should not be used in the user's program DPLl DPL2 DBL WBL SPECIAL The following statement must appear in the user's program COMMON/PRESS/PATM and a value must be assigned to PATM. 21 MODEL VALIDATION No attempt was made to check the error limits reported by Brooker (1970) for individual equations. Instead, a validation program was written (Appendix E) which computes all remaining properties given DB and HA. These computed prOperties were then used, one at a time, with the given DB to recompute all remaining prOperties. The resulting output is in table form for easy comparison of all computed values. Appendix F contains a number of results from the validation program. The values of DB and HA were chosen so that comparison with the psychro- metric chart would be easy. Only a few points show variation to the accuracy printed and this is never more than one or two in the least significant digit. Comparison of the computed values with the chart shows the model to be at least as accurate as the chart. 22 RECOMMENDATIONS FOR FURTHER STUDY Although the inclusion of a two-dimensional search enables this model to completely describe the psychrometric chart, it is not recommended as a standard. Differences in computer hardware and soft- ware preclude this. Ideally, a model of the psychrometric chart should need neither a one-dimensional nor two-dimensional search. Work should be done to derive theoretical or empirical relationships among more of the properties. A true standard would need to specify accuracy and number of significant digits to be carried through the calculations. In addition, the effect of changes in atmospheric pressure may not be properly reflected by this model. More experimental data at various pressures is needed. 23 SUMMARY This package of subprograms was written to facilitate the use of psychrometric parameters and the processing of psychrometric data. The included Tables and examples are intended to assist the user in the application of these routines and are not to be considered the only possibilities. Much thought was, however, was put into the lists of suggested calls and calling orders and it is believed that they are the most efficient. The routines have been thoroughly tested. They gave excellent results when compared to the thermodynamic values read from the standard psychrometric charts. 24 B IBL IOGRAPHY BIBLIOGRAPHY ASAE (1971) Agricultural Engineers Yearbook. American Society of Agricultural Engineers. St. Joseph, Michigan ASHRAE (1963) ASHRAE Guide and Data Book: Fundamentals and Equipment. American Society of Heating, Refrigerating, and Air Conditioning Engineers. New York, New York. Bakker-Arkema, F. W., Lerew, L. E., and DeBoer, S. F. (1972) Grain dryer simulation. Research Report - Michigan Agricultural Experiment Station, Michigan State University, East Lansing, Michigan. Brooker, D. B. (1967) Mathematical model of the psychrometric chart. Transactions of the ASAE (10) 558. Brooker, D. B. (1970) Modeling of the psychrometric chart. Proceedings of the Institute for Simulation of Cooling and Drying Beds of Agricultural Products, Michigan State University, Nov. 2-4, 1970. Brunt, David (1941) Physical and dynamical meteorology. Cambridge University Press, pp. 85-87. Dekker, T. J. (1967) Finding a zero by means of successive linear interpolation. In Conference of Constructive Aspects of the Fundamental Theorem in Algebra, Zurich, Switzerland. Wiley and Sons, Inc., New York, New York. Farmer, D. M. (1971) Personal communication. Keenan, J. H. and F. G. Keyes (1936). Thermodynamic prOperties of steam. John Wiley and Sons, Inc. Powell, M. J. D. (1964) An efficient method for finding the minimum of a function of several variables without calculating derivatives. Computer Journal 7(2)155. Rosenbrock, H. H. and C. Story (1966). Computational Techniques for Chemical Engineers. Pergamon Press. Schluckbier, G. W. (1972) Personal communication. Smith, L. B., Keyes, F. G. and H. T. Gerry (1934) Proc. American Academy of Arts and Science. (69)l37. 26 APPENDICES APPENDIX A Original References References for equations given by Brooker (1970). Equation Number Reference 1 Brooker (1967) 2 Smith, et a1. (1934) 3 Brooker (1967) 4 .L219- 5 Brooker (1970) 6 Brunt (1941) 7 Brooker (1970) 8 $239: 9 Ibid. 10 Ibid. 11 27 l‘tl .ll‘ll‘ lr. -{Tls APPENDIX B: Model Listing C""‘ P S Y C H R O H E I R I C N O O E L C”"’ L.E. LEREN, PROGRAMMER ' COOOO. COOOOO C""' DESCRIPTION Cfffff GROUP OF HIGHLY INTERACTIVE FUNCTION SUBPROGRAMS BASED C””* ON EQUATIONS DY BROOKER (1967,1970) AND A ROOT FINDING SUBI Cfffff ROUTINE BASED ON AN ALGORITHM BY DEKKER (1967) USED TO FIND C""‘ AND CONVERT PSYCHROMETRIC PROPERTIES COJOO. C"“‘ FUNCTION SUBPROGRAMS CONTAINED C““’ DBPSS C“"’ DBL (USED INTERNALLY ONLY) C""' DPDBENS C’tttf DPHAS C""‘ DPLI (USED INTERNALLY ONLY) C”“' DPLZ (USED INTERNALLY ONLY) Cffttf ENDBJP thfff HADBRH COOOOO “ADP C‘OO‘. HAPV COCO.‘ “L08 0900.. P508 COOOOO PVHA C"“‘ PVOBVS C"“‘ PVOJHB 0"‘2'. RHPs=v (ENTRY-IRHDBHA) C""’ VSDBiA C““’ NBOBiAS C"“‘ HJL (USED INTERNALLY ONLY) COO... C""‘ SUBROUTINES CONTAINED Cfttff ZEROIN COOOOO FUNCTION DBPSS(PS,GI,GZ,EPS) COMMON /SPECIAL/PV,TB,XTRA COMMON/ICHECK/ICHK EXTERNAL DBL A=Gl 3 8:62 5 XTRAIPS CALL ZEROIN (A,B,EPS,DBL) UBFSS=(A+B)/Zo IF(ICHK.EQ.1) PRINT ZIJ,GI,GZ,DBPSS 210 FORMAT(1H0,130(iHS)/10X54HN A R N I N G-INO 0R MULTIPLE ROOTS DETE ICTED BY ZEROIN/30X33HCALLEO FROM DBFSS O ZNUPGZd.14/30X17HRETURNED VALUE 159621.14) RETURN END FUNCTION DBL(T) COMMON lSPECIAL/PV,TB,PS OBL=PSIPSDB(T) RETURN END GUESSES ARE,GZI.14,4M A sssss sssss sssss sssss sssss sssss sssss sssss sssss sssss sssss sssss sssss sssss sssss sssss sssss sssss sssss sssss sssss sssss sssss sssss sssss sssss sssss sssss sssss sssss sssss sssss sssss sssss 29 FUNCTION DPDJENS(OB.EN.GI,GZ,EPS) COMMON /SPECIAL/PV,TB,XTRA COMMON/ICHECK/ICHK EXTERNAL DPLZ A851 S B=GZ XTRA=.2h05'(DB-459.69)-EN ~TB=DB CALL ZEROIN (A,B,EPS,DPLZ) DPDBENS=(A+B)IZ. IF(ICHK.EQ.1) PRINT 210,61,GZ,DPDBENS 210 FORMAT(1HO,13J(1H3)/10XSHHH A R N I N G--NO OR MULTIPLE ROOTS DETE 10TED BY ZEROIN/3OX33HCALLED FROM DPUBENS. GUESSES ARE,GZi.1h,hH A 2N0’62201“/3Ux1IHRETURNED VALUE IS,521o1h) RETURN END FUNCTION OPHAS(HA,GI,GZ,EPS) COMMON lSPECIAL/PV.OB,XTRA COHHON/ICHEC(/ICHK EXTERNAL OPLI PV=PVHA(HA) A=51 3 8:52 CALL ZEROIN (A.B,EPS,OPL1) DPHAS=(A+B)/2o IF(ICHK.EQ.1) PRINT 210,61,G2,DPHAS . 21o FORMAT(1HO,130(1H5)/10X5AHH A a N I N G--NO OR MULTIPLE ROOTS oere 10150 av ZEROIN/3OX33HCALLED FROM DPHAS . GUESSES ARE.621.1«,4N A 2ND,GZZ.1h/3OX17HRETURNED VALUE 15,521.1u» RETURN END FUNCTION-OPL1(TDP) COMMON /SPECIAL/PV,DB,XTRA DPL1=PV-PSDB(TDP) RETURN END FUNCTION DPLZ (TDP) COMMON lSPECIAL/PV,DB,XTRA HA=HAOP(TDP) T2=HA’HLDB(TDP)+.Qh8'HA‘(D3-TDP) IF(TDP-Q91.69) 1,2,2 1 DPL2=XTRA+T2 -HA‘(193o35+.Q85‘(h91.69-TDP)) RETURN Z DPL2=XTRA+T2§HA’(TDP-991.69) RETURN END 3O FUNCTION ENDBOP (DB,DP) COHHON/PRESS/ PATH HASHADP(DP) T1=.2QOS’(DB-Q59.69)+HLDB(OP)‘HA+.4“8‘HA‘(DB-DP) IF(DP-h91.69) 1,2,2 ENDBDP=T1-HA’(lb3.35+.485‘(h91.69-DP)) RETURN ENDBDP=T1+HA‘(DP-h91.69) RETURN END FUNCTION HADBRH(DB,RH) HADBRH=HAPV(RH’PSDB(DB)) RETURN END FUNCTION HADP(DP) HADP=HAPV(PSDB(DP)) RETURN END FUNCTION HAPV(PV) COMMON/PRESS/PATN HAPV=.6219‘PV/(PATM-PV) RETURN END FUNCTION HLDB (DB) IF(DB-991.69) 192,2 HLDB=1220o8QQ-.OSC77’(DB-Q59.69) RETURN IF(DB-609.69) 39“,“ HLDB=IU75o6965-.56983‘(DB-459o69) RETURN HLDB=SQRT(1359673.21h-.9125275587‘OB'OB) RETURN END I It [In .1 J 31 FUNCTION P503 (03) ' DATA R,A,8.CpOgEgF,G/.3206182232E049-o27h055258361h26E05y.55189607 A6328951E029-o95137038Q112655E-1go215321191636356E‘99-o“62026656819 8982E-8,.2k1612720987hE01,.1215“6516706055E°2/ IF(DB-h91.69) 1,2,2 1 PSOB=EXP(23.392h-11286.6489/DB-.QGOST’ALOG(DB)) RETURN 2 PSDB=R’EXP((AODB’(B+DB'(C+DB’(0+DB‘E))))/(DB‘(F-G’DB))) RETURN END , FUNCTION PVHA (HA) COMMON/PRESS/PATH PVHA=HA’PATMI(.6219+HA) RETURN END FUNCTION PVOBVS(DB,VS) COMMON/PRESS/PATH PVDBVS=PATM~53o35’DB/VS/1hho RETURN END FUNCTION PVUBHB(DB,HB) COMMON/PRESS/PATN A=FSDB(HB) 3 B=.6219Q‘HLDB(HB)’PATM S C=.2905‘(A-PATM)‘(HB-DB) PVDBNB=(A’B-S‘PATM)/(B+.15577‘C) RETURN END FUNCTION RHFSPV(DI,02) A=01 S 8:02 5 GO TO 1 ENTRY RHDBHA A=PSDB(D1) S B=PVHA(02) 1 RHPSPV=BIA RETURN END 32 FUNCTION VSDBHA (DapMA) COMMON IPRESS/PATM VSDBHA853o35’DB’(.62199HA)/1hh./.6219/PATM RETURN END FUNCTION HBDBHAS(DB,HA,01,GZ,EPS) ' COMMON /SPEOIAL/PV,TB,XTRA COMMON/ICHEC