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 "PorousFlow2PhaseHysPP.h"
      11             : 
      12             : registerMooseObject("PorousFlowApp", PorousFlow2PhaseHysPP);
      13             : 
      14             : InputParameters
      15         364 : PorousFlow2PhaseHysPP::validParams()
      16             : {
      17         364 :   InputParameters params = PorousFlowHystereticCapillaryPressure::validParams();
      18         728 :   params.addRequiredCoupledVar("phase0_porepressure",
      19             :                                "Variable that is the porepressure of phase 0.  This is assumed to "
      20             :                                "be the liquid phase.  It will be <= phase1_porepressure.");
      21         728 :   params.addRequiredCoupledVar("phase1_porepressure",
      22             :                                "Variable that is the porepressure of phase 1 (the gas phase)");
      23         728 :   params.addParam<unsigned>(
      24         728 :       "liquid_phase", 0, "Phase number of the liquid phase (more precisely, the phase of phase0)");
      25         364 :   params.addClassDescription(
      26             :       "This Material is used for 2-phase situations.  It calculates the 2 saturations given the 2 "
      27             :       "porepressures, assuming the capillary pressure is hysteretic.  Derivatives of these "
      28             :       "quantities are also computed");
      29         364 :   return params;
      30           0 : }
      31             : 
      32         284 : PorousFlow2PhaseHysPP::PorousFlow2PhaseHysPP(const InputParameters & parameters)
      33             :   : PorousFlowHystereticCapillaryPressure(parameters),
      34         143 :     _pc(_nodal_material ? declareProperty<Real>("PorousFlow_hysteretic_capillary_pressure_nodal")
      35         425 :                         : declareProperty<Real>("PorousFlow_hysteretic_capillary_pressure_qp")),
      36         284 :     _phase0_porepressure(_nodal_material ? coupledDofValues("phase0_porepressure")
      37         425 :                                          : coupledValue("phase0_porepressure")),
      38         284 :     _phase0_gradp_qp(coupledGradient("phase0_porepressure")),
      39         284 :     _phase0_porepressure_varnum(coupled("phase0_porepressure")),
      40         284 :     _p0var(_dictator.isPorousFlowVariable(_phase0_porepressure_varnum)
      41         284 :                ? _dictator.porousFlowVariableNum(_phase0_porepressure_varnum)
      42             :                : 0),
      43             : 
      44         284 :     _phase1_porepressure(_nodal_material ? coupledDofValues("phase1_porepressure")
      45         425 :                                          : coupledValue("phase1_porepressure")),
      46         284 :     _phase1_gradp_qp(coupledGradient("phase1_porepressure")),
      47         284 :     _phase1_porepressure_varnum(coupled("phase1_porepressure")),
      48         284 :     _p1var(_dictator.isPorousFlowVariable(_phase1_porepressure_varnum)
      49         284 :                ? _dictator.porousFlowVariableNum(_phase1_porepressure_varnum)
      50         284 :                : 0)
      51             : {
      52         284 :   if (_num_phases != 2)
      53           4 :     mooseError("The Dictator announces that the number of phases is ",
      54           2 :                _dictator.numPhases(),
      55             :                " whereas PorousFlow2PhaseHysPP can only be used for 2-phase simulation.  When you "
      56             :                "have efficient government, you have dictatorship.");
      57         282 : }
      58             : 
      59             : void
      60         184 : PorousFlow2PhaseHysPP::initQpStatefulProperties()
      61             : {
      62         184 :   PorousFlowHystereticCapillaryPressure::initQpStatefulProperties();
      63         184 :   buildQpPPSS();
      64         184 : }
      65             : 
      66             : void
      67       30768 : PorousFlow2PhaseHysPP::computeQpProperties()
      68             : {
      69             :   // size stuff correctly and prepare the derivative matrices with zeroes
      70       30768 :   PorousFlowHystereticCapillaryPressure::computeQpProperties();
      71             : 
      72       30768 :   buildQpPPSS();
      73       30768 :   const Real pc = _pc[_qp];                // >= 0
      74       30768 :   const Real ds = dliquidSaturationQp(pc); // dS/d(pc)
      75             : 
      76       30768 :   if (!_nodal_material)
      77             :   {
      78       15384 :     (*_gradp_qp)[_qp][0] = _phase0_gradp_qp[_qp];
      79       15384 :     (*_gradp_qp)[_qp][1] = _phase1_gradp_qp[_qp];
      80       15384 :     (*_grads_qp)[_qp][0] = ds * ((*_gradp_qp)[_qp][1] - (*_gradp_qp)[_qp][0]);
      81       15384 :     (*_grads_qp)[_qp][1] = -(*_grads_qp)[_qp][0];
      82             :   }
      83             : 
      84             :   // the derivatives of porepressure with respect to porepressure
      85             :   // remain fixed (at unity) throughout the simulation
      86       30768 :   if (_dictator.isPorousFlowVariable(_phase0_porepressure_varnum))
      87             :   {
      88       30768 :     (*_dporepressure_dvar)[_qp][0][_p0var] = 1.0;
      89       30768 :     if (!_nodal_material)
      90       15384 :       (*_dgradp_qp_dgradv)[_qp][0][_p0var] = 1.0;
      91             :   }
      92       30768 :   if (_dictator.isPorousFlowVariable(_phase1_porepressure_varnum))
      93             :   {
      94       30768 :     (*_dporepressure_dvar)[_qp][1][_p1var] = 1.0;
      95       30768 :     if (!_nodal_material)
      96       15384 :       (*_dgradp_qp_dgradv)[_qp][1][_p1var] = 1.0;
      97             :   }
      98             : 
      99       30768 :   if (_dictator.isPorousFlowVariable(_phase0_porepressure_varnum))
     100             :   {
     101       30768 :     (*_dsaturation_dvar)[_qp][0][_p0var] = -ds;
     102       30768 :     (*_dsaturation_dvar)[_qp][1][_p0var] = ds;
     103             :   }
     104       30768 :   if (_dictator.isPorousFlowVariable(_phase1_porepressure_varnum))
     105             :   {
     106       30768 :     (*_dsaturation_dvar)[_qp][0][_p1var] = ds;
     107       30768 :     (*_dsaturation_dvar)[_qp][1][_p1var] = -ds;
     108             :   }
     109             : 
     110       30768 :   _pc[_qp] = _phase1_porepressure[_qp] - _phase0_porepressure[_qp]; // this is >= 0
     111       30768 :   if (!_nodal_material)
     112             :   {
     113       15384 :     const Real d2s_qp = d2liquidSaturationQp(pc); // d^2(S_qp)/d(pc_qp)^2
     114       15384 :     if (_dictator.isPorousFlowVariable(_phase0_porepressure_varnum))
     115             :     {
     116       15384 :       (*_dgrads_qp_dgradv)[_qp][0][_p0var] = -ds;
     117       15384 :       (*_dgrads_qp_dv)[_qp][0][_p0var] = -d2s_qp * (_phase1_gradp_qp[_qp] - _phase0_gradp_qp[_qp]);
     118       15384 :       (*_dgrads_qp_dgradv)[_qp][1][_p0var] = ds;
     119       15384 :       (*_dgrads_qp_dv)[_qp][1][_p0var] = d2s_qp * (_phase1_gradp_qp[_qp] - _phase0_gradp_qp[_qp]);
     120             :     }
     121       15384 :     if (_dictator.isPorousFlowVariable(_phase1_porepressure_varnum))
     122             :     {
     123       15384 :       (*_dgrads_qp_dgradv)[_qp][0][_p1var] = ds;
     124       15384 :       (*_dgrads_qp_dv)[_qp][0][_p1var] = d2s_qp * (_phase1_gradp_qp[_qp] - _phase0_gradp_qp[_qp]);
     125       15384 :       (*_dgrads_qp_dgradv)[_qp][1][_p1var] = -ds;
     126       15384 :       (*_dgrads_qp_dv)[_qp][1][_p1var] = -d2s_qp * (_phase1_gradp_qp[_qp] - _phase0_gradp_qp[_qp]);
     127             :     }
     128             :   }
     129       30768 : }
     130             : 
     131             : void
     132       30952 : PorousFlow2PhaseHysPP::buildQpPPSS()
     133             : {
     134       30952 :   _porepressure[_qp][0] = _phase0_porepressure[_qp];
     135       30952 :   _porepressure[_qp][1] = _phase1_porepressure[_qp];
     136       30952 :   _pc[_qp] = _phase1_porepressure[_qp] - _phase0_porepressure[_qp]; // this is >= 0
     137       30952 :   const Real sat = liquidSaturationQp(_pc[_qp]);
     138       30952 :   _saturation[_qp][0] = sat;
     139       30952 :   _saturation[_qp][1] = 1.0 - sat;
     140       30952 : }