Media.Water.IF97_Utilities and Media.Air.ReferenceAir: h_pT_der / h_dT_der have the wrong sign in the pressure/density term #4823

Describe the bug

The derivative functions attached to h_pT and h_dT return the negative of the true pressure/density derivative of the enthalpy. The temperature terms are correct.

For any equation of state, h = u + p/d and (∂u/∂d)_T = (p − T·(∂p/∂T)_d)/d², so

(∂h/∂d)_T = (d·p_d − T·p_T)/d² and (∂h/∂p)_T = (∂h/∂d)_T / p_d,

with p_d = aux.pd and p_T = aux.pt. The code has the numerator with the opposite sign:

• IF97_Utilities.h_pT_der, region 3 (L7584-L7587): (-aux.rho*aux.pd + T*aux.pt)/(aux.rho*aux.rho*aux.pd)
• IF97_Utilities.h_dT_der (L8028-L8037), every branch:
  • region 3: (-d*aux.pd + T*aux.pt)/(d*d)
  • region 4: T*aux.dpT/(d*d), while Clapeyron gives (∂h/∂d)_T = −T·(dp/dT)_sat/d² in the two-phase region
  • regions 1, 2, 5: -(-1/d + T*aux.vt)/(d*d*aux.vp), while (∂h/∂p)_T = 1/d − T·v_T and (∂d/∂p)_T = −d²·v_p give (∂h/∂d)_T = −(1/d − T·v_T)/(d²·v_p)
• ReferenceAir.Air_Utilities.h_pT_der (L1549-L1561) and h_dT_der (L1801-L1812): the same negated numerator

hvl_p_der in the same file uses the correct (bpro.d*bpro.pd - bpro.T*bpro.pt) (L1017), and h_pT_der regions 1, 2, 5 (1/aux.rho - aux.T*aux.vt) are also correct.

These functions are live derivative annotations (derivative(noDerivative=aux) = h_pT_der / h_dT_der at IF97_Utilities.mo L7552, L7994, ReferenceAir.mo L1532, L1784). A tool uses them whenever it differentiates h with p,T or d,T states, for example WaterIF97_pT in region 3 (above 16.5 MPa, 623 to 863 K) and every ReferenceAir.Air_pT / Air_dT state. For an ideal gas (∂h/∂d)_T = 0, so ideal-gas tests cannot catch this.

In a closed rigid tank heated with 10 kJ, WaterIF97_pT at 25 MPa / 650 K loses 8.75 kJ from the energy balance, and ReferenceAir.Air_pT at 100 bar / 200 K loses 515 J of 1 kJ. The same tanks with p,h states, which do not call these functions, balance to the solver tolerance.

Expected behavior

h_dT_der and h_pT_der equal the total derivatives of h_dT and h_pT, so err in the model below stays at solver tolerance for all four media.

Proposed fix (density/pressure terms only; the temperature terms stay as they are):

 // IF97_Utilities.h_pT_der, region 3, and ReferenceAir.Air_Utilities.h_pT_der
-h_der := ((-aux.rho*aux.pd + T*aux.pt)/(aux.rho*aux.rho*aux.pd))*p_der + ...
+h_der := ((aux.rho*aux.pd - T*aux.pt)/(aux.rho*aux.rho*aux.pd))*p_der + ...

 // IF97_Utilities.h_dT_der, region 3, and ReferenceAir.Air_Utilities.h_dT_der
-h_der := ((-d*aux.pd + T*aux.pt)/(d*d))*d_der + ...
+h_der := ((d*aux.pd - T*aux.pt)/(d*d))*d_der + ...

 // IF97_Utilities.h_dT_der, region 4
-h_der := T*aux.dpT/(d*d)*d_der + ...
+h_der := -T*aux.dpT/(d*d)*d_der + ...

 // IF97_Utilities.h_dT_der, regions 1, 2, 5
-h_der := (-(-1/d + T*aux.vt)/(d*d*aux.vp))*d_der + ...
+h_der := ((-1/d + T*aux.vt)/(d*d*aux.vp))*d_der + ...

Minimal reproducible example

A closed, rigid tank with heat input Q has constant mass and dU/dt = Q, so err = U - U0 - Q*time must stay 0. U is computed from the medium state.

partial model TankBase
  replaceable package Medium = Modelica.Media.Interfaces.PartialMedium;
  parameter Modelica.Units.SI.Volume V = 1e-3;
  parameter Modelica.Units.SI.HeatFlowRate Q = 1000;
  parameter Modelica.Units.SI.Pressure p0;
  parameter Modelica.Units.SI.Temperature T0;
  Medium.BaseProperties medium(preferredMediumStates = true);
  Modelica.Units.SI.Mass m = medium.d*V;
  Modelica.Units.SI.Energy U = m*medium.u;
  Modelica.Units.SI.Energy E_in(start = 0, fixed = true) "Heat added";
  Real err(unit = "J") = U - U0 - E_in "Energy balance defect";
  parameter Modelica.Units.SI.Energy U0(fixed = false);
initial equation
  medium.p = p0;
  medium.T = T0;
  U0 = U;
equation
  der(E_in) = Q;
  der(m) = 0;
  der(U) = Q;
end TankBase;

model Tank_IF97_pT "IF97, p,T states, region 3"
  extends TankBase(redeclare package Medium = Modelica.Media.Water.WaterIF97_pT, p0 = 25e6, T0 = 650);
end Tank_IF97_pT;

model Tank_IF97_ph "Same tank with p,h states"
  extends TankBase(redeclare package Medium = Modelica.Media.Water.WaterIF97_ph, p0 = 25e6, T0 = 650);
end Tank_IF97_ph;

model Tank_RefAir_pT "ReferenceAir, p,T states"
  extends TankBase(redeclare package Medium = Modelica.Media.Air.ReferenceAir.Air_pT, p0 = 100e5, T0 = 200, Q = 100);
end Tank_RefAir_pT;

model Tank_RefAir_ph "Same tank with p,h states"
  extends TankBase(redeclare package Medium = Modelica.Media.Air.ReferenceAir.Air_ph, p0 = 100e5, T0 = 200, Q = 100);
end Tank_RefAir_ph;

Simulated with stopTime = 10, tolerance = 1e-8:

Tank_IF97_pT:    medium.T = 650.7867431   err = -8754.468567
Tank_IF97_ph:    medium.T = 656.3902494   err = -0.007367730259
Tank_RefAir_pT:  medium.T = 202.8687364   err = -515.0483899
Tank_RefAir_ph:  medium.T = 205.9248615   err = 0.0001940232177

Environment

OpenModelica 1.25.5, MSL master at 4c40388. Not tested in Dymola.

Additional context

We found this by checking each derivative annotation against the derivative of its function with a theorem prover.