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 "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("SolidMechanicsApp", FiniteStrainUObasedCP);
19 :
20 : InputParameters
21 312 : FiniteStrainUObasedCP::validParams()
22 : {
23 312 : InputParameters params = ComputeStressBase::validParams();
24 312 : params.addClassDescription("UserObject based Crystal Plasticity system.");
25 624 : params.addParam<Real>("rtol", 1e-6, "Constitutive stress residue relative tolerance");
26 624 : params.addParam<Real>("abs_tol", 1e-6, "Constitutive stress residue absolute tolerance");
27 624 : params.addParam<Real>(
28 624 : "stol", 1e-2, "Constitutive slip system resistance relative residual tolerance");
29 624 : params.addParam<Real>(
30 624 : "zero_tol", 1e-12, "Tolerance for residual check when variable value is zero");
31 624 : params.addParam<unsigned int>("maxiter", 100, "Maximum number of iterations for stress update");
32 624 : params.addParam<unsigned int>(
33 624 : "maxiter_state_variable", 100, "Maximum number of iterations for state variable update");
34 624 : MooseEnum tan_mod_options("exact none", "none"); // Type of read
35 624 : params.addParam<MooseEnum>("tan_mod_type",
36 : tan_mod_options,
37 : "Type of tangent moduli for preconditioner: default elastic");
38 624 : params.addParam<unsigned int>(
39 624 : "maximum_substep_iteration", 1, "Maximum number of substep iteration");
40 624 : params.addParam<bool>("use_line_search", false, "Use line search in constitutive update");
41 624 : params.addParam<Real>("min_line_search_step_size", 0.01, "Minimum line search step size");
42 624 : params.addParam<Real>("line_search_tol", 0.5, "Line search bisection method tolerance");
43 624 : params.addParam<unsigned int>(
44 624 : "line_search_maxiter", 20, "Line search bisection method maximum number of iteration");
45 624 : MooseEnum line_search_method("CUT_HALF BISECTION", "CUT_HALF");
46 624 : params.addParam<MooseEnum>(
47 : "line_search_method", line_search_method, "The method used in line search");
48 624 : 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 624 : 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 624 : 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 624 : 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 312 : return params;
62 312 : }
63 :
64 234 : FiniteStrainUObasedCP::FiniteStrainUObasedCP(const InputParameters & parameters)
65 : : ComputeStressBase(parameters),
66 234 : _num_uo_slip_rates(parameters.get<std::vector<UserObjectName>>("uo_slip_rates").size()),
67 234 : _num_uo_slip_resistances(
68 234 : parameters.get<std::vector<UserObjectName>>("uo_slip_resistances").size()),
69 234 : _num_uo_state_vars(parameters.get<std::vector<UserObjectName>>("uo_state_vars").size()),
70 234 : _num_uo_state_var_evol_rate_comps(
71 234 : parameters.get<std::vector<UserObjectName>>("uo_state_var_evol_rate_comps").size()),
72 468 : _rtol(getParam<Real>("rtol")),
73 468 : _abs_tol(getParam<Real>("abs_tol")),
74 468 : _stol(getParam<Real>("stol")),
75 468 : _zero_tol(getParam<Real>("zero_tol")),
76 468 : _maxiter(getParam<unsigned int>("maxiter")),
77 468 : _maxiterg(getParam<unsigned int>("maxiter_state_variable")),
78 468 : _tan_mod_type(getParam<MooseEnum>("tan_mod_type")),
79 468 : _max_substep_iter(getParam<unsigned int>("maximum_substep_iteration")),
80 468 : _use_line_search(getParam<bool>("use_line_search")),
81 468 : _min_lsrch_step(getParam<Real>("min_line_search_step_size")),
82 468 : _lsrch_tol(getParam<Real>("line_search_tol")),
83 468 : _lsrch_max_iter(getParam<unsigned int>("line_search_maxiter")),
84 468 : _lsrch_method(getParam<MooseEnum>("line_search_method")),
85 234 : _fp(declareProperty<RankTwoTensor>("fp")), // Plastic deformation gradient
86 468 : _fp_old(getMaterialPropertyOld<RankTwoTensor>(
87 : "fp")), // Plastic deformation gradient of previous increment
88 234 : _pk2(declareProperty<RankTwoTensor>("pk2")), // 2nd Piola-Kirchoff Stress
89 468 : _pk2_old(getMaterialPropertyOld<RankTwoTensor>(
90 : "pk2")), // 2nd Piola Kirchoff Stress of previous increment
91 234 : _lag_e(declareProperty<RankTwoTensor>("lage")), // Lagrangian strain
92 234 : _update_rot(declareProperty<RankTwoTensor>(
93 : "update_rot")), // Rotation tensor considering material rotation and crystal orientation
94 468 : _update_rot_old(getMaterialPropertyOld<RankTwoTensor>("update_rot")),
95 234 : _elasticity_tensor_name(_base_name + "elasticity_tensor"),
96 234 : _elasticity_tensor(getMaterialPropertyByName<RankFourTensor>(_elasticity_tensor_name)),
97 468 : _deformation_gradient(getMaterialProperty<RankTwoTensor>("deformation_gradient")),
98 468 : _deformation_gradient_old(getMaterialPropertyOld<RankTwoTensor>("deformation_gradient")),
99 1170 : _crysrot(getMaterialProperty<RankTwoTensor>("crysrot"))
100 : {
101 234 : _err_tol = false;
102 :
103 : _delta_dfgrd.zero();
104 :
105 : // resize the material properties for each userobject
106 234 : _mat_prop_slip_rates.resize(_num_uo_slip_rates);
107 234 : _mat_prop_slip_resistances.resize(_num_uo_slip_resistances);
108 234 : _mat_prop_state_vars.resize(_num_uo_state_vars);
109 234 : _mat_prop_state_vars_old.resize(_num_uo_state_vars);
110 234 : _mat_prop_state_var_evol_rate_comps.resize(_num_uo_state_var_evol_rate_comps);
111 :
112 : // resize the flow direction
113 234 : _flow_direction.resize(_num_uo_slip_rates);
114 :
115 : // resize local state variables
116 234 : _state_vars_old.resize(_num_uo_state_vars);
117 234 : _state_vars_old_stored.resize(_num_uo_state_vars);
118 234 : _state_vars_prev.resize(_num_uo_state_vars);
119 :
120 : // resize user objects
121 234 : _uo_slip_rates.resize(_num_uo_slip_rates);
122 234 : _uo_slip_resistances.resize(_num_uo_slip_resistances);
123 234 : _uo_state_vars.resize(_num_uo_state_vars);
124 234 : _uo_state_var_evol_rate_comps.resize(_num_uo_state_var_evol_rate_comps);
125 :
126 : // assign the user objects
127 468 : for (unsigned int i = 0; i < _num_uo_slip_rates; ++i)
128 : {
129 234 : _uo_slip_rates[i] = &getUserObjectByName<CrystalPlasticitySlipRate>(
130 234 : parameters.get<std::vector<UserObjectName>>("uo_slip_rates")[i]);
131 234 : _mat_prop_slip_rates[i] = &declareProperty<std::vector<Real>>(
132 468 : parameters.get<std::vector<UserObjectName>>("uo_slip_rates")[i]);
133 234 : _flow_direction[i] = &declareProperty<std::vector<RankTwoTensor>>(
134 468 : parameters.get<std::vector<UserObjectName>>("uo_slip_rates")[i] + "_flow_direction");
135 : }
136 :
137 468 : for (unsigned int i = 0; i < _num_uo_slip_resistances; ++i)
138 : {
139 234 : _uo_slip_resistances[i] = &getUserObjectByName<CrystalPlasticitySlipResistance>(
140 234 : parameters.get<std::vector<UserObjectName>>("uo_slip_resistances")[i]);
141 234 : _mat_prop_slip_resistances[i] = &declareProperty<std::vector<Real>>(
142 468 : parameters.get<std::vector<UserObjectName>>("uo_slip_resistances")[i]);
143 : }
144 :
145 468 : for (unsigned int i = 0; i < _num_uo_state_vars; ++i)
146 : {
147 234 : _uo_state_vars[i] = &getUserObjectByName<CrystalPlasticityStateVariable>(
148 234 : parameters.get<std::vector<UserObjectName>>("uo_state_vars")[i]);
149 234 : _mat_prop_state_vars[i] = &declareProperty<std::vector<Real>>(
150 468 : parameters.get<std::vector<UserObjectName>>("uo_state_vars")[i]);
151 234 : _mat_prop_state_vars_old[i] = &getMaterialPropertyOld<std::vector<Real>>(
152 468 : parameters.get<std::vector<UserObjectName>>("uo_state_vars")[i]);
153 : }
154 :
155 468 : for (unsigned int i = 0; i < _num_uo_state_var_evol_rate_comps; ++i)
156 : {
157 234 : _uo_state_var_evol_rate_comps[i] = &getUserObjectByName<CrystalPlasticityStateVarRateComponent>(
158 234 : parameters.get<std::vector<UserObjectName>>("uo_state_var_evol_rate_comps")[i]);
159 234 : _mat_prop_state_var_evol_rate_comps[i] = &declareProperty<std::vector<Real>>(
160 468 : parameters.get<std::vector<UserObjectName>>("uo_state_var_evol_rate_comps")[i]);
161 : }
162 :
163 234 : _substep_dt = 0.0;
164 234 : }
165 :
166 : void
167 6208 : FiniteStrainUObasedCP::initQpStatefulProperties()
168 : {
169 12416 : for (unsigned int i = 0; i < _num_uo_slip_rates; ++i)
170 : {
171 6208 : (*_mat_prop_slip_rates[i])[_qp].resize(_uo_slip_rates[i]->variableSize());
172 6208 : (*_flow_direction[i])[_qp].resize(_uo_slip_rates[i]->variableSize());
173 : }
174 :
175 12416 : for (unsigned int i = 0; i < _num_uo_slip_resistances; ++i)
176 6208 : (*_mat_prop_slip_resistances[i])[_qp].resize(_uo_slip_resistances[i]->variableSize());
177 :
178 12416 : for (unsigned int i = 0; i < _num_uo_state_vars; ++i)
179 : {
180 6208 : (*_mat_prop_state_vars[i])[_qp].resize(_uo_state_vars[i]->variableSize());
181 6208 : _state_vars_old[i].resize(_uo_state_vars[i]->variableSize());
182 6208 : _state_vars_old_stored[i].resize(_uo_state_vars[i]->variableSize());
183 6208 : _state_vars_prev[i].resize(_uo_state_vars[i]->variableSize());
184 : }
185 :
186 12416 : for (unsigned int i = 0; i < _num_uo_state_var_evol_rate_comps; ++i)
187 6208 : (*_mat_prop_state_var_evol_rate_comps[i])[_qp].resize(
188 6208 : _uo_state_var_evol_rate_comps[i]->variableSize());
189 :
190 6208 : _stress[_qp].zero();
191 6208 : _pk2[_qp].zero();
192 6208 : _lag_e[_qp].zero();
193 :
194 6208 : _fp[_qp].setToIdentity();
195 6208 : _update_rot[_qp].setToIdentity();
196 :
197 12416 : for (unsigned int i = 0; i < _num_uo_state_vars; ++i)
198 : // Initializes slip system related properties
199 6208 : _uo_state_vars[i]->initSlipSysProps((*_mat_prop_state_vars[i])[_qp], _q_point[_qp]);
200 6208 : }
201 :
202 : /**
203 : * Solves stress residual equation using NR.
204 : * Updates slip system resistances iteratively.
205 : */
206 : void
207 312366 : FiniteStrainUObasedCP::computeQpStress()
208 : {
209 : // Userobject based crystal plasticity does not support face/boundary material property
210 : // calculation.
211 312366 : 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 309806 : _dfgrd_tmp_old = _deformation_gradient_old[_qp];
219 309806 : if (_dfgrd_tmp_old.det() == 0)
220 0 : _dfgrd_tmp_old.addIa(1.0);
221 :
222 309806 : _delta_dfgrd = _deformation_gradient[_qp] - _dfgrd_tmp_old;
223 :
224 : // Saves the old stateful properties that are modified during sub stepping
225 619612 : for (unsigned int i = 0; i < _num_uo_state_vars; ++i)
226 309806 : _state_vars_old[i] = (*_mat_prop_state_vars_old[i])[_qp];
227 :
228 619612 : for (unsigned int i = 0; i < _num_uo_slip_rates; ++i)
229 309806 : _uo_slip_rates[i]->calcFlowDirection(_qp, (*_flow_direction[i])[_qp]);
230 :
231 : do
232 : {
233 311246 : _err_tol = false;
234 :
235 311246 : preSolveQp();
236 :
237 311246 : _substep_dt = _dt / num_substep;
238 :
239 624366 : for (unsigned int istep = 0; istep < num_substep; ++istep)
240 : {
241 314606 : _dfgrd_tmp = (static_cast<Real>(istep) + 1) / num_substep * _delta_dfgrd + _dfgrd_tmp_old;
242 :
243 314606 : solveQp();
244 :
245 314606 : if (_err_tol)
246 : {
247 1486 : substep_iter++;
248 1486 : num_substep *= 2;
249 1486 : break;
250 : }
251 : }
252 :
253 311246 : if (substep_iter > _max_substep_iter && _err_tol)
254 46 : throw MooseException("FiniteStrainUObasedCP: Constitutive failure.");
255 311200 : } while (_err_tol);
256 :
257 309760 : postSolveQp();
258 : }
259 :
260 : void
261 311246 : FiniteStrainUObasedCP::preSolveQp()
262 : {
263 622492 : for (unsigned int i = 0; i < _num_uo_state_vars; ++i)
264 311246 : (*_mat_prop_state_vars[i])[_qp] = _state_vars_old_stored[i] = _state_vars_old[i];
265 :
266 311246 : _pk2[_qp] = _pk2_old[_qp];
267 311246 : _fp_old_inv = _fp_old[_qp].inverse();
268 311246 : }
269 :
270 : void
271 314606 : FiniteStrainUObasedCP::solveQp()
272 : {
273 314606 : preSolveStatevar();
274 314606 : solveStatevar();
275 314606 : if (_err_tol)
276 : return;
277 :
278 313120 : postSolveStatevar();
279 : }
280 :
281 : void
282 309760 : FiniteStrainUObasedCP::postSolveQp()
283 : {
284 309760 : _stress[_qp] = _fe * _pk2[_qp] * _fe.transpose() / _fe.det();
285 :
286 : // Calculate jacobian for preconditioner
287 309760 : calcTangentModuli();
288 :
289 309760 : RankTwoTensor iden(RankTwoTensor::initIdentity);
290 :
291 309760 : _lag_e[_qp] = _deformation_gradient[_qp].transpose() * _deformation_gradient[_qp] - iden;
292 309760 : _lag_e[_qp] = _lag_e[_qp] * 0.5;
293 :
294 309760 : RankTwoTensor rot;
295 : // Calculate material rotation
296 309760 : _deformation_gradient[_qp].getRUDecompositionRotation(rot);
297 309760 : _update_rot[_qp] = rot * _crysrot[_qp];
298 309760 : }
299 :
300 : void
301 314606 : FiniteStrainUObasedCP::preSolveStatevar()
302 : {
303 629212 : for (unsigned int i = 0; i < _num_uo_state_vars; ++i)
304 314606 : (*_mat_prop_state_vars[i])[_qp] = _state_vars_old_stored[i];
305 :
306 629212 : for (unsigned int i = 0; i < _num_uo_slip_resistances; ++i)
307 314606 : _uo_slip_resistances[i]->calcSlipResistance(_qp, (*_mat_prop_slip_resistances[i])[_qp]);
308 :
309 314606 : _fp_inv = _fp_old_inv;
310 314606 : }
311 :
312 : void
313 314606 : 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 634478 : while (iter_flag && iterg < _maxiterg)
321 : {
322 321358 : preSolveStress();
323 321358 : solveStress();
324 321358 : if (_err_tol)
325 : return;
326 :
327 319872 : postSolveStress();
328 :
329 : // Update slip system resistance and state variable
330 319872 : updateSlipSystemResistanceAndStateVariable();
331 :
332 319872 : if (_err_tol)
333 : return;
334 :
335 319872 : iter_flag = isStateVariablesConverged();
336 319872 : iterg++;
337 : }
338 :
339 313120 : 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 319872 : FiniteStrainUObasedCP::isStateVariablesConverged()
350 : {
351 : Real diff;
352 :
353 632992 : for (unsigned int i = 0; i < _num_uo_state_vars; ++i)
354 : {
355 319872 : unsigned int n = (*_mat_prop_state_vars[i])[_qp].size();
356 4644096 : for (unsigned j = 0; j < n; j++)
357 : {
358 4330976 : diff = std::abs((*_mat_prop_state_vars[i])[_qp][j] - _state_vars_prev[i][j]);
359 :
360 4330976 : if (std::abs(_state_vars_old_stored[i][j]) < _zero_tol && diff > _zero_tol)
361 : return true;
362 4330976 : if (std::abs(_state_vars_old_stored[i][j]) > _zero_tol &&
363 4330976 : diff > _stol * std::abs(_state_vars_old_stored[i][j]))
364 : return true;
365 : }
366 : }
367 : return false;
368 : }
369 :
370 : void
371 313120 : FiniteStrainUObasedCP::postSolveStatevar()
372 : {
373 626240 : for (unsigned int i = 0; i < _num_uo_state_vars; ++i)
374 313120 : _state_vars_old_stored[i] = (*_mat_prop_state_vars[i])[_qp];
375 :
376 313120 : _fp_old_inv = _fp_inv;
377 313120 : }
378 :
379 : void
380 321358 : FiniteStrainUObasedCP::preSolveStress()
381 : {
382 321358 : }
383 :
384 : void
385 321358 : FiniteStrainUObasedCP::solveStress()
386 : {
387 : unsigned int iter = 0;
388 321358 : RankTwoTensor dpk2;
389 : Real rnorm, rnorm0, rnorm_prev;
390 :
391 : // Calculate stress residual
392 321358 : calcResidJacob();
393 321358 : 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 321358 : rnorm = _resid.L2norm();
405 : rnorm0 = rnorm;
406 :
407 : // Check for stress residual tolerance
408 1158230 : while (rnorm > _rtol * rnorm0 && rnorm0 > _abs_tol && iter < _maxiter)
409 : {
410 : // Calculate stress increment
411 838358 : dpk2 = -_jac.invSymm() * _resid;
412 838358 : _pk2[_qp] = _pk2[_qp] + dpk2;
413 :
414 838358 : calcResidJacob();
415 :
416 838358 : 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 836872 : rnorm = _resid.L2norm();
429 :
430 836872 : 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 836872 : if (_use_line_search)
440 12320 : rnorm = _resid.L2norm();
441 :
442 836872 : iter++;
443 : }
444 :
445 319872 : 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 319872 : FiniteStrainUObasedCP::postSolveStress()
465 : {
466 319872 : _fp[_qp] = _fp_inv.inverse();
467 319872 : }
468 :
469 : void
470 319872 : FiniteStrainUObasedCP::updateSlipSystemResistanceAndStateVariable()
471 : {
472 639744 : for (unsigned int i = 0; i < _num_uo_state_vars; ++i)
473 319872 : _state_vars_prev[i] = (*_mat_prop_state_vars[i])[_qp];
474 :
475 639744 : for (unsigned int i = 0; i < _num_uo_state_var_evol_rate_comps; ++i)
476 319872 : _uo_state_var_evol_rate_comps[i]->calcStateVariableEvolutionRateComponent(
477 319872 : _qp, (*_mat_prop_state_var_evol_rate_comps[i])[_qp]);
478 :
479 639744 : for (unsigned int i = 0; i < _num_uo_state_vars; ++i)
480 : {
481 319872 : if (!_uo_state_vars[i]->updateStateVariable(
482 319872 : _qp, _substep_dt, (*_mat_prop_state_vars[i])[_qp], _state_vars_old_stored[i]))
483 0 : _err_tol = true;
484 : }
485 :
486 639744 : for (unsigned int i = 0; i < _num_uo_slip_resistances; ++i)
487 319872 : _uo_slip_resistances[i]->calcSlipResistance(_qp, (*_mat_prop_slip_resistances[i])[_qp]);
488 319872 : }
489 :
490 : // Calculates stress residual equation and jacobian
491 : void
492 1159716 : FiniteStrainUObasedCP::calcResidJacob()
493 : {
494 1159716 : calcResidual();
495 1159716 : if (_err_tol)
496 : return;
497 1158230 : calcJacobian();
498 : }
499 :
500 : void
501 1159716 : FiniteStrainUObasedCP::getSlipRates()
502 : {
503 2317946 : for (unsigned int i = 0; i < _num_uo_slip_rates; ++i)
504 : {
505 1159716 : if (!_uo_slip_rates[i]->calcSlipRate(_qp, _substep_dt, (*_mat_prop_slip_rates[i])[_qp]))
506 : {
507 1486 : _err_tol = true;
508 1486 : return;
509 : }
510 : }
511 : }
512 :
513 : void
514 1159716 : FiniteStrainUObasedCP::calcResidual()
515 : {
516 1159716 : RankTwoTensor iden(RankTwoTensor::initIdentity), ce, ee, ce_pk2, eqv_slip_incr, pk2_new;
517 :
518 1159716 : getSlipRates();
519 1159716 : if (_err_tol)
520 : return;
521 :
522 2316460 : for (unsigned int i = 0; i < _num_uo_slip_rates; ++i)
523 16906742 : for (unsigned int j = 0; j < _uo_slip_rates[i]->variableSize(); ++j)
524 : eqv_slip_incr +=
525 15748512 : (*_flow_direction[i])[_qp][j] * (*_mat_prop_slip_rates[i])[_qp][j] * _substep_dt;
526 :
527 1158230 : eqv_slip_incr = iden - eqv_slip_incr;
528 1158230 : _fp_inv = _fp_old_inv * eqv_slip_incr;
529 1158230 : _fe = _dfgrd_tmp * _fp_inv;
530 :
531 1158230 : ce = _fe.transpose() * _fe;
532 1158230 : ee = ce - iden;
533 1158230 : ee *= 0.5;
534 :
535 1158230 : pk2_new = _elasticity_tensor[_qp] * ee;
536 :
537 1158230 : _resid = _pk2[_qp] - pk2_new;
538 : }
539 :
540 : void
541 1158230 : FiniteStrainUObasedCP::calcJacobian()
542 : {
543 1158230 : RankFourTensor dfedfpinv, deedfe, dfpinvdpk2;
544 :
545 4632920 : for (const auto i : make_range(Moose::dim))
546 13898760 : for (const auto j : make_range(Moose::dim))
547 41696280 : for (const auto k : make_range(Moose::dim))
548 31272210 : dfedfpinv(i, j, k, j) = _dfgrd_tmp(i, k);
549 :
550 4632920 : for (const auto i : make_range(Moose::dim))
551 13898760 : for (const auto j : make_range(Moose::dim))
552 41696280 : for (const auto k : make_range(Moose::dim))
553 : {
554 31272210 : deedfe(i, j, k, i) = deedfe(i, j, k, i) + _fe(k, j) * 0.5;
555 31272210 : deedfe(i, j, k, j) = deedfe(i, j, k, j) + _fe(k, i) * 0.5;
556 : }
557 :
558 2316460 : for (unsigned int i = 0; i < _num_uo_slip_rates; ++i)
559 : {
560 1158230 : unsigned int nss = _uo_slip_rates[i]->variableSize();
561 1158230 : std::vector<RankTwoTensor> dtaudpk2(nss), dfpinvdslip(nss);
562 : std::vector<Real> dslipdtau;
563 1158230 : dslipdtau.resize(nss);
564 1158230 : _uo_slip_rates[i]->calcSlipRateDerivative(_qp, _substep_dt, dslipdtau);
565 16906742 : for (unsigned int j = 0; j < nss; j++)
566 : {
567 15748512 : dtaudpk2[j] = (*_flow_direction[i])[_qp][j];
568 15748512 : dfpinvdslip[j] = -_fp_old_inv * (*_flow_direction[i])[_qp][j];
569 15748512 : dfpinvdpk2 += (dfpinvdslip[j] * dslipdtau[j] * _substep_dt).outerProduct(dtaudpk2[j]);
570 : }
571 : }
572 : _jac =
573 2316460 : RankFourTensor::IdentityFour() - (_elasticity_tensor[_qp] * deedfe * dfedfpinv * dfpinvdpk2);
574 1158230 : }
575 :
576 : void
577 309760 : FiniteStrainUObasedCP::calcTangentModuli()
578 : {
579 309760 : switch (_tan_mod_type)
580 : {
581 309760 : case 0:
582 309760 : elastoPlasticTangentModuli();
583 309760 : break;
584 0 : default:
585 0 : elasticTangentModuli();
586 : }
587 309760 : }
588 :
589 : void
590 309760 : FiniteStrainUObasedCP::elastoPlasticTangentModuli()
591 : {
592 309760 : RankFourTensor tan_mod;
593 309760 : RankTwoTensor pk2fet, fepk2;
594 309760 : RankFourTensor deedfe, dsigdpk2dfe, dfedf;
595 :
596 : // Fill in the matrix stiffness material property
597 1239040 : for (const auto i : make_range(Moose::dim))
598 3717120 : for (const auto j : make_range(Moose::dim))
599 11151360 : for (const auto k : make_range(Moose::dim))
600 : {
601 8363520 : deedfe(i, j, k, i) = deedfe(i, j, k, i) + _fe(k, j) * 0.5;
602 8363520 : deedfe(i, j, k, j) = deedfe(i, j, k, j) + _fe(k, i) * 0.5;
603 : }
604 :
605 : usingTensorIndices(i_, j_, k_, l_);
606 309760 : dsigdpk2dfe = _fe.times<i_, k_, j_, l_>(_fe) * _elasticity_tensor[_qp] * deedfe;
607 :
608 309760 : pk2fet = _pk2[_qp] * _fe.transpose();
609 309760 : fepk2 = _fe * _pk2[_qp];
610 :
611 1239040 : for (const auto i : make_range(Moose::dim))
612 3717120 : for (const auto j : make_range(Moose::dim))
613 11151360 : for (const auto l : make_range(Moose::dim))
614 : {
615 8363520 : tan_mod(i, j, i, l) += pk2fet(l, j);
616 8363520 : tan_mod(i, j, j, l) += fepk2(i, l);
617 : }
618 :
619 309760 : tan_mod += dsigdpk2dfe;
620 :
621 309760 : Real je = _fe.det();
622 309760 : if (je > 0.0)
623 309760 : tan_mod /= je;
624 :
625 1239040 : for (const auto i : make_range(Moose::dim))
626 3717120 : for (const auto j : make_range(Moose::dim))
627 11151360 : for (const auto l : make_range(Moose::dim))
628 8363520 : dfedf(i, j, i, l) = _fp_inv(l, j);
629 :
630 309760 : _Jacobian_mult[_qp] = tan_mod * dfedf;
631 309760 : }
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 : }
|
