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 "HenckyIsotropicElasticity.h"
       6             : #include "RaccoonUtils.h"
       7             : 
       8             : registerMooseObject("raccoonApp", HenckyIsotropicElasticity);
       9             : 
      10             : InputParameters
      11         198 : HenckyIsotropicElasticity::validParams()
      12             : {
      13         198 :   InputParameters params = LargeDeformationElasticityModel::validParams();
      14         198 :   params.addClassDescription("Isotropic Hencky-type hyperelasticity model. The logarithmic right "
      15             :                              "Cauchy-Green strain tensor is used as the strain measure.");
      16             : 
      17         396 :   params.addRequiredParam<MaterialPropertyName>("bulk_modulus", "The bulk modulus $K$");
      18         396 :   params.addRequiredParam<MaterialPropertyName>("shear_modulus", "The shear modulus $G$");
      19             : 
      20         396 :   params.addRequiredCoupledVar("phase_field", "Name of the phase-field (damage) variable");
      21         396 :   params.addParam<MaterialPropertyName>(
      22             :       "strain_energy_density",
      23             :       "psie",
      24             :       "Name of the strain energy density computed by this material model");
      25         396 :   params.addParam<MaterialPropertyName>("degradation_function", "g", "The degradation function");
      26         396 :   params.addParam<MooseEnum>(
      27         594 :       "decomposition", MooseEnum("NONE SPECTRAL VOLDEV", "NONE"), "The decomposition method");
      28             : 
      29         198 :   return params;
      30           0 : }
      31             : 
      32          45 : HenckyIsotropicElasticity::HenckyIsotropicElasticity(const InputParameters & parameters)
      33             :   : LargeDeformationElasticityModel(parameters),
      34             :     DerivativeMaterialPropertyNameInterface(),
      35         135 :     _K(getADMaterialPropertyByName<Real>(prependBaseName("bulk_modulus", true))),
      36          90 :     _G(getADMaterialPropertyByName<Real>(prependBaseName("shear_modulus", true))),
      37             : 
      38          90 :     _d_name(getVar("phase_field", 0)->name()),
      39             : 
      40             :     // The strain energy density and its derivatives
      41          45 :     _psie_name(prependBaseName("strain_energy_density", true)),
      42          45 :     _psie(declareADProperty<Real>(_psie_name)),
      43          45 :     _psie_active(declareADProperty<Real>(_psie_name + "_active")),
      44          90 :     _dpsie_dd(declareADProperty<Real>(derivativePropertyName(_psie_name, {_d_name}))),
      45             : 
      46             :     // The degradation function and its derivatives
      47          45 :     _g_name(prependBaseName("degradation_function", true)),
      48          45 :     _g(getADMaterialProperty<Real>(_g_name)),
      49         135 :     _dg_dd(getADMaterialProperty<Real>(derivativePropertyName(_g_name, {_d_name}))),
      50             : 
      51         135 :     _decomposition(getParam<MooseEnum>("decomposition").getEnum<Decomposition>())
      52             : {
      53         135 : }
      54             : 
      55             : ADRankTwoTensor
      56     1648109 : HenckyIsotropicElasticity::computeMandelStress(const ADRankTwoTensor & Fe,
      57             :                                                const bool plasticity_update)
      58             : {
      59     1648109 :   ADRankTwoTensor stress;
      60             : 
      61     1648109 :   if (_decomposition == Decomposition::none)
      62     3296218 :     stress = computeMandelStressNoDecomposition(Fe, plasticity_update);
      63           0 :   else if (_decomposition == Decomposition::voldev)
      64           0 :     stress = computeMandelStressVolDevDecomposition(Fe, plasticity_update);
      65           0 :   else if (_decomposition == Decomposition::spectral)
      66           0 :     stress = computeMandelStressSpectralDecomposition(Fe, plasticity_update);
      67             :   else
      68           0 :     paramError("decomposition", "Unsupported decomposition type.");
      69             : 
      70     1648109 :   return stress;
      71             : }
      72             : 
      73             : ADRankTwoTensor
      74     1648109 : HenckyIsotropicElasticity::computeMandelStressNoDecomposition(const ADRankTwoTensor & Fe,
      75             :                                                               const bool plasticity_update)
      76             : {
      77             :   ADRankTwoTensor strain = Fe;
      78             :   // If this is called during a plasticity update, we need to first exponentiate Fe, where Fe should
      79             :   // be some plastic flow. The following operations cancel out with an exponentiation of Fe, so we
      80             :   // only do this in the case of exponentiate == false
      81     1648109 :   if (!plasticity_update)
      82     1038128 :     strain = 0.5 * RaccoonUtils::log(Fe.transpose() * Fe);
      83             : 
      84     1648109 :   const ADRankTwoTensor I2(ADRankTwoTensor::initIdentity);
      85             :   // Here, we keep the volumetric part no matter what. But ideally, in the case of J2 plasticity,
      86             :   // the volumetric part of the flow should be zero, and we could save some operations.
      87     6592436 :   ADRankTwoTensor stress_intact = _K[_qp] * strain.trace() * I2 + 2 * _G[_qp] * strain.deviatoric();
      88     1648109 :   ADRankTwoTensor stress = _g[_qp] * stress_intact;
      89             : 
      90             :   // If plasticity_update == false, then we are not in the middle of a plasticity update, hence we
      91             :   // compute the strain energy density
      92     1648109 :   if (!plasticity_update)
      93             :   {
      94             :     // It is convenient that the Mandel stress is conjugate to the log strain
      95     1038128 :     _psie_active[_qp] = 0.5 * stress_intact.doubleContraction(strain);
      96     1038128 :     _psie[_qp] = _g[_qp] * _psie_active[_qp];
      97     1038128 :     _dpsie_dd[_qp] = _dg_dd[_qp] * _psie_active[_qp];
      98             :   }
      99             : 
     100     1648109 :   return stress;
     101             : }
     102             : 
     103             : ADRankTwoTensor
     104           0 : HenckyIsotropicElasticity::computeMandelStressVolDevDecomposition(const ADRankTwoTensor & Fe,
     105             :                                                                   const bool plasticity_update)
     106             : {
     107             :   ADRankTwoTensor strain = Fe;
     108             :   // If this is called during a plasticity update, we need to first exponentiate Fe, where Fe should
     109             :   // be some plastic flow. The following operations cancel out with an exponentiation of Fe, so we
     110             :   // only do this in the case of plasticity_update == false
     111           0 :   if (!plasticity_update)
     112           0 :     strain = 0.5 * RaccoonUtils::log(Fe.transpose() * Fe);
     113             : 
     114           0 :   const ADRankTwoTensor I2(ADRankTwoTensor::initIdentity);
     115             :   // Here, we keep the volumetric part no matter what. But ideally, in the case of J2 plasticity,
     116             :   // the volumetric part of the flow should be zero, and we could save some operations.
     117           0 :   ADReal strain_tr = strain.trace();
     118           0 :   ADReal strain_tr_pos = RaccoonUtils::Macaulay(strain_tr);
     119             :   ADReal strain_tr_neg = strain_tr - strain_tr_pos;
     120           0 :   ADRankTwoTensor strain_dev = strain.deviatoric();
     121             : 
     122           0 :   ADRankTwoTensor stress_intact = _K[_qp] * strain_tr * I2 + 2 * _G[_qp] * strain_dev;
     123           0 :   ADRankTwoTensor stress_neg = _K[_qp] * strain_tr_neg * I2;
     124           0 :   ADRankTwoTensor stress_pos = stress_intact - stress_neg;
     125           0 :   ADRankTwoTensor stress = _g[_qp] * stress_pos + stress_neg;
     126             : 
     127             :   // If plasticity_update == false, then we are not in the middle of a plasticity update, hence we
     128             :   // compute the strain energy density
     129           0 :   if (!plasticity_update)
     130             :   {
     131             :     ADReal psie_intact =
     132           0 :         0.5 * _K[_qp] * strain_tr * strain_tr + _G[_qp] * strain_dev.doubleContraction(strain_dev);
     133           0 :     ADReal psie_inactive = 0.5 * _K[_qp] * strain_tr_neg * strain_tr_neg;
     134           0 :     _psie_active[_qp] = psie_intact - psie_inactive;
     135           0 :     _psie[_qp] = _g[_qp] * _psie_active[_qp] + psie_inactive;
     136           0 :     _dpsie_dd[_qp] = _dg_dd[_qp] * _psie_active[_qp];
     137             :   }
     138             : 
     139           0 :   return stress;
     140             : }
     141             : 
     142             : ADRankTwoTensor
     143           0 : HenckyIsotropicElasticity::computeMandelStressSpectralDecomposition(const ADRankTwoTensor & Fe,
     144             :                                                                     const bool plasticity_update)
     145             : {
     146             :   ADRankTwoTensor strain = Fe;
     147             :   // If this is called during a plasticity update, we need to first exponentiate Fe, where Fe should
     148             :   // be some plastic flow. The following operations cancel out with an exponentiation of Fe, so we
     149             :   // only do this in the case of plasticity_update == false
     150           0 :   if (!plasticity_update)
     151           0 :     strain = 0.5 * RaccoonUtils::log(Fe.transpose() * Fe);
     152             : 
     153           0 :   const ADReal lambda = _K[_qp] - 2 * _G[_qp] / LIBMESH_DIM;
     154           0 :   const ADRankTwoTensor I2(ADRankTwoTensor::initIdentity);
     155           0 :   ADReal strain_tr = strain.trace();
     156           0 :   ADReal strain_tr_pos = RaccoonUtils::Macaulay(strain_tr);
     157             : 
     158             :   // Spectral decomposition
     159           0 :   ADRankTwoTensor strain_pos = RaccoonUtils::spectralDecomposition(strain);
     160             : 
     161             :   // Stress
     162           0 :   ADRankTwoTensor stress_intact = _K[_qp] * strain.trace() * I2 + 2 * _G[_qp] * strain.deviatoric();
     163           0 :   ADRankTwoTensor stress_pos = lambda * strain_tr_pos * I2 + 2 * _G[_qp] * strain_pos;
     164           0 :   ADRankTwoTensor stress_neg = stress_intact - stress_pos;
     165           0 :   ADRankTwoTensor stress = _g[_qp] * stress_pos + stress_neg;
     166             :   // If plasticity_update == false, then we are not in the middle of a plasticity update, hence we
     167             :   // compute the strain energy density
     168           0 :   if (!plasticity_update)
     169             :   {
     170             :     ADReal psie_intact =
     171           0 :         0.5 * lambda * strain_tr * strain_tr + _G[_qp] * strain.doubleContraction(strain);
     172           0 :     _psie_active[_qp] = 0.5 * lambda * strain_tr_pos * strain_tr_pos +
     173           0 :                         _G[_qp] * strain_pos.doubleContraction(strain_pos);
     174             :     ADReal psie_inactive = psie_intact - _psie_active[_qp];
     175           0 :     _psie[_qp] = _g[_qp] * _psie_active[_qp] + psie_inactive;
     176           0 :     _dpsie_dd[_qp] = _dg_dd[_qp] * _psie_active[_qp];
     177             :   }
     178           0 :   return stress;
     179             : }