Relative Humidity is not calculated properly in iapws.humidAir. (added more info)

[che0815]

When using iapws.humidAir it returns incorrect values for RH which where used to initialise the state in the first place, check code below.
Here the relative humidity is used to calculate the molar fraction for the initialisation of hum_air_state.
However hum_air_state.RH does not return the proper phi.

import iapws

T_degC = 14
T_K    = T_degC + 273.15
  
p_mbar = 1000# 013.25
p_bar  = p_mbar * 1e-3
p_Pa   = p_mbar * 1e+2
p_MPa  = p_Pa   * 1e-6

phi_p  = 60
phi_r  = phi_p * 1e-2

p_sat_w = iapws.iapws97._PSat_T(T_K) #Temperature, [K]; Returns:Pressure, [MPa]

p_w = phi_r * p_sat_w
xa = (p_MPa - p_w)/p_MPa

hum_air_state = iapws.humidAir.HumidAir(T=T_K, P=p_MPa, xa=xa)

print(phi_r) #OUT 0.6
print(hum_air_state.RH) #16.210589116125586

[jjgomera]

Hi,
Your input calculation has a concept error:

p_sat_w = iapws.iapws97._PSat_T(T_K) #Temperature, [K]; Returns:Pressure, [MPa]

This line calculate the saturation pressure for steam water over liquid water without air, the resulted value is very high and because that the RH value returned in calculation is 16 (1600% relative humidity)

dry_air = iapws.humidAir.HumidAir(T=287.15, P=0.1, xa=1)
print(xw, 1-dry_air.xa_sat)
0.009593664339180234 0.000591814663271

As you can see the air at 14ºC and 0.1MPa can be 0.0591% of water at saturation state, but you are forcing a water composition of 0.9%. so you're imposing a humid air with a 1600% of relative humidity

[jjgomera]

If you are interested in psychrometric calculation I do recommend you this tools in my other project pychemqt: https://github.com/jjgomera/pychemqt/blob/master/tools/UI_psychrometry.py

Documentation: https://pychemqt.readthedocs.io/en/latest/tools.UI_psychrometry.html

[che0815]

sorry, I did not realise (until now) where I need to look for a response to this/my issue.
You are certainly right I am making a conceptual error. What I actually intended to do is to calculate the fluid properties of humid air (R, gamma, cp). I was therefore looking for the partial pressure of the water vapour. And in order to calculate the partial pressure of water vapour in humid air, I needed the saturation pressure of water.

Am I reading it correctly that you are saying I cannot use p_w = phi_r * p_sat_w in order to calculate the partial pressure of water vapour in humid air?

[che0815]

using coolprop and magnus in comparison gives similar values:
1027.6441580839664 < partial pwv from coolprop
1023.4469246354241 < partial pwv from iapws
1021.0032144355836 < partial pwv from magnus

so I think this is not where the error is.

[che0815]

so the only mistake might be the calculation of xa then:

xa = (p_MPa - p_wv)/p_MPa

I expected xa to be the mole fraction of dry air in humid air, and calculated it via the partial pressures. Not sure what I should have done differently.

[che0815]

using above equation and xw = 1-xa, which aligns nicely with coolprop again
0.010142059295178548 < psi_w from coolprop, HAPropsSI('psi_w','P',p_Pa,'H',h,'R',phi_r)
0.010100635821716497 < psi_w via pwv from iapws

so to me it still seems the problem (might as well be mine) is located in the definition of the state.
iapws.humidAir.HumidAir(T=T_K, P=p_MPa, xa=xa)

[che0815]

I think I did a very bad job the first time to phrase what I am doing. Therefore, I am posting the code to reproduce the error here with comments.

TL;DR
While Saturation Pressure, Partial Pressure and Mole Fraction of IAPWS and CoolProp show the same results. I am not able to reproduce Relative Humidity (RH) as a cross check.

import modules

import iapws
from CoolProp.HumidAirProp import HAPropsSI

thermodynamic boundary conditions

T_degC = 15
T_K = T_degC + 273.15

p_mbar = 1013.25
p_bar = p_mbar * 1e-3
p_Pa = p_mbar * 1e+2
p_MPa = p_Pa * 1e-6

phi_p = 60
phi_r = phi_p * 1e-2

IAPWS

p_sat_w_iapws = iapws.iapws97._PSat_T(T_K) #Temperature, [K]; Returns:Pressure, [MPa]

p_w_iapws = phi_r * p_sat_w
xa = (p_MPa - p_w)/p_MPa
xw = (p_w)/p_MPa

Coolprop

h = HAPropsSI('H','T',T_K,'P',p_Pa,'R',phi_r)
h_sat = HAPropsSI('H','T',T_K,'P',p_Pa,'R',1.0)

