Modelica.Media.R134a problem #3163

beutlich · 2019-10-21T11:09:13Z

Reported by CTG on 15 Oct 2019 07:56 UTC

Although it is not an OpenModelica related problem, I would like to share the following experience.
While checking my own media libraries compatible with Modelica.Media, I have found what I think is a bug in the Modelica.Media.R134a package. The calculation fail is shown in the attached file. Coming from a density between bubble and dew densities, it fails in the calculation of the vapor quality. The origin of the problem seems a bad implementation of the setState_dTX function: the calculation of Helmholtz derivatives and specific enthalpy is done before to check the phase fractions, as if the given density would be a one phase density. This produces a complete fail for calculations in the two phase region coming from density and temperature.
I do not know if this is a known bug, or I am wrong in some point.

ModelicaTests.mo

Originally reported at https://openmodelica.org/forum/default-topic/2801-modelica-media-r134a-problem#p9417

HansOlsson · 2019-10-21T11:39:33Z

The bug-report lacks two parts:

The actual error messages.
The desired result. In this case the goal seems to be to get a specific density and temperature, and I don't know if it is achievable.

beutlich · 2019-10-21T11:54:49Z

The bug-report lacks two parts

True. Hopefully @wischhusen has enough information for an analysis.

HansOlsson · 2019-10-21T12:37:13Z

As for the desired result I believe that it can be computed as follows:

model ModelicaR134ATest
  package Medium = Modelica.Media.R134a.R134a_ph;
  parameter Real T=20 + 273.15;
  Medium.SaturationProperties sat;
  Medium.ThermodynamicState StateBub;
  Medium.ThermodynamicState StateDew;
  Medium.ThermodynamicState StateBi
    "a biphasic state with density halfway between bubble and dew";
  Medium.Density dBi "density half way between bubble and dew";
  Real gasFraction;
  Real h(start=3e5);
equation 
  sat =  Medium.setSat_T(T);
  StateBub =  Medium.setBubbleState(sat);
  StateDew =  Medium.setDewState(sat);
  dBi =  (StateBub.d + StateDew.d)/2;
    // system of equations:
    StateBi =  Medium.setState_phX(StateBub.p, h);
    StateBi.d=dBi;
    // end system
  gasFraction =  Medium.vapourQuality(StateBi);
end ModelicaR134ATest;

beutlich · 2019-12-04T20:08:04Z

@casella Maybe you can help out here. Thanks.

casella · 2019-12-04T22:44:54Z

@beutlich this model was implemented by @wischhusen, if I am not mistaken, I hope he or someone from his company can have a look and decide whether this is a quick fix or not. I don't think I can manage to do this unless I take three days off to study the problem thoroughly, which is unlikely to happen for a while.

wischhusen · 2019-12-17T10:20:17Z

I resolved the problem like this:

redeclare function extends setState_dTX
  "Set state for density and temperature (X not used since single substance)"
protected 
  Modelica.Media.Common.HelmholtzDerivs f "Helmholtz derivatives";
  Modelica.SIunits.SpecificHeatCapacity R "Specific gas constant";
  SaturationProperties sat "Saturation temperature and pressure";
  Modelica.SIunits.Density dl "Liquid density";
  Modelica.SIunits.Density dv "Vapor density";

