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 "MazarsDamage.h"
      16             : 
      17             : #include "RankTwoTensor.h"
      18             : #include "RankFourTensor.h"
      19             : #include "ElasticityTensorTools.h"
      20             : #include "MooseUtils.h"
      21             : 
      22             : #include "libmesh/utility.h"
      23             : 
      24             : registerMooseObject("BlackBearApp", MazarsDamage);
      25             : 
      26             : InputParameters
      27         818 : MazarsDamage::validParams()
      28             : {
      29         818 :   InputParameters params = ScalarDamageBase::validParams();
      30         818 :   params.addClassDescription("Mazars scalar damage model");
      31        1636 :   params.addRequiredCoupledVar("tensile_strength",
      32             :                                "Tensile stress threshold for damage initiation");
      33        1636 :   params.addRequiredParam<Real>("a_t",
      34             :                                 "A_t parameter that controls the shape of the response in tension");
      35        1636 :   params.addRequiredParam<Real>("b_t",
      36             :                                 "B_t parameter that controls the shape of the response in tension");
      37        1636 :   params.addRequiredParam<Real>(
      38             :       "a_c", "A_c parameter that controls the shape of the response in compression");
      39        1636 :   params.addRequiredParam<Real>(
      40             :       "b_c", "B_c parameter that controls the shape of the response in compression");
      41         818 :   return params;
      42           0 : }
      43             : 
      44         630 : MazarsDamage::MazarsDamage(const InputParameters & parameters)
      45             :   : ScalarDamageBase(parameters),
      46             :     GuaranteeConsumer(this),
      47         630 :     _tensile_strength(coupledValue("tensile_strength")),
      48        1260 :     _a_t(getParam<Real>("a_t")),
      49        1260 :     _a_c(getParam<Real>("a_c")),
      50        1260 :     _b_t(getParam<Real>("b_t")),
      51        1260 :     _b_c(getParam<Real>("b_c")),
      52         630 :     _kappa(declareProperty<Real>(_base_name + "kappa")),
      53        1260 :     _kappa_old(getMaterialPropertyOld<Real>(_base_name + "kappa")),
      54        1260 :     _stress(getMaterialProperty<RankTwoTensor>(_base_name + "stress")),
      55        1260 :     _mechanical_strain(getMaterialProperty<RankTwoTensor>(_base_name + "mechanical_strain")),
      56         630 :     _elasticity_tensor_name(_base_name + "elasticity_tensor"),
      57         630 :     _elasticity_tensor(getMaterialPropertyByName<RankFourTensor>(_elasticity_tensor_name)),
      58         630 :     _eigval(3, 0.0),
      59         630 :     _positive_stress(3, 0.0),
      60        1260 :     _positive_strain(3, 0.0)
      61             : {
      62         630 : }
      63             : 
      64             : void
      65        4928 : MazarsDamage::initQpStatefulProperties()
      66             : {
      67        4928 :   ScalarDamageBase::initQpStatefulProperties();
      68        4928 :   _kappa[_qp] = 0.0;
      69        4928 : }
      70             : 
      71             : void
      72         129 : MazarsDamage::initialSetup()
      73             : {
      74         258 :   if (!hasGuaranteedMaterialProperty(_elasticity_tensor_name, Guarantee::ISOTROPIC))
      75           3 :     mooseError("MazarsDamage requires that the elasticity tensor be guaranteed isotropic");
      76         126 : }
      77             : 
      78             : void
      79     1775264 : MazarsDamage::updateQpDamageIndex()
      80             : {
      81     1775264 :   _stress[_qp].symmetricEigenvalues(_eigval);
      82     7101056 :   for (unsigned int i = 0; i < 3; ++i)
      83     7075209 :     _positive_stress[i] = std::max(_eigval[i], 0.0);
      84             : 
      85     1775264 :   _mechanical_strain[_qp].symmetricEigenvalues(_eigval);
      86     7101056 :   for (unsigned int i = 0; i < 3; ++i)
      87     8796823 :     _positive_strain[i] = std::max(_eigval[i], 0.0);
      88             : 
      89             :   const Real equiv_tensile_strain =
      90     1775264 :       std::sqrt(Utility::pow<2>(_positive_strain[0]) + Utility::pow<2>(_positive_strain[1]) +
      91             :                 Utility::pow<2>(_positive_strain[2]));
      92             : 
      93             :   Real alpha_t = 0.0;
      94             :   Real alpha_c = 0.0;
      95             : 
      96     1775264 :   const Real e = ElasticityTensorTools::getIsotropicYoungsModulus(_elasticity_tensor[_qp]);
      97             :   const Real nu = ElasticityTensorTools::getIsotropicPoissonsRatio(_elasticity_tensor[_qp]);
      98             : 
      99     1775264 :   const Real sum_pos_stress = _positive_stress[0] + _positive_stress[1] + _positive_stress[2];
     100             : 
     101     7101056 :   for (unsigned int i = 0; i < 3; ++i)
     102             :   {
     103     5325792 :     const Real epsilon_ti = ((1.0 + nu) * _positive_stress[i] - nu * sum_pos_stress) / e;
     104     5325792 :     const Real epsilon_ci = _eigval[i] - epsilon_ti; //_eigval contains the principal strains
     105     5325792 :     if (!MooseUtils::absoluteFuzzyEqual(equiv_tensile_strain, 0.0))
     106             :     {
     107             :       const Real ets2 = Utility::pow<2>(equiv_tensile_strain);
     108     5294190 :       alpha_t += epsilon_ti * _positive_strain[i] / ets2;
     109     5294190 :       alpha_c += epsilon_ci * _positive_strain[i] / ets2;
     110             :     }
     111             :   }
     112             : 
     113     1775264 :   const Real kappa_0 = _tensile_strength[_qp] / e; // Threshold strain for damage
     114     1775264 :   Real & kappa = _kappa[_qp];
     115     1775264 :   kappa = std::max(kappa_0, _kappa_old[_qp]);
     116             : 
     117     1775264 :   _damage_index[_qp] = 0.0;
     118     1775264 :   if (equiv_tensile_strain > kappa)
     119             :   {
     120      404214 :     kappa = equiv_tensile_strain;
     121             :     const Real tensile_damage =
     122      404214 :         1.0 - kappa_0 * (1.0 - _a_t) / kappa - _a_t * std::exp(-_b_t * (kappa - kappa_0));
     123             :     const Real compressive_damage =
     124      404214 :         1.0 - kappa_0 * (1.0 - _a_c) / kappa - _a_c * std::exp(-_b_c * (kappa - kappa_0));
     125      404214 :     _damage_index[_qp] = std::max(alpha_t * tensile_damage + alpha_c * compressive_damage, 0.0);
     126             :   }
     127             : 
     128             :   // Prevent damage index from decreasing or from being greater than 1
     129     1775264 :   _damage_index[_qp] = std::max(_damage_index[_qp], _damage_index_old[_qp]);
     130     1775264 :   _damage_index[_qp] = std::min(_damage_index[_qp], 1.0);
     131     1775264 : }