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         140 : KLRNucleationMicroForce::validParams()
      14             : {
      15         140 :   InputParameters params = NucleationMicroForceBase::validParams();
      16             : 
      17         140 :   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         280 :   params.addRequiredParam<MaterialPropertyName>(
      21             :       "tensile_strength", "The tensile strength of the material beyond which the material fails.");
      22         280 :   params.addRequiredParam<MaterialPropertyName>(
      23             :       "compressive_strength",
      24             :       "The compressive strength of the material beyond which the material fails.");
      25         280 :   params.addRequiredParam<MaterialPropertyName>("delta", "delta");
      26         140 :   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             :   using std::pow;
      42             :   using std::sqrt;
      43             :   // The bulk modulus
      44      785760 :   ADReal K = _lambda[_qp] + 2 * _mu[_qp] / 3;
      45             : 
      46             :   // The mobility
      47      392880 :   ADReal M = _Gc[_qp] / _L[_qp] / _c0[_qp];
      48             : 
      49             :   // Invariants of the stress
      50      392880 :   ADReal I1 = _stress[_qp].trace();
      51      392880 :   ADRankTwoTensor stress_dev = _stress[_qp].deviatoric();
      52      785760 :   ADReal J2 = 0.5 * stress_dev.doubleContraction(stress_dev);
      53             : 
      54             :   // Just to be extra careful... J2 is for sure non-negative but descritization and interpolation
      55             :   // might bring surprise
      56             :   mooseAssert(J2 >= 0, "Negative J2");
      57             : 
      58             :   // define zero J2's derivative
      59      392880 :   if (MooseUtils::absoluteFuzzyEqual(J2, 0))
      60        1200 :     J2.value() = libMesh::TOLERANCE * libMesh::TOLERANCE;
      61             : 
      62      392880 :   if (MooseUtils::absoluteFuzzyEqual(I1, 0))
      63        1200 :     I1.value() = libMesh::TOLERANCE;
      64             : 
      65             :   // Parameters in the strength surface
      66     1178640 :   ADReal gamma_0 = _sigma_ts[_qp] / 6.0 / (3.0 * _lambda[_qp] + 2.0 * _mu[_qp]) +
      67      392880 :                    _sigma_ts[_qp] / 6.0 / _mu[_qp];
      68      785760 :   ADReal gamma_1 = (1.0 + _delta[_qp]) / (2.0 * _sigma_ts[_qp] * _sigma_cs[_qp]);
      69      392880 :   ADReal beta_0 = _delta[_qp] * M;
      70      785760 :   ADReal beta_1 = -gamma_1 * M * (_sigma_cs[_qp] - _sigma_ts[_qp]) + gamma_0;
      71     1571520 :   ADReal beta_2 = sqrt(3.0) * (-gamma_1 * M * (_sigma_cs[_qp] + _sigma_ts[_qp]) + gamma_0);
      72      392880 :   ADReal beta_3 = _L[_qp] * _sigma_ts[_qp] / _mu[_qp] / K / _Gc[_qp];
      73             : 
      74             :   // Compute the external driving force required to recover the desired strength envelope.
      75      392880 :   _ex_driving[_qp] =
      76      392880 :       (beta_2 * sqrt(J2) + beta_1 * I1 + beta_0) +
      77     1178640 :       (1.0 - sqrt(I1 * I1) / I1) / pow(_g[_qp], 1.5) *
      78     1964400 :           (J2 / 2.0 / _mu[_qp] + I1 * I1 / 6.0 / (3.0 * _lambda[_qp] + 2.0 * _mu[_qp]));
      79             : 
      80     1178640 :   _stress_balance[_qp] = J2 / _mu[_qp] + pow(I1, 2) / 9.0 / K - _ex_driving[_qp] - M;
      81             : 
      82      392880 :   _druck_prager_balance[_qp] =
      83      392880 :       sqrt(J2) +
      84      785760 :       (_sigma_cs[_qp] - _sigma_ts[_qp]) / sqrt(3.0) / (_sigma_cs[_qp] + _sigma_ts[_qp]) * I1 -
      85     1178640 :       2.0 * _sigma_ts[_qp] * _sigma_cs[_qp] / sqrt(3.0) / (_sigma_cs[_qp] + _sigma_ts[_qp]);
      86      392880 : }