Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://www.mooseframework.org
       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 "HEVPFlowRatePowerLawJ2.h"
      11             : 
      12             : registerMooseObject("TensorMechanicsApp", HEVPFlowRatePowerLawJ2);
      13             : 
      14             : InputParameters
      15          30 : HEVPFlowRatePowerLawJ2::validParams()
      16             : {
      17          30 :   InputParameters params = HEVPFlowRateUOBase::validParams();
      18          60 :   params.addParam<Real>(
      19          60 :       "reference_flow_rate", 0.001, "Reference flow rate for rate dependent flow");
      20          60 :   params.addParam<Real>("flow_rate_exponent", 10.0, "Power law exponent in flow rate equation");
      21          60 :   params.addParam<Real>("flow_rate_tol", 1e3, "Tolerance for flow rate");
      22          30 :   params.addClassDescription(
      23             :       "User object to evaluate power law flow rate and flow direction based on J2");
      24             : 
      25          30 :   return params;
      26           0 : }
      27             : 
      28          15 : HEVPFlowRatePowerLawJ2::HEVPFlowRatePowerLawJ2(const InputParameters & parameters)
      29             :   : HEVPFlowRateUOBase(parameters),
      30          15 :     _ref_flow_rate(getParam<Real>("reference_flow_rate")),
      31          30 :     _flow_rate_exponent(getParam<Real>("flow_rate_exponent")),
      32          45 :     _flow_rate_tol(getParam<Real>("flow_rate_tol"))
      33             : {
      34          15 : }
      35             : 
      36             : bool
      37      533696 : HEVPFlowRatePowerLawJ2::computeValue(unsigned int qp, Real & val) const
      38             : {
      39      533696 :   RankTwoTensor pk2_dev = computePK2Deviatoric(_pk2[qp], _ce[qp]);
      40      533696 :   Real eqv_stress = computeEqvStress(pk2_dev, _ce[qp]);
      41      533696 :   val = std::pow(eqv_stress / _strength[qp], _flow_rate_exponent) * _ref_flow_rate;
      42             : 
      43      533696 :   if (val > _flow_rate_tol)
      44             :   {
      45             : #ifdef DEBUG
      46             :     mooseWarning(
      47             :         "Flow rate greater than ", _flow_rate_tol, " ", val, " ", eqv_stress, " ", _strength[qp]);
      48             : #endif
      49       19776 :     return false;
      50             :   }
      51             :   return true;
      52             : }
      53             : 
      54             : bool
      55      513920 : HEVPFlowRatePowerLawJ2::computeDirection(unsigned int qp, RankTwoTensor & val) const
      56             : {
      57      513920 :   RankTwoTensor pk2_dev = computePK2Deviatoric(_pk2[qp], _ce[qp]);
      58      513920 :   Real eqv_stress = computeEqvStress(pk2_dev, _ce[qp]);
      59             : 
      60             :   val.zero();
      61      513920 :   if (eqv_stress > 0.0)
      62      513840 :     val = 1.5 / eqv_stress * _ce[qp] * pk2_dev * _ce[qp];
      63             : 
      64      513920 :   return true;
      65             : }
      66             : 
      67             : bool
      68      630128 : HEVPFlowRatePowerLawJ2::computeDerivative(unsigned int qp,
      69             :                                           const std::string & coupled_var_name,
      70             :                                           Real & val) const
      71             : {
      72      630128 :   val = 0.0;
      73             : 
      74      630128 :   if (_strength_prop_name == coupled_var_name)
      75             :   {
      76      447920 :     RankTwoTensor pk2_dev = computePK2Deviatoric(_pk2[qp], _ce[qp]);
      77      447920 :     Real eqv_stress = computeEqvStress(pk2_dev, _ce[qp]);
      78      447920 :     val = -_ref_flow_rate * _flow_rate_exponent *
      79      447920 :           std::pow(eqv_stress / _strength[qp], _flow_rate_exponent) / _strength[qp];
      80             :   }
      81             : 
      82      630128 :   return true;
      83             : }
      84             : 
      85             : bool
      86      447920 : HEVPFlowRatePowerLawJ2::computeTensorDerivative(unsigned int qp,
      87             :                                                 const std::string & coupled_var_name,
      88             :                                                 RankTwoTensor & val) const
      89             : {
      90             :   val.zero();
      91             : 
      92      447920 :   if (_pk2_prop_name == coupled_var_name)
      93             :   {
      94      447920 :     RankTwoTensor pk2_dev = computePK2Deviatoric(_pk2[qp], _ce[qp]);
      95      447920 :     Real eqv_stress = computeEqvStress(pk2_dev, _ce[qp]);
      96      447920 :     Real dflowrate_dseqv = _ref_flow_rate * _flow_rate_exponent *
      97      447920 :                            std::pow(eqv_stress / _strength[qp], _flow_rate_exponent - 1.0) /
      98             :                            _strength[qp];
      99             : 
     100      447920 :     RankTwoTensor tau = pk2_dev * _ce[qp];
     101      447920 :     RankTwoTensor dseqv_dpk2dev;
     102             : 
     103             :     dseqv_dpk2dev.zero();
     104      447920 :     if (eqv_stress > 0.0)
     105      447920 :       dseqv_dpk2dev = 1.5 / eqv_stress * tau * _ce[qp];
     106             : 
     107      447920 :     RankTwoTensor ce_inv = _ce[qp].inverse();
     108             : 
     109      447920 :     RankFourTensor dpk2dev_dpk2;
     110     1791680 :     for (const auto i : make_range(Moose::dim))
     111     5375040 :       for (const auto j : make_range(Moose::dim))
     112    16125120 :         for (const auto k : make_range(Moose::dim))
     113    48375360 :           for (const auto l : make_range(Moose::dim))
     114             :           {
     115    36281520 :             dpk2dev_dpk2(i, j, k, l) = 0.0;
     116    36281520 :             if (i == k && j == l)
     117     4031280 :               dpk2dev_dpk2(i, j, k, l) = 1.0;
     118    36281520 :             dpk2dev_dpk2(i, j, k, l) -= ce_inv(i, j) * _ce[qp](k, l) / 3.0;
     119             :           }
     120      895840 :     val = dflowrate_dseqv * dpk2dev_dpk2.transposeMajor() * dseqv_dpk2dev;
     121             :   }
     122      447920 :   return true;
     123             : }
     124             : 
     125             : RankTwoTensor
     126     1943456 : HEVPFlowRatePowerLawJ2::computePK2Deviatoric(const RankTwoTensor & pk2,
     127             :                                              const RankTwoTensor & ce) const
     128             : {
     129     1943456 :   return pk2 - (pk2.doubleContraction(ce) * ce.inverse()) / 3.0;
     130             : }
     131             : 
     132             : Real
     133     1943456 : HEVPFlowRatePowerLawJ2::computeEqvStress(const RankTwoTensor & pk2_dev,
     134             :                                          const RankTwoTensor & ce) const
     135             : {
     136     1943456 :   RankTwoTensor sdev = pk2_dev * ce;
     137     1943456 :   Real val = sdev.doubleContraction(sdev.transpose());
     138     1943456 :   return std::sqrt(1.5 * val);
     139             : }