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 "KLRNucleationMicroForce.h"
       6             : 
       7             : registerMooseObjectReplaced("raccoonApp",
       8             :                             KLRNucleationMicroForce,
       9             :                             "12/31/2024 23:59",
      10             :                             LDLNucleationMicroForce);
      11             : 
      12             : InputParameters
      13         142 : KLRNucleationMicroForce::validParams()
      14             : {
      15         142 :   InputParameters params = NucleationMicroForceBase::validParams();
      16             : 
      17         142 :   params.addClassDescription("This class computes the external driving force for nucleation given "
      18             :                              "a Drucker-Prager strength envelope developed by Kumar et al. (2022)");
      19             : 
      20         284 :   params.addRequiredParam<MaterialPropertyName>(
      21             :       "tensile_strength", "The tensile strength of the material beyond which the material fails.");
      22         284 :   params.addRequiredParam<MaterialPropertyName>(
      23             :       "compressive_strength",
      24             :       "The compressive strength of the material beyond which the material fails.");
      25         284 :   params.addRequiredParam<MaterialPropertyName>("delta", "delta");
      26         142 :   return params;
      27           0 : }
      28             : 
      29           3 : KLRNucleationMicroForce::KLRNucleationMicroForce(const InputParameters & parameters)
      30             :   : NucleationMicroForceBase(parameters),
      31           9 :     _sigma_ts(getADMaterialProperty<Real>(prependBaseName("tensile_strength", true))),
      32           9 :     _sigma_cs(getADMaterialProperty<Real>(prependBaseName("compressive_strength", true))),
      33           6 :     _delta(getADMaterialProperty<Real>(prependBaseName("delta", true))),
      34           6 :     _druck_prager_balance(declareADProperty<Real>("druck_prager_balance"))
      35             : {
      36           3 : }
      37             : 
      38             : void
      39      392880 : KLRNucleationMicroForce::computeQpProperties()
      40             : {
      41             :   // The bulk modulus
      42      785760 :   ADReal K = _lambda[_qp] + 2 * _mu[_qp] / 3;
      43             : 
      44             :   // The mobility
      45      392880 :   ADReal M = _Gc[_qp] / _L[_qp] / _c0[_qp];
      46             : 
      47             :   // Invariants of the stress
      48      392880 :   ADReal I1 = _stress[_qp].trace();
      49      392880 :   ADRankTwoTensor stress_dev = _stress[_qp].deviatoric();
      50      785760 :   ADReal J2 = 0.5 * stress_dev.doubleContraction(stress_dev);
      51             : 
      52             :   // Just to be extra careful... J2 is for sure non-negative but descritization and interpolation
      53             :   // might bring surprise
      54             :   mooseAssert(J2 >= 0, "Negative J2");
      55             : 
      56             :   // define zero J2's derivative
      57      392880 :   if (MooseUtils::absoluteFuzzyEqual(J2, 0))
      58        1200 :     J2.value() = libMesh::TOLERANCE * libMesh::TOLERANCE;
      59             : 
      60      392880 :   if (MooseUtils::absoluteFuzzyEqual(I1, 0))
      61        1200 :     I1.value() = libMesh::TOLERANCE;
      62             : 
      63             :   // Parameters in the strength surface
      64     1178640 :   ADReal gamma_0 = _sigma_ts[_qp] / 6.0 / (3.0 * _lambda[_qp] + 2.0 * _mu[_qp]) +
      65      392880 :                    _sigma_ts[_qp] / 6.0 / _mu[_qp];
      66      785760 :   ADReal gamma_1 = (1.0 + _delta[_qp]) / (2.0 * _sigma_ts[_qp] * _sigma_cs[_qp]);
      67      392880 :   ADReal beta_0 = _delta[_qp] * M;
      68      785760 :   ADReal beta_1 = -gamma_1 * M * (_sigma_cs[_qp] - _sigma_ts[_qp]) + gamma_0;
      69     1571520 :   ADReal beta_2 = std::sqrt(3.0) * (-gamma_1 * M * (_sigma_cs[_qp] + _sigma_ts[_qp]) + gamma_0);
      70      392880 :   ADReal beta_3 = _L[_qp] * _sigma_ts[_qp] / _mu[_qp] / K / _Gc[_qp];
      71             : 
      72             :   // Compute the external driving force required to recover the desired strength envelope.
      73      392880 :   _ex_driving[_qp] =
      74      392880 :       (beta_2 * std::sqrt(J2) + beta_1 * I1 + beta_0) +
      75     1178640 :       (1.0 - std::sqrt(I1 * I1) / I1) / std::pow(_g[_qp], 1.5) *
      76     1964400 :           (J2 / 2.0 / _mu[_qp] + I1 * I1 / 6.0 / (3.0 * _lambda[_qp] + 2.0 * _mu[_qp]));
      77             : 
      78     1178640 :   _stress_balance[_qp] = J2 / _mu[_qp] + std::pow(I1, 2) / 9.0 / K - _ex_driving[_qp] - M;
      79             : 
      80      392880 :   _druck_prager_balance[_qp] =
      81      392880 :       std::sqrt(J2) +
      82      785760 :       (_sigma_cs[_qp] - _sigma_ts[_qp]) / std::sqrt(3.0) / (_sigma_cs[_qp] + _sigma_ts[_qp]) * I1 -
      83     1178640 :       2.0 * _sigma_ts[_qp] * _sigma_cs[_qp] / std::sqrt(3.0) / (_sigma_cs[_qp] + _sigma_ts[_qp]);
      84      392880 : }