PorousFlowVanGenuchten.h

Go to the documentation of this file.

//* 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 "MooseTypes.h"
#include "libmesh/utility.h"

namespace PorousFlowVanGenuchten
{
Real effectiveSaturation(Real p, Real alpha, Real m);

Real dEffectiveSaturation(Real p, Real alpha, Real m);

Real d2EffectiveSaturation(Real p, Real alpha, Real m);

Real capillaryPressure(Real seff, Real alpha, Real m, Real pc_max);

Real dCapillaryPressure(Real seff, Real alpha, Real m, Real pc_max);

Real d2CapillaryPressure(Real seff, Real alpha, Real m, Real pc_max);

template <typename T>
T
relativePermeability(const T & seff, Real m)
{
  if (MetaPhysicL::raw_value(seff) <= 0.0)
    return 0.0;
  else if (MetaPhysicL::raw_value(seff) >= 1.0)
    return 1.0;

  const T a = 1.0 - std::pow(seff, 1.0 / m);
  const T b = 1.0 - std::pow(a, m);

  return std::sqrt(seff) * Utility::pow<2>(b);
}

Real dRelativePermeability(Real seff, Real m);

Real d2RelativePermeability(Real seff, Real m);

template <typename T>
T
relativePermeabilityNW(const T & seff, Real m)
{
  if (MetaPhysicL::raw_value(seff) <= 0.0)
    return 0.0;
  else if (MetaPhysicL::raw_value(seff) >= 1.0)
    return 1.0;

  const T a = std::pow(1.0 - seff, 1.0 / m);
  const T b = std::pow(1.0 - a, 2.0 * m);

  return std::sqrt(seff) * b;
}

Real dRelativePermeabilityNW(Real seff, Real m);

Real d2RelativePermeabilityNW(Real seff, Real m);

struct LowCapillaryPressureExtension
{
  enum ExtensionStrategy
  {
    NONE,
    QUADRATIC,
    EXPONENTIAL
  };
  ExtensionStrategy strategy;
  Real S;
  Real Pc;
  Real dPc;

  LowCapillaryPressureExtension()
    : strategy(LowCapillaryPressureExtension::NONE),
      S(0.0),
      Pc(std::numeric_limits<Real>::max()),
      dPc(std::numeric_limits<Real>::lowest()){};

  LowCapillaryPressureExtension(const ExtensionStrategy & strategy, Real S, Real Pc, Real dPc)
    : strategy(strategy), S(S), Pc(Pc), dPc(dPc){};
};

struct HighCapillaryPressureExtension
{
  enum ExtensionStrategy
  {
    NONE,
    POWER
  };
  ExtensionStrategy strategy;
  Real S;
  Real Pc;
  Real dPc;

  HighCapillaryPressureExtension()
    : strategy(HighCapillaryPressureExtension::NONE),
      S(1.0),
      Pc(0.0),
      dPc(std::numeric_limits<Real>::lowest()){};

  HighCapillaryPressureExtension(const ExtensionStrategy & strategy, Real S, Real Pc, Real dPc)
    : strategy(strategy), S(S), Pc(Pc), dPc(dPc){};
};

Real capillaryPressureHys(
    Real sl,
    Real slmin,
    Real sgrdel,
    Real alpha,
    Real n,
    const LowCapillaryPressureExtension & low_ext = LowCapillaryPressureExtension(),
    const HighCapillaryPressureExtension & high_ext = HighCapillaryPressureExtension());

Real dcapillaryPressureHys(
    Real sl,
    Real slmin,
    Real sgrdel,
    Real alpha,
    Real n,
    const LowCapillaryPressureExtension & low_ext = LowCapillaryPressureExtension(),
    const HighCapillaryPressureExtension & high_ext = HighCapillaryPressureExtension());

Real d2capillaryPressureHys(
    Real sl,
    Real slmin,
    Real sgrdel,
    Real alpha,
    Real n,
    const LowCapillaryPressureExtension & low_ext = LowCapillaryPressureExtension(),
    const HighCapillaryPressureExtension & high_ext = HighCapillaryPressureExtension());

Real
saturationHys(Real pc,
              Real slmin,
              Real sgrdel,
              Real alpha,
              Real n,
              const LowCapillaryPressureExtension & low_ext = LowCapillaryPressureExtension(),
              const HighCapillaryPressureExtension & high_ext = HighCapillaryPressureExtension());

Real
dsaturationHys(Real pc,
               Real slmin,
               Real sgrdel,
               Real alpha,
               Real n,
               const LowCapillaryPressureExtension & low_ext = LowCapillaryPressureExtension(),
               const HighCapillaryPressureExtension & high_ext = HighCapillaryPressureExtension());

Real
d2saturationHys(Real pc,
                Real slmin,
                Real sgrdel,
                Real alpha,
                Real n,
                const LowCapillaryPressureExtension & low_ext = LowCapillaryPressureExtension(),
                const HighCapillaryPressureExtension & high_ext = HighCapillaryPressureExtension());

Real relativePermeabilityHys(Real sl,
                             Real slr,
                             Real sgrdel,
                             Real sgrmax,
                             Real sldel,
                             Real m,
                             Real upper_liquid_param,
                             Real y0,
                             Real y0p,
                             Real y1,
                             Real y1p);

Real drelativePermeabilityHys(Real sl,
                              Real slr,
                              Real sgrdel,
                              Real sgrmax,
                              Real sldel,
                              Real m,
                              Real upper_liquid_param,
                              Real y0,
                              Real y0p,
                              Real y1,
                              Real y1p);

Real relativePermeabilityNWHys(Real sl,
                               Real slr,
                               Real sgrdel,
                               Real sgrmax,
                               Real sldel,
                               Real m,
                               Real gamma,
                               Real k_rg_max,
                               Real y0p);

Real drelativePermeabilityNWHys(Real sl,
                                Real slr,
                                Real sgrdel,
                                Real sgrmax,
                                Real sldel,
                                Real m,
                                Real gamma,
                                Real k_rg_max,
                                Real y0p);
}