Various mixture-related fluid properties improvements

modules/fluid_properties/include/fluidproperties/VaporMixtureFluidProperties.h

@@ -116,6 +116,7 @@ class VaporMixtureFluidProperties : public FluidProperties
   propfunc(c, v, e)
   propfunc(rho, p, T)
   propfunc(e, p, T)
+  propfunc(s, p, T)
   propfunc(c, p, T)
   propfunc(cp, p, T)
   propfunc(cv, p, T)

modules/fluid_properties/include/userobjects/FluidPropertiesInterrogator.h

@@ -173,6 +173,8 @@ class FluidPropertiesInterrogator : public GeneralUserObject
   const bool _has_2phase;
   /// flag that user provided 2-phase NCG fluid properties
   const bool _has_2phase_ncg;
+  /// flag that user provided 2-phase NCG partial pressure fluid properties
+  const bool _has_2phase_ncg_partial_pressure;
 
   /// pointer to liquid fluid properties object (if provided 2-phase object)
   const SinglePhaseFluidProperties * const _fp_liquid;

modules/fluid_properties/include/utils/BrentsMethod.h

@@ -16,6 +16,13 @@
 // C++ includes
 #include <functional>
 
+/**
+ * Brent's method is used to find the root of a function f(x), i.e., find
+ * x such that f(x) = 0.
+ *
+ * First, brackets x1 and x2 are found such that f(x) changes sign between
+ * x1 and x2, implying that there is a root between the points.
+ */
 namespace BrentsMethod
 {
 /**

modules/fluid_properties/src/fluidproperties/SinglePhaseFluidProperties.C

@@ -309,13 +309,13 @@ SinglePhaseFluidProperties::criticalTemperature() const
 Real
 SinglePhaseFluidProperties::criticalDensity() const
 {
-  mooseError(__PRETTY_FUNCTION__, " not implemented.");
+  return rho_from_p_T(criticalPressure(), criticalTemperature());
 }
 
 Real
 SinglePhaseFluidProperties::criticalInternalEnergy() const
 {
-  mooseError(__PRETTY_FUNCTION__, " not implemented.");
+  return e_from_p_rho(criticalPressure(), criticalDensity());
 }
 
 Real

modules/fluid_properties/src/userobjects/FluidPropertiesInterrogator.C

@@ -53,6 +53,7 @@ FluidPropertiesInterrogator::FluidPropertiesInterrogator(const InputParameters &
     _has_vapor_mixture(dynamic_cast<const VaporMixtureFluidProperties * const>(_fp)),
     _has_2phase(_fp_2phase),
     _has_2phase_ncg(_fp_2phase_ncg),
+    _has_2phase_ncg_partial_pressure(_fp_2phase_ncg_partial_pressure),
     _fp_liquid(_has_2phase
                    ? &getUserObjectByName<SinglePhaseFluidProperties>(_fp_2phase->getLiquidName())
                    : nullptr),
@@ -387,7 +388,9 @@ FluidPropertiesInterrogator::computeVaporMixture(bool throw_error_if_no_match)
 
   // determine how state is specified
   std::vector<std::vector<std::string>> parameter_sets = {
-      {"p", "T", "x_ncg"}, {"rho", "e", "x_ncg"}, {"rho", "rhou", "rhoE", "x_ncg"}, {"p", "T"}};
+      {"p", "T", "x_ncg"}, {"rho", "e", "x_ncg"}, {"rho", "rhou", "rhoE", "x_ncg"}};
+  if (_has_2phase_ncg_partial_pressure)
+    parameter_sets.push_back({"p", "T"});
   auto specified = getSpecifiedSetMap(parameter_sets, "vapor mixture", throw_error_if_no_match);
 
   // compute/determine rho, e, p, T, vel
@@ -420,7 +423,7 @@ FluidPropertiesInterrogator::computeVaporMixture(bool throw_error_if_no_match)
 
     specified_a_set = true;
   }
-  else if (_fp_2phase_ncg_partial_pressure && specified["p,T"])
+  else if (_has_2phase_ncg_partial_pressure && specified["p,T"])
   {
     p = getParam<Real>("p");
     T = getParam<Real>("T");
@@ -456,6 +459,7 @@ FluidPropertiesInterrogator::computeVaporMixture(bool throw_error_if_no_match)
 
     const Real v = 1.0 / rho;
     const Real h = e + p / rho;
+    const Real s = _fp_vapor_mixture->s_from_p_T(p, T, x_ncg);
     const Real c = _fp_vapor_mixture->c_from_p_T(p, T, x_ncg);
     const Real mu = _fp_vapor_mixture->mu_from_p_T(p, T, x_ncg);
     const Real cp = _fp_vapor_mixture->cp_from_p_T(p, T, x_ncg);
@@ -469,6 +473,7 @@ FluidPropertiesInterrogator::computeVaporMixture(bool throw_error_if_no_match)
     params.set<Real>("e") = e;
     params.set<Real>("v") = v;
     params.set<Real>("h") = h;
+    params.set<Real>("s") = s;
     params.set<Real>("c") = c;
     params.set<Real>("mu") = mu;
     params.set<Real>("cp") = cp;
@@ -534,7 +539,7 @@ void
 FluidPropertiesInterrogator::buildJSONVaporMixture(nlohmann::json & json,
                                                    const InputParameters & params)
 {
-  for (auto p : {"p", "T", "rho", "e", "v", "h", "c", "mu", "cp", "cv", "k"})
+  for (auto p : {"p", "T", "rho", "e", "v", "h", "s", "c", "mu", "cp", "cv", "k"})
     if (params.have_parameter<Real>(p))
       json["static"][p] = params.get<Real>(p);
 
@@ -761,6 +766,7 @@ FluidPropertiesInterrogator::outputVaporMixtureStaticProperties(const InputParam
   outputProperty("Specific volume", params.get<Real>("v"), "m^3/kg");
   outputProperty("Specific internal energy", params.get<Real>("e"), "J/kg");
   outputProperty("Specific enthalpy", params.get<Real>("h"), "J/kg");
+  outputProperty("Specific entropy", params.get<Real>("s"), "J/kg");
   _console << std::endl;
   outputProperty("Sound speed", params.get<Real>("c"), "m/s");
   outputProperty("Dynamic viscosity", params.get<Real>("mu"), "Pa-s");

modules/fluid_properties/src/utils/BrentsMethod.C

@@ -34,8 +34,11 @@ bracket(std::function<Real(Real)> const & f, Real & x1, Real & x2)
   if (f1 * f2 > 0.0)
   {
     unsigned int iter = 0;
+    std::stringstream debug_ss;
     while (f1 * f2 > 0.0)
     {
+      debug_ss << "  iteration " << iter << ": (x1,x2) = (" << x1 << "," << x2 << "), (f1,f2) = ("
+               << f1 << "," << f2 << ")\n";
       if (std::abs(f1) < std::abs(f2))
       {
         x1 += factor * (x1 - x2);
@@ -52,7 +55,7 @@ bracket(std::function<Real(Real)> const & f, Real & x1, Real & x2)
       iter++;
       if (iter >= n)
         throw MooseException("No bracketing interval found by BrentsMethod::bracket after " +
-                             Moose::stringify(n) + " iterations");
+                             Moose::stringify(n) + " iterations:\n" + debug_ss.str());
     }
   }
 }
@@ -71,8 +74,11 @@ root(std::function<Real(Real)> const & f, Real x1, Real x2, Real tol)
     throw MooseException("Root must be bracketed in BrentsMethod::root");
 
   fc = fb;
+  std::stringstream debug_ss;
   for (unsigned int i = 1; i <= iter_max; ++i)
   {
+    debug_ss << "  iteration " << i << ": dx = " << xm << ", x = " << b << ", f(x) = " << fb
+             << "\n";
     if (fb * fc > 0.0)
     {
       // Rename a,b and c and adjust bounding interval d
@@ -154,7 +160,8 @@ root(std::function<Real(Real)> const & f, Real x1, Real x2, Real tol)
     fb = f(b);
   }
 
-  throw MooseException("Maximum number of iterations exceeded in BrentsMethod::root");
+  throw MooseException("Maximum number of iterations exceeded in BrentsMethod::root:\n" +
+                       debug_ss.str());
   return 0.0; // Should never get here
 }
 } // namespace BrentsMethod

modules/fluid_properties/test/include/utils/FluidPropertiesTestUtils.h

@@ -28,7 +28,7 @@
 // difference approximation
 #define REL_TOL_DERIVATIVE 1e-6
 
-// Macro for computing relative error
+// Macro for performing relative error test
 #define REL_TEST(value, ref_value, tol)                                                            \
   if (std::abs(ref_value) < 1e-15)                                                                 \
     ABS_TEST(value, ref_value, tol);                                                               \
@@ -37,6 +37,15 @@
                        (MetaPhysicL::raw_value(ref_value))),                                       \
               tol);
 
+// Macro for performing relative or absolute error test
+#define REL_ABS_TEST(value, ref_value, rel_tol, abs_tol)                                           \
+  if (std::abs(ref_value - value) < abs_tol)                                                       \
+    EXPECT_TRUE(true);                                                                             \
+  else                                                                                             \
+    EXPECT_LE(std::abs(((MetaPhysicL::raw_value(value)) - (MetaPhysicL::raw_value(ref_value))) /   \
+                       (MetaPhysicL::raw_value(ref_value))),                                       \
+              rel_tol);
+
 // Macro for computing absolute error
 #define ABS_TEST(value, ref_value, tol)                                                            \
   EXPECT_LE(std::abs(((MetaPhysicL::raw_value(value)) - (MetaPhysicL::raw_value(ref_value)))),     \

modules/fluid_properties/test/include/utils/VaporMixtureFluidPropertiesUtils.h

@@ -12,8 +12,9 @@
 #include "SinglePhaseFluidPropertiesTestUtils.h"
 
 // Macro for performing a derivative test
-#define VAPOR_MIX_DERIV_TEST(f, a, b, x, tol)                                                      \
+#define VAPOR_MIX_DERIV_TEST(f, a, b, x, rel_tol)                                                  \
   {                                                                                                \
+    const Real abs_tol = 1e-10;                                                                    \
     const Real da = REL_PERTURBATION * a;                                                          \
     const Real db = REL_PERTURBATION * b;                                                          \
     std::vector<Real> dx(x.size());                                                                \
@@ -33,9 +34,9 @@
     Real f_value, df_da, df_db;                                                                    \
     std::vector<Real> df_dx(x.size());                                                             \
     f(a, b, x, f_value, df_da, df_db, df_dx);                                                      \
-    REL_TEST(f(a, b, x), f_value, REL_TOL_CONSISTENCY);                                            \
-    REL_TEST(df_da, df_da_fd, tol);                                                                \
-    REL_TEST(df_db, df_db_fd, tol);                                                                \
+    REL_ABS_TEST(f(a, b, x), f_value, REL_TOL_CONSISTENCY, abs_tol);                               \
+    REL_ABS_TEST(df_da, df_da_fd, rel_tol, abs_tol);                                               \
+    REL_ABS_TEST(df_db, df_db_fd, rel_tol, abs_tol);                                               \
     for (unsigned int i = 0; i < x.size(); ++i)                                                    \
-      REL_TEST(df_dx[i], df_dx_fd[i], tol);                                                        \
+      REL_ABS_TEST(df_dx[i], df_dx_fd[i], rel_tol, abs_tol);                                       \
   }

modules/fluid_properties/test/tests/fp_interrogator/tests

@@ -329,6 +329,7 @@ Density:                                   0.8972309616  kg/m\^3
 Specific volume:                            1.114540227  m\^3/kg
 Specific internal energy:              1.0334658129e\+06  J/kg
 Specific enthalpy:                     1.1449198356e\+06  J/kg
+Specific entropy:                      1.4838960040e\+04  J/kg
 
 Sound speed:                                351.3883482  m/s
 Dynamic viscosity:                     1.8230000000e-05  Pa-s
@@ -374,7 +375,7 @@ Thermal conductivity:                           0.02568  W/\(m-K\)'
                  '"k":0.[0-9]+,"mu":1.823e-05,"p":100000.0,"rho":0.[0-9]+,"s":2486.[0-9]+,'
                  '"v":1.[0-9]+}},"vapor-mixture":{"static":{"T":372.[0-9]+,"c":351.[0-9]+,'
                  '"cp":3071.[0-9]+,"cv":2772.[0-9]+,"e":1033465.[0-9]+,"h":1144919.[0-9]+,'
-                 '"k":0.[0-9]+,"mu":1.[0-9]+e-05,"p":100000.0,"rho":0.[0-9]+,"v":1.[0-9]+}}}'
+                 '"k":0.[0-9]+,"mu":1.[0-9]+e-05,"p":100000.0,"rho":0.[0-9]+,"s":14838.[0-9]+,"v":1.[0-9]+}}}'
 
     requirement = "The fluid properties interrogator shall output two-phase "
       "and static-state, single-phase fluid properties for (p, T) input with "
@@ -399,6 +400,7 @@ Density:                                    1.187005237  kg/m\^3
 Specific volume:                           0.8424562661  m\^3/kg
 Specific internal energy:              2.4771659033e\+06  J/kg
 Specific enthalpy:                     3.2387829780e\+06  J/kg
+Specific entropy:                           5966.595696  J/kg
 
 Sound speed:                                997.8878004  m/s
 Dynamic viscosity:                     1.8230000000e-05  Pa-s
@@ -422,7 +424,7 @@ Thermal conductivity:                           0.02568  W/\(m-K\)'
     threading = '!pthreads'
     expect_out = '{"static":{"T":2547.[0-9]+,"c":997.[0-9]+,"cp":1271.[0-9]+,'
                  '"cv":972.[0-9]+,"e":2477165.[0-9]+,"h":3238782.[0-9]+,"k":0.[0-9]+,'
-                 '"mu":1.[0-9]+e-05,"p":904043.[0-9]+,"rho":1.[0-9]+,"v":0.[0-9]+}}'
+                 '"mu":1.[0-9]+e-05,"p":904043.[0-9]+,"rho":1.[0-9]+,"s":5966.[0-9]+,"v":0.[0-9]+}}'
 
     requirement = "The fluid properties interrogator shall output static-state, "
       "single-phase fluid properties for (rho, e) input with a vapor mixture "