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 64 : HillElastoPlasticityStressUpdateTempl<is_ad>::validParams()
20 : {
21 64 : InputParameters params = AnisotropicReturnPlasticityStressUpdateBaseTempl<is_ad>::validParams();
22 64 : 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 128 : params.addRequiredParam<Real>("hardening_constant",
28 : "Hardening constant (H) for anisotropic plasticity");
29 128 : params.addParam<Real>(
30 128 : "hardening_exponent", 1.0, "Hardening exponent (n) for anisotropic plasticity");
31 128 : params.addRequiredParam<Real>("yield_stress",
32 : "Yield stress (constant value) for anisotropic plasticity");
33 128 : params.addParam<bool>(
34 : "local_cylindrical_csys",
35 128 : false,
36 : "Compute inelastic strain increment in local cylindrical coordinate system");
37 :
38 128 : params.addParam<MooseEnum>("axis",
39 192 : 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 64 : return params;
44 0 : }
45 :
46 : template <bool is_ad>
47 48 : HillElastoPlasticityStressUpdateTempl<is_ad>::HillElastoPlasticityStressUpdateTempl(
48 : const InputParameters & parameters)
49 : : AnisotropicReturnPlasticityStressUpdateBaseTempl<is_ad>(parameters),
50 96 : _qsigma(0.0),
51 48 : _eigenvalues_hill(6),
52 48 : _eigenvectors_hill(6, 6),
53 48 : _anisotropic_elastic_tensor(6, 6),
54 48 : _elasticity_tensor_name(this->_base_name + "elasticity_tensor"),
55 48 : _elasticity_tensor(
56 48 : this->template getGenericMaterialProperty<RankFourTensor, is_ad>(_elasticity_tensor_name)),
57 96 : _hardening_constant(this->template getParam<Real>("hardening_constant")),
58 96 : _hardening_exponent(this->template getParam<Real>("hardening_exponent")),
59 96 : _hardening_variable(this->template declareGenericProperty<Real, is_ad>(this->_base_name +
60 : "hardening_variable")),
61 48 : _hardening_variable_old(
62 48 : this->template getMaterialPropertyOld<Real>(this->_base_name + "hardening_variable")),
63 96 : _elasticity_eigenvectors(this->template declareGenericProperty<DenseMatrix<Real>, is_ad>(
64 : this->_base_name + "elasticity_eigenvectors")),
65 96 : _elasticity_eigenvalues(this->template declareGenericProperty<DenseVector<Real>, is_ad>(
66 : this->_base_name + "elasticity_eigenvalues")),
67 96 : _b_eigenvectors(this->template declareGenericProperty<DenseMatrix<Real>, is_ad>(
68 : this->_base_name + "b_eigenvectors")),
69 96 : _b_eigenvalues(this->template declareGenericProperty<DenseVector<Real>, is_ad>(
70 : this->_base_name + "b_eigenvalues")),
71 96 : _alpha_matrix(this->template declareGenericProperty<DenseMatrix<Real>, is_ad>(this->_base_name +
72 : "alpha_matrix")),
73 96 : _sigma_tilde(this->template declareGenericProperty<DenseVector<Real>, is_ad>(this->_base_name +
74 : "sigma_tilde")),
75 96 : _sigma_tilde_rotated(this->template declareGenericProperty<DenseVector<Real>, is_ad>(
76 : this->_base_name + "sigma_tilde_rotated")),
77 48 : _hardening_derivative(0.0),
78 48 : _yield_condition(1.0),
79 96 : _yield_stress(this->template getParam<Real>("yield_stress")),
80 96 : _hill_tensor(this->template getMaterialPropertyByName<DenseMatrix<Real>>(this->_base_name +
81 : "hill_tensor")),
82 48 : _rotation_matrix(
83 48 : this->template declareProperty<RankTwoTensor>(this->_base_name + "rotation_matrix")),
84 48 : _rotation_matrix_transpose(this->template declareProperty<RankTwoTensor>(
85 : this->_base_name + "rotation_matrix_transpose")),
86 48 : _rotation_matrix_old(
87 48 : this->template getMaterialPropertyOld<RankTwoTensor>(this->_base_name + "rotation_matrix")),
88 96 : _rotation_matrix_transpose_old(this->template getMaterialPropertyOld<RankTwoTensor>(
89 : this->_base_name + "rotation_matrix_transpose")),
90 96 : _local_cylindrical_csys(this->template getParam<bool>("local_cylindrical_csys")),
91 144 : _axis(this->template getParam<MooseEnum>("axis").template getEnum<Axis>())
92 : {
93 120 : 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 60 : 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 48 : }
108 :
109 : template <bool is_ad>
110 : void
111 4080 : HillElastoPlasticityStressUpdateTempl<is_ad>::initQpStatefulProperties()
112 : {
113 : RealVectorValue normal_vector(0, 0, 0);
114 : RealVectorValue axial_vector(0, 0, 0);
115 :
116 4080 : 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 4080 : 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 4080 : AnisotropicReturnPlasticityStressUpdateBaseTempl<is_ad>::initQpStatefulProperties();
170 4080 : }
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 86840 : 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 86840 : _hardening_variable[_qp] = _hardening_variable_old[_qp];
196 86840 : _plasticity_strain[_qp] = _plasticity_strain_old[_qp];
197 86840 : _effective_inelastic_strain[_qp] = _effective_inelastic_strain_old[_qp];
198 :
199 86840 : _elasticity_eigenvectors[_qp].resize(6, 6);
200 86840 : _elasticity_eigenvalues[_qp].resize(6);
201 :
202 86840 : _b_eigenvectors[_qp].resize(6, 6);
203 86840 : _b_eigenvalues[_qp].resize(6);
204 :
205 86840 : _alpha_matrix[_qp].resize(6, 6);
206 86840 : _sigma_tilde[_qp].resize(6);
207 86840 : _sigma_tilde_rotated[_qp].resize(6);
208 :
209 : // Algebra needed for the generalized return mapping of anisotropic elastoplasticity
210 86840 : computeElasticityTensorEigenDecomposition();
211 :
212 86840 : _yield_condition = 1.0; // Some positive value
213 86840 : _yield_condition = -computeResidual(stress_dev, stress, 0.0);
214 86840 : }
215 :
216 : template <bool is_ad>
217 : void
218 86840 : HillElastoPlasticityStressUpdateTempl<is_ad>::computeElasticityTensorEigenDecomposition()
219 : {
220 : // Helper method to compute qp matrices required for the elasto-plasticity algorithm.
221 : using std::sqrt;
222 :
223 86840 : GenericRankFourTensor<is_ad> elasticity_tensor = _elasticity_tensor[_qp];
224 :
225 86840 : if (_local_cylindrical_csys)
226 : {
227 39480 : elasticity_tensor.rotate(_rotation_matrix_transpose[_qp]);
228 : }
229 :
230 86840 : ElasticityTensorTools::toVoigtNotation<is_ad>(_anisotropic_elastic_tensor, elasticity_tensor);
231 :
232 : const unsigned int dimension = _anisotropic_elastic_tensor.n();
233 :
234 : AnisotropyMatrixReal A = AnisotropyMatrixReal::Zero();
235 607880 : for (unsigned int index_i = 0; index_i < dimension; index_i++)
236 3647280 : for (unsigned int index_j = 0; index_j < dimension; index_j++)
237 3126240 : A(index_i, index_j) = MetaPhysicL::raw_value(_anisotropic_elastic_tensor(index_i, index_j));
238 :
239 : Eigen::SelfAdjointEigenSolver<AnisotropyMatrixReal> es(A);
240 :
241 : auto lambda = es.eigenvalues();
242 : auto v = es.eigenvectors();
243 :
244 607880 : for (unsigned int index_i = 0; index_i < dimension; index_i++)
245 521040 : _elasticity_eigenvalues[_qp](index_i) = lambda(index_i);
246 :
247 607880 : for (unsigned int index_i = 0; index_i < dimension; index_i++)
248 3647280 : for (unsigned int index_j = 0; index_j < dimension; index_j++)
249 3126240 : _elasticity_eigenvectors[_qp](index_i, index_j) = v(index_i, index_j);
250 :
251 : // Compute sqrt(Delta_c) * QcT * A * Qc * sqrt(Delta_c)
252 86840 : GenericDenseMatrix<is_ad> sqrt_Delta(6, 6);
253 607880 : for (unsigned int i = 0; i < 6; i++)
254 : {
255 521040 : sqrt_Delta(i, i) = sqrt(_elasticity_eigenvalues[_qp](i));
256 : }
257 :
258 86840 : GenericDenseMatrix<is_ad> eigenvectors_elasticity_transpose(6, 6);
259 86840 : _elasticity_eigenvectors[_qp].get_transpose(eigenvectors_elasticity_transpose);
260 :
261 86840 : GenericDenseMatrix<is_ad> b_matrix(6, 6);
262 86840 : b_matrix = _hill_tensor[_qp];
263 :
264 : // Right multiply by matrix of eigenvectors transpose
265 86840 : b_matrix.right_multiply(_elasticity_eigenvectors[_qp]);
266 86840 : b_matrix.right_multiply(sqrt_Delta);
267 86840 : b_matrix.left_multiply(eigenvectors_elasticity_transpose);
268 86840 : b_matrix.left_multiply(sqrt_Delta);
269 :
270 : AnisotropyMatrixReal B_eigen = AnisotropyMatrixReal::Zero();
271 607880 : for (unsigned int index_i = 0; index_i < dimension; index_i++)
272 3647280 : for (unsigned int index_j = 0; index_j < dimension; index_j++)
273 3126240 : B_eigen(index_i, index_j) = MetaPhysicL::raw_value(b_matrix(index_i, index_j));
274 :
275 86840 : if (isBlockDiagonal(B_eigen))
276 : {
277 3840 : Eigen::SelfAdjointEigenSolver<AnisotropyMatrixRealBlock> es_b_left(B_eigen.block<3, 3>(0, 0));
278 3840 : Eigen::SelfAdjointEigenSolver<AnisotropyMatrixRealBlock> es_b_right(B_eigen.block<3, 3>(3, 3));
279 :
280 : auto lambda_b_left = es_b_left.eigenvalues();
281 : auto v_b_left = es_b_left.eigenvectors();
282 :
283 : auto lambda_b_right = es_b_right.eigenvalues();
284 : auto v_b_right = es_b_right.eigenvectors();
285 :
286 3840 : _b_eigenvalues[_qp](0) = lambda_b_left(0);
287 3840 : _b_eigenvalues[_qp](1) = lambda_b_left(1);
288 3840 : _b_eigenvalues[_qp](2) = lambda_b_left(2);
289 3840 : _b_eigenvalues[_qp](3) = lambda_b_right(0);
290 3840 : _b_eigenvalues[_qp](4) = lambda_b_right(1);
291 3840 : _b_eigenvalues[_qp](5) = lambda_b_right(2);
292 :
293 3840 : _b_eigenvectors[_qp](0, 0) = v_b_left(0, 0);
294 3840 : _b_eigenvectors[_qp](0, 1) = v_b_left(0, 1);
295 3840 : _b_eigenvectors[_qp](0, 2) = v_b_left(0, 2);
296 3840 : _b_eigenvectors[_qp](1, 0) = v_b_left(1, 0);
297 3840 : _b_eigenvectors[_qp](1, 1) = v_b_left(1, 1);
298 3840 : _b_eigenvectors[_qp](1, 2) = v_b_left(1, 2);
299 3840 : _b_eigenvectors[_qp](2, 0) = v_b_left(2, 0);
300 3840 : _b_eigenvectors[_qp](2, 1) = v_b_left(2, 1);
301 3840 : _b_eigenvectors[_qp](2, 2) = v_b_left(2, 2);
302 3840 : _b_eigenvectors[_qp](3, 3) = v_b_right(0, 0);
303 3840 : _b_eigenvectors[_qp](3, 4) = v_b_right(0, 1);
304 3840 : _b_eigenvectors[_qp](3, 5) = v_b_right(0, 2);
305 3840 : _b_eigenvectors[_qp](4, 3) = v_b_right(1, 0);
306 3840 : _b_eigenvectors[_qp](4, 4) = v_b_right(1, 1);
307 3840 : _b_eigenvectors[_qp](4, 5) = v_b_right(1, 2);
308 3840 : _b_eigenvectors[_qp](5, 3) = v_b_right(2, 0);
309 3840 : _b_eigenvectors[_qp](5, 4) = v_b_right(2, 1);
310 3840 : _b_eigenvectors[_qp](5, 5) = v_b_right(2, 2);
311 : }
312 : else
313 : {
314 : Eigen::SelfAdjointEigenSolver<AnisotropyMatrixReal> es_b(B_eigen);
315 :
316 : auto lambda_b = es_b.eigenvalues();
317 : auto v_b = es_b.eigenvectors();
318 581000 : for (unsigned int index_i = 0; index_i < dimension; index_i++)
319 498000 : _b_eigenvalues[_qp](index_i) = lambda_b(index_i);
320 :
321 581000 : for (unsigned int index_i = 0; index_i < dimension; index_i++)
322 3486000 : for (unsigned int index_j = 0; index_j < dimension; index_j++)
323 2988000 : _b_eigenvectors[_qp](index_i, index_j) = v_b(index_i, index_j);
324 : }
325 :
326 86840 : _alpha_matrix[_qp] = sqrt_Delta;
327 86840 : _alpha_matrix[_qp].right_multiply(_b_eigenvectors[_qp]);
328 86840 : _alpha_matrix[_qp].left_multiply(_elasticity_eigenvectors[_qp]);
329 86840 : }
330 :
331 : template <bool is_ad>
332 : GenericReal<is_ad>
333 255092 : HillElastoPlasticityStressUpdateTempl<is_ad>::computeOmega(
334 : const GenericReal<is_ad> & delta_gamma, const GenericDenseVector<is_ad> & /*stress_trial*/)
335 : {
336 255092 : GenericDenseVector<is_ad> K(6);
337 255092 : GenericReal<is_ad> omega = 0.0;
338 :
339 1785644 : for (unsigned int i = 0; i < 6; i++)
340 : {
341 6122208 : K(i) = _b_eigenvalues[_qp](i) / Utility::pow<2>(1 + delta_gamma * _b_eigenvalues[_qp](i));
342 :
343 1530552 : if (_local_cylindrical_csys)
344 2109600 : omega += K(i) * _sigma_tilde_rotated[_qp](i) * _sigma_tilde_rotated[_qp](i);
345 : else
346 951504 : omega += K(i) * _sigma_tilde[_qp](i) * _sigma_tilde[_qp](i);
347 : }
348 255092 : omega *= 0.5;
349 :
350 255092 : if (omega == 0.0)
351 9240 : omega = 1.0e-36;
352 :
353 : using std::sqrt;
354 510184 : return sqrt(omega);
355 : }
356 :
357 : template <bool is_ad>
358 : Real
359 121120 : HillElastoPlasticityStressUpdateTempl<is_ad>::computeReferenceResidual(
360 : const GenericDenseVector<is_ad> & /*effective_trial_stress*/,
361 : const GenericDenseVector<is_ad> & /*stress_new*/,
362 : const GenericReal<is_ad> & /*residual*/,
363 : const GenericReal<is_ad> & /*scalar_effective_inelastic_strain*/)
364 : {
365 : // Equation is already normalized.
366 121120 : return 1.0;
367 : }
368 :
369 : template <bool is_ad>
370 : GenericReal<is_ad>
371 207960 : HillElastoPlasticityStressUpdateTempl<is_ad>::computeResidual(
372 : const GenericDenseVector<is_ad> & /*stress_dev*/,
373 : const GenericDenseVector<is_ad> & stress_sigma,
374 : const GenericReal<is_ad> & delta_gamma)
375 : {
376 : using std::pow;
377 :
378 : // If in elastic regime, just return
379 207960 : if (_yield_condition <= 0.0)
380 39908 : return 0.0;
381 :
382 : GenericReal<is_ad> omega;
383 :
384 168052 : if (_local_cylindrical_csys)
385 : {
386 : // Get stress_tilde_rotated
387 107640 : GenericRankTwoTensor<is_ad> stress;
388 107640 : RankTwoScalarTools::VoigtVectorToRankTwoTensor<is_ad>(stress_sigma, stress);
389 :
390 107640 : stress.rotate(_rotation_matrix_transpose[_qp]);
391 :
392 107640 : GenericDenseVector<is_ad> stress_sigma_rotated(6);
393 107640 : RankTwoScalarTools::RankTwoTensorToVoigtVector<is_ad>(stress, stress_sigma_rotated);
394 :
395 107640 : GenericDenseVector<is_ad> stress_tilde_rotated(6);
396 107640 : GenericDenseMatrix<is_ad> alpha_temp_rotated(_alpha_matrix[_qp]);
397 107640 : alpha_temp_rotated.lu_solve(stress_sigma_rotated, stress_tilde_rotated);
398 :
399 : // Material property used in computeStressFinalize
400 107640 : _sigma_tilde_rotated[_qp] = stress_tilde_rotated;
401 107640 : omega = computeOmega(delta_gamma, stress_tilde_rotated);
402 107640 : }
403 :
404 : else
405 : {
406 : // Get stress_tilde
407 60412 : GenericDenseVector<is_ad> stress_tilde(6);
408 60412 : GenericDenseMatrix<is_ad> alpha_temp(_alpha_matrix[_qp]);
409 60412 : alpha_temp.lu_solve(stress_sigma, stress_tilde);
410 :
411 : // Material property used in computeStressFinalize
412 60412 : _sigma_tilde[_qp] = stress_tilde;
413 60412 : omega = computeOmega(delta_gamma, stress_tilde);
414 60412 : }
415 :
416 : // Hardening variable is \alpha isotropic hardening for now.
417 168052 : _hardening_variable[_qp] = computeHardeningValue(delta_gamma, omega);
418 :
419 : // A small value of 1.0e-30 is added to the hardening variable to avoid numerical issues
420 : // related to hardening variable becoming negative early in the iteration leading to non-
421 : // convergence. Note that pow(x,n) requires x to be positive if n is less than 1.
422 :
423 168052 : GenericReal<is_ad> s_y =
424 168052 : _hardening_constant * pow(_hardening_variable[_qp] + 1.0e-30, _hardening_exponent) +
425 168052 : _yield_stress;
426 :
427 168052 : GenericReal<is_ad> residual = 0.0;
428 168052 : residual = s_y / omega - 1.0;
429 :
430 168052 : return residual;
431 : }
432 :
433 : template <bool is_ad>
434 : GenericReal<is_ad>
435 43520 : HillElastoPlasticityStressUpdateTempl<is_ad>::computeDerivative(
436 : const GenericDenseVector<is_ad> & /*stress_dev*/,
437 : const GenericDenseVector<is_ad> & stress_sigma,
438 : const GenericReal<is_ad> & delta_gamma)
439 : {
440 : using std::pow;
441 :
442 : // If in elastic regime, return unit derivative
443 43520 : if (_yield_condition <= 0.0)
444 0 : return 1.0;
445 :
446 43520 : GenericReal<is_ad> omega = computeOmega(delta_gamma, stress_sigma);
447 43520 : _hardening_derivative = computeHardeningDerivative();
448 :
449 43520 : GenericReal<is_ad> hardeningVariable = computeHardeningValue(delta_gamma, omega);
450 43520 : GenericReal<is_ad> sy =
451 87040 : _hardening_constant * pow(hardeningVariable + 1.0e-30, _hardening_exponent) + _yield_stress;
452 43520 : GenericReal<is_ad> sy_alpha = _hardening_derivative;
453 :
454 : GenericReal<is_ad> omega_gamma;
455 : GenericReal<is_ad> sy_gamma;
456 :
457 43520 : computeDeltaDerivatives(delta_gamma, stress_sigma, sy_alpha, omega, omega_gamma, sy_gamma);
458 87040 : GenericReal<is_ad> residual_derivative = 1 / omega * (sy_gamma - 1 / omega * omega_gamma * sy);
459 :
460 43520 : return residual_derivative;
461 : }
462 :
463 : template <bool is_ad>
464 : void
465 43520 : HillElastoPlasticityStressUpdateTempl<is_ad>::computeDeltaDerivatives(
466 : const GenericReal<is_ad> & delta_gamma,
467 : const GenericDenseVector<is_ad> & /*stress_trial*/,
468 : const GenericReal<is_ad> & sy_alpha,
469 : GenericReal<is_ad> & omega,
470 : GenericReal<is_ad> & omega_gamma,
471 : GenericReal<is_ad> & sy_gamma)
472 : {
473 43520 : omega_gamma = 0.0;
474 43520 : sy_gamma = 0.0;
475 :
476 43520 : GenericDenseVector<is_ad> K_deltaGamma(6);
477 :
478 43520 : if (_local_cylindrical_csys)
479 34080 : omega = computeOmega(delta_gamma, _sigma_tilde_rotated[_qp]);
480 : else
481 9440 : omega = computeOmega(delta_gamma, _sigma_tilde[_qp]);
482 :
483 43520 : GenericDenseVector<is_ad> K(6);
484 304640 : for (unsigned int i = 0; i < 6; i++)
485 1044480 : K(i) = _b_eigenvalues[_qp](i) / Utility::pow<2>(1 + delta_gamma * _b_eigenvalues[_qp](i));
486 :
487 304640 : for (unsigned int i = 0; i < 6; i++)
488 261120 : K_deltaGamma(i) =
489 1044480 : -2.0 * _b_eigenvalues[_qp](i) * K(i) / (1 + _b_eigenvalues[_qp](i) * delta_gamma);
490 :
491 304640 : for (unsigned int i = 0; i < 6; i++)
492 : {
493 261120 : if (_local_cylindrical_csys)
494 408960 : omega_gamma += K_deltaGamma(i) * _sigma_tilde_rotated[_qp](i) * _sigma_tilde_rotated[_qp](i);
495 : else
496 113280 : omega_gamma += K_deltaGamma(i) * _sigma_tilde[_qp](i) * _sigma_tilde[_qp](i);
497 : }
498 :
499 87040 : omega_gamma /= 4.0 * omega;
500 87040 : sy_gamma = 2.0 * sy_alpha * (omega + delta_gamma * omega_gamma);
501 43520 : }
502 :
503 : template <bool is_ad>
504 : GenericReal<is_ad>
505 211572 : HillElastoPlasticityStressUpdateTempl<is_ad>::computeHardeningValue(
506 : const GenericReal<is_ad> & delta_gamma, const GenericReal<is_ad> & omega)
507 : {
508 423144 : return _hardening_variable_old[_qp] + 2.0 * delta_gamma * omega;
509 : }
510 :
511 : template <bool is_ad>
512 : Real
513 43520 : HillElastoPlasticityStressUpdateTempl<is_ad>::computeHardeningDerivative()
514 : {
515 : using std::pow;
516 43520 : return _hardening_constant * _hardening_exponent *
517 43520 : MetaPhysicL::raw_value(pow(_hardening_variable[_qp] + 1.0e-30, _hardening_exponent - 1));
518 : }
519 :
520 : template <bool is_ad>
521 : void
522 77600 : HillElastoPlasticityStressUpdateTempl<is_ad>::computeStrainFinalize(
523 : GenericRankTwoTensor<is_ad> & inelasticStrainIncrement,
524 : const GenericRankTwoTensor<is_ad> & stress,
525 : const GenericDenseVector<is_ad> & stress_dev,
526 : const GenericReal<is_ad> & delta_gamma)
527 : {
528 : using std::sqrt;
529 :
530 : // e^P = delta_gamma * hill_tensor * stress
531 77600 : GenericDenseVector<is_ad> inelasticStrainIncrement_vector(6);
532 77600 : GenericDenseVector<is_ad> hill_stress(6);
533 77600 : GenericDenseVector<is_ad> stress_vector(6);
534 :
535 77600 : if (_local_cylindrical_csys)
536 : {
537 0 : GenericRankTwoTensor<is_ad> stress_rotated(stress);
538 34080 : stress_rotated.rotate(_rotation_matrix_transpose[_qp]);
539 :
540 34080 : RankTwoScalarTools::RankTwoTensorToVoigtVector<is_ad>(stress_rotated, stress_vector);
541 : }
542 : else
543 : {
544 43520 : RankTwoScalarTools::RankTwoTensorToVoigtVector<is_ad>(stress, stress_vector);
545 : }
546 :
547 77600 : _hill_tensor[_qp].vector_mult(hill_stress, stress_vector);
548 77600 : hill_stress.scale(delta_gamma);
549 : inelasticStrainIncrement_vector = hill_stress;
550 :
551 77600 : inelasticStrainIncrement(0, 0) = inelasticStrainIncrement_vector(0);
552 77600 : inelasticStrainIncrement(1, 1) = inelasticStrainIncrement_vector(1);
553 77600 : inelasticStrainIncrement(2, 2) = inelasticStrainIncrement_vector(2);
554 77600 : inelasticStrainIncrement(0, 1) = inelasticStrainIncrement(1, 0) =
555 77600 : inelasticStrainIncrement_vector(3) / 2.0;
556 77600 : inelasticStrainIncrement(1, 2) = inelasticStrainIncrement(2, 1) =
557 77600 : inelasticStrainIncrement_vector(4) / 2.0;
558 77600 : inelasticStrainIncrement(0, 2) = inelasticStrainIncrement(2, 0) =
559 77600 : inelasticStrainIncrement_vector(5) / 2.0;
560 :
561 77600 : if (_local_cylindrical_csys)
562 34080 : inelasticStrainIncrement.rotate(_rotation_matrix[_qp]);
563 :
564 : // Calculate equivalent plastic strain
565 77600 : GenericDenseVector<is_ad> Mepsilon(6);
566 77600 : _hill_tensor[_qp].vector_mult(Mepsilon, inelasticStrainIncrement_vector);
567 77600 : GenericReal<is_ad> eq_plastic_strain_inc = Mepsilon.dot(inelasticStrainIncrement_vector);
568 :
569 77600 : if (eq_plastic_strain_inc > 0.0)
570 37692 : eq_plastic_strain_inc = sqrt(eq_plastic_strain_inc);
571 :
572 155200 : _effective_inelastic_strain[_qp] = _effective_inelastic_strain_old[_qp] + eq_plastic_strain_inc;
573 :
574 77600 : AnisotropicReturnPlasticityStressUpdateBaseTempl<is_ad>::computeStrainFinalize(
575 : inelasticStrainIncrement, stress, stress_dev, delta_gamma);
576 77600 : }
577 :
578 : template <bool is_ad>
579 : void
580 86840 : HillElastoPlasticityStressUpdateTempl<is_ad>::computeStressFinalize(
581 : const GenericRankTwoTensor<is_ad> & /*plastic_strain_increment*/,
582 : const GenericReal<is_ad> & delta_gamma,
583 : GenericRankTwoTensor<is_ad> & stress_new,
584 : const GenericDenseVector<is_ad> & /*stress_dev*/,
585 : const GenericRankTwoTensor<is_ad> & /*stress_old*/,
586 : const GenericRankFourTensor<is_ad> & /*elasticity_tensor*/)
587 : {
588 : // Need to compute this iteration's stress tensor based on the scalar variable
589 : // For deviatoric
590 : // sigma(n+1) = {Alpha [I + delta_gamma*Delta_b]^(-1) A^-1} sigma(trial)
591 :
592 86840 : if (_yield_condition <= 0.0)
593 43748 : return;
594 43092 : GenericDenseMatrix<is_ad> inv_matrix(6, 6);
595 43092 : inv_matrix.zero();
596 :
597 301644 : for (unsigned int i = 0; i < 6; i++)
598 1034208 : inv_matrix(i, i) = 1 / (1 + delta_gamma * _b_eigenvalues[_qp](i));
599 :
600 43092 : _alpha_matrix[_qp].right_multiply(inv_matrix);
601 :
602 43092 : GenericDenseVector<is_ad> stress_output(6);
603 43092 : if (_local_cylindrical_csys)
604 39480 : _alpha_matrix[_qp].vector_mult(stress_output, _sigma_tilde_rotated[_qp]);
605 : else
606 3612 : _alpha_matrix[_qp].vector_mult(stress_output, _sigma_tilde[_qp]);
607 :
608 43092 : RankTwoScalarTools::VoigtVectorToRankTwoTensor<is_ad>(stress_output, stress_new);
609 :
610 43092 : if (_local_cylindrical_csys)
611 39480 : stress_new.rotate(_rotation_matrix[_qp]);
612 43092 : }
613 :
614 : template <bool is_ad>
615 : Real
616 0 : HillElastoPlasticityStressUpdateTempl<is_ad>::computeStrainEnergyRateDensity(
617 : const GenericMaterialProperty<RankTwoTensor, is_ad> & /*stress*/,
618 : const GenericMaterialProperty<RankTwoTensor, is_ad> & /*strain_rate*/)
619 : {
620 0 : mooseError("computeStrainEnergyRateDensity not implemented for anisotropic plasticity.");
621 : return 0.0;
622 : }
623 :
624 : template class HillElastoPlasticityStressUpdateTempl<false>;
625 : template class HillElastoPlasticityStressUpdateTempl<true>;
|
