ReferenceMoistAir, substances mass conservation issue

[GaelEnee]

Hello,

I would like to report an issue with the Modelica.Media.Air.ReferenceMoistAir media.

From unexpected simulation results, I implemented the following simple model (cf. enclosed model): a constant mass flow source (Modelica.Fluid.Sources.MassFlowSource_T) with X = {0, 1} connected to a constant volume tank (Modelica.Fluid.Vessels.ClosedVolume) with X = {0.5, 0.5} @ t = 0. The tank water mass is expected to remain the same and the total mass was expected to increase linearly. Instead, the tank water mass increased and the total fluid mass increasing beyond expectation and beyond the water mass increase.

XBalance.zip

The issue has been observed with Dymola 2020 using Modelica Standard Library - Version 3.2.3 and confirmed using Modelica Standard Library - Version 4.0.0. As a comparison, using Modelica.Media.IdealGases.MixtureGases.AirSteam, the observed error is more likely due to computing small error propagation.

Quick investigations pointed to ReferenceMoistAir density and derivated functions.

Please let me know if problem if the problem has already been reported and if solution has been developed.

Thank you !

[GaelEnee]

Hello,

After investigations, the problem seem to be solvable implementing the followings:

1/ Correction of the derivative within Utilities.pd_pTX_der():

    pd_der := (xw_der*(Modelica.Media.Air.ReferenceMoistAir.k_mair + xw) -
      xw*xw_der)*p + xw/(Modelica.Media.Air.ReferenceMoistAir.k_mair + xw)*
      p_der;

becomes

    pd_der := xw_der*Modelica.Media.Air.ReferenceMoistAir.k_mair/
      (Modelica.Media.Air.ReferenceMoistAir.k_mair + xw)^2*p
      + xw/(Modelica.Media.Air.ReferenceMoistAir.k_mair + xw)*p_der;

2/ Reformulating of Utilities.rho_pTX() and Utilities.rho_pTX_der() using X[1] and Xsat rather than xw and xws for numerical considerations (might hide errors in the current code). I can provide an example upon request.

Globally, the whole media code might be subject to improvements and simplifications.

[wischhusen]

I will check as soon as possible and revert back. Many thanks for the work.

[wischhusen]

Sorry, I made two changes for issues in my master. Is someone else able to put

    pd_der := xw_der*Modelica.Media.Air.ReferenceMoistAir.k_mair/
      (Modelica.Media.Air.ReferenceMoistAir.k_mair + xw)^2*p
      + xw/(Modelica.Media.Air.ReferenceMoistAir.k_mair + xw)*p_der;

in function pd_pTX_der since the other PR is pending and I do not won't to roll back changes? Many thanks.

[beutlich]

Is someone else able

Done. See #3834.

[wischhusen]

Many thanks!

[thorade]

Here is a test model for this function:

model pd_der_test
  replaceable package Medium = Modelica.Media.Air.ReferenceMoistAir;

  Real der_T = 300;
  Real T = 73.15 + der_T * time;
  Real der_p = -1e5;
  Real p = 1e7 + der_p * time;
  Real[Medium.nS] X = Medium.reference_X;

  Real var_from_function = Medium.Utilities.pd_pTX(p, T, X);
  Real timederiv = Medium.Utilities.pd_pTX_der(p, T, X, der_p, der_T, {0,0});
  Real var_from_timederiv;

  Real error_abs =  var_from_function - var_from_timederiv;
  Real error_rel =  abs(error_abs / var_from_function);

initial equation 
  var_from_timederiv = var_from_function;

equation 
  // var_from_timederiv = integral(timederiv);
  // but integral() operator does not exist
  // formulate it the other way round  
  der(var_from_timederiv) = timederiv;

  annotation (experiment(StopTime=100));
end pd_der_test;

[beutlich]

Here is a test model for this function:

Thanks. If you like you can add it to the corresponding PR #3834.

[beutlich]

Here is a test model for this function:

Edit: I do not see how #3834 fixes the two error variables of test model pd_der_test.

[hubertus65]

To me, the proposed derivative looks wrong. If the function is:
//pd := xw/(Modelica.Media.Air.ReferenceMoistAir.k_mair + xw)*p;
And for better readability, I propose a shorthand and minor re-write:
Real k = Modelica.Media.Air.ReferenceMoistAir.k_mair; //pd := xw/(k + xw)*p;
I get:
pd_der = (p*xw_der + xw*p_der)/(k + xw) - p*xw*xw_der/(k + xw)^2
Could somebody please verify this, also with the test function from Matthis?

[beutlich]

Could somebody please verify this, also with the test function from Matthis?

That is also correct, i.e. for pd = xw/(k + xw)*p the derivative is pd_der = (p*xw_der + xw*p_der)/(k + xw) - p*xw*xw_der/(k + xw)^2 or simplified pd_der = xw_der*k/(k + xw)^2*p + xw/(k + xw)*p_der as proposed by @GaelEnee or @wischhusen in #3666 (comment).

Still, my previous concern holds:

Edit: I do not see how #3834 fixes the two error variables of test model pd_der_test.

[HansOlsson]

I believe there's a better solution will comment in #3834