Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://www.mooseframework.org
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("TensorMechanicsApp", ADHillElastoPlasticityStressUpdate);
15 : registerMooseObject("TensorMechanicsApp", HillElastoPlasticityStressUpdate);
16 :
17 : template <bool is_ad>
18 : InputParameters
19 48 : HillElastoPlasticityStressUpdateTempl<is_ad>::validParams()
20 : {
21 48 : InputParameters params = AnisotropicReturnPlasticityStressUpdateBaseTempl<is_ad>::validParams();
22 48 : 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 96 : params.addRequiredParam<Real>("hardening_constant",
28 : "Hardening constant (H) for anisotropic plasticity");
29 96 : params.addParam<Real>(
30 96 : "hardening_exponent", 1.0, "Hardening exponent (n) for anisotropic plasticity");
31 96 : params.addRequiredParam<Real>("yield_stress",
32 : "Yield stress (constant value) for anisotropic plasticity");
33 96 : params.addParam<bool>(
34 : "local_cylindrical_csys",
35 96 : false,
36 : "Compute inelastic strain increment in local cylindrical coordinate system");
37 :
38 96 : params.addParam<MooseEnum>("axis",
39 144 : 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 48 : return params;
44 0 : }
45 :
46 : template <bool is_ad>
47 36 : HillElastoPlasticityStressUpdateTempl<is_ad>::HillElastoPlasticityStressUpdateTempl(
48 : const InputParameters & parameters)
49 : : AnisotropicReturnPlasticityStressUpdateBaseTempl<is_ad>(parameters),
50 72 : _qsigma(0.0),
51 36 : _eigenvalues_hill(6),
52 36 : _eigenvectors_hill(6, 6),
53 36 : _anisotropic_elastic_tensor(6, 6),
54 36 : _elasticity_tensor_name(this->_base_name + "elasticity_tensor"),
55 36 : _elasticity_tensor(
56 36 : this->template getGenericMaterialProperty<RankFourTensor, is_ad>(_elasticity_tensor_name)),
57 72 : _hardening_constant(this->template getParam<Real>("hardening_constant")),
58 72 : _hardening_exponent(this->template getParam<Real>("hardening_exponent")),
59 72 : _hardening_variable(this->template declareGenericProperty<Real, is_ad>(this->_base_name +
60 : "hardening_variable")),
61 36 : _hardening_variable_old(
62 36 : this->template getMaterialPropertyOld<Real>(this->_base_name + "hardening_variable")),
63 72 : _elasticity_eigenvectors(this->template declareGenericProperty<DenseMatrix<Real>, is_ad>(
64 : this->_base_name + "elasticity_eigenvectors")),
65 72 : _elasticity_eigenvalues(this->template declareGenericProperty<DenseVector<Real>, is_ad>(
66 : this->_base_name + "elasticity_eigenvalues")),
67 72 : _b_eigenvectors(this->template declareGenericProperty<DenseMatrix<Real>, is_ad>(
68 : this->_base_name + "b_eigenvectors")),
69 72 : _b_eigenvalues(this->template declareGenericProperty<DenseVector<Real>, is_ad>(
70 : this->_base_name + "b_eigenvalues")),
71 72 : _alpha_matrix(this->template declareGenericProperty<DenseMatrix<Real>, is_ad>(this->_base_name +
72 : "alpha_matrix")),
73 72 : _sigma_tilde(this->template declareGenericProperty<DenseVector<Real>, is_ad>(this->_base_name +
74 : "sigma_tilde")),
75 72 : _sigma_tilde_rotated(this->template declareGenericProperty<DenseVector<Real>, is_ad>(
76 : this->_base_name + "sigma_tilde_rotated")),
77 36 : _hardening_derivative(0.0),
78 36 : _yield_condition(1.0),
79 72 : _yield_stress(this->template getParam<Real>("yield_stress")),
80 72 : _hill_tensor(this->template getMaterialPropertyByName<DenseMatrix<Real>>(this->_base_name +
81 : "hill_tensor")),
82 36 : _rotation_matrix(
83 36 : this->template declareProperty<RankTwoTensor>(this->_base_name + "rotation_matrix")),
84 36 : _rotation_matrix_transpose(this->template declareProperty<RankTwoTensor>(
85 : this->_base_name + "rotation_matrix_transpose")),
86 36 : _rotation_matrix_old(
87 36 : this->template getMaterialPropertyOld<RankTwoTensor>(this->_base_name + "rotation_matrix")),
88 72 : _rotation_matrix_transpose_old(this->template getMaterialPropertyOld<RankTwoTensor>(
89 : this->_base_name + "rotation_matrix_transpose")),
90 72 : _local_cylindrical_csys(this->template getParam<bool>("local_cylindrical_csys")),
91 108 : _axis(this->template getParam<MooseEnum>("axis").template getEnum<Axis>())
92 : {
93 90 : 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 45 : 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 36 : }
108 :
109 : template <bool is_ad>
110 : void
111 3504 : HillElastoPlasticityStressUpdateTempl<is_ad>::initQpStatefulProperties()
112 : {
113 : RealVectorValue normal_vector(0, 0, 0);
114 : RealVectorValue axial_vector(0, 0, 0);
115 :
116 3504 : if (_local_cylindrical_csys)
117 : {
118 2160 : if (_axis == Axis::X)
119 : {
120 : normal_vector(0) = 0.0;
121 720 : normal_vector(1) = _q_point[_qp](1);
122 720 : normal_vector(2) = _q_point[_qp](2);
123 720 : axial_vector(0) = 1.0;
124 : }
125 1440 : else if (_axis == Axis::Y)
126 : {
127 720 : normal_vector(0) = _q_point[_qp](0);
128 : normal_vector(1) = 0.0;
129 720 : normal_vector(2) = _q_point[_qp](2);
130 720 : axial_vector(1) = 1.0;
131 : }
132 720 : else if (_axis == Axis::Z)
133 : {
134 720 : normal_vector(0) = _q_point[_qp](0);
135 720 : normal_vector(1) = _q_point[_qp](1);
136 : normal_vector(2) = 0.0;
137 720 : axial_vector(2) = 1.0;
138 : }
139 : else
140 0 : mooseError("\nInvalid value for axis parameter provided in HillElastoPlasticityStressUpdate");
141 : }
142 :
143 3504 : if (_local_cylindrical_csys)
144 : {
145 2160 : Real nv_norm = normal_vector.norm();
146 2160 : Real av_norm = axial_vector.norm();
147 :
148 2160 : 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 2160 : 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 2160 : RankTwoScalarTools::setRotationMatrix(
165 2160 : normal_vector, axial_vector, _rotation_matrix[_qp], false);
166 2160 : _rotation_matrix_transpose[_qp] = _rotation_matrix[_qp].transpose();
167 : }
168 :
169 3504 : AnisotropicReturnPlasticityStressUpdateBaseTempl<is_ad>::initQpStatefulProperties();
170 3504 : }
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 60024 : 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 60024 : _hardening_variable[_qp] = _hardening_variable_old[_qp];
196 60024 : _plasticity_strain[_qp] = _plasticity_strain_old[_qp];
197 60024 : _effective_inelastic_strain[_qp] = _effective_inelastic_strain_old[_qp];
198 :
199 60024 : _elasticity_eigenvectors[_qp].resize(6, 6);
200 60024 : _elasticity_eigenvalues[_qp].resize(6);
201 :
202 60024 : _b_eigenvectors[_qp].resize(6, 6);
203 60024 : _b_eigenvalues[_qp].resize(6);
204 :
205 60024 : _alpha_matrix[_qp].resize(6, 6);
206 60024 : _sigma_tilde[_qp].resize(6);
207 60024 : _sigma_tilde_rotated[_qp].resize(6);
208 :
209 : // Algebra needed for the generalized return mapping of anisotropic elastoplasticity
210 60024 : computeElasticityTensorEigenDecomposition();
211 :
212 60024 : _yield_condition = 1.0; // Some positive value
213 60024 : _yield_condition = -computeResidual(stress_dev, stress, 0.0);
214 60024 : }
215 :
216 : template <bool is_ad>
217 : void
218 60024 : HillElastoPlasticityStressUpdateTempl<is_ad>::computeElasticityTensorEigenDecomposition()
219 : {
220 : // Helper method to compute qp matrices required for the elasto-plasticity algorithm.
221 :
222 60024 : GenericRankFourTensor<is_ad> elasticity_tensor = _elasticity_tensor[_qp];
223 :
224 60024 : if (_local_cylindrical_csys)
225 : {
226 26520 : elasticity_tensor.rotate(_rotation_matrix_transpose[_qp]);
227 : }
228 :
229 60024 : 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 420168 : for (unsigned int index_i = 0; index_i < dimension; index_i++)
235 2521008 : for (unsigned int index_j = 0; index_j < dimension; index_j++)
236 2160864 : 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 420168 : for (unsigned int index_i = 0; index_i < dimension; index_i++)
244 360144 : _elasticity_eigenvalues[_qp](index_i) = lambda(index_i);
245 :
246 420168 : for (unsigned int index_i = 0; index_i < dimension; index_i++)
247 2521008 : for (unsigned int index_j = 0; index_j < dimension; index_j++)
248 2160864 : _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 60024 : GenericDenseMatrix<is_ad> sqrt_Delta(6, 6);
252 420168 : for (unsigned int i = 0; i < 6; i++)
253 : {
254 360144 : sqrt_Delta(i, i) = std::sqrt(_elasticity_eigenvalues[_qp](i));
255 : }
256 :
257 60024 : GenericDenseMatrix<is_ad> eigenvectors_elasticity_transpose(6, 6);
258 60024 : _elasticity_eigenvectors[_qp].get_transpose(eigenvectors_elasticity_transpose);
259 :
260 60024 : GenericDenseMatrix<is_ad> b_matrix(6, 6);
261 60024 : b_matrix = _hill_tensor[_qp];
262 :
263 : // Right multiply by matrix of eigenvectors transpose
264 60024 : b_matrix.right_multiply(_elasticity_eigenvectors[_qp]);
265 60024 : b_matrix.right_multiply(sqrt_Delta);
266 60024 : b_matrix.left_multiply(eigenvectors_elasticity_transpose);
267 60024 : b_matrix.left_multiply(sqrt_Delta);
268 :
269 : AnisotropyMatrixReal B_eigen = AnisotropyMatrixReal::Zero();
270 420168 : for (unsigned int index_i = 0; index_i < dimension; index_i++)
271 2521008 : for (unsigned int index_j = 0; index_j < dimension; index_j++)
272 2160864 : B_eigen(index_i, index_j) = MetaPhysicL::raw_value(b_matrix(index_i, index_j));
273 :
274 60024 : if (isBlockDiagonal(B_eigen))
275 : {
276 3936 : Eigen::SelfAdjointEigenSolver<AnisotropyMatrixRealBlock> es_b_left(B_eigen.block<3, 3>(0, 0));
277 3936 : 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 3936 : _b_eigenvalues[_qp](0) = lambda_b_left(0);
286 3936 : _b_eigenvalues[_qp](1) = lambda_b_left(1);
287 3936 : _b_eigenvalues[_qp](2) = lambda_b_left(2);
288 3936 : _b_eigenvalues[_qp](3) = lambda_b_right(0);
289 3936 : _b_eigenvalues[_qp](4) = lambda_b_right(1);
290 3936 : _b_eigenvalues[_qp](5) = lambda_b_right(2);
291 :
292 3936 : _b_eigenvectors[_qp](0, 0) = v_b_left(0, 0);
293 3936 : _b_eigenvectors[_qp](0, 1) = v_b_left(0, 1);
294 3936 : _b_eigenvectors[_qp](0, 2) = v_b_left(0, 2);
295 3936 : _b_eigenvectors[_qp](1, 0) = v_b_left(1, 0);
296 3936 : _b_eigenvectors[_qp](1, 1) = v_b_left(1, 1);
297 3936 : _b_eigenvectors[_qp](1, 2) = v_b_left(1, 2);
298 3936 : _b_eigenvectors[_qp](2, 0) = v_b_left(2, 0);
299 3936 : _b_eigenvectors[_qp](2, 1) = v_b_left(2, 1);
300 3936 : _b_eigenvectors[_qp](2, 2) = v_b_left(2, 2);
301 3936 : _b_eigenvectors[_qp](3, 3) = v_b_right(0, 0);
302 3936 : _b_eigenvectors[_qp](3, 4) = v_b_right(0, 1);
303 3936 : _b_eigenvectors[_qp](3, 5) = v_b_right(0, 2);
304 3936 : _b_eigenvectors[_qp](4, 3) = v_b_right(1, 0);
305 3936 : _b_eigenvectors[_qp](4, 4) = v_b_right(1, 1);
306 3936 : _b_eigenvectors[_qp](4, 5) = v_b_right(1, 2);
307 3936 : _b_eigenvectors[_qp](5, 3) = v_b_right(2, 0);
308 3936 : _b_eigenvectors[_qp](5, 4) = v_b_right(2, 1);
309 3936 : _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 392616 : for (unsigned int index_i = 0; index_i < dimension; index_i++)
318 336528 : _b_eigenvalues[_qp](index_i) = lambda_b(index_i);
319 :
320 392616 : for (unsigned int index_i = 0; index_i < dimension; index_i++)
321 2355696 : for (unsigned int index_j = 0; index_j < dimension; index_j++)
322 2019168 : _b_eigenvectors[_qp](index_i, index_j) = v_b(index_i, index_j);
323 : }
324 :
325 60024 : _alpha_matrix[_qp] = sqrt_Delta;
326 60024 : _alpha_matrix[_qp].right_multiply(_b_eigenvectors[_qp]);
327 60024 : _alpha_matrix[_qp].left_multiply(_elasticity_eigenvectors[_qp]);
328 60024 : }
329 :
330 : template <bool is_ad>
331 : GenericReal<is_ad>
332 164080 : HillElastoPlasticityStressUpdateTempl<is_ad>::computeOmega(
333 : const GenericReal<is_ad> & delta_gamma, const GenericDenseVector<is_ad> & /*stress_trial*/)
334 : {
335 164080 : GenericDenseVector<is_ad> K(6);
336 164080 : GenericReal<is_ad> omega = 0.0;
337 :
338 1148560 : for (unsigned int i = 0; i < 6; i++)
339 : {
340 3937920 : K(i) = _b_eigenvalues[_qp](i) / Utility::pow<2>(1 + delta_gamma * _b_eigenvalues[_qp](i));
341 :
342 984480 : if (_local_cylindrical_csys)
343 1332000 : omega += K(i) * _sigma_tilde_rotated[_qp](i) * _sigma_tilde_rotated[_qp](i);
344 : else
345 636960 : omega += K(i) * _sigma_tilde[_qp](i) * _sigma_tilde[_qp](i);
346 : }
347 164080 : omega *= 0.5;
348 :
349 164080 : if (omega == 0.0)
350 9336 : omega = 1.0e-36;
351 :
352 328160 : return std::sqrt(omega);
353 : }
354 :
355 : template <bool is_ad>
356 : Real
357 77600 : 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 77600 : return 1.0;
365 : }
366 :
367 : template <bool is_ad>
368 : GenericReal<is_ad>
369 137624 : 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 137624 : if (_yield_condition <= 0.0)
377 27368 : return 0.0;
378 :
379 : GenericReal<is_ad> omega;
380 :
381 110256 : if (_local_cylindrical_csys)
382 : {
383 : // Get stress_tilde_rotated
384 68760 : GenericRankTwoTensor<is_ad> stress;
385 68760 : RankTwoScalarTools::VoigtVectorToRankTwoTensor<is_ad>(stress_sigma, stress);
386 :
387 68760 : stress.rotate(_rotation_matrix_transpose[_qp]);
388 :
389 68760 : GenericDenseVector<is_ad> stress_sigma_rotated(6);
390 68760 : RankTwoScalarTools::RankTwoTensorToVoigtVector<is_ad>(stress, stress_sigma_rotated);
391 :
392 68760 : GenericDenseVector<is_ad> stress_tilde_rotated(6);
393 68760 : GenericDenseMatrix<is_ad> alpha_temp_rotated(_alpha_matrix[_qp]);
394 68760 : alpha_temp_rotated.lu_solve(stress_sigma_rotated, stress_tilde_rotated);
395 :
396 : // Material property used in computeStressFinalize
397 68760 : _sigma_tilde_rotated[_qp] = stress_tilde_rotated;
398 68760 : omega = computeOmega(delta_gamma, stress_tilde_rotated);
399 68760 : }
400 :
401 : else
402 : {
403 : // Get stress_tilde
404 41496 : GenericDenseVector<is_ad> stress_tilde(6);
405 41496 : GenericDenseMatrix<is_ad> alpha_temp(_alpha_matrix[_qp]);
406 41496 : alpha_temp.lu_solve(stress_sigma, stress_tilde);
407 :
408 : // Material property used in computeStressFinalize
409 41496 : _sigma_tilde[_qp] = stress_tilde;
410 41496 : omega = computeOmega(delta_gamma, stress_tilde);
411 41496 : }
412 :
413 : // Hardening variable is \alpha isotropic hardening for now.
414 110256 : _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 110256 : GenericReal<is_ad> s_y =
421 110256 : _hardening_constant * std::pow(_hardening_variable[_qp] + 1.0e-30, _hardening_exponent) +
422 110256 : _yield_stress;
423 :
424 110256 : GenericReal<is_ad> residual = 0.0;
425 110256 : residual = s_y / omega - 1.0;
426 :
427 110256 : return residual;
428 : }
429 :
430 : template <bool is_ad>
431 : GenericReal<is_ad>
432 26912 : 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 26912 : if (_yield_condition <= 0.0)
439 0 : return 1.0;
440 :
441 26912 : GenericReal<is_ad> omega = computeOmega(delta_gamma, stress_sigma);
442 26912 : _hardening_derivative = computeHardeningDerivative();
443 :
444 26912 : GenericReal<is_ad> hardeningVariable = computeHardeningValue(delta_gamma, omega);
445 26912 : GenericReal<is_ad> sy =
446 26912 : _hardening_constant * std::pow(hardeningVariable + 1.0e-30, _hardening_exponent) +
447 26912 : _yield_stress;
448 26912 : GenericReal<is_ad> sy_alpha = _hardening_derivative;
449 :
450 : GenericReal<is_ad> omega_gamma;
451 : GenericReal<is_ad> sy_gamma;
452 :
453 26912 : computeDeltaDerivatives(delta_gamma, stress_sigma, sy_alpha, omega, omega_gamma, sy_gamma);
454 53824 : GenericReal<is_ad> residual_derivative = 1 / omega * (sy_gamma - 1 / omega * omega_gamma * sy);
455 :
456 26912 : return residual_derivative;
457 : }
458 :
459 : template <bool is_ad>
460 : void
461 26912 : 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 26912 : omega_gamma = 0.0;
470 26912 : sy_gamma = 0.0;
471 :
472 26912 : GenericDenseVector<is_ad> K_deltaGamma(6);
473 :
474 26912 : if (_local_cylindrical_csys)
475 21120 : omega = computeOmega(delta_gamma, _sigma_tilde_rotated[_qp]);
476 : else
477 5792 : omega = computeOmega(delta_gamma, _sigma_tilde[_qp]);
478 :
479 26912 : GenericDenseVector<is_ad> K(6);
480 188384 : for (unsigned int i = 0; i < 6; i++)
481 645888 : K(i) = _b_eigenvalues[_qp](i) / Utility::pow<2>(1 + delta_gamma * _b_eigenvalues[_qp](i));
482 :
483 188384 : for (unsigned int i = 0; i < 6; i++)
484 161472 : K_deltaGamma(i) =
485 645888 : -2.0 * _b_eigenvalues[_qp](i) * K(i) / (1 + _b_eigenvalues[_qp](i) * delta_gamma);
486 :
487 188384 : for (unsigned int i = 0; i < 6; i++)
488 : {
489 161472 : if (_local_cylindrical_csys)
490 253440 : omega_gamma += K_deltaGamma(i) * _sigma_tilde_rotated[_qp](i) * _sigma_tilde_rotated[_qp](i);
491 : else
492 69504 : omega_gamma += K_deltaGamma(i) * _sigma_tilde[_qp](i) * _sigma_tilde[_qp](i);
493 : }
494 :
495 53824 : omega_gamma /= 4.0 * omega;
496 53824 : sy_gamma = 2.0 * sy_alpha * (omega + delta_gamma * omega_gamma);
497 26912 : }
498 :
499 : template <bool is_ad>
500 : GenericReal<is_ad>
501 137168 : HillElastoPlasticityStressUpdateTempl<is_ad>::computeHardeningValue(
502 : const GenericReal<is_ad> & delta_gamma, const GenericReal<is_ad> & omega)
503 : {
504 274336 : return _hardening_variable_old[_qp] + 2.0 * delta_gamma * omega;
505 : }
506 :
507 : template <bool is_ad>
508 : Real
509 26912 : HillElastoPlasticityStressUpdateTempl<is_ad>::computeHardeningDerivative()
510 : {
511 26912 : return _hardening_constant * _hardening_exponent *
512 26912 : MetaPhysicL::raw_value(
513 26912 : std::pow(_hardening_variable[_qp] + 1.0e-30, _hardening_exponent - 1));
514 : }
515 :
516 : template <bool is_ad>
517 : void
518 50688 : 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 50688 : GenericDenseVector<is_ad> inelasticStrainIncrement_vector(6);
526 50688 : GenericDenseVector<is_ad> hill_stress(6);
527 50688 : GenericDenseVector<is_ad> stress_vector(6);
528 :
529 50688 : if (_local_cylindrical_csys)
530 : {
531 0 : GenericRankTwoTensor<is_ad> stress_rotated(stress);
532 21120 : stress_rotated.rotate(_rotation_matrix_transpose[_qp]);
533 :
534 21120 : RankTwoScalarTools::RankTwoTensorToVoigtVector<is_ad>(stress_rotated, stress_vector);
535 : }
536 : else
537 : {
538 29568 : RankTwoScalarTools::RankTwoTensorToVoigtVector<is_ad>(stress, stress_vector);
539 : }
540 :
541 50688 : _hill_tensor[_qp].vector_mult(hill_stress, stress_vector);
542 50688 : hill_stress.scale(delta_gamma);
543 : inelasticStrainIncrement_vector = hill_stress;
544 :
545 50688 : inelasticStrainIncrement(0, 0) = inelasticStrainIncrement_vector(0);
546 50688 : inelasticStrainIncrement(1, 1) = inelasticStrainIncrement_vector(1);
547 50688 : inelasticStrainIncrement(2, 2) = inelasticStrainIncrement_vector(2);
548 50688 : inelasticStrainIncrement(0, 1) = inelasticStrainIncrement(1, 0) =
549 50688 : inelasticStrainIncrement_vector(3) / 2.0;
550 50688 : inelasticStrainIncrement(1, 2) = inelasticStrainIncrement(2, 1) =
551 50688 : inelasticStrainIncrement_vector(4) / 2.0;
552 50688 : inelasticStrainIncrement(0, 2) = inelasticStrainIncrement(2, 0) =
553 50688 : inelasticStrainIncrement_vector(5) / 2.0;
554 :
555 50688 : if (_local_cylindrical_csys)
556 21120 : inelasticStrainIncrement.rotate(_rotation_matrix[_qp]);
557 :
558 : // Calculate equivalent plastic strain
559 50688 : GenericDenseVector<is_ad> Mepsilon(6);
560 50688 : _hill_tensor[_qp].vector_mult(Mepsilon, inelasticStrainIncrement_vector);
561 50688 : GenericReal<is_ad> eq_plastic_strain_inc = Mepsilon.dot(inelasticStrainIncrement_vector);
562 :
563 50688 : if (eq_plastic_strain_inc > 0.0)
564 23320 : eq_plastic_strain_inc = std::sqrt(eq_plastic_strain_inc);
565 :
566 101376 : _effective_inelastic_strain[_qp] = _effective_inelastic_strain_old[_qp] + eq_plastic_strain_inc;
567 :
568 50688 : AnisotropicReturnPlasticityStressUpdateBaseTempl<is_ad>::computeStrainFinalize(
569 : inelasticStrainIncrement, stress, stress_dev, delta_gamma);
570 50688 : }
571 :
572 : template <bool is_ad>
573 : void
574 60024 : 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 60024 : if (_yield_condition <= 0.0)
587 31304 : return;
588 28720 : GenericDenseMatrix<is_ad> inv_matrix(6, 6);
589 28720 : inv_matrix.zero();
590 :
591 201040 : for (unsigned int i = 0; i < 6; i++)
592 689280 : inv_matrix(i, i) = 1 / (1 + delta_gamma * _b_eigenvalues[_qp](i));
593 :
594 28720 : _alpha_matrix[_qp].right_multiply(inv_matrix);
595 :
596 28720 : GenericDenseVector<is_ad> stress_output(6);
597 28720 : if (_local_cylindrical_csys)
598 26520 : _alpha_matrix[_qp].vector_mult(stress_output, _sigma_tilde_rotated[_qp]);
599 : else
600 2200 : _alpha_matrix[_qp].vector_mult(stress_output, _sigma_tilde[_qp]);
601 :
602 28720 : RankTwoScalarTools::VoigtVectorToRankTwoTensor<is_ad>(stress_output, stress_new);
603 :
604 28720 : if (_local_cylindrical_csys)
605 26520 : stress_new.rotate(_rotation_matrix[_qp]);
606 28720 : }
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>;
|
