Made it possible to use MoltenSalt_ph as incompressible medium

[casella]

This fix makes it possible to also use MoltenSalt_ph with singleState=true, i.e. as incompressible medium. It is 100% backwards compatible, since it doesn't change the default setting of singleState, which is still false.

[jdpipe]

Thanks Francesco for providing this input! Obviously the change you propose could alter the results obtained by our existing models, and so I wanted to ask whether you can provide some reference to a preferred theoretical basis for modelling of incompressible fluids which is consistent with the available for molten salt and related media, eg c_p(T) and rho(T). Our conclusion, possibly incorrect, was that, if c_p(T) and rho(T) are both arbitrary functions independent of p, then in general u = h - p v = h - p / rho must hold at all points, with h = integ(c_p(T), T) in general. But is there some thermodynamic relation that allows you to replace p with p_ref as you indicated, without loss of correctness?

[casella]

@jdpipe, first of all this PR won't change your results at all, as long as you don't change singleState when you use the medium model, and keep it to its default of false, as it is today.

The rationale of singleState is as follows.

In general you have a dynamic mass balance and a dynamic energy balance equation for each volume. If you consider the fluid compressibility, i.e. you assume drho/dp <> 0, and you make a change of state variable to (p,T) or (p,h), you end up with two dynamic equations in the two states p,T or p,h. If the fluid is a liquid, the resulting system is usually (very) stiff, because the hydraulic RC time constant turns out to be many orders of magnitude smaller than the residence time, which roughly determines the thermal time constant.

If you are mainly interested in the thermal phenomena and you want to get rid of the fast time constant, you have to end up with a single state variable, i.e. either T or h. Unfortunately, if you set drho/dp = 0 but still consider the fact that both T = T(p,h) and rho = rho(p,h), when you make the automatic state variable change from der(M) and der(E) to der(p) and der(h), you still end up with a second order system, that may still be stiff.

However, incompressible fluids such as water or molten salts are often handled at relatively moderate pressures, e.g. a few bars in building HVAC systems, and have relatively high density, usually around 1000 kg/m3. In that case, considering the effect of pressure onto enthalpy is theoretically correct, but in practice it amounts to hair splitting. If you have a fluid with 1000 kg/m3 density, a change of pressure of 1 bar, i.e. 1e5 Pa, gives a contribution of 100 J/kg, which corresponds to 0.025 degC of temperature change in the case of water.

Hence, in this case one can consider rho = rho(h), cp = cp(h), and e = h - pref/rho(h), so that everything is just a function of a single thermal state variable. Ditto for the case of T, which is not covered by this PR.

If you start from models written as der(M) = ... and der(E) = ...., and set singleState=true, using these state functions and stateSelect will end up having a first-order system, that takes into account the thermal energy storage in the fluid and the thermal expansion effects with more than reasonable accuracy, avoiding the burden of the extra pressure state.

Of course it is not mandatory to do so, but that's the point of my PR, to make it possible optionally. In fact, you could then run two different simulations, one with singleState = true and the other one with singleState = false, and check that the results in terms of temperatures and power flows are practically the same.