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         140 : CNHIsotropicElasticity::validParams()
      12             : {
      13         140 :   InputParameters params = LargeDeformationElasticityModel::validParams();
      14         140 :   params.addClassDescription("Isotropic Compressible Neo-Hookean hyperelasticity model.");
      15             : 
      16         280 :   params.addRequiredParam<MaterialPropertyName>("bulk_modulus", "The bulk modulus $K$");
      17         280 :   params.addRequiredParam<MaterialPropertyName>("shear_modulus", "The shear modulus $G$");
      18             : 
      19         280 :   params.addRequiredCoupledVar("phase_field", "Name of the phase-field (damage) variable");
      20         280 :   params.addParam<MaterialPropertyName>(
      21             :       "strain_energy_density",
      22             :       "psie",
      23             :       "Name of the strain energy density computed by this material model");
      24         280 :   params.addParam<MaterialPropertyName>("degradation_function", "g", "The degradation function");
      25         280 :   params.addParam<MooseEnum>(
      26         420 :       "decomposition", MooseEnum("NONE SPECTRAL VOLDEV", "NONE"), "The decomposition method");
      27             : 
      28         140 :   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             :   using std::sqrt;
      75             :   // We use the left Cauchy-Green strain
      76        3144 :   ADRankTwoTensor strain;
      77        3144 :   if (plasticity_update)
      78             :   {
      79        1048 :     ADRankTwoTensor expFe = RaccoonUtils::exp(Fe);
      80        2096 :     strain = expFe * expFe.transpose();
      81             :   }
      82             :   else
      83        4192 :     strain = Fe * Fe.transpose();
      84             : 
      85        3144 :   ADReal J = sqrt(strain.det());
      86             : 
      87        3144 :   const ADRankTwoTensor I2(ADRankTwoTensor::initIdentity);
      88             :   // Here, we keep the volumetric part no matter what. But ideally, in the case of J2 plasticity,
      89             :   // the volumetric part of the flow should be zero, and we could save some operations.
      90       12576 :   ADRankTwoTensor stress_intact = 0.5 * _K[_qp] * (J * J - 1) * I2 + _G[_qp] * strain.deviatoric();
      91        3144 :   ADRankTwoTensor stress = _g[_qp] * stress_intact;
      92             : 
      93             :   // If plasticity_update == false, then we are not in the middle of a plasticity update, hence we
      94             :   // compute the strain energy density
      95        3144 :   if (!plasticity_update)
      96             :   {
      97             :     using std::log;
      98             :     using std::pow;
      99        2096 :     ADRankTwoTensor strain_bar = pow(J, -2. / 3.) * strain;
     100        6288 :     ADReal U = 0.5 * _K[_qp] * (0.5 * (J * J - 1) - log(J));
     101        6288 :     ADReal W = 0.5 * _G[_qp] * (strain_bar.trace() - 3.0);
     102        2096 :     _psie_active[_qp] = U + W;
     103        4192 :     _psie[_qp] = _g[_qp] * _psie_active[_qp];
     104        4192 :     _dpsie_dd[_qp] = _dg_dd[_qp] * _psie_active[_qp];
     105             :   }
     106             : 
     107        3144 :   return stress;
     108             : }
     109             : 
     110             : ADRankTwoTensor
     111           0 : CNHIsotropicElasticity::computeMandelStressVolDevDecomposition(const ADRankTwoTensor & Fe,
     112             :                                                                const bool plasticity_update)
     113             : {
     114             :   using std::sqrt;
     115             :   // We use the left Cauchy-Green strain
     116           0 :   ADRankTwoTensor strain;
     117           0 :   if (plasticity_update)
     118             :   {
     119           0 :     ADRankTwoTensor expFe = RaccoonUtils::exp(Fe);
     120           0 :     strain = expFe * expFe.transpose();
     121             :   }
     122             :   else
     123           0 :     strain = Fe * Fe.transpose();
     124             : 
     125           0 :   ADReal J = sqrt(strain.det());
     126             : 
     127           0 :   const ADRankTwoTensor I2(ADRankTwoTensor::initIdentity);
     128             :   // Here, we keep the volumetric part no matter what. But ideally, in the case of J2 plasticity,
     129             :   // the volumetric part of the flow should be zero, and we could save some operations.
     130           0 :   ADRankTwoTensor stress_vol = 0.5 * _K[_qp] * (J * J - 1) * I2;
     131           0 :   ADRankTwoTensor stress_dev = _G[_qp] * strain.deviatoric();
     132             :   ADRankTwoTensor stress =
     133           0 :       J > 1 ? _g[_qp] * (stress_vol + stress_dev) : stress_vol + _g[_qp] * stress_dev;
     134             : 
     135             :   // If plasticity_update == false, then we are not in the middle of a plasticity update, hence we
     136             :   // compute the strain energy density
     137           0 :   if (!plasticity_update)
     138             :   {
     139             :     using std::log;
     140             :     using std::pow;
     141           0 :     ADRankTwoTensor strain_bar = pow(J, -2. / 3.) * strain;
     142           0 :     ADReal U = 0.5 * _K[_qp] * (0.5 * (J * J - 1) - log(J));
     143           0 :     ADReal W = 0.5 * _G[_qp] * (strain_bar.trace() - 3.0);
     144           0 :     _psie_active[_qp] = J > 1 ? U + W : W;
     145           0 :     _psie[_qp] = _g[_qp] * _psie_active[_qp];
     146           0 :     _dpsie_dd[_qp] = _dg_dd[_qp] * _psie_active[_qp];
     147             :   }
     148             : 
     149           0 :   return stress;
     150             : }