www.mooseframework.org
ComputeSmearedCrackingStress.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 "ElasticityTensorTools.h"
12 #include "StressUpdateBase.h"
13 #include "Conversion.h"
14 
16 
17 template <>
18 InputParameters
20 {
21  InputParameters params = validParams<ComputeMultipleInelasticStress>();
22  params.addClassDescription("Compute stress using a fixed smeared cracking model");
23  MooseEnum cracking_release("abrupt exponential power", "abrupt");
24  params.addDeprecatedParam<MooseEnum>(
25  "cracking_release",
26  cracking_release,
27  "The cracking release type. 'abrupt' (default) gives an abrupt "
28  "stress release, 'exponential' uses an exponential softening model, "
29  "and 'power' uses a power law",
30  "This is replaced by the use of 'softening_models' together with a separate block defining "
31  "a softening model");
32  params.addParam<std::vector<MaterialName>>(
33  "softening_models",
34  "The material objects used to compute softening behavior for loading a crack."
35  "Either 1 or 3 models must be specified. If a single model is specified, it is"
36  "used for all directions. If 3 models are specified, they will be used for the"
37  "3 crack directions in sequence");
38  params.addDeprecatedParam<Real>(
39  "cracking_residual_stress",
40  0.0,
41  "The fraction of the cracking stress allowed to be maintained following a crack.",
42  "This is replaced by the use of 'softening_models' together with a separate block defining "
43  "a softening model");
44  params.addRequiredCoupledVar(
45  "cracking_stress",
46  "The stress threshold beyond which cracking occurs. Negative values prevent cracking.");
47  MultiMooseEnum direction("x y z");
48  params.addParam<MultiMooseEnum>(
49  "prescribed_crack_directions", direction, "Prescribed directions of first cracks");
50  params.addParam<unsigned int>(
51  "max_cracks", 3, "The maximum number of cracks allowed at a material point.");
52  params.addRangeCheckedParam<Real>("cracking_neg_fraction",
53  0,
54  "cracking_neg_fraction <= 1 & cracking_neg_fraction >= 0",
55  "The fraction of the cracking strain at which "
56  "a transitition begins during decreasing "
57  "strain to the original stiffness.");
58  params.addDeprecatedParam<Real>(
59  "cracking_beta",
60  1.0,
61  "Coefficient used to control the softening in the exponential model. "
62  "When set to 1, the initial softening slope is equal to the negative "
63  "of the Young's modulus. Smaller numbers scale down that slope.",
64  "This is replaced by the use of 'softening_models' together with a separate block defining "
65  "a softening model");
66  params.addParam<Real>(
67  "max_stress_correction",
68  1.0,
69  "Maximum permitted correction to the predicted stress as a ratio of the "
70  "stress change to the predicted stress from the previous step's damage level. "
71  "Values less than 1 will improve robustness, but not be as accurate.");
72 
73  params.addRangeCheckedParam<Real>(
74  "shear_retention_factor",
75  0.0,
76  "shear_retention_factor>=0 & shear_retention_factor<=1.0",
77  "Fraction of original shear stiffness to be retained after cracking");
78  params.set<std::vector<MaterialName>>("inelastic_models") = {};
79 
80  return params;
81 }
82 
84  : ComputeMultipleInelasticStress(parameters),
85  _cracking_release(getParam<MooseEnum>("cracking_release").getEnum<CrackingRelease>()),
86  _cracking_residual_stress(getParam<Real>("cracking_residual_stress")),
87  _cracking_stress(coupledValue("cracking_stress")),
88  _max_cracks(getParam<unsigned int>("max_cracks")),
89  _cracking_neg_fraction(getParam<Real>("cracking_neg_fraction")),
90  _cracking_beta(getParam<Real>("cracking_beta")),
91  _shear_retention_factor(getParam<Real>("shear_retention_factor")),
92  _max_stress_correction(getParam<Real>("max_stress_correction")),
93  _crack_damage(declareProperty<RealVectorValue>(_base_name + "crack_damage")),
94  _crack_damage_old(getMaterialPropertyOld<RealVectorValue>(_base_name + "crack_damage")),
95  _crack_flags(declareProperty<RealVectorValue>(_base_name + "crack_flags")),
96  _crack_rotation(declareProperty<RankTwoTensor>(_base_name + "crack_rotation")),
97  _crack_rotation_old(getMaterialPropertyOld<RankTwoTensor>(_base_name + "crack_rotation")),
98  _crack_initiation_strain(
99  declareProperty<RealVectorValue>(_base_name + "crack_initiation_strain")),
100  _crack_initiation_strain_old(
101  getMaterialPropertyOld<RealVectorValue>(_base_name + "crack_initiation_strain")),
102  _crack_max_strain(declareProperty<RealVectorValue>(_base_name + "crack_max_strain")),
103  _crack_max_strain_old(getMaterialPropertyOld<RealVectorValue>(_base_name + "crack_max_strain"))
104 {
105  MultiMooseEnum prescribed_crack_directions =
106  getParam<MultiMooseEnum>("prescribed_crack_directions");
107  if (prescribed_crack_directions.size() > 0)
108  {
109  if (prescribed_crack_directions.size() > 3)
110  mooseError("A maximum of three crack directions may be specified");
111  for (unsigned int i = 0; i < prescribed_crack_directions.size(); ++i)
112  {
113  for (unsigned int j = 0; j < i; ++j)
114  if (prescribed_crack_directions[i] == prescribed_crack_directions[j])
115  mooseError("Entries in 'prescribed_crack_directions' cannot be repeated");
117  static_cast<unsigned int>(prescribed_crack_directions.get(i)));
118  }
119 
120  // Fill in the last remaining direction if 2 are specified
121  if (_prescribed_crack_directions.size() == 2)
122  {
123  std::set<unsigned int> available_dirs = {0, 1, 2};
124  for (auto dir : _prescribed_crack_directions)
125  if (available_dirs.erase(dir) != 1)
126  mooseError("Invalid prescribed crack direction:" + Moose::stringify(dir));
127  if (available_dirs.size() != 1)
128  mooseError("Error in finding remaining available crack direction");
129  _prescribed_crack_directions.push_back(*available_dirs.begin());
130  }
131  }
132 
133  if (parameters.isParamSetByUser("softening_models"))
134  {
135  if (parameters.isParamSetByUser("cracking_release"))
136  mooseError("In ComputeSmearedCrackingStress cannot specify both 'cracking_release' and "
137  "'softening_models'");
138  if (parameters.isParamSetByUser("cracking_residual_stress"))
139  mooseError("In ComputeSmearedCrackingStress cannot specify both 'cracking_residual_stress' "
140  "and 'softening_models'");
141  if (parameters.isParamSetByUser("cracking_beta"))
142  mooseError("In ComputeSmearedCrackingStress cannot specify both 'cracking_beta' and "
143  "'softening_models'");
144  }
145 }
146 
147 void
149 {
151 
152  _crack_damage[_qp] = 0.0;
153 
154  _crack_initiation_strain[_qp] = 0.0;
155  _crack_max_strain[_qp](0) = 0.0;
156 
157  switch (_prescribed_crack_directions.size())
158  {
159  case 0:
160  {
161  _crack_rotation[_qp] = RankTwoTensor::Identity();
162  break;
163  }
164  case 1:
165  {
166  _crack_rotation[_qp].zero();
167  switch (_prescribed_crack_directions[0])
168  {
169  case 0:
170  {
171  _crack_rotation[_qp](0, 0) = 1.0;
172  _crack_rotation[_qp](1, 1) = 1.0;
173  _crack_rotation[_qp](2, 2) = 1.0;
174  break;
175  }
176  case 1:
177  {
178  _crack_rotation[_qp](1, 0) = 1.0;
179  _crack_rotation[_qp](0, 1) = 1.0;
180  _crack_rotation[_qp](2, 2) = 1.0;
181  break;
182  }
183  case 2:
184  {
185  _crack_rotation[_qp](2, 0) = 1.0;
186  _crack_rotation[_qp](0, 1) = 1.0;
187  _crack_rotation[_qp](1, 2) = 1.0;
188  break;
189  }
190  }
191  break;
192  }
193  case 2:
194  {
195  mooseError("Number of prescribed crack directions cannot be 2");
196  break;
197  }
198  case 3:
199  {
200  for (unsigned int i = 0; i < _prescribed_crack_directions.size(); ++i)
201  {
202  RealVectorValue crack_dir_vec;
203  crack_dir_vec(_prescribed_crack_directions[i]) = 1.0;
204  _crack_rotation[_qp].fillColumn(i, crack_dir_vec);
205  }
206  }
207  }
208 }
209 
210 void
212 {
214 
216  mooseError("ComputeSmearedCrackingStress requires that the elasticity tensor be "
217  "guaranteed isotropic");
218 
219  std::vector<MaterialName> soft_matls = getParam<std::vector<MaterialName>>("softening_models");
220  if (soft_matls.size() != 0)
221  {
222  for (auto soft_matl : soft_matls)
223  {
225  dynamic_cast<SmearedCrackSofteningBase *>(&getMaterialByName(soft_matl));
226  if (scsb)
227  _softening_models.push_back(scsb);
228  else
229  mooseError("Model " + soft_matl +
230  " is not a softening model that can be used with ComputeSmearedCrackingStress");
231  }
232  if (_softening_models.size() == 1)
233  {
234  // Reuse the same model in all 3 directions
235  _softening_models.push_back(_softening_models[0]);
236  _softening_models.push_back(_softening_models[0]);
237  }
238  else if (_softening_models.size() != 3)
239  mooseError("If 'softening_models' is specified in ComputeSmearedCrackingStress, either 1 or "
240  "3 models must be provided");
241  }
242 }
243 
244 void
246 {
247  bool force_elasticity_rotation = false;
248 
249  if (!previouslyCracked())
251  else
252  {
254 
255  // Propagate behavior from the (now inactive) inelastic models
257  for (auto model : _models)
258  {
259  model->setQp(_qp);
260  model->propagateQpStatefulProperties();
261  }
262 
263  // Since the other inelastic models are inactive, they will not limit the time step
264  _matl_timestep_limit[_qp] = std::numeric_limits<Real>::max();
265 
266  // update _local_elasticity_tensor based on cracking state in previous time step
268 
269  // Calculate stress in intermediate configuration
271 
273  force_elasticity_rotation = true;
274  }
275 
276  // compute crack status and adjust stress
278 
280  {
281  finiteStrainRotation(force_elasticity_rotation);
283  }
284 }
285 
286 void
288 {
289  const Real youngs_modulus =
291 
292  bool cracking_locally_active = false;
293 
294  const Real cracking_stress = _cracking_stress[_qp];
295 
296  if (cracking_stress > 0)
297  {
298  RealVectorValue stiffness_ratio_local(1.0, 1.0, 1.0);
299  const RankTwoTensor & R = _crack_rotation_old[_qp];
300  RankTwoTensor ePrime(_elastic_strain_old[_qp]);
301  ePrime.rotate(R.transpose());
302 
303  for (unsigned int i = 0; i < 3; ++i)
304  {
305  // Update elasticity tensor based on crack status of the end of last time step
306  if (_crack_damage_old[_qp](i) > 0.0)
307  {
308  if (_cracking_neg_fraction == 0.0 && MooseUtils::absoluteFuzzyLessThan(ePrime(i, i), 0.0))
309  stiffness_ratio_local(i) = 1.0;
310  else if (_cracking_neg_fraction > 0.0 &&
311  ePrime(i, i) < _crack_initiation_strain_old[_qp](i) * _cracking_neg_fraction &&
312  ePrime(i, i) > -_crack_initiation_strain_old[_qp](i) * _cracking_neg_fraction)
313  {
314  const Real etr = _cracking_neg_fraction * _crack_initiation_strain_old[_qp](i);
315  const Real Eo = cracking_stress / _crack_initiation_strain_old[_qp](i);
316  const Real Ec = Eo * (1.0 - _crack_damage_old[_qp](i));
317  const Real a = (Ec - Eo) / (4 * etr);
318  const Real b = (Ec + Eo) / 2;
319  // Compute the ratio of the current transition stiffness to the original stiffness
320  stiffness_ratio_local(i) = (2.0 * a * etr + b) / Eo;
321  cracking_locally_active = true;
322  }
323  else
324  {
325  stiffness_ratio_local(i) = (1.0 - _crack_damage_old[_qp](i));
326  cracking_locally_active = true;
327  }
328  }
329  }
330 
331  if (cracking_locally_active)
332  {
333  // Update the elasticity tensor in the crack coordinate system
334  const bool c0_coupled = MooseUtils::absoluteFuzzyEqual(stiffness_ratio_local(0), 1.0);
335  const bool c1_coupled = MooseUtils::absoluteFuzzyEqual(stiffness_ratio_local(1), 1.0);
336  const bool c2_coupled = MooseUtils::absoluteFuzzyEqual(stiffness_ratio_local(2), 1.0);
337 
338  const Real c01 = (c0_coupled && c1_coupled ? 1.0 : 0.0);
339  const Real c02 = (c0_coupled && c2_coupled ? 1.0 : 0.0);
340  const Real c12 = (c1_coupled && c2_coupled ? 1.0 : 0.0);
341 
342  const Real c01_shear_retention = (c0_coupled && c1_coupled ? 1.0 : _shear_retention_factor);
343  const Real c02_shear_retention = (c0_coupled && c2_coupled ? 1.0 : _shear_retention_factor);
344  const Real c12_shear_retention = (c1_coupled && c2_coupled ? 1.0 : _shear_retention_factor);
345 
346  // Filling with 9 components is sufficient because these are the only nonzero entries
347  // for isotropic or orthotropic materials.
348  std::vector<Real> local_elastic(9);
349 
350  local_elastic[0] = (c0_coupled ? _elasticity_tensor[_qp](0, 0, 0, 0)
351  : stiffness_ratio_local(0) * youngs_modulus);
352  local_elastic[1] = _elasticity_tensor[_qp](0, 0, 1, 1) * c01;
353  local_elastic[2] = _elasticity_tensor[_qp](0, 0, 2, 2) * c02;
354  local_elastic[3] = (c1_coupled ? _elasticity_tensor[_qp](1, 1, 1, 1)
355  : stiffness_ratio_local(1) * youngs_modulus);
356  local_elastic[4] = _elasticity_tensor[_qp](1, 1, 2, 2) * c12;
357  local_elastic[5] = (c2_coupled ? _elasticity_tensor[_qp](2, 2, 2, 2)
358  : stiffness_ratio_local(2) * youngs_modulus);
359  local_elastic[6] = _elasticity_tensor[_qp](1, 2, 1, 2) * c12_shear_retention;
360  local_elastic[7] = _elasticity_tensor[_qp](0, 2, 0, 2) * c02_shear_retention;
361  local_elastic[8] = _elasticity_tensor[_qp](0, 1, 0, 1) * c01_shear_retention;
362 
363  _local_elasticity_tensor.fillFromInputVector(local_elastic, RankFourTensor::symmetric9);
364 
365  // Rotate the modified elasticity tensor back into global coordinates
366  _local_elasticity_tensor.rotate(R);
367  }
368  }
369  if (!cracking_locally_active)
371 }
372 
373 void
375 {
376  const Real youngs_modulus =
378  const Real cracking_alpha = -youngs_modulus;
379 
380  Real cracking_stress = _cracking_stress[_qp];
381 
382  if (cracking_stress > 0)
383  {
384  // Initializing crack states
386 
387  for (unsigned i = 0; i < 3; ++i)
388  {
389  _crack_max_strain[_qp](i) = _crack_max_strain_old[_qp](i);
391  _crack_damage[_qp](i) = _crack_damage_old[_qp](i);
392  }
393 
394  // Compute crack orientations: updated _crack_rotation[_qp] based on current strain
395  RealVectorValue strain_in_crack_dir;
396  computeCrackStrainAndOrientation(strain_in_crack_dir);
397 
398  for (unsigned i = 0; i < 3; ++i)
399  {
400  if (strain_in_crack_dir(i) > _crack_max_strain[_qp](i))
401  _crack_max_strain[_qp](i) = strain_in_crack_dir(i);
402  }
403 
404  // Check for new cracks.
405  // Rotate stress to cracked orientation.
406  const RankTwoTensor & R = _crack_rotation[_qp];
407  RankTwoTensor sigmaPrime(_stress[_qp]);
408  sigmaPrime.rotate(R.transpose()); // stress in crack coordinates
409 
410  unsigned int num_cracks = 0;
411  for (unsigned int i = 0; i < 3; ++i)
412  {
413  if (_crack_damage_old[_qp](i) > 0.0)
414  ++num_cracks;
415  }
416 
417  bool cracked(false);
418  RealVectorValue sigma;
419  for (unsigned int i = 0; i < 3; ++i)
420  {
421  sigma(i) = sigmaPrime(i, i);
422 
423  Real stiffness_ratio = 1.0 - _crack_damage[_qp](i);
424 
425  const bool pre_existing_crack = (_crack_damage_old[_qp](i) > 0.0);
426  const bool met_stress_criterion = (sigma(i) > cracking_stress);
427  const bool loading_existing_crack = (strain_in_crack_dir(i) >= _crack_max_strain[_qp](i));
428  const bool allowed_to_crack = (pre_existing_crack || num_cracks < _max_cracks);
429  bool new_crack = false;
430 
431  cracked |= pre_existing_crack;
432 
433  // Adjustments for newly created cracks
434  if (met_stress_criterion && !pre_existing_crack && allowed_to_crack)
435  {
436  new_crack = true;
437  ++num_cracks;
438 
439  // Assume Poisson's ratio drops to zero for this direction. Stiffness is then Young's
440  // modulus.
441  _crack_initiation_strain[_qp](i) = cracking_stress / youngs_modulus;
442 
443  if (_crack_max_strain[_qp](i) < _crack_initiation_strain[_qp](i))
445  }
446 
447  // Update stress and stiffness ratio according to specified crack release model
448  if (new_crack || (pre_existing_crack && loading_existing_crack))
449  {
450  cracked = true;
451  if (_softening_models.size() != 0)
453  stiffness_ratio,
454  strain_in_crack_dir(i),
455  _crack_initiation_strain[_qp](i),
456  _crack_max_strain[_qp](i),
457  cracking_stress,
458  youngs_modulus);
459  else
461  sigma(i),
462  stiffness_ratio,
463  strain_in_crack_dir(i),
464  cracking_stress,
465  cracking_alpha,
466  youngs_modulus);
467  _crack_damage[_qp](i) = 1.0 - stiffness_ratio;
468  }
469 
470  else if (cracked && _cracking_neg_fraction > 0 &&
471  _crack_initiation_strain[_qp](i) * _cracking_neg_fraction > strain_in_crack_dir(i) &&
472  -_crack_initiation_strain[_qp](i) * _cracking_neg_fraction < strain_in_crack_dir(i))
473  {
474  const Real etr = _cracking_neg_fraction * _crack_initiation_strain[_qp](i);
475  const Real Eo = cracking_stress / _crack_initiation_strain[_qp](i);
476  const Real Ec = Eo * (1.0 - _crack_damage_old[_qp](i));
477  const Real a = (Ec - Eo) / (4.0 * etr);
478  const Real b = 0.5 * (Ec + Eo);
479  const Real c = 0.25 * (Ec - Eo) * etr;
480  sigma(i) = (a * strain_in_crack_dir(i) + b) * strain_in_crack_dir(i) + c;
481  }
482  }
483 
484  if (cracked)
486  }
487 
488  _crack_flags[_qp](0) = 1.0 - _crack_damage[_qp](2);
489  _crack_flags[_qp](1) = 1.0 - _crack_damage[_qp](1);
490  _crack_flags[_qp](2) = 1.0 - _crack_damage[_qp](0);
491 }
492 
493 void
495  RealVectorValue & strain_in_crack_dir)
496 {
497  // The rotation tensor is ordered such that directions for pre-existing cracks appear first
498  // in the list of columns. For example, if there is one existing crack, its direction is in the
499  // first column in the rotation tensor.
500  const unsigned int num_known_dirs = getNumKnownCrackDirs();
501 
502  if (num_known_dirs == 0)
503  {
504  std::vector<Real> eigval(3, 0.0);
505  RankTwoTensor eigvec;
506 
507  _elastic_strain[_qp].symmetricEigenvaluesEigenvectors(eigval, eigvec);
508 
509  // If the elastic strain is beyond the cracking strain, save the eigen vectors as
510  // the rotation tensor. Reverse their order so that the third principal strain
511  // (most tensile) will correspond to the first crack.
512  _crack_rotation[_qp].fillColumn(0, eigvec.column(2));
513  _crack_rotation[_qp].fillColumn(1, eigvec.column(1));
514  _crack_rotation[_qp].fillColumn(2, eigvec.column(0));
515 
516  strain_in_crack_dir(0) = eigval[2];
517  strain_in_crack_dir(1) = eigval[1];
518  strain_in_crack_dir(2) = eigval[0];
519  }
520  else if (num_known_dirs == 1)
521  {
522  // This is easily the most complicated case.
523  // 1. Rotate the elastic strain to the orientation associated with the known
524  // crack.
525  // 2. Extract the lower 2x2 block into a separate tensor.
526  // 3. Run the eigen solver on the result.
527  // 4. Update the rotation tensor to reflect the effect of the 2 eigenvectors.
528 
529  // 1.
530  const RankTwoTensor & R = _crack_rotation[_qp];
531  RankTwoTensor ePrime(_elastic_strain[_qp]);
532  ePrime.rotate(R.transpose()); // elastic strain in crack coordinates
533 
534  // 2.
535  ColumnMajorMatrix e2x2(2, 2);
536  e2x2(0, 0) = ePrime(1, 1);
537  e2x2(1, 0) = ePrime(2, 1);
538  e2x2(0, 1) = ePrime(1, 2);
539  e2x2(1, 1) = ePrime(2, 2);
540 
541  // 3.
542  ColumnMajorMatrix e_val2x1(2, 1);
543  ColumnMajorMatrix e_vec2x2(2, 2);
544  e2x2.eigen(e_val2x1, e_vec2x2);
545 
546  // 4.
547  RankTwoTensor eigvec(
548  1.0, 0.0, 0.0, 0.0, e_vec2x2(0, 1), e_vec2x2(1, 1), 0.0, e_vec2x2(0, 0), e_vec2x2(1, 0));
549 
550  _crack_rotation[_qp] = _crack_rotation_old[_qp] * eigvec; // Roe implementation
551 
552  strain_in_crack_dir(0) = ePrime(0, 0);
553  strain_in_crack_dir(1) = e_val2x1(1, 0);
554  strain_in_crack_dir(2) = e_val2x1(0, 0);
555  }
556  else if (num_known_dirs == 2 || num_known_dirs == 3)
557  {
558  // Rotate to cracked orientation and pick off the strains in the rotated
559  // coordinate directions.
560  const RankTwoTensor & R = _crack_rotation[_qp];
561  RankTwoTensor ePrime(_elastic_strain[_qp]);
562  ePrime.rotate(R.transpose()); // elastic strain in crack coordinates
563 
564  strain_in_crack_dir(0) = ePrime(0, 0);
565  strain_in_crack_dir(1) = ePrime(1, 1);
566  strain_in_crack_dir(2) = ePrime(2, 2);
567  }
568  else
569  mooseError("Invalid number of known crack directions");
570 }
571 
572 unsigned int
574 {
575  unsigned int num_known_dirs = 0;
576  for (unsigned int i = 0; i < 3; ++i)
577  {
578  if (_crack_damage_old[_qp](i) > 0.0 || _prescribed_crack_directions.size() >= i + 1)
579  ++num_known_dirs;
580  }
581  return num_known_dirs;
582 }
583 
584 void
586  Real & sigma,
587  Real & stiffness_ratio,
588  const Real strain_in_crack_dir,
589  const Real cracking_stress,
590  const Real cracking_alpha,
591  const Real youngs_modulus)
592 {
593  switch (_cracking_release)
594  {
596  {
597  if (sigma > cracking_stress)
598  {
599  stiffness_ratio /= 3.0;
600  sigma = stiffness_ratio * youngs_modulus * strain_in_crack_dir;
601  }
602  break;
603  }
605  {
606  const Real crack_max_strain = _crack_max_strain[_qp](i);
607  mooseAssert(crack_max_strain >= _crack_initiation_strain[_qp](i),
608  "crack_max_strain must be >= crack_initiation_strain");
609 
610  // Compute stress that follows exponental curve
611  sigma =
612  cracking_stress * (_cracking_residual_stress +
613  (1.0 - _cracking_residual_stress) *
614  std::exp(cracking_alpha * _cracking_beta / cracking_stress *
615  (crack_max_strain - _crack_initiation_strain[_qp](i))));
616  // Compute ratio of current stiffness to original stiffness
617  stiffness_ratio =
618  sigma * _crack_initiation_strain[_qp](i) / (crack_max_strain * cracking_stress);
619  break;
620  }
622  {
623  if (_cracking_residual_stress == 0)
624  {
625  const Real tiny = 1e-16;
626  stiffness_ratio = tiny;
627  sigma = tiny * _crack_initiation_strain[_qp](i) * youngs_modulus;
628  }
629  else
630  {
631  sigma = _cracking_residual_stress * cracking_stress;
632  stiffness_ratio = sigma / (_crack_max_strain[_qp](i) * youngs_modulus);
633  }
634  break;
635  }
636  }
637 
638  if (stiffness_ratio < 0)
639  {
640  std::stringstream err;
641  err << "Negative stiffness ratio: " << i << " " << stiffness_ratio << ", "
642  << _crack_max_strain[_qp](i) << ", " << _crack_initiation_strain[_qp](i) << ", "
643  << std::endl;
644  mooseError(err.str());
645  }
646 }
647 
648 void
650  const RealVectorValue & sigma)
651 {
652  // Get transformation matrix
653  const RankTwoTensor & R = _crack_rotation[_qp];
654  // Rotate to crack frame
655  tensor.rotate(R.transpose());
656 
657  // Reset stress if cracked
658  for (unsigned int i = 0; i < 3; ++i)
659  if (_crack_damage[_qp](i) > 0.0)
660  {
661  const Real stress_correction_ratio = (tensor(i, i) - sigma(i)) / tensor(i, i);
662  if (stress_correction_ratio > _max_stress_correction)
663  tensor(i, i) *= (1.0 - _max_stress_correction);
664  else if (stress_correction_ratio < -_max_stress_correction)
665  tensor(i, i) *= (1.0 + _max_stress_correction);
666  else
667  tensor(i, i) = sigma(i);
668  }
669 
670  // Rotate back to global frame
671  tensor.rotate(R);
672 }
673 
674 bool
676 {
677  for (unsigned int i = 0; i < 3; ++i)
678  if (_crack_damage_old[_qp](i) > 0.0)
679  return true;
680  return false;
681 }
ComputeSmearedCrackingStress computes the stress for a finite strain material with smeared cracking...
MaterialProperty< RankFourTensor > & _Jacobian_mult
derivative of stress w.r.t. strain (_dstress_dstrain)
virtual void computeQpStress() override
Compute the stress and store it in the _stress material property for the current quadrature point...
void updateStressTensorForCracking(RankTwoTensor &tensor, const RealVectorValue &sigma)
Updates the full stress tensor to account for the effect of cracking using the provided stresses in t...
ComputeMultipleInelasticStress computes the stress, the consistent tangent operator (or an approximat...
const MaterialProperty< RankTwoTensor > & _rotation_increment
Rotation increment material property.
MaterialProperty< RealVectorValue > & _crack_max_strain
MaterialProperty< RealVectorValue > & _crack_initiation_strain
const unsigned int _max_cracks
Maximum number of cracks permitted at a material point.
const MaterialProperty< RankTwoTensor > & _elastic_strain_old
Strain tensors.
const MaterialProperty< RankTwoTensor > & _crack_rotation_old
MaterialProperty< Real > & _matl_timestep_limit
const MaterialProperty< RankTwoTensor > & _strain_increment
void updateLocalElasticityTensor()
Update the local elasticity tensor (_local_elasticity_tensor) due to the effects of cracking...
const MaterialProperty< RealVectorValue > & _crack_max_strain_old
MaterialProperty< RankTwoTensor > & _crack_rotation
const MaterialProperty< RankFourTensor > & _elasticity_tensor
Elasticity tensor material property.
bool previouslyCracked()
Check to see whether there was cracking in any diretion in the previous time step.
std::vector< SmearedCrackSofteningBase * > _softening_models
The user-supplied list of softening models to be used in the 3 crack directions.
MaterialProperty< RankTwoTensor > & _stress
Stress material property.
ComputeSmearedCrackingStress(const InputParameters &parameters)
const MaterialProperty< RankTwoTensor > & _inelastic_strain_old
old value of inelastic strain
virtual void computeQpStressIntermediateConfiguration()
Compute the stress for the current QP, but do not rotate tensors from the intermediate configuration ...
const std::string _elasticity_tensor_name
Name of the elasticity tensor material property.
const Real _cracking_beta
Controls slope of exponential softening curve.
void computeCrackStrainAndOrientation(RealVectorValue &strain_in_crack_dir)
Compute the crack strain in the crack coordinate system.
InputParameters validParams< ComputeSmearedCrackingStress >()
enum ComputeSmearedCrackingStress::CrackingRelease _cracking_release
const VariableValue & _cracking_stress
Threshold at which cracking initiates if tensile stress exceeds it.
const Real _max_stress_correction
Controls the maximum amount that the damaged elastic stress is corrected to folow the release model d...
virtual void initQpStatefulProperties() override
MaterialProperty< RealVectorValue > & _crack_flags
Vector of values going from 1 to 0 as crack damage accumulates.
virtual unsigned int getNumKnownCrackDirs() const
Get the number of known crack directions.
SmearedCrackSofteningBase is the base class for a set of models that define the softening behavior of...
Real getIsotropicYoungsModulus(const RankFourTensor &elasticity_tensor)
Get the Young&#39;s modulus for an isotropic elasticity tensor param elasticity_tensor the tensor (must b...
InputParameters validParams< ComputeMultipleInelasticStress >()
const Real _cracking_neg_fraction
Defines transition to changed stiffness during unloading.
MaterialProperty< RealVectorValue > & _crack_damage
const Real _shear_retention_factor
Controls the amount of shear retained.
const MaterialProperty< RealVectorValue > & _crack_damage_old
std::vector< unsigned int > _prescribed_crack_directions
User-prescribed cracking directions.
const MaterialProperty< RealVectorValue > & _crack_initiation_strain_old
virtual void finiteStrainRotation(const bool force_elasticity_rotation=false)
Rotate _elastic_strain, _stress, _inelastic_strain, and _Jacobian_mult to the new configuration...
std::vector< StressUpdateBase * > _models
The user supplied list of inelastic models to use in the simulation.
const bool _perform_finite_strain_rotations
after updateQpState, rotate the stress, elastic_strain, inelastic_strain and Jacobian_mult using _rot...
virtual void computeCrackingRelease(int i, Real &sigma, Real &stiffness_ratio, const Real strain_in_crack_dir, const Real cracking_stress, const Real cracking_alpha, const Real youngs_modulus)
Compute the effect of the cracking release model on the stress and stiffness in the direction of a si...
MaterialProperty< RankTwoTensor > & _elastic_strain
Elastic strain material property.
registerMooseObject("TensorMechanicsApp", ComputeSmearedCrackingStress)
bool hasGuaranteedMaterialProperty(const MaterialPropertyName &prop, Guarantee guarantee)
virtual void updateCrackingStateAndStress()
Update all cracking-related state variables and the stress tensor due to cracking in all directions...
const Real _cracking_residual_stress
Input parameters for smeared crack models.
MaterialProperty< RankTwoTensor > & _inelastic_strain
The sum of the inelastic strains that come from the plastic models.
CrackingRelease
Enum defining the crack release model.