diff --git a/modules/fluid_properties/doc/content/source/fluidproperties/TemperaturePressureFunctionFluidProperties.md b/modules/fluid_properties/doc/content/source/fluidproperties/TemperaturePressureFunctionFluidProperties.md index 943edc755514..b27db4546e4e 100644 --- a/modules/fluid_properties/doc/content/source/fluidproperties/TemperaturePressureFunctionFluidProperties.md +++ b/modules/fluid_properties/doc/content/source/fluidproperties/TemperaturePressureFunctionFluidProperties.md @@ -9,15 +9,36 @@ 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)) \\ - c_v = c_v^{user} \\ - e = c_v * T + \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} -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. +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} +\begin{aligned} + c_v = c_v^{u} \\ + 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^{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 $^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. @@ -31,9 +52,10 @@ 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, 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 5901ae2a244f..a69b615f326f 100644 --- a/modules/fluid_properties/include/fluidproperties/TemperaturePressureFunctionFluidProperties.h +++ b/modules/fluid_properties/include/fluidproperties/TemperaturePressureFunctionFluidProperties.h @@ -348,6 +348,17 @@ 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; + /// 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/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/src/fluidproperties/TemperaturePressureFunctionFluidProperties.C b/modules/fluid_properties/src/fluidproperties/TemperaturePressureFunctionFluidProperties.C index 3a1696a02d51..a91db7e2252c 100644 --- a/modules/fluid_properties/src/fluidproperties/TemperaturePressureFunctionFluidProperties.C +++ b/modules/fluid_properties/src/fluidproperties/TemperaturePressureFunctionFluidProperties.C @@ -23,8 +23,16 @@ 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 internal energy at the reference temperature"); + params.addParam("T_ref", 0, "Reference temperature for the specific internal energy"); + 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 " "conductivity, density, and viscosity as functions of temperature and pressure."); @@ -33,8 +41,16 @@ 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")), + _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."); } void @@ -43,6 +59,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 +70,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 +127,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 +153,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 +161,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 +243,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 @@ -234,30 +354,62 @@ 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; } 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 + { + 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_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)) + // to current V. The dT term is the largest one though + 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); + 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); + } } 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 @@ -271,54 +423,81 @@ TemperaturePressureFunctionFluidProperties::beta_from_p_T(Real pressure, Real te 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); + // 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 + { + 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); + // 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; + } } 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 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..093de73b16cb 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,22 @@ 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 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("dT_integration_intervals") = 0.01; + + _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/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); diff --git a/modules/fluid_properties/unit/src/TemperaturePressureFunctionFluidPropertiesTest.C b/modules/fluid_properties/unit/src/TemperaturePressureFunctionFluidPropertiesTest.C index 9a20103139ab..f306d53c52dc 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,90 @@ 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; + // 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.46051060781; + 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, 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, 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, newton_tol); + 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.6748718396; + 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, 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, 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, newton_tol); + 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 +192,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 +208,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.0458189499477; + 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); + 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 +240,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 = 5374151.12329934724; + 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); + 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); } /**