#!/usr/bin/env python2 # -*- coding: utf-8 -*- """ psychro.py Utility routines for calculations related to psychrometrics. See: http://en.wikipedia.org/wiki/Arden_Buck_equation and Buck, A. L. (1981), "New equations for computing vapor pressure and enhancement factor", J. Appl. Meteorol. 20: 1527–1532 linked from the Wikipedia article. """ from calculable import Calculable, CalculableAttribute import math def absolute(T): """ Convert temperature in °C to K. """ return T + 273.15 def wexler(T, supercooled=False, pressure=100e3): """ Calculate the equilibrium ("saturated") vapour pressure of water at given temperature T. T - Degrees Celcius. supercooled - Pass True to use the over-water equation even when temperature is below 0°C. pressure - Air pressure in Pa. Pass None for pure water vapour (no nitrogen, oxygen, etc). Result is in pascals. (Note, above reference gives mb.) """ # Check in applicable temperature range. if T > 50.0: raise ValueError('Too hot (greater than 50°C): %s°C' % (str(T),)) if T < -80.0: raise ValueError('Too cold (less than -80°C): %s°C' % (str(T),)) if (T < -40.0) and supercooled: raise ValueError('Too cold for supercooled water (less than -40°C): %s°C' % (str(T),)) # Calculate pressure for pure water vapour (no air). T0 = absolute(T) poly = [T0 ** i for i in range(-2, 5)] + [math.log(T0)] if (T >= 0.0) or supercooled: # Over water coeff = [-2991.2729, -6017.0128, 18.87643854, -0.028354721, 0.17838301e-4, -0.84150417e-9, 0.44412543e-12, 2.858487] else: # Over ice coeff = [0.0, -5865.36960, 22.241033, 0.013749042, -0.34031775e-4, 0.26967687e-7, 0.0, 0.6918651] e = math.exp(sum(map(float.__mul__, coeff, poly))) # Apply small correction (enhancement factor) for air pressure, if required. if pressure != None: if (T >= 0.0) or supercooled: class f: # Arden Buck table 3, fw5 A = 4.1e-4 B = 3.48e-6 C = 7.4e-10 D = 30.6 E = -3.8e-2 else: class f: # Arden Buck table 3, fi5 A = 4.8e-4 B = 3.47e-6 C = 5.9e-10 D = 23.8 E = -3.1e-2 P = pressure / 100.0 ft = T + f.D + (f.E * P) f = 1 + f.A + (P * (f.B + f.C * (ft * ft))) e = e * f return e def testWexler(): """ Test function wexler against data in Arden Buck table 1. """ data = { -80: { 'ei': 0.0005481 }, -70: { 'ei': 0.0026189 }, -60: { 'ei': 0.010820, 'emi': 0.01089 }, -50: { 'ei': 0.039402, 'emi': 0.03963 }, -40: { 'ew': 0.19047, 'ei': 0.12849, 'emw': 0.19146, 'emi': 0.12915 }, -30: { 'ew': 0.5106, 'ei': 0.38024, 'emw': 0.5130, 'emi': 0.38203 }, -20: { 'ew': 1.2563, 'ei': 1.0328, 'emw': 1.2618, 'emi': 1.0373 }, -10: { 'ew': 2.8657, 'ei': 2.5992, 'emw': 2.8774, 'emi': 2.6098 }, 0: { 'ew': 6.1121, 'ei': 6.1115, 'emw': 6.1360, 'emi': 6.1360 }, 10: { 'ew': 12.279, 'emw': 12.327 }, 20: { 'ew': 23.385, 'emw': 23.479 }, 30: { 'ew': 42.452, 'emw': 42.633 }, 40: { 'ew': 73.813, 'emw': 74.157 }, 50: { 'ew': 123.45, 'emw': 124.09 } } for t, d in data.items(): for k, v in d.items(): if k == 'ew': p = wexler(t, supercooled=True, pressure=None) elif k == 'ei': p = wexler(t, supercooled=False, pressure=None) elif k == 'emw': p = wexler(t, supercooled=True, pressure=100e3) elif k == 'emi': p = wexler(t, supercooled=False, pressure=100e3) else: raise ValueError('testWexler: bad key') v *= 100.0 if (p < (v * 0.9999)) or ((v * 1.001) < p): raise ValueError('wexler fail: %f, %s, %f, %f' % (t, k, v, p)) class GasMass(Calculable): """ Representation of a blob of ideal gas which obeys the gas law: pV = nRT where: p = pressure V = volume n = quantity (number of moles) R = ideal gas constant T = absolute temperature The mass of the gas is the quantity times the molar mass. Molar mass is defined by sub-classes. """ R = 8.314472 # Ideal gas constant, J·mol⁻¹·K⁻¹ _ = CalculableAttribute attributes = ( _('p', 'Pressure', 'Pa') (lambda s: s.n * s.R * s.T / s.V), _('V', 'Volume', 'm³') (lambda s: s.n * s.R * s.T / s.p), _('n', 'Quantity', 'mol') (lambda s: s.p * s.V / (s.R * s.T)) (lambda s: s.m / s.molarMass), _('t', 'Temperature', '°C') (lambda s: s.T - 273.15), _('T', 'Absolute temperature', 'K') (lambda s: s.t + 273.15), _('m', 'Mass', 'kg') (lambda s: s.V * s.den) (lambda s: s.n * s.molarMass), _('den', 'Density', 'kg/m³') (lambda s: s.m / s.V), _('H', 'Enthalpy', 'J') (lambda s: (s.m * s.t * s.specificHeatCapacity) + (s.m * s.latentHeatOfVaporization)) ) del _ class DryAirMass(GasMass): """ Air (nitrogen, oxygen, etc) without water vapour. """ molarMass = 0.02897 # kg/mol specificHeatCapacity = 1005.0 # J/(kg·K), constant pressure, just above freezing, # http://www.engineeringtoolbox.com/air-properties-d_156.html latentHeatOfVaporization = 0.0 # If we're condensing oxygen, nitrogen, etc, in the house... class WaterVapourMass(GasMass): """ Gaseous water on its own. """ molarMass = 0.018 # kg/mol specificHeatCapacity = 1859.0 # J/(kg·K), constant pressure, just above freezing, # http://www.engineeringtoolbox.com/water-vapor-d_979.html latentHeatOfVaporization = 2470e3 # J/kg, approx value for around 10 °C class AirMass(Calculable): """ Mixture of dry air and water vapour. """ _ = CalculableAttribute attributes = ( _('p', 'Pressure', 'Pa') (lambda s: s.da.p + s.wv.p), _('V', 'Volume', 'm³') (lambda s: s.da.V) (lambda s: s.wv.V) (lambda s: s.n * s.R * s.T / s.p), _('t', 'Temperature', '°C') (lambda s: s.da.t), _('m', 'Mass', 'kg') (lambda s: s.da.m + s.wv.m) (lambda s: s.V * s.den) (lambda s: s.n * s.molarMass), _('den', 'Density', 'kg/m³') (lambda s: s.m / s.V), _('svp', 'Saturated Vapour Pressure', 'Pa') (lambda s: wexler(s.wv.t, s.p)), _('rh', 'Relative humidity') (lambda s: s.wv.p / s.svp), _('H', 'Enthalpy', 'J') (lambda s: s.da.H + s.wv.H), _('da', 'Dry Air', cls=DryAirMass). var('p') (lambda s: s.p - s.wv.p). var('V') (lambda s: s.V). var('t') (lambda s: s.t) (), _('wv', 'Water Vapour', cls=WaterVapourMass). var('p') (lambda s: s.rh * s.svp). var('V') (lambda s: s.V). var('t') (lambda s: s.t) () ) del _ if __name__ == '__main__': testWexler()