SinglePhaseFluidProperties
SinglePhaseFluidProperties is a base class for all single-phase fluid properties objects. Its main function is to provide interfaces for computing various properties from different combinations of other properties.
Properties
The following properties are considered in this class, where the "Name" column gives the identifier for the property used in the available interfaces:
| Name | Symbol | Description |
|---|---|---|
beta | Volumetric expansion coefficient | |
c | Speed of sound | |
cp | Isobaric specific heat capacity | |
cv | Isochoric specific heat capacity | |
e | Specific internal energy | |
g | Gibbs free energy | |
gamma | Ratio of specific heats, | |
h | Specific enthalpy | |
k | Thermal conductivity | |
mu | Dynamic viscosity | |
p | Pressure | |
rho | Density | |
s | Specific entropy | |
T | Temperature | |
v | Specific volume |
Because two independent, intensive thermodynamic properties define a thermodynamic state of a pure fluid, interfaces in these objects are of the form where is the desired thermodynamic property and and are independent, intensive thermodynamic properties that define the thermodynamic state. The corresponding function name is fname_from_aname_bname, where fname, aname, and bname are the names in the table above, corresponding to , , and , respectively. The following table lists which properties are available from various combinations of properties (e.g., "Yes" in the column for the row denotes that the interface fname_from_aname_bname is available):
| Name | |||||||||
|---|---|---|---|---|---|---|---|---|---|
| Yes | |||||||||
| Yes | Yes | ||||||||
| Yes | Yes | ||||||||
| Yes | Yes | Yes | |||||||
| Yes | Yes | Yes | Yes | ||||||
| Yes | |||||||||
| Yes | Yes | ||||||||
| Yes | Yes | ||||||||
| Yes | Yes | Yes | |||||||
| Yes | Yes | Yes | |||||||
| Yes | Yes | Yes | |||||||
| Yes | Yes | ||||||||
| Yes | Yes | Yes | Yes | ||||||
| Yes | Yes | ||||||||
| Yes |
Interfaces are also provided for getting derivatives of fluid properties with respect to the input arguments. These interfaces are named the same as their non-derivative counterparts, but have no return value but 3 additional (output) arguments, corresponding to the property value and then the derivatives of each of the two input arguments. For example, has the interface rho_from_p_T(p, T, rho, drho_dp, drho_dT), where drho_dp and drho_dT correspond to and , respectively.
commentnote:Automatic Differentiation
Fluid properties objects have interfaces for taking advantage of MOOSE's Automatic Differentiation capability. See the example in the next section.
Additionally, the following interfaces are available:
fluidName(): The fluid name.molarMass(): The fluid's molar mass (kg/mol).
The full list of available methods can be found in either the source code or the Modules Doxygen page for each FluidProperties class.
Default Analytical Fluid Properties Relations
SinglePhaseFluidProperties provides a number of default implementations for some fluid properties where analytical relations hold for all single phase fluid properties. Some of these fluid properties are also implemented along with their derivatives with regards to the input variables, when these derivatives can also be analytically described. Relevant automatic differentiation (AD) implementations are also provided through a macro to avoid duplicated code.
The full list of available methods can be found in either the source code or the Doxygen page.
Variable Set Conversions
Different fluid applications may require different variable sets, such as (pressure, temperature) or (specific volume, specific internal energy), depending on the flow regimes of interest and relatedly the numerical discretization. Fluid properties are not necessarily implemented or known for all variable sets, so conversions from one variable set to another can be helpful.
For many fluids, analytical closures for these conversions are not known, so SinglePhaseFluidProperties defines several routines for iteratively converting from one variable set to another. This leverages the numerical inversion methods utilities. Notably, the following routines are provided:
void p_T_from_v_e(const CppType & v,
const CppType & e,
Real p0,
Real T0,
CppType & p,
CppType & T,
bool & conversion_succeeded) const;
void p_T_from_v_h(const T & v,
const T & h,
Real p0,
Real T0,
T & pressure,
T & temperature,
bool & conversion_succeeded) const;
void p_T_from_h_s(const T & h,
const T & s,
Real p0,
Real T0,
T & pressure,
T & temperature,
bool & conversion_succeeded) const;
These routines may then be used to convert from one variable set to another before obtaining the desired fluid property. For example, this routine converts (pressure, temperature) to (specific volume, specific energy) to compute entropy.
SinglePhaseFluidProperties::s_from_p_T(const Real pressure, const Real temperature) const
{
Real v, e;
v_e_from_p_T(pressure, temperature, v, e);
return s_from_v_e(v, e);
}
(moose/modules/fluid_properties/include/fluidproperties/SinglePhaseFluidProperties.h)
// This file is part of the MOOSE framework
// https://www.mooseframework.org
//
// All rights reserved, see COPYRIGHT for full restrictions
// https://github.com/idaholab/moose/blob/master/COPYRIGHT
//
// Licensed under LGPL 2.1, please see LICENSE for details
// https://www.gnu.org/licenses/lgpl-2.1.html
#pragma once
#include "FluidProperties.h"
#include "NewtonInversion.h"
#include "metaphysicl/dualnumberarray.h"
/**
* Adds AD versions of each fluid property. These functions use the Real versions of these methods
* to compute the AD variables complete with derivatives. Typically, these do not need to be
* overriden in derived classes.
*/
#define propfuncAD(want, prop1, prop2) \
virtual ADReal want##_from_##prop1##_##prop2(const ADReal & p1, const ADReal & p2) const \
{ \
Real x = 0; \
Real raw1 = p1.value(); \
Real raw2 = p2.value(); \
Real dxd1 = 0; \
Real dxd2 = 0; \
want##_from_##prop1##_##prop2(raw1, raw2, x, dxd1, dxd2); \
\
ADReal result = x; \
result.derivatives() = p1.derivatives() * dxd1 + p2.derivatives() * dxd2; \
return result; \
} \
\
virtual void want##_from_##prop1##_##prop2(const ADReal & prop1, \
const ADReal & prop2, \
ADReal & val, \
ADReal & d##want##d1, \
ADReal & d##want##d2) const \
{ \
unimplementedDerivativeMethod(__PRETTY_FUNCTION__); \
Real dummy, tmp1, tmp2; \
val = want##_from_##prop1##_##prop2(prop1, prop2); \
want##_from_##prop1##_##prop2(prop1.value(), prop2.value(), dummy, tmp1, tmp2); \
d##want##d1 = tmp1; \
d##want##d2 = tmp2; \
}
/**
* Adds function definitions with not implemented error. These functions should be overriden in
* derived classes where required. AD versions are constructed automatically using propfuncAD.
*/
#define propfunc(want, prop1, prop2) \
virtual Real want##_from_##prop1##_##prop2(Real, Real) const \
{ \
mooseError( \
"The fluid properties class '", \
type(), \
"' has not implemented the method below. If your application requires this method, you " \
"must either implement it or use a different fluid properties class.\n\n", \
__PRETTY_FUNCTION__); \
} \
\
virtual void want##_from_##prop1##_##prop2( \
Real prop1, Real prop2, Real & val, Real & d##want##d1, Real & d##want##d2) const \
{ \
unimplementedDerivativeMethod(__PRETTY_FUNCTION__); \
d##want##d1 = 0; \
d##want##d2 = 0; \
val = want##_from_##prop1##_##prop2(prop1, prop2); \
} \
\
propfuncAD(want, prop1, prop2)
/**
* Adds Real declarations of functions that have a default implementation.
* Important: properties declared using this macro must be defined in SinglePhaseFluidProperties.C.
* AD versions are constructed automatically using propfuncAD.
*/
#define propfuncWithDefault(want, prop1, prop2) \
virtual Real want##_from_##prop1##_##prop2(Real, Real) const; \
virtual void want##_from_##prop1##_##prop2( \
Real prop1, Real prop2, Real & val, Real & d##want##d1, Real & d##want##d2) const; \
\
propfuncAD(want, prop1, prop2)
/**
* Adds Real and ADReal declarations of functions that have an implementation.
*/
#define propfuncWithDefinitionOverride(want, prop1, prop2) \
Real want##_from_##prop1##_##prop2(Real, Real) const override; \
void want##_from_##prop1##_##prop2( \
Real prop1, Real prop2, Real & val, Real & d##want##d1, Real & d##want##d2) const override; \
ADReal want##_from_##prop1##_##prop2(const ADReal &, const ADReal &) const override; \
void want##_from_##prop1##_##prop2(const ADReal & prop1, \
const ADReal & prop2, \
ADReal & val, \
ADReal & d##want##d1, \
ADReal & d##want##d2) const override; \
template <typename CppType> \
CppType want##_from_##prop1##_##prop2##_template(const CppType & prop1, const CppType & prop2) \
const; \
template <typename CppType> \
void want##_from_##prop1##_##prop2##_template(const CppType & prop1, \
const CppType & prop2, \
CppType & val, \
CppType & d##want##d1, \
CppType & d##want##d2) const
/**
* Common class for single phase fluid properties
*/
class SinglePhaseFluidProperties : public FluidProperties
{
public:
static InputParameters validParams();
SinglePhaseFluidProperties(const InputParameters & parameters);
virtual ~SinglePhaseFluidProperties();
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Woverloaded-virtual"
// clang-format off
/**
* @brief Compute a fluid property given for the state defined by two given properties.
*
* For all functions, the first two arguments are the given properties that define the fluid
* state. For the two-argument variants, the desired property is the return value.
* The five-argument variants also provide partial derivatives dx/da and dx/db where x is the
* desired property being computed, a is the first given property, and b is the second given
* property. The desired property, dx/da, and dx/db are stored into the 3rd, 4th, and 5th
* arguments respectively.
*
* Properties/parameters used in these function are listed below with their units:
*
* @begincode
* p pressure [Pa]
* T temperature [K]
* e specific internal energy [J/kg]
* v specific volume [m^3/kg]
* rho density [kg/m^3]
* h specific enthalpy [J/kg]
* s specific entropy [J/(kg*K)]
* mu viscosity [Pa*s]
* k thermal conductivity [W/(m*K)]
* c speed of sound [m/s]
* cp constant-pressure specific heat [J/K]
* cv constant-volume specific heat [J/K]
* beta volumetric thermal expansion coefficient [1/K]
* g Gibbs free energy [J]
* pp_sat partial pressure at saturation [Pa]
* gamma Adiabatic ratio (cp/cv) [-]
* @endcode
*
* As an example:
*
* @begincode
* // calculate pressure given specific vol and energy:
* auto pressure = your_fluid_properties_object.p_from_v_e(specific_vol, specific_energy);
*
* // or use the derivative variant:
* Real dp_dv = 0; // derivative will be stored into here
* Real dp_de = 0; // derivative will be stored into here
* your_fluid_properties_object.p_from_v_e(specific_vol, specific_energy, pressure, dp_dv, dp_de);
* @endcode
*
* Automatic differentiation (AD) support is provided through x_from_a_b(ADReal a, ADReal b) and
* x_from_a_b(ADReal a, ADReal b, ADReal x, ADReal dx_da, ADReal dx_db) versions of the
* functions where a and b must be ADReal/DualNumber's calculated using all AD-supporting values:
*
* @begincode
* auto v = 1/rho; // rho must be an AD non-linear variable.
* auto e = rhoE/rho - vel_energy; // rhoE and vel_energy must be AD variables/numbers also.
* auto pressure = your_fluid_properties_object.p_from_v_e(v, e);
* // pressure now contains partial derivatives w.r.t. all degrees of freedom
* @endcode
*/
///@{
propfunc(p, v, e)
propfunc(T, v, e)
propfunc(c, v, e)
propfunc(cp, v, e)
propfunc(cv, v, e)
propfunc(mu, v, e)
propfunc(k, v, e)
propfuncWithDefault(s, v, e)
propfunc(s, h, p)
propfunc(rho, p, s)
propfunc(e, v, h)
propfuncWithDefault(s, p, T)
propfunc(pp_sat, p, T)
propfunc(mu, rho, T)
propfunc(k, rho, T)
propfuncWithDefault(c, p, T)
propfuncWithDefault(cp, p, T)
propfuncWithDefault(cv, p, T)
propfuncWithDefault(mu, p, T)
propfuncWithDefault(k, p, T)
propfunc(rho, p, T)
propfunc(e, p, rho)
propfunc(e, T, v)
propfunc(p, T, v)
propfunc(h, T, v)
propfunc(s, T, v)
propfunc(cv, T, v)
propfunc(h, p, T)
propfuncWithDefault(h, v, e)
propfunc(g, v, e)
propfuncWithDefault(p, h, s)
propfunc(T, h, p) // temporary, until uniformization
propfuncWithDefault(T, p, h)
propfuncWithDefault(beta, p, T)
propfuncWithDefault(v, p, T)
propfuncWithDefault(e, p, T)
propfuncWithDefault(gamma, v, e)
propfuncWithDefault(gamma, p, T)
///@}
// clang-format on
#undef propfunc
#undef propfuncWithDefault
#undef propfuncAD
/**
* Fluid name
* @return string representing fluid name
*/
virtual std::string fluidName() const;
/**
* Molar mass [kg/mol]
* @return molar mass
*/
virtual Real molarMass() const;
/**
* Critical pressure
* @return critical pressure (Pa)
*/
virtual Real criticalPressure() const;
/**
* Critical temperature
* @return critical temperature (K)
*/
virtual Real criticalTemperature() const;
/**
* Critical density
* @return critical density (kg/m^3)
*/
virtual Real criticalDensity() const;
/**
* Critical specific internal energy
* @return specific internal energy (J/kg)
*/
virtual Real criticalInternalEnergy() const;
/**
* Triple point pressure
* @return triple point pressure (Pa)
*/
virtual Real triplePointPressure() const;
/**
* Triple point temperature
* @return triple point temperature (K)
*/
virtual Real triplePointTemperature() const;
/**
* Specific internal energy from temperature and specific volume
*
* @param[in] T temperature
* @param[in] v specific volume
*/
virtual Real e_spndl_from_v(Real v) const;
/**
* Specific internal energy from temperature and specific volume
*
* @param[in] T temperature
* @param[in] v specific volume
*/
virtual void v_e_spndl_from_T(Real T, Real & v, Real & e) const;
/**
* Vapor pressure. Used to delineate liquid and gas phases.
* Valid for temperatures between the triple point temperature
* and the critical temperature
*
* @param T fluid temperature (K)
* @param[out] saturation pressure (Pa)
* @param[out] derivative of saturation pressure wrt temperature (Pa/K)
*/
virtual Real vaporPressure(Real T) const;
virtual void vaporPressure(Real T, Real & psat, Real & dpsat_dT) const;
virtual ADReal vaporPressure(const ADReal & T) const;
/**
* Vapor temperature. Used to delineate liquid and gas phases.
* Valid for pressures between the triple point pressure
* and the critical pressure
*
* @param p fluid pressure (Pa)
* @param[out] saturation temperature (K)
* @param[out] derivative of saturation temperature wrt pressure
*/
virtual Real vaporTemperature(Real p) const;
virtual void vaporTemperature(Real p, Real & Tsat, Real & dTsat_dp) const;
virtual ADReal vaporTemperature(const ADReal & p) const;
/**
* Henry's law coefficients for dissolution in water
* @return Henry's constant coefficients
*/
virtual std::vector<Real> henryCoefficients() const;
template <typename CppType>
void v_e_from_p_T(const CppType & p, const CppType & T, CppType & v, CppType & e) const;
template <typename CppType>
void v_e_from_p_T(const CppType & p,
const CppType & T,
CppType & v,
CppType & dv_dp,
CppType & dv_dT,
CppType & e,
CppType & de_dp,
CppType & de_dT) const;
/**
* Combined methods. These methods are particularly useful for the PorousFlow
* module, where density and viscosity are typically both computed everywhere.
* The combined methods allow the most efficient means of calculating both
* properties, especially where rho(p, T) and mu(rho, T). In this case, an
* extra density calculation would be required to calculate mu(p, T). All
* property names are described above.
*/
virtual void rho_mu_from_p_T(Real p, Real T, Real & rho, Real & mu) const;
virtual void rho_mu_from_p_T(Real p,
Real T,
Real & rho,
Real & drho_dp,
Real & drho_dT,
Real & mu,
Real & dmu_dp,
Real & dmu_dT) const;
virtual void rho_mu_from_p_T(const ADReal & p, const ADReal & T, ADReal & rho, ADReal & mu) const;
virtual void rho_e_from_p_T(Real p,
Real T,
Real & rho,
Real & drho_dp,
Real & drho_dT,
Real & e,
Real & de_dp,
Real & de_dT) const;
/**
* Determines (p,T) from (v,e) using Newton Solve in 2D
* Useful for conversion between different sets of state variables
*
* @param[in] v specific volume (m^3 / kg)
* @param[in] e specific internal energy (J / kg)
* @param[in] p0 initial guess for pressure (Pa / kg)
* @param[in] T0 initial guess for temperature (K)
* @param[out] fluid pressure (Pa / kg)
* @param[out] Temperature (K)
*/
template <typename CppType>
void p_T_from_v_e(const CppType & v,
const CppType & e,
Real p0,
Real T0,
CppType & p,
CppType & T,
bool & conversion_succeeded) const;
/**
* Determines (p,T) from (v,h) using Newton Solve in 2D
* Useful for conversion between different sets of state variables
*
* @param[in] v specific volume (m^3 / kg)
* @param[in] h specific enthalpy (J / kg)
* @param[in] p0 initial guess for pressure (Pa / kg)
* @param[in] T0 initial guess for temperature (K)
* @param[out] fluid pressure (Pa / kg)
* @param[out] Temperature (K)
*/
template <typename T>
void p_T_from_v_h(const T & v,
const T & h,
Real p0,
Real T0,
T & pressure,
T & temperature,
bool & conversion_succeeded) const;
/**
* Determines (p,T) from (h,s) using Newton Solve in 2D
* Useful for conversion between different sets of state variables
*
* @param[in] h specific enthalpy (J / kg)
* @param[in] s specific entropy (J/K*kg)
* @param[in] p0 initial guess for pressure (Pa / kg)
* @param[in] T0 initial guess for temperature (K)
* @param[out] fluid pressure (Pa / kg)
* @param[out] Temperature (K)
*/
template <typename T>
void p_T_from_h_s(const T & h,
const T & s,
Real p0,
Real T0,
T & pressure,
T & temperature,
bool & conversion_succeeded) const;
protected:
/**
* Computes the dependent variable z and its derivatives with respect to the independent
* variables x and y using the simple two parameter \p z_from_x_y functor. The derivatives are
* computed using a compound automatic differentiation type
*/
template <typename T, typename Functor>
static void
xyDerivatives(const T x, const T & y, T & z, T & dz_dx, T & dz_dy, const Functor & z_from_x_y);
/**
* Given a type example, this method returns zero and unity representations of that type (first
* and second members of returned pair respectively)
*/
template <typename T>
static std::pair<T, T> makeZeroAndOne(const T &);
/**
* Newton's method may be used to convert between variable sets
* _tolerance, _T_initial_guess, and _p_initial_guess are the parameters for these
* iterative solves
*/
const Real _tolerance;
const Real _T_initial_guess;
const Real _p_initial_guess;
private:
void unimplementedDerivativeMethod(const std::string & property_function_name) const
{
const std::string message =
"The fluid properties class '" + type() +
"' has not implemented the method below, which computes derivatives of fluid properties "
"with regards to the flow variables. If your application requires this "
"method, you must either implement it or use a different fluid properties "
" class.\n\n" +
property_function_name;
if (_allow_imperfect_jacobians)
mooseDoOnce(mooseWarning(message + "\nThe unimplemented derivatives for this fluid property "
"are currently neglected, set to 0."));
else
mooseError(message + "\n\nYou can avoid this error by neglecting the "
"unimplemented derivatives of fluid properties by setting the "
"'allow_imperfect_jacobians' parameter");
}
};
#pragma GCC diagnostic pop
template <typename T>
std::pair<T, T>
SinglePhaseFluidProperties::makeZeroAndOne(const T & /*ex*/)
{
return {T{0, 0}, T{1, 0}};
}
template <>
inline std::pair<Real, Real>
SinglePhaseFluidProperties::makeZeroAndOne(const Real & /*ex*/)
{
return {Real{0}, Real{1}};
}
template <typename T, typename Functor>
void
SinglePhaseFluidProperties::xyDerivatives(
const T x, const T & y, T & z, T & dz_dx, T & dz_dy, const Functor & z_from_x_y)
{
typedef MetaPhysicL::DualNumber<T, MetaPhysicL::NumberArray<2, T>> CompoundType;
const auto [zero, one] = makeZeroAndOne(x);
CompoundType x_c(x, zero);
auto & x_cd = x_c.derivatives();
x_cd[0] = one;
CompoundType y_c(y, zero);
auto & y_cd = y_c.derivatives();
y_cd[1] = one;
const auto z_c = z_from_x_y(x_c, y_c);
z = z_c.value();
dz_dx = z_c.derivatives()[0];
dz_dy = z_c.derivatives()[1];
}
template <typename CppType>
void
SinglePhaseFluidProperties::p_T_from_v_e(const CppType & v, // v value
const CppType & e, // e value
const Real p0, // initial guess
const Real T0, // initial guess
CppType & p, // returned pressure
CppType & T, // returned temperature
bool & conversion_succeeded) const
{
auto v_lambda = [&](const CppType & pressure,
const CppType & temperature,
CppType & new_v,
CppType & dv_dp,
CppType & dv_dT) { v_from_p_T(pressure, temperature, new_v, dv_dp, dv_dT); };
auto e_lambda = [&](const CppType & pressure,
const CppType & temperature,
CppType & new_e,
CppType & de_dp,
CppType & de_dT) { e_from_p_T(pressure, temperature, new_e, de_dp, de_dT); };
try
{
FluidPropertiesUtils::NewtonSolve2D(
v, e, p0, T0, p, T, _tolerance, _tolerance, v_lambda, e_lambda);
conversion_succeeded = true;
}
catch (MooseException &)
{
conversion_succeeded = false;
}
if (!conversion_succeeded)
mooseDoOnce(mooseWarning("Conversion from (v, e)=(", v, ", ", e, ") to (p, T) failed"));
}
template <typename T>
void
SinglePhaseFluidProperties::p_T_from_v_h(const T & v, // v value
const T & h, // e value
const Real p0, // initial guess
const Real T0, // initial guess
T & pressure, // returned pressure
T & temperature, // returned temperature
bool & conversion_succeeded) const
{
auto v_lambda = [&](const T & pressure, const T & temperature, T & new_v, T & dv_dp, T & dv_dT)
{ v_from_p_T(pressure, temperature, new_v, dv_dp, dv_dT); };
auto h_lambda = [&](const T & pressure, const T & temperature, T & new_h, T & dh_dp, T & dh_dT)
{ h_from_p_T(pressure, temperature, new_h, dh_dp, dh_dT); };
try
{
FluidPropertiesUtils::NewtonSolve2D(
v, h, p0, T0, pressure, temperature, _tolerance, _tolerance, v_lambda, h_lambda);
conversion_succeeded = true;
}
catch (MooseException &)
{
conversion_succeeded = false;
}
if (!conversion_succeeded)
mooseDoOnce(mooseWarning("Conversion from (v, h)=(", v, ", ", h, ") to (p, T) failed"));
}
template <typename T>
void
SinglePhaseFluidProperties::p_T_from_h_s(const T & h, // h value
const T & s, // s value
const Real p0, // initial guess
const Real T0, // initial guess
T & pressure, // returned pressure
T & temperature, // returned temperature
bool & conversion_succeeded) const
{
auto h_lambda = [&](const T & pressure, const T & temperature, T & new_h, T & dh_dp, T & dh_dT)
{ h_from_p_T(pressure, temperature, new_h, dh_dp, dh_dT); };
auto s_lambda = [&](const T & pressure, const T & temperature, T & new_s, T & ds_dp, T & ds_dT)
{ s_from_p_T(pressure, temperature, new_s, ds_dp, ds_dT); };
try
{
FluidPropertiesUtils::NewtonSolve2D(
h, s, p0, T0, pressure, temperature, _tolerance, _tolerance, h_lambda, s_lambda);
conversion_succeeded = true;
}
catch (MooseException &)
{
conversion_succeeded = false;
}
if (!conversion_succeeded)
mooseDoOnce(mooseWarning("Conversion from (h, s)=(", h, ", ", s, ") to (p, T) failed"));
}
template <typename CppType>
void
SinglePhaseFluidProperties::v_e_from_p_T(const CppType & p,
const CppType & T,
CppType & v,
CppType & e) const
{
const CppType rho = rho_from_p_T(p, T);
v = 1.0 / rho;
try
{
// more likely to not involve a Newton search
e = e_from_p_T(p, T);
}
catch (...)
{
e = e_from_p_rho(p, rho);
}
}
template <typename CppType>
void
SinglePhaseFluidProperties::v_e_from_p_T(const CppType & p,
const CppType & T,
CppType & v,
CppType & dv_dp,
CppType & dv_dT,
CppType & e,
CppType & de_dp,
CppType & de_dT) const
{
CppType rho, drho_dp, drho_dT;
rho_from_p_T(p, T, rho, drho_dp, drho_dT);
v = 1.0 / rho;
const CppType dv_drho = -1.0 / (rho * rho);
dv_dp = dv_drho * drho_dp;
dv_dT = dv_drho * drho_dT;
CppType de_dp_partial, de_drho;
e_from_p_rho(p, rho, e, de_dp_partial, de_drho);
de_dp = de_dp_partial + de_drho * drho_dp;
de_dT = de_drho * drho_dT;
}
(moose/modules/fluid_properties/include/fluidproperties/SinglePhaseFluidProperties.h)
// This file is part of the MOOSE framework
// https://www.mooseframework.org
//
// All rights reserved, see COPYRIGHT for full restrictions
// https://github.com/idaholab/moose/blob/master/COPYRIGHT
//
// Licensed under LGPL 2.1, please see LICENSE for details
// https://www.gnu.org/licenses/lgpl-2.1.html
#pragma once
#include "FluidProperties.h"
#include "NewtonInversion.h"
#include "metaphysicl/dualnumberarray.h"
/**
* Adds AD versions of each fluid property. These functions use the Real versions of these methods
* to compute the AD variables complete with derivatives. Typically, these do not need to be
* overriden in derived classes.
*/
#define propfuncAD(want, prop1, prop2) \
virtual ADReal want##_from_##prop1##_##prop2(const ADReal & p1, const ADReal & p2) const \
{ \
Real x = 0; \
Real raw1 = p1.value(); \
Real raw2 = p2.value(); \
Real dxd1 = 0; \
Real dxd2 = 0; \
want##_from_##prop1##_##prop2(raw1, raw2, x, dxd1, dxd2); \
\
ADReal result = x; \
result.derivatives() = p1.derivatives() * dxd1 + p2.derivatives() * dxd2; \
return result; \
} \
\
virtual void want##_from_##prop1##_##prop2(const ADReal & prop1, \
const ADReal & prop2, \
ADReal & val, \
ADReal & d##want##d1, \
ADReal & d##want##d2) const \
{ \
unimplementedDerivativeMethod(__PRETTY_FUNCTION__); \
Real dummy, tmp1, tmp2; \
val = want##_from_##prop1##_##prop2(prop1, prop2); \
want##_from_##prop1##_##prop2(prop1.value(), prop2.value(), dummy, tmp1, tmp2); \
d##want##d1 = tmp1; \
d##want##d2 = tmp2; \
}
/**
* Adds function definitions with not implemented error. These functions should be overriden in
* derived classes where required. AD versions are constructed automatically using propfuncAD.
*/
#define propfunc(want, prop1, prop2) \
virtual Real want##_from_##prop1##_##prop2(Real, Real) const \
{ \
mooseError( \
"The fluid properties class '", \
type(), \
"' has not implemented the method below. If your application requires this method, you " \
"must either implement it or use a different fluid properties class.\n\n", \
__PRETTY_FUNCTION__); \
} \
\
virtual void want##_from_##prop1##_##prop2( \
Real prop1, Real prop2, Real & val, Real & d##want##d1, Real & d##want##d2) const \
{ \
unimplementedDerivativeMethod(__PRETTY_FUNCTION__); \
d##want##d1 = 0; \
d##want##d2 = 0; \
val = want##_from_##prop1##_##prop2(prop1, prop2); \
} \
\
propfuncAD(want, prop1, prop2)
/**
* Adds Real declarations of functions that have a default implementation.
* Important: properties declared using this macro must be defined in SinglePhaseFluidProperties.C.
* AD versions are constructed automatically using propfuncAD.
*/
#define propfuncWithDefault(want, prop1, prop2) \
virtual Real want##_from_##prop1##_##prop2(Real, Real) const; \
virtual void want##_from_##prop1##_##prop2( \
Real prop1, Real prop2, Real & val, Real & d##want##d1, Real & d##want##d2) const; \
\
propfuncAD(want, prop1, prop2)
/**
* Adds Real and ADReal declarations of functions that have an implementation.
*/
#define propfuncWithDefinitionOverride(want, prop1, prop2) \
Real want##_from_##prop1##_##prop2(Real, Real) const override; \
void want##_from_##prop1##_##prop2( \
Real prop1, Real prop2, Real & val, Real & d##want##d1, Real & d##want##d2) const override; \
ADReal want##_from_##prop1##_##prop2(const ADReal &, const ADReal &) const override; \
void want##_from_##prop1##_##prop2(const ADReal & prop1, \
const ADReal & prop2, \
ADReal & val, \
ADReal & d##want##d1, \
ADReal & d##want##d2) const override; \
template <typename CppType> \
CppType want##_from_##prop1##_##prop2##_template(const CppType & prop1, const CppType & prop2) \
const; \
template <typename CppType> \
void want##_from_##prop1##_##prop2##_template(const CppType & prop1, \
const CppType & prop2, \
CppType & val, \
CppType & d##want##d1, \
CppType & d##want##d2) const
/**
* Common class for single phase fluid properties
*/
class SinglePhaseFluidProperties : public FluidProperties
{
public:
static InputParameters validParams();
SinglePhaseFluidProperties(const InputParameters & parameters);
virtual ~SinglePhaseFluidProperties();
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Woverloaded-virtual"
// clang-format off
/**
* @brief Compute a fluid property given for the state defined by two given properties.
*
* For all functions, the first two arguments are the given properties that define the fluid
* state. For the two-argument variants, the desired property is the return value.
* The five-argument variants also provide partial derivatives dx/da and dx/db where x is the
* desired property being computed, a is the first given property, and b is the second given
* property. The desired property, dx/da, and dx/db are stored into the 3rd, 4th, and 5th
* arguments respectively.
*
* Properties/parameters used in these function are listed below with their units:
*
* @begincode
* p pressure [Pa]
* T temperature [K]
* e specific internal energy [J/kg]
* v specific volume [m^3/kg]
* rho density [kg/m^3]
* h specific enthalpy [J/kg]
* s specific entropy [J/(kg*K)]
* mu viscosity [Pa*s]
* k thermal conductivity [W/(m*K)]
* c speed of sound [m/s]
* cp constant-pressure specific heat [J/K]
* cv constant-volume specific heat [J/K]
* beta volumetric thermal expansion coefficient [1/K]
* g Gibbs free energy [J]
* pp_sat partial pressure at saturation [Pa]
* gamma Adiabatic ratio (cp/cv) [-]
* @endcode
*
* As an example:
*
* @begincode
* // calculate pressure given specific vol and energy:
* auto pressure = your_fluid_properties_object.p_from_v_e(specific_vol, specific_energy);
*
* // or use the derivative variant:
* Real dp_dv = 0; // derivative will be stored into here
* Real dp_de = 0; // derivative will be stored into here
* your_fluid_properties_object.p_from_v_e(specific_vol, specific_energy, pressure, dp_dv, dp_de);
* @endcode
*
* Automatic differentiation (AD) support is provided through x_from_a_b(ADReal a, ADReal b) and
* x_from_a_b(ADReal a, ADReal b, ADReal x, ADReal dx_da, ADReal dx_db) versions of the
* functions where a and b must be ADReal/DualNumber's calculated using all AD-supporting values:
*
* @begincode
* auto v = 1/rho; // rho must be an AD non-linear variable.
* auto e = rhoE/rho - vel_energy; // rhoE and vel_energy must be AD variables/numbers also.
* auto pressure = your_fluid_properties_object.p_from_v_e(v, e);
* // pressure now contains partial derivatives w.r.t. all degrees of freedom
* @endcode
*/
///@{
propfunc(p, v, e)
propfunc(T, v, e)
propfunc(c, v, e)
propfunc(cp, v, e)
propfunc(cv, v, e)
propfunc(mu, v, e)
propfunc(k, v, e)
propfuncWithDefault(s, v, e)
propfunc(s, h, p)
propfunc(rho, p, s)
propfunc(e, v, h)
propfuncWithDefault(s, p, T)
propfunc(pp_sat, p, T)
propfunc(mu, rho, T)
propfunc(k, rho, T)
propfuncWithDefault(c, p, T)
propfuncWithDefault(cp, p, T)
propfuncWithDefault(cv, p, T)
propfuncWithDefault(mu, p, T)
propfuncWithDefault(k, p, T)
propfunc(rho, p, T)
propfunc(e, p, rho)
propfunc(e, T, v)
propfunc(p, T, v)
propfunc(h, T, v)
propfunc(s, T, v)
propfunc(cv, T, v)
propfunc(h, p, T)
propfuncWithDefault(h, v, e)
propfunc(g, v, e)
propfuncWithDefault(p, h, s)
propfunc(T, h, p) // temporary, until uniformization
propfuncWithDefault(T, p, h)
propfuncWithDefault(beta, p, T)
propfuncWithDefault(v, p, T)
propfuncWithDefault(e, p, T)
propfuncWithDefault(gamma, v, e)
propfuncWithDefault(gamma, p, T)
///@}
// clang-format on
#undef propfunc
#undef propfuncWithDefault
#undef propfuncAD
/**
* Fluid name
* @return string representing fluid name
*/
virtual std::string fluidName() const;
/**
* Molar mass [kg/mol]
* @return molar mass
*/
virtual Real molarMass() const;
/**
* Critical pressure
* @return critical pressure (Pa)
*/
virtual Real criticalPressure() const;
/**
* Critical temperature
* @return critical temperature (K)
*/
virtual Real criticalTemperature() const;
/**
* Critical density
* @return critical density (kg/m^3)
*/
virtual Real criticalDensity() const;
/**
* Critical specific internal energy
* @return specific internal energy (J/kg)
*/
virtual Real criticalInternalEnergy() const;
/**
* Triple point pressure
* @return triple point pressure (Pa)
*/
virtual Real triplePointPressure() const;
/**
* Triple point temperature
* @return triple point temperature (K)
*/
virtual Real triplePointTemperature() const;
/**
* Specific internal energy from temperature and specific volume
*
* @param[in] T temperature
* @param[in] v specific volume
*/
virtual Real e_spndl_from_v(Real v) const;
/**
* Specific internal energy from temperature and specific volume
*
* @param[in] T temperature
* @param[in] v specific volume
*/
virtual void v_e_spndl_from_T(Real T, Real & v, Real & e) const;
/**
* Vapor pressure. Used to delineate liquid and gas phases.
* Valid for temperatures between the triple point temperature
* and the critical temperature
*
* @param T fluid temperature (K)
* @param[out] saturation pressure (Pa)
* @param[out] derivative of saturation pressure wrt temperature (Pa/K)
*/
virtual Real vaporPressure(Real T) const;
virtual void vaporPressure(Real T, Real & psat, Real & dpsat_dT) const;
virtual ADReal vaporPressure(const ADReal & T) const;
/**
* Vapor temperature. Used to delineate liquid and gas phases.
* Valid for pressures between the triple point pressure
* and the critical pressure
*
* @param p fluid pressure (Pa)
* @param[out] saturation temperature (K)
* @param[out] derivative of saturation temperature wrt pressure
*/
virtual Real vaporTemperature(Real p) const;
virtual void vaporTemperature(Real p, Real & Tsat, Real & dTsat_dp) const;
virtual ADReal vaporTemperature(const ADReal & p) const;
/**
* Henry's law coefficients for dissolution in water
* @return Henry's constant coefficients
*/
virtual std::vector<Real> henryCoefficients() const;
template <typename CppType>
void v_e_from_p_T(const CppType & p, const CppType & T, CppType & v, CppType & e) const;
template <typename CppType>
void v_e_from_p_T(const CppType & p,
const CppType & T,
CppType & v,
CppType & dv_dp,
CppType & dv_dT,
CppType & e,
CppType & de_dp,
CppType & de_dT) const;
/**
* Combined methods. These methods are particularly useful for the PorousFlow
* module, where density and viscosity are typically both computed everywhere.
* The combined methods allow the most efficient means of calculating both
* properties, especially where rho(p, T) and mu(rho, T). In this case, an
* extra density calculation would be required to calculate mu(p, T). All
* property names are described above.
*/
virtual void rho_mu_from_p_T(Real p, Real T, Real & rho, Real & mu) const;
virtual void rho_mu_from_p_T(Real p,
Real T,
Real & rho,
Real & drho_dp,
Real & drho_dT,
Real & mu,
Real & dmu_dp,
Real & dmu_dT) const;
virtual void rho_mu_from_p_T(const ADReal & p, const ADReal & T, ADReal & rho, ADReal & mu) const;
virtual void rho_e_from_p_T(Real p,
Real T,
Real & rho,
Real & drho_dp,
Real & drho_dT,
Real & e,
Real & de_dp,
Real & de_dT) const;
/**
* Determines (p,T) from (v,e) using Newton Solve in 2D
* Useful for conversion between different sets of state variables
*
* @param[in] v specific volume (m^3 / kg)
* @param[in] e specific internal energy (J / kg)
* @param[in] p0 initial guess for pressure (Pa / kg)
* @param[in] T0 initial guess for temperature (K)
* @param[out] fluid pressure (Pa / kg)
* @param[out] Temperature (K)
*/
template <typename CppType>
void p_T_from_v_e(const CppType & v,
const CppType & e,
Real p0,
Real T0,
CppType & p,
CppType & T,
bool & conversion_succeeded) const;
/**
* Determines (p,T) from (v,h) using Newton Solve in 2D
* Useful for conversion between different sets of state variables
*
* @param[in] v specific volume (m^3 / kg)
* @param[in] h specific enthalpy (J / kg)
* @param[in] p0 initial guess for pressure (Pa / kg)
* @param[in] T0 initial guess for temperature (K)
* @param[out] fluid pressure (Pa / kg)
* @param[out] Temperature (K)
*/
template <typename T>
void p_T_from_v_h(const T & v,
const T & h,
Real p0,
Real T0,
T & pressure,
T & temperature,
bool & conversion_succeeded) const;
/**
* Determines (p,T) from (h,s) using Newton Solve in 2D
* Useful for conversion between different sets of state variables
*
* @param[in] h specific enthalpy (J / kg)
* @param[in] s specific entropy (J/K*kg)
* @param[in] p0 initial guess for pressure (Pa / kg)
* @param[in] T0 initial guess for temperature (K)
* @param[out] fluid pressure (Pa / kg)
* @param[out] Temperature (K)
*/
template <typename T>
void p_T_from_h_s(const T & h,
const T & s,
Real p0,
Real T0,
T & pressure,
T & temperature,
bool & conversion_succeeded) const;
protected:
/**
* Computes the dependent variable z and its derivatives with respect to the independent
* variables x and y using the simple two parameter \p z_from_x_y functor. The derivatives are
* computed using a compound automatic differentiation type
*/
template <typename T, typename Functor>
static void
xyDerivatives(const T x, const T & y, T & z, T & dz_dx, T & dz_dy, const Functor & z_from_x_y);
/**
* Given a type example, this method returns zero and unity representations of that type (first
* and second members of returned pair respectively)
*/
template <typename T>
static std::pair<T, T> makeZeroAndOne(const T &);
/**
* Newton's method may be used to convert between variable sets
* _tolerance, _T_initial_guess, and _p_initial_guess are the parameters for these
* iterative solves
*/
const Real _tolerance;
const Real _T_initial_guess;
const Real _p_initial_guess;
private:
void unimplementedDerivativeMethod(const std::string & property_function_name) const
{
const std::string message =
"The fluid properties class '" + type() +
"' has not implemented the method below, which computes derivatives of fluid properties "
"with regards to the flow variables. If your application requires this "
"method, you must either implement it or use a different fluid properties "
" class.\n\n" +
property_function_name;
if (_allow_imperfect_jacobians)
mooseDoOnce(mooseWarning(message + "\nThe unimplemented derivatives for this fluid property "
"are currently neglected, set to 0."));
else
mooseError(message + "\n\nYou can avoid this error by neglecting the "
"unimplemented derivatives of fluid properties by setting the "
"'allow_imperfect_jacobians' parameter");
}
};
#pragma GCC diagnostic pop
template <typename T>
std::pair<T, T>
SinglePhaseFluidProperties::makeZeroAndOne(const T & /*ex*/)
{
return {T{0, 0}, T{1, 0}};
}
template <>
inline std::pair<Real, Real>
SinglePhaseFluidProperties::makeZeroAndOne(const Real & /*ex*/)
{
return {Real{0}, Real{1}};
}
template <typename T, typename Functor>
void
SinglePhaseFluidProperties::xyDerivatives(
const T x, const T & y, T & z, T & dz_dx, T & dz_dy, const Functor & z_from_x_y)
{
typedef MetaPhysicL::DualNumber<T, MetaPhysicL::NumberArray<2, T>> CompoundType;
const auto [zero, one] = makeZeroAndOne(x);
CompoundType x_c(x, zero);
auto & x_cd = x_c.derivatives();
x_cd[0] = one;
CompoundType y_c(y, zero);
auto & y_cd = y_c.derivatives();
y_cd[1] = one;
const auto z_c = z_from_x_y(x_c, y_c);
z = z_c.value();
dz_dx = z_c.derivatives()[0];
dz_dy = z_c.derivatives()[1];
}
template <typename CppType>
void
SinglePhaseFluidProperties::p_T_from_v_e(const CppType & v, // v value
const CppType & e, // e value
const Real p0, // initial guess
const Real T0, // initial guess
CppType & p, // returned pressure
CppType & T, // returned temperature
bool & conversion_succeeded) const
{
auto v_lambda = [&](const CppType & pressure,
const CppType & temperature,
CppType & new_v,
CppType & dv_dp,
CppType & dv_dT) { v_from_p_T(pressure, temperature, new_v, dv_dp, dv_dT); };
auto e_lambda = [&](const CppType & pressure,
const CppType & temperature,
CppType & new_e,
CppType & de_dp,
CppType & de_dT) { e_from_p_T(pressure, temperature, new_e, de_dp, de_dT); };
try
{
FluidPropertiesUtils::NewtonSolve2D(
v, e, p0, T0, p, T, _tolerance, _tolerance, v_lambda, e_lambda);
conversion_succeeded = true;
}
catch (MooseException &)
{
conversion_succeeded = false;
}
if (!conversion_succeeded)
mooseDoOnce(mooseWarning("Conversion from (v, e)=(", v, ", ", e, ") to (p, T) failed"));
}
template <typename T>
void
SinglePhaseFluidProperties::p_T_from_v_h(const T & v, // v value
const T & h, // e value
const Real p0, // initial guess
const Real T0, // initial guess
T & pressure, // returned pressure
T & temperature, // returned temperature
bool & conversion_succeeded) const
{
auto v_lambda = [&](const T & pressure, const T & temperature, T & new_v, T & dv_dp, T & dv_dT)
{ v_from_p_T(pressure, temperature, new_v, dv_dp, dv_dT); };
auto h_lambda = [&](const T & pressure, const T & temperature, T & new_h, T & dh_dp, T & dh_dT)
{ h_from_p_T(pressure, temperature, new_h, dh_dp, dh_dT); };
try
{
FluidPropertiesUtils::NewtonSolve2D(
v, h, p0, T0, pressure, temperature, _tolerance, _tolerance, v_lambda, h_lambda);
conversion_succeeded = true;
}
catch (MooseException &)
{
conversion_succeeded = false;
}
if (!conversion_succeeded)
mooseDoOnce(mooseWarning("Conversion from (v, h)=(", v, ", ", h, ") to (p, T) failed"));
}
template <typename T>
void
SinglePhaseFluidProperties::p_T_from_h_s(const T & h, // h value
const T & s, // s value
const Real p0, // initial guess
const Real T0, // initial guess
T & pressure, // returned pressure
T & temperature, // returned temperature
bool & conversion_succeeded) const
{
auto h_lambda = [&](const T & pressure, const T & temperature, T & new_h, T & dh_dp, T & dh_dT)
{ h_from_p_T(pressure, temperature, new_h, dh_dp, dh_dT); };
auto s_lambda = [&](const T & pressure, const T & temperature, T & new_s, T & ds_dp, T & ds_dT)
{ s_from_p_T(pressure, temperature, new_s, ds_dp, ds_dT); };
try
{
FluidPropertiesUtils::NewtonSolve2D(
h, s, p0, T0, pressure, temperature, _tolerance, _tolerance, h_lambda, s_lambda);
conversion_succeeded = true;
}
catch (MooseException &)
{
conversion_succeeded = false;
}
if (!conversion_succeeded)
mooseDoOnce(mooseWarning("Conversion from (h, s)=(", h, ", ", s, ") to (p, T) failed"));
}
template <typename CppType>
void
SinglePhaseFluidProperties::v_e_from_p_T(const CppType & p,
const CppType & T,
CppType & v,
CppType & e) const
{
const CppType rho = rho_from_p_T(p, T);
v = 1.0 / rho;
try
{
// more likely to not involve a Newton search
e = e_from_p_T(p, T);
}
catch (...)
{
e = e_from_p_rho(p, rho);
}
}
template <typename CppType>
void
SinglePhaseFluidProperties::v_e_from_p_T(const CppType & p,
const CppType & T,
CppType & v,
CppType & dv_dp,
CppType & dv_dT,
CppType & e,
CppType & de_dp,
CppType & de_dT) const
{
CppType rho, drho_dp, drho_dT;
rho_from_p_T(p, T, rho, drho_dp, drho_dT);
v = 1.0 / rho;
const CppType dv_drho = -1.0 / (rho * rho);
dv_dp = dv_drho * drho_dp;
dv_dT = dv_drho * drho_dT;
CppType de_dp_partial, de_drho;
e_from_p_rho(p, rho, e, de_dp_partial, de_drho);
de_dp = de_dp_partial + de_drho * drho_dp;
de_dT = de_drho * drho_dT;
}
(moose/modules/fluid_properties/include/fluidproperties/SinglePhaseFluidProperties.h)
// This file is part of the MOOSE framework
// https://www.mooseframework.org
//
// All rights reserved, see COPYRIGHT for full restrictions
// https://github.com/idaholab/moose/blob/master/COPYRIGHT
//
// Licensed under LGPL 2.1, please see LICENSE for details
// https://www.gnu.org/licenses/lgpl-2.1.html
#pragma once
#include "FluidProperties.h"
#include "NewtonInversion.h"
#include "metaphysicl/dualnumberarray.h"
/**
* Adds AD versions of each fluid property. These functions use the Real versions of these methods
* to compute the AD variables complete with derivatives. Typically, these do not need to be
* overriden in derived classes.
*/
#define propfuncAD(want, prop1, prop2) \
virtual ADReal want##_from_##prop1##_##prop2(const ADReal & p1, const ADReal & p2) const \
{ \
Real x = 0; \
Real raw1 = p1.value(); \
Real raw2 = p2.value(); \
Real dxd1 = 0; \
Real dxd2 = 0; \
want##_from_##prop1##_##prop2(raw1, raw2, x, dxd1, dxd2); \
\
ADReal result = x; \
result.derivatives() = p1.derivatives() * dxd1 + p2.derivatives() * dxd2; \
return result; \
} \
\
virtual void want##_from_##prop1##_##prop2(const ADReal & prop1, \
const ADReal & prop2, \
ADReal & val, \
ADReal & d##want##d1, \
ADReal & d##want##d2) const \
{ \
unimplementedDerivativeMethod(__PRETTY_FUNCTION__); \
Real dummy, tmp1, tmp2; \
val = want##_from_##prop1##_##prop2(prop1, prop2); \
want##_from_##prop1##_##prop2(prop1.value(), prop2.value(), dummy, tmp1, tmp2); \
d##want##d1 = tmp1; \
d##want##d2 = tmp2; \
}
/**
* Adds function definitions with not implemented error. These functions should be overriden in
* derived classes where required. AD versions are constructed automatically using propfuncAD.
*/
#define propfunc(want, prop1, prop2) \
virtual Real want##_from_##prop1##_##prop2(Real, Real) const \
{ \
mooseError( \
"The fluid properties class '", \
type(), \
"' has not implemented the method below. If your application requires this method, you " \
"must either implement it or use a different fluid properties class.\n\n", \
__PRETTY_FUNCTION__); \
} \
\
virtual void want##_from_##prop1##_##prop2( \
Real prop1, Real prop2, Real & val, Real & d##want##d1, Real & d##want##d2) const \
{ \
unimplementedDerivativeMethod(__PRETTY_FUNCTION__); \
d##want##d1 = 0; \
d##want##d2 = 0; \
val = want##_from_##prop1##_##prop2(prop1, prop2); \
} \
\
propfuncAD(want, prop1, prop2)
/**
* Adds Real declarations of functions that have a default implementation.
* Important: properties declared using this macro must be defined in SinglePhaseFluidProperties.C.
* AD versions are constructed automatically using propfuncAD.
*/
#define propfuncWithDefault(want, prop1, prop2) \
virtual Real want##_from_##prop1##_##prop2(Real, Real) const; \
virtual void want##_from_##prop1##_##prop2( \
Real prop1, Real prop2, Real & val, Real & d##want##d1, Real & d##want##d2) const; \
\
propfuncAD(want, prop1, prop2)
/**
* Adds Real and ADReal declarations of functions that have an implementation.
*/
#define propfuncWithDefinitionOverride(want, prop1, prop2) \
Real want##_from_##prop1##_##prop2(Real, Real) const override; \
void want##_from_##prop1##_##prop2( \
Real prop1, Real prop2, Real & val, Real & d##want##d1, Real & d##want##d2) const override; \
ADReal want##_from_##prop1##_##prop2(const ADReal &, const ADReal &) const override; \
void want##_from_##prop1##_##prop2(const ADReal & prop1, \
const ADReal & prop2, \
ADReal & val, \
ADReal & d##want##d1, \
ADReal & d##want##d2) const override; \
template <typename CppType> \
CppType want##_from_##prop1##_##prop2##_template(const CppType & prop1, const CppType & prop2) \
const; \
template <typename CppType> \
void want##_from_##prop1##_##prop2##_template(const CppType & prop1, \
const CppType & prop2, \
CppType & val, \
CppType & d##want##d1, \
CppType & d##want##d2) const
/**
* Common class for single phase fluid properties
*/
class SinglePhaseFluidProperties : public FluidProperties
{
public:
static InputParameters validParams();
SinglePhaseFluidProperties(const InputParameters & parameters);
virtual ~SinglePhaseFluidProperties();
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Woverloaded-virtual"
// clang-format off
/**
* @brief Compute a fluid property given for the state defined by two given properties.
*
* For all functions, the first two arguments are the given properties that define the fluid
* state. For the two-argument variants, the desired property is the return value.
* The five-argument variants also provide partial derivatives dx/da and dx/db where x is the
* desired property being computed, a is the first given property, and b is the second given
* property. The desired property, dx/da, and dx/db are stored into the 3rd, 4th, and 5th
* arguments respectively.
*
* Properties/parameters used in these function are listed below with their units:
*
* @begincode
* p pressure [Pa]
* T temperature [K]
* e specific internal energy [J/kg]
* v specific volume [m^3/kg]
* rho density [kg/m^3]
* h specific enthalpy [J/kg]
* s specific entropy [J/(kg*K)]
* mu viscosity [Pa*s]
* k thermal conductivity [W/(m*K)]
* c speed of sound [m/s]
* cp constant-pressure specific heat [J/K]
* cv constant-volume specific heat [J/K]
* beta volumetric thermal expansion coefficient [1/K]
* g Gibbs free energy [J]
* pp_sat partial pressure at saturation [Pa]
* gamma Adiabatic ratio (cp/cv) [-]
* @endcode
*
* As an example:
*
* @begincode
* // calculate pressure given specific vol and energy:
* auto pressure = your_fluid_properties_object.p_from_v_e(specific_vol, specific_energy);
*
* // or use the derivative variant:
* Real dp_dv = 0; // derivative will be stored into here
* Real dp_de = 0; // derivative will be stored into here
* your_fluid_properties_object.p_from_v_e(specific_vol, specific_energy, pressure, dp_dv, dp_de);
* @endcode
*
* Automatic differentiation (AD) support is provided through x_from_a_b(ADReal a, ADReal b) and
* x_from_a_b(ADReal a, ADReal b, ADReal x, ADReal dx_da, ADReal dx_db) versions of the
* functions where a and b must be ADReal/DualNumber's calculated using all AD-supporting values:
*
* @begincode
* auto v = 1/rho; // rho must be an AD non-linear variable.
* auto e = rhoE/rho - vel_energy; // rhoE and vel_energy must be AD variables/numbers also.
* auto pressure = your_fluid_properties_object.p_from_v_e(v, e);
* // pressure now contains partial derivatives w.r.t. all degrees of freedom
* @endcode
*/
///@{
propfunc(p, v, e)
propfunc(T, v, e)
propfunc(c, v, e)
propfunc(cp, v, e)
propfunc(cv, v, e)
propfunc(mu, v, e)
propfunc(k, v, e)
propfuncWithDefault(s, v, e)
propfunc(s, h, p)
propfunc(rho, p, s)
propfunc(e, v, h)
propfuncWithDefault(s, p, T)
propfunc(pp_sat, p, T)
propfunc(mu, rho, T)
propfunc(k, rho, T)
propfuncWithDefault(c, p, T)
propfuncWithDefault(cp, p, T)
propfuncWithDefault(cv, p, T)
propfuncWithDefault(mu, p, T)
propfuncWithDefault(k, p, T)
propfunc(rho, p, T)
propfunc(e, p, rho)
propfunc(e, T, v)
propfunc(p, T, v)
propfunc(h, T, v)
propfunc(s, T, v)
propfunc(cv, T, v)
propfunc(h, p, T)
propfuncWithDefault(h, v, e)
propfunc(g, v, e)
propfuncWithDefault(p, h, s)
propfunc(T, h, p) // temporary, until uniformization
propfuncWithDefault(T, p, h)
propfuncWithDefault(beta, p, T)
propfuncWithDefault(v, p, T)
propfuncWithDefault(e, p, T)
propfuncWithDefault(gamma, v, e)
propfuncWithDefault(gamma, p, T)
///@}
// clang-format on
#undef propfunc
#undef propfuncWithDefault
#undef propfuncAD
/**
* Fluid name
* @return string representing fluid name
*/
virtual std::string fluidName() const;
/**
* Molar mass [kg/mol]
* @return molar mass
*/
virtual Real molarMass() const;
/**
* Critical pressure
* @return critical pressure (Pa)
*/
virtual Real criticalPressure() const;
/**
* Critical temperature
* @return critical temperature (K)
*/
virtual Real criticalTemperature() const;
/**
* Critical density
* @return critical density (kg/m^3)
*/
virtual Real criticalDensity() const;
/**
* Critical specific internal energy
* @return specific internal energy (J/kg)
*/
virtual Real criticalInternalEnergy() const;
/**
* Triple point pressure
* @return triple point pressure (Pa)
*/
virtual Real triplePointPressure() const;
/**
* Triple point temperature
* @return triple point temperature (K)
*/
virtual Real triplePointTemperature() const;
/**
* Specific internal energy from temperature and specific volume
*
* @param[in] T temperature
* @param[in] v specific volume
*/
virtual Real e_spndl_from_v(Real v) const;
/**
* Specific internal energy from temperature and specific volume
*
* @param[in] T temperature
* @param[in] v specific volume
*/
virtual void v_e_spndl_from_T(Real T, Real & v, Real & e) const;
/**
* Vapor pressure. Used to delineate liquid and gas phases.
* Valid for temperatures between the triple point temperature
* and the critical temperature
*
* @param T fluid temperature (K)
* @param[out] saturation pressure (Pa)
* @param[out] derivative of saturation pressure wrt temperature (Pa/K)
*/
virtual Real vaporPressure(Real T) const;
virtual void vaporPressure(Real T, Real & psat, Real & dpsat_dT) const;
virtual ADReal vaporPressure(const ADReal & T) const;
/**
* Vapor temperature. Used to delineate liquid and gas phases.
* Valid for pressures between the triple point pressure
* and the critical pressure
*
* @param p fluid pressure (Pa)
* @param[out] saturation temperature (K)
* @param[out] derivative of saturation temperature wrt pressure
*/
virtual Real vaporTemperature(Real p) const;
virtual void vaporTemperature(Real p, Real & Tsat, Real & dTsat_dp) const;
virtual ADReal vaporTemperature(const ADReal & p) const;
/**
* Henry's law coefficients for dissolution in water
* @return Henry's constant coefficients
*/
virtual std::vector<Real> henryCoefficients() const;
template <typename CppType>
void v_e_from_p_T(const CppType & p, const CppType & T, CppType & v, CppType & e) const;
template <typename CppType>
void v_e_from_p_T(const CppType & p,
const CppType & T,
CppType & v,
CppType & dv_dp,
CppType & dv_dT,
CppType & e,
CppType & de_dp,
CppType & de_dT) const;
/**
* Combined methods. These methods are particularly useful for the PorousFlow
* module, where density and viscosity are typically both computed everywhere.
* The combined methods allow the most efficient means of calculating both
* properties, especially where rho(p, T) and mu(rho, T). In this case, an
* extra density calculation would be required to calculate mu(p, T). All
* property names are described above.
*/
virtual void rho_mu_from_p_T(Real p, Real T, Real & rho, Real & mu) const;
virtual void rho_mu_from_p_T(Real p,
Real T,
Real & rho,
Real & drho_dp,
Real & drho_dT,
Real & mu,
Real & dmu_dp,
Real & dmu_dT) const;
virtual void rho_mu_from_p_T(const ADReal & p, const ADReal & T, ADReal & rho, ADReal & mu) const;
virtual void rho_e_from_p_T(Real p,
Real T,
Real & rho,
Real & drho_dp,
Real & drho_dT,
Real & e,
Real & de_dp,
Real & de_dT) const;
/**
* Determines (p,T) from (v,e) using Newton Solve in 2D
* Useful for conversion between different sets of state variables
*
* @param[in] v specific volume (m^3 / kg)
* @param[in] e specific internal energy (J / kg)
* @param[in] p0 initial guess for pressure (Pa / kg)
* @param[in] T0 initial guess for temperature (K)
* @param[out] fluid pressure (Pa / kg)
* @param[out] Temperature (K)
*/
template <typename CppType>
void p_T_from_v_e(const CppType & v,
const CppType & e,
Real p0,
Real T0,
CppType & p,
CppType & T,
bool & conversion_succeeded) const;
/**
* Determines (p,T) from (v,h) using Newton Solve in 2D
* Useful for conversion between different sets of state variables
*
* @param[in] v specific volume (m^3 / kg)
* @param[in] h specific enthalpy (J / kg)
* @param[in] p0 initial guess for pressure (Pa / kg)
* @param[in] T0 initial guess for temperature (K)
* @param[out] fluid pressure (Pa / kg)
* @param[out] Temperature (K)
*/
template <typename T>
void p_T_from_v_h(const T & v,
const T & h,
Real p0,
Real T0,
T & pressure,
T & temperature,
bool & conversion_succeeded) const;
/**
* Determines (p,T) from (h,s) using Newton Solve in 2D
* Useful for conversion between different sets of state variables
*
* @param[in] h specific enthalpy (J / kg)
* @param[in] s specific entropy (J/K*kg)
* @param[in] p0 initial guess for pressure (Pa / kg)
* @param[in] T0 initial guess for temperature (K)
* @param[out] fluid pressure (Pa / kg)
* @param[out] Temperature (K)
*/
template <typename T>
void p_T_from_h_s(const T & h,
const T & s,
Real p0,
Real T0,
T & pressure,
T & temperature,
bool & conversion_succeeded) const;
protected:
/**
* Computes the dependent variable z and its derivatives with respect to the independent
* variables x and y using the simple two parameter \p z_from_x_y functor. The derivatives are
* computed using a compound automatic differentiation type
*/
template <typename T, typename Functor>
static void
xyDerivatives(const T x, const T & y, T & z, T & dz_dx, T & dz_dy, const Functor & z_from_x_y);
/**
* Given a type example, this method returns zero and unity representations of that type (first
* and second members of returned pair respectively)
*/
template <typename T>
static std::pair<T, T> makeZeroAndOne(const T &);
/**
* Newton's method may be used to convert between variable sets
* _tolerance, _T_initial_guess, and _p_initial_guess are the parameters for these
* iterative solves
*/
const Real _tolerance;
const Real _T_initial_guess;
const Real _p_initial_guess;
private:
void unimplementedDerivativeMethod(const std::string & property_function_name) const
{
const std::string message =
"The fluid properties class '" + type() +
"' has not implemented the method below, which computes derivatives of fluid properties "
"with regards to the flow variables. If your application requires this "
"method, you must either implement it or use a different fluid properties "
" class.\n\n" +
property_function_name;
if (_allow_imperfect_jacobians)
mooseDoOnce(mooseWarning(message + "\nThe unimplemented derivatives for this fluid property "
"are currently neglected, set to 0."));
else
mooseError(message + "\n\nYou can avoid this error by neglecting the "
"unimplemented derivatives of fluid properties by setting the "
"'allow_imperfect_jacobians' parameter");
}
};
#pragma GCC diagnostic pop
template <typename T>
std::pair<T, T>
SinglePhaseFluidProperties::makeZeroAndOne(const T & /*ex*/)
{
return {T{0, 0}, T{1, 0}};
}
template <>
inline std::pair<Real, Real>
SinglePhaseFluidProperties::makeZeroAndOne(const Real & /*ex*/)
{
return {Real{0}, Real{1}};
}
template <typename T, typename Functor>
void
SinglePhaseFluidProperties::xyDerivatives(
const T x, const T & y, T & z, T & dz_dx, T & dz_dy, const Functor & z_from_x_y)
{
typedef MetaPhysicL::DualNumber<T, MetaPhysicL::NumberArray<2, T>> CompoundType;
const auto [zero, one] = makeZeroAndOne(x);
CompoundType x_c(x, zero);
auto & x_cd = x_c.derivatives();
x_cd[0] = one;
CompoundType y_c(y, zero);
auto & y_cd = y_c.derivatives();
y_cd[1] = one;
const auto z_c = z_from_x_y(x_c, y_c);
z = z_c.value();
dz_dx = z_c.derivatives()[0];
dz_dy = z_c.derivatives()[1];
}
template <typename CppType>
void
SinglePhaseFluidProperties::p_T_from_v_e(const CppType & v, // v value
const CppType & e, // e value
const Real p0, // initial guess
const Real T0, // initial guess
CppType & p, // returned pressure
CppType & T, // returned temperature
bool & conversion_succeeded) const
{
auto v_lambda = [&](const CppType & pressure,
const CppType & temperature,
CppType & new_v,
CppType & dv_dp,
CppType & dv_dT) { v_from_p_T(pressure, temperature, new_v, dv_dp, dv_dT); };
auto e_lambda = [&](const CppType & pressure,
const CppType & temperature,
CppType & new_e,
CppType & de_dp,
CppType & de_dT) { e_from_p_T(pressure, temperature, new_e, de_dp, de_dT); };
try
{
FluidPropertiesUtils::NewtonSolve2D(
v, e, p0, T0, p, T, _tolerance, _tolerance, v_lambda, e_lambda);
conversion_succeeded = true;
}
catch (MooseException &)
{
conversion_succeeded = false;
}
if (!conversion_succeeded)
mooseDoOnce(mooseWarning("Conversion from (v, e)=(", v, ", ", e, ") to (p, T) failed"));
}
template <typename T>
void
SinglePhaseFluidProperties::p_T_from_v_h(const T & v, // v value
const T & h, // e value
const Real p0, // initial guess
const Real T0, // initial guess
T & pressure, // returned pressure
T & temperature, // returned temperature
bool & conversion_succeeded) const
{
auto v_lambda = [&](const T & pressure, const T & temperature, T & new_v, T & dv_dp, T & dv_dT)
{ v_from_p_T(pressure, temperature, new_v, dv_dp, dv_dT); };
auto h_lambda = [&](const T & pressure, const T & temperature, T & new_h, T & dh_dp, T & dh_dT)
{ h_from_p_T(pressure, temperature, new_h, dh_dp, dh_dT); };
try
{
FluidPropertiesUtils::NewtonSolve2D(
v, h, p0, T0, pressure, temperature, _tolerance, _tolerance, v_lambda, h_lambda);
conversion_succeeded = true;
}
catch (MooseException &)
{
conversion_succeeded = false;
}
if (!conversion_succeeded)
mooseDoOnce(mooseWarning("Conversion from (v, h)=(", v, ", ", h, ") to (p, T) failed"));
}
template <typename T>
void
SinglePhaseFluidProperties::p_T_from_h_s(const T & h, // h value
const T & s, // s value
const Real p0, // initial guess
const Real T0, // initial guess
T & pressure, // returned pressure
T & temperature, // returned temperature
bool & conversion_succeeded) const
{
auto h_lambda = [&](const T & pressure, const T & temperature, T & new_h, T & dh_dp, T & dh_dT)
{ h_from_p_T(pressure, temperature, new_h, dh_dp, dh_dT); };
auto s_lambda = [&](const T & pressure, const T & temperature, T & new_s, T & ds_dp, T & ds_dT)
{ s_from_p_T(pressure, temperature, new_s, ds_dp, ds_dT); };
try
{
FluidPropertiesUtils::NewtonSolve2D(
h, s, p0, T0, pressure, temperature, _tolerance, _tolerance, h_lambda, s_lambda);
conversion_succeeded = true;
}
catch (MooseException &)
{
conversion_succeeded = false;
}
if (!conversion_succeeded)
mooseDoOnce(mooseWarning("Conversion from (h, s)=(", h, ", ", s, ") to (p, T) failed"));
}
template <typename CppType>
void
SinglePhaseFluidProperties::v_e_from_p_T(const CppType & p,
const CppType & T,
CppType & v,
CppType & e) const
{
const CppType rho = rho_from_p_T(p, T);
v = 1.0 / rho;
try
{
// more likely to not involve a Newton search
e = e_from_p_T(p, T);
}
catch (...)
{
e = e_from_p_rho(p, rho);
}
}
template <typename CppType>
void
SinglePhaseFluidProperties::v_e_from_p_T(const CppType & p,
const CppType & T,
CppType & v,
CppType & dv_dp,
CppType & dv_dT,
CppType & e,
CppType & de_dp,
CppType & de_dT) const
{
CppType rho, drho_dp, drho_dT;
rho_from_p_T(p, T, rho, drho_dp, drho_dT);
v = 1.0 / rho;
const CppType dv_drho = -1.0 / (rho * rho);
dv_dp = dv_drho * drho_dp;
dv_dT = dv_drho * drho_dT;
CppType de_dp_partial, de_drho;
e_from_p_rho(p, rho, e, de_dp_partial, de_drho);
de_dp = de_dp_partial + de_drho * drho_dp;
de_dT = de_drho * drho_dT;
}
(moose/modules/fluid_properties/src/fluidproperties/SinglePhaseFluidProperties.C)
// This file is part of the MOOSE framework
// https://www.mooseframework.org
//
// All rights reserved, see COPYRIGHT for full restrictions
// https://github.com/idaholab/moose/blob/master/COPYRIGHT
//
// Licensed under LGPL 2.1, please see LICENSE for details
// https://www.gnu.org/licenses/lgpl-2.1.html
#include "SinglePhaseFluidProperties.h"
InputParameters
SinglePhaseFluidProperties::validParams()
{
InputParameters params = FluidProperties::validParams();
params.set<std::string>("fp_type") = "single-phase-fp";
// Variable set conversion parameters
params.addRangeCheckedParam<Real>(
"tolerance", 1e-8, "tolerance > 0", "Tolerance for 2D Newton variable set conversion");
params.addRangeCheckedParam<Real>(
"T_initial_guess",
400,
"T_initial_guess > 0",
"Temperature initial guess for Newton Method variable set conversion");
params.addRangeCheckedParam<Real>(
"p_initial_guess",
2e5,
"p_initial_guess > 0",
"Pressure initial guess for Newton Method variable set conversion");
params.addParamNamesToGroup("tolerance T_initial_guess p_initial_guess",
"Variable set conversions Newton solve");
return params;
}
SinglePhaseFluidProperties::SinglePhaseFluidProperties(const InputParameters & parameters)
: FluidProperties(parameters),
// downstream apps are creating fluid properties without their parameters, hence the workaround
_tolerance(isParamValid("tolerance") ? getParam<Real>("tolerance") : 1e-8),
_T_initial_guess(isParamValid("T_initial_guess") ? getParam<Real>("T_initial_guess") : 400),
_p_initial_guess(isParamValid("p_initial_guess") ? getParam<Real>("p_initial_guess") : 2e5)
{
}
SinglePhaseFluidProperties::~SinglePhaseFluidProperties() {}
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Woverloaded-virtual"
Real
SinglePhaseFluidProperties::s_from_p_T(const Real pressure, const Real temperature) const
{
Real v, e;
v_e_from_p_T(pressure, temperature, v, e);
return s_from_v_e(v, e);
}
void
SinglePhaseFluidProperties::s_from_p_T(
const Real pressure, const Real temperature, Real & s, Real & ds_dp, Real & ds_dT) const
{
Real v, e, dv_dp, dv_dT, de_dp, de_dT;
v_e_from_p_T(pressure, temperature, v, dv_dp, dv_dT, e, de_dp, de_dT);
Real ds_dv, ds_de;
s_from_v_e(v, e, s, ds_dv, ds_de);
ds_dp = ds_dv * dv_dp + ds_de * de_dp;
ds_dT = ds_dv * dv_dT + ds_de * de_dT;
}
Real
SinglePhaseFluidProperties::s_from_v_e(const Real v, const Real e) const
{
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);
const Real s = s_from_p_T(p, T);
return s;
}
void
SinglePhaseFluidProperties::s_from_v_e(
const Real v, const Real e, Real & s, Real & ds_dv, Real & ds_de) const
{
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);
s = s_from_p_T(p, T);
ds_dv = p / T;
ds_de = 1 / T;
}
Real
SinglePhaseFluidProperties::c_from_p_T(Real p, Real T) const
{
Real v, e;
v_e_from_p_T(p, T, v, e);
return c_from_v_e(v, e);
}
void
SinglePhaseFluidProperties::c_from_p_T(Real p, Real T, Real & c, Real & dc_dp, Real & dc_dT) const
{
Real v, e, dv_dp, dv_dT, de_dp, de_dT;
v_e_from_p_T(p, T, v, dv_dp, dv_dT, e, de_dp, de_dT);
Real dc_dv, dc_de;
c_from_v_e(v, e, c, dc_dv, dc_de);
dc_dp = dc_dv * dv_dp + dc_de * de_dp;
dc_dT = dc_dv * dv_dT + dc_de * de_dT;
}
Real
SinglePhaseFluidProperties::mu_from_p_T(Real p, Real T) const
{
Real v, e;
v_e_from_p_T(p, T, v, e);
return mu_from_v_e(v, e);
}
void
SinglePhaseFluidProperties::mu_from_p_T(
Real p, Real T, Real & mu, Real & dmu_dp, Real & dmu_dT) const
{
Real v, e, dv_dp, dv_dT, de_dp, de_dT;
v_e_from_p_T(p, T, v, dv_dp, dv_dT, e, de_dp, de_dT);
Real dmu_dv, dmu_de;
mu_from_v_e(v, e, mu, dmu_dv, dmu_de);
dmu_dp = dmu_dv * dv_dp + dmu_de * de_dp;
dmu_dT = dmu_dv * dv_dT + dmu_de * de_dT;
}
Real
SinglePhaseFluidProperties::cv_from_p_T(Real p, Real T) const
{
Real v, e;
v_e_from_p_T(p, T, v, e);
return cv_from_v_e(v, e);
}
void
SinglePhaseFluidProperties::cv_from_p_T(
Real p, Real T, Real & cv, Real & dcv_dp, Real & dcv_dT) const
{
Real v, e, dv_dp, dv_dT, de_dp, de_dT;
v_e_from_p_T(p, T, v, dv_dp, dv_dT, e, de_dp, de_dT);
Real dcv_dv, dcv_de;
cv_from_v_e(v, e, cv, dcv_dv, dcv_de);
dcv_dp = dcv_dv * dv_dp + dcv_de * de_dp;
dcv_dT = dcv_dv * dv_dT + dcv_de * de_dT;
}
Real
SinglePhaseFluidProperties::cp_from_p_T(Real p, Real T) const
{
Real v, e;
v_e_from_p_T(p, T, v, e);
return cp_from_v_e(v, e);
}
void
SinglePhaseFluidProperties::cp_from_p_T(
Real p, Real T, Real & cp, Real & dcp_dp, Real & dcp_dT) const
{
Real v, e, dv_dp, dv_dT, de_dp, de_dT;
v_e_from_p_T(p, T, v, dv_dp, dv_dT, e, de_dp, de_dT);
Real dcp_dv, dcp_de;
cp_from_v_e(v, e, cp, dcp_dv, dcp_de);
dcp_dp = dcp_dv * dv_dp + dcp_de * de_dp;
dcp_dT = dcp_dv * dv_dT + dcp_de * de_dT;
}
Real
SinglePhaseFluidProperties::k_from_p_T(Real p, Real T) const
{
Real v, e;
v_e_from_p_T(p, T, v, e);
return k_from_v_e(v, e);
}
void
SinglePhaseFluidProperties::k_from_p_T(Real p, Real T, Real & k, Real & dk_dp, Real & dk_dT) const
{
Real v, e, dv_dp, dv_dT, de_dp, de_dT;
v_e_from_p_T(p, T, v, dv_dp, dv_dT, e, de_dp, de_dT);
Real dk_dv, dk_de;
k_from_v_e(v, e, k, dk_dv, dk_de);
dk_dp = dk_dv * dv_dp + dk_de * de_dp;
dk_dT = dk_dv * dv_dT + dk_de * de_dT;
}
Real
SinglePhaseFluidProperties::h_from_v_e(Real v, Real e) const
{
return e + v * p_from_v_e(v, e);
}
void
SinglePhaseFluidProperties::h_from_v_e(Real v, Real e, Real & h, Real & dh_dv, Real & dh_de) const
{
Real p, dp_dv, dp_de;
p_from_v_e(v, e, p, dp_dv, dp_de);
h = e + v * p;
dh_dv = p + v * dp_dv;
dh_de = 1 + v * dp_de;
}
Real
SinglePhaseFluidProperties::e_from_p_T(Real p, Real T) const
{
const Real rho = rho_from_p_T(p, T);
return e_from_p_rho(p, rho);
}
void
SinglePhaseFluidProperties::e_from_p_T(Real p, Real T, Real & e, Real & de_dp, Real & de_dT) const
{
// From rho(p,T), compute: drho(p,T)/dp, drho(p,T)/dT
Real rho = 0., drho_dp = 0., drho_dT = 0.;
rho_from_p_T(p, T, rho, drho_dp, drho_dT);
// From e(p, rho), compute: de(p,rho)/dp, de(p,rho)/drho
Real depr_dp = 0., depr_drho = 0.;
e_from_p_rho(p, rho, e, depr_dp, depr_drho);
// Using partial derivative rules, we have:
// de(p,T)/dp = de(p,rho)/dp * dp/dp + de(p,rho)/drho * drho(p,T)/dp, (dp/dp == 1)
// de(p,T)/dT = de(p,rho)/dp * dp/dT + de(p,rho)/drho * drho(p,T)/dT, (dp/dT == 0)
de_dp = depr_dp + depr_drho * drho_dp;
de_dT = depr_drho * drho_dT;
}
Real
SinglePhaseFluidProperties::v_from_p_T(Real p, Real T) const
{
const Real rho = rho_from_p_T(p, T);
return 1.0 / rho;
}
void
SinglePhaseFluidProperties::v_from_p_T(Real p, Real T, Real & v, Real & dv_dp, Real & dv_dT) const
{
Real rho, drho_dp, drho_dT;
rho_from_p_T(p, T, rho, drho_dp, drho_dT);
v = 1.0 / rho;
const Real dv_drho = -1.0 / (rho * rho);
dv_dp = dv_drho * drho_dp;
dv_dT = dv_drho * drho_dT;
}
void
SinglePhaseFluidProperties::beta_from_p_T(Real, Real, Real &, Real &, Real &) const
{
mooseError(__PRETTY_FUNCTION__, " is not implemented.");
}
Real
SinglePhaseFluidProperties::beta_from_p_T(Real p, Real T) const
{
// The volumetric thermal expansion coefficient is defined as
// 1/v dv/dT)_p
// It is the fractional change rate of volume with respect to temperature change
// at constant pressure. Here it is coded as
// - 1/rho drho/dT)_p
// using chain rule with v = v(rho)
Real rho, drho_dp, drho_dT;
rho_from_p_T(p, T, rho, drho_dp, drho_dT);
return -drho_dT / rho;
}
Real
SinglePhaseFluidProperties::molarMass() const
{
mooseError(__PRETTY_FUNCTION__, " not implemented.");
}
std::string
SinglePhaseFluidProperties::fluidName() const
{
return std::string("");
}
Real
SinglePhaseFluidProperties::criticalPressure() const
{
mooseError(__PRETTY_FUNCTION__, " not implemented.");
}
Real
SinglePhaseFluidProperties::criticalTemperature() const
{
mooseError(__PRETTY_FUNCTION__, " not implemented.");
}
Real
SinglePhaseFluidProperties::criticalDensity() const
{
mooseError(__PRETTY_FUNCTION__, " not implemented.");
}
Real
SinglePhaseFluidProperties::criticalInternalEnergy() const
{
mooseError(__PRETTY_FUNCTION__, " not implemented.");
}
Real
SinglePhaseFluidProperties::triplePointPressure() const
{
mooseError(__PRETTY_FUNCTION__, " not implemented.");
}
Real
SinglePhaseFluidProperties::triplePointTemperature() const
{
mooseError(__PRETTY_FUNCTION__, " not implemented.");
}
Real
SinglePhaseFluidProperties::gamma_from_v_e(Real v, Real e) const
{
return cp_from_v_e(v, e) / cv_from_v_e(v, e);
}
void
SinglePhaseFluidProperties::gamma_from_v_e(
Real v, Real e, Real & gamma, Real & dgamma_dv, Real & dgamma_de) const
{
unimplementedDerivativeMethod(__PRETTY_FUNCTION__);
dgamma_dv = 0.0;
dgamma_de = 0.0;
gamma = gamma_from_v_e(v, e);
}
Real
SinglePhaseFluidProperties::gamma_from_p_T(Real p, Real T) const
{
return cp_from_p_T(p, T) / cv_from_p_T(p, T);
}
void
SinglePhaseFluidProperties::gamma_from_p_T(
Real p, Real T, Real & gamma, Real & dgamma_dp, Real & dgamma_dT) const
{
unimplementedDerivativeMethod(__PRETTY_FUNCTION__);
dgamma_dp = 0.0;
dgamma_dT = 0.0;
gamma = gamma_from_p_T(p, T);
}
Real
SinglePhaseFluidProperties::vaporPressure(Real) const
{
mooseError(__PRETTY_FUNCTION__, " not implemented.");
}
std::vector<Real>
SinglePhaseFluidProperties::henryCoefficients() const
{
mooseError(__PRETTY_FUNCTION__, " not implemented.");
}
void
SinglePhaseFluidProperties::vaporPressure(Real T, Real & p, Real & dp_dT) const
{
unimplementedDerivativeMethod(__PRETTY_FUNCTION__);
dp_dT = 0.0;
p = vaporPressure(T);
}
ADReal
SinglePhaseFluidProperties::vaporPressure(const ADReal & T) const
{
Real p = 0.0;
Real temperature = T.value();
Real dpdT = 0.0;
vaporPressure(temperature, p, dpdT);
ADReal result = p;
result.derivatives() = T.derivatives() * dpdT;
return result;
}
Real
SinglePhaseFluidProperties::vaporTemperature(Real) const
{
mooseError(__PRETTY_FUNCTION__, " not implemented.");
}
void
SinglePhaseFluidProperties::vaporTemperature(Real p, Real & T, Real & dT_dp) const
{
unimplementedDerivativeMethod(__PRETTY_FUNCTION__);
dT_dp = 0.0;
T = vaporTemperature(p);
}
ADReal
SinglePhaseFluidProperties::vaporTemperature(const ADReal & p) const
{
Real T = 0.0;
Real pressure = p.value();
Real dTdp = 0.0;
vaporTemperature(pressure, T, dTdp);
ADReal result = T;
result.derivatives() = p.derivatives() * dTdp;
return result;
}
void
SinglePhaseFluidProperties::rho_e_from_p_T(Real p,
Real T,
Real & rho,
Real & drho_dp,
Real & drho_dT,
Real & e,
Real & de_dp,
Real & de_dT) const
{
rho_from_p_T(p, T, rho, drho_dp, drho_dT);
e_from_p_T(p, T, e, de_dp, de_dT);
}
void
SinglePhaseFluidProperties::rho_mu_from_p_T(Real p, Real T, Real & rho, Real & mu) const
{
rho = rho_from_p_T(p, T);
mu = mu_from_p_T(p, T);
}
void
SinglePhaseFluidProperties::rho_mu_from_p_T(Real p,
Real T,
Real & rho,
Real & drho_dp,
Real & drho_dT,
Real & mu,
Real & dmu_dp,
Real & dmu_dT) const
{
rho_from_p_T(p, T, rho, drho_dp, drho_dT);
mu_from_p_T(p, T, mu, dmu_dp, dmu_dT);
}
void
SinglePhaseFluidProperties::rho_mu_from_p_T(const ADReal & p,
const ADReal & T,
ADReal & rho,
ADReal & mu) const
{
rho = rho_from_p_T(p, T);
mu = mu_from_p_T(p, T);
}
Real
SinglePhaseFluidProperties::e_spndl_from_v(Real) const
{
mooseError(__PRETTY_FUNCTION__, " not implemented.");
}
void
SinglePhaseFluidProperties::v_e_spndl_from_T(Real, Real &, Real &) const
{
mooseError(__PRETTY_FUNCTION__, " not implemented.");
}
Real
SinglePhaseFluidProperties::T_from_p_h(Real p, Real h) const
{
const Real s = s_from_h_p(h, p);
const Real rho = rho_from_p_s(p, s);
const Real v = 1. / rho;
const Real e = e_from_v_h(v, h);
return T_from_v_e(v, e);
}
Real
SinglePhaseFluidProperties::p_from_h_s(Real h, Real s) const
{
Real p0 = _p_initial_guess;
Real T0 = _T_initial_guess;
Real p, T;
bool conversion_succeeded = true;
p_T_from_h_s(h, s, p0, T0, p, T, conversion_succeeded);
return p;
}
void
SinglePhaseFluidProperties::p_from_h_s(Real h, Real s, Real & p, Real & dp_dh, Real & dp_ds) const
{
Real p0 = _p_initial_guess;
Real T0 = _T_initial_guess;
Real T;
bool conversion_succeeded = true;
p_T_from_h_s(h, s, p0, T0, p, T, conversion_succeeded);
dp_dh = rho_from_p_T(p, T);
dp_ds = -T * rho_from_p_T(p, T);
}
void
SinglePhaseFluidProperties::T_from_p_h(Real p, Real h, Real & T, Real & dT_dp, Real & dT_dh) const
{
Real s, ds_dh, ds_dp;
s_from_h_p(h, p, s, ds_dh, ds_dp);
Real rho, drho_dp_partial, drho_ds;
rho_from_p_s(p, s, rho, drho_dp_partial, drho_ds);
const Real drho_dp = drho_dp_partial + drho_ds * ds_dp;
const Real drho_dh = drho_ds * ds_dh;
const Real v = 1.0 / rho;
const Real dv_drho = -1.0 / (rho * rho);
const Real dv_dp = dv_drho * drho_dp;
const Real dv_dh = dv_drho * drho_dh;
Real e, de_dv, de_dh_partial;
e_from_v_h(v, h, e, de_dv, de_dh_partial);
const Real de_dp = de_dv * dv_dp;
const Real de_dh = de_dh_partial + de_dv * dv_dh;
Real dT_dv, dT_de;
T_from_v_e(v, e, T, dT_dv, dT_de);
dT_dp = dT_dv * dv_dp + dT_de * de_dp;
dT_dh = dT_dv * dv_dh + dT_de * de_dh;
}
#pragma GCC diagnostic pop