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 "CompositePowerLawCreepStressUpdate.h"
      11             : #include <cmath>
      12             : 
      13             : registerMooseObject("SolidMechanicsApp", CompositePowerLawCreepStressUpdate);
      14             : registerMooseObject("SolidMechanicsApp", ADCompositePowerLawCreepStressUpdate);
      15             : 
      16             : template <bool is_ad>
      17             : InputParameters
      18         146 : CompositePowerLawCreepStressUpdateTempl<is_ad>::validParams()
      19             : {
      20         146 :   InputParameters params = RadialReturnCreepStressUpdateBaseTempl<is_ad>::validParams();
      21         146 :   params.addClassDescription(
      22             :       "This class uses the stress update material in a radial return isotropic power law creep "
      23             :       "model. This class can be used in conjunction with other creep and plasticity materials "
      24             :       "for more complex simulations. This class is an extension to include multi-phase "
      25             :       "capability.");
      26             : 
      27             :   // Linear strain hardening parameters
      28         292 :   params.addRequiredCoupledVar("temperature", "Coupled temperature");
      29         292 :   params.addRequiredParam<std::vector<Real>>(
      30             :       "coefficient",
      31             :       "a vector of leading coefficient / Dorn Constant in power-law equation for each material.");
      32         292 :   params.addRequiredParam<std::vector<Real>>(
      33             :       "n_exponent",
      34             :       "a vector of Exponent on effective stress in power-law equation for each material");
      35         292 :   params.addParam<Real>("m_exponent", 0.0, "Exponent on time in power-law equation");
      36         292 :   params.addRequiredParam<std::vector<Real>>("activation_energy",
      37             :                                              "a vector of Activation energy for Arrhenius-type "
      38             :                                              "equation of the Dorn Constant for each material");
      39         292 :   params.addParam<Real>("gas_constant", 8.3143, "Universal gas constant");
      40         292 :   params.addParam<Real>("start_time", 0.0, "Start time (if not zero)");
      41         292 :   params.addRequiredParam<std::vector<MaterialPropertyName>>(
      42             :       "switching_functions", "a vector of switching functions for each material");
      43         146 :   return params;
      44           0 : }
      45             : 
      46             : template <bool is_ad>
      47         111 : CompositePowerLawCreepStressUpdateTempl<is_ad>::CompositePowerLawCreepStressUpdateTempl(
      48             :     const InputParameters & parameters)
      49             :   : RadialReturnCreepStressUpdateBaseTempl<is_ad>(parameters),
      50         333 :     _temperature(this->isParamValid("temperature")
      51         222 :                      ? &this->template coupledGenericValue<is_ad>("temperature")
      52             :                      : nullptr),
      53         222 :     _coefficient(this->template getParam<std::vector<Real>>("coefficient")),
      54         222 :     _n_exponent(this->template getParam<std::vector<Real>>("n_exponent")),
      55         222 :     _m_exponent(this->template getParam<Real>("m_exponent")),
      56         222 :     _activation_energy(this->template getParam<std::vector<Real>>("activation_energy")),
      57         222 :     _gas_constant(this->template getParam<Real>("gas_constant")),
      58         222 :     _start_time(this->template getParam<Real>("start_time")),
      59         333 :     _switching_func_names(
      60         111 :         this->template getParam<std::vector<MaterialPropertyName>>("switching_functions"))
      61             : 
      62             : {
      63         111 :   _num_materials = _switching_func_names.size();
      64         111 :   if (_n_exponent.size() != _num_materials)
      65           2 :     this->paramError("n_exponent", "n exponent must be equal to the number of switching functions");
      66             : 
      67         109 :   if (_coefficient.size() != _num_materials)
      68           2 :     this->paramError("coefficient",
      69             :                      "number of Dorn constant must be equal to the number of switching functions");
      70             : 
      71         107 :   if (_activation_energy.size() != _num_materials)
      72           2 :     this->paramError("activation_energy",
      73             :                      "activation energy must be equal to the number of swithing functions");
      74             : 
      75         105 :   if (_start_time < this->_app.getStartTime() && (std::trunc(_m_exponent) != _m_exponent))
      76           0 :     this->paramError("start_time",
      77             :                      "Start time must be equal to or greater than the Executioner start_time if a "
      78             :                      "non-integer m_exponent is used");
      79         105 :   _switchingFunc.resize(_num_materials);
      80             :   // set switching functions material properties for each phase
      81         315 :   for (unsigned int i = 0; i < _num_materials; ++i)
      82             :   {
      83         210 :     _switchingFunc[i] =
      84         210 :         &this->template getGenericMaterialProperty<Real, is_ad>(_switching_func_names[i]);
      85             :   }
      86         105 : }
      87             : 
      88             : template <bool is_ad>
      89             : void
      90     1786464 : CompositePowerLawCreepStressUpdateTempl<is_ad>::computeStressInitialize(
      91             :     const GenericReal<is_ad> & effective_trial_stress,
      92             :     const GenericRankFourTensor<is_ad> & elasticity_tensor)
      93             : {
      94             :   using std::pow;
      95             : 
      96           0 :   RadialReturnStressUpdateTempl<is_ad>::computeStressInitialize(effective_trial_stress,
      97             :                                                                 elasticity_tensor);
      98     1786464 :   _exp_time = pow(_t - _start_time, _m_exponent);
      99     1786464 : }
     100             : 
     101             : template <bool is_ad>
     102             : template <typename ScalarType>
     103             : ScalarType
     104    19604789 : CompositePowerLawCreepStressUpdateTempl<is_ad>::computeResidualInternal(
     105             :     const GenericReal<is_ad> & effective_trial_stress, const ScalarType & scalar)
     106             : {
     107             :   using std::pow, std::exp;
     108             : 
     109    19604789 :   const ScalarType stress_delta = effective_trial_stress - _three_shear_modulus * scalar;
     110           0 :   ScalarType creep_rate = 0.0;
     111    58814367 :   for (const auto n : make_range(_num_materials))
     112             :   {
     113    39209578 :     creep_rate += _coefficient[n] * pow(stress_delta, _n_exponent[n]) * (*_switchingFunc[n])[_qp] *
     114    39209578 :                   exp(-_activation_energy[n] / (_gas_constant * (*_temperature)[_qp])) * _exp_time;
     115             :   }
     116    19604789 :   return creep_rate * _dt - scalar;
     117             : }
     118             : 
     119             : template <bool is_ad>
     120             : GenericReal<is_ad>
     121    19604789 : CompositePowerLawCreepStressUpdateTempl<is_ad>::computeDerivative(
     122             :     const GenericReal<is_ad> & effective_trial_stress, const GenericReal<is_ad> & scalar)
     123             : {
     124             :   using std::pow, std::exp;
     125             : 
     126    19604789 :   const GenericReal<is_ad> stress_delta = effective_trial_stress - _three_shear_modulus * scalar;
     127           0 :   GenericReal<is_ad> creep_rate_derivative = 0.0;
     128    58814367 :   for (const auto n : make_range(_num_materials))
     129             :   {
     130    39209578 :     creep_rate_derivative += -_coefficient[n] * _three_shear_modulus * _n_exponent[n] *
     131    39209578 :                              pow(stress_delta, _n_exponent[n] - 1.0) * (*_switchingFunc[n])[_qp] *
     132    39209578 :                              exp(-_activation_energy[n] / (_gas_constant * (*_temperature)[_qp])) *
     133    39209578 :                              _exp_time;
     134             :   }
     135    19604789 :   return creep_rate_derivative * _dt - 1.0;
     136             : }
     137             : 
     138             : template <bool is_ad>
     139             : Real
     140           0 : CompositePowerLawCreepStressUpdateTempl<is_ad>::computeStrainEnergyRateDensity(
     141             :     const GenericMaterialProperty<RankTwoTensor, is_ad> & stress,
     142             :     const GenericMaterialProperty<RankTwoTensor, is_ad> & strain_rate)
     143             : {
     144           0 :   GenericReal<is_ad> interpolated_exponent = 0.0;
     145           0 :   for (unsigned int n = 0; n < _num_materials; ++n)
     146             :   {
     147           0 :     interpolated_exponent += (_n_exponent[n] / (_n_exponent[n] + 1)) * (*_switchingFunc[n])[_qp];
     148             :   }
     149           0 :   return MetaPhysicL::raw_value(interpolated_exponent *
     150           0 :                                 stress[_qp].doubleContraction((strain_rate)[_qp]));
     151             : }
     152             : 
     153             : template <bool is_ad>
     154             : void
     155     1786464 : CompositePowerLawCreepStressUpdateTempl<is_ad>::computeStressFinalize(
     156             :     const GenericRankTwoTensor<is_ad> & plastic_strain_increment)
     157             : {
     158     1786464 :   _creep_strain[_qp] += plastic_strain_increment;
     159     1786464 : }
     160             : 
     161             : template <bool is_ad>
     162             : void
     163     1786464 : CompositePowerLawCreepStressUpdateTempl<is_ad>::resetIncrementalMaterialProperties()
     164             : {
     165     1786464 :   _creep_strain[_qp] = _creep_strain_old[_qp];
     166     1786464 : }
     167             : 
     168             : template <bool is_ad>
     169             : bool
     170     1786464 : CompositePowerLawCreepStressUpdateTempl<is_ad>::substeppingCapabilityEnabled()
     171             : {
     172     1786464 :   return this->_use_substepping != RadialReturnStressUpdateTempl<is_ad>::SubsteppingType::NONE;
     173             : }
     174             : 
     175             : template class CompositePowerLawCreepStressUpdateTempl<false>;
     176             : template class CompositePowerLawCreepStressUpdateTempl<true>;
     177             : template Real
     178             : CompositePowerLawCreepStressUpdateTempl<false>::computeResidualInternal<Real>(const Real &,
     179             :                                                                               const Real &);
     180             : template ADReal
     181             : CompositePowerLawCreepStressUpdateTempl<true>::computeResidualInternal<ADReal>(const ADReal &,
     182             :                                                                                const ADReal &);
     183             : template ChainedReal
     184             : CompositePowerLawCreepStressUpdateTempl<false>::computeResidualInternal<ChainedReal>(
     185             :     const Real &, const ChainedReal &);
     186             : template ChainedADReal
     187             : CompositePowerLawCreepStressUpdateTempl<true>::computeResidualInternal<ChainedADReal>(
     188             :     const ADReal &, const ChainedADReal &);