www.mooseframework.org
FiniteStrainCrystalPlasticity.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
3 /* */
4 /* All contents are licensed under LGPL V2.1 */
5 /* See LICENSE for full restrictions */
6 /****************************************************************/
8 #include "petscblaslapack.h"
9 #include "libmesh/utility.h"
10 
11 #include <fstream>
12 
13 template <>
14 InputParameters
16 {
17  InputParameters params = validParams<ComputeStressBase>();
18  params.addClassDescription(
19  "Crystal Plasticity base class: FCC system with power law flow rule implemented");
20  params.addRequiredParam<int>("nss", "Number of slip systems");
21  params.addParam<std::vector<Real>>("gprops", "Initial values of slip system resistances");
22  params.addParam<std::vector<Real>>("hprops", "Hardening properties");
23  params.addParam<std::vector<Real>>("flowprops", "Parameters used in slip rate equations");
24  params.addRequiredParam<FileName>("slip_sys_file_name",
25  "Name of the file containing the slip system");
26  params.addParam<FileName>(
27  "slip_sys_res_prop_file_name",
28  "",
29  "Name of the file containing the initial values of slip system resistances");
30  params.addParam<FileName>(
31  "slip_sys_flow_prop_file_name",
32  "",
33  "Name of the file containing the values of slip rate equation parameters");
34  params.addParam<FileName>(
35  "slip_sys_hard_prop_file_name",
36  "",
37  "Name of the file containing the values of hardness evolution parameters");
38  params.addParam<Real>("rtol", 1e-6, "Constitutive stress residue relative tolerance");
39  params.addParam<Real>("abs_tol", 1e-6, "Constitutive stress residue absolute tolerance");
40  params.addParam<Real>("gtol", 1e2, "Constitutive slip system resistance residual tolerance");
41  params.addParam<Real>("slip_incr_tol", 2e-2, "Maximum allowable slip in an increment");
42  params.addParam<unsigned int>("maxiter", 100, "Maximum number of iterations for stress update");
43  params.addParam<unsigned int>(
44  "maxitergss", 100, "Maximum number of iterations for slip system resistance update");
45  params.addParam<unsigned int>(
46  "num_slip_sys_flowrate_props",
47  2,
48  "Number of flow rate properties for a slip system"); // Used for reading flow rate parameters
49  params.addParam<UserObjectName>("read_prop_user_object",
50  "The ElementReadPropertyFile "
51  "GeneralUserObject to read element "
52  "specific property values from file");
53  MooseEnum tan_mod_options("exact none", "none"); // Type of read
54  params.addParam<MooseEnum>("tan_mod_type",
55  tan_mod_options,
56  "Type of tangent moduli for preconditioner: default elastic");
57  MooseEnum intvar_read_options("slip_sys_file slip_sys_res_file none", "none");
58  params.addParam<MooseEnum>(
59  "intvar_read_type",
60  intvar_read_options,
61  "Read from options for initial value of internal variables: Default from .i file");
62  params.addParam<unsigned int>("num_slip_sys_props",
63  0,
64  "Number of slip system specific properties provided in the file "
65  "containing slip system normals and directions");
66  params.addParam<bool>(
67  "gen_random_stress_flag",
68  false,
69  "Flag to generate random stress to perform time cutback on constitutive failure");
70  params.addParam<bool>("input_random_scaling_var",
71  false,
72  "Flag to input scaling variable: _Cijkl(0,0,0,0) when false");
73  params.addParam<Real>("random_scaling_var",
74  1e9,
75  "Random scaling variable: Large value can cause non-positive definiteness");
76  params.addParam<unsigned int>(
77  "random_seed",
78  2000,
79  "Random integer used to generate random stress when constitutive failure occurs");
80  params.addParam<unsigned int>(
81  "maximum_substep_iteration", 1, "Maximum number of substep iteration");
82  params.addParam<bool>("use_line_search", false, "Use line search in constitutive update");
83  params.addParam<Real>("min_line_search_step_size", 0.01, "Minimum line search step size");
84  params.addParam<Real>("line_search_tol", 0.5, "Line search bisection method tolerance");
85  params.addParam<unsigned int>(
86  "line_search_maxiter", 20, "Line search bisection method maximum number of iteration");
87  MooseEnum line_search_method("CUT_HALF BISECTION", "CUT_HALF");
88  params.addParam<MooseEnum>(
89  "line_search_method", line_search_method, "The method used in line search");
90 
91  return params;
92 }
93 
95  : ComputeStressBase(parameters),
96  _nss(getParam<int>("nss")),
97  _gprops(getParam<std::vector<Real>>("gprops")),
98  _hprops(getParam<std::vector<Real>>("hprops")),
99  _flowprops(getParam<std::vector<Real>>("flowprops")),
100  _slip_sys_file_name(getParam<FileName>("slip_sys_file_name")),
101  _slip_sys_res_prop_file_name(getParam<FileName>("slip_sys_res_prop_file_name")),
102  _slip_sys_flow_prop_file_name(getParam<FileName>("slip_sys_flow_prop_file_name")),
103  _slip_sys_hard_prop_file_name(getParam<FileName>("slip_sys_hard_prop_file_name")),
104  _rtol(getParam<Real>("rtol")),
105  _abs_tol(getParam<Real>("abs_tol")),
106  _gtol(getParam<Real>("gtol")),
107  _slip_incr_tol(getParam<Real>("slip_incr_tol")),
108  _maxiter(getParam<unsigned int>("maxiter")),
109  _maxiterg(getParam<unsigned int>("maxitergss")),
110  _num_slip_sys_flowrate_props(getParam<unsigned int>("num_slip_sys_flowrate_props")),
111  _tan_mod_type(getParam<MooseEnum>("tan_mod_type")),
112  _intvar_read_type(getParam<MooseEnum>("intvar_read_type")),
113  _num_slip_sys_props(getParam<unsigned int>("num_slip_sys_props")),
114  _gen_rndm_stress_flag(getParam<bool>("gen_random_stress_flag")),
115  _input_rndm_scale_var(getParam<bool>("input_random_scaling_var")),
116  _rndm_scale_var(getParam<Real>("random_scaling_var")),
117  _rndm_seed(getParam<unsigned int>("random_seed")),
118  _max_substep_iter(getParam<unsigned int>("maximum_substep_iteration")),
119  _use_line_search(getParam<bool>("use_line_search")),
120  _min_lsrch_step(getParam<Real>("min_line_search_step_size")),
121  _lsrch_tol(getParam<Real>("line_search_tol")),
122  _lsrch_max_iter(getParam<unsigned int>("line_search_maxiter")),
123  _lsrch_method(getParam<MooseEnum>("line_search_method")),
124  _fp(declareProperty<RankTwoTensor>("fp")), // Plastic deformation gradient
125  _fp_old(getMaterialPropertyOld<RankTwoTensor>(
126  "fp")), // Plastic deformation gradient of previous increment
127  _pk2(declareProperty<RankTwoTensor>("pk2")), // 2nd Piola Kirchoff Stress
128  _pk2_old(getMaterialPropertyOld<RankTwoTensor>(
129  "pk2")), // 2nd Piola Kirchoff Stress of previous increment
130  _lag_e(declareProperty<RankTwoTensor>("lage")), // Lagrangian strain
131  _lag_e_old(
132  getMaterialPropertyOld<RankTwoTensor>("lage")), // Lagrangian strain of previous increment
133  _gss(declareProperty<std::vector<Real>>("gss")), // Slip system resistances
134  _gss_old(getMaterialPropertyOld<std::vector<Real>>(
135  "gss")), // Slip system resistances of previous increment
136  _acc_slip(declareProperty<Real>("acc_slip")), // Accumulated slip
137  _acc_slip_old(
138  getMaterialPropertyOld<Real>("acc_slip")), // Accumulated alip of previous increment
139  _update_rot(declareProperty<RankTwoTensor>(
140  "update_rot")), // Rotation tensor considering material rotation and crystal orientation
141  _deformation_gradient(getMaterialProperty<RankTwoTensor>("deformation_gradient")),
142  _deformation_gradient_old(getMaterialPropertyOld<RankTwoTensor>("deformation_gradient")),
143  _elasticity_tensor(getMaterialProperty<RankFourTensor>("elasticity_tensor")),
144  _crysrot(getMaterialProperty<RankTwoTensor>("crysrot")),
145  _mo(_nss * LIBMESH_DIM),
146  _no(_nss * LIBMESH_DIM),
147  _slip_incr(_nss),
148  _tau(_nss),
149  _dslipdtau(_nss),
150  _s0(_nss),
151  _gss_tmp(_nss),
152  _gss_tmp_old(_nss),
153  _dgss_dsliprate(_nss, _nss)
154 {
155  _err_tol = false;
156 
157  if (_num_slip_sys_props > 0)
159 
160  _pk2_tmp.zero();
161  _delta_dfgrd.zero();
162 
163  _first_step_iter = false;
164  _last_step_iter = false;
165  // Initialize variables in the first iteration of substepping
166  _first_substep = true;
167 
168  _read_from_slip_sys_file = false;
169  if (_intvar_read_type == "slip_sys_file")
171 
172  if (_read_from_slip_sys_file && !(_num_slip_sys_props > 0))
173  mooseError("Crystal Plasticity Error: Specify number of internal variable's initial values to "
174  "be read from slip system file");
175 
176  getSlipSystems();
177 
178  RankTwoTensor::initRandom(_rndm_seed);
179 }
180 
181 void
183 {
184  _stress[_qp].zero();
185 
186  _fp[_qp].zero();
187  _fp[_qp].addIa(1.0);
188 
189  _pk2[_qp].zero();
190  _acc_slip[_qp] = 0.0;
191  _lag_e[_qp].zero();
192 
193  _update_rot[_qp].zero();
194  _update_rot[_qp].addIa(1.0);
195 
196  initSlipSysProps(); // Initializes slip system related properties
198 }
199 
200 void
202 {
203  switch (_intvar_read_type)
204  {
205  case 0:
207  break;
208  case 1:
210  break;
211  default:
213  }
214 
215  if (_slip_sys_flow_prop_file_name.length() != 0)
217  else
219 
220  if (_slip_sys_hard_prop_file_name.length() != 0)
222  else
224 }
225 
226 void
228 {
229  _gss[_qp].resize(_nss);
230 
231  for (unsigned int i = 0; i < _nss; ++i)
232  _gss[_qp][i] = _slip_sys_props(i);
233 }
234 
235 // Read initial slip system resistances from .txt file. See test.
236 void
238 {
239  _gss[_qp].resize(_nss);
240 
241  MooseUtils::checkFileReadable(_slip_sys_res_prop_file_name);
242 
243  std::ifstream file;
244  file.open(_slip_sys_res_prop_file_name.c_str());
245 
246  for (unsigned int i = 0; i < _nss; ++i)
247  if (!(file >> _gss[_qp][i]))
248  mooseError("Error FiniteStrainCrystalPlasticity: Premature end of slip_sys_res_prop file");
249 
250  file.close();
251 }
252 
253 // Read initial slip system resistances from .i file
254 void
256 {
257  if (_gprops.size() <= 0)
258  mooseError("FiniteStrainCrystalPLasticity: Error in reading slip system resistance properties: "
259  "Specify input in .i file or in slip_sys_res_prop_file or in slip_sys_file");
260 
261  _gss[_qp].resize(_nss, 0.0);
262 
263  unsigned int num_data_grp = 3; // Number of data per group e.g. start_slip_sys, end_slip_sys,
264  // value
265 
266  for (unsigned int i = 0; i < _gprops.size() / num_data_grp; ++i)
267  {
268  Real vs, ve;
269  unsigned int is, ie;
270 
271  vs = _gprops[i * num_data_grp];
272  ve = _gprops[i * num_data_grp + 1];
273 
274  if (vs <= 0 || ve <= 0)
275  mooseError("FiniteStrainCrystalPLasticity: Indices in gss property read must be positive "
276  "integers: is = ",
277  vs,
278  " ie = ",
279  ve);
280 
281  if (vs != floor(vs) || ve != floor(ve))
282  mooseError("FiniteStrainCrystalPLasticity: Error in reading slip system resistances: Values "
283  "specifying start and end number of slip system groups should be integer");
284 
285  is = static_cast<unsigned int>(vs);
286  ie = static_cast<unsigned int>(ve);
287 
288  if (is > ie)
289  mooseError("FiniteStrainCrystalPLasticity: Start index is = ",
290  is,
291  " should be greater than end index ie = ",
292  ie,
293  " in slip system resistance property read");
294 
295  for (unsigned int j = is; j <= ie; ++j)
296  _gss[_qp][j - 1] = _gprops[i * num_data_grp + 2];
297  }
298 
299  for (unsigned int i = 0; i < _nss; ++i)
300  if (_gss[_qp][i] <= 0.0)
301  mooseError("FiniteStrainCrystalPLasticity: Value of resistance for slip system ",
302  i + 1,
303  " non positive");
304 }
305 
306 // Read flow rate parameters from .txt file. See test.
307 void
309 {
310  _a0.resize(_nss);
311  _xm.resize(_nss);
312 
313  MooseUtils::checkFileReadable(_slip_sys_flow_prop_file_name);
314 
315  std::ifstream file;
316  file.open(_slip_sys_flow_prop_file_name.c_str());
317 
318  std::vector<Real> vec;
319  vec.resize(_num_slip_sys_flowrate_props);
320 
321  for (unsigned int i = 0; i < _nss; ++i)
322  {
323  for (unsigned int j = 0; j < _num_slip_sys_flowrate_props; ++j)
324  if (!(file >> vec[j]))
325  mooseError(
326  "Error FiniteStrainCrystalPlasticity: Premature end of slip_sys_flow_rate_param file");
327 
328  _a0(i) = vec[0];
329  _xm(i) = vec[1];
330  }
331 
332  file.close();
333 }
334 
335 // Read flow rate parameters from .i file
336 void
338 {
339  if (_flowprops.size() <= 0)
340  mooseError("FiniteStrainCrystalPLasticity: Error in reading flow rate properties: Specify "
341  "input in .i file or a slip_sys_flow_prop_file_name");
342 
343  _a0.resize(_nss);
344  _xm.resize(_nss);
345 
346  unsigned int num_data_grp = 2 + _num_slip_sys_flowrate_props; // Number of data per group e.g.
347  // start_slip_sys, end_slip_sys,
348  // value1, value2, ..
349 
350  for (unsigned int i = 0; i < _flowprops.size() / num_data_grp; ++i)
351  {
352  Real vs, ve;
353  unsigned int is, ie;
354 
355  vs = _flowprops[i * num_data_grp];
356  ve = _flowprops[i * num_data_grp + 1];
357 
358  if (vs <= 0 || ve <= 0)
359  mooseError("FiniteStrainCrystalPLasticity: Indices in flow rate parameter read must be "
360  "positive integers: is = ",
361  vs,
362  " ie = ",
363  ve);
364 
365  if (vs != floor(vs) || ve != floor(ve))
366  mooseError("FiniteStrainCrystalPLasticity: Error in reading flow props: Values specifying "
367  "start and end number of slip system groups should be integer");
368 
369  is = static_cast<unsigned int>(vs);
370  ie = static_cast<unsigned int>(ve);
371 
372  if (is > ie)
373  mooseError("FiniteStrainCrystalPLasticity: Start index is = ",
374  is,
375  " should be greater than end index ie = ",
376  ie,
377  " in flow rate parameter read");
378 
379  for (unsigned int j = is; j <= ie; ++j)
380  {
381  _a0(j - 1) = _flowprops[i * num_data_grp + 2];
382  _xm(j - 1) = _flowprops[i * num_data_grp + 3];
383  }
384  }
385 
386  for (unsigned int i = 0; i < _nss; ++i)
387  {
388  if (!(_a0(i) > 0.0 && _xm(i) > 0.0))
389  {
390  mooseWarning(
391  "FiniteStrainCrystalPlasticity: Non-positive flow rate parameters ", _a0(i), ",", _xm(i));
392  break;
393  }
394  }
395 }
396 
397 // Read hardness parameters from .txt file
398 void
400 {
401 }
402 
403 // Read hardness parameters from .i file
404 void
406 {
407  if (_hprops.size() <= 0)
408  mooseError("FiniteStrainCrystalPLasticity: Error in reading hardness properties: Specify input "
409  "in .i file or a slip_sys_hard_prop_file_name");
410 
411  _r = _hprops[0];
412  _h0 = _hprops[1];
413  _tau_init = _hprops[2];
414  _tau_sat = _hprops[3];
415 }
416 
417 // Read slip systems from file
418 void
420 {
421  Real vec[LIBMESH_DIM];
422  std::ifstream fileslipsys;
423 
424  MooseUtils::checkFileReadable(_slip_sys_file_name);
425 
426  fileslipsys.open(_slip_sys_file_name.c_str());
427 
428  for (unsigned int i = 0; i < _nss; ++i)
429  {
430  // Read the slip normal
431  for (unsigned int j = 0; j < LIBMESH_DIM; ++j)
432  if (!(fileslipsys >> vec[j]))
433  mooseError("Crystal Plasticity Error: Premature end of file reading slip system file \n");
434 
435  // Normalize the vectors
436  Real mag;
437  mag = Utility::pow<2>(vec[0]) + Utility::pow<2>(vec[1]) + Utility::pow<2>(vec[2]);
438  mag = std::sqrt(mag);
439 
440  for (unsigned j = 0; j < LIBMESH_DIM; ++j)
441  _no(i * LIBMESH_DIM + j) = vec[j] / mag;
442 
443  // Read the slip direction
444  for (unsigned int j = 0; j < LIBMESH_DIM; ++j)
445  if (!(fileslipsys >> vec[j]))
446  mooseError("Crystal Plasticity Error: Premature end of file reading slip system file \n");
447 
448  // Normalize the vectors
449  mag = Utility::pow<2>(vec[0]) + Utility::pow<2>(vec[1]) + Utility::pow<2>(vec[2]);
450  mag = std::sqrt(mag);
451 
452  for (unsigned int j = 0; j < LIBMESH_DIM; ++j)
453  _mo(i * LIBMESH_DIM + j) = vec[j] / mag;
454 
455  mag = 0.0;
456  for (unsigned int j = 0; j < LIBMESH_DIM; ++j)
457  mag += _mo(i * LIBMESH_DIM + j) * _no(i * LIBMESH_DIM + j);
458 
459  if (std::abs(mag) > 1e-8)
460  mooseError(
461  "Crystal Plasicity Error: Slip direction and normal not orthonormal, System number = ",
462  i,
463  "\n");
464 
466  for (unsigned int j = 0; j < _num_slip_sys_props; ++j)
467  if (!(fileslipsys >> _slip_sys_props(i * _num_slip_sys_props + j)))
468  mooseError("Crystal Plasticity Error: Premature end of file reading slip system file - "
469  "check in slip system file read input options/values\n");
470  }
471 
472  fileslipsys.close();
473 }
474 
475 // Initialize addtional stateful material properties
476 void
478 {
479 }
480 
485 void
487 {
488  unsigned int substep_iter = 1; // Depth of substepping; Limited to maximum substep iteration
489  unsigned int num_substep = 1; // Calculated from substep_iter as 2^substep_iter
490  Real dt_original = _dt; // Stores original _dt; Reset at the end of solve
491  _first_substep = true; // Initialize variables at substep_iter = 1
492 
493  if (_max_substep_iter > 1)
494  {
496  if (_dfgrd_tmp_old.det() == 0)
497  _dfgrd_tmp_old.addIa(1.0);
498 
500  _err_tol = true; // Indicator to continue substepping
501  }
502 
503  // Substepping loop
504  while (_err_tol && _max_substep_iter > 1)
505  {
506  _dt = dt_original / num_substep;
507 
508  for (unsigned int istep = 0; istep < num_substep; ++istep)
509  {
510  _first_step_iter = false;
511  if (istep == 0)
512  _first_step_iter = true;
513 
514  _last_step_iter = false;
515  if (istep == num_substep - 1)
516  _last_step_iter = true;
517 
518  _dfgrd_scale_factor = (static_cast<Real>(istep) + 1) / num_substep;
519 
520  preSolveQp();
521  solveQp();
522 
523  if (_err_tol)
524  {
525  substep_iter++;
526  num_substep *= 2;
527  break;
528  }
529  }
530 
531  _first_substep = false; // Prevents reinitialization
532  _dt = dt_original; // Resets dt
533 
534 #ifdef DEBUG
535  if (substep_iter > _max_substep_iter)
536  mooseWarning("FiniteStrainCrystalPlasticity: Failure with substepping");
537 #endif
538 
539  if (!_err_tol || substep_iter > _max_substep_iter)
540  postSolveQp(); // Evaluate variables after successful solve or indicate failure
541  }
542 
543  // No substepping
544  if (_max_substep_iter == 1)
545  {
546  preSolveQp();
547  solveQp();
548  postSolveQp();
549  }
550 }
551 
552 void
554 {
555  // Initialize variable
556  if (_first_substep)
557  {
558  _Jacobian_mult[_qp].zero(); // Initializes jacobian for preconditioner
560  }
561 
562  if (_max_substep_iter == 1)
563  _dfgrd_tmp = _deformation_gradient[_qp]; // Without substepping
564  else
566 
567  _err_tol = false;
568 }
569 
570 void
572 {
574  solveStatevar();
575  if (_err_tol)
576  return;
578 }
579 
580 void
582 {
583  if (_err_tol)
584  {
585  _err_tol = false;
587  {
589  _rndm_scale_var = _elasticity_tensor[_qp](0, 0, 0, 0);
590 
591  _stress[_qp] = RankTwoTensor::genRandomSymmTensor(_rndm_scale_var, 1.0);
592  }
593  else
594  mooseError("FiniteStrainCrystalPlasticity: Constitutive failure");
595  }
596  else
597  {
598  _stress[_qp] = _fe * _pk2[_qp] * _fe.transpose() / _fe.det();
599 
600  _Jacobian_mult[_qp] += calcTangentModuli(); // Calculate jacobian for preconditioner
601 
602  RankTwoTensor iden;
603  iden.addIa(1.0);
604 
605  _lag_e[_qp] = _deformation_gradient[_qp].transpose() * _deformation_gradient[_qp] - iden;
606  _lag_e[_qp] = _lag_e[_qp] * 0.5;
607 
608  RankTwoTensor rot;
609  rot = get_current_rotation(_deformation_gradient[_qp]); // Calculate material rotation
610  _update_rot[_qp] = rot * _crysrot[_qp];
611  }
612 }
613 
614 void
616 {
617  if (_max_substep_iter == 1) // No substepping
618  {
619  _gss_tmp = _gss_old[_qp];
621  }
622  else
623  {
624  if (_first_step_iter)
625  {
626  _gss_tmp = _gss_tmp_old = _gss_old[_qp];
628  }
629  else
631  }
632 }
633 
634 void
636 {
637  Real gmax, gdiff;
638  unsigned int iterg;
639  std::vector<Real> gss_prev(_nss);
640 
641  gmax = 1.1 * _gtol;
642  iterg = 0;
643 
644  while (gmax > _gtol && iterg < _maxiterg) // Check for slip system resistance update tolerance
645  {
646  preSolveStress();
647  solveStress();
648  if (_err_tol)
649  return;
650  postSolveStress();
651 
652  gss_prev = _gss_tmp;
653 
654  update_slip_system_resistance(); // Update slip system resistance
655 
656  gmax = 0.0;
657  for (unsigned i = 0; i < _nss; ++i)
658  {
659  gdiff = std::abs(gss_prev[i] - _gss_tmp[i]); // Calculate increment size
660 
661  if (gdiff > gmax)
662  gmax = gdiff;
663  }
664  iterg++;
665  }
666 
667  if (iterg == _maxiterg)
668  {
669 #ifdef DEBUG
670  mooseWarning("FiniteStrainCrystalPLasticity: Hardness Integration error gmax", gmax, "\n");
671 #endif
672  _err_tol = true;
673  }
674 }
675 
676 void
678 {
679  if (_max_substep_iter == 1) // No substepping
680  {
681  _gss[_qp] = _gss_tmp;
682  _acc_slip[_qp] = _accslip_tmp;
683  }
684  else
685  {
686  if (_last_step_iter)
687  {
688  _gss[_qp] = _gss_tmp;
689  _acc_slip[_qp] = _accslip_tmp;
690  }
691  else
692  {
695  }
696  }
697 }
698 
699 void
701 {
702  if (_max_substep_iter == 1) // No substepping
703  {
704  _pk2_tmp = _pk2_old[_qp];
705  _fp_old_inv = _fp_old[_qp].inverse();
708  }
709  else
710  {
711  if (_first_step_iter)
712  {
713  _pk2_tmp = _pk2_tmp_old = _pk2_old[_qp];
714  _fp_old_inv = _fp_old[_qp].inverse();
715  }
716  else
718 
721  }
722 }
723 
724 void
726 {
727  unsigned int iter = 0;
728  RankTwoTensor resid, dpk2;
729  RankFourTensor jac;
730  Real rnorm, rnorm0, rnorm_prev;
731 
732  calc_resid_jacob(resid, jac); // Calculate stress residual
733  if (_err_tol)
734  {
735 #ifdef DEBUG
736  mooseWarning(
737  "FiniteStrainCrystalPLasticity: Slip increment exceeds tolerance - Element number ",
738  _current_elem->id(),
739  " Gauss point = ",
740  _qp);
741 #endif
742  return;
743  }
744 
745  rnorm = resid.L2norm();
746  rnorm0 = rnorm;
747 
748  while (rnorm > _rtol * rnorm0 && rnorm0 > _abs_tol &&
749  iter < _maxiter) // Check for stress residual tolerance
750  {
751  dpk2 = -jac.invSymm() * resid; // Calculate stress increment
752  _pk2_tmp = _pk2_tmp + dpk2; // Update stress
753  calc_resid_jacob(resid, jac);
754  internalVariableUpdateNRiteration(); // update _fp_prev_inv
755 
756  if (_err_tol)
757  {
758 #ifdef DEBUG
759  mooseWarning(
760  "FiniteStrainCrystalPLasticity: Slip increment exceeds tolerance - Element number ",
761  _current_elem->id(),
762  " Gauss point = ",
763  _qp);
764 #endif
765  return;
766  }
767 
768  rnorm_prev = rnorm;
769  rnorm = resid.L2norm();
770 
771  if (_use_line_search && rnorm > rnorm_prev && !line_search_update(rnorm_prev, dpk2))
772  {
773 #ifdef DEBUG
774  mooseWarning("FiniteStrainCrystalPLasticity: Failed with line search");
775 #endif
776  _err_tol = true;
777  return;
778  }
779 
780  if (_use_line_search)
781  rnorm = resid.L2norm();
782 
783  iter++;
784  }
785 
786  if (iter >= _maxiter)
787  {
788 #ifdef DEBUG
789  mooseWarning("FiniteStrainCrystalPLasticity: Stress Integration error rmax = ", rnorm);
790 #endif
791  _err_tol = true;
792  }
793 }
794 
795 void
797 {
798  if (_max_substep_iter == 1) // No substepping
799  {
800  _fp[_qp] = _fp_inv.inverse();
801  _pk2[_qp] = _pk2_tmp;
802  }
803  else
804  {
805  if (_last_step_iter)
806  {
807  _fp[_qp] = _fp_inv.inverse();
808  _pk2[_qp] = _pk2_tmp;
809  }
810  else
811  {
814  }
815  }
816 }
817 
818 // Update slip system resistance. Overide to incorporate new slip system resistance laws
819 void
821 {
822  updateGss();
823 }
824 
829 void
831 {
832  DenseVector<Real> hb(_nss);
833  Real qab;
834 
835  Real a = _hprops[4]; // Kalidindi
836 
838  for (unsigned int i = 0; i < _nss; ++i)
839  _accslip_tmp += std::abs(_slip_incr(i));
840 
841  // Real val = std::cosh(_h0 * _accslip_tmp / (_tau_sat - _tau_init)); // Karthik
842  // val = _h0 * std::pow(1.0/val,2.0); // Kalidindi
843 
844  for (unsigned int i = 0; i < _nss; ++i)
845  // hb(i)=val;
846  hb(i) = _h0 * std::pow(std::abs(1.0 - _gss_tmp[i] / _tau_sat), a) *
847  copysign(1.0, 1.0 - _gss_tmp[i] / _tau_sat);
848 
849  for (unsigned int i = 0; i < _nss; ++i)
850  {
851  if (_max_substep_iter == 1) // No substepping
852  _gss_tmp[i] = _gss_old[_qp][i];
853  else
854  _gss_tmp[i] = _gss_tmp_old[i];
855 
856  for (unsigned int j = 0; j < _nss; ++j)
857  {
858  unsigned int iplane, jplane;
859  iplane = i / 3;
860  jplane = j / 3;
861 
862  if (iplane == jplane) // Kalidindi
863  qab = 1.0;
864  else
865  qab = _r;
866 
867  _gss_tmp[i] += qab * hb(j) * std::abs(_slip_incr(j));
868  _dgss_dsliprate(i, j) = qab * hb(j) * copysign(1.0, _slip_incr(j)) * _dt;
869  }
870  }
871 }
872 
873 // Calculates stress residual equation and jacobian
874 void
875 FiniteStrainCrystalPlasticity::calc_resid_jacob(RankTwoTensor & resid, RankFourTensor & jac)
876 {
877  calcResidual(resid);
878  if (_err_tol)
879  return;
880  calcJacobian(jac);
881 }
882 
883 void
885 {
886  RankTwoTensor iden, ce, ee, ce_pk2, eqv_slip_incr, pk2_new;
887 
888  iden.zero();
889  iden.addIa(1.0);
890 
891  _fe = _dfgrd_tmp * _fp_prev_inv; // _fp_inv ==> _fp_prev_inv
892 
893  ce = _fe.transpose() * _fe;
894  ce_pk2 = ce * _pk2_tmp;
895  ce_pk2 = ce_pk2 / _fe.det();
896 
897  // Calculate Schmid tensor and resolved shear stresses
898  for (unsigned int i = 0; i < _nss; ++i)
899  _tau(i) = ce_pk2.doubleContraction(_s0[i]);
900 
901  getSlipIncrements(); // Calculate dslip,dslipdtau
902 
903  if (_err_tol)
904  return;
905 
906  eqv_slip_incr.zero();
907  for (unsigned int i = 0; i < _nss; ++i)
908  eqv_slip_incr += _s0[i] * _slip_incr(i);
909 
910  eqv_slip_incr = iden - eqv_slip_incr;
911  _fp_inv = _fp_old_inv * eqv_slip_incr;
912  _fe = _dfgrd_tmp * _fp_inv;
913 
914  ce = _fe.transpose() * _fe;
915  ee = ce - iden;
916  ee *= 0.5;
917 
918  pk2_new = _elasticity_tensor[_qp] * ee;
919 
920  resid = _pk2_tmp - pk2_new;
921 }
922 
923 void
925 {
926  RankFourTensor dfedfpinv, deedfe, dfpinvdpk2;
927 
928  std::vector<RankTwoTensor> dtaudpk2(_nss), dfpinvdslip(_nss);
929 
930  for (unsigned int i = 0; i < _nss; ++i)
931  {
932  dtaudpk2[i] = _s0[i];
933  dfpinvdslip[i] = -_fp_old_inv * _s0[i];
934  }
935 
936  for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
937  for (unsigned int j = 0; j < LIBMESH_DIM; ++j)
938  for (unsigned int k = 0; k < LIBMESH_DIM; ++k)
939  dfedfpinv(i, j, k, j) = _dfgrd_tmp(i, k);
940 
941  for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
942  for (unsigned int j = 0; j < LIBMESH_DIM; ++j)
943  for (unsigned int k = 0; k < LIBMESH_DIM; ++k)
944  {
945  deedfe(i, j, k, i) = deedfe(i, j, k, i) + _fe(k, j) * 0.5;
946  deedfe(i, j, k, j) = deedfe(i, j, k, j) + _fe(k, i) * 0.5;
947  }
948 
949  for (unsigned int i = 0; i < _nss; ++i)
950  dfpinvdpk2 += (dfpinvdslip[i] * _dslipdtau(i)).outerProduct(dtaudpk2[i]);
951 
952  jac =
953  RankFourTensor::IdentityFour() - (_elasticity_tensor[_qp] * deedfe * dfedfpinv * dfpinvdpk2);
954 }
955 
956 // Calculate slip increment,dslipdtau. Override to modify.
957 void
959 {
960  for (unsigned int i = 0; i < _nss; ++i)
961  {
962  _slip_incr(i) = _a0(i) * std::pow(std::abs(_tau(i) / _gss_tmp[i]), 1.0 / _xm(i)) *
963  copysign(1.0, _tau(i)) * _dt;
964  if (std::abs(_slip_incr(i)) > _slip_incr_tol)
965  {
966  _err_tol = true;
967 #ifdef DEBUG
968  mooseWarning("Maximum allowable slip increment exceeded ", std::abs(_slip_incr(i)));
969 #endif
970  return;
971  }
972  }
973 
974  for (unsigned int i = 0; i < _nss; ++i)
975  _dslipdtau(i) = _a0(i) / _xm(i) *
976  std::pow(std::abs(_tau(i) / _gss_tmp[i]), 1.0 / _xm(i) - 1.0) / _gss_tmp[i] *
977  _dt;
978 }
979 
980 // Calls getMatRot to perform RU factorization of a tensor.
981 RankTwoTensor
983 {
984  return getMatRot(a);
985 }
986 
987 // Performs RU factorization of a tensor
988 RankTwoTensor
990 {
991  RankTwoTensor rot;
992  RankTwoTensor c, diag, evec;
993  PetscScalar cmat[LIBMESH_DIM][LIBMESH_DIM], work[10];
994  PetscReal w[LIBMESH_DIM];
995  PetscBLASInt nd = LIBMESH_DIM, lwork = 10, info;
996 
997  c = a.transpose() * a;
998 
999  for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
1000  for (unsigned int j = 0; j < LIBMESH_DIM; ++j)
1001  cmat[i][j] = c(i, j);
1002 
1003  LAPACKsyev_("V", "U", &nd, &cmat[0][0], &nd, w, work, &lwork, &info);
1004 
1005  if (info != 0)
1006  mooseError("FiniteStrainCrystalPLasticity: DSYEV function call in getMatRot function failed");
1007 
1008  diag.zero();
1009 
1010  for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
1011  diag(i, i) = std::sqrt(w[i]);
1012 
1013  for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
1014  for (unsigned int j = 0; j < LIBMESH_DIM; ++j)
1015  evec(i, j) = cmat[i][j];
1016 
1017  rot = a * ((evec.transpose() * diag * evec).inverse());
1018 
1019  return rot;
1020 }
1021 
1022 // Calculates tangent moduli which is used for global solve
1023 void
1025 {
1026 }
1027 
1028 RankFourTensor
1030 {
1031  RankFourTensor tan_mod;
1032 
1033  switch (_tan_mod_type)
1034  {
1035  case 0:
1036  tan_mod = elastoPlasticTangentModuli();
1037  break;
1038  default:
1039  tan_mod = elasticTangentModuli();
1040  }
1041 
1042  return tan_mod;
1043 }
1044 
1045 void
1047 {
1048  DenseVector<Real> mo(LIBMESH_DIM * _nss), no(LIBMESH_DIM * _nss);
1049 
1050  // Update slip direction and normal with crystal orientation
1051  for (unsigned int i = 0; i < _nss; ++i)
1052  {
1053  for (unsigned int j = 0; j < LIBMESH_DIM; ++j)
1054  {
1055  mo(i * LIBMESH_DIM + j) = 0.0;
1056  for (unsigned int k = 0; k < LIBMESH_DIM; ++k)
1057  mo(i * LIBMESH_DIM + j) =
1058  mo(i * LIBMESH_DIM + j) + _crysrot[_qp](j, k) * _mo(i * LIBMESH_DIM + k);
1059  }
1060 
1061  for (unsigned int j = 0; j < LIBMESH_DIM; ++j)
1062  {
1063  no(i * LIBMESH_DIM + j) = 0.0;
1064  for (unsigned int k = 0; k < LIBMESH_DIM; ++k)
1065  no(i * LIBMESH_DIM + j) =
1066  no(i * LIBMESH_DIM + j) + _crysrot[_qp](j, k) * _no(i * LIBMESH_DIM + k);
1067  }
1068  }
1069 
1070  // Calculate Schmid tensor and resolved shear stresses
1071  for (unsigned int i = 0; i < _nss; ++i)
1072  for (unsigned int j = 0; j < LIBMESH_DIM; ++j)
1073  for (unsigned int k = 0; k < LIBMESH_DIM; ++k)
1074  _s0[i](j, k) = mo(i * LIBMESH_DIM + j) * no(i * LIBMESH_DIM + k);
1075 }
1076 
1077 RankFourTensor
1079 {
1080  RankFourTensor tan_mod;
1081  RankTwoTensor pk2fet, fepk2;
1082  RankFourTensor deedfe, dsigdpk2dfe;
1083 
1084  // Fill in the matrix stiffness material property
1085 
1086  for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
1087  for (unsigned int j = 0; j < LIBMESH_DIM; ++j)
1088  for (unsigned int k = 0; k < LIBMESH_DIM; ++k)
1089  {
1090  deedfe(i, j, k, i) = deedfe(i, j, k, i) + _fe(k, j) * 0.5;
1091  deedfe(i, j, k, j) = deedfe(i, j, k, j) + _fe(k, i) * 0.5;
1092  }
1093 
1094  dsigdpk2dfe = _fe.mixedProductIkJl(_fe) * _elasticity_tensor[_qp] * deedfe;
1095 
1096  pk2fet = _pk2_tmp * _fe.transpose();
1097  fepk2 = _fe * _pk2_tmp;
1098 
1099  for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
1100  for (unsigned int j = 0; j < LIBMESH_DIM; ++j)
1101  for (unsigned int l = 0; l < LIBMESH_DIM; ++l)
1102  {
1103  tan_mod(i, j, i, l) = tan_mod(i, j, i, l) + pk2fet(l, j);
1104  tan_mod(i, j, j, l) = tan_mod(i, j, j, l) + fepk2(i, l);
1105  }
1106 
1107  tan_mod += dsigdpk2dfe;
1108 
1109  Real je = _fe.det();
1110  if (je > 0.0)
1111  tan_mod /= je;
1112 
1113  return tan_mod;
1114 }
1115 
1116 RankFourTensor
1118 {
1119  return _elasticity_tensor[_qp]; // update jacobian_mult
1120 }
1121 
1122 bool
1123 FiniteStrainCrystalPlasticity::line_search_update(const Real rnorm_prev, const RankTwoTensor dpk2)
1124 {
1125  if (_lsrch_method == "CUT_HALF")
1126  {
1127  Real rnorm;
1128  RankTwoTensor resid;
1129  Real step = 1.0;
1130 
1131  do
1132  {
1133  _pk2_tmp = _pk2_tmp - step * dpk2;
1134  step /= 2.0;
1135  _pk2_tmp = _pk2_tmp + step * dpk2;
1136 
1137  calcResidual(resid);
1138  rnorm = resid.L2norm();
1139  } while (rnorm > rnorm_prev && step > _min_lsrch_step);
1140 
1141  if (rnorm > rnorm_prev && step <= _min_lsrch_step)
1142  return false;
1143 
1144  return true;
1145  }
1146  else if (_lsrch_method == "BISECTION")
1147  {
1148  unsigned int count = 0;
1149  Real step_a = 0.0;
1150  Real step_b = 1.0;
1151  Real step = 1.0;
1152  Real s_m = 1000.0;
1153  Real rnorm = 1000.0;
1154 
1155  RankTwoTensor resid;
1156  calcResidual(resid);
1157  Real s_b = resid.doubleContraction(dpk2);
1158  Real rnorm1 = resid.L2norm();
1159  _pk2_tmp = _pk2_tmp - dpk2;
1160  calcResidual(resid);
1161  Real s_a = resid.doubleContraction(dpk2);
1162  Real rnorm0 = resid.L2norm();
1163  _pk2_tmp = _pk2_tmp + dpk2;
1164 
1165  if ((rnorm1 / rnorm0) < _lsrch_tol || s_a * s_b > 0)
1166  {
1167  calcResidual(resid);
1168  return true;
1169  }
1170 
1171  while ((rnorm / rnorm0) > _lsrch_tol && count < _lsrch_max_iter)
1172  {
1173  _pk2_tmp = _pk2_tmp - step * dpk2;
1174  step = 0.5 * (step_b + step_a);
1175  _pk2_tmp = _pk2_tmp + step * dpk2;
1176  calcResidual(resid);
1177  s_m = resid.doubleContraction(dpk2);
1178  rnorm = resid.L2norm();
1179 
1180  if (s_m * s_a < 0.0)
1181  {
1182  step_b = step;
1183  s_b = s_m;
1184  }
1185  if (s_m * s_b < 0.0)
1186  {
1187  step_a = step;
1188  s_a = s_m;
1189  }
1190  count++;
1191  }
1192 
1193  if ((rnorm / rnorm0) < _lsrch_tol && count < _lsrch_max_iter)
1194  return true;
1195 
1196  return false;
1197  }
1198  else
1199  {
1200  mooseError("Line search meothod is not provided.");
1201  return false;
1202  }
1203 }
1204 
1205 void
1207 {
1208  _fp_prev_inv = _fp_inv; // update _fp_prev_inv
1209 }
const MaterialProperty< RankTwoTensor > & _pk2_old
MaterialProperty< RankFourTensor > & _Jacobian_mult
derivative of stress w.r.t. strain (_dstress_dstrain)
void internalVariableUpdateNRiteration()
This function updates internal variables after each NewTon Raphson iteration (_fp_inv) ...
MaterialProperty< RankTwoTensor > & _lag_e
InputParameters validParams< FiniteStrainCrystalPlasticity >()
const MaterialProperty< RankTwoTensor > & _deformation_gradient
ComputeStressBase is the base class for stress tensors.
bool _input_rndm_scale_var
Input option for scaling variable to generate random stress when convergence fails.
virtual void getHardnessParams()
This function assign flow rate parameters from .i file - see test.
bool _use_line_search
Flag to activate line serach.
const MaterialProperty< std::vector< Real > > & _gss_old
const MaterialProperty< RankFourTensor > & _elasticity_tensor
virtual void calcJacobian(RankFourTensor &)
This function calculate jacobian.
virtual void computeQpStress()
This function updates the stress at a quadrature point.
MaterialProperty< RankTwoTensor > & _pk2
virtual void initQpStatefulProperties()
This function initializes the stateful properties such as stress, plastic deformation gradient...
virtual void solveStress()
This function solves for stress, updates plastic deformation gradient.
std::string _slip_sys_flow_prop_file_name
File should contain values of the flow rate equation parameters.
virtual void solveStatevar()
This function solves internal variables.
FiniteStrainCrystalPlasticity(const InputParameters &parameters)
MaterialProperty< RankTwoTensor > & _stress
virtual void update_slip_system_resistance()
This function updates the slip system resistances.
unsigned int _num_slip_sys_flowrate_props
Number of slip system flow rate parameters.
virtual void calcResidual(RankTwoTensor &)
This function calculate stress residual.
const MaterialProperty< RankTwoTensor > & _fp_old
RankTwoTensor getMatRot(const RankTwoTensor &a)
This function perform RU decomposition to obtain the rotation tensor.
virtual void readFileHardnessParams()
This function read hardness parameters from file.
virtual void updateGss()
This function updates the slip system resistances.
Real _rtol
Stress residual equation relative tolerance.
std::string _slip_sys_hard_prop_file_name
The hardening parameters in this class are read from .i file. The user can override to read from file...
unsigned int _maxiter
Maximum number of iterations for stress update.
virtual void assignSlipSysRes()
This function assign initial values of slip system resistances/internal variables read from getSlipSy...
Real _dfgrd_scale_factor
Scales the substepping increment to obtain deformation gradient at a substep iteration.
const unsigned int _nss
Number of slip system resistance.
virtual RankFourTensor elastoPlasticTangentModuli()
This function calculate the exact tangent moduli for preconditioner.
RankTwoTensor _delta_dfgrd
Flag to check whether convergence is achieved.
virtual RankFourTensor elasticTangentModuli()
This function calculate the elastic tangent moduli for preconditioner.
MaterialProperty< std::vector< Real > > & _gss
virtual RankFourTensor calcTangentModuli()
This function calculate the tangent moduli for preconditioner.
bool line_search_update(const Real rnorm_prev, const RankTwoTensor)
This function performs the line search update.
const MaterialProperty< RankTwoTensor > & _deformation_gradient_old
const MaterialProperty< RankTwoTensor > & _crysrot
void calc_schmid_tensor()
This function calculate the Schmid tensor.
MaterialProperty< RankTwoTensor > & _fp
virtual void postSolveStatevar()
This function update internal variable after solve.
unsigned int _lsrch_max_iter
Line search bisection method maximum iteration number.
virtual void getFlowRateParams()
This function assign flow rate parameters - see test.
Real _gtol
Internal variable update equation tolerance.
std::string _slip_sys_res_prop_file_name
File should contain initial values of the slip system resistances.
virtual void postSolveQp()
This function update stress and internal variable after solve.
virtual void preSolveStress()
This function set variables for stress solve.
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)
InputParameters validParams< ComputeStressBase >()
virtual void postSolveStress()
This function update stress and plastic deformation gradient after solve.
virtual void calc_resid_jacob(RankTwoTensor &, RankFourTensor &)
This function calls the residual and jacobian functions used in the stress update algorithm...
MooseEnum _tan_mod_type
Type of tangent moduli calculation.
Real _lsrch_tol
Line search bisection method tolerance.
virtual void computeQpElasticityTensor()
This function updates the elasticity tensor at a quadrature point.
virtual void getSlipSystems()
This function reads slip system from file - see test.
unsigned int _num_slip_sys_props
Number of slip system specific properties provided in the file containing slip system normals and dir...
virtual void readFileFlowRateParams()
This function read flow rate parameters from file - see test.
unsigned int _max_substep_iter
Maximum number of substep iterations.
bool _first_step_iter
Flags to reset variables and reinitialize variables.
virtual void initSlipSysProps()
This function initializes slip system resistances.
virtual void getInitSlipSysRes()
This function assign slip system resistances - see test.
virtual void initAdditionalProps()
This function initializes additional parameters.
virtual void preSolveStatevar()
This function set variables for internal variable solve.
MaterialProperty< RankTwoTensor > & _update_rot
Real _abs_tol
Stress residual equation absolute tolerance.
virtual void readFileInitSlipSysRes()
This function read slip system resistances from file - see test.
unsigned int _maxiterg
Maximum number of iterations for internal variable update.
virtual void getSlipIncrements()
This function updates the slip increments.
const MaterialProperty< Real > & _acc_slip_old
RankTwoTensor get_current_rotation(const RankTwoTensor &a)
This function perform RU decomposition to obtain the rotation tensor.
virtual void solveQp()
This function solves stress and internal variables.
MooseEnum _intvar_read_type
Read from options for initial values of internal variables.
Real _min_lsrch_step
Minimum line search step size.
virtual void preSolveQp()
This function set variables for stress and internal variable solve.
std::string _slip_sys_file_name
File should contain slip plane normal and direction. See test.
Real _slip_incr_tol
Slip increment tolerance.