Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://mooseframework.inl.gov
       3             : //*
       4             : //* All rights reserved, see COPYRIGHT for full restrictions
       5             : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
       6             : //*
       7             : //* Licensed under LGPL 2.1, please see LICENSE for details
       8             : //* https://www.gnu.org/licenses/lgpl-2.1.html
       9             : 
      10             : #include "RichardsRelPermVG.h"
      11             : #include "libmesh/utility.h"
      12             : 
      13             : registerMooseObject("RichardsApp", RichardsRelPermVG);
      14             : 
      15             : InputParameters
      16          24 : RichardsRelPermVG::validParams()
      17             : {
      18          24 :   InputParameters params = RichardsRelPerm::validParams();
      19          48 :   params.addRequiredRangeCheckedParam<Real>(
      20             :       "simm",
      21             :       "simm >= 0 & simm < 1",
      22             :       "Immobile saturation.  Must be between 0 and 1.  Define s = "
      23             :       "(seff - simm)/(1 - simm).  Then relperm = s^(1/2) * (1 - (1 "
      24             :       "- s^(1/m))^m)^2");
      25          48 :   params.addRequiredRangeCheckedParam<Real>(
      26             :       "m",
      27             :       "m > 0 & m < 1",
      28             :       "van-Genuchten m parameter.  Must be between 0 and 1, and optimally "
      29             :       "should be set >0.5.  Define s = (seff - simm)/(1 - simm).  Then "
      30             :       "relperm = s^(1/2) * (1 - (1 - s^(1/m))^m)^2");
      31          24 :   params.addClassDescription("VG form of relative permeability.  Define s = (seff - simm)/(1 - "
      32             :                              "simm).  Then relperm = s^(1/2) * (1 - (1 - s^(1/m))^m)^2, if s>0, "
      33             :                              "and relperm=0 otherwise");
      34          24 :   return params;
      35           0 : }
      36             : 
      37          10 : RichardsRelPermVG::RichardsRelPermVG(const InputParameters & parameters)
      38          30 :   : RichardsRelPerm(parameters), _simm(getParam<Real>("simm")), _m(getParam<Real>("m"))
      39             : {
      40          10 : }
      41             : 
      42             : Real
      43       96798 : RichardsRelPermVG::relperm(Real seff) const
      44             : {
      45       96798 :   if (seff >= 1.0)
      46             :     return 1.0;
      47             : 
      48       96798 :   if (seff <= _simm)
      49             :     return 0.0;
      50             : 
      51       96798 :   Real s_internal = (seff - _simm) / (1.0 - _simm);
      52       96798 :   Real krel = std::sqrt(s_internal) *
      53       96798 :               Utility::pow<2>(1.0 - std::pow(1.0 - std::pow(s_internal, 1.0 / _m), _m));
      54             : 
      55             :   // bound, just in case
      56       96798 :   if (krel < 0.0)
      57             :     krel = 0.0;
      58       96798 :   if (krel > 1.0)
      59             :     krel = 1.0;
      60             : 
      61             :   return krel;
      62             : }
      63             : 
      64             : Real
      65      192984 : RichardsRelPermVG::drelperm(Real seff) const
      66             : {
      67      192984 :   if (seff >= 1.0)
      68             :     return 0.0;
      69             : 
      70      192984 :   if (seff <= _simm)
      71             :     return 0.0;
      72             : 
      73      192984 :   Real s_internal = (seff - _simm) / (1.0 - _simm);
      74      192984 :   Real tmp = 1.0 - std::pow(s_internal, 1.0 / _m);
      75      192984 :   Real tmpp = -1.0 / _m * std::pow(s_internal, 1.0 / _m - 1.0);
      76      192984 :   Real tmp2 = 1.0 - std::pow(tmp, _m);
      77      192984 :   Real tmp2p = -_m * std::pow(tmp, _m - 1.0) * tmpp;
      78             :   // Real krel = std::sqrt(s_internal)*std::pow(tmp2, 2);
      79             :   Real krelp =
      80      192984 :       0.5 * std::pow(s_internal, -0.5) * tmp2 * tmp2 + 2.0 * std::sqrt(s_internal) * tmp2 * tmp2p;
      81      192984 :   return krelp / (1.0 - _simm);
      82             : }
      83             : 
      84             : Real
      85       96798 : RichardsRelPermVG::d2relperm(Real seff) const
      86             : {
      87       96798 :   if (seff >= 1.0)
      88             :     return 0.0;
      89             : 
      90       96798 :   if (seff <= _simm)
      91             :     return 0.0;
      92             : 
      93       96798 :   Real s_internal = (seff - _simm) / (1.0 - _simm);
      94       96798 :   Real tmp = 1.0 - std::pow(s_internal, 1.0 / _m);
      95       96798 :   Real tmpp = -1.0 / _m * std::pow(s_internal, 1.0 / _m - 1.0);
      96       96798 :   Real tmppp = -1.0 / _m * (1.0 / _m - 1.0) * std::pow(s_internal, 1.0 / _m - 2);
      97       96798 :   Real tmp2 = 1.0 - std::pow(tmp, _m);
      98       96798 :   Real tmp2p = -_m * std::pow(tmp, _m - 1.0) * tmpp;
      99       96798 :   Real tmp2pp = -_m * (_m - 1.0) * std::pow(tmp, _m - 2.0) * tmpp * tmpp -
     100       96798 :                 _m * std::pow(tmp, _m - 1.0) * tmppp;
     101             :   // Real krel = std::sqrt(s_internal)*std::pow(tmp2, 2);
     102             :   // Real krelp = 0.5 * std::pow(s_internal, -0.5)*std::pow(tmp2, 2) +
     103             :   // 2*std::sqrt(s_internal)*tmp2*tmp2p;
     104       96798 :   Real krelpp = -0.25 * std::pow(s_internal, -1.5) * tmp2 * tmp2 +
     105       96798 :                 2.0 * 0.5 * std::pow(s_internal, -0.5) * 2.0 * tmp2 * tmp2p +
     106       96798 :                 2.0 * std::sqrt(s_internal) * (tmp2p * tmp2p + tmp2 * tmp2pp);
     107             : 
     108       96798 :   return krelpp / Utility::pow<2>(1.0 - _simm);
     109             : }