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 "KLRNucleationMicroForce.h" 6 : 7 : registerMooseObjectReplaced("raccoonApp", 8 : KLRNucleationMicroForce, 9 : "12/31/2024 23:59", 10 : LDLNucleationMicroForce); 11 : 12 : InputParameters 13 142 : KLRNucleationMicroForce::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. (2022)"); 19 : 20 284 : params.addRequiredParam<MaterialPropertyName>( 21 : "tensile_strength", "The tensile strength of the material beyond which the material fails."); 22 284 : params.addRequiredParam<MaterialPropertyName>( 23 : "compressive_strength", 24 : "The compressive strength of the material beyond which the material fails."); 25 284 : params.addRequiredParam<MaterialPropertyName>("delta", "delta"); 26 142 : return params; 27 0 : } 28 : 29 3 : KLRNucleationMicroForce::KLRNucleationMicroForce(const InputParameters & parameters) 30 : : NucleationMicroForceBase(parameters), 31 9 : _sigma_ts(getADMaterialProperty<Real>(prependBaseName("tensile_strength", true))), 32 9 : _sigma_cs(getADMaterialProperty<Real>(prependBaseName("compressive_strength", true))), 33 6 : _delta(getADMaterialProperty<Real>(prependBaseName("delta", true))), 34 6 : _druck_prager_balance(declareADProperty<Real>("druck_prager_balance")) 35 : { 36 3 : } 37 : 38 : void 39 392880 : KLRNucleationMicroForce::computeQpProperties() 40 : { 41 : // The bulk modulus 42 785760 : ADReal K = _lambda[_qp] + 2 * _mu[_qp] / 3; 43 : 44 : // The mobility 45 392880 : ADReal M = _Gc[_qp] / _L[_qp] / _c0[_qp]; 46 : 47 : // Invariants of the stress 48 392880 : ADReal I1 = _stress[_qp].trace(); 49 392880 : ADRankTwoTensor stress_dev = _stress[_qp].deviatoric(); 50 785760 : ADReal J2 = 0.5 * stress_dev.doubleContraction(stress_dev); 51 : 52 : // Just to be extra careful... J2 is for sure non-negative but descritization and interpolation 53 : // might bring surprise 54 : mooseAssert(J2 >= 0, "Negative J2"); 55 : 56 : // define zero J2's derivative 57 392880 : if (MooseUtils::absoluteFuzzyEqual(J2, 0)) 58 1200 : J2.value() = libMesh::TOLERANCE * libMesh::TOLERANCE; 59 : 60 392880 : if (MooseUtils::absoluteFuzzyEqual(I1, 0)) 61 1200 : I1.value() = libMesh::TOLERANCE; 62 : 63 : // Parameters in the strength surface 64 1178640 : ADReal gamma_0 = _sigma_ts[_qp] / 6.0 / (3.0 * _lambda[_qp] + 2.0 * _mu[_qp]) + 65 392880 : _sigma_ts[_qp] / 6.0 / _mu[_qp]; 66 785760 : ADReal gamma_1 = (1.0 + _delta[_qp]) / (2.0 * _sigma_ts[_qp] * _sigma_cs[_qp]); 67 392880 : ADReal beta_0 = _delta[_qp] * M; 68 785760 : ADReal beta_1 = -gamma_1 * M * (_sigma_cs[_qp] - _sigma_ts[_qp]) + gamma_0; 69 1571520 : ADReal beta_2 = std::sqrt(3.0) * (-gamma_1 * M * (_sigma_cs[_qp] + _sigma_ts[_qp]) + gamma_0); 70 392880 : ADReal beta_3 = _L[_qp] * _sigma_ts[_qp] / _mu[_qp] / K / _Gc[_qp]; 71 : 72 : // Compute the external driving force required to recover the desired strength envelope. 73 392880 : _ex_driving[_qp] = 74 392880 : (beta_2 * std::sqrt(J2) + beta_1 * I1 + beta_0) + 75 1178640 : (1.0 - std::sqrt(I1 * I1) / I1) / std::pow(_g[_qp], 1.5) * 76 1964400 : (J2 / 2.0 / _mu[_qp] + I1 * I1 / 6.0 / (3.0 * _lambda[_qp] + 2.0 * _mu[_qp])); 77 : 78 1178640 : _stress_balance[_qp] = J2 / _mu[_qp] + std::pow(I1, 2) / 9.0 / K - _ex_driving[_qp] - M; 79 : 80 392880 : _druck_prager_balance[_qp] = 81 392880 : std::sqrt(J2) + 82 785760 : (_sigma_cs[_qp] - _sigma_ts[_qp]) / std::sqrt(3.0) / (_sigma_cs[_qp] + _sigma_ts[_qp]) * I1 - 83 1178640 : 2.0 * _sigma_ts[_qp] * _sigma_cs[_qp] / std::sqrt(3.0) / (_sigma_cs[_qp] + _sigma_ts[_qp]); 84 392880 : }