algorithm 
   R := R134aData.R;
   sat := setSat_T(T);
   dl := bubbleDensity(sat);
   dv := dewDensity(sat);
   if d < dl and d > dv and T<R134aData.data.FTCRIT then
    f := Common.HelmholtzDerivs(d=d, T=T, R=R, delta=0, tau=0, f=0, fdelta=0,
          fdeltadelta=0, ftau=0, ftautau=0, fdeltatau=0);
   state.p := saturationPressure(T);
   state.h := 1/(dl/dv + 1)*(dewEnthalpy(sat) - bubbleEnthalpy(sat))
               + bubbleEnthalpy(sat);
   state.phase := 2;
   else
   f := f_R134a(d=d, T=T);
   state.p := d*R*T*f.delta*f.fdelta;
   state.h := R*T*(f.tau*f.ftau + f.delta*f.fdelta);
   state.phase := 1;
   end if;
   state.T := T;
   state.d := d;

   annotation (Documentation(revisions="<html>
<p>2019-12-17  Stefan Wischhusen: Two-phase calculation corrected.</p>
<p>2012-08-01  Stefan Wischhusen: Corrected passing-error of inputs.</p>
</html>",
      info="<html>
<p>Although the medium package is explicit for pressure and specific enthalpy, this function may be used in order to calculate the thermodynamic state record used as input by many functions. It will calculate the missing states:</p>
<ul>
<li>pressure</li>
<li>specific enthalpy</li>
</ul>
<p>
Example:
</p>
<pre>
     parameter Medium.Density d = 4;
     parameter Medium.Temperature T = 298;

     Medium.SpecficEntropy s;

     <strong>equation</strong>

     s = Medium.specificEntropy(setState_dTX(d, T, fill(0, Medium.nX)));
</pre>

</html>"));
end setState_dTX;

@beutlich could you be so kind to create a pull request before I mess the whole thing again with my outdated fork?

HansOlsson · 2019-12-17T11:22:11Z

When I managed to mess up my fork too much I found
http://stackoverflow.com/questions/42332769/how-do-i-reset-the-git-master-branch-to-the-upstream-branch-in-a-forked-reposito useful, i.e.:

git checkout master
git reset upstream/master
git pull --rebase upstream master
git push origin master --force

Obviously don't do it if your fork contains anything useful.

beutlich · 2019-12-18T20:23:46Z

@beutlich could you be so kind to create a pull request before I mess the whole thing again with my outdated fork?

Done by #3275. Please note that R was renamed to R_s, thus your proposal was suited for maint/3.2.3 only.

wischhusen · 2019-12-20T14:22:59Z

Implemented last line. Result is somehow unchanged for gas fraction but never mind.

Thanks @HansOlsson for the tipp on GIT but the reset process got stuck again. I will try again next year.
Here is the solution by @casella and me:

redeclare function extends setState_dTX
  "Set state for density and temperature (X not used since single substance)"

protected 
      Modelica.Media.Common.HelmholtzDerivs f "Helmholtz derivatives";
      Modelica.SIunits.SpecificHeatCapacity R_s "Specific gas constant";
      SaturationProperties sat "Saturation temperature and pressure";
      Modelica.SIunits.Density dl "Liquid density";
      Modelica.SIunits.Density dv "Vapor density";

algorithm 
      R_s := R134aData.data.R_s;
      sat := setSat_T(T);
      dl := bubbleDensity(sat);
      dv := dewDensity(sat);
      if d < dl and d > dv and T < R134aData.data.FTCRIT then
        f := Modelica.Media.Common.HelmholtzDerivs(
          d=d, T=T, R_s=R_s, delta=0, tau=0, f=0, fdelta=0,
          fdeltadelta=0, ftau=0, ftautau=0, fdeltatau=0);
        state.p := saturationPressure(T);
        state.h := (dl*dv/d*(dewEnthalpy(sat) - bubbleEnthalpy(sat))
                  - dv*dewEnthalpy(sat) + dl*bubbleEnthalpy(sat))/(dl - dv);
        state.phase := 2;
      else
        f := f_R134a(d=d, T=T);
        state.p := d*R_s*T*f.delta*f.fdelta;
        state.h := R_s*T*(f.tau*f.ftau + f.delta*f.fdelta);
        state.phase := 1;
      end if;
      state.T := T;
      state.d := d;

  annotation (Documentation(revisions="<html>
<p>2012-08-01  Stefan Wischhusen: Corrected passing-error of inputs.</p>
</html>",
      info="<html>
<p>Although the medium package is explicit for pressure and specific enthalpy, this function may be used in order to calculate the thermodynamic state record used as input by many functions. It will calculate the missing states:</p>
<ul>
<li>pressure</li>
<li>specific enthalpy</li>
</ul>
<p>
Example:
</p>
<pre>
     parameter Medium.Density d = 4;
     parameter Medium.Temperature T = 298;

     Medium.SpecficEntropy s;

     <strong>equation</strong>

     s = Medium.specificEntropy(setState_dTX(d, T, fill(0, Medium.nX)));
</pre>

</html>"));
end setState_dTX;

casella · 2019-12-20T15:14:48Z

@wischhusen, I would suggest that we add some test to cross check that results with setState_ph and setState_dT are consistent. I can do that if you don't have time

wischhusen · 2019-12-20T15:35:09Z

@casella many thanks! I am leaving for christmas vacation. I could do it from Jan. 6th. If you have time I would really appreciate!

casella · 2019-12-20T15:46:09Z

I'll see what I can do.

* refs #3163: Fix two-phase calculation of setState_dTX * refs #3163: Fix specific enthalpy in two-phase * Fix typo * refs #3163: Use calculated saturation pressure Co-authored-by: Stefan Wischhusen <wischhusen@xrg-simulation.de>

beutlich added the L: Media Issue addresses Modelica.Media label Oct 21, 2019

beutlich added this to the MSL4.0.0 milestone Oct 21, 2019

beutlich assigned wischhusen Oct 21, 2019

beutlich added the question Unexplained or undecided issue label Oct 21, 2019

modelica deleted a comment from wischhusen Oct 23, 2019

beutlich assigned casella Dec 4, 2019

beutlich added bug Critical/severe issue and removed question Unexplained or undecided issue labels Dec 18, 2019

beutlich pushed a commit to beutlich/ModelicaStandardLibrary that referenced this issue Dec 18, 2019

refs modelica#3163: Fix two-phase calculation of setState_dTX

a35e5f0

beutlich pushed a commit to beutlich/ModelicaStandardLibrary that referenced this issue Dec 18, 2019

refs modelica#3163: Fix two-phase calculation of setState_dTX

100cffc

beutlich mentioned this issue Dec 18, 2019

Fix two-phase calculation of setState_dTX for R134a #3275

Merged

beutlich added a commit to beutlich/ModelicaStandardLibrary that referenced this issue Dec 20, 2019

refs modelica#3163: Fix specific enthalpy in two-phase

cf08313

beutlich pushed a commit to beutlich/ModelicaStandardLibrary that referenced this issue Jan 8, 2020

refs modelica#3163: Fix two-phase calculation of setState_dTX

0c0c915

beutlich added a commit to beutlich/ModelicaStandardLibrary that referenced this issue Jan 8, 2020

refs modelica#3163: Fix specific enthalpy in two-phase

84fb108

beutlich added a commit to beutlich/ModelicaStandardLibrary that referenced this issue Jan 8, 2020

refs modelica#3163: Use calculated saturation pressure

8759bef

dietmarw closed this as completed in #3275 Jan 9, 2020

beutlich added a commit to beutlich/ModelicaStandardLibrary that referenced this issue Jan 9, 2020

refs modelica#3163: Fix two-phase calculation of setState_dTX for R134a

0bd93e4

beutlich mentioned this issue Jan 9, 2020

Fix two-phase calculation of setState_dTX for R134a #3298

Merged

beutlich added a commit that referenced this issue Jan 10, 2020

refs #3163: Fix two-phase calculation of setState_dTX for R134a (#3298)

1e613c3

beutlich mentioned this issue Jan 14, 2020

Erroneous two-phase calculation of function setState_psX for R134a medium #3320

Closed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Modelica.Media.R134a problem #3163

Modelica.Media.R134a problem #3163

beutlich commented Oct 21, 2019

HansOlsson commented Oct 21, 2019

beutlich commented Oct 21, 2019

HansOlsson commented Oct 21, 2019

beutlich commented Dec 4, 2019

casella commented Dec 4, 2019

wischhusen commented Dec 17, 2019 •

edited by beutlich

HansOlsson commented Dec 17, 2019

beutlich commented Dec 18, 2019

wischhusen commented Dec 20, 2019 •

edited by beutlich

casella commented Dec 20, 2019

wischhusen commented Dec 20, 2019

casella commented Dec 20, 2019

Modelica.Media.R134a problem #3163

Modelica.Media.R134a problem #3163

Comments

beutlich commented Oct 21, 2019

HansOlsson commented Oct 21, 2019

beutlich commented Oct 21, 2019

HansOlsson commented Oct 21, 2019

beutlich commented Dec 4, 2019

casella commented Dec 4, 2019

wischhusen commented Dec 17, 2019 • edited by beutlich

HansOlsson commented Dec 17, 2019

beutlich commented Dec 18, 2019

wischhusen commented Dec 20, 2019 • edited by beutlich

casella commented Dec 20, 2019

wischhusen commented Dec 20, 2019

casella commented Dec 20, 2019

wischhusen commented Dec 17, 2019 •

edited by beutlich

wischhusen commented Dec 20, 2019 •

edited by beutlich