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 "HenckyIsotropicElasticity.h" 6 : #include "RaccoonUtils.h" 7 : 8 : registerMooseObject("raccoonApp", HenckyIsotropicElasticity); 9 : 10 : InputParameters 11 198 : HenckyIsotropicElasticity::validParams() 12 : { 13 198 : InputParameters params = LargeDeformationElasticityModel::validParams(); 14 198 : params.addClassDescription("Isotropic Hencky-type hyperelasticity model. The logarithmic right " 15 : "Cauchy-Green strain tensor is used as the strain measure."); 16 : 17 396 : params.addRequiredParam<MaterialPropertyName>("bulk_modulus", "The bulk modulus $K$"); 18 396 : params.addRequiredParam<MaterialPropertyName>("shear_modulus", "The shear modulus $G$"); 19 : 20 396 : params.addRequiredCoupledVar("phase_field", "Name of the phase-field (damage) variable"); 21 396 : params.addParam<MaterialPropertyName>( 22 : "strain_energy_density", 23 : "psie", 24 : "Name of the strain energy density computed by this material model"); 25 396 : params.addParam<MaterialPropertyName>("degradation_function", "g", "The degradation function"); 26 396 : params.addParam<MooseEnum>( 27 594 : "decomposition", MooseEnum("NONE SPECTRAL VOLDEV", "NONE"), "The decomposition method"); 28 : 29 198 : return params; 30 0 : } 31 : 32 45 : HenckyIsotropicElasticity::HenckyIsotropicElasticity(const InputParameters & parameters) 33 : : LargeDeformationElasticityModel(parameters), 34 : DerivativeMaterialPropertyNameInterface(), 35 135 : _K(getADMaterialPropertyByName<Real>(prependBaseName("bulk_modulus", true))), 36 90 : _G(getADMaterialPropertyByName<Real>(prependBaseName("shear_modulus", true))), 37 : 38 90 : _d_name(getVar("phase_field", 0)->name()), 39 : 40 : // The strain energy density and its derivatives 41 45 : _psie_name(prependBaseName("strain_energy_density", true)), 42 45 : _psie(declareADProperty<Real>(_psie_name)), 43 45 : _psie_active(declareADProperty<Real>(_psie_name + "_active")), 44 90 : _dpsie_dd(declareADProperty<Real>(derivativePropertyName(_psie_name, {_d_name}))), 45 : 46 : // The degradation function and its derivatives 47 45 : _g_name(prependBaseName("degradation_function", true)), 48 45 : _g(getADMaterialProperty<Real>(_g_name)), 49 135 : _dg_dd(getADMaterialProperty<Real>(derivativePropertyName(_g_name, {_d_name}))), 50 : 51 135 : _decomposition(getParam<MooseEnum>("decomposition").getEnum<Decomposition>()) 52 : { 53 135 : } 54 : 55 : ADRankTwoTensor 56 1648109 : HenckyIsotropicElasticity::computeMandelStress(const ADRankTwoTensor & Fe, 57 : const bool plasticity_update) 58 : { 59 1648109 : ADRankTwoTensor stress; 60 : 61 1648109 : if (_decomposition == Decomposition::none) 62 3296218 : stress = computeMandelStressNoDecomposition(Fe, plasticity_update); 63 0 : else if (_decomposition == Decomposition::voldev) 64 0 : stress = computeMandelStressVolDevDecomposition(Fe, plasticity_update); 65 0 : else if (_decomposition == Decomposition::spectral) 66 0 : stress = computeMandelStressSpectralDecomposition(Fe, plasticity_update); 67 : else 68 0 : paramError("decomposition", "Unsupported decomposition type."); 69 : 70 1648109 : return stress; 71 : } 72 : 73 : ADRankTwoTensor 74 1648109 : HenckyIsotropicElasticity::computeMandelStressNoDecomposition(const ADRankTwoTensor & Fe, 75 : const bool plasticity_update) 76 : { 77 : ADRankTwoTensor strain = Fe; 78 : // If this is called during a plasticity update, we need to first exponentiate Fe, where Fe should 79 : // be some plastic flow. The following operations cancel out with an exponentiation of Fe, so we 80 : // only do this in the case of exponentiate == false 81 1648109 : if (!plasticity_update) 82 1038128 : strain = 0.5 * RaccoonUtils::log(Fe.transpose() * Fe); 83 : 84 1648109 : const ADRankTwoTensor I2(ADRankTwoTensor::initIdentity); 85 : // Here, we keep the volumetric part no matter what. But ideally, in the case of J2 plasticity, 86 : // the volumetric part of the flow should be zero, and we could save some operations. 87 6592436 : ADRankTwoTensor stress_intact = _K[_qp] * strain.trace() * I2 + 2 * _G[_qp] * strain.deviatoric(); 88 1648109 : ADRankTwoTensor stress = _g[_qp] * stress_intact; 89 : 90 : // If plasticity_update == false, then we are not in the middle of a plasticity update, hence we 91 : // compute the strain energy density 92 1648109 : if (!plasticity_update) 93 : { 94 : // It is convenient that the Mandel stress is conjugate to the log strain 95 1038128 : _psie_active[_qp] = 0.5 * stress_intact.doubleContraction(strain); 96 1038128 : _psie[_qp] = _g[_qp] * _psie_active[_qp]; 97 1038128 : _dpsie_dd[_qp] = _dg_dd[_qp] * _psie_active[_qp]; 98 : } 99 : 100 1648109 : return stress; 101 : } 102 : 103 : ADRankTwoTensor 104 0 : HenckyIsotropicElasticity::computeMandelStressVolDevDecomposition(const ADRankTwoTensor & Fe, 105 : const bool plasticity_update) 106 : { 107 : ADRankTwoTensor strain = Fe; 108 : // If this is called during a plasticity update, we need to first exponentiate Fe, where Fe should 109 : // be some plastic flow. The following operations cancel out with an exponentiation of Fe, so we 110 : // only do this in the case of plasticity_update == false 111 0 : if (!plasticity_update) 112 0 : strain = 0.5 * RaccoonUtils::log(Fe.transpose() * Fe); 113 : 114 0 : const ADRankTwoTensor I2(ADRankTwoTensor::initIdentity); 115 : // Here, we keep the volumetric part no matter what. But ideally, in the case of J2 plasticity, 116 : // the volumetric part of the flow should be zero, and we could save some operations. 117 0 : ADReal strain_tr = strain.trace(); 118 0 : ADReal strain_tr_pos = RaccoonUtils::Macaulay(strain_tr); 119 : ADReal strain_tr_neg = strain_tr - strain_tr_pos; 120 0 : ADRankTwoTensor strain_dev = strain.deviatoric(); 121 : 122 0 : ADRankTwoTensor stress_intact = _K[_qp] * strain_tr * I2 + 2 * _G[_qp] * strain_dev; 123 0 : ADRankTwoTensor stress_neg = _K[_qp] * strain_tr_neg * I2; 124 0 : ADRankTwoTensor stress_pos = stress_intact - stress_neg; 125 0 : ADRankTwoTensor stress = _g[_qp] * stress_pos + stress_neg; 126 : 127 : // If plasticity_update == false, then we are not in the middle of a plasticity update, hence we 128 : // compute the strain energy density 129 0 : if (!plasticity_update) 130 : { 131 : ADReal psie_intact = 132 0 : 0.5 * _K[_qp] * strain_tr * strain_tr + _G[_qp] * strain_dev.doubleContraction(strain_dev); 133 0 : ADReal psie_inactive = 0.5 * _K[_qp] * strain_tr_neg * strain_tr_neg; 134 0 : _psie_active[_qp] = psie_intact - psie_inactive; 135 0 : _psie[_qp] = _g[_qp] * _psie_active[_qp] + psie_inactive; 136 0 : _dpsie_dd[_qp] = _dg_dd[_qp] * _psie_active[_qp]; 137 : } 138 : 139 0 : return stress; 140 : } 141 : 142 : ADRankTwoTensor 143 0 : HenckyIsotropicElasticity::computeMandelStressSpectralDecomposition(const ADRankTwoTensor & Fe, 144 : const bool plasticity_update) 145 : { 146 : ADRankTwoTensor strain = Fe; 147 : // If this is called during a plasticity update, we need to first exponentiate Fe, where Fe should 148 : // be some plastic flow. The following operations cancel out with an exponentiation of Fe, so we 149 : // only do this in the case of plasticity_update == false 150 0 : if (!plasticity_update) 151 0 : strain = 0.5 * RaccoonUtils::log(Fe.transpose() * Fe); 152 : 153 0 : const ADReal lambda = _K[_qp] - 2 * _G[_qp] / LIBMESH_DIM; 154 0 : const ADRankTwoTensor I2(ADRankTwoTensor::initIdentity); 155 0 : ADReal strain_tr = strain.trace(); 156 0 : ADReal strain_tr_pos = RaccoonUtils::Macaulay(strain_tr); 157 : 158 : // Spectral decomposition 159 0 : ADRankTwoTensor strain_pos = RaccoonUtils::spectralDecomposition(strain); 160 : 161 : // Stress 162 0 : ADRankTwoTensor stress_intact = _K[_qp] * strain.trace() * I2 + 2 * _G[_qp] * strain.deviatoric(); 163 0 : ADRankTwoTensor stress_pos = lambda * strain_tr_pos * I2 + 2 * _G[_qp] * strain_pos; 164 0 : ADRankTwoTensor stress_neg = stress_intact - stress_pos; 165 0 : ADRankTwoTensor stress = _g[_qp] * stress_pos + stress_neg; 166 : // If plasticity_update == false, then we are not in the middle of a plasticity update, hence we 167 : // compute the strain energy density 168 0 : if (!plasticity_update) 169 : { 170 : ADReal psie_intact = 171 0 : 0.5 * lambda * strain_tr * strain_tr + _G[_qp] * strain.doubleContraction(strain); 172 0 : _psie_active[_qp] = 0.5 * lambda * strain_tr_pos * strain_tr_pos + 173 0 : _G[_qp] * strain_pos.doubleContraction(strain_pos); 174 : ADReal psie_inactive = psie_intact - _psie_active[_qp]; 175 0 : _psie[_qp] = _g[_qp] * _psie_active[_qp] + psie_inactive; 176 0 : _dpsie_dd[_qp] = _dg_dd[_qp] * _psie_active[_qp]; 177 : } 178 0 : return stress; 179 : }