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 "CNHIsotropicElasticity.h" 6 : #include "RaccoonUtils.h" 7 : 8 : registerMooseObject("raccoonApp", CNHIsotropicElasticity); 9 : 10 : InputParameters 11 142 : CNHIsotropicElasticity::validParams() 12 : { 13 142 : InputParameters params = LargeDeformationElasticityModel::validParams(); 14 142 : params.addClassDescription("Isotropic Compressible Neo-Hookean hyperelasticity model."); 15 : 16 284 : params.addRequiredParam<MaterialPropertyName>("bulk_modulus", "The bulk modulus $K$"); 17 284 : params.addRequiredParam<MaterialPropertyName>("shear_modulus", "The shear modulus $G$"); 18 : 19 284 : params.addRequiredCoupledVar("phase_field", "Name of the phase-field (damage) variable"); 20 284 : params.addParam<MaterialPropertyName>( 21 : "strain_energy_density", 22 : "psie", 23 : "Name of the strain energy density computed by this material model"); 24 284 : params.addParam<MaterialPropertyName>("degradation_function", "g", "The degradation function"); 25 284 : params.addParam<MooseEnum>( 26 426 : "decomposition", MooseEnum("NONE SPECTRAL VOLDEV", "NONE"), "The decomposition method"); 27 : 28 142 : return params; 29 0 : } 30 : 31 3 : CNHIsotropicElasticity::CNHIsotropicElasticity(const InputParameters & parameters) 32 : : LargeDeformationElasticityModel(parameters), 33 : DerivativeMaterialPropertyNameInterface(), 34 9 : _K(getADMaterialPropertyByName<Real>(prependBaseName("bulk_modulus", true))), 35 6 : _G(getADMaterialPropertyByName<Real>(prependBaseName("shear_modulus", true))), 36 : 37 6 : _d_name(getVar("phase_field", 0)->name()), 38 : 39 : // The strain energy density and its derivatives 40 3 : _psie_name(prependBaseName("strain_energy_density", true)), 41 3 : _psie(declareADProperty<Real>(_psie_name)), 42 3 : _psie_active(declareADProperty<Real>(_psie_name + "_active")), 43 6 : _dpsie_dd(declareADProperty<Real>(derivativePropertyName(_psie_name, {_d_name}))), 44 : 45 : // The degradation function and its derivatives 46 3 : _g_name(prependBaseName("degradation_function", true)), 47 3 : _g(getADMaterialProperty<Real>(_g_name)), 48 9 : _dg_dd(getADMaterialProperty<Real>(derivativePropertyName(_g_name, {_d_name}))), 49 : 50 9 : _decomposition(getParam<MooseEnum>("decomposition").getEnum<Decomposition>()) 51 : { 52 9 : } 53 : 54 : ADRankTwoTensor 55 3144 : CNHIsotropicElasticity::computeMandelStress(const ADRankTwoTensor & Fe, 56 : const bool plasticity_update) 57 : { 58 3144 : ADRankTwoTensor stress; 59 : 60 3144 : if (_decomposition == Decomposition::none) 61 6288 : stress = computeMandelStressNoDecomposition(Fe, plasticity_update); 62 0 : else if (_decomposition == Decomposition::voldev) 63 0 : stress = computeMandelStressVolDevDecomposition(Fe, plasticity_update); 64 : else 65 0 : paramError("decomposition", "Unsupported decomposition type."); 66 : 67 3144 : return stress; 68 : } 69 : 70 : ADRankTwoTensor 71 3144 : CNHIsotropicElasticity::computeMandelStressNoDecomposition(const ADRankTwoTensor & Fe, 72 : const bool plasticity_update) 73 : { 74 : // We use the left Cauchy-Green strain 75 3144 : ADRankTwoTensor strain; 76 3144 : if (plasticity_update) 77 : { 78 1048 : ADRankTwoTensor expFe = RaccoonUtils::exp(Fe); 79 1048 : strain = expFe * expFe.transpose(); 80 : } 81 : else 82 4192 : strain = Fe * Fe.transpose(); 83 : 84 3144 : ADReal J = std::sqrt(strain.det()); 85 : 86 3144 : const ADRankTwoTensor I2(ADRankTwoTensor::initIdentity); 87 : // Here, we keep the volumetric part no matter what. But ideally, in the case of J2 plasticity, 88 : // the volumetric part of the flow should be zero, and we could save some operations. 89 12576 : ADRankTwoTensor stress_intact = 0.5 * _K[_qp] * (J * J - 1) * I2 + _G[_qp] * strain.deviatoric(); 90 3144 : ADRankTwoTensor stress = _g[_qp] * stress_intact; 91 : 92 : // If plasticity_update == false, then we are not in the middle of a plasticity update, hence we 93 : // compute the strain energy density 94 3144 : if (!plasticity_update) 95 : { 96 2096 : ADRankTwoTensor strain_bar = std::pow(J, -2. / 3.) * strain; 97 6288 : ADReal U = 0.5 * _K[_qp] * (0.5 * (J * J - 1) - std::log(J)); 98 6288 : ADReal W = 0.5 * _G[_qp] * (strain_bar.trace() - 3.0); 99 2096 : _psie_active[_qp] = U + W; 100 4192 : _psie[_qp] = _g[_qp] * _psie_active[_qp]; 101 4192 : _dpsie_dd[_qp] = _dg_dd[_qp] * _psie_active[_qp]; 102 : } 103 : 104 3144 : return stress; 105 : } 106 : 107 : ADRankTwoTensor 108 0 : CNHIsotropicElasticity::computeMandelStressVolDevDecomposition(const ADRankTwoTensor & Fe, 109 : const bool plasticity_update) 110 : { 111 : // We use the left Cauchy-Green strain 112 0 : ADRankTwoTensor strain; 113 0 : if (plasticity_update) 114 : { 115 0 : ADRankTwoTensor expFe = RaccoonUtils::exp(Fe); 116 0 : strain = expFe * expFe.transpose(); 117 : } 118 : else 119 0 : strain = Fe * Fe.transpose(); 120 : 121 0 : ADReal J = std::sqrt(strain.det()); 122 : 123 0 : const ADRankTwoTensor I2(ADRankTwoTensor::initIdentity); 124 : // Here, we keep the volumetric part no matter what. But ideally, in the case of J2 plasticity, 125 : // the volumetric part of the flow should be zero, and we could save some operations. 126 0 : ADRankTwoTensor stress_vol = 0.5 * _K[_qp] * (J * J - 1) * I2; 127 0 : ADRankTwoTensor stress_dev = _G[_qp] * strain.deviatoric(); 128 : ADRankTwoTensor stress = 129 0 : J > 1 ? _g[_qp] * (stress_vol + stress_dev) : stress_vol + _g[_qp] * stress_dev; 130 : 131 : // If plasticity_update == false, then we are not in the middle of a plasticity update, hence we 132 : // compute the strain energy density 133 0 : if (!plasticity_update) 134 : { 135 0 : ADRankTwoTensor strain_bar = std::pow(J, -2. / 3.) * strain; 136 0 : ADReal U = 0.5 * _K[_qp] * (0.5 * (J * J - 1) - std::log(J)); 137 0 : ADReal W = 0.5 * _G[_qp] * (strain_bar.trace() - 3.0); 138 0 : _psie_active[_qp] = J > 1 ? U + W : W; 139 0 : _psie[_qp] = _g[_qp] * _psie_active[_qp]; 140 0 : _dpsie_dd[_qp] = _dg_dd[_qp] * _psie_active[_qp]; 141 : } 142 : 143 0 : return stress; 144 : }