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 "ComputeMultipleCrystalPlasticityStress.h"
11 :
12 : #include "CrystalPlasticityStressUpdateBase.h"
13 : #include "libmesh/utility.h"
14 : #include "Conversion.h"
15 : #include "MooseException.h"
16 :
17 : registerMooseObject("TensorMechanicsApp", ComputeMultipleCrystalPlasticityStress);
18 :
19 : InputParameters
20 585 : ComputeMultipleCrystalPlasticityStress::validParams()
21 : {
22 585 : InputParameters params = ComputeFiniteStrainElasticStress::validParams();
23 :
24 585 : params.addClassDescription(
25 : "Crystal Plasticity base class: handles the Newton iteration over the stress residual and "
26 : "calculates the Jacobian based on constitutive laws from multiple material classes "
27 : "that are inherited from CrystalPlasticityStressUpdateBase");
28 1170 : params.addParam<std::string>(
29 : "base_name",
30 : "Optional parameter that allows the user to define multiple mechanics material systems on "
31 : "the same block, i.e. for multiple phases");
32 :
33 1170 : params.addRequiredParam<std::vector<MaterialName>>(
34 : "crystal_plasticity_models",
35 : "The material objects to use to calculate crystal plasticity stress and strains.");
36 1170 : params.addParam<std::vector<MaterialName>>(
37 : "eigenstrain_names", {}, "The material objects to calculate eigenstrains.");
38 1170 : params.addParam<MooseEnum>("tan_mod_type",
39 1755 : MooseEnum("exact none", "none"),
40 : "Type of tangent moduli for preconditioner: default elastic");
41 1170 : params.addParam<Real>("rtol", 1e-6, "Constitutive stress residual relative tolerance");
42 1170 : params.addParam<Real>("abs_tol", 1e-6, "Constitutive stress residual absolute tolerance");
43 1170 : params.addParam<unsigned int>("maxiter", 100, "Maximum number of iterations for stress update");
44 1170 : params.addParam<unsigned int>(
45 1170 : "maxiter_state_variable", 100, "Maximum number of iterations for state variable update");
46 1170 : params.addParam<unsigned int>(
47 1170 : "maximum_substep_iteration", 1, "Maximum number of substep iteration");
48 1170 : params.addParam<bool>("use_line_search", false, "Use line search in constitutive update");
49 1170 : params.addParam<Real>("min_line_search_step_size", 0.01, "Minimum line search step size");
50 1170 : params.addParam<Real>("line_search_tol", 0.5, "Line search bisection method tolerance");
51 1170 : params.addParam<unsigned int>(
52 1170 : "line_search_maxiter", 20, "Line search bisection method maximum number of iteration");
53 1170 : params.addParam<MooseEnum>("line_search_method",
54 1755 : MooseEnum("CUT_HALF BISECTION", "CUT_HALF"),
55 : "The method used in line search");
56 1170 : params.addParam<bool>(
57 : "print_state_variable_convergence_error_messages",
58 1170 : false,
59 : "Whether or not to print warning messages from the crystal plasticity specific convergence "
60 : "checks on the stress measure and general constitutive model quantinties.");
61 585 : return params;
62 0 : }
63 :
64 438 : ComputeMultipleCrystalPlasticityStress::ComputeMultipleCrystalPlasticityStress(
65 438 : const InputParameters & parameters)
66 : : ComputeFiniteStrainElasticStress(parameters),
67 438 : _num_models(getParam<std::vector<MaterialName>>("crystal_plasticity_models").size()),
68 876 : _num_eigenstrains(getParam<std::vector<MaterialName>>("eigenstrain_names").size()),
69 912 : _base_name(isParamValid("base_name") ? getParam<std::string>("base_name") + "_" : ""),
70 876 : _elasticity_tensor(getMaterialPropertyByName<RankFourTensor>(_base_name + "elasticity_tensor")),
71 876 : _rtol(getParam<Real>("rtol")),
72 876 : _abs_tol(getParam<Real>("abs_tol")),
73 876 : _maxiter(getParam<unsigned int>("maxiter")),
74 876 : _maxiterg(getParam<unsigned int>("maxiter_state_variable")),
75 876 : _tan_mod_type(getParam<MooseEnum>("tan_mod_type").getEnum<TangentModuliType>()),
76 876 : _max_substep_iter(getParam<unsigned int>("maximum_substep_iteration")),
77 876 : _use_line_search(getParam<bool>("use_line_search")),
78 876 : _min_line_search_step_size(getParam<Real>("min_line_search_step_size")),
79 876 : _line_search_tolerance(getParam<Real>("line_search_tol")),
80 876 : _line_search_max_iterations(getParam<unsigned int>("line_search_maxiter")),
81 876 : _line_search_method(getParam<MooseEnum>("line_search_method").getEnum<LineSearchMethod>()),
82 438 : _plastic_deformation_gradient(declareProperty<RankTwoTensor>("plastic_deformation_gradient")),
83 438 : _plastic_deformation_gradient_old(
84 438 : getMaterialPropertyOld<RankTwoTensor>("plastic_deformation_gradient")),
85 438 : _eigenstrain_deformation_gradient(
86 438 : _num_eigenstrains ? &declareProperty<RankTwoTensor>("eigenstrain_deformation_gradient")
87 : : nullptr),
88 438 : _eigenstrain_deformation_gradient_old(
89 438 : _num_eigenstrains
90 507 : ? &getMaterialPropertyOld<RankTwoTensor>("eigenstrain_deformation_gradient")
91 : : nullptr),
92 876 : _deformation_gradient(getMaterialProperty<RankTwoTensor>(_base_name + "deformation_gradient")),
93 438 : _deformation_gradient_old(
94 438 : getMaterialPropertyOld<RankTwoTensor>(_base_name + "deformation_gradient")),
95 438 : _pk2(declareProperty<RankTwoTensor>("second_piola_kirchhoff_stress")),
96 876 : _pk2_old(getMaterialPropertyOld<RankTwoTensor>("second_piola_kirchhoff_stress")),
97 438 : _total_lagrangian_strain(
98 438 : declareProperty<RankTwoTensor>("total_lagrangian_strain")), // Lagrangian strain
99 438 : _updated_rotation(declareProperty<RankTwoTensor>("updated_rotation")),
100 438 : _crysrot(getMaterialProperty<RankTwoTensor>(
101 438 : _base_name + "crysrot")), // defined in the elasticity tensor classes for crystal plasticity
102 1752 : _print_convergence_message(getParam<bool>("print_state_variable_convergence_error_messages"))
103 : {
104 438 : _convergence_failed = false;
105 438 : }
106 :
107 : void
108 11392 : ComputeMultipleCrystalPlasticityStress::initQpStatefulProperties()
109 : {
110 11392 : _plastic_deformation_gradient[_qp].zero();
111 11392 : _plastic_deformation_gradient[_qp].addIa(1.0);
112 :
113 11392 : if (_num_eigenstrains)
114 : {
115 3200 : (*_eigenstrain_deformation_gradient)[_qp].zero();
116 3200 : (*_eigenstrain_deformation_gradient)[_qp].addIa(1.0);
117 : }
118 : else
119 : {
120 : // set to identity if no eigenstrain is added
121 : _inverse_eigenstrain_deformation_grad.zero();
122 8192 : _inverse_eigenstrain_deformation_grad.addIa(1.0);
123 : }
124 :
125 11392 : _pk2[_qp].zero();
126 :
127 11392 : _total_lagrangian_strain[_qp].zero();
128 :
129 11392 : _updated_rotation[_qp].zero();
130 11392 : _updated_rotation[_qp].addIa(1.0);
131 :
132 24064 : for (unsigned int i = 0; i < _num_models; ++i)
133 : {
134 12672 : _models[i]->setQp(_qp);
135 12672 : _models[i]->initQpStatefulProperties();
136 : }
137 :
138 14624 : for (unsigned int i = 0; i < _num_eigenstrains; ++i)
139 : {
140 3232 : _eigenstrains[i]->setQp(_qp);
141 3232 : _eigenstrains[i]->initQpStatefulProperties();
142 : }
143 11392 : }
144 :
145 : void
146 396 : ComputeMultipleCrystalPlasticityStress::initialSetup()
147 : {
148 : // get crystal plasticity models
149 : std::vector<MaterialName> model_names =
150 792 : getParam<std::vector<MaterialName>>("crystal_plasticity_models");
151 :
152 837 : for (unsigned int i = 0; i < _num_models; ++i)
153 : {
154 : CrystalPlasticityStressUpdateBase * model =
155 441 : dynamic_cast<CrystalPlasticityStressUpdateBase *>(&getMaterialByName(model_names[i]));
156 :
157 441 : if (model)
158 : {
159 441 : _models.push_back(model);
160 : // TODO: check to make sure that the material model is compatible with this class
161 : }
162 : else
163 0 : mooseError("Model " + model_names[i] +
164 : " is not compatible with ComputeMultipleCrystalPlasticityStress");
165 : }
166 :
167 : // get crystal plasticity eigenstrains
168 : std::vector<MaterialName> eigenstrain_names =
169 792 : getParam<std::vector<MaterialName>>("eigenstrain_names");
170 :
171 474 : for (unsigned int i = 0; i < _num_eigenstrains; ++i)
172 : {
173 : ComputeCrystalPlasticityEigenstrainBase * eigenstrain =
174 78 : dynamic_cast<ComputeCrystalPlasticityEigenstrainBase *>(
175 78 : &getMaterialByName(eigenstrain_names[i]));
176 :
177 78 : if (eigenstrain)
178 78 : _eigenstrains.push_back(eigenstrain);
179 : else
180 0 : mooseError("Eigenstrain" + eigenstrain_names[i] +
181 : " is not compatible with ComputeMultipleCrystalPlasticityStress");
182 : }
183 396 : }
184 :
185 : void
186 846917 : ComputeMultipleCrystalPlasticityStress::computeQpStress()
187 : {
188 1784074 : for (unsigned int i = 0; i < _num_models; ++i)
189 937157 : _models[i]->setQp(_qp);
190 :
191 1146415 : for (unsigned int i = 0; i < _num_eigenstrains; ++i)
192 299498 : _eigenstrains[i]->setQp(_qp);
193 :
194 846917 : updateStress(_stress[_qp], _Jacobian_mult[_qp]); // This is NOT the exact jacobian
195 846785 : }
196 :
197 : void
198 846917 : ComputeMultipleCrystalPlasticityStress::updateStress(RankTwoTensor & cauchy_stress,
199 : RankFourTensor & jacobian_mult)
200 : {
201 : // Does not support face/boundary material property calculation
202 846917 : if (isBoundaryMaterial())
203 : return;
204 :
205 : // Initialize substepping variables
206 : unsigned int substep_iter = 1;
207 : unsigned int num_substep = 1;
208 :
209 845637 : _temporary_deformation_gradient_old = _deformation_gradient_old[_qp];
210 845637 : if (_temporary_deformation_gradient_old.det() == 0)
211 0 : _temporary_deformation_gradient_old.addIa(1.0);
212 :
213 845637 : _delta_deformation_gradient = _deformation_gradient[_qp] - _temporary_deformation_gradient_old;
214 :
215 : // Loop through all models and calculate the schmid tensor for the current state of the crystal
216 : // lattice
217 : // Not sure if we should pass in the updated or the original rotation here
218 : // If not, then we should not need to compute the flow direction every iteration here
219 1781514 : for (unsigned int i = 0; i < _num_models; ++i)
220 935877 : _models[i]->calculateFlowDirection(_crysrot[_qp]);
221 :
222 : do
223 : {
224 935027 : _convergence_failed = false;
225 935027 : preSolveQp();
226 :
227 935027 : _substep_dt = _dt / num_substep;
228 1960294 : for (unsigned int i = 0; i < _num_models; ++i)
229 1025267 : _models[i]->setSubstepDt(_substep_dt);
230 :
231 : // calculate F^{eigen} only when we have eigenstrain
232 935027 : if (_num_eigenstrains)
233 341410 : calculateEigenstrainDeformationGrad();
234 :
235 2007074 : for (unsigned int istep = 0; istep < num_substep; ++istep)
236 : {
237 : _temporary_deformation_gradient =
238 1161569 : (static_cast<Real>(istep) + 1) / num_substep * _delta_deformation_gradient;
239 1161569 : _temporary_deformation_gradient += _temporary_deformation_gradient_old;
240 :
241 1161569 : solveQp();
242 :
243 1161565 : if (_convergence_failed)
244 : {
245 89508 : if (_print_convergence_message)
246 0 : mooseWarning(
247 : "The crystal plasticity constitutive model has failed to converge. Increasing "
248 : "the number of substeps.");
249 :
250 89508 : substep_iter++;
251 89508 : num_substep *= 2;
252 89508 : break;
253 : }
254 : }
255 :
256 935013 : if (substep_iter > _max_substep_iter && _convergence_failed)
257 118 : mooseException("ComputeMultipleCrystalPlasticityStress: Constitutive failure");
258 934895 : } while (_convergence_failed);
259 :
260 845505 : postSolveQp(cauchy_stress, jacobian_mult);
261 : }
262 :
263 : void
264 935027 : ComputeMultipleCrystalPlasticityStress::preSolveQp()
265 : {
266 1960294 : for (unsigned int i = 0; i < _num_models; ++i)
267 1025267 : _models[i]->setInitialConstitutiveVariableValues();
268 :
269 935027 : _pk2[_qp] = _pk2_old[_qp];
270 935027 : _inverse_plastic_deformation_grad_old = _plastic_deformation_gradient_old[_qp].inverse();
271 935027 : }
272 :
273 : void
274 1161569 : ComputeMultipleCrystalPlasticityStress::solveQp()
275 : {
276 2413378 : for (unsigned int i = 0; i < _num_models; ++i)
277 : {
278 1251809 : _models[i]->setSubstepConstitutiveVariableValues();
279 1251809 : _models[i]->calculateSlipResistance();
280 : }
281 :
282 1161569 : _inverse_plastic_deformation_grad = _inverse_plastic_deformation_grad_old;
283 :
284 1161569 : solveStateVariables();
285 1161565 : if (_convergence_failed)
286 : return; // pop back up and take a smaller substep
287 :
288 2234354 : for (unsigned int i = 0; i < _num_models; ++i)
289 1162297 : _models[i]->updateSubstepConstitutiveVariableValues();
290 :
291 : // save off the old F^{p} inverse now that have converged on the stress and state variables
292 1072057 : _inverse_plastic_deformation_grad_old = _inverse_plastic_deformation_grad;
293 : }
294 :
295 : void
296 845505 : ComputeMultipleCrystalPlasticityStress::postSolveQp(RankTwoTensor & cauchy_stress,
297 : RankFourTensor & jacobian_mult)
298 : {
299 845505 : cauchy_stress = _elastic_deformation_gradient * _pk2[_qp] *
300 845505 : _elastic_deformation_gradient.transpose() / _elastic_deformation_gradient.det();
301 :
302 845505 : calcTangentModuli(jacobian_mult);
303 :
304 845505 : _total_lagrangian_strain[_qp] =
305 845505 : _deformation_gradient[_qp].transpose() * _deformation_gradient[_qp] -
306 845505 : RankTwoTensor::Identity();
307 845505 : _total_lagrangian_strain[_qp] = _total_lagrangian_strain[_qp] * 0.5;
308 :
309 : // Calculate crystal rotation to track separately
310 845505 : RankTwoTensor rot;
311 845505 : _elastic_deformation_gradient.getRUDecompositionRotation(rot);
312 845505 : _updated_rotation[_qp] = rot * _crysrot[_qp];
313 845505 : }
314 :
315 : void
316 1161569 : ComputeMultipleCrystalPlasticityStress::solveStateVariables()
317 : {
318 : unsigned int iteration;
319 : bool iter_flag = true;
320 :
321 : iteration = 0;
322 : // Check for slip system resistance update tolerance
323 : do
324 : {
325 1498836 : solveStress();
326 1498833 : if (_convergence_failed)
327 : return;
328 :
329 1409937 : _plastic_deformation_gradient[_qp] =
330 1409937 : _inverse_plastic_deformation_grad.inverse(); // the postSoveStress
331 :
332 : // Update slip system resistance and state variable after the stress has been finalized
333 : // We loop through all the models for each calculation
334 : // in order to make sure that when coupling appears, the state variables are updated based on
335 : // the same components
336 2973408 : for (unsigned int i = 0; i < _num_models; ++i)
337 1563471 : _models[i]->cacheStateVariablesBeforeUpdate();
338 :
339 2973408 : for (unsigned int i = 0; i < _num_models; ++i)
340 1563471 : _models[i]->calculateStateVariableEvolutionRateComponent();
341 :
342 2973407 : for (unsigned int i = 0; i < _num_models; ++i)
343 1563471 : if (!_models[i]->updateStateVariables())
344 58 : _convergence_failed = true;
345 :
346 2973406 : for (unsigned int i = 0; i < _num_models; ++i)
347 1563470 : _models[i]->calculateSlipResistance();
348 :
349 1409936 : if (_convergence_failed)
350 : return;
351 :
352 2572943 : for (unsigned int i = 0; i < _num_models; ++i)
353 : {
354 : // iter_flag = false, stop iteration if all values are converged and good to go
355 : // iter_flag = true, continue iteration if any value is not converged
356 1500886 : if (!_models[i]->areConstitutiveStateVariablesConverged())
357 : {
358 : iter_flag = true;
359 : break;
360 : }
361 : else
362 : iter_flag = false; // iter_flag = false, stop iteration only when all models returns true
363 : }
364 :
365 1409878 : if (iter_flag)
366 : {
367 337821 : if (_print_convergence_message)
368 0 : mooseWarning("ComputeMultipleCrystalPlasticityStress: State variables (or the system "
369 : "resistance) did not converge at element ",
370 0 : _current_elem->id(),
371 : " and qp ",
372 0 : _qp,
373 : "\n");
374 : }
375 1072057 : iteration++;
376 337821 : } while (iter_flag && iteration < _maxiterg);
377 :
378 1072611 : if (iteration == _maxiterg)
379 : {
380 554 : if (_print_convergence_message)
381 0 : mooseWarning(
382 : "ComputeMultipleCrystalPlasticityStress: Hardness Integration error. Reached the "
383 : "maximum number of iterations to solve for the state variables at element ",
384 0 : _current_elem->id(),
385 : " and qp ",
386 0 : _qp,
387 : "\n");
388 :
389 554 : _convergence_failed = true;
390 : }
391 : }
392 :
393 : void
394 1498836 : ComputeMultipleCrystalPlasticityStress::solveStress()
395 : {
396 : unsigned int iteration = 0;
397 1498836 : RankTwoTensor dpk2;
398 : Real rnorm, rnorm0, rnorm_prev;
399 :
400 : // Calculate stress residual
401 1498836 : calculateResidualAndJacobian();
402 1498836 : if (_convergence_failed)
403 : {
404 0 : if (_print_convergence_message)
405 0 : mooseWarning("ComputeMultipleCrystalPlasticityStress: the slip increment exceeds tolerance "
406 : "at element ",
407 0 : _current_elem->id(),
408 : " and Gauss point ",
409 0 : _qp);
410 :
411 0 : return;
412 : }
413 :
414 1498836 : rnorm = _residual_tensor.L2norm();
415 1498836 : rnorm0 = rnorm;
416 :
417 : // Check for stress residual tolerance; different from user object version which
418 : // compares the absolute tolerance of only the original rnorm value
419 6470911 : while (rnorm > _rtol * rnorm0 && rnorm > _abs_tol && iteration < _maxiter)
420 : {
421 : // Calculate stress increment
422 5057938 : dpk2 = -_jacobian.invSymm() * _residual_tensor;
423 5057938 : _pk2[_qp] = _pk2[_qp] + dpk2;
424 :
425 5057938 : calculateResidualAndJacobian();
426 :
427 5057936 : if (_convergence_failed)
428 : {
429 85861 : if (_print_convergence_message)
430 0 : mooseWarning("ComputeMultipleCrystalPlasticityStress: the slip increment exceeds tolerance "
431 : "at element ",
432 0 : _current_elem->id(),
433 : " and Gauss point ",
434 0 : _qp);
435 :
436 85861 : return;
437 : }
438 :
439 4972075 : rnorm_prev = rnorm;
440 4972075 : rnorm = _residual_tensor.L2norm();
441 :
442 4972075 : if (_use_line_search && rnorm > rnorm_prev && !lineSearchUpdate(rnorm_prev, dpk2))
443 : {
444 0 : if (_print_convergence_message)
445 0 : mooseWarning("ComputeMultipleCrystalPlasticityStress: Failed with line search");
446 :
447 0 : _convergence_failed = true;
448 0 : return;
449 : }
450 :
451 4972075 : if (_use_line_search)
452 339648 : rnorm = _residual_tensor.L2norm();
453 :
454 4972075 : iteration++;
455 : }
456 :
457 1412973 : if (iteration >= _maxiter)
458 : {
459 3036 : if (_print_convergence_message)
460 1 : mooseWarning("ComputeMultipleCrystalPlasticityStress: Stress Integration error rmax = ",
461 : rnorm,
462 : " and the tolerance is ",
463 1 : _rtol * rnorm0,
464 : " when the rnorm0 value is ",
465 : rnorm0,
466 : " for element ",
467 0 : _current_elem->id(),
468 : " and qp ",
469 1 : _qp);
470 :
471 3035 : _convergence_failed = true;
472 : }
473 : }
474 :
475 : // Calculates stress residual equation and jacobian
476 : void
477 6556774 : ComputeMultipleCrystalPlasticityStress::calculateResidualAndJacobian()
478 : {
479 6556774 : calculateResidual();
480 6556772 : if (_convergence_failed)
481 : return;
482 6470911 : calculateJacobian();
483 : }
484 :
485 : void
486 6582694 : ComputeMultipleCrystalPlasticityStress::calculateResidual()
487 : {
488 6582694 : RankTwoTensor ce, elastic_strain, ce_pk2, equivalent_slip_increment_per_model,
489 6582694 : equivalent_slip_increment, pk2_new;
490 :
491 : equivalent_slip_increment.zero();
492 :
493 : // calculate slip rate in order to compute F^{p-1}
494 13609042 : for (unsigned int i = 0; i < _num_models; ++i)
495 : {
496 : equivalent_slip_increment_per_model.zero();
497 :
498 : // calculate shear stress with consideration of contribution from other physics
499 7112211 : _models[i]->calculateShearStress(
500 7112211 : _pk2[_qp], _inverse_eigenstrain_deformation_grad, _num_eigenstrains);
501 :
502 7112211 : _convergence_failed = !_models[i]->calculateSlipRate();
503 :
504 7112209 : if (_convergence_failed)
505 : return;
506 :
507 7026348 : _models[i]->calculateEquivalentSlipIncrement(equivalent_slip_increment_per_model);
508 7026348 : equivalent_slip_increment += equivalent_slip_increment_per_model;
509 : }
510 :
511 : RankTwoTensor residual_equivalent_slip_increment =
512 6496831 : RankTwoTensor::Identity() - equivalent_slip_increment;
513 : _inverse_plastic_deformation_grad =
514 6496831 : _inverse_plastic_deformation_grad_old * residual_equivalent_slip_increment;
515 :
516 12993662 : _elastic_deformation_gradient = _temporary_deformation_gradient *
517 12993662 : _inverse_eigenstrain_deformation_grad *
518 6496831 : _inverse_plastic_deformation_grad;
519 :
520 6496831 : ce = _elastic_deformation_gradient.transpose() * _elastic_deformation_gradient;
521 :
522 6496831 : elastic_strain = ce - RankTwoTensor::Identity();
523 6496831 : elastic_strain *= 0.5;
524 :
525 6496831 : pk2_new = _elasticity_tensor[_qp] * elastic_strain;
526 6496831 : _residual_tensor = _pk2[_qp] - pk2_new;
527 : }
528 :
529 : void
530 6470911 : ComputeMultipleCrystalPlasticityStress::calculateJacobian()
531 : {
532 : // may not need to cache the dfpinvdpk2 here. need to double check
533 6470911 : RankFourTensor dfedfpinv, deedfe, dfpinvdpk2, dfpinvdpk2_per_model;
534 :
535 6470911 : RankTwoTensor ffeiginv = _temporary_deformation_gradient * _inverse_eigenstrain_deformation_grad;
536 :
537 25883644 : for (const auto i : make_range(Moose::dim))
538 77650932 : for (const auto j : make_range(Moose::dim))
539 232952796 : for (const auto k : make_range(Moose::dim))
540 174714597 : dfedfpinv(i, j, k, j) = ffeiginv(i, k);
541 :
542 25883644 : for (const auto i : make_range(Moose::dim))
543 77650932 : for (const auto j : make_range(Moose::dim))
544 232952796 : for (const auto k : make_range(Moose::dim))
545 : {
546 174714597 : deedfe(i, j, k, i) = deedfe(i, j, k, i) + _elastic_deformation_gradient(k, j) * 0.5;
547 174714597 : deedfe(i, j, k, j) = deedfe(i, j, k, j) + _elastic_deformation_gradient(k, i) * 0.5;
548 : }
549 :
550 13471339 : for (unsigned int i = 0; i < _num_models; ++i)
551 : {
552 7000428 : _models[i]->calculateTotalPlasticDeformationGradientDerivative(
553 : dfpinvdpk2_per_model,
554 7000428 : _inverse_plastic_deformation_grad_old,
555 7000428 : _inverse_eigenstrain_deformation_grad,
556 7000428 : _num_eigenstrains);
557 7000428 : dfpinvdpk2 += dfpinvdpk2_per_model;
558 : }
559 :
560 : _jacobian =
561 12941822 : RankFourTensor::IdentityFour() - (_elasticity_tensor[_qp] * deedfe * dfedfpinv * dfpinvdpk2);
562 6470911 : }
563 :
564 : void
565 845505 : ComputeMultipleCrystalPlasticityStress::calcTangentModuli(RankFourTensor & jacobian_mult)
566 : {
567 845505 : switch (_tan_mod_type)
568 : {
569 845505 : case TangentModuliType::EXACT:
570 845505 : elastoPlasticTangentModuli(jacobian_mult);
571 845505 : break;
572 0 : default:
573 0 : elasticTangentModuli(jacobian_mult);
574 : }
575 845505 : }
576 :
577 : void
578 845505 : ComputeMultipleCrystalPlasticityStress::elastoPlasticTangentModuli(RankFourTensor & jacobian_mult)
579 : {
580 845505 : RankFourTensor tan_mod;
581 845505 : RankTwoTensor pk2fet, fepk2, feiginvfpinv;
582 845505 : RankFourTensor deedfe, dsigdpk2dfe, dfedf;
583 :
584 : // Fill in the matrix stiffness material property
585 3382020 : for (const auto i : make_range(Moose::dim))
586 10146060 : for (const auto j : make_range(Moose::dim))
587 30438180 : for (const auto k : make_range(Moose::dim))
588 : {
589 22828635 : deedfe(i, j, k, i) = deedfe(i, j, k, i) + _elastic_deformation_gradient(k, j) * 0.5;
590 22828635 : deedfe(i, j, k, j) = deedfe(i, j, k, j) + _elastic_deformation_gradient(k, i) * 0.5;
591 : }
592 :
593 : usingTensorIndices(i_, j_, k_, l_);
594 845505 : dsigdpk2dfe = _elastic_deformation_gradient.times<i_, k_, j_, l_>(_elastic_deformation_gradient) *
595 845505 : _elasticity_tensor[_qp] * deedfe;
596 :
597 845505 : pk2fet = _pk2[_qp] * _elastic_deformation_gradient.transpose();
598 845505 : fepk2 = _elastic_deformation_gradient * _pk2[_qp];
599 :
600 3382020 : for (const auto i : make_range(Moose::dim))
601 10146060 : for (const auto j : make_range(Moose::dim))
602 30438180 : for (const auto l : make_range(Moose::dim))
603 : {
604 22828635 : tan_mod(i, j, i, l) += pk2fet(l, j);
605 22828635 : tan_mod(i, j, j, l) += fepk2(i, l);
606 : }
607 :
608 845505 : tan_mod += dsigdpk2dfe;
609 :
610 845505 : const auto je = _elastic_deformation_gradient.det();
611 845505 : if (je > 0.0)
612 845505 : tan_mod /= je;
613 :
614 845505 : feiginvfpinv = _inverse_eigenstrain_deformation_grad * _inverse_plastic_deformation_grad;
615 3382020 : for (const auto i : make_range(Moose::dim))
616 10146060 : for (const auto j : make_range(Moose::dim))
617 30438180 : for (const auto l : make_range(Moose::dim))
618 22828635 : dfedf(i, j, i, l) = feiginvfpinv(l, j);
619 :
620 845505 : jacobian_mult = tan_mod * dfedf;
621 845505 : }
622 :
623 : void
624 0 : ComputeMultipleCrystalPlasticityStress::elasticTangentModuli(RankFourTensor & jacobian_mult)
625 : {
626 : // update jacobian_mult
627 0 : jacobian_mult = _elasticity_tensor[_qp];
628 0 : }
629 :
630 : bool
631 25920 : ComputeMultipleCrystalPlasticityStress::lineSearchUpdate(const Real & rnorm_prev,
632 : const RankTwoTensor & dpk2)
633 : {
634 25920 : if (_line_search_method == LineSearchMethod::CutHalf)
635 : {
636 : Real rnorm;
637 25920 : Real step = 1.0;
638 :
639 : do
640 : {
641 25920 : _pk2[_qp] = _pk2[_qp] - step * dpk2;
642 25920 : step /= 2.0;
643 25920 : _pk2[_qp] = _pk2[_qp] + step * dpk2;
644 :
645 25920 : calculateResidual();
646 25920 : rnorm = _residual_tensor.L2norm();
647 25920 : } while (rnorm > rnorm_prev && step > _min_line_search_step_size);
648 :
649 : // has norm improved or is the step still above minumum search step size?
650 25920 : return (rnorm <= rnorm_prev || step > _min_line_search_step_size);
651 : }
652 0 : else if (_line_search_method == LineSearchMethod::Bisection)
653 : {
654 : unsigned int count = 0;
655 : Real step_a = 0.0;
656 : Real step_b = 1.0;
657 0 : Real step = 1.0;
658 : Real s_m = 1000.0;
659 : Real rnorm = 1000.0;
660 :
661 0 : calculateResidual();
662 0 : auto s_b = _residual_tensor.doubleContraction(dpk2);
663 0 : const auto rnorm1 = _residual_tensor.L2norm();
664 0 : _pk2[_qp] = _pk2[_qp] - dpk2;
665 0 : calculateResidual();
666 0 : auto s_a = _residual_tensor.doubleContraction(dpk2);
667 0 : const auto rnorm0 = _residual_tensor.L2norm();
668 0 : _pk2[_qp] = _pk2[_qp] + dpk2;
669 :
670 0 : if ((rnorm1 / rnorm0) < _line_search_tolerance || s_a * s_b > 0)
671 : {
672 0 : calculateResidual();
673 0 : return true;
674 : }
675 :
676 0 : while ((rnorm / rnorm0) > _line_search_tolerance && count < _line_search_max_iterations)
677 : {
678 0 : _pk2[_qp] = _pk2[_qp] - step * dpk2;
679 0 : step = 0.5 * (step_b + step_a);
680 0 : _pk2[_qp] = _pk2[_qp] + step * dpk2;
681 0 : calculateResidual();
682 0 : s_m = _residual_tensor.doubleContraction(dpk2);
683 0 : rnorm = _residual_tensor.L2norm();
684 :
685 0 : if (s_m * s_a < 0.0)
686 : {
687 : step_b = step;
688 : s_b = s_m;
689 : }
690 0 : if (s_m * s_b < 0.0)
691 : {
692 : step_a = step;
693 : s_a = s_m;
694 : }
695 0 : count++;
696 : }
697 :
698 : // below tolerance and max iterations?
699 0 : return ((rnorm / rnorm0) < _line_search_tolerance && count < _line_search_max_iterations);
700 : }
701 : else
702 0 : mooseError("Line search method is not provided.");
703 : }
704 :
705 : void
706 341410 : ComputeMultipleCrystalPlasticityStress::calculateEigenstrainDeformationGrad()
707 : {
708 : _inverse_eigenstrain_deformation_grad.zero();
709 341410 : _inverse_eigenstrain_deformation_grad.addIa(1.0);
710 :
711 729114 : for (unsigned int i = 0; i < _num_eigenstrains; ++i)
712 : {
713 387714 : _eigenstrains[i]->setSubstepDt(_substep_dt);
714 387714 : _eigenstrains[i]->computeQpProperties();
715 387704 : _inverse_eigenstrain_deformation_grad *= _eigenstrains[i]->getDeformationGradientInverse();
716 : }
717 341400 : (*_eigenstrain_deformation_gradient)[_qp] = _inverse_eigenstrain_deformation_grad.inverse();
718 341400 : }
|
