Saturation fallback

[vegardjervell]

I've been a bit annoyed sometimes that safe_bubP and friends would sometimes fail, even though the envelope traces worked (specifically for mixtures). I also noticed that safe_bubP and friends sometimes ended up returning points well inside the two-phase region without indicating an error.

To remedy this, I've implemented a fallback solver that is used if either

• safe_[dew/bub][P/T] fails, or
• The value computed is in the two phase region (within a tolerance) on the phase fractions

The solver is a naive bracket solver that uses twoPhaseTPflash to locate the saturation line, so in order to make sure the bracket solver didn't force termination if twoPhaseTPflash failed, I also had to make twoPhaseTPflash_safe, which is the old twoPhaseTPflash, except it returns with an error flag instead of terminating if it hits an error. The subroutine twoPhaseTPflash calls twoPhaseTPflash_safe and terminates exactly as it used to if it gets an error.

On the positive side, safe_bubP and friends will now never silently return points deep in the two-phase region, and will successfully find the saturation line in a lot of situations where it previously failed. The caveat is that all calls to safe_bubP and friends will incur a call to twoPhaseTPflash for multicomponent mixtures, which costs a little runtime. This cost could maybe be reduced by using something from the stability module to check the solution instead of a full TP-flash.

I've attached some figures to illustrate the issue as it was, and that the fallback solver handles them. Note that the fallback solver is never called if the existing saturation functions successfully located the saturation point, or if the system is a pure component.

TP_no fallback.pdf
TP_t_no fallback.pdf
TP_t_with fallback.pdf
TP_with fallback.pdf
TXY_no fallback.pdf
TXY_with fallback.pdf

[vegardjervell]

@morteham I've added some checks that ensures the bracket solver will terminate either if no two-phase region is found or if no single-phase region is found. Tested that the PXY and TXY envelopes for N2/CO2 can be produced using the saturation point functions (i.e. that the solvers don't get stuck for too long in the cases where the envelope doesn't close).

As far as I'm concerned this can be merged now. Just let me know if there's something that you think should be modified or tested more extensively.