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             : //  "VG1" form of relative permeability
      11             : //
      12             : #include "RichardsRelPermVG1.h"
      13             : 
      14             : registerMooseObject("RichardsApp", RichardsRelPermVG1);
      15             : 
      16             : InputParameters
      17          14 : RichardsRelPermVG1::validParams()
      18             : {
      19          14 :   InputParameters params = RichardsRelPermVG::validParams();
      20          28 :   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 = s^(1/2) * (1 - (1 "
      25             :       "- s^(1/m))^m)^2");
      26          28 :   params.addRequiredRangeCheckedParam<Real>(
      27             :       "m",
      28             :       "m > 0 & m < 1",
      29             :       "van-Genuchten m parameter.  Must be between 0 and 1, and optimally "
      30             :       "should be set >0.5.  Define s = (seff - simm)/(1 - simm).  Then "
      31             :       "relperm = s^(1/2) * (1 - (1 - s^(1/m))^m)^2");
      32          28 :   params.addRequiredRangeCheckedParam<Real>(
      33             :       "scut", "scut > 0 & scut < 1", "cutoff in effective saturation.");
      34          14 :   params.addClassDescription("VG1 form of relative permeability.  Define s = (seff - simm)/(1 - "
      35             :                              "simm).  Then relperm = s^(1/2) * (1 - (1 - s^(1/m))^m)^2, if s>0, "
      36             :                              "and relperm=0 otherwise");
      37          14 :   return params;
      38           0 : }
      39             : 
      40           6 : RichardsRelPermVG1::RichardsRelPermVG1(const InputParameters & parameters)
      41             :   : RichardsRelPermVG(parameters),
      42           6 :     _simm(getParam<Real>("simm")),
      43          12 :     _m(getParam<Real>("m")),
      44          12 :     _scut(getParam<Real>("scut")),
      45           6 :     _vg1_const(0),
      46           6 :     _vg1_linear(0),
      47           6 :     _vg1_quad(0),
      48           6 :     _vg1_cub(0)
      49             : {
      50           6 :   _vg1_const = RichardsRelPermVG::relperm(_scut);
      51           6 :   _vg1_linear = RichardsRelPermVG::drelperm(_scut);
      52           6 :   _vg1_quad = RichardsRelPermVG::d2relperm(_scut);
      53           6 :   _vg1_cub = (1 - _vg1_const - _vg1_linear * (1 - _scut) - _vg1_quad * std::pow(1 - _scut, 2)) /
      54             :              std::pow(1 - _scut, 3);
      55           6 : }
      56             : 
      57             : void
      58           4 : RichardsRelPermVG1::initialSetup()
      59             : {
      60           4 :   _console << "Relative permeability of VG1 type has cubic coefficients " << _vg1_const << " "
      61           4 :            << _vg1_linear << " " << _vg1_quad << " " << _vg1_cub << std::endl;
      62           4 : }
      63             : 
      64             : Real
      65       98046 : RichardsRelPermVG1::relperm(Real seff) const
      66             : {
      67       98046 :   if (seff >= 1.0)
      68             :     return 1.0;
      69             : 
      70       98046 :   if (seff <= _simm)
      71             :     return 0.0;
      72             : 
      73       98046 :   Real s_internal = (seff - _simm) / (1.0 - _simm);
      74             : 
      75       98046 :   if (s_internal < _scut)
      76       96186 :     return RichardsRelPermVG::relperm(seff);
      77             : 
      78        1860 :   Real krel = _vg1_const + _vg1_linear * (s_internal - _scut) +
      79        1860 :               _vg1_quad * std::pow(s_internal - _scut, 2) +
      80        1860 :               _vg1_cub * std::pow(s_internal - _scut, 3);
      81             : 
      82             :   // bound, just in case
      83        1860 :   if (krel < 0)
      84             :   {
      85             :     krel = 0;
      86             :   }
      87        1860 :   if (krel > 1)
      88             :   {
      89             :     krel = 1;
      90             :   }
      91             :   return krel;
      92             : }
      93             : 
      94             : Real
      95      195486 : RichardsRelPermVG1::drelperm(Real seff) const
      96             : {
      97      195486 :   if (seff >= 1.0)
      98             :     return 0.0;
      99             : 
     100      195486 :   if (seff <= _simm)
     101             :     return 0.0;
     102             : 
     103      195486 :   Real s_internal = (seff - _simm) / (1.0 - _simm);
     104             : 
     105      195486 :   if (s_internal < _scut)
     106      192372 :     return RichardsRelPermVG::drelperm(seff);
     107             : 
     108        3114 :   Real krelp = _vg1_linear + 2 * _vg1_quad * (s_internal - _scut) +
     109        3114 :                3 * _vg1_cub * std::pow(s_internal - _scut, 2);
     110        3114 :   return krelp / (1.0 - _simm);
     111             : }
     112             : 
     113             : Real
     114       98046 : RichardsRelPermVG1::d2relperm(Real seff) const
     115             : {
     116       98046 :   if (seff >= 1.0)
     117             :     return 0.0;
     118             : 
     119       98046 :   if (seff <= _simm)
     120             :     return 0.0;
     121             : 
     122       98046 :   Real s_internal = (seff - _simm) / (1.0 - _simm);
     123             : 
     124       98046 :   if (s_internal < _scut)
     125       96186 :     return RichardsRelPermVG::d2relperm(seff);
     126             : 
     127        1860 :   Real krelpp = 2 * _vg1_quad + 6 * _vg1_cub * (s_internal - _scut);
     128        1860 :   return krelpp / std::pow(1.0 - _simm, 2);
     129             : }