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