Different result using substance.p() and substance.ps() in multiphase substance

[juanjpnv]

Hello! While working with multi-phase water, I observed that the saturated pressure methods substance.ps() yield different results compared to the regular pressure method substance.p() when the quality is set to zero or one.

Exemple:

water = pm.get('mp.H2O')

water.ps(T=400)       # return 2.457652635418225 bar
water.p(T=400, x=0) # return 2.45743614 bar

The discrepancy is quite small, appearing in the 4th decimal place. However, other thermodynamic properties do not exhibit similar divergence.

Regarding the pressure case, is this variation expected? Could it be related to floating-point arithmetic or Pyromat's internal calculations?

[chmarti1]

That's a really good question, and thanks for your patience - it took me quite a while to get you an answer. Yes, this variation is expected. Here's why...

The model used to calculate the saturation points is actually a different equation than the one used to calculate all other points. See the original source for all the details, but that's quite a bit of reading. Here's the short version:

There is a global model that is used to calculate all data points that are not saturated. It is formulated as A(T,d) (helmholtz free energy as a function of temperature and density). It is possible to calculate the equilibrium saturation points, but that is a computationally expensive operation. Instead, the original sources also supply polynomial models for the saturation line. See the PYroMat handbook Section 7.2.7 page 112.

These two models agree with one another very closely, but as you point out, it is not numerically possible for them to be exactly the same. This is tolerated because they agree to within better than the reported model precision.

Does this cause a problem? Not usually, but it is something we have to be aware of sometimes. For example, I need to be aware that the saturation line is not precise when I am writing numerical iteration routines. I sometimes need to allow the algorithm to wander just a little to either side of the saturation line even though I know that isn't strictly where the solution will lie.

Thanks for your question!

[chmarti1]

Now that Gibbs energy is a property supported by the multi-phase back end, it is practical to implement the Maxwell criteria directly to polish the saturation properties to be in more precise agreement with the equation of state. The Maxwell criteria are that the Gibbs energy and pressure for the vapor and liquid phases will be equal in equilibrium.

$$ g(T, \rho_L) = g(T, \rho_V) $$

$$ p(T, \rho_L) = p(T, \rho_V) $$

Especially if the initial guess is very close to the correct solution, root polishing using Newton iteration becomes a practical possibility.

In the next release of PYroMat, this feature is planned in the mp1 class, and work is underway to implement a mp2 class, which does not require the original source to include saturation ancillary equations.