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             : //  "Power" form of relative permeability
      11             : //
      12             : #include "RichardsRelPermPower.h"
      13             : 
      14             : registerMooseObject("RichardsApp", RichardsRelPermPower);
      15             : 
      16             : InputParameters
      17         843 : RichardsRelPermPower::validParams()
      18             : {
      19         843 :   InputParameters params = RichardsRelPerm::validParams();
      20        1686 :   params.addRequiredRangeCheckedParam<Real>(
      21             :       "simm",
      22             :       "simm >= 0 & simm < 1",
      23             :       "Immobile saturation.  Must be between 0 and 1.   Define s = "
      24             :       "(seff - simm)/(1 - simm).  Then relperm = (n+1)s^n - "
      25             :       "ns^(n+1)");
      26        1686 :   params.addRequiredRangeCheckedParam<Real>("n",
      27             :                                             "n >= 2",
      28             :                                             "Exponent.  Must be >= 2.   Define s = "
      29             :                                             "(seff - simm)/(1 - simm).  Then "
      30             :                                             "relperm = (n+1)s^n - ns^(n+1)");
      31         843 :   params.addClassDescription("Power form of relative permeability.  Define s = (seff - simm)/(1 - "
      32             :                              "simm).  Then relperm = (n+1)s^n - ns^(n+1) if s<simm, otherwise "
      33             :                              "relperm=1");
      34         843 :   return params;
      35           0 : }
      36             : 
      37         416 : RichardsRelPermPower::RichardsRelPermPower(const InputParameters & parameters)
      38        1248 :   : RichardsRelPerm(parameters), _simm(getParam<Real>("simm")), _n(getParam<Real>("n"))
      39             : {
      40         416 : }
      41             : 
      42             : Real
      43     6756526 : RichardsRelPermPower::relperm(Real seff) const
      44             : {
      45     6756526 :   if (seff >= 1.0)
      46             :     return 1.0;
      47             : 
      48     2855696 :   if (seff <= _simm)
      49             :     return 0.0;
      50             : 
      51     2718004 :   Real s_internal = (seff - _simm) / (1.0 - _simm);
      52     2718004 :   Real krel = (_n + 1) * std::pow(s_internal, _n) - _n * std::pow(s_internal, _n + 1);
      53             : 
      54             :   // bound, just in case
      55     2718004 :   if (krel < 0)
      56             :   {
      57             :     krel = 0;
      58             :   }
      59     2718004 :   if (krel > 1)
      60             :   {
      61             :     krel = 1;
      62             :   }
      63             :   return krel;
      64             : }
      65             : 
      66             : Real
      67    11160762 : RichardsRelPermPower::drelperm(Real seff) const
      68             : {
      69    11160762 :   if (seff >= 1.0)
      70             :     return 0.0;
      71             : 
      72     4440077 :   if (seff <= _simm)
      73             :     return 0.0;
      74             : 
      75     4213322 :   Real s_internal = (seff - _simm) / (1.0 - _simm);
      76             :   Real krelp =
      77     4213322 :       (_n + 1) * _n * std::pow(s_internal, _n - 1) - _n * (_n + 1) * std::pow(s_internal, _n);
      78     4213322 :   return krelp / (1.0 - _simm);
      79             : }
      80             : 
      81             : Real
      82     4552486 : RichardsRelPermPower::d2relperm(Real seff) const
      83             : {
      84     4552486 :   if (seff >= 1.0)
      85             :     return 0.0;
      86             : 
      87     1732487 :   if (seff <= _simm)
      88             :     return 0.0;
      89             : 
      90     1643332 :   Real s_internal = (seff - _simm) / (1.0 - _simm);
      91     1643332 :   Real krelpp = (_n + 1) * _n * (_n - 1) * std::pow(s_internal, _n - 2) -
      92     1643332 :                 _n * (_n + 1) * _n * std::pow(s_internal, _n - 1);
      93     1643332 :   return krelpp / std::pow(1.0 - _simm, 2);
      94             : }