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         436 : PorousFlowHystereticRelativePermeabilityBase::validParams()
      14             : {
      15         436 :   InputParameters params = PorousFlowMaterialBase::validParams();
      16        1308 :   params.addRangeCheckedParam<Real>(
      17             :       "S_lr",
      18         872 :       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         872 :   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         872 :   params.addRequiredRangeCheckedParam<Real>(
      33             :       "m", "m > 0 & m < 1", "van Genuchten m parameter.  Suggested value is around 0.9");
      34         872 :   params.addPrivateParam<std::string>("pf_material_type", "relative_permeability");
      35         436 :   params.addPrivateParam<bool>("is_ad", false);
      36         436 :   params.addClassDescription("PorousFlow material that computes relative permeability for 1-phase "
      37             :                              "or 2-phase models that include hysteresis");
      38         436 :   return params;
      39           0 : }
      40             : 
      41         338 : PorousFlowHystereticRelativePermeabilityBase::PorousFlowHystereticRelativePermeabilityBase(
      42         338 :     const InputParameters & parameters)
      43             :   : PorousFlowMaterialBase(parameters),
      44         338 :     _s_lr(getParam<Real>("S_lr")),
      45         676 :     _s_gr_max(getParam<Real>("S_gr_max")),
      46         676 :     _m(getParam<Real>("m")),
      47         676 :     _saturation(_nodal_material
      48         338 :                     ? getMaterialProperty<std::vector<Real>>("PorousFlow_saturation_nodal")
      49         942 :                     : getMaterialProperty<std::vector<Real>>("PorousFlow_saturation_qp")),
      50         338 :     _hys_order(_nodal_material ? getMaterialProperty<unsigned>("PorousFlow_hysteresis_order_nodal")
      51         942 :                                : getMaterialProperty<unsigned>("PorousFlow_hysteresis_order_qp")),
      52         676 :     _hys_order_old(_nodal_material
      53         338 :                        ? getMaterialPropertyOld<unsigned>("PorousFlow_hysteresis_order_nodal")
      54         942 :                        : getMaterialPropertyOld<unsigned>("PorousFlow_hysteresis_order_qp")),
      55         338 :     _hys_sat_tps(
      56         338 :         _nodal_material
      57         338 :             ? getMaterialProperty<std::array<Real, PorousFlowConstants::MAX_HYSTERESIS_ORDER>>(
      58             :                   "PorousFlow_hysteresis_saturation_tps_nodal")
      59         942 :             : getMaterialProperty<std::array<Real, PorousFlowConstants::MAX_HYSTERESIS_ORDER>>(
      60             :                   "PorousFlow_hysteresis_saturation_tps_qp")),
      61         338 :     _relative_permeability(
      62         338 :         _nodal_material ? declareProperty<Real>("PorousFlow_relative_permeability_nodal" + _phase)
      63         640 :                         : declareProperty<Real>("PorousFlow_relative_permeability_qp" + _phase)),
      64         676 :     _drelative_permeability_ds(
      65         338 :         _nodal_material
      66         410 :             ? declarePropertyDerivative<Real>("PorousFlow_relative_permeability_nodal" + _phase,
      67             :                                               _saturation_variable_name)
      68        1244 :             : declarePropertyDerivative<Real>("PorousFlow_relative_permeability_qp" + _phase,
      69             :                                               _saturation_variable_name)),
      70         676 :     _s_gr_tp0(_nodal_material
      71         338 :                   ? declareProperty<Real>("PorousFlow_hysteresis_s_gr_tp0_nodal" + _phase)
      72         978 :                   : 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         338 :   if (_s_gr_max >= 1 - _s_lr)
      79           2 :     paramError("S_gr_max", "Must be less than 1 - S_lr");
      80         336 : }
      81             : 
      82             : void
      83         292 : PorousFlowHystereticRelativePermeabilityBase::initQpStatefulProperties()
      84             : {
      85         292 :   PorousFlowMaterialBase::initQpStatefulProperties();
      86         292 :   if (_hys_order[_qp] >= 1)
      87          16 :     computeTurningPoint0Info(_hys_sat_tps[_qp].at(0));
      88             : 
      89         292 :   computeRelPermQp();
      90         292 : }
      91             : 
      92             : void
      93       43962 : PorousFlowHystereticRelativePermeabilityBase::computeQpProperties()
      94             : {
      95       43962 :   PorousFlowMaterialBase::computeQpProperties();
      96             : 
      97       43962 :   if (_hys_order[_qp] != _hys_order_old[_qp] && _hys_order[_qp] == 1)
      98        1908 :     computeTurningPoint0Info(_hys_sat_tps[_qp].at(0));
      99             : 
     100       43962 :   computeRelPermQp();
     101       43962 : }
     102             : 
     103             : Real
     104        1924 : PorousFlowHystereticRelativePermeabilityBase::landSat(Real slDel) const
     105             : {
     106        1924 :   const Real a = 1.0 / _s_gr_max - 1.0 / (1.0 - _s_lr);
     107        1924 :   return (1.0 - slDel) / (1.0 + a * (1.0 - slDel));
     108             : }
     109             : 
     110             : void
     111        1924 : PorousFlowHystereticRelativePermeabilityBase::computeTurningPoint0Info(Real tp_sat)
     112             : {
     113             :   // s_gr_tp0 is the Land expression as a function of the zeroth turning point saturation
     114        1924 :   _s_gr_tp0[_qp] = landSat(tp_sat);
     115        1924 : }