diff --git a/src/c/psychrolib.c b/src/c/psychrolib.c index 2060e561..f3b37df3 100644 --- a/src/c/psychrolib.c +++ b/src/c/psychrolib.c @@ -53,12 +53,16 @@ * Global constants *****************************************************************************************************/ -#define R_DA_IP 53.350 // Universal gas constant for dry air (IP version) in ft∙lbf/lb_da/R - // Reference: ASHRAE Handbook - Fundamentals (2017) ch. 1 -#define R_DA_SI 287.042 // Universal gas constant for dry air (SI version) in J/kg_da/K - // Reference: ASHRAE Handbook - Fundamentals (2017) ch. 1 +#define R_DA_IP 53.350 // Universal gas constant for dry air (IP version) in ft∙lbf/lb_da/R + // Reference: ASHRAE Handbook - Fundamentals (2017) ch. 1 +#define R_DA_SI 287.042 // Universal gas constant for dry air (SI version) in J/kg_da/K + // Reference: ASHRAE Handbook - Fundamentals (2017) ch. 1 +#define INVALID -99999 // Invalid value -#define INVALID -99999 // Invalid value +#define MAX_ITER_COUNT 100 // Maximum number of iterations before exiting while loops. + +#define MIN_HUM_RATIO 1e-7 // Minimum acceptable humidity ratio used/returned by any functions. + // Any value above 0 or below the MIN_HUM_RATIO will be reset to this value. /****************************************************************************************************** @@ -286,6 +290,41 @@ double GetRelHumFromVapPres // (o) Relative humidity [0-1] return VapPres/GetSatVapPres(TDryBulb); } +// Helper function returning the derivative of the natural log of the saturation vapor pressure +// as a function of dry-bulb temperature. +// Reference: ASHRAE Handbook - Fundamentals (2017) ch. 1 eqn. 5 & 6 +double dLnPws_ // (o) Derivative of natural log of vapor pressure of saturated air in Psi [IP] or Pa [SI] + ( double TDryBulb // (i) Dry bulb temperature in °F [IP] or °C [SI] + ) +{ + double dLnPws, T; + + if (isIP()) + { + T = GetTRankineFromTFahrenheit(TDryBulb); + + if (TDryBulb <= 32.) + dLnPws = 1.0214165E+04 / pow(T, 2) - 5.3765794E-03 + 2 * 1.9202377E-07 * T + + 2 * 3.5575832E-10 * pow(T, 2) - 4 * 9.0344688E-14 * pow(T, 3) + 4.1635019 / T; + else + dLnPws = 1.0440397E+04 / pow(T, 2) - 2.7022355E-02 + 2 * 1.2890360E-05 * T + - 3 * 2.4780681E-09 * pow(T, 2) + 6.5459673 / T; + } + else + { + T = GetTKelvinFromTCelsius(TDryBulb); + + if (TDryBulb <= 0.) + dLnPws = 5.6745359E+03 / pow(T, 2) - 9.677843E-03 + 2 * 6.2215701E-07 * T + + 3 * 2.0747825E-09 * pow(T, 2) - 4 * 9.484024E-13 * pow(T, 3) + 4.1635019 / T; + else + dLnPws = 5.8002206E+03 / pow(T, 2) - 4.8640239E-02 + 2 * 4.1764768E-05 * T + - 3 * 1.4452093E-08 * pow(T, 2) + 6.5459673 / T; + } + + return dLnPws; +} + // Return dew-point temperature given dry-bulb temperature and vapor pressure. // Reference: ASHRAE Handbook - Fundamentals (2017) ch. 1 eqn. 5 and 6 // Notes: the dew point temperature is solved by inverting the equation giving water vapor pressure @@ -301,55 +340,71 @@ double GetTDewPointFromVapPres // (o) Dew Point temperature in °F [IP] or °C , double VapPres // (i) Partial pressure of water vapor in moist air in Psi [IP] or Pa [SI] ) { - // Bounds and step size as a function of the system of units - double _BOUNDS[2]; // Domain of validity of the equations - double _STEPSIZE; // Temperature step used for the calculation of numerical derivatives + // Bounds function of the system of units + double BOUNDS[2]; // Domain of validity of the equations + double T_WATER_FREEZE, T_WATER_FREEZE_LOW, T_WATER_FREEZE_HIGH, PWS_FREEZE_LOW, PWS_FREEZE_HIGH; + if (isIP()) { - _BOUNDS[0] = -148.; - _BOUNDS[1] = 392.; - _STEPSIZE = 0.01 * 9. / 5.; + BOUNDS[0] = -148.; + BOUNDS[1] = 392.; + T_WATER_FREEZE = 32.; } else { - _BOUNDS[0] = -100.; - _BOUNDS[1] = 200.; - _STEPSIZE = 0.01; + BOUNDS[0] = -100.; + BOUNDS[1] = 200.; + T_WATER_FREEZE = 0.; } - double TMidPoint = (_BOUNDS[0] + _BOUNDS[1]) / 2.; // Midpoint of domain of validity - // Bounds outside which a solution cannot be found - ASSERT (VapPres >= GetSatVapPres(_BOUNDS[0]) && VapPres <= GetSatVapPres(_BOUNDS[1]), "Partial pressure of water vapor is outside range of validity of equations") + ASSERT (VapPres >= GetSatVapPres(BOUNDS[0]) && VapPres <= GetSatVapPres(BOUNDS[1]), "Partial pressure of water vapor is outside range of validity of equations") + + // Vapor pressure contained within the discontinuity of the Pws function: return temperature of freezing + T_WATER_FREEZE_LOW = T_WATER_FREEZE - PSYCHROLIB_TOLERANCE / 10.; // Temperature just below freezing + T_WATER_FREEZE_HIGH = T_WATER_FREEZE + PSYCHROLIB_TOLERANCE / 10.; // Temperature just above freezing + PWS_FREEZE_LOW = GetSatVapPres(T_WATER_FREEZE_LOW); + PWS_FREEZE_HIGH = GetSatVapPres(T_WATER_FREEZE_HIGH); + + // Restrict iteration to either left or right part of the saturation vapor pressure curve + // to avoid iterating back and forth across the discontinuity of the curve at the freezing point + // When the partial pressure of water vapor is within the discontinuity of GetSatVapPres, + // simply return the freezing point of water. + if (VapPres < PWS_FREEZE_LOW) + BOUNDS[1] = T_WATER_FREEZE_LOW; + else if (VapPres > PWS_FREEZE_LOW) + BOUNDS[0] = T_WATER_FREEZE_HIGH; + else + return T_WATER_FREEZE; + // We use NR to approximate the solution. // First guess - double Tdp = TDryBulb; // Calculated value of dew point temperatures, solved for iteratively in °F [IP] or °C [SI] - double lnVP = log(VapPres); // Natural logarithm of partial pressure of water vapor pressure in moist air + double TDewPoint = TDryBulb; // Calculated value of dew point temperatures, solved for iteratively in °F [IP] or °C [SI] + double lnVP = log(VapPres); // Natural logarithm of partial pressure of water vapor pressure in moist air + + double TDewPoint_iter; // Value of TDewPoint used in NR calculation + double lnVP_iter; // Value of log of vapor water pressure used in NR calculation + int index = 1; - double Tdp_c; // Value of Tdp used in NR calculation - double lnVP_c; // Value of log of vapor water pressure used in NR calculation - double d_Tdp; // Value of temperature step used in NR calculation do { - // Current point - Tdp_c = Tdp; - lnVP_c = log(GetSatVapPres(Tdp_c)); - // Step - negative in the right part of the curve, positive in the left part - // to avoid going past the domain of validity of eqn. 5 and 6 - // when Tdp_c is close to its bounds - if (Tdp_c > TMidPoint) - d_Tdp = -_STEPSIZE; - else - d_Tdp = _STEPSIZE; - // Derivative of function, calculated numerically - double d_lnVP = (log(GetSatVapPres(Tdp_c + d_Tdp)) - lnVP_c) / d_Tdp; + TDewPoint_iter = TDewPoint; // TDewPoint used in NR calculation + lnVP_iter = log(GetSatVapPres(TDewPoint_iter)); + + // Derivative of function, calculated analytically + double d_lnVP = dLnPws_(TDewPoint_iter); + // New estimate, bounded by domain of validity of eqn. 5 and 6 - Tdp = Tdp_c - (lnVP_c - lnVP) / d_lnVP; - Tdp = max(Tdp, _BOUNDS[0]); - Tdp = min(Tdp, _BOUNDS[1]); + TDewPoint = TDewPoint_iter - (lnVP_iter - lnVP) / d_lnVP; + TDewPoint = max(TDewPoint, BOUNDS[0]); + TDewPoint = min(TDewPoint, BOUNDS[1]); + + ASSERT (index <= MAX_ITER_COUNT, "Convergence not reached in GetTDewPointFromVapPres. Stopping.") + + index++; } - while (fabs(Tdp - Tdp_c) > PSYCHROLIB_TOLERANCE); - return min(Tdp, TDryBulb); + while (fabs(TDewPoint - TDewPoint_iter) > PSYCHROLIB_TOLERANCE); + return min(TDewPoint, TDryBulb); } // Return vapor pressure given dew point temperature. @@ -376,11 +431,13 @@ double GetTWetBulbFromHumRatio // (o) Wet bulb temperature in °F [IP] or °C [ { // Declarations double Wstar; - double TDewPoint, TWetBulb, TWetBulbSup, TWetBulbInf; + double TDewPoint, TWetBulb, TWetBulbSup, TWetBulbInf, BoundedHumRatio; + int index = 1; ASSERT (HumRatio >= 0., "Humidity ratio is negative") + BoundedHumRatio = max(HumRatio, MIN_HUM_RATIO); - TDewPoint = GetTDewPointFromHumRatio(TDryBulb, HumRatio, Pressure); + TDewPoint = GetTDewPointFromHumRatio(TDryBulb, BoundedHumRatio, Pressure); // Initial guesses TWetBulbSup = TDryBulb; @@ -388,19 +445,23 @@ double GetTWetBulbFromHumRatio // (o) Wet bulb temperature in °F [IP] or °C [ TWetBulb = (TWetBulbInf + TWetBulbSup) / 2.; // Bisection loop - while(TWetBulbSup - TWetBulbInf > PSYCHROLIB_TOLERANCE) + while ((TWetBulbSup - TWetBulbInf) > PSYCHROLIB_TOLERANCE) { // Compute humidity ratio at temperature Tstar Wstar = GetHumRatioFromTWetBulb(TDryBulb, TWetBulb, Pressure); // Get new bounds - if (Wstar > HumRatio) + if (Wstar > BoundedHumRatio) TWetBulbSup = TWetBulb; else TWetBulbInf = TWetBulb; // New guess of wet bulb temperature TWetBulb = (TWetBulbSup+TWetBulbInf) / 2.; + + ASSERT (index <= MAX_ITER_COUNT, "Convergence not reached in GetTWetBulbFromHumRatio. Stopping.") + + index++; } return TWetBulb; @@ -439,8 +500,8 @@ double GetHumRatioFromTWetBulb // (o) Humidity Ratio in lb_H₂O lb_Air⁻¹ [I HumRatio = ((2830. - 0.24 * TWetBulb) * Wsstar - 1.006 * (TDryBulb - TWetBulb)) / (2830. + 1.86 * TDryBulb - 2.1 * TWetBulb); } - - return HumRatio; + // Validity check. + return max(HumRatio, MIN_HUM_RATIO); } // Return humidity ratio given dry-bulb temperature, relative humidity, and pressure. @@ -516,9 +577,14 @@ double GetHumRatioFromVapPres // (o) Humidity Ratio in lb_H₂O lb_Air⁻¹ [I , double Pressure // (i) Atmospheric pressure in Psi [IP] or Pa [SI] ) { + double HumRatio; + ASSERT (VapPres >= 0., "Partial pressure of water vapor in moist air is negative") - return 0.621945 * VapPres / (Pressure - VapPres); + HumRatio = 0.621945 * VapPres / (Pressure - VapPres); + + // Validity check. + return max(HumRatio, MIN_HUM_RATIO); } // Return vapor pressure given humidity ratio and pressure. @@ -528,11 +594,12 @@ double GetVapPresFromHumRatio // (o) Partial pressure of water vapor in moist , double Pressure // (i) Atmospheric pressure in Psi [IP] or Pa [SI] ) { - double VapPres; + double VapPres, BoundedHumRatio; ASSERT (HumRatio >= 0., "Humidity ratio is negative") + BoundedHumRatio = max(HumRatio, MIN_HUM_RATIO); - VapPres = Pressure * HumRatio/(0.621945 + HumRatio); + VapPres = Pressure * BoundedHumRatio / (0.621945 + BoundedHumRatio); return VapPres; } @@ -547,8 +614,12 @@ double GetSpecificHumFromHumRatio // (o) Specific humidity ratio in lb_H₂O lb_ ( double HumRatio // (i) Humidity ratio in lb_H₂O lb_Dry_Air⁻¹ [IP] or kg_H₂O kg_Dry_Air⁻¹ [SI] ) { + double BoundedHumRatio; + ASSERT (HumRatio >= 0., "Humidity ratio is negative") - return HumRatio / (1.0 + HumRatio); + BoundedHumRatio = max(HumRatio, MIN_HUM_RATIO); + + return BoundedHumRatio / (1.0 + BoundedHumRatio); } // Return the humidity ratio (aka mixing ratio) from specific humidity @@ -557,8 +628,14 @@ double GetHumRatioFromSpecificHum // (o) Humidity ratio in lb_H₂O lb_Dry_Air ( double SpecificHum // (i) Specific humidity ratio in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI] ) { + double HumRatio; + ASSERT (SpecificHum >= 0.0 && SpecificHum < 1.0, "Specific humidity is outside range [0,1[") - return SpecificHum / (1.0 - SpecificHum); + + HumRatio = SpecificHum / (1.0 - SpecificHum); + + // Validity check + return max(HumRatio, MIN_HUM_RATIO); } @@ -618,12 +695,15 @@ double GetTDryBulbFromEnthalpyAndHumRatio // (o) Dry-bulb temperature in °F [I , double HumRatio // (i) Humidity ratio in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI] ) { + double BoundedHumRatio; + ASSERT (HumRatio >= 0., "Humidity ratio is negative") + BoundedHumRatio = max(HumRatio, MIN_HUM_RATIO); if (isIP()) - return (MoistAirEnthalpy - 1061.0 * HumRatio) / (0.240 + 0.444 * HumRatio); + return (MoistAirEnthalpy - 1061.0 * BoundedHumRatio) / (0.240 + 0.444 * BoundedHumRatio); else - return (MoistAirEnthalpy / 1000.0 - 2501.0 * HumRatio) / (1.006 + 1.86 * HumRatio); + return (MoistAirEnthalpy / 1000.0 - 2501.0 * BoundedHumRatio) / (1.006 + 1.86 * BoundedHumRatio); } // Return humidity ratio from enthalpy and dry-bulb temperature. @@ -634,10 +714,14 @@ double GetHumRatioFromEnthalpyAndTDryBulb // (o) Humidity ratio in lb_H₂O lb_ , double TDryBulb // (i) Dry-bulb temperature in °F [IP] or °C [SI] ) { + double HumRatio; if (isIP()) - return (MoistAirEnthalpy - 0.240 * TDryBulb) / (1061.0 + 0.444 * TDryBulb); + HumRatio = (MoistAirEnthalpy - 0.240 * TDryBulb) / (1061.0 + 0.444 * TDryBulb); else - return (MoistAirEnthalpy / 1000.0 - 1.006 * TDryBulb) / (2501.0 + 1.86 * TDryBulb); + HumRatio = (MoistAirEnthalpy / 1000.0 - 1.006 * TDryBulb) / (2501.0 + 1.86 * TDryBulb); + + // Validity check. + return max(HumRatio, MIN_HUM_RATIO); } @@ -692,10 +776,13 @@ double GetSatHumRatio // (o) Humidity ratio of saturated air in lb_H , double Pressure // (i) Atmospheric pressure in Psi [IP] or Pa [SI] ) { - double SatVaporPres; + double SatVaporPres, SatHumRatio; SatVaporPres = GetSatVapPres(TDryBulb); - return 0.621945 * SatVaporPres / (Pressure - SatVaporPres); + SatHumRatio = 0.621945 * SatVaporPres / (Pressure - SatVaporPres); + + // Validity check. + return max(SatHumRatio, MIN_HUM_RATIO); } // Return saturated air enthalpy given dry-bulb temperature and pressure. @@ -738,9 +825,12 @@ double GetDegreeOfSaturation // (o) Degree of saturation [] , double Pressure // (i) Atmospheric pressure in Psi [IP] or Pa [SI] ) { + double BoundedHumRatio; + ASSERT (HumRatio >= 0., "Humidity ratio is negative") + BoundedHumRatio = max(HumRatio, MIN_HUM_RATIO); - return HumRatio / GetSatHumRatio(TDryBulb, Pressure); + return BoundedHumRatio / GetSatHumRatio(TDryBulb, Pressure); } // Return moist air enthalpy given dry-bulb temperature and humidity ratio. @@ -750,12 +840,15 @@ double GetMoistAirEnthalpy // (o) Moist Air Enthalpy in Btu lb⁻¹ [IP] or , double HumRatio // (i) Humidity ratio in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI] ) { + double BoundedHumRatio; + ASSERT (HumRatio >= 0., "Humidity ratio is negative") + BoundedHumRatio = max(HumRatio, MIN_HUM_RATIO); if (isIP()) - return 0.240 * TDryBulb + HumRatio*(1061. + 0.444 * TDryBulb); + return 0.240 * TDryBulb + BoundedHumRatio*(1061. + 0.444 * TDryBulb); else - return (1.006 * TDryBulb + HumRatio*(2501. + 1.86 * TDryBulb)) * 1000.; + return (1.006 * TDryBulb + BoundedHumRatio*(2501. + 1.86 * TDryBulb)) * 1000.; } // Return moist air specific volume given dry-bulb temperature, humidity ratio, and pressure. @@ -768,12 +861,15 @@ double GetMoistAirVolume // (o) Specific Volume ft³ lb⁻¹ [IP] or in m , double Pressure // (i) Atmospheric pressure in Psi [IP] or Pa [SI] ) { - ASSERT(HumRatio >= 0., "Humidity ratio is negative") + double BoundedHumRatio; + + ASSERT (HumRatio >= 0., "Humidity ratio is negative") + BoundedHumRatio = max(HumRatio, MIN_HUM_RATIO); if (isIP()) - return R_DA_IP * GetTRankineFromTFahrenheit(TDryBulb) * (1. + 1.607858*HumRatio) / (144. * Pressure); + return R_DA_IP * GetTRankineFromTFahrenheit(TDryBulb) * (1. + 1.607858 * BoundedHumRatio) / (144. * Pressure); else - return R_DA_SI * GetTKelvinFromTCelsius(TDryBulb) * (1. + 1.607858*HumRatio) / Pressure; + return R_DA_SI * GetTKelvinFromTCelsius(TDryBulb) * (1. + 1.607858 * BoundedHumRatio) / Pressure; } // Return moist air density given humidity ratio, dry bulb temperature, and pressure. @@ -784,9 +880,12 @@ double GetMoistAirDensity // (o) Moist air density in lb ft⁻³ [IP] or k , double Pressure // (i) Atmospheric pressure in Psi [IP] or Pa [SI] ) { + double BoundedHumRatio; + ASSERT (HumRatio >= 0., "Humidity ratio is negative") + BoundedHumRatio = max(HumRatio, MIN_HUM_RATIO); - return (1. + HumRatio) / GetMoistAirVolume(TDryBulb, HumRatio, Pressure); + return (1. + BoundedHumRatio) / GetMoistAirVolume(TDryBulb, BoundedHumRatio, Pressure); } diff --git a/src/fortran/psychrolib.f90 b/src/fortran/psychrolib.f90 index f7293bc0..8b6bf88d 100644 --- a/src/fortran/psychrolib.f90 +++ b/src/fortran/psychrolib.f90 @@ -88,6 +88,7 @@ module psychrolib public :: CalcPsychrometricsFromTWetBulb public :: CalcPsychrometricsFromTDewPoint public :: CalcPsychrometricsFromRelHum + public :: dLnPws_ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -111,7 +112,14 @@ module psychrolib !+ Unit system to use. real :: PSYCHROLIB_TOLERANCE = 1.0 - !+ Tolerance of temperature calculations + !+ Tolerance of temperature calculations. + + integer, parameter :: MAX_ITER_COUNT = 100 + !+ Maximum number of iterations before exiting while loops. + + real, parameter :: MIN_HUM_RATIO = 1e-7 + !+ Minimum acceptable humidity ratio used/returned by any functions. + !+ Any value above 0 or below the MIN_HUM_RATIO will be reset to this value. contains @@ -383,6 +391,45 @@ function GetRelHumFromVapPres(TDryBulb, VapPres) result(RelHum) RelHum = VapPres / GetSatVapPres(TDryBulb) end function GetRelHumFromVapPres + function dLnPws_(TDryBulb) result(dLnPws) + !+ Helper function returning the derivative of the natural log of the saturation vapor pressure + !+ as a function of dry-bulb temperature. + !+ Reference: + !+ ASHRAE Handbook - Fundamentals (2017) ch. 1 eqn 5 + + real, intent(in) :: TDryBulb + !+ Dry-bulb temperature in °F [IP] or °C [SI] + real :: dLnPws + !+ Derivative of natural log of vapor pressure of saturated air in Psi [IP] or Pa [SI] + real :: T + !+ Dry bulb temperature in R [IP] or K [SI] + + if (isIP()) then + + T = GetTRankineFromTFahrenheit(TDryBulb) + + if (TDryBulb <= 32.) then + dLnPws = 1.0214165E+04 / T**2 - 5.3765794E-03 + 2 * 1.9202377E-07 * T & + + 2 * 3.5575832E-10 * T**2 - 4 * 9.0344688E-14 * T**3 + 4.1635019 / T + else + dLnPws = 1.0440397E+04 / T**2 - 2.7022355E-02 + 2 * 1.2890360E-05 * T & + - 3 * 2.4780681E-09 * T**2 + 6.5459673 / T + end if + + else + + T = GetTKelvinFromTCelsius(TDryBulb) + + if (TDryBulb <= 0) then + dLnPws = 5.6745359E+03 / T**2 - 9.677843E-03 + 2 * 6.2215701E-07 * T & + + 3 * 2.0747825E-09 * T**2 - 4 * 9.484024E-13 * T**3 + 4.1635019 / T + else + dLnPws = 5.8002206E+03 / T**2 - 4.8640239E-02 + 2 * 4.1764768E-05 * T & + - 3 * 1.4452093E-08 * T**2 + 6.5459673 / T + end if + end if + end function dLnPws_ + function GetTDewPointFromVapPres(TDryBulb, VapPres) result(TDewPoint) !+ Return dew-point temperature given dry-bulb temperature and vapor pressure. !+ References: @@ -403,73 +450,85 @@ function GetTDewPointFromVapPres(TDryBulb, VapPres) result(TDewPoint) !+ Partial pressure of water vapor in moist air in Psi [IP] or Pa [SI] real :: TDewPoint !+ Dew-point temperature in °F [IP] or °C [SI] - real :: Tdp - !+ Calculated value of dew point temperatures, solved for iteratively in °F [IP] or °C [SI] real :: lnVP !+ Natural logarithm of partial pressure of water vapor pressure in moist air real :: d_lnVP !+ Derivative of function, calculated numerically - real :: lnVP_c + real :: lnVP_iter !+ Value of log of vapor water pressure used in NR calculation - real :: Tdp_c - !+ Value of Tdp used in NR calculation - real :: d_Tdp - !+ Value of temperature step used in NR calculation - real :: STEPSIZE - !+ Size of timestep (dimensionless) - real :: TMidPoint - !+ Midpoint of domain of validity + real :: TDewPoint_iter + !+ Value of TDewPoint used in NR calculation real, dimension(2) :: BOUNDS !+ Valid temperature range in °F [IP] or °C [SI] integer :: index !+ Index used in the calculation + real :: T_WATER_FREEZE, T_WATER_FREEZE_LOW, T_WATER_FREEZE_HIGH, PWS_FREEZE_LOW, PWS_FREEZE_HIGH ! Bounds and step size as a function of the system of units if (isIP()) then BOUNDS(1) = -148.0 BOUNDS(2) = 392.0 - STEPSIZE = 0.01 * 9.0 / 5.0 + T_WATER_FREEZE = 32. else BOUNDS(1) = -100.0 BOUNDS(2) = 200.0 - STEPSIZE = 0.01 + T_WATER_FREEZE = 0. end if - TMidPoint = (BOUNDS(2) + BOUNDS(2)) / 2.0 - - ! Bounds outside which a solution cannot be found + ! Validity check -- bounds outside which a solution cannot be found if (VapPres < GetSatVapPres(BOUNDS(1)) .or. VapPres > GetSatVapPres(BOUNDS(2))) then error stop "Error: partial pressure of water vapor is outside range of validity of equations" end if - ! First guess - Tdp = TDryBulb + ! Vapor pressure contained within the discontinuity of the Pws function: return temperature of freezing + T_WATER_FREEZE_LOW = T_WATER_FREEZE - PSYCHROLIB_TOLERANCE / 10. ! Temperature just below freezing + T_WATER_FREEZE_LOW = T_WATER_FREEZE - PSYCHROLIB_TOLERANCE / 10. ! Temperature just below freezing + T_WATER_FREEZE_HIGH = T_WATER_FREEZE + PSYCHROLIB_TOLERANCE ! Temperature just above freezing + PWS_FREEZE_LOW = GetSatVapPres(T_WATER_FREEZE_LOW) + PWS_FREEZE_HIGH = GetSatVapPres(T_WATER_FREEZE_HIGH) + + ! Restrict iteration to either left or right part of the saturation vapor pressure curve + ! to avoid iterating back and forth across the discontinuity of the curve at the freezing point + ! When the partial pressure of water vapor is within the discontinuity of GetSatVapPres, + ! simply return the freezing point of water. + if (VapPres < PWS_FREEZE_LOW) then + BOUNDS(2) = T_WATER_FREEZE_LOW + else if (VapPres > PWS_FREEZE_HIGH) then + BOUNDS(1) = T_WATER_FREEZE_HIGH + else + TDewPoint = T_WATER_FREEZE + return + end if + + ! We use NR to approximate the solution. + TDewPoint = TDryBulb lnVP = log(VapPres) - Tdp_c = Tdp index = 1 - do while ((abs(Tdp - Tdp_c) > PSYCHROLIB_TOLERANCE) .or. (index < 2)) - ! Current point - Tdp_c = Tdp - lnVP_c = log(GetSatVapPres(Tdp_c)) - ! Step - negative in the right part of the curve, positive in the left part - ! to avoid going past the domain of validity of eqn. 5 and 6 - ! when Tdp_c is close to its bounds - if (Tdp_c > TMidPoint) then - d_Tdp = -STEPSIZE - else - d_Tdp = STEPSIZE + do while (.true.) + TDewPoint_iter = TDewPoint ! TDewPoint_iter used in NR calculation + lnVP_iter = log(GetSatVapPres(TDewPoint_iter)) + + ! Derivative of function, calculated analytically + d_lnVP = dLnPws_(TDewPoint_iter) + + ! New estimate, bounded by the search domain defined above + TDewPoint = TDewPoint_iter - (lnVP_iter - lnVP) / d_lnVP + TDewPoint = max(TDewPoint, BOUNDS(1)) + TDewPoint = min(TDewPoint, BOUNDS(2)) + + if (abs(TDewPoint - TDewPoint_iter) <= PSYCHROLIB_TOLERANCE) then + exit + end if + + if (index > MAX_ITER_COUNT) then + error stop "Convergence not reached in GetTDewPointFromVapPres. Stopping." end if - ! Derivative of function, calculated numerically - d_lnVP = (log(GetSatVapPres(Tdp_c + d_Tdp)) - lnVP_c) / d_Tdp - ! New estimate, bounded by domain of validity of eqn. 5 and 6 - Tdp = Tdp_c - (lnVP_c - lnVP) / d_lnVP - Tdp = max(Tdp, BOUNDS(1)) - Tdp = min(Tdp, BOUNDS(2)) + index = index + 1 end do - TDewPoint = min(Tdp, TDryBulb) + TDewPoint = min(TDewPoint, TDryBulb) end function GetTDewPointFromVapPres function GetVapPresFromTDewPoint(TDewPoint) result(VapPres) @@ -511,26 +570,32 @@ function GetTWetBulbFromHumRatio(TDryBulb, HumRatio, Pressure) result(TWetBulb) !+ Lower value of wet bulb temperature in bissection method (initial guess is from dew point temperature) in °F [IP] or °C [SI] real :: Wstar !+ Humidity ratio at temperature Tstar in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI] + real :: BoundedHumRatio + !+ Humidity ratio bounded to MIN_HUM_RATIO + integer :: index + !+ index used in iteration if (HumRatio < 0.0) then error stop "Error: humidity ratio cannot be negative" end if + BoundedHumRatio = max(HumRatio, MIN_HUM_RATIO) - TDewPoint = GetTDewPointFromHumRatio(TDryBulb, HumRatio, Pressure) + TDewPoint = GetTDewPointFromHumRatio(TDryBulb, BoundedHumRatio, Pressure) ! Initial guesses TWetBulbSup = TDryBulb TWetBulbInf = TDewPoint TWetBulb = (TWetBulbInf + TWetBulbSup) / 2.0 + index = 1 ! Bisection loop - do while(TWetBulbSup - TWetBulbInf > PSYCHROLIB_TOLERANCE) + do while ((TWetBulbSup - TWetBulbInf) > PSYCHROLIB_TOLERANCE) ! Compute humidity ratio at temperature Tstar Wstar = GetHumRatioFromTWetBulb(TDryBulb, TWetBulb, Pressure) ! Get new bounds - if (Wstar > HumRatio) then + if (Wstar > BoundedHumRatio) then TWetBulbSup = TWetBulb else TWetBulbInf = TWetBulb @@ -538,6 +603,12 @@ function GetTWetBulbFromHumRatio(TDryBulb, HumRatio, Pressure) result(TWetBulb) ! New guess of wet bulb temperature TWetBulb = (TWetBulbSup + TWetBulbInf) / 2.0 + + if (index > MAX_ITER_COUNT) then + error stop "Convergence not reached in GetTWetBulbFromHumRatio. Stopping." + end if + + index = index + 1 end do end function GetTWetBulbFromHumRatio @@ -580,6 +651,9 @@ function GetHumRatioFromTWetBulb(TDryBulb, TWetBulb, Pressure) result(HumRatio) / (2830.0 + 1.86 * TDryBulb - 2.1 * TWetBulb) end if end if + + ! Validity check. + HumRatio = max(HumRatio, MIN_HUM_RATIO) end function GetHumRatioFromTWetBulb function GetHumRatioFromRelHum(TDryBulb, RelHum, Pressure) result(HumRatio) @@ -695,6 +769,9 @@ function GetHumRatioFromVapPres(VapPres, Pressure) result(HumRatio) end if HumRatio = 0.621945 * VapPres / (Pressure-VapPres) + + ! Validity check. + HumRatio = max(HumRatio, MIN_HUM_RATIO) end function GetHumRatioFromVapPres function GetVapPresFromHumRatio(HumRatio, Pressure) result(VapPres) @@ -708,12 +785,15 @@ function GetVapPresFromHumRatio(HumRatio, Pressure) result(VapPres) !+ Atmospheric pressure in Psi [IP] or Pa [SI] real :: VapPres !+ Partial pressure of water vapor in moist air in Psi [IP] or Pa [SI] + real :: BoundedHumRatio + !+ Humidity ratio bounded to MIN_HUM_RATIO if (HumRatio < 0.0) then error stop "Error: humidity ratio is negative" end if + BoundedHumRatio = max(HumRatio, MIN_HUM_RATIO) - VapPres = Pressure * HumRatio / (0.621945 + HumRatio) + VapPres = Pressure * BoundedHumRatio / (0.621945 + BoundedHumRatio) end function GetVapPresFromHumRatio @@ -730,12 +810,15 @@ function GetSpecificHumFromHumRatio(HumRatio) result(SpecificHum) !+ Humidity ratio in lb_H₂O lb_Dry_Air⁻¹ [IP] or kg_H₂O kg_Dry_Air⁻¹ [SI] real :: SpecificHum !+ Specific humidity in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI] + real :: BoundedHumRatio + !+ Humidity ratio bounded to MIN_HUM_RATIO if (HumRatio < 0.0) then error stop "Error: humidity ratio cannot be negative" end if + BoundedHumRatio = max(HumRatio, MIN_HUM_RATIO) - SpecificHum = HumRatio / (1.0 + HumRatio) + SpecificHum = BoundedHumRatio / (1.0 + BoundedHumRatio) end function GetSpecificHumFromHumRatio function GetHumRatioFromSpecificHum(SpecificHum) result(HumRatio) @@ -753,6 +836,9 @@ function GetHumRatioFromSpecificHum(SpecificHum) result(HumRatio) end if HumRatio = SpecificHum / (1.0 - SpecificHum) + + ! Validity check. + HumRatio = max(HumRatio, MIN_HUM_RATIO) end function GetHumRatioFromSpecificHum @@ -836,15 +922,18 @@ function GetTDryBulbFromEnthalpyAndHumRatio(MoistAirEnthalpy, HumRatio) result(T !+ Humidity ratio in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI] real :: TDryBulb !+ Dry-bulb temperature in °F [IP] or °C [SI] + real :: BoundedHumRatio + !+ Humidity ratio bounded to MIN_HUM_RATIO if (HumRatio < 0.0) then error stop "Error: humidity ratio is negative" end if + BoundedHumRatio = max(HumRatio, MIN_HUM_RATIO) if (isIP()) then - TDryBulb = (MoistAirEnthalpy - 1061.0 * HumRatio) / (0.240 + 0.444 * HumRatio) + TDryBulb = (MoistAirEnthalpy - 1061.0 * BoundedHumRatio) / (0.240 + 0.444 * BoundedHumRatio) else - TDryBulb = (MoistAirEnthalpy / 1000.0 - 2501.0 * HumRatio) / (1.006 + 1.86 * HumRatio) + TDryBulb = (MoistAirEnthalpy / 1000.0 - 2501.0 * BoundedHumRatio) / (1.006 + 1.86 * BoundedHumRatio) end if end function GetTDryBulbFromEnthalpyAndHumRatio @@ -867,6 +956,9 @@ function GetHumRatioFromEnthalpyAndTDryBulb(MoistAirEnthalpy, TDryBulb) result(H else HumRatio = (MoistAirEnthalpy / 1000.0 - 1.006 * TDryBulb) / (2501.0 + 1.86 * TDryBulb) end if + + ! Validity check. + HumRatio = max(HumRatio, MIN_HUM_RATIO) end function GetHumRatioFromEnthalpyAndTDryBulb @@ -938,6 +1030,9 @@ function GetSatHumRatio(TDryBulb, Pressure) result(SatHumRatio) SatVaporPres = GetSatVapPres(TDryBulb) SatHumRatio = 0.621945 * SatVaporPres / (Pressure-SatVaporPres) + + ! Validity check. + SatHumRatio = max(SatHumRatio, MIN_HUM_RATIO) end function GetSatHumRatio function GetSatAirEnthalpy(TDryBulb, Pressure) result(SatAirEnthalpy) @@ -1000,12 +1095,15 @@ function GetDegreeOfSaturation(TDryBulb, HumRatio, Pressure) result(DegreeOfSatu !+ Atmospheric pressure in Psi [IP] or Pa [SI] real :: DegreeOfSaturation !+ Degree of saturation in arbitrary unit + real :: BoundedHumRatio + !+ Humidity ratio bounded to MIN_HUM_RATIO if (HumRatio < 0.0) then error stop "Error: humidity ratio is negative" end if + BoundedHumRatio = max(HumRatio, MIN_HUM_RATIO) - DegreeOfSaturation = HumRatio/GetSatHumRatio(TDryBulb, Pressure) + DegreeOfSaturation = BoundedHumRatio / GetSatHumRatio(TDryBulb, Pressure) end function GetDegreeOfSaturation function GetMoistAirEnthalpy(TDryBulb, HumRatio) result(MoistAirEnthalpy) @@ -1019,15 +1117,18 @@ function GetMoistAirEnthalpy(TDryBulb, HumRatio) result(MoistAirEnthalpy) !+ Humidity ratio in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI] real :: MoistAirEnthalpy !+ Moist air enthalpy in Btu lb⁻¹ [IP] or J kg⁻¹ + real :: BoundedHumRatio + !+ Humidity ratio bounded to MIN_HUM_RATIO if (HumRatio < 0.0) then error stop "Error: humidity ratio is negative" end if + BoundedHumRatio = max(HumRatio, MIN_HUM_RATIO) if (isIP()) then - MoistAirEnthalpy = 0.240 * TDryBulb + HumRatio * (1061.0 + 0.444 * TDryBulb) + MoistAirEnthalpy = 0.240 * TDryBulb + BoundedHumRatio * (1061.0 + 0.444 * TDryBulb) else - MoistAirEnthalpy = (1.006 * TDryBulb + HumRatio * (2501.0 + 1.86 * TDryBulb)) * 1000.0 + MoistAirEnthalpy = (1.006 * TDryBulb + BoundedHumRatio * (2501.0 + 1.86 * TDryBulb)) * 1000.0 end if end function GetMoistAirEnthalpy @@ -1047,15 +1148,18 @@ function GetMoistAirVolume(TDryBulb, HumRatio, Pressure) result(MoistAirVolume) !+ Atmospheric pressure in Psi [IP] or Pa [SI] real :: MoistAirVolume !+ Specific volume of moist air in ft³ lb⁻¹ of dry air [IP] or in m³ kg⁻¹ of dry air [SI] + real :: BoundedHumRatio + !+ Humidity ratio bounded to MIN_HUM_RATIO if (HumRatio < 0.0) then error stop "Error: humidity ratio is negative" end if + BoundedHumRatio = max(HumRatio, MIN_HUM_RATIO) if (isIP()) then - MoistAirVolume = R_DA_IP * GetTRankineFromTFahrenheit(TDryBulb) * (1.0 + 1.607858 * HumRatio) / (144.0 * Pressure) + MoistAirVolume = R_DA_IP * GetTRankineFromTFahrenheit(TDryBulb) * (1.0 + 1.607858 * BoundedHumRatio) / (144.0 * Pressure) else - MoistAirVolume = R_DA_SI * GetTKelvinFromTCelsius(TDryBulb) * (1.0 + 1.607858 * HumRatio) / Pressure + MoistAirVolume = R_DA_SI * GetTKelvinFromTCelsius(TDryBulb) * (1.0 + 1.607858 * BoundedHumRatio) / Pressure end if end function GetMoistAirVolume @@ -1072,12 +1176,15 @@ function GetMoistAirDensity(TDryBulb, HumRatio, Pressure) result(MoistAirDensity !+ Atmospheric pressure in Psi [IP] or Pa [SI] real :: MoistAirDensity !+ Moist air density in lb ft⁻³ [IP] or kg m⁻³ [SI] + real :: BoundedHumRatio + !+ Humidity ratio bounded to MIN_HUM_RATIO if (HumRatio < 0.0) then error stop "Error: humidity ratio is negative" end if + BoundedHumRatio = max(HumRatio, MIN_HUM_RATIO) - MoistAirDensity = (1.0 + HumRatio) / GetMoistAirVolume(TDryBulb, HumRatio, Pressure) + MoistAirDensity = (1.0 + BoundedHumRatio) / GetMoistAirVolume(TDryBulb, BoundedHumRatio, Pressure) end function GetMoistAirDensity diff --git a/src/js/psychrolib.js b/src/js/psychrolib.js index 757dc622..d2b35cee 100644 --- a/src/js/psychrolib.js +++ b/src/js/psychrolib.js @@ -57,12 +57,16 @@ function Psychrometrics() { * Global constants *****************************************************************************************************/ - var R_DA_IP = 53.350; // Universal gas constant for dry air (IP version) in ft lb_Force lb_DryAir⁻¹ R⁻¹ - // Reference: ASHRAE Handbook - Fundamentals (2017) ch. 1 - var R_DA_SI = 287.042; // Universal gas constant for dry air (SI version) in J kg_DryAir⁻¹ K⁻¹ - // Reference: ASHRAE Handbook - Fundamentals (2017) ch. 1 + var R_DA_IP = 53.350; // Universal gas constant for dry air (IP version) in ft lb_Force lb_DryAir⁻¹ R⁻¹ + // Reference: ASHRAE Handbook - Fundamentals (2017) ch. 1 + var R_DA_SI = 287.042; // Universal gas constant for dry air (SI version) in J kg_DryAir⁻¹ K⁻¹ + // Reference: ASHRAE Handbook - Fundamentals (2017) ch. 1 + var INVALID = -99999; // Invalid value (dimensionless) - var INVALID = -99999; // Invalid value (dimensionless) + var MAX_ITER_COUNT = 100 // Maximum number of iterations before exiting while loops. + + var MIN_HUM_RATIO = 1e-7 // Minimum acceptable humidity ratio used/returned by any functions. + // Any value above 0 or below the MIN_HUM_RATIO will be reset to this value. /****************************************************************************************************** @@ -252,6 +256,40 @@ function Psychrometrics() { return VapPres / this.GetSatVapPres(TDryBulb); } + // Helper function returning the derivative of the natural log of the saturation vapor pressure + // as a function of dry-bulb temperature. + // Reference: ASHRAE Handbook - Fundamentals (2017) ch. 1 eqn. 5 & 6 + this.dLnPws_ = function // (o) Derivative of natural log of vapor pressure of saturated air in Psi [IP] or Pa [SI] + ( TDryBulb // (i) Dry bulb temperature in °F [IP] or °C [SI] + ) { + var dLnPws, T; + + if (this.isIP()) + { + T = this.GetTRankineFromTFahrenheit(TDryBulb); + + if (TDryBulb <= 32.) + dLnPws = 1.0214165E+04 / pow(T, 2) - 5.3765794E-03 + 2 * 1.9202377E-07 * T + + 2 * 3.5575832E-10 * pow(T, 2) - 4 * 9.0344688E-14 * pow(T, 3) + 4.1635019 / T; + else + dLnPws = 1.0440397E+04 / pow(T, 2) - 2.7022355E-02 + 2 * 1.2890360E-05 * T + - 3 * 2.4780681E-09 * pow(T, 2) + 6.5459673 / T; + } + else + { + T = this.GetTKelvinFromTCelsius(TDryBulb); + + if (TDryBulb <= 0.) + dLnPws = 5.6745359E+03 / pow(T, 2) - 9.677843E-03 + 2 * 6.2215701E-07 * T + + 3 * 2.0747825E-09 * pow(T, 2) - 4 * 9.484024E-13 * pow(T, 3) + 4.1635019 / T; + else + dLnPws = 5.8002206E+03 / pow(T, 2) - 4.8640239E-02 + 2 * 4.1764768E-05 * T + - 3 * 1.4452093E-08 * pow(T, 2) + 6.5459673 / T; + } + + return dLnPws; + } + // Return dew-point temperature given dry-bulb temperature and vapor pressure. // Reference: ASHRAE Handbook - Fundamentals (2017) ch. 1 eqn. 5 and 6 // Notes: the dew point temperature is solved by inverting the equation giving water vapor pressure @@ -266,53 +304,71 @@ function Psychrometrics() { ( TDryBulb // (i) Dry bulb temperature in °F [IP] or °C [SI] , VapPres // (i) Partial pressure of water vapor in moist air in Psi [IP] or Pa [SI] ) { - // Bounds and step size as a function of the system of units - var _STEPSIZE; // Temperature step used for the calculation of numerical derivatives + // Bounds function of the system of units + var BOUNDS // Domain of validity of the equations + var T_WATER_FREEZE, T_WATER_FREEZE_LOW, T_WATER_FREEZE_HIGH, PWS_FREEZE_LOW, PWS_FREEZE_HIGH; + if (this.isIP()) { - var _BOUNDS = [-148., 392.] // Domain of validity of the equations - _STEPSIZE = 0.01 * 9. / 5. + BOUNDS = [-148., 392.]; // Domain of validity of the equations + T_WATER_FREEZE = 32.; } else { - var _BOUNDS = [-100., 200.] // Domain of validity of the equations - _STEPSIZE = 0.01; + BOUNDS = [-100., 200.]; // Domain of validity of the equations + T_WATER_FREEZE = 0.; } - var TMidPoint = (_BOUNDS[0] + _BOUNDS[1]) / 2.; // Midpoint of domain of validity - // Bounds outside which a solution cannot be found - if (VapPres < this.GetSatVapPres(_BOUNDS[0]) || VapPres > this.GetSatVapPres(_BOUNDS[1])) + if (VapPres < this.GetSatVapPres(BOUNDS[0]) || VapPres > this.GetSatVapPres(BOUNDS[1])) throw new Error("Partial pressure of water vapor is outside range of validity of equations"); + // Vapor pressure contained within the discontinuity of the Pws function: return temperature of freezing + T_WATER_FREEZE_LOW = T_WATER_FREEZE - PSYCHROLIB_TOLERANCE / 10.; // Temperature just below freezing + T_WATER_FREEZE_HIGH = T_WATER_FREEZE + PSYCHROLIB_TOLERANCE / 10.; // Temperature just above freezing + PWS_FREEZE_LOW = this.GetSatVapPres(T_WATER_FREEZE_LOW); + PWS_FREEZE_HIGH = this.GetSatVapPres(T_WATER_FREEZE_HIGH); + + // Restrict iteration to either left or right part of the saturation vapor pressure curve + // to avoid iterating back and forth across the discontinuity of the curve at the freezing point + // When the partial pressure of water vapor is within the discontinuity of GetSatVapPres, + // simply return the freezing point of water. + if (VapPres < PWS_FREEZE_LOW) + BOUNDS[1] = T_WATER_FREEZE_LOW; + else if (VapPres > PWS_FREEZE_HIGH) + BOUNDS[0] = T_WATER_FREEZE_HIGH; + else + return T_WATER_FREEZE; + + // We use NR to approximate the solution. // First guess - var Tdp = TDryBulb; // Calculated value of dew point temperatures, solved for iteratively in °F [IP] or °C [SI] - var lnVP = log(VapPres); // Natural logarithm of partial pressure of water vapor pressure in moist air + var TDewPoint = TDryBulb; // Calculated value of dew point temperatures, solved for iteratively in °F [IP] or °C [SI] + var lnVP = log(VapPres); // Natural logarithm of partial pressure of water vapor pressure in moist air - var Tdp_c; // Value of Tdp used in NR calculation - var lnVP_c; // Value of log of vapor water pressure used in NR calculation - var d_Tdp; // Value of temperature step used in NR calculation + var TDewPoint_iter; // Value of TDewPoint used in NR calculation + var lnVP_iter; // Value of log of vapor water pressure used in NR calculation + var index = 1; do { // Current point - Tdp_c = Tdp; - lnVP_c = log(this.GetSatVapPres(Tdp_c)); - // Step - negative in the right part of the curve, positive in the left part - // to avoid going past the domain of validity of eqn. 5 and 6 - // when Tdp_c is close to its bounds - if (Tdp_c > TMidPoint) - d_Tdp = -_STEPSIZE; - else - d_Tdp = _STEPSIZE; - // Derivative of function, calculated numerically - var d_lnVP = (log(this.GetSatVapPres(Tdp_c + d_Tdp)) - lnVP_c) / d_Tdp; + TDewPoint_iter = TDewPoint; + lnVP_iter = log(this.GetSatVapPres(TDewPoint_iter)); + + // Derivative of function, calculated analytically + var d_lnVP = this.dLnPws_(TDewPoint_iter); + // New estimate, bounded by domain of validity of eqn. 5 and 6 - Tdp = Tdp_c - (lnVP_c - lnVP) / d_lnVP; - Tdp = max(Tdp, _BOUNDS[0]); - Tdp = min(Tdp, _BOUNDS[1]); + TDewPoint = TDewPoint_iter - (lnVP_iter - lnVP) / d_lnVP; + TDewPoint = max(TDewPoint, BOUNDS[0]); + TDewPoint = min(TDewPoint, BOUNDS[1]); + + if (index > MAX_ITER_COUNT) + throw new Error("Convergence not reached in GetTDewPointFromVapPres. Stopping."); + + index++; } - while (abs(Tdp - Tdp_c) > PSYCHROLIB_TOLERANCE); - return min(Tdp, TDryBulb); + while (abs(TDewPoint - TDewPoint_iter) > PSYCHROLIB_TOLERANCE); + return min(TDewPoint, TDryBulb); } // Return vapor pressure given dew point temperature. @@ -337,12 +393,14 @@ function Psychrometrics() { ) { // Declarations var Wstar; - var TDewPoint, TWetBulb, TWetBulbSup, TWetBulbInf; + var TDewPoint, TWetBulb, TWetBulbSup, TWetBulbInf, BoundedHumRatio; + var index = 1; if (!(HumRatio >= 0.)) throw new Error("Humidity ratio is negative"); + BoundedHumRatio = max(HumRatio, MIN_HUM_RATIO); - TDewPoint = this.GetTDewPointFromHumRatio(TDryBulb, HumRatio, Pressure); + TDewPoint = this.GetTDewPointFromHumRatio(TDryBulb, BoundedHumRatio, Pressure); // Initial guesses TWetBulbSup = TDryBulb; @@ -350,18 +408,23 @@ function Psychrometrics() { TWetBulb = (TWetBulbInf + TWetBulbSup) / 2.; // Bisection loop - while (TWetBulbSup - TWetBulbInf > PSYCHROLIB_TOLERANCE) { + while ((TWetBulbSup - TWetBulbInf) > PSYCHROLIB_TOLERANCE) { // Compute humidity ratio at temperature Tstar Wstar = this.GetHumRatioFromTWetBulb(TDryBulb, TWetBulb, Pressure); // Get new bounds - if (Wstar > HumRatio) + if (Wstar > BoundedHumRatio) TWetBulbSup = TWetBulb; else TWetBulbInf = TWetBulb; // New guess of wet bulb temperature TWetBulb = (TWetBulbSup + TWetBulbInf) / 2.; + + if (index > MAX_ITER_COUNT) + throw new Error("Convergence not reached in GetTWetBulbFromHumRatio. Stopping."); + + index++; } return TWetBulb; @@ -400,8 +463,8 @@ function Psychrometrics() { HumRatio = ((2830. - 0.24 * TWetBulb) * Wsstar - 1.006 * (TDryBulb - TWetBulb)) / (2830. + 1.86 * TDryBulb - 2.1 * TWetBulb); } - - return HumRatio; + // Validity check. + return max(HumRatio, MIN_HUM_RATIO); } // Return humidity ratio given dry-bulb temperature, relative humidity, and pressure. @@ -475,10 +538,15 @@ function Psychrometrics() { ( VapPres // (i) Partial pressure of water vapor in moist air in Psi [IP] or Pa [SI] , Pressure // (i) Atmospheric pressure in Psi [IP] or Pa [SI] ) { + var HumRatio; + if (!(VapPres >= 0.)) throw new Error("Partial pressure of water vapor in moist air is negative"); - return 0.621945 * VapPres / (Pressure - VapPres); + HumRatio = 0.621945 * VapPres / (Pressure - VapPres); + + // Validity check. + return max(HumRatio, MIN_HUM_RATIO); } // Return vapor pressure given humidity ratio and pressure. @@ -487,12 +555,13 @@ function Psychrometrics() { ( HumRatio // (i) Humidity ratio in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI] , Pressure // (i) Atmospheric pressure in Psi [IP] or Pa [SI] ) { - var VapPres; + var VapPres, BoundedHumRatio; if (!(HumRatio >= 0.)) throw new Error("Humidity ratio is negative"); + BoundedHumRatio = max(HumRatio, MIN_HUM_RATIO); - VapPres = Pressure * HumRatio / (0.621945 + HumRatio); + VapPres = Pressure * BoundedHumRatio / (0.621945 + BoundedHumRatio); return VapPres; } @@ -506,10 +575,12 @@ function Psychrometrics() { this.GetSpecificHumFromHumRatio = function // (o) Specific humidity ratio in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI] ( HumRatio // (i) Humidity ratio in lb_H₂O lb_Dry_Air⁻¹ [IP] or kg_H₂O kg_Dry_Air⁻¹ [SI] ) { + var BoundedHumRatio; if (!(HumRatio >= 0.)) throw new Error("Humidity ratio is negative"); + BoundedHumRatio = max(HumRatio, MIN_HUM_RATIO); - return HumRatio / (1.0 + HumRatio); + return BoundedHumRatio / (1.0 + BoundedHumRatio); } // Return the humidity ratio (aka mixing ratio) from specific humidity @@ -517,10 +588,15 @@ function Psychrometrics() { this.GetHumRatioFromSpecificHum = function // (o) Humidity ratio in lb_H₂O lb_Dry_Air⁻¹ [IP] or kg_H₂O kg_Dry_Air⁻¹ [SI] ( SpecificHum // (i) Specific humidity ratio in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI] ) { + var HumRatio; + if (!(SpecificHum >= 0.0 && SpecificHum < 1.0)) throw new Error("Specific humidity is outside range [0, 1["); - return SpecificHum / (1.0 - SpecificHum); + HumRatio = SpecificHum / (1.0 - SpecificHum); + + // Validity check + return max(HumRatio, MIN_HUM_RATIO); } @@ -576,13 +652,15 @@ function Psychrometrics() { ( MoistAirEnthalpy // (i) Moist air enthalpy in Btu lb⁻¹ [IP] or J kg⁻¹ , HumRatio // (i) Humidity ratio in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI] ) { + var BoundedHumRatio; if (!(HumRatio >= 0.)) throw new Error("Humidity ratio is negative"); + BoundedHumRatio = max(HumRatio, MIN_HUM_RATIO); if (this.isIP()) - return (MoistAirEnthalpy - 1061.0 * HumRatio) / (0.240 + 0.444 * HumRatio); + return (MoistAirEnthalpy - 1061.0 * BoundedHumRatio) / (0.240 + 0.444 * BoundedHumRatio); else - return (MoistAirEnthalpy / 1000.0 - 2501.0 * HumRatio) / (1.006 + 1.86 * HumRatio); + return (MoistAirEnthalpy / 1000.0 - 2501.0 * BoundedHumRatio) / (1.006 + 1.86 * BoundedHumRatio); } // Return humidity ratio from enthalpy and dry-bulb temperature. @@ -592,10 +670,14 @@ function Psychrometrics() { ( MoistAirEnthalpy // (i) Moist air enthalpy in Btu lb⁻¹ [IP] or J kg⁻¹ , TDryBulb // (i) Dry-bulb temperature in °F [IP] or °C [SI] ) { + var HumRatio; if (this.isIP()) - return (MoistAirEnthalpy - 0.240 * TDryBulb) / (1061.0 + 0.444 * TDryBulb); + HumRatio = (MoistAirEnthalpy - 0.240 * TDryBulb) / (1061.0 + 0.444 * TDryBulb); else - return (MoistAirEnthalpy / 1000.0 - 1.006 * TDryBulb) / (2501.0 + 1.86 * TDryBulb); + HumRatio = (MoistAirEnthalpy / 1000.0 - 1.006 * TDryBulb) / (2501.0 + 1.86 * TDryBulb); + + // Validity check. + return max(HumRatio, MIN_HUM_RATIO); } @@ -650,10 +732,13 @@ function Psychrometrics() { ( TDryBulb // (i) Dry bulb temperature in °F [IP] or °C [SI] , Pressure // (i) Atmospheric pressure in Psi [IP] or Pa [SI] ) { - var SatVaporPres; + var SatVaporPres, SatHumRatio; SatVaporPres = this.GetSatVapPres(TDryBulb); - return 0.621945 * SatVaporPres / (Pressure - SatVaporPres); + SatHumRatio = 0.621945 * SatVaporPres / (Pressure - SatVaporPres); + + // Validity check. + return max(SatHumRatio, MIN_HUM_RATIO); } // Return saturated air enthalpy given dry-bulb temperature and pressure. @@ -695,11 +780,13 @@ function Psychrometrics() { , HumRatio // (i) Humidity ratio in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI] , Pressure // (i) Atmospheric pressure in Psi [IP] or Pa [SI] ) { + var BoundedHumRatio; if (!(HumRatio >= 0.)) throw new Error("Humidity ratio is negative"); + BoundedHumRatio = max(HumRatio, MIN_HUM_RATIO); - return HumRatio / this.GetSatHumRatio(TDryBulb, Pressure); + return BoundedHumRatio / this.GetSatHumRatio(TDryBulb, Pressure); } // Return moist air enthalpy given dry-bulb temperature and humidity ratio. @@ -708,11 +795,16 @@ function Psychrometrics() { ( TDryBulb // (i) Dry bulb temperature in °F [IP] or °C [SI] , HumRatio // (i) Humidity ratio in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI] ) { + var BoundedHumRatio; + + if (!(HumRatio >= 0.)) + throw new Error("Humidity ratio is negative"); + BoundedHumRatio = max(HumRatio, MIN_HUM_RATIO); if (this.isIP()) - return 0.240 * TDryBulb + HumRatio * (1061. + 0.444 * TDryBulb); + return 0.240 * TDryBulb + BoundedHumRatio * (1061. + 0.444 * TDryBulb); else - return (1.006 * TDryBulb + HumRatio * (2501. + 1.86 * TDryBulb)) * 1000.; + return (1.006 * TDryBulb + BoundedHumRatio * (2501. + 1.86 * TDryBulb)) * 1000.; } // Return moist air specific volume given dry-bulb temperature, humidity ratio, and pressure. @@ -724,14 +816,16 @@ function Psychrometrics() { , HumRatio // (i) Humidity ratio in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI] , Pressure // (i) Atmospheric pressure in Psi [IP] or Pa [SI] ) { + var BoundedHumRatio; if (!(HumRatio >= 0.)) throw new Error("Humidity ratio is negative"); + BoundedHumRatio = max(HumRatio, MIN_HUM_RATIO); if (this.isIP()) - return R_DA_IP * this.GetTRankineFromTFahrenheit(TDryBulb) * (1. + 1.607858 * HumRatio) / (144. * Pressure); + return R_DA_IP * this.GetTRankineFromTFahrenheit(TDryBulb) * (1. + 1.607858 * BoundedHumRatio) / (144. * Pressure); else - return R_DA_SI * this.GetTKelvinFromTCelsius(TDryBulb) * (1. + 1.607858 * HumRatio) / Pressure; + return R_DA_SI * this.GetTKelvinFromTCelsius(TDryBulb) * (1. + 1.607858 * BoundedHumRatio) / Pressure; } // Return moist air density given humidity ratio, dry bulb temperature, and pressure. @@ -741,11 +835,13 @@ function Psychrometrics() { , HumRatio // (i) Humidity ratio in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI] , Pressure // (i) Atmospheric pressure in Psi [IP] or Pa [SI] ) { + var BoundedHumRatio; if (!(HumRatio >= 0.)) throw new Error("Humidity ratio is negative"); + BoundedHumRatio = max(HumRatio, MIN_HUM_RATIO); - return (1. + HumRatio) / this.GetMoistAirVolume(TDryBulb, HumRatio, Pressure); + return (1. + BoundedHumRatio) / this.GetMoistAirVolume(TDryBulb, BoundedHumRatio, Pressure); } diff --git a/src/python/psychrolib.py b/src/python/psychrolib.py index 0a9ff266..cd3a0181 100644 --- a/src/python/psychrolib.py +++ b/src/python/psychrolib.py @@ -73,6 +73,17 @@ """ +MAX_ITER_COUNT = 100 +"""int: Maximum number of iterations before exiting while loops. + +""" + +MIN_HUM_RATIO = 1e-7 +"""float: Minimum acceptable humidity ratio used/returned by any functions. + Any value above 0 or below the MIN_HUM_RATIO will be reset to this value. + +""" + ####################################################################################################### # Helper functions @@ -116,7 +127,7 @@ def SetUnitSystem(Units: UnitSystem) -> None: # Define tolerance on temperature calculations # The tolerance is the same in IP and SI if Units == IP: - PSYCHROLIB_TOLERANCE = 0.001 * 9 / 5 + PSYCHROLIB_TOLERANCE = 0.001 * 9. / 5. else: PSYCHROLIB_TOLERANCE = 0.001 @@ -375,6 +386,40 @@ def GetRelHumFromVapPres(TDryBulb: float, VapPres: float) -> float: RelHum = VapPres / GetSatVapPres(TDryBulb) return RelHum +def dLnPws_(TDryBulb: float) -> float: + """ + Helper function returning the derivative of the natural log of the saturation vapor pressure + as a function of dry-bulb temperature. + + Args: + TDryBulb : Dry-bulb temperature in °F [IP] or °C [SI] + + Returns: + Derivative of natural log of vapor pressure of saturated air in Psi [IP] or Pa [SI] + + Reference: + ASHRAE Handbook - Fundamentals (2017) ch. 1 eqn 5 & 6 + + """ + if isIP(): + T = GetTRankineFromTFahrenheit(TDryBulb) + if TDryBulb <= 32.: + dLnPws = 1.0214165E+04 / math.pow(T, 2) - 5.3765794E-03 + 2 * 1.9202377E-07 * T \ + + 2 * 3.5575832E-10 * math.pow(T, 2) - 4 * 9.0344688E-14 * math.pow(T, 3) + 4.1635019 / T + else: + dLnPws = 1.0440397E+04 / math.pow(T, 2) - 2.7022355E-02 + 2 * 1.2890360E-05 * T \ + - 3 * 2.4780681E-09 * math.pow(T, 2) + 6.5459673 / T + else: + T = GetTKelvinFromTCelsius(TDryBulb) + if TDryBulb <= 0.: + dLnPws = 5.6745359E+03 / math.pow(T, 2) - 9.677843E-03 + 2 * 6.2215701E-07 * T \ + + 3 * 2.0747825E-09 * math.pow(T, 2) - 4 * 9.484024E-13 * math.pow(T, 3) + 4.1635019 / T + else: + dLnPws = 5.8002206E+03 / math.pow(T, 2) - 4.8640239E-02 + 2 * 4.1764768E-05 * T \ + - 3 * 1.4452093E-08 * math.pow(T, 2) + 6.5459673 / T + + return dLnPws + def GetTDewPointFromVapPres(TDryBulb: float, VapPres: float) -> float: """ Return dew-point temperature given dry-bulb temperature and vapor pressure. @@ -401,42 +446,59 @@ def GetTDewPointFromVapPres(TDryBulb: float, VapPres: float) -> float: """ if isIP(): - _BOUNDS = -148, 392 - _STEPSIZE = 0.01 * 9 / 5 + BOUNDS = [-148, 392] + T_WATER_FREEZE = 32. else: - _BOUNDS = -100, 200 - _STEPSIZE = 0.01 - - TMidPoint = (_BOUNDS[0] + _BOUNDS[1]) / 2. # Midpoint of domain of validity + BOUNDS = [-100, 200] + T_WATER_FREEZE = 0. - if VapPres < GetSatVapPres(_BOUNDS[0]) or VapPres > GetSatVapPres(_BOUNDS[1]): + # Validity check -- bounds outside which a solution cannot be found + if VapPres < GetSatVapPres(BOUNDS[0]) or VapPres > GetSatVapPres(BOUNDS[1]): raise ValueError("Partial pressure of water vapor is outside range of validity of equations") + # Vapor pressure contained within the discontinuity of the Pws function: return temperature of freezing + T_WATER_FREEZE_LOW = T_WATER_FREEZE - PSYCHROLIB_TOLERANCE / 10. # Temperature just below freezing + T_WATER_FREEZE_HIGH = T_WATER_FREEZE + PSYCHROLIB_TOLERANCE / 10. # Temperature just above freezing + PWS_FREEZE_LOW = GetSatVapPres(T_WATER_FREEZE_LOW) + PWS_FREEZE_HIGH = GetSatVapPres(T_WATER_FREEZE_HIGH) + + # Restrict iteration to either left or right part of the saturation vapor pressure curve + # to avoid iterating back and forth across the discontinuity of the curve at the freezing point + # When the partial pressure of water vapor is within the discontinuity of GetSatVapPres, + # simply return the freezing point of water. + if (VapPres < PWS_FREEZE_LOW): + BOUNDS[1] = T_WATER_FREEZE_LOW + elif (VapPres > PWS_FREEZE_HIGH): + BOUNDS[0] = T_WATER_FREEZE_HIGH + else: + return T_WATER_FREEZE + + # We use NR to approximate the solution. # First guess TDewPoint = TDryBulb # Calculated value of dew point temperatures, solved for iteratively - lnVP = math.log(VapPres) # Partial pressure of water vapor in moist air - while True: - TDewPoint_iter = TDewPoint # Value of Tdp used in NR calculation - # Step - negative in the right part of the curve, positive in the left part - # to avoid going past the domain of validity of eqn. 5 and 6 - # when TDewPoint_iter is close to its bounds - if TDewPoint_iter > TMidPoint: - StepSize = -_STEPSIZE - else: - StepSize = _STEPSIZE + index = 1 + while True: + TDewPoint_iter = TDewPoint # TDewPoint used in NR calculation lnVP_iter = math.log(GetSatVapPres(TDewPoint_iter)) - # Derivative of function, calculated numerically - d_lnVP = (math.log(GetSatVapPres(TDewPoint_iter + StepSize)) - lnVP_iter) / StepSize - # New estimate, bounded by domain of validity of eqn. 5 and 6 + + # Derivative of function, calculated analytically + d_lnVP = dLnPws_(TDewPoint_iter) + + # New estimate, bounded by the search domain defined above TDewPoint = TDewPoint_iter - (lnVP_iter - lnVP) / d_lnVP - TDewPoint = max(TDewPoint, _BOUNDS[0]) - TDewPoint = min(TDewPoint, _BOUNDS[1]) + TDewPoint = max(TDewPoint, BOUNDS[0]) + TDewPoint = min(TDewPoint, BOUNDS[1]) + + if ((math.fabs(TDewPoint - TDewPoint_iter) <= PSYCHROLIB_TOLERANCE)): + break + + if (index > MAX_ITER_COUNT): + raise ValueError("Convergence not reached in GetTDewPointFromVapPres. Stopping.") - if math.fabs(TDewPoint - TDewPoint_iter) <= PSYCHROLIB_TOLERANCE: - break + index = index + 1 TDewPoint = min(TDewPoint, TDryBulb) return TDewPoint @@ -481,28 +543,35 @@ def GetTWetBulbFromHumRatio(TDryBulb: float, HumRatio: float, Pressure: float) - """ if HumRatio < 0: raise ValueError("Humidity ratio cannot be negative") + BoundedHumRatio = max(HumRatio, MIN_HUM_RATIO) - TDewPoint = GetTDewPointFromHumRatio(TDryBulb, HumRatio, Pressure) + TDewPoint = GetTDewPointFromHumRatio(TDryBulb, BoundedHumRatio, Pressure) # Initial guesses TWetBulbSup = TDryBulb TWetBulbInf = TDewPoint TWetBulb = (TWetBulbInf + TWetBulbSup) / 2 + index = 1 # Bisection loop - while (TWetBulbSup - TWetBulbInf > PSYCHROLIB_TOLERANCE): + while ((TWetBulbSup - TWetBulbInf) > PSYCHROLIB_TOLERANCE): # Compute humidity ratio at temperature Tstar Wstar = GetHumRatioFromTWetBulb(TDryBulb, TWetBulb, Pressure) # Get new bounds - if Wstar > HumRatio: + if Wstar > BoundedHumRatio: TWetBulbSup = TWetBulb else: TWetBulbInf = TWetBulb # New guess of wet bulb temperature TWetBulb = (TWetBulbSup + TWetBulbInf) / 2 + + if (index >= MAX_ITER_COUNT): + raise ValueError("Convergence not reached in GetTWetBulbFromHumRatio. Stopping.") + + index = index + 1 return TWetBulb def GetHumRatioFromTWetBulb(TDryBulb: float, TWetBulb: float, Pressure: float) -> float: @@ -540,7 +609,8 @@ def GetHumRatioFromTWetBulb(TDryBulb: float, TWetBulb: float, Pressure: float) - else: HumRatio = ((2830. - 0.24 * TWetBulb) * Wsstar - 1.006 * (TDryBulb - TWetBulb)) \ / (2830. + 1.86 * TDryBulb - 2.1 * TWetBulb) - return HumRatio + # Validity check. + return max(HumRatio, MIN_HUM_RATIO) def GetHumRatioFromRelHum(TDryBulb: float, RelHum: float, Pressure: float) -> float: """ @@ -654,7 +724,9 @@ def GetHumRatioFromVapPres(VapPres: float, Pressure: float) -> float: raise ValueError("Partial pressure of water vapor in moist air cannot be negative") HumRatio = 0.621945 * VapPres / (Pressure - VapPres) - return HumRatio + + # Validity check. + return max(HumRatio, MIN_HUM_RATIO) def GetVapPresFromHumRatio(HumRatio: float, Pressure: float) -> float: """ @@ -673,8 +745,9 @@ def GetVapPresFromHumRatio(HumRatio: float, Pressure: float) -> float: """ if HumRatio < 0: raise ValueError("Humidity ratio is negative") + BoundedHumRatio = max(HumRatio, MIN_HUM_RATIO) - VapPres = Pressure * HumRatio / (0.621945 + HumRatio) + VapPres = Pressure * BoundedHumRatio / (0.621945 + BoundedHumRatio) return VapPres @@ -698,8 +771,9 @@ def GetSpecificHumFromHumRatio(HumRatio: float) -> float: """ if HumRatio < 0: raise ValueError("Humidity ratio cannot be negative") + BoundedHumRatio = max(HumRatio, MIN_HUM_RATIO) - SpecificHum = HumRatio / (1.0 + HumRatio) + SpecificHum = BoundedHumRatio / (1.0 + BoundedHumRatio) return SpecificHum def GetHumRatioFromSpecificHum(SpecificHum: float) -> float: @@ -720,7 +794,9 @@ def GetHumRatioFromSpecificHum(SpecificHum: float) -> float: raise ValueError("Specific humidity is outside range [0, 1[") HumRatio = SpecificHum / (1.0 - SpecificHum) - return SpecificHum + + # Validity check. + return max(HumRatio, MIN_HUM_RATIO) ####################################################################################################### @@ -821,11 +897,12 @@ def GetTDryBulbFromEnthalpyAndHumRatio(MoistAirEnthalpy: float, HumRatio: float) """ if HumRatio < 0: raise ValueError("Humidity ratio is negative") + BoundedHumRatio = max(HumRatio, MIN_HUM_RATIO) if isIP(): - TDryBulb = (MoistAirEnthalpy - 1061.0 * HumRatio) / (0.240 + 0.444 * HumRatio) + TDryBulb = (MoistAirEnthalpy - 1061.0 * BoundedHumRatio) / (0.240 + 0.444 * BoundedHumRatio) else: - TDryBulb = (MoistAirEnthalpy / 1000.0 - 2501.0 * HumRatio) / (1.006 + 1.86 * HumRatio) + TDryBulb = (MoistAirEnthalpy / 1000.0 - 2501.0 * BoundedHumRatio) / (1.006 + 1.86 * BoundedHumRatio) return TDryBulb def GetHumRatioFromEnthalpyAndTDryBulb(MoistAirEnthalpy: float, TDryBulb: float) -> float: @@ -851,7 +928,9 @@ def GetHumRatioFromEnthalpyAndTDryBulb(MoistAirEnthalpy: float, TDryBulb: float) HumRatio = (MoistAirEnthalpy - 0.240 * TDryBulb) / (1061.0 + 0.444 * TDryBulb) else: HumRatio = (MoistAirEnthalpy / 1000.0 - 1.006 * TDryBulb) / (2501.0 + 1.86 * TDryBulb) - return HumRatio + + # Validity check. + return max(HumRatio, MIN_HUM_RATIO) ####################################################################################################### @@ -917,7 +996,9 @@ def GetSatHumRatio(TDryBulb: float, Pressure: float) -> float: """ SatVaporPres = GetSatVapPres(TDryBulb) SatHumRatio = 0.621945 * SatVaporPres / (Pressure - SatVaporPres) - return SatHumRatio + + # Validity check. + return max(SatHumRatio, MIN_HUM_RATIO) def GetSatAirEnthalpy(TDryBulb: float, Pressure: float) -> float: """ @@ -988,9 +1069,10 @@ def GetDegreeOfSaturation(TDryBulb: float, HumRatio: float, Pressure: float) -> """ if HumRatio < 0: raise ValueError("Humidity ratio is negative") + BoundedHumRatio = max(HumRatio, MIN_HUM_RATIO) SatHumRatio = GetSatHumRatio(TDryBulb, Pressure) - DegreeOfSaturation = HumRatio / SatHumRatio + DegreeOfSaturation = BoundedHumRatio / SatHumRatio return DegreeOfSaturation def GetMoistAirEnthalpy(TDryBulb: float, HumRatio: float) -> float: @@ -1010,11 +1092,12 @@ def GetMoistAirEnthalpy(TDryBulb: float, HumRatio: float) -> float: """ if HumRatio < 0: raise ValueError("Humidity ratio is negative") + BoundedHumRatio = max(HumRatio, MIN_HUM_RATIO) if isIP(): - MoistAirEnthalpy = 0.240 * TDryBulb + HumRatio * (1061 + 0.444 * TDryBulb) + MoistAirEnthalpy = 0.240 * TDryBulb + BoundedHumRatio * (1061 + 0.444 * TDryBulb) else: - MoistAirEnthalpy = (1.006 * TDryBulb + HumRatio * (2501. + 1.86 * TDryBulb)) * 1000 + MoistAirEnthalpy = (1.006 * TDryBulb + BoundedHumRatio * (2501. + 1.86 * TDryBulb)) * 1000 return MoistAirEnthalpy def GetMoistAirVolume(TDryBulb: float, HumRatio: float, Pressure: float) -> float: @@ -1039,11 +1122,12 @@ def GetMoistAirVolume(TDryBulb: float, HumRatio: float, Pressure: float) -> floa """ if HumRatio < 0: raise ValueError("Humidity ratio is negative") + BoundedHumRatio = max(HumRatio, MIN_HUM_RATIO) if isIP(): - MoistAirVolume = R_DA_IP * GetTRankineFromTFahrenheit(TDryBulb) * (1 + 1.607858 * HumRatio) / (144 * Pressure) + MoistAirVolume = R_DA_IP * GetTRankineFromTFahrenheit(TDryBulb) * (1 + 1.607858 * BoundedHumRatio) / (144 * Pressure) else: - MoistAirVolume = R_DA_SI * GetTKelvinFromTCelsius(TDryBulb) * (1 + 1.607858 * HumRatio) / Pressure + MoistAirVolume = R_DA_SI * GetTKelvinFromTCelsius(TDryBulb) * (1 + 1.607858 * BoundedHumRatio) / Pressure return MoistAirVolume def GetMoistAirDensity(TDryBulb: float, HumRatio: float, Pressure:float) -> float: @@ -1064,9 +1148,10 @@ def GetMoistAirDensity(TDryBulb: float, HumRatio: float, Pressure:float) -> floa """ if HumRatio < 0: raise ValueError("Humidity ratio is negative") + BoundedHumRatio = max(HumRatio, MIN_HUM_RATIO) - MoistAirVolume = GetMoistAirVolume(TDryBulb, HumRatio, Pressure) - MoistAirDensity = (1 + HumRatio) / MoistAirVolume + MoistAirVolume = GetMoistAirVolume(TDryBulb, BoundedHumRatio, Pressure) + MoistAirDensity = (1 + BoundedHumRatio) / MoistAirVolume return MoistAirDensity diff --git a/src/vba/psychrolib.bas b/src/vba/psychrolib.bas index c75ffd62..15012ab1 100644 --- a/src/vba/psychrolib.bas +++ b/src/vba/psychrolib.bas @@ -64,6 +64,9 @@ End Enum Private Const R_DA_IP = 53.35 ' Universal gas constant for dry air (IP version) in ft lbf/lb_DryAir/R Private Const R_DA_SI = 287.042 ' Universal gas constant for dry air (SI version) in J/kg_DryAir/K +Private Const MAX_ITER_COUNT = 100 ' Maximum number of iterations before exiting while loops. +Private Const MIN_HUM_RATIO = 1e-7 ' Minimum acceptable humidity ratio used/returned by any functions. + ' Any value above 0 or below the MIN_HUM_RATIO will be reset to this value. '****************************************************************************************************** @@ -481,6 +484,43 @@ ErrHandler: End Function + +Private Function dLnPws_(TDryBulb As Variant) As Variant +' +' Helper function returning the derivative of the natural log of the saturation vapor pressure +' as a function of dry-bulb temperature. +' +' Args: +' TDryBulb : Dry-bulb temperature in °F [IP] or °C [SI] +' +' Returns: +' Derivative of natural log of vapor pressure of saturated air in Psi [IP] or Pa [SI] +' +' Reference: +' ASHRAE Handbook - Fundamentals (2017) ch. 1 eqn 5 & 6 +' + Dim T As Variant + If (isIP()) Then + T = GetTRankineFromTFahrenheit(TDryBulb) + If (TDryBulb <= 32.) Then + dLnPws_ = 10214.165 / T ^ 2 - 0.0053765794 + 2 * 0.00000019202377 * T _ + + 2 * 3.5575832E-10 * T ^ 2 - 4 * 9.0344688E-14 * T ^ 3 + 4.1635019 / T + Else + dLnPws_ = 10440.397 / T ^ 2 - 0.027022355 + 2 * 0.00001289036 * T _ + - 3 * 2.4780681E-09 * T ^ 2 + 6.5459673 / T + End If + Else + T = GetTKelvinFromTCelsius(TDryBulb) + If (TDryBulb <= 0.) Then + dLnPws_ = 5674.5359 / T ^ 2 - 0.009677843 + 2 * 0.00000062215701 * T _ + + 3 * 2.0747825E-09 * T ^ 2 - 4 * 9.484024E-13 * T ^ 3 + 4.1635019 / T + Else + dLnPws_ = 5800.2206 / T ^ 2 - 0.048640239 + 2 * 0.000041764768 * T _ + - 3 * 0.000000014452093 * T ^ 2 + 6.5459673 / T + End If + End If +End Function + Function GetTDewPointFromVapPres(ByVal TDryBulb As Variant, ByVal VapPres As Variant) As Variant ' ' Return dew-point temperature given dry-bulb temperature and vapor pressure. @@ -505,62 +545,81 @@ Function GetTDewPointFromVapPres(ByVal TDryBulb As Variant, ByVal VapPres As Var ' Convergence is usually achieved in 3 to 5 iterations. ' TDryBulb is not really needed here, just used for convenience. ' - Dim BOUNDS_(2) As Variant - Dim STEPSIZE_ As Variant + Dim BOUNDS(2) As Variant + Dim T_WATER_FREEZE As Variant, T_WATER_FREEZE_LOW As Variant, T_WATER_FREEZE_HIGH As Variant + Dim PWS_FREEZE_LOW As Variant, PWS_FREEZE_HIGH As Variant + Dim PSYCHROLIB_TOLERANCE As Variant + If (isIP()) Then - BOUNDS_(1) = -148 - BOUNDS_(2) = 392 - STEPSIZE_ = 0.01 * 9 / 5 + BOUNDS(1) = -148. + BOUNDS(2) = 392. + T_WATER_FREEZE = 32. Else - BOUNDS_(1) = -100 - BOUNDS_(2) = 200 - STEPSIZE_ = 0.01 + BOUNDS(1) = -100. + BOUNDS(2) = 200. + T_WATER_FREEZE = 0. End If - Dim TMidPoint As Variant - TMidPoint = (BOUNDS_(1) + BOUNDS_(2)) / 2# ' Midpoint of domain of validity - On Error GoTo ErrHandler - If ((VapPres < GetSatVapPres(BOUNDS_(1))) Or (VapPres > GetSatVapPres(BOUNDS_(2)))) Then + If ((VapPres < GetSatVapPres(BOUNDS(1))) Or (VapPres > GetSatVapPres(BOUNDS(2)))) Then MyMsgBox ("Partial pressure of water vapor is outside range of validity of equations") GoTo ErrHandler End If - ' First guess + PSYCHROLIB_TOLERANCE = GetTol() + ' Vapor pressure contained within the discontinuity of the Pws function: return temperature of freezing + T_WATER_FREEZE_LOW = T_WATER_FREEZE - PSYCHROLIB_TOLERANCE / 10. ' Temperature just below freezing + T_WATER_FREEZE_HIGH = T_WATER_FREEZE + PSYCHROLIB_TOLERANCE / 10. ' Temperature just above freezing + PWS_FREEZE_LOW = GetSatVapPres(T_WATER_FREEZE_LOW) + PWS_FREEZE_HIGH = GetSatVapPres(T_WATER_FREEZE_HIGH) + + ' Restrict iteration to either left or right part of the saturation vapor pressure curve + ' to avoid iterating back and forth across the discontinuity of the curve at the freezing point + ' When the partial pressure of water vapor is within the discontinuity of GetSatVapPres, + ' simply return the freezing point of water. + If (VapPres < PWS_FREEZE_LOW) Then + BOUNDS(2) = T_WATER_FREEZE_LOW + ElseIf (VapPres > PWS_FREEZE_HIGH) Then + BOUNDS(1) = T_WATER_FREEZE_HIGH + Else + GetTDewPointFromVapPres = T_WATER_FREEZE + Exit Function + End If + Dim TDewPoint As Variant Dim lnVP As Variant Dim d_lnVP As Variant Dim TDewPoint_iter As Variant - Dim StepSize As Variant Dim lnVP_iter + Dim index As Variant + index = 1 + + ' We use NR to approximate the solution. + ' First guess TDewPoint = TDryBulb ' Calculated value of dew point temperatures, solved for iteratively - Dim Tol As Variant + lnVP = Log(VapPres) ' Partial pressure of water vapor in moist air - lnVP = Log(VapPres) ' Partial pressure of water vapor in moist air - Tol = GetTol() Do TDewPoint_iter = TDewPoint ' Value of Tdp used in NR calculation + lnVP_iter = Log(GetSatVapPres(TDewPoint_iter)) - ' Step - negative in the right part of the curve, positive in the left part - ' to avoid going past the domain of validity of eqn. 5 and 6 - ' when TDewPoint_iter is close to its bounds - If (TDewPoint_iter > TMidPoint) Then - StepSize = -STEPSIZE_ - Else - StepSize = STEPSIZE_ - End If + ' Derivative of function, calculated analytically + d_lnVP = dLnPws_(TDewPoint_iter) - lnVP_iter = Log(GetSatVapPres(TDewPoint_iter)) - ' Derivative of function, calculated numerically - d_lnVP = (Log(GetSatVapPres(TDewPoint_iter + StepSize)) - lnVP_iter) / StepSize - ' New estimate, bounded by domain of validity of eqn. 5 and 6 + ' New estimate, bounded by domain of validity of eqn. 5 and 6 and by the freezing point TDewPoint = TDewPoint_iter - (lnVP_iter - lnVP) / d_lnVP - TDewPoint = Max(TDewPoint, BOUNDS_(1)) - TDewPoint = Min(TDewPoint, BOUNDS_(2)) + TDewPoint = Max(TDewPoint, BOUNDS(1)) + TDewPoint = Min(TDewPoint, BOUNDS(2)) + + If (index > MAX_ITER_COUNT) Then + GoTo ErrHandler + End If + + index = index + 1 - Loop While (Abs(TDewPoint - TDewPoint_iter) > Tol) + Loop While (Abs(TDewPoint - TDewPoint_iter) > PSYCHROLIB_TOLERANCE) TDewPoint = Min(TDewPoint, TDryBulb) GetTDewPointFromVapPres = TDewPoint @@ -617,7 +676,7 @@ Function GetTWetBulbFromHumRatio(ByVal TDryBulb As Variant, ByVal HumRatio As Va ' Declarations Dim Wstar As Variant Dim TDewPoint As Variant, TWetBulb As Variant, TWetBulbSup As Variant, TWetBulbInf As Variant - Dim Tol As Variant + Dim Tol As Variant, BoundedHumRatio As Variant, index As Variant On Error GoTo ErrHandler @@ -625,8 +684,9 @@ Function GetTWetBulbFromHumRatio(ByVal TDryBulb As Variant, ByVal HumRatio As Va MyMsgBox ("Humidity ratio cannot be negative") GoTo ErrHandler End If + BoundedHumRatio = max(HumRatio, MIN_HUM_RATIO) - TDewPoint = GetTDewPointFromHumRatio(TDryBulb, HumRatio, Pressure) + TDewPoint = GetTDewPointFromHumRatio(TDryBulb, BoundedHumRatio, Pressure) ' Initial guesses TWetBulbSup = TDryBulb @@ -635,21 +695,27 @@ Function GetTWetBulbFromHumRatio(ByVal TDryBulb As Variant, ByVal HumRatio As Va ' Bisection loop Tol = GetTol() - While (TWetBulbSup - TWetBulbInf > Tol) + index = 0 + While ((TWetBulbSup - TWetBulbInf) > Tol) - ' Compute humidity ratio at temperature Tstar - Wstar = GetHumRatioFromTWetBulb(TDryBulb, TWetBulb, Pressure) + ' Compute humidity ratio at temperature Tstar + Wstar = GetHumRatioFromTWetBulb(TDryBulb, TWetBulb, Pressure) - ' Get new bounds - If (Wstar > HumRatio) Then - TWetBulbSup = TWetBulb - Else - TWetBulbInf = TWetBulb - End If + ' Get new bounds + If (Wstar > BoundedHumRatio) Then + TWetBulbSup = TWetBulb + Else + TWetBulbInf = TWetBulb + End If + + ' New guess of wet bulb temperature + TWetBulb = (TWetBulbSup + TWetBulbInf) / 2 - ' New guess of wet bulb temperature - TWetBulb = (TWetBulbSup + TWetBulbInf) / 2 + If (index > MAX_ITER_COUNT) Then + GoTo ErrHandler + End If + index = index + 1 Wend GetTWetBulbFromHumRatio = TWetBulb @@ -675,7 +741,7 @@ Function GetHumRatioFromTWetBulb(ByVal TDryBulb As Variant, ByVal TWetBulb As Va ' Reference: ' ASHRAE Handbook - Fundamentals (2017) ch. 1 eqn 33 and 35 - Dim Wsstar As Variant + Dim Wsstar As Variant, HumRatio As Variant Wsstar = GetSatHumRatio(TWetBulb, Pressure) On Error GoTo ErrHandler @@ -687,17 +753,19 @@ Function GetHumRatioFromTWetBulb(ByVal TDryBulb As Variant, ByVal TWetBulb As Va If isIP() Then If (TWetBulb >= 32) Then - GetHumRatioFromTWetBulb = ((1093 - 0.556 * TWetBulb) * Wsstar - 0.24 * (TDryBulb - TWetBulb)) / (1093 + 0.444 * TDryBulb - TWetBulb) + HumRatio = ((1093 - 0.556 * TWetBulb) * Wsstar - 0.24 * (TDryBulb - TWetBulb)) / (1093 + 0.444 * TDryBulb - TWetBulb) Else - GetHumRatioFromTWetBulb = ((1220 - 0.04 * TWetBulb) * Wsstar - 0.24 * (TDryBulb - TWetBulb)) / (1220 + 0.444 * TDryBulb - 0.48 * TWetBulb) + HumRatio = ((1220 - 0.04 * TWetBulb) * Wsstar - 0.24 * (TDryBulb - TWetBulb)) / (1220 + 0.444 * TDryBulb - 0.48 * TWetBulb) End If Else If (TWetBulb >= 0) Then - GetHumRatioFromTWetBulb = ((2501 - 2.326 * TWetBulb) * Wsstar - 1.006 * (TDryBulb - TWetBulb)) / (2501 + 1.86 * TDryBulb - 4.186 * TWetBulb) + HumRatio = ((2501 - 2.326 * TWetBulb) * Wsstar - 1.006 * (TDryBulb - TWetBulb)) / (2501 + 1.86 * TDryBulb - 4.186 * TWetBulb) Else - GetHumRatioFromTWetBulb = ((2830# - 0.24 * TWetBulb) * Wsstar - 1.006 * (TDryBulb - TWetBulb)) / (2830# + 1.86 * TDryBulb - 2.1 * TWetBulb) + HumRatio = ((2830 - 0.24 * TWetBulb) * Wsstar - 1.006 * (TDryBulb - TWetBulb)) / (2830 + 1.86 * TDryBulb - 2.1 * TWetBulb) End If End If + ' Validity check. + GetHumRatioFromTWetBulb = max(HumRatio, MIN_HUM_RATIO) Exit Function ErrHandler: @@ -850,6 +918,8 @@ Function GetHumRatioFromVapPres(ByVal VapPres As Variant, ByVal Pressure As Vari ' Reference: ' ASHRAE Handbook - Fundamentals (2017) ch. 1 eqn 20 ' + Dim HumRatio As Variant + On Error GoTo ErrHandler If VapPres < 0 Then @@ -857,7 +927,9 @@ Function GetHumRatioFromVapPres(ByVal VapPres As Variant, ByVal Pressure As Vari GoTo ErrHandler End If - GetHumRatioFromVapPres = 0.621945 * VapPres / (Pressure - VapPres) + HumRatio = 0.621945 * VapPres / (Pressure - VapPres) + ' Validity check. + GetHumRatioFromVapPres = max(HumRatio, MIN_HUM_RATIO) Exit Function ErrHandler: @@ -880,7 +952,7 @@ Function GetVapPresFromHumRatio(ByVal HumRatio As Variant, ByVal Pressure As Var ' ASHRAE Handbook - Fundamentals (2017) ch. 1 eqn 20 solved for pw ' - Dim VapPres As Variant + Dim VapPres As Variant, BoundedHumRatio As Variant On Error GoTo ErrHandler @@ -888,8 +960,9 @@ Function GetVapPresFromHumRatio(ByVal HumRatio As Variant, ByVal Pressure As Var MyMsgBox ("Humidity ratio is negative") GoTo ErrHandler End If + BoundedHumRatio = max(HumRatio, MIN_HUM_RATIO) - VapPres = Pressure * HumRatio / (0.621945 + HumRatio) + VapPres = Pressure * BoundedHumRatio / (0.621945 + BoundedHumRatio) GetVapPresFromHumRatio = VapPres Exit Function @@ -959,7 +1032,7 @@ Function GetHumRatioFromSpecificHum(ByVal SpecificHum As Variant) As Variant End If HumRatio = SpecificHum / (1.0 - SpecificHum) - GetHumRatioFromSpecificHum = HumRatio + GetHumRatioFromSpecificHum = max(HumRatio, MIN_HUM_RATIO) Exit Function ErrHandler: @@ -1213,12 +1286,13 @@ Function GetSatHumRatio(ByVal TDryBulb As Variant, ByVal Pressure As Variant) As ' Reference: ' ASHRAE Handbook - Fundamentals (2017) ch. 1 eqn 36, solved for W ' - Dim SatVaporPres As Variant + Dim SatVaporPres As Variant, SatHumRatio As Variant On Error GoTo ErrHandler SatVaporPres = GetSatVapPres(TDryBulb) - GetSatHumRatio = 0.621945 * SatVaporPres / (Pressure - SatVaporPres) + SatHumRatio = 0.621945 * SatVaporPres / (Pressure - SatVaporPres) + GetSatHumRatio = max(SatHumRatio, MIN_HUM_RATIO) Exit Function ErrHandler: @@ -1307,6 +1381,8 @@ Function GetDegreeOfSaturation(ByVal TDryBulb As Variant, ByVal HumRatio As Vari ' ' Notes: ' This definition is absent from the 2017 Handbook. Using 2009 version instead. +' + Dim BoundedHumRatio As Variant On Error GoTo ErrHandler @@ -1314,8 +1390,9 @@ Function GetDegreeOfSaturation(ByVal TDryBulb As Variant, ByVal HumRatio As Vari MyMsgBox ("Humidity ratio is negative") GoTo ErrHandler End If + BoundedHumRatio = max(HumRatio, MIN_HUM_RATIO) - GetDegreeOfSaturation = HumRatio / GetSatHumRatio(TDryBulb, Pressure) + GetDegreeOfSaturation = BoundedHumRatio / GetSatHumRatio(TDryBulb, Pressure) Exit Function ErrHandler: @@ -1337,17 +1414,20 @@ Function GetMoistAirEnthalpy(ByVal TDryBulb As Variant, ByVal HumRatio As Varian ' Reference: ' ASHRAE Handbook - Fundamentals (2017) ch. 1 eqn 30 ' + Dim BoundedHumRatio As Variant + On Error GoTo ErrHandler If (HumRatio < 0) Then MyMsgBox ("Humidity ratio is negative") GoTo ErrHandler End If + BoundedHumRatio = max(HumRatio, MIN_HUM_RATIO) If (isIP()) Then - GetMoistAirEnthalpy = 0.24 * TDryBulb + HumRatio * (1061 + 0.444 * TDryBulb) + GetMoistAirEnthalpy = 0.24 * TDryBulb + BoundedHumRatio * (1061 + 0.444 * TDryBulb) Else - GetMoistAirEnthalpy = (1.006 * TDryBulb + HumRatio * (2501 + 1.86 * TDryBulb)) * 1000 + GetMoistAirEnthalpy = (1.006 * TDryBulb + BoundedHumRatio * (2501 + 1.86 * TDryBulb)) * 1000 End If Exit Function @@ -1375,17 +1455,20 @@ Function GetMoistAirVolume(ByVal TDryBulb As Variant, ByVal HumRatio As Variant, ' In IP units, R_DA_IP / 144 equals 0.370486 which is the coefficient appearing in eqn 26 ' The factor 144 is for the conversion of Psi = lb/in² to lb/ft². ' + Dim BoundedHumRatio As Variant + On Error GoTo ErrHandler If (HumRatio < 0) Then MyMsgBox ("Humidity ratio is negative") GoTo ErrHandler End If + BoundedHumRatio = max(HumRatio, MIN_HUM_RATIO) If (isIP()) Then - GetMoistAirVolume = R_DA_IP * GetTRankineFromTFahrenheit(TDryBulb) * (1 + 1.607858 * HumRatio) / (144 * Pressure) + GetMoistAirVolume = R_DA_IP * GetTRankineFromTFahrenheit(TDryBulb) * (1 + 1.607858 * BoundedHumRatio) / (144 * Pressure) Else - GetMoistAirVolume = R_DA_SI * GetTKelvinFromTCelsius(TDryBulb) * (1 + 1.607858 * HumRatio) / Pressure + GetMoistAirVolume = R_DA_SI * GetTKelvinFromTCelsius(TDryBulb) * (1 + 1.607858 * BoundedHumRatio) / Pressure End If Exit Function @@ -1409,7 +1492,7 @@ Function GetMoistAirDensity(ByVal TDryBulb As Variant, ByVal HumRatio As Variant ' Reference: ' ASHRAE Handbook - Fundamentals (2017) ch. 1 eqn 11 ' - Dim MoistAirVolume As Variant + Dim MoistAirVolume As Variant, BoundedHumRatio As Variant On Error GoTo ErrHandler @@ -1417,9 +1500,10 @@ Function GetMoistAirDensity(ByVal TDryBulb As Variant, ByVal HumRatio As Variant MyMsgBox ("Humidity ratio is negative") GoTo ErrHandler End If + BoundedHumRatio = max(HumRatio, MIN_HUM_RATIO) - MoistAirVolume = GetMoistAirVolume(TDryBulb, HumRatio, Pressure) - GetMoistAirDensity = (1 + HumRatio) / MoistAirVolume + MoistAirVolume = GetMoistAirVolume(TDryBulb, BoundedHumRatio, Pressure) + GetMoistAirDensity = (1 + BoundedHumRatio) / MoistAirVolume Exit Function ErrHandler: diff --git a/tests/js/test_psychrolib_ip.js b/tests/js/test_psychrolib_ip.js index 439c560f..d798bc01 100644 --- a/tests/js/test_psychrolib_ip.js +++ b/tests/js/test_psychrolib_ip.js @@ -91,6 +91,15 @@ it('test_VapPres_TDewPoint', function () { expect(psyjs.GetTDewPointFromVapPres(140.0, VapPres)).to.be.closeTo(122.0, 0.001) }); +// Test that the NR in GetTDewPointFromVapPres converges. +// This test was known problem in versions of PsychroLib <= 2.0.0 +it('test_GetTDewPointFromVapPres_convergence', function () { + for (var TDryBulb = -148; TDryBulb <= 392; TDryBulb += 1) + for (var RelHum = 0; RelHum <= 1; RelHum += 0.1) + for (var Pressure = 8.6; Pressure <= 17.4; Pressure += 1) + psyjs.GetTWetBulbFromRelHum(TDryBulb, RelHum, Pressure) +}); + // Test of relationships between humidity ratio and vapour pressure it('test_HumRatio_VapPres', function () { var HumRatio = psyjs.GetHumRatioFromVapPres(0.45973, 14.175) // conditions at 77 F, std atm pressure at 1000 ft @@ -121,6 +130,9 @@ it('test_HumRatio_TWetBulb', function () { checkRelDiff(HumRatio,0.00114657481090184, 0.0003) var TWetBulb = psyjs.GetTWetBulbFromHumRatio(30.2, HumRatio, 14.1751) expect(TWetBulb).to.be.closeTo(23.0, 0.001) + + // Low HumRatio -- this should evaluate true as we clamp the HumRation to 1e-07. + expect(psyjs.GetTWetBulbFromHumRatio(23,1e-09,95461)).to.equal(psyjs.GetTWetBulbFromHumRatio(23,1e-07,95461)) }); /** diff --git a/tests/js/test_psychrolib_si.js b/tests/js/test_psychrolib_si.js index c24e5cca..9a637a6e 100644 --- a/tests/js/test_psychrolib_si.js +++ b/tests/js/test_psychrolib_si.js @@ -91,6 +91,23 @@ it('test_VapPres_TDewPoint', function () { expect(psyjs.GetTDewPointFromVapPres(60.0, VapPres)).to.be.closeTo(50.0, 0.001) }); +// Test of relationships between wet bulb temperature and relative humidity +// This test was known to cause a convergence issue in GetTDewPointFromVapPres +// in versions of PsychroLib <= 2.0.0 +it('test_TWetBulb_RelHum', function () { + var TWetBulb = psyjs.GetTWetBulbFromRelHum(7, 0.61, 100000) + checkRelDiff(TWetBulb, 3.92667433781955, 0.001) +}); + +// Test that the NR in GetTDewPointFromVapPres converges. +// This test was known problem in versions of PsychroLib <= 2.0.0 +it('test_GetTDewPointFromVapPres_convergence', function () { + for (var TDryBulb = -100; TDryBulb <= 200; TDryBulb += 1) + for (var RelHum = 0; RelHum <= 1; RelHum += 0.1) + for (var Pressure = 60000; Pressure <= 120000; Pressure += 10000) + psyjs.GetTWetBulbFromRelHum(TDryBulb, RelHum, Pressure) +}); + // Test of relationships between humidity ratio and vapour pressure it('test_HumRatio_VapPres', function () { var HumRatio = psyjs.GetHumRatioFromVapPres(3169.7, 95461) // conditions at 25 C, std atm pressure at 500 m @@ -121,6 +138,9 @@ it('test_HumRatio_TWetBulb', function () { checkRelDiff(HumRatio, 0.00120399819933844, 0.0003) var TWetBulb = psyjs.GetTWetBulbFromHumRatio(-1, HumRatio, 95461) expect(TWetBulb).to.be.closeTo(-5, 0.001) + + // Low HumRatio -- this should evaluate true as we clamp the HumRation to 1e-07. + expect(psyjs.GetTWetBulbFromHumRatio(-5,1e-09,95461)).to.equal(psyjs.GetTWetBulbFromHumRatio(-5,1e-07,95461)) }); /** diff --git a/tests/test_psychrolib_ip.py b/tests/test_psychrolib_ip.py index aed5efc3..39319e0e 100644 --- a/tests/test_psychrolib_ip.py +++ b/tests/test_psychrolib_ip.py @@ -1,6 +1,7 @@ # Copyright (c) 2018 D. Thevenard and D. Meyer. Licensed under the MIT License. # Test of PsychroLib in IP units for Python, C, and Fortran. +import numpy as np import pytest pytestmark = pytest.mark.usefixtures('SetUnitSystem_IP') @@ -28,6 +29,18 @@ def test_GetSatVapPres(psy): assert psy.GetSatVapPres(212) == pytest.approx(14.7094, rel = 0.0003) assert psy.GetSatVapPres(300) == pytest.approx(67.0206, rel = 0.0003) +# Test that the NR in GetTDewPointFromVapPres converges. +# This test was known problem in versions of PsychroLib <= 2.0.0 +def test_GetTDewPointFromVapPres_convergence(psy): + TDryBulb = np.arange(-148, 392, 1) + RelHum = np.arange(0, 1, 0.1) + Pressure = np.arange(8.6, 17.4, 1) + for T in TDryBulb: + for RH in RelHum: + for p in Pressure: + psy.GetTWetBulbFromRelHum(T, RH, p) + print('GetTDewPointFromVapPres converged') + # Test saturation humidity ratio # The values are tested against those published in Table 2 of ch. 1 of the 2017 ASHRAE Handbook - Fundamentals # Agreement is not terrific - up to 2% difference with the values published in the table @@ -97,6 +110,8 @@ def test_HumRatio_TWetBulb(psy): assert HumRatio == pytest.approx(0.00114657481090184, rel = 0.0003) TWetBulb = psy.GetTWetBulbFromHumRatio(30.2, HumRatio, 14.1751) assert TWetBulb == pytest.approx(23.0, abs = 0.001) + # Low HumRatio -- this should evaluate true as we clamp the HumRation to 1e-07. + assert psy.GetTWetBulbFromHumRatio(25,1e-09,95461) == psy.GetTWetBulbFromHumRatio(25,1e-07,95461) ############################################################################### diff --git a/tests/test_psychrolib_si.py b/tests/test_psychrolib_si.py index cf461ad2..6759b710 100644 --- a/tests/test_psychrolib_si.py +++ b/tests/test_psychrolib_si.py @@ -1,6 +1,7 @@ # Copyright (c) 2018 D. Thevenard and D. Meyer. Licensed under the MIT License. # Test of PsychroLib in SI units for Python, C, and Fortran. +import numpy as np import pytest pytestmark = pytest.mark.usefixtures('SetUnitSystem_SI') @@ -68,6 +69,25 @@ def test_VapPres_TDewPoint(psy): VapPres = psy.GetVapPresFromTDewPoint(50.0) assert psy.GetTDewPointFromVapPres(60.0, VapPres) == pytest.approx(50.0, abs = 0.001) +# Test of relationships between wet bulb temperature and relative humidity +# This test was known to cause a convergence issue in GetTDewPointFromVapPres +# in versions of PsychroLib <= 2.0.0 +def test_TWetBulb_RelHum(psy): + TWetBulb = psy.GetTWetBulbFromRelHum(7, 0.61, 100000) + assert TWetBulb == pytest.approx(3.92667433781955, rel = 0.001) + +# Test that the NR in GetTDewPointFromVapPres converges. +# This test was known problem in versions of PsychroLib <= 2.0.0 +def test_GetTDewPointFromVapPres_convergence(psy): + TDryBulb = np.arange(-100, 200, 1) + RelHum = np.arange(0, 1, 0.1) + Pressure = np.arange(60000, 120000, 10000) + for T in TDryBulb: + for RH in RelHum: + for p in Pressure: + psy.GetTWetBulbFromRelHum(T, RH, p) + print('GetTDewPointFromVapPres converged') + # Test of relationships between humidity ratio and vapour pressure # Humidity ratio values to test against are calculated with Excel def test_HumRatio_VapPres(psy): @@ -97,6 +117,8 @@ def test_HumRatio_TWetBulb(psy): assert HumRatio == pytest.approx(0.00120399819933844, rel = 0.0003) TWetBulb = psy.GetTWetBulbFromHumRatio(-1, HumRatio, 95461) assert TWetBulb == pytest.approx(-5, abs = 0.001) + # Low HumRatio -- this should evaluate true as we clamp the HumRation to 1e-07. + assert psy.GetTWetBulbFromHumRatio(-5,1e-09,95461) == psy.GetTWetBulbFromHumRatio(-5,1e-07,95461) ############################################################################### diff --git a/tests/vba/test_psychrolib_ip.bas b/tests/vba/test_psychrolib_ip.bas index 3d51093f..7cf6417f 100644 --- a/tests/vba/test_psychrolib_ip.bas +++ b/tests/vba/test_psychrolib_ip.bas @@ -27,6 +27,7 @@ Sub RunAllTests() Call test_GetStandardAtmTemperature Call test_SeaLevel_Station_Pressure Call test_AllPsychrometrics + Call test_GetTDewPointFromVapPres_convergence Debug.Print "# of tests run :", TestCount Debug.Print "# of issues found:", IssueCount End Sub @@ -261,3 +262,20 @@ Sub test_AllPsychrometrics() Call TestExpression("CalcPsychrometricsFromRelHum", TWetBulb, 65, abst:=0.1) End Sub + +'############################################################################## +' Test of the convergence of the NR method in GetTDewPointFromVapPres +' over a wide range of inputs +' This test was known problem in versions of PsychroLib <= 2.0.0 +'############################################################################## +' +Sub test_GetTDewPointFromVapPres_convergence() + For TDryBulb = -148 To 392 Step 2 + For RelHum = 0 To 1 Step 0.1 + For Pressure = 8.6 To 17.4 Step 0.8 + TWetBulb = GetTWetBulbFromRelHum(TDryBulb, RelHum, Pressure) + Next Pressure + Next RelHum + Next TDryBulb + Call TestExpression("GetTDewPointFromVapPres convergence test", 1, 1, abst:=0.1) +End Sub diff --git a/tests/vba/test_psychrolib_ip.xlsm b/tests/vba/test_psychrolib_ip.xlsm index 5628e945..1044a328 100755 Binary files a/tests/vba/test_psychrolib_ip.xlsm and b/tests/vba/test_psychrolib_ip.xlsm differ diff --git a/tests/vba/test_psychrolib_si.bas b/tests/vba/test_psychrolib_si.bas index 8d1efe7f..660a064b 100644 --- a/tests/vba/test_psychrolib_si.bas +++ b/tests/vba/test_psychrolib_si.bas @@ -27,6 +27,7 @@ Sub RunAllTests() Call test_GetStandardAtmTemperature Call test_SeaLevel_Station_Pressure Call test_AllPsychrometrics + Call test_GetTDewPointFromVapPres_convergence Debug.Print "# of tests run :", TestCount Debug.Print "# of issues found:", IssueCount End Sub @@ -126,6 +127,14 @@ Sub test_VapPres_TDewPoint() Call TestExpression("GetTDewPointFromVapPres", GetTDewPointFromVapPres(60#, VapPres), 50#, abst:=0.001) End Sub +' Test of relationships between wet bulb temperature and relative humidity +' This test was known to cause a convergence issue in GetTDewPointFromVapPres +' in versions of PsychroLib <= 2.0.0 +Sub test_TWetBulb_RelHum() +TWetBulb = GetTWetBulbFromRelHum(7, 0.61, 100000) +Call TestExpression("GetTWetBulbFromRelHum", TWetBulb, 3.92667433781955, relt:=0.001) +End Sub + ' Test of relationships between humidity ratio and vapour pressure ' Humidity ratio values to test against are calculated with Excel Sub test_HumRatio_VapPres() @@ -261,3 +270,21 @@ Sub test_AllPsychrometrics() Call TestExpression("CalcPsychrometricsFromRelHum", TWetBulb, 20, abst:=0.1) End Sub +'############################################################################## +' Test of the convergence of the NR method in GetTDewPointFromVapPres +' over a wide range of inputs +' This test was known problem in versions of PsychroLib <= 2.0.0 +'############################################################################## +' +Sub test_GetTDewPointFromVapPres_convergence() + For TDryBulb = -100 To 200 Step 1 + For RelHum = 0 To 1 Step 0.1 + For Pressure = 60000 To 120000 Step 10000 + TWetBulb = GetTWetBulbFromRelHum(TDryBulb, RelHum, Pressure) + Next Pressure + Next RelHum + Next TDryBulb + Call TestExpression("GetTDewPointFromVapPres convergence test", 1, 1, abst:=0.1) +End Sub + + diff --git a/tests/vba/test_psychrolib_si.xlsm b/tests/vba/test_psychrolib_si.xlsm index 38229995..15528456 100755 Binary files a/tests/vba/test_psychrolib_si.xlsm and b/tests/vba/test_psychrolib_si.xlsm differ