SUBROUTINE TYPE70 (TIME,XIN,OUT,T,DTDT,PAR,INFO,ICNTRL,*) C*********************************************************************** C* Copyright ASHRAE. Toolkit for HVAC System Energy Calculations C*********************************************************************** C* C* ****************WARNING******************* C* THE SAMPLE DATA INCLUDED BELOW WAS GENERATED USING THE TRNSYS C* DEBUG PROGRAM. THEY WERE NOT INCLUDED WITH THE DOCUMENTATION!! C* C********************************************************************** C* C* SUBROUTINE: DRYWETCOIL C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: Calculate the performance of a cooling C* coil when the external fin surface is C* part wet and part dry. Results include C* outlet air temperature and humidity, C* outlet liquid temperature, sensible and C* total cooling capacities, and the wet C* fraction of the air-side surface area. C*********************************************************************** C* INPUT VARIABLES DESCRIPTION(UNITS) SAMPLE VALUES C* C* XIN(1) MLiq Liquid mass flow rate(kg/s) 1.60 C* XIN(2) TLiqEnt Entering liquid temperature(C) 6.50 C* XIN(3) MAir Dry air mass flow rate(kg/s) 2.887 C* XIN(4) TAirEnt Entering air dry bulb temperature(C) 25.60 C* XIN(5) HAirEnt Entering air enthalpy(J/kg) 16.0 C* XIN(6) WAirEnt Entering air humidity ratio(-) .008 C* XIN(7) TDewPt Entering air dew point(C) 10.5 C* XIN(8) ATot External (air-side) surface area(m2) 4.35 C* XIN(9) UIntTot Int. overall heat transfer coeff.(W/m2C) 1100.0 C* XIN(10) DryUExtTot Ext. overall heat transfer coeff.(W/m2C) 400.0 C* XIN(11) WetUExtTot Ext. overall heat transfer coeff.(W/m2C) 400.0 C* C* OUTPUT VARIABLES C* OUT(1) TLiqLvg Leaving liquid temperature(C) 9.277 C* OUT(2) TAirLvg Leaving air dry bulb temperature(C) 19.33 C* OUT(3) WAirLvg Leaving air humidity ratio(-) .007982 C* OUT(4) QTot Total heat transfer rate(W) 18600.0 C* OUT(5) QSen Sensible heat transfer rate(W) 18470.0 C* OUT(6) FWet Fraction of surface area wet(-) .1358 C* OUT(7) ErrStat Error status indicator,0=ok,1=error(-) 0.0 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: DRYCOIL C WETCOIL C FUNCTIONS CALLED: XITERATE 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 qDry Heat transfer rate for dry coil (W) C qTotWet Total heat transfer rate for wet coil (W) C qSenWet Sensible heat transfer rate for wet coil (W) C aWet Air-side area of wet coil (m2) C aDry Air-side area of dry coil (m2) C dryUA Overall heat transfer coefficient for dry coil (W/C) C wetUA Overall heat transfer coefficient for wet coil (W/C) C tLiqBnd Liquid temperature at wet/dry boundary (C) C tAirBnd Air temperature at wet/dry boundary (C) C tSurfBnd Surface temperature at wet/dry boundary (C) C tNewLiqBnd Estimated liquid temperature at wet/dry boundary (C) 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*********************************************************************** DOUBLE PRECISION XIN(11),OUT(7) INTEGER ErrStat,INFO(15),IOPT,NI,NP,ND REAL MAir,MLiq,PROP(16) CHARACTER*3 YCHECK(11),OCHECK(7) INTEGER Patm,CpAir,CpWat,CpLiq,CpVap,DViscAir, & DViscLiq,KAir,KLiq,RhoLiq,Hfg,RAir, & TKelMult,TAbsAdd,PaMult,PAbsAdd COMMON /LUNITS/LUR,LUW,IFORM,LUK DATA itmax/50/,configHX/1./ DATA YCHECK/'MF2','TE1','MF2','TE1','NAV','DM1','TE1','AR1', & 'NAV','NAV','NAV'/ DATA OCHECK/'TE1','TE1','DM1','PW2','PW2','DM1','DM1'/ 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 IOPT = -1 NI = 11 !CORRECT NUMBER OF INPUTS NP = 0 !CORRECT NUMBER OF PARAMETERS ND = 0 !CORRECT NUMBER OF DERIVATIVES MLIQ = XIN(1) TLIQENT = XIN(2) MAIR = XIN(3) TAIRENT = XIN(4) HAIRENT = XIN(5) WAIRENT = XIN(6) TDEWPT = XIN(7) ATOT = XIN(8) UINTTOT = XIN(9) DRYUEXTTOT = XIN(10) WETUEXTTOT = XIN(11) IF (INFO(7).EQ.-1) THEN CALL TYPECK(IOPT,INFO,NI,NP,ND) C CHECKS TO SEE IF USER'S INPUT MATCHES CORRECT NUMBERS CALL RCHECK(INFO,YCHECK,OCHECK) C CHECKS TO SEE IF INPUT AND OUTPUT UNITS MATCH ENDIF C1*** Iterate on FWet to converge on surface temperature equal to C1*** entering air dewpoint at wet/dry boundary. C1*** Preliminary estimates of coil performance to begin iteration TLiqLvg = TAirEnt qDry = 0. qTotWet = 0. qSenWet = 0. C2*** Estimate liquid temperature at boundary as entering air dew point tLiqBnd = TDewPt C2*** Estimate fraction wet surface area based on liquid temperatures FWet = (tLiqBnd-TLiqEnt)/(TLiqLvg-TLiqEnt) C1*** BEGIN LOOP to converge on FWet DO 100 iter = 1,itmax aWet = FWet*ATot aDry = ATot-aWet dryUA = aDry/(1./UIntTot+1./DryUExtTot) wetUAExt = WetUExtTot*aWet wetUAInt = UIntTot*aWet tLiqBnd = TLiqEnt+FWet*(TLiqLvg-TliqEnt) C1*** BEGIN LOOP to converge on liquid temperature at wet/dry boundary DO 50 itT = 1,itmax C1*** Calculate dry coil performance with estimated liquid temperature C1*** at the boundary. CALL DRYCOIL(Prop,MLiq,tLiqBnd,MAir,TAirEnt,WAirEnt, & dryUA,configHX, & TLiqLvg,tAirBnd,wAirBnd,qDry, & ErrStat) C1*** Calculate wet coil performance with calculated air temperature C1*** at the boundary. CALL WETCOIL (Prop,MLiq,TLiqEnt,MAir,tAirBnd,wAirBnd, & wetUAInt,wetUAExt,configHX, & tNewLiqBnd,TAirLvg,WAirLvg,qTotWet,qSenWet, & dumWet,tSurfBnd,ErrStat) errorT = tNewLiqBnd - tLiqBnd tLiqBnd = XITERATE(tLiqBnd,errorT,X1T,F1T,X2T,F2T,itT,icvgT) IF(icvgT .EQ. 1) GO TO 60 50 CONTINUE C2*** Boundary temperature not converged after maximum specified iterations. C2*** Print error message, set return error flag WRITE(LUW,1001) itmax 1001 FORMAT(/1X,'*** ERROR IN SUBROUTINE PARTWET ***'/ & 1X,' Liquid temperature not converged at boundary ' & 'after ',I2,' iterations'/) ErrStat = 1 C1*** Estimate new value for fraction wet surface area 60 CONTINUE C1*** If surface is dry, calculate dry coil performance and return IF(FWet .LE. 0.0 .AND. tSurfBnd .GE. TDewPt) THEN dryUA = aTot/(1./UIntTot+1./DryUExtTot) CALL DRYCOIL(Prop,MLiq,TLiqEnt,MAir,TAirEnt,WAirEnt, & dryUA,configHX, & TLiqLvg,TAirLvg,WAirLvg,QTot, & ErrStat) QSen = QTot FWet = 0. GO TO 999 ENDIF error = tSurfBnd - TDewPt FWet = XITERATE(FWet,error,X1,F1,X2,F2,iter,icvg) C1*** If converged, leave iteration loop IF (icvg .EQ. 1) GO TO 110 C2*** Surface temperature not converged. Repeat calculations with new C2*** estimate of fraction wet surface area. IF(FWet .GT. 1) FWet = 1. IF(FWet .LT. 0) FWet = 0. 100 CONTINUE C2*** Surface temperature not converged after maximum specified iterations. C2*** Print error message, set return error flag WRITE(LUW,1002) itmax 1002 FORMAT(/1X,'*** ERROR IN SUBROUTINE PARTWET ***'/ & 1X,' Wet/Dry boundary surface temperature not ' & 'converged after ',I2,' iterations'/) ErrStat = 1 110 CONTINUE C1*** Calculate outlet air temperature and humidity from enthalpies and C1*** surface conditions. QTot = qDry+qTotWet QSen = qDry+qSenWet 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 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,ERRSTAT REAL Prop(16),MAIR,MLIQ 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) 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,ERRSTAT REAL Prop(16),MAIR,MLIQ,INTRESIST 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/ 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 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/ 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*** 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./ 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*** 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 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 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 END