Fix Modelica.Media.R134a.R134a_ph.setState_pTX (#4626)

* Fixed R134a dofpT() to always converge in all cases

* Added test of fixed R134a to ModelicaTest. The test is based on asserts, no need to generate reference results.

* Remove trailing whitespace.

* Add empty comparisonSignals.txt for test intended to be tested with asserts only.

---------

Co-authored-by: Malte Lenz <maltel@wolfram.com>

Author: casella

Modelica/Media/R134a.mo

@@ -2473,6 +2473,12 @@ This function integrates the derivative of temperature w.r.t. time in order to a
       output SI.Density d "Density";
 
     protected
+      constant Real p_breaks[:]=R134aData.pbreaks
+        "Grid points of reduced pressure";
+      constant Real dl_coef[:, 4]=R134aData.dlcoef
+        "Coefficients of cubic spline for rho_liq(p)";
+      constant Real dv_coef[:, 4]=R134aData.dvcoef
+        "Coefficients of cubic spline for rho_vap(p)";
       constant Real T_breaks[:]=R134aData.Tbreaks
         "Grid points of reduced temperature";
       constant Real dlt_coef[:, 4]=R134aData.dltcoef
@@ -2481,8 +2487,10 @@ This function integrates the derivative of temperature w.r.t. time in order to a
         "Coefficients of cubic spline for rho_vap(T)";
 
       Boolean liquid "Is liquid";
-      Boolean supercritical "Is supercritcal";
+      Boolean supercritical "Is supercritical (p > pcrit)";
+      Boolean highT "= true for high temperature";
       Integer int "Interval number";
+      Real pred "Reduced pressure";
       Real Tred "Reduced temperature";
       Real localx "Abscissa of local spline";
       Integer i "Newton iteration number";
@@ -2500,9 +2508,18 @@ This function integrates the derivative of temperature w.r.t. time in order to a
       i := 0;
       error := 0;
       found := false;
+      pred := p/R134aData.data.FPCRIT;
       Tred := T/R134aData.data.FTCRIT;
-      (int,error) := Common.FindInterval(Tred, T_breaks);
-      localx := Tred - T_breaks[int];
+      highT := (Tred > 0.95);
+      // Guess values for density computed as function of p for high temperature, as function of T for low Temperature
+      // this improves convergence properties (see discussion in #3695)
+      if highT then
+        (int,error) := Common.FindInterval(pred, p_breaks);
+        localx := pred - p_breaks[int];
+      else
+        (int,error) := Common.FindInterval(Tred, T_breaks);
+        localx := Tred - T_breaks[int];
+      end if;
       // set decent initial guesses for d and T
       supercritical := p > R134aData.data.FPCRIT;
       if supercritical then
@@ -2511,11 +2528,11 @@ This function integrates the derivative of temperature w.r.t. time in order to a
       else
         liquid := T <= Modelica.Media.R134a.R134a_ph.saturationTemperature(p);
         if liquid then
-          d := R134aData.data.FDCRIT*Common.CubicSplineEval(localx, dlt_coef[int,
-            1:4])*1.02;
+          d := if highT then R134aData.data.FDCRIT*Common.CubicSplineEval(localx, dl_coef[int,1:4])*1.02
+                        else R134aData.data.FDCRIT*Common.CubicSplineEval(localx, dlt_coef[int,1:4])*1.02;
         else
-          d := R134aData.data.FDCRIT*Common.CubicSplineEval(localx, dvt_coef[int,
-            1:4])*0.95;
+          d := if highT then R134aData.data.FDCRIT*Common.CubicSplineEval(localx, dv_coef[int,1:4])*0.95
+                        else R134aData.data.FDCRIT*Common.CubicSplineEval(localx, dvt_coef[int,1:4])*0.95;
         end if;
       end if;
 

ModelicaTest/Media.mo

@@ -971,6 +971,49 @@ is given to compare the approximation.
       annotation (experiment(StopTime=1));
     end R134a_setState_phX;
 
+    model R134a_pTX_phX_all "Tests the R134a_ph.setState_pTX and setState_phX functions for consistency over the entire valid range of p and T"
+      extends Modelica.Icons.Example;
+      package Medium = Modelica.Media.R134a.R134a_ph;
+      parameter Medium.AbsolutePressure pmin = 0.0042e5;
+      parameter Medium.AbsolutePressure pmax = 700e5;
+      parameter Modelica.Units.SI.PerUnit p_increase = 2;
+      parameter Medium.Temperature Tmin = 190;
+      parameter Medium.Temperature Tmax = 450;
+
+      Medium.ThermodynamicState state = Medium.setState_pTX(p,T);
+      Medium.Temperature Tsat = Medium.saturationTemperature(p);
+      Medium.AbsolutePressure p;
+      Medium.Temperature T;
+      Medium.Density d = Medium.density(state);
+      Medium.SpecificEnthalpy h = Medium.specificEnthalpy(state);
+      Medium.ThermalConductivity k = Medium.thermalConductivity(state);
+
+      Medium.ThermodynamicState state_check = Medium.setState_phX(p,h);
+      Medium.Temperature T_check = Medium.temperature(state_check);
+      Medium.Density d_check = Medium.density(state_check);
+
+      Modelica.Units.SI.PerUnit dT;
+    initial equation
+      p = pmin;
+      T = Tmin;
+      dT = 100*45;
+    equation
+      der(T) = dT;
+      when abs(T - Tsat) < 2 then
+        // skip region close to Tsat
+        reinit(T, Tsat+4*sign(dT));
+      end when;
+      when {T > Tmax, T < Tmin} then
+        dT = -pre(dT);
+        p = pre(p)*p_increase;
+      end when;
+      assert(0.7e5 < h and h < 6e5, "Specific enthalpy h out of bounds");
+      assert(d > 0 and d < 1600, "Density out of bounds");
+      assert(abs(T - T_check)/T < 1e-3 and abs(d - d_check)/d < 1e-2, "Convergence error");
+    annotation(
+        experiment(StartTime = 0, StopTime = 1, Tolerance = 1e-06, Interval = 0.0001));
+    end R134a_pTX_phX_all;
+
     model WaterIF97_dewEnthalpy "Test dewEnthalpy of WaterIF97"
       extends Modelica.Icons.Example;
       replaceable package Medium = Modelica.Media.Water.StandardWater "Medium model";

ModelicaTest/Resources/Reference/ModelicaTest/Media/TestOnly/R134a_pTX_phX_all/comparisonSignals.txt

@@ -0,0 +1 @@
+time