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 "LDLNucleationMicroForce.h"
       6             : 
       7             : registerADMooseObject("raccoonApp", LDLNucleationMicroForce);
       8             : 
       9             : InputParameters
      10         142 : LDLNucleationMicroForce::validParams()
      11             : {
      12         142 :   InputParameters params = NucleationMicroForceBase::validParams();
      13             : 
      14         142 :   params.addClassDescription(
      15             :       "This class computes the external driving force for nucleation given "
      16             :       "a Drucker-Prager strength envelope developed by Larsen et al. (2024)");
      17         284 :   params.addRequiredParam<MaterialPropertyName>(
      18             :       "tensile_strength", "The tensile strength of the material beyond which the material fails.");
      19         284 :   params.addRequiredParam<MaterialPropertyName>(
      20             :       "hydrostatic_strength",
      21             :       "The hydrostatic strength of the material beyond which the material fails.");
      22         284 :   params.addParam<MaterialPropertyName>("delta", "delta", "Name of the unitless coefficient delta");
      23         284 :   params.addParam<bool>("h_correction", false, "Whether to use h correction formula for delta");
      24         142 :   return params;
      25           0 : }
      26             : 
      27           3 : LDLNucleationMicroForce::LDLNucleationMicroForce(const InputParameters & parameters)
      28             :   : NucleationMicroForceBase(parameters),
      29           9 :     _sigma_ts(getADMaterialProperty<Real>(prependBaseName("tensile_strength", true))),
      30           6 :     _sigma_hs(getADMaterialProperty<Real>(prependBaseName("hydrostatic_strength", true))),
      31           3 :     _delta(declareADProperty<Real>(prependBaseName("delta", true))),
      32           9 :     _h_correction(getParam<bool>("h_correction"))
      33             : {
      34           3 : }
      35             : 
      36             : void
      37      486960 : LDLNucleationMicroForce::computeQpProperties()
      38             : {
      39             :   // The bulk modulus
      40      973920 :   ADReal K = _lambda[_qp] + 2.0 * _mu[_qp] / 3.0;
      41             : 
      42             :   // The Young's modulus
      43     1460880 :   ADReal E = 9.0 * _mu[_qp] * K / (_mu[_qp] + 3.0 * K);
      44             : 
      45             :   // The mobility
      46      486960 :   ADReal M = _Gc[_qp] / _L[_qp] / _c0[_qp];
      47             : 
      48             :   // Invariants of the stress
      49      486960 :   ADReal I1 = _stress[_qp].trace();
      50      486960 :   ADRankTwoTensor stress_dev = _stress[_qp].deviatoric();
      51      973920 :   ADReal J2 = 0.5 * stress_dev.doubleContraction(stress_dev);
      52             : 
      53             :   // Just to be extra careful... J2 is for sure non-negative but discretization and interpolation
      54             :   // might bring surprise
      55      486960 :   if (J2 < 0)
      56           0 :     mooseException("Negative J2");
      57             : 
      58             :   // define zero J2's derivative
      59      486960 :   if (MooseUtils::absoluteFuzzyEqual(J2, 0))
      60        1200 :     J2.value() = libMesh::TOLERANCE * libMesh::TOLERANCE;
      61             : 
      62             :   // Compute critical energy
      63      486960 :   ADReal W_ts = _sigma_ts[_qp] * _sigma_ts[_qp] / 2.0 / E;
      64      486960 :   ADReal W_hs = _sigma_hs[_qp] * _sigma_hs[_qp] / 2.0 / K;
      65             : 
      66             :   // Compute delta
      67      486960 :   if (!_h_correction)
      68             :   {
      69             :     // Use formula without h correction
      70           0 :     _delta[_qp] = (_sigma_ts[_qp] + (1 + 2 * std::sqrt(3)) * _sigma_hs[_qp]) /
      71           0 :                       (8 + 3 * std::sqrt(3)) / _sigma_hs[_qp] * 3.0 / 16.0 *
      72           0 :                       (_Gc[_qp] / W_ts / _L[_qp]) +
      73           0 :                   3.0 / 8.0;
      74             :   }
      75             :   else
      76             :   {
      77             :     // Get mesh size of current element
      78      486960 :     ADReal h = _current_elem->hmin();
      79             : 
      80             :     // Use formula with h correction
      81     1947840 :     _delta[_qp] = std::pow(1 + 3.0 / 8.0 * h / _L[_qp], -2) *
      82      973920 :                       (_sigma_ts[_qp] + (1 + 2 * std::sqrt(3.0)) * _sigma_hs[_qp]) /
      83      486960 :                       (8 + 3 * std::sqrt(3.0)) / _sigma_hs[_qp] * 3 / 16 *
      84      486960 :                       (_Gc[_qp] / W_ts / _L[_qp]) +
      85     1947840 :                   std::pow(1 + 3.0 / 8.0 * h / _L[_qp], -1) * 2 / 5;
      86             :   }
      87             : 
      88             :   // Parameters in the strength surface
      89             :   ADReal alpha_1 =
      90      973920 :       -_delta[_qp] * _Gc[_qp] / 8.0 / _sigma_hs[_qp] / _L[_qp] + 2.0 / 3.0 * W_hs / _sigma_hs[_qp];
      91      973920 :   ADReal alpha_2 = -(std::sqrt(3.0) / 8.0 * _delta[_qp] * (3.0 * _sigma_hs[_qp] - _sigma_ts[_qp]) /
      92      486960 :                          (_sigma_hs[_qp] * _sigma_ts[_qp]) * _Gc[_qp] / _L[_qp] +
      93      486960 :                      2.0 / std::sqrt(3.0) * W_hs / _sigma_hs[_qp] -
      94      973920 :                      2.0 * std::sqrt(3.0) * W_ts / _sigma_ts[_qp]);
      95             : 
      96             :   // Compute the external driving force required to recover the desired strength envelope.
      97      486960 :   _ex_driving[_qp] =
      98      486960 :       alpha_2 * std::sqrt(J2) + alpha_1 * I1 +
      99     1460880 :       (1.0 - std::sqrt(I1 * I1) / I1) / std::pow(_g[_qp], 1.5) *
     100     2434800 :           (J2 / 2.0 / _mu[_qp] + I1 * I1 / 6.0 / (3.0 * _lambda[_qp] + 2.0 * _mu[_qp]));
     101             : 
     102      486960 :   _stress_balance[_qp] =
     103     1947840 :       J2 / _mu[_qp] + std::pow(I1, 2) / 9.0 / K - _ex_driving[_qp] - M * _delta[_qp];
     104      486960 : }