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 "KLBFNucleationMicroForce.h"
       6             : 
       7             : registerMooseObjectReplaced("raccoonApp",
       8             :                             KLBFNucleationMicroForce,
       9             :                             "12/31/2024 23:59",
      10             :                             LDLNucleationMicroForce);
      11             : 
      12             : InputParameters
      13         140 : KLBFNucleationMicroForce::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. (2020)");
      19         280 :   params.addRequiredParam<MaterialPropertyName>(
      20             :       "tensile_strength", "The tensile strength of the material beyond which the material fails.");
      21         280 :   params.addRequiredParam<MaterialPropertyName>(
      22             :       "compressive_strength",
      23             :       "The compressive strength of the material beyond which the material fails.");
      24         280 :   params.addRequiredParam<MaterialPropertyName>("delta", "delta");
      25         140 :   return params;
      26           0 : }
      27             : 
      28           3 : KLBFNucleationMicroForce::KLBFNucleationMicroForce(const InputParameters & parameters)
      29             :   : NucleationMicroForceBase(parameters),
      30           9 :     _sigma_ts(getADMaterialProperty<Real>(prependBaseName("tensile_strength", true))),
      31           9 :     _sigma_cs(getADMaterialProperty<Real>(prependBaseName("compressive_strength", true))),
      32           9 :     _delta(getADMaterialProperty<Real>(prependBaseName("delta", true)))
      33             : {
      34           3 : }
      35             : 
      36             : void
      37       30000 : KLBFNucleationMicroForce::computeQpProperties()
      38             : {
      39             :   using std::pow;
      40             :   using std::sqrt;
      41             :   // The bulk modulus
      42       60000 :   ADReal K = _lambda[_qp] + 2 * _mu[_qp] / 3;
      43             : 
      44             :   // The mobility
      45       30000 :   ADReal M = _Gc[_qp] / _L[_qp] / _c0[_qp];
      46             : 
      47             :   // Invariants of the stress
      48       30000 :   ADReal I1 = _stress[_qp].trace();
      49       30000 :   ADRankTwoTensor stress_dev = _stress[_qp].deviatoric();
      50       60000 :   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       30000 :   if (MooseUtils::absoluteFuzzyEqual(J2, 0))
      58        1200 :     J2.value() = libMesh::TOLERANCE * libMesh::TOLERANCE;
      59             : 
      60             :   // Parameters in the strength surface
      61             :   ADReal gamma_0 =
      62       60000 :       (_mu[_qp] + 3 * K) * _sigma_ts[_qp] * _L[_qp] / _Gc[_qp] / 18 / _mu[_qp] / _mu[_qp] / K / K;
      63       60000 :   ADReal gamma_1 = (1.0 + _delta[_qp]) / (2.0 * _sigma_ts[_qp] * _sigma_cs[_qp]);
      64       90000 :   ADReal gamma_2 = (8 * _mu[_qp] + 24 * K - 27 * _sigma_ts[_qp]) / 144 / _mu[_qp] / K;
      65       30000 :   ADReal beta_0 = _delta[_qp] * M;
      66       60000 :   ADReal beta_1 = (-gamma_1 * M - gamma_2) * (_sigma_cs[_qp] - _sigma_ts[_qp]) -
      67       30000 :                   gamma_0 * (pow(_sigma_cs[_qp], 3) - pow(_sigma_ts[_qp], 3));
      68      120000 :   ADReal beta_2 = sqrt(3.0) * ((-gamma_1 * M + gamma_2) * (_sigma_cs[_qp] + _sigma_ts[_qp]) +
      69       30000 :                                gamma_0 * (pow(_sigma_cs[_qp], 3) + pow(_sigma_ts[_qp], 3)));
      70       60000 :   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       90000 :   _ex_driving[_qp] = (beta_2 * sqrt(J2) + beta_1 * I1 + beta_0) / (1 + beta_3 * I1 * I1);
      74             : 
      75       90000 :   _stress_balance[_qp] = J2 / _mu[_qp] + pow(I1, 2) / 9.0 / K - _ex_driving[_qp] - M;
      76       30000 : }