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 "PorousFlowPorosityHMBiotModulus.h"
      11             : #include "libmesh/utility.h"
      12             : 
      13             : registerMooseObject("PorousFlowApp", PorousFlowPorosityHMBiotModulus);
      14             : 
      15             : InputParameters
      16         312 : PorousFlowPorosityHMBiotModulus::validParams()
      17             : {
      18         312 :   InputParameters params = PorousFlowPorosity::validParams();
      19         312 :   params.set<bool>("mechanical") = true;
      20         312 :   params.set<bool>("fluid") = true;
      21         624 :   params.addRequiredRangeCheckedParam<Real>("constant_biot_modulus",
      22             :                                             "constant_biot_modulus>0",
      23             :                                             "Biot modulus, which is constant for this Material");
      24         624 :   params.addRequiredRangeCheckedParam<Real>(
      25             :       "constant_fluid_bulk_modulus",
      26             :       "constant_fluid_bulk_modulus>0",
      27             :       "Fluid bulk modulus, which is constant for this Material");
      28         312 :   params.addClassDescription(
      29             :       "This Material calculates the porosity for hydro-mechanical simulations, assuming that the "
      30             :       "Biot modulus and the fluid bulk modulus are both constant.  This is useful for comparing "
      31             :       "with solutions from poroelasticity theory, but is less accurate than PorousFlowPorosity");
      32         312 :   return params;
      33           0 : }
      34             : 
      35         240 : PorousFlowPorosityHMBiotModulus::PorousFlowPorosityHMBiotModulus(const InputParameters & parameters)
      36             :   : PorousFlowPorosity(parameters),
      37         120 :     _porosity_old(_nodal_material ? getMaterialPropertyOld<Real>("PorousFlow_porosity_nodal")
      38         480 :                                   : getMaterialPropertyOld<Real>("PorousFlow_porosity_qp")),
      39         480 :     _biot_modulus(getParam<Real>("constant_biot_modulus")),
      40         480 :     _fluid_bulk_modulus(getParam<Real>("constant_fluid_bulk_modulus")),
      41         480 :     _pf_old(_nodal_material
      42         240 :                 ? getMaterialPropertyOld<Real>("PorousFlow_effective_fluid_pressure_nodal")
      43         480 :                 : getMaterialPropertyOld<Real>("PorousFlow_effective_fluid_pressure_qp")),
      44         480 :     _vol_strain_qp_old(getMaterialPropertyOld<Real>("PorousFlow_total_volumetric_strain_qp")),
      45         480 :     _vol_strain_rate_qp(getMaterialProperty<Real>("PorousFlow_volumetric_strain_rate_qp")),
      46         480 :     _dvol_strain_rate_qp_dvar(getMaterialProperty<std::vector<RealGradient>>(
      47         240 :         "dPorousFlow_volumetric_strain_rate_qp_dvar"))
      48             : {
      49         240 : }
      50             : 
      51             : void
      52     1713600 : PorousFlowPorosityHMBiotModulus::computeQpProperties()
      53             : {
      54             :   // Note that in the following _strain[_qp] is evaluated at q quadpoint
      55             :   // So _porosity_nodal[_qp], which should be the nodal value of porosity
      56             :   // actually uses the strain at a quadpoint.  This
      57             :   // is OK for LINEAR elements, as strain is constant over the element anyway.
      58             : 
      59             :   const unsigned qp_to_use =
      60     1713600 :       (_nodal_material && (_bnd || _strain_at_nearest_qp) ? nearestQP(_qp) : _qp);
      61             : 
      62     1713600 :   const Real denom = (1.0 + _vol_strain_rate_qp[qp_to_use] * _dt + _vol_strain_qp_old[qp_to_use]);
      63     1713600 :   const Real s_p = _porosity_old[_qp] * (1 + _vol_strain_qp_old[qp_to_use]);
      64     1713600 :   _porosity[_qp] = (s_p * std::exp(-((*_pf)[_qp] - _pf_old[_qp]) / _fluid_bulk_modulus) +
      65     1713600 :                     ((*_pf)[_qp] - _pf_old[_qp]) / _biot_modulus +
      66     1713600 :                     _biot * ((*_vol_strain_qp)[qp_to_use] - _vol_strain_qp_old[qp_to_use])) /
      67             :                    denom;
      68             : 
      69     1713600 :   (*_dporosity_dvar)[_qp].resize(_num_var);
      70     8568000 :   for (unsigned int v = 0; v < _num_var; ++v)
      71     6854400 :     (*_dporosity_dvar)[_qp][v] =
      72     6854400 :         (*_dpf_dvar)[_qp][v] *
      73     6854400 :         (-s_p * std::exp(-((*_pf)[_qp] - _pf_old[_qp]) / _fluid_bulk_modulus) /
      74     6854400 :              _fluid_bulk_modulus +
      75     6854400 :          1.0 / _biot_modulus) /
      76             :         denom;
      77             : 
      78     1713600 :   (*_dporosity_dgradvar)[_qp].resize(_num_var);
      79     8568000 :   for (unsigned int v = 0; v < _num_var; ++v)
      80     6854400 :     (*_dporosity_dgradvar)[_qp][v] =
      81     6854400 :         _biot * (*_dvol_strain_qp_dvar)[qp_to_use][v] / denom -
      82     6854400 :         _porosity[_qp] / denom * _dvol_strain_rate_qp_dvar[qp_to_use][v] * _dt;
      83     1713600 : }