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