Pressure-dependency of enthalpy is not accounted for

[BenPortner]

I am currently working on an isothermal compressor unit for biosteam. I realized that the power demand is not calculated correctly at high pressures. This is because the current implementation of Mixture.H does not account for pressure. This is okay for ideal gases, but not for real gases. E.g. compressing hydrogen from 20 bar to 350 bar:

current behavior

import thermosteam as tmo
tmo.settings.set_thermo(["H2"])
T = 298.15
feed = tmo.Stream(H2=1, T=T, P=20e5, phase='g')
S0 = feed.S
H0 = feed.H
feed.P = 350e5
feed.T = T
S = feed.S
H = feed.H

print(f"T*dS: {T*(S-S0):.2f}, dH: {H-H0:.2f}") # T*dS: -7095.27, dH: 0.00
power_tmo = H-H0 - T*(S-S0) # 7095 kJ/hr

expected behavior

import CoolProp.CoolProp as cp

feed = cp.AbstractState("HEOS", 'hydrogen')
feed.update(cp.PT_INPUTS, 20e5, T)
S0 = feed.smolar()
H0 = feed.hmolar()
feed.update(cp.PT_INPUTS, 350e5, T)
S = feed.smolar()
H = feed.hmolar()

print(f"T*dS: {T*(S-S0):.2f}, dH: {H-H0:.2f}") # T*dS: -7214.74, dH: 381.30
power_cp = H-H0 - T*(S-S0) # 7596 kJ/hr
pass

Note that the entropy solution is also slightly off (~100 kJ/mol) but the difference is much smaller than the error introduced by dH (~400 kJ/mol).

[yoelcortes]

Hi Ben,

To include pure component excess energies, set thermo.mixture.include_excess_energies = True. Unfortunately, it looks like thermo fails to calculate excess enthalpy/entropy at P=350e5, T=298.15.

import thermosteam as tmo
thermo = tmo.Thermo(['H2'])
thermo.mixture.include_excess_energies = True
tmo.settings.set_thermo(thermo)
T = 298.15
feed = tmo.Stream(H2=1, T=T, P=20e5, phase='g')
S0 = feed.S # This works
H0 = feed.H # This works
feed.P = 350e5
feed.T = T
S = feed.S # This fails 
H = feed.H # This also fails

I think Caleb might be able to help you out further on this. You can post an issue in thermo that the following fails:

from thermo.eos import IG, PR 
PR(T=298.15, P=350e5, Tc=33.2, Pc=1296960.0, omega=-0.22).S_dep_g # AttributeError: 'PR' object has no attribute 'S_dep_g'

If you'd like to use coolprop for all thermodyamic properties, one potential solution is to create your own Mixture object that interfaces with coolprop (instead of the default one created by thermosteam):

https://biosteam.readthedocs.io/en/latest/API/thermosteam/mixture/Mixture.html

Hope this helps, Thanks!

[yoelcortes]

Hi Ben,

I saw the issue on CalebBell/thermo#118, thanks for posting! Here are the new results with the fallback when phase identification suggests another phase (which doesn't actually matter since it is supercritical).

import thermosteam as tmo
import thermo as tm
thermo = tmo.Thermo([tmo.Chemical('H2', eos=tm.SRK)])
thermo.mixture.include_excess_energies = True
tmo.settings.set_thermo(thermo)
T = 298.15
feed = tmo.Stream(H2=1, T=T, P=20e5, phase='g')
S0 = feed.S # This works
H0 = feed.H # This works
feed.P = 350e5
feed.T = T
S = feed.S 
H = feed.H 
power_cp = H-H0 - T*(S-S0) # 7596 kJ/hr in coolprop
print(f"T*dS: {T*(S-S0):.2f}, dH: {H-H0:.2f}, power_cp: {power_cp:.2f}") # T*dS: -7214.74, dH: 381.30 in coolprop

Output

T*dS: -7261.40, dH: 307.44, power_cp: 7568.85

Note that I'm using SRK, which is the closest we can get to coolprop's HEOS. I'm not sure which one is more accurate (you may need to check literature), but both seem quite close. EDIT: The calculation of S and H also depends on the Cp model used, so that might also play a role if you see any deviations from literature.

Hope this helps!

[BenPortner]

Hi @yoelcortes,

I think they are close enough. Thank you very much for your help! 😃
I will now close this issue.