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 "HillElastoPlasticityStressUpdate.h"
11 : #include "RankTwoScalarTools.h"
12 : #include "ElasticityTensorTools.h"
13 :
14 : registerMooseObject("SolidMechanicsApp", ADHillElastoPlasticityStressUpdate);
15 : registerMooseObject("SolidMechanicsApp", HillElastoPlasticityStressUpdate);
16 :
17 : template <bool is_ad>
18 : InputParameters
19 96 : HillElastoPlasticityStressUpdateTempl<is_ad>::validParams()
20 : {
21 96 : InputParameters params = AnisotropicReturnPlasticityStressUpdateBaseTempl<is_ad>::validParams();
22 96 : params.addClassDescription(
23 : "This class uses the generalized radial return for anisotropic elasto-plasticity model."
24 : "This class can be used in conjunction with other creep and plasticity materials for "
25 : "more complex simulations.");
26 :
27 192 : params.addRequiredParam<Real>("hardening_constant",
28 : "Hardening constant (H) for anisotropic plasticity");
29 192 : params.addParam<Real>(
30 192 : "hardening_exponent", 1.0, "Hardening exponent (n) for anisotropic plasticity");
31 192 : params.addRequiredParam<Real>("yield_stress",
32 : "Yield stress (constant value) for anisotropic plasticity");
33 192 : params.addParam<bool>(
34 : "local_cylindrical_csys",
35 192 : false,
36 : "Compute inelastic strain increment in local cylindrical coordinate system");
37 :
38 192 : params.addParam<MooseEnum>("axis",
39 288 : MooseEnum("x y z", "y"),
40 : "The axis of cylindrical component about which to rotate when "
41 : "computing inelastic strain increment in local coordinate system");
42 :
43 96 : return params;
44 0 : }
45 :
46 : template <bool is_ad>
47 72 : HillElastoPlasticityStressUpdateTempl<is_ad>::HillElastoPlasticityStressUpdateTempl(
48 : const InputParameters & parameters)
49 : : AnisotropicReturnPlasticityStressUpdateBaseTempl<is_ad>(parameters),
50 144 : _qsigma(0.0),
51 72 : _eigenvalues_hill(6),
52 72 : _eigenvectors_hill(6, 6),
53 72 : _anisotropic_elastic_tensor(6, 6),
54 72 : _elasticity_tensor_name(this->_base_name + "elasticity_tensor"),
55 72 : _elasticity_tensor(
56 72 : this->template getGenericMaterialProperty<RankFourTensor, is_ad>(_elasticity_tensor_name)),
57 144 : _hardening_constant(this->template getParam<Real>("hardening_constant")),
58 144 : _hardening_exponent(this->template getParam<Real>("hardening_exponent")),
59 144 : _hardening_variable(this->template declareGenericProperty<Real, is_ad>(this->_base_name +
60 : "hardening_variable")),
61 72 : _hardening_variable_old(
62 72 : this->template getMaterialPropertyOld<Real>(this->_base_name + "hardening_variable")),
63 144 : _elasticity_eigenvectors(this->template declareGenericProperty<DenseMatrix<Real>, is_ad>(
64 : this->_base_name + "elasticity_eigenvectors")),
65 144 : _elasticity_eigenvalues(this->template declareGenericProperty<DenseVector<Real>, is_ad>(
66 : this->_base_name + "elasticity_eigenvalues")),
67 144 : _b_eigenvectors(this->template declareGenericProperty<DenseMatrix<Real>, is_ad>(
68 : this->_base_name + "b_eigenvectors")),
69 144 : _b_eigenvalues(this->template declareGenericProperty<DenseVector<Real>, is_ad>(
70 : this->_base_name + "b_eigenvalues")),
71 144 : _alpha_matrix(this->template declareGenericProperty<DenseMatrix<Real>, is_ad>(this->_base_name +
72 : "alpha_matrix")),
73 144 : _sigma_tilde(this->template declareGenericProperty<DenseVector<Real>, is_ad>(this->_base_name +
74 : "sigma_tilde")),
75 144 : _sigma_tilde_rotated(this->template declareGenericProperty<DenseVector<Real>, is_ad>(
76 : this->_base_name + "sigma_tilde_rotated")),
77 72 : _hardening_derivative(0.0),
78 72 : _yield_condition(1.0),
79 144 : _yield_stress(this->template getParam<Real>("yield_stress")),
80 144 : _hill_tensor(this->template getMaterialPropertyByName<DenseMatrix<Real>>(this->_base_name +
81 : "hill_tensor")),
82 72 : _rotation_matrix(
83 72 : this->template declareProperty<RankTwoTensor>(this->_base_name + "rotation_matrix")),
84 72 : _rotation_matrix_transpose(this->template declareProperty<RankTwoTensor>(
85 : this->_base_name + "rotation_matrix_transpose")),
86 72 : _rotation_matrix_old(
87 72 : this->template getMaterialPropertyOld<RankTwoTensor>(this->_base_name + "rotation_matrix")),
88 144 : _rotation_matrix_transpose_old(this->template getMaterialPropertyOld<RankTwoTensor>(
89 : this->_base_name + "rotation_matrix_transpose")),
90 144 : _local_cylindrical_csys(this->template getParam<bool>("local_cylindrical_csys")),
91 216 : _axis(this->template getParam<MooseEnum>("axis").template getEnum<Axis>())
92 : {
93 180 : if (_local_cylindrical_csys && !parameters.isParamSetByUser("axis"))
94 : {
95 0 : this->paramError(
96 : "local_cylindrical_csys",
97 : "\nIf parameter local_cylindrical_csys is provided then parameter axis should be also "
98 : "provided.");
99 : }
100 :
101 90 : if (!_local_cylindrical_csys && parameters.isParamSetByUser("axis"))
102 : {
103 0 : this->paramError("axis",
104 : "\nIf parameter axis is provided then the parameter local_cylindrical_csys "
105 : "should be set to true");
106 : }
107 72 : }
108 :
109 : template <bool is_ad>
110 : void
111 7008 : HillElastoPlasticityStressUpdateTempl<is_ad>::initQpStatefulProperties()
112 : {
113 : RealVectorValue normal_vector(0, 0, 0);
114 : RealVectorValue axial_vector(0, 0, 0);
115 :
116 7008 : if (_local_cylindrical_csys)
117 : {
118 4320 : if (_axis == Axis::X)
119 : {
120 : normal_vector(0) = 0.0;
121 1440 : normal_vector(1) = _q_point[_qp](1);
122 1440 : normal_vector(2) = _q_point[_qp](2);
123 1440 : axial_vector(0) = 1.0;
124 : }
125 2880 : else if (_axis == Axis::Y)
126 : {
127 1440 : normal_vector(0) = _q_point[_qp](0);
128 : normal_vector(1) = 0.0;
129 1440 : normal_vector(2) = _q_point[_qp](2);
130 1440 : axial_vector(1) = 1.0;
131 : }
132 1440 : else if (_axis == Axis::Z)
133 : {
134 1440 : normal_vector(0) = _q_point[_qp](0);
135 1440 : normal_vector(1) = _q_point[_qp](1);
136 : normal_vector(2) = 0.0;
137 1440 : axial_vector(2) = 1.0;
138 : }
139 : else
140 0 : mooseError("\nInvalid value for axis parameter provided in HillElastoPlasticityStressUpdate");
141 : }
142 :
143 7008 : if (_local_cylindrical_csys)
144 : {
145 4320 : Real nv_norm = normal_vector.norm();
146 4320 : Real av_norm = axial_vector.norm();
147 :
148 4320 : if (!(MooseUtils::absoluteFuzzyEqual(nv_norm, 0.0)))
149 : normal_vector /= nv_norm;
150 : else
151 : {
152 0 : mooseError("The normal vector cannot be a zero vector in "
153 : "HillElastoPlasticityStressUpdate");
154 : }
155 :
156 4320 : if (!(MooseUtils::absoluteFuzzyEqual(av_norm, 0.0)))
157 : axial_vector /= av_norm;
158 : else
159 : {
160 0 : mooseError("The axial vector cannot be a zero vector in "
161 : "HillElastoPlasticityStressUpdate");
162 : }
163 :
164 4320 : RankTwoScalarTools::setRotationMatrix(
165 4320 : normal_vector, axial_vector, _rotation_matrix[_qp], false);
166 4320 : _rotation_matrix_transpose[_qp] = _rotation_matrix[_qp].transpose();
167 : }
168 :
169 7008 : AnisotropicReturnPlasticityStressUpdateBaseTempl<is_ad>::initQpStatefulProperties();
170 7008 : }
171 :
172 : template <bool is_ad>
173 : void
174 0 : HillElastoPlasticityStressUpdateTempl<is_ad>::propagateQpStatefulProperties()
175 : {
176 0 : _hardening_variable[_qp] = _hardening_variable_old[_qp];
177 0 : _plasticity_strain[_qp] = _plasticity_strain_old[_qp];
178 :
179 0 : if (_local_cylindrical_csys)
180 : {
181 0 : _rotation_matrix[_qp] = _rotation_matrix_old[_qp];
182 0 : _rotation_matrix_transpose[_qp] = _rotation_matrix_transpose_old[_qp];
183 : }
184 :
185 0 : AnisotropicReturnPlasticityStressUpdateBaseTempl<is_ad>::propagateQpStatefulProperties();
186 0 : }
187 :
188 : template <bool is_ad>
189 : void
190 108576 : HillElastoPlasticityStressUpdateTempl<is_ad>::computeStressInitialize(
191 : const GenericDenseVector<is_ad> & stress_dev,
192 : const GenericDenseVector<is_ad> & stress,
193 : const GenericRankFourTensor<is_ad> & /*elasticity_tensor*/)
194 : {
195 108576 : _hardening_variable[_qp] = _hardening_variable_old[_qp];
196 108576 : _plasticity_strain[_qp] = _plasticity_strain_old[_qp];
197 108576 : _effective_inelastic_strain[_qp] = _effective_inelastic_strain_old[_qp];
198 :
199 108576 : _elasticity_eigenvectors[_qp].resize(6, 6);
200 108576 : _elasticity_eigenvalues[_qp].resize(6);
201 :
202 108576 : _b_eigenvectors[_qp].resize(6, 6);
203 108576 : _b_eigenvalues[_qp].resize(6);
204 :
205 108576 : _alpha_matrix[_qp].resize(6, 6);
206 108576 : _sigma_tilde[_qp].resize(6);
207 108576 : _sigma_tilde_rotated[_qp].resize(6);
208 :
209 : // Algebra needed for the generalized return mapping of anisotropic elastoplasticity
210 108576 : computeElasticityTensorEigenDecomposition();
211 :
212 108576 : _yield_condition = 1.0; // Some positive value
213 108576 : _yield_condition = -computeResidual(stress_dev, stress, 0.0);
214 108576 : }
215 :
216 : template <bool is_ad>
217 : void
218 108576 : HillElastoPlasticityStressUpdateTempl<is_ad>::computeElasticityTensorEigenDecomposition()
219 : {
220 : // Helper method to compute qp matrices required for the elasto-plasticity algorithm.
221 :
222 108576 : GenericRankFourTensor<is_ad> elasticity_tensor = _elasticity_tensor[_qp];
223 :
224 108576 : if (_local_cylindrical_csys)
225 : {
226 49440 : elasticity_tensor.rotate(_rotation_matrix_transpose[_qp]);
227 : }
228 :
229 108576 : ElasticityTensorTools::toVoigtNotation<is_ad>(_anisotropic_elastic_tensor, elasticity_tensor);
230 :
231 : const unsigned int dimension = _anisotropic_elastic_tensor.n();
232 :
233 : AnisotropyMatrixReal A = AnisotropyMatrixReal::Zero();
234 760032 : for (unsigned int index_i = 0; index_i < dimension; index_i++)
235 4560192 : for (unsigned int index_j = 0; index_j < dimension; index_j++)
236 3908736 : A(index_i, index_j) = MetaPhysicL::raw_value(_anisotropic_elastic_tensor(index_i, index_j));
237 :
238 : Eigen::SelfAdjointEigenSolver<AnisotropyMatrixReal> es(A);
239 :
240 : auto lambda = es.eigenvalues();
241 : auto v = es.eigenvectors();
242 :
243 760032 : for (unsigned int index_i = 0; index_i < dimension; index_i++)
244 651456 : _elasticity_eigenvalues[_qp](index_i) = lambda(index_i);
245 :
246 760032 : for (unsigned int index_i = 0; index_i < dimension; index_i++)
247 4560192 : for (unsigned int index_j = 0; index_j < dimension; index_j++)
248 3908736 : _elasticity_eigenvectors[_qp](index_i, index_j) = v(index_i, index_j);
249 :
250 : // Compute sqrt(Delta_c) * QcT * A * Qc * sqrt(Delta_c)
251 108576 : GenericDenseMatrix<is_ad> sqrt_Delta(6, 6);
252 760032 : for (unsigned int i = 0; i < 6; i++)
253 : {
254 651456 : sqrt_Delta(i, i) = std::sqrt(_elasticity_eigenvalues[_qp](i));
255 : }
256 :
257 108576 : GenericDenseMatrix<is_ad> eigenvectors_elasticity_transpose(6, 6);
258 108576 : _elasticity_eigenvectors[_qp].get_transpose(eigenvectors_elasticity_transpose);
259 :
260 108576 : GenericDenseMatrix<is_ad> b_matrix(6, 6);
261 108576 : b_matrix = _hill_tensor[_qp];
262 :
263 : // Right multiply by matrix of eigenvectors transpose
264 108576 : b_matrix.right_multiply(_elasticity_eigenvectors[_qp]);
265 108576 : b_matrix.right_multiply(sqrt_Delta);
266 108576 : b_matrix.left_multiply(eigenvectors_elasticity_transpose);
267 108576 : b_matrix.left_multiply(sqrt_Delta);
268 :
269 : AnisotropyMatrixReal B_eigen = AnisotropyMatrixReal::Zero();
270 760032 : for (unsigned int index_i = 0; index_i < dimension; index_i++)
271 4560192 : for (unsigned int index_j = 0; index_j < dimension; index_j++)
272 3908736 : B_eigen(index_i, index_j) = MetaPhysicL::raw_value(b_matrix(index_i, index_j));
273 :
274 108576 : if (isBlockDiagonal(B_eigen))
275 : {
276 5248 : Eigen::SelfAdjointEigenSolver<AnisotropyMatrixRealBlock> es_b_left(B_eigen.block<3, 3>(0, 0));
277 5248 : Eigen::SelfAdjointEigenSolver<AnisotropyMatrixRealBlock> es_b_right(B_eigen.block<3, 3>(3, 3));
278 :
279 : auto lambda_b_left = es_b_left.eigenvalues();
280 : auto v_b_left = es_b_left.eigenvectors();
281 :
282 : auto lambda_b_right = es_b_right.eigenvalues();
283 : auto v_b_right = es_b_right.eigenvectors();
284 :
285 5248 : _b_eigenvalues[_qp](0) = lambda_b_left(0);
286 5248 : _b_eigenvalues[_qp](1) = lambda_b_left(1);
287 5248 : _b_eigenvalues[_qp](2) = lambda_b_left(2);
288 5248 : _b_eigenvalues[_qp](3) = lambda_b_right(0);
289 5248 : _b_eigenvalues[_qp](4) = lambda_b_right(1);
290 5248 : _b_eigenvalues[_qp](5) = lambda_b_right(2);
291 :
292 5248 : _b_eigenvectors[_qp](0, 0) = v_b_left(0, 0);
293 5248 : _b_eigenvectors[_qp](0, 1) = v_b_left(0, 1);
294 5248 : _b_eigenvectors[_qp](0, 2) = v_b_left(0, 2);
295 5248 : _b_eigenvectors[_qp](1, 0) = v_b_left(1, 0);
296 5248 : _b_eigenvectors[_qp](1, 1) = v_b_left(1, 1);
297 5248 : _b_eigenvectors[_qp](1, 2) = v_b_left(1, 2);
298 5248 : _b_eigenvectors[_qp](2, 0) = v_b_left(2, 0);
299 5248 : _b_eigenvectors[_qp](2, 1) = v_b_left(2, 1);
300 5248 : _b_eigenvectors[_qp](2, 2) = v_b_left(2, 2);
301 5248 : _b_eigenvectors[_qp](3, 3) = v_b_right(0, 0);
302 5248 : _b_eigenvectors[_qp](3, 4) = v_b_right(0, 1);
303 5248 : _b_eigenvectors[_qp](3, 5) = v_b_right(0, 2);
304 5248 : _b_eigenvectors[_qp](4, 3) = v_b_right(1, 0);
305 5248 : _b_eigenvectors[_qp](4, 4) = v_b_right(1, 1);
306 5248 : _b_eigenvectors[_qp](4, 5) = v_b_right(1, 2);
307 5248 : _b_eigenvectors[_qp](5, 3) = v_b_right(2, 0);
308 5248 : _b_eigenvectors[_qp](5, 4) = v_b_right(2, 1);
309 5248 : _b_eigenvectors[_qp](5, 5) = v_b_right(2, 2);
310 : }
311 : else
312 : {
313 : Eigen::SelfAdjointEigenSolver<AnisotropyMatrixReal> es_b(B_eigen);
314 :
315 : auto lambda_b = es_b.eigenvalues();
316 : auto v_b = es_b.eigenvectors();
317 723296 : for (unsigned int index_i = 0; index_i < dimension; index_i++)
318 619968 : _b_eigenvalues[_qp](index_i) = lambda_b(index_i);
319 :
320 723296 : for (unsigned int index_i = 0; index_i < dimension; index_i++)
321 4339776 : for (unsigned int index_j = 0; index_j < dimension; index_j++)
322 3719808 : _b_eigenvectors[_qp](index_i, index_j) = v_b(index_i, index_j);
323 : }
324 :
325 108576 : _alpha_matrix[_qp] = sqrt_Delta;
326 108576 : _alpha_matrix[_qp].right_multiply(_b_eigenvectors[_qp]);
327 108576 : _alpha_matrix[_qp].left_multiply(_elasticity_eigenvectors[_qp]);
328 108576 : }
329 :
330 : template <bool is_ad>
331 : GenericReal<is_ad>
332 316688 : HillElastoPlasticityStressUpdateTempl<is_ad>::computeOmega(
333 : const GenericReal<is_ad> & delta_gamma, const GenericDenseVector<is_ad> & /*stress_trial*/)
334 : {
335 316688 : GenericDenseVector<is_ad> K(6);
336 316688 : GenericReal<is_ad> omega = 0.0;
337 :
338 2216816 : for (unsigned int i = 0; i < 6; i++)
339 : {
340 7600512 : K(i) = _b_eigenvalues[_qp](i) / Utility::pow<2>(1 + delta_gamma * _b_eigenvalues[_qp](i));
341 :
342 1900128 : if (_local_cylindrical_csys)
343 2620800 : omega += K(i) * _sigma_tilde_rotated[_qp](i) * _sigma_tilde_rotated[_qp](i);
344 : else
345 1179456 : omega += K(i) * _sigma_tilde[_qp](i) * _sigma_tilde[_qp](i);
346 : }
347 316688 : omega *= 0.5;
348 :
349 316688 : if (omega == 0.0)
350 12448 : omega = 1.0e-36;
351 :
352 633376 : return std::sqrt(omega);
353 : }
354 :
355 : template <bool is_ad>
356 : Real
357 149952 : HillElastoPlasticityStressUpdateTempl<is_ad>::computeReferenceResidual(
358 : const GenericDenseVector<is_ad> & /*effective_trial_stress*/,
359 : const GenericDenseVector<is_ad> & /*stress_new*/,
360 : const GenericReal<is_ad> & /*residual*/,
361 : const GenericReal<is_ad> & /*scalar_effective_inelastic_strain*/)
362 : {
363 : // Equation is already normalized.
364 149952 : return 1.0;
365 : }
366 :
367 : template <bool is_ad>
368 : GenericReal<is_ad>
369 258528 : HillElastoPlasticityStressUpdateTempl<is_ad>::computeResidual(
370 : const GenericDenseVector<is_ad> & /*stress_dev*/,
371 : const GenericDenseVector<is_ad> & stress_sigma,
372 : const GenericReal<is_ad> & delta_gamma)
373 : {
374 :
375 : // If in elastic regime, just return
376 258528 : if (_yield_condition <= 0.0)
377 49488 : return 0.0;
378 :
379 : GenericReal<is_ad> omega;
380 :
381 209040 : if (_local_cylindrical_csys)
382 : {
383 : // Get stress_tilde_rotated
384 133920 : GenericRankTwoTensor<is_ad> stress;
385 133920 : RankTwoScalarTools::VoigtVectorToRankTwoTensor<is_ad>(stress_sigma, stress);
386 :
387 133920 : stress.rotate(_rotation_matrix_transpose[_qp]);
388 :
389 133920 : GenericDenseVector<is_ad> stress_sigma_rotated(6);
390 133920 : RankTwoScalarTools::RankTwoTensorToVoigtVector<is_ad>(stress, stress_sigma_rotated);
391 :
392 133920 : GenericDenseVector<is_ad> stress_tilde_rotated(6);
393 133920 : GenericDenseMatrix<is_ad> alpha_temp_rotated(_alpha_matrix[_qp]);
394 133920 : alpha_temp_rotated.lu_solve(stress_sigma_rotated, stress_tilde_rotated);
395 :
396 : // Material property used in computeStressFinalize
397 133920 : _sigma_tilde_rotated[_qp] = stress_tilde_rotated;
398 133920 : omega = computeOmega(delta_gamma, stress_tilde_rotated);
399 133920 : }
400 :
401 : else
402 : {
403 : // Get stress_tilde
404 75120 : GenericDenseVector<is_ad> stress_tilde(6);
405 75120 : GenericDenseMatrix<is_ad> alpha_temp(_alpha_matrix[_qp]);
406 75120 : alpha_temp.lu_solve(stress_sigma, stress_tilde);
407 :
408 : // Material property used in computeStressFinalize
409 75120 : _sigma_tilde[_qp] = stress_tilde;
410 75120 : omega = computeOmega(delta_gamma, stress_tilde);
411 75120 : }
412 :
413 : // Hardening variable is \alpha isotropic hardening for now.
414 209040 : _hardening_variable[_qp] = computeHardeningValue(delta_gamma, omega);
415 :
416 : // A small value of 1.0e-30 is added to the hardening variable to avoid numerical issues
417 : // related to hardening variable becoming negative early in the iteration leading to non-
418 : // convergence. Note that std::pow(x,n) requires x to be positive if n is less than 1.
419 :
420 209040 : GenericReal<is_ad> s_y =
421 209040 : _hardening_constant * std::pow(_hardening_variable[_qp] + 1.0e-30, _hardening_exponent) +
422 209040 : _yield_stress;
423 :
424 209040 : GenericReal<is_ad> residual = 0.0;
425 209040 : residual = s_y / omega - 1.0;
426 :
427 209040 : return residual;
428 : }
429 :
430 : template <bool is_ad>
431 : GenericReal<is_ad>
432 53824 : HillElastoPlasticityStressUpdateTempl<is_ad>::computeDerivative(
433 : const GenericDenseVector<is_ad> & /*stress_dev*/,
434 : const GenericDenseVector<is_ad> & stress_sigma,
435 : const GenericReal<is_ad> & delta_gamma)
436 : {
437 : // If in elastic regime, return unit derivative
438 53824 : if (_yield_condition <= 0.0)
439 0 : return 1.0;
440 :
441 53824 : GenericReal<is_ad> omega = computeOmega(delta_gamma, stress_sigma);
442 53824 : _hardening_derivative = computeHardeningDerivative();
443 :
444 53824 : GenericReal<is_ad> hardeningVariable = computeHardeningValue(delta_gamma, omega);
445 53824 : GenericReal<is_ad> sy =
446 53824 : _hardening_constant * std::pow(hardeningVariable + 1.0e-30, _hardening_exponent) +
447 53824 : _yield_stress;
448 53824 : GenericReal<is_ad> sy_alpha = _hardening_derivative;
449 :
450 : GenericReal<is_ad> omega_gamma;
451 : GenericReal<is_ad> sy_gamma;
452 :
453 53824 : computeDeltaDerivatives(delta_gamma, stress_sigma, sy_alpha, omega, omega_gamma, sy_gamma);
454 107648 : GenericReal<is_ad> residual_derivative = 1 / omega * (sy_gamma - 1 / omega * omega_gamma * sy);
455 :
456 53824 : return residual_derivative;
457 : }
458 :
459 : template <bool is_ad>
460 : void
461 53824 : HillElastoPlasticityStressUpdateTempl<is_ad>::computeDeltaDerivatives(
462 : const GenericReal<is_ad> & delta_gamma,
463 : const GenericDenseVector<is_ad> & /*stress_trial*/,
464 : const GenericReal<is_ad> & sy_alpha,
465 : GenericReal<is_ad> & omega,
466 : GenericReal<is_ad> & omega_gamma,
467 : GenericReal<is_ad> & sy_gamma)
468 : {
469 53824 : omega_gamma = 0.0;
470 53824 : sy_gamma = 0.0;
471 :
472 53824 : GenericDenseVector<is_ad> K_deltaGamma(6);
473 :
474 53824 : if (_local_cylindrical_csys)
475 42240 : omega = computeOmega(delta_gamma, _sigma_tilde_rotated[_qp]);
476 : else
477 11584 : omega = computeOmega(delta_gamma, _sigma_tilde[_qp]);
478 :
479 53824 : GenericDenseVector<is_ad> K(6);
480 376768 : for (unsigned int i = 0; i < 6; i++)
481 1291776 : K(i) = _b_eigenvalues[_qp](i) / Utility::pow<2>(1 + delta_gamma * _b_eigenvalues[_qp](i));
482 :
483 376768 : for (unsigned int i = 0; i < 6; i++)
484 322944 : K_deltaGamma(i) =
485 1291776 : -2.0 * _b_eigenvalues[_qp](i) * K(i) / (1 + _b_eigenvalues[_qp](i) * delta_gamma);
486 :
487 376768 : for (unsigned int i = 0; i < 6; i++)
488 : {
489 322944 : if (_local_cylindrical_csys)
490 506880 : omega_gamma += K_deltaGamma(i) * _sigma_tilde_rotated[_qp](i) * _sigma_tilde_rotated[_qp](i);
491 : else
492 139008 : omega_gamma += K_deltaGamma(i) * _sigma_tilde[_qp](i) * _sigma_tilde[_qp](i);
493 : }
494 :
495 107648 : omega_gamma /= 4.0 * omega;
496 107648 : sy_gamma = 2.0 * sy_alpha * (omega + delta_gamma * omega_gamma);
497 53824 : }
498 :
499 : template <bool is_ad>
500 : GenericReal<is_ad>
501 262864 : HillElastoPlasticityStressUpdateTempl<is_ad>::computeHardeningValue(
502 : const GenericReal<is_ad> & delta_gamma, const GenericReal<is_ad> & omega)
503 : {
504 525728 : return _hardening_variable_old[_qp] + 2.0 * delta_gamma * omega;
505 : }
506 :
507 : template <bool is_ad>
508 : Real
509 53824 : HillElastoPlasticityStressUpdateTempl<is_ad>::computeHardeningDerivative()
510 : {
511 53824 : return _hardening_constant * _hardening_exponent *
512 53824 : MetaPhysicL::raw_value(
513 53824 : std::pow(_hardening_variable[_qp] + 1.0e-30, _hardening_exponent - 1));
514 : }
515 :
516 : template <bool is_ad>
517 : void
518 96128 : HillElastoPlasticityStressUpdateTempl<is_ad>::computeStrainFinalize(
519 : GenericRankTwoTensor<is_ad> & inelasticStrainIncrement,
520 : const GenericRankTwoTensor<is_ad> & stress,
521 : const GenericDenseVector<is_ad> & stress_dev,
522 : const GenericReal<is_ad> & delta_gamma)
523 : {
524 : // e^P = delta_gamma * hill_tensor * stress
525 96128 : GenericDenseVector<is_ad> inelasticStrainIncrement_vector(6);
526 96128 : GenericDenseVector<is_ad> hill_stress(6);
527 96128 : GenericDenseVector<is_ad> stress_vector(6);
528 :
529 96128 : if (_local_cylindrical_csys)
530 : {
531 0 : GenericRankTwoTensor<is_ad> stress_rotated(stress);
532 42240 : stress_rotated.rotate(_rotation_matrix_transpose[_qp]);
533 :
534 42240 : RankTwoScalarTools::RankTwoTensorToVoigtVector<is_ad>(stress_rotated, stress_vector);
535 : }
536 : else
537 : {
538 53888 : RankTwoScalarTools::RankTwoTensorToVoigtVector<is_ad>(stress, stress_vector);
539 : }
540 :
541 96128 : _hill_tensor[_qp].vector_mult(hill_stress, stress_vector);
542 96128 : hill_stress.scale(delta_gamma);
543 : inelasticStrainIncrement_vector = hill_stress;
544 :
545 96128 : inelasticStrainIncrement(0, 0) = inelasticStrainIncrement_vector(0);
546 96128 : inelasticStrainIncrement(1, 1) = inelasticStrainIncrement_vector(1);
547 96128 : inelasticStrainIncrement(2, 2) = inelasticStrainIncrement_vector(2);
548 96128 : inelasticStrainIncrement(0, 1) = inelasticStrainIncrement(1, 0) =
549 96128 : inelasticStrainIncrement_vector(3) / 2.0;
550 96128 : inelasticStrainIncrement(1, 2) = inelasticStrainIncrement(2, 1) =
551 96128 : inelasticStrainIncrement_vector(4) / 2.0;
552 96128 : inelasticStrainIncrement(0, 2) = inelasticStrainIncrement(2, 0) =
553 96128 : inelasticStrainIncrement_vector(5) / 2.0;
554 :
555 96128 : if (_local_cylindrical_csys)
556 42240 : inelasticStrainIncrement.rotate(_rotation_matrix[_qp]);
557 :
558 : // Calculate equivalent plastic strain
559 96128 : GenericDenseVector<is_ad> Mepsilon(6);
560 96128 : _hill_tensor[_qp].vector_mult(Mepsilon, inelasticStrainIncrement_vector);
561 96128 : GenericReal<is_ad> eq_plastic_strain_inc = Mepsilon.dot(inelasticStrainIncrement_vector);
562 :
563 96128 : if (eq_plastic_strain_inc > 0.0)
564 46640 : eq_plastic_strain_inc = std::sqrt(eq_plastic_strain_inc);
565 :
566 192256 : _effective_inelastic_strain[_qp] = _effective_inelastic_strain_old[_qp] + eq_plastic_strain_inc;
567 :
568 96128 : AnisotropicReturnPlasticityStressUpdateBaseTempl<is_ad>::computeStrainFinalize(
569 : inelasticStrainIncrement, stress, stress_dev, delta_gamma);
570 96128 : }
571 :
572 : template <bool is_ad>
573 : void
574 108576 : HillElastoPlasticityStressUpdateTempl<is_ad>::computeStressFinalize(
575 : const GenericRankTwoTensor<is_ad> & /*plastic_strain_increment*/,
576 : const GenericReal<is_ad> & delta_gamma,
577 : GenericRankTwoTensor<is_ad> & stress_new,
578 : const GenericDenseVector<is_ad> & /*stress_dev*/,
579 : const GenericRankTwoTensor<is_ad> & /*stress_old*/,
580 : const GenericRankFourTensor<is_ad> & /*elasticity_tensor*/)
581 : {
582 : // Need to compute this iteration's stress tensor based on the scalar variable
583 : // For deviatoric
584 : // sigma(n+1) = {Alpha [I + delta_gamma*Delta_b]^(-1) A^-1} sigma(trial)
585 :
586 108576 : if (_yield_condition <= 0.0)
587 54736 : return;
588 53840 : GenericDenseMatrix<is_ad> inv_matrix(6, 6);
589 53840 : inv_matrix.zero();
590 :
591 376880 : for (unsigned int i = 0; i < 6; i++)
592 1292160 : inv_matrix(i, i) = 1 / (1 + delta_gamma * _b_eigenvalues[_qp](i));
593 :
594 53840 : _alpha_matrix[_qp].right_multiply(inv_matrix);
595 :
596 53840 : GenericDenseVector<is_ad> stress_output(6);
597 53840 : if (_local_cylindrical_csys)
598 49440 : _alpha_matrix[_qp].vector_mult(stress_output, _sigma_tilde_rotated[_qp]);
599 : else
600 4400 : _alpha_matrix[_qp].vector_mult(stress_output, _sigma_tilde[_qp]);
601 :
602 53840 : RankTwoScalarTools::VoigtVectorToRankTwoTensor<is_ad>(stress_output, stress_new);
603 :
604 53840 : if (_local_cylindrical_csys)
605 49440 : stress_new.rotate(_rotation_matrix[_qp]);
606 53840 : }
607 :
608 : template <bool is_ad>
609 : Real
610 0 : HillElastoPlasticityStressUpdateTempl<is_ad>::computeStrainEnergyRateDensity(
611 : const GenericMaterialProperty<RankTwoTensor, is_ad> & /*stress*/,
612 : const GenericMaterialProperty<RankTwoTensor, is_ad> & /*strain_rate*/)
613 : {
614 0 : mooseError("computeStrainEnergyRateDensity not implemented for anisotropic plasticity.");
615 : return 0.0;
616 : }
617 :
618 : template class HillElastoPlasticityStressUpdateTempl<false>;
619 : template class HillElastoPlasticityStressUpdateTempl<true>;
|
