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