www.mooseframework.org
ReturnMappingModel.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 #include "ReturnMappingModel.h"
8 
10 #include "Conversion.h"
11 
12 template <>
13 InputParameters
15 {
16  InputParameters params = validParams<ConstitutiveModel>();
18  params.addParam<Real>("max_inelastic_increment",
19  1e-4,
20  "The maximum inelastic strain increment allowed in a time step");
21  return params;
22 }
23 
24 ReturnMappingModel::ReturnMappingModel(const InputParameters & parameters,
25  const std::string inelastic_strain_name)
26  : ConstitutiveModel(parameters),
28  _effective_strain_increment(0),
29  _effective_inelastic_strain(
30  declareProperty<Real>("effective_" + inelastic_strain_name + "_strain")),
31  _effective_inelastic_strain_old(
32  getMaterialPropertyOld<Real>("effective_" + inelastic_strain_name + "_strain")),
33  _max_inelastic_increment(parameters.get<Real>("max_inelastic_increment"))
34 {
35 }
36 
37 void
39 {
40  _effective_inelastic_strain[_qp] = 0.0;
41 }
42 
43 void
44 ReturnMappingModel::computeStress(const Elem & current_elem,
45  const SymmElasticityTensor & elasticityTensor,
46  const SymmTensor & stress_old,
47  SymmTensor & strain_increment,
48  SymmTensor & stress_new)
49 {
50  // Given the stretching, compute the stress increment and add it to the old stress. Also update
51  // the creep strain
52  // stress = stressOld + stressIncrement
53  if (_t_step == 0 && !_app.isRestarting())
54  return;
55 
56  stress_new = elasticityTensor * strain_increment;
57  stress_new += stress_old;
58 
59  SymmTensor inelastic_strain_increment;
60  computeStress(current_elem,
61  elasticityTensor,
62  stress_old,
63  strain_increment,
64  stress_new,
65  inelastic_strain_increment);
66 }
67 
68 void
69 ReturnMappingModel::computeStress(const Elem & /*current_elem*/,
70  const SymmElasticityTensor & elasticityTensor,
71  const SymmTensor & stress_old,
72  SymmTensor & strain_increment,
73  SymmTensor & stress_new,
74  SymmTensor & inelastic_strain_increment)
75 {
76  // compute deviatoric trial stress
77  SymmTensor dev_trial_stress(stress_new);
78  dev_trial_stress.addDiag(-dev_trial_stress.trace() / 3.0);
79 
80  // compute effective trial stress
81  Real dts_squared = dev_trial_stress.doubleContraction(dev_trial_stress);
82  Real effective_trial_stress = std::sqrt(1.5 * dts_squared);
83 
84  // compute effective strain increment
85  SymmTensor dev_strain_increment(strain_increment);
86  dev_strain_increment.addDiag(-strain_increment.trace() / 3.0);
87  _effective_strain_increment = dev_strain_increment.doubleContraction(dev_strain_increment);
89 
90  const SymmIsotropicElasticityTensor * iso_e_t =
91  dynamic_cast<const SymmIsotropicElasticityTensor *>(&elasticityTensor);
92  if (!iso_e_t)
93  mooseError("Models derived from ReturnMappingModel require a SymmIsotropicElasticityTensor");
94  _three_shear_modulus = 3.0 * iso_e_t->shearModulus();
95 
96  computeStressInitialize(effective_trial_stress, elasticityTensor);
97 
98  Real scalar;
99  returnMappingSolve(effective_trial_stress, scalar, _console);
100 
101  // compute inelastic and elastic strain increments
103  {
104  if (effective_trial_stress < 0.01)
105  effective_trial_stress = 0.01;
106 
107  inelastic_strain_increment = dev_trial_stress;
108  inelastic_strain_increment *= (1.5 * scalar / effective_trial_stress);
109  }
110  else
111  {
112  if (scalar != 0.0)
113  inelastic_strain_increment = dev_trial_stress * (1.5 * scalar / effective_trial_stress);
114  else
115  inelastic_strain_increment = 0.0;
116  }
117 
118  strain_increment -= inelastic_strain_increment;
120 
121  // compute stress increment
122  stress_new = elasticityTensor * strain_increment;
123 
124  // update stress
125  stress_new += stress_old;
126 
127  computeStressFinalize(inelastic_strain_increment);
128 }
129 
130 Real
131 ReturnMappingModel::computeReferenceResidual(const Real effective_trial_stress, const Real scalar)
132 {
133  return effective_trial_stress / _three_shear_modulus - scalar;
134 }
135 
136 Real
138 {
139  Real scalar_inelastic_strain_incr;
140 
141  scalar_inelastic_strain_incr =
143  if (MooseUtils::absoluteFuzzyEqual(scalar_inelastic_strain_incr, 0.0))
144  return std::numeric_limits<Real>::max();
145 
146  return _dt * _max_inelastic_increment / scalar_inelastic_strain_incr;
147 }
This class defines a basic set of capabilities any elasticity tensor should have. ...
MaterialProperty< Real > & _effective_inelastic_strain
bool _legacy_return_mapping
Whether to use the legacy return mapping algorithm and compute residuals in the legacy manner...
virtual void initQpStatefulProperties() override
Base class that provides capability for Newton return mapping iterations on a single variable...
virtual void computeStress(const Elem &current_elem, const SymmElasticityTensor &elasticityTensor, const SymmTensor &stress_old, SymmTensor &strain_increment, SymmTensor &stress_new) override
Real _three_shear_modulus
3 * shear modulus
void returnMappingSolve(const Real effective_trial_stress, Real &scalar, const ConsoleStream &console)
Perform the return mapping iterations.
virtual Real computeReferenceResidual(const Real effective_trial_stress, const Real scalar) override
Compute a reference quantity to be used for checking relative convergence.
Real trace() const
Definition: SymmTensor.h:95
virtual void computeStressInitialize(Real, const SymmElasticityTensor &)
Perform any necessary initialization before return mapping iterations.
ReturnMappingModel(const InputParameters &parameters, const std::string inelastic_strain_name="")
const MaterialProperty< Real > & _effective_inelastic_strain_old
virtual void computeStressFinalize(const SymmTensor &)
Perform any necessary steps to finalize state after return mapping iterations.
Real computeTimeStepLimit()
Compute the limiting value of the time step for this material.
Real shearModulus() const
Return the shear modulus...
void addDiag(Real value)
Definition: SymmTensor.h:279
InputParameters validParams< ConstitutiveModel >()
Real doubleContraction(const SymmTensor &rhs) const
Definition: SymmTensor.h:257
InputParameters validParams< SingleVariableReturnMappingSolution >()
Defines an Isotropic Elasticity Tensor.
InputParameters validParams< ReturnMappingModel >()