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         152 : LDLNucleationMicroForce::validParams()
      11             : {
      12         152 :   InputParameters params = NucleationMicroForceBase::validParams();
      13             : 
      14         152 :   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         304 :   params.addRequiredParam<MaterialPropertyName>(
      18             :       "tensile_strength", "The tensile strength of the material beyond which the material fails.");
      19         304 :   params.addRequiredParam<MaterialPropertyName>(
      20             :       "hydrostatic_strength",
      21             :       "The hydrostatic strength of the material beyond which the material fails.");
      22         304 :   params.addParam<MaterialPropertyName>("delta", "delta", "Name of the unitless coefficient delta");
      23         304 :   params.addParam<bool>("h_correction", false, "Whether to use h correction formula for delta");
      24         152 :   return params;
      25           0 : }
      26             : 
      27          12 : LDLNucleationMicroForce::LDLNucleationMicroForce(const InputParameters & parameters)
      28             :   : NucleationMicroForceBase(parameters),
      29          36 :     _sigma_ts(getADMaterialProperty<Real>(prependBaseName("tensile_strength", true))),
      30          24 :     _sigma_hs(getADMaterialProperty<Real>(prependBaseName("hydrostatic_strength", true))),
      31          12 :     _delta(declareADProperty<Real>(prependBaseName("delta", true))),
      32          36 :     _h_correction(getParam<bool>("h_correction"))
      33             : {
      34          12 : }
      35             : 
      36             : void
      37     1035840 : LDLNucleationMicroForce::computeQpProperties()
      38             : {
      39             :   using std::pow;
      40             :   using std::sqrt;
      41             :   // The bulk modulus
      42     2071680 :   ADReal K = _lambda[_qp] + 2.0 * _mu[_qp] / 3.0;
      43             : 
      44             :   // The Young's modulus
      45     3107520 :   ADReal E = 9.0 * _mu[_qp] * K / (_mu[_qp] + 3.0 * K);
      46             : 
      47             :   // The mobility
      48     1035840 :   ADReal M = _Gc[_qp] / _L[_qp] / _c0[_qp];
      49             : 
      50     1035840 :   ADReal l_ch = E * _Gc[_qp] / _sigma_ts[_qp] / _sigma_ts[_qp] * 3.0 / 8.0;
      51     2071680 :   if (_L[_qp] > l_ch * 0.8)
      52       35761 :     flagSolutionWarning(
      53             :         "The reg length might be too large to construct a reasonable numerical strength surface.");
      54             : 
      55             :   // Invariants of the stress
      56     1035840 :   ADReal I1 = _stress[_qp].trace();
      57     1035840 :   ADRankTwoTensor stress_dev = _stress[_qp].deviatoric();
      58     2071680 :   ADReal J2 = 0.5 * stress_dev.doubleContraction(stress_dev);
      59             : 
      60             :   // Just to be extra careful... J2 is for sure non-negative but discretization and interpolation
      61             :   // might bring surprise
      62     1035840 :   if (J2 < 0)
      63           0 :     mooseException("Negative J2");
      64             : 
      65             :   // define zero J2's derivative
      66     1035840 :   if (MooseUtils::absoluteFuzzyEqual(J2, 0))
      67        4800 :     J2.value() = libMesh::TOLERANCE * libMesh::TOLERANCE;
      68             : 
      69             :   // Compute critical energy
      70     1035840 :   ADReal W_ts = _sigma_ts[_qp] * _sigma_ts[_qp] / 2.0 / E;
      71     1035840 :   ADReal W_hs = _sigma_hs[_qp] * _sigma_hs[_qp] / 2.0 / K;
      72             : 
      73             :   // Compute delta
      74     1035840 :   if (!_h_correction)
      75             :   {
      76             :     // Use formula without h correction
      77           0 :     _delta[_qp] = (_sigma_ts[_qp] + (1 + 2 * sqrt(3)) * _sigma_hs[_qp]) / (8 + 3 * sqrt(3)) /
      78           0 :                       _sigma_hs[_qp] * 3.0 / 16.0 * (_Gc[_qp] / W_ts / _L[_qp]) +
      79           0 :                   3.0 / 8.0;
      80             :   }
      81             :   else
      82             :   {
      83             :     using std::pow;
      84             :     // Get mesh size of current element
      85     1035840 :     ADReal h = _current_elem->hmin();
      86             :     // if (_L[_qp]/h *_L[_qp]/l_ch > 0.16)
      87     4071840 :     if ((h > _L[_qp] * 0.5) && (_L[_qp] > l_ch * 0.7))
      88     1000083 :       flagSolutionWarning(
      89             :           "The mesh size might be too coarse for a valid h_correction in the complete "
      90             :           "model and lead to unexpected numerical strength surface. You may either refine "
      91             :           "the mesh, reduce the reg length, or turn h_correction=false.");
      92             : 
      93             :     // Use formula with h correction
      94     4143360 :     _delta[_qp] = pow(1 + 3.0 / 8.0 * h / _L[_qp], -2) *
      95     2071680 :                       (_sigma_ts[_qp] + (1 + 2 * sqrt(3.0)) * _sigma_hs[_qp]) /
      96     2071680 :                       (8 + 3 * sqrt(3.0)) / _sigma_hs[_qp] * 3 / 16 * (_Gc[_qp] / W_ts / _L[_qp]) +
      97     4143360 :                   pow(1 + 3.0 / 8.0 * h / _L[_qp], -1) * 2 / 5;
      98             :   }
      99             : 
     100             :   // Parameters in the strength surface
     101             :   ADReal alpha_1 =
     102     2071680 :       -_delta[_qp] * _Gc[_qp] / 8.0 / _sigma_hs[_qp] / _L[_qp] + 2.0 / 3.0 * W_hs / _sigma_hs[_qp];
     103             :   ADReal alpha_2 =
     104     2071680 :       -(sqrt(3.0) / 8.0 * _delta[_qp] * (3.0 * _sigma_hs[_qp] - _sigma_ts[_qp]) /
     105     1035840 :             (_sigma_hs[_qp] * _sigma_ts[_qp]) * _Gc[_qp] / _L[_qp] +
     106     3107520 :         2.0 / sqrt(3.0) * W_hs / _sigma_hs[_qp] - 2.0 * sqrt(3.0) * W_ts / _sigma_ts[_qp]);
     107             : 
     108             :   // Compute the external driving force required to recover the desired strength envelope.
     109     1035840 :   _ex_driving[_qp] =
     110     1035840 :       alpha_2 * sqrt(J2) + alpha_1 * I1 +
     111     1040640 :       (I1 > 0 ? 0 : 2) / pow(_g[_qp], 1.5) *
     112     5179200 :           (J2 / 2.0 / _mu[_qp] + I1 * I1 / 6.0 / (3.0 * _lambda[_qp] + 2.0 * _mu[_qp]));
     113             : 
     114     4143360 :   _stress_balance[_qp] = J2 / _mu[_qp] + pow(I1, 2) / 9.0 / K - _ex_driving[_qp] - M * _delta[_qp];
     115     1035840 : }