Wrong enthalpies after system simulation

[BenPortner]

It seems that system simulations, which include feedback streams somehow interfer with enthalpy caching. E.g. the mixer outlet enthalpy of the following system comes out wrong:

# setup thermodynamic backend
bst.settings.set_thermo(["CH4"])
bst.settings.mixture.include_excess_energies = True

# prevent "no cooling agent that can cool under ..." error
bst.HeatUtility.cooling_agents[-1].T = 100

# inlet
inlet = bst.Stream("inlet", CH4=1, units="kg/s", T=6.85+273.15, P=1e5, phase="g")

# recycle mixer
mixer = bst.Mixer("mixer", ins=(inlet, "recylce"), outs="C1_in", rigorous=True)
mixer_out = mixer.outs[0]
mixer_recycle_in = mixer.ins[1]

# compressor first stage
C1 = bst.units.IsentropicCompressor("C1", ins=mixer_out, outs="C1_out", P=5e5, eta=1)
C1_out = C1.outs[0]

# intercooling first stage
HX1 = bst.units.HXutility("HX1", ins=C1_out, outs="HX1_out", T=inlet.T)
HX1_out = HX1.outs[0]

# compressor second stage
C2 = bst.units.IsentropicCompressor("C2", ins=HX1_out, outs="C2_out", P=25e5, eta=1)
C2_out = C2.outs[0]

# intercooling second stage
HX2 = bst.units.HXutility("HX2", ins=C2_out, outs="HX2_out", T=inlet.T)
HX2_out = HX2.outs[0]

# compressor third stage
C3 = bst.units.IsentropicCompressor("C3", ins=HX2_out, outs="C3_out", P=100e5, eta=1)
C3_out = C3.outs[0]

# cooler
cooler = bst.units.HXutility("cooler", ins=C3_out, outs="cooler_out", T=-63.15+273.15)
cooler_out = cooler.outs[0]

# recycle HX
regenerator = bst.units.HXprocess("regenerator", ins=(cooler_out, "flash_gas"), outs=("throttle_in", mixer_recycle_in), dT=14.78)
throttle_in = regenerator.outs[0]

# throttling
valve = bst.units.IsenthalpicValve("expansion", ins=throttle_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 regenerator
regenerator.ins[1] = gas_out

# run simulation
sys.simulate()

print(f"Wrong mixer outlet enthalpy: {mixer.outs[0].H}") # -3.8e5
print(f"Correct mixer outlet enthalpy: {mixer.outs[0].copy().H}") # -1.4e6

[yoelcortes]

@BenPortner, thanks again for submitting this bug. I previously fixed the bug in Stream objects, but forgot about multistream objects :(. It is fixed now and has tests so it doesn't repeat.

Thanks!

[BenPortner]

@yoelcortes I'm reopening this because the initial problem is fixed, but something is still not right. When doing Stream.vle() repeatedly on valve_out, the vapor fraction keeps changing:

valve_out
"""
Out[2]: MultiStream: s.3
 phases: ('g', 'l'), T: 111.47 K, P: 100000 Pa
 flow (kmol/hr): (g) CH4  366
                 (l) CH4  227
"""
valve_out.vle(H=valve_out.H, P=valve_out.P)
valve_out
"""
Out[4]: MultiStream: s.3
 phases: ('g', 'l'), T: 111.47 K, P: 100000 Pa
 flow (kmol/hr): (g) CH4  368
                 (l) CH4  224
"""
valve_out.vle(H=valve_out.H, P=valve_out.P)
valve_out
"""
Out[6]: MultiStream: s.3
 phases: ('g', 'l'), T: 111.47 K, P: 100000 Pa
 flow (kmol/hr): (g) CH4  371
                 (l) CH4  222
"""
...

[yoelcortes]

@BenPortner

I believe this issue might be a bug in Python/Numpy. For some reason, passing a zip object returns the wrong result, but if you convert the values to a list it works as expected. I tried to recreate a simpler version of this issue so that I can submit it to numpy or python, but could not...

In any case, I just converted the zip object to a list and it fixed the problem (https://github.com/BioSTEAMDevelopmentGroup/thermosteam/commit/ef8a080f0b50c85a7bf25f80698390644c41d62a). I might post this picture in a couple of weeks to python/numpy and see if they can help me recreate the issue for them.

Thanks again for posting this bug!

[yoelcortes]

@BenPortner

How silly of me. This was not a bug in Python/Numpy. I just realized that phase_mol is iterated twice if include_excess_energies is True. But zip objects are only iterated through once (second iteration of a zip object does nothing). So if phase_mol is a zip object, excess energies won't be accounted for (and phase_mol must therefore be a list or tuple).

[BenPortner]

Good job catching this one @yoelcortes! I also had trouble narrowing this one down to its source. I'm glad that you could 😃

Btw. a nice side effect of fixing this bug is that the solver got more stable, too. The following system did not converge with the old thermosteam version. After your fixes it converges without problems:

import biosteam as bst

# setup thermodynamic backend
bst.settings.set_thermo(["N2"])
bst.settings.mixture.include_excess_energies = True
bst.settings.chemicals["N2"].reset_free_energies()

# inlet
inlet = bst.Stream("inlet", N2=2.75, units="kg/s", T=6.85+273.15, P=200e5, phase="g")

# recycle HX
regenerator = bst.units.HXprocess(
    "regenerator", ins=(inlet, "flash_gas"),
    outs=("throttle_in", 'out'), dT=9.55,
    # phase0='g', phase1='g'  this is not necessary anymore
)
throttle_in = regenerator.outs[0]

# throttling
valve = bst.units.IsenthalpicValve("expansion", ins=throttle_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 regenerator
regenerator.ins[1] = gas_out

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