Add an option for a functionalized cp for TemperaturePressureFunctionFP

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
 

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;
 };

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;