Line data    Source code

       1             : //* This file is part of the RACCOON application
       2             : //* being developed at Dolbow lab at Duke University
       3             : //* http://dolbow.pratt.duke.edu
       4             : 
       5             : #include "SmallDeformationIsotropicElasticity.h"
       6             : #include "RaccoonUtils.h"
       7             : 
       8             : registerMooseObject("raccoonApp", SmallDeformationIsotropicElasticity);
       9             : 
      10             : InputParameters
      11         214 : SmallDeformationIsotropicElasticity::validParams()
      12             : {
      13         214 :   InputParameters params = SmallDeformationElasticityModel::validParams();
      14         214 :   params.addClassDescription("Isotropic elasticity under small strain asumptions.");
      15             : 
      16         428 :   params.addRequiredParam<MaterialPropertyName>("bulk_modulus", "The bulk modulus $K$");
      17         428 :   params.addRequiredParam<MaterialPropertyName>("shear_modulus", "The shear modulus $G$");
      18             : 
      19         428 :   params.addRequiredCoupledVar("phase_field", "Name of the phase-field (damage) variable");
      20         428 :   params.addParam<MaterialPropertyName>(
      21             :       "strain_energy_density",
      22             :       "psie",
      23             :       "Name of the strain energy density computed by this material model");
      24         428 :   params.addParam<MaterialPropertyName>("degradation_function", "g", "The degradation function");
      25         428 :   params.addParam<MooseEnum>(
      26         642 :       "decomposition", MooseEnum("NONE SPECTRAL VOLDEV", "NONE"), "The decomposition method");
      27             : 
      28         214 :   return params;
      29           0 : }
      30             : 
      31          57 : SmallDeformationIsotropicElasticity::SmallDeformationIsotropicElasticity(
      32          57 :     const InputParameters & parameters)
      33             :   : SmallDeformationElasticityModel(parameters),
      34             :     DerivativeMaterialPropertyNameInterface(),
      35         171 :     _K(getADMaterialPropertyByName<Real>(prependBaseName("bulk_modulus", true))),
      36         114 :     _G(getADMaterialPropertyByName<Real>(prependBaseName("shear_modulus", true))),
      37             : 
      38         114 :     _d_name(getVar("phase_field", 0)->name()),
      39             : 
      40             :     // The strain energy density and its derivatives
      41          57 :     _psie_name(prependBaseName("strain_energy_density", true)),
      42          57 :     _psie(declareADProperty<Real>(_psie_name)),
      43          57 :     _psie_active(declareADProperty<Real>(_psie_name + "_active")),
      44         114 :     _dpsie_dd(declareADProperty<Real>(derivativePropertyName(_psie_name, {_d_name}))),
      45             : 
      46             :     // The degradation function and its derivatives
      47          57 :     _g_name(prependBaseName("degradation_function", true)),
      48          57 :     _g(getADMaterialProperty<Real>(_g_name)),
      49         171 :     _dg_dd(getADMaterialProperty<Real>(derivativePropertyName(_g_name, {_d_name}))),
      50             : 
      51         171 :     _decomposition(getParam<MooseEnum>("decomposition").getEnum<Decomposition>())
      52             : {
      53         171 : }
      54             : 
      55             : ADRankTwoTensor
      56    12320102 : SmallDeformationIsotropicElasticity::computeStress(const ADRankTwoTensor & strain)
      57             : {
      58    12320102 :   ADRankTwoTensor stress;
      59             : 
      60    12320102 :   if (_decomposition == Decomposition::none)
      61    20079156 :     stress = computeStressNoDecomposition(strain);
      62     2280524 :   else if (_decomposition == Decomposition::spectral)
      63     3657320 :     stress = computeStressSpectralDecomposition(strain);
      64      451864 :   else if (_decomposition == Decomposition::voldev)
      65      903728 :     stress = computeStressVolDevDecomposition(strain);
      66             :   else
      67           0 :     paramError("decomposition", "Unsupported decomposition type.");
      68             : 
      69    12320102 :   return stress;
      70             : }
      71             : 
      72             : ADRankTwoTensor
      73    10039578 : SmallDeformationIsotropicElasticity::computeStressNoDecomposition(const ADRankTwoTensor & strain)
      74             : {
      75    10039578 :   const ADRankTwoTensor I2(ADRankTwoTensor::initIdentity);
      76    40158312 :   ADRankTwoTensor stress_intact = _K[_qp] * strain.trace() * I2 + 2 * _G[_qp] * strain.deviatoric();
      77    10039578 :   ADRankTwoTensor stress = _g[_qp] * stress_intact;
      78             : 
      79    20079156 :   _psie_active[_qp] = 0.5 * stress_intact.doubleContraction(strain);
      80    20079156 :   _psie[_qp] = _g[_qp] * _psie_active[_qp];
      81    20079156 :   _dpsie_dd[_qp] = _dg_dd[_qp] * _psie_active[_qp];
      82             : 
      83    10039578 :   return stress;
      84             : }
      85             : 
      86             : ADRankTwoTensor
      87     1828660 : SmallDeformationIsotropicElasticity::computeStressSpectralDecomposition(
      88             :     const ADRankTwoTensor & strain)
      89             : {
      90     3657320 :   const ADReal lambda = _K[_qp] - 2 * _G[_qp] / LIBMESH_DIM;
      91     1828660 :   const ADRankTwoTensor I2(ADRankTwoTensor::initIdentity);
      92     1828660 :   ADReal strain_tr = strain.trace();
      93     1828660 :   ADReal strain_tr_pos = RaccoonUtils::Macaulay(strain_tr);
      94             : 
      95             :   // Spectral decomposition
      96     1828660 :   ADRankTwoTensor strain_pos = RaccoonUtils::spectralDecomposition(strain);
      97             : 
      98             :   // Stress
      99     7314640 :   ADRankTwoTensor stress_intact = _K[_qp] * strain.trace() * I2 + 2 * _G[_qp] * strain.deviatoric();
     100     3657320 :   ADRankTwoTensor stress_pos = lambda * strain_tr_pos * I2 + 2 * _G[_qp] * strain_pos;
     101     1828660 :   ADRankTwoTensor stress_neg = stress_intact - stress_pos;
     102     1828660 :   ADRankTwoTensor stress = _g[_qp] * stress_pos + stress_neg;
     103             : 
     104             :   // Strain energy density
     105             :   ADReal psie_intact =
     106     3657320 :       0.5 * lambda * strain_tr * strain_tr + _G[_qp] * strain.doubleContraction(strain);
     107     3657320 :   _psie_active[_qp] = 0.5 * lambda * strain_tr_pos * strain_tr_pos +
     108     3657320 :                       _G[_qp] * strain_pos.doubleContraction(strain_pos);
     109             :   ADReal psie_inactive = psie_intact - _psie_active[_qp];
     110     3657320 :   _psie[_qp] = _g[_qp] * _psie_active[_qp] + psie_inactive;
     111     3657320 :   _dpsie_dd[_qp] = _dg_dd[_qp] * _psie_active[_qp];
     112             : 
     113     1828660 :   return stress;
     114             : }
     115             : 
     116             : ADRankTwoTensor
     117      451864 : SmallDeformationIsotropicElasticity::computeStressVolDevDecomposition(
     118             :     const ADRankTwoTensor & strain)
     119             : {
     120      451864 :   const ADRankTwoTensor I2(ADRankTwoTensor::initIdentity);
     121             : 
     122             :   // Volumetric-deviatoric decomposition
     123      451864 :   ADReal strain_tr = strain.trace();
     124      451864 :   ADReal strain_tr_pos = RaccoonUtils::Macaulay(strain_tr);
     125             :   ADReal strain_tr_neg = strain_tr - strain_tr_pos;
     126      451864 :   ADRankTwoTensor strain_dev = strain.deviatoric();
     127             : 
     128             :   // Stress
     129     1807456 :   ADRankTwoTensor stress_intact = _K[_qp] * strain.trace() * I2 + 2 * _G[_qp] * strain.deviatoric();
     130      451864 :   ADRankTwoTensor stress_neg = _K[_qp] * strain_tr_neg * I2;
     131      451864 :   ADRankTwoTensor stress_pos = stress_intact - stress_neg;
     132      451864 :   ADRankTwoTensor stress = _g[_qp] * stress_pos + stress_neg;
     133             : 
     134             :   // Strain energy density
     135             :   ADReal psie_intact =
     136      903728 :       0.5 * _K[_qp] * strain_tr * strain_tr + _G[_qp] * strain_dev.doubleContraction(strain_dev);
     137      903728 :   ADReal psie_inactive = 0.5 * _K[_qp] * strain_tr_neg * strain_tr_neg;
     138      451864 :   _psie_active[_qp] = psie_intact - psie_inactive;
     139      903728 :   _psie[_qp] = _g[_qp] * _psie_active[_qp] + psie_inactive;
     140      903728 :   _dpsie_dd[_qp] = _dg_dd[_qp] * _psie_active[_qp];
     141             : 
     142      451864 :   return stress;
     143             : }