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 "FiniteStrainUObasedCP.h"
11 : #include "petscblaslapack.h"
12 : #include "MooseException.h"
13 : #include "CrystalPlasticitySlipRate.h"
14 : #include "CrystalPlasticitySlipResistance.h"
15 : #include "CrystalPlasticityStateVariable.h"
16 : #include "CrystalPlasticityStateVarRateComponent.h"
17 :
18 : registerMooseObject("TensorMechanicsApp", FiniteStrainUObasedCP);
19 :
20 : InputParameters
21 156 : FiniteStrainUObasedCP::validParams()
22 : {
23 156 : InputParameters params = ComputeStressBase::validParams();
24 156 : params.addClassDescription("UserObject based Crystal Plasticity system.");
25 312 : params.addParam<Real>("rtol", 1e-6, "Constitutive stress residue relative tolerance");
26 312 : params.addParam<Real>("abs_tol", 1e-6, "Constitutive stress residue absolute tolerance");
27 312 : params.addParam<Real>(
28 312 : "stol", 1e-2, "Constitutive slip system resistance relative residual tolerance");
29 312 : params.addParam<Real>(
30 312 : "zero_tol", 1e-12, "Tolerance for residual check when variable value is zero");
31 312 : params.addParam<unsigned int>("maxiter", 100, "Maximum number of iterations for stress update");
32 312 : params.addParam<unsigned int>(
33 312 : "maxiter_state_variable", 100, "Maximum number of iterations for state variable update");
34 312 : MooseEnum tan_mod_options("exact none", "none"); // Type of read
35 312 : params.addParam<MooseEnum>("tan_mod_type",
36 : tan_mod_options,
37 : "Type of tangent moduli for preconditioner: default elastic");
38 312 : params.addParam<unsigned int>(
39 312 : "maximum_substep_iteration", 1, "Maximum number of substep iteration");
40 312 : params.addParam<bool>("use_line_search", false, "Use line search in constitutive update");
41 312 : params.addParam<Real>("min_line_search_step_size", 0.01, "Minimum line search step size");
42 312 : params.addParam<Real>("line_search_tol", 0.5, "Line search bisection method tolerance");
43 312 : params.addParam<unsigned int>(
44 312 : "line_search_maxiter", 20, "Line search bisection method maximum number of iteration");
45 312 : MooseEnum line_search_method("CUT_HALF BISECTION", "CUT_HALF");
46 312 : params.addParam<MooseEnum>(
47 : "line_search_method", line_search_method, "The method used in line search");
48 312 : params.addRequiredParam<std::vector<UserObjectName>>(
49 : "uo_slip_rates",
50 : "List of names of user objects that define the slip rates for this material.");
51 312 : params.addRequiredParam<std::vector<UserObjectName>>(
52 : "uo_slip_resistances",
53 : "List of names of user objects that define the slip resistances for this material.");
54 312 : params.addRequiredParam<std::vector<UserObjectName>>(
55 : "uo_state_vars",
56 : "List of names of user objects that define the state variable for this material.");
57 312 : params.addRequiredParam<std::vector<UserObjectName>>(
58 : "uo_state_var_evol_rate_comps",
59 : "List of names of user objects that define the state "
60 : "variable evolution rate components for this material.");
61 156 : return params;
62 156 : }
63 :
64 117 : FiniteStrainUObasedCP::FiniteStrainUObasedCP(const InputParameters & parameters)
65 : : ComputeStressBase(parameters),
66 117 : _num_uo_slip_rates(parameters.get<std::vector<UserObjectName>>("uo_slip_rates").size()),
67 117 : _num_uo_slip_resistances(
68 117 : parameters.get<std::vector<UserObjectName>>("uo_slip_resistances").size()),
69 117 : _num_uo_state_vars(parameters.get<std::vector<UserObjectName>>("uo_state_vars").size()),
70 117 : _num_uo_state_var_evol_rate_comps(
71 117 : parameters.get<std::vector<UserObjectName>>("uo_state_var_evol_rate_comps").size()),
72 234 : _rtol(getParam<Real>("rtol")),
73 234 : _abs_tol(getParam<Real>("abs_tol")),
74 234 : _stol(getParam<Real>("stol")),
75 234 : _zero_tol(getParam<Real>("zero_tol")),
76 234 : _maxiter(getParam<unsigned int>("maxiter")),
77 234 : _maxiterg(getParam<unsigned int>("maxiter_state_variable")),
78 234 : _tan_mod_type(getParam<MooseEnum>("tan_mod_type")),
79 234 : _max_substep_iter(getParam<unsigned int>("maximum_substep_iteration")),
80 234 : _use_line_search(getParam<bool>("use_line_search")),
81 234 : _min_lsrch_step(getParam<Real>("min_line_search_step_size")),
82 234 : _lsrch_tol(getParam<Real>("line_search_tol")),
83 234 : _lsrch_max_iter(getParam<unsigned int>("line_search_maxiter")),
84 234 : _lsrch_method(getParam<MooseEnum>("line_search_method")),
85 117 : _fp(declareProperty<RankTwoTensor>("fp")), // Plastic deformation gradient
86 234 : _fp_old(getMaterialPropertyOld<RankTwoTensor>(
87 : "fp")), // Plastic deformation gradient of previous increment
88 117 : _pk2(declareProperty<RankTwoTensor>("pk2")), // 2nd Piola-Kirchoff Stress
89 234 : _pk2_old(getMaterialPropertyOld<RankTwoTensor>(
90 : "pk2")), // 2nd Piola Kirchoff Stress of previous increment
91 117 : _lag_e(declareProperty<RankTwoTensor>("lage")), // Lagrangian strain
92 117 : _update_rot(declareProperty<RankTwoTensor>(
93 : "update_rot")), // Rotation tensor considering material rotation and crystal orientation
94 234 : _update_rot_old(getMaterialPropertyOld<RankTwoTensor>("update_rot")),
95 117 : _elasticity_tensor_name(_base_name + "elasticity_tensor"),
96 117 : _elasticity_tensor(getMaterialPropertyByName<RankFourTensor>(_elasticity_tensor_name)),
97 234 : _deformation_gradient(getMaterialProperty<RankTwoTensor>("deformation_gradient")),
98 234 : _deformation_gradient_old(getMaterialPropertyOld<RankTwoTensor>("deformation_gradient")),
99 585 : _crysrot(getMaterialProperty<RankTwoTensor>("crysrot"))
100 : {
101 117 : _err_tol = false;
102 :
103 : _delta_dfgrd.zero();
104 :
105 : // resize the material properties for each userobject
106 117 : _mat_prop_slip_rates.resize(_num_uo_slip_rates);
107 117 : _mat_prop_slip_resistances.resize(_num_uo_slip_resistances);
108 117 : _mat_prop_state_vars.resize(_num_uo_state_vars);
109 117 : _mat_prop_state_vars_old.resize(_num_uo_state_vars);
110 117 : _mat_prop_state_var_evol_rate_comps.resize(_num_uo_state_var_evol_rate_comps);
111 :
112 : // resize the flow direction
113 117 : _flow_direction.resize(_num_uo_slip_rates);
114 :
115 : // resize local state variables
116 117 : _state_vars_old.resize(_num_uo_state_vars);
117 117 : _state_vars_old_stored.resize(_num_uo_state_vars);
118 117 : _state_vars_prev.resize(_num_uo_state_vars);
119 :
120 : // resize user objects
121 117 : _uo_slip_rates.resize(_num_uo_slip_rates);
122 117 : _uo_slip_resistances.resize(_num_uo_slip_resistances);
123 117 : _uo_state_vars.resize(_num_uo_state_vars);
124 117 : _uo_state_var_evol_rate_comps.resize(_num_uo_state_var_evol_rate_comps);
125 :
126 : // assign the user objects
127 234 : for (unsigned int i = 0; i < _num_uo_slip_rates; ++i)
128 : {
129 117 : _uo_slip_rates[i] = &getUserObjectByName<CrystalPlasticitySlipRate>(
130 117 : parameters.get<std::vector<UserObjectName>>("uo_slip_rates")[i]);
131 117 : _mat_prop_slip_rates[i] = &declareProperty<std::vector<Real>>(
132 234 : parameters.get<std::vector<UserObjectName>>("uo_slip_rates")[i]);
133 117 : _flow_direction[i] = &declareProperty<std::vector<RankTwoTensor>>(
134 234 : parameters.get<std::vector<UserObjectName>>("uo_slip_rates")[i] + "_flow_direction");
135 : }
136 :
137 234 : for (unsigned int i = 0; i < _num_uo_slip_resistances; ++i)
138 : {
139 117 : _uo_slip_resistances[i] = &getUserObjectByName<CrystalPlasticitySlipResistance>(
140 117 : parameters.get<std::vector<UserObjectName>>("uo_slip_resistances")[i]);
141 117 : _mat_prop_slip_resistances[i] = &declareProperty<std::vector<Real>>(
142 234 : parameters.get<std::vector<UserObjectName>>("uo_slip_resistances")[i]);
143 : }
144 :
145 234 : for (unsigned int i = 0; i < _num_uo_state_vars; ++i)
146 : {
147 117 : _uo_state_vars[i] = &getUserObjectByName<CrystalPlasticityStateVariable>(
148 117 : parameters.get<std::vector<UserObjectName>>("uo_state_vars")[i]);
149 117 : _mat_prop_state_vars[i] = &declareProperty<std::vector<Real>>(
150 234 : parameters.get<std::vector<UserObjectName>>("uo_state_vars")[i]);
151 117 : _mat_prop_state_vars_old[i] = &getMaterialPropertyOld<std::vector<Real>>(
152 234 : parameters.get<std::vector<UserObjectName>>("uo_state_vars")[i]);
153 : }
154 :
155 234 : for (unsigned int i = 0; i < _num_uo_state_var_evol_rate_comps; ++i)
156 : {
157 117 : _uo_state_var_evol_rate_comps[i] = &getUserObjectByName<CrystalPlasticityStateVarRateComponent>(
158 117 : parameters.get<std::vector<UserObjectName>>("uo_state_var_evol_rate_comps")[i]);
159 117 : _mat_prop_state_var_evol_rate_comps[i] = &declareProperty<std::vector<Real>>(
160 234 : parameters.get<std::vector<UserObjectName>>("uo_state_var_evol_rate_comps")[i]);
161 : }
162 :
163 117 : _substep_dt = 0.0;
164 117 : }
165 :
166 : void
167 3104 : FiniteStrainUObasedCP::initQpStatefulProperties()
168 : {
169 6208 : for (unsigned int i = 0; i < _num_uo_slip_rates; ++i)
170 : {
171 3104 : (*_mat_prop_slip_rates[i])[_qp].resize(_uo_slip_rates[i]->variableSize());
172 3104 : (*_flow_direction[i])[_qp].resize(_uo_slip_rates[i]->variableSize());
173 : }
174 :
175 6208 : for (unsigned int i = 0; i < _num_uo_slip_resistances; ++i)
176 3104 : (*_mat_prop_slip_resistances[i])[_qp].resize(_uo_slip_resistances[i]->variableSize());
177 :
178 6208 : for (unsigned int i = 0; i < _num_uo_state_vars; ++i)
179 : {
180 3104 : (*_mat_prop_state_vars[i])[_qp].resize(_uo_state_vars[i]->variableSize());
181 3104 : _state_vars_old[i].resize(_uo_state_vars[i]->variableSize());
182 3104 : _state_vars_old_stored[i].resize(_uo_state_vars[i]->variableSize());
183 3104 : _state_vars_prev[i].resize(_uo_state_vars[i]->variableSize());
184 : }
185 :
186 6208 : for (unsigned int i = 0; i < _num_uo_state_var_evol_rate_comps; ++i)
187 3104 : (*_mat_prop_state_var_evol_rate_comps[i])[_qp].resize(
188 3104 : _uo_state_var_evol_rate_comps[i]->variableSize());
189 :
190 3104 : _stress[_qp].zero();
191 3104 : _pk2[_qp].zero();
192 3104 : _lag_e[_qp].zero();
193 :
194 3104 : _fp[_qp].setToIdentity();
195 3104 : _update_rot[_qp].setToIdentity();
196 :
197 6208 : for (unsigned int i = 0; i < _num_uo_state_vars; ++i)
198 : // Initializes slip system related properties
199 3104 : _uo_state_vars[i]->initSlipSysProps((*_mat_prop_state_vars[i])[_qp], _q_point[_qp]);
200 3104 : }
201 :
202 : /**
203 : * Solves stress residual equation using NR.
204 : * Updates slip system resistances iteratively.
205 : */
206 : void
207 172295 : FiniteStrainUObasedCP::computeQpStress()
208 : {
209 : // Userobject based crystal plasticity does not support face/boundary material property
210 : // calculation.
211 172295 : if (isBoundaryMaterial())
212 : return;
213 : // Depth of substepping; Limited to maximum substep iteration
214 : unsigned int substep_iter = 1;
215 : // Calculated from substep_iter as 2^substep_iter
216 : unsigned int num_substep = 1;
217 :
218 171015 : _dfgrd_tmp_old = _deformation_gradient_old[_qp];
219 171015 : if (_dfgrd_tmp_old.det() == 0)
220 0 : _dfgrd_tmp_old.addIa(1.0);
221 :
222 171015 : _delta_dfgrd = _deformation_gradient[_qp] - _dfgrd_tmp_old;
223 :
224 : // Saves the old stateful properties that are modified during sub stepping
225 342030 : for (unsigned int i = 0; i < _num_uo_state_vars; ++i)
226 171015 : _state_vars_old[i] = (*_mat_prop_state_vars_old[i])[_qp];
227 :
228 342030 : for (unsigned int i = 0; i < _num_uo_slip_rates; ++i)
229 171015 : _uo_slip_rates[i]->calcFlowDirection(_qp, (*_flow_direction[i])[_qp]);
230 :
231 : do
232 : {
233 171735 : _err_tol = false;
234 :
235 171735 : preSolveQp();
236 :
237 171735 : _substep_dt = _dt / num_substep;
238 :
239 344407 : for (unsigned int istep = 0; istep < num_substep; ++istep)
240 : {
241 173415 : _dfgrd_tmp = (static_cast<Real>(istep) + 1) / num_substep * _delta_dfgrd + _dfgrd_tmp_old;
242 :
243 173415 : solveQp();
244 :
245 173415 : if (_err_tol)
246 : {
247 743 : substep_iter++;
248 743 : num_substep *= 2;
249 743 : break;
250 : }
251 : }
252 :
253 171735 : if (substep_iter > _max_substep_iter && _err_tol)
254 23 : throw MooseException("FiniteStrainUObasedCP: Constitutive failure.");
255 171712 : } while (_err_tol);
256 :
257 170992 : postSolveQp();
258 : }
259 :
260 : void
261 171735 : FiniteStrainUObasedCP::preSolveQp()
262 : {
263 343470 : for (unsigned int i = 0; i < _num_uo_state_vars; ++i)
264 171735 : (*_mat_prop_state_vars[i])[_qp] = _state_vars_old_stored[i] = _state_vars_old[i];
265 :
266 171735 : _pk2[_qp] = _pk2_old[_qp];
267 171735 : _fp_old_inv = _fp_old[_qp].inverse();
268 171735 : }
269 :
270 : void
271 173415 : FiniteStrainUObasedCP::solveQp()
272 : {
273 173415 : preSolveStatevar();
274 173415 : solveStatevar();
275 173415 : if (_err_tol)
276 : return;
277 :
278 172672 : postSolveStatevar();
279 : }
280 :
281 : void
282 170992 : FiniteStrainUObasedCP::postSolveQp()
283 : {
284 170992 : _stress[_qp] = _fe * _pk2[_qp] * _fe.transpose() / _fe.det();
285 :
286 : // Calculate jacobian for preconditioner
287 170992 : calcTangentModuli();
288 :
289 170992 : RankTwoTensor iden(RankTwoTensor::initIdentity);
290 :
291 170992 : _lag_e[_qp] = _deformation_gradient[_qp].transpose() * _deformation_gradient[_qp] - iden;
292 170992 : _lag_e[_qp] = _lag_e[_qp] * 0.5;
293 :
294 170992 : RankTwoTensor rot;
295 : // Calculate material rotation
296 170992 : _deformation_gradient[_qp].getRUDecompositionRotation(rot);
297 170992 : _update_rot[_qp] = rot * _crysrot[_qp];
298 170992 : }
299 :
300 : void
301 173415 : FiniteStrainUObasedCP::preSolveStatevar()
302 : {
303 346830 : for (unsigned int i = 0; i < _num_uo_state_vars; ++i)
304 173415 : (*_mat_prop_state_vars[i])[_qp] = _state_vars_old_stored[i];
305 :
306 346830 : for (unsigned int i = 0; i < _num_uo_slip_resistances; ++i)
307 173415 : _uo_slip_resistances[i]->calcSlipResistance(_qp, (*_mat_prop_slip_resistances[i])[_qp]);
308 :
309 173415 : _fp_inv = _fp_old_inv;
310 173415 : }
311 :
312 : void
313 173415 : FiniteStrainUObasedCP::solveStatevar()
314 : {
315 : unsigned int iterg;
316 : bool iter_flag = true;
317 :
318 : iterg = 0;
319 : // Check for slip system resistance update tolerance
320 349511 : while (iter_flag && iterg < _maxiterg)
321 : {
322 176839 : preSolveStress();
323 176839 : solveStress();
324 176839 : if (_err_tol)
325 : return;
326 :
327 176096 : postSolveStress();
328 :
329 : // Update slip system resistance and state variable
330 176096 : updateSlipSystemResistanceAndStateVariable();
331 :
332 176096 : if (_err_tol)
333 : return;
334 :
335 176096 : iter_flag = isStateVariablesConverged();
336 176096 : iterg++;
337 : }
338 :
339 172672 : if (iterg == _maxiterg)
340 : {
341 : #ifdef DEBUG
342 : mooseWarning("FiniteStrainUObasedCP: Hardness Integration error\n");
343 : #endif
344 0 : _err_tol = true;
345 : }
346 : }
347 :
348 : bool
349 176096 : FiniteStrainUObasedCP::isStateVariablesConverged()
350 : {
351 : Real diff;
352 :
353 348768 : for (unsigned int i = 0; i < _num_uo_state_vars; ++i)
354 : {
355 176096 : unsigned int n = (*_mat_prop_state_vars[i])[_qp].size();
356 2543072 : for (unsigned j = 0; j < n; j++)
357 : {
358 2370400 : diff = std::abs((*_mat_prop_state_vars[i])[_qp][j] - _state_vars_prev[i][j]);
359 :
360 2370400 : if (std::abs(_state_vars_old_stored[i][j]) < _zero_tol && diff > _zero_tol)
361 : return true;
362 2370400 : if (std::abs(_state_vars_old_stored[i][j]) > _zero_tol &&
363 2370400 : diff > _stol * std::abs(_state_vars_old_stored[i][j]))
364 : return true;
365 : }
366 : }
367 : return false;
368 : }
369 :
370 : void
371 172672 : FiniteStrainUObasedCP::postSolveStatevar()
372 : {
373 345344 : for (unsigned int i = 0; i < _num_uo_state_vars; ++i)
374 172672 : _state_vars_old_stored[i] = (*_mat_prop_state_vars[i])[_qp];
375 :
376 172672 : _fp_old_inv = _fp_inv;
377 172672 : }
378 :
379 : void
380 176839 : FiniteStrainUObasedCP::preSolveStress()
381 : {
382 176839 : }
383 :
384 : void
385 176839 : FiniteStrainUObasedCP::solveStress()
386 : {
387 : unsigned int iter = 0;
388 176839 : RankTwoTensor dpk2;
389 : Real rnorm, rnorm0, rnorm_prev;
390 :
391 : // Calculate stress residual
392 176839 : calcResidJacob();
393 176839 : if (_err_tol)
394 : {
395 : #ifdef DEBUG
396 : mooseWarning("FiniteStrainUObasedCP: Slip increment exceeds tolerance - Element number ",
397 : _current_elem->id(),
398 : " Gauss point = ",
399 : _qp);
400 : #endif
401 : return;
402 : }
403 :
404 176839 : rnorm = _resid.L2norm();
405 : rnorm0 = rnorm;
406 :
407 : // Check for stress residual tolerance
408 619157 : while (rnorm > _rtol * rnorm0 && rnorm0 > _abs_tol && iter < _maxiter)
409 : {
410 : // Calculate stress increment
411 443061 : dpk2 = -_jac.invSymm() * _resid;
412 443061 : _pk2[_qp] = _pk2[_qp] + dpk2;
413 :
414 443061 : calcResidJacob();
415 :
416 443061 : if (_err_tol)
417 : {
418 : #ifdef DEBUG
419 : mooseWarning("FiniteStrainUObasedCP: Slip increment exceeds tolerance - Element number ",
420 : _current_elem->id(),
421 : " Gauss point = ",
422 : _qp);
423 : #endif
424 : return;
425 : }
426 :
427 : rnorm_prev = rnorm;
428 442318 : rnorm = _resid.L2norm();
429 :
430 442318 : if (_use_line_search && rnorm > rnorm_prev && !lineSearchUpdate(rnorm_prev, dpk2))
431 : {
432 : #ifdef DEBUG
433 : mooseWarning("FiniteStrainUObasedCP: Failed with line search");
434 : #endif
435 0 : _err_tol = true;
436 0 : return;
437 : }
438 :
439 442318 : if (_use_line_search)
440 6592 : rnorm = _resid.L2norm();
441 :
442 442318 : iter++;
443 : }
444 :
445 176096 : if (iter >= _maxiter)
446 : {
447 : #ifdef DEBUG
448 : mooseWarning("FiniteStrainUObasedCP: Stress Integration error rmax = ",
449 : rnorm,
450 : " and the tolerance is ",
451 : _rtol * rnorm0,
452 : " when the rnorm0 value is ",
453 : rnorm0,
454 : " for element ",
455 : _current_elem->id(),
456 : " and qp ",
457 : _qp);
458 : #endif
459 0 : _err_tol = true;
460 : }
461 : }
462 :
463 : void
464 176096 : FiniteStrainUObasedCP::postSolveStress()
465 : {
466 176096 : _fp[_qp] = _fp_inv.inverse();
467 176096 : }
468 :
469 : void
470 176096 : FiniteStrainUObasedCP::updateSlipSystemResistanceAndStateVariable()
471 : {
472 352192 : for (unsigned int i = 0; i < _num_uo_state_vars; ++i)
473 176096 : _state_vars_prev[i] = (*_mat_prop_state_vars[i])[_qp];
474 :
475 352192 : for (unsigned int i = 0; i < _num_uo_state_var_evol_rate_comps; ++i)
476 176096 : _uo_state_var_evol_rate_comps[i]->calcStateVariableEvolutionRateComponent(
477 176096 : _qp, (*_mat_prop_state_var_evol_rate_comps[i])[_qp]);
478 :
479 352192 : for (unsigned int i = 0; i < _num_uo_state_vars; ++i)
480 : {
481 176096 : if (!_uo_state_vars[i]->updateStateVariable(
482 176096 : _qp, _substep_dt, (*_mat_prop_state_vars[i])[_qp], _state_vars_old_stored[i]))
483 0 : _err_tol = true;
484 : }
485 :
486 352192 : for (unsigned int i = 0; i < _num_uo_slip_resistances; ++i)
487 176096 : _uo_slip_resistances[i]->calcSlipResistance(_qp, (*_mat_prop_slip_resistances[i])[_qp]);
488 176096 : }
489 :
490 : // Calculates stress residual equation and jacobian
491 : void
492 619900 : FiniteStrainUObasedCP::calcResidJacob()
493 : {
494 619900 : calcResidual();
495 619900 : if (_err_tol)
496 : return;
497 619157 : calcJacobian();
498 : }
499 :
500 : void
501 619900 : FiniteStrainUObasedCP::getSlipRates()
502 : {
503 1239057 : for (unsigned int i = 0; i < _num_uo_slip_rates; ++i)
504 : {
505 619900 : if (!_uo_slip_rates[i]->calcSlipRate(_qp, _substep_dt, (*_mat_prop_slip_rates[i])[_qp]))
506 : {
507 743 : _err_tol = true;
508 743 : return;
509 : }
510 : }
511 : }
512 :
513 : void
514 619900 : FiniteStrainUObasedCP::calcResidual()
515 : {
516 619900 : RankTwoTensor iden(RankTwoTensor::initIdentity), ce, ee, ce_pk2, eqv_slip_incr, pk2_new;
517 :
518 619900 : getSlipRates();
519 619900 : if (_err_tol)
520 : return;
521 :
522 1238314 : for (unsigned int i = 0; i < _num_uo_slip_rates; ++i)
523 9002717 : for (unsigned int j = 0; j < _uo_slip_rates[i]->variableSize(); ++j)
524 : eqv_slip_incr +=
525 8383560 : (*_flow_direction[i])[_qp][j] * (*_mat_prop_slip_rates[i])[_qp][j] * _substep_dt;
526 :
527 619157 : eqv_slip_incr = iden - eqv_slip_incr;
528 619157 : _fp_inv = _fp_old_inv * eqv_slip_incr;
529 619157 : _fe = _dfgrd_tmp * _fp_inv;
530 :
531 619157 : ce = _fe.transpose() * _fe;
532 619157 : ee = ce - iden;
533 619157 : ee *= 0.5;
534 :
535 619157 : pk2_new = _elasticity_tensor[_qp] * ee;
536 :
537 619157 : _resid = _pk2[_qp] - pk2_new;
538 : }
539 :
540 : void
541 619157 : FiniteStrainUObasedCP::calcJacobian()
542 : {
543 619157 : RankFourTensor dfedfpinv, deedfe, dfpinvdpk2;
544 :
545 2476628 : for (const auto i : make_range(Moose::dim))
546 7429884 : for (const auto j : make_range(Moose::dim))
547 22289652 : for (const auto k : make_range(Moose::dim))
548 16717239 : dfedfpinv(i, j, k, j) = _dfgrd_tmp(i, k);
549 :
550 2476628 : for (const auto i : make_range(Moose::dim))
551 7429884 : for (const auto j : make_range(Moose::dim))
552 22289652 : for (const auto k : make_range(Moose::dim))
553 : {
554 16717239 : deedfe(i, j, k, i) = deedfe(i, j, k, i) + _fe(k, j) * 0.5;
555 16717239 : deedfe(i, j, k, j) = deedfe(i, j, k, j) + _fe(k, i) * 0.5;
556 : }
557 :
558 1238314 : for (unsigned int i = 0; i < _num_uo_slip_rates; ++i)
559 : {
560 619157 : unsigned int nss = _uo_slip_rates[i]->variableSize();
561 619157 : std::vector<RankTwoTensor> dtaudpk2(nss), dfpinvdslip(nss);
562 : std::vector<Real> dslipdtau;
563 619157 : dslipdtau.resize(nss);
564 619157 : _uo_slip_rates[i]->calcSlipRateDerivative(_qp, _substep_dt, dslipdtau);
565 9002717 : for (unsigned int j = 0; j < nss; j++)
566 : {
567 8383560 : dtaudpk2[j] = (*_flow_direction[i])[_qp][j];
568 8383560 : dfpinvdslip[j] = -_fp_old_inv * (*_flow_direction[i])[_qp][j];
569 8383560 : dfpinvdpk2 += (dfpinvdslip[j] * dslipdtau[j] * _substep_dt).outerProduct(dtaudpk2[j]);
570 : }
571 : }
572 : _jac =
573 1238314 : RankFourTensor::IdentityFour() - (_elasticity_tensor[_qp] * deedfe * dfedfpinv * dfpinvdpk2);
574 619157 : }
575 :
576 : void
577 170992 : FiniteStrainUObasedCP::calcTangentModuli()
578 : {
579 170992 : switch (_tan_mod_type)
580 : {
581 170992 : case 0:
582 170992 : elastoPlasticTangentModuli();
583 170992 : break;
584 0 : default:
585 0 : elasticTangentModuli();
586 : }
587 170992 : }
588 :
589 : void
590 170992 : FiniteStrainUObasedCP::elastoPlasticTangentModuli()
591 : {
592 170992 : RankFourTensor tan_mod;
593 170992 : RankTwoTensor pk2fet, fepk2;
594 170992 : RankFourTensor deedfe, dsigdpk2dfe, dfedf;
595 :
596 : // Fill in the matrix stiffness material property
597 683968 : for (const auto i : make_range(Moose::dim))
598 2051904 : for (const auto j : make_range(Moose::dim))
599 6155712 : for (const auto k : make_range(Moose::dim))
600 : {
601 4616784 : deedfe(i, j, k, i) = deedfe(i, j, k, i) + _fe(k, j) * 0.5;
602 4616784 : deedfe(i, j, k, j) = deedfe(i, j, k, j) + _fe(k, i) * 0.5;
603 : }
604 :
605 : usingTensorIndices(i_, j_, k_, l_);
606 170992 : dsigdpk2dfe = _fe.times<i_, k_, j_, l_>(_fe) * _elasticity_tensor[_qp] * deedfe;
607 :
608 170992 : pk2fet = _pk2[_qp] * _fe.transpose();
609 170992 : fepk2 = _fe * _pk2[_qp];
610 :
611 683968 : for (const auto i : make_range(Moose::dim))
612 2051904 : for (const auto j : make_range(Moose::dim))
613 6155712 : for (const auto l : make_range(Moose::dim))
614 : {
615 4616784 : tan_mod(i, j, i, l) += pk2fet(l, j);
616 4616784 : tan_mod(i, j, j, l) += fepk2(i, l);
617 : }
618 :
619 170992 : tan_mod += dsigdpk2dfe;
620 :
621 170992 : Real je = _fe.det();
622 170992 : if (je > 0.0)
623 170992 : tan_mod /= je;
624 :
625 683968 : for (const auto i : make_range(Moose::dim))
626 2051904 : for (const auto j : make_range(Moose::dim))
627 6155712 : for (const auto l : make_range(Moose::dim))
628 4616784 : dfedf(i, j, i, l) = _fp_inv(l, j);
629 :
630 170992 : _Jacobian_mult[_qp] = tan_mod * dfedf;
631 170992 : }
632 :
633 : void
634 0 : FiniteStrainUObasedCP::elasticTangentModuli()
635 : {
636 : // update jacobian_mult
637 0 : _Jacobian_mult[_qp] = _elasticity_tensor[_qp];
638 0 : }
639 :
640 : bool
641 0 : FiniteStrainUObasedCP::lineSearchUpdate(const Real rnorm_prev, const RankTwoTensor dpk2)
642 : {
643 0 : switch (_lsrch_method)
644 : {
645 0 : case 0: // CUT_HALF
646 : {
647 : Real rnorm;
648 0 : Real step = 1.0;
649 :
650 : do
651 : {
652 0 : _pk2[_qp] = _pk2[_qp] - step * dpk2;
653 0 : step /= 2.0;
654 0 : _pk2[_qp] = _pk2[_qp] + step * dpk2;
655 :
656 0 : calcResidual();
657 0 : rnorm = _resid.L2norm();
658 0 : } while (rnorm > rnorm_prev && step > _min_lsrch_step);
659 :
660 : // has norm improved or is the step still above minumum search step size?
661 0 : return (rnorm <= rnorm_prev || step > _min_lsrch_step);
662 : }
663 :
664 0 : case 1: // BISECTION
665 : {
666 : unsigned int count = 0;
667 : Real step_a = 0.0;
668 : Real step_b = 1.0;
669 0 : Real step = 1.0;
670 : Real s_m = 1000.0;
671 : Real rnorm = 1000.0;
672 :
673 0 : calcResidual();
674 0 : Real s_b = _resid.doubleContraction(dpk2);
675 0 : Real rnorm1 = _resid.L2norm();
676 0 : _pk2[_qp] = _pk2[_qp] - dpk2;
677 0 : calcResidual();
678 0 : Real s_a = _resid.doubleContraction(dpk2);
679 0 : Real rnorm0 = _resid.L2norm();
680 0 : _pk2[_qp] = _pk2[_qp] + dpk2;
681 :
682 0 : if ((rnorm1 / rnorm0) < _lsrch_tol || s_a * s_b > 0)
683 : {
684 0 : calcResidual();
685 0 : return true;
686 : }
687 :
688 0 : while ((rnorm / rnorm0) > _lsrch_tol && count < _lsrch_max_iter)
689 : {
690 0 : _pk2[_qp] = _pk2[_qp] - step * dpk2;
691 0 : step = 0.5 * (step_b + step_a);
692 0 : _pk2[_qp] = _pk2[_qp] + step * dpk2;
693 0 : calcResidual();
694 0 : s_m = _resid.doubleContraction(dpk2);
695 0 : rnorm = _resid.L2norm();
696 :
697 0 : if (s_m * s_a < 0.0)
698 : {
699 : step_b = step;
700 : s_b = s_m;
701 : }
702 0 : if (s_m * s_b < 0.0)
703 : {
704 : step_a = step;
705 : s_a = s_m;
706 : }
707 0 : count++;
708 : }
709 :
710 : // below tolerance and max iterations?
711 0 : return ((rnorm / rnorm0) < _lsrch_tol && count < _lsrch_max_iter);
712 : }
713 :
714 0 : default:
715 0 : mooseError("Line search method is not provided.");
716 : }
717 : }
|
