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 1630 : ComputeMultipleCrystalPlasticityStress::validParams()
21 : {
22 1630 : InputParameters params = ComputeFiniteStrainElasticStress::validParams();
23 :
24 1630 : 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 3260 : 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 3260 : params.addRequiredParam<std::vector<MaterialName>>(
34 : "crystal_plasticity_models",
35 : "The material objects to use to calculate crystal plasticity stress and strains.");
36 3260 : params.addParam<std::vector<MaterialName>>(
37 : "eigenstrain_names", {}, "The material objects to calculate eigenstrains.");
38 3260 : params.addParam<MooseEnum>("tan_mod_type",
39 4890 : MooseEnum("exact none", "none"),
40 : "Type of tangent moduli for preconditioner: default elastic");
41 3260 : params.addParam<Real>("rtol", 1e-6, "Constitutive stress residual relative tolerance");
42 3260 : params.addParam<Real>("abs_tol", 1e-6, "Constitutive stress residual absolute tolerance");
43 3260 : params.addParam<unsigned int>("maxiter", 100, "Maximum number of iterations for stress update");
44 3260 : params.addParam<unsigned int>(
45 3260 : "maxiter_state_variable", 100, "Maximum number of iterations for state variable update");
46 3260 : params.addParam<unsigned int>(
47 3260 : "maximum_substep_iteration", 1, "Maximum number of substep iteration");
48 3260 : params.addParam<bool>("use_line_search", false, "Use line search in constitutive update");
49 3260 : params.addParam<Real>("min_line_search_step_size", 0.01, "Minimum line search step size");
50 3260 : params.addParam<Real>("line_search_tol", 0.5, "Line search bisection method tolerance");
51 3260 : params.addParam<unsigned int>(
52 3260 : "line_search_maxiter", 20, "Line search bisection method maximum number of iteration");
53 3260 : params.addParam<MooseEnum>("line_search_method",
54 4890 : MooseEnum("CUT_HALF BISECTION", "CUT_HALF"),
55 : "The method used in line search");
56 3260 : params.addParam<bool>(
57 : "print_state_variable_convergence_error_messages",
58 3260 : 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 1630 : return params;
62 0 : }
63 :
64 1218 : ComputeMultipleCrystalPlasticityStress::ComputeMultipleCrystalPlasticityStress(
65 1218 : const InputParameters & parameters)
66 : : ComputeFiniteStrainElasticStress(parameters),
67 1218 : _num_models(getParam<std::vector<MaterialName>>("crystal_plasticity_models").size()),
68 2436 : _num_eigenstrains(getParam<std::vector<MaterialName>>("eigenstrain_names").size()),
69 2604 : _base_name(isParamValid("base_name") ? getParam<std::string>("base_name") + "_" : ""),
70 2436 : _elasticity_tensor(getMaterialPropertyByName<RankFourTensor>(_base_name + "elasticity_tensor")),
71 2436 : _rtol(getParam<Real>("rtol")),
72 2436 : _abs_tol(getParam<Real>("abs_tol")),
73 2436 : _maxiter(getParam<unsigned int>("maxiter")),
74 2436 : _maxiterg(getParam<unsigned int>("maxiter_state_variable")),
75 2436 : _tan_mod_type(getParam<MooseEnum>("tan_mod_type").getEnum<TangentModuliType>()),
76 2436 : _max_substep_iter(getParam<unsigned int>("maximum_substep_iteration")),
77 2436 : _use_line_search(getParam<bool>("use_line_search")),
78 2436 : _min_line_search_step_size(getParam<Real>("min_line_search_step_size")),
79 2436 : _line_search_tolerance(getParam<Real>("line_search_tol")),
80 2436 : _line_search_max_iterations(getParam<unsigned int>("line_search_maxiter")),
81 2436 : _line_search_method(getParam<MooseEnum>("line_search_method").getEnum<LineSearchMethod>()),
82 1218 : _plastic_deformation_gradient(declareProperty<RankTwoTensor>("plastic_deformation_gradient")),
83 1218 : _plastic_deformation_gradient_old(
84 1218 : getMaterialPropertyOld<RankTwoTensor>("plastic_deformation_gradient")),
85 1218 : _eigenstrain_deformation_gradient(
86 1218 : _num_eigenstrains ? &declareProperty<RankTwoTensor>("eigenstrain_deformation_gradient")
87 : : nullptr),
88 1218 : _eigenstrain_deformation_gradient_old(
89 1218 : _num_eigenstrains
90 1413 : ? &getMaterialPropertyOld<RankTwoTensor>("eigenstrain_deformation_gradient")
91 : : nullptr),
92 2436 : _deformation_gradient(getMaterialProperty<RankTwoTensor>(_base_name + "deformation_gradient")),
93 1218 : _deformation_gradient_old(
94 1218 : getMaterialPropertyOld<RankTwoTensor>(_base_name + "deformation_gradient")),
95 1218 : _pk2(declareProperty<RankTwoTensor>("second_piola_kirchhoff_stress")),
96 2436 : _pk2_old(getMaterialPropertyOld<RankTwoTensor>("second_piola_kirchhoff_stress")),
97 1218 : _total_lagrangian_strain(
98 1218 : declareProperty<RankTwoTensor>("total_lagrangian_strain")), // Lagrangian strain
99 1218 : _updated_rotation(declareProperty<RankTwoTensor>("updated_rotation")),
100 1218 : _crysrot(getMaterialProperty<RankTwoTensor>(
101 1218 : _base_name + "crysrot")), // defined in the elasticity tensor classes for crystal plasticity
102 4872 : _print_convergence_message(getParam<bool>("print_state_variable_convergence_error_messages"))
103 : {
104 1218 : _convergence_failed = false;
105 1218 : }
106 :
107 : void
108 27984 : ComputeMultipleCrystalPlasticityStress::initQpStatefulProperties()
109 : {
110 27984 : _plastic_deformation_gradient[_qp].zero();
111 27984 : _plastic_deformation_gradient[_qp].addIa(1.0);
112 :
113 27984 : if (_num_eigenstrains)
114 : {
115 7568 : (*_eigenstrain_deformation_gradient)[_qp].zero();
116 7568 : (*_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 20416 : _inverse_eigenstrain_deformation_grad.addIa(1.0);
123 : }
124 :
125 27984 : _pk2[_qp].zero();
126 :
127 27984 : _total_lagrangian_strain[_qp].zero();
128 :
129 27984 : _updated_rotation[_qp].zero();
130 27984 : _updated_rotation[_qp].addIa(1.0);
131 :
132 59168 : for (unsigned int i = 0; i < _num_models; ++i)
133 : {
134 31184 : _models[i]->setQp(_qp);
135 31184 : _models[i]->initQpStatefulProperties();
136 : }
137 :
138 35632 : for (unsigned int i = 0; i < _num_eigenstrains; ++i)
139 : {
140 7648 : _eigenstrains[i]->setQp(_qp);
141 7648 : _eigenstrains[i]->initQpStatefulProperties();
142 : }
143 27984 : }
144 :
145 : void
146 1134 : ComputeMultipleCrystalPlasticityStress::initialSetup()
147 : {
148 : // get crystal plasticity models
149 : std::vector<MaterialName> model_names =
150 2268 : getParam<std::vector<MaterialName>>("crystal_plasticity_models");
151 :
152 2412 : for (unsigned int i = 0; i < _num_models; ++i)
153 : {
154 : CrystalPlasticityStressUpdateBase * model =
155 1278 : dynamic_cast<CrystalPlasticityStressUpdateBase *>(&getMaterialByName(model_names[i]));
156 :
157 1278 : if (model)
158 : {
159 1278 : _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 2268 : getParam<std::vector<MaterialName>>("eigenstrain_names");
170 :
171 1350 : for (unsigned int i = 0; i < _num_eigenstrains; ++i)
172 : {
173 : ComputeCrystalPlasticityEigenstrainBase * eigenstrain =
174 216 : dynamic_cast<ComputeCrystalPlasticityEigenstrainBase *>(
175 216 : &getMaterialByName(eigenstrain_names[i]));
176 :
177 216 : if (eigenstrain)
178 216 : _eigenstrains.push_back(eigenstrain);
179 : else
180 0 : mooseError("Eigenstrain" + eigenstrain_names[i] +
181 : " is not compatible with ComputeMultipleCrystalPlasticityStress");
182 : }
183 1134 : }
184 :
185 : void
186 2255745 : ComputeMultipleCrystalPlasticityStress::computeQpStress()
187 : {
188 4797442 : for (unsigned int i = 0; i < _num_models; ++i)
189 : {
190 2541697 : _models[i]->setQp(_qp);
191 2541697 : _models[i]->setMaterialVectorSize();
192 : }
193 :
194 3061781 : for (unsigned int i = 0; i < _num_eigenstrains; ++i)
195 806036 : _eigenstrains[i]->setQp(_qp);
196 :
197 2255745 : updateStress(_stress[_qp], _Jacobian_mult[_qp]); // This is NOT the exact jacobian
198 2255425 : }
199 :
200 : void
201 2255745 : ComputeMultipleCrystalPlasticityStress::updateStress(RankTwoTensor & cauchy_stress,
202 : RankFourTensor & jacobian_mult)
203 : {
204 : // Does not support face/boundary material property calculation
205 2255745 : if (isBoundaryMaterial())
206 : return;
207 :
208 : // Initialize substepping variables
209 : unsigned int substep_iter = 1;
210 : unsigned int num_substep = 1;
211 :
212 2252545 : _temporary_deformation_gradient_old = _deformation_gradient_old[_qp];
213 2252545 : if (_temporary_deformation_gradient_old.det() == 0)
214 0 : _temporary_deformation_gradient_old.addIa(1.0);
215 :
216 2252545 : _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 4791042 : for (unsigned int i = 0; i < _num_models; ++i)
223 2538497 : _models[i]->calculateFlowDirection(_crysrot[_qp]);
224 :
225 : do
226 : {
227 2463776 : _convergence_failed = false;
228 2463776 : preSolveQp();
229 :
230 2463776 : _substep_dt = _dt / num_substep;
231 5213504 : for (unsigned int i = 0; i < _num_models; ++i)
232 2749728 : _models[i]->setSubstepDt(_substep_dt);
233 :
234 : // calculate F^{eigen} only when we have eigenstrain
235 : _inverse_eigenstrain_deformation_grad.zero();
236 2463776 : _inverse_eigenstrain_deformation_grad.addIa(1.0);
237 2463776 : if (_num_eigenstrains)
238 899824 : calculateEigenstrainDeformationGrad();
239 :
240 5235133 : for (unsigned int istep = 0; istep < num_substep; ++istep)
241 : {
242 2982908 : _temporary_deformation_gradient =
243 2982908 : (static_cast<Real>(istep) + 1) / num_substep * _delta_deformation_gradient;
244 2982908 : _temporary_deformation_gradient += _temporary_deformation_gradient_old;
245 :
246 2982908 : solveQp();
247 :
248 2982900 : if (_convergence_failed)
249 : {
250 211523 : 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 211523 : substep_iter++;
256 211523 : num_substep *= 2;
257 211523 : break;
258 : }
259 : }
260 :
261 2463748 : if (substep_iter > _max_substep_iter && _convergence_failed)
262 292 : mooseException("ComputeMultipleCrystalPlasticityStress: Constitutive failure");
263 2463456 : } while (_convergence_failed);
264 :
265 2252225 : postSolveQp(cauchy_stress, jacobian_mult);
266 : }
267 :
268 : void
269 2463776 : ComputeMultipleCrystalPlasticityStress::preSolveQp()
270 : {
271 5213504 : for (unsigned int i = 0; i < _num_models; ++i)
272 2749728 : _models[i]->setInitialConstitutiveVariableValues();
273 :
274 2463776 : _pk2[_qp] = _pk2_old[_qp];
275 2463776 : _inverse_plastic_deformation_grad_old = _plastic_deformation_gradient_old[_qp].inverse();
276 2463776 : }
277 :
278 : void
279 2982908 : ComputeMultipleCrystalPlasticityStress::solveQp()
280 : {
281 6251768 : for (unsigned int i = 0; i < _num_models; ++i)
282 : {
283 3268860 : _models[i]->setSubstepConstitutiveVariableValues();
284 3268860 : _models[i]->calculateSlipResistance();
285 : }
286 :
287 2982908 : _inverse_plastic_deformation_grad = _inverse_plastic_deformation_grad_old;
288 :
289 2982908 : solveStateVariables();
290 2982900 : if (_convergence_failed)
291 : return; // pop back up and take a smaller substep
292 :
293 5828706 : for (unsigned int i = 0; i < _num_models; ++i)
294 3057329 : _models[i]->updateSubstepConstitutiveVariableValues();
295 :
296 : // save off the old F^{p} inverse now that have converged on the stress and state variables
297 2771377 : _inverse_plastic_deformation_grad_old = _inverse_plastic_deformation_grad;
298 : }
299 :
300 : void
301 2252225 : ComputeMultipleCrystalPlasticityStress::postSolveQp(RankTwoTensor & cauchy_stress,
302 : RankFourTensor & jacobian_mult)
303 : {
304 2252225 : cauchy_stress = _elastic_deformation_gradient * _pk2[_qp] *
305 2252225 : _elastic_deformation_gradient.transpose() / _elastic_deformation_gradient.det();
306 :
307 2252225 : calcTangentModuli(jacobian_mult);
308 :
309 2252225 : _total_lagrangian_strain[_qp] =
310 2252225 : _deformation_gradient[_qp].transpose() * _deformation_gradient[_qp] -
311 2252225 : RankTwoTensor::Identity();
312 2252225 : _total_lagrangian_strain[_qp] = _total_lagrangian_strain[_qp] * 0.5;
313 :
314 : // Calculate crystal rotation to track separately
315 2252225 : RankTwoTensor rot;
316 2252225 : _elastic_deformation_gradient.getRUDecompositionRotation(rot);
317 2252225 : _updated_rotation[_qp] = rot * _crysrot[_qp];
318 2252225 : }
319 :
320 : void
321 2982908 : 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 4463042 : solveStress();
331 4463036 : if (_convergence_failed)
332 : return;
333 :
334 4254619 : _plastic_deformation_gradient[_qp] =
335 4254619 : _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 9073759 : for (unsigned int i = 0; i < _num_models; ++i)
342 4819140 : _models[i]->cacheStateVariablesBeforeUpdate();
343 :
344 9073759 : for (unsigned int i = 0; i < _num_models; ++i)
345 4819140 : _models[i]->calculateStateVariableEvolutionRateComponent();
346 :
347 9073757 : for (unsigned int i = 0; i < _num_models; ++i)
348 4819140 : if (!_models[i]->updateStateVariables())
349 146 : _convergence_failed = true;
350 :
351 9073755 : for (unsigned int i = 0; i < _num_models; ++i)
352 4819138 : _models[i]->calculateSlipResistance();
353 :
354 4254617 : if (_convergence_failed)
355 : return;
356 :
357 7314543 : 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 4543111 : 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 4254471 : if (iter_flag)
371 : {
372 1483039 : 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 2771432 : iteration++;
381 1483039 : } while (iter_flag && iteration < _maxiterg);
382 :
383 2774337 : if (iteration == _maxiterg)
384 : {
385 2960 : 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 2960 : _convergence_failed = true;
395 : }
396 : }
397 :
398 : void
399 4463042 : ComputeMultipleCrystalPlasticityStress::solveStress()
400 : {
401 : unsigned int iteration = 0;
402 4463042 : RankTwoTensor dpk2;
403 : Real rnorm, rnorm0, rnorm_prev;
404 :
405 : // Calculate stress residual
406 4463042 : calculateResidualAndJacobian();
407 4463042 : 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 4463042 : rnorm = _residual_tensor.L2norm();
420 4463042 : 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 22005899 : while (rnorm > _rtol * rnorm0 && rnorm > _abs_tol && iteration < _maxiter)
425 : {
426 : // Calculate stress increment
427 17744616 : dpk2 = -_jacobian.invSymm() * _residual_tensor;
428 17744616 : _pk2[_qp] = _pk2[_qp] + dpk2;
429 :
430 17744616 : calculateResidualAndJacobian();
431 :
432 17744612 : if (_convergence_failed)
433 : {
434 201755 : 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 201755 : return;
442 : }
443 :
444 17542857 : rnorm_prev = rnorm;
445 17542857 : rnorm = _residual_tensor.L2norm();
446 :
447 17542857 : 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 17542857 : if (_use_line_search)
457 945728 : rnorm = _residual_tensor.L2norm();
458 :
459 17542857 : iteration++;
460 : }
461 :
462 4261283 : if (iteration >= _maxiter)
463 : {
464 6664 : 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 6662 : _convergence_failed = true;
477 : }
478 : }
479 :
480 : // Calculates stress residual equation and jacobian
481 : void
482 22207658 : ComputeMultipleCrystalPlasticityStress::calculateResidualAndJacobian()
483 : {
484 22207658 : calculateResidual();
485 22207654 : if (_convergence_failed)
486 : return;
487 22005899 : calculateJacobian();
488 : }
489 :
490 : void
491 22272458 : ComputeMultipleCrystalPlasticityStress::calculateResidual()
492 : {
493 22272458 : RankTwoTensor ce, elastic_strain, ce_pk2, equivalent_slip_increment_per_model,
494 22272458 : equivalent_slip_increment, pk2_new;
495 :
496 : equivalent_slip_increment.zero();
497 :
498 : // calculate slip rate in order to compute F^{p-1}
499 47318982 : 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 25248283 : _models[i]->calculateShearStress(
505 25248283 : _pk2[_qp], _inverse_eigenstrain_deformation_grad, _num_eigenstrains);
506 :
507 25248283 : _convergence_failed = !_models[i]->calculateSlipRate();
508 :
509 25248279 : if (_convergence_failed)
510 : return;
511 :
512 25046524 : _models[i]->calculateEquivalentSlipIncrement(equivalent_slip_increment_per_model);
513 25046524 : equivalent_slip_increment += equivalent_slip_increment_per_model;
514 : }
515 :
516 : RankTwoTensor residual_equivalent_slip_increment =
517 22070699 : RankTwoTensor::Identity() - equivalent_slip_increment;
518 22070699 : _inverse_plastic_deformation_grad =
519 22070699 : _inverse_plastic_deformation_grad_old * residual_equivalent_slip_increment;
520 :
521 22070699 : _elastic_deformation_gradient = _temporary_deformation_gradient *
522 22070699 : _inverse_eigenstrain_deformation_grad *
523 : _inverse_plastic_deformation_grad;
524 :
525 22070699 : ce = _elastic_deformation_gradient.transpose() * _elastic_deformation_gradient;
526 :
527 22070699 : elastic_strain = ce - RankTwoTensor::Identity();
528 22070699 : elastic_strain *= 0.5;
529 :
530 22070699 : pk2_new = _elasticity_tensor[_qp] * elastic_strain;
531 22070699 : _residual_tensor = _pk2[_qp] - pk2_new;
532 : }
533 :
534 : void
535 22005899 : ComputeMultipleCrystalPlasticityStress::calculateJacobian()
536 : {
537 : // may not need to cache the dfpinvdpk2 here. need to double check
538 22005899 : RankFourTensor dfedfpinv, deedfe, dfpinvdpk2, dfpinvdpk2_per_model;
539 :
540 22005899 : RankTwoTensor ffeiginv = _temporary_deformation_gradient * _inverse_eigenstrain_deformation_grad;
541 :
542 88023596 : for (const auto i : make_range(Moose::dim))
543 264070788 : for (const auto j : make_range(Moose::dim))
544 792212364 : for (const auto k : make_range(Moose::dim))
545 594159273 : dfedfpinv(i, j, k, j) = ffeiginv(i, k);
546 :
547 88023596 : for (const auto i : make_range(Moose::dim))
548 264070788 : for (const auto j : make_range(Moose::dim))
549 792212364 : for (const auto k : make_range(Moose::dim))
550 : {
551 594159273 : deedfe(i, j, k, i) = deedfe(i, j, k, i) + _elastic_deformation_gradient(k, j) * 0.5;
552 594159273 : deedfe(i, j, k, j) = deedfe(i, j, k, j) + _elastic_deformation_gradient(k, i) * 0.5;
553 : }
554 :
555 46987623 : for (unsigned int i = 0; i < _num_models; ++i)
556 : {
557 24981724 : _models[i]->calculateTotalPlasticDeformationGradientDerivative(
558 : dfpinvdpk2_per_model,
559 24981724 : _inverse_plastic_deformation_grad_old,
560 24981724 : _inverse_eigenstrain_deformation_grad,
561 24981724 : _num_eigenstrains);
562 24981724 : dfpinvdpk2 += dfpinvdpk2_per_model;
563 : }
564 :
565 : _jacobian =
566 44011798 : RankFourTensor::IdentityFour() - (_elasticity_tensor[_qp] * deedfe * dfedfpinv * dfpinvdpk2);
567 22005899 : }
568 :
569 : void
570 2252225 : ComputeMultipleCrystalPlasticityStress::calcTangentModuli(RankFourTensor & jacobian_mult)
571 : {
572 2252225 : switch (_tan_mod_type)
573 : {
574 2252225 : case TangentModuliType::EXACT:
575 2252225 : elastoPlasticTangentModuli(jacobian_mult);
576 2252225 : break;
577 0 : default:
578 0 : elasticTangentModuli(jacobian_mult);
579 : }
580 2252225 : }
581 :
582 : void
583 2252225 : ComputeMultipleCrystalPlasticityStress::elastoPlasticTangentModuli(RankFourTensor & jacobian_mult)
584 : {
585 2252225 : RankFourTensor tan_mod;
586 2252225 : RankTwoTensor pk2fet, fepk2, feiginvfpinv;
587 2252225 : RankFourTensor deedfe, dsigdpk2dfe, dfedf;
588 :
589 : // Fill in the matrix stiffness material property
590 9008900 : for (const auto i : make_range(Moose::dim))
591 27026700 : for (const auto j : make_range(Moose::dim))
592 81080100 : for (const auto k : make_range(Moose::dim))
593 : {
594 60810075 : deedfe(i, j, k, i) = deedfe(i, j, k, i) + _elastic_deformation_gradient(k, j) * 0.5;
595 60810075 : 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 2252225 : dsigdpk2dfe = _elastic_deformation_gradient.times<i_, k_, j_, l_>(_elastic_deformation_gradient) *
600 2252225 : _elasticity_tensor[_qp] * deedfe;
601 :
602 2252225 : pk2fet = _pk2[_qp] * _elastic_deformation_gradient.transpose();
603 2252225 : fepk2 = _elastic_deformation_gradient * _pk2[_qp];
604 :
605 9008900 : for (const auto i : make_range(Moose::dim))
606 27026700 : for (const auto j : make_range(Moose::dim))
607 81080100 : for (const auto l : make_range(Moose::dim))
608 : {
609 60810075 : tan_mod(i, j, i, l) += pk2fet(l, j);
610 60810075 : tan_mod(i, j, j, l) += fepk2(i, l);
611 : }
612 :
613 2252225 : tan_mod += dsigdpk2dfe;
614 :
615 2252225 : const auto je = _elastic_deformation_gradient.det();
616 2252225 : if (je > 0.0)
617 2252225 : tan_mod /= je;
618 :
619 2252225 : feiginvfpinv = _inverse_eigenstrain_deformation_grad * _inverse_plastic_deformation_grad;
620 9008900 : for (const auto i : make_range(Moose::dim))
621 27026700 : for (const auto j : make_range(Moose::dim))
622 81080100 : for (const auto l : make_range(Moose::dim))
623 60810075 : dfedf(i, j, i, l) = feiginvfpinv(l, j);
624 :
625 2252225 : jacobian_mult = tan_mod * dfedf;
626 2252225 : }
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 64800 : ComputeMultipleCrystalPlasticityStress::lineSearchUpdate(const Real & rnorm_prev,
637 : const RankTwoTensor & dpk2)
638 : {
639 64800 : if (_line_search_method == LineSearchMethod::CutHalf)
640 : {
641 : Real rnorm;
642 64800 : Real step = 1.0;
643 :
644 : do
645 : {
646 64800 : _pk2[_qp] = _pk2[_qp] - step * dpk2;
647 64800 : step /= 2.0;
648 64800 : _pk2[_qp] = _pk2[_qp] + step * dpk2;
649 :
650 64800 : calculateResidual();
651 64800 : rnorm = _residual_tensor.L2norm();
652 64800 : } 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 64800 : 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 899824 : ComputeMultipleCrystalPlasticityStress::calculateEigenstrainDeformationGrad()
712 : {
713 : _inverse_eigenstrain_deformation_grad.zero();
714 899824 : _inverse_eigenstrain_deformation_grad.addIa(1.0);
715 :
716 1911612 : for (unsigned int i = 0; i < _num_eigenstrains; ++i)
717 : {
718 1011808 : _eigenstrains[i]->setSubstepDt(_substep_dt);
719 1011808 : _eigenstrains[i]->computeQpProperties();
720 1011788 : _inverse_eigenstrain_deformation_grad *= _eigenstrains[i]->getDeformationGradientInverse();
721 : }
722 899804 : (*_eigenstrain_deformation_gradient)[_qp] = _inverse_eigenstrain_deformation_grad.inverse();
723 899804 : }
|
