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 "PorousFlow2PhaseHysPS.h"
      11             : 
      12             : registerMooseObject("PorousFlowApp", PorousFlow2PhaseHysPS);
      13             : 
      14             : InputParameters
      15         364 : PorousFlow2PhaseHysPS::validParams()
      16             : {
      17         364 :   InputParameters params = PorousFlowHystereticCapillaryPressure::validParams();
      18         728 :   params.addRequiredCoupledVar(
      19             :       "phase0_porepressure",
      20             :       "Variable that is the porepressure of phase 0, which is the liquid phase");
      21         728 :   params.addRequiredCoupledVar(
      22             :       "phase1_saturation", "Variable that is the saturation of phase 1, which is the gas phase");
      23         364 :   params.addClassDescription("This Material is used for 2-phase situations.  It calculates the 2 "
      24             :                              "saturations and 2 porepressures, assuming the capillary pressure is "
      25             :                              "hysteretic.  Derivatives of these quantities are also computed");
      26         364 :   return params;
      27           0 : }
      28             : 
      29         284 : PorousFlow2PhaseHysPS::PorousFlow2PhaseHysPS(const InputParameters & parameters)
      30             :   : PorousFlowHystereticCapillaryPressure(parameters),
      31         143 :     _pc(_nodal_material ? declareProperty<Real>("PorousFlow_hysteretic_capillary_pressure_nodal")
      32         425 :                         : declareProperty<Real>("PorousFlow_hysteretic_capillary_pressure_qp")),
      33         284 :     _phase0_porepressure(_nodal_material ? coupledDofValues("phase0_porepressure")
      34         425 :                                          : coupledValue("phase0_porepressure")),
      35         284 :     _phase0_gradp_qp(coupledGradient("phase0_porepressure")),
      36         284 :     _phase0_porepressure_varnum(coupled("phase0_porepressure")),
      37         284 :     _pvar(_dictator.isPorousFlowVariable(_phase0_porepressure_varnum)
      38         284 :               ? _dictator.porousFlowVariableNum(_phase0_porepressure_varnum)
      39             :               : 0),
      40             : 
      41         284 :     _phase1_saturation(_nodal_material ? coupledDofValues("phase1_saturation")
      42         425 :                                        : coupledValue("phase1_saturation")),
      43         284 :     _phase1_grads_qp(coupledGradient("phase1_saturation")),
      44         284 :     _phase1_saturation_varnum(coupled("phase1_saturation")),
      45         284 :     _svar(_dictator.isPorousFlowVariable(_phase1_saturation_varnum)
      46         284 :               ? _dictator.porousFlowVariableNum(_phase1_saturation_varnum)
      47         284 :               : 0)
      48             : {
      49         284 :   if (_num_phases != 2)
      50           4 :     mooseError("The Dictator announces that the number of phases is ",
      51           2 :                _dictator.numPhases(),
      52             :                " whereas PorousFlow2PhaseHysPS can only be used for 2-phase simulation.  The "
      53             :                "Dictator is always watching.");
      54         282 : }
      55             : 
      56             : void
      57         184 : PorousFlow2PhaseHysPS::initQpStatefulProperties()
      58             : {
      59         184 :   PorousFlowHystereticCapillaryPressure::initQpStatefulProperties();
      60         184 :   buildQpPPSS();
      61         184 : }
      62             : 
      63             : void
      64       27468 : PorousFlow2PhaseHysPS::computeQpProperties()
      65             : {
      66             :   // size stuff correctly and prepare the derivative matrices with zeroes
      67       27468 :   PorousFlowHystereticCapillaryPressure::computeQpProperties();
      68             : 
      69       27468 :   buildQpPPSS();
      70       27468 :   const Real dpc = dcapillaryPressureQp(1.0 - _phase1_saturation[_qp]); // d(Pc)/d(S0)
      71             : 
      72       27468 :   if (!_nodal_material)
      73             :   {
      74       13734 :     (*_grads_qp)[_qp][0] = -_phase1_grads_qp[_qp];
      75       13734 :     (*_grads_qp)[_qp][1] = _phase1_grads_qp[_qp];
      76       13734 :     (*_gradp_qp)[_qp][0] = _phase0_gradp_qp[_qp];
      77       13734 :     (*_gradp_qp)[_qp][1] = _phase0_gradp_qp[_qp] - dpc * (*_grads_qp)[_qp][1];
      78             :   }
      79             : 
      80             :   // _porepressure depends on _phase0_porepressure, and its derivative is 1
      81       27468 :   if (_dictator.isPorousFlowVariable(_phase0_porepressure_varnum))
      82             :   {
      83             :     // _phase0_porepressure is a PorousFlow variable
      84       82404 :     for (unsigned phase = 0; phase < _num_phases; ++phase)
      85             :     {
      86       54936 :       (*_dporepressure_dvar)[_qp][phase][_pvar] = 1.0;
      87       54936 :       if (!_nodal_material)
      88       27468 :         (*_dgradp_qp_dgradv)[_qp][phase][_pvar] = 1.0;
      89             :     }
      90             :   }
      91             : 
      92             :   // _saturation is only dependent on _phase1_saturation, and its derivative is +/- 1
      93       27468 :   if (_dictator.isPorousFlowVariable(_phase1_saturation_varnum))
      94             :   {
      95             :     // _phase1_saturation is a PorousFlow variable
      96             :     // _phase1_porepressure depends on saturation through the capillary pressure function
      97       27468 :     (*_dsaturation_dvar)[_qp][0][_svar] = -1.0;
      98       27468 :     (*_dsaturation_dvar)[_qp][1][_svar] = 1.0;
      99       27468 :     (*_dporepressure_dvar)[_qp][1][_svar] = -dpc;
     100             : 
     101       27468 :     if (!_nodal_material)
     102             :     {
     103       13734 :       (*_dgrads_qp_dgradv)[_qp][0][_svar] = -1.0;
     104       13734 :       (*_dgrads_qp_dgradv)[_qp][1][_svar] = 1.0;
     105       13734 :       const Real d2pc_qp = d2capillaryPressureQp(1.0 - _phase1_saturation[_qp]); // d^2(Pc)/dS0^2
     106       13734 :       (*_dgradp_qp_dv)[_qp][1][_svar] = d2pc_qp * (*_grads_qp)[_qp][1];
     107       13734 :       (*_dgradp_qp_dgradv)[_qp][1][_svar] = -dpc;
     108             :     }
     109             :   }
     110       27468 : }
     111             : 
     112             : void
     113       27652 : PorousFlow2PhaseHysPS::buildQpPPSS()
     114             : {
     115       27652 :   _saturation[_qp][0] = 1.0 - _phase1_saturation[_qp];
     116       27652 :   _saturation[_qp][1] = _phase1_saturation[_qp];
     117       27652 :   _pc[_qp] = capillaryPressureQp(1.0 - _phase1_saturation[_qp]);
     118       27652 :   _porepressure[_qp][0] = _phase0_porepressure[_qp];
     119       27652 :   _porepressure[_qp][1] = _phase0_porepressure[_qp] + _pc[_qp];
     120       27652 : }