Some Wet-bulb related calculations Rutger Kassies, 06-12-2012¶
In [1]:
# function declaration
def vp(t): # vapor pressure, t in celcius
vp = 6.108*np.exp((17.27*t)/(237.3+t))
return vp
def ea2td(ea): # eact to dewpoint, ea in hPa
td = (237.3*np.log(ea/6.108))/(17.27-np.log(ea/6.108))
return td
def psy(lv, p): # psychometric constant
# 1004 = specific heat
# 0.622 = ratio of molecular weight
psy = (1004*p)/(lv*0.622)
return psy
def lv(t): # latent heat of vaporization, t in celcius
lv = 1000*(2501-2.361*t)
return lv
def wb(ta, tw, psyc, td):
# function to calculate the wet-bulb by iterating.
ea_ta = vp(td)
i = 1
while i < 50: # would be better to check for convergence
# calculate the eact from the wet-bult, this should be equal to ea_ta
ea_tw = vp(tw) - psyc * ( ta - tw)
# so use the difference to adjust twet
dea = ea_ta - ea_tw
# emperical adjustment, might never converge at extreme values!
tw = tw + (dea * (0.5 - (ea_ta/450)))
i = i + 1
return tw
In [2]:
# declare arrays for evaluation, reshape to allow breadcasting
tair = np.linspace(-10,5,151).reshape(151,1)
RH = np.linspace(1,100,100).reshape(1,100)
RH_grd, tair_grd = np.meshgrid(RH.ravel(), tair.ravel())
# calc actual and saturated vapor pressure
esat = vp(tair)
eact = esat * ( RH / 100.)
# calc dew-point
td = ea2td(eact)
# calc psychometric constant at standard pressure
psyc = psy(lv(tair), 1013.25)
# calculate initial tw to start iteration
tw_ini = ( tair + td ) / 2.
# calc tw
tw = wb(tair, tw_ini, psyc, td)
In [3]:
fig, ax = plt.subplots(2, 2, figsize=(14,14))
levels=np.arange(-20,22,2)
mpl.rcParams['contour.negative_linestyle'] = 'solid'
loc=(0,0)
cs1 = ax[loc].contourf(tair_grd, td, tw, levels=levels, cmap=mpl.cm.Blues)
cs1_line = ax[loc].contour(tair_grd, td, tw, levels=levels, colors='k')
ax[loc].clabel(cs1_line, inline=1, fontsize=11, fmt='%1.0f')
loc=(0,1)
cs1 = ax[loc].contourf(tair_grd,RH_grd, tw, levels=levels, cmap=mpl.cm.Blues)
cs1_line = ax[loc].contour(tair_grd,RH_grd, tw, levels=levels, colors='k')
ax[loc].clabel(cs1_line, inline=1, fontsize=11, fmt='%1.0f')
loc=(1,0)
cs1 = ax[loc].contourf(tair_grd,RH_grd, (tair_grd + td) /2. , levels=levels, cmap=mpl.cm.Blues)
cs1_line = ax[loc].contour(tair_grd,RH_grd, (tair_grd + td) /2., levels=levels, colors='k')
ax[loc].clabel(cs1_line, inline=1, fontsize=11, fmt='%1.0f')
errlevels=[-20,-15,-10,-8,-6,-4,-2,-1,-.5,-.25,0]
loc=(1,1)
cs1 = ax[loc].contourf(tair_grd,RH_grd, ((tair_grd + td) /2.)-tw , levels=errlevels, cmap=mpl.cm.Reds_r)
cs1_line = ax[loc].contour(tair_grd,RH_grd, ((tair_grd + td) /2.)-tw, levels=errlevels, colors='k')
ax[loc].clabel(cs1_line, inline=1, fontsize=11, fmt='%1.2f')
graph_props = [(True, (-10,5), (-40,5), 'Ta [$^\circ$C]', 'Td [$^\circ$C]','Twet' ),\
(True, (-10,5), (0,100), 'Ta [$^\circ$C]', 'RH [%]','Twet' ),\
(True, (-10,5), (0,100), 'Ta [$^\circ$C]', 'RH [%]','(Ta + Td) / 2' ),\
(True, (-10,5), (0,100), 'Ta [$^\circ$C]', 'RH [%]','error of (Ta + Td) / 2' )]
for n, axs in enumerate(ax.ravel()):
grid, xlim, ylim, xlabel, ylabel, title = graph_props[n]
axs.grid(grid)
axs.set_xlim(xlim)
axs.set_ylim(ylim)
axs.set_xlabel(xlabel, fontsize=10)
axs.set_ylabel(ylabel, fontsize=10)
axs.set_title(title, fontsize=15)
axs.tick_params(labelsize=10, axis='both')
