Numerical issue in water saturation routine

[jranalli]

Found a case where the numbers go a bit crazy. This is nowhere near the critical point (which is 220 bar). Expected value is around 592 K, which is also not especially near the critical point.

import pyromat as pm
h2o = pm.get('mp.H2O')
h2o.Ts(p=110) # yields nan

[chmarti1]

Confirmed. This only occurs in a tight numerical range around 110bar. I can recreate the failure for 110, 111, and 112 bar. I get valid values everywhere else.

[jranalli]

MORE INFO!

At that particular pressure, the iteration just randomly happens to result in an iteration guess around 648 for temperature that turns out to be above the critical point. That's leading _ps() to yield imaginary values.

So I believe the error comes in because you're allowing the temperature limits in _Ts() to go above the critical point on line 1952. I think the imaginary values come from the fact that the polynomial likely can't ever actually exceed the critical point, since the slope becomes horizontal there. Likely if you make Tc a hard cap it will make it harder to converge as you approach the critical point, but will at least not yield A+Bj's.

[chmarti1]

Confirmed. The _Ts() algorithm allows iteration to wander 1% above critical and 1% below triple temperature. The original thinking was that even if the underlying polynomial is meaningless outside of those bounds, it should still be continuous. That approach is flawed, because the formulation uses fractional powers of a term that becomes negative above the critical point. I can correct this issue by enforcing a strict limit against guesses outside the critical point. In my tests, that works fine.

The only downside is that the saturation temperature will fail when precisely the critical pressure is given. Before I release the fix, I'll see if it makes sense to check for those values and coerce them to the known critical temperature. I think that will work fine.

[chmarti1]

This issue is corrected in version 2.2.1. I eliminated the 1% expansion on the iteration boundaries on the high side. It completely mitigates the problem. I have left the 1% expansion on the low side, so the triple point is still inside of the valid iteration domain even with numerical error.