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.setDocUnit("epsilon_f", "unitless");
      39         330 :   params.addRequiredParam<Real>("creep_law_exponent", "Exponent of the creep power law.");
      40         330 :   params.setDocUnit("creep_law_exponent", "unitless");
      41         330 :   params.addParam<Real>(
      42         330 :       "reduction_factor", 1000.0, "Reduction on the load-carrying capacity (stress).");
      43         330 :   params.setDocUnit("reduction_factor", "unitless");
      44         495 :   params.addRangeCheckedParam<Real>(
      45             :       "reduction_damage_threshold",
      46         330 :       0.9,
      47             :       "reduction_damage_threshold <= 1.0 & reduction_damage_threshold >= 0.0",
      48             :       "Starting value of damage that will trigger linear reduction of "
      49             :       "quadrature point's load-carrying capacity.");
      50         330 :   params.setDocUnit("reduction_damage_threshold", "unitless");
      51         330 :   params.addRequiredParam<std::vector<std::string>>(
      52             :       "creep_strain_names", "Names of the creep strains driving the damage.");
      53         165 :   return params;
      54           0 : }
      55             : 
      56             : template <bool is_ad>
      57         126 : SteelCreepDamageOhTempl<is_ad>::SteelCreepDamageOhTempl(const InputParameters & parameters)
      58             :   : ScalarDamageBaseTempl<is_ad>(parameters),
      59             :     GuaranteeConsumer(this),
      60         126 :     _creep_strain_names(this->template getParam<std::vector<std::string>>("creep_strain_names")),
      61         126 :     _creep_model(_creep_strain_names.size()),
      62         126 :     _combined_creep_strain(this->template declareGenericProperty<RankTwoTensor, is_ad>(
      63         126 :         _base_name + "combined_creep_strain")),
      64         126 :     _combined_creep_strain_old(
      65         126 :         this->template getMaterialPropertyOld<RankTwoTensor>(_base_name + "combined_creep_strain")),
      66         252 :     _epsilon_f(this->template getParam<Real>("epsilon_f")),
      67         252 :     _reduction_factor(this->template getParam<Real>("reduction_factor")),
      68         252 :     _reduction_damage_threshold(this->template getParam<Real>("reduction_damage_threshold")),
      69         252 :     _creep_law_exponent(this->template getParam<Real>("creep_law_exponent")),
      70         126 :     _stress(this->template getGenericMaterialProperty<RankTwoTensor, is_ad>(_base_name + "stress")),
      71         252 :     _omega(this->template declareGenericProperty<Real, is_ad>(_base_name + "omega")),
      72         378 :     _omega_old(this->template getMaterialPropertyOld<Real>(_base_name + "omega"))
      73             : {
      74         252 :   for (unsigned int i = 0; i < _creep_strain_names.size(); ++i)
      75             :   {
      76         126 :     _creep_model[i] = &this->template getGenericMaterialProperty<RankTwoTensor, is_ad>(
      77             :         _base_name + _creep_strain_names[i]);
      78             :   }
      79         126 :   if (MooseUtils::absoluteFuzzyEqual(_creep_law_exponent, -0.5, TOLERANCE))
      80           0 :     this->template paramError(
      81             :         "creep_law_exponent",
      82             :         "creep_law_exponent cannot be -0.5 due to singularities in the multiaxial update of "
      83             :         "the uniaxial ductility value.");
      84         126 : }
      85             : 
      86             : template <bool is_ad>
      87             : void
      88         840 : SteelCreepDamageOhTempl<is_ad>::initQpStatefulProperties()
      89             : {
      90         840 :   ScalarDamageBaseTempl<is_ad>::initQpStatefulProperties();
      91         840 :   _omega[_qp] = 0.0;
      92         840 :   _combined_creep_strain[_qp].zero();
      93         840 : }
      94             : 
      95             : template <bool is_ad>
      96             : void
      97      234696 : SteelCreepDamageOhTempl<is_ad>::updateQpDamageIndex()
      98             : {
      99             :   Real epsilon_f_star;
     100             : 
     101      234696 :   const auto & stress = MetaPhysicL::raw_value(_stress[_qp]);
     102      234696 :   const auto & vonMises = RankTwoScalarTools::vonMisesStress(stress);
     103             : 
     104      234696 :   if (vonMises > TOLERANCE)
     105             :   {
     106      233422 :     const auto & h = RankTwoScalarTools::triaxialityStress(stress);
     107             : 
     108             :     // Let's only modify axial ductility for significant values of stress triaxiality
     109      233422 :     if (h > TOLERANCE)
     110             :       // Update creep ductility due to triaxiliaty effects
     111      200256 :       epsilon_f_star =
     112      200256 :           _epsilon_f *
     113      200256 :           std::sinh(2.0 / 3.0 * (_creep_law_exponent - 0.5) / (_creep_law_exponent + 0.5)) /
     114      200256 :           std::sinh(2.0 * (_creep_law_exponent - 0.5) / (_creep_law_exponent + 0.5) * h);
     115             :     else
     116       33166 :       epsilon_f_star = _epsilon_f;
     117             :   }
     118             :   else
     119        1274 :     epsilon_f_star = _epsilon_f;
     120             : 
     121             :   // If the value obtained from multiaxiality equation isn't reasonable, do not update damage
     122             :   // This would also cause an FPE a few lines below
     123      234696 :   if (std::abs(epsilon_f_star) < TOLERANCE)
     124             :   {
     125       40467 :     _damage_index[_qp] = _damage_index_old[_qp];
     126       40467 :     _omega[_qp] = _omega_old[_qp];
     127       40467 :     return;
     128             :   }
     129             : 
     130      194229 :   _combined_creep_strain[_qp].zero();
     131      388458 :   for (unsigned int i = 0; i < _creep_strain_names.size(); ++i)
     132      194229 :     _combined_creep_strain[_qp] += (*_creep_model[i])[_qp];
     133             : 
     134      194229 :   GenericRankTwoTensor<is_ad> creep_increment =
     135      194229 :       _combined_creep_strain[_qp] - _combined_creep_strain_old[_qp];
     136             : 
     137             :   // Avoid derivative's divide by zero error
     138             :   Real epsilon_ad = 1.0e-14;
     139             : 
     140             :   // Get equivalent creep strain increment
     141      194229 :   GenericReal<is_ad> equivalent_creep_increment =
     142      129360 :       1 / std::sqrt(2) *
     143      129549 :       std::sqrt(Utility::pow<2>(creep_increment(0, 0) - creep_increment(1, 1)) +
     144      129549 :                 Utility::pow<2>(creep_increment(1, 1) - creep_increment(2, 2)) +
     145      129549 :                 Utility::pow<2>(creep_increment(2, 2) - creep_increment(0, 0)) +
     146      194229 :                 1.5 * creep_increment(0, 1) * creep_increment(0, 1) +
     147      194229 :                 1.5 * creep_increment(1, 2) * creep_increment(1, 2) +
     148      258909 :                 1.5 * creep_increment(2, 0) * creep_increment(2, 0) + epsilon_ad);
     149             : 
     150      258909 :   _omega[_qp] = _omega_old[_qp] + equivalent_creep_increment / epsilon_f_star;
     151             : 
     152      194229 :   if (_omega[_qp] > 1.0)
     153      141729 :     _omega[_qp] = 1.0;
     154             : 
     155             :   Real threshold = MetaPhysicL::raw_value(_omega[_qp]);
     156             : 
     157      194229 :   if (threshold < _reduction_damage_threshold)
     158             :   {
     159             :     // If threshold is not reached, there is no damage.
     160       52500 :     _damage_index[_qp] = 0.0;
     161             :     return;
     162             :   }
     163             : 
     164      141729 :   if (threshold > 1.0)
     165             :     threshold = 1;
     166             : 
     167             :   // Cast damage index from the omega model to the traditional damage index.
     168      141729 :   Real factor = (threshold - _reduction_damage_threshold) / (1 - _reduction_damage_threshold) *
     169      141729 :                 _reduction_factor;
     170      141729 :   _damage_index[_qp] = 1.0 - 1.0 / factor;
     171             : 
     172             :   // Account for a corner case where threshold == _reduction_damage_threshold
     173      141729 :   if (_damage_index[_qp] < 0.0)
     174           0 :     _damage_index[_qp] = 0.0;
     175             : }
     176             : 
     177             : template class SteelCreepDamageOhTempl<false>;
     178             : template class SteelCreepDamageOhTempl<true>;