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         142 : KLBFNucleationMicroForce::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. (2020)");
      19         284 :   params.addRequiredParam<MaterialPropertyName>(
      20             :       "tensile_strength", "The tensile strength of the material beyond which the material fails.");
      21         284 :   params.addRequiredParam<MaterialPropertyName>(
      22             :       "compressive_strength",
      23             :       "The compressive strength of the material beyond which the material fails.");
      24         284 :   params.addRequiredParam<MaterialPropertyName>("delta", "delta");
      25         142 :   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             :   // The bulk modulus
      40       60000 :   ADReal K = _lambda[_qp] + 2 * _mu[_qp] / 3;
      41             : 
      42             :   // The mobility
      43       30000 :   ADReal M = _Gc[_qp] / _L[_qp] / _c0[_qp];
      44             : 
      45             :   // Invariants of the stress
      46       30000 :   ADReal I1 = _stress[_qp].trace();
      47       30000 :   ADRankTwoTensor stress_dev = _stress[_qp].deviatoric();
      48       60000 :   ADReal J2 = 0.5 * stress_dev.doubleContraction(stress_dev);
      49             : 
      50             :   // Just to be extra careful... J2 is for sure non-negative but descritization and interpolation
      51             :   // might bring surprise
      52             :   mooseAssert(J2 >= 0, "Negative J2");
      53             : 
      54             :   // define zero J2's derivative
      55       30000 :   if (MooseUtils::absoluteFuzzyEqual(J2, 0))
      56        1200 :     J2.value() = libMesh::TOLERANCE * libMesh::TOLERANCE;
      57             : 
      58             :   // Parameters in the strength surface
      59             :   ADReal gamma_0 =
      60       60000 :       (_mu[_qp] + 3 * K) * _sigma_ts[_qp] * _L[_qp] / _Gc[_qp] / 18 / _mu[_qp] / _mu[_qp] / K / K;
      61       60000 :   ADReal gamma_1 = (1.0 + _delta[_qp]) / (2.0 * _sigma_ts[_qp] * _sigma_cs[_qp]);
      62       90000 :   ADReal gamma_2 = (8 * _mu[_qp] + 24 * K - 27 * _sigma_ts[_qp]) / 144 / _mu[_qp] / K;
      63       30000 :   ADReal beta_0 = _delta[_qp] * M;
      64       60000 :   ADReal beta_1 = (-gamma_1 * M - gamma_2) * (_sigma_cs[_qp] - _sigma_ts[_qp]) -
      65       30000 :                   gamma_0 * (std::pow(_sigma_cs[_qp], 3) - std::pow(_sigma_ts[_qp], 3));
      66             :   ADReal beta_2 =
      67      120000 :       std::sqrt(3.0) * ((-gamma_1 * M + gamma_2) * (_sigma_cs[_qp] + _sigma_ts[_qp]) +
      68       30000 :                         gamma_0 * (std::pow(_sigma_cs[_qp], 3) + std::pow(_sigma_ts[_qp], 3)));
      69       60000 :   ADReal beta_3 = _L[_qp] * _sigma_ts[_qp] / _mu[_qp] / K / _Gc[_qp];
      70             : 
      71             :   // Compute the external driving force required to recover the desired strength envelope.
      72       90000 :   _ex_driving[_qp] = (beta_2 * std::sqrt(J2) + beta_1 * I1 + beta_0) / (1 + beta_3 * I1 * I1);
      73             : 
      74       90000 :   _stress_balance[_qp] = J2 / _mu[_qp] + std::pow(I1, 2) / 9.0 / K - _ex_driving[_qp] - M;
      75       30000 : }