eConstThermo Eref, Tref, and Hc with reactinTwoPhaseEulerFoam
 

Hi,

I've been playing around with the wallBoiling tutorial case for reactingTwoPhaseEulerFoam located in $FOAM_TUTORIALS/multiphase/reactingTwoPhaseEulerFoam/RAS/wallBoiling/
Specifically, the thermophysicalProperties.(gas | liquid) - there are a few properties which are relevant to this thread:

• Eref [J/kg] - the eConstThermo.H file says "Reference sensible enthalpy around which to linearise", which is not what appears to be the implementation in eConstThermoI.H, where it is treated as a reference sensible internal energy.
• Tref [K] - the reference temperature used for Eref.
• Hf [J/kg] - the enthalpy of formation of the compound at Pstd, Tstd.
• Cv [J/kg K] - the heat capacity at constant volume.

Looking at the code of eConstThermoI.H, we see that Es is defined as you might expect:

So (ignoring the internal energy separation term from EquationOfState) the sensible internal energy is just a linear extrapolation from the point (Tref_, Esref_) with gradient Cv_, which are all constructed from the thermophysicalProperties file. All good here.

Now, my points of confusion:

1. Why does changing Tref.gas and Eref.gas to different (but equivalent) values produce different results? E.g. for the gas in the wallBoiling case, Cv.gas = 12078.4 J/kg K. But (Tref, Eref) = (373.55, 2675500) (image centre) and (273.55, 1467660) (image right) give different results. This also occurs when using hConstThermo.
2. In what universe does water vapour have a specific heat capacity of 12kJ/kg K at 100 deg C?!
3. What even is Eref? And for that matter, Href in hConstThermo? My understanding is that there is only one way to define sensible enthalpy, as , where ? But Eref causes a non-zero Es at Tstd. My only explanation is that the phase change model uses Es or Hs only, and doesn't consider Hf. So Eref / Href is used a workaround in order to get a correct enthalpy difference between liquid and gas phases in the phase change code. From what I can tell, this appears to be the case (specifically looking at ThermalPhaseChangeSystem.C). Can anyone confirm/deny this?
4. Why is Hf completely irrelevant to the results of the wallBoiling tutorial case? Setting Hf.gas=1e6 (image left) has no impact on the results. Seems to confirm the theory of Hf being completely ignored.
5. Why is Hf used in the laminar/bubbleColumnEvaporating tutorial case? I assume it's to do with the gas phase being a multi-component mixture of H2O and air, but they have no in-phase reactions? Is Hf also ignored there?

Point number 1 is most confusing to me - I see no-where else in the code where Eref and Tref are used (except in the Gibb's Free Energy at standard pressure, but that function shouldn't be affected by the example I gave, and regardless, this isn't used anywhere for this case). And yet, the results change unexpectedly.

I have had a look at a bunch of different theory books (Fluent theory guide, Versteeg, OpenFoam "Little User Guide"), but I haven't yet found an adequate explanation - so I have given up resolving these questions on my own.

I would really appreciate any help with this!

#1: After going through the solver line-by-line with a debugger, and comparing between the two cases, I found that the most likely culprit is a tiny floating-point rounding error which accumulates over time.

The only difference I found between the original Tref & Eref combination and the modified (but equivalent) version is a ~2e-26 difference between the off-diagonal matrix coefficients, which I tracked down to floating-point rounding. This is possibly due to the definition of sensible internal energy in eConstThermoI.H amplifying these rounding errors.

It seems I just happened to stumble across a tutorial case that was particularly sensitive to these tiny errors accumulating over time.

#2: Still not sure. Maybe a typo that was never corrected?

#3 and #4: After looking through the evaporation code again, I haven't seen any reference to Hf or Hc - so Hf is irrelevant, and therefore a non-zero sensible enthalpy must be defined at Tstd in order to capture the latent heat of vaporisation.

Any hints on #2 or #5?