Line data    Source code

       1             : /****************************************************************/
       2             : /*               DO NOT MODIFY THIS HEADER                      */
       3             : /*                       BlackBear                              */
       4             : /*                                                              */
       5             : /*           (c) 2017 Battelle Energy Alliance, LLC             */
       6             : /*                   ALL RIGHTS RESERVED                        */
       7             : /*                                                              */
       8             : /*          Prepared by Battelle Energy Alliance, LLC           */
       9             : /*            Under Contract No. DE-AC07-05ID14517              */
      10             : /*            With the U. S. Department of Energy               */
      11             : /*                                                              */
      12             : /*            See COPYRIGHT for full restrictions               */
      13             : /****************************************************************/
      14             : 
      15             : #include "SteelCreepDamageOh.h"
      16             : 
      17             : #include "RankTwoTensor.h"
      18             : #include "RankFourTensor.h"
      19             : #include "ElasticityTensorTools.h"
      20             : #include "MooseUtils.h"
      21             : #include <cmath>
      22             : #include "RankTwoScalarTools.h"
      23             : 
      24             : #include "libmesh/utility.h"
      25             : 
      26             : registerMooseObject("BlackBearApp", SteelCreepDamageOh);
      27             : registerMooseObject("BlackBearApp", ADSteelCreepDamageOh);
      28             : 
      29             : template <bool is_ad>
      30             : InputParameters
      31         165 : SteelCreepDamageOhTempl<is_ad>::validParams()
      32             : {
      33         165 :   InputParameters params = ScalarDamageBaseTempl<is_ad>::validParams();
      34         165 :   params.addClassDescription(
      35             :       "Steel scalar damage model: Material 'suddenly' loses load-carrying "
      36             :       "capacity. This can model, e.g., 316H creep failure. See Oh et al (2011).");
      37         330 :   params.addRequiredParam<Real>("epsilon_f", "Uniaxial creep fracture strain (creep ductility).");
      38         330 :   params.addRequiredParam<Real>("creep_law_exponent", "Exponent of the creep power law.");
      39         330 :   params.addParam<Real>(
      40         330 :       "reduction_factor", 1000.0, "Reduction on the load-carrying capacity (stress).");
      41         495 :   params.addRangeCheckedParam<Real>(
      42             :       "reduction_damage_threshold",
      43         330 :       0.9,
      44             :       "reduction_damage_threshold <= 1.0 & reduction_damage_threshold >= 0.0",
      45             :       "Starting value of damage that will trigger linear reduction of "
      46             :       "quadrature point's load-carrying capacity.");
      47         330 :   params.addRequiredParam<std::vector<std::string>>(
      48             :       "creep_strain_names", "Names of the creep strains driving the damage.");
      49         165 :   return params;
      50           0 : }
      51             : 
      52             : template <bool is_ad>
      53         126 : SteelCreepDamageOhTempl<is_ad>::SteelCreepDamageOhTempl(const InputParameters & parameters)
      54             :   : ScalarDamageBaseTempl<is_ad>(parameters),
      55             :     GuaranteeConsumer(this),
      56         126 :     _creep_strain_names(this->template getParam<std::vector<std::string>>("creep_strain_names")),
      57         126 :     _creep_model(_creep_strain_names.size()),
      58         126 :     _combined_creep_strain(this->template declareGenericProperty<RankTwoTensor, is_ad>(
      59         126 :         _base_name + "combined_creep_strain")),
      60         126 :     _combined_creep_strain_old(
      61         126 :         this->template getMaterialPropertyOld<RankTwoTensor>(_base_name + "combined_creep_strain")),
      62         252 :     _epsilon_f(this->template getParam<Real>("epsilon_f")),
      63         252 :     _reduction_factor(this->template getParam<Real>("reduction_factor")),
      64         252 :     _reduction_damage_threshold(this->template getParam<Real>("reduction_damage_threshold")),
      65         252 :     _creep_law_exponent(this->template getParam<Real>("creep_law_exponent")),
      66         126 :     _stress(this->template getGenericMaterialProperty<RankTwoTensor, is_ad>(_base_name + "stress")),
      67         252 :     _omega(this->template declareGenericProperty<Real, is_ad>(_base_name + "omega")),
      68         378 :     _omega_old(this->template getMaterialPropertyOld<Real>(_base_name + "omega"))
      69             : {
      70         252 :   for (unsigned int i = 0; i < _creep_strain_names.size(); ++i)
      71             :   {
      72         126 :     _creep_model[i] = &this->template getGenericMaterialProperty<RankTwoTensor, is_ad>(
      73             :         _base_name + _creep_strain_names[i]);
      74             :   }
      75         126 :   if (MooseUtils::absoluteFuzzyEqual(_creep_law_exponent, -0.5, TOLERANCE))
      76           0 :     this->template paramError(
      77             :         "creep_law_exponent",
      78             :         "creep_law_exponent cannot be -0.5 due to singularities in the multiaxial update of "
      79             :         "the uniaxial ductility value.");
      80         126 : }
      81             : 
      82             : template <bool is_ad>
      83             : void
      84         840 : SteelCreepDamageOhTempl<is_ad>::initQpStatefulProperties()
      85             : {
      86         840 :   ScalarDamageBaseTempl<is_ad>::initQpStatefulProperties();
      87         840 :   _omega[_qp] = 0.0;
      88         840 :   _combined_creep_strain[_qp].zero();
      89         840 : }
      90             : 
      91             : template <bool is_ad>
      92             : void
      93      234696 : SteelCreepDamageOhTempl<is_ad>::updateQpDamageIndex()
      94             : {
      95             :   Real epsilon_f_star;
      96             : 
      97      234696 :   const auto & stress = MetaPhysicL::raw_value(_stress[_qp]);
      98      234696 :   const auto & vonMises = RankTwoScalarTools::vonMisesStress(stress);
      99             : 
     100      234696 :   if (vonMises > TOLERANCE)
     101             :   {
     102      233422 :     const auto & h = RankTwoScalarTools::triaxialityStress(stress);
     103             : 
     104             :     // Let's only modify axial ductility for significant values of stress triaxiality
     105      233422 :     if (h > TOLERANCE)
     106             :       // Update creep ductility due to triaxiliaty effects
     107      200256 :       epsilon_f_star =
     108      200256 :           _epsilon_f *
     109      200256 :           std::sinh(2.0 / 3.0 * (_creep_law_exponent - 0.5) / (_creep_law_exponent + 0.5)) /
     110      200256 :           std::sinh(2.0 * (_creep_law_exponent - 0.5) / (_creep_law_exponent + 0.5) * h);
     111             :     else
     112       33166 :       epsilon_f_star = _epsilon_f;
     113             :   }
     114             :   else
     115        1274 :     epsilon_f_star = _epsilon_f;
     116             : 
     117             :   // If the value obtained from multiaxiality equation isn't reasonable, do not update damage
     118             :   // This would also cause an FPE a few lines below
     119      234696 :   if (std::abs(epsilon_f_star) < TOLERANCE)
     120             :   {
     121       40467 :     _damage_index[_qp] = _damage_index_old[_qp];
     122       40467 :     _omega[_qp] = _omega_old[_qp];
     123       40467 :     return;
     124             :   }
     125             : 
     126      194229 :   _combined_creep_strain[_qp].zero();
     127      388458 :   for (unsigned int i = 0; i < _creep_strain_names.size(); ++i)
     128      194229 :     _combined_creep_strain[_qp] += (*_creep_model[i])[_qp];
     129             : 
     130      194229 :   GenericRankTwoTensor<is_ad> creep_increment =
     131      194229 :       _combined_creep_strain[_qp] - _combined_creep_strain_old[_qp];
     132             : 
     133             :   // Avoid derivative's divide by zero error
     134             :   Real epsilon_ad = 1.0e-14;
     135             : 
     136             :   // Get equivalent creep strain increment
     137      194229 :   GenericReal<is_ad> equivalent_creep_increment =
     138      129360 :       1 / std::sqrt(2) *
     139      129549 :       std::sqrt(Utility::pow<2>(creep_increment(0, 0) - creep_increment(1, 1)) +
     140      129549 :                 Utility::pow<2>(creep_increment(1, 1) - creep_increment(2, 2)) +
     141      129549 :                 Utility::pow<2>(creep_increment(2, 2) - creep_increment(0, 0)) +
     142      194229 :                 1.5 * creep_increment(0, 1) * creep_increment(0, 1) +
     143      194229 :                 1.5 * creep_increment(1, 2) * creep_increment(1, 2) +
     144      258909 :                 1.5 * creep_increment(2, 0) * creep_increment(2, 0) + epsilon_ad);
     145             : 
     146      258909 :   _omega[_qp] = _omega_old[_qp] + equivalent_creep_increment / epsilon_f_star;
     147             : 
     148      194229 :   if (_omega[_qp] > 1.0)
     149      141729 :     _omega[_qp] = 1.0;
     150             : 
     151             :   Real threshold = MetaPhysicL::raw_value(_omega[_qp]);
     152             : 
     153      194229 :   if (threshold < _reduction_damage_threshold)
     154             :   {
     155             :     // If threshold is not reached, there is no damage.
     156       52500 :     _damage_index[_qp] = 0.0;
     157             :     return;
     158             :   }
     159             : 
     160      141729 :   if (threshold > 1.0)
     161             :     threshold = 1;
     162             : 
     163             :   // Cast damage index from the omega model to the traditional damage index.
     164      141729 :   Real factor = (threshold - _reduction_damage_threshold) / (1 - _reduction_damage_threshold) *
     165      141729 :                 _reduction_factor;
     166      141729 :   _damage_index[_qp] = 1.0 - 1.0 / factor;
     167             : 
     168             :   // Account for a corner case where threshold == _reduction_damage_threshold
     169      141729 :   if (_damage_index[_qp] < 0.0)
     170           0 :     _damage_index[_qp] = 0.0;
     171             : }
     172             : 
     173             : template class SteelCreepDamageOhTempl<false>;
     174             : template class SteelCreepDamageOhTempl<true>;