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