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 "CNHIsotropicElasticity.h"
       6             : #include "RaccoonUtils.h"
       7             : 
       8             : registerMooseObject("raccoonApp", CNHIsotropicElasticity);
       9             : 
      10             : InputParameters
      11         142 : CNHIsotropicElasticity::validParams()
      12             : {
      13         142 :   InputParameters params = LargeDeformationElasticityModel::validParams();
      14         142 :   params.addClassDescription("Isotropic Compressible Neo-Hookean hyperelasticity model.");
      15             : 
      16         284 :   params.addRequiredParam<MaterialPropertyName>("bulk_modulus", "The bulk modulus $K$");
      17         284 :   params.addRequiredParam<MaterialPropertyName>("shear_modulus", "The shear modulus $G$");
      18             : 
      19         284 :   params.addRequiredCoupledVar("phase_field", "Name of the phase-field (damage) variable");
      20         284 :   params.addParam<MaterialPropertyName>(
      21             :       "strain_energy_density",
      22             :       "psie",
      23             :       "Name of the strain energy density computed by this material model");
      24         284 :   params.addParam<MaterialPropertyName>("degradation_function", "g", "The degradation function");
      25         284 :   params.addParam<MooseEnum>(
      26         426 :       "decomposition", MooseEnum("NONE SPECTRAL VOLDEV", "NONE"), "The decomposition method");
      27             : 
      28         142 :   return params;
      29           0 : }
      30             : 
      31           3 : CNHIsotropicElasticity::CNHIsotropicElasticity(const InputParameters & parameters)
      32             :   : LargeDeformationElasticityModel(parameters),
      33             :     DerivativeMaterialPropertyNameInterface(),
      34           9 :     _K(getADMaterialPropertyByName<Real>(prependBaseName("bulk_modulus", true))),
      35           6 :     _G(getADMaterialPropertyByName<Real>(prependBaseName("shear_modulus", true))),
      36             : 
      37           6 :     _d_name(getVar("phase_field", 0)->name()),
      38             : 
      39             :     // The strain energy density and its derivatives
      40           3 :     _psie_name(prependBaseName("strain_energy_density", true)),
      41           3 :     _psie(declareADProperty<Real>(_psie_name)),
      42           3 :     _psie_active(declareADProperty<Real>(_psie_name + "_active")),
      43           6 :     _dpsie_dd(declareADProperty<Real>(derivativePropertyName(_psie_name, {_d_name}))),
      44             : 
      45             :     // The degradation function and its derivatives
      46           3 :     _g_name(prependBaseName("degradation_function", true)),
      47           3 :     _g(getADMaterialProperty<Real>(_g_name)),
      48           9 :     _dg_dd(getADMaterialProperty<Real>(derivativePropertyName(_g_name, {_d_name}))),
      49             : 
      50           9 :     _decomposition(getParam<MooseEnum>("decomposition").getEnum<Decomposition>())
      51             : {
      52           9 : }
      53             : 
      54             : ADRankTwoTensor
      55        3144 : CNHIsotropicElasticity::computeMandelStress(const ADRankTwoTensor & Fe,
      56             :                                             const bool plasticity_update)
      57             : {
      58        3144 :   ADRankTwoTensor stress;
      59             : 
      60        3144 :   if (_decomposition == Decomposition::none)
      61        6288 :     stress = computeMandelStressNoDecomposition(Fe, plasticity_update);
      62           0 :   else if (_decomposition == Decomposition::voldev)
      63           0 :     stress = computeMandelStressVolDevDecomposition(Fe, plasticity_update);
      64             :   else
      65           0 :     paramError("decomposition", "Unsupported decomposition type.");
      66             : 
      67        3144 :   return stress;
      68             : }
      69             : 
      70             : ADRankTwoTensor
      71        3144 : CNHIsotropicElasticity::computeMandelStressNoDecomposition(const ADRankTwoTensor & Fe,
      72             :                                                            const bool plasticity_update)
      73             : {
      74             :   // We use the left Cauchy-Green strain
      75        3144 :   ADRankTwoTensor strain;
      76        3144 :   if (plasticity_update)
      77             :   {
      78        1048 :     ADRankTwoTensor expFe = RaccoonUtils::exp(Fe);
      79        1048 :     strain = expFe * expFe.transpose();
      80             :   }
      81             :   else
      82        4192 :     strain = Fe * Fe.transpose();
      83             : 
      84        3144 :   ADReal J = std::sqrt(strain.det());
      85             : 
      86        3144 :   const ADRankTwoTensor I2(ADRankTwoTensor::initIdentity);
      87             :   // Here, we keep the volumetric part no matter what. But ideally, in the case of J2 plasticity,
      88             :   // the volumetric part of the flow should be zero, and we could save some operations.
      89       12576 :   ADRankTwoTensor stress_intact = 0.5 * _K[_qp] * (J * J - 1) * I2 + _G[_qp] * strain.deviatoric();
      90        3144 :   ADRankTwoTensor stress = _g[_qp] * stress_intact;
      91             : 
      92             :   // If plasticity_update == false, then we are not in the middle of a plasticity update, hence we
      93             :   // compute the strain energy density
      94        3144 :   if (!plasticity_update)
      95             :   {
      96        2096 :     ADRankTwoTensor strain_bar = std::pow(J, -2. / 3.) * strain;
      97        6288 :     ADReal U = 0.5 * _K[_qp] * (0.5 * (J * J - 1) - std::log(J));
      98        6288 :     ADReal W = 0.5 * _G[_qp] * (strain_bar.trace() - 3.0);
      99        2096 :     _psie_active[_qp] = U + W;
     100        4192 :     _psie[_qp] = _g[_qp] * _psie_active[_qp];
     101        4192 :     _dpsie_dd[_qp] = _dg_dd[_qp] * _psie_active[_qp];
     102             :   }
     103             : 
     104        3144 :   return stress;
     105             : }
     106             : 
     107             : ADRankTwoTensor
     108           0 : CNHIsotropicElasticity::computeMandelStressVolDevDecomposition(const ADRankTwoTensor & Fe,
     109             :                                                                const bool plasticity_update)
     110             : {
     111             :   // We use the left Cauchy-Green strain
     112           0 :   ADRankTwoTensor strain;
     113           0 :   if (plasticity_update)
     114             :   {
     115           0 :     ADRankTwoTensor expFe = RaccoonUtils::exp(Fe);
     116           0 :     strain = expFe * expFe.transpose();
     117             :   }
     118             :   else
     119           0 :     strain = Fe * Fe.transpose();
     120             : 
     121           0 :   ADReal J = std::sqrt(strain.det());
     122             : 
     123           0 :   const ADRankTwoTensor I2(ADRankTwoTensor::initIdentity);
     124             :   // Here, we keep the volumetric part no matter what. But ideally, in the case of J2 plasticity,
     125             :   // the volumetric part of the flow should be zero, and we could save some operations.
     126           0 :   ADRankTwoTensor stress_vol = 0.5 * _K[_qp] * (J * J - 1) * I2;
     127           0 :   ADRankTwoTensor stress_dev = _G[_qp] * strain.deviatoric();
     128             :   ADRankTwoTensor stress =
     129           0 :       J > 1 ? _g[_qp] * (stress_vol + stress_dev) : stress_vol + _g[_qp] * stress_dev;
     130             : 
     131             :   // If plasticity_update == false, then we are not in the middle of a plasticity update, hence we
     132             :   // compute the strain energy density
     133           0 :   if (!plasticity_update)
     134             :   {
     135           0 :     ADRankTwoTensor strain_bar = std::pow(J, -2. / 3.) * strain;
     136           0 :     ADReal U = 0.5 * _K[_qp] * (0.5 * (J * J - 1) - std::log(J));
     137           0 :     ADReal W = 0.5 * _G[_qp] * (strain_bar.trace() - 3.0);
     138           0 :     _psie_active[_qp] = J > 1 ? U + W : W;
     139           0 :     _psie[_qp] = _g[_qp] * _psie_active[_qp];
     140           0 :     _dpsie_dd[_qp] = _dg_dd[_qp] * _psie_active[_qp];
     141             :   }
     142             : 
     143           0 :   return stress;
     144             : }