Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "HillPlasticityStressUpdate.h"
11 : #include "ElasticityTensorTools.h"
12 :
13 : registerMooseObject("SolidMechanicsApp", ADHillPlasticityStressUpdate);
14 : registerMooseObject("SolidMechanicsApp", HillPlasticityStressUpdate);
15 :
16 : template <bool is_ad>
17 : InputParameters
18 28 : HillPlasticityStressUpdateTempl<is_ad>::validParams()
19 : {
20 28 : InputParameters params = AnisotropicReturnPlasticityStressUpdateBaseTempl<is_ad>::validParams();
21 28 : params.addClassDescription(
22 : "This class uses the generalized radial return for anisotropic plasticity model."
23 : "This class can be used in conjunction with other creep and plasticity materials for "
24 : "more complex simulations.");
25 :
26 56 : params.addRequiredParam<Real>("hardening_constant",
27 : "Hardening constant (H) for anisotropic plasticity");
28 56 : params.addParam<Real>(
29 56 : "hardening_exponent", 1.0, "Hardening exponent (n) for anisotropic plasticity");
30 56 : params.addRequiredParam<Real>("yield_stress",
31 : "Yield stress (constant value) for anisotropic plasticity");
32 :
33 28 : return params;
34 0 : }
35 :
36 : template <bool is_ad>
37 21 : HillPlasticityStressUpdateTempl<is_ad>::HillPlasticityStressUpdateTempl(
38 : const InputParameters & parameters)
39 : : AnisotropicReturnPlasticityStressUpdateBaseTempl<is_ad>(parameters),
40 42 : _qsigma(0.0),
41 21 : _eigenvalues_hill(6),
42 21 : _eigenvectors_hill(6, 6),
43 42 : _hardening_constant(this->template getParam<Real>("hardening_constant")),
44 42 : _hardening_exponent(this->template getParam<Real>("hardening_exponent")),
45 42 : _hardening_variable(this->template declareGenericProperty<Real, is_ad>(this->_base_name +
46 : "hardening_variable")),
47 21 : _hardening_variable_old(
48 21 : this->template getMaterialPropertyOld<Real>(this->_base_name + "hardening_variable")),
49 21 : _hardening_derivative(0.0),
50 21 : _yield_condition(1.0),
51 42 : _yield_stress(this->template getParam<Real>("yield_stress")),
52 42 : _hill_tensor(this->template getMaterialPropertyByName<DenseMatrix<Real>>(this->_base_name +
53 : "hill_tensor")),
54 42 : _stress_np1(6)
55 : {
56 21 : }
57 :
58 : template <bool is_ad>
59 : void
60 0 : HillPlasticityStressUpdateTempl<is_ad>::propagateQpStatefulProperties()
61 : {
62 0 : _hardening_variable[_qp] = _hardening_variable_old[_qp];
63 0 : _plasticity_strain[_qp] = _plasticity_strain_old[_qp];
64 0 : AnisotropicReturnPlasticityStressUpdateBaseTempl<is_ad>::propagateQpStatefulProperties();
65 0 : }
66 :
67 : template <bool is_ad>
68 : void
69 11320 : HillPlasticityStressUpdateTempl<is_ad>::computeStressInitialize(
70 : const GenericDenseVector<is_ad> & stress_dev,
71 : const GenericDenseVector<is_ad> & /*stress*/,
72 : const GenericRankFourTensor<is_ad> & elasticity_tensor)
73 : {
74 11320 : _hardening_variable[_qp] = _hardening_variable_old[_qp];
75 11320 : _plasticity_strain[_qp] = _plasticity_strain_old[_qp];
76 11320 : _effective_inelastic_strain[_qp] = _effective_inelastic_strain_old[_qp];
77 :
78 22640 : _two_shear_modulus = 2.0 * ElasticityTensorTools::getIsotropicShearModulus(elasticity_tensor);
79 :
80 : // Hill constants: We use directly the transformation tensor, which won't be updated if not
81 : // necessary in the Hill tensor material.
82 11320 : computeHillTensorEigenDecomposition(_hill_tensor[_qp]);
83 :
84 11320 : _yield_condition = 1.0; // Some positive value
85 11320 : _yield_condition = -computeResidual(stress_dev, stress_dev, 0.0);
86 11320 : }
87 :
88 : template <bool is_ad>
89 : GenericReal<is_ad>
90 11320 : HillPlasticityStressUpdateTempl<is_ad>::computeOmega(const GenericReal<is_ad> & delta_gamma,
91 : const GenericDenseVector<is_ad> & stress_trial)
92 : {
93 11320 : GenericDenseVector<is_ad> K(6);
94 11320 : GenericReal<is_ad> omega = 0.0;
95 :
96 79240 : for (unsigned int i = 0; i < 6; i++)
97 : {
98 67920 : K(i) = _eigenvalues_hill(i) /
99 203760 : (Utility::pow<2>(1 + _two_shear_modulus * delta_gamma * _eigenvalues_hill(i)));
100 67920 : omega += K(i) * stress_trial(i) * stress_trial(i);
101 : }
102 11320 : omega *= 0.5;
103 :
104 11320 : if (omega == 0.0)
105 160 : omega = 1.0e-36;
106 :
107 : using std::sqrt;
108 22640 : return sqrt(omega);
109 : }
110 :
111 : template <bool is_ad>
112 : void
113 0 : HillPlasticityStressUpdateTempl<is_ad>::computeDeltaDerivatives(
114 : const GenericReal<is_ad> & delta_gamma,
115 : const GenericDenseVector<is_ad> & stress_trial,
116 : const GenericReal<is_ad> & sy_alpha,
117 : GenericReal<is_ad> & omega,
118 : GenericReal<is_ad> & omega_gamma,
119 : GenericReal<is_ad> & sy_gamma)
120 : {
121 0 : omega_gamma = 0.0;
122 0 : sy_gamma = 0.0;
123 :
124 0 : GenericDenseVector<is_ad> K_deltaGamma(6);
125 0 : omega = computeOmega(delta_gamma, stress_trial);
126 :
127 0 : GenericDenseVector<is_ad> K(6);
128 0 : for (unsigned int i = 0; i < 6; i++)
129 0 : K(i) = _eigenvalues_hill(i) /
130 0 : (Utility::pow<2>(1 + _two_shear_modulus * delta_gamma * _eigenvalues_hill(i)));
131 :
132 0 : for (unsigned int i = 0; i < 6; i++)
133 0 : K_deltaGamma(i) = -2.0 * _two_shear_modulus * _eigenvalues_hill(i) * K(i) /
134 0 : (1 + _two_shear_modulus * delta_gamma * _eigenvalues_hill(i));
135 :
136 0 : for (unsigned int i = 0; i < 6; i++)
137 0 : omega_gamma += K_deltaGamma(i) * stress_trial(i) * stress_trial(i);
138 :
139 0 : omega_gamma /= 4.0 * omega;
140 0 : sy_gamma = 2.0 * sy_alpha * (omega + delta_gamma * omega_gamma);
141 0 : }
142 :
143 : template <bool is_ad>
144 : Real
145 11160 : HillPlasticityStressUpdateTempl<is_ad>::computeReferenceResidual(
146 : const GenericDenseVector<is_ad> & /*effective_trial_stress*/,
147 : const GenericDenseVector<is_ad> & /*stress_new*/,
148 : const GenericReal<is_ad> & /*residual*/,
149 : const GenericReal<is_ad> & /*scalar_effective_inelastic_strain*/)
150 : {
151 : // Equation is already normalized.
152 11160 : return 1.0;
153 : }
154 :
155 : template <bool is_ad>
156 : GenericReal<is_ad>
157 22480 : HillPlasticityStressUpdateTempl<is_ad>::computeResidual(
158 : const GenericDenseVector<is_ad> & stress_dev,
159 : const GenericDenseVector<is_ad> & /*stress_sigma*/,
160 : const GenericReal<is_ad> & delta_gamma)
161 : {
162 : using std::pow;
163 :
164 : // If in elastic regime, just return
165 22480 : if (_yield_condition <= 0.0)
166 11160 : return 0.0;
167 :
168 11320 : GenericDenseMatrix<is_ad> eigenvectors_hill_transpose(6, 6);
169 :
170 11320 : _eigenvectors_hill.get_transpose(eigenvectors_hill_transpose);
171 11320 : eigenvectors_hill_transpose.vector_mult(_stress_np1, stress_dev);
172 :
173 11320 : GenericReal<is_ad> omega = computeOmega(delta_gamma, _stress_np1);
174 :
175 : // Hardening variable is \alpha isotropic hardening for now.
176 11320 : _hardening_variable[_qp] = computeHardeningValue(delta_gamma, omega);
177 11320 : GenericReal<is_ad> s_y =
178 11320 : _hardening_constant * pow(_hardening_variable[_qp] + 1.0e-30, _hardening_exponent) +
179 11320 : _yield_stress;
180 :
181 11320 : GenericReal<is_ad> residual = 0.0;
182 11320 : residual = s_y / omega - 1.0;
183 :
184 11320 : return residual;
185 11320 : }
186 :
187 : template <bool is_ad>
188 : GenericReal<is_ad>
189 0 : HillPlasticityStressUpdateTempl<is_ad>::computeDerivative(
190 : const GenericDenseVector<is_ad> & /*stress_dev*/,
191 : const GenericDenseVector<is_ad> & /*stress_sigma*/,
192 : const GenericReal<is_ad> & delta_gamma)
193 : {
194 : // If in elastic regime, return unit derivative
195 0 : if (_yield_condition <= 0.0)
196 0 : return 1.0;
197 :
198 0 : GenericReal<is_ad> omega = computeOmega(delta_gamma, _stress_np1);
199 0 : _hardening_derivative = computeHardeningDerivative();
200 :
201 0 : GenericReal<is_ad> sy =
202 0 : _hardening_derivative * computeHardeningValue(delta_gamma, omega) + _yield_stress;
203 0 : GenericReal<is_ad> sy_alpha = _hardening_derivative;
204 :
205 : GenericReal<is_ad> omega_gamma;
206 : GenericReal<is_ad> sy_gamma;
207 :
208 0 : computeDeltaDerivatives(delta_gamma, _stress_np1, sy_alpha, omega, omega_gamma, sy_gamma);
209 0 : GenericReal<is_ad> residual_derivative = 1 / omega * (sy_gamma - 1 / omega * omega_gamma * sy);
210 :
211 0 : return residual_derivative;
212 : }
213 :
214 : template <bool is_ad>
215 : void
216 11320 : HillPlasticityStressUpdateTempl<is_ad>::computeHillTensorEigenDecomposition(
217 : const DenseMatrix<Real> & hill_tensor)
218 : {
219 : const unsigned int dimension = hill_tensor.n();
220 :
221 : AnisotropyMatrixReal A;
222 79240 : for (unsigned int index_i = 0; index_i < dimension; index_i++)
223 475440 : for (unsigned int index_j = 0; index_j < dimension; index_j++)
224 407520 : A(index_i, index_j) = MetaPhysicL::raw_value(hill_tensor(index_i, index_j));
225 :
226 11320 : if (isBlockDiagonal(A))
227 : {
228 160 : Eigen::SelfAdjointEigenSolver<AnisotropyMatrixRealBlock> es(A.block<3, 3>(0, 0));
229 :
230 : auto lambda = es.eigenvalues();
231 : auto v = es.eigenvectors();
232 :
233 160 : _eigenvalues_hill(0) = lambda(0);
234 160 : _eigenvalues_hill(1) = lambda(1);
235 160 : _eigenvalues_hill(2) = lambda(2);
236 160 : _eigenvalues_hill(3) = A(3, 3);
237 160 : _eigenvalues_hill(4) = A(4, 4);
238 160 : _eigenvalues_hill(5) = A(5, 5);
239 :
240 160 : _eigenvectors_hill(0, 0) = v(0, 0);
241 160 : _eigenvectors_hill(0, 1) = v(0, 1);
242 160 : _eigenvectors_hill(0, 2) = v(0, 2);
243 160 : _eigenvectors_hill(1, 0) = v(1, 0);
244 160 : _eigenvectors_hill(1, 1) = v(1, 1);
245 160 : _eigenvectors_hill(1, 2) = v(1, 2);
246 160 : _eigenvectors_hill(2, 0) = v(2, 0);
247 160 : _eigenvectors_hill(2, 1) = v(2, 1);
248 160 : _eigenvectors_hill(2, 2) = v(2, 2);
249 160 : _eigenvectors_hill(3, 3) = 1.0;
250 160 : _eigenvectors_hill(4, 4) = 1.0;
251 160 : _eigenvectors_hill(5, 5) = 1.0;
252 : }
253 : else
254 : {
255 : Eigen::SelfAdjointEigenSolver<AnisotropyMatrixReal> es_b(A);
256 :
257 : auto lambda_b = es_b.eigenvalues();
258 : auto v_b = es_b.eigenvectors();
259 78120 : for (unsigned int index_i = 0; index_i < dimension; index_i++)
260 66960 : _eigenvalues_hill(index_i) = lambda_b(index_i);
261 :
262 78120 : for (unsigned int index_i = 0; index_i < dimension; index_i++)
263 468720 : for (unsigned int index_j = 0; index_j < dimension; index_j++)
264 401760 : _eigenvectors_hill(index_i, index_j) = v_b(index_i, index_j);
265 : }
266 11320 : }
267 :
268 : template <bool is_ad>
269 : GenericReal<is_ad>
270 11320 : HillPlasticityStressUpdateTempl<is_ad>::computeHardeningValue(
271 : const GenericReal<is_ad> & delta_gamma, const GenericReal<is_ad> & omega)
272 : {
273 22640 : return _hardening_variable_old[_qp] + 2.0 * delta_gamma * omega;
274 : }
275 :
276 : template <bool is_ad>
277 : Real
278 0 : HillPlasticityStressUpdateTempl<is_ad>::computeHardeningDerivative()
279 : {
280 : using std::pow;
281 :
282 0 : return _hardening_constant * _hardening_exponent *
283 0 : MetaPhysicL::raw_value(pow(_hardening_variable[_qp] + 1.0e-30, _hardening_exponent - 1));
284 : }
285 :
286 : template <bool is_ad>
287 : void
288 11160 : HillPlasticityStressUpdateTempl<is_ad>::computeStrainFinalize(
289 : GenericRankTwoTensor<is_ad> & inelasticStrainIncrement,
290 : const GenericRankTwoTensor<is_ad> & stress,
291 : const GenericDenseVector<is_ad> & stress_dev,
292 : const GenericReal<is_ad> & delta_gamma)
293 : {
294 : using std::sqrt;
295 :
296 : // e^P = delta_gamma * hill_tensor * stress
297 11160 : GenericDenseVector<is_ad> inelasticStrainIncrement_vector(6);
298 11160 : GenericDenseVector<is_ad> hill_stress(6);
299 11160 : _hill_tensor[_qp].vector_mult(hill_stress, stress_dev);
300 11160 : hill_stress.scale(delta_gamma);
301 : inelasticStrainIncrement_vector = hill_stress;
302 :
303 11160 : inelasticStrainIncrement(0, 0) = inelasticStrainIncrement_vector(0);
304 11160 : inelasticStrainIncrement(1, 1) = inelasticStrainIncrement_vector(1);
305 11160 : inelasticStrainIncrement(2, 2) = inelasticStrainIncrement_vector(2);
306 11160 : inelasticStrainIncrement(0, 1) = inelasticStrainIncrement(1, 0) =
307 11160 : inelasticStrainIncrement_vector(3) / 2.0;
308 11160 : inelasticStrainIncrement(1, 2) = inelasticStrainIncrement(2, 1) =
309 11160 : inelasticStrainIncrement_vector(4) / 2.0;
310 11160 : inelasticStrainIncrement(0, 2) = inelasticStrainIncrement(2, 0) =
311 11160 : inelasticStrainIncrement_vector(5) / 2.0;
312 :
313 : // Calculate equivalent plastic strain
314 11160 : GenericDenseVector<is_ad> Mepsilon(6);
315 11160 : _hill_tensor[_qp].vector_mult(Mepsilon, inelasticStrainIncrement_vector);
316 11160 : GenericReal<is_ad> eq_plastic_strain_inc = Mepsilon.dot(inelasticStrainIncrement_vector);
317 :
318 11160 : if (eq_plastic_strain_inc > 0.0)
319 0 : eq_plastic_strain_inc = sqrt(eq_plastic_strain_inc);
320 :
321 22320 : _effective_inelastic_strain[_qp] = _effective_inelastic_strain_old[_qp] + eq_plastic_strain_inc;
322 :
323 11160 : AnisotropicReturnPlasticityStressUpdateBaseTempl<is_ad>::computeStrainFinalize(
324 : inelasticStrainIncrement, stress, stress_dev, delta_gamma);
325 11160 : }
326 :
327 : template <bool is_ad>
328 : void
329 11320 : HillPlasticityStressUpdateTempl<is_ad>::computeStressFinalize(
330 : const GenericRankTwoTensor<is_ad> & /*plastic_strain_increment*/,
331 : const GenericReal<is_ad> & delta_gamma,
332 : GenericRankTwoTensor<is_ad> & stress_new,
333 : const GenericDenseVector<is_ad> & stress_dev,
334 : const GenericRankTwoTensor<is_ad> & /*sstress_old*/,
335 : const GenericRankFourTensor<is_ad> & /*elasticity_tensor*/)
336 : {
337 : // Need to compute this iteration's stress tensor based on the scalar variable for deviatoric
338 : // s(n+1) = {Q [I + 2*nu*delta_gamma*Delta]^(-1) Q^T} s(trial)
339 :
340 11320 : if (_yield_condition <= 0.0)
341 11320 : return;
342 :
343 0 : GenericDenseMatrix<is_ad> inv_matrix(6, 6);
344 0 : for (unsigned int i = 0; i < 6; i++)
345 0 : inv_matrix(i, i) = 1 / (1 + _two_shear_modulus * delta_gamma * _eigenvalues_hill(i));
346 :
347 0 : GenericDenseMatrix<is_ad> eigenvectors_hill_transpose(6, 6);
348 :
349 0 : _eigenvectors_hill.get_transpose(eigenvectors_hill_transpose);
350 0 : GenericDenseMatrix<is_ad> eigenvectors_hill_copy(_eigenvectors_hill);
351 :
352 : // Right multiply by matrix of eigenvectors transpose
353 0 : inv_matrix.right_multiply(eigenvectors_hill_transpose);
354 : // Right multiply eigenvector matrix by [I + 2*nu*delta_gamma*Delta]^(-1) Q^T
355 0 : eigenvectors_hill_copy.right_multiply(inv_matrix);
356 :
357 0 : GenericDenseVector<is_ad> stress_np1(6);
358 0 : eigenvectors_hill_copy.vector_mult(stress_np1, stress_dev);
359 :
360 0 : GenericRankTwoTensor<is_ad> stress_new_volumetric = stress_new - stress_new.deviatoric();
361 :
362 0 : stress_new(0, 0) = stress_new_volumetric(0, 0) + stress_np1(0);
363 0 : stress_new(1, 1) = stress_new_volumetric(1, 1) + stress_np1(1);
364 0 : stress_new(2, 2) = stress_new_volumetric(2, 2) + stress_np1(2);
365 0 : stress_new(0, 1) = stress_new(1, 0) = stress_np1(3);
366 0 : stress_new(1, 2) = stress_new(2, 1) = stress_np1(4);
367 0 : stress_new(0, 2) = stress_new(2, 0) = stress_np1(5);
368 :
369 0 : GenericReal<is_ad> omega = computeOmega(delta_gamma, _stress_np1);
370 0 : _hardening_variable[_qp] = computeHardeningValue(delta_gamma, omega);
371 0 : }
372 :
373 : template class HillPlasticityStressUpdateTempl<false>;
374 : template class HillPlasticityStressUpdateTempl<true>;
|
