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