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             :   // The bulk modulus
      40     2071680 :   ADReal K = _lambda[_qp] + 2.0 * _mu[_qp] / 3.0;
      41             : 
      42             :   // The Young's modulus
      43     3107520 :   ADReal E = 9.0 * _mu[_qp] * K / (_mu[_qp] + 3.0 * K);
      44             : 
      45             :   // The mobility
      46     1035840 :   ADReal M = _Gc[_qp] / _L[_qp] / _c0[_qp];
      47             : 
      48     1035840 :   ADReal l_ch = E * _Gc[_qp] / _sigma_ts[_qp] / _sigma_ts[_qp] * 3.0 / 8.0;
      49     2071680 :   if (_L[_qp] > l_ch * 0.8)
      50       35761 :     flagSolutionWarning(
      51             :         "The reg length might be too large to construct a reasonable numerical strength surface.");
      52             : 
      53             :   // Invariants of the stress
      54     1035840 :   ADReal I1 = _stress[_qp].trace();
      55     1035840 :   ADRankTwoTensor stress_dev = _stress[_qp].deviatoric();
      56     2071680 :   ADReal J2 = 0.5 * stress_dev.doubleContraction(stress_dev);
      57             : 
      58             :   // Just to be extra careful... J2 is for sure non-negative but discretization and interpolation
      59             :   // might bring surprise
      60     1035840 :   if (J2 < 0)
      61           0 :     mooseException("Negative J2");
      62             : 
      63             :   // define zero J2's derivative
      64     1035840 :   if (MooseUtils::absoluteFuzzyEqual(J2, 0))
      65        4800 :     J2.value() = libMesh::TOLERANCE * libMesh::TOLERANCE;
      66             : 
      67             :   // Compute critical energy
      68     1035840 :   ADReal W_ts = _sigma_ts[_qp] * _sigma_ts[_qp] / 2.0 / E;
      69     1035840 :   ADReal W_hs = _sigma_hs[_qp] * _sigma_hs[_qp] / 2.0 / K;
      70             : 
      71             :   // Compute delta
      72     1035840 :   if (!_h_correction)
      73             :   {
      74             :     // Use formula without h correction
      75           0 :     _delta[_qp] = (_sigma_ts[_qp] + (1 + 2 * std::sqrt(3)) * _sigma_hs[_qp]) /
      76           0 :                       (8 + 3 * std::sqrt(3)) / _sigma_hs[_qp] * 3.0 / 16.0 *
      77           0 :                       (_Gc[_qp] / W_ts / _L[_qp]) +
      78           0 :                   3.0 / 8.0;
      79             :   }
      80             :   else
      81             :   {
      82             :     // Get mesh size of current element
      83     1035840 :     ADReal h = _current_elem->hmin();
      84             :     // if (_L[_qp]/h *_L[_qp]/l_ch > 0.16)
      85     4071840 :     if ((h > _L[_qp] * 0.5) && (_L[_qp] > l_ch * 0.7))
      86     1000083 :       flagSolutionWarning(
      87             :           "The mesh size might be too coarse for a valid h_correction in the complete "
      88             :           "model and lead to unexpected numerical strength surface. You may either refine "
      89             :           "the mesh, reduce the reg length, or turn h_correction=false.");
      90             : 
      91             :     // Use formula with h correction
      92     4143360 :     _delta[_qp] = std::pow(1 + 3.0 / 8.0 * h / _L[_qp], -2) *
      93     2071680 :                       (_sigma_ts[_qp] + (1 + 2 * std::sqrt(3.0)) * _sigma_hs[_qp]) /
      94     1035840 :                       (8 + 3 * std::sqrt(3.0)) / _sigma_hs[_qp] * 3 / 16 *
      95     1035840 :                       (_Gc[_qp] / W_ts / _L[_qp]) +
      96     4143360 :                   std::pow(1 + 3.0 / 8.0 * h / _L[_qp], -1) * 2 / 5;
      97             :   }
      98             : 
      99             :   // Parameters in the strength surface
     100             :   ADReal alpha_1 =
     101     2071680 :       -_delta[_qp] * _Gc[_qp] / 8.0 / _sigma_hs[_qp] / _L[_qp] + 2.0 / 3.0 * W_hs / _sigma_hs[_qp];
     102     2071680 :   ADReal alpha_2 = -(std::sqrt(3.0) / 8.0 * _delta[_qp] * (3.0 * _sigma_hs[_qp] - _sigma_ts[_qp]) /
     103     1035840 :                          (_sigma_hs[_qp] * _sigma_ts[_qp]) * _Gc[_qp] / _L[_qp] +
     104     1035840 :                      2.0 / std::sqrt(3.0) * W_hs / _sigma_hs[_qp] -
     105     2071680 :                      2.0 * std::sqrt(3.0) * W_ts / _sigma_ts[_qp]);
     106             : 
     107             :   // Compute the external driving force required to recover the desired strength envelope.
     108     1035840 :   _ex_driving[_qp] =
     109     1035840 :       alpha_2 * std::sqrt(J2) + alpha_1 * I1 +
     110     3107520 :       (1.0 - std::sqrt(I1 * I1) / I1) / std::pow(_g[_qp], 1.5) *
     111     5179200 :           (J2 / 2.0 / _mu[_qp] + I1 * I1 / 6.0 / (3.0 * _lambda[_qp] + 2.0 * _mu[_qp]));
     112             : 
     113     1035840 :   _stress_balance[_qp] =
     114     4143360 :       J2 / _mu[_qp] + std::pow(I1, 2) / 9.0 / K - _ex_driving[_qp] - M * _delta[_qp];
     115     1035840 : }