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