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 "PorousFlowHystereticRelativePermeabilityBase.h"
      11             : 
      12             : InputParameters
      13         896 : PorousFlowHystereticRelativePermeabilityBase::validParams()
      14             : {
      15         896 :   InputParameters params = PorousFlowMaterialBase::validParams();
      16        2688 :   params.addRangeCheckedParam<Real>(
      17             :       "S_lr",
      18        1792 :       0.0,
      19             :       "S_lr >= 0 & S_lr < 1",
      20             :       "Liquid residual saturation.  At liquid saturation = S_lr, the liquid relative permeability "
      21             :       "is zero and the gas relative permeability is k_rg_max.  Almost definitely you need to set "
      22             :       "S_lr equal to the quantity used for your hysteretic capillary-pressure curve, if you are "
      23             :       "using one.");
      24        1792 :   params.addRequiredRangeCheckedParam<Real>(
      25             :       "S_gr_max",
      26             :       "S_gr_max >= 0 & S_gr_max < 1",
      27             :       "Residual gas saturation.  For liquid saturation = 1 - S_gr_max, the gas wetting "
      28             :       "relative-permeability is zero, if the turning point was <= S_lr.  1 - S_gr_max is the "
      29             :       "maximum liquid saturation for which the van Genuchten expression is valid for the liquid "
      30             :       "relative-permeability wetting curve.  You must ensure S_gr_max < 1 - S_lr.  Often S_gr_max "
      31             :       "= -0.3136 * ln(porosity) - 0.1334 is used");
      32        1792 :   params.addRequiredRangeCheckedParam<Real>(
      33             :       "m", "m > 0 & m < 1", "van Genuchten m parameter.  Suggested value is around 0.9");
      34        1792 :   params.addPrivateParam<std::string>("pf_material_type", "relative_permeability");
      35         896 :   params.addPrivateParam<bool>("is_ad", false);
      36         896 :   params.addClassDescription("PorousFlow material that computes relative permeability for 1-phase "
      37             :                              "or 2-phase models that include hysteresis");
      38         896 :   return params;
      39           0 : }
      40             : 
      41         698 : PorousFlowHystereticRelativePermeabilityBase::PorousFlowHystereticRelativePermeabilityBase(
      42         698 :     const InputParameters & parameters)
      43             :   : PorousFlowMaterialBase(parameters),
      44         698 :     _s_lr(getParam<Real>("S_lr")),
      45        1396 :     _s_gr_max(getParam<Real>("S_gr_max")),
      46        1396 :     _m(getParam<Real>("m")),
      47        1396 :     _saturation(_nodal_material
      48         698 :                     ? getMaterialProperty<std::vector<Real>>("PorousFlow_saturation_nodal")
      49        2022 :                     : getMaterialProperty<std::vector<Real>>("PorousFlow_saturation_qp")),
      50         698 :     _hys_order(_nodal_material ? getMaterialProperty<unsigned>("PorousFlow_hysteresis_order_nodal")
      51        2022 :                                : getMaterialProperty<unsigned>("PorousFlow_hysteresis_order_qp")),
      52        1396 :     _hys_order_old(_nodal_material
      53         698 :                        ? getMaterialPropertyOld<unsigned>("PorousFlow_hysteresis_order_nodal")
      54        2022 :                        : getMaterialPropertyOld<unsigned>("PorousFlow_hysteresis_order_qp")),
      55         698 :     _hys_sat_tps(
      56         698 :         _nodal_material
      57         698 :             ? getMaterialProperty<std::array<Real, PorousFlowConstants::MAX_HYSTERESIS_ORDER>>(
      58             :                   "PorousFlow_hysteresis_saturation_tps_nodal")
      59        2022 :             : getMaterialProperty<std::array<Real, PorousFlowConstants::MAX_HYSTERESIS_ORDER>>(
      60             :                   "PorousFlow_hysteresis_saturation_tps_qp")),
      61         698 :     _relative_permeability(
      62         698 :         _nodal_material ? declareProperty<Real>("PorousFlow_relative_permeability_nodal" + _phase)
      63        1360 :                         : declareProperty<Real>("PorousFlow_relative_permeability_qp" + _phase)),
      64        1396 :     _drelative_permeability_ds(
      65         698 :         _nodal_material
      66         770 :             ? declarePropertyDerivative<Real>("PorousFlow_relative_permeability_nodal" + _phase,
      67             :                                               _saturation_variable_name)
      68        2684 :             : declarePropertyDerivative<Real>("PorousFlow_relative_permeability_qp" + _phase,
      69             :                                               _saturation_variable_name)),
      70        1396 :     _s_gr_tp0(_nodal_material
      71         698 :                   ? declareProperty<Real>("PorousFlow_hysteresis_s_gr_tp0_nodal" + _phase)
      72        2058 :                   : declareProperty<Real>("PorousFlow_hysteresis_s_gr_tp0_qp" + _phase))
      73             : {
      74             :   mooseAssert(_dictator.numPhases() <= 2,
      75             :               "PorousFlowHystereticRelativePermeability Materials can only be used for "
      76             :               "1-phase or 2-phase models.  The Dictator proclaims that the number of phases "
      77             :               "exceeds 2, and dictators are never incorrect.");
      78         698 :   if (_s_gr_max >= 1 - _s_lr)
      79           2 :     paramError("S_gr_max", "Must be less than 1 - S_lr");
      80         696 : }
      81             : 
      82             : void
      83         452 : PorousFlowHystereticRelativePermeabilityBase::initQpStatefulProperties()
      84             : {
      85         452 :   PorousFlowMaterialBase::initQpStatefulProperties();
      86         452 :   if (_hys_order[_qp] >= 1)
      87          16 :     computeTurningPoint0Info(_hys_sat_tps[_qp].at(0));
      88             : 
      89         452 :   computeRelPermQp();
      90         452 : }
      91             : 
      92             : void
      93       62028 : PorousFlowHystereticRelativePermeabilityBase::computeQpProperties()
      94             : {
      95       62028 :   PorousFlowMaterialBase::computeQpProperties();
      96             : 
      97       62028 :   if (_hys_order[_qp] != _hys_order_old[_qp] && _hys_order[_qp] == 1)
      98        2862 :     computeTurningPoint0Info(_hys_sat_tps[_qp].at(0));
      99             : 
     100       62028 :   computeRelPermQp();
     101       62028 : }
     102             : 
     103             : Real
     104        2878 : PorousFlowHystereticRelativePermeabilityBase::landSat(Real slDel) const
     105             : {
     106        2878 :   const Real a = 1.0 / _s_gr_max - 1.0 / (1.0 - _s_lr);
     107        2878 :   return (1.0 - slDel) / (1.0 + a * (1.0 - slDel));
     108             : }
     109             : 
     110             : void
     111        2878 : PorousFlowHystereticRelativePermeabilityBase::computeTurningPoint0Info(Real tp_sat)
     112             : {
     113             :   // s_gr_tp0 is the Land expression as a function of the zeroth turning point saturation
     114        2878 :   _s_gr_tp0[_qp] = landSat(tp_sat);
     115        2878 : }