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