Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add option to advect enthalpy instead of cp T for non constant cp fluids #29569

Open
wants to merge 8 commits into
base: next
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from 7 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Original file line number Diff line number Diff line change
Expand Up @@ -69,7 +69,8 @@ class NaClFluidProperties : public SinglePhaseFluidProperties

virtual Real cp_from_p_T(Real pressure, Real temperature) const override;

using SinglePhaseFluidProperties::cp_from_p_T;
virtual void cp_from_p_T(
Real pressure, Real temperature, Real & cp, Real & dcp_dp, Real & dcp_dT) const override;

virtual Real cv_from_p_T(Real pressure, Real temperature) const override;

Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -79,7 +79,7 @@ NaClFluidProperties::rho_from_p_T(Real pressure, Real temperature) const
{
// Correlation needs pressure in bar
Real pbar = pressure * 1.0e-5;
// Correlation requires temperature in Celcius
// Correlation requires temperature in Celsius
Real Tc = temperature - _T_c2k;

// Halite density at 0 Pa
Expand All @@ -99,7 +99,7 @@ NaClFluidProperties::rho_from_p_T(

// Correlation needs pressure in bar
Real pbar = pressure * 1.0e-5;
// Correlation requires temperature in Celcius
// Correlation requires temperature in Celsius
Real Tc = temperature - _T_c2k;

// Halite density at 0 Pa
Expand Down Expand Up @@ -137,19 +137,42 @@ NaClFluidProperties::cp_from_p_T(Real pressure, Real temperature) const
{
// Correlation needs pressure in bar
Real pbar = pressure * 10.0e-5;
// Correlation requires temperature in Celcius
// Correlation requires temperature in Celsius
Real Tc = temperature - _T_c2k;
// Triple point temperature of NaCl (in C)
Real Tt = _T_triple - _T_c2k;
// Coefficients used in the correlation
Real r3 = -1.7099e-3 - 3.82734e-6 * Tc - 8.65455e-9 * Tc * Tc;
Real r4 = 5.29063e-8 - 9.63084e-11 * Tc + 6.50745e-13 * Tc * Tc;

// Halite isobaric heat capapcity
// Halite isobaric heat capacity
return 1148.81 + 0.551548 * (Tc - Tt) + 2.64309e-4 * (Tc - Tt) * (Tc - Tt) + r3 * pbar +
r4 * pbar * pbar;
}

void
NaClFluidProperties::cp_from_p_T(
Real pressure, Real temperature, Real & cp, Real & dcp_dp, Real & dcp_dT) const
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could we please add unit test to compare the computed cp against tabulated values for a few points in the input sapce?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

they already exist, the dependence on pressure is just so small it did not catch it

{
// Correlation needs pressure in bar
Real pbar = pressure * 10.0e-5;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Conversion to bar should be 1e-5 (1 bar = 1e5 Pa)

Copy link
Contributor Author

@GiudGiud GiudGiud Jan 12, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

very true. this is like this in the original property :O ( I copy pasted that line)

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

fixed

// Correlation requires temperature in Celsius
Real Tc = temperature - _T_c2k;
// Triple point temperature of NaCl (in C)
Real Tt = _T_triple - _T_c2k;
// Coefficients used in the correlation
Real r3 = -1.7099e-3 - 3.82734e-6 * Tc - 8.65455e-9 * Tc * Tc;
Real r4 = 5.29063e-8 - 9.63084e-11 * Tc + 6.50745e-13 * Tc * Tc;
Real dr3_dT = -3.82734e-6 - 2 * 8.65455e-9 * Tc;
Real dr4_dT = -9.63084e-11 + 2 * 6.50745e-13 * Tc;

// Halite isobaric heat capacity
cp = 1148.81 + 0.551548 * (Tc - Tt) + 2.64309e-4 * (Tc - Tt) * (Tc - Tt) + r3 * pbar +
r4 * pbar * pbar;
dcp_dp = r3 * 10.e-5 + 2 * r4 * pbar * 10.e-5;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Conversion to bar may also be wrong here

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

fixed

dcp_dT = 0.551548 + 2 * 2.64309e-4 * (Tc - Tt) + dr3_dT * pbar + dr4_dT * pbar * pbar;
}

Real
NaClFluidProperties::cv_from_p_T(Real pressure, Real temperature) const
{
Expand All @@ -159,7 +182,7 @@ NaClFluidProperties::cv_from_p_T(Real pressure, Real temperature) const
Real
NaClFluidProperties::k_from_p_T(Real /*pressure*/, Real temperature) const
{
// Correlation requires temperature in Celcius
// Correlation requires temperature in Celsius
Real Tc = temperature - _T_c2k;

return 6.82793 - 3.16584e-2 * Tc + 1.03451e-4 * Tc * Tc - 1.48207e-7 * Tc * Tc * Tc;
Expand All @@ -169,7 +192,7 @@ void
NaClFluidProperties::k_from_p_T(
Real /*pressure*/, Real temperature, Real & k, Real & dk_dp, Real & dk_dT) const
{
// Correlation requires temperature in Celcius
// Correlation requires temperature in Celsius
Real Tc = temperature - _T_c2k;

k = 6.82793 - 3.16584e-2 * Tc + 1.03451e-4 * Tc * Tc - 1.48207e-7 * Tc * Tc * Tc;
Expand All @@ -182,7 +205,7 @@ NaClFluidProperties::h_from_p_T(Real pressure, Real temperature) const
{
// Correlation needs pressure in bar
Real pbar = pressure * 1.0e-5;
// Correlation requires temperature in Celcius
// Correlation requires temperature in Celsius
Real Tc = temperature - _T_c2k;
// Triple point temperature of water (in C)
Real Tt = 273.16 - _T_c2k;
Expand All @@ -200,7 +223,7 @@ NaClFluidProperties::h_from_p_T(
{
// Correlation needs pressure in bar
Real pbar = pressure * 1.0e-5;
// Correlation requires temperature in Celcius
// Correlation requires temperature in Celsius
Real Tc = temperature - _T_c2k;
// Triple point temperature of water (in C)
Real Tt = 273.16 - _T_c2k;
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -104,6 +104,7 @@ TEST_F(NaClFluidPropertiesTest, 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, tol)
}

/**
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,8 @@

#include "FunctorMaterial.h"

class SinglePhaseFluidProperties;

/**
* This is the material class used to compute enthalpy for the incompressible/weakly-compressible
* finite-volume implementation of the Navier-Stokes equations
Expand All @@ -23,12 +25,24 @@ class INSFVEnthalpyFunctorMaterial : public FunctorMaterial
INSFVEnthalpyFunctorMaterial(const InputParameters & parameters);

protected:
/// whether we can use a constant cp as a shortcut to compute enthalpy
bool _assume_constant_cp;
GiudGiud marked this conversation as resolved.
Show resolved Hide resolved

/// A fluid properties user object to compute enthalpy
const SinglePhaseFluidProperties * _fp;

/// density
const Moose::Functor<ADReal> & _rho;

/// the temperature
const Moose::Functor<ADReal> & _temperature;

/// the pressure
const Moose::Functor<ADReal> * _pressure;

/// the specific heat capacity
const Moose::Functor<ADReal> & _cp;

/// the specific enthalpy
const Moose::Functor<ADReal> * _h;
};
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,7 @@

#include "INSFVEnthalpyFunctorMaterial.h"
#include "NS.h"
#include "SinglePhaseFluidProperties.h"

registerMooseObjectRenamed("NavierStokesApp",
INSFVEnthalpyMaterial,
Expand All @@ -23,7 +24,7 @@ INSFVEnthalpyFunctorMaterial::validParams()
params.addClassDescription(
"This is the material class used to compute enthalpy for the "
"incompressible/weakly-compressible finite-volume implementation of the Navier-Stokes "
"equations. Note that this class assumes that cp is a constant");
"equations.");
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

With these changes this will look similar to @freiler 's functor material. Do you see value in having only one? Or should we keep two? I think if it can be done cleanly, we should have only one.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I havent been keeping track with his PR but I think his functor material might be focusing on computing T from h rather than h from T

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If the class assumes that cp is constant, please add it to the name of the functormaterial, e.g., INSFVConstantCpEnthalpyFunctorMaterial or something. Otherwise, this is can lead to errors

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

the assumptions depend on the parameter

params.addRequiredParam<MooseFunctorName>(NS::density, "The value for the density");
params.addRequiredParam<MooseFunctorName>("temperature", "the temperature");
params.addParam<MooseFunctorName>(
Expand All @@ -32,31 +33,92 @@ INSFVEnthalpyFunctorMaterial::validParams()
NS::enthalpy_density, NS::enthalpy_density, "the name of the (extensive) enthalpy");
params.addParam<MooseFunctorName>(
NS::specific_enthalpy, NS::specific_enthalpy, "the name of the specific enthalpy");

// To handle non constant cp
params.addParam<bool>("assume_constant_cp", true, "Whether to assume cp is constant");
GiudGiud marked this conversation as resolved.
Show resolved Hide resolved
params.addParam<UserObjectName>(
NS::fluid, "Fluid properties, to be used when cp is not constant to compute enthalpy");
params.addParam<MooseFunctorName>(
NS::pressure, "Pressure functor, to be used when cp is not constant to compute enthalpy");
params.addParam<MooseFunctorName>(
NS::specific_enthalpy + "_in",
"Specific enthalpy functor, to be used when cp is not constant to compute the enthalpy, as "
"an alternative to using a 'fp' FluidProperties object");

return params;
}

INSFVEnthalpyFunctorMaterial::INSFVEnthalpyFunctorMaterial(const InputParameters & parameters)
: FunctorMaterial(parameters),
_assume_constant_cp(getParam<bool>("assume_constant_cp")),
GiudGiud marked this conversation as resolved.
Show resolved Hide resolved
_fp(isParamValid(NS::fluid)
? &UserObjectInterface::getUserObject<SinglePhaseFluidProperties>(NS::fluid)
: nullptr),
_rho(getFunctor<ADReal>(NS::density)),
_temperature(getFunctor<ADReal>("temperature")),
_cp(getFunctor<ADReal>(NS::cp))
_pressure(isParamValid("pressure") ? &getFunctor<ADReal>("pressure") : nullptr),
_cp(getFunctor<ADReal>(NS::cp)),
_h(isParamValid(NS::specific_enthalpy + "_in")
? &getFunctor<ADReal>(NS::specific_enthalpy + "_in")
: nullptr)
{
const auto & rho_h =
addFunctorProperty<ADReal>(NS::enthalpy_density,
[this](const auto & r, const auto & t)
{ return _rho(r, t) * _cp(r, t) * _temperature(r, t); });
// We have to use a warning because fp is often in the global parameters
if (_assume_constant_cp && _fp)
GiudGiud marked this conversation as resolved.
Show resolved Hide resolved
paramWarning(
"fp", "No need to specify fluid properties if assuming the specific enthalpy is constant");
if (!_assume_constant_cp && ((!_fp || !_pressure) && !_h))
paramError("fp",
"Must specify both fluid properties and pressure or an enthalpy functor if not "
"assuming the specific enthalpy is constant");

if (_assume_constant_cp)
GiudGiud marked this conversation as resolved.
Show resolved Hide resolved
{
const auto & rho_h =
addFunctorProperty<ADReal>(NS::enthalpy_density,
[this](const auto & r, const auto & t)
{ return _rho(r, t) * _cp(r, t) * _temperature(r, t); });

const auto & h = addFunctorProperty<ADReal>(NS::specific_enthalpy,
[this](const auto & r, const auto & t)
{ return _cp(r, t) * _temperature(r, t); });

addFunctorProperty<ADReal>(NS::time_deriv(getParam<MooseFunctorName>(NS::specific_enthalpy)),
[this](const auto & r, const auto & t)
{ return _cp(r, t) * _temperature.dot(r, t); });

addFunctorProperty<ADReal>(
"rho_cp_temp", [&rho_h](const auto & r, const auto & t) -> ADReal { return rho_h(r, t); });

const auto & h = addFunctorProperty<ADReal>(NS::specific_enthalpy,
[this](const auto & r, const auto & t)
{ return _cp(r, t) * _temperature(r, t); });
addFunctorProperty<ADReal>("cp_temp",
[&h](const auto & r, const auto & t) -> ADReal { return h(r, t); });
}
else if (_h)
{
addFunctorProperty<ADReal>(NS::enthalpy_density,
[this](const auto & r, const auto & t)
{ return _rho(r, t) * (*_h)(r, t); });

addFunctorProperty<ADReal>(NS::time_deriv(getParam<MooseFunctorName>(NS::specific_enthalpy)),
[this](const auto & r, const auto & t)
{ return _cp(r, t) * _temperature.dot(r, t); });
addFunctorProperty<ADReal>(NS::time_deriv(getParam<MooseFunctorName>(NS::specific_enthalpy)),
[this](const auto & r, const auto & t) { return _h->dot(r, t); });
}
else
{
addFunctorProperty<ADReal>(
NS::enthalpy_density,
[this](const auto & r, const auto & t)
{ return _rho(r, t) * _fp->h_from_p_T((*_pressure)(r, t), _temperature(r, t)); });

addFunctorProperty<ADReal>(
"rho_cp_temp", [&rho_h](const auto & r, const auto & t) -> ADReal { return rho_h(r, t); });
addFunctorProperty<ADReal>(NS::specific_enthalpy,
[this](const auto & r, const auto & t)
{ return _fp->h_from_p_T((*_pressure)(r, t), _temperature(r, t)); });

addFunctorProperty<ADReal>("cp_temp",
[&h](const auto & r, const auto & t) -> ADReal { return h(r, t); });
addFunctorProperty<ADReal>(
NS::time_deriv(getParam<MooseFunctorName>(NS::specific_enthalpy)),
[this](const auto & r, const auto & t)
{
Real h, dh_dp, dh_dT;
_fp->h_from_p_T((*_pressure)(r, t).value(), _temperature(r, t).value(), h, dh_dp, dh_dT);
return dh_dT * _temperature.dot(r, t) + dh_dp * _pressure->dot(r, t);
});
}
}
Original file line number Diff line number Diff line change
Expand Up @@ -10,5 +10,7 @@
design = 'FunctorErgunDragCoefficients.md'
valgrind = 'none'
method = '!dbg'
# GlobalParams is setting fp to an object not using it at the time
allow_warnings = true
[]
[]
Loading