www.mooseframework.org
SolidModel.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 /****************************************************************/
7 
8 #include "SolidModel.h"
9 #include "AxisymmetricRZ.h"
10 #include "NonlinearRZ.h"
11 #include "SphericalR.h"
12 #include "Linear.h"
13 #include "Nonlinear3D.h"
14 #include "PlaneStrain.h"
15 #include "NonlinearPlaneStrain.h"
16 #include "VolumetricModel.h"
17 #include "ConstitutiveModel.h"
19 #include "MooseApp.h"
20 #include "Problem.h"
21 #include "PiecewiseLinear.h"
22 
23 #include "libmesh/quadrature.h"
24 
25 template <>
26 InputParameters
28 {
29  MooseEnum formulation(
30  "Nonlinear3D NonlinearRZ AxisymmetricRZ SphericalR Linear PlaneStrain NonlinearPlaneStrain");
31  MooseEnum compute_method("NoShearRetention ShearRetention");
32 
33  InputParameters params = validParams<Material>();
34  params.addParam<std::string>(
35  "appended_property_name", "", "Name appended to material properties to make them unique");
36  params.addParam<Real>("bulk_modulus", "The bulk modulus for the material.");
37  params.addParam<Real>("lambda", "Lame's first parameter for the material.");
38  params.addParam<Real>("poissons_ratio", "Poisson's ratio for the material.");
39  params.addParam<FunctionName>("poissons_ratio_function",
40  "Poisson's ratio as a function of temperature.");
41  params.addParam<Real>("shear_modulus", "The shear modulus of the material.");
42  params.addParam<Real>("youngs_modulus", "Young's modulus of the material.");
43  params.addParam<FunctionName>("youngs_modulus_function",
44  "Young's modulus as a function of temperature.");
45  params.addParam<Real>("thermal_expansion", "The thermal expansion coefficient.");
46  params.addParam<FunctionName>("thermal_expansion_function",
47  "Thermal expansion coefficient as a function of temperature.");
48  params.addCoupledVar("temp", "Coupled Temperature");
49  params.addParam<Real>(
50  "stress_free_temperature",
51  "The stress-free temperature. If not specified, the initial temperature is used.");
52  params.addParam<Real>("thermal_expansion_reference_temperature",
53  "Reference temperature for mean thermal expansion function.");
54  MooseEnum cte_function_type("instantaneous mean");
55  params.addParam<MooseEnum>("thermal_expansion_function_type",
56  cte_function_type,
57  "Type of thermal expansion function. Choices are: " +
58  cte_function_type.getRawNames());
59  params.addParam<std::vector<Real>>("initial_stress",
60  "The initial stress tensor (xx, yy, zz, xy, yz, zx)");
61  params.addParam<std::string>(
62  "cracking_release",
63  "abrupt",
64  "The cracking release type. Choices are abrupt (default) and exponential.");
65  params.addParam<Real>("cracking_stress",
66  0.0,
67  "The stress threshold beyond which cracking occurs. Must be positive.");
68  params.addParam<Real>(
69  "cracking_residual_stress",
70  0.0,
71  "The fraction of the cracking stress allowed to be maintained following a crack.");
72  params.addParam<Real>("cracking_beta", 1.0, "The coefficient used in the exponetional model.");
73  params.addParam<MooseEnum>(
74  "compute_method", compute_method, "The method used in the stress calculation.");
75  params.addParam<FunctionName>(
76  "cracking_stress_function", "", "The cracking stress as a function of time and location");
77  params.addParam<std::vector<unsigned int>>(
78  "active_crack_planes", "Planes on which cracks are allowed (0,1,2 -> x,z,theta in RZ)");
79  params.addParam<unsigned int>(
80  "max_cracks", 3, "The maximum number of cracks allowed at a material point.");
81  params.addParam<Real>("cracking_neg_fraction",
82  "The fraction of the cracking strain at which a "
83  "transitition begins during decreasing strain to "
84  "the original stiffness.");
85  params.addParam<MooseEnum>("formulation",
86  formulation,
87  "Element formulation. Choices are: " + formulation.getRawNames());
88  params.addParam<std::string>("increment_calculation",
89  "RashidApprox",
90  "The algorithm to use when computing the "
91  "incremental strain and rotation (RashidApprox or "
92  "Eigen). For use with Nonlinear3D/RZ formulation.");
93  params.addParam<bool>("large_strain",
94  false,
95  "Whether to include large strain terms in "
96  "AxisymmetricRZ, SphericalR, and PlaneStrain "
97  "formulations.");
98  params.addParam<bool>("compute_JIntegral", false, "Whether to compute the J Integral.");
99  params.addParam<bool>(
100  "compute_InteractionIntegral", false, "Whether to compute the Interaction Integral.");
101  params.addParam<bool>("store_stress_older",
102  false,
103  "Parameter which indicates whether the older "
104  "stress state, required for HHT time "
105  "integration, needs to be stored");
106  params.addCoupledVar("disp_r", "The r displacement");
107  params.addCoupledVar("disp_x", "The x displacement");
108  params.addCoupledVar("disp_y", "The y displacement");
109  params.addCoupledVar("disp_z", "The z displacement");
110  params.addCoupledVar("strain_zz", "The zz strain");
111  params.addCoupledVar("scalar_strain_zz", "The zz strain (scalar variable)");
112  params.addParam<std::vector<std::string>>(
113  "dep_matl_props", "Names of material properties this material depends on.");
114  params.addParam<std::string>("constitutive_model", "ConstitutiveModel to use (optional)");
115  params.addParam<bool>("volumetric_locking_correction",
116  true,
117  "Set to false to turn off volumetric locking correction");
118  return params;
119 }
120 
121 namespace
122 {
124 getCrackingModel(const std::string & name)
125 {
126  std::string n(name);
127  std::transform(n.begin(), n.end(), n.begin(), ::tolower);
129  if (n == "abrupt")
131  else if (n == "exponential")
133  else if (n == "power")
135  if (cm == SolidModel::CR_UNKNOWN)
136  mooseError("Unknown cracking model");
137  return cm;
138 }
139 }
140 
141 SolidModel::SolidModel(const InputParameters & parameters)
142  : DerivativeMaterialInterface<Material>(parameters),
143  _appended_property_name(getParam<std::string>("appended_property_name")),
144  _bulk_modulus_set(parameters.isParamValid("bulk_modulus")),
145  _lambda_set(parameters.isParamValid("lambda")),
146  _poissons_ratio_set(parameters.isParamValid("poissons_ratio")),
147  _shear_modulus_set(parameters.isParamValid("shear_modulus")),
148  _youngs_modulus_set(parameters.isParamValid("youngs_modulus")),
149  _bulk_modulus(_bulk_modulus_set ? getParam<Real>("bulk_modulus") : -1),
150  _lambda(_lambda_set ? getParam<Real>("lambda") : -1),
151  _poissons_ratio(_poissons_ratio_set ? getParam<Real>("poissons_ratio") : -1),
152  _shear_modulus(_shear_modulus_set ? getParam<Real>("shear_modulus") : -1),
153  _youngs_modulus(_youngs_modulus_set ? getParam<Real>("youngs_modulus") : -1),
154  _youngs_modulus_function(
155  isParamValid("youngs_modulus_function") ? &getFunction("youngs_modulus_function") : NULL),
156  _poissons_ratio_function(
157  isParamValid("poissons_ratio_function") ? &getFunction("poissons_ratio_function") : NULL),
158  _cracking_release(getCrackingModel(getParam<std::string>("cracking_release"))),
159  _cracking_stress(
160  parameters.isParamValid("cracking_stress")
161  ? (getParam<Real>("cracking_stress") > 0 ? getParam<Real>("cracking_stress") : -1)
162  : -1),
163  _cracking_residual_stress(getParam<Real>("cracking_residual_stress")),
164  _cracking_beta(getParam<Real>("cracking_beta")),
165  _compute_method(getParam<MooseEnum>("compute_method")),
166  _cracking_stress_function(getParam<FunctionName>("cracking_stress_function") != ""
167  ? &getFunction("cracking_stress_function")
168  : NULL),
169  _cracking_alpha(0),
170  _active_crack_planes(3, 1),
171  _max_cracks(getParam<unsigned int>("max_cracks")),
172  _cracking_neg_fraction(
173  isParamValid("cracking_neg_fraction") ? getParam<Real>("cracking_neg_fraction") : 0),
174  _has_temp(isCoupled("temp")),
175  _temperature(_has_temp ? coupledValue("temp") : _zero),
176  _temperature_old(_has_temp ? coupledValueOld("temp") : _zero),
177  _temp_grad(coupledGradient("temp")),
178  _alpha(parameters.isParamValid("thermal_expansion") ? getParam<Real>("thermal_expansion") : 0.),
179  _alpha_function(parameters.isParamValid("thermal_expansion_function")
180  ? &getFunction("thermal_expansion_function")
181  : NULL),
182  _piecewise_linear_alpha_function(NULL),
183  _has_stress_free_temp(false),
184  _stress_free_temp(0.0),
185  _ref_temp(0.0),
186  _volumetric_models(),
187  _dep_matl_props(),
188  _stress(createProperty<SymmTensor>("stress")),
189  _stress_old_prop(getPropertyOld<SymmTensor>("stress")),
190  _stress_old(0),
191  _total_strain(createProperty<SymmTensor>("total_strain")),
192  _total_strain_old(getPropertyOld<SymmTensor>("total_strain")),
193  _elastic_strain(createProperty<SymmTensor>("elastic_strain")),
194  _elastic_strain_old(getPropertyOld<SymmTensor>("elastic_strain")),
195  _crack_flags(NULL),
196  _crack_flags_old(NULL),
197  _crack_flags_local(),
198  _crack_count(NULL),
199  _crack_count_old(NULL),
200  _crack_rotation(NULL),
201  _crack_rotation_old(NULL),
202  _crack_strain(NULL),
203  _crack_strain_old(NULL),
204  _crack_max_strain(NULL),
205  _crack_max_strain_old(NULL),
206  _principal_strain(3, 1),
207  _elasticity_tensor(createProperty<SymmElasticityTensor>("elasticity_tensor")),
208  _Jacobian_mult(createProperty<SymmElasticityTensor>("Jacobian_mult")),
209  _d_strain_dT(),
210  _d_stress_dT(createProperty<SymmTensor>("d_stress_dT")),
211  _total_strain_increment(0),
212  _strain_increment(0),
213  _compute_JIntegral(getParam<bool>("compute_JIntegral")),
214  _compute_InteractionIntegral(getParam<bool>("compute_InteractionIntegral")),
215  _store_stress_older(getParam<bool>("store_stress_older")),
216  _SED(NULL),
217  _SED_old(NULL),
218  _Eshelby_tensor(NULL),
219  _J_thermal_term_vec(NULL),
220  _current_instantaneous_thermal_expansion_coef(NULL),
221  _block_id(std::vector<SubdomainID>(blockIDs().begin(), blockIDs().end())),
222  _constitutive_active(false),
223  _step_zero(declareRestartableData<bool>("step_zero", true)),
224  _step_one(declareRestartableData<bool>("step_one", true)),
225  _element(NULL),
226  _local_elasticity_tensor(NULL)
227 {
228  bool same_coord_type = true;
229 
230  for (unsigned int i = 1; i < _block_id.size(); ++i)
231  same_coord_type &=
232  (_subproblem.getCoordSystem(_block_id[0]) == _subproblem.getCoordSystem(_block_id[i]));
233  if (!same_coord_type)
234  mooseError("Material '",
235  name(),
236  "' was specified on multiple blocks that do not have the same coordinate system");
237  // Use the first block to figure out the coordinate system (the above check ensures that they are
238  // the same)
239  _coord_type = _subproblem.getCoordSystem(_block_id[0]);
241 
242  const std::vector<std::string> & dmp = getParam<std::vector<std::string>>("dep_matl_props");
243  _dep_matl_props.insert(dmp.begin(), dmp.end());
244  for (std::set<std::string>::const_iterator i = _dep_matl_props.begin();
245  i != _dep_matl_props.end();
246  ++i)
247  {
248  // Tell MOOSE that we need this MaterialProperty. This enables dependency checking.
249  getMaterialProperty<Real>(*i);
250  }
251 
253 
254  if (_cracking_stress > 0)
255  {
256  _crack_flags = &createProperty<RealVectorValue>("crack_flags");
257  _crack_flags_old = &getPropertyOld<RealVectorValue>("crack_flags");
259  {
260  _crack_count = &createProperty<RealVectorValue>("crack_count");
261  _crack_count_old = &getPropertyOld<RealVectorValue>("crack_count");
262  }
263  _crack_rotation = &createProperty<ColumnMajorMatrix>("crack_rotation");
264  _crack_rotation_old = &getPropertyOld<ColumnMajorMatrix>("crack_rotation");
265  _crack_max_strain = &createProperty<RealVectorValue>("crack_max_strain");
266  _crack_max_strain_old = &getPropertyOld<RealVectorValue>("crack_max_strain");
267  _crack_strain = &createProperty<RealVectorValue>("crack_strain");
268  _crack_strain_old = &getPropertyOld<RealVectorValue>("crack_strain");
269 
270  if (parameters.isParamValid("active_crack_planes"))
271  {
272  const std::vector<unsigned int> & planes =
273  getParam<std::vector<unsigned>>("active_crack_planes");
274  for (unsigned i(0); i < 3; ++i)
275  _active_crack_planes[i] = 0;
276 
277  for (unsigned i(0); i < planes.size(); ++i)
278  {
279  if (planes[i] > 2)
280  mooseError("Active planes must be 0, 1, or 2");
281  _active_crack_planes[planes[i]] = 1;
282  }
283  }
284  if (_cracking_residual_stress < 0 || _cracking_residual_stress > 1)
285  {
286  mooseError("cracking_residual_stress must be between 0 and 1");
287  }
288  if (isParamValid("cracking_neg_fraction") &&
289  (_cracking_neg_fraction <= 0 || _cracking_neg_fraction > 1))
290  {
291  mooseError("cracking_neg_fraction must be > zero and <= 1");
292  }
293  }
294 
295  if (parameters.isParamValid("stress_free_temperature"))
296  {
297  _has_stress_free_temp = true;
298  _stress_free_temp = getParam<Real>("stress_free_temperature");
299  if (!_has_temp)
300  mooseError("Cannot specify stress_free_temperature without coupling to temperature");
301  }
302 
303  if (parameters.isParamValid("thermal_expansion_function_type"))
304  {
305  if (!_alpha_function)
306  mooseError("thermal_expansion_function_type can only be set when thermal_expansion_function "
307  "is used");
308  MooseEnum tec = getParam<MooseEnum>("thermal_expansion_function_type");
309  if (tec == "mean")
310  _mean_alpha_function = true;
311  else if (tec == "instantaneous")
312  _mean_alpha_function = false;
313  else
314  mooseError("Invalid option for thermal_expansion_function_type");
315  }
316  else
317  _mean_alpha_function = false;
318 
319  if (parameters.isParamValid("thermal_expansion_reference_temperature"))
320  {
321  if (!_alpha_function)
322  mooseError("thermal_expansion_reference_temperature can only be set when "
323  "thermal_expansion_function is used");
325  mooseError("thermal_expansion_reference_temperature can only be set when "
326  "thermal_expansion_function_type = mean");
327  _ref_temp = getParam<Real>("thermal_expansion_reference_temperature");
328  if (!_has_temp)
329  mooseError(
330  "Cannot specify thermal_expansion_reference_temperature without coupling to temperature");
331  }
332 
334  {
335  if (!parameters.isParamValid("thermal_expansion_reference_temperature") ||
337  mooseError(
338  "Must specify both stress_free_temperature and thermal_expansion_reference_temperature "
339  "if thermal_expansion_function_type = mean");
340  }
341 
342  if (parameters.isParamValid("thermal_expansion") &&
343  parameters.isParamValid("thermal_expansion_function"))
344  mooseError("Cannot specify both thermal_expansion and thermal_expansion_function");
345 
346  if (_compute_JIntegral)
347  {
348  _SED = &declareProperty<Real>("strain_energy_density");
349  _SED_old = &getMaterialPropertyOld<Real>("strain_energy_density");
350  _Eshelby_tensor = &declareProperty<RankTwoTensor>("Eshelby_tensor");
351  _J_thermal_term_vec = &declareProperty<RealVectorValue>("J_thermal_term_vec");
353  &declareProperty<Real>("current_instantaneous_thermal_expansion_coef");
354  }
355 
357  !hasMaterialProperty<Real>("current_instantaneous_thermal_expansion_coef"))
359  &declareProperty<Real>("current_instantaneous_thermal_expansion_coef");
360 }
361 
363 
365 {
367  delete _element;
368 }
369 
371 
372 void
374 {
375  int num_elastic_constants = _bulk_modulus_set + _lambda_set + _poissons_ratio_set +
377 
378  if (num_elastic_constants != 2)
379  {
380  std::string err("Exactly two elastic constants must be defined for material '");
381  err += name();
382  err += "'.";
383  mooseError(err);
384  }
385 
386  if (_bulk_modulus_set && _bulk_modulus <= 0)
387  {
388  std::string err("Bulk modulus must be positive in material '");
389  err += name();
390  err += "'.";
391  mooseError(err);
392  }
393  if (_poissons_ratio_set && (_poissons_ratio <= -1.0 || _poissons_ratio >= 0.5))
394  {
395  std::string err("Poissons ratio must be greater than -1 and less than 0.5 in material '");
396  err += name();
397  err += "'.";
398  mooseError(err);
399  }
401  {
402  std::string err("Shear modulus must not be negative in material '");
403  err += name();
404  err += "'.";
405  mooseError(err);
406  }
407  if (_youngs_modulus_set && _youngs_modulus <= 0)
408  {
409  std::string err("Youngs modulus must be positive in material '");
410  err += name();
411  err += "'.";
412  mooseError(err);
413  }
414 
415  // Calculate lambda, the shear modulus, and Young's modulus
416  if (_lambda_set && _shear_modulus_set) // First and second Lame
417  {
421  }
422  else if (_lambda_set && _poissons_ratio_set)
423  {
424  _shear_modulus = (_lambda * (1.0 - 2.0 * _poissons_ratio)) / (2.0 * _poissons_ratio);
427  }
428  else if (_lambda_set && _bulk_modulus_set)
429  {
430  _shear_modulus = 3.0 * (_bulk_modulus - _lambda) / 2.0;
434  }
435  else if (_lambda_set && _youngs_modulus_set)
436  {
438  ((_youngs_modulus - 3.0 * _lambda) / 4.0) +
439  (std::sqrt((_youngs_modulus - 3.0 * _lambda) * (_youngs_modulus - 3.0 * _lambda) +
440  8.0 * _lambda * _youngs_modulus) /
441  4.0);
443  }
445  {
446  _lambda = (2.0 * _shear_modulus * _poissons_ratio) / (1.0 - 2.0 * _poissons_ratio);
449  }
451  {
452  _lambda = _bulk_modulus - 2.0 * _shear_modulus / 3.0;
456  (3 * _bulk_modulus - 2 * _shear_modulus) / (2 * (3 * _bulk_modulus + _shear_modulus));
457  }
458  else if (_shear_modulus_set && _youngs_modulus_set)
459  {
463  }
465  {
468  (3.0 * _bulk_modulus * (1.0 - 2.0 * _poissons_ratio)) / (2.0 * (1.0 + _poissons_ratio));
471  }
472  else if (_youngs_modulus_set && _poissons_ratio_set) // Young's Modulus and Poisson's Ratio
473  {
475  ((1.0 + _poissons_ratio) * (1 - 2.0 * _poissons_ratio));
476  _shear_modulus = _youngs_modulus / (2.0 * (1.0 + _poissons_ratio));
477  }
478  else if (_youngs_modulus_set && _bulk_modulus_set)
479  {
480  _lambda = 3.0 * _bulk_modulus * (3.0 * _bulk_modulus - _youngs_modulus) /
481  (9.0 * _bulk_modulus - _youngs_modulus);
485  }
486 
487  _lambda_set = true;
488  _shear_modulus_set = true;
489  _youngs_modulus_set = true;
490  _poissons_ratio_set = true;
491 }
492 
494 
495 void
497 {
498  bool constant(true);
499 
502  {
503  constant = false;
504  }
505 
507  mooseAssert(_youngs_modulus_set, "Internal error: Youngs modulus not set");
508  mooseAssert(_poissons_ratio_set, "Internal error: Poissons ratio not set");
511  iso->calculate(0);
512  elasticityTensor(iso);
513 }
514 
516 
517 void
519 {
520  // if (_cracking_stress > 0)
521  // {
522  // _cracked_this_step_count.clear();
523  // }
524 }
525 
527 
528 void
530 {
531  // if (_cracking_stress > 0)
532  // {
533  // for (std::map<Point, unsigned>::iterator i = _cracked_this_step.begin();
534  // i != _cracked_this_step.end(); ++i)
535  // {
536  // if (i->second)
537  // {
538  // ++_cracked_this_step_count[i->first];
539  // }
540  // }
541  // }
542 }
543 
545 
546 void
548 {
551 }
552 
554 
555 void
557 {
558  bool modified = false;
559  _d_strain_dT.zero();
560 
561  const SubdomainID current_block = _current_elem->subdomain_id();
563  {
564  MooseSharedPointer<ConstitutiveModel> cm = _constitutive_model[current_block];
565 
566  // Let's be a little careful and check for a non-existent
567  // ConstitutiveModel, which could be returned as a default value
568  // from std::map::operator[]
569  if (!cm)
570  mooseError("ConstitutiveModel not available for block ", current_block);
571 
572  cm->setQp(_qp);
573  modified |= cm->modifyStrainIncrement(*_current_elem, _strain_increment, _d_strain_dT);
574  }
575 
576  if (!modified)
577  {
579  }
580 
582 }
583 
585 
586 void
588 {
589  if (_has_temp && !_step_zero)
590  {
591  Real inc_thermal_strain;
592  Real d_thermal_strain_d_temp;
593 
594  Real old_temp;
596  old_temp = _stress_free_temp;
597  else
598  old_temp = _temperature_old[_qp];
599 
600  Real current_temp = _temperature[_qp];
601 
602  Real delta_t = current_temp - old_temp;
603 
604  Real alpha = _alpha;
605 
606  if (_alpha_function)
607  {
608  Point p;
609  Real alpha_current_temp = _alpha_function->value(current_temp, p);
610  Real alpha_old_temp = _alpha_function->value(old_temp, p);
611 
613  {
614  Real alpha_stress_free_temperature = _alpha_function->value(_stress_free_temp, p);
615  Real small(1e-6);
616 
617  Real numerator = alpha_current_temp * (current_temp - _ref_temp) -
618  alpha_old_temp * (old_temp - _ref_temp);
619  Real denominator = 1.0 + alpha_stress_free_temperature * (_stress_free_temp - _ref_temp);
620  if (denominator < small)
621  mooseError("Denominator too small in thermal strain calculation");
622  inc_thermal_strain = numerator / denominator;
623  d_thermal_strain_d_temp = alpha_current_temp * (current_temp - _ref_temp);
624  }
625  else
626  {
627  inc_thermal_strain = delta_t * 0.5 * (alpha_current_temp + alpha_old_temp);
628  d_thermal_strain_d_temp = alpha_current_temp;
629  }
630  }
631  else
632  {
633  inc_thermal_strain = delta_t * alpha;
634  d_thermal_strain_d_temp = alpha;
635  }
636 
637  _strain_increment.addDiag(-inc_thermal_strain);
638  _d_strain_dT.addDiag(-d_thermal_strain_d_temp);
639  }
640 }
641 
643 
644 void
646 {
647  const Real V0Vold = 1 / _element->volumeRatioOld(_qp);
648  const SubdomainID current_block = _current_elem->subdomain_id();
649  const std::vector<MooseSharedPointer<VolumetricModel>> & vm(_volumetric_models[current_block]);
650  for (unsigned int i(0); i < vm.size(); ++i)
651  {
652  vm[i]->modifyStrain(_qp, V0Vold, _strain_increment, _d_strain_dT);
653  }
654 }
655 
657 
658 void
659 SolidModel::rotateSymmetricTensor(const ColumnMajorMatrix & R,
660  const SymmTensor & T,
661  SymmTensor & result)
662 {
663 
664  // R T Rt
665  // 00 01 02 00 01 02 00 10 20
666  // 10 11 12 * 10 11 12 * 01 11 21
667  // 20 21 22 20 21 22 02 12 22
668  //
669  const Real T00 = R(0, 0) * T.xx() + R(0, 1) * T.xy() + R(0, 2) * T.zx();
670  const Real T01 = R(0, 0) * T.xy() + R(0, 1) * T.yy() + R(0, 2) * T.yz();
671  const Real T02 = R(0, 0) * T.zx() + R(0, 1) * T.yz() + R(0, 2) * T.zz();
672 
673  const Real T10 = R(1, 0) * T.xx() + R(1, 1) * T.xy() + R(1, 2) * T.zx();
674  const Real T11 = R(1, 0) * T.xy() + R(1, 1) * T.yy() + R(1, 2) * T.yz();
675  const Real T12 = R(1, 0) * T.zx() + R(1, 1) * T.yz() + R(1, 2) * T.zz();
676 
677  const Real T20 = R(2, 0) * T.xx() + R(2, 1) * T.xy() + R(2, 2) * T.zx();
678  const Real T21 = R(2, 0) * T.xy() + R(2, 1) * T.yy() + R(2, 2) * T.yz();
679  const Real T22 = R(2, 0) * T.zx() + R(2, 1) * T.yz() + R(2, 2) * T.zz();
680 
681  result.xx(T00 * R(0, 0) + T01 * R(0, 1) + T02 * R(0, 2));
682  result.yy(T10 * R(1, 0) + T11 * R(1, 1) + T12 * R(1, 2));
683  result.zz(T20 * R(2, 0) + T21 * R(2, 1) + T22 * R(2, 2));
684  result.xy(T00 * R(1, 0) + T01 * R(1, 1) + T02 * R(1, 2));
685  result.yz(T10 * R(2, 0) + T11 * R(2, 1) + T12 * R(2, 2));
686  result.zx(T00 * R(2, 0) + T01 * R(2, 1) + T02 * R(2, 2));
687 }
688 
690 
691 void
693 {
694  if (isParamValid("initial_stress"))
695  {
696  const std::vector<Real> & s = getParam<std::vector<Real>>("initial_stress");
697  if (6 != s.size())
698  {
699  mooseError("initial_stress must give six values");
700  }
701  _stress[_qp].fillFromInputVector(s);
702  }
703 
704  if (_cracking_stress_function != NULL)
705  {
706  _cracking_stress = _cracking_stress_function->value(_t, _q_point[_qp]);
707  }
708  if (_cracking_stress > 0)
709  {
710  (*_crack_flags)[_qp](0) = (*_crack_flags)[_qp](1) = (*_crack_flags)[_qp](2) = 1;
711  if (_crack_count)
712  {
713  (*_crack_count)[_qp](0) = (*_crack_count)[_qp](1) = (*_crack_count)[_qp](2) = 0;
714  }
715 
716  (*_crack_rotation)[_qp].identity();
717  }
718  if (_SED)
719  (*_SED)[_qp] = 0;
720 }
721 
723 
724 void
726 {
727  if (_t_step >= 1)
728  _step_zero = false;
729 
730  if (_t_step >= 2)
731  _step_one = false;
732 
733  elementInit();
734  _element->init();
735 
736  for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
737  {
740 
742 
744 
746  computeStress();
747  else
749 
750  if (_compute_JIntegral)
752 
754 
756 
757  finalizeStress();
758 
759  if (_compute_JIntegral)
760  computeEshelby();
761 
764 
767 
769  }
770 }
771 
773 
774 void
776 {
777  mooseAssert(_SED, "_SED not initialized");
778  mooseAssert(_SED_old, "_SED_old not initialized");
779  (*_SED)[_qp] = (*_SED_old)[_qp] + _stress[_qp].doubleContraction(_strain_increment) / 2 +
780  _stress_old_prop[_qp].doubleContraction(_strain_increment) / 2;
781 }
782 
784 
785 void
787 {
788  mooseAssert(_SED, "_SED not initialized");
789  mooseAssert(_Eshelby_tensor, "_Eshelby_tensor not initialized");
790  // Cauchy stress (sigma) in a colum major matrix:
791  ColumnMajorMatrix stress_CMM;
792  stress_CMM(0, 0) = _stress[_qp].xx();
793  stress_CMM(0, 1) = _stress[_qp].xy();
794  stress_CMM(0, 2) = _stress[_qp].xz();
795  stress_CMM(1, 0) = _stress[_qp].xy();
796  stress_CMM(1, 1) = _stress[_qp].yy();
797  stress_CMM(1, 2) = _stress[_qp].yz();
798  stress_CMM(2, 0) = _stress[_qp].xz();
799  stress_CMM(2, 1) = _stress[_qp].yz();
800  stress_CMM(2, 2) = _stress[_qp].zz();
801 
802  // Deformation gradient (F):
803  ColumnMajorMatrix F;
805  // Displacement gradient (H):
806  ColumnMajorMatrix H(F);
807  H.addDiag(-1.0);
808  Real detF = _element->detMatrix(F);
809  ColumnMajorMatrix Finv;
810  _element->invertMatrix(F, Finv);
811  ColumnMajorMatrix FinvT;
812  FinvT = Finv.transpose();
813  ColumnMajorMatrix HT;
814  HT = H.transpose();
815 
816  // 1st Piola-Kirchoff Stress (P):
817  ColumnMajorMatrix piola;
818  piola = stress_CMM * FinvT;
819  piola *= detF;
820 
821  // HTP = H^T * P = H^T * detF * sigma * FinvT;
822  ColumnMajorMatrix HTP;
823  HTP = HT * piola;
824 
825  ColumnMajorMatrix WI;
826  WI.identity();
827  WI *= (*_SED)[_qp];
828  WI *= detF;
829  (*_Eshelby_tensor)[_qp] = WI - HTP;
830 }
831 
833 
834 void
836 {
837  // Given the stretching, compute the stress increment and add it to the old stress. Also update
838  // the creep strain
839  // stress = stressOld + stressIncrement
840 
841  if (_step_zero)
842  return;
843 
844  const SubdomainID current_block = _current_elem->subdomain_id();
845  MooseSharedPointer<ConstitutiveModel> cm = _constitutive_model[current_block];
846 
847  mooseAssert(_constitutive_active, "Logic error. ConstitutiveModel not active.");
848 
849  // Let's be a little careful and check for a non-existent
850  // ConstitutiveModel, which could be returned as a default value
851  // from std::map::operator[]
852  if (!cm)
853  mooseError("Logic error. No ConstitutiveModel for current_block=", current_block, ".");
854 
855  cm->setQp(_qp);
856  cm->computeStress(
857  *_current_elem, *elasticityTensor(), _stress_old, _strain_increment, _stress[_qp]);
858 }
859 
861 void
863 {
864  if (_cracking_stress_function != NULL)
865  {
866  _cracking_stress = _cracking_stress_function->value(_t, _q_point[_qp]);
867  }
868 
870 
872 
874 
876 
878 
879  if (changed || _cracking_stress > 0)
880  {
882  }
883 }
884 
886 
887 bool
889 {
890  bool changed(false);
892  {
893  const SubdomainID current_block = _current_elem->subdomain_id();
894  MooseSharedPointer<ConstitutiveModel> cm = _constitutive_model[current_block];
895 
896  // Let's be a little careful and check for a non-existent
897  // ConstitutiveModel, which could be returned as a default value
898  // from std::map::operator[]
899  if (!cm)
900  mooseError("ConstitutiveModel not available for block ", current_block);
901 
902  cm->setQp(_qp);
903  changed |= cm->updateElasticityTensor(tensor);
904  }
905 
907  {
908  SymmIsotropicElasticityTensor * t = dynamic_cast<SymmIsotropicElasticityTensor *>(&tensor);
909  if (!t)
910  {
911  mooseError("Cannot use Youngs modulus or Poissons ratio functions");
912  }
913  t->unsetConstants();
914  Point p;
916  ? _youngs_modulus_function->value(_temperature[_qp], p)
917  : _youngs_modulus));
919  ? _poissons_ratio_function->value(_temperature[_qp], p)
920  : _poissons_ratio));
921  changed = true;
922  }
923  return changed;
924 }
925 
927 
928 void
930 {
931  std::vector<SymmTensor *> t(3);
932  t[0] = &_elastic_strain[_qp];
933  t[1] = &_total_strain[_qp];
934  t[2] = &_stress[_qp];
936 }
937 
939 
940 void
942 {
943  mooseAssert(_local_elasticity_tensor, "null elasticity tensor");
944 
945  // _Jacobian_mult[_qp] = *_local_elasticity_tensor;
946  // _d_stress_dT[_qp] = *_local_elasticity_tensor * _d_strain_dT;
949 }
950 
952 
953 void
955 {
956 
958 
960 
961  // Load in the volumetric models and constitutive models
962  bool set_constitutive_active = false;
963  for (unsigned i(0); i < _block_id.size(); ++i)
964  {
965 
966  // const std::vector<Material*> * mats_p;
967  std::vector<MooseSharedPointer<Material>> const * mats_p;
968  std::string suffix;
969  if (_bnd)
970  {
971  mats_p = &_fe_problem.getMaterialWarehouse()[Moose::FACE_MATERIAL_DATA].getActiveBlockObjects(
972  _block_id[i], _tid);
973  suffix = "_face";
974  }
975  else
976  mats_p = &_fe_problem.getMaterialWarehouse().getActiveBlockObjects(_block_id[i], _tid);
977 
978  const std::vector<MooseSharedPointer<Material>> & mats = *mats_p;
979 
980  for (unsigned int j = 0; j < mats.size(); ++j)
981  {
982  MooseSharedPointer<VolumetricModel> vm =
983  MooseSharedNamespace::dynamic_pointer_cast<VolumetricModel>(mats[j]);
984  if (vm)
985  {
986  const std::vector<std::string> & dep_matl_props = vm->getDependentMaterialProperties();
987  for (unsigned k = 0; k < dep_matl_props.size(); ++k)
988  {
989  if ("" != dep_matl_props[k] &&
990  _dep_matl_props.find(dep_matl_props[k]) == _dep_matl_props.end())
991  {
992  mooseError("A VolumetricModel depends on " + dep_matl_props[k] +
993  ", but that material property was not given in the dep_matl_props line.");
994  }
995  }
996  _volumetric_models[_block_id[i]].push_back(vm);
997  }
998  }
999 
1000  for (std::map<SubdomainID, MooseSharedPointer<ConstitutiveModel>>::iterator iter =
1001  _constitutive_model.begin();
1002  iter != _constitutive_model.end();
1003  ++iter)
1004  {
1005  iter->second->initialSetup();
1006  }
1007 
1008  if (isParamValid("constitutive_model") && !_constitutive_active)
1009  {
1010  // User-defined name of the constitutive model (a Material object)
1011  std::string constitutive_model = getParam<std::string>("constitutive_model") + suffix;
1012 
1013  for (unsigned int j = 0; j < mats.size(); ++j)
1014  {
1015  MooseSharedPointer<ConstitutiveModel> cm =
1016  MooseSharedNamespace::dynamic_pointer_cast<ConstitutiveModel>(mats[j]);
1017 
1018  if (cm && cm->name() == constitutive_model)
1019  {
1020  _constitutive_model[_block_id[i]] = cm;
1021  set_constitutive_active = true;
1022  break;
1023  }
1024  }
1025 
1026  if (!set_constitutive_active)
1027  mooseError("Unable to find constitutive model " + constitutive_model);
1028  }
1029  }
1030  if (set_constitutive_active)
1031  _constitutive_active = true;
1032 
1034  {
1035  // Make sure that timeDerivative is supported for _alpha_function. If not, it will error out.
1036  Point dummy_point;
1037  Real dummy_temp = 0;
1038  _alpha_function->timeDerivative(dummy_temp, dummy_point);
1039  }
1040 }
1041 
1043 
1044 void
1046 {
1047  bool cracking_locally_active(false);
1048  if (_cracking_stress > 0)
1049  {
1050  // Compute whether cracking has occurred
1051  (*_crack_rotation)[_qp] = (*_crack_rotation_old)[_qp];
1052 
1053  ColumnMajorMatrix RT((*_crack_rotation)[_qp].transpose());
1054  SymmTensor ePrime;
1055  rotateSymmetricTensor(RT, _elastic_strain[_qp], ePrime);
1056 
1057  for (unsigned int i(0); i < 3; ++i)
1058  {
1059  (*_crack_max_strain)[_qp](i) = (*_crack_max_strain_old)[_qp](i);
1060 
1061  if (_cracking_neg_fraction == 0 && ePrime(i, i) < 0)
1062  {
1063  _crack_flags_local(i) = 1;
1064  }
1065  else if (_cracking_neg_fraction > 0 &&
1066  (*_crack_strain)[_qp](i) * _cracking_neg_fraction > ePrime(i, i))
1067  {
1068  if (-(*_crack_strain)[_qp](i) * _cracking_neg_fraction > ePrime(i, i))
1069  {
1070  _crack_flags_local(i) = 1;
1071  }
1072  else
1073  {
1074  // s = a*e^2 + b*e + c
1075  // a = (Ec-Eo)/(4etr)
1076  // b = (Ec+Eo)/2
1077  // c = (Ec-Eo)*etr/4
1078  // etr = _cracking_neg_fraction * strain when crack occurred
1079  const Real etr = _cracking_neg_fraction * (*_crack_strain)[_qp](i);
1080  const Real Eo = _cracking_stress / (*_crack_strain)[_qp](i);
1081  const Real Ec = Eo * (*_crack_flags_old)[_qp](i);
1082  const Real a = (Ec - Eo) / (4 * etr);
1083  const Real b = (Ec + Eo) / 2;
1084  // Compute the ratio of the current transition stiffness to the original stiffness
1085  _crack_flags_local(i) = (2 * a * etr + b) / Eo;
1086  cracking_locally_active = true;
1087  }
1088  }
1089  else
1090  {
1091  _crack_flags_local(i) = (*_crack_flags_old)[_qp](i);
1092  if (_crack_flags_local(i) < 1)
1093  {
1094  cracking_locally_active = true;
1095  }
1096  }
1097  }
1098  }
1099  if (cracking_locally_active)
1100  {
1101  // Adjust the elasticity matrix for cracking. This must be used by the
1102  // constitutive law.
1103  if (_compute_method == "ShearRetention")
1105  else
1107 
1108  ColumnMajorMatrix R_9x9(9, 9);
1109  const ColumnMajorMatrix & R((*_crack_rotation)[_qp]);
1112 
1114  }
1115 }
1116 
1118 
1119 void
1120 SolidModel::applyCracksToTensor(SymmTensor & tensor, const RealVectorValue & sigma)
1121 {
1122  // Form transformation matrix R*E*R^T
1123  const ColumnMajorMatrix & R((*_crack_rotation)[_qp]);
1124 
1125  // Rotate to crack frame
1126  rotateSymmetricTensor(R.transpose(), tensor, tensor);
1127 
1128  // Reset stress if cracked
1129  if ((*_crack_flags)[_qp](0) < 1)
1130  {
1131  tensor(0, 0) = sigma(0);
1132  }
1133  if ((*_crack_flags)[_qp](1) < 1)
1134  {
1135  tensor(1, 1) = sigma(1);
1136  }
1137  if ((*_crack_flags)[_qp](2) < 1)
1138  {
1139  tensor(2, 2) = sigma(2);
1140  }
1141 
1142  // Rotate back to global frame
1143  rotateSymmetricTensor(R, tensor, tensor);
1144 }
1145 
1147 
1148 void
1149 SolidModel::computeCrackStrainAndOrientation(ColumnMajorMatrix & principal_strain)
1150 {
1151  // The rotation tensor is ordered such that known dirs appear last in the list of
1152  // columns. So, if one dir is known, it corresponds with the last column in the
1153  // rotation tensor.
1154  //
1155  // This convention is based on the eigen routine returning eigen values in
1156  // ascending order.
1157  const unsigned int numKnownDirs = getNumKnownCrackDirs();
1158 
1160 
1161  (*_crack_rotation)[_qp] = (*_crack_rotation_old)[_qp];
1162 
1163  if (numKnownDirs == 0)
1164  {
1165  ColumnMajorMatrix e_vec(3, 3);
1166  _elastic_strain[_qp].columnMajorMatrix().eigen(principal_strain, e_vec);
1167  // If the elastic strain is beyond the cracking strain, save the eigen vectors as
1168  // the rotation tensor.
1169  (*_crack_rotation)[_qp] = e_vec;
1170  }
1171  else if (numKnownDirs == 1)
1172  {
1173  // This is easily the most complicated case.
1174  // 1. Rotate the elastic strain to the orientation associated with the known
1175  // crack.
1176  // 2. Extract the upper 2x2 diagonal block into a separate tensor.
1177  // 3. Run the eigen solver on the result.
1178  // 4. Update the rotation tensor to reflect the effect of the 2 eigenvectors.
1179 
1180  // 1.
1181  ColumnMajorMatrix RT((*_crack_rotation)[_qp].transpose());
1182  SymmTensor ePrime;
1183  rotateSymmetricTensor(RT, _elastic_strain[_qp], ePrime);
1184 
1185  // 2.
1186  ColumnMajorMatrix e2x2(2, 2);
1187  e2x2(0, 0) = ePrime(0, 0);
1188  e2x2(1, 0) = ePrime(1, 0);
1189  e2x2(0, 1) = ePrime(0, 1);
1190  e2x2(1, 1) = ePrime(1, 1);
1191 
1192  // 3.
1193  ColumnMajorMatrix e_val2x1(2, 1);
1194  ColumnMajorMatrix e_vec2x2(2, 2);
1195  e2x2.eigen(e_val2x1, e_vec2x2);
1196 
1197  // 4.
1198  ColumnMajorMatrix e_vec(3, 3);
1199  e_vec(0, 0) = e_vec2x2(0, 0);
1200  e_vec(1, 0) = e_vec2x2(1, 0);
1201  e_vec(2, 0) = 0;
1202  e_vec(0, 1) = e_vec2x2(0, 1);
1203  e_vec(1, 1) = e_vec2x2(1, 1);
1204  e_vec(2, 1) = 0;
1205  e_vec(2, 0) = 0;
1206  e_vec(2, 1) = 0;
1207  e_vec(2, 2) = 1;
1208  (*_crack_rotation)[_qp] = (*_crack_rotation_old)[_qp] * e_vec;
1209 
1210  principal_strain(0, 0) = e_val2x1(0, 0);
1211  principal_strain(1, 0) = e_val2x1(1, 0);
1212  principal_strain(2, 0) = ePrime(2, 2);
1213  }
1214  else if (numKnownDirs == 2 || numKnownDirs == 3)
1215  {
1216  // Rotate to cracked orientation and pick off the strains in the rotated
1217  // coordinate directions.
1218  ColumnMajorMatrix RT((*_crack_rotation)[_qp].transpose());
1219  SymmTensor ePrime;
1220  rotateSymmetricTensor(RT, _elastic_strain[_qp], ePrime);
1221  principal_strain(0, 0) = ePrime.xx();
1222  principal_strain(1, 0) = ePrime.yy();
1223  principal_strain(2, 0) = ePrime.zz();
1224  }
1225  else
1226  {
1227  mooseError("Invalid number of known crack directions");
1228  }
1229 }
1230 
1232 
1233 void
1235 {
1236  if (_cracking_stress_function != NULL)
1237  {
1238  _cracking_stress = _cracking_stress_function->value(_t, _q_point[_qp]);
1239  }
1240 
1241  if (_cracking_stress > 0)
1242  {
1243 
1245 
1246  for (unsigned i(0); i < 3; ++i)
1247  {
1248  if (_principal_strain(i, 0) > (*_crack_max_strain_old)[_qp](i))
1249  {
1250  (*_crack_max_strain)[_qp](i) = _principal_strain(i, 0);
1251  }
1252  }
1253 
1254  // Check for new cracks.
1255  // This must be done in the crack-local coordinate frame.
1256 
1257  // Rotate stress to cracked orientation.
1258  ColumnMajorMatrix R((*_crack_rotation)[_qp]);
1259  ColumnMajorMatrix RT((*_crack_rotation)[_qp].transpose());
1260  SymmTensor sigmaPrime;
1261  rotateSymmetricTensor(RT, _stress[_qp], sigmaPrime);
1262 
1263  unsigned int num_cracks(0);
1264  for (unsigned i(0); i < 3; ++i)
1265  {
1266  _crack_flags_local(i) = 1;
1267  (*_crack_strain)[_qp](i) = (*_crack_strain_old)[_qp](i);
1268  if ((*_crack_flags_old)[_qp](i) < 1)
1269  {
1270  ++num_cracks;
1271  }
1272  }
1273 
1274  bool new_crack(false);
1275  bool cracked(false);
1276  RealVectorValue sigma;
1277  for (unsigned i(0); i < 3; ++i)
1278  {
1279  sigma(i) = sigmaPrime(i, i);
1280  (*_crack_flags)[_qp](i) = (*_crack_flags_old)[_qp](i);
1281  if (sigma(i) <= 1e-4)
1282  {
1283  if ((*_crack_flags)[_qp](i) == 1)
1284  {
1285  (*_crack_max_strain)[_qp](i) = _principal_strain(i, 0);
1286  }
1287  }
1288 
1289  // _cracked_this_step[_q_point[_qp]] = 0;
1290  Real crackFactor(1);
1291  if (_cracking_release == CR_POWER)
1292  {
1293  (*_crack_count)[_qp](i) = (*_crack_count_old)[_qp](i);
1294  }
1295  if ((_cracking_release == CR_POWER && sigma(i) > _cracking_stress &&
1296  _active_crack_planes[i] == 1
1297  // && (*_crack_count)[_qp](i) == 0
1298  )
1299  // || _cracked_this_step_count[_q_point[_qp]] > 5
1300  )
1301  {
1302  cracked = true;
1303  ++((*_crack_count)[_qp](i));
1304  // _cracked_this_step[_q_point[_qp]] = 1;
1305  // Assume Poisson's ratio drops to zero for this direction. Stiffness is then Young's
1306  // modulus.
1307  const Real stiff = _youngs_modulus_function
1308  ? _youngs_modulus_function->value(_temperature[_qp], Point())
1309  : _youngs_modulus;
1310 
1311  if ((*_crack_count_old)[_qp](i) == 0)
1312  {
1313  new_crack = true;
1314  ++num_cracks;
1315 
1316  (*_crack_strain)[_qp](i) = _cracking_stress / stiff;
1317  }
1318  // Compute stress, factor....
1319  (*_crack_flags)[_qp](i) *= 1. / 3.;
1320 
1321  if ((*_crack_max_strain)[_qp](i) < (*_crack_strain)[_qp](i))
1322  {
1323  (*_crack_max_strain)[_qp](i) = (*_crack_strain)[_qp](i);
1324  }
1325  sigma(i) = (*_crack_flags)[_qp](i) * stiff * _principal_strain(i, 0);
1326  }
1327  else if ((_cracking_release != CR_POWER && (*_crack_flags_old)[_qp](i) == 1 &&
1328  sigma(i) > _cracking_stress && num_cracks < _max_cracks &&
1329  _active_crack_planes[i] == 1)
1330  // || _cracked_this_step_count[_q_point[_qp]] > 5
1331  )
1332  {
1333  // A new crack
1334  // _cracked_this_step[_q_point[_qp]] = 1;
1335 
1336  cracked = true;
1337  new_crack = true;
1338  ++num_cracks;
1339 
1340  // Assume Poisson's ratio drops to zero for this direction. Stiffness is then Young's
1341  // modulus.
1342  const Real stiff = _youngs_modulus_function
1343  ? _youngs_modulus_function->value(_temperature[_qp], Point())
1344  : _youngs_modulus;
1345 
1346  (*_crack_strain)[_qp](i) = _cracking_stress / stiff;
1347  if ((*_crack_max_strain)[_qp](i) < (*_crack_strain)[_qp](i))
1348  {
1349  (*_crack_max_strain)[_qp](i) = (*_crack_strain)[_qp](i);
1350  }
1351 
1352  crackFactor = computeCrackFactor(i, sigma(i), (*_crack_flags)[_qp](i));
1353 
1354  (*_crack_flags)[_qp](i) = crackFactor;
1355  _crack_flags_local(i) = crackFactor;
1356 
1357  // May want to set the old value. This may help with the nonlinear solve
1358  // since the stress cannot bounce between just below the critical stress and
1359  // effectively zero. However, this may set allow cracking prematurely.
1360  // (*_crack_flags_old)[_qp](i) = crackFactor;
1361  // (*_crack_strain_old)[_qp](i) = _principal_strain(i,0);
1362  }
1363  else if (_cracking_release != CR_POWER && (*_crack_flags_old)[_qp](i) < 1 &&
1364  std::abs(_principal_strain(i, 0) - (*_crack_max_strain)[_qp](i)) < 1e-10)
1365  {
1366  // Previously cracked,
1367  // Crack opening
1368  cracked = true;
1369  crackFactor = computeCrackFactor(i, sigma(i), (*_crack_flags)[_qp](i));
1370  (*_crack_flags)[_qp](i) = crackFactor;
1371  _crack_flags_local(i) = crackFactor;
1372  }
1373  else if (_cracking_neg_fraction > 0 &&
1376  {
1377  // s = a*e^2 + b*e + c
1378  // a = (Ec-Eo)/(4etr)
1379  // b = (Ec+Eo)/2
1380  // c = (Ec-Eo)*etr/4
1381  // etr = _cracking_neg_fraction * strain when crack occurred
1382  cracked = true;
1383  const Real etr = _cracking_neg_fraction * (*_crack_strain)[_qp](i);
1384  const Real Eo = _cracking_stress / (*_crack_strain)[_qp](i);
1385  const Real Ec = Eo * (*_crack_flags_old)[_qp](i);
1386  const Real a = (Ec - Eo) / (4 * etr);
1387  const Real b = (Ec + Eo) / 2;
1388  const Real c = (Ec - Eo) * etr / 4;
1389  sigma(i) = (a * _principal_strain(i, 0) + b) * _principal_strain(i, 0) + c;
1390  }
1391  }
1392 
1393  if (!new_crack)
1394  {
1395  (*_crack_rotation)[_qp] = (*_crack_rotation_old)[_qp];
1396  }
1397  if (cracked)
1398  {
1399  applyCracksToTensor(_stress[_qp], sigma);
1400  }
1401  }
1402 }
1403 
1404 Real
1405 SolidModel::computeCrackFactor(int i, Real & sigma, Real & flagVal)
1406 {
1408  {
1409  if ((*_crack_max_strain)[_qp](i) < (*_crack_strain)[_qp](i))
1410  {
1411  std::stringstream err;
1412  err << "Max strain less than crack strain: " << i << " " << sigma << ", "
1413  << (*_crack_max_strain)[_qp](i) << ", " << (*_crack_strain)[_qp](i) << ", "
1414  << _principal_strain(0, 0) << ", " << _principal_strain(1, 0) << ", "
1415  << _principal_strain(2, 0) << _elastic_strain[_qp] << std::endl;
1416  mooseError(err.str());
1417  }
1418  const Real crackMaxStrain((*_crack_max_strain)[_qp](i));
1419  // Compute stress that follows exponental curve
1423  (crackMaxStrain - (*_crack_strain)[_qp](i))));
1424  // Compute ratio of current stiffness to original stiffness
1425  flagVal = sigma * (*_crack_strain)[_qp](i) / (crackMaxStrain * _cracking_stress);
1426  }
1427  else
1428  {
1429  if (_cracking_residual_stress == 0)
1430  {
1431  const Real tiny(1e-16);
1432  flagVal = tiny;
1433  sigma = tiny * (*_crack_strain)[_qp](i) * _youngs_modulus;
1434  }
1435  else
1436  {
1438  flagVal = sigma / ((*_crack_max_strain)[_qp](i) * _youngs_modulus);
1439  }
1440  }
1441  if (flagVal < 0)
1442  {
1443  std::stringstream err;
1444  err << "Negative crack flag found: " << i << " " << flagVal << ", "
1445  << (*_crack_max_strain)[_qp](i) << ", " << (*_crack_strain)[_qp](i) << ", " << std::endl;
1446  mooseError(err.str());
1447  }
1448  return flagVal;
1449 }
1450 
1451 unsigned int
1453 {
1454  const unsigned fromElement = _element->getNumKnownCrackDirs();
1455  unsigned int retVal(0);
1456  for (unsigned int i(0); i < 3 - fromElement; ++i)
1457  {
1458  retVal += ((*_crack_flags_old)[_qp](i) < 1);
1459  }
1460  return retVal + fromElement;
1461 }
1462 
1465 {
1466  std::string mat_name = name();
1467  InputParameters parameters = emptyInputParameters();
1468  parameters += this->parameters();
1469 
1471 
1472  std::string formulation = getParam<MooseEnum>("formulation");
1473  std::transform(formulation.begin(), formulation.end(), formulation.begin(), ::tolower);
1474  if (formulation == "nonlinear3d")
1475  {
1476  if (!isCoupled("disp_x") || !isCoupled("disp_y") || !isCoupled("disp_z"))
1477  mooseError("Nonlinear3D requires all three displacements");
1478 
1479  if (isCoupled("disp_r"))
1480  mooseError("Linear must not define disp_r");
1481 
1482  if (_coord_type == Moose::COORD_RZ)
1483  mooseError("Nonlinear3D formulation requested for coord_type = RZ problem");
1484 
1485  element = new SolidMechanics::Nonlinear3D(*this, mat_name, parameters);
1486  }
1487  else if (formulation == "nonlinearrz")
1488  {
1489  if (!isCoupled("disp_r") || !isCoupled("disp_z"))
1490  mooseError("NonlinearRZ must define disp_r and disp_z");
1491 
1492  element = new SolidMechanics::NonlinearRZ(*this, mat_name, parameters);
1493  }
1494  else if (formulation == "axisymmetricrz")
1495  {
1496  if (!isCoupled("disp_r") || !isCoupled("disp_z"))
1497  mooseError("AxisymmetricRZ must define disp_r and disp_z");
1498  element = new SolidMechanics::AxisymmetricRZ(*this, mat_name, parameters);
1499  }
1500  else if (formulation == "sphericalr")
1501  {
1502  if (!isCoupled("disp_r"))
1503  mooseError("SphericalR must define disp_r");
1504  element = new SolidMechanics::SphericalR(*this, mat_name, parameters);
1505  }
1506  else if (formulation == "planestrain")
1507  {
1508  if (!isCoupled("disp_x") || !isCoupled("disp_y"))
1509  mooseError("PlaneStrain must define disp_x and disp_y");
1510  element = new SolidMechanics::PlaneStrain(*this, mat_name, parameters);
1511  }
1512  else if (formulation == "nonlinearplanestrain")
1513  {
1514  if (!isCoupled("disp_x") || !isCoupled("disp_y"))
1515  mooseError("NonlinearPlaneStrain must define disp_x and disp_y");
1516  element = new SolidMechanics::NonlinearPlaneStrain(*this, mat_name, parameters);
1517  }
1518  else if (formulation == "linear")
1519  {
1520  if (isCoupled("disp_r"))
1521  mooseError("Linear must not define disp_r");
1522  if (_coord_type == Moose::COORD_RZ)
1523  mooseError("Linear formulation requested for coord_type = RZ problem");
1524  element = new SolidMechanics::Linear(*this, mat_name, parameters);
1525  }
1526  else if (formulation != "")
1527  mooseError("Unknown formulation: " + formulation);
1528 
1529  if (!element && _coord_type == Moose::COORD_RZ)
1530  {
1531  if (!isCoupled("disp_r") || !isCoupled("disp_z"))
1532  {
1533  std::string err(name());
1534  err += ": RZ coord sys requires disp_r and disp_z for AxisymmetricRZ formulation";
1535  mooseError(err);
1536  }
1537  element = new SolidMechanics::AxisymmetricRZ(*this, mat_name, parameters);
1538  }
1539  else if (!element && _coord_type == Moose::COORD_RSPHERICAL)
1540  {
1541  if (!isCoupled("disp_r"))
1542  {
1543  std::string err(name());
1544  err += ": RSPHERICAL coord sys requires disp_r for SphericalR formulation";
1545  mooseError(err);
1546  }
1547  element = new SolidMechanics::SphericalR(*this, mat_name, parameters);
1548  }
1549 
1550  if (!element)
1551  {
1552  if (isCoupled("disp_x") && isCoupled("disp_y") && isCoupled("disp_z"))
1553  {
1554  if (isCoupled("disp_r"))
1555  mooseError("Error with displacement specification in material " + mat_name);
1556  element = new SolidMechanics::Nonlinear3D(*this, mat_name, parameters);
1557  }
1558  else if (isCoupled("disp_x") && isCoupled("disp_y"))
1559  {
1560  if (isCoupled("disp_r"))
1561  mooseError("Error with displacement specification in material " + mat_name);
1562  element = new SolidMechanics::PlaneStrain(*this, mat_name, parameters);
1563  }
1564  else if (isCoupled("disp_r") && isCoupled("disp_z"))
1565  {
1566  if (_coord_type != Moose::COORD_RZ)
1567  mooseError("RZ coord system not specified, but disp_r and disp_z are");
1568  element = new SolidMechanics::AxisymmetricRZ(*this, mat_name, parameters);
1569  }
1570  else if (isCoupled("disp_r"))
1571  {
1572  if (_coord_type != Moose::COORD_RSPHERICAL)
1573  mooseError("RSPHERICAL coord system not specified, but disp_r is");
1574  element = new SolidMechanics::SphericalR(*this, mat_name, parameters);
1575  }
1576  else if (isCoupled("disp_x"))
1577  element = new SolidMechanics::Linear(*this, mat_name, parameters);
1578  else
1579  mooseError("Unable to determine formulation for material " + mat_name);
1580  }
1581 
1582  mooseAssert(element, "No Element created for material " + mat_name);
1583 
1584  return element;
1585 }
1586 
1587 void
1588 SolidModel::createConstitutiveModel(const std::string & cm_name)
1589 {
1590 
1591  Factory & factory = _app.getFactory();
1592  InputParameters params = factory.getValidParams(cm_name);
1593  // These set_attributes calls are to make isParamSetByUser() work correctly on
1594  // these parameters in the ConstitutiveModel class, and are needed only for the
1595  // legacy_return_mapping option.
1596  params.set_attributes("absolute_tolerance", false);
1597  params.set_attributes("relative_tolerance", false);
1598  params.set_attributes("max_its", false);
1599  params += parameters();
1600  MooseSharedPointer<ConstitutiveModel> cm =
1601  factory.create<ConstitutiveModel>(cm_name, name() + "Model", params, _tid);
1602 
1603  _models_to_free.insert(
1604  cm); // Keep track of the dynamic memory that is created internally to this object
1605 
1606  _constitutive_active = true;
1607  for (unsigned i(0); i < _block_id.size(); ++i)
1608  {
1609  _constitutive_model[_block_id[i]] = cm;
1610  }
1611 }
1612 
1613 void
1615 {
1616  for (_qp = 0; _qp < n_points; ++_qp)
1617  {
1619  }
1621  {
1622  const SubdomainID current_block = _current_elem->subdomain_id();
1623  MooseSharedPointer<ConstitutiveModel> cm = _constitutive_model[current_block];
1624  cm->initStatefulProperties(n_points);
1625  }
1626 }
1627 
1628 void
1630 {
1631  mooseAssert(_J_thermal_term_vec, "_J_thermal_term_vec not initialized");
1632 
1633  Real stress_trace = _stress[_qp].xx() + _stress[_qp].yy() + _stress[_qp].zz();
1634 
1636  for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
1637  {
1638  Real dthermstrain_dx =
1640  (*_J_thermal_term_vec)[_qp](i) = stress_trace * dthermstrain_dx;
1641  }
1642 }
1643 
1644 void
1646 {
1648  "_current_instantaneous_thermal_expansion_coef not initialized");
1649 
1650  (*_current_instantaneous_thermal_expansion_coef)[_qp] = 0.0;
1651 
1652  if (_alpha_function)
1653  {
1654  Point p;
1655  Real current_temp = _temperature[_qp];
1656 
1657  if (!_mean_alpha_function)
1658  {
1659  Real alpha = _alpha_function->value(current_temp, p);
1660  (*_current_instantaneous_thermal_expansion_coef)[_qp] = alpha;
1661  }
1662  else
1663  {
1664  Real small(1e-6);
1665  Real dalphabar_dT = _alpha_function->timeDerivative(current_temp, p);
1666  Real alphabar_Tsf = _alpha_function->value(_stress_free_temp, p);
1667  Real alphabar = _alpha_function->value(current_temp, p);
1668  Real numerator = dalphabar_dT * (current_temp - _ref_temp) + alphabar;
1669  Real denominator = 1.0 + alphabar_Tsf * (_stress_free_temp - _ref_temp);
1670  if (denominator < small)
1671  mooseError("Denominator too small in thermal strain calculation");
1672  (*_current_instantaneous_thermal_expansion_coef)[_qp] = numerator / denominator;
1673  }
1674  }
1675  else
1677 }
virtual std::vector< std::string > getDependentMaterialProperties() const
This class defines a basic set of capabilities any elasticity tensor should have. ...
SymmElasticityTensor * elasticityTensor() const
Definition: SolidModel.h:217
std::set< std::string > _dep_matl_props
Definition: SolidModel.h:104
MaterialProperty< ColumnMajorMatrix > * _crack_rotation
Definition: SolidModel.h:125
const MaterialProperty< RealVectorValue > * _crack_count_old
Definition: SolidModel.h:124
Real _poissons_ratio
Definition: SolidModel.h:70
const MaterialProperty< Real > * _SED_old
Definition: SolidModel.h:151
MaterialProperty< Real > * _SED
Definition: SolidModel.h:150
const bool _compute_JIntegral
Definition: SolidModel.h:145
const unsigned int _max_cracks
Definition: SolidModel.h:86
Real yy() const
Definition: SymmTensor.h:130
const Real _cracking_residual_stress
Definition: SolidModel.h:79
Real _cracking_stress
Definition: SolidModel.h:78
const Real _alpha
Definition: SolidModel.h:95
SymmElasticityTensor * _local_elasticity_tensor
Definition: SolidModel.h:265
bool _constitutive_active
Definition: SolidModel.h:246
SymmTensor _total_strain_increment
Definition: SolidModel.h:142
Element is the base class for all of this module&#39;s solid mechanics element formulations.
Definition: Element.h:23
MaterialProperty< SymmElasticityTensor > & _Jacobian_mult
Definition: SolidModel.h:134
virtual void computeDeformationGradient(unsigned int, ColumnMajorMatrix &)
Definition: Element.h:48
virtual Real volumeRatioOld(unsigned) const
Definition: Element.h:58
virtual void modifyStrainIncrement()
Modify increment for things like thermal strain.
Definition: SolidModel.C:556
MaterialProperty< RealVectorValue > * _crack_count
Definition: SolidModel.h:123
const bool _compute_InteractionIntegral
Definition: SolidModel.h:146
Nonlinear3D is the base class for all 3D nonlinear solid mechanics material models.
Definition: Nonlinear3D.h:22
virtual void computeStrain(const unsigned qp, const SymmTensor &total_strain_old, SymmTensor &total_strain_new, SymmTensor &strain_increment)=0
const VariableValue & _temperature_old
Definition: SolidModel.h:93
virtual bool updateElasticityTensor(SymmElasticityTensor &tensor)
Return true if the elasticity tensor changed.
Definition: SolidModel.C:888
void applyCracksToTensor(SymmTensor &tensor, const RealVectorValue &sigma)
Definition: SolidModel.C:1120
virtual void checkElasticConstants()
Definition: SolidModel.C:373
Real xx() const
Definition: SymmTensor.h:129
Real _stress_free_temp
Definition: SolidModel.h:99
Real _bulk_modulus
Definition: SolidModel.h:68
virtual void computeEshelby()
Definition: SolidModel.C:786
bool _youngs_modulus_set
Definition: SolidModel.h:66
bool _poissons_ratio_set
Definition: SolidModel.h:64
virtual Real computeCrackFactor(int i, Real &sigma, Real &flagVal)
Definition: SolidModel.C:1405
virtual void computePreconditioning()
Definition: SolidModel.C:941
virtual void initQpStatefulProperties()
Definition: SolidModel.C:692
Function * _youngs_modulus_function
Definition: SolidModel.h:74
const CRACKING_RELEASE _cracking_release
Definition: SolidModel.h:77
Function *const _cracking_stress_function
Definition: SolidModel.h:82
virtual void createElasticityTensor()
Definition: SolidModel.C:496
MaterialProperty< SymmTensor > & _total_strain
Definition: SolidModel.h:114
std::vector< SubdomainID > _block_id
Definition: SolidModel.h:241
const MaterialProperty< SymmTensor > & _total_strain_old
Definition: SolidModel.h:115
const MaterialProperty< SymmTensor > & _stress_old_prop
Definition: SolidModel.h:109
virtual unsigned int getNumKnownCrackDirs() const
Definition: Element.h:63
virtual void crackingStressRotation()
Definition: SolidModel.C:1234
virtual void computeCurrentInstantaneousThermalExpansionCoefficient()
Definition: SolidModel.C:1645
virtual void adjustForCrackingWithShearRetention(const RealVectorValue &crack_flags)
MaterialProperty< SymmTensor > & _d_stress_dT
Definition: SolidModel.h:140
virtual void adjustForCracking(const RealVectorValue &crack_flags)
MaterialProperty< SymmElasticityTensor > & _elasticity_tensor
Definition: SolidModel.h:133
const std::string _compute_method
Definition: SolidModel.h:81
virtual void elementInit()
Definition: SolidModel.h:172
Real zz() const
Definition: SymmTensor.h:131
void zero()
Definition: SymmTensor.h:273
bool & _step_zero
Restartable data to check for the zeroth and first time steps for thermal calculations.
Definition: SolidModel.h:254
const VariableGradient & _temp_grad
Definition: SolidModel.h:94
MaterialProperty< RealVectorValue > * _J_thermal_term_vec
Definition: SolidModel.h:153
bool & _step_one
Definition: SolidModel.h:255
bool _mean_alpha_function
Definition: SolidModel.h:100
bool _shear_modulus_set
Definition: SolidModel.h:65
const MaterialProperty< RealVectorValue > * _crack_flags_old
Definition: SolidModel.h:121
void computeCrackStrainAndOrientation(ColumnMajorMatrix &principal_strain)
Definition: SolidModel.C:1149
virtual void applyThermalStrain()
Definition: SolidModel.C:587
Real _cracking_alpha
Definition: SolidModel.h:84
Function * _poissons_ratio_function
Definition: SolidModel.h:75
static Real detMatrix(const ColumnMajorMatrix &A)
Definition: Element.C:28
const Real _cracking_neg_fraction
Definition: SolidModel.h:87
const bool _has_temp
Definition: SolidModel.h:91
SymmTensor _d_strain_dT
Definition: SolidModel.h:137
std::set< MooseSharedPointer< ConstitutiveModel > > _models_to_free
Definition: SolidModel.h:245
bool _bulk_modulus_set
Definition: SolidModel.h:62
InputParameters validParams< SolidModel >()
Definition: SolidModel.C:27
virtual void computeProperties()
Definition: SolidModel.C:725
Moose::CoordinateSystemType _coord_type
Definition: SolidModel.h:58
virtual void applyVolumetricStrain()
Definition: SolidModel.C:645
void calculate(unsigned int qp)
Public function that will be called whenever the values for this matrix need to be filled in...
std::map< SubdomainID, MooseSharedPointer< ConstitutiveModel > > _constitutive_model
Definition: SolidModel.h:243
Real _youngs_modulus
Definition: SolidModel.h:72
virtual void initStatefulProperties(unsigned n_points)
Definition: SolidModel.C:1614
MaterialProperty< SymmTensor > & _elastic_strain
Definition: SolidModel.h:117
void rotateFromLocalToGlobal(const ColumnMajorMatrix &R)
MaterialProperty< RealVectorValue > * _crack_flags
Definition: SolidModel.h:120
ColumnMajorMatrix _principal_strain
Definition: SolidModel.h:131
virtual void computeConstitutiveModelStress()
Compute the stress (sigma += deltaSigma)
Definition: SolidModel.C:835
SymmTensor _strain_increment
Definition: SolidModel.h:143
SymmTensor _stress_old
Definition: SolidModel.h:112
NonlinearPlaneStrain is a class for large deformation plane strain.
void form9x9Rotation(const ColumnMajorMatrix &R_3x3, ColumnMajorMatrix &R_9x9) const
const MaterialProperty< ColumnMajorMatrix > * _crack_rotation_old
Definition: SolidModel.h:126
Real xy() const
Definition: SymmTensor.h:132
const MaterialProperty< RealVectorValue > * _crack_max_strain_old
Definition: SolidModel.h:130
const VariableValue & _temperature
Definition: SolidModel.h:92
const MaterialProperty< SymmTensor > & _elastic_strain_old
Definition: SolidModel.h:118
MaterialProperty< RealVectorValue > * _crack_max_strain
Definition: SolidModel.h:129
const MaterialProperty< RealVectorValue > * _crack_strain_old
Definition: SolidModel.h:128
bool _has_stress_free_temp
Definition: SolidModel.h:98
SolidMechanics::Element * _element
Definition: SolidModel.h:263
virtual void finalizeStress(std::vector< SymmTensor * > &)
Rotate stress to current configuration.
Definition: Element.h:61
Real yz() const
Definition: SymmTensor.h:133
NonlinearRZ is the base class for all RZ nonlinear solid mechanics material models.
Definition: NonlinearRZ.h:18
Real zx() const
Definition: SymmTensor.h:134
MaterialProperty< RealVectorValue > * _crack_strain
Definition: SolidModel.h:127
SolidMechanics::Element * createElement()
Definition: SolidModel.C:1464
MaterialProperty< SymmTensor > & _stress
Definition: SolidModel.h:106
MaterialProperty< Real > * _current_instantaneous_thermal_expansion_coef
Definition: SolidModel.h:156
virtual void computeStrainEnergyDensity()
Definition: SolidModel.C:775
virtual void finalizeStress()
Rotate stress to current configuration.
Definition: SolidModel.C:929
virtual void initialSetup()
Definition: SolidModel.C:954
virtual void computeThermalJvec()
Definition: SolidModel.C:1629
Function * _alpha_function
Definition: SolidModel.h:96
Real _ref_temp
Definition: SolidModel.h:101
virtual ~SolidModel()
Definition: SolidModel.C:364
void addDiag(Real value)
Definition: SymmTensor.h:279
virtual void timestepSetup()
Definition: SolidModel.C:518
virtual void crackingStrainDirections()
Determine cracking directions. Rotate elasticity tensor.
Definition: SolidModel.C:1045
virtual unsigned int getNumKnownCrackDirs() const
Definition: SolidModel.C:1452
bool _lambda_set
Definition: SolidModel.h:63
static void invertMatrix(const ColumnMajorMatrix &A, ColumnMajorMatrix &Ainv)
Definition: Element.C:49
void setYoungsModulus(const Real E)
Set the Young&#39;s Modulus.
virtual void jacobianSetup()
Definition: SolidModel.C:529
virtual void init()
Definition: Element.h:46
Real _lambda
Definition: SolidModel.h:69
virtual void computeStress()
Compute the stress (sigma += deltaSigma)
Definition: SolidModel.h:183
std::map< SubdomainID, std::vector< MooseSharedPointer< VolumetricModel > > > _volumetric_models
Definition: SolidModel.h:103
std::vector< unsigned int > _active_crack_planes
Definition: SolidModel.h:85
Real _shear_modulus
Definition: SolidModel.h:71
SolidModel(const InputParameters &parameters)
Definition: SolidModel.C:141
void computeElasticityTensor()
Definition: SolidModel.C:862
void createConstitutiveModel(const std::string &cm_name)
Definition: SolidModel.C:1588
Defines an Isotropic Elasticity Tensor.
const SolidMechanics::Element * element() const
Definition: SolidModel.h:219
MaterialProperty< RankTwoTensor > * _Eshelby_tensor
Definition: SolidModel.h:152
RealVectorValue _crack_flags_local
Definition: SolidModel.h:122
const Real _cracking_beta
Definition: SolidModel.h:80
static void rotateSymmetricTensor(const ColumnMajorMatrix &R, const SymmTensor &T, SymmTensor &result)
Definition: SolidModel.C:659
void setPoissonsRatio(const Real nu)
Set Poissons Ratio.