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