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         148 : ADMultiplePowerLawCreepStressUpdate::validParams()
      17             : {
      18         148 :   InputParameters params = ADRadialReturnCreepStressUpdateBase::validParams();
      19         148 :   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         296 :   params.addCoupledVar("temperature", "Coupled temperature");
      26         296 :   params.addRequiredParam<std::vector<Real>>("coefficient",
      27             :                                              "Leading coefficient in power-law equation");
      28         296 :   params.addRequiredParam<std::vector<Real>>("n_exponent",
      29             :                                              "Exponent on effective stress in power-law equation");
      30         296 :   params.addRequiredParam<std::vector<Real>>("m_exponent",
      31             :                                              "Exponent on time in power-law equation");
      32         296 :   params.addRequiredParam<std::vector<Real>>("stress_thresholds",
      33             :                                              "Stress intervals to switch creep behavior");
      34         296 :   params.addRequiredParam<std::vector<Real>>("activation_energy", "Activation energy");
      35         296 :   params.addParam<Real>("gas_constant", 8.3143, "Universal gas constant");
      36         296 :   params.addParam<Real>("start_time", 0.0, "Start time (if not zero)");
      37         148 :   return params;
      38           0 : }
      39             : 
      40         110 : ADMultiplePowerLawCreepStressUpdate::ADMultiplePowerLawCreepStressUpdate(
      41         110 :     const InputParameters & parameters)
      42             :   : ADRadialReturnCreepStressUpdateBase(parameters),
      43         110 :     _temperature(isParamValid("temperature") ? &adCoupledValue("temperature") : nullptr),
      44         220 :     _coefficient(getParam<std::vector<Real>>("coefficient")),
      45         220 :     _n_exponent(getParam<std::vector<Real>>("n_exponent")),
      46         220 :     _m_exponent(getParam<std::vector<Real>>("m_exponent")),
      47         220 :     _stress_thresholds(getParam<std::vector<Real>>("stress_thresholds")),
      48         220 :     _activation_energy(getParam<std::vector<Real>>("activation_energy")),
      49         220 :     _gas_constant(getParam<Real>("gas_constant")),
      50         220 :     _start_time(getParam<Real>("start_time")),
      51         110 :     _exponential(1),
      52         110 :     _exp_time(0),
      53         110 :     _stress_index(0),
      54         220 :     _number_of_models(_m_exponent.size())
      55             : {
      56             : 
      57         110 :   if (!std::is_sorted(_stress_thresholds.begin(), _stress_thresholds.end()))
      58           2 :     paramError("stress_thresholds",
      59             :                "Stress thresholds input must be provided in increasing ordered");
      60             : 
      61         108 :   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         108 :   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         288 :   for (unsigned int i = 0; i < _number_of_models; i++)
      73         180 :     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         108 : }
      78             : 
      79             : void
      80     4495896 : ADMultiplePowerLawCreepStressUpdate::computeStressInitialize(
      81             :     const ADReal & effective_trial_stress, const ADRankFourTensor & elasticity_tensor)
      82             : {
      83     4495896 :   ADRadialReturnStressUpdate::computeStressInitialize(effective_trial_stress, elasticity_tensor);
      84             : 
      85     4495896 :   _stress_index = stressIndex(effective_trial_stress);
      86             : 
      87     4495896 :   if (_temperature)
      88             :     _exponential =
      89           0 :         std::exp(-_activation_energy[_stress_index] / (_gas_constant * (*_temperature)[_qp]));
      90             : 
      91     4495896 :   _exp_time = std::pow(_t - _start_time, _m_exponent[_stress_index]);
      92     4495896 : }
      93             : 
      94             : ADReal
      95    27641968 : ADMultiplePowerLawCreepStressUpdate::computeResidual(const ADReal & effective_trial_stress,
      96             :                                                      const ADReal & scalar)
      97             : {
      98    27641968 :   const ADReal stress_delta = effective_trial_stress - _three_shear_modulus * scalar;
      99    27641968 :   const ADReal creep_rate = _coefficient[_stress_index] *
     100    55283936 :                             std::pow(stress_delta, _n_exponent[_stress_index]) * _exponential *
     101    27641968 :                             _exp_time;
     102    27641968 :   return creep_rate * _dt - scalar;
     103             : }
     104             : 
     105             : ADReal
     106    27641968 : ADMultiplePowerLawCreepStressUpdate::computeDerivative(const ADReal & effective_trial_stress,
     107             :                                                        const ADReal & scalar)
     108             : {
     109    27641968 :   const ADReal stress_delta = effective_trial_stress - _three_shear_modulus * scalar;
     110             :   const ADReal creep_rate_derivative =
     111    27641968 :       -_coefficient[_stress_index] * _three_shear_modulus * _n_exponent[_stress_index] *
     112    27641968 :       std::pow(stress_delta, _n_exponent[_stress_index] - 1.0) * _exponential * _exp_time;
     113    55283936 :   return creep_rate_derivative * _dt - 1.0;
     114             : }
     115             : 
     116             : Real
     117      225216 : ADMultiplePowerLawCreepStressUpdate::computeStrainEnergyRateDensity(
     118             :     const ADMaterialProperty<RankTwoTensor> & stress,
     119             :     const ADMaterialProperty<RankTwoTensor> & strain_rate)
     120             : {
     121      225216 :   if (_n_exponent[_stress_index] <= 1)
     122             :     return 0.0;
     123             : 
     124      225216 :   Real creep_factor = _n_exponent[_stress_index] / (_n_exponent[_stress_index] + 1);
     125             : 
     126      450432 :   return MetaPhysicL::raw_value(creep_factor * stress[_qp].doubleContraction((strain_rate)[_qp]));
     127             : }
     128             : 
     129             : bool
     130     4495896 : ADMultiplePowerLawCreepStressUpdate::substeppingCapabilityEnabled()
     131             : {
     132     4495896 :   return this->_use_substepping != ADRadialReturnStressUpdate::SubsteppingType::NONE;
     133             : }
     134             : 
     135             : std::size_t
     136     4495896 : ADMultiplePowerLawCreepStressUpdate::stressIndex(const ADReal & effective_trial_stress)
     137             : {
     138             :   // Check the correct model for *this* radial return
     139             :   std::size_t i = 0;
     140     6631580 :   while (i < _number_of_models - 1 && effective_trial_stress > _stress_thresholds[i + 1])
     141             :     ++i;
     142             : 
     143     4495896 :   return i;
     144             : }