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