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