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 "ADMultiplePowerLawCreepStressUpdate.h"
      11             : #include <algorithm>
      12             : 
      13             : registerMooseObject("SolidMechanicsApp", ADMultiplePowerLawCreepStressUpdate);
      14             : 
      15             : InputParameters
      16          98 : ADMultiplePowerLawCreepStressUpdate::validParams()
      17             : {
      18          98 :   InputParameters params = ADRadialReturnCreepStressUpdateBase::validParams();
      19          98 :   params.addClassDescription(
      20             :       "This class uses the stress update material in a radial return isotropic power law creep "
      21             :       "model. This class can be used in conjunction with other creep and plasticity materials "
      22             :       "for more complex simulations.");
      23             : 
      24             :   // Linear strain hardening parameters
      25         196 :   params.addCoupledVar("temperature", "Coupled temperature");
      26         196 :   params.addRequiredParam<std::vector<Real>>("coefficient",
      27             :                                              "Leading coefficient in power-law equation");
      28         196 :   params.addRequiredParam<std::vector<Real>>("n_exponent",
      29             :                                              "Exponent on effective stress in power-law equation");
      30         196 :   params.addRequiredParam<std::vector<Real>>("m_exponent",
      31             :                                              "Exponent on time in power-law equation");
      32         196 :   params.addRequiredParam<std::vector<Real>>("stress_thresholds",
      33             :                                              "Stress intervals to switch creep behavior");
      34         196 :   params.addRequiredParam<std::vector<Real>>("activation_energy", "Activation energy");
      35         196 :   params.addParam<Real>("gas_constant", 8.3143, "Universal gas constant");
      36         196 :   params.addParam<Real>("start_time", 0.0, "Start time (if not zero)");
      37          98 :   return params;
      38           0 : }
      39             : 
      40          73 : ADMultiplePowerLawCreepStressUpdate::ADMultiplePowerLawCreepStressUpdate(
      41          73 :     const InputParameters & parameters)
      42             :   : ADRadialReturnCreepStressUpdateBase(parameters),
      43          73 :     _temperature(isParamValid("temperature") ? &adCoupledValue("temperature") : nullptr),
      44         146 :     _coefficient(getParam<std::vector<Real>>("coefficient")),
      45         146 :     _n_exponent(getParam<std::vector<Real>>("n_exponent")),
      46         146 :     _m_exponent(getParam<std::vector<Real>>("m_exponent")),
      47         146 :     _stress_thresholds(getParam<std::vector<Real>>("stress_thresholds")),
      48         146 :     _activation_energy(getParam<std::vector<Real>>("activation_energy")),
      49         146 :     _gas_constant(getParam<Real>("gas_constant")),
      50         146 :     _start_time(getParam<Real>("start_time")),
      51          73 :     _exponential(1),
      52          73 :     _exp_time(0),
      53          73 :     _stress_index(0),
      54         146 :     _number_of_models(_m_exponent.size())
      55             : {
      56             : 
      57          73 :   if (!std::is_sorted(_stress_thresholds.begin(), _stress_thresholds.end()))
      58           1 :     paramError("stress_thresholds",
      59             :                "Stress thresholds input must be provided in increasing ordered");
      60             : 
      61          72 :   if (_coefficient.size() != _n_exponent.size() || _coefficient.size() != _m_exponent.size() ||
      62             :       _coefficient.size() != _activation_energy.size())
      63           0 :     paramError(
      64             :         "n_exponent",
      65             :         "Inputs to ADMultiplePowerLawCreepStressUpdate creep models must have the same size");
      66             : 
      67          72 :   if (_number_of_models != _stress_thresholds.size() - 1)
      68           0 :     paramError("stress_thresholds",
      69             :                "The number of creep models must be the number of stress thresholds minute. That "
      70             :                "is, one creep model will be valid between two stress thresholds");
      71             : 
      72         192 :   for (unsigned int i = 0; i < _number_of_models; i++)
      73         120 :     if (_start_time < this->_app.getStartTime() && (std::trunc(_m_exponent[i]) != _m_exponent[i]))
      74           0 :       paramError("start_time",
      75             :                  "Start time must be equal to or greater than the Executioner start_time if a "
      76             :                  "non-integer m_exponent is used");
      77          72 : }
      78             : 
      79             : void
      80     3461712 : ADMultiplePowerLawCreepStressUpdate::computeStressInitialize(
      81             :     const ADReal & effective_trial_stress, const ADRankFourTensor & elasticity_tensor)
      82             : {
      83             :   using std::exp;
      84             : 
      85     3461712 :   ADRadialReturnStressUpdate::computeStressInitialize(effective_trial_stress, elasticity_tensor);
      86             : 
      87     3461712 :   _stress_index = stressIndex(effective_trial_stress);
      88             : 
      89     3461712 :   if (_temperature)
      90           0 :     _exponential = exp(-_activation_energy[_stress_index] / (_gas_constant * (*_temperature)[_qp]));
      91             : 
      92     3461712 :   _exp_time = pow(_t - _start_time, _m_exponent[_stress_index]);
      93     3461712 : }
      94             : 
      95             : ADReal
      96    21295546 : ADMultiplePowerLawCreepStressUpdate::computeResidual(const ADReal & effective_trial_stress,
      97             :                                                      const ADReal & scalar)
      98             : {
      99             :   using std::pow;
     100             : 
     101    21295546 :   const ADReal stress_delta = effective_trial_stress - _three_shear_modulus * scalar;
     102    21295546 :   const ADReal creep_rate = _coefficient[_stress_index] *
     103    42591092 :                             pow(stress_delta, _n_exponent[_stress_index]) * _exponential *
     104    21295546 :                             _exp_time;
     105    21295546 :   return creep_rate * _dt - scalar;
     106             : }
     107             : 
     108             : ADReal
     109    21295546 : ADMultiplePowerLawCreepStressUpdate::computeDerivative(const ADReal & effective_trial_stress,
     110             :                                                        const ADReal & scalar)
     111             : {
     112             :   using std::pow;
     113             : 
     114    21295546 :   const ADReal stress_delta = effective_trial_stress - _three_shear_modulus * scalar;
     115             :   const ADReal creep_rate_derivative =
     116    21295546 :       -_coefficient[_stress_index] * _three_shear_modulus * _n_exponent[_stress_index] *
     117    21295546 :       pow(stress_delta, _n_exponent[_stress_index] - 1.0) * _exponential * _exp_time;
     118    42591092 :   return creep_rate_derivative * _dt - 1.0;
     119             : }
     120             : 
     121             : Real
     122      172692 : ADMultiplePowerLawCreepStressUpdate::computeStrainEnergyRateDensity(
     123             :     const ADMaterialProperty<RankTwoTensor> & stress,
     124             :     const ADMaterialProperty<RankTwoTensor> & strain_rate)
     125             : {
     126      172692 :   if (_n_exponent[_stress_index] <= 1)
     127             :     return 0.0;
     128             : 
     129      172692 :   Real creep_factor = _n_exponent[_stress_index] / (_n_exponent[_stress_index] + 1);
     130             : 
     131      345384 :   return MetaPhysicL::raw_value(creep_factor * stress[_qp].doubleContraction((strain_rate)[_qp]));
     132             : }
     133             : 
     134             : bool
     135     3461712 : ADMultiplePowerLawCreepStressUpdate::substeppingCapabilityEnabled()
     136             : {
     137     3461712 :   return this->_use_substepping != ADRadialReturnStressUpdate::SubsteppingType::NONE;
     138             : }
     139             : 
     140             : std::size_t
     141     3461712 : ADMultiplePowerLawCreepStressUpdate::stressIndex(const ADReal & effective_trial_stress)
     142             : {
     143             :   // Check the correct model for *this* radial return
     144             :   std::size_t i = 0;
     145     5105930 :   while (i < _number_of_models - 1 && effective_trial_stress > _stress_thresholds[i + 1])
     146             :     ++i;
     147             : 
     148     3461712 :   return i;
     149             : }