Consider phases in reaction enthalpy calculation

[BenPortner]

Description

The documentation clearly states that Reaction.dH does not currently account for the heat of vaporization. However, this leads to faulty reaction enthalpies in cases where e.g. water is not in liquid state at the end of a reaction. Fixing this should be easy. There is already a function thermosteam.chemicals.reaction.Hf, which takes CASRN and phase as attributes and returns the correct heat of formation. Using this function instead of the STP-heats of formation in thermosteam.reaction._reaction.Reaction.dH should solve the issue.

Minimum error example

import thermosteam as tmo

tmo.settings.set_thermo(["CH4", "O2", "CO2", "H2O"])
combustion = tmo.Reaction(
    'CH4,g + O2,g -> CO2,g + H2O,g',
    correct_atomic_balance=True,
    reactant='CH4',
    X=1,
)
print(combustion.dH) # expected: -800e3, actual: -890e3

[yoelcortes]

Hi Ben,

Thanks for the issue! Heats of vaporization are temperature dependent and I consider temperature dependence out of the scope of a reaction object. We can compute heat of reactions with phase changes using Stream objects together with Reaction objects:

import thermosteam as tmo
tmo.settings.set_thermo(['Water', 'Methane', 'CO2', 'O2'])
combustion = tmo.Reaction('Methane + O2 -> Water + CO2',
                          reactant='Methane', X=1,
                          correct_atomic_balance=True)
feed = tmo.Stream(Methane=1, O2=2, T=298.15, phase='g')
H0 = feed.Hnet
combustion(feed)
Hf = feed.Hnet
print(f'Lower heating value at 298.15 K reference state: {Hf - H0:.2f} kJ/hr')

feed = tmo.Stream(Methane=1, O2=2, T=400, phase='g')
H0 = feed.Hnet
combustion(feed)
Hf = feed.Hnet
print(f'Lower heating value at 400 K reference state: {Hf - H0:.2f} kJ/hr')

Output

Lower heating value at 298.15 K reference state: -802735.67 kJ/hr
Lower heating value at 400 K reference state: -801756.16 kJ/hr

I think this issue is related to #29. You can check it out for additional context. At the time the code had many bugs, but the same ideas apply.

Thanks!

[BenPortner]

Hi @yoelcortes,

thanks for the link to the related issue. It seems I didn't search carefully, sorry! I hope it is fine if I go on in this thread nevertheless:

I agree that the temperature dependence of the heat of vaporization (dH_vap) is out of the scope of the reaction object. However, I do believe it is in the scope of the reaction object to calculate the correct heat of reaction - and that is currently not the case. Moreover, I believe that it is not necessary to model the temperature dependence of dH_vap to get the reaction enthalpy right. The enthalpies of formation (Hf) are defined at STP. Hence, the difference between Hf_gas and Hf_liquid is the heat of vaporization at STP. Hf is thus needed at STP, not at any other temperature. Am I wrong?

[yoelcortes]

Yeah, adding latent heats at 298.15 K makes sense. I think it all checks out, so I'll go ahead and add it for reactions where phases are specified. I might not get to it today, but surely by this weekend.

Thanks!

[yoelcortes]

Done! I added tests so I'm pretty confident of results, but I'll let you close this issue so you can get a chance to test it too.

Thanks!

[BenPortner]

Amazing work @yoelcortes! Thanks for the swift implementation.

I tested it for the methane combustion example and it worked perfectly. However, there seems to be a problem when involving liquid hydrogen:

import thermosteam as tmo

tmo.settings.set_thermo(["H2O", "H2", "O2"])
electro_liquef = tmo.Reaction(
    'H2O,l -> H2,l + O2,g',
    correct_atomic_balance=True,
    reactant='H2O',
    X=1,
)
print(electro_liquef.dH / tmo.Chemical("H2O").MW)

"""
RuntimeError: Failed to extrapolate enthalpy of vaporization method 'DIPPR_PERRY_8E' at T=298.15 K for component with CASRN '1333-74-0 (H2)'
"""

More generally, it seems the heat of vaporization cannot be calculated for hydrogen:

tmo.Chemical("H2").Hvap(298) # same error

Not sure if I should open a new issue for this?

Thanks for the amazing work! Ben

[yoelcortes]

Thanks, Ben! Hydrogen cannot be liquefied above it's critical temperature, no matter how much pressure is applied (the gas and liquid phases are the same phase). So I think it's best just to use 'H2O,l -> H2,g + O2,g' instead since Hvap does not exist.

Thanks!

[BenPortner]

Silly me, hydrogen cannot be liquified above 33K, so there is no formation enthalpy of liquid H2 at 298K. Makes sense! 😅

It works fine with gaseous H2, as you suggest. Looks like everything is good! Closing the issue now 😄

Thanks. Ben