From 0ad51993616f2e97bdcc99cfd5e56ae9e1a091c0 Mon Sep 17 00:00:00 2001 From: Guillaume Giudicelli Date: Sat, 4 Jan 2025 15:23:08 +0100 Subject: [PATCH 1/6] Add an option for variable cp for TP function-fp Use a different definition to compute cp from cv refs #21245 --- ...peraturePressureFunctionFluidProperties.md | 29 +- ...mperaturePressureFunctionFluidProperties.h | 11 +- ...mperaturePressureFunctionFluidProperties.C | 306 ++++++++++++++---- 3 files changed, 277 insertions(+), 69 deletions(-) diff --git a/modules/fluid_properties/doc/content/source/fluidproperties/TemperaturePressureFunctionFluidProperties.md b/modules/fluid_properties/doc/content/source/fluidproperties/TemperaturePressureFunctionFluidProperties.md index 943edc755514..9ee80b27c3fb 100644 --- a/modules/fluid_properties/doc/content/source/fluidproperties/TemperaturePressureFunctionFluidProperties.md +++ b/modules/fluid_properties/doc/content/source/fluidproperties/TemperaturePressureFunctionFluidProperties.md @@ -12,12 +12,31 @@ functions. The following relations hold true: \rho = \rho^{user}(time=0, (x=T, y=P, z=0)) \\ \mu = \mu^{user}(time=0, (x=T, y=P, z=0)) \\ k = k^{user}(time=0, (x=T, y=P, z=0)) \\ +\end{aligned} +\end{equation} + +There are two options for specific heat. Either the user sets a constant specific isochoric heat capacity + +\begin{equation} +\begin{aligned} c_v = c_v^{user} \\ - e = c_v * T + e = e_{ref} + c_v * (T - T_{ref}) +\end{aligned} +\end{equation} + +Or, the user uses a function of temperature and pressure (same arguments as for density) + +\begin{equation} +\begin{aligned} + x = T \\ + y = P \\ + cp = cp^{user}(time=0, (x=T, y=P, z=0)) + cv = cp - \dfrac{\alpha ^ 2 T}{\rho \beta_T} \end{aligned} \end{equation} -with $T$ the temperature, $P$ the pressure, $\rho$ the density, $\mu$ the dynamic viscosity, $k$ the thermal conductivity, $c_v$ the specific isochoric heat capacity, $e$ the specific energy and the $user$ exponent indicating a user-passed parameter. +with $T$ the temperature, $P$ the pressure, $\rho$ the density, $\mu$ the dynamic viscosity, $k$ the thermal conductivity, $c_v$ the specific isochoric heat capacity, $e$ the specific energy, $\alpha$ the coefficient of thermal expansion, $\beta_T$ the isothermal +compressibility, and the $user$ exponent indicating a user-passed parameter. The derivatives of the fluid properties are obtained using the `Function`(s) gradient components and the appropriate derivative chaining for derived properties. @@ -31,9 +50,9 @@ Support for the conservative (specific volume, internal energy) variable set is partial. Notable missing implementations are routines for entropy, the speed of sound, and some conversions between specific enthalpy and specific energy. -!alert warning -Due to the approximations made when computing the isobaric heat capacity from the constant -isochoric heat capacity, this material should only be used for nearly-incompressible fluids. +!alert note +When using a function for the isobaric specific heat capacity, Simpson's rule is used to compute +the specific energy from the isochoric specific heat capacity. This limits the accuracy of the calculation. ## Example Input File Syntax diff --git a/modules/fluid_properties/include/fluidproperties/TemperaturePressureFunctionFluidProperties.h b/modules/fluid_properties/include/fluidproperties/TemperaturePressureFunctionFluidProperties.h index 5901ae2a244f..c9343be517b4 100644 --- a/modules/fluid_properties/include/fluidproperties/TemperaturePressureFunctionFluidProperties.h +++ b/modules/fluid_properties/include/fluidproperties/TemperaturePressureFunctionFluidProperties.h @@ -348,6 +348,15 @@ class TemperaturePressureFunctionFluidProperties : public SinglePhaseFluidProper /// function defining dynamic viscosity as a function of temperature and pressure const Function * _mu_function; + /// function defining specific heat as a function of temperature and pressure + const Function * _cp_function; + /// constant isochoric specific heat - const Real & _cv; + const Real _cv; + /// whether a constant isochoric specific heat is used + const bool _cv_is_constant; + /// Reference specific energy + const Real _e_ref; + /// Reference temperature for the reference specific energy + const Real _T_ref; }; diff --git a/modules/fluid_properties/src/fluidproperties/TemperaturePressureFunctionFluidProperties.C b/modules/fluid_properties/src/fluidproperties/TemperaturePressureFunctionFluidProperties.C index 3a1696a02d51..d4146b96d166 100644 --- a/modules/fluid_properties/src/fluidproperties/TemperaturePressureFunctionFluidProperties.C +++ b/modules/fluid_properties/src/fluidproperties/TemperaturePressureFunctionFluidProperties.C @@ -23,8 +23,13 @@ TemperaturePressureFunctionFluidProperties::validParams() params.addRequiredParam( "mu", "Dynamic viscosity function of temperature and pressure [Pa-s]"); - params.addRequiredRangeCheckedParam( - "cv", "cv > 0", "Constant isochoric specific heat [J/(kg-K)]"); + params.addParam( + "cp", "Isobaric specific heat function of temperature and pressure [J/(kg-K)]"); + params.addRangeCheckedParam( + "cv", 0, "cv >= 0", "Constant isochoric specific heat [J/(kg-K)]"); + params.addParam("e_ref", 0, "Specific energy at the reference temperature"); + params.addParam("T_ref", 0, "Reference temperature for the specific energy"); + params.addClassDescription( "Single-phase fluid properties that allows to provide thermal " "conductivity, density, and viscosity as functions of temperature and pressure."); @@ -33,8 +38,15 @@ TemperaturePressureFunctionFluidProperties::validParams() TemperaturePressureFunctionFluidProperties::TemperaturePressureFunctionFluidProperties( const InputParameters & parameters) - : SinglePhaseFluidProperties(parameters), _initialized(false), _cv(getParam("cv")) + : SinglePhaseFluidProperties(parameters), + _initialized(false), + _cv(getParam("cv")), + _cv_is_constant(_cv != 0), + _e_ref(getParam("e_ref")), + _T_ref(getParam("T_ref")) { + if (isParamValid("cp") && _cv_is_constant) + paramError("cp", "Can only set a function for cp if a constant is not specified for cv."); } void @@ -43,6 +55,7 @@ TemperaturePressureFunctionFluidProperties::initialSetup() _k_function = &getFunction("k"); _rho_function = &getFunction("rho"); _mu_function = &getFunction("mu"); + _cp_function = isParamValid("cp") ? &getFunction("cp") : nullptr; _initialized = true; } @@ -53,9 +66,22 @@ TemperaturePressureFunctionFluidProperties::fluidName() const } Real -TemperaturePressureFunctionFluidProperties::T_from_v_e(Real /* v */, Real e) const -{ - return e / _cv; +TemperaturePressureFunctionFluidProperties::T_from_v_e(Real v, Real e) const +{ + if (_cv_is_constant) + return _T_ref + (e - _e_ref) / _cv; + else + { + const Real p0 = _p_initial_guess; + const Real T0 = _T_initial_guess; + Real p, T; + bool conversion_succeeded = true; + p_T_from_v_e(v, e, p0, T0, p, T, conversion_succeeded); + if (conversion_succeeded) + return T; + else + mooseError("T_from_v_e calculation failed."); + } } Real @@ -97,9 +123,24 @@ TemperaturePressureFunctionFluidProperties::T_from_p_rho(Real p, Real rho) const Real TemperaturePressureFunctionFluidProperties::cp_from_v_e(Real v, Real e) const { - Real p = p_from_v_e(v, e); - Real T = T_from_v_e(v, e); - return cp_from_p_T(p, T); + if (_cv_is_constant) + { + Real p = p_from_v_e(v, e); + Real T = T_from_v_e(v, e); + return cp_from_p_T(p, T); + } + else + { + const Real p0 = _p_initial_guess; + const Real T0 = _T_initial_guess; + Real p, T; + bool conversion_succeeded = true; + p_T_from_v_e(v, e, p0, T0, p, T, conversion_succeeded); + if (conversion_succeeded) + return cp_from_p_T(p, T); + else + mooseError("cp_from_v_e calculation failed. p= ", p, " T = ", T); + } } void @@ -108,7 +149,7 @@ TemperaturePressureFunctionFluidProperties::cp_from_v_e( { cp = cp_from_v_e(v, e); // Using finite difference to get around difficulty of implementation - Real eps = 1e-8; + Real eps = 1e-10; Real cp_pert = cp_from_v_e(v * (1 + eps), e); dcp_dv = (cp_pert - cp) / eps / v; cp_pert = cp_from_v_e(v, e * (1 + eps)); @@ -116,18 +157,63 @@ TemperaturePressureFunctionFluidProperties::cp_from_v_e( } Real -TemperaturePressureFunctionFluidProperties::cv_from_v_e(Real /* v */, Real /* e */) const -{ - return _cv; +TemperaturePressureFunctionFluidProperties::cv_from_v_e(Real v, Real e) const +{ + if (_cv_is_constant) + return _cv; + else + { + const Real p0 = _p_initial_guess; + const Real T0 = _T_initial_guess; + Real p, T; + bool conversion_succeeded = true; + p_T_from_v_e(v, e, p0, T0, p, T, conversion_succeeded); + if (conversion_succeeded) + return cv_from_p_T(p, T); + else + mooseError("cp_from_v_e calculation failed."); + } } void TemperaturePressureFunctionFluidProperties::cv_from_v_e( Real v, Real e, Real & cv, Real & dcv_dv, Real & dcv_de) const { - cv = cv_from_v_e(v, e); - dcv_dv = 0.0; - dcv_de = 0.0; + if (_cv_is_constant) + { + cv = cv_from_v_e(v, e); + dcv_dv = 0.0; + dcv_de = 0.0; + } + else + { + const Real p0 = _p_initial_guess; + const Real T0 = _T_initial_guess; + Real p, T; + bool conversion_succeeded = true; + p_T_from_v_e(v, e, p0, T0, p, T, conversion_succeeded); + Real dcv_dp, dcv_dT; + cv_from_p_T(p, T, cv, dcv_dp, dcv_dT); + if (!conversion_succeeded) + mooseError("cp_from_v_e and derivatives calculation failed."); + + Real p1, T1; + p_T_from_v_e(v * (1 + 1e-6), e, p0, T0, p1, T1, conversion_succeeded); + Real dp_dv = (p1 - p) / (v * 1e-6); + Real dT_dv = (T1 - T) / (v * 1e-6); + if (!conversion_succeeded) + mooseError("cp_from_v_e and derivatives calculation failed."); + + Real p2, T2; + p_T_from_v_e(v, e * (1 + 1e-6), p0, T0, p2, T2, conversion_succeeded); + Real dp_de = (p2 - p) / (e * 1e-6); + Real dT_de = (T2 - T) / (e * 1e-6); + if (!conversion_succeeded) + mooseError("cp_from_v_e and derivatives calculation failed."); + + dcv_dv = dcv_dp * dp_dv + dcv_dT * dT_dv; + dcv_de = dcv_dp * dp_de + dcv_dT * dT_de; + } } Real @@ -153,17 +239,47 @@ TemperaturePressureFunctionFluidProperties::p_from_v_e(Real v, Real e) const Real TemperaturePressureFunctionFluidProperties::mu_from_v_e(Real v, Real e) const { - Real temperature = T_from_v_e(v, e); - Real pressure = p_from_v_e(v, e); - return mu_from_p_T(pressure, temperature); + if (_cv_is_constant) + { + Real temperature = T_from_v_e(v, e); + Real pressure = p_from_v_e(v, e); + return mu_from_p_T(pressure, temperature); + } + else + { + const Real p0 = _p_initial_guess; + const Real T0 = _T_initial_guess; + Real p, T; + bool conversion_succeeded = true; + p_T_from_v_e(v, e, p0, T0, p, T, conversion_succeeded); + if (conversion_succeeded) + return mu_from_p_T(p, T); + else + mooseError("mu_from_v_e calculation failed."); + } } Real TemperaturePressureFunctionFluidProperties::k_from_v_e(Real v, Real e) const { - Real temperature = T_from_v_e(v, e); - Real pressure = p_from_v_e(v, e); - return k_from_p_T(pressure, temperature); + if (_cv_is_constant) + { + Real temperature = T_from_v_e(v, e); + Real pressure = p_from_v_e(v, e); + return k_from_p_T(pressure, temperature); + } + else + { + const Real p0 = _p_initial_guess; + const Real T0 = _T_initial_guess; + Real p, T; + bool conversion_succeeded = true; + p_T_from_v_e(v, e, p0, T0, p, T, conversion_succeeded); + if (conversion_succeeded) + return k_from_p_T(p, T); + else + mooseError("k_from_v_e calculation failed."); + } } Real @@ -239,25 +355,65 @@ TemperaturePressureFunctionFluidProperties::h_from_p_T( } Real -TemperaturePressureFunctionFluidProperties::e_from_p_T(Real /* pressure */, Real temperature) const -{ - return _cv * temperature; +TemperaturePressureFunctionFluidProperties::e_from_p_T(Real pressure, Real temperature) const +{ + if (_cv_is_constant) + return _e_ref + _cv * (temperature - _T_ref); + else + { + // Use Simpson's 3/8 rule to integrate cv + Real h = (temperature - _T_ref) / 3; + Real integral = + 3. / 8. * h * + (cv_from_p_T(pressure, _T_ref) + 3 * cv_from_p_T(pressure, _T_ref + h) + + 3 * cv_from_p_T(pressure, _T_ref + 2 * h) + cv_from_p_T(pressure, temperature)); + return _e_ref + integral; + } } void TemperaturePressureFunctionFluidProperties::e_from_p_T( Real pressure, Real temperature, Real & e, Real & de_dp, Real & de_dT) const { - e = e_from_p_T(pressure, temperature); - de_dp = 0.0; - de_dT = _cv; + if (_cv_is_constant) + { + e = e_from_p_T(pressure, temperature); + de_dp = 0.0; + de_dT = _cv; + } + else + { + e = e_from_p_T(pressure, temperature); + // Use Simpson's 3/8 rule to integrate cv then take the derivative + Real h = (temperature - _T_ref) / 3; + Real cv0, dcv_dp0, dcv_dT0; + cv_from_p_T(pressure, _T_ref, cv0, dcv_dp0, dcv_dT0); + Real cv1, dcv_dp1, dcv_dT1; + cv_from_p_T(pressure, _T_ref + h, cv1, dcv_dp1, dcv_dT1); + Real cv2, dcv_dp2, dcv_dT2; + cv_from_p_T(pressure, _T_ref + 2 * h, cv2, dcv_dp2, dcv_dT2); + Real cv3, dcv_dp3, dcv_dT3; + cv_from_p_T(pressure, temperature, cv3, dcv_dp3, dcv_dT3); + de_dp = 3. / 8 * h * (dcv_dp0 + 3 * dcv_dp1 + 3 * dcv_dp2 + dcv_dp3); + // Consistency with how the energy is computed is more important here + de_dT = + 1. / 8. * (cv0 + 3 * cv1 + 3 * cv2 + cv3) + 3. / 8 * h * (dcv_dT1 + 2 * dcv_dT2 + dcv_dT3); + } } Real TemperaturePressureFunctionFluidProperties::e_from_p_rho(Real p, Real rho) const { - Real temperature = T_from_p_rho(p, rho); - return _cv * temperature; + if (_cv_is_constant) + { + Real temperature = T_from_p_rho(p, rho); + return _e_ref + _cv * (temperature - _T_ref); + } + else + { + Real temperature = T_from_p_rho(p, rho); + return e_from_p_T(p, temperature); + } } Real @@ -265,60 +421,84 @@ TemperaturePressureFunctionFluidProperties::beta_from_p_T(Real pressure, Real te { Real rho, drho_dp, drho_dT; rho_from_p_T(pressure, temperature, rho, drho_dp, drho_dT); - return -drho_dT / rho; + return -drho_dp / rho; } Real TemperaturePressureFunctionFluidProperties::cp_from_p_T(Real p, Real T) const { - Real rho, drho_dp, drho_dT; - rho_from_p_T(p, T, rho, drho_dp, drho_dT); - Real dv_dT = -1 / rho / rho * drho_dT; - // neglecting dp_dT term due to difficulty in computing it in the general case - // this is not OK for gases, but should be ok for nearly incompressible fluids - return _cv + p * dv_dT; - // an alternative would be to use finite differencing for the p * v term in its entirety + if (_cv_is_constant) + { + Real rho, drho_dp, drho_dT; + rho_from_p_T(p, T, rho, drho_dp, drho_dT); + Real alpha = -1. / rho * drho_dT; + return _cv + MathUtils::pow(alpha, 2) * T / rho / beta_from_p_T(p, T); + } + else + { + if (!_initialized) + const_cast(this)->initialSetup(); + return _cp_function->value(0, Point(T, p, 0)); + } } void TemperaturePressureFunctionFluidProperties::cp_from_p_T( Real pressure, Real temperature, Real & cp, Real & dcp_dp, Real & dcp_dT) const { - Real e, de_dp, de_dT; - e_from_p_T(pressure, temperature, e, de_dp, de_dT); - Real rho, drho_dp, drho_dT; - rho_from_p_T(pressure, temperature, rho, drho_dp, drho_dT); - Real dv_dT = -1 / rho / rho * drho_dT; - cp = _cv + pressure * dv_dT; - - // use finite difference for second derivative to avoid having to define hessian for the function - Real eps = 1e-8; - Real rho_pert, drho_dp_pert, drho_dT_pert; - rho_from_p_T(pressure, temperature * (1 + eps), rho_pert, drho_dp_pert, drho_dT_pert); - Real d2v_dT2 = - -(-1 / rho / rho * drho_dT + 1 / rho_pert / rho_pert * drho_dT_pert) / eps / temperature; - rho_from_p_T(pressure * (1 + eps), temperature, rho_pert, drho_dp_pert, drho_dT_pert); - Real d2v_dTdp = - -(-1 / rho / rho * drho_dT + 1 / rho_pert / rho_pert * drho_dT_pert) / eps / pressure; - - dcp_dp = dv_dT + pressure * d2v_dTdp; - dcp_dT = pressure * d2v_dT2; + if (_cv_is_constant) + { + cp = cp_from_p_T(pressure, temperature); + Real eps = 1e-8; + Real cp_1p = cp_from_p_T(pressure * (1 + eps), temperature); + Real cp_1T = cp_from_p_T(pressure, temperature * (1 + eps)); + dcp_dp = (cp_1p - cp) / (pressure * eps); + dcp_dT = (cp_1T - cp) / (temperature * eps); + } + else + { + cp = cp_from_p_T(pressure, temperature); + const RealVectorValue grad_function = + _cp_function->gradient(0, Point(temperature, pressure, 0)); + dcp_dT = grad_function(0); + dcp_dp = grad_function(1); + } } Real -TemperaturePressureFunctionFluidProperties::cv_from_p_T(Real /* pressure */, - Real /* temperature */) const +TemperaturePressureFunctionFluidProperties::cv_from_p_T(Real pressure, Real temperature) const { - return _cv; + if (_cv_is_constant) + return _cv; + else + { + Real rho, drho_dp, drho_dT; + rho_from_p_T(pressure, temperature, rho, drho_dp, drho_dT); + Real alpha = -1. / rho * drho_dT; + return cp_from_p_T(pressure, temperature) - + MathUtils::pow(alpha, 2) * temperature / rho / beta_from_p_T(pressure, temperature); + } } void TemperaturePressureFunctionFluidProperties::cv_from_p_T( Real pressure, Real temperature, Real & cv, Real & dcv_dp, Real & dcv_dT) const { - cv = cv_from_p_T(pressure, temperature); - dcv_dp = 0.0; - dcv_dT = 0.0; + if (_cv_is_constant) + { + cv = cv_from_p_T(pressure, temperature); + dcv_dp = 0.0; + dcv_dT = 0.0; + } + else + { + cv = cv_from_p_T(pressure, temperature); + Real eps = 1e-10; + Real cv_1p = cv_from_p_T(pressure * (1 + eps), temperature); + Real cv_1T = cv_from_p_T(pressure, temperature * (1 + eps)); + dcv_dp = (cv_1p - cv) / (pressure * eps); + dcv_dT = (cv_1T - cv) / (temperature * eps); + } } Real From b1509ef5c479e36e439b96645dee55bd9a2d0e8d Mon Sep 17 00:00:00 2001 From: Guillaume Giudicelli Date: Sat, 4 Jan 2025 15:23:32 +0100 Subject: [PATCH 2/6] Add unit tests for the new variable cp option refs #21245 Update regression test as cp expression changed --- .../gold/example_out.csv | 2 +- ...aturePressureFunctionFluidPropertiesTest.h | 19 +++ ...aturePressureFunctionFluidPropertiesTest.C | 146 ++++++++++++++++-- 3 files changed, 154 insertions(+), 13 deletions(-) diff --git a/modules/fluid_properties/test/tests/temperature_pressure_function/gold/example_out.csv b/modules/fluid_properties/test/tests/temperature_pressure_function/gold/example_out.csv index cbdc65b63141..3c04e9431b1b 100644 --- a/modules/fluid_properties/test/tests/temperature_pressure_function/gold/example_out.csv +++ b/modules/fluid_properties/test/tests/temperature_pressure_function/gold/example_out.csv @@ -1,3 +1,3 @@ time,cp_avg,cv_diff,e_diff,h_avg,k_diff,mu_diff,rho_diff 0,0,0,0,0,0,0,0 -1,3999.9944970494,0,0,1600065.0618087,0,0,0 +1,4000.0190768956,0,0,1600065.0618087,0,0,0 diff --git a/modules/fluid_properties/unit/include/TemperaturePressureFunctionFluidPropertiesTest.h b/modules/fluid_properties/unit/include/TemperaturePressureFunctionFluidPropertiesTest.h index cd2b9fc1eb50..8811cfaf48fe 100644 --- a/modules/fluid_properties/unit/include/TemperaturePressureFunctionFluidPropertiesTest.h +++ b/modules/fluid_properties/unit/include/TemperaturePressureFunctionFluidPropertiesTest.h @@ -35,9 +35,14 @@ class TemperaturePressureFunctionFluidPropertiesTest : public MooseObjectUnitTes params.set("_object_name") = "mu_function"; _fe_problem->addFunction("ParsedFunction", "mu_function", params); + params.set("expression") = "3000. + 3*x + 5e-4*y"; + params.set("_object_name") = "cp_function"; + _fe_problem->addFunction("ParsedFunction", "cp_function", params); + _fe_problem->getFunction("k_function").initialSetup(); _fe_problem->getFunction("rho_function").initialSetup(); _fe_problem->getFunction("mu_function").initialSetup(); + _fe_problem->getFunction("cp_function").initialSetup(); buildObjects(); } @@ -55,7 +60,21 @@ class TemperaturePressureFunctionFluidPropertiesTest : public MooseObjectUnitTes _fe_problem->addUserObject("TemperaturePressureFunctionFluidProperties", "fp", uo_params); _fp = &_fe_problem->getUserObject("fp"); + + InputParameters uo_params2 = + _factory.getValidParams("TemperaturePressureFunctionFluidProperties"); + // Set the three functions and the specific isobaric heat capacity + uo_params2.set("k") = "k_function"; + uo_params2.set("rho") = "rho_function"; + uo_params2.set("mu") = "mu_function"; + uo_params2.set("cp") = "cp_function"; + uo_params2.set("T_initial_guess") = 250; + uo_params2.set("p_initial_guess") = 1e7; + + _fe_problem->addUserObject("TemperaturePressureFunctionFluidProperties", "fp_cp", uo_params2); + _fp_cp = &_fe_problem->getUserObject("fp_cp"); } TemperaturePressureFunctionFluidProperties * _fp; + TemperaturePressureFunctionFluidProperties * _fp_cp; }; diff --git a/modules/fluid_properties/unit/src/TemperaturePressureFunctionFluidPropertiesTest.C b/modules/fluid_properties/unit/src/TemperaturePressureFunctionFluidPropertiesTest.C index 9a20103139ab..a549b8f407e7 100644 --- a/modules/fluid_properties/unit/src/TemperaturePressureFunctionFluidPropertiesTest.C +++ b/modules/fluid_properties/unit/src/TemperaturePressureFunctionFluidPropertiesTest.C @@ -19,9 +19,9 @@ TEST_F(TemperaturePressureFunctionFluidPropertiesTest, fluidName) } /** - * Verify calculation of the fluid properties + * Verify calculation of the fluid properties with a constant cv */ -TEST_F(TemperaturePressureFunctionFluidPropertiesTest, properties) +TEST_F(TemperaturePressureFunctionFluidPropertiesTest, properties_constant_cv) { const Real cv = 4186.0; @@ -38,9 +38,11 @@ TEST_F(TemperaturePressureFunctionFluidPropertiesTest, properties) Real thermal_cond = 14 + 2e-2 * T + 3e-5 * p; Real visc = 1e-3 + 1e-5 * T - 3e-9 * p; Real density = 1 / v; - Real dv_dT = -1 / density / density * 2.5; - Real cp = cv + p * dv_dT; + Real alpha = -1. / density * 2.5; + Real beta = -1. / density * 32e-5; + Real cp = cv + MathUtils::pow(alpha, 2) * T / density / beta; + // Testing the properties with a constant cv ABS_TEST(_fp->cp_from_p_T(p, T), cp, tol); ABS_TEST(_fp->cv_from_p_T(p, T), cv, tol); ABS_TEST(_fp->k_from_p_T(p, T), thermal_cond, tol); @@ -73,8 +75,9 @@ TEST_F(TemperaturePressureFunctionFluidPropertiesTest, properties) thermal_cond = 14 + 2e-2 * T + 3e-5 * p; visc = 1e-3 + 1e-5 * T - 3e-9 * p; density = 1 / v; - dv_dT = -1 / density / density * 2.5; - cp = cv + p * dv_dT; + alpha = -1. / density * 2.5; + beta = -1. / density * 32e-5; + cp = cv + MathUtils::pow(alpha, 2) * T / density / beta; ABS_TEST(_fp->cp_from_p_T(p, T), cp, tol); ABS_TEST(_fp->cv_from_p_T(p, T), cv, tol); @@ -98,6 +101,94 @@ TEST_F(TemperaturePressureFunctionFluidPropertiesTest, properties) ABS_TEST(_fp->mu_from_v_e(v, e), visc, tol); } +/** + * Verify calculation of the fluid properties with a function cp + */ +TEST_F(TemperaturePressureFunctionFluidPropertiesTest, properties_function_cp) +{ + const Real tol = REL_TOL_CONSISTENCY; + // Simpson integration only does 1e-5 + const Real simpson_tol = 2e5 * tol; + // any property post using Simpson to compute 'e' + const Real large_tol = 100 * tol; + + { + Real p = 8.56E7; + Real T = 200.0; + Real cp = 3000. + 3 * T + 5e-4 * p; + Real e = 9220000.4605; + Real v = 1 / (1400 + 2.5 * T + 32e-5 * p); + Real h = e + p * v; + + // See header for expressions + Real thermal_cond = 14 + 2e-2 * T + 3e-5 * p; + Real visc = 1e-3 + 1e-5 * T - 3e-9 * p; + Real density = 1 / v; + Real alpha = -1. / density * 2.5; + Real beta = -1. / density * 32e-5; + Real cv = cp - MathUtils::pow(alpha, 2) * T / density / beta; + + ABS_TEST(_fp_cp->cp_from_p_T(p, T), cp, tol); + ABS_TEST(_fp_cp->cv_from_p_T(p, T), cv, tol); + ABS_TEST(_fp_cp->k_from_p_T(p, T), thermal_cond, tol); + ABS_TEST(_fp_cp->k_from_p_T(p, T), thermal_cond, tol); + ABS_TEST(_fp_cp->rho_from_p_T(p, T), density, tol); + ABS_TEST(_fp_cp->v_from_p_T(p, T), 1 / density, tol); + ABS_TEST(_fp_cp->e_from_p_T(p, T), e, simpson_tol); + ABS_TEST(_fp_cp->e_from_p_rho(p, 1. / v), e, simpson_tol); + ABS_TEST(_fp_cp->mu_from_p_T(p, T), visc, tol); + ABS_TEST(_fp_cp->h_from_p_T(p, T), h, simpson_tol); + ABS_TEST(_fp_cp->cp_from_v_e(v, e), cp, large_tol); + ABS_TEST(_fp_cp->cv_from_v_e(v, e), cv, large_tol); + ABS_TEST(_fp_cp->k_from_v_e(v, e), thermal_cond, large_tol); + ABS_TEST(_fp_cp->T_from_v_e(v, e), T, large_tol); + ABS_TEST(_fp_cp->T_from_p_h(p, _fp_cp->h_from_p_T(p, T)), T, large_tol); + ABS_TEST(_fp_cp->T_from_p_rho(p, 1. / v), T, large_tol); + ABS_TEST(_fp_cp->p_from_v_e(v, e), + p, + 1.5 * simpson_tol); // uses a Newton solve for variable set inversion + ABS_TEST(_fp_cp->mu_from_v_e(v, e), visc, large_tol); + } + + { + Real p = 1.06841E8; + Real T = 300.0; + Real cp = 3000. + 3 * T + 5e-4 * p; + Real e = 17061150.67487172; + Real v = 1 / (1400 + 2.5 * T + 32e-5 * p); + Real h = e + p * v; + + // See header for expressions + Real thermal_cond = 14 + 2e-2 * T + 3e-5 * p; + Real visc = 1e-3 + 1e-5 * T - 3e-9 * p; + Real density = 1 / v; + Real alpha = -1. / density * 2.5; + Real beta = -1. / density * 32e-5; + Real cv = cp - MathUtils::pow(alpha, 2) * T / density / beta; + + ABS_TEST(_fp_cp->cp_from_p_T(p, T), cp, tol); + ABS_TEST(_fp_cp->cv_from_p_T(p, T), cv, tol); + ABS_TEST(_fp_cp->k_from_p_T(p, T), thermal_cond, tol); + ABS_TEST(_fp_cp->k_from_p_T(p, T), thermal_cond, tol); + ABS_TEST(_fp_cp->rho_from_p_T(p, T), density, tol); + ABS_TEST(_fp_cp->v_from_p_T(p, T), 1 / density, tol); + ABS_TEST(_fp_cp->e_from_p_T(p, T), e, simpson_tol); + ABS_TEST(_fp_cp->e_from_p_rho(p, 1. / v), e, simpson_tol); + ABS_TEST(_fp_cp->mu_from_p_T(p, T), visc, tol); + ABS_TEST(_fp_cp->h_from_p_T(p, T), h, simpson_tol); + ABS_TEST(_fp_cp->cp_from_v_e(v, e), cp, 1.5 * large_tol); + ABS_TEST(_fp_cp->cv_from_v_e(v, e), cv, 1.5 * large_tol); + ABS_TEST(_fp_cp->k_from_v_e(v, e), thermal_cond, large_tol); + ABS_TEST(_fp_cp->T_from_v_e(v, e), T, 1.5 * large_tol); + ABS_TEST(_fp_cp->T_from_p_h(p, _fp_cp->h_from_p_T(p, T)), T, large_tol); + ABS_TEST(_fp_cp->T_from_p_rho(p, 1. / v), T, large_tol); + ABS_TEST(_fp_cp->p_from_v_e(v, e), + p, + 10 * simpson_tol); // uses a Newton solve for variable set inversion + ABS_TEST(_fp_cp->mu_from_v_e(v, e), visc, large_tol); + } +} + /** * Verify calculation of the derivatives by comparing with finite * differences @@ -105,9 +196,10 @@ TEST_F(TemperaturePressureFunctionFluidPropertiesTest, properties) TEST_F(TemperaturePressureFunctionFluidPropertiesTest, derivatives) { const Real tol = REL_TOL_DERIVATIVE; - const Real large_tol = 30 * tol; + // Finite difference just does not get much + const Real large_tol = 1000 * tol; - // p, T and v,e are not consistent here but that's ok + // p, T and v, e are not consistent here but that's ok Real p = 1.0E7; Real T = 10.0; Real e = 8.372E5; @@ -120,13 +212,28 @@ TEST_F(TemperaturePressureFunctionFluidPropertiesTest, derivatives) DERIV_TEST(_fp->e_from_p_T, p, T, tol); DERIV_TEST(_fp->h_from_p_T, p, T, tol); DERIV_TEST(_fp->k_from_p_T, p, T, tol); - DERIV_TEST(_fp->cp_from_p_T, p, T, large_tol); // uses finite differencing + DERIV_TEST(_fp->cp_from_p_T, p, T, 4 * large_tol); // uses finite differencing DERIV_TEST(_fp->cv_from_p_T, p, T, tol); DERIV_TEST(_fp->cp_from_v_e, v, e, large_tol); // uses finite differencing DERIV_TEST(_fp->cv_from_v_e, v, e, tol); + e = 80150; + v = 2.16216E-4; + + AD_DERIV_TEST(_fp_cp->rho_from_p_T, p, T, tol); + DERIV_TEST(_fp_cp->rho_from_p_T, p, T, tol); + DERIV_TEST(_fp_cp->v_from_p_T, p, T, tol); + DERIV_TEST(_fp_cp->mu_from_p_T, p, T, tol); + DERIV_TEST(_fp_cp->e_from_p_T, p, T, tol); + DERIV_TEST(_fp_cp->h_from_p_T, p, T, tol); + DERIV_TEST(_fp_cp->k_from_p_T, p, T, tol); + DERIV_TEST(_fp_cp->cp_from_p_T, p, T, large_tol); // uses finite differencing + DERIV_TEST(_fp_cp->cv_from_p_T, p, T, large_tol); + DERIV_TEST(_fp_cp->cp_from_v_e, v, e, large_tol); // uses finite differencing + DERIV_TEST(_fp_cp->cv_from_v_e, v, e, large_tol); + p = 5.0E7; - T = 90.0; + T = 190.0; e = 1.6744E6; v = 6.25E-4; @@ -137,10 +244,25 @@ TEST_F(TemperaturePressureFunctionFluidPropertiesTest, derivatives) DERIV_TEST(_fp->e_from_p_T, p, T, tol); DERIV_TEST(_fp->h_from_p_T, p, T, tol); DERIV_TEST(_fp->k_from_p_T, p, T, tol); - DERIV_TEST(_fp->cp_from_p_T, p, T, large_tol); // uses finite differencing + DERIV_TEST(_fp->cp_from_p_T, p, T, 2 * large_tol); // uses finite differencing DERIV_TEST(_fp->cv_from_p_T, p, T, tol); - DERIV_TEST(_fp->cp_from_v_e, v, e, large_tol); // uses finite differencing + DERIV_TEST(_fp->cp_from_v_e, v, e, 2 * large_tol); // uses finite differencing DERIV_TEST(_fp->cv_from_v_e, v, e, tol); + + e = 5.37415e+06; + v = 5.59441e-05; + + AD_DERIV_TEST(_fp_cp->rho_from_p_T, p, T, tol); + DERIV_TEST(_fp_cp->rho_from_p_T, p, T, tol); + DERIV_TEST(_fp_cp->v_from_p_T, p, T, tol); + DERIV_TEST(_fp_cp->mu_from_p_T, p, T, tol); + DERIV_TEST(_fp_cp->e_from_p_T, p, T, tol); + DERIV_TEST(_fp_cp->h_from_p_T, p, T, tol); + DERIV_TEST(_fp_cp->k_from_p_T, p, T, tol); + DERIV_TEST(_fp_cp->cp_from_p_T, p, T, large_tol); // uses finite differencing + DERIV_TEST(_fp_cp->cv_from_p_T, p, T, large_tol); + DERIV_TEST(_fp_cp->cp_from_v_e, v, e, large_tol); // uses finite differencing + DERIV_TEST(_fp_cp->cv_from_v_e, v, e, large_tol); } /** From 25194473a389aeacecd95cf5eabe1017b9bc8db6 Mon Sep 17 00:00:00 2001 From: Guillaume Giudicelli Date: Tue, 7 Jan 2025 01:32:25 +0100 Subject: [PATCH 3/6] Address Josh's review - better doco - fix beta (different notation than wikipedia) and cp-cv conversions - use a parameterized integral for the e(T) = int(cv(T)) approximation - adapt documentation - adapt unit test to new values from numerical integration - clean up tolerances in unit test + document remaining approximation in computation of e Co-authored-by: joshuahansel --- ...peraturePressureFunctionFluidProperties.md | 17 +++--- ...mperaturePressureFunctionFluidProperties.h | 2 + ...mperaturePressureFunctionFluidProperties.C | 61 +++++++++---------- ...aturePressureFunctionFluidPropertiesTest.h | 3 +- ...aturePressureFunctionFluidPropertiesTest.C | 40 ++++++------ 5 files changed, 62 insertions(+), 61 deletions(-) diff --git a/modules/fluid_properties/doc/content/source/fluidproperties/TemperaturePressureFunctionFluidProperties.md b/modules/fluid_properties/doc/content/source/fluidproperties/TemperaturePressureFunctionFluidProperties.md index 9ee80b27c3fb..7b088d146b19 100644 --- a/modules/fluid_properties/doc/content/source/fluidproperties/TemperaturePressureFunctionFluidProperties.md +++ b/modules/fluid_properties/doc/content/source/fluidproperties/TemperaturePressureFunctionFluidProperties.md @@ -9,12 +9,14 @@ functions. The following relations hold true: \begin{aligned} x = T \\ y = P \\ - \rho = \rho^{user}(time=0, (x=T, y=P, z=0)) \\ - \mu = \mu^{user}(time=0, (x=T, y=P, z=0)) \\ - k = k^{user}(time=0, (x=T, y=P, z=0)) \\ + \rho = \rho^{user}(t=0, (x=T, y=P, z=0)) \\ + \mu = \mu^{user}(t=0, (x=T, y=P, z=0)) \\ + k = k^{user}(t=0, (x=T, y=P, z=0)) \\ \end{aligned} \end{equation} +Both the time (`t`) and Z-axis dimension are not used here. A fluid property made to depend on time will +not be properly updated by this `FluidProperties` object. There are two options for specific heat. Either the user sets a constant specific isochoric heat capacity \begin{equation} @@ -30,12 +32,12 @@ Or, the user uses a function of temperature and pressure (same arguments as for \begin{aligned} x = T \\ y = P \\ - cp = cp^{user}(time=0, (x=T, y=P, z=0)) + cp = cp^{user}(t=0, (x=T, y=P, z=0)) cv = cp - \dfrac{\alpha ^ 2 T}{\rho \beta_T} \end{aligned} \end{equation} -with $T$ the temperature, $P$ the pressure, $\rho$ the density, $\mu$ the dynamic viscosity, $k$ the thermal conductivity, $c_v$ the specific isochoric heat capacity, $e$ the specific energy, $\alpha$ the coefficient of thermal expansion, $\beta_T$ the isothermal +with $T$ the temperature, $P$ the pressure, $\rho$ the density, $\mu$ the dynamic viscosity, $k$ the thermal conductivity, $c_v$ the specific isochoric heat capacity, $e$ the specific internal energy, $T_{ref}$ a reference temperature at which the specific internal energy is equal to a reference energy $e_{ref}$, $\alpha$ the coefficient of thermal expansion, $\beta_T$ the isothermal compressibility, and the $user$ exponent indicating a user-passed parameter. The derivatives of the fluid properties are obtained using the `Function`(s) gradient components @@ -51,8 +53,9 @@ partial. Notable missing implementations are routines for entropy, the speed of conversions between specific enthalpy and specific energy. !alert note -When using a function for the isobaric specific heat capacity, Simpson's rule is used to compute -the specific energy from the isochoric specific heat capacity. This limits the accuracy of the calculation. +When using a function for the isobaric specific heat capacity, a numerical integration is performed to compute +$e(p,T)$ as $e_{ref} + \int_{T_{ref}}^T c_v(p,T) dT$. Note that this neglects the $dV$ term. This is exact +for incompressible fluids and ideal gases. ## Example Input File Syntax diff --git a/modules/fluid_properties/include/fluidproperties/TemperaturePressureFunctionFluidProperties.h b/modules/fluid_properties/include/fluidproperties/TemperaturePressureFunctionFluidProperties.h index c9343be517b4..4ae185145431 100644 --- a/modules/fluid_properties/include/fluidproperties/TemperaturePressureFunctionFluidProperties.h +++ b/modules/fluid_properties/include/fluidproperties/TemperaturePressureFunctionFluidProperties.h @@ -359,4 +359,6 @@ class TemperaturePressureFunctionFluidProperties : public SinglePhaseFluidProper const Real _e_ref; /// Reference temperature for the reference specific energy const Real _T_ref; + /// Number of steps to take when integrating the specific heat to compute the specific energy + const unsigned int _n_integration; }; diff --git a/modules/fluid_properties/src/fluidproperties/TemperaturePressureFunctionFluidProperties.C b/modules/fluid_properties/src/fluidproperties/TemperaturePressureFunctionFluidProperties.C index d4146b96d166..c856d2772829 100644 --- a/modules/fluid_properties/src/fluidproperties/TemperaturePressureFunctionFluidProperties.C +++ b/modules/fluid_properties/src/fluidproperties/TemperaturePressureFunctionFluidProperties.C @@ -27,8 +27,12 @@ TemperaturePressureFunctionFluidProperties::validParams() "cp", "Isobaric specific heat function of temperature and pressure [J/(kg-K)]"); params.addRangeCheckedParam( "cv", 0, "cv >= 0", "Constant isochoric specific heat [J/(kg-K)]"); - params.addParam("e_ref", 0, "Specific energy at the reference temperature"); - params.addParam("T_ref", 0, "Reference temperature for the specific energy"); + params.addParam("e_ref", 0, "Specific internal energy at the reference temperature"); + params.addParam("T_ref", 0, "Reference temperature for the specific internal energy"); + params.addParam( + "n_integration_intervals", + 10, + "Number of intervals for integrating cv(T) to compute e(T) from e(T_ref)"); params.addClassDescription( "Single-phase fluid properties that allows to provide thermal " @@ -43,10 +47,11 @@ TemperaturePressureFunctionFluidProperties::TemperaturePressureFunctionFluidProp _cv(getParam("cv")), _cv_is_constant(_cv != 0), _e_ref(getParam("e_ref")), - _T_ref(getParam("T_ref")) + _T_ref(getParam("T_ref")), + _n_integration(getParam("n_integration_intervals")) { if (isParamValid("cp") && _cv_is_constant) - paramError("cp", "Can only set a function for cp if a constant is not specified for cv."); + paramError("cp", "The parameter 'cp' may only be specified if 'cv' is unspecified or is zero."); } void @@ -350,7 +355,7 @@ TemperaturePressureFunctionFluidProperties::h_from_p_T( Real rho, drho_dp, drho_dT; rho_from_p_T(pressure, temperature, rho, drho_dp, drho_dT); h = e + pressure / rho; - dh_dp = de_dp + 1 / rho - pressure / rho / rho * drho_dp; + dh_dp = de_dp + 1. / rho - pressure / rho / rho * drho_dp; dh_dT = de_dT - pressure * drho_dT / rho / rho; } @@ -361,12 +366,14 @@ TemperaturePressureFunctionFluidProperties::e_from_p_T(Real pressure, Real tempe return _e_ref + _cv * (temperature - _T_ref); else { - // Use Simpson's 3/8 rule to integrate cv - Real h = (temperature - _T_ref) / 3; - Real integral = - 3. / 8. * h * - (cv_from_p_T(pressure, _T_ref) + 3 * cv_from_p_T(pressure, _T_ref + h) + - 3 * cv_from_p_T(pressure, _T_ref + 2 * h) + cv_from_p_T(pressure, temperature)); + const Real h = (temperature - _T_ref) / _n_integration; + Real integral = 0; + // Centered step integration is second-order + for (const auto i : make_range(_n_integration)) + integral += cv_from_p_T(pressure, _T_ref + (i + 0.5) * h); + integral *= h; + // we are still missing the dV or dP term to go from V/P_ref (e_ref = e(T_ref, V/P_ref)) + // to current V. The dT term is the largest one though return _e_ref + integral; } } @@ -384,20 +391,9 @@ TemperaturePressureFunctionFluidProperties::e_from_p_T( else { e = e_from_p_T(pressure, temperature); - // Use Simpson's 3/8 rule to integrate cv then take the derivative - Real h = (temperature - _T_ref) / 3; - Real cv0, dcv_dp0, dcv_dT0; - cv_from_p_T(pressure, _T_ref, cv0, dcv_dp0, dcv_dT0); - Real cv1, dcv_dp1, dcv_dT1; - cv_from_p_T(pressure, _T_ref + h, cv1, dcv_dp1, dcv_dT1); - Real cv2, dcv_dp2, dcv_dT2; - cv_from_p_T(pressure, _T_ref + 2 * h, cv2, dcv_dp2, dcv_dT2); - Real cv3, dcv_dp3, dcv_dT3; - cv_from_p_T(pressure, temperature, cv3, dcv_dp3, dcv_dT3); - de_dp = 3. / 8 * h * (dcv_dp0 + 3 * dcv_dp1 + 3 * dcv_dp2 + dcv_dp3); - // Consistency with how the energy is computed is more important here - de_dT = - 1. / 8. * (cv0 + 3 * cv1 + 3 * cv2 + cv3) + 3. / 8 * h * (dcv_dT1 + 2 * dcv_dT2 + dcv_dT3); + Real ep = e_from_p_T(pressure * (1 + 1e-8), temperature); + de_dp = (ep - e) / (pressure * 1e-8); + de_dT = cv_from_p_T(pressure, temperature); } } @@ -421,7 +417,7 @@ TemperaturePressureFunctionFluidProperties::beta_from_p_T(Real pressure, Real te { Real rho, drho_dp, drho_dT; rho_from_p_T(pressure, temperature, rho, drho_dp, drho_dT); - return -drho_dp / rho; + return -drho_dT / rho; } Real @@ -431,8 +427,10 @@ TemperaturePressureFunctionFluidProperties::cp_from_p_T(Real p, Real T) const { Real rho, drho_dp, drho_dT; rho_from_p_T(p, T, rho, drho_dp, drho_dT); - Real alpha = -1. / rho * drho_dT; - return _cv + MathUtils::pow(alpha, 2) * T / rho / beta_from_p_T(p, T); + // Wikipedia notation for thermal expansion / compressibility coefficients + Real alpha = -drho_dT / rho; + Real beta = -drho_dp / rho; + return _cv + MathUtils::pow(alpha, 2) * T / rho / beta; } else { @@ -474,9 +472,10 @@ TemperaturePressureFunctionFluidProperties::cv_from_p_T(Real pressure, Real temp { Real rho, drho_dp, drho_dT; rho_from_p_T(pressure, temperature, rho, drho_dp, drho_dT); - Real alpha = -1. / rho * drho_dT; - return cp_from_p_T(pressure, temperature) - - MathUtils::pow(alpha, 2) * temperature / rho / beta_from_p_T(pressure, temperature); + // Wikipedia notation for thermal expansion / compressibility coefficients + Real alpha = -drho_dT / rho; + Real beta = -drho_dp / rho; + return cp_from_p_T(pressure, temperature) - MathUtils::pow(alpha, 2) * temperature / rho / beta; } } diff --git a/modules/fluid_properties/unit/include/TemperaturePressureFunctionFluidPropertiesTest.h b/modules/fluid_properties/unit/include/TemperaturePressureFunctionFluidPropertiesTest.h index 8811cfaf48fe..c14cbb25aa19 100644 --- a/modules/fluid_properties/unit/include/TemperaturePressureFunctionFluidPropertiesTest.h +++ b/modules/fluid_properties/unit/include/TemperaturePressureFunctionFluidPropertiesTest.h @@ -63,13 +63,14 @@ class TemperaturePressureFunctionFluidPropertiesTest : public MooseObjectUnitTes InputParameters uo_params2 = _factory.getValidParams("TemperaturePressureFunctionFluidProperties"); - // Set the three functions and the specific isobaric heat capacity + // Set the three functions and the specific heat capacity uo_params2.set("k") = "k_function"; uo_params2.set("rho") = "rho_function"; uo_params2.set("mu") = "mu_function"; uo_params2.set("cp") = "cp_function"; uo_params2.set("T_initial_guess") = 250; uo_params2.set("p_initial_guess") = 1e7; + uo_params2.set("n_integration_intervals") = 2000; _fe_problem->addUserObject("TemperaturePressureFunctionFluidProperties", "fp_cp", uo_params2); _fp_cp = &_fe_problem->getUserObject("fp_cp"); diff --git a/modules/fluid_properties/unit/src/TemperaturePressureFunctionFluidPropertiesTest.C b/modules/fluid_properties/unit/src/TemperaturePressureFunctionFluidPropertiesTest.C index a549b8f407e7..ef992e470c62 100644 --- a/modules/fluid_properties/unit/src/TemperaturePressureFunctionFluidPropertiesTest.C +++ b/modules/fluid_properties/unit/src/TemperaturePressureFunctionFluidPropertiesTest.C @@ -107,16 +107,16 @@ TEST_F(TemperaturePressureFunctionFluidPropertiesTest, properties_constant_cv) TEST_F(TemperaturePressureFunctionFluidPropertiesTest, properties_function_cp) { const Real tol = REL_TOL_CONSISTENCY; - // Simpson integration only does 1e-5 - const Real simpson_tol = 2e5 * tol; - // any property post using Simpson to compute 'e' - const Real large_tol = 100 * tol; + // properties computed using numerical integral evaluation + const Real large_tol = 1e-8; + // properties using both numerical integration and Newton's method + const Real newton_tol = 1e-4; { Real p = 8.56E7; Real T = 200.0; Real cp = 3000. + 3 * T + 5e-4 * p; - Real e = 9220000.4605; + Real e = 9220000.46051058359; Real v = 1 / (1400 + 2.5 * T + 32e-5 * p); Real h = e + p * v; @@ -134,19 +134,17 @@ TEST_F(TemperaturePressureFunctionFluidPropertiesTest, properties_function_cp) ABS_TEST(_fp_cp->k_from_p_T(p, T), thermal_cond, tol); ABS_TEST(_fp_cp->rho_from_p_T(p, T), density, tol); ABS_TEST(_fp_cp->v_from_p_T(p, T), 1 / density, tol); - ABS_TEST(_fp_cp->e_from_p_T(p, T), e, simpson_tol); - ABS_TEST(_fp_cp->e_from_p_rho(p, 1. / v), e, simpson_tol); + ABS_TEST(_fp_cp->e_from_p_T(p, T), e, tol); + ABS_TEST(_fp_cp->e_from_p_rho(p, 1. / v), e, 10 * large_tol); ABS_TEST(_fp_cp->mu_from_p_T(p, T), visc, tol); - ABS_TEST(_fp_cp->h_from_p_T(p, T), h, simpson_tol); + ABS_TEST(_fp_cp->h_from_p_T(p, T), h, tol); ABS_TEST(_fp_cp->cp_from_v_e(v, e), cp, large_tol); ABS_TEST(_fp_cp->cv_from_v_e(v, e), cv, large_tol); ABS_TEST(_fp_cp->k_from_v_e(v, e), thermal_cond, large_tol); ABS_TEST(_fp_cp->T_from_v_e(v, e), T, large_tol); ABS_TEST(_fp_cp->T_from_p_h(p, _fp_cp->h_from_p_T(p, T)), T, large_tol); ABS_TEST(_fp_cp->T_from_p_rho(p, 1. / v), T, large_tol); - ABS_TEST(_fp_cp->p_from_v_e(v, e), - p, - 1.5 * simpson_tol); // uses a Newton solve for variable set inversion + ABS_TEST(_fp_cp->p_from_v_e(v, e), p, newton_tol); ABS_TEST(_fp_cp->mu_from_v_e(v, e), visc, large_tol); } @@ -154,7 +152,7 @@ TEST_F(TemperaturePressureFunctionFluidPropertiesTest, properties_function_cp) Real p = 1.06841E8; Real T = 300.0; Real cp = 3000. + 3 * T + 5e-4 * p; - Real e = 17061150.67487172; + Real e = 17061150.6748719364; Real v = 1 / (1400 + 2.5 * T + 32e-5 * p); Real h = e + p * v; @@ -172,19 +170,17 @@ TEST_F(TemperaturePressureFunctionFluidPropertiesTest, properties_function_cp) ABS_TEST(_fp_cp->k_from_p_T(p, T), thermal_cond, tol); ABS_TEST(_fp_cp->rho_from_p_T(p, T), density, tol); ABS_TEST(_fp_cp->v_from_p_T(p, T), 1 / density, tol); - ABS_TEST(_fp_cp->e_from_p_T(p, T), e, simpson_tol); - ABS_TEST(_fp_cp->e_from_p_rho(p, 1. / v), e, simpson_tol); + ABS_TEST(_fp_cp->e_from_p_T(p, T), e, tol); + ABS_TEST(_fp_cp->e_from_p_rho(p, 1. / v), e, tol); ABS_TEST(_fp_cp->mu_from_p_T(p, T), visc, tol); - ABS_TEST(_fp_cp->h_from_p_T(p, T), h, simpson_tol); + ABS_TEST(_fp_cp->h_from_p_T(p, T), h, tol); ABS_TEST(_fp_cp->cp_from_v_e(v, e), cp, 1.5 * large_tol); ABS_TEST(_fp_cp->cv_from_v_e(v, e), cv, 1.5 * large_tol); ABS_TEST(_fp_cp->k_from_v_e(v, e), thermal_cond, large_tol); ABS_TEST(_fp_cp->T_from_v_e(v, e), T, 1.5 * large_tol); ABS_TEST(_fp_cp->T_from_p_h(p, _fp_cp->h_from_p_T(p, T)), T, large_tol); ABS_TEST(_fp_cp->T_from_p_rho(p, 1. / v), T, large_tol); - ABS_TEST(_fp_cp->p_from_v_e(v, e), - p, - 10 * simpson_tol); // uses a Newton solve for variable set inversion + ABS_TEST(_fp_cp->p_from_v_e(v, e), p, newton_tol); ABS_TEST(_fp_cp->mu_from_v_e(v, e), visc, large_tol); } } @@ -217,8 +213,8 @@ TEST_F(TemperaturePressureFunctionFluidPropertiesTest, derivatives) DERIV_TEST(_fp->cp_from_v_e, v, e, large_tol); // uses finite differencing DERIV_TEST(_fp->cv_from_v_e, v, e, tol); - e = 80150; - v = 2.16216E-4; + e = 80150.045818949831; + v = 0.00021621621621621621; AD_DERIV_TEST(_fp_cp->rho_from_p_T, p, T, tol); DERIV_TEST(_fp_cp->rho_from_p_T, p, T, tol); @@ -249,8 +245,8 @@ TEST_F(TemperaturePressureFunctionFluidPropertiesTest, derivatives) DERIV_TEST(_fp->cp_from_v_e, v, e, 2 * large_tol); // uses finite differencing DERIV_TEST(_fp->cv_from_v_e, v, e, tol); - e = 5.37415e+06; - v = 5.59441e-05; + e = 5374151.1232993276; + v = 5.5944055944055945e-05; AD_DERIV_TEST(_fp_cp->rho_from_p_T, p, T, tol); DERIV_TEST(_fp_cp->rho_from_p_T, p, T, tol); From 68479ec37a63621a3b9955b095d15147af12eeaf Mon Sep 17 00:00:00 2001 From: Guillaume Giudicelli Date: Sun, 12 Jan 2025 15:57:30 -0700 Subject: [PATCH 4/6] Address FP unit test failures on M1 Mac --- .../src/fluidproperties/SimpleFluidProperties.C | 13 ++++--------- .../unit/src/BrineFluidPropertiesTest.C | 4 ++-- 2 files changed, 6 insertions(+), 11 deletions(-) diff --git a/modules/fluid_properties/src/fluidproperties/SimpleFluidProperties.C b/modules/fluid_properties/src/fluidproperties/SimpleFluidProperties.C index a6697617c020..66ad87b448c9 100644 --- a/modules/fluid_properties/src/fluidproperties/SimpleFluidProperties.C +++ b/modules/fluid_properties/src/fluidproperties/SimpleFluidProperties.C @@ -163,20 +163,15 @@ SimpleFluidProperties::c_from_p_T( } Real -SimpleFluidProperties::c_from_v_e(Real v, Real e) const +SimpleFluidProperties::c_from_v_e(Real v, Real /*e*/) const { - Real T = T_from_v_e(v, e); - Real p = p_from_v_e(v, e); - return std::sqrt(_bulk_modulus / rho_from_p_T(p, T)); + return std::sqrt(_bulk_modulus * v); } void -SimpleFluidProperties::c_from_v_e(Real v, Real e, Real & c, Real & dc_dv, Real & dc_de) const +SimpleFluidProperties::c_from_v_e(Real v, Real /*e*/, Real & c, Real & dc_dv, Real & dc_de) const { - Real T = T_from_v_e(v, e); - Real p = p_from_v_e(v, e); - - c = std::sqrt(_bulk_modulus / rho_from_p_T(p, T)); + c = std::sqrt(_bulk_modulus * v); dc_dv = 0.5 * std::sqrt(_bulk_modulus / v); dc_de = 0.0; diff --git a/modules/fluid_properties/unit/src/BrineFluidPropertiesTest.C b/modules/fluid_properties/unit/src/BrineFluidPropertiesTest.C index 8d68ea1a3aec..50cb3cf2a0c5 100644 --- a/modules/fluid_properties/unit/src/BrineFluidPropertiesTest.C +++ b/modules/fluid_properties/unit/src/BrineFluidPropertiesTest.C @@ -170,7 +170,7 @@ TEST_F(BrineFluidPropertiesTest, derivatives) Real h = 0.0, dh_dp = 0.0, dh_dT = 0.0, dh_dx = 0.0; _fp->h_from_p_T_X(p, T, x, h, dh_dp, dh_dT, dh_dx); - ABS_TEST(h, _fp->h_from_p_T_X(p, T, x), REL_TOL_CONSISTENCY); + ABS_TEST(h, _fp->h_from_p_T_X(p, T, x), 1.2 * REL_TOL_CONSISTENCY); REL_TEST(dh_dp, dh_dp_fd, 1.0e-4); REL_TEST(dh_dT, dh_dT_fd, 1.0e-6); REL_TEST(dh_dx, dh_dx_fd, 1.0e-6); @@ -183,7 +183,7 @@ TEST_F(BrineFluidPropertiesTest, derivatives) Real e = 0.0, de_dp = 0.0, de_dT = 0.0, de_dx = 0.0; _fp->e_from_p_T_X(p, T, x, e, de_dp, de_dT, de_dx); - ABS_TEST(e, _fp->e_from_p_T_X(p, T, x), REL_TOL_CONSISTENCY); + ABS_TEST(e, _fp->e_from_p_T_X(p, T, x), 1.2 * REL_TOL_CONSISTENCY); REL_TEST(de_dp, de_dp_fd, 1.0e-3); REL_TEST(de_dT, de_dT_fd, 1.0e-6); REL_TEST(de_dx, de_dx_fd, 1.0e-6); From 829d04ffbc1c7c1408b4632de4a9c7f1a6cd6cdf Mon Sep 17 00:00:00 2001 From: Guillaume Giudicelli Date: Mon, 13 Jan 2025 12:27:56 -0700 Subject: [PATCH 5/6] Address Josh's review - doco improvement - add T to integration param name --- .../TemperaturePressureFunctionFluidProperties.md | 14 +++++++------- .../TemperaturePressureFunctionFluidProperties.h | 2 +- .../TemperaturePressureFunctionFluidProperties.C | 8 ++++---- ...emperaturePressureFunctionFluidPropertiesTest.h | 2 +- 4 files changed, 13 insertions(+), 13 deletions(-) diff --git a/modules/fluid_properties/doc/content/source/fluidproperties/TemperaturePressureFunctionFluidProperties.md b/modules/fluid_properties/doc/content/source/fluidproperties/TemperaturePressureFunctionFluidProperties.md index 7b088d146b19..b27db4546e4e 100644 --- a/modules/fluid_properties/doc/content/source/fluidproperties/TemperaturePressureFunctionFluidProperties.md +++ b/modules/fluid_properties/doc/content/source/fluidproperties/TemperaturePressureFunctionFluidProperties.md @@ -9,9 +9,9 @@ functions. The following relations hold true: \begin{aligned} x = T \\ y = P \\ - \rho = \rho^{user}(t=0, (x=T, y=P, z=0)) \\ - \mu = \mu^{user}(t=0, (x=T, y=P, z=0)) \\ - k = k^{user}(t=0, (x=T, y=P, z=0)) \\ + \rho = \rho^{u}(t=0, (x=T, y=P, z=0)) \\ + \mu = \mu^{u}(t=0, (x=T, y=P, z=0)) \\ + k = k^{u}(t=0, (x=T, y=P, z=0)) \\ \end{aligned} \end{equation} @@ -21,7 +21,7 @@ There are two options for specific heat. Either the user sets a constant specifi \begin{equation} \begin{aligned} - c_v = c_v^{user} \\ + c_v = c_v^{u} \\ e = e_{ref} + c_v * (T - T_{ref}) \end{aligned} \end{equation} @@ -32,13 +32,13 @@ Or, the user uses a function of temperature and pressure (same arguments as for \begin{aligned} x = T \\ y = P \\ - cp = cp^{user}(t=0, (x=T, y=P, z=0)) - cv = cp - \dfrac{\alpha ^ 2 T}{\rho \beta_T} + cp = cp^{u}(t=0, (x=T, y=P, z=0)) + cv = cp - \dfrac{\alpha^2 T}{\rho \beta_T} \end{aligned} \end{equation} with $T$ the temperature, $P$ the pressure, $\rho$ the density, $\mu$ the dynamic viscosity, $k$ the thermal conductivity, $c_v$ the specific isochoric heat capacity, $e$ the specific internal energy, $T_{ref}$ a reference temperature at which the specific internal energy is equal to a reference energy $e_{ref}$, $\alpha$ the coefficient of thermal expansion, $\beta_T$ the isothermal -compressibility, and the $user$ exponent indicating a user-passed parameter. +compressibility, and the $^u$ exponent indicating a user-passed parameter. The derivatives of the fluid properties are obtained using the `Function`(s) gradient components and the appropriate derivative chaining for derived properties. diff --git a/modules/fluid_properties/include/fluidproperties/TemperaturePressureFunctionFluidProperties.h b/modules/fluid_properties/include/fluidproperties/TemperaturePressureFunctionFluidProperties.h index 4ae185145431..f34db5a8cc22 100644 --- a/modules/fluid_properties/include/fluidproperties/TemperaturePressureFunctionFluidProperties.h +++ b/modules/fluid_properties/include/fluidproperties/TemperaturePressureFunctionFluidProperties.h @@ -360,5 +360,5 @@ class TemperaturePressureFunctionFluidProperties : public SinglePhaseFluidProper /// Reference temperature for the reference specific energy const Real _T_ref; /// Number of steps to take when integrating the specific heat to compute the specific energy - const unsigned int _n_integration; + const unsigned int _n_integration_dT; }; diff --git a/modules/fluid_properties/src/fluidproperties/TemperaturePressureFunctionFluidProperties.C b/modules/fluid_properties/src/fluidproperties/TemperaturePressureFunctionFluidProperties.C index c856d2772829..46f7290adc0d 100644 --- a/modules/fluid_properties/src/fluidproperties/TemperaturePressureFunctionFluidProperties.C +++ b/modules/fluid_properties/src/fluidproperties/TemperaturePressureFunctionFluidProperties.C @@ -30,7 +30,7 @@ TemperaturePressureFunctionFluidProperties::validParams() params.addParam("e_ref", 0, "Specific internal energy at the reference temperature"); params.addParam("T_ref", 0, "Reference temperature for the specific internal energy"); params.addParam( - "n_integration_intervals", + "dT_integration_intervals", 10, "Number of intervals for integrating cv(T) to compute e(T) from e(T_ref)"); @@ -48,7 +48,7 @@ TemperaturePressureFunctionFluidProperties::TemperaturePressureFunctionFluidProp _cv_is_constant(_cv != 0), _e_ref(getParam("e_ref")), _T_ref(getParam("T_ref")), - _n_integration(getParam("n_integration_intervals")) + _n_integration_dT(getParam("dT_integration_intervals")) { if (isParamValid("cp") && _cv_is_constant) paramError("cp", "The parameter 'cp' may only be specified if 'cv' is unspecified or is zero."); @@ -366,10 +366,10 @@ TemperaturePressureFunctionFluidProperties::e_from_p_T(Real pressure, Real tempe return _e_ref + _cv * (temperature - _T_ref); else { - const Real h = (temperature - _T_ref) / _n_integration; + const Real h = (temperature - _T_ref) / _n_integration_dT; Real integral = 0; // Centered step integration is second-order - for (const auto i : make_range(_n_integration)) + for (const auto i : make_range(_n_integration_dT)) integral += cv_from_p_T(pressure, _T_ref + (i + 0.5) * h); integral *= h; // we are still missing the dV or dP term to go from V/P_ref (e_ref = e(T_ref, V/P_ref)) diff --git a/modules/fluid_properties/unit/include/TemperaturePressureFunctionFluidPropertiesTest.h b/modules/fluid_properties/unit/include/TemperaturePressureFunctionFluidPropertiesTest.h index c14cbb25aa19..e49e65776cf3 100644 --- a/modules/fluid_properties/unit/include/TemperaturePressureFunctionFluidPropertiesTest.h +++ b/modules/fluid_properties/unit/include/TemperaturePressureFunctionFluidPropertiesTest.h @@ -70,7 +70,7 @@ class TemperaturePressureFunctionFluidPropertiesTest : public MooseObjectUnitTes uo_params2.set("cp") = "cp_function"; uo_params2.set("T_initial_guess") = 250; uo_params2.set("p_initial_guess") = 1e7; - uo_params2.set("n_integration_intervals") = 2000; + uo_params2.set("dT_integration_intervals") = 2000; _fe_problem->addUserObject("TemperaturePressureFunctionFluidProperties", "fp_cp", uo_params2); _fp_cp = &_fe_problem->getUserObject("fp_cp"); From 683486e3307baee9275fb51f28b201040eecefe2 Mon Sep 17 00:00:00 2001 From: Guillaume Giudicelli Date: Wed, 22 Jan 2025 08:53:44 -0700 Subject: [PATCH 6/6] Switch to an interval size based precision parameter for cv(T) integration to compute e --- .../TemperaturePressureFunctionFluidProperties.h | 4 ++-- .../TemperaturePressureFunctionFluidProperties.C | 14 +++++++------- ...emperaturePressureFunctionFluidPropertiesTest.h | 2 +- ...emperaturePressureFunctionFluidPropertiesTest.C | 8 ++++---- 4 files changed, 14 insertions(+), 14 deletions(-) diff --git a/modules/fluid_properties/include/fluidproperties/TemperaturePressureFunctionFluidProperties.h b/modules/fluid_properties/include/fluidproperties/TemperaturePressureFunctionFluidProperties.h index f34db5a8cc22..a69b615f326f 100644 --- a/modules/fluid_properties/include/fluidproperties/TemperaturePressureFunctionFluidProperties.h +++ b/modules/fluid_properties/include/fluidproperties/TemperaturePressureFunctionFluidProperties.h @@ -359,6 +359,6 @@ class TemperaturePressureFunctionFluidProperties : public SinglePhaseFluidProper const Real _e_ref; /// Reference temperature for the reference specific energy const Real _T_ref; - /// Number of steps to take when integrating the specific heat to compute the specific energy - const unsigned int _n_integration_dT; + /// Size of temperature intervals when integrating the specific heat to compute the specific energy + const Real _integration_dT; }; diff --git a/modules/fluid_properties/src/fluidproperties/TemperaturePressureFunctionFluidProperties.C b/modules/fluid_properties/src/fluidproperties/TemperaturePressureFunctionFluidProperties.C index 46f7290adc0d..a91db7e2252c 100644 --- a/modules/fluid_properties/src/fluidproperties/TemperaturePressureFunctionFluidProperties.C +++ b/modules/fluid_properties/src/fluidproperties/TemperaturePressureFunctionFluidProperties.C @@ -29,10 +29,9 @@ TemperaturePressureFunctionFluidProperties::validParams() "cv", 0, "cv >= 0", "Constant isochoric specific heat [J/(kg-K)]"); params.addParam("e_ref", 0, "Specific internal energy at the reference temperature"); params.addParam("T_ref", 0, "Reference temperature for the specific internal energy"); - params.addParam( - "dT_integration_intervals", - 10, - "Number of intervals for integrating cv(T) to compute e(T) from e(T_ref)"); + params.addParam("dT_integration_intervals", + 10, + "Size of intervals for integrating cv(T) to compute e(T) from e(T_ref)"); params.addClassDescription( "Single-phase fluid properties that allows to provide thermal " @@ -48,7 +47,7 @@ TemperaturePressureFunctionFluidProperties::TemperaturePressureFunctionFluidProp _cv_is_constant(_cv != 0), _e_ref(getParam("e_ref")), _T_ref(getParam("T_ref")), - _n_integration_dT(getParam("dT_integration_intervals")) + _integration_dT(getParam("dT_integration_intervals")) { if (isParamValid("cp") && _cv_is_constant) paramError("cp", "The parameter 'cp' may only be specified if 'cv' is unspecified or is zero."); @@ -366,10 +365,11 @@ TemperaturePressureFunctionFluidProperties::e_from_p_T(Real pressure, Real tempe return _e_ref + _cv * (temperature - _T_ref); else { - const Real h = (temperature - _T_ref) / _n_integration_dT; + const int n_intervals = std::ceil(std::abs(temperature - _T_ref) / _integration_dT); + const auto h = (temperature - _T_ref) / n_intervals; Real integral = 0; // Centered step integration is second-order - for (const auto i : make_range(_n_integration_dT)) + for (const auto i : make_range(n_intervals)) integral += cv_from_p_T(pressure, _T_ref + (i + 0.5) * h); integral *= h; // we are still missing the dV or dP term to go from V/P_ref (e_ref = e(T_ref, V/P_ref)) diff --git a/modules/fluid_properties/unit/include/TemperaturePressureFunctionFluidPropertiesTest.h b/modules/fluid_properties/unit/include/TemperaturePressureFunctionFluidPropertiesTest.h index e49e65776cf3..093de73b16cb 100644 --- a/modules/fluid_properties/unit/include/TemperaturePressureFunctionFluidPropertiesTest.h +++ b/modules/fluid_properties/unit/include/TemperaturePressureFunctionFluidPropertiesTest.h @@ -70,7 +70,7 @@ class TemperaturePressureFunctionFluidPropertiesTest : public MooseObjectUnitTes uo_params2.set("cp") = "cp_function"; uo_params2.set("T_initial_guess") = 250; uo_params2.set("p_initial_guess") = 1e7; - uo_params2.set("dT_integration_intervals") = 2000; + uo_params2.set("dT_integration_intervals") = 0.01; _fe_problem->addUserObject("TemperaturePressureFunctionFluidProperties", "fp_cp", uo_params2); _fp_cp = &_fe_problem->getUserObject("fp_cp"); diff --git a/modules/fluid_properties/unit/src/TemperaturePressureFunctionFluidPropertiesTest.C b/modules/fluid_properties/unit/src/TemperaturePressureFunctionFluidPropertiesTest.C index ef992e470c62..f306d53c52dc 100644 --- a/modules/fluid_properties/unit/src/TemperaturePressureFunctionFluidPropertiesTest.C +++ b/modules/fluid_properties/unit/src/TemperaturePressureFunctionFluidPropertiesTest.C @@ -116,7 +116,7 @@ TEST_F(TemperaturePressureFunctionFluidPropertiesTest, properties_function_cp) Real p = 8.56E7; Real T = 200.0; Real cp = 3000. + 3 * T + 5e-4 * p; - Real e = 9220000.46051058359; + Real e = 9220000.46051060781; Real v = 1 / (1400 + 2.5 * T + 32e-5 * p); Real h = e + p * v; @@ -152,7 +152,7 @@ TEST_F(TemperaturePressureFunctionFluidPropertiesTest, properties_function_cp) Real p = 1.06841E8; Real T = 300.0; Real cp = 3000. + 3 * T + 5e-4 * p; - Real e = 17061150.6748719364; + Real e = 17061150.6748718396; Real v = 1 / (1400 + 2.5 * T + 32e-5 * p); Real h = e + p * v; @@ -213,7 +213,7 @@ TEST_F(TemperaturePressureFunctionFluidPropertiesTest, derivatives) DERIV_TEST(_fp->cp_from_v_e, v, e, large_tol); // uses finite differencing DERIV_TEST(_fp->cv_from_v_e, v, e, tol); - e = 80150.045818949831; + e = 80150.0458189499477; v = 0.00021621621621621621; AD_DERIV_TEST(_fp_cp->rho_from_p_T, p, T, tol); @@ -245,7 +245,7 @@ TEST_F(TemperaturePressureFunctionFluidPropertiesTest, derivatives) DERIV_TEST(_fp->cp_from_v_e, v, e, 2 * large_tol); // uses finite differencing DERIV_TEST(_fp->cv_from_v_e, v, e, tol); - e = 5374151.1232993276; + e = 5374151.12329934724; v = 5.5944055944055945e-05; AD_DERIV_TEST(_fp_cp->rho_from_p_T, p, T, tol);