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 : }