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 "EigenADReal.h"
       6             : #include "LargeDeformationJ2PlasticityBeBar.h"
       7             : #include "CNHIsotropicElasticity.h"
       8             : #include "RaccoonUtils.h"
       9             : 
      10             : registerMooseObject("raccoonApp", LargeDeformationJ2PlasticityBeBar);
      11             : 
      12             : InputParameters
      13         144 : LargeDeformationJ2PlasticityBeBar::validParams()
      14             : {
      15         144 :   InputParameters params = LargeDeformationPlasticityModel::validParams();
      16         144 :   params.addClassDescription("Large deformation $J_2$ plasticity using the bebar (modified "
      17             :                              "left Cauchy-Green) update. Elastic parameters are sourced from "
      18             :                              "the associated CNHIsotropicElasticity model.");
      19         288 :   params.addRequiredCoupledVar("phase_field", "Name of the phase-field (damage) variable");
      20         288 :   params.addParam<MaterialPropertyName>(
      21             :       "strain_energy_density_corr",
      22             :       "psie_corr",
      23             :       "Name of the strain energy density computed by this material model");
      24         288 :   params.addParam<bool>("apply_strain_energy_split",
      25         288 :                         true,
      26             :                         "Whether to apply volumetric/deviatoric split to the active strain "
      27             :                         "energy (psie_active). If true (default), only the positive part of the "
      28             :                         "energy is assigned to psie_active. If false, the full unsplit energy U + "
      29             :                         "W is used for psie_active.");
      30         144 :   return params;
      31           0 : }
      32             : 
      33           6 : LargeDeformationJ2PlasticityBeBar::LargeDeformationJ2PlasticityBeBar(
      34           6 :     const InputParameters & parameters)
      35             :   : LargeDeformationPlasticityModel(parameters),
      36             :     DerivativeMaterialPropertyNameInterface(),
      37           6 :     _bebar(declareADProperty<RankTwoTensor>(prependBaseName("be_bar"))),
      38          18 :     _bebar_old(getMaterialPropertyOldByName<RankTwoTensor>(prependBaseName("be_bar"))),
      39           6 :     _bebar_det(declareADProperty<Real>(prependBaseName("be_bar_det"))),
      40          18 :     _F_old(getMaterialPropertyOldByName<RankTwoTensor>(prependBaseName("deformation_gradient"))),
      41          12 :     _F(getADMaterialPropertyByName<RankTwoTensor>(prependBaseName("deformation_gradient"))),
      42          12 :     _d_name(getVar("phase_field", 0)->name()),
      43           6 :     _psie_name(prependBaseName("strain_energy_density_corr", true)),
      44           6 :     _psie_corr(declareADProperty<Real>(_psie_name)),
      45           6 :     _psie_active_corr(declareADProperty<Real>(_psie_name + "_active")),
      46          12 :     _dpsie_dd_corr(declareADProperty<Real>(derivativePropertyName(_psie_name, {_d_name}))),
      47           6 :     _psie_unsplit(declareADProperty<Real>(_psie_name + "_unsplit")),
      48          24 :     _apply_strain_energy_split(getParam<bool>("apply_strain_energy_split"))
      49             : {
      50           6 :   _check_range = true;
      51          12 : }
      52             : 
      53             : void
      54           6 : LargeDeformationJ2PlasticityBeBar::setElasticityModel(
      55             :     LargeDeformationElasticityModel * elasticity_model)
      56             : {
      57           6 :   auto * cnh = dynamic_cast<CNHIsotropicElasticity *>(elasticity_model);
      58           6 :   if (!cnh)
      59           0 :     mooseError(type(), ": requires a CNHIsotropicElasticity model; got ", elasticity_model->type());
      60           6 :   _K = &cnh->getK();
      61           6 :   _G = &cnh->getG();
      62           6 :   _ge = &cnh->getDegradation();
      63           6 :   _dge_dd = &cnh->getDegradationDerivative();
      64           6 :   _psie_cnh = &cnh->getPsie();
      65           6 :   _psie_active_cnh = &cnh->getPsieActive();
      66           6 :   _dpsie_dd_cnh = &cnh->getDpsieDD();
      67           6 :   LargeDeformationPlasticityModel::setElasticityModel(elasticity_model);
      68           6 : }
      69             : 
      70             : void
      71          32 : LargeDeformationJ2PlasticityBeBar::initQpStatefulProperties()
      72             : {
      73          32 :   _Fp[_qp].setToIdentity();
      74          32 :   _ep[_qp] = 0;
      75          32 :   _bebar[_qp].setToIdentity();
      76          32 : }
      77             : 
      78             : void
      79        1176 : LargeDeformationJ2PlasticityBeBar::updateState(ADRankTwoTensor & stress, ADRankTwoTensor & /*Fe*/)
      80             : {
      81             :   using std::cbrt;
      82             :   using std::sqrt;
      83             : 
      84        1176 :   ADRankTwoTensor I2;
      85        1176 :   I2.setToIdentity();
      86             : 
      87        1176 :   ADReal delta_ep = 0;
      88             : 
      89             :   // Compute incremental fbar
      90        1176 :   ADRankTwoTensor f = _F[_qp] * _F_old[_qp].inverse();
      91        1176 :   ADRankTwoTensor fbar = f / cbrt(f.det());
      92             : 
      93             :   // Compute bebar_trial
      94        1176 :   _bebar[_qp] = fbar * _bebar_old[_qp] * fbar.transpose();
      95             : 
      96             :   // Compute trial deviatoric stress and its norm
      97        2352 :   ADRankTwoTensor s_trial = (*_G)[_qp] * (*_ge)[_qp] * _bebar[_qp].deviatoric();
      98        1176 :   ADReal s_trial_norm = s_trial.doubleContraction(s_trial);
      99        1176 :   if (MooseUtils::absoluteFuzzyEqual(s_trial_norm, 0))
     100           8 :     s_trial_norm.value() = libMesh::TOLERANCE * libMesh::TOLERANCE;
     101        1176 :   s_trial_norm = sqrt(s_trial_norm);
     102             : 
     103        2352 :   _Np[_qp] = sqrt(1.5) * s_trial / s_trial_norm;
     104             : 
     105             :   // Return mapping
     106        1176 :   ADReal phi = computeResidual(s_trial_norm, delta_ep);
     107        1176 :   if (phi > 0)
     108             :   {
     109        1064 :     returnMappingSolve(s_trial_norm, delta_ep, _console);
     110             : 
     111        2128 :     _ep[_qp] = _ep_old[_qp] + delta_ep;
     112        1064 :     ADRankTwoTensor s = s_trial - (2.0 / 3.0) * (*_ge)[_qp] * (*_G)[_qp] * delta_ep *
     113        2128 :                                       _bebar[_qp].trace() * _Np[_qp];
     114             : 
     115        1064 :     ADReal J = _F[_qp].det();
     116        2128 :     ADReal p = 0.5 * (*_K)[_qp] * (J * J - 1);
     117        4256 :     ADRankTwoTensor tau = J >= 1.0 ? (*_ge)[_qp] * p * I2 + s : p * I2 + s;
     118        1064 :     stress = tau / J;
     119             : 
     120        1064 :     computeCorrectionTerm(s / (*_ge)[_qp] / (*_G)[_qp]);
     121             :   }
     122             :   else
     123             :   {
     124         112 :     _ep[_qp] = _ep_old[_qp];
     125             : 
     126         112 :     ADReal J = _F[_qp].det();
     127         224 :     ADReal p = 0.5 * (*_K)[_qp] * (J * J - 1);
     128         448 :     ADRankTwoTensor tau = J >= 1.0 ? (*_ge)[_qp] * p * I2 + s_trial : p * I2 + s_trial;
     129         224 :     stress = tau / J;
     130             :   }
     131             : 
     132        1176 :   computeStrainEnergyDensity();
     133        1176 :   _hardening_model->plasticEnergy(_ep[_qp]);
     134        1176 :   _hardening_model->plasticDissipation(delta_ep, _ep[_qp], 0);
     135             : 
     136        1176 :   if (_t_step > 0)
     137             :   {
     138        2336 :     _heat[_qp] = _hardening_model->plasticDissipation(delta_ep, _ep[_qp], 1) * delta_ep / _dt;
     139        2336 :     _heat[_qp] += _hardening_model->thermalConjugate(_ep[_qp]) * delta_ep / _dt;
     140             :   }
     141             :   else
     142           8 :     _heat[_qp] = 0;
     143             : 
     144        1176 :   _bebar_det[_qp] = _bebar[_qp].det();
     145        1176 : }
     146             : 
     147             : Real
     148        2808 : LargeDeformationJ2PlasticityBeBar::computeReferenceResidual(const ADReal & effective_trial_stress,
     149             :                                                             const ADReal & delta_ep)
     150             : {
     151        2808 :   return MetaPhysicL::raw_value(effective_trial_stress - std::sqrt(2.0 / 3.0) * (*_ge)[_qp] *
     152        2808 :                                                              (*_G)[_qp] * delta_ep *
     153        5616 :                                                              _bebar[_qp].trace());
     154             : }
     155             : 
     156             : ADReal
     157        3984 : LargeDeformationJ2PlasticityBeBar::computeResidual(const ADReal & effective_trial_stress,
     158             :                                                    const ADReal & delta_ep)
     159             : {
     160        3984 :   ADReal ep = _ep_old[_qp] + delta_ep;
     161        3984 :   return effective_trial_stress -
     162        7968 :          std::sqrt(2.0 / 3.0) * (_hardening_model->plasticEnergy(ep, 1) +
     163        3984 :                                  _hardening_model->plasticDissipation(delta_ep, ep, 1)) -
     164        7968 :          std::sqrt(2.0 / 3.0) * (*_ge)[_qp] * (*_G)[_qp] * delta_ep * _bebar[_qp].trace();
     165             : }
     166             : 
     167             : ADReal
     168        2808 : LargeDeformationJ2PlasticityBeBar::computeDerivative(const ADReal & /*effective_trial_stress*/,
     169             :                                                      const ADReal & delta_ep)
     170             : {
     171        2808 :   ADReal ep = _ep_old[_qp] + delta_ep;
     172        5616 :   return -std::sqrt(2.0 / 3.0) * (_hardening_model->plasticEnergy(ep, 2) +
     173        2808 :                                   _hardening_model->plasticDissipation(delta_ep, ep, 2)) -
     174        5616 :          std::sqrt(2.0 / 3.0) * (*_ge)[_qp] * (*_G)[_qp] * _bebar[_qp].trace();
     175             : }
     176             : 
     177             : void
     178        1064 : LargeDeformationJ2PlasticityBeBar::computeCorrectionTerm(const ADRankTwoTensor & devbebar)
     179             : {
     180             :   using std::cbrt;
     181             :   using std::sqrt;
     182             : 
     183        1064 :   ADRankTwoTensor I2(RankTwoTensorTempl<ADReal>::initIdentity);
     184             : 
     185        1064 :   ADReal a = devbebar(0, 0);
     186        1064 :   ADReal b = devbebar(1, 1);
     187        1064 :   ADReal c = devbebar(2, 2);
     188        1064 :   ADReal d = devbebar(1, 2);
     189        1064 :   ADReal e = devbebar(0, 2);
     190        1064 :   ADReal h = devbebar(0, 1);
     191             : 
     192        1064 :   ADReal A = a + b + c;
     193        1064 :   ADReal B = a * b + a * c + b * c - d * d - e * e - h * h;
     194        1064 :   ADReal C = a * b * c + 2.0 * d * e * h - a * d * d - b * e * e - c * h * h - 1.0;
     195             : 
     196             :   ADReal D = cbrt(
     197        1064 :       -2 * A * A * A +
     198        1064 :       3 * sqrt(3) *
     199        5320 :           sqrt(4 * A * A * A * C - A * A * B * B - 18 * A * B * C + 4 * B * B * B + 27 * C * C) +
     200        3192 :       9 * A * B - 27 * C);
     201             : 
     202        6384 :   ADReal Ie_bar = D / 3 / cbrt(2) - cbrt(2) * (3 * B - A * A) / 3 / D - A / 3;
     203        1064 :   _bebar[_qp] = devbebar + Ie_bar * I2;
     204        1064 : }
     205             : 
     206             : void
     207        1176 : LargeDeformationJ2PlasticityBeBar::computeStrainEnergyDensity()
     208             : {
     209             :   using std::log;
     210        1176 :   ADReal J = _F[_qp].det();
     211        3528 :   ADReal U = 0.5 * (*_K)[_qp] * (0.5 * (J * J - 1) - log(J));
     212        3528 :   ADReal W = 0.5 * (*_G)[_qp] * (_bebar[_qp].trace() - 3.0);
     213             : 
     214        1176 :   _psie_unsplit[_qp] = U + W;
     215             : 
     216        1176 :   ADReal E_el_pos = J >= 1.0 ? U + W : W;
     217        1176 :   ADReal E_el_neg = J >= 1.0 ? 0.0 : U;
     218             : 
     219        1176 :   _psie_active_corr[_qp] = _apply_strain_energy_split ? E_el_pos : _psie_unsplit[_qp];
     220        2352 :   _psie_corr[_qp] = (*_ge)[_qp] * E_el_pos + E_el_neg;
     221        2352 :   _dpsie_dd_corr[_qp] = (*_dge_dd)[_qp] * _psie_active_corr[_qp];
     222             : 
     223             :   // Overwrite CNH's energy properties so that psie_active is correct for fracture driving.
     224             :   // CNH computes these from Fe, which BeBar does not update; the bebar-based values are correct.
     225        1176 :   (*_psie_cnh)[_qp] = _psie_corr[_qp];
     226        1176 :   (*_psie_active_cnh)[_qp] = _psie_active_corr[_qp];
     227        1176 :   (*_dpsie_dd_cnh)[_qp] = _dpsie_dd_corr[_qp];
     228        1176 : }