Numerics.h

Go to the documentation of this file.

//* This file is part of the MOOSE framework
//* https://mooseframework.inl.gov
//*
//* 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 "libmesh/libmesh_common.h"
#include "libmesh/vector_value.h"
#include "libmesh/dense_vector.h"

#include "ADReal.h"
#include "NavierStokesMethods.h"
#include "HeatTransferUtils.h"

using namespace libMesh;

namespace THM
{

// Default value for magnitude of acceleration due to gravity
static const Real gravity_const = 9.81;

// Default value for gravitational acceleration vector
static VectorValue<Real> default_gravity_vector = VectorValue<Real>(0.0, 0.0, -gravity_const);

// Stefan-Boltzman constant, in [W/m^2-K]
static const Real Stefan_Boltzman_const = 5.670e-8;

template <typename T>
int
sgn(T val)
{
  return (T(0) < val) - (val < T(0));
}

bool absoluteFuzzyEqualVectors(const RealVectorValue & a,
                               const RealVectorValue & b,
                               const Real & tol = libMesh::TOLERANCE * libMesh::TOLERANCE);

bool areParallelVectors(const RealVectorValue & a,
                        const RealVectorValue & b,
                        const Real & tol = libMesh::TOLERANCE * libMesh::TOLERANCE);

bool haveSameDirection(const RealVectorValue & a,
                       const RealVectorValue & b,
                       const Real & tol = libMesh::TOLERANCE * libMesh::TOLERANCE);

Real
applyQuotientRule(const Real & num, const Real & den, const Real & dnum_dy, const Real & dden_dy);

DenseVector<Real> applyQuotientRule(const Real & num,
                                    const Real & den,
                                    const DenseVector<Real> & dnum_dy,
                                    const DenseVector<Real> & dden_dy);

template <typename T1, typename T2, typename T3, typename T4, typename T5>
auto
Reynolds(const T1 & volume_fraction, const T2 & rho, const T3 & vel, const T4 & D_h, const T5 & mu)
{
  return volume_fraction * HeatTransferUtils::reynolds(rho, vel, D_h, mu);
}

template <typename T1, typename T2, typename T3>
auto
Prandtl(const T1 & cp, const T2 & mu, const T3 & k)
{
  return HeatTransferUtils::prandtl(cp, mu, k);
}

template <typename T1, typename T2, typename T3, typename T4, typename T5, typename T6>
auto
Peclet(const T1 & volume_fraction,
       const T2 & cp,
       const T3 & rho,
       const T4 & vel,
       const T5 & D_h,
       const T6 & k)
{
  const auto diffusivity = HeatTransferUtils::thermalDiffusivity(k, rho, cp);
  return volume_fraction * HeatTransferUtils::peclet(vel, D_h, diffusivity);
}

template <typename T1, typename T2, typename T3, typename T4, typename T5>
auto
Grashof(const T1 & beta,
        const T2 & dT,
        const T3 & D_h,
        const T4 & rho_liquid,
        const T5 & mu_liquid,
        const Real & gravity_magnitude)
{
  using std::pow;
  return gravity_magnitude * beta * dT * pow(D_h, 3) * (rho_liquid * rho_liquid) /
         (mu_liquid * mu_liquid);
}

template <typename T1, typename T2>
auto
Laplace(const T1 & surf_tension, const T2 & delta_rho, const Real & gravity_magnitude)
{
  using std::sqrt;
  return sqrt(surf_tension / (gravity_magnitude * delta_rho));
}

template <typename T1, typename T2, typename T3, typename T4>
auto
viscosityNumber(const T1 & viscosity,
                const T2 & surf_tension,
                const T3 & rho_k,
                const T4 & delta_rho,
                const Real & gravity_magnitude)
{
  using std::sqrt;
  return viscosity /
         sqrt(rho_k * surf_tension * sqrt(surf_tension / gravity_magnitude / delta_rho));
}

using NS::wallHeatTransferCoefficient;

template <typename T1, typename T2>
auto
Dean(const T1 & Re, const T2 & doD)
{
  using std::sqrt;
  return Re * sqrt(doD);
}

void
vel_from_arhoA_arhouA(Real arhoA, Real arhouA, Real & vel, Real & dvel_darhoA, Real & dvel_darhouA);

ADReal vel_from_arhoA_arhouA(ADReal arhoA, ADReal arhouA);

Real dvel_darhoA(Real arhoA, Real arhouA);

Real dvel_darhouA(Real arhoA);

void rho_from_arhoA_alpha_A(
    Real arhoA, Real alpha, Real A, Real & rho, Real & drho_darhoA, Real & drho_dalpha);

ADReal rho_from_arhoA_alpha_A(ADReal arhoA, ADReal alpha, ADReal A);

void v_from_rhoA_A(Real rhoA, Real A, Real & v, Real & dv_drhoA);

ADReal v_from_rhoA_A(ADReal rhoA, ADReal A);

void
v_from_arhoA_alpha_A(Real arhoA, Real alpha, Real A, Real & v, Real & dv_darhoA, Real & dv_dalpha);

ADReal v_from_arhoA_alpha_A(ADReal arhoA, ADReal alpha, ADReal A);

void v_from_rho(Real rho, Real & v, Real & dv_drho);

Real dv_dalpha_liquid(Real area, Real arhoA, bool is_liquid);

Real dv_darhoA(Real area, Real arhoA);

void e_from_arhoA_arhouA_arhoEA(Real arhoA,
                                Real arhouA,
                                Real arhoEA,
                                Real & e,
                                Real & de_darhoA,
                                Real & de_darhouA,
                                Real & de_darhoEA);

ADReal e_from_arhoA_arhouA_arhoEA(ADReal arhoA, ADReal arhouA, ADReal arhoEA);

void e_from_E_vel(Real E, Real vel, Real & e, Real & de_dE, Real & de_dvel);

ADReal e_from_E_vel(ADReal E, ADReal vel);

Real de_darhoA(Real arhoA, Real arhouA, Real arhoEA);

Real de_darhouA(Real arhoA, Real arhouA);

Real de_darhoEA(Real arhoA);

void E_from_arhoA_arhoEA(Real arhoA, Real arhoEA, Real & E, Real & dE_darhoA, Real & dE_darhoEA);

ADReal E_from_arhoA_arhoEA(ADReal arhoA, ADReal arhoEA);

void E_from_e_vel(Real e, Real vel, Real & E, Real & dE_de, Real & dE_dvel);

void h_from_e_p_rho(Real e, Real p, Real rho, Real & h, Real & dh_de, Real & dh_dp, Real & dh_drho);

ADReal h_from_e_p_rho(ADReal e, ADReal p, ADReal rho);

bool isInlet(Real vel, Real normal);
bool isInlet(ADReal vel, Real normal);

bool isOutlet(Real vel, Real normal);
bool isOutlet(ADReal vel, Real normal);
}