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         126 : CompositePowerLawCreepStressUpdateTempl<is_ad>::validParams()
      19             : {
      20         126 :   InputParameters params = RadialReturnCreepStressUpdateBaseTempl<is_ad>::validParams();
      21         126 :   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         252 :   params.addRequiredCoupledVar("temperature", "Coupled temperature");
      29         252 :   params.addRequiredParam<std::vector<Real>>(
      30             :       "coefficient",
      31             :       "a vector of leading coefficient / Dorn Constant in power-law equation for each material.");
      32         252 :   params.addRequiredParam<std::vector<Real>>(
      33             :       "n_exponent",
      34             :       "a vector of Exponent on effective stress in power-law equation for each material");
      35         252 :   params.addParam<Real>("m_exponent", 0.0, "Exponent on time in power-law equation");
      36         252 :   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         252 :   params.addParam<Real>("gas_constant", 8.3143, "Universal gas constant");
      40         252 :   params.addParam<Real>("start_time", 0.0, "Start time (if not zero)");
      41         252 :   params.addRequiredParam<std::vector<MaterialPropertyName>>(
      42             :       "switching_functions", "a vector of switching functions for each material");
      43         126 :   return params;
      44           0 : }
      45             : 
      46             : template <bool is_ad>
      47          96 : CompositePowerLawCreepStressUpdateTempl<is_ad>::CompositePowerLawCreepStressUpdateTempl(
      48             :     const InputParameters & parameters)
      49             :   : RadialReturnCreepStressUpdateBaseTempl<is_ad>(parameters),
      50         288 :     _temperature(this->isParamValid("temperature")
      51         192 :                      ? &this->template coupledGenericValue<is_ad>("temperature")
      52             :                      : nullptr),
      53         192 :     _coefficient(this->template getParam<std::vector<Real>>("coefficient")),
      54         192 :     _n_exponent(this->template getParam<std::vector<Real>>("n_exponent")),
      55         192 :     _m_exponent(this->template getParam<Real>("m_exponent")),
      56         192 :     _activation_energy(this->template getParam<std::vector<Real>>("activation_energy")),
      57         192 :     _gas_constant(this->template getParam<Real>("gas_constant")),
      58         192 :     _start_time(this->template getParam<Real>("start_time")),
      59         192 :     _switching_func_names(
      60         192 :         this->template getParam<std::vector<MaterialPropertyName>>("switching_functions"))
      61             : 
      62             : {
      63          96 :   _num_materials = _switching_func_names.size();
      64          96 :   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          94 :   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          92 :   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          90 :   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          90 :   _switchingFunc.resize(_num_materials);
      80             :   // set switching functions material properties for each phase
      81         270 :   for (unsigned int i = 0; i < _num_materials; ++i)
      82             :   {
      83         180 :     _switchingFunc[i] =
      84         180 :         &this->template getGenericMaterialProperty<Real, is_ad>(_switching_func_names[i]);
      85             :   }
      86          90 : }
      87             : 
      88             : template <bool is_ad>
      89             : void
      90     2176512 : CompositePowerLawCreepStressUpdateTempl<is_ad>::computeStressInitialize(
      91             :     const GenericReal<is_ad> & effective_trial_stress,
      92             :     const GenericRankFourTensor<is_ad> & elasticity_tensor)
      93             : {
      94           0 :   RadialReturnStressUpdateTempl<is_ad>::computeStressInitialize(effective_trial_stress,
      95             :                                                                 elasticity_tensor);
      96     2176512 :   _exp_time = std::pow(_t - _start_time, _m_exponent);
      97     2176512 : }
      98             : 
      99             : template <bool is_ad>
     100             : template <typename ScalarType>
     101             : ScalarType
     102    24006024 : CompositePowerLawCreepStressUpdateTempl<is_ad>::computeResidualInternal(
     103             :     const GenericReal<is_ad> & effective_trial_stress, const ScalarType & scalar)
     104             : {
     105    24006024 :   const ScalarType stress_delta = effective_trial_stress - _three_shear_modulus * scalar;
     106           0 :   ScalarType creep_rate = 0.0;
     107    72018072 :   for (const auto n : make_range(_num_materials))
     108             :   {
     109    48012048 :     creep_rate +=
     110    48012048 :         _coefficient[n] * std::pow(stress_delta, _n_exponent[n]) * (*_switchingFunc[n])[_qp] *
     111    48012048 :         std::exp(-_activation_energy[n] / (_gas_constant * (*_temperature)[_qp])) * _exp_time;
     112             :   }
     113    24006024 :   return creep_rate * _dt - scalar;
     114             : }
     115             : 
     116             : template <bool is_ad>
     117             : GenericReal<is_ad>
     118    24006024 : CompositePowerLawCreepStressUpdateTempl<is_ad>::computeDerivative(
     119             :     const GenericReal<is_ad> & effective_trial_stress, const GenericReal<is_ad> & scalar)
     120             : {
     121    24006024 :   const GenericReal<is_ad> stress_delta = effective_trial_stress - _three_shear_modulus * scalar;
     122           0 :   GenericReal<is_ad> creep_rate_derivative = 0.0;
     123    72018072 :   for (const auto n : make_range(_num_materials))
     124             :   {
     125    48012048 :     creep_rate_derivative +=
     126    48012048 :         -_coefficient[n] * _three_shear_modulus * _n_exponent[n] *
     127    48012048 :         std::pow(stress_delta, _n_exponent[n] - 1.0) * (*_switchingFunc[n])[_qp] *
     128    48012048 :         std::exp(-_activation_energy[n] / (_gas_constant * (*_temperature)[_qp])) * _exp_time;
     129             :   }
     130    24006024 :   return creep_rate_derivative * _dt - 1.0;
     131             : }
     132             : 
     133             : template <bool is_ad>
     134             : Real
     135           0 : CompositePowerLawCreepStressUpdateTempl<is_ad>::computeStrainEnergyRateDensity(
     136             :     const GenericMaterialProperty<RankTwoTensor, is_ad> & stress,
     137             :     const GenericMaterialProperty<RankTwoTensor, is_ad> & strain_rate)
     138             : {
     139           0 :   GenericReal<is_ad> interpolated_exponent = 0.0;
     140           0 :   for (unsigned int n = 0; n < _num_materials; ++n)
     141             :   {
     142           0 :     interpolated_exponent += (_n_exponent[n] / (_n_exponent[n] + 1)) * (*_switchingFunc[n])[_qp];
     143             :   }
     144           0 :   return MetaPhysicL::raw_value(interpolated_exponent *
     145           0 :                                 stress[_qp].doubleContraction((strain_rate)[_qp]));
     146             : }
     147             : 
     148             : template <bool is_ad>
     149             : void
     150     2176512 : CompositePowerLawCreepStressUpdateTempl<is_ad>::computeStressFinalize(
     151             :     const GenericRankTwoTensor<is_ad> & plastic_strain_increment)
     152             : {
     153     2176512 :   _creep_strain[_qp] += plastic_strain_increment;
     154     2176512 : }
     155             : 
     156             : template <bool is_ad>
     157             : void
     158     2176512 : CompositePowerLawCreepStressUpdateTempl<is_ad>::resetIncrementalMaterialProperties()
     159             : {
     160     2176512 :   _creep_strain[_qp] = _creep_strain_old[_qp];
     161     2176512 : }
     162             : 
     163             : template <bool is_ad>
     164             : bool
     165     2176512 : CompositePowerLawCreepStressUpdateTempl<is_ad>::substeppingCapabilityEnabled()
     166             : {
     167     2176512 :   return this->_use_substepping != RadialReturnStressUpdateTempl<is_ad>::SubsteppingType::NONE;
     168             : }
     169             : 
     170             : template class CompositePowerLawCreepStressUpdateTempl<false>;
     171             : template class CompositePowerLawCreepStressUpdateTempl<true>;
     172             : template Real
     173             : CompositePowerLawCreepStressUpdateTempl<false>::computeResidualInternal<Real>(const Real &,
     174             :                                                                               const Real &);
     175             : template ADReal
     176             : CompositePowerLawCreepStressUpdateTempl<true>::computeResidualInternal<ADReal>(const ADReal &,
     177             :                                                                                const ADReal &);
     178             : template ChainedReal
     179             : CompositePowerLawCreepStressUpdateTempl<false>::computeResidualInternal<ChainedReal>(
     180             :     const Real &, const ChainedReal &);
     181             : template ChainedADReal
     182             : CompositePowerLawCreepStressUpdateTempl<true>::computeResidualInternal<ChainedADReal>(
     183             :     const ADReal &, const ChainedADReal &);