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             : //  "Broadbridge-White" form of relative permeability (P Broadbridge and I White ``Constant rate
      11             : //  rainfall infiltration: A versatile nonlinear model 1. Analytic Solution'', Water Resources
      12             : //  Research 24 (1988) 145-154)
      13             : //
      14             : #include "RichardsRelPermBW.h"
      15             : #include "libmesh/utility.h"
      16             : 
      17             : registerMooseObject("RichardsApp", RichardsRelPermBW);
      18             : 
      19             : InputParameters
      20          20 : RichardsRelPermBW::validParams()
      21             : {
      22          20 :   InputParameters params = RichardsRelPerm::validParams();
      23          40 :   params.addRequiredRangeCheckedParam<Real>(
      24             :       "Sn",
      25             :       "Sn >= 0",
      26             :       "Low saturation.  This must be < Ss, and non-negative.  This is BW's "
      27             :       "initial effective saturation, below which effective saturation never goes "
      28             :       "in their simulations/models.  If Kn=0 then Sn is the immobile saturation.");
      29          60 :   params.addRangeCheckedParam<Real>(
      30             :       "Ss",
      31          40 :       1.0,
      32             :       "Ss <= 1",
      33             :       "High saturation.  This must be > Sn and <= 1.  Effective saturation "
      34             :       "where porepressure = 0.  Effective saturation never exceeds this "
      35             :       "value in BW's simulations/models.");
      36          60 :   params.addRangeCheckedParam<Real>(
      37          40 :       "Kn", 0.0, "Kn >= 0", "Relative permeability at Seff = Sn.  Must be < Ks");
      38          60 :   params.addRangeCheckedParam<Real>(
      39          40 :       "Ks", 1.0, "Ks <= 1", "Relative permeability at Seff = Ss.  Must be > Kn");
      40          40 :   params.addRequiredRangeCheckedParam<Real>(
      41             :       "C",
      42             :       "C > 1",
      43             :       "BW's C parameter.  Must be > 1.   Define s = (seff - Sn)/(Ss - Sn).  Then "
      44             :       "relperm = Kn + s^2(c-1)(Kn-Ks)/(c-s) if 0<s<1, otherwise relperm = Kn if "
      45             :       "s<=0, otherwise relperm = Ks if s>=1.");
      46          20 :   params.addClassDescription("Broadbridge-White form of relative permeability.  Define s = (seff - "
      47             :                              "Sn)/(Ss - Sn).  Then relperm = Kn + s^2(c-1)(Kn-Ks)/(c-s) if 0<s<1, "
      48             :                              "otherwise relperm = Kn if s<=0, otherwise relperm = Ks if s>=1.");
      49          20 :   return params;
      50           0 : }
      51             : 
      52          10 : RichardsRelPermBW::RichardsRelPermBW(const InputParameters & parameters)
      53             :   : RichardsRelPerm(parameters),
      54          10 :     _sn(getParam<Real>("Sn")),
      55          20 :     _ss(getParam<Real>("Ss")),
      56          20 :     _kn(getParam<Real>("Kn")),
      57          20 :     _ks(getParam<Real>("Ks")),
      58          30 :     _c(getParam<Real>("C"))
      59             : {
      60          10 :   if (_ss <= _sn)
      61           1 :     mooseError("In BW relative permeability Sn set to ",
      62           1 :                _sn,
      63             :                " and Ss set to ",
      64           1 :                _ss,
      65             :                " but these must obey Ss > Sn");
      66           9 :   if (_ks <= _kn)
      67           1 :     mooseError("In BW relative permeability Kn set to ",
      68           1 :                _kn,
      69             :                " and Ks set to ",
      70           1 :                _ks,
      71             :                " but these must obey Ks > Kn");
      72           8 :   _coef = (_ks - _kn) * (_c - 1); // shorthand coefficient
      73           8 : }
      74             : 
      75             : Real
      76      696382 : RichardsRelPermBW::relperm(Real seff) const
      77             : {
      78      696382 :   if (seff <= _sn)
      79          30 :     return _kn;
      80             : 
      81      696352 :   if (seff >= _ss)
      82        5442 :     return _ks;
      83             : 
      84      690910 :   const Real s_internal = (seff - _sn) / (_ss - _sn);
      85      690910 :   const Real krel = _kn + _coef * Utility::pow<2>(s_internal) / (_c - s_internal);
      86      690910 :   return krel;
      87             : }
      88             : 
      89             : Real
      90     1392158 : RichardsRelPermBW::drelperm(Real seff) const
      91             : {
      92     1392158 :   if (seff <= _sn)
      93             :     return 0.0;
      94             : 
      95     1392128 :   if (seff >= _ss)
      96             :     return 0.0;
      97             : 
      98     1381274 :   const Real s_internal = (seff - _sn) / (_ss - _sn);
      99     1381274 :   const Real krelp = _coef * (2.0 * s_internal / (_c - s_internal) +
     100     1381274 :                               Utility::pow<2>(s_internal) / Utility::pow<2>(_c - s_internal));
     101     1381274 :   return krelp / (_ss - _sn);
     102             : }
     103             : 
     104             : Real
     105      696382 : RichardsRelPermBW::d2relperm(Real seff) const
     106             : {
     107      696382 :   if (seff <= _sn)
     108             :     return 0.0;
     109             : 
     110      696352 :   if (seff >= _ss)
     111             :     return 0.0;
     112             : 
     113      690910 :   const Real s_internal = (seff - _sn) / (_ss - _sn);
     114             :   const Real krelpp =
     115      690910 :       _coef * (2.0 / (_c - s_internal) + 4.0 * s_internal / Utility::pow<2>(_c - s_internal) +
     116      690910 :                2.0 * Utility::pow<2>(s_internal) / Utility::pow<3>(_c - s_internal));
     117      690910 :   return krelpp / Utility::pow<2>(_ss - _sn);
     118             : }