p_sat_w_coolprop = HAPropsSI('P_w','P',p_Pa,'H',h_sat,'R',1.0)
p_w_coolprop = HAPropsSI('P_w','P',p_Pa,'H',h,'R',phi_r)
psi_w_coolprop = HAPropsSI('psi_w','P',p_Pa,'H',h,'R',phi_r)

#Comparison
print('Saturation Pressure')
print('{} < IAPWS'.format(p_sat_w_iapws*1e6))
print('{} < Coolprop'.format(p_sat_w_coolprop))

print('Partial Pressure')
print('{} < IAPWS'.format(p_sat_w_iapws1e6phi_r))
print('{} < Coolprop (direct)'.format(p_w_coolprop))
print('{} < Coolprop (via phi)'.format(p_sat_w_coolprop*phi_r))

print('Mole Fraction of Water Vapour')
print('{} < IAPWS (via phi)'.format((p_sat_w_iapws*phi_r)/p_MPa))
print('{} < Coolprop (direct)'.format(psi_w_coolprop))

print('Cross-Check Relative Humidity')
hum_air_state = iapws.humidAir.HumidAir(T=T_K, P=p_MPa, xw=psi_w_coolprop)
print('{} < relative Humidity (Thermodynamic Boundary Conditions)'.format(phi_r))
print('{} < hum_air_state relative Humidity (.RH)'.format(hum_air_state.RH))

[jjgomera]

Thanks for reporting, and excuse me for my late response.
You're right, the saturation value had a bug so the RH returned fail too. Try now with last commit

[che0815]

Hi, thanks for following up (and sorry for my late response as well).
I downloaded version 1.4 (conda install -c conda-forge iapws). And it seems something has changed however the sanity check still fails, or maybe I am doing something wrong.

import iapws

# setting boundary conditions
T_degC = 15
T_K    = T_degC + 273.15
  
p_mbar = 1013.25
p_bar  = p_mbar * 1e-3
p_Pa   = p_mbar * 1e+2
p_MPa  = p_Pa   * 1e-6

phi_p  = 60
phi_r  = phi_p * 1e-2

# computing saturation and partial pressure
p_sat_w_iapws = iapws.iapws97._PSat_T(T_K) #Temperature, [K]; Returns:Pressure, [MPa]

p_w_iapws = phi_r * p_sat_w_iapws
xa_iapws = (p_MPa - p_w_iapws)/p_MPa
xw_iapws = (p_w_iapws)/p_MPa

# sanity check
hum_air_state_xw = iapws.humidAir.HumidAir(T=T_K, P=p_MPa, xw=xw_iapws)
hum_air_state_xa = iapws.humidAir.HumidAir(T=T_K, P=p_MPa, xa=xa_iapws)
print('{} < relative Humidity (Thermodynamic Boundary Conditions)'.format(phi_r))
print('{} < hum_air_state relative Humidity via xw (.RH)'.format(hum_air_state_xw.RH))
print('{} < hum_air_state relative Humidity via xa (.RH)'.format(hum_air_state_xa.RH))

I get:

0.6 < relative Humidity (Thermodynamic Boundary Conditions)
-0.3775459391856669 < hum_air_state relative Humidity via xw (.RH)
-0.3775459391856669 < hum_air_state relative Humidity via xa (.RH)

[schymans]

Probably related to this, in Version 1.5.4, I get:

from iapws import HumidAir
HA = HumidAir()
HA._eq(290, 0.101325)

1.0188116920478836
It seems that any values of temperature below 297.35 return values >1, which is impossible. Seems like a bug to me. Or am I misunderstanding what _eq should calculate? It says "Saturation mass fraction of dry air in humid air [kg/kg]".

[jjgomera]

Hi schymans,
your problem is different, it is a problem with the convergence in the method, solved, I think with the last commit, try it

By the way, thanks for remembering this mistake, I totally forgot it

[schymans]

Thank you @jjgomera! This was ultra-fast, I'm glad you found and fixed the bug so quickly. The results now pass my sanity test, so we can keep going from here.

[jjgomera]

By the way, I review original bug and the sanity check work now:

0.6 < relative Humidity (Thermodynamic Boundary Conditions)
0.5974895833971026 < hum_air_state relative Humidity via xw (.RH)
0.5974895833971026 < hum_air_state relative Humidity via xa (.RH)

Really it was same bug, thanks.

[schymans]

Awesome! Just one thing, iapws.__version__ now states 1.5.2, whereas I had 1.5.4 before. Could it be that some changes done for 1.5.4 have not been merged back to master?

[jjgomera]

Are you sure you had 1.5.4? Maybe it was 1.4? The version is unchanged from 19 jun 2020 as 1.5.2

https://github.com/jjgomera/iapws/commits/master/iapws/VERSION

and this is the last used, maybe is necessary 1.5.3 for add last commits including this bug fix to pipy package