improve solver stability close to critical point

[BenPortner]

One reason for solver instabilities is that liquid enthalpies skyrocket for T>Tc. I opened an issue in CalebBell/thermo to address this. Until a more elegant solution is found, I propose to set phase='g' if T>Tc while getting stream properties like enthalpy, viscosity, etc. Using the proposed fix, the nitrogen liquefaction example in https://github.com/BioSTEAMDevelopmentGroup/biosteam/pull/121 works out of the box with CoolProp methods, without the need for a second convergence cycle.

Note: I tried to put the Tcmax calculation into Stream.__init__ to reduce the overhead but since there are 1283212837 other ways to create Streams in thermosteam (perceived reality 😜) and the Tcmax calculation has to be added to all of them, I could not get it to work.

[yoelcortes]

@BenPortner, I see. The proposed solution raises some issues with multiple components:

>>> from thermosteam import *
>>> settings.set_thermo(['Water', 'CO2'])
>>> s = Stream(Water=2, CO2=0.00000001, T=320)
>>> s.rho # Proposed solution would give 0.6860779610037084
987.4305968948237

Instead, I added a small hook in thermosteam's PhaseHandle class to use the gas phase when the temperature is above critical. I think this is consistent with how supercritical pure component properties should be the same for either 'l' or 'g' phases. Your code should also work now without having to re-simulate :).

Thanks!

[BenPortner]

@yoelcortes Reopening this because there seems to be an issue still. With my original version, using COOLPROP in the nitrogen liquefaction example (copied again below) works straight away. With your version it still fails. Not sure why.

import biosteam as bst

# setup thermodynamic backend
bst.settings.set_thermo(["N2"])
bst.settings.chemicals["N2"].Cn.l.method = "COOLPROP"
bst.settings.chemicals["N2"].reset_free_energies()
bst.settings.mixture.include_excess_energies = True # this will ensure pressure dependent enthalpies are used

# inlet
inlet = bst.Stream("inlet", N2=1, T=280, P=200e5, phase="g")

# recycle HX
regenerator = bst.units.HXprocess("regenerator", ins=(inlet, "flash_gas"), outs=("valve_in", 'out'), dT=10)
valve_in = regenerator.outs[0]

# expansion valve
valve = bst.units.IsenthalpicValve("expansion", ins=valve_in, outs="valve_out", P=1e5, vle=True)
valve_out = valve.outs[0]

# flash drum
flash = bst.units.Flash("flash", ins=valve_out, outs=("flash_gas", "flash_liquid"), Q=0, P=1e5)
gas_out = flash.outs[0]
liquid_out = flash.outs[1]

# connect flash gas to recycle HX
regenerator.ins[1] = gas_out

# run simulation
sys = bst.main_flowsheet.create_system('nitrogen_liquefaction')
sys.simulate()

Furthermore, I cannot reproduce the s.rho # 0.686... error you mention. It works fine on my machine.

tmo.settings.set_thermo(['Water', 'CO2'])
s = tmo.Stream(Water=2, CO2=0.00000001, T=320)
s.rho # Proposed solution would give 0.6860779610037084
Out[2]: 987.4305968948237

[yoelcortes]

@BenPortner, thanks for reopening. Ahh, I misread your code... I thought it used the min of Tcs (instead of max). This presents a different issue: as soon as another chemical is defined, you would get different results. For example, if you use bst.settings.set_thermo(["N2", "H2O"]) (H2O has a higher critical point), then you would go back to using the liquid heat capacity of N2 because H2O has the highest critical point. This time I tried it locally to confirm jaja (I'll keep this practice too).

I fixed the bug with the solution I provided earlier. Now hooks from phase handles will be used for ideal mixtures too. I tested it and the system converges too.

Thanks!