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 112 : HillElastoPlasticityStressUpdateTempl<is_ad>::validParams()
20 : {
21 112 : InputParameters params = AnisotropicReturnPlasticityStressUpdateBaseTempl<is_ad>::validParams();
22 112 : 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 224 : params.addRequiredParam<Real>("hardening_constant",
28 : "Hardening constant (H) for anisotropic plasticity");
29 224 : params.addParam<Real>(
30 224 : "hardening_exponent", 1.0, "Hardening exponent (n) for anisotropic plasticity");
31 224 : params.addRequiredParam<Real>("yield_stress",
32 : "Yield stress (constant value) for anisotropic plasticity");
33 224 : params.addParam<bool>(
34 : "local_cylindrical_csys",
35 224 : false,
36 : "Compute inelastic strain increment in local cylindrical coordinate system");
37 :
38 224 : params.addParam<MooseEnum>("axis",
39 336 : 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 112 : return params;
44 0 : }
45 :
46 : template <bool is_ad>
47 84 : HillElastoPlasticityStressUpdateTempl<is_ad>::HillElastoPlasticityStressUpdateTempl(
48 : const InputParameters & parameters)
49 : : AnisotropicReturnPlasticityStressUpdateBaseTempl<is_ad>(parameters),
50 168 : _qsigma(0.0),
51 84 : _eigenvalues_hill(6),
52 84 : _eigenvectors_hill(6, 6),
53 84 : _anisotropic_elastic_tensor(6, 6),
54 84 : _elasticity_tensor_name(this->_base_name + "elasticity_tensor"),
55 84 : _elasticity_tensor(
56 84 : this->template getGenericMaterialProperty<RankFourTensor, is_ad>(_elasticity_tensor_name)),
57 168 : _hardening_constant(this->template getParam<Real>("hardening_constant")),
58 168 : _hardening_exponent(this->template getParam<Real>("hardening_exponent")),
59 168 : _hardening_variable(this->template declareGenericProperty<Real, is_ad>(this->_base_name +
60 : "hardening_variable")),
61 84 : _hardening_variable_old(
62 84 : this->template getMaterialPropertyOld<Real>(this->_base_name + "hardening_variable")),
63 168 : _elasticity_eigenvectors(this->template declareGenericProperty<DenseMatrix<Real>, is_ad>(
64 : this->_base_name + "elasticity_eigenvectors")),
65 168 : _elasticity_eigenvalues(this->template declareGenericProperty<DenseVector<Real>, is_ad>(
66 : this->_base_name + "elasticity_eigenvalues")),
67 168 : _b_eigenvectors(this->template declareGenericProperty<DenseMatrix<Real>, is_ad>(
68 : this->_base_name + "b_eigenvectors")),
69 168 : _b_eigenvalues(this->template declareGenericProperty<DenseVector<Real>, is_ad>(
70 : this->_base_name + "b_eigenvalues")),
71 168 : _alpha_matrix(this->template declareGenericProperty<DenseMatrix<Real>, is_ad>(this->_base_name +
72 : "alpha_matrix")),
73 168 : _sigma_tilde(this->template declareGenericProperty<DenseVector<Real>, is_ad>(this->_base_name +
74 : "sigma_tilde")),
75 168 : _sigma_tilde_rotated(this->template declareGenericProperty<DenseVector<Real>, is_ad>(
76 : this->_base_name + "sigma_tilde_rotated")),
77 84 : _hardening_derivative(0.0),
78 84 : _yield_condition(1.0),
79 168 : _yield_stress(this->template getParam<Real>("yield_stress")),
80 168 : _hill_tensor(this->template getMaterialPropertyByName<DenseMatrix<Real>>(this->_base_name +
81 : "hill_tensor")),
82 84 : _rotation_matrix(
83 84 : this->template declareProperty<RankTwoTensor>(this->_base_name + "rotation_matrix")),
84 84 : _rotation_matrix_transpose(this->template declareProperty<RankTwoTensor>(
85 : this->_base_name + "rotation_matrix_transpose")),
86 84 : _rotation_matrix_old(
87 84 : this->template getMaterialPropertyOld<RankTwoTensor>(this->_base_name + "rotation_matrix")),
88 168 : _rotation_matrix_transpose_old(this->template getMaterialPropertyOld<RankTwoTensor>(
89 : this->_base_name + "rotation_matrix_transpose")),
90 168 : _local_cylindrical_csys(this->template getParam<bool>("local_cylindrical_csys")),
91 252 : _axis(this->template getParam<MooseEnum>("axis").template getEnum<Axis>())
92 : {
93 210 : 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 105 : 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 84 : }
108 :
109 : template <bool is_ad>
110 : void
111 6800 : HillElastoPlasticityStressUpdateTempl<is_ad>::initQpStatefulProperties()
112 : {
113 : RealVectorValue normal_vector(0, 0, 0);
114 : RealVectorValue axial_vector(0, 0, 0);
115 :
116 6800 : if (_local_cylindrical_csys)
117 : {
118 3600 : if (_axis == Axis::X)
119 : {
120 : normal_vector(0) = 0.0;
121 1200 : normal_vector(1) = _q_point[_qp](1);
122 1200 : normal_vector(2) = _q_point[_qp](2);
123 1200 : axial_vector(0) = 1.0;
124 : }
125 2400 : else if (_axis == Axis::Y)
126 : {
127 1200 : normal_vector(0) = _q_point[_qp](0);
128 : normal_vector(1) = 0.0;
129 1200 : normal_vector(2) = _q_point[_qp](2);
130 1200 : axial_vector(1) = 1.0;
131 : }
132 1200 : else if (_axis == Axis::Z)
133 : {
134 1200 : normal_vector(0) = _q_point[_qp](0);
135 1200 : normal_vector(1) = _q_point[_qp](1);
136 : normal_vector(2) = 0.0;
137 1200 : axial_vector(2) = 1.0;
138 : }
139 : else
140 0 : mooseError("\nInvalid value for axis parameter provided in HillElastoPlasticityStressUpdate");
141 : }
142 :
143 6800 : if (_local_cylindrical_csys)
144 : {
145 3600 : Real nv_norm = normal_vector.norm();
146 3600 : Real av_norm = axial_vector.norm();
147 :
148 3600 : 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 3600 : 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 3600 : RankTwoScalarTools::setRotationMatrix(
165 3600 : normal_vector, axial_vector, _rotation_matrix[_qp], false);
166 3600 : _rotation_matrix_transpose[_qp] = _rotation_matrix[_qp].transpose();
167 : }
168 :
169 6800 : AnisotropicReturnPlasticityStressUpdateBaseTempl<is_ad>::initQpStatefulProperties();
170 6800 : }
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 135560 : 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 135560 : _hardening_variable[_qp] = _hardening_variable_old[_qp];
196 135560 : _plasticity_strain[_qp] = _plasticity_strain_old[_qp];
197 135560 : _effective_inelastic_strain[_qp] = _effective_inelastic_strain_old[_qp];
198 :
199 135560 : _elasticity_eigenvectors[_qp].resize(6, 6);
200 135560 : _elasticity_eigenvalues[_qp].resize(6);
201 :
202 135560 : _b_eigenvectors[_qp].resize(6, 6);
203 135560 : _b_eigenvalues[_qp].resize(6);
204 :
205 135560 : _alpha_matrix[_qp].resize(6, 6);
206 135560 : _sigma_tilde[_qp].resize(6);
207 135560 : _sigma_tilde_rotated[_qp].resize(6);
208 :
209 : // Algebra needed for the generalized return mapping of anisotropic elastoplasticity
210 135560 : computeElasticityTensorEigenDecomposition();
211 :
212 135560 : _yield_condition = 1.0; // Some positive value
213 135560 : _yield_condition = -computeResidual(stress_dev, stress, 0.0);
214 135560 : }
215 :
216 : template <bool is_ad>
217 : void
218 135560 : HillElastoPlasticityStressUpdateTempl<is_ad>::computeElasticityTensorEigenDecomposition()
219 : {
220 : // Helper method to compute qp matrices required for the elasto-plasticity algorithm.
221 : using std::sqrt;
222 :
223 135560 : GenericRankFourTensor<is_ad> elasticity_tensor = _elasticity_tensor[_qp];
224 :
225 135560 : if (_local_cylindrical_csys)
226 : {
227 59400 : elasticity_tensor.rotate(_rotation_matrix_transpose[_qp]);
228 : }
229 :
230 135560 : 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 948920 : for (unsigned int index_i = 0; index_i < dimension; index_i++)
236 5693520 : for (unsigned int index_j = 0; index_j < dimension; index_j++)
237 4880160 : 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 948920 : for (unsigned int index_i = 0; index_i < dimension; index_i++)
245 813360 : _elasticity_eigenvalues[_qp](index_i) = lambda(index_i);
246 :
247 948920 : for (unsigned int index_i = 0; index_i < dimension; index_i++)
248 5693520 : for (unsigned int index_j = 0; index_j < dimension; index_j++)
249 4880160 : _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 135560 : GenericDenseMatrix<is_ad> sqrt_Delta(6, 6);
253 948920 : for (unsigned int i = 0; i < 6; i++)
254 : {
255 813360 : sqrt_Delta(i, i) = sqrt(_elasticity_eigenvalues[_qp](i));
256 : }
257 :
258 135560 : GenericDenseMatrix<is_ad> eigenvectors_elasticity_transpose(6, 6);
259 135560 : _elasticity_eigenvectors[_qp].get_transpose(eigenvectors_elasticity_transpose);
260 :
261 135560 : GenericDenseMatrix<is_ad> b_matrix(6, 6);
262 135560 : b_matrix = _hill_tensor[_qp];
263 :
264 : // Right multiply by matrix of eigenvectors transpose
265 135560 : b_matrix.right_multiply(_elasticity_eigenvectors[_qp]);
266 135560 : b_matrix.right_multiply(sqrt_Delta);
267 135560 : b_matrix.left_multiply(eigenvectors_elasticity_transpose);
268 135560 : b_matrix.left_multiply(sqrt_Delta);
269 :
270 : AnisotropyMatrixReal B_eigen = AnisotropyMatrixReal::Zero();
271 948920 : for (unsigned int index_i = 0; index_i < dimension; index_i++)
272 5693520 : for (unsigned int index_j = 0; index_j < dimension; index_j++)
273 4880160 : B_eigen(index_i, index_j) = MetaPhysicL::raw_value(b_matrix(index_i, index_j));
274 :
275 135560 : if (isBlockDiagonal(B_eigen))
276 : {
277 6400 : Eigen::SelfAdjointEigenSolver<AnisotropyMatrixRealBlock> es_b_left(B_eigen.block<3, 3>(0, 0));
278 6400 : 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 6400 : _b_eigenvalues[_qp](0) = lambda_b_left(0);
287 6400 : _b_eigenvalues[_qp](1) = lambda_b_left(1);
288 6400 : _b_eigenvalues[_qp](2) = lambda_b_left(2);
289 6400 : _b_eigenvalues[_qp](3) = lambda_b_right(0);
290 6400 : _b_eigenvalues[_qp](4) = lambda_b_right(1);
291 6400 : _b_eigenvalues[_qp](5) = lambda_b_right(2);
292 :
293 6400 : _b_eigenvectors[_qp](0, 0) = v_b_left(0, 0);
294 6400 : _b_eigenvectors[_qp](0, 1) = v_b_left(0, 1);
295 6400 : _b_eigenvectors[_qp](0, 2) = v_b_left(0, 2);
296 6400 : _b_eigenvectors[_qp](1, 0) = v_b_left(1, 0);
297 6400 : _b_eigenvectors[_qp](1, 1) = v_b_left(1, 1);
298 6400 : _b_eigenvectors[_qp](1, 2) = v_b_left(1, 2);
299 6400 : _b_eigenvectors[_qp](2, 0) = v_b_left(2, 0);
300 6400 : _b_eigenvectors[_qp](2, 1) = v_b_left(2, 1);
301 6400 : _b_eigenvectors[_qp](2, 2) = v_b_left(2, 2);
302 6400 : _b_eigenvectors[_qp](3, 3) = v_b_right(0, 0);
303 6400 : _b_eigenvectors[_qp](3, 4) = v_b_right(0, 1);
304 6400 : _b_eigenvectors[_qp](3, 5) = v_b_right(0, 2);
305 6400 : _b_eigenvectors[_qp](4, 3) = v_b_right(1, 0);
306 6400 : _b_eigenvectors[_qp](4, 4) = v_b_right(1, 1);
307 6400 : _b_eigenvectors[_qp](4, 5) = v_b_right(1, 2);
308 6400 : _b_eigenvectors[_qp](5, 3) = v_b_right(2, 0);
309 6400 : _b_eigenvectors[_qp](5, 4) = v_b_right(2, 1);
310 6400 : _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 904120 : for (unsigned int index_i = 0; index_i < dimension; index_i++)
319 774960 : _b_eigenvalues[_qp](index_i) = lambda_b(index_i);
320 :
321 904120 : for (unsigned int index_i = 0; index_i < dimension; index_i++)
322 5424720 : for (unsigned int index_j = 0; index_j < dimension; index_j++)
323 4649760 : _b_eigenvectors[_qp](index_i, index_j) = v_b(index_i, index_j);
324 : }
325 :
326 135560 : _alpha_matrix[_qp] = sqrt_Delta;
327 135560 : _alpha_matrix[_qp].right_multiply(_b_eigenvectors[_qp]);
328 135560 : _alpha_matrix[_qp].left_multiply(_elasticity_eigenvectors[_qp]);
329 135560 : }
330 :
331 : template <bool is_ad>
332 : GenericReal<is_ad>
333 391548 : HillElastoPlasticityStressUpdateTempl<is_ad>::computeOmega(
334 : const GenericReal<is_ad> & delta_gamma, const GenericDenseVector<is_ad> & /*stress_trial*/)
335 : {
336 391548 : GenericDenseVector<is_ad> K(6);
337 391548 : GenericReal<is_ad> omega = 0.0;
338 :
339 2740836 : for (unsigned int i = 0; i < 6; i++)
340 : {
341 9397152 : K(i) = _b_eigenvalues[_qp](i) / Utility::pow<2>(1 + delta_gamma * _b_eigenvalues[_qp](i));
342 :
343 2349288 : if (_local_cylindrical_csys)
344 3166560 : omega += K(i) * _sigma_tilde_rotated[_qp](i) * _sigma_tilde_rotated[_qp](i);
345 : else
346 1532016 : omega += K(i) * _sigma_tilde[_qp](i) * _sigma_tilde[_qp](i);
347 : }
348 391548 : omega *= 0.5;
349 :
350 391548 : if (omega == 0.0)
351 14680 : omega = 1.0e-36;
352 :
353 : using std::sqrt;
354 783096 : return sqrt(omega);
355 : }
356 :
357 : template <bool is_ad>
358 : Real
359 187232 : 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 187232 : return 1.0;
367 : }
368 :
369 : template <bool is_ad>
370 : GenericReal<is_ad>
371 322792 : 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 322792 : if (_yield_condition <= 0.0)
380 63948 : return 0.0;
381 :
382 : GenericReal<is_ad> omega;
383 :
384 258844 : if (_local_cylindrical_csys)
385 : {
386 : // Get stress_tilde_rotated
387 161640 : GenericRankTwoTensor<is_ad> stress;
388 161640 : RankTwoScalarTools::VoigtVectorToRankTwoTensor<is_ad>(stress_sigma, stress);
389 :
390 161640 : stress.rotate(_rotation_matrix_transpose[_qp]);
391 :
392 161640 : GenericDenseVector<is_ad> stress_sigma_rotated(6);
393 161640 : RankTwoScalarTools::RankTwoTensorToVoigtVector<is_ad>(stress, stress_sigma_rotated);
394 :
395 161640 : GenericDenseVector<is_ad> stress_tilde_rotated(6);
396 161640 : GenericDenseMatrix<is_ad> alpha_temp_rotated(_alpha_matrix[_qp]);
397 161640 : alpha_temp_rotated.lu_solve(stress_sigma_rotated, stress_tilde_rotated);
398 :
399 : // Material property used in computeStressFinalize
400 161640 : _sigma_tilde_rotated[_qp] = stress_tilde_rotated;
401 161640 : omega = computeOmega(delta_gamma, stress_tilde_rotated);
402 161640 : }
403 :
404 : else
405 : {
406 : // Get stress_tilde
407 97204 : GenericDenseVector<is_ad> stress_tilde(6);
408 97204 : GenericDenseMatrix<is_ad> alpha_temp(_alpha_matrix[_qp]);
409 97204 : alpha_temp.lu_solve(stress_sigma, stress_tilde);
410 :
411 : // Material property used in computeStressFinalize
412 97204 : _sigma_tilde[_qp] = stress_tilde;
413 97204 : omega = computeOmega(delta_gamma, stress_tilde);
414 97204 : }
415 :
416 : // Hardening variable is \alpha isotropic hardening for now.
417 258844 : _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 258844 : GenericReal<is_ad> s_y =
424 258844 : _hardening_constant * pow(_hardening_variable[_qp] + 1.0e-30, _hardening_exponent) +
425 258844 : _yield_stress;
426 :
427 258844 : GenericReal<is_ad> residual = 0.0;
428 258844 : residual = s_y / omega - 1.0;
429 :
430 258844 : return residual;
431 : }
432 :
433 : template <bool is_ad>
434 : GenericReal<is_ad>
435 66352 : 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 66352 : if (_yield_condition <= 0.0)
444 0 : return 1.0;
445 :
446 66352 : GenericReal<is_ad> omega = computeOmega(delta_gamma, stress_sigma);
447 66352 : _hardening_derivative = computeHardeningDerivative();
448 :
449 66352 : GenericReal<is_ad> hardeningVariable = computeHardeningValue(delta_gamma, omega);
450 66352 : GenericReal<is_ad> sy =
451 132704 : _hardening_constant * pow(hardeningVariable + 1.0e-30, _hardening_exponent) + _yield_stress;
452 66352 : GenericReal<is_ad> sy_alpha = _hardening_derivative;
453 :
454 : GenericReal<is_ad> omega_gamma;
455 : GenericReal<is_ad> sy_gamma;
456 :
457 66352 : computeDeltaDerivatives(delta_gamma, stress_sigma, sy_alpha, omega, omega_gamma, sy_gamma);
458 132704 : GenericReal<is_ad> residual_derivative = 1 / omega * (sy_gamma - 1 / omega * omega_gamma * sy);
459 :
460 66352 : return residual_derivative;
461 : }
462 :
463 : template <bool is_ad>
464 : void
465 66352 : 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 66352 : omega_gamma = 0.0;
474 66352 : sy_gamma = 0.0;
475 :
476 66352 : GenericDenseVector<is_ad> K_deltaGamma(6);
477 :
478 66352 : if (_local_cylindrical_csys)
479 51120 : omega = computeOmega(delta_gamma, _sigma_tilde_rotated[_qp]);
480 : else
481 15232 : omega = computeOmega(delta_gamma, _sigma_tilde[_qp]);
482 :
483 66352 : GenericDenseVector<is_ad> K(6);
484 464464 : for (unsigned int i = 0; i < 6; i++)
485 1592448 : K(i) = _b_eigenvalues[_qp](i) / Utility::pow<2>(1 + delta_gamma * _b_eigenvalues[_qp](i));
486 :
487 464464 : for (unsigned int i = 0; i < 6; i++)
488 398112 : K_deltaGamma(i) =
489 1592448 : -2.0 * _b_eigenvalues[_qp](i) * K(i) / (1 + _b_eigenvalues[_qp](i) * delta_gamma);
490 :
491 464464 : for (unsigned int i = 0; i < 6; i++)
492 : {
493 398112 : if (_local_cylindrical_csys)
494 613440 : omega_gamma += K_deltaGamma(i) * _sigma_tilde_rotated[_qp](i) * _sigma_tilde_rotated[_qp](i);
495 : else
496 182784 : omega_gamma += K_deltaGamma(i) * _sigma_tilde[_qp](i) * _sigma_tilde[_qp](i);
497 : }
498 :
499 132704 : omega_gamma /= 4.0 * omega;
500 132704 : sy_gamma = 2.0 * sy_alpha * (omega + delta_gamma * omega_gamma);
501 66352 : }
502 :
503 : template <bool is_ad>
504 : GenericReal<is_ad>
505 325196 : HillElastoPlasticityStressUpdateTempl<is_ad>::computeHardeningValue(
506 : const GenericReal<is_ad> & delta_gamma, const GenericReal<is_ad> & omega)
507 : {
508 650392 : return _hardening_variable_old[_qp] + 2.0 * delta_gamma * omega;
509 : }
510 :
511 : template <bool is_ad>
512 : Real
513 66352 : HillElastoPlasticityStressUpdateTempl<is_ad>::computeHardeningDerivative()
514 : {
515 : using std::pow;
516 66352 : return _hardening_constant * _hardening_exponent *
517 66352 : MetaPhysicL::raw_value(pow(_hardening_variable[_qp] + 1.0e-30, _hardening_exponent - 1));
518 : }
519 :
520 : template <bool is_ad>
521 : void
522 120880 : 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 120880 : GenericDenseVector<is_ad> inelasticStrainIncrement_vector(6);
532 120880 : GenericDenseVector<is_ad> hill_stress(6);
533 120880 : GenericDenseVector<is_ad> stress_vector(6);
534 :
535 120880 : if (_local_cylindrical_csys)
536 : {
537 0 : GenericRankTwoTensor<is_ad> stress_rotated(stress);
538 51120 : stress_rotated.rotate(_rotation_matrix_transpose[_qp]);
539 :
540 51120 : RankTwoScalarTools::RankTwoTensorToVoigtVector<is_ad>(stress_rotated, stress_vector);
541 : }
542 : else
543 : {
544 69760 : RankTwoScalarTools::RankTwoTensorToVoigtVector<is_ad>(stress, stress_vector);
545 : }
546 :
547 120880 : _hill_tensor[_qp].vector_mult(hill_stress, stress_vector);
548 120880 : hill_stress.scale(delta_gamma);
549 : inelasticStrainIncrement_vector = hill_stress;
550 :
551 120880 : inelasticStrainIncrement(0, 0) = inelasticStrainIncrement_vector(0);
552 120880 : inelasticStrainIncrement(1, 1) = inelasticStrainIncrement_vector(1);
553 120880 : inelasticStrainIncrement(2, 2) = inelasticStrainIncrement_vector(2);
554 120880 : inelasticStrainIncrement(0, 1) = inelasticStrainIncrement(1, 0) =
555 120880 : inelasticStrainIncrement_vector(3) / 2.0;
556 120880 : inelasticStrainIncrement(1, 2) = inelasticStrainIncrement(2, 1) =
557 120880 : inelasticStrainIncrement_vector(4) / 2.0;
558 120880 : inelasticStrainIncrement(0, 2) = inelasticStrainIncrement(2, 0) =
559 120880 : inelasticStrainIncrement_vector(5) / 2.0;
560 :
561 120880 : if (_local_cylindrical_csys)
562 51120 : inelasticStrainIncrement.rotate(_rotation_matrix[_qp]);
563 :
564 : // Calculate equivalent plastic strain
565 120880 : GenericDenseVector<is_ad> Mepsilon(6);
566 120880 : _hill_tensor[_qp].vector_mult(Mepsilon, inelasticStrainIncrement_vector);
567 120880 : GenericReal<is_ad> eq_plastic_strain_inc = Mepsilon.dot(inelasticStrainIncrement_vector);
568 :
569 120880 : if (eq_plastic_strain_inc > 0.0)
570 56932 : eq_plastic_strain_inc = sqrt(eq_plastic_strain_inc);
571 :
572 241760 : _effective_inelastic_strain[_qp] = _effective_inelastic_strain_old[_qp] + eq_plastic_strain_inc;
573 :
574 120880 : AnisotropicReturnPlasticityStressUpdateBaseTempl<is_ad>::computeStrainFinalize(
575 : inelasticStrainIncrement, stress, stress_dev, delta_gamma);
576 120880 : }
577 :
578 : template <bool is_ad>
579 : void
580 135560 : 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 135560 : if (_yield_condition <= 0.0)
593 70348 : return;
594 65212 : GenericDenseMatrix<is_ad> inv_matrix(6, 6);
595 65212 : inv_matrix.zero();
596 :
597 456484 : for (unsigned int i = 0; i < 6; i++)
598 1565088 : inv_matrix(i, i) = 1 / (1 + delta_gamma * _b_eigenvalues[_qp](i));
599 :
600 65212 : _alpha_matrix[_qp].right_multiply(inv_matrix);
601 :
602 65212 : GenericDenseVector<is_ad> stress_output(6);
603 65212 : if (_local_cylindrical_csys)
604 59400 : _alpha_matrix[_qp].vector_mult(stress_output, _sigma_tilde_rotated[_qp]);
605 : else
606 5812 : _alpha_matrix[_qp].vector_mult(stress_output, _sigma_tilde[_qp]);
607 :
608 65212 : RankTwoScalarTools::VoigtVectorToRankTwoTensor<is_ad>(stress_output, stress_new);
609 :
610 65212 : if (_local_cylindrical_csys)
611 59400 : stress_new.rotate(_rotation_matrix[_qp]);
612 65212 : }
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>;
|
