SUBROUTINE TYPE74 (TIME,XIN,OUT,T,DTDT,PAR,INFO,ICNTRL,*) C*********************************************************************** C* Copyright ASHRAE. Toolkit for HVAC System Energy Calculations C************************************************************************ C* SUBROUTINE: CCSIM C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: Model the performance of a counterflow C* crossflow cooling coil. The model accounts C* for condensation on the outside surface. C* Three conditions are possible: all wet, C* partially wet or all dry. Output includes C* outlet air temperature and humidity, outlet C* water temperature, sensible and total C* cooling capacities and the wet fraction of C* air-side surface area. C*********************************************************************** C* INPUT VARIABLES DESCRIPTION(UNITS) SAMPLE VALUES C* XIN(1) MLiq Liquid mass flow rate(kg/s) 4.2 C* XIN(2) TLiqEnt Entering water temperature(C) 5.5556 C* XIN(3) MAir Dry air mass flow rate(kg/s) 3.2 C* XIN(4) TAirEnt Entering air dry bulb temperature(C) 25.0 C* XIN(5) WAirEnt Entering air humidity ratio(-) .01 C* C* OUTPUT VARIABLES C* OUT(1) TLiqLvg Leaving water temperature(C) 9.16554 C* OUT(2) TAirLvg Leaving air dry bulb temperature(C) 11.0299 C* OUT(3) WAirLvg Leaving air humidity ratio(-) .0078074 C* OUT(4) QTot Total heat transfer rate(W) 63467.1 C* OUT(5) QSen Sensible heat transfer rate(W) 45779.5 C* OUT(6) FWet Fraction of surface area wet(-) 1.0 C* OUT(7) ErrStat Error status indicator,0=ok,1=error(-) 0.0 C* C* PARAMETERS C* PAR(1) MLiqRat Liquid mass flow rate at rating(kg/s) 4.2 C* PAR(2) TLiqRat Entering water temperature at rating(C) 5.5556 C* PAR(3) MAirRat Dry air mass flow rate at rating(kg/s) 6.4 C* PAR(4) TAirRat Entering air dry bulb temperature at rating(C) 26.6667 C* PAR(5) WAirRat Entering air humidity ratio at rating(-) .0112 C* PAR(6) QTotRat Total heat transfer rate at rating(W) 88000.0 C* PAR(7) QSenRat Sensible heat transfer rate at rating(W) 66000.0 C* C* PROPERTIES C* CpAir Dry air specific heat (J/kg C) C* CpVap Water vapor specific heat (J/kg C) C* CpLiq Liquid specific heat (J/kg C) C*********************************************************************** C MAJOR RESTRICTIONS: General application is for heat exchanger C with four or more rows in a counterflow C configuration. C Approximates part-wet operation as C either fully wet or fully dry. C Constant UA. C C DEVELOPER: Michael J. Brandemuehl, PhD, PE C University of Colorado at Boulder C C DATE: January 1, 1992 C C INCLUDE FILES: coilsim.inc C prop.inc C SUBROUTINES CALLED: DRYCOIL C WETCOIL C HEATEX C UAHX C BYPASS C FUNCTIONS CALLED: DEWPOINT C DRYBULB C ENTHALPY C ENTHSAT C HUMTH C C REVISION HISTORY: None C C REFERENCE: TRNSYS. 1990. A Transient System C Simulation Program: Reference Manual. C Solar Energy Laboratory, Univ. Wisconsin- C Madison, pp. 4.6.8-1 - 4.6.8-12. C C Threlkeld, J.L. 1970. Thermal C Environmental Engineering, 2nd Edition, C Englewood Cliffs: Prentice-Hall,Inc. C pp. 254-270. C*********************************************************************** C INTERNAL VARIABLES C P(UAExt) Overall external dry UA/total external area (W/C) C P(UAInt) Overall internal UA/total external area (W/C) C P(UATot) Overall heat transfer coefficient (W/C) C uaH Enthalpy-based overall transfer coefficient (kg/s) C configHX Heat exchanger configuration (-) C hAirRat Entering air enthalpy at rating (J/kg) C hAirLvgRat Leaving air enthalpy at rating (J/kg) C hAdpRat Air enthalpy at apparatus dew point at rating (J/kg) C hLiqRatSat Saturated enthalpy at entering liquid temp (J/kg) C tAirLvgRat Leaving air temperature at rating (C) C tDewRat Entering air dewpoint at rating (C) C tSurfEnt Coil surface temperature at air entrance (C) C capAir Air-side capacity rate (W/C) C capAirH Enthalpy-based air-side capacity rate (kg/s) C capLiqH Enthalpy-based liquid-side capacity rate (kg/s) C small Small number in place of zero C large Large number in place of infinity C*********************************************************************** DOUBLE PRECISION XIN,OUT DIMENSION XIN(5),OUT(7),PAR(7),INFO(15) INTEGER Patm,CpAir,CpWat,CpLiq,CpVap,DViscAir, & DViscLiq,KAir,KLiq,RhoLiq,Hfg,RAir, & TKelMult,TAbsAdd,PaMult,PAbsAdd,ERRSTAT REAL Prop(16),LARGE,UATOT,UAINT,UAEXT,PAR INTEGER INFO,IOPT,NI,NP,ND CHARACTER*3 YCHECK(5),OCHECK(7) COMMON /LUNITS/LUR,LUW,IFORM,LUK PARAMETER (Patm = 1) PARAMETER (CpAir = 2) PARAMETER (CpWat = 3) PARAMETER (CpVap = 4) PARAMETER (CpLiq = 5) PARAMETER (DViscAir = 6) PARAMETER (DViscLiq = 7) PARAMETER (KAir = 8) PARAMETER (KLiq = 9) PARAMETER (RhoLiq = 10) PARAMETER (Hfg = 11) PARAMETER (RAir = 12) PARAMETER (TKelMult = 13) PARAMETER (TAbsAdd = 14) PARAMETER (PaMult = 15) PARAMETER (PAbsAdd = 16) PROP(PATM) = 101325.0 PROP(CPAIR) = 1006.0 PROP(CPWAT) = 4186.0 PROP(CPVAP) = 1805.0 PROP(CPLIQ) = 4186.0 PROP(DVISCAIR) = .0000182 PROP(DVISCLIQ) = .00144 PROP(KAIR) = .026 PROP(KLIQ) = .604 PROP(RHOLIQ) = 998.0 PROP(HFG) = 2501000.0 PROP(RAIR) = 287.055 PROP(TKELMULT) = 1.0 PROP(TABSADD) = 273.15 PROP(PAMULT) = 1.0 PROP(PABSADD) = 0.0 DATA small/1.E-9/, large /1.E20/, configHX /1./ DATA YCHECK/'MF2','TE1','MF2','TE1','DM1'/ DATA OCHECK/'TE1','TE1','DM1','PW2','PW2','DM1','DM1'/ ErrStat = 0 IOPT = -1 NI = 5 !CORRECT NUMBER OF INPUTS NP = 7 !CORRECT NUMBER OF PARAMETERS ND = 0 !CORRECT NUMBER OF DERIVATIVES MLIQ = XIN(1) TLIQENT = XIN(2) MAIR = XIN(3) TAIRENT = XIN(4) WAIRENT = XIN(5) MLIQRAT = PAR(1) TLIQRAT = PAR(2) MAIRRAT = PAR(3) TAIRRAT = PAR(4) WAIRRAT = PAR(5) QTOTRAT = PAR(6) QSENRAT = PAR(7) IF (INFO(7).EQ.-1) THEN CALL TYPECK(IOPT,INFO,NI,NP,ND) C CHECKS TO SEE IF USER'S INFO MATCHES CORRECT NUMBERS CALL RCHECK(INFO,YCHECK,OCHECK) C CHECKS TO SEE IF INPUT AND OUTPUT UNITS MATCH ENDIF C2********************************************************************** C2 The code between these bars of asterisks is used to set internal C2 parameters and is independent of component input values. In an C2 hourly simulation, this block of code may be skipped after the C2 first call. C1*** Calculate properties of air and liquid at rating point hAirRat = ENTHALPY(Prop,TAirRat,WAirRat) hLiqRatSat = ENTHSAT(Prop,TLiqRat) C1*** Calculate leaving air states at rating point hAirLvgRat = hAirRat - QTotRat/MAirRat hDummy = hAirRat - (QTotRat-QSenRat)/MAirRat wAirLvgRat = HUMTH(Prop,TAirRat,hDummy) tAirLvgRat = DRYBULB(Prop,hAirLvgRat,wAirLvgRat) C1*** Calculate coil UA assuming wet coil at rating C2*** Heat transfer in a wet coil is calculated based on enthalpy C2*** rather than temperature to include latent effects. Air enthalpies C2*** are evaluated using conventional psychrometric equations. The C2*** corresponding enthalpies of the coil and water are related to C2*** that of the air through "fictitious enthalpies," defined as the C2*** enthalpy of saturated air at the temperature of the coil or water. C2 C2*** While heat transfer rates are commonly expressed as the product C2*** of an overall heat transfer coefficient, UA, and a temperature C2*** difference, the use of enthalpy-based heat transfer calculations C2*** requires an enthalpy-based heat transfer coefficient, UAH. C2 C2*** q = UAH * (H1-H2) C2 C2*** where UAH = UA / cp C2*** UA = conventional heat transfer coefficient C2*** cp = specific heat across enthalpy difference C2 C2*** When using fictitious enthalpies, a corresponding fictitious C2*** specific heat must be defined. C2 C2*** EnthSat1-EnthSat2 = cpSat * (Temp1-Temp2) C2 C2*** UAH can be calculated from a combination of series or parallel C2*** enthalpy resistances, similar to thermal resistances modified for C2*** enthalpy as above. Enthalpy capacity rates relate heat transfer C2*** to the enthalpy change of a fluid between inlet and outlet. C2 C2*** q = CapH * (HAirLvg - HAirEnt) C2 C2*** On the air side, enthalpy capacity rate is the air mass flow rate. C2*** On the water side, the enthalpy capacity rate is based on the C2*** enthalpy of saturated air at the water temperature. C1*** Estimate cpSat using entering air dewpoint and water temperature tDewRat = DEWPOINT(Prop,WAirRat) cpSat = (ENTHSAT(Prop,tDewRat)-hLiqRatSat) & /(tDewRat-TLiqRat) C1*** Calculate overall heat transfer coefficient from fluid states C1*** and known total heat transfer capAirH = MAirRat capLiqH = MLiqRat * (Prop(CpLiq)/cpSat) uaH = UAHX(capAirH,hAirRat,capLiqH,hLiqRatSat,QTotRat, & configHX,ErrStat) C1*** Determine air-side coefficient, UAExt, assuming that the C1*** surface temperature is at the apparatus dewpoint temperature CALL BYPASS(Prop,TairRat,WAirRat,tAirLvgRat,wAirLvgRat, & tAdpRat,wAdpRat,bfRat,ErrStat) hAdpRat = ENTHALPY(Prop,tAdpRat,wAdpRat) IF (hAdpRat .LE. hLiqRatSat) THEN UAExt = uaH*Prop(CpAir) ELSE capAir = MAirRat*(Prop(CpAir)+WAirRat*Prop(CpVap)) UAExt = -LOG(bfRat)*capAir ENDIF C1*** Calculate liquid-side coefficient, UAInt, from enthalpy-based C1*** overall coefficient and air-side coefficient UAInt = cpSat/MIN((1./uaH - Prop(CpAir)/UAExt),large) UATot = 1./(1./UAExt+1./UAInt) C2********************************************************************** C1*** If both flows are zero, set outputs to inputs and return IF (ABS(MAir) .LT. small .AND. ABS(MLiq) .LT. small) THEN TLiqLvg = TLiqEnt TAirLvg = TAirEnt WAirLvg = WAirEnt GO TO 999 ENDIF C1*** IF coil is completely dry THEN tDewPt = DEWPOINT (Prop,WAirEnt) IF (tDewPt .LE. TLiqEnt) THEN C1*** Calculate the leaving conditions and performance of dry coil CALL DRYCOIL (Prop,MLiq,TLiqEnt,MAir,TAirEnt,WAirEnt, & UATot,configHX, & TLiqLvg,TAirLvg,WAirLvg,QTot,ErrStat) QSen = QTot FWet = 0. ELSE C1*** ELSE Assume external surface of coil is completely wet C1*** Calculate the leaving conditions and performance of wet coil CALL WETCOIL (Prop,MLiq,TLiqEnt,MAir,TAirEnt,WAirEnt, & UAInt,UAExt,configHX, & TLiqLvg,TAirLvg,WAirLvg,QTot,QSen,FWet, & tSurfEnt,ErrStat) C1*** IF coil is only partially wet THEN IF (tDewPt .LT. tSurfEnt) THEN C1*** Calculate the leaving conditions and performance of dry coil CALL DRYCOIL (Prop,MLiq,TLiqEnt,MAir,TAirEnt,WAirEnt, & UATot,configHX, & dryTLiqLvg,dryTAirLvg,dryWAirLvg,dryQTot,ErrStat) C1*** IF heat transfer from drycoil calculations is greater than that C1*** from wetcoil calculations THEN approximate the coil as dry. IF (dryQTot .GT. QTot) THEN TLiqLvg = dryTLiqLvg TAirLvg = dryTAirLvg WAirLvg = dryWAirLvg QTot = dryQTot QSen = QTot FWet = 0. ENDIF ENDIF ENDIF 999 CONTINUE OUT(1) = TLIQLVG OUT(2) = TAIRLVG OUT(3) = WAIRLVG OUT(4) = QTOT OUT(5) = QSEN OUT(6) = FWET OUT(7) = ERRSTAT RETURN 1 END C*********************************************************************** C* FILE: PROP.INC C* C* This file assigns a numbers to air and water property names to be C* used in the "Prop" array. C*********************************************************************** C DEVELOPER: Inger Andresen C Michael J. Brandemuehl, PhD, PE C C DATE: July 1, 1991 C C FILES REQUIRED: None C*********************************************************************** C INTERNAL VARIABLES: C Patm Atmospheric pressure (Pa) C CpAir Specific heat of dry air (J/kg C) C CpLiq Specific heat of liquid water (J/kg C) C CpVap Specific heat of saturated water vapor (J/kg C) C DViscAir Air dynamic viscosity (kg/m s) C DViscLiq Liquid dynamic viscosity (kg/m s) C KAir Air thermal conductivity (W/m C) C KLiq Liquid thermal conductivity (W/m C) C RhoLiq Liquid density (kg/m3) C Hfg Latent heat of vaporization of water (J/kg) C RAir Gas constant for air (J/kg C) C TKelMult Multiplying factor to convert user T to Kelvin C TAbsAdd Additive factor to convert user P to Kelvin C tKel = Prop(TKelMult)*T + Prop(TKelAdd) C PaMult Multiplying factor to convert user P to Pascals C PAbsAdd Additive factor to convert user P to Pascals C Pa = Prop(PaMult)*P + Prop(PaAdd) C*********************************************************************** C C INTEGER Patm,CpAir,CpWat,CpLiq,CpVap,DViscAir, C & DViscLiq,KAir,KLiq,RhoLiq,Hfg,RAir, C & TKelMult,TAbsAdd,PaMult,PAbsAdd C REAL Prop(16) C C PARAMETER (Patm = 1) C PARAMETER (CpAir = 2) C PARAMETER (CpWat = 3) C PARAMETER (CpVap = 4) C PARAMETER (CpLiq = 5) C PARAMETER (DViscAir = 6) C PARAMETER (DViscLiq = 7) C PARAMETER (KAir = 8) C PARAMETER (KLiq = 9) C PARAMETER (RhoLiq = 10) C PARAMETER (Hfg = 11) C PARAMETER (RAir = 12) C PARAMETER (TKelMult = 13) C PARAMETER (TAbsAdd = 14) C PARAMETER (PaMult = 15) C PARAMETER (PAbsAdd = 16) SUBROUTINE DRYCOIL (Prop,MLiq,TLiqEnt,MAir,TAirEnt,WAirEnt, & UA,ConfigHX, & TLiqLvg,TAirLvg,WAirLvg,Q, & ErrStat) C*********************************************************************** C* Copyright ASHRAE. Toolkit for HVAC System Energy Calculations C*********************************************************************** C* SUBROUTINE: DRYCOIL C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: Calculate the performance of a sensible C* air-liquid heat exchanger. Calculated C* results include outlet air temperature C* and humidity, outlet water temperature, C* and heat transfer rate. C*********************************************************************** C* INPUT VARIABLES C* MLiq Liquid mass flow rate (kg/s) C* TLiqEnt Entering water temperature (C) C* MAir Dry air mass flow rate (kg/s) C* TAirEnt Entering air dry bulb temperature (C) C* WAirEnt Entering air humidity ratio (-) C* C* UA Overall heat transfer coefficient (W/C) C* ConfigHX Heat exchanger configuration (-) C* 1 - Counterflow C* 2 - Parallel flow C* 3 - Cross flow, both streams unmixed C* 4 - Cross flow, both streams mixed C* 5 - Cross flow, stream 1 unmixed C* 6 - Cross flow, stream 2 unmixed C* C* OUTPUT VARIABLES C* TLiqLvg Leaving water temperature (C) C* TAirLvg Leaving air dry bulb temperature (C) C* WAirLvg Leaving air humidity ratio (-) C* Q Heat transfer rate (W) C* ErrStat Error status indicator, 0 = ok, 1 = error (-) C* C* PROPERTIES C* CpAir Specific heat of air (J/kg C) C* CpVap Specific heat of water vapor (J/kg C) C* CpLiq Specific heat of liquid (J/kg C) C*********************************************************************** C MAJOR RESTRICTIONS: Models coil using effectiveness-Ntu model. C C DEVELOPER: Shauna Gabel C Michael J. Brandemuehl, PhD, PE C University of Colorado at Boulder C C DATE: January 1, 1992 C C INCLUDE FILES: prop.inc C SUBROUTINES CALLED: HEATEX C FUNCTIONS CALLED: None C C REVISION HISTORY: None C C REFERENCE: Kays, W.M. and A.L. London. 1964. C Compact Heat Exchangers, 2nd Edition, C New York: McGraw-Hill. C C Threlkeld, J.L. 1970. Thermal C Environmental Engineering, 2nd Edition, C Englewood Cliffs: Prentice-Hall, Inc. C pp. 254-270. C*********************************************************************** C INTERNAL VARIABLES: C capAir Air-side capacity rate (W/C) C capLiq Water-side capacity rate (W/C) C*********************************************************************** INTEGER Patm,CpAir,CpWat,CpLiq,CpVap,DViscAir, & DViscLiq,KAir,KLiq,RhoLiq,Hfg,RAir, & TKelMult,TAbsAdd,PaMult,PAbsAdd REAL Prop(16) PARAMETER (Patm = 1) PARAMETER (CpAir = 2) PARAMETER (CpWat = 3) PARAMETER (CpVap = 4) PARAMETER (CpLiq = 5) PARAMETER (DViscAir = 6) PARAMETER (DViscLiq = 7) PARAMETER (KAir = 8) PARAMETER (KLiq = 9) PARAMETER (RhoLiq = 10) PARAMETER (Hfg = 11) PARAMETER (RAir = 12) PARAMETER (TKelMult = 13) PARAMETER (TAbsAdd = 14) PARAMETER (PaMult = 15) PARAMETER (PAbsAdd = 16) REAL MAir,MLiq INTEGER Errstat Errstat = 0 C2*** Calculate air and water capacity rates capAir = MAir*(Prop(CpAir)+WAirEnt*Prop(CpVap)) capLiq = MLiq*Prop(CpLiq) C1*** Determine the air and water outlet conditions CALL HEATEX (capLiq,TLiqEnt,capAir,TAirEnt,UA,ConfigHX, & TLiqLvg,TAirLvg) C1*** Calculate the total and sensible heat transfer rate Q = capAir*(TAirEnt-TAirLvg) WAirLvg = WAirEnt RETURN END SUBROUTINE WETCOIL (Prop,MLiq,TLiqEnt,MAir,TAirEnt,WAirEnt, & UAIntTot,UAExtTot,ConfigHX, & TLiqLvg,TAirLvg,WAirLvg,QTot,QSen,FWet, & TSurfEnt,ErrStat) C*********************************************************************** C* Copyright ASHRAE. Toolkit for HVAC System Energy Calculations C*********************************************************************** C* SUBROUTINE: WETCOIL C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: Calculate the performance of a cooling C* coil when the external fin surface is C* complete wet. Results include C* outlet air temperature and humidity, C* outlet water temperature, sensible and C* total cooling capacities, and the wet C* fraction of the air-side surface area. C*********************************************************************** C* INPUT VARIABLES C* MLiq Liquid mass flow rate (kg/s) C* TLiqEnt Entering water temperature (C) C* MAir Dry air mass flow rate (kg/s) C* TAirEnt Entering air dry bulb temperature (C) C* WAirEnt Entering air humidity ratio (-) C* C* UAIntTot Internal overall heat transfer coefficient (W/m2 C) C* UAExtTot External overall heat transfer coefficient (W/m2 C) C* ConfigHX Heat exchanger configuration (-) C* 1 - Counterflow C* 2 - Parallel flow C* 3 - Cross flow, both streams unmixed C* 4 - Cross flow, both streams mixed C* 5 - Cross flow, stream 1 unmixed C* 6 - Cross flow, stream 2 unmixed C* C* OUTPUT VARIABLES C* TLiqLvg Leaving water temperature (C) C* TAirLvg Leaving air dry bulb temperature (C) C* WAirLvg Leaving air humidity ratio (-) C* QTot Total heat transfer rate (W) C* QSen Sensible heat transfer rate (W) C* FWet Fraction of surface area wet (-) C* TSurfEnt Surface temperature at air entrance (C) C* ErrStat Error status indicator, 0 = ok, 1 = error (-) C* C* PROPERTIES C* CpLiq Specific heat of liquid (J/kg C) C* CpAir Specific heat of dry air (J/kg C) C*********************************************************************** C MAJOR RESTRICTIONS: Models coil as counterflow heat exchanger C Approximates saturated air enthalpy as C a linear function of temperature C C DEVELOPER: Shauna Gabel C Michael J. Brandemuehl, PhD, PE C University of Colorado at Boulder C C DATE: January 1, 1992 C C INCLUDE FILES: prop.inc C SUBROUTINES CALLED: HEATEX C WCOILOUT C FUNCTIONS CALLED: ENTHALPY C ENTHSAT C TAIRSAT C C REVISION HISTORY: None C C REFERENCE: Elmahdy, A.H. and Mitalas, G.P. 1977. C "A Simple Model for Cooling and C Dehumidifying Coils for Use In Calculating C Energy Requirements for Buildings," C ASHRAE Transactions,Vol.83 Part 2, C pp. 103-117. C C TRNSYS. 1990. A Transient System C Simulation Program: Reference Manual. C Solar Energy Laboratory, Univ. Wisconsin- C Madison, pp. 4.6.8-1 - 4.6.8-12. C C Threlkeld, J.L. 1970. Thermal C Environmental Engineering, 2nd Edition, C Englewood Cliffs: Prentice-Hall,Inc. C pp. 254-270. C*********************************************************************** C INTERNAL VARIABLES: C extResist Air-side resistance to heat transfer (m2 C/W) C intResist Liquid-side resistance to heat transfer (m2 C/W) C tDewEnt Entering air dew point (C) C uaH Overall enthalpy heat transfer coefficient (kg/s) C capAirWet Air-side capacity rate (kg/s) C capLiqWet Liquid-side capacity rate (kg/s) C resistRatio Ratio of resistances (-) C hAirLvg Outlet air enthalpy C hLiqEntSat Saturated enthalpy of air at (J/kg) C entering water temperature C hLiqLvgSat Saturated enthalpy of air at exit (J/kg) C water temperature C hSurfEntSat Saturated enthalpy of air at (J/kg) C entering surface temperature C hSurfLvgSat Saturated enthalpy of air at exit (J/kg) C surface temperature C cpSat Coefficient for equation below (J/kg C) C EnthSat1-EnthSat2 = cpSat*(TSat1-TSat2) C (all water and surface temperatures are C related to saturated air enthalpies for C wet surface heat transfer calculations) C************************************************************************ INTEGER Patm,CpAir,CpWat,CpLiq,CpVap,DViscAir, & DViscLiq,KAir,KLiq,RhoLiq,Hfg,RAir, & TKelMult,TAbsAdd,PaMult,PAbsAdd REAL Prop(16) PARAMETER (Patm = 1) PARAMETER (CpAir = 2) PARAMETER (CpWat = 3) PARAMETER (CpVap = 4) PARAMETER (CpLiq = 5) PARAMETER (DViscAir = 6) PARAMETER (DViscLiq = 7) PARAMETER (KAir = 8) PARAMETER (KLiq = 9) PARAMETER (RhoLiq = 10) PARAMETER (Hfg = 11) PARAMETER (RAir = 12) PARAMETER (TKelMult = 13) PARAMETER (TAbsAdd = 14) PARAMETER (PaMult = 15) PARAMETER (PAbsAdd = 16) REAL MAir,MLiq,intResist INTEGER ErrStat DATA small/1.E-9/ FWet = 1. extResist = 1./UAExtTot intResist = 1./UAIntTot C1*** Calculate enthalpies of entering air and water hAirEnt = ENTHALPY(Prop,TAirEnt,WAirEnt) hLiqEntSat = ENTHSAT(Prop,TLiqEnt) C1*** Estimate cpSat using entering air dewpoint and water temperature tDewEnt = DEWPOINT(Prop,WAirEnt) cpSat = (ENTHSAT(Prop,tDewEnt)-hLiqEntSat) & /(tDewEnt-TLiqEnt) C1*** Enthalpy-based heat transfer calculations C2*** Heat transfer in a wet coil is calculated based on enthalpy C2*** rather than temperature to include latent effects. Air enthalpies C2*** are evaluated using conventional psychrometric equations. The C2*** corresponding enthalpies of the coil and water are related to C2*** that of the air through "fictitious enthalpies," defined as the C2*** enthalpy of saturated air at the temperature of the coil or water. C2 C2*** While heat transfer rates are commonly expressed as the product C2*** of an overall heat transfer coefficient, UA, and a temperature C2*** difference, the use of enthalpy-based heat transfer calculations C2*** requires an enthalpy-based heat transfer coefficient, UAH. C2 C2*** q = UAH * (H1-H2) C2 C2*** where UAH = UA / cp C2*** UA = conventional heat transfer coefficient C2*** cp = specific heat across enthalpy difference C2 C2*** When using fictitious enthalpies, a corresponding fictitious C2*** specific heat must be defined. C2 C2*** EnthSat1-EnthSat2 = cpSat * (Temp1-Temp2) C2 C2*** UAH can be calculated from a combination of series or parallel C2*** enthalpy resistances, similar to thermal resistances modified for C2*** enthalpy as above. Enthalpy capacity rates relate heat transfer C2*** to the enthalpy change of a fluid between inlet and outlet. C2 C2*** q = CapH * (HAirLvg - HAirEnt) C2 C2*** On the air side, enthalpy capacity rate is the air mass flow rate. C2*** On the water side, the enthalpy capacity rate is based on the C2*** enthalpy of saturated air at the water temperature. C1*** Determine air and water enthalpy outlet conditions by modeling C1*** coil as counterflow enthalpy heat exchanger uaH = 1./(cpSat*intResist+Prop(CpAir)*extResist) capAirWet = MAir capLiqWet = MLiq * (Prop(CpLiq)/cpSat) CALL HEATEX(capAirWet,hAirEnt,capLiqWet,hLiqEntSat,uaH, & ConfigHX,hAirLvg,hLiqLvgSat) C1*** Calculate entering and leaving external surface conditions from C1*** air and water conditions and the ratio of resistances resistRatio = (intResist)/(intResist + & Prop(CpAir)/cpSat*extResist) hSurfEntSat = hLiqLvgSat + resistRatio*(hAirEnt-hLiqLvgSat) hSurfLvgSat = hLiqEntSat + resistRatio*(hAirLvg-hLiqEntSat) TSurfEnt = TAIRSAT(Prop,hSurfEntSat) C1*** Calculate outlet air temperature and humidity from enthalpies and C1*** surface conditions. QTot = MAir*(hAirEnt-hAirLvg) TLiqLvg = TLiqEnt+QTot/MAX(MLiq,small)/Prop(CpLiq) CALL WCOILOUT (Prop,MAir,TAirEnt,WAirEnt,hAirEnt,hAirLvg, & UAExtTot,TAirLvg,WAirLvg,QSen,ErrStat) 999 RETURN END SUBROUTINE HEATEX (Cap1,In1,Cap2,In2,UA,ConfigHX,Out1,Out2) C*********************************************************************** C* Copyright ASHRAE. Toolkit for HVAC System Energy Calculations C*********************************************************************** C* SUBROUTINE: HEATEX C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: Calculate the outlet states of a simple C* heat exchanger using the effectiveness-Ntu C* method of analysis. C*********************************************************************** C* INPUT VARIABLES C* Cap1 Capacity rate of stream 1 (W/C) C* In1 Inlet state of stream 1 (C) C* Cap2 Capacity rate of stream 2 (W/C) C* In2 Inlet state of stream 2 (C) C* UA Overall heat transfer coefficient (W/C) C* ConfigHX Heat exchanger configuration (-) C* 1 - Counterflow C* 2 - Parallel flow C* 3 - Cross flow, both streams unmixed C* 4 - Cross flow, both streams mixed C* 5 - Cross flow, stream 1 unmixed C* 6 - Cross flow, stream 2 unmixed C* C* OUTPUT VARIABLES C* Out1 Outlet state of stream 1 (C) C* Out2 Outlet state of stream 2 (C) C*********************************************************************** C MAJOR RESTRICTIONS: None C C DEVELOPER: Shauna Gabel C Michael J. Brandemuehl, PhD, PE C University of Colorado at Boulder C C DATE: January 1, 1992 C C INCLUDE FILES: None C SUBROUTINES CALLED: None C FUNCTIONS CALLED: None C C REVISION HISTORY: None C C REFERENCE: Kays, W.M. and A.L. London. 1964. C Compact Heat Exchangers, 2nd Ed., McGraw- C Hill: New York. C*********************************************************************** C* INTERNAL VARIABLES: C* cMin Minimum capacity rate of the streams (W/C) C* cMax Maximum capacity rate of the streams (W/C) C* cRatio Ratio of minimum to maximum capacity rate C* ntu Number of transfer units (-) C* effectiveness Heat exchanger effectiveness (-) C* qMax Maximum heat transfer possible (W) C*********************************************************************** REAL ntu,qMax,In1,In2,large DATA small/1.E-15/, large/1.E15/ C1*** Ntu and Cmin/Cmax (cRatio) calculations cMin = MIN(Cap1,Cap2) cMax = MAX(Cap1,Cap2) IF( cMax .EQ. 0.) THEN cRatio = 1. ELSE cRatio = cMin/cMax ENDIF IF( cMin .EQ. 0.) THEN ntu = large ELSE ntu = ua/cMin ENDIF C1*** Calculate effectiveness for special limiting cases mode = NINT(ConfigHX) IF(ntu .LE. 0) THEN effectiveness = 0. ELSE IF(cRatio .LT. small) THEN C2*** Cmin/Cmax = 0 and effectiveness is independent of configuration effectiveness = 1 - EXP(-ntu) C1*** Calculate effectiveness depending on heat exchanger configuration ELSE IF (mode .EQ. 1) THEN C2*** Counterflow IF (ABS(cRatio-1.) .LT. small) THEN effectiveness = ntu/(ntu+1.) ELSE e=EXP(-ntu*(1-cRatio)) effectiveness = (1-e)/(1-cRatio*e) ENDIF ELSE IF (mode .EQ. 2) THEN C2*** Parallel flow effectiveness = (1-EXP(-ntu*(1+cRatio)))/(1+cRatio) ELSE IF (mode .EQ. 3) THEN C2*** Cross flow, both streams unmixed eta = ntu**(-0.22) effectiveness = 1 - EXP((EXP(-ntu*cRatio*eta)-1)/(cRatio*eta)) ELSE IF (mode .EQ. 4) THEN C2*** Cross flow, both streams mixed effectiveness = ((1/(1-EXP(-ntu)))+ & (cRatio/(1-EXP(-ntu*cRatio)))-(1/(-ntu)))**(-1) ELSE C2*** One stream is mixed and one is unmixed. Determine whether the C2*** minimum or maximum capacity rate stream is mixed. IF ( (ABS(Cap1-cMin).LT.small .AND. mode.EQ.5) .OR. & (ABS(Cap2-cMin).LT.small .AND. mode.EQ.6) ) THEN C2*** Cross flow, stream with minimum capacity rate unmixed effectiveness = (1-EXP(-cRatio*(1-EXP(-ntu))))/cRatio ELSE C2*** Cross flow, stream with maximum capacity rate unmixed effectiveness = 1-EXP(-(1-EXP(-ntu*cRatio))/cRatio) ENDIF ENDIF C1*** Determine leaving conditions for the two streams qMax = MAX(cMin,small)*(In1-In2) Out1 = In1 - effectiveness*qMax/MAX(Cap1,small) Out2 = In2 + effectiveness*qMax/MAX(Cap2,small) RETURN END REAL FUNCTION UAHX (Cap1,In1,Cap2,In2,Q,ConfigHX,ErrStat) C*********************************************************************** C* Copyright ASHRAE. Toolkit for HVAC System Energy Calculations C*********************************************************************** C* FUNCTION: UAHX C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: Calculate the UA of a heat exchanger C* using the effectiveness-Ntu relationships C* given the entering capacity rate and C* temperature of each flow stream, the C* heat transfer rate under these conditions C* and the heat exchanger configuration. C*********************************************************************** C* INPUT VARIABLES C* Cap1 Capacity rate of stream 1 (W/C) C* In1 Inlet state of stream 1 (C) C* Cap2 Capacity rate of stream 2 (W/C) C* In2 Inlet state of stream 2 (C) C* Q Heat transfer rate (W) C* ConfigHX Heat exchanger configuration (-) C* 1 - Counterflow C* 2 - Parallel flow C* 3 - Cross flow, both streams unmixed C* 4 - Cross flow, both streams mixed C* 5 - Cross flow, stream 1 unmixed C* 6 - Cross flow, stream 2 unmixed C* C* OUTPUT VARIABLES C* UAHX Overall heat transfer coefficient (W/C) C* ErrStat Error status indicator, 0 = ok, 1 = error (-) C*********************************************************************** C MAJOR RESTRICTIONS: Models coil using effectiveness Ntu model C C DEVELOPER: Michael J. Brandemuehl, PhD, PE C University of Colorado at Boulder C C DATE: January 1, 1992 C C INCLUDE FILES: None C SUBROUTINES CALLED: HEATEX C FUNCTIONS CALLED: XITERATE C C REVISION HISTORY: None C C REFERENCE: None C*********************************************************************** C INTERNAL VARIABLES C qEstimate Estimated heat transfer in iteration (W) C ua Estimated heat transfer coefficient (W/C) C error Deviation of dependent variable in iteration C icvg Iteration convergence flag C iter Iteration index C itmax Maximum number of iterations C F1,F2 Previous values of error in iteration C X1,X2 Previous values of independent variable in iteration C*********************************************************************** REAL In1,In2 INTEGER ErrStat DATA itmax/20/ ErrStat = 0 C1*** Check for Q out of range (effectiveness > 1) IF(ABS(Q) .GT. ABS(MIN(Cap1,Cap2)*(In1-In2))) THEN WRITE(LUW,1001) 1001 FORMAT(/1X,'*** ERROR IN SUBROUTINE UAHX ***'/ & 1X,' Given Q is impossible for given inlet states'/) ErrStat = 1 ENDIF C1*** Estimate UAHX ua = ABS(Q/(In1-In2)) C1*** BEGIN LOOP to iteratively calculate UAHX DO 100 iter = 1,itmax C1*** Calculate heat transfer rate for estimated UAHX CALL HEATEX (Cap1,In1,Cap2,In2,ua,ConfigHx,out1,out2) qEstimate = Cap1*(In1-out1) C1*** Calculate new estimate for UAHX error = ABS(qEstimate) - ABS(Q) ua = XITERATE(ua,error,X1,F1,X2,F2,iter,icvg) C1*** If converged, leave loop IF (icvg .EQ. 1) GO TO 110 100 CONTINUE C1*** If not converged after itmax iterations, return error code WRITE(LUW,1005) itmax 1005 FORMAT(/1X,'*** ERROR IN SUBROUTINE UAHX ***'/ & 1X,' UA has not converged after',I2, & ' iterations'/) ErrStat = 1 110 CONTINUE UAHX = ua RETURN END SUBROUTINE BYPASS(Prop,TEnt,WEnt,TLvg,WLvg, & TAdp,WAdp,BF,ErrStat) C*********************************************************************** C* Copyright ASHRAE. Toolkit for HVAC System Energy Calculations C************************************************************************ C* SUBROUTINE: BYPASS C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: Calculate apparatus dew point and bypass C* factor given entering and leaving moist C* air conditions of cooling coil. C*********************************************************************** C* INPUT VARIABLES C* TEnt Entering air temperature (C) C* WEnt Entering air humidity ratio (-) C* TLvg Leaving air temperature (C) C* WLvg Leaving air humidity ratio (-) C* C* OUTPUT VARIABLES C* TAdp Apparatus dewpoint temperature (C) C* WAdp Apparatus dewpoint humidity ratio (-) C* BF Bypass factor (-) C* ErrStat Error status indicator, 0 = ok, 1 = error (-) C*********************************************************************** C MAJOR RESTRICTIONS: None C C DEVELOPER: Michael J. Brandemuehl, PhD, PE C University of Colorado at Boulder C C Hugh Henderson & Kannan Rengarajan C Florida Solar Energy Center C C DATE: January 1, 1992 C C INCLUDE FILES: prop.inc C SUBROUTINES CALLED: None C FUNCTIONS CALLED: DEWPOINT C ENTHALPY C HUMRATIO C SATPRESS C XITERATE C C REVISION HISTORY: None C C REFERENCE: 1989 ASHRAE Handbook - Fundamentals C*********************************************************************** C INTERNAL VARIABLES: C hEnt Entering air enthalpy C hLvg Leaving air enthalpy C hAdp Air enthalpy at apparatus dew point C slope Ratio temperature difference to humidity difference C between entering and leaving air states C tAdpEst Estimate of TAdp from slope C error Deviation of dependent variable in iteration C iter Iteration counter C icvg Iteration convergence flag C F1,F2 Previous values of dependent variable in XITERATE C X1,X2 Previous values of independent variable in XITERATE C************************************************************************ INTEGER Patm,CpAir,CpWat,CpLiq,CpVap,DViscAir, & DViscLiq,KAir,KLiq,RhoLiq,Hfg,RAir, & TKelMult,TAbsAdd,PaMult,PAbsAdd,ERRSTAT REAL Prop(16) PARAMETER (Patm = 1) PARAMETER (CpAir = 2) PARAMETER (CpWat = 3) PARAMETER (CpVap = 4) PARAMETER (CpLiq = 5) PARAMETER (DViscAir = 6) PARAMETER (DViscLiq = 7) PARAMETER (KAir = 8) PARAMETER (KLiq = 9) PARAMETER (RhoLiq = 10) PARAMETER (Hfg = 11) PARAMETER (RAir = 12) PARAMETER (TKelMult = 13) PARAMETER (TAbsAdd = 14) PARAMETER (PaMult = 15) PARAMETER (PAbsAdd = 16) DATA itmax/20/ C1*** Iterate to determine apparatus dewpoint at which the ADP C1*** equals the temperature calculated by extending the line between C1*** entering and leaving conditions to the saturation curve C1*** Calculate "slope" of temperature vs. humidity ratio between C1*** entering and leaving states slope = (TEnt-TLvg)/(WEnt-WLvg) C1*** Initialize iteration parameters TAdp = DEWPOINT(Prop,WLvg) DO 100 iter=1,itmax C1*** Calculate apparatus dewpoint and compare with predicted value C1*** using entering conditions and slope WAdp = HUMRATIO(Prop(Patm),SATPRESS(Prop,TAdp)) TAdpEst = TEnt - slope*(WEnt-WAdp) error = TAdp-TAdpEst TAdp = XITERATE(TAdp,error,X1,F1,X2,F2,iter,icvg) C1*** If converged, exit loop IF (icvg .EQ. 1) GO TO 110 100 CONTINUE C1*** Apparatus dewpoint has not converged after maximum iterations. C1*** Print error message, set return error flag, and RETURN WRITE(LUW,1001) itmax 1001 FORMAT(/1X,'*** ERROR IN SUBOUTINE BYPASS ***'/ & 1X,' Apparatus dewpoint has not ' & 'converged after ',I2,' iterations'/) ErrStat = 1 110 CONTINUE C1*** Calculate bypass factor from enthalpies hLvg = ENTHALPY(Prop,TLvg,WLvg) hEnt = ENTHALPY(Prop,TEnt,WEnt) hAdp = ENTHALPY(Prop,TAdp,WAdp) BF = (hLvg-hAdp)/(hEnt-hAdp) 999 RETURN END REAL FUNCTION DEWPOINT (Prop,W) C*********************************************************************** C* Copyright ASHRAE. Toolkit for HVAC System Energy Calculations C*********************************************************************** C* FUNCTION: DEWPOINT C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: Calculate the dewpoint temperature given C* humidity ratio C*********************************************************************** C* INPUT VARIABLES C* W Humidity ratio (-) C* C* OUTPUT VARIABLES C* DewPoint Dew point temperature of air (C) C* C* PROPERTIES C* Patm Atmospheric pressure (Pa) C*********************************************************************** C MAJOR RESTRICTIONS: None C C DEVELOPER: Michael J. Brandemuehl, PhD, PE C University of Colorado at Boulder C C DATE: January 1, 1992 C C INCLUDE FILES: None C SUBROUTINES CALLED: None C FUNCTIONS CALLED: None C C REVISION HISTORY: None C C REFERENCE: 1989 ASHRAE Handbook - Fundamentals C*********************************************************************** C INTERNAL VARIABLES: C pw Partial water vapor pressure (Pa) C small Small number C*********************************************************************** INTEGER Patm,CpAir,CpWat,CpLiq,CpVap,DViscAir, & DViscLiq,KAir,KLiq,RhoLiq,Hfg,RAir, & TKelMult,TAbsAdd,PaMult,PAbsAdd REAL Prop(16) PARAMETER (Patm = 1) PARAMETER (CpAir = 2) PARAMETER (CpWat = 3) PARAMETER (CpVap = 4) PARAMETER (CpLiq = 5) PARAMETER (DViscAir = 6) PARAMETER (DViscLiq = 7) PARAMETER (KAir = 8) PARAMETER (KLiq = 9) PARAMETER (RhoLiq = 10) PARAMETER (Hfg = 11) PARAMETER (RAir = 12) PARAMETER (TKelMult = 13) PARAMETER (TAbsAdd = 14) PARAMETER (PaMult = 15) PARAMETER (PAbsAdd = 16) DATA small/1.E-9/ C1*** Test for "dry" air IF (W .LT. small) THEN DewPoint = -999 ELSE C1*** Calculate the partial water vapor pressure as a function of C1*** humidity ratio. pw= Prop(Patm)*W/(.62198+W) C1*** Calculate dewpoint as saturation temperature at water vapor C1*** partial pressure DewPoint = SATTEMP(Prop,pw) ENDIF 999 RETURN END REAL FUNCTION DRYBULB (Prop,H,W) C*********************************************************************** C* Copyright ASHRAE. Toolkit for HVAC System Energy Calculations C*********************************************************************** C* FUNCTION: DRYBULB C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: Calculate the dry bulb temperature of C* moist air from enthalpy and humidity. C*********************************************************************** C* INPUT VARIABLES: C* H Enthalpy (J/kg) C* W Humidity ratio (-) C* C* OUTPUT VARIABLES: C* Drybulb Dry bulb temperature (C) C* C* PROPERTIES: C* CpAir Specific heat of air (J/kg C) C* CpVap Specific heat of water vapor (J/kg C) C* Hfg Reference heat of vaporization of water (J/kg) C*********************************************************************** C MAJOR RESTRICTIONS: Uses perfect gas relationships C Fit for enthalpy of saturated water vapor C C DEVELOPER: Shauna Gabel C Michael J. Brandemuehl, PhD, PE C University of Colorado at Boulder C C DATE: January 1, 1992 C C INCLUDE FILES: PROP.INC C SUBROUTINES CALLED: None C FUNCTIONS CALLED: None C C REVISION HISTORY: None C C REFERENCE: 1989 ASHRAE Handbook - Fundamentals C*********************************************************************** INTEGER Patm,CpAir,CpWat,CpLiq,CpVap,DViscAir, & DViscLiq,KAir,KLiq,RhoLiq,Hfg,RAir, & TKelMult,TAbsAdd,PaMult,PAbsAdd REAL Prop(16) PARAMETER (Patm = 1) PARAMETER (CpAir = 2) PARAMETER (CpWat = 3) PARAMETER (CpVap = 4) PARAMETER (CpLiq = 5) PARAMETER (DViscAir = 6) PARAMETER (DViscLiq = 7) PARAMETER (KAir = 8) PARAMETER (KLiq = 9) PARAMETER (RhoLiq = 10) PARAMETER (Hfg = 11) PARAMETER (RAir = 12) PARAMETER (TKelMult = 13) PARAMETER (TAbsAdd = 14) PARAMETER (PaMult = 15) PARAMETER (PAbsAdd = 16) C1*** Calculate the dry bulb temperature as a function of enthalpy and C1*** humidity ratio. C2*** hDryAir = Prop(CpAir)*TDB C2*** hSatVap = Prop(Hfg) + Prop(CpVap)*TDB C2*** Enthalpy = hDryAir + W*hSatVap Drybulb = (H-Prop(Hfg)*W)/(Prop(CpAir)+Prop(CpVap)*W) RETURN END REAL FUNCTION ENTHALPY (Prop,TDB,W) C*********************************************************************** C* Copyright ASHRAE. Toolkit for HVAC System Energy Calculations C*********************************************************************** C* FUNCTION: ENTHALPY C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: Calculate the enthalpy of moist air. C*********************************************************************** C* INPUT VARIABLES: C* TDB Dry bulb temperature (C) C* W Humidity ratio (-) C* C* OUTPUT VARIABLES: C* Enthalpy Enthalpy of moist air (J/kg) C* C* PROPERTIES: C* CpAir Specific heat of air (J/kg C) C* CpVap Specific heat of water vapor (J/kg C) C* Hfg Reference heat of vaporization of water (J/kg) C*********************************************************************** C MAJOR RESTRICTIONS Uses perfect gas relationships C Fit for enthalpy of saturated water vapor C C DEVELOPER: Shauna Gabel C Michael J. Brandemuehl, PhD, PE C University of Colorado at Boulder C C DATE: January 1, 1992 C C INCLUDE FILES: PROP.INC C SUBROUTINES CALLED: None C FUNCTIONS CALLED: None C C REVISION HISTORY: None C C REFERENCE: 1989 ASHRAE Handbook - Fundamentals C*********************************************************************** INTEGER Patm,CpAir,CpWat,CpLiq,CpVap,DViscAir, & DViscLiq,KAir,KLiq,RhoLiq,Hfg,RAir, & TKelMult,TAbsAdd,PaMult,PAbsAdd REAL Prop(16) PARAMETER (Patm = 1) PARAMETER (CpAir = 2) PARAMETER (CpWat = 3) PARAMETER (CpVap = 4) PARAMETER (CpLiq = 5) PARAMETER (DViscAir = 6) PARAMETER (DViscLiq = 7) PARAMETER (KAir = 8) PARAMETER (KLiq = 9) PARAMETER (RhoLiq = 10) PARAMETER (Hfg = 11) PARAMETER (RAir = 12) PARAMETER (TKelMult = 13) PARAMETER (TAbsAdd = 14) PARAMETER (PaMult = 15) PARAMETER (PAbsAdd = 16) C1*** Calculate the enthalpy as a function of dry bulb temperature and C1*** humidity ratio. hDryAir = Prop(CpAir)*TDB hSatVap = Prop(Hfg) + Prop(CpVap)*TDB Enthalpy = hDryAir + W*hSatVap RETURN END REAL FUNCTION ENTHSAT (Prop,TDB) C*********************************************************************** C* Copyright ASHRAE. Toolkit for HVAC System Energy Calculations C*********************************************************************** C* FUNCTION: ENTHSAT C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: Calculate the enthalpy at saturation C* for given dry bulb temperature C*********************************************************************** C* INPUT VARIABLES C* TDB Dry bulb temperature (C) C* C* OUTPUT VARIABLES C* EnthSat Enthalpy at saturation (J/kg) C* C* PROPERTIES C* Patm Atmospheric pressure (Pa) C*********************************************************************** C MAJOR RESTRICTIONS: None C C DEVELOPER: Shauna Gabel C Michael J. Brandemuehl, PhD, PE C University of Colorado at Boulder C C DATE: January 1, 1992 C C INCLUDE FILES: PROP.INC C SUBROUTINES CALLED: None C FUNCTIONS CALLED: SATPRESS C HUMRATIO C ENTHALPY C C REVISION HISTORY: None C C REFERENCE: 1989 ASHRAE Handbook - Fundamentals C*********************************************************************** C INTERNAL VARIABLES: C psat Saturated water vapor pressure (Pa) C w Humidity ratio (-) C*********************************************************************** INTEGER Patm,CpAir,CpWat,CpLiq,CpVap,DViscAir, & DViscLiq,KAir,KLiq,RhoLiq,Hfg,RAir, & TKelMult,TAbsAdd,PaMult,PAbsAdd REAL Prop(16) PARAMETER (Patm = 1) PARAMETER (CpAir = 2) PARAMETER (CpWat = 3) PARAMETER (CpVap = 4) PARAMETER (CpLiq = 5) PARAMETER (DViscAir = 6) PARAMETER (DViscLiq = 7) PARAMETER (KAir = 8) PARAMETER (KLiq = 9) PARAMETER (RhoLiq = 10) PARAMETER (Hfg = 11) PARAMETER (RAir = 12) PARAMETER (TKelMult = 13) PARAMETER (TAbsAdd = 14) PARAMETER (PaMult = 15) PARAMETER (PAbsAdd = 16) C1*** Calculate the saturation pressure at the given temperature. psat = SATPRESS (Prop,TDB) C1*** Calculate the humidity ratio from the saturation pressure w = HUMRATIO (Prop(Patm),psat) C1*** Calculate the enthalpy as a function of dry bulb temperature C1*** and humidity ratio. ENTHSAT = ENTHALPY (Prop,TDB,w) RETURN END REAL FUNCTION HUMRATIO (Patm,Pw) C*********************************************************************** C* Copyright ASHRAE. Toolkit for HVAC System Energy Calculations C*********************************************************************** C* FUNCTION: HUMRATIO C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: Calculate the humidity ratio from water C* vapor pressure and atmospheric pressure C*********************************************************************** C* INPUT VARIABLES C* Patm Atmospheric pressure (Pa) C* Pw Partial water vapor pressure (Pa) C* C* OUTPUT VARIABLES C* HumRatio Humidity ratio (-) C*********************************************************************** C MAJOR RESRICTIONS: None C C DEVELOPER: Shauna Gabel C Michael J. Brandemuehl, PhD, PE C University of Colorado at Boulder C C DATE: January 1, 1992 C C INCLUDE FILES: None C SUBROUTINES CALLED: None C FUNCTIONS CALLED: None C C REVISION HISTORY: None C C REFERENCE: 1989 ASHRAE Handbook - Fundamentals C************************************************************************ C1*** Calculate the humidity ratio. HumRatio = 0.62198*Pw/(Patm-Pw) RETURN END REAL FUNCTION HUMTH (Prop,TDB,H) C*********************************************************************** C* Copyright ASHRAE. Toolkit for HVAC System Energy Calculations C*********************************************************************** C* FUNCTION: HUMTH C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: Calculate the humidity ratio of moist air C* from dry bulb temperature and enthalpy. C*********************************************************************** C* INPUT VARIABLES: C* H Enthalpy (J/kg) C* TDB Dry bulb temperature (C) C* C* OUTPUT VARIABLES: C* HumTH Humidity ratio (-) C* C* PROPERTIES: C* CpAir Specific heat of air (J/kg C) C* CpVap Specific heat of water vapor (J/kg C) C* Hfg Reference heat of vaporization of water (J/kg) C*********************************************************************** C MAJOR RESTRICTIONS: Uses perfect gas relationships C Fit for enthalpy of saturated water vapor C C DEVELOPER: Shauna Gabel C Michael J. Brandemuehl, PhD, PE C University of Colorado at Boulder C C DATE: January 1, 1992 C C INCLUDE FILES: prop.inc C SUBROUTINES CALLED: None C FUNCTIONS CALLED: None C C REVISION HISTORY: None C C REFERENCE: 1989 ASHRAE Handbook - Fundamentals C*********************************************************************** INTEGER Patm,CpAir,CpWat,CpLiq,CpVap,DViscAir, & DViscLiq,KAir,KLiq,RhoLiq,Hfg,RAir, & TKelMult,TAbsAdd,PaMult,PAbsAdd REAL Prop(16) PARAMETER (Patm = 1) PARAMETER (CpAir = 2) PARAMETER (CpWat = 3) PARAMETER (CpVap = 4) PARAMETER (CpLiq = 5) PARAMETER (DViscAir = 6) PARAMETER (DViscLiq = 7) PARAMETER (KAir = 8) PARAMETER (KLiq = 9) PARAMETER (RhoLiq = 10) PARAMETER (Hfg = 11) PARAMETER (RAir = 12) PARAMETER (TKelMult = 13) PARAMETER (TAbsAdd = 14) PARAMETER (PaMult = 15) PARAMETER (PAbsAdd = 16) C1*** Calculate humidity ratio from dry bulb temperature and enthalpy C2*** hDryAir = Prop(CpAir)*TDB C2*** hSatVap = Prop(Hfg) + Prop(CpVap)*TDB C2*** Enthalpy = hDryAir + W*hSatVap HumTH = (H-Prop(CpAir)*TDB)/(Prop(Hfg)+Prop(CpVap)*TDB) RETURN END REAL FUNCTION RELHUM (Patm,Psat,HumRatio) C*********************************************************************** C* Copyright ASHRAE. Toolkit for HVAC System Energy Calculations C*********************************************************************** C* FUNCTION: RELHUM C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: Calculate the relative humidity from C* saturation and atmospheric pressures C*********************************************************************** C* INPUT VARIABLES C* Patm Atmospheric pressure (Pa) C* Psat Saturation pressure (Pa) C* HumRatio Humidity ratio (-) C* C* OUTPUT VARIABLES C* RelHum Relative humidity (-) C****** ****************************************************************** C MAJOR RESTRICTIONS: None C C DEVELOPER: Shauna Gabel C Michael J. Brandemuehl, PhD, PE C University of Colorado at Boulder C C DATE: January 1, 1992 C C INCLUDE FILES: None C SUBROUTINES CALLED: None C FUNCTIONS CALLED: None C C REVISION HISTORY: None C C REFERENCE: 1989 ASHRAE Handbook - Fundamentals C*********************************************************************** C INTERNAL VARIABLES: C pw Partial water vapor pressure (Pa) C*********************************************************************** C1*** Calculate the partial water vapor pressure as a function of C1*** humidity ratio. pw = Patm*HumRatio/(.62198+HumRatio) C1*** Calculate the relative humidity as a function of partial water C1*** vapor pressure and water vapor pressure at saturation. RelHum = pw/Psat RETURN END REAL FUNCTION RHODRY (Prop,TDB,W) C*********************************************************************** C* Copyright ASHRAE. Toolkit for HVAC System Energy Calculations C*********************************************************************** C* FUNCTION: RHODRY C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: Calculate dry air density. C*********************************************************************** C* INPUT VARIABLES C* TDB Dry bulb temperature (C) C* W Humidity ratio (-) C* C* OUTPUT VARIABLES C* RhoDry Density of dry air (kg/m3) C* C* PROPERTIES C* Patm Atmospheric pressure (Pa) C* RAir Gas constant for air (J/kg C) C* TAbsAdd Additive constant to convert user T to absolute T C*********************************************************************** C MAJOR RESTRICTIONS: Perfect gas relationships C C DEVELOPER: Shauna Gabel C Michael J. Brandemuehl, PhD, PE C University of Colorado at Boulder C C DATE: January 1, 1992 C C INCLUDE FILES: prop.inc C SUBROUTINES CALLED: None C FUNCTIONS CALLED: None C C REVISION HISTORY: None C C REFERENCE: 1989 ASHRAE Handbook - Fundamentals C*********************************************************************** C INTERNAL VARIABLES: C pAir Partial pressure of dry air (Pa) C*********************************************************************** INTEGER Patm,CpAir,CpWat,CpLiq,CpVap,DViscAir, & DViscLiq,KAir,KLiq,RhoLiq,Hfg,RAir, & TKelMult,TAbsAdd,PaMult,PAbsAdd REAL Prop(16) PARAMETER (Patm = 1) PARAMETER (CpAir = 2) PARAMETER (CpWat = 3) PARAMETER (CpVap = 4) PARAMETER (CpLiq = 5) PARAMETER (DViscAir = 6) PARAMETER (DViscLiq = 7) PARAMETER (KAir = 8) PARAMETER (KLiq = 9) PARAMETER (RhoLiq = 10) PARAMETER (Hfg = 11) PARAMETER (RAir = 12) PARAMETER (TKelMult = 13) PARAMETER (TAbsAdd = 14) PARAMETER (PaMult = 15) PARAMETER (PAbsAdd = 16) C1*** Calculate the dry air density from perfect gas laws. pAir = 0.62198*Prop(Patm)/(0.62198+W) RhoDry = pAir/Prop(RAir)/(TDB+Prop(TAbsAdd)) RETURN END REAL FUNCTION RHOMOIST (RhoDry,W) C*********************************************************************** C* Copyright ASHRAE. Toolkit for HVAC System Energy Calculations C*********************************************************************** C* FUNCTION: RHOMOIST C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: Calculate moist air density from dry air C* density and humidity ratio C*********************************************************************** C* INPUT VARIABLES: C* RhoDry Dry air density (kg/m3) C* W Humidity ratio (-) C* C* OUTPUT VARIABLES: C* RhoMoist Density of dry air (kg/m3) C*********************************************************************** C MAJOR RESTRICTIONS: None C C DEVELOPER: Shauna Gabel C Michael J. Brandemuehl, PhD, PE C University of Colorado at Boulder C C DATE: January 1, 1992 C C INCLUDE FILES: None C SUBROUTINES CALLED: None C FUNCTIONS CALLED: None C C REVISION HISTORY: None C C REFERENCE: 1989 ASHRAE Handbook - Fundamentals C************************************************************************ C1*** Calculate the moist air density RhoMoist = RhoDry*(1.+W) RETURN END REAL FUNCTION SATPRESS (Prop,T) C*********************************************************************** C* Copyright ASHRAE. Toolkit for HVAC System Energy Calculations C*********************************************************************** C* SUBROUTINE: SATPRESS C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: Calculate saturation pressure of water C* vapor as a function of temperature C*********************************************************************** C* INPUT VARIABLES C* T Temperature (C) C* C* OUTPUT VARIABLES C* SatPress Saturation pressure (Pa) C* C* PROPERTIES C* TKelMult Multiplying factor to convert user T to Kelvin C* TAbsAdd Additive factor to convert user T to absolute T C* tKel = Prop(TKelMult) * (T + Prop(TAbsAdd)) C* PaMult Multiplying factor to convert user P to Pascals C* PAbsAdd Additive factor to convert user P to absolute P C* Pa = Prop(PaMult) * (P + Prop(PAbsAdd)) C*********************************************************************** C MAJOR RESTRICTIONS: 173.16 K <= Temp <= 473.15 K C C DEVELOPER: Shauna Gabel C Michael J. Brandemuehl, PhD, PE C University of Colorado at Boulder C C DATE: January 1, 1992 C C INCLUDE FILES: prop.inc C SUBROUTINES CALLED: None C FUNCTIONS CALLED: None C C REVISION HISTORY: None C C REFERENCE: 1989 ASHRAE Handbook - Fundamentals C C Hyland, R.W., and A. Wexler. 1983. C Formulations for the thermodynamic C properties of the saturated phases of H2O C from 173.15 K to 473.15 K. ASHRAE C Transactions, Vol. 89, No. 2A, pp. 500-519 C*********************************************************************** C INTERNAL VARIABLES: C tKel Temperature in Kelvin (K) C pascals Saturation pressure (Pa) C*********************************************************************** INTEGER Patm,CpAir,CpWat,CpLiq,CpVap,DViscAir, & DViscLiq,KAir,KLiq,RhoLiq,Hfg,RAir, & TKelMult,TAbsAdd,PaMult,PAbsAdd REAL Prop(16) PARAMETER (Patm = 1) PARAMETER (CpAir = 2) PARAMETER (CpWat = 3) PARAMETER (CpVap = 4) PARAMETER (CpLiq = 5) PARAMETER (DViscAir = 6) PARAMETER (DViscLiq = 7) PARAMETER (KAir = 8) PARAMETER (KLiq = 9) PARAMETER (RhoLiq = 10) PARAMETER (Hfg = 11) PARAMETER (RAir = 12) PARAMETER (TKelMult = 13) PARAMETER (TAbsAdd = 14) PARAMETER (PaMult = 15) PARAMETER (PAbsAdd = 16) DATA C1/-5674.5359/,C2/6.3925247/,C3/-0.9677843E-2/ DATA C4/0.62215701E-6/,C5/0.20747825E-8/,C6/-0.9484024E-12/ DATA C7/4.1635019/,C8/-5800.2206/,C9/1.3914993/,C10/-0.048640239/ DATA C11/0.41764768E-4/,C12/-0.14452093E-7/,C13/6.5459673/ C1*** Convert temperature from user units to Kelvin. tKel = Prop(TKelMult)*(T+Prop(TAbsAdd)) C1*** If below freezing, calculate saturation pressure over ice. IF (tKel .LT. 273.15) THEN pascals = EXP(C1/tKel+C2+C3*tKel+C4*tKel**2+C5*tKel**3+C6* & tKel**4+C7*ALOG(tKel)) C1*** If above freezing, calculate saturation pressure over liquid water. ELSE IF (tKel .GE. 273.15) THEN pascals = EXP(C8/tKel+C9+C10*tKel+C11*tKel**2+C12*tKel**3+C13 & *ALOG(tKel)) ENDIF C1*** Convert pressure from Pascals to user units SatPress = pascals/Prop(PaMult) - Prop(PAbsAdd) RETURN END REAL FUNCTION SATTEMP (Prop,P) C*********************************************************************** C* Copyright ASHRAE. Toolkit for HVAC System Energy Calculations C*********************************************************************** C* FUNCTION: SATTEMP C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: Calculate the saturation (boiling) C* temperature of water given pressure C*********************************************************************** C* INPUT VARIABLES C* P Pressure (Pa) C* C* OUTPUT VARIABLES C* SatTemp Saturation temperature of water vapor (C) C*********************************************************************** C MAJOR RESTRICTIONS: None C C DEVELOPER: Shauna Gabel C Michael J. Brandemuehl, PhD, PE C University of Colorado at Boulder C C DATE: January 1, 1992 C C INCLUDE FILES: prop.inc C SUBROUTINES CALLED: None C FUNCTIONS CALLED: SATPRESS C XITERATE C C REVISION HISTORY: None C C REFERENCE: 1989 ASHRAE Handbook - Fundamentals C*********************************************************************** C INTERNAL VARIABLES: C tSat Water temperature guess (C) C pSat Pressure corresponding to temp. guess (Pa) C error Deviation of dependent variable in iteration C iter Iteration counter C icvg Iteration convergence flag C F1,F2 Previous values of dependent variable in XITERATE C X1,X2 Previous values of independent variable in XITERATE C*********************************************************************** DATA itmax/50/ C1*** Use an iterative process to determine the saturation temperature C1*** at a given pressure using a correlation of saturated water vapor C1*** pressure as a function of temperature C1*** Initial guess of boiling temperature tSat = 100. C1*** Iterate to find the saturation temperature C1*** of water given the total pressure C2*** Set iteration loop parameters DO 100 iter = 1,itmax C1*** Calculate saturation pressure for estimated boiling temperature pSat = SATPRESS(Prop,tSat) C1*** Compare with specified pressure and update estimate of temperature error = P - pSat tSat = XITERATE (tSat,error,X1,F1,X2,F2,iter,icvg) C2*** If converged leave loop iteration IF (icvg .EQ. 1) GO TO 110 C2*** Water temperature not converged, repeat calculations with new C2*** estimate of water temperature 100 CONTINUE C1*** Saturation temperature has not converged after maximum specified C1*** iterations. Print error message, set return error flag, and RETURN WRITE(LUW,1001) itmax 1001 FORMAT(/1X,'*** ERROR IN FUNCTION SatTemp ***'/ & 1X,' Saturation temperature has not ' & 'converged after ',I2,' iterations'/) 110 SatTemp = tSat RETURN END REAL FUNCTION TAIRSAT (Prop,HSat) C*********************************************************************** C* Copyright ASHRAE. Toolkit for HVAC System Energy Calculations C*********************************************************************** C* FUNCTION: ENTHSAT C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: Calculate the dry bulb temperature given C* enthalpy at saturation. C*********************************************************************** C* INPUT VARIABLES: C* HSat Enthalpy at saturation (J/kg) C* C* OUTPUT VARIABLES: C* TAirSat Dry bulb temperature (C) C*********************************************************************** C MAJOR RESTRICTIONS: None C C DEVELOPER: Michael J. Brandemuehl, PhD, PE C University of Colorado at Boulder C C DATE: January 1, 1992 C C INCLUDE FILES: prop.inc C SUBROUTINES CALLED: None C FUNCTIONS CALLED: ENTHSAT C C REVISION HISTORY: None C C REFERENCE: 1989 ASHRAE Handbook - Fundamentals C*********************************************************************** C INTERNAL VARIABLES: C error Deviation of dependent variable in iteration C iter Iteration counter C icvg Iteration convergence flag C F1,F2 Previous values of dependent variable in XITERATE C X1,X2 Previous values of independent variable in XITERATE C*********************************************************************** DATA itmax/20/,tSat/50./ C1*** Estimate saturation temperature if reasonable value not available IF(tSat .LT. -200. .OR. tSat .GT. 1000.) tSat = 50. C1*** Calculate saturation temperature by iteration using function to C1*** calculate saturation enthalpy from temperature DO 100 iter=1,itmax error = HSat - ENTHSAT(Prop,tSat) tSat = XITERATE(tSat,error,X1,F1,X2,F2,iter,icvg) C1*** If converged, leave iteration loop. IF (icvg .EQ. 1) GO TO 110 C1*** Temperature not converged, repeat calculation with new C1*** estimate of temperature. 100 CONTINUE C1*** Temperature has not converged after maximum specified C1*** iterations. Print error message and RETURN WRITE(LUW,1001) itmax 1001 FORMAT(/1X,'*** ERROR IN FUNCTION TAIRSAT ***'/ & 1X,' Temperature has not ' & 'converged after ',I2,' iterations'/) 110 CONTINUE TAirSat = tSat RETURN END REAL FUNCTION WETBULB (Prop,TDB,W) C*********************************************************************** C* Copyright ASHRAE. Toolkit for HVAC System Energy Calculations C*********************************************************************** C* FUNCTION: WETBULB C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: Calculate wet bulb temperature from dry C* bulb temperature and humidity ratio C*********************************************************************** C* INPUT VARIABLES C* TDB Dry bulb temperature (C) C* W Humidity ratio of air (-) C* C* OUTPUT VARIABLES C* WetBulb Wet bulb temperature (C) C* C* PROPERTIES: C* Patm Atmospheric pressure (Pa) C* Hfg Latent heat of vaporization of water (J/kg) C* CpAir Specific heat of air (J/kg C) C* CpVap Specific heat of water vapor (J/kg C) C* CpWat Specific heat of water (J/kg C) C*********************************************************************** C MAJOR RESTRICTIONS: None C C DEVELOPER: Shauna Gabel C Michael J. Brandemuehl, PhD, PE C University of Colorado at Boulder C C DATE: January 1, 1992 C C INCLUDE FILES: prop.inc C SUBROUTINES CALLED: None C FUNCTIONS CALLED: SATPRESS C HUMRATIO C SATTEMP C XITERATE C C REVISION HISTORY: None C C REFERENCE: 1989 ASHRAE Handbook - Fundamentals C*********************************************************************** C INTERNAL VARIABLES: C tBoil Boiling temperature of water at given pressure (C) C psatStar Saturation pressure at wet bulb temperature (C) C wStar Humidity ratio as a function of PsatStar (-) C newW Humidity ratio calculated with wet bulb guess (-) C error Deviation of dependent variable in iteration C iter Iteration counter C icvg Iteration convergence flag C F1,F2 Previous values of dependent variable in XITERATE C X1,X2 Previous values of independent variable in XITERATE C*********************************************************************** INTEGER Patm,CpAir,CpWat,CpLiq,CpVap,DViscAir, & DViscLiq,KAir,KLiq,RhoLiq,Hfg,RAir, & TKelMult,TAbsAdd,PaMult,PAbsAdd REAL Prop(16) PARAMETER (Patm = 1) PARAMETER (CpAir = 2) PARAMETER (CpWat = 3) PARAMETER (CpVap = 4) PARAMETER (CpLiq = 5) PARAMETER (DViscAir = 6) PARAMETER (DViscLiq = 7) PARAMETER (KAir = 8) PARAMETER (KLiq = 9) PARAMETER (RhoLiq = 10) PARAMETER (Hfg = 11) PARAMETER (RAir = 12) PARAMETER (TKelMult = 13) PARAMETER (TAbsAdd = 14) PARAMETER (PaMult = 15) PARAMETER (PAbsAdd = 16) REAL newW DATA itmax/20/ C1*** Initial temperature guess tBoil = SATTEMP (Prop,Prop(Patm)) WetBulb = MAX( MIN(WetBulb,TDB,(tBoil-0.1)), 0.) C1*** Begin iteration loop DO 100 iter = 1,itmax IF (WetBulb .GE. (tBoil-0.09) ) WETBULB = tBoil-0.1 C1*** Determine the saturation pressure for wet bulb temperature psatStar = SATPRESS (Prop,WetBulb) C1*** Determine humidity ratio for given saturation pressure wStar = HUMRATIO (Prop(Patm),psatStar) C1*** Calculate new humidity ratio and determine difference from known C1*** humidity ratio newW = ((Prop(Hfg)-(Prop(CpWat)-Prop(CpVap))*WetBulb)*wStar- & Prop(CpAir)*(TDB-WetBulb))/(Prop(Hfg)+Prop(CpVap)*TDB & -Prop(CpWat)*WetBulb) C1*** Check error, if not satisfied, calculate new guess and iterate error = W-newW WetBulb = XITERATE(WetBulb,error,X1,F1,X2,F2,iter,icvg) C1*** If converged, leave iteration loop. IF (icvg .EQ. 1) GO TO 900 C1*** Wet bulb temperature not converged, repeat calculation with new C1*** estimate of wet bulb temperature. 100 CONTINUE C1*** Wet bulb temperature has not converged after maximum specified C1*** iterations. Print error message, set return error flag, and RETURN WRITE(LUW,1009) itmax 1009 FORMAT(/1X,'*** ERROR IN FUNCTION WetBulb ***'/ & 1X,' Wet bulb temperature has not ' & 'converged after ',I2,' iterations'/) 900 IF (WetBulb .GT. TDB) WetBulb = TDB 999 RETURN END  REAL FUNCTION XITERATE (X0,F0,X1,F1,X2,F2,ICount,ICvg) C*********************************************************************** C* Copyright ASHRAE. Toolkit for HVAC System Energy Calculations C*********************************************************************** C* C* SUBROUTINE: XITERATE C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: Iterately solves for the value of X which C* satisfies F(X)=0. Given Xi,F(Xi) pairs, C* the subroutine tests for convergence and C* provides a new guess for the value of the C* independent variable X. C*********************************************************************** C* INPUT VARIABLES C* F0 Current value of the function F(X) C* X0 Current value of X C* F1,F2 Two previous values of F(Xi) C* X1,X2 Two previous values of X C* C* NOTE: F1,X1,F2,X2 MUST BE STORED AND SAVED IN CALLING C* ROUTINE. THEY NEED NO INITIALIZATION C* C* ICount Number of iterations C* C* OUTPUT VARIABLES C* XIterate New estimate of X for F(X)=0 C* ICvg Convergence flag ICvg = 0: Not converged C* ICvg = 1: Converged C*********************************************************************** C DEVELOPER: Michael J. Brandemuehl, PhD, PE C University of Colorado at Boulder C C DATE: January 1, 1992 C C INCLUDE FILES: None C SUBROUTINES CALLED: None C FUNCTIONS CALLED: None C C REVISION HISTORY: None C C REFERENCE: None C*********************************************************************** C INTERNAL VARIABLES C small Small number used in place of zero C mode Number of points used in fit C mode = 1: Use XPerburb to get new X C mode = 2: Linear equation to get new X C mode > 2: Quadratic equation to get new X C coef(i) Coefficients for quadratic fit C F(X) = coef(1) + coef(2)*X + coef(3)*X*X C check Term under radical in quadratic solution C FiQ,XiQ Double precision values of Fi,Xi C slope Slope for linear fit C tolRel Relative error tolerance C xPerturb Perturbation applied to X to initialize iteration C*********************************************************************** DOUBLE PRECISION coef(3),check,F0Q,F1Q,F2Q,X0Q,X1Q,X2Q DATA tolRel/1.E-5/,xPerturb/0.1/,small/1.E-9/ C1*** Check for convergence by comparing change in X IF ((ABS(X0-X1) .LT. tolRel*MAX(ABS(X0),small) .AND. & ICount .NE. 1) .OR. F0 .EQ. 0.) THEN XIterate = X0 ICvg=1 RETURN ENDIF C1*** Not converged. C2*** If after the second iteration there are enough previous points to C2 fit a quadratic for the new X. If the quadratic fit is not C2 applicable, mode will be set to 1 or 2 and a new X will be C2 determined by incrementing X from xPerturb or from a linear fit. ICvg=0 mode=ICount 10 IF (mode .EQ. 1) THEN C1*** New guess is specified by xPerturb IF (ABS(X0) .GT. small) THEN XIterate = X0*(1.+xPerturb) ELSE XIterate = xPerturb ENDIF ELSEIF (mode .EQ. 2) THEN C1*** New guess calculated from LINEAR FIT of most recent two points SLOPE=(F1-F0)/(X1-X0) IF(slope.EQ.0) THEN mode=1 GO TO 10 ENDIF XIterate=X0-F0/SLOPE ELSE C1*** New guess calculated from QUADRATIC FIT C1*** If two Xi are equal, set mode for linear fit and return to top IF (X0 .EQ. X1) THEN X1=X2 F1=F2 mode=2 GO TO 10 ELSEIF (X0 .EQ. X2) THEN mode=2 GO TO 10 ENDIF C1*** Determine quadratic coefficients from the three data points C1*** using double precision. F2Q=F2 F1Q=F1 F0Q=F0 X2Q=X2 X1Q=X1 X0Q=X0 coef(3)=((F2Q-F0Q)/(X2Q-X0Q)-(F1Q-F0Q)/(X1Q-X0Q))/(X2Q-X1Q) coef(2)=(F1Q-F0Q)/(X1Q-X0Q)-(X1Q+X0Q)*coef(3) coef(1)=F0-(coef(2)+coef(3)*X0Q)*X0Q C1*** If points are colinear, set mode for linear fit and return to top IF (ABS(coef(3)) .LT. 1.D-10) THEN mode=2 GO TO 10 ENDIF C1*** Check for precision. If the coefficients do not accurately C1*** predict the given data points due to round-off errors, set C1*** mode for a linear fit and return to top. IF (ABS((coef(1)+(coef(2)+coef(3)*X1Q)*X1Q-F1Q)/F1Q) .GT. & 1.D-4) THEN mode=2 GO TO 10 ENDIF C1*** Check for imaginary roots. If no real roots, set mode to C1*** estimate new X by simply incrementing by xPerturb check=coef(2)**2-4*coef(1)*coef(3) IF (check .LT. 0) THEN C1*** Imaginary roots -- go back to linear fit mode=2 GO TO 10 ELSEIF (check .GT. 0) THEN C1*** Real unequal roots -- determine root nearest to most recent guess XIterate=(-coef(2)+SQRT(check))/coef(3)/2 xOther=-XIterate-coef(2)/coef(3) IF (ABS(XIterate-X0) .GT. ABS(xOther-X0)) XIterate=xOther ELSE C1*** Real Equal Roots -- one solution XIterate=-coef(2)/coef(3)/2 ENDIF ENDIF C1*** Set previous variable values for the next iteration IF (mode .LT. 3) THEN C1*** No valid previous points to eliminate. X2=X1 F2=F1 X1=X0 F1=F0 ELSE C1*** Eliminate one previous point based on sign and magnitude of F(X) C2*** Keep the current point and eliminate one of the previous ones. IF (F1*F0 .GT. 0 .AND. F2*F0 .GT. 0) THEN C2*** All previous points of same sign. Eliminate one with biggest F(X) IF (ABS(F2) .GT. ABS(F1)) THEN X2=X1 F2=F1 ENDIF ELSE C1*** Points of different sign. C1*** Eliminate the previous one with the same sign as current F(X). IF (F2*F0 .GT. 0) THEN X2=X1 F2=F1 ENDIF ENDIF X1=X0 F1=F0 ENDIF RETURN END SUBROUTINE BESSEL (X,I0,I1,K0,K1) C*********************************************************************** C* SUBROUTINE: BESSEL C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: Determine modified Bessel functions of the C* first kind I0 and I1 and of the second C* kind K0 and K1 using a polynomial C* approximation. C* C* REFERENCE: Press, W.H. et al. 1986. Numerical C* Recipies, Cambridge University Press: C* Cambridge. C*********************************************************************** REAL X,I0,I1,K0,K1 DATA T1,T2,T3,T4,T5,T6,T7/1.0,3.5156229,3.0899424, & 1.2067492,0.2659732,0.0360768,0.0045813/ DATA U1,U2,U3,U4,U5,U6,U7,U8,U9/0.39894228,0.01328592, & 0.00225319,-0.00157565,0.00916281,-0.02057706, & 0.02635537,-0.01647633,0.00392377/ DATA V1,V2,V3,V4,V5,V6,V7/0.5,0.87890594,0.51498869, & 0.15084934,0.02658733,0.00301532,0.00032411/ DATA O1,O2,O3,O4,O5,O6,O7,O8,O9/0.39894228,-0.03988024, & -0.00362018,0.00163801,-0.01031555,0.02282967, & -0.02895312,0.01787654,-0.00420059/ DATA P1,P2,P3,P4,P5,P6,P7/-0.57721566,0.42278420, & 0.23069756,0.03488590,0.00262698,0.00010750,0.0000074/ DATA Q1,Q2,Q3,Q4,Q5,Q6,Q7/1.25331414,-0.07832358,0.02189568, & -0.01062446,0.00587872,-0.00251540,0.00053208/ DATA R1,R2,R3,R4,R5,R6,R7/1.0,0.15443144,-0.67278579, & -0.18156897,-0.01919402,-0.00110404,-0.00004686/ DATA S1,S2,S3,S4,S5,S6,S7/1.25331414,0.23498619,-0.03655620, & 0.01504268,-0.00780353,0.00325614,-0.00068245/ C1*** Returns I0 IF (ABS(X) .LT. 3.75) THEN Y = (X/3.75)**2 I0 = T1+Y*(T2+Y*(T3+Y*(T4+Y*(T5+Y*(T6+Y*T7))))) ELSE AX = ABS(X) Y = 3.75/AX I0 = (EXP(AX)/SQRT(AX))*(U1+Y*(U2+Y*(U3+Y*(U4+Y* & (U5+Y*(U6+Y*(U7+Y*(U8+Y*U9)))))))) ENDIF C1*** Returns I1 IF (ABS(X) .LT. 3.75) THEN Y = (X/3.75)**2 I1 = X*(V1+Y*(V2+Y*(V3+Y*(V4+Y*(V5+Y*(V6+Y*V7)))))) ELSE AX = ABS(X) Y=3.75/AX I1 = (EXP(AX)/SQRT(AX))*(O1+Y*(O2+Y*(O3+Y*(O4+Y* & (O5+Y*(O6+Y*(O7+Y*(O8+Y*O9)))))))) IF (X .LT. 0) I1 = -I1 ENDIF C1*** Returns K0 IF (X .LE. 2) THEN Y = X*X/4.0 K0 = (-LOG(X/2.0)*I0)+(P1+Y*(P2+Y*(P3+Y*(P4+Y*(P5+Y* & (P6+Y*P7)))))) ELSE Y = (2.0/X) K0 = (EXP(-X)/SQRT(X))*(Q1+Y*(Q2+Y*(Q3+Y*(Q4+Y*(Q5+Y* & (Q6+Y*Q7)))))) ENDIF C1*** Returns K1 IF (X .LE. 2.0) THEN Y = X*X/4.0 K1 = (LOG(X/2.0)*I1)+(1.0/X)*(R1+Y*(R2+Y*(R3+Y* & (R4+Y*(R5+Y*(R6+Y*R7)))))) ELSE Y = (2.0/X) K1 = (EXP(-X)/SQRT(X))*(S1+Y*(S2+Y*(S3+Y*(S4+Y* & (S5+Y*(S6+Y*S7)))))) ENDIF RETURN 1 END SUBROUTINE WCOILOUT (Prop,MAir,TAirEnt,WAirEnt,HAirEnt,HAirLvg, & UAExt,TAirLvg,WAirLvg,QSen,ErrStat) C*********************************************************************** C* Copyright ASHRAE. Toolkit for HVAC System Energy Calculations C*********************************************************************** C* SUBROUTINE: WCOILOUT C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: Calculate the leaving air temperature, C* the leaving air humidity ratio and the C* sensible cooling capacity of wet cooling C* coil. C*********************************************************************** C* INPUT VARIABLES C* MAir Dry air mass flow rate (kg/s) C* TAirEnt Entering air dry bulb temperature (C) C* WAirEnt Entering air humidity ratio (-) C* HAirEnt Entering air enthalpy (J/kg) C* HAirLvg Leaving air enthalpy (J/kg) C* UAExt Heat transfer coefficient for external surface (W/C) C* C* OUTPUT VARIABLES C* TAirLvg Leaving air dry bulb temperature (C) C* WAirLvg Leaving air humidity ratio (-) C* Qsen Sensible heat transfer rate (W) C* ErrStat Error status indicator, 0 = ok, 1 = error (-) C* C* PROPERTIES C* Patm Atmospheric pressure (-) C* CpAir Specific heat of air (J/kg C) C* CpVap Specific heat of water vapor (J/kg C) C*********************************************************************** C MAJOR RESTRICTIONS: Assumes condensate at uniform temperature. C C DEVELOPER: Shauna Gabel C Michael J. Brandemuehl, PhD, PE C University of Colorado at Boulder C C DATE: January 1, 1992 C C INCLUDE FILES: prop.inc C SUBROUTINES CALLED: None C FUNCTIONS CALLED: TAIRSAT C TDB_H C C REVISION HISTORY: None C C REFERENCE: Elmahdy, A.H. and Mitalas, G.P. 1977. C "A Simple Model for Cooling and C Dehumidifying Coils for Use In Calculating C Energy Requirements for Buildings," C ASHRAE Transactions,Vol.83 Part 2, C pp. 103-117. C*********************************************************************** C INTERNAL VARIABLES: C capAir Air capacity rate (W/C) C ntu Number of heat transfer units (-) C effectiveness Heat exchanger effectiveness (-) C hCondSat Saturated air enthalpy at temperature of C condensate (J/kg) C tempCond Temperature of condensate (C) C*********************************************************************** INTEGER Patm,CpAir,CpWat,CpLiq,CpVap,DViscAir, & DViscLiq,KAir,KLiq,RhoLiq,Hfg,RAir, & TKelMult,TAbsAdd,PaMult,PAbsAdd,ERRSTAT REAL Prop(16),NTU,MAIR PARAMETER (Patm = 1) PARAMETER (CpAir = 2) PARAMETER (CpWat = 3) PARAMETER (CpVap = 4) PARAMETER (CpLiq = 5) PARAMETER (DViscAir = 6) PARAMETER (DViscLiq = 7) PARAMETER (KAir = 8) PARAMETER (KLiq = 9) PARAMETER (RhoLiq = 10) PARAMETER (Hfg = 11) PARAMETER (RAir = 12) PARAMETER (TKelMult = 13) PARAMETER (TAbsAdd = 14) PARAMETER (PaMult = 15) PARAMETER (PAbsAdd = 16) DATA small/1.E-9/ ErrStat = 0 C1*** Determine the temperature effectiveness, assuming the temperature C1*** of the condensate is constant (Cmin/Cmax = 0) and the specific heat C1*** of moist air is constant capAir = MAir*(Prop(CpAir)+WAirEnt*Prop(CpVap)) ntu = UAExt/MAX(capAir,small) effectiveness = 1 - EXP(-ntu) C1*** Calculate coil surface enthalpy and temperature at the exit C1*** of the wet part of the coil using the effectiveness relation effectiveness = MAX(effectiveness,small) hCondSat = HAirEnt-(HAirEnt-HAirLvg)/effectiveness C1*** Calculate condensate temperature as the saturation temperature C1*** at given saturation enthalpy tempCond = TAIRSAT(Prop,hCondSat) C1*** Calculate exit air conditions and sensible heat transfer IF (tempCond .LT. DEWPOINT(Prop,WAirEnt)) THEN TAirLvg = TAirEnt-(TAirEnt-tempCond)*effectiveness WAirLvg = HUMTH(Prop,TAirLvg,HAirLvg) ELSE WAirLvg = WAirEnt TAirLvg = DRYBULB(Prop,HAirLvg,WAirLvg) ENDIF Qsen = capAir*(TAirEnt-TAirLvg) RETURN 1